v0.14.0
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  */
16 struct GenericSliding {
17 
18  GenericSliding() = default;
19  virtual ~GenericSliding() = default;
20 
21  struct OpGetActiveDofsLambda
23  boost::shared_ptr<VectorDouble> activeVariablesPtr;
24  OpGetActiveDofsLambda(const std::string field_name,
25  boost::shared_ptr<VectorDouble> &active_variables_ptr)
28  activeVariablesPtr(active_variables_ptr) {}
29  MoFEMErrorCode doWork(int side, EntityType type,
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,
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,
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>
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  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  */
334  struct DriverElementOrientation {
335 
336  int elementOrientation;
337 
338  virtual MoFEMErrorCode
340  const FEMethod *fe_method_ptr) {
342  elementOrientation = 1;
344  }
346  const EntityHandle face) {
348  elementOrientation = 1;
350  }
351  };
352 
354 
355  struct MyTriangleFE : public MoFEM::FaceElementForcesAndSourcesCore {
356 
357  Mat B;
358  Vec F;
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 
384  MyTriangleFE &feRhs;
386  MyTriangleFE &feLhs;
388 
389  DriverElementOrientation &crackFrontOrientation;
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 
417  struct OpJacobian
419 
420  const int tAg;
421  boost::shared_ptr<VectorDouble> activeVariablesPtr;
422  boost::shared_ptr<VectorDouble> resultsPtr;
423  boost::shared_ptr<MatrixDouble> jacobianPtr;
424  DriverElementOrientation &oRientation;
425  bool evaluateJacobian;
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 
442  MoFEMErrorCode doWork(int side, EntityType type,
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>();
477  FTensor::Tensor1<adouble, 3> t_position;
478  FTensor::Tensor1<adouble, 3> t_position_ksi;
479  FTensor::Tensor1<adouble, 3> t_position_eta;
480  adouble lambda;
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));
601  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
602  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
603  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
604  material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
605  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
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));
632  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
633  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
634 
636  }
637 };
638 
639 struct EdgeSlidingConstrains : public GenericSliding {
640 
641  struct CalculateEdgeBase {
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 
876  struct MyEdgeFE : public MoFEM::EdgeElementForcesAndSourcesCore {
877 
878  Mat B;
879  Vec F;
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 
905  MyEdgeFE &feRhs;
906  MyEdgeFE &getLoopFeRhs() { return feRhs; }
907  MyEdgeFE &feLhs;
908  MyEdgeFE &getLoopFeLhs() { return feLhs; }
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 
933  struct OpJacobian
935 
936  const int tAg;
937  boost::shared_ptr<VectorDouble> activeVariablesPtr;
938  boost::shared_ptr<VectorDouble> resultsPtr;
939  boost::shared_ptr<MatrixDouble> jacobianPtr;
940  bool evaluateJacobian;
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 
955  MoFEMErrorCode doWork(int side, EntityType type,
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();
1011  FTensor::Tensor1<adouble, 3> t_position;
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 
1092  MoFEMErrorCode setOperators(int tag, Range edges, Range faces,
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();
1120  feRhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
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();
1134  feLhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
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));
1141  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1142  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1143  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1144  material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
1145  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1146  material_field_name, material_field_name, jacobian_ptr));
1147 
1149  }
1150 
1152  setOperatorsConstrainOnly(int tag, Range edges, Range faces,
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();
1178  feLhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
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));
1185  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
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__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::Types::VectorDouble6
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:95
SurfaceSlidingConstrains::OpJacobian::resultsPtr
boost::shared_ptr< VectorDouble > resultsPtr
Definition: SurfaceSlidingConstrains.hpp:422
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
EdgeSlidingConstrains::setOperatorsConstrainOnly
MoFEMErrorCode setOperatorsConstrainOnly(int tag, Range edges, Range faces, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
Definition: SurfaceSlidingConstrains.hpp:1152
SurfaceSlidingConstrains::DriverElementOrientation::getElementOrientation
virtual MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field, const FEMethod *fe_method_ptr)
Definition: SurfaceSlidingConstrains.hpp:339
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1644
EdgeSlidingConstrains::setOperators
MoFEMErrorCode setOperators(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name, const double *alpha=nullptr)
Definition: SurfaceSlidingConstrains.hpp:1103
EdgeSlidingConstrains::OpJacobian::activeVariablesPtr
boost::shared_ptr< VectorDouble > activeVariablesPtr
Definition: SurfaceSlidingConstrains.hpp:937
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
EdgeElementForcesAndSourcesCore
FTensor::Tensor1< adouble, 3 >
EntityHandle
FTensor::cross
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
EdgeSlidingConstrains::mField
MoFEM::Interface & mField
Definition: SurfaceSlidingConstrains.hpp:874
SurfaceSlidingConstrains::OpJacobian::OpJacobian
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)
Definition: SurfaceSlidingConstrains.hpp:429
EdgeSlidingConstrains::MyEdgeFE::B
Mat B
Definition: SurfaceSlidingConstrains.hpp:878
EdgeSlidingConstrains::getOptions
MoFEMErrorCode getOptions()
Definition: SurfaceSlidingConstrains.hpp:912
EdgeSlidingConstrains::EdgeSlidingConstrains
EdgeSlidingConstrains(MoFEM::Interface &m_field)
Definition: SurfaceSlidingConstrains.hpp:925
GenericSliding::OpGetActiveDofsPositions::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfaceSlidingConstrains.hpp:50
SurfaceSlidingConstrains::feRhsPtr
boost::shared_ptr< MyTriangleFE > feRhsPtr
Definition: SurfaceSlidingConstrains.hpp:382
MoFEM::ForcesAndSourcesCore::UserDataOperator::getSNESf
Vec getSNESf() const
Definition: ForcesAndSourcesCore.hpp:1112
EdgeSlidingConstrains::CalculateEdgeBase::setTags
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)
Definition: SurfaceSlidingConstrains.hpp:716
SurfaceSlidingConstrains::OpJacobian::tAg
const int tAg
Definition: SurfaceSlidingConstrains.hpp:420
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
EdgeSlidingConstrains::MyEdgeFE::getRule
int getRule(int order)
Definition: SurfaceSlidingConstrains.hpp:884
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::SnesMethod::snes_B
Mat & snes_B
preconditioner of jacobian matrix
Definition: LoopMethods.hpp:124
GenericSliding::OpGetActiveDofsLambda::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfaceSlidingConstrains.hpp:29
GenericSliding::OpAssembleRhs::resultsPtr
boost::shared_ptr< VectorDouble > resultsPtr
Definition: SurfaceSlidingConstrains.hpp:64
SurfaceSlidingConstrains::feRhs
MyTriangleFE & feRhs
Definition: SurfaceSlidingConstrains.hpp:384
SurfaceSlidingConstrains::mField
MoFEM::Interface & mField
Definition: SurfaceSlidingConstrains.hpp:353
GenericSliding
Implementation of surface sliding constrains.
Definition: SurfaceSlidingConstrains.hpp:16
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator::getEdgeFE
const EdgeElementForcesAndSourcesCore * getEdgeFE()
get pointer to this finite element
Definition: EdgeElementForcesAndSourcesCore.hpp:210
GenericSliding::OpGetActiveDofsPositions::activeVariablesPtr
boost::shared_ptr< VectorDouble > activeVariablesPtr
Definition: SurfaceSlidingConstrains.hpp:43
EdgeSlidingConstrains::getLoopFeRhs
MyEdgeFE & getLoopFeRhs()
Definition: SurfaceSlidingConstrains.hpp:906
EdgeSlidingConstrains::OpJacobian
Definition: SurfaceSlidingConstrains.hpp:933
SurfaceSlidingConstrains::MyTriangleFE::F
Vec F
Definition: SurfaceSlidingConstrains.hpp:358
SurfaceSlidingConstrains::OpJacobian::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfaceSlidingConstrains.hpp:442
MoFEM::ForcesAndSourcesCore::preProcess
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: ForcesAndSourcesCore.cpp:1993
SurfaceSlidingConstrains::aLpha
double aLpha
Definition: SurfaceSlidingConstrains.hpp:391
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition: ForcesAndSourcesCore.hpp:569
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
sdf.r
int r
Definition: sdf.py:8
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
SurfaceSlidingConstrains::getOptions
MoFEMErrorCode getOptions()
Definition: SurfaceSlidingConstrains.hpp:392
GenericSliding::OpAssembleRhs::OpAssembleRhs
OpAssembleRhs(const std::string field_name, boost::shared_ptr< VectorDouble > &results_ptr)
Definition: SurfaceSlidingConstrains.hpp:66
SurfaceSlidingConstrains::setOperators
MoFEMErrorCode setOperators(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name, const double *alpha=nullptr)
Definition: SurfaceSlidingConstrains.hpp:563
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1236
GenericSliding::OpGetActiveDofsLambda::OpGetActiveDofsLambda
OpGetActiveDofsLambda(const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr)
Definition: SurfaceSlidingConstrains.hpp:24
MoFEM::ForcesAndSourcesCore::UserDataOperator::ForcesAndSourcesCore
friend class ForcesAndSourcesCore
Definition: ForcesAndSourcesCore.hpp:990
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEName
std::string getFEName() const
Get name of the element.
Definition: ForcesAndSourcesCore.hpp:1064
SurfaceSlidingConstrains::OpJacobian
Definition: SurfaceSlidingConstrains.hpp:417
EdgeSlidingConstrains::MyEdgeFE::F
Vec F
Definition: SurfaceSlidingConstrains.hpp:879
FEMethod
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
SurfaceSlidingConstrains::OpJacobian::oRientation
DriverElementOrientation & oRientation
Definition: SurfaceSlidingConstrains.hpp:424
GenericSliding::OpGetActiveDofsPositions::OpGetActiveDofsPositions
OpGetActiveDofsPositions(const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr)
Definition: SurfaceSlidingConstrains.hpp:44
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
SurfaceSlidingConstrains::~SurfaceSlidingConstrains
virtual ~SurfaceSlidingConstrains()=default
EdgeSlidingConstrains::CalculateEdgeBase::saveEdges
static MoFEMErrorCode saveEdges(moab::Interface &moab, std::string name, Range edges, Range *faces=nullptr)
Definition: SurfaceSlidingConstrains.hpp:857
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SurfaceSlidingConstrains::getLoopFeLhs
MyTriangleFE & getLoopFeLhs()
Definition: SurfaceSlidingConstrains.hpp:387
EdgeSlidingConstrains::aLpha
double aLpha
Definition: SurfaceSlidingConstrains.hpp:910
SurfaceSlidingConstrains::feLhs
MyTriangleFE & feLhs
Definition: SurfaceSlidingConstrains.hpp:386
SurfaceSlidingConstrains::getLoopFeRhs
MyTriangleFE & getLoopFeRhs()
Definition: SurfaceSlidingConstrains.hpp:385
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
SurfaceSlidingConstrains::OpJacobian::jacobianPtr
boost::shared_ptr< MatrixDouble > jacobianPtr
Definition: SurfaceSlidingConstrains.hpp:423
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
SurfaceSlidingConstrains::OpJacobian::aLpha
const double & aLpha
Definition: SurfaceSlidingConstrains.hpp:427
EdgeSlidingConstrains::setOperatorsConstrainOnly
MoFEMErrorCode setOperatorsConstrainOnly(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
Definition: SurfaceSlidingConstrains.hpp:1165
convert.type
type
Definition: convert.py:64
EdgeSlidingConstrains::MyEdgeFE::MyEdgeFE
MyEdgeFE(MoFEM::Interface &m_field)
Definition: SurfaceSlidingConstrains.hpp:881
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPCOL
@ OPCOL
operator doWork function is executed on FE columns
Definition: ForcesAndSourcesCore.hpp:568
GenericSliding::OpGetActiveDofsLambda
Definition: SurfaceSlidingConstrains.hpp:21
SurfaceSlidingConstrains::DriverElementOrientation::getElementOrientation
virtual MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field, const EntityHandle face)
Definition: SurfaceSlidingConstrains.hpp:345
SurfaceSlidingConstrains::SurfaceSlidingConstrains
SurfaceSlidingConstrains(MoFEM::Interface &m_field, DriverElementOrientation &orientation)
Definition: SurfaceSlidingConstrains.hpp:405
SurfaceSlidingConstrains::feLhsPtr
boost::shared_ptr< MyTriangleFE > feLhsPtr
Definition: SurfaceSlidingConstrains.hpp:382
GenericSliding::~GenericSliding
virtual ~GenericSliding()=default
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
EdgeSlidingConstrains::setOperators
MoFEMErrorCode setOperators(int tag, Range edges, Range faces, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
Definition: SurfaceSlidingConstrains.hpp:1092
SurfaceSlidingConstrains
Shape preserving constrains, i.e. nodes sliding on body surface.
Definition: SurfaceSlidingConstrains.hpp:323
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
SurfaceSlidingConstrains::MyTriangleFE::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: SurfaceSlidingConstrains.hpp:365
EdgeSlidingConstrains::OpJacobian::tAg
const int tAg
Definition: SurfaceSlidingConstrains.hpp:936
SurfaceSlidingConstrains::OpJacobian::activeVariablesPtr
boost::shared_ptr< VectorDouble > activeVariablesPtr
Definition: SurfaceSlidingConstrains.hpp:421
EdgeSlidingConstrains::feRhsPtr
boost::shared_ptr< MyEdgeFE > feRhsPtr
Definition: SurfaceSlidingConstrains.hpp:903
EdgeSlidingConstrains::OpJacobian::resultsPtr
boost::shared_ptr< VectorDouble > resultsPtr
Definition: SurfaceSlidingConstrains.hpp:938
EdgeSlidingConstrains::feRhs
MyEdgeFE & feRhs
Definition: SurfaceSlidingConstrains.hpp:905
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
GenericSliding::OpAssembleLhs::jacobianPtr
boost::shared_ptr< MatrixDouble > jacobianPtr
Definition: SurfaceSlidingConstrains.hpp:99
EdgeSlidingConstrains::CalculateEdgeBase::createTag
static MoFEMErrorCode createTag(moab::Interface &moab, Tag &th0, Tag &th1, Tag &th2, Tag &th3)
Definition: SurfaceSlidingConstrains.hpp:643
GenericSliding::OpAssembleRhs
Definition: SurfaceSlidingConstrains.hpp:62
EdgeSlidingConstrains::OpJacobian::evaluateJacobian
bool evaluateJacobian
Definition: SurfaceSlidingConstrains.hpp:940
EdgeSlidingConstrains::MyEdgeFE::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: SurfaceSlidingConstrains.hpp:886
MoFEM::ForcesAndSourcesCore::UserDataOperator::getSNESB
Mat getSNESB() const
Definition: ForcesAndSourcesCore.hpp:1136
FaceElementForcesAndSourcesCore
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:1042
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator
default operator for EDGE element
Definition: EdgeElementForcesAndSourcesCore.hpp:68
GenericSliding::OpAssembleLhs
Definition: SurfaceSlidingConstrains.hpp:97
convert.n
n
Definition: convert.py:82
EdgeSlidingConstrains::OpJacobian::OpJacobian
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)
Definition: SurfaceSlidingConstrains.hpp:944
SurfaceSlidingConstrains::setOperatorsConstrainOnly
MoFEMErrorCode setOperatorsConstrainOnly(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
Definition: SurfaceSlidingConstrains.hpp:612
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
EdgeSlidingConstrains::MyEdgeFE
Definition: SurfaceSlidingConstrains.hpp:876
adouble
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:122
FTensor::dd
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
SurfaceSlidingConstrains::crackFrontOrientation
DriverElementOrientation & crackFrontOrientation
Definition: SurfaceSlidingConstrains.hpp:389
FTensor::Tensor0
Definition: Tensor0.hpp:16
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
Definition: ForcesAndSourcesCore.hpp:1004
EdgeSlidingConstrains::CalculateEdgeBase::numberSurfaces
static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges, Range tris)
Definition: SurfaceSlidingConstrains.hpp:662
GenericSliding::OpAssembleRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfaceSlidingConstrains.hpp:72
MoFEM::Types::VectorInt
UBlasVector< int > VectorInt
Definition: Types.hpp:67
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
GenericSliding::OpAssembleLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: SurfaceSlidingConstrains.hpp:110
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
SurfaceSlidingConstrains::MyTriangleFE
Definition: SurfaceSlidingConstrains.hpp:355
SurfaceSlidingConstrains::MyTriangleFE::getRule
int getRule(int order)
Definition: SurfaceSlidingConstrains.hpp:363
EdgeSlidingConstrains::OpJacobian::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: SurfaceSlidingConstrains.hpp:955
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
GenericSliding::OpGetActiveDofsPositions
Definition: SurfaceSlidingConstrains.hpp:41
EdgeSlidingConstrains::OpJacobian::aLpha
const double & aLpha
Definition: SurfaceSlidingConstrains.hpp:942
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
GenericSliding::OpGetActiveDofsLambda::activeVariablesPtr
boost::shared_ptr< VectorDouble > activeVariablesPtr
Definition: SurfaceSlidingConstrains.hpp:23
GenericSliding::GenericSliding
GenericSliding()=default
EdgeSlidingConstrains::feLhsPtr
boost::shared_ptr< MyEdgeFE > feLhsPtr
Definition: SurfaceSlidingConstrains.hpp:903
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
EdgeSlidingConstrains::feLhs
MyEdgeFE & feLhs
Definition: SurfaceSlidingConstrains.hpp:907
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:82
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getFaceFE
FaceElementForcesAndSourcesCore * getFaceFE()
return pointer to Generic Triangle Finite Element object
Definition: FaceElementForcesAndSourcesCore.hpp:337
SurfaceSlidingConstrains::OpJacobian::evaluateJacobian
bool evaluateJacobian
Definition: SurfaceSlidingConstrains.hpp:425
EdgeSlidingConstrains::getLoopFeLhs
MyEdgeFE & getLoopFeLhs()
Definition: SurfaceSlidingConstrains.hpp:908
SurfaceSlidingConstrains::DriverElementOrientation
Class implemented by user to detect face orientation.
Definition: SurfaceSlidingConstrains.hpp:334
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor1CoordsAtGaussPts
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
Definition: ForcesAndSourcesCore.hpp:1269
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
SurfaceSlidingConstrains::DriverElementOrientation::elementOrientation
int elementOrientation
Definition: SurfaceSlidingConstrains.hpp:336
tol
double tol
Definition: mesh_smoothing.cpp:26
GenericSliding::OpAssembleLhs::OpAssembleLhs
OpAssembleLhs(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > &jacobian_ptr)
Definition: SurfaceSlidingConstrains.hpp:101
EdgeSlidingConstrains
Definition: SurfaceSlidingConstrains.hpp:639
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
SurfaceSlidingConstrains::MyTriangleFE::B
Mat B
Definition: SurfaceSlidingConstrains.hpp:357
EdgeSlidingConstrains::OpJacobian::jacobianPtr
boost::shared_ptr< MatrixDouble > jacobianPtr
Definition: SurfaceSlidingConstrains.hpp:939
SurfaceSlidingConstrains::MyTriangleFE::MyTriangleFE
MyTriangleFE(MoFEM::Interface &m_field)
Definition: SurfaceSlidingConstrains.hpp:360