v0.13.1
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) {}
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>
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) {}
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>
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
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 SETERRQ2(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>
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,
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 SETERRQ1(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
340 const FEMethod *fe_method_ptr) {
344 }
346 const EntityHandle face) {
350 }
351 };
352
354
356
357 Mat B;
359
361 : MoFEM::FaceElementForcesAndSourcesCore(m_field), B(PETSC_NULL),
362 F(PETSC_NULL) {}
363 int getRule(int order) { return 2 * order; };
364
367
369
370 if (B != PETSC_NULL) {
371 snes_B = B;
372 }
373
374 if (F != PETSC_NULL) {
375 snes_f = F;
376 }
377
379 }
380 };
381
382 boost::shared_ptr<MyTriangleFE> feRhsPtr, feLhsPtr;
383
388
390
391 double aLpha;
394 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
395 "Get surface sliding constrains element scaling",
396 "none");
397 CHKERRQ(ierr);
398 CHKERR PetscOptionsScalar("-surface_sliding_alpha", "scaling parameter", "",
399 aLpha, &aLpha, PETSC_NULL);
400 ierr = PetscOptionsEnd();
401 CHKERRQ(ierr);
403 }
404
406 DriverElementOrientation &orientation)
407 : mField(m_field), feRhsPtr(new MyTriangleFE(m_field)),
408 feLhsPtr(new MyTriangleFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr),
409 crackFrontOrientation(orientation), aLpha(1) {
410
411 ierr = getOptions();
412 CHKERRABORT(PETSC_COMM_WORLD, ierr);
413 }
414
415 virtual ~SurfaceSlidingConstrains() = default;
416
419
420 const int tAg;
421 boost::shared_ptr<VectorDouble> activeVariablesPtr;
422 boost::shared_ptr<VectorDouble> resultsPtr;
423 boost::shared_ptr<MatrixDouble> jacobianPtr;
426
427 const double &aLpha;
428
429 OpJacobian(int tag, const std::string field_name,
430 boost::shared_ptr<VectorDouble> &active_variables_ptr,
431 boost::shared_ptr<VectorDouble> &results_ptr,
432 boost::shared_ptr<MatrixDouble> &jacobian_ptr,
433 DriverElementOrientation &orientation, bool evaluate_jacobian,
434 const double &alpha)
437 tAg(tag), activeVariablesPtr(active_variables_ptr),
438 resultsPtr(results_ptr), jacobianPtr(jacobian_ptr),
439 oRientation(orientation), evaluateJacobian(evaluate_jacobian),
440 aLpha(alpha) {}
441
445 if (type != MBVERTEX)
447
448 VectorInt &indices = data.getIndices();
449
450 trace_on(tAg);
451
452 ublas::vector<adouble> lambda_dofs(3);
453 for (int dd = 0; dd != 3; ++dd) {
454 lambda_dofs[dd] <<= (*activeVariablesPtr)[dd];
455 }
456 ublas::vector<adouble> position_dofs(9);
457 for (int dd = 0; dd != 9; ++dd) {
458 position_dofs[dd] <<= (*activeVariablesPtr)[3 + dd];
459 }
460
467
469 getFEMethod());
471
472 int nb_gauss_pts = data.getN().size1();
473 int nb_base_functions = data.getN().size2();
474 auto t_base1 = data.getFTensor0N();
475 auto t_base2 = data.getFTensor0N();
476 auto t_diff_base = data.getFTensor1DiffN<2>();
478 FTensor::Tensor1<adouble, 3> t_position_ksi;
479 FTensor::Tensor1<adouble, 3> t_position_eta;
483
484 ublas::vector<adouble> c_vec(3);
485 ublas::vector<adouble> f_vec(9);
486 c_vec.clear();
487 f_vec.clear();
488
489 auto t_coord_ref = getFTensor1CoordsAtGaussPts();
490
491 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
492
493 FTensor::Tensor1<adouble *, 3> t_position_dofs(
494 &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
495 FTensor::Tensor0<adouble *> t_lambda_dof(&lambda_dofs[0]);
496
497 t_position(i) = 0;
498 t_position_ksi(i) = 0;
499 t_position_eta(i) = 0;
500 lambda = 0;
501
502 for (int bb = 0; bb != nb_base_functions; ++bb) {
503 t_position(i) += t_base1 * t_position_dofs(i);
504 t_position_ksi(i) += t_diff_base(N0) * t_position_dofs(i);
505 t_position_eta(i) += t_diff_base(N1) * t_position_dofs(i);
506 lambda += t_base1 * t_lambda_dof;
507 ++t_base1;
508 ++t_position_dofs;
509 ++t_lambda_dof;
510 ++t_diff_base;
511 }
512
513 t_delta(i) = t_position(i) - t_coord_ref(i);
514 t_normal(k) = FTensor::cross(t_position_ksi(i), t_position_eta(j), k);
515
516 double w = getGaussPts()(2, gg) * 0.5 * aLpha;
517 adouble val;
518 FTensor::Tensor0<adouble *> t_c(&c_vec[0]);
519 FTensor::Tensor1<adouble *, 3> t_f(&f_vec[0], &f_vec[1], &f_vec[2], 3);
520
521 adouble area = sqrt(t_normal(i) * t_normal(i));
522
523 for (int bb = 0; bb != nb_base_functions; ++bb) {
524 if (indices[bb] != -1) {
525 val = w * eo * t_base2;
526 t_c += val * t_normal(i) * t_delta(i);
527 val *= lambda;
528 t_f(i) += val * t_normal(i);
529 }
530 ++t_c;
531 ++t_f;
532 ++t_base2;
533 }
534
535 ++t_coord_ref;
536 }
537
538 for (int rr = 0; rr != 3; ++rr)
539 c_vec[rr] >>= (*resultsPtr)[rr];
540
541 for (int rr = 0; rr != 9; ++rr)
542 f_vec(rr) >>= (*resultsPtr)[3 + rr];
543
544 trace_off();
545
546 if (evaluateJacobian) {
547 double *jac_ptr[3 + 9];
548 for (int rr = 0; rr != 3 + 9; ++rr)
549 jac_ptr[rr] = &(*jacobianPtr)(rr, 0);
550
551 // play recorder for jacobians
552 int r =
553 ::jacobian(tAg, 3 + 9, 3 + 9, &(*activeVariablesPtr)[0], jac_ptr);
554 if (r < 0)
555 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
556 "ADOL-C function evaluation with error");
557 }
558
560 }
561 };
562
564 const std::string lagrange_multipliers_field_name,
565 const std::string material_field_name,
566 const double *alpha = nullptr) {
568
569 if (alpha != nullptr) {
570 aLpha = *alpha;
571 }
572
573 boost::shared_ptr<VectorDouble> active_variables_ptr(
574 new VectorDouble(3 + 9));
575 boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
576 boost::shared_ptr<MatrixDouble> jacobian_ptr(
577 new MatrixDouble(3 + 9, 3 + 9));
578
579 feRhs.getOpPtrVector().clear();
581 lagrange_multipliers_field_name, active_variables_ptr));
583 material_field_name, active_variables_ptr));
584 feRhs.getOpPtrVector().push_back(new OpJacobian(
585 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
586 jacobian_ptr, crackFrontOrientation, false, aLpha));
587 feRhs.getOpPtrVector().push_back(
588 new OpAssembleRhs<3, 9>(lagrange_multipliers_field_name, results_ptr));
589 feRhs.getOpPtrVector().push_back(
590 new OpAssembleRhs<3, 9>(material_field_name, results_ptr));
591
592 // Adding operators to calculate the left hand side
593 feLhs.getOpPtrVector().clear();
595 lagrange_multipliers_field_name, active_variables_ptr));
597 material_field_name, active_variables_ptr));
598 feLhs.getOpPtrVector().push_back(new OpJacobian(
599 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
600 jacobian_ptr, crackFrontOrientation, true, aLpha));
602 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
604 material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
606 material_field_name, material_field_name, jacobian_ptr));
607
609 }
610
613 const std::string lagrange_multipliers_field_name,
614 const std::string material_field_name) {
616
617 boost::shared_ptr<VectorDouble> active_variables_ptr(
618 new VectorDouble(3 + 9));
619 boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
620 boost::shared_ptr<MatrixDouble> jacobian_ptr(
621 new MatrixDouble(3 + 9, 3 + 9));
622
623 // Adding operators to calculate the left hand side
624 feLhs.getOpPtrVector().clear();
626 lagrange_multipliers_field_name, active_variables_ptr));
628 material_field_name, active_variables_ptr));
629 feLhs.getOpPtrVector().push_back(new OpJacobian(
630 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
631 jacobian_ptr, crackFrontOrientation, true, aLpha));
633 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
634
636 }
637};
638
640
642
643 static MoFEMErrorCode createTag(moab::Interface &moab, Tag &th0, Tag &th1,
644 Tag &th2, Tag &th3) {
646 double def_val[] = {0, 0, 0};
647 CHKERR moab.tag_get_handle("EDGE_BASE0", 3, MB_TYPE_DOUBLE, th0,
648 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
649 CHKERR moab.tag_get_handle("EDGE_BASE1", 3, MB_TYPE_DOUBLE, th1,
650 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
651 int def_numb_val[] = {-1};
652 CHKERR moab.tag_get_handle("PATCH_NUMBER", 1, MB_TYPE_INTEGER, th2,
653 MB_TAG_CREAT | MB_TAG_SPARSE, def_numb_val);
654 int def_orientation_val[] = {1};
655 CHKERR moab.tag_get_handle("PATCH_ORIENTATION", 1, MB_TYPE_INTEGER, th3,
656 MB_TAG_CREAT | MB_TAG_SPARSE,
657 def_orientation_val);
658
660 }
661
663 Range tris) {
665
666 auto get_edges = [&](const Range &ents) {
667 Range edges;
668 CHKERR moab.get_adjacencies(ents, 1, false, edges,
669 moab::Interface::UNION);
670 return edges;
671 };
672
673 auto get_face_adj = [&, get_edges](const Range &faces) {
674 Range adj_faces;
675 CHKERR moab.get_adjacencies(subtract(get_edges(faces), edges), 2, false,
676 adj_faces, moab::Interface::UNION);
677 return intersect(adj_faces, tris);
678 };
679
680 auto get_patch = [&, get_face_adj](const EntityHandle face) {
681 Range patch_ents;
682 patch_ents.insert(face);
683 unsigned int nb0;
684 do {
685 nb0 = patch_ents.size();
686 patch_ents.merge(get_face_adj(patch_ents));
687 } while (nb0 != patch_ents.size());
688 return patch_ents;
689 };
690
691 auto get_patches = [&]() {
692 std::vector<Range> patches;
693 while (!tris.empty()) {
694 patches.push_back(get_patch(tris[0]));
695 tris = subtract(tris, patches.back());
696 }
697 return patches;
698 };
699
700 Tag th0, th1, th2, th3;
701 CHKERR createTag(moab, th0, th1, th2, th3);
702
703 auto patches = get_patches();
704 int pp = 0;
705 for (auto patch : patches) {
706 // cerr << "pp: " << pp << endl;
707 // cerr << patch << endl;
708 std::vector<int> tags_vals(patch.size(), pp);
709 CHKERR moab.tag_set_data(th2, patch, &*tags_vals.begin());
710 ++pp;
711 }
712
714 }
715
717 moab::Interface &moab, Range edges, Range tris,
718 bool number_pathes = true,
719 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
720 surface_orientation = nullptr,
721 MoFEM::Interface *m_field_ptr = nullptr) {
723 Tag th0, th1, th2, th3;
724 CHKERR createTag(moab, th0, th1, th2, th3);
725 if (number_pathes) {
726 CHKERR numberSurfaces(moab, edges, tris);
727 }
728 for (auto edge : edges) {
729 Range adj_faces;
730 CHKERR moab.get_adjacencies(&edge, 1, 2, false, adj_faces);
731 adj_faces = intersect(adj_faces, tris);
732 if (adj_faces.size() != 2) {
733 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
734 "Should be 2 faces adjacent to edge but is %d",
735 adj_faces.size());
736 }
738 auto get_tensor_from_vec = [](VectorDouble3 &v) {
740 &v[0], &v[1], &v[2]);
741 };
742
746
747 std::array<double, 3> areas;
748 auto calculate_normals = [&, get_tensor_from_vec]() {
749 int ff = 0;
750 for (auto face : adj_faces) {
751 double &x = (v[ff])[0];
752 double &y = (v[ff])[1];
753 double &z = (v[ff])[2];
754 moab::Util::normal(&moab, face, x, y, z);
755 auto t_n = get_tensor_from_vec(v[ff]);
756 areas[ff] = sqrt(t_n(i) * t_n(i));
757 t_n(i) /= areas[ff];
758 int orientation;
759 // FIXME: this tag is empty
760 CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
761 if (orientation == -1) {
762 t_n(i) *= -1;
763 }
764 if (surface_orientation && m_field_ptr) {
765 CHKERR surface_orientation->getElementOrientation(*m_field_ptr,
766 face);
767 int eo = surface_orientation->elementOrientation;
768 if (eo == -1) {
769 t_n(i) *= -1;
770 }
771 }
772 ++ff;
773 }
774 };
775 calculate_normals();
776
777 auto get_patch_number = [&]() {
778 std::vector<int> p = {0, 0};
779 CHKERR moab.tag_get_data(th2, adj_faces, &*p.begin());
780 return p;
781 };
782
783 std::array<EntityHandle, 2> faces_handles = {adj_faces[0],
784 adj_faces[1]};
785 auto order_normals = [&, get_patch_number]() {
786 auto p = get_patch_number();
787 if (p[0] < p[1]) {
788 v[0].swap(v[1]);
789 faces_handles[0] = adj_faces[1];
790 faces_handles[1] = adj_faces[0];
791 }
792 };
793 order_normals();
794
796 auto t_n0 = get_tensor_from_vec(v[0]);
797 auto t_n1 = get_tensor_from_vec(v[1]);
798 t_cross(k) = FTensor::cross(t_n0(i), t_n1(j), k);
799
800 auto get_base_for_coplanar_faces = [&]() {
802
803 std::vector<EntityHandle> face_conn;
804 CHKERR moab.get_connectivity(&faces_handles[1], 1, face_conn, true);
805 std::vector<EntityHandle> edge_conn;
806 CHKERR moab.get_connectivity(&edge, 1, edge_conn, true);
807 std::sort(face_conn.begin(), face_conn.end());
808 std::sort(edge_conn.begin(), edge_conn.end());
809 int n = 0;
810 for (; n != 2; ++n)
811 if (face_conn[n] != edge_conn[n])
812 break;
813 VectorDouble3 coords_face(3);
814 CHKERR moab.get_coords(&face_conn[n], 1,
815 &*coords_face.data().begin());
816 VectorDouble6 coords_edge(6);
817 CHKERR moab.get_coords(&edge_conn[0], 2,
818 &*coords_edge.data().begin());
819 auto t_edge_n0 = FTensor::Tensor1<double, 3>(
820 coords_edge[0], coords_edge[1], coords_edge[2]);
821 auto t_edge_n1 = FTensor::Tensor1<double, 3>(
822 coords_edge[3], coords_edge[4], coords_edge[5]);
823 auto t_face_n = get_tensor_from_vec(coords_face);
824 t_face_n(i) -= t_edge_n0(i);
826 t_delta(i) = t_edge_n1(i) - t_edge_n0(i);
827 t_delta(i) /= sqrt(t_delta(i) * t_delta(i));
828 t_n1(k) = FTensor::cross(t_n0(i), t_delta(j), k);
829 if (t_n1(i) * t_face_n(i) < 0)
830 t_n1(i) *= -1;
831
833 };
834
835 auto get_base_for_not_planar_faces = [&]() {
837 t_n1(k) = FTensor::cross(t_n0(i), t_cross(j), k);
838 t_n1(i) /= sqrt(t_n1(i) * t_n1(i));
840 };
841
842 // This a case when faces adjacent to edge are coplanar !!!
843 constexpr double tol = 1e-6;
844 if ((t_cross(i) * t_cross(i)) < tol)
845 CHKERR get_base_for_coplanar_faces();
846 else
847 CHKERR get_base_for_not_planar_faces();
848
849 VectorDouble3 &v0 = v[0];
850 VectorDouble3 &v1 = v[1];
851 CHKERR moab.tag_set_data(th0, &edge, 1, &v0[0]);
852 CHKERR moab.tag_set_data(th1, &edge, 1, &v1[0]);
853 }
855 }
856
857 static MoFEMErrorCode saveEdges(moab::Interface &moab, std::string name,
858 Range edges, Range *faces = nullptr) {
860 EntityHandle meshset;
861 Tag ths[4];
862 CHKERR createTag(moab, ths[0], ths[1], ths[2], ths[3]);
863 CHKERR moab.create_meshset(MESHSET_SET, meshset);
864 CHKERR moab.add_entities(meshset, edges);
865 if (faces != nullptr) {
866 CHKERR moab.add_entities(meshset, *faces);
867 }
868 CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1, ths, 4);
869 CHKERR moab.delete_entities(&meshset, 1);
871 }
872 };
873
875
877
878 Mat B;
880
882 : MoFEM::EdgeElementForcesAndSourcesCore(m_field), B(PETSC_NULL),
883 F(PETSC_NULL) {}
884 int getRule(int order) { return 2 * order; };
885
888
890
891 if (B != PETSC_NULL) {
892 snes_B = B;
893 }
894
895 if (F != PETSC_NULL) {
896 snes_f = F;
897 }
898
900 }
901 };
902
903 boost::shared_ptr<MyEdgeFE> feRhsPtr, feLhsPtr;
904
909
910 double aLpha;
911
914 ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
915 "Get edge sliding constrains element scaling",
916 "none");
917 CHKERRQ(ierr);
918 CHKERR PetscOptionsScalar("-edge_sliding_alpha", "scaling parameter", "",
919 aLpha, &aLpha, PETSC_NULL);
920 ierr = PetscOptionsEnd();
921 CHKERRQ(ierr);
923 }
924
926 : mField(m_field), feRhsPtr(new MyEdgeFE(m_field)),
927 feLhsPtr(new MyEdgeFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr),
928 aLpha(1) {
929 ierr = getOptions();
930 CHKERRABORT(PETSC_COMM_WORLD, ierr);
931 }
932
935
936 const int tAg;
937 boost::shared_ptr<VectorDouble> activeVariablesPtr;
938 boost::shared_ptr<VectorDouble> resultsPtr;
939 boost::shared_ptr<MatrixDouble> jacobianPtr;
941
942 const double &aLpha;
943
944 OpJacobian(int tag, const std::string field_name,
945 boost::shared_ptr<VectorDouble> &active_variables_ptr,
946 boost::shared_ptr<VectorDouble> &results_ptr,
947 boost::shared_ptr<MatrixDouble> &jacobian_ptr,
948 bool evaluate_jacobian, const double &alpha)
951 tAg(tag), activeVariablesPtr(active_variables_ptr),
952 resultsPtr(results_ptr), jacobianPtr(jacobian_ptr),
953 evaluateJacobian(evaluate_jacobian), aLpha(alpha) {}
954
958 if (type != MBVERTEX)
960
965
966 Tag th0, th1, th2, th3;
968 th1, th2, th3);
969 FTensor::Tensor1<double, 3> t_edge_base0, t_edge_base1;
970 EntityHandle fe_ent = getFEEntityHandle();
971 CHKERR getEdgeFE()->mField.get_moab().tag_get_data(th0, &fe_ent, 1,
972 &t_edge_base0(0));
973 CHKERR getEdgeFE()->mField.get_moab().tag_get_data(th1, &fe_ent, 1,
974 &t_edge_base1(0));
975
976 VectorInt &indices = data.getIndices();
977
978 trace_on(tAg);
979
980 ublas::vector<adouble> lambda_dofs(4);
981 for (int dd = 0; dd != 4; ++dd) {
982 lambda_dofs[dd] <<= (*activeVariablesPtr)[dd];
983 }
984 ublas::vector<adouble> position_dofs(6);
985 for (int dd = 0; dd != 6; ++dd) {
986 position_dofs[dd] <<= (*activeVariablesPtr)[4 + dd];
987 }
988
990 &position_dofs[0], &position_dofs[1], &position_dofs[2]);
992 &position_dofs[3], &position_dofs[4], &position_dofs[5]);
993
995 t_tangent(i) = t_node1(i) - t_node0(i);
996 adouble l = sqrt(t_tangent(i) * t_tangent(i));
997 t_tangent(i) /= l;
998
999 adouble t_dot0, t_dot1;
1000 t_dot0 = t_edge_base0(i) * t_tangent(i);
1001 t_dot1 = t_edge_base1(i) * t_tangent(i);
1002
1003 FTensor::Tensor1<adouble, 3> t_base0, t_base1;
1004 t_base0(i) = t_edge_base0(i) - t_dot0 * t_tangent(i);
1005 t_base1(i) = t_edge_base1(i) - t_dot1 * t_tangent(i);
1006 t_base0(i) /= sqrt(t_base0(i) * t_base0(i));
1007 t_base1(i) /= sqrt(t_base1(i) * t_base1(i));
1008
1009 auto t_base_fun1 = data.getFTensor0N();
1010 auto t_base_fun2 = data.getFTensor0N();
1014 auto t_coord_ref = getFTensor1CoordsAtGaussPts();
1015
1016 ublas::vector<adouble> c_vec(4);
1017 ublas::vector<adouble> f_vec(6);
1018 c_vec.clear();
1019 f_vec.clear();
1020
1021 int nb_gauss_pts = data.getN().size1();
1022 int nb_base_functions = data.getN().size2();
1023 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1024
1025 FTensor::Tensor1<adouble *, 3> t_position_dofs(
1026 &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
1027 FTensor::Tensor1<adouble *, 2> t_lambda_dof(&lambda_dofs[0],
1028 &lambda_dofs[1], 2);
1029
1030 t_position(i) = 0;
1031 t_lambda(j) = 0;
1032 for (int bb = 0; bb != nb_base_functions; ++bb) {
1033 t_position(i) += t_base_fun1 * t_position_dofs(i);
1034 t_lambda(j) += t_base_fun1 * t_lambda_dof(j);
1035 ++t_base_fun1;
1036 ++t_position_dofs;
1037 ++t_lambda_dof;
1038 }
1039
1040 t_delta(i) = t_position(i) - t_coord_ref(i);
1041 adouble dot0 = t_base0(i) * t_delta(i);
1042 adouble dot1 = t_base1(i) * t_delta(i);
1043
1044 adouble w = getGaussPts()(1, gg) * l * aLpha;
1045 adouble val, val1, val2;
1046 FTensor::Tensor1<adouble *, 2> t_c(&c_vec[0], &c_vec[1], 2);
1047 FTensor::Tensor1<adouble *, 3> t_f(&f_vec[0], &f_vec[1], &f_vec[2], 3);
1048 for (int bb = 0; bb != nb_base_functions; ++bb) {
1049 if (indices[2 * bb] != -1) {
1050 val = w * t_base_fun2;
1051 t_c(N0) += val * dot0;
1052 t_c(N1) += val * dot1;
1053 val1 = val * t_lambda(N0);
1054 val2 = val * t_lambda(N1);
1055 t_f(i) += val1 * t_base0(i) + val2 * t_base1(i);
1056 }
1057 ++t_c;
1058 ++t_f;
1059 ++t_base_fun2;
1060 }
1061
1062 ++t_coord_ref;
1063 }
1064
1065 for (int rr = 0; rr != 4; ++rr) {
1066 c_vec[rr] >>= (*resultsPtr)[rr];
1067 }
1068 for (int rr = 0; rr != 6; ++rr) {
1069 f_vec(rr) >>= (*resultsPtr)[4 + rr];
1070 }
1071
1072 trace_off();
1073
1074 if (evaluateJacobian) {
1075 double *jac_ptr[4 + 6];
1076 for (int rr = 0; rr != 4 + 6; ++rr) {
1077 jac_ptr[rr] = &(*jacobianPtr)(rr, 0);
1078 }
1079 // play recorder for jacobians
1080 int r =
1081 ::jacobian(tAg, 4 + 6, 4 + 6, &(*activeVariablesPtr)[0], jac_ptr);
1082 if (r < 0) {
1083 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1084 "ADOL-C function evaluation with error");
1085 }
1086 }
1087
1089 }
1090 };
1091
1093 const std::string lagrange_multipliers_field_name,
1094 const std::string material_field_name) {
1097 edges, faces);
1098 CHKERR setOperators(tag, lagrange_multipliers_field_name,
1099 material_field_name);
1101 }
1102
1104 const std::string lagrange_multipliers_field_name,
1105 const std::string material_field_name,
1106 const double *alpha = nullptr) {
1108
1109 if (alpha) {
1110 aLpha = *alpha;
1111 }
1112
1113 boost::shared_ptr<VectorDouble> active_variables_ptr(
1114 new VectorDouble(4 + 6));
1115 boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(4 + 6));
1116 boost::shared_ptr<MatrixDouble> jacobian_ptr(
1117 new MatrixDouble(4 + 6, 4 + 6));
1118
1119 feRhs.getOpPtrVector().clear();
1121 lagrange_multipliers_field_name, active_variables_ptr));
1123 material_field_name, active_variables_ptr));
1124 feRhs.getOpPtrVector().push_back(new OpJacobian(
1125 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1126 jacobian_ptr, false, aLpha));
1127 feRhs.getOpPtrVector().push_back(
1128 new OpAssembleRhs<4, 6>(lagrange_multipliers_field_name, results_ptr));
1129 feRhs.getOpPtrVector().push_back(
1130 new OpAssembleRhs<4, 6>(material_field_name, results_ptr));
1131
1132 // Adding operators to calculate the left hand side
1133 feLhs.getOpPtrVector().clear();
1135 lagrange_multipliers_field_name, active_variables_ptr));
1137 material_field_name, active_variables_ptr));
1138 feLhs.getOpPtrVector().push_back(new OpJacobian(
1139 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1140 jacobian_ptr, true, aLpha));
1142 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1144 material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
1146 material_field_name, material_field_name, jacobian_ptr));
1147
1149 }
1150
1153 const std::string lagrange_multipliers_field_name,
1154 const std::string material_field_name) {
1156
1158 edges, faces);
1159 CHKERR setOperatorsConstrainOnly(tag, lagrange_multipliers_field_name,
1160 material_field_name);
1162 }
1163
1166 const std::string lagrange_multipliers_field_name,
1167 const std::string material_field_name) {
1169
1170 boost::shared_ptr<VectorDouble> active_variables_ptr(
1171 new VectorDouble(4 + 6));
1172 boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(4 + 6));
1173 boost::shared_ptr<MatrixDouble> jacobian_ptr(
1174 new MatrixDouble(4 + 6, 4 + 6));
1175
1176 // Adding operators to calculate the left hand side
1177 feLhs.getOpPtrVector().clear();
1179 lagrange_multipliers_field_name, active_variables_ptr));
1181 material_field_name, active_variables_ptr));
1182 feLhs.getOpPtrVector().push_back(new OpJacobian(
1183 tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1184 jacobian_ptr, true, aLpha));
1186 lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1187
1189 }
1190};
1191
1192#endif // WITH_ADOL_C
1193
1194#endif // __SURFACE_SLIDING_CONSTRAINS_HPP__
static Number< 2 > N2
static Index< 'p', 3 > p
static Number< 1 > N1
static Number< 0 > N0
ForcesAndSourcesCore::UserDataOperator UserDataOperator
EntitiesFieldData::EntData EntData
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ 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()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'n', SPACE_DIM > n
constexpr double lambda
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double tol
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:65
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:95
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
double w(const double g, const double t)
const double r
rate factor
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< MatrixDouble > jacobianPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VectorDouble > resultsPtr
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)
MoFEMErrorCode setOperators(int tag, Range edges, Range faces, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
boost::shared_ptr< MyEdgeFE > feRhsPtr
EdgeSlidingConstrains(MoFEM::Interface &m_field)
boost::shared_ptr< MyEdgeFE > feLhsPtr
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, 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)
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)
boost::shared_ptr< MatrixDouble > jacobianPtr
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
boost::shared_ptr< VectorDouble > activeVariablesPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpGetActiveDofsLambda(const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr)
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
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
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 > activeVariablesPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< VectorDouble > resultsPtr
boost::shared_ptr< MatrixDouble > jacobianPtr
Shape preserving constrains, i.e. nodes sliding on body surface.
virtual ~SurfaceSlidingConstrains()=default
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)
boost::shared_ptr< MyTriangleFE > feRhsPtr
DriverElementOrientation & crackFrontOrientation