v0.10.0
SurfaceSlidingConstrains.hpp
Go to the documentation of this file.
1 /** \file SurfaceSlidingConstrains.hpp
2  * \brief Implementing surface sliding constrains
3  */
4 
5 /* This file is part of MoFEM.
6  * MoFEM is free software: you can redistribute it and/or modify it under
7  * the terms of the GNU Lesser General Public License as published by the
8  * Free Software Foundation, either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14  * License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
18 
19 #ifndef __SURFACE_SLIDING_CONSTRAINS_HPP__
20 #define __SURFACE_SLIDING_CONSTRAINS_HPP__
21 
22 #ifdef WITH_ADOL_C
23 
24 /**
25  * @brief Implementation of surface sliding constrains
26  *
27  */
29 
30  GenericSliding() = default;
31  virtual ~GenericSliding() = default;
32 
35  boost::shared_ptr<VectorDouble> activeVariablesPtr;
36  OpGetActiveDofsLambda(const std::string field_name,
37  boost::shared_ptr<VectorDouble> &active_variables_ptr)
39  field_name, UserDataOperator::OPCOL),
40  activeVariablesPtr(active_variables_ptr) {}
41  MoFEMErrorCode doWork(int side, EntityType type,
44  if (type == MBVERTEX) {
45  for (unsigned int dd = 0; dd != data.getFieldData().size(); ++dd)
46  (*activeVariablesPtr)[dd] = data.getFieldData()[dd];
47  }
49  }
50  };
51 
52  template <int SizeLambda>
55  boost::shared_ptr<VectorDouble> activeVariablesPtr;
57  const std::string field_name,
58  boost::shared_ptr<VectorDouble> &active_variables_ptr)
60  field_name, UserDataOperator::OPCOL),
61  activeVariablesPtr(active_variables_ptr) {}
62  MoFEMErrorCode doWork(int side, EntityType type,
65  if (type == MBVERTEX) {
66  for (unsigned int dd = 0; dd != data.getFieldData().size(); ++dd)
67  (*activeVariablesPtr)[SizeLambda + dd] = data.getFieldData()[dd];
68  }
70  }
71  };
72 
73  template <int SizeLambda, int SizePositions>
75 
76  boost::shared_ptr<VectorDouble> resultsPtr;
77 
78  OpAssembleRhs(const std::string field_name,
79  boost::shared_ptr<VectorDouble> &results_ptr)
81  field_name, UserDataOperator::OPROW),
82  resultsPtr(results_ptr) {}
83 
84  MoFEMErrorCode doWork(int side, EntityType type,
87  if (type != MBVERTEX)
89  VectorInt &indices = data.getIndices();
90  int shift = 0;
91  if (indices.empty())
93  else if (indices.size() == SizeLambda)
94  shift = 0;
95  else if (indices.size() == SizePositions)
96  shift = SizeLambda;
97  else
98  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
99  "Element %s: Data inconsistency nb of indices %d",
100  getFEName().c_str(), indices.size());
101 
102  CHKERR VecSetOption(getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
103  CHKERR VecSetValues(getSNESf(), indices.size(), &indices[0],
104  &(*resultsPtr)[shift], ADD_VALUES);
106  }
107  };
108  template <int SizeLambda, int SizePositions>
110 
111  boost::shared_ptr<MatrixDouble> jacobianPtr;
112 
113  OpAssembleLhs(const std::string field_name_row,
114  const std::string field_name_col,
115  boost::shared_ptr<MatrixDouble> &jacobian_ptr)
117  field_name_row, field_name_col, UserDataOperator::OPROWCOL),
118  jacobianPtr(jacobian_ptr) {
119  sYmm = false;
120  }
121 
122  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
123  EntityType col_type,
127  if (row_type != MBVERTEX)
129  if (col_type != MBVERTEX)
131  VectorInt &row_indices = row_data.getIndices();
132  VectorInt &col_indices = col_data.getIndices();
133  if (row_indices.empty() || col_indices.empty())
135 
136  int shift_row = 0;
137  if (row_indices.size() == SizeLambda)
138  shift_row = 0;
139  else if (row_indices.size() == SizePositions)
140  shift_row = SizeLambda;
141  else
142  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
143  "Data inconsistency");
144 
145  int shift_col = 0;
146  if (col_indices.size() == SizeLambda)
147  shift_col = 0;
148  else if (col_indices.size() == SizePositions)
149  shift_col = SizeLambda;
150  else
151  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
152  "Data inconsistency");
153 
154  MatrixDouble jac(row_indices.size(), col_indices.size());
155  for (unsigned int rr = 0; rr != row_indices.size(); ++rr) {
156  for (unsigned int cc = 0; cc != col_indices.size(); ++cc) {
157  jac(rr, cc) = (*jacobianPtr)(shift_row + rr, shift_col + cc);
158  if (jac(rr, cc) != jac(rr, cc))
159  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
160  "Jacobian assemble not a number jac(rr,cc) = %3.4f",
161  jac(rr, cc));
162  }
163  }
164  CHKERR MatSetValues(getSNESB(), row_indices.size(), &row_indices[0],
165  col_indices.size(), &col_indices[0], &jac(0, 0),
166  ADD_VALUES);
168  }
169  };
170 };
171 
172 /** \brief Shape preserving constrains, i.e. nodes sliding on body surface.
173 
174  Derivation and implementation of constrains preserving body surface,
175  i.e. body shape and volume.
176 
177  The idea starts form observation that body shape can be globally characterized
178  by constant calculated as volume over its area
179  \f[
180  \frac{V}{A} = C
181  \f]
182  Above equation expressed in integral form is
183  \f[
184  \int_\Omega \textrm{d}V = C \int_\Gamma \textrm{d}S
185  \f]
186  where notting that,
187  \f[
188  \frac{1}{3}
189  \int_\Omega \textrm{div}[\mathbf{X}] \textrm{d}\Omega
190  =
191  C \int_\Gamma \textrm{d}S
192  \f]
193  and applying Gauss theorem we get
194  \f[
195  \int_\Gamma
196  \mathbf{X}\cdot \frac{\mathbf{N}}{\|\mathbf{N}\|}
197  \textrm{d}\Gamma
198  =
199  3C \int_\Gamma \textrm{d}S.
200  \f]
201  Drooping integrals on both sides, and linearizing equation, we get
202  \f[
203  \frac{\mathbf{N}}{\|\mathbf{N}\|} \cdot \delta \mathbf{X}
204  =
205  3C - \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot \mathbf{X}
206  \f]
207  where \f$\delta \mathbf{X}\f$ is displacement sub-inctrement. Above equation
208  is a constrain if satisfied in body shape and volume is conserved. Final form
209  of constrain equation is \f[ \mathcal{r} =
210  \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot \mathbf{X}
211  -
212  \frac{\mathbf{N_0}}{\|\mathbf{N_0}\|}\cdot \mathbf{X}_0 =
213  \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot (\mathbf{X}-\mathbf{X}_0)
214  \f]
215 
216  In the spirit of finite element method the constrain equation is multiplied
217  by shape functions and enforce using Lagrange multiplier method
218  \f[
219  \int_\Gamma \mathbf{N}^\mathsf{T}_\lambda
220  \left(
221  \frac{\mathbf{N}}{\|\mathbf{N}\|}\mathbf{N}_\mathbf{X}\cdot
222  (\overline{\mathbf{X}}-\overline{\mathbf{X}}_0)
223  \right)
224  \|\mathbf{N}\|
225  \textrm{d}\Gamma
226  =
227  \mathbf{0}.
228  \f]
229  Above equation is nonlinear, applying to it Taylor expansion, we can get form
230  which can be used with Newton interactive method \f[ \begin{split}
231  &\int_\Gamma \mathbf{N}^\mathsf{T}_\lambda
232  \left\{
233  \mathbf{N}\mathbf{N}_\mathbf{X}
234  +
235  \left(\mathbf{X}-\mathbf{X}_0\right) \cdot
236  \left(
237  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
238  -
239  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
240  \right)
241  \right\}
242  \textrm{d}\Gamma
243  \cdot
244  \delta
245  \overline{\mathbf{X}}\\
246  =
247  &\int_\Gamma \mathbf{N}^\mathsf{T}_\lambda
248  \mathbf{N}\cdot(\mathbf{X}-\mathbf{X}_0)
249  \textrm{d}\Gamma
250  \end{split}.
251  \f]
252  Equation expressing forces on shape as result of constrains, as result
253  Lagrange multiplier method have following form \f[ \begin{split}
254  &\int_\Gamma
255  \mathbf{N}^\mathsf{T}_\mathbf{X} \cdot \mathbf{N}
256  \mathbf{N}_\lambda
257  \textrm{d}\Gamma
258  \cdot
259  \delta\overline{\lambda}\\
260  +
261  &\int_\Gamma
262  \lambda
263  \mathbf{N}^\mathsf{T}_\mathbf{X}
264  \left(
265  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
266  -
267  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
268  \right)
269  \textrm{d}\Gamma
270  \delta\overline{\mathbf{X}}\\
271  =
272  &\int_\Gamma
273  \lambda
274  \mathbf{N}^\mathsf{T}_\mathbf{X} \cdot \mathbf{N}
275  \textrm{d}\Gamma
276  \end{split}
277  \f]
278 
279  Above equations are assembled into global system of equations as following
280  \f[
281  \left[
282  \begin{array}{cc}
283  \mathbf{K} + \mathbf{B} & \mathbf{C}^\mathsf{T} \\
284  \mathbf{C} + \mathbf{A} & 0
285  \end{array}
286  \right]
287  \left\{
288  \begin{array}{c}
289  \delta \overline{\mathbf{X}} \\
290  \delta \overline{\lambda}
291  \end{array}
292  \right\}=
293  \left[
294  \begin{array}{c}
295  \mathbf{f} - \mathbf{C}^\mathsf{T}\overline{\lambda} \\
296  \overline{\mathbf{r}}
297  \end{array}
298  \right]
299  \f]
300  where
301  \f[
302  \mathbf{C}=
303  \int_\Gamma
304  \mathbf{N}_\lambda^\mathsf{T}
305  \mathbf{N} \cdot
306  \mathbf{N}_\mathbf{X}
307  \textrm{d}\Gamma,
308  \f]
309  \f[
310  \mathbf{B}=
311  \int_\Gamma
312  \lambda
313  (\mathbf{X}-\mathbf{X}_0)\cdot
314  \left(
315  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
316  -
317  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
318  \right)
319  \textrm{d}\Gamma
320  \f]
321  and
322  \f[
323  \mathbf{A}=
324  \int_\Gamma
325  \mathbf{N}^\mathsf{T}_\lambda
326  \left(\mathbf{X}-\mathbf{X}_0\right) \cdot
327  \left(
328  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
329  -
330  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
331  \right).
332  \f]
333 
334 */
336 
337  /** \brief Class implemented by user to detect face orientation
338 
339  If mesh generated is with surface mesher, usually you don't have to do
340  nothing, all elements on the surface have consistent orientation. In case of
341  internal faces or if you do something with mesh connectivity which breaks
342  orientation on the face, you have to implement method which will set
343  orientation to face.
344 
345  */
347 
349 
350  virtual MoFEMErrorCode
352  const FEMethod *fe_method_ptr) {
354  elementOrientation = 1;
356  }
357  };
358 
360 
362 
363  Mat B;
364  Vec F;
365 
367  : MoFEM::FaceElementForcesAndSourcesCore(m_field), B(PETSC_NULL),
368  F(PETSC_NULL) {}
369  int getRule(int order) { return 2 * order; };
370 
373 
374  CHKERR MoFEM::FaceElementForcesAndSourcesCore::preProcess();
375 
376  if (B != PETSC_NULL) {
377  snes_B = B;
378  }
379 
380  if (F != PETSC_NULL) {
381  snes_f = F;
382  }
383 
385  }
386  };
387 
388  boost::shared_ptr<MyTriangleFE> feRhsPtr, feLhsPtr;
389 
394 
396 
397  double aLpha;
400  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
401  "Get surface sliding constrains element scaling",
402  "none");
403  CHKERRQ(ierr);
404  CHKERR PetscOptionsScalar("-surface_sliding_alpha", "scaling parameter", "",
405  aLpha, &aLpha, PETSC_NULL);
406  ierr = PetscOptionsEnd();
407  CHKERRQ(ierr);
409  }
410 
412  DriverElementOrientation &orientation)
413  : mField(m_field), feRhsPtr(new MyTriangleFE(m_field)),
414  feLhsPtr(new MyTriangleFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr),
415  crackFrontOrientation(orientation), aLpha(1) {
416 
417  ierr = getOptions();
418  CHKERRABORT(PETSC_COMM_WORLD, ierr);
419  }
420 
421  virtual ~SurfaceSlidingConstrains() = default;
422 
423  struct OpJacobian
425 
426  const int tAg;
427  boost::shared_ptr<VectorDouble> activeVariablesPtr;
428  boost::shared_ptr<VectorDouble> resultsPtr;
429  boost::shared_ptr<MatrixDouble> jacobianPtr;
432 
433  const double &aLpha;
434 
435  OpJacobian(int tag, const std::string field_name,
436  boost::shared_ptr<VectorDouble> &active_variables_ptr,
437  boost::shared_ptr<VectorDouble> &results_ptr,
438  boost::shared_ptr<MatrixDouble> &jacobian_ptr,
439  DriverElementOrientation &orientation, bool evaluate_jacobian,
440  const double &alpha)
442  field_name, UserDataOperator::OPCOL),
443  tAg(tag), activeVariablesPtr(active_variables_ptr),
444  resultsPtr(results_ptr), jacobianPtr(jacobian_ptr),
445  oRientation(orientation), evaluateJacobian(evaluate_jacobian),
446  aLpha(alpha) {}
447 
448  MoFEMErrorCode doWork(int side, EntityType type,
451  if (type != MBVERTEX)
453 
454  VectorInt &indices = data.getIndices();
455 
456  trace_on(tAg);
457 
458  ublas::vector<adouble> lambda_dofs(3);
459  for (int dd = 0; dd != 3; ++dd) {
460  lambda_dofs[dd] <<= (*activeVariablesPtr)[dd];
461  }
462  ublas::vector<adouble> position_dofs(9);
463  for (int dd = 0; dd != 9; ++dd) {
464  position_dofs[dd] <<= (*activeVariablesPtr)[3 + dd];
465  }
466 
473 
475  getFEMethod());
477 
478  int nb_gauss_pts = data.getN().size1();
479  int nb_base_functions = data.getN().size2();
480  auto t_base1 = data.getFTensor0N();
481  auto t_base2 = data.getFTensor0N();
482  auto t_diff_base = data.getFTensor1DiffN<2>();
483  FTensor::Tensor1<adouble, 3> t_position;
484  FTensor::Tensor1<adouble, 3> t_position_ksi;
485  FTensor::Tensor1<adouble, 3> t_position_eta;
486  adouble lambda;
489 
490  ublas::vector<adouble> c_vec(3);
491  ublas::vector<adouble> f_vec(9);
492  c_vec.clear();
493  f_vec.clear();
494 
495  auto t_coord_ref = getFTensor1CoordsAtGaussPts();
496 
497  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
498 
499  FTensor::Tensor1<adouble *, 3> t_position_dofs(
500  &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
501  FTensor::Tensor0<adouble *> t_lambda_dof(&lambda_dofs[0]);
502 
503  t_position(i) = 0;
504  t_position_ksi(i) = 0;
505  t_position_eta(i) = 0;
506  lambda = 0;
507 
508  for (int bb = 0; bb != nb_base_functions; ++bb) {
509  t_position(i) += t_base1 * t_position_dofs(i);
510  t_position_ksi(i) += t_diff_base(N0) * t_position_dofs(i);
511  t_position_eta(i) += t_diff_base(N1) * t_position_dofs(i);
512  lambda += t_base1 * t_lambda_dof;
513  ++t_base1;
514  ++t_position_dofs;
515  ++t_lambda_dof;
516  ++t_diff_base;
517  }
518 
519  t_delta(i) = t_position(i) - t_coord_ref(i);
520  t_normal(k) = FTensor::cross(t_position_ksi(i), t_position_eta(j), k);
521 
522  double w = getGaussPts()(2, gg) * 0.5 * aLpha;
523  adouble val;
524  FTensor::Tensor0<adouble *> t_c(&c_vec[0]);
525  FTensor::Tensor1<adouble *, 3> t_f(&f_vec[0], &f_vec[1], &f_vec[2], 3);
526 
527  adouble area = sqrt(t_normal(i) * t_normal(i));
528 
529  for (int bb = 0; bb != nb_base_functions; ++bb) {
530  if (indices[bb] != -1) {
531  val = w * eo * t_base2;
532  t_c += val * t_normal(i) * t_delta(i);
533  val *= lambda;
534  t_f(i) += val * t_normal(i);
535  }
536  ++t_c;
537  ++t_f;
538  ++t_base2;
539  }
540 
541  ++t_coord_ref;
542  }
543 
544  for (int rr = 0; rr != 3; ++rr)
545  c_vec[rr] >>= (*resultsPtr)[rr];
546 
547  for (int rr = 0; rr != 9; ++rr)
548  f_vec(rr) >>= (*resultsPtr)[3 + rr];
549 
550 
551  trace_off();
552 
553  if (evaluateJacobian) {
554  double *jac_ptr[3 + 9];
555  for (int rr = 0; rr != 3 + 9; ++rr)
556  jac_ptr[rr] = &(*jacobianPtr)(rr, 0);
557 
558  // play recorder for jacobians
559  int r =
560  ::jacobian(tAg, 3 + 9, 3 + 9, &(*activeVariablesPtr)[0], jac_ptr);
561  if (r < 0)
562  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
563  "ADOL-C function evaluation with error");
564 
565  }
566 
568  }
569  };
570 
572  const std::string lagrange_multipliers_field_name,
573  const std::string material_field_name,
574  const double *alpha = nullptr) {
576 
577  if (alpha != nullptr) {
578  aLpha = *alpha;
579  }
580 
581  boost::shared_ptr<VectorDouble> active_variables_ptr(
582  new VectorDouble(3 + 9));
583  boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
584  boost::shared_ptr<MatrixDouble> jacobian_ptr(
585  new MatrixDouble(3 + 9, 3 + 9));
586 
587  feRhs.getOpPtrVector().clear();
589  lagrange_multipliers_field_name, active_variables_ptr));
591  material_field_name, active_variables_ptr));
592  feRhs.getOpPtrVector().push_back(new OpJacobian(
593  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
594  jacobian_ptr, crackFrontOrientation, false, aLpha));
595  feRhs.getOpPtrVector().push_back(
596  new OpAssembleRhs<3, 9>(lagrange_multipliers_field_name, results_ptr));
597  feRhs.getOpPtrVector().push_back(
598  new OpAssembleRhs<3, 9>(material_field_name, results_ptr));
599 
600  // Adding operators to calculate the left hand side
601  feLhs.getOpPtrVector().clear();
603  lagrange_multipliers_field_name, active_variables_ptr));
605  material_field_name, active_variables_ptr));
606  feLhs.getOpPtrVector().push_back(new OpJacobian(
607  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
608  jacobian_ptr, crackFrontOrientation, true, aLpha));
609  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
610  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
611  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
612  material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
613  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
614  material_field_name, material_field_name, jacobian_ptr));
615 
617  }
618 
621  const std::string lagrange_multipliers_field_name,
622  const std::string material_field_name) {
624 
625  boost::shared_ptr<VectorDouble> active_variables_ptr(
626  new VectorDouble(3 + 9));
627  boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
628  boost::shared_ptr<MatrixDouble> jacobian_ptr(
629  new MatrixDouble(3 + 9, 3 + 9));
630 
631  // Adding operators to calculate the left hand side
632  feLhs.getOpPtrVector().clear();
634  lagrange_multipliers_field_name, active_variables_ptr));
636  material_field_name, active_variables_ptr));
637  feLhs.getOpPtrVector().push_back(new OpJacobian(
638  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
639  jacobian_ptr, crackFrontOrientation, true, aLpha));
640  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
641  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
642 
644  }
645 };
646 
648 
650 
651  static MoFEMErrorCode createTag(moab::Interface &moab, Tag &th0, Tag &th1,
652  Tag &th2, Tag &th3) {
654  double def_val[] = {0, 0, 0};
655  CHKERR moab.tag_get_handle("EDGE_BASE0", 3, MB_TYPE_DOUBLE, th0,
656  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
657  CHKERR moab.tag_get_handle("EDGE_BASE1", 3, MB_TYPE_DOUBLE, th1,
658  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
659  int def_numb_val[] = {-1};
660  CHKERR moab.tag_get_handle("PATCH_NUMBER", 1, MB_TYPE_INTEGER, th2,
661  MB_TAG_CREAT | MB_TAG_SPARSE, def_numb_val);
662  int def_orientation_val[] = {1};
663  CHKERR moab.tag_get_handle("PATCH_ORIENTATION", 1, MB_TYPE_INTEGER, th3,
664  MB_TAG_CREAT | MB_TAG_SPARSE,
665  def_orientation_val);
666 
668  }
669 
670  static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges,
671  Range tris) {
673 
674  auto get_edges = [&](const Range &ents) {
675  Range edges;
676  CHKERR moab.get_adjacencies(ents, 1, false, edges,
677  moab::Interface::UNION);
678  return edges;
679  };
680 
681  auto get_face_adj = [&, get_edges](const Range &faces) {
682  Range adj_faces;
683  CHKERR moab.get_adjacencies(subtract(get_edges(faces), edges), 2, false,
684  adj_faces, moab::Interface::UNION);
685  return intersect(adj_faces, tris);
686  };
687 
688  auto get_patch = [&, get_face_adj](const EntityHandle face) {
689  Range patch_ents;
690  patch_ents.insert(face);
691  unsigned int nb0;
692  do {
693  nb0 = patch_ents.size();
694  patch_ents.merge(get_face_adj(patch_ents));
695  } while (nb0 != patch_ents.size());
696  return patch_ents;
697  };
698 
699  auto get_patches = [&]() {
700  std::vector<Range> patches;
701  while (!tris.empty()) {
702  patches.push_back(get_patch(tris[0]));
703  tris = subtract(tris, patches.back());
704  }
705  return patches;
706  };
707 
708  Tag th0, th1, th2, th3;
709  CHKERR createTag(moab, th0, th1, th2, th3);
710 
711  auto patches = get_patches();
712  int pp = 0;
713  for (auto patch : patches) {
714  // cerr << "pp: " << pp << endl;
715  // cerr << patch << endl;
716  std::vector<int> tags_vals(patch.size(), pp);
717  CHKERR moab.tag_set_data(th2, patch, &*tags_vals.begin());
718  ++pp;
719  }
720 
722  }
723 
724  static MoFEMErrorCode setTags(moab::Interface &moab, Range edges,
725  Range tris, bool number_pathes = true) {
727  Tag th0, th1, th2, th3;
728  CHKERR createTag(moab, th0, th1, th2, th3);
729  if (number_pathes) {
730  CHKERR numberSurfaces(moab, edges, tris);
731  }
732  for (auto edge : edges) {
733  Range adj_faces;
734  CHKERR moab.get_adjacencies(&edge, 1, 2, false, adj_faces);
735  adj_faces = intersect(adj_faces, tris);
736  if (adj_faces.size() != 2) {
737  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
738  "Should be 2 faces adjacent to edge but is %d",
739  adj_faces.size());
740  }
741  VectorDouble3 v[2] = {VectorDouble3(3), VectorDouble3(3)};
742  auto get_tensor_from_vec = [](VectorDouble3 &v) {
744  &v[0], &v[1], &v[2]);
745  };
746 
750 
751  std::array<double, 3> areas;
752  auto calculate_normals = [&, get_tensor_from_vec]() {
753  int ff = 0;
754  for (auto face : adj_faces) {
755  double &x = (v[ff])[0];
756  double &y = (v[ff])[1];
757  double &z = (v[ff])[2];
758  moab::Util::normal(&moab, face, x, y, z);
759  auto t_n = get_tensor_from_vec(v[ff]);
760  areas[ff] = sqrt(t_n(i) * t_n(i));
761  t_n(i) /= areas[ff];
762  int orientation;
763  CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
764  if (orientation == -1) {
765  t_n(i) *= -1;
766  }
767  ++ff;
768  }
769  };
770  calculate_normals();
771 
772  auto get_patch_number = [&]() {
773  std::vector<int> p = {0, 0};
774  CHKERR moab.tag_get_data(th2, adj_faces, &*p.begin());
775  return p;
776  };
777 
778  std::array<EntityHandle, 2> faces_handles = {adj_faces[0],
779  adj_faces[1]};
780  auto order_normals = [&, get_patch_number]() {
781  auto p = get_patch_number();
782  if (p[0] < p[1]) {
783  v[0].swap(v[1]);
784  faces_handles[0] = adj_faces[1];
785  faces_handles[1] = adj_faces[0];
786  }
787  };
788  order_normals();
789 
791  auto t_n0 = get_tensor_from_vec(v[0]);
792  auto t_n1 = get_tensor_from_vec(v[1]);
793  t_cross(k) = FTensor::cross(t_n0(i), t_n1(j), k);
794 
795  auto get_base_for_coplanar_faces = [&]() {
797 
798  std::vector<EntityHandle> face_conn;
799  CHKERR moab.get_connectivity(&faces_handles[1], 1, face_conn, true);
800  std::vector<EntityHandle> edge_conn;
801  CHKERR moab.get_connectivity(&edge, 1, edge_conn, true);
802  std::sort(face_conn.begin(), face_conn.end());
803  std::sort(edge_conn.begin(), edge_conn.end());
804  int n = 0;
805  for (; n != 2; ++n)
806  if (face_conn[n] != edge_conn[n])
807  break;
808  VectorDouble3 coords_face(3);
809  CHKERR moab.get_coords(&face_conn[n], 1,
810  &*coords_face.data().begin());
811  VectorDouble6 coords_edge(6);
812  CHKERR moab.get_coords(&edge_conn[0], 2,
813  &*coords_edge.data().begin());
814  auto t_edge_n0 = FTensor::Tensor1<double, 3>(
815  coords_edge[0], coords_edge[1], coords_edge[2]);
816  auto t_edge_n1 = FTensor::Tensor1<double, 3>(
817  coords_edge[3], coords_edge[4], coords_edge[5]);
818  auto t_face_n = get_tensor_from_vec(coords_face);
819  t_face_n(i) -= t_edge_n0(i);
821  t_delta(i) = t_edge_n1(i) - t_edge_n0(i);
822  t_delta(i) /= sqrt(t_delta(i) * t_delta(i));
823  t_n1(k) = FTensor::cross(t_n0(i), t_delta(j), k);
824  if (t_n1(i) * t_face_n(i) < 0)
825  t_n1(i) *= -1;
826 
828  };
829 
830  auto get_base_for_not_planar_faces = [&]() {
832  t_n1(k) = FTensor::cross(t_n0(i), t_cross(j), k);
833  t_n1(i) /= sqrt(t_n1(i) * t_n1(i));
835  };
836 
837  // This a case when faces adjacent to edge are coplanar !!!
838  constexpr double tol = 1e-6;
839  if ((t_cross(i) * t_cross(i)) < tol)
840  CHKERR get_base_for_coplanar_faces();
841  else
842  CHKERR get_base_for_not_planar_faces();
843 
844  VectorDouble3 &v0 = v[0];
845  VectorDouble3 &v1 = v[1];
846  CHKERR moab.tag_set_data(th0, &edge, 1, &v0[0]);
847  CHKERR moab.tag_set_data(th1, &edge, 1, &v1[0]);
848  }
850  }
851 
852  static MoFEMErrorCode saveEdges(moab::Interface &moab, std::string name,
853  Range edges, Range *faces = nullptr) {
855  EntityHandle meshset;
856  Tag ths[4];
857  CHKERR createTag(moab, ths[0], ths[1], ths[2], ths[3]);
858  CHKERR moab.create_meshset(MESHSET_SET, meshset);
859  CHKERR moab.add_entities(meshset, edges);
860  if (faces != nullptr) {
861  CHKERR moab.add_entities(meshset, *faces);
862  }
863  CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1, ths, 4);
864  CHKERR moab.delete_entities(&meshset, 1);
866  }
867  };
868 
870 
872 
873  Mat B;
874  Vec F;
875 
877  : MoFEM::EdgeElementForcesAndSourcesCore(m_field), B(PETSC_NULL),
878  F(PETSC_NULL) {}
879  int getRule(int order) { return 2 * order; };
880 
883 
885 
886  if (B != PETSC_NULL) {
887  snes_B = B;
888  }
889 
890  if (F != PETSC_NULL) {
891  snes_f = F;
892  }
893 
895  }
896  };
897 
898  boost::shared_ptr<MyEdgeFE> feRhsPtr, feLhsPtr;
899 
901  MyEdgeFE &getLoopFeRhs() { return feRhs; }
903  MyEdgeFE &getLoopFeLhs() { return feLhs; }
904 
905  double aLpha;
906 
909  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
910  "Get edge sliding constrains element scaling",
911  "none");
912  CHKERRQ(ierr);
913  CHKERR PetscOptionsScalar("-edge_sliding_alpha", "scaling parameter", "",
914  aLpha, &aLpha, PETSC_NULL);
915  ierr = PetscOptionsEnd();
916  CHKERRQ(ierr);
918  }
919 
921  : mField(m_field), feRhsPtr(new MyEdgeFE(m_field)),
922  feLhsPtr(new MyEdgeFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr),
923  aLpha(1) {
924  ierr = getOptions();
925  CHKERRABORT(PETSC_COMM_WORLD, ierr);
926  }
927 
928  struct OpJacobian
930 
931  const int tAg;
932  boost::shared_ptr<VectorDouble> activeVariablesPtr;
933  boost::shared_ptr<VectorDouble> resultsPtr;
934  boost::shared_ptr<MatrixDouble> jacobianPtr;
936 
937  const double &aLpha;
938 
939  OpJacobian(int tag, const std::string field_name,
940  boost::shared_ptr<VectorDouble> &active_variables_ptr,
941  boost::shared_ptr<VectorDouble> &results_ptr,
942  boost::shared_ptr<MatrixDouble> &jacobian_ptr,
943  bool evaluate_jacobian, const double &alpha)
945  field_name, UserDataOperator::OPCOL),
946  tAg(tag), activeVariablesPtr(active_variables_ptr),
947  resultsPtr(results_ptr), jacobianPtr(jacobian_ptr),
948  evaluateJacobian(evaluate_jacobian), aLpha(alpha) {}
949 
950  MoFEMErrorCode doWork(int side, EntityType type,
953  if (type != MBVERTEX)
955 
960 
961  Tag th0, th1, th2, th3;
963  th1, th2, th3);
964  FTensor::Tensor1<double, 3> t_edge_base0, t_edge_base1;
965  EntityHandle fe_ent = getFEEntityHandle();
966  CHKERR getEdgeFE()->mField.get_moab().tag_get_data(th0, &fe_ent, 1,
967  &t_edge_base0(0));
968  CHKERR getEdgeFE()->mField.get_moab().tag_get_data(th1, &fe_ent, 1,
969  &t_edge_base1(0));
970 
971  VectorInt &indices = data.getIndices();
972 
973  trace_on(tAg);
974 
975  ublas::vector<adouble> lambda_dofs(4);
976  for (int dd = 0; dd != 4; ++dd) {
977  lambda_dofs[dd] <<= (*activeVariablesPtr)[dd];
978  }
979  ublas::vector<adouble> position_dofs(6);
980  for (int dd = 0; dd != 6; ++dd) {
981  position_dofs[dd] <<= (*activeVariablesPtr)[4 + dd];
982  }
983 
985  &position_dofs[0], &position_dofs[1], &position_dofs[2]);
987  &position_dofs[3], &position_dofs[4], &position_dofs[5]);
988 
990  t_tangent(i) = t_node1(i) - t_node0(i);
991  adouble l = sqrt(t_tangent(i) * t_tangent(i));
992  t_tangent(i) /= l;
993 
994  adouble t_dot0, t_dot1;
995  t_dot0 = t_edge_base0(i) * t_tangent(i);
996  t_dot1 = t_edge_base1(i) * t_tangent(i);
997 
998  FTensor::Tensor1<adouble, 3> t_base0, t_base1;
999  t_base0(i) = t_edge_base0(i) - t_dot0 * t_tangent(i);
1000  t_base1(i) = t_edge_base1(i) - t_dot1 * t_tangent(i);
1001  t_base0(i) /= sqrt(t_base0(i) * t_base0(i));
1002  t_base1(i) /= sqrt(t_base1(i) * t_base1(i));
1003 
1004  auto t_base_fun1 = data.getFTensor0N();
1005  auto t_base_fun2 = data.getFTensor0N();
1006  FTensor::Tensor1<adouble, 3> t_position;
1009  auto t_coord_ref = getFTensor1CoordsAtGaussPts();
1010 
1011  ublas::vector<adouble> c_vec(4);
1012  ublas::vector<adouble> f_vec(6);
1013  c_vec.clear();
1014  f_vec.clear();
1015 
1016  int nb_gauss_pts = data.getN().size1();
1017  int nb_base_functions = data.getN().size2();
1018  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1019 
1020  FTensor::Tensor1<adouble *, 3> t_position_dofs(
1021  &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
1022  FTensor::Tensor1<adouble *, 2> t_lambda_dof(&lambda_dofs[0],
1023  &lambda_dofs[1], 2);
1024 
1025  t_position(i) = 0;
1026  t_lambda(j) = 0;
1027  for (int bb = 0; bb != nb_base_functions; ++bb) {
1028  t_position(i) += t_base_fun1 * t_position_dofs(i);
1029  t_lambda(j) += t_base_fun1 * t_lambda_dof(j);
1030  ++t_base_fun1;
1031  ++t_position_dofs;
1032  ++t_lambda_dof;
1033  }
1034 
1035  t_delta(i) = t_position(i) - t_coord_ref(i);
1036  adouble dot0 = t_base0(i) * t_delta(i);
1037  adouble dot1 = t_base1(i) * t_delta(i);
1038 
1039  adouble w = getGaussPts()(1, gg) * l * aLpha;
1040  adouble val, val1, val2;
1041  FTensor::Tensor1<adouble *, 2> t_c(&c_vec[0], &c_vec[1], 2);
1042  FTensor::Tensor1<adouble *, 3> t_f(&f_vec[0], &f_vec[1], &f_vec[2], 3);
1043  for (int bb = 0; bb != nb_base_functions; ++bb) {
1044  if (indices[2 * bb] != -1) {
1045  val = w * t_base_fun2;
1046  t_c(N0) += val * dot0;
1047  t_c(N1) += val * dot1;
1048  val1 = val * t_lambda(N0);
1049  val2 = val * t_lambda(N1);
1050  t_f(i) += val1 * t_base0(i) + val2 * t_base1(i);
1051  }
1052  ++t_c;
1053  ++t_f;
1054  ++t_base_fun2;
1055  }
1056 
1057  ++t_coord_ref;
1058  }
1059 
1060  for (int rr = 0; rr != 4; ++rr) {
1061  c_vec[rr] >>= (*resultsPtr)[rr];
1062  }
1063  for (int rr = 0; rr != 6; ++rr) {
1064  f_vec(rr) >>= (*resultsPtr)[4 + rr];
1065  }
1066 
1067  trace_off();
1068 
1069  if (evaluateJacobian) {
1070  double *jac_ptr[4 + 6];
1071  for (int rr = 0; rr != 4 + 6; ++rr) {
1072  jac_ptr[rr] = &(*jacobianPtr)(rr, 0);
1073  }
1074  // play recorder for jacobians
1075  int r =
1076  ::jacobian(tAg, 4 + 6, 4 + 6, &(*activeVariablesPtr)[0], jac_ptr);
1077  if (r < 0) {
1078  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1079  "ADOL-C function evaluation with error");
1080  }
1081  }
1082 
1084  }
1085  };
1086 
1087  MoFEMErrorCode setOperators(int tag, Range edges, Range faces,
1088  const std::string lagrange_multipliers_field_name,
1089  const std::string material_field_name) {
1092  edges, faces);
1093  CHKERR setOperators(tag, lagrange_multipliers_field_name,
1094  material_field_name);
1096  }
1097 
1099  const std::string lagrange_multipliers_field_name,
1100  const std::string material_field_name,
1101  const double *alpha = nullptr) {
1103 
1104  if (alpha) {
1105  aLpha = *alpha;
1106  }
1107 
1108  boost::shared_ptr<VectorDouble> active_variables_ptr(
1109  new VectorDouble(4 + 6));
1110  boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(4 + 6));
1111  boost::shared_ptr<MatrixDouble> jacobian_ptr(
1112  new MatrixDouble(4 + 6, 4 + 6));
1113 
1114  feRhs.getOpPtrVector().clear();
1115  feRhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
1116  lagrange_multipliers_field_name, active_variables_ptr));
1118  material_field_name, active_variables_ptr));
1119  feRhs.getOpPtrVector().push_back(new OpJacobian(
1120  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1121  jacobian_ptr, false, aLpha));
1122  feRhs.getOpPtrVector().push_back(
1123  new OpAssembleRhs<4, 6>(lagrange_multipliers_field_name, results_ptr));
1124  feRhs.getOpPtrVector().push_back(
1125  new OpAssembleRhs<4, 6>(material_field_name, results_ptr));
1126 
1127  // Adding operators to calculate the left hand side
1128  feLhs.getOpPtrVector().clear();
1129  feLhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
1130  lagrange_multipliers_field_name, active_variables_ptr));
1132  material_field_name, active_variables_ptr));
1133  feLhs.getOpPtrVector().push_back(new OpJacobian(
1134  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1135  jacobian_ptr, true, aLpha));
1136  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1137  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1138  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1139  material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
1140  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1141  material_field_name, material_field_name, jacobian_ptr));
1142 
1144  }
1145 
1147  setOperatorsConstrainOnly(int tag, Range edges, Range faces,
1148  const std::string lagrange_multipliers_field_name,
1149  const std::string material_field_name) {
1151 
1153  edges, faces);
1154  CHKERR setOperatorsConstrainOnly(tag, lagrange_multipliers_field_name,
1155  material_field_name);
1157  }
1158 
1161  const std::string lagrange_multipliers_field_name,
1162  const std::string material_field_name) {
1164 
1165  boost::shared_ptr<VectorDouble> active_variables_ptr(
1166  new VectorDouble(4 + 6));
1167  boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(4 + 6));
1168  boost::shared_ptr<MatrixDouble> jacobian_ptr(
1169  new MatrixDouble(4 + 6, 4 + 6));
1170 
1171  // Adding operators to calculate the left hand side
1172  feLhs.getOpPtrVector().clear();
1173  feLhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
1174  lagrange_multipliers_field_name, active_variables_ptr));
1176  material_field_name, active_variables_ptr));
1177  feLhs.getOpPtrVector().push_back(new OpJacobian(
1178  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1179  jacobian_ptr, true, aLpha));
1180  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1181  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1182 
1184  }
1185 };
1186 
1187 #endif // WITH_ADOL_C
1188 
1189 #endif // __SURFACE_SLIDING_CONSTRAINS_HPP__
MoFEMErrorCode setOperators(int tag, Range edges, Range faces, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
OpAssembleRhs(const std::string field_name, boost::shared_ptr< VectorDouble > &results_ptr)
Deprecated interface functions.
const std::string & getFEName() const
Get name of the element.
virtual moab::Interface & get_moab()=0
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
SurfaceSlidingConstrains(MoFEM::Interface &m_field, DriverElementOrientation &orientation)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
static Index< 'p', 3 > p
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
FTensor::Tensor1< double, 2 > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, 2 > &t_disp)
Definition: ContactOps.hpp:137
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
boost::shared_ptr< VectorDouble > activeVariablesPtr
virtual ~SurfaceSlidingConstrains()=default
Class implemented by user to detect face orientation.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
static Index< 'l', 3 > l
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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)
Implementation of surface sliding constrains.
static Index< 'n', 3 > n
boost::shared_ptr< MyEdgeFE > feRhsPtr
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:95
Vec & snes_f
residual
static MoFEMErrorCode createTag(moab::Interface &moab, Tag &th0, Tag &th1, Tag &th2, Tag &th3)
MoFEMErrorCode setOperatorsConstrainOnly(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
DriverElementOrientation & crackFrontOrientation
bool sYmm
If true assume that matrix is symmetric structure.
Mat & snes_B
preconditioner of jacobian matrix
double tol
GenericSliding()=default
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
static Index< 'i', 3 > i
static MoFEMErrorCode saveEdges(moab::Interface &moab, std::string name, Range edges, Range *faces=nullptr)
boost::shared_ptr< VectorDouble > activeVariablesPtr
boost::shared_ptr< MyEdgeFE > feLhsPtr
CHKERRQ(ierr)
static Number< 2 > N2
DataForcesAndSourcesCore::EntData EntData
MoFEMErrorCode setOperatorsConstrainOnly(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
boost::shared_ptr< VectorDouble > resultsPtr
MoFEMErrorCode setOperators(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name, const double *alpha=nullptr)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEMErrorCode setOperators(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name, const double *alpha=nullptr)
static Index< 'j', 3 > j
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
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
boost::shared_ptr< VectorDouble > resultsPtr
boost::shared_ptr< VectorDouble > activeVariablesPtr
EdgeSlidingConstrains(MoFEM::Interface &m_field)
OpGetActiveDofsPositions(const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr)
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
#define CHKERR
Inline error check.
Definition: definitions.h:604
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
virtual MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field, const FEMethod *fe_method_ptr)
const double r
rate factor
static Number< 1 > N1
boost::shared_ptr< MatrixDouble > jacobianPtr
MoFEMErrorCode preProcess()
function is run at the beginning of loop
constexpr int order
virtual ~GenericSliding()=default
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
boost::shared_ptr< VectorDouble > activeVariablesPtr
boost::shared_ptr< MyTriangleFE > feLhsPtr
Shape preserving constrains, i.e. nodes sliding on body surface.
boost::shared_ptr< MyTriangleFE > feRhsPtr
boost::shared_ptr< MatrixDouble > jacobianPtr
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)
static Index< 'k', 3 > k
MoFEMErrorCode setOperatorsConstrainOnly(int tag, Range edges, Range faces, const std::string lagrange_multipliers_field_name, const std::string material_field_name)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
OpAssembleLhs(const std::string field_name_row, const std::string field_name_col, boost::shared_ptr< MatrixDouble > &jacobian_ptr)
boost::shared_ptr< VectorDouble > resultsPtr
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
static Number< 0 > N0
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
static MoFEMErrorCode setTags(moab::Interface &moab, Range edges, Range tris, bool number_pathes=true)
static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges, Range tris)
boost::shared_ptr< MatrixDouble > jacobianPtr
double w(const double g, const double t)
Definition: ContactOps.hpp:169
OpGetActiveDofsLambda(const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr)