v0.15.0
Loading...
Searching...
No Matches
SurfaceSlidingConstrains.hpp
Go to the documentation of this file.
1/** \file SurfaceSlidingConstrains.hpp
2 * \brief Implementing surface sliding constrains
3 */
4
5
6
7#ifndef __SURFACE_SLIDING_CONSTRAINS_HPP__
8#define __SURFACE_SLIDING_CONSTRAINS_HPP__
9
10#ifdef WITH_ADOL_C
11
12/**
13 * @brief Implementation of surface sliding constrains
14 *
15 */
17
18 GenericSliding() = default;
19 virtual ~GenericSliding() = default;
20
23 boost::shared_ptr<VectorDouble> activeVariablesPtr;
25 boost::shared_ptr<VectorDouble> &active_variables_ptr)
28 activeVariablesPtr(active_variables_ptr) {}
29 MoFEMErrorCode doWork(int side, EntityType type,
30 EntitiesFieldData::EntData &data) {
32 if (type == MBVERTEX) {
33 for (unsigned int dd = 0; dd != data.getFieldData().size(); ++dd)
34 (*activeVariablesPtr)[dd] = data.getFieldData()[dd];
35 }
37 }
38 };
39
40 template <int SizeLambda>
41 struct OpGetActiveDofsPositions
43 boost::shared_ptr<VectorDouble> activeVariablesPtr;
45 const std::string field_name,
46 boost::shared_ptr<VectorDouble> &active_variables_ptr)
49 activeVariablesPtr(active_variables_ptr) {}
50 MoFEMErrorCode doWork(int side, EntityType type,
51 EntitiesFieldData::EntData &data) {
53 if (type == MBVERTEX) {
54 for (unsigned int dd = 0; dd != data.getFieldData().size(); ++dd)
55 (*activeVariablesPtr)[SizeLambda + dd] = data.getFieldData()[dd];
56 }
58 }
59 };
60
61 template <int SizeLambda, int SizePositions>
62 struct OpAssembleRhs : public MoFEM::ForcesAndSourcesCore::UserDataOperator {
63
64 boost::shared_ptr<VectorDouble> resultsPtr;
65
66 OpAssembleRhs(const std::string field_name,
67 boost::shared_ptr<VectorDouble> &results_ptr)
70 resultsPtr(results_ptr) {}
71
72 MoFEMErrorCode doWork(int side, EntityType type,
73 EntitiesFieldData::EntData &data) {
75 if (type != MBVERTEX)
77 VectorInt &indices = data.getIndices();
78 int shift = 0;
79 if (indices.empty())
81 else if (indices.size() == SizeLambda)
82 shift = 0;
83 else if (indices.size() == SizePositions)
84 shift = SizeLambda;
85 else
86 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
87 "Element %s: Data inconsistency nb of indices %d",
88 getFEName().c_str(), indices.size());
89
90 CHKERR VecSetOption(getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
91 CHKERR VecSetValues(getSNESf(), indices.size(), &indices[0],
92 &(*resultsPtr)[shift], ADD_VALUES);
94 }
95 };
96 template <int SizeLambda, int SizePositions>
97 struct OpAssembleLhs : public MoFEM::ForcesAndSourcesCore::UserDataOperator {
98
99 boost::shared_ptr<MatrixDouble> jacobianPtr;
100
101 OpAssembleLhs(const std::string field_name_row,
102 const std::string field_name_col,
103 boost::shared_ptr<MatrixDouble> &jacobian_ptr)
105 field_name_row, field_name_col, UserDataOperator::OPROWCOL),
106 jacobianPtr(jacobian_ptr) {
107 sYmm = false;
108 }
109
110 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
111 EntityType col_type,
112 EntitiesFieldData::EntData &row_data,
113 EntitiesFieldData::EntData &col_data) {
115 if (row_type != MBVERTEX)
117 if (col_type != MBVERTEX)
119 VectorInt &row_indices = row_data.getIndices();
120 VectorInt &col_indices = col_data.getIndices();
121 if (row_indices.empty() || col_indices.empty())
123
124 int shift_row = 0;
125 if (row_indices.size() == SizeLambda)
126 shift_row = 0;
127 else if (row_indices.size() == SizePositions)
128 shift_row = SizeLambda;
129 else
130 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
131 "Data inconsistency");
132
133 int shift_col = 0;
134 if (col_indices.size() == SizeLambda)
135 shift_col = 0;
136 else if (col_indices.size() == SizePositions)
137 shift_col = SizeLambda;
138 else
139 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
140 "Data inconsistency");
141
142 MatrixDouble jac(row_indices.size(), col_indices.size());
143 for (unsigned int rr = 0; rr != row_indices.size(); ++rr) {
144 for (unsigned int cc = 0; cc != col_indices.size(); ++cc) {
145 jac(rr, cc) = (*jacobianPtr)(shift_row + rr, shift_col + cc);
146 if (jac(rr, cc) != jac(rr, cc))
147 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
148 "Jacobian assemble not a number jac(rr,cc) = %3.4f",
149 jac(rr, cc));
150 }
151 }
152 CHKERR MatSetValues(getSNESB(), row_indices.size(), &row_indices[0],
153 col_indices.size(), &col_indices[0], &jac(0, 0),
154 ADD_VALUES);
156 }
157 };
158};
159
160/** \brief Shape preserving constrains, i.e. nodes sliding on body surface.
161
162 Derivation and implementation of constrains preserving body surface,
163 i.e. body shape and volume.
164
165 The idea starts form observation that body shape can be globally characterized
166 by constant calculated as volume over its area
167 \f[
168 \frac{V}{A} = C
169 \f]
170 Above equation expressed in integral form is
171 \f[
172 \int_\Omega \textrm{d}V = C \int_\Gamma \textrm{d}S
173 \f]
174 where notting that,
175 \f[
176 \frac{1}{3}
177 \int_\Omega \textrm{div}[\mathbf{X}] \textrm{d}\Omega
178 =
179 C \int_\Gamma \textrm{d}S
180 \f]
181 and applying Gauss theorem we get
182 \f[
183 \int_\Gamma
184 \mathbf{X}\cdot \frac{\mathbf{N}}{\|\mathbf{N}\|}
185 \textrm{d}\Gamma
186 =
187 3C \int_\Gamma \textrm{d}S.
188 \f]
189 Drooping integrals on both sides, and linearizing equation, we get
190 \f[
191 \frac{\mathbf{N}}{\|\mathbf{N}\|} \cdot \delta \mathbf{X}
192 =
193 3C - \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot \mathbf{X}
194 \f]
195 where \f$\delta \mathbf{X}\f$ is displacement sub-inctrement. Above equation
196 is a constrain if satisfied in body shape and volume is conserved. Final form
197 of constrain equation is \f[ \mathcal{r} =
198 \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot \mathbf{X}
199 -
200 \frac{\mathbf{N_0}}{\|\mathbf{N_0}\|}\cdot \mathbf{X}_0 =
201 \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot (\mathbf{X}-\mathbf{X}_0)
202 \f]
203
204 In the spirit of finite element method the constrain equation is multiplied
205 by shape functions and enforce using Lagrange multiplier method
206 \f[
207 \int_\Gamma \mathbf{N}^\mathsf{T}_\lambda
208 \left(
209 \frac{\mathbf{N}}{\|\mathbf{N}\|}\mathbf{N}_\mathbf{X}\cdot
210 (\overline{\mathbf{X}}-\overline{\mathbf{X}}_0)
211 \right)
212 \|\mathbf{N}\|
213 \textrm{d}\Gamma
214 =
215 \mathbf{0}.
216 \f]
217 Above equation is nonlinear, applying to it Taylor expansion, we can get form
218 which can be used with Newton interactive method \f[ \begin{split}
219 &\int_\Gamma \mathbf{N}^\mathsf{T}_\lambda
220 \left\{
221 \mathbf{N}\mathbf{N}_\mathbf{X}
222 +
223 \left(\mathbf{X}-\mathbf{X}_0\right) \cdot
224 \left(
225 \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
226 -
227 \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
228 \right)
229 \right\}
230 \textrm{d}\Gamma
231 \cdot
232 \delta
233 \overline{\mathbf{X}}\\
234 =
235 &\int_\Gamma \mathbf{N}^\mathsf{T}_\lambda
236 \mathbf{N}\cdot(\mathbf{X}-\mathbf{X}_0)
237 \textrm{d}\Gamma
238 \end{split}.
239 \f]
240 Equation expressing forces on shape as result of constrains, as result
241 Lagrange multiplier method have following form \f[ \begin{split}
242 &\int_\Gamma
243 \mathbf{N}^\mathsf{T}_\mathbf{X} \cdot \mathbf{N}
244 \mathbf{N}_\lambda
245 \textrm{d}\Gamma
246 \cdot
247 \delta\overline{\lambda}\\
248 +
249 &\int_\Gamma
250 \lambda
251 \mathbf{N}^\mathsf{T}_\mathbf{X}
252 \left(
253 \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
254 -
255 \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
256 \right)
257 \textrm{d}\Gamma
258 \delta\overline{\mathbf{X}}\\
259 =
260 &\int_\Gamma
261 \lambda
262 \mathbf{N}^\mathsf{T}_\mathbf{X} \cdot \mathbf{N}
263 \textrm{d}\Gamma
264 \end{split}
265 \f]
266
267 Above equations are assembled into global system of equations as following
268 \f[
269 \left[
270 \begin{array}{cc}
271 \mathbf{K} + \mathbf{B} & \mathbf{C}^\mathsf{T} \\
272 \mathbf{C} + \mathbf{A} & 0
273 \end{array}
274 \right]
275 \left\{
276 \begin{array}{c}
277 \delta \overline{\mathbf{X}} \\
278 \delta \overline{\lambda}
279 \end{array}
280 \right\}=
281 \left[
282 \begin{array}{c}
283 \mathbf{f} - \mathbf{C}^\mathsf{T}\overline{\lambda} \\
284 \overline{\mathbf{r}}
285 \end{array}
286 \right]
287 \f]
288 where
289 \f[
290 \mathbf{C}=
291 \int_\Gamma
292 \mathbf{N}_\lambda^\mathsf{T}
293 \mathbf{N} \cdot
294 \mathbf{N}_\mathbf{X}
295 \textrm{d}\Gamma,
296 \f]
297 \f[
298 \mathbf{B}=
299 \int_\Gamma
300 \lambda
301 (\mathbf{X}-\mathbf{X}_0)\cdot
302 \left(
303 \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
304 -
305 \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
306 \right)
307 \textrm{d}\Gamma
308 \f]
309 and
310 \f[
311 \mathbf{A}=
312 \int_\Gamma
313 \mathbf{N}^\mathsf{T}_\lambda
314 \left(\mathbf{X}-\mathbf{X}_0\right) \cdot
315 \left(
316 \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
317 -
318 \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
319 \right).
320 \f]
321
322*/
324
325 /** \brief Class implemented by user to detect face orientation
326
327 If mesh generated is with surface mesher, usually you don't have to do
328 nothing, all elements on the surface have consistent orientation. In case of
329 internal faces or if you do something with mesh connectivity which breaks
330 orientation on the face, you have to implement method which will set
331 orientation to face.
332
333 */
335
337
338 virtual MoFEMErrorCode
345 virtual MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field,
346 const EntityHandle face) {
350 }
351 };
352
354
356
357 Mat B;
358 Vec F;
359
361 : MoFEM::FaceElementForcesAndSourcesCore(m_field), B(PETSC_NULLPTR),
362 F(PETSC_NULLPTR) {}
363 int getRule(int order) { return 2 * order; };
364
365 MoFEMErrorCode preProcess() {
367
369
370 if (B != PETSC_NULLPTR) {
371 snes_B = B;
372 }
373
374 if (F != PETSC_NULLPTR) {
375 snes_f = F;
376 }
377
379 }
380 };
381
382 boost::shared_ptr<MyTriangleFE> feRhsPtr, feLhsPtr;
383
388
390
391 double aLpha;
392 MoFEMErrorCode getOptions() {
394 PetscOptionsBegin(PETSC_COMM_WORLD, "",
395 "Get surface sliding constrains element scaling",
396 "none");
397 CHKERR PetscOptionsScalar("-surface_sliding_alpha", "scaling parameter", "",
398 aLpha, &aLpha, PETSC_NULLPTR);
399 PetscOptionsEnd();
401 }
402
404 DriverElementOrientation &orientation)
405 : mField(m_field), feRhsPtr(new MyTriangleFE(m_field)),
406 feLhsPtr(new MyTriangleFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr),
407 crackFrontOrientation(orientation), aLpha(1) {
408
409 ierr = getOptions();
410 CHKERRABORT(PETSC_COMM_WORLD, ierr);
411 }
412
413 virtual ~SurfaceSlidingConstrains() = default;
414
415 struct OpJacobian
417
418 const int tAg;
419 boost::shared_ptr<VectorDouble> activeVariablesPtr;
420 boost::shared_ptr<VectorDouble> resultsPtr;
421 boost::shared_ptr<MatrixDouble> jacobianPtr;
424
425 const double &aLpha;
426
427 OpJacobian(int tag, const std::string field_name,
428 boost::shared_ptr<VectorDouble> &active_variables_ptr,
429 boost::shared_ptr<VectorDouble> &results_ptr,
430 boost::shared_ptr<MatrixDouble> &jacobian_ptr,
431 DriverElementOrientation &orientation, bool evaluate_jacobian,
432 const double &alpha)
435 tAg(tag), activeVariablesPtr(active_variables_ptr),
436 resultsPtr(results_ptr), jacobianPtr(jacobian_ptr),
437 oRientation(orientation), evaluateJacobian(evaluate_jacobian),
438 aLpha(alpha) {}
439
440 MoFEMErrorCode doWork(int side, EntityType type,
441 EntitiesFieldData::EntData &data) {
443 if (type != MBVERTEX)
445
446 VectorInt &indices = data.getIndices();
447
448 trace_on(tAg);
449
450 ublas::vector<adouble> lambda_dofs(3);
451 for (int dd = 0; dd != 3; ++dd) {
452 lambda_dofs[dd] <<= (*activeVariablesPtr)[dd];
453 }
454 ublas::vector<adouble> position_dofs(9);
455 for (int dd = 0; dd != 9; ++dd) {
456 position_dofs[dd] <<= (*activeVariablesPtr)[3 + dd];
457 }
458
459 FTensor::Index<'i', 3> i;
460 FTensor::Index<'j', 3> j;
461 FTensor::Index<'k', 3> k;
465
467 getFEMethod());
469
470 int nb_gauss_pts = data.getN().size1();
471 int nb_base_functions = data.getN().size2();
472 auto t_base1 = data.getFTensor0N();
473 auto t_base2 = data.getFTensor0N();
474 auto t_diff_base = data.getFTensor1DiffN<2>();
476 FTensor::Tensor1<adouble, 3> t_position_ksi;
477 FTensor::Tensor1<adouble, 3> t_position_eta;
481
482 ublas::vector<adouble> c_vec(3);
483 ublas::vector<adouble> f_vec(9);
484 c_vec.clear();
485 f_vec.clear();
486
487 auto t_coord_ref = getFTensor1CoordsAtGaussPts();
488
489 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
490
491 FTensor::Tensor1<adouble *, 3> t_position_dofs(
492 &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
493 FTensor::Tensor0<adouble *> t_lambda_dof(&lambda_dofs[0]);
494
495 t_position(i) = 0;
496 t_position_ksi(i) = 0;
497 t_position_eta(i) = 0;
498 lambda = 0;
499
500 for (int bb = 0; bb != nb_base_functions; ++bb) {
501 t_position(i) += t_base1 * t_position_dofs(i);
502 t_position_ksi(i) += t_diff_base(N0) * t_position_dofs(i);
503 t_position_eta(i) += t_diff_base(N1) * t_position_dofs(i);
504 lambda += t_base1 * t_lambda_dof;
505 ++t_base1;
506 ++t_position_dofs;
507 ++t_lambda_dof;
508 ++t_diff_base;
509 }
510
511 t_delta(i) = t_position(i) - t_coord_ref(i);
512 t_normal(k) = FTensor::cross(t_position_ksi(i), t_position_eta(j), k);
513
514 double w = getGaussPts()(2, gg) * 0.5 * aLpha;
515 adouble val;
516 FTensor::Tensor0<adouble *> t_c(&c_vec[0]);
517 FTensor::Tensor1<adouble *, 3> t_f(&f_vec[0], &f_vec[1], &f_vec[2], 3);
518
519 adouble area = sqrt(t_normal(i) * t_normal(i));
520
521 for (int bb = 0; bb != nb_base_functions; ++bb) {
522 if (indices[bb] != -1) {
523 val = w * eo * t_base2;
524 t_c += val * t_normal(i) * t_delta(i);
525 val *= lambda;
526 t_f(i) += val * t_normal(i);
527 }
528 ++t_c;
529 ++t_f;
530 ++t_base2;
531 }
532
533 ++t_coord_ref;
534 }
535
536 for (int rr = 0; rr != 3; ++rr)
537 c_vec[rr] >>= (*resultsPtr)[rr];
538
539 for (int rr = 0; rr != 9; ++rr)
540 f_vec(rr) >>= (*resultsPtr)[3 + rr];
541
542 trace_off();
543
544 if (evaluateJacobian) {
545 double *jac_ptr[3 + 9];
546 for (int rr = 0; rr != 3 + 9; ++rr)
547 jac_ptr[rr] = &(*jacobianPtr)(rr, 0);
548
549 // play recorder for jacobians
550 int r =
551 ::jacobian(tAg, 3 + 9, 3 + 9, &(*activeVariablesPtr)[0], jac_ptr);
552 if (r < 0)
553 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
554 "ADOL-C function evaluation with error");
555 }
556
558 }
559 };
560
561 MoFEMErrorCode setOperators(int tag,
562 const std::string lagrange_multipliers_field_name,
563 const std::string material_field_name,
564 const double *alpha = nullptr) {
566
567 if (alpha != nullptr) {
568 aLpha = *alpha;
569 }
570
571 boost::shared_ptr<VectorDouble> active_variables_ptr(
572 new VectorDouble(3 + 9));
573 boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
574 boost::shared_ptr<MatrixDouble> jacobian_ptr(
575 new MatrixDouble(3 + 9, 3 + 9));
576
577 feRhs.getOpPtrVector().clear();
579 lagrange_multipliers_field_name, active_variables_ptr));
581 material_field_name, active_variables_ptr));
582 feRhs.getOpPtrVector().push_back(new OpJacobian(
583 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
584 jacobian_ptr, crackFrontOrientation, false, aLpha));
585 feRhs.getOpPtrVector().push_back(
586 new OpAssembleRhs<3, 9>(lagrange_multipliers_field_name, results_ptr));
587 feRhs.getOpPtrVector().push_back(
588 new OpAssembleRhs<3, 9>(material_field_name, results_ptr));
589
590 // Adding operators to calculate the left hand side
591 feLhs.getOpPtrVector().clear();
593 lagrange_multipliers_field_name, active_variables_ptr));
595 material_field_name, active_variables_ptr));
596 feLhs.getOpPtrVector().push_back(new OpJacobian(
597 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
598 jacobian_ptr, crackFrontOrientation, true, aLpha));
600 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
602 material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
604 material_field_name, material_field_name, jacobian_ptr));
605
607 }
608
609 MoFEMErrorCode
611 const std::string lagrange_multipliers_field_name,
612 const std::string material_field_name) {
614
615 boost::shared_ptr<VectorDouble> active_variables_ptr(
616 new VectorDouble(3 + 9));
617 boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
618 boost::shared_ptr<MatrixDouble> jacobian_ptr(
619 new MatrixDouble(3 + 9, 3 + 9));
620
621 // Adding operators to calculate the left hand side
622 feLhs.getOpPtrVector().clear();
624 lagrange_multipliers_field_name, active_variables_ptr));
626 material_field_name, active_variables_ptr));
627 feLhs.getOpPtrVector().push_back(new OpJacobian(
628 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
629 jacobian_ptr, crackFrontOrientation, true, aLpha));
631 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
632
634 }
635};
636
638
639 struct CalculateEdgeBase {
640
641 static MoFEMErrorCode createTag(moab::Interface &moab, Tag &th0, Tag &th1,
642 Tag &th2, Tag &th3) {
644 double def_val[] = {0, 0, 0};
645 CHKERR moab.tag_get_handle("EDGE_BASE0", 3, MB_TYPE_DOUBLE, th0,
646 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
647 CHKERR moab.tag_get_handle("EDGE_BASE1", 3, MB_TYPE_DOUBLE, th1,
648 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
649 int def_numb_val[] = {-1};
650 CHKERR moab.tag_get_handle("PATCH_NUMBER", 1, MB_TYPE_INTEGER, th2,
651 MB_TAG_CREAT | MB_TAG_SPARSE, def_numb_val);
652 int def_orientation_val[] = {1};
653 CHKERR moab.tag_get_handle("PATCH_ORIENTATION", 1, MB_TYPE_INTEGER, th3,
654 MB_TAG_CREAT | MB_TAG_SPARSE,
655 def_orientation_val);
656
658 }
659
660 static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges,
661 Range tris) {
663
664 auto get_edges = [&](const Range &ents) {
665 Range edges;
666 CHKERR moab.get_adjacencies(ents, 1, false, edges,
667 moab::Interface::UNION);
668 return edges;
669 };
670
671 auto get_face_adj = [&, get_edges](const Range &faces) {
672 Range adj_faces;
673 CHKERR moab.get_adjacencies(subtract(get_edges(faces), edges), 2, false,
674 adj_faces, moab::Interface::UNION);
675 return intersect(adj_faces, tris);
676 };
677
678 auto get_patch = [&, get_face_adj](const EntityHandle face) {
679 Range patch_ents;
680 patch_ents.insert(face);
681 unsigned int nb0;
682 do {
683 nb0 = patch_ents.size();
684 patch_ents.merge(get_face_adj(patch_ents));
685 } while (nb0 != patch_ents.size());
686 return patch_ents;
687 };
688
689 auto get_patches = [&]() {
690 std::vector<Range> patches;
691 while (!tris.empty()) {
692 patches.push_back(get_patch(tris[0]));
693 tris = subtract(tris, patches.back());
694 }
695 return patches;
696 };
697
698 Tag th0, th1, th2, th3;
699 CHKERR createTag(moab, th0, th1, th2, th3);
700
701 auto patches = get_patches();
702 int pp = 0;
703 for (auto patch : patches) {
704 // cerr << "pp: " << pp << endl;
705 // cerr << patch << endl;
706 std::vector<int> tags_vals(patch.size(), pp);
707 CHKERR moab.tag_set_data(th2, patch, &*tags_vals.begin());
708 ++pp;
709 }
710
712 }
713
714 static MoFEMErrorCode setTags(
715 moab::Interface &moab, Range edges, Range tris,
716 bool number_pathes = true,
717 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
718 surface_orientation = nullptr,
719 MoFEM::Interface *m_field_ptr = nullptr) {
721 Tag th0, th1, th2, th3;
722 CHKERR createTag(moab, th0, th1, th2, th3);
723 if (number_pathes) {
724 CHKERR numberSurfaces(moab, edges, tris);
725 }
726 for (auto edge : edges) {
727 Range adj_faces;
728 CHKERR moab.get_adjacencies(&edge, 1, 2, false, adj_faces);
729 adj_faces = intersect(adj_faces, tris);
730 if (adj_faces.size() != 2) {
731 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
732 "Should be 2 faces adjacent to edge but is %d",
733 adj_faces.size());
734 }
735 VectorDouble3 v[2] = {VectorDouble3(3), VectorDouble3(3)};
736 auto get_tensor_from_vec = [](VectorDouble3 &v) {
738 &v[0], &v[1], &v[2]);
739 };
740
741 FTensor::Index<'i', 3> i;
742 FTensor::Index<'j', 3> j;
743 FTensor::Index<'k', 3> k;
744
745 std::array<double, 3> areas;
746 auto calculate_normals = [&, get_tensor_from_vec]() {
747 int ff = 0;
748 for (auto face : adj_faces) {
749 double &x = (v[ff])[0];
750 double &y = (v[ff])[1];
751 double &z = (v[ff])[2];
752 moab::Util::normal(&moab, face, x, y, z);
753 auto t_n = get_tensor_from_vec(v[ff]);
754 areas[ff] = sqrt(t_n(i) * t_n(i));
755 t_n(i) /= areas[ff];
756 int orientation;
757 // FIXME: this tag is empty
758 CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
759 if (orientation == -1) {
760 t_n(i) *= -1;
761 }
762 if (surface_orientation && m_field_ptr) {
763 CHKERR surface_orientation->getElementOrientation(*m_field_ptr,
764 face);
765 int eo = surface_orientation->elementOrientation;
766 if (eo == -1) {
767 t_n(i) *= -1;
768 }
769 }
770 ++ff;
771 }
772 };
773 calculate_normals();
774
775 auto get_patch_number = [&]() {
776 std::vector<int> p = {0, 0};
777 CHKERR moab.tag_get_data(th2, adj_faces, &*p.begin());
778 return p;
779 };
780
781 std::array<EntityHandle, 2> faces_handles = {adj_faces[0],
782 adj_faces[1]};
783 auto order_normals = [&, get_patch_number]() {
784 auto p = get_patch_number();
785 if (p[0] < p[1]) {
786 v[0].swap(v[1]);
787 faces_handles[0] = adj_faces[1];
788 faces_handles[1] = adj_faces[0];
789 }
790 };
791 order_normals();
792
794 auto t_n0 = get_tensor_from_vec(v[0]);
795 auto t_n1 = get_tensor_from_vec(v[1]);
796 t_cross(k) = FTensor::cross(t_n0(i), t_n1(j), k);
797
798 auto get_base_for_coplanar_faces = [&]() {
800
801 std::vector<EntityHandle> face_conn;
802 CHKERR moab.get_connectivity(&faces_handles[1], 1, face_conn, true);
803 std::vector<EntityHandle> edge_conn;
804 CHKERR moab.get_connectivity(&edge, 1, edge_conn, true);
805 std::sort(face_conn.begin(), face_conn.end());
806 std::sort(edge_conn.begin(), edge_conn.end());
807 int n = 0;
808 for (; n != 2; ++n)
809 if (face_conn[n] != edge_conn[n])
810 break;
811 VectorDouble3 coords_face(3);
812 CHKERR moab.get_coords(&face_conn[n], 1,
813 &*coords_face.data().begin());
814 VectorDouble6 coords_edge(6);
815 CHKERR moab.get_coords(&edge_conn[0], 2,
816 &*coords_edge.data().begin());
817 auto t_edge_n0 = FTensor::Tensor1<double, 3>(
818 coords_edge[0], coords_edge[1], coords_edge[2]);
819 auto t_edge_n1 = FTensor::Tensor1<double, 3>(
820 coords_edge[3], coords_edge[4], coords_edge[5]);
821 auto t_face_n = get_tensor_from_vec(coords_face);
822 t_face_n(i) -= t_edge_n0(i);
824 t_delta(i) = t_edge_n1(i) - t_edge_n0(i);
825 t_delta(i) /= sqrt(t_delta(i) * t_delta(i));
826 t_n1(k) = FTensor::cross(t_n0(i), t_delta(j), k);
827 if (t_n1(i) * t_face_n(i) < 0)
828 t_n1(i) *= -1;
829
831 };
832
833 auto get_base_for_not_planar_faces = [&]() {
835 t_n1(k) = FTensor::cross(t_n0(i), t_cross(j), k);
836 t_n1(i) /= sqrt(t_n1(i) * t_n1(i));
838 };
839
840 // This a case when faces adjacent to edge are coplanar !!!
841 constexpr double tol = 1e-6;
842 if ((t_cross(i) * t_cross(i)) < tol)
843 CHKERR get_base_for_coplanar_faces();
844 else
845 CHKERR get_base_for_not_planar_faces();
846
847 VectorDouble3 &v0 = v[0];
848 VectorDouble3 &v1 = v[1];
849 CHKERR moab.tag_set_data(th0, &edge, 1, &v0[0]);
850 CHKERR moab.tag_set_data(th1, &edge, 1, &v1[0]);
851 }
853 }
854
855 static MoFEMErrorCode saveEdges(moab::Interface &moab, std::string name,
856 Range edges, Range *faces = nullptr) {
858 EntityHandle meshset;
859 Tag ths[4];
860 CHKERR createTag(moab, ths[0], ths[1], ths[2], ths[3]);
861 CHKERR moab.create_meshset(MESHSET_SET, meshset);
862 CHKERR moab.add_entities(meshset, edges);
863 if (faces != nullptr) {
864 CHKERR moab.add_entities(meshset, *faces);
865 }
866 CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1, ths, 4);
867 CHKERR moab.delete_entities(&meshset, 1);
869 }
870 };
871
873
875
876 Mat B;
877 Vec F;
878
880 : MoFEM::EdgeElementForcesAndSourcesCore(m_field), B(PETSC_NULLPTR),
881 F(PETSC_NULLPTR) {}
882 int getRule(int order) { return 2 * order; };
883
884 MoFEMErrorCode preProcess() {
886
888
889 if (B != PETSC_NULLPTR) {
890 snes_B = B;
891 }
892
893 if (F != PETSC_NULLPTR) {
894 snes_f = F;
895 }
896
898 }
899 };
900
901 boost::shared_ptr<MyEdgeFE> feRhsPtr, feLhsPtr;
902
907
908 double aLpha;
909
910 MoFEMErrorCode getOptions() {
912 PetscOptionsBegin(PETSC_COMM_WORLD, "",
913 "Get edge sliding constrains element scaling",
914 "none");
915 CHKERR PetscOptionsScalar("-edge_sliding_alpha", "scaling parameter", "",
916 aLpha, &aLpha, PETSC_NULLPTR);
917 PetscOptionsEnd();
919 }
920
922 : mField(m_field), feRhsPtr(new MyEdgeFE(m_field)),
923 feLhsPtr(new MyEdgeFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr),
924 aLpha(1) {
925 ierr = getOptions();
926 CHKERRABORT(PETSC_COMM_WORLD, ierr);
927 }
928
929 struct OpJacobian
931
932 const int tAg;
933 boost::shared_ptr<VectorDouble> activeVariablesPtr;
934 boost::shared_ptr<VectorDouble> resultsPtr;
935 boost::shared_ptr<MatrixDouble> jacobianPtr;
937
938 const double &aLpha;
939
940 OpJacobian(int tag, const std::string field_name,
941 boost::shared_ptr<VectorDouble> &active_variables_ptr,
942 boost::shared_ptr<VectorDouble> &results_ptr,
943 boost::shared_ptr<MatrixDouble> &jacobian_ptr,
944 bool evaluate_jacobian, const double &alpha)
947 tAg(tag), activeVariablesPtr(active_variables_ptr),
948 resultsPtr(results_ptr), jacobianPtr(jacobian_ptr),
949 evaluateJacobian(evaluate_jacobian), aLpha(alpha) {}
950
951 MoFEMErrorCode doWork(int side, EntityType type,
952 EntitiesFieldData::EntData &data) {
954 if (type != MBVERTEX)
956
957 FTensor::Index<'i', 3> i;
958 FTensor::Index<'j', 2> j;
961
962 Tag th0, th1, th2, th3;
964 th1, th2, th3);
965 FTensor::Tensor1<double, 3> t_edge_base0, t_edge_base1;
967 CHKERR getEdgeFE()->mField.get_moab().tag_get_data(th0, &fe_ent, 1,
968 &t_edge_base0(0));
969 CHKERR getEdgeFE()->mField.get_moab().tag_get_data(th1, &fe_ent, 1,
970 &t_edge_base1(0));
971
972 VectorInt &indices = data.getIndices();
973
974 trace_on(tAg);
975
976 ublas::vector<adouble> lambda_dofs(4);
977 for (int dd = 0; dd != 4; ++dd) {
978 lambda_dofs[dd] <<= (*activeVariablesPtr)[dd];
979 }
980 ublas::vector<adouble> position_dofs(6);
981 for (int dd = 0; dd != 6; ++dd) {
982 position_dofs[dd] <<= (*activeVariablesPtr)[4 + dd];
983 }
984
986 &position_dofs[0], &position_dofs[1], &position_dofs[2]);
988 &position_dofs[3], &position_dofs[4], &position_dofs[5]);
989
991 t_tangent(i) = t_node1(i) - t_node0(i);
992 adouble l = sqrt(t_tangent(i) * t_tangent(i));
993 t_tangent(i) /= l;
994
995 adouble t_dot0, t_dot1;
996 t_dot0 = t_edge_base0(i) * t_tangent(i);
997 t_dot1 = t_edge_base1(i) * t_tangent(i);
998
999 FTensor::Tensor1<adouble, 3> t_base0, t_base1;
1000 t_base0(i) = t_edge_base0(i) - t_dot0 * t_tangent(i);
1001 t_base1(i) = t_edge_base1(i) - t_dot1 * t_tangent(i);
1002 t_base0(i) /= sqrt(t_base0(i) * t_base0(i));
1003 t_base1(i) /= sqrt(t_base1(i) * t_base1(i));
1004
1005 auto t_base_fun1 = data.getFTensor0N();
1006 auto t_base_fun2 = data.getFTensor0N();
1010 auto t_coord_ref = getFTensor1CoordsAtGaussPts();
1011
1012 ublas::vector<adouble> c_vec(4);
1013 ublas::vector<adouble> f_vec(6);
1014 c_vec.clear();
1015 f_vec.clear();
1016
1017 int nb_gauss_pts = data.getN().size1();
1018 int nb_base_functions = data.getN().size2();
1019 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1020
1021 FTensor::Tensor1<adouble *, 3> t_position_dofs(
1022 &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
1023 FTensor::Tensor1<adouble *, 2> t_lambda_dof(&lambda_dofs[0],
1024 &lambda_dofs[1], 2);
1025
1026 t_position(i) = 0;
1027 t_lambda(j) = 0;
1028 for (int bb = 0; bb != nb_base_functions; ++bb) {
1029 t_position(i) += t_base_fun1 * t_position_dofs(i);
1030 t_lambda(j) += t_base_fun1 * t_lambda_dof(j);
1031 ++t_base_fun1;
1032 ++t_position_dofs;
1033 ++t_lambda_dof;
1034 }
1035
1036 t_delta(i) = t_position(i) - t_coord_ref(i);
1037 adouble dot0 = t_base0(i) * t_delta(i);
1038 adouble dot1 = t_base1(i) * t_delta(i);
1039
1040 adouble w = getGaussPts()(1, gg) * l * aLpha;
1041 adouble val, val1, val2;
1042 FTensor::Tensor1<adouble *, 2> t_c(&c_vec[0], &c_vec[1], 2);
1043 FTensor::Tensor1<adouble *, 3> t_f(&f_vec[0], &f_vec[1], &f_vec[2], 3);
1044 for (int bb = 0; bb != nb_base_functions; ++bb) {
1045 if (indices[2 * bb] != -1) {
1046 val = w * t_base_fun2;
1047 t_c(N0) += val * dot0;
1048 t_c(N1) += val * dot1;
1049 val1 = val * t_lambda(N0);
1050 val2 = val * t_lambda(N1);
1051 t_f(i) += val1 * t_base0(i) + val2 * t_base1(i);
1052 }
1053 ++t_c;
1054 ++t_f;
1055 ++t_base_fun2;
1056 }
1057
1058 ++t_coord_ref;
1059 }
1060
1061 for (int rr = 0; rr != 4; ++rr) {
1062 c_vec[rr] >>= (*resultsPtr)[rr];
1063 }
1064 for (int rr = 0; rr != 6; ++rr) {
1065 f_vec(rr) >>= (*resultsPtr)[4 + rr];
1066 }
1067
1068 trace_off();
1069
1070 if (evaluateJacobian) {
1071 double *jac_ptr[4 + 6];
1072 for (int rr = 0; rr != 4 + 6; ++rr) {
1073 jac_ptr[rr] = &(*jacobianPtr)(rr, 0);
1074 }
1075 // play recorder for jacobians
1076 int r =
1077 ::jacobian(tAg, 4 + 6, 4 + 6, &(*activeVariablesPtr)[0], jac_ptr);
1078 if (r < 0) {
1079 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1080 "ADOL-C function evaluation with error");
1081 }
1082 }
1083
1085 }
1086 };
1087
1088 MoFEMErrorCode setOperators(int tag, Range edges, Range faces,
1089 const std::string lagrange_multipliers_field_name,
1090 const std::string material_field_name) {
1093 edges, faces);
1094 CHKERR setOperators(tag, lagrange_multipliers_field_name,
1095 material_field_name);
1097 }
1098
1099 MoFEMErrorCode setOperators(int tag,
1100 const std::string lagrange_multipliers_field_name,
1101 const std::string material_field_name,
1102 const double *alpha = nullptr) {
1104
1105 if (alpha) {
1106 aLpha = *alpha;
1107 }
1108
1109 boost::shared_ptr<VectorDouble> active_variables_ptr(
1110 new VectorDouble(4 + 6));
1111 boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(4 + 6));
1112 boost::shared_ptr<MatrixDouble> jacobian_ptr(
1113 new MatrixDouble(4 + 6, 4 + 6));
1114
1115 feRhs.getOpPtrVector().clear();
1117 lagrange_multipliers_field_name, active_variables_ptr));
1119 material_field_name, active_variables_ptr));
1120 feRhs.getOpPtrVector().push_back(new OpJacobian(
1121 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1122 jacobian_ptr, false, aLpha));
1123 feRhs.getOpPtrVector().push_back(
1124 new OpAssembleRhs<4, 6>(lagrange_multipliers_field_name, results_ptr));
1125 feRhs.getOpPtrVector().push_back(
1126 new OpAssembleRhs<4, 6>(material_field_name, results_ptr));
1127
1128 // Adding operators to calculate the left hand side
1129 feLhs.getOpPtrVector().clear();
1131 lagrange_multipliers_field_name, active_variables_ptr));
1133 material_field_name, active_variables_ptr));
1134 feLhs.getOpPtrVector().push_back(new OpJacobian(
1135 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1136 jacobian_ptr, true, aLpha));
1138 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1140 material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
1142 material_field_name, material_field_name, jacobian_ptr));
1143
1145 }
1146
1147 MoFEMErrorCode
1149 const std::string lagrange_multipliers_field_name,
1150 const std::string material_field_name) {
1152
1154 edges, faces);
1155 CHKERR setOperatorsConstrainOnly(tag, lagrange_multipliers_field_name,
1156 material_field_name);
1158 }
1159
1160 MoFEMErrorCode
1162 const std::string lagrange_multipliers_field_name,
1163 const std::string material_field_name) {
1165
1166 boost::shared_ptr<VectorDouble> active_variables_ptr(
1167 new VectorDouble(4 + 6));
1168 boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(4 + 6));
1169 boost::shared_ptr<MatrixDouble> jacobian_ptr(
1170 new MatrixDouble(4 + 6, 4 + 6));
1171
1172 // Adding operators to calculate the left hand side
1173 feLhs.getOpPtrVector().clear();
1175 lagrange_multipliers_field_name, active_variables_ptr));
1177 material_field_name, active_variables_ptr));
1178 feLhs.getOpPtrVector().push_back(new OpJacobian(
1179 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1180 jacobian_ptr, true, aLpha));
1182 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1183
1185 }
1186};
1187
1188#endif // WITH_ADOL_C
1189
1190#endif // __SURFACE_SLIDING_CONSTRAINS_HPP__
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static PetscErrorCode ierr
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
constexpr int order
FTensor::Index< 'i', SPACE_DIM > i
static double lambda
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double tol
auto cross(const Tensor1_Expr< A, T, 3, i > &a, const Tensor1_Expr< B, U, 3, j > &b, const Index< k, 3 > &)
Definition cross.hpp:10
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
constexpr auto field_name
static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges, Range tris)
static MoFEMErrorCode createTag(moab::Interface &moab, Tag &th0, Tag &th1, Tag &th2, Tag &th3)
static MoFEMErrorCode setTags(moab::Interface &moab, Range edges, Range tris, bool number_pathes=true, boost::shared_ptr< SurfaceSlidingConstrains::DriverElementOrientation > surface_orientation=nullptr, MoFEM::Interface *m_field_ptr=nullptr)
static MoFEMErrorCode saveEdges(moab::Interface &moab, std::string name, Range edges, Range *faces=nullptr)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::shared_ptr< VectorDouble > activeVariablesPtr
boost::shared_ptr< VectorDouble > resultsPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpJacobian(int tag, const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr, boost::shared_ptr< VectorDouble > &results_ptr, boost::shared_ptr< MatrixDouble > &jacobian_ptr, bool evaluate_jacobian, const double &alpha)
boost::shared_ptr< MatrixDouble > jacobianPtr
MoFEMErrorCode setOperators(int tag, Range edges, Range faces, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
EdgeSlidingConstrains(MoFEM::Interface &m_field)
boost::shared_ptr< MyEdgeFE > feRhsPtr
MoFEMErrorCode setOperators(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name, const double *alpha=nullptr)
boost::shared_ptr< MyEdgeFE > feLhsPtr
MoFEMErrorCode setOperatorsConstrainOnly(int tag, Range edges, Range faces, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
MoFEMErrorCode setOperatorsConstrainOnly(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
boost::shared_ptr< MatrixDouble > jacobianPtr
OpAssembleLhs(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > &jacobian_ptr)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpAssembleRhs(const std::string field_name, boost::shared_ptr< VectorDouble > &results_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VectorDouble > resultsPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpGetActiveDofsLambda(const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr)
boost::shared_ptr< VectorDouble > activeVariablesPtr
OpGetActiveDofsPositions(const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VectorDouble > activeVariablesPtr
Implementation of surface sliding constrains.
virtual ~GenericSliding()=default
GenericSliding()=default
virtual moab::Interface & get_moab()=0
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
const EdgeElementForcesAndSourcesCore * getEdgeFE()
get pointer to this finite element
FaceElementForcesAndSourcesCore * getFaceFE()
return pointer to Generic Triangle Finite Element object
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
std::string getFEName() const
Get name of the element.
@ OPCOL
operator doWork function is executed on FE columns
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Vec & snes_f
residual
Mat & snes_B
preconditioner of jacobian matrix
Class implemented by user to detect face orientation.
virtual MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field, const EntityHandle face)
virtual MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field, const FEMethod *fe_method_ptr)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
OpJacobian(int tag, const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr, boost::shared_ptr< VectorDouble > &results_ptr, boost::shared_ptr< MatrixDouble > &jacobian_ptr, DriverElementOrientation &orientation, bool evaluate_jacobian, const double &alpha)
boost::shared_ptr< VectorDouble > resultsPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VectorDouble > activeVariablesPtr
boost::shared_ptr< MatrixDouble > jacobianPtr
Shape preserving constrains, i.e. nodes sliding on body surface.
boost::shared_ptr< MyTriangleFE > feRhsPtr
virtual ~SurfaceSlidingConstrains()=default
DriverElementOrientation & crackFrontOrientation
boost::shared_ptr< MyTriangleFE > feLhsPtr
SurfaceSlidingConstrains(MoFEM::Interface &m_field, DriverElementOrientation &orientation)
MoFEMErrorCode setOperators(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name, const double *alpha=nullptr)
MoFEMErrorCode setOperatorsConstrainOnly(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name)