v0.9.1
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 
32  boost::shared_ptr<VectorDouble> activeVariablesPtr;
33  OpGetActiveDofsLambda(const std::string field_name,
34  boost::shared_ptr<VectorDouble> &active_variables_ptr)
36  field_name, UserDataOperator::OPCOL),
37  activeVariablesPtr(active_variables_ptr) {}
38  MoFEMErrorCode doWork(int side, EntityType type,
41  if (type == MBVERTEX) {
42  for (unsigned int dd = 0; dd != data.getFieldData().size(); ++dd)
43  (*activeVariablesPtr)[dd] = data.getFieldData()[dd];
44  }
46  }
47  };
48 
49  template <int SizeLambda>
52  boost::shared_ptr<VectorDouble> activeVariablesPtr;
54  const std::string field_name,
55  boost::shared_ptr<VectorDouble> &active_variables_ptr)
57  field_name, UserDataOperator::OPCOL),
58  activeVariablesPtr(active_variables_ptr) {}
59  MoFEMErrorCode doWork(int side, EntityType type,
62  if (type == MBVERTEX) {
63  for (unsigned int dd = 0; dd != data.getFieldData().size(); ++dd)
64  (*activeVariablesPtr)[SizeLambda + dd] = data.getFieldData()[dd];
65  }
67  }
68  };
69 
70  template <int SizeLambda, int SizePositions>
72 
73  boost::shared_ptr<VectorDouble> resultsPtr;
74 
75  OpAssembleRhs(const std::string field_name,
76  boost::shared_ptr<VectorDouble> &results_ptr)
78  field_name, UserDataOperator::OPROW),
79  resultsPtr(results_ptr) {}
80 
81  MoFEMErrorCode doWork(int side, EntityType type,
84  if (type != MBVERTEX)
86  VectorInt &indices = data.getIndices();
87  int shift = 0;
88  if (indices.empty())
90  else if (indices.size() == SizeLambda)
91  shift = 0;
92  else if (indices.size() == SizePositions)
93  shift = SizeLambda;
94  else
95  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
96  "Element %s: Data inconsistency nb of indices %d",
97  getFEName().c_str(), indices.size());
98 
99  CHKERR VecSetOption(getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
100  CHKERR VecSetValues(getSNESf(), indices.size(), &indices[0],
101  &(*resultsPtr)[shift], ADD_VALUES);
103  }
104  };
105  template <int SizeLambda, int SizePositions>
107 
108  boost::shared_ptr<MatrixDouble> jacobianPtr;
109 
110  OpAssembleLhs(const std::string field_name_row,
111  const std::string field_name_col,
112  boost::shared_ptr<MatrixDouble> &jacobian_ptr)
114  field_name_row, field_name_col, UserDataOperator::OPROWCOL),
115  jacobianPtr(jacobian_ptr) {
116  sYmm = false;
117  }
118 
119  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
120  EntityType col_type,
124  if (row_type != MBVERTEX)
126  if (col_type != MBVERTEX)
128  VectorInt &row_indices = row_data.getIndices();
129  VectorInt &col_indices = col_data.getIndices();
130  if (row_indices.empty() || col_indices.empty())
132 
133  int shift_row = 0;
134  if (row_indices.size() == SizeLambda)
135  shift_row = 0;
136  else if (row_indices.size() == SizePositions)
137  shift_row = SizeLambda;
138  else
139  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
140  "Data inconsistency");
141 
142  int shift_col = 0;
143  if (col_indices.size() == SizeLambda)
144  shift_col = 0;
145  else if (col_indices.size() == SizePositions)
146  shift_col = SizeLambda;
147  else
148  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
149  "Data inconsistency");
150 
151  MatrixDouble jac(row_indices.size(), col_indices.size());
152  for (unsigned int rr = 0; rr != row_indices.size(); ++rr) {
153  for (unsigned int cc = 0; cc != col_indices.size(); ++cc) {
154  jac(rr, cc) = (*jacobianPtr)(shift_row + rr, shift_col + cc);
155  if (jac(rr, cc) != jac(rr, cc))
156  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
157  "Jacobian assemble not a number jac(rr,cc) = %3.4f",
158  jac(rr, cc));
159  }
160  }
161  CHKERR MatSetValues(getSNESB(), row_indices.size(), &row_indices[0],
162  col_indices.size(), &col_indices[0], &jac(0, 0),
163  ADD_VALUES);
165  }
166  };
167 };
168 
169 /** \brief Shape preserving constrains, i.e. nodes sliding on body surface.
170 
171  Derivation and implementation of constrains preserving body surface,
172  i.e. body shape and volume.
173 
174  The idea starts form observation that body shape can be globally characterized
175  by constant calculated as volume over its area
176  \f[
177  \frac{V}{A} = C
178  \f]
179  Above equation expressed in integral form is
180  \f[
181  \int_\Omega \textrm{d}V = C \int_\Gamma \textrm{d}S
182  \f]
183  where notting that,
184  \f[
185  \frac{1}{3}
186  \int_\Omega \textrm{div}[\mathbf{X}] \textrm{d}\Omega
187  =
188  C \int_\Gamma \textrm{d}S
189  \f]
190  and applying Gauss theorem we get
191  \f[
192  \int_\Gamma
193  \mathbf{X}\cdot \frac{\mathbf{N}}{\|\mathbf{N}\|}
194  \textrm{d}\Gamma
195  =
196  3C \int_\Gamma \textrm{d}S.
197  \f]
198  Drooping integrals on both sides, and linearizing equation, we get
199  \f[
200  \frac{\mathbf{N}}{\|\mathbf{N}\|} \cdot \delta \mathbf{X}
201  =
202  3C - \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot \mathbf{X}
203  \f]
204  where \f$\delta \mathbf{X}\f$ is displacement sub-inctrement. Above equation
205  is a constrain if satisfied in body shape and volume is conserved. Final form
206  of constrain equation is \f[ \mathcal{r} =
207  \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot \mathbf{X}
208  -
209  \frac{\mathbf{N_0}}{\|\mathbf{N_0}\|}\cdot \mathbf{X}_0 =
210  \frac{\mathbf{N}}{\|\mathbf{N}\|}\cdot (\mathbf{X}-\mathbf{X}_0)
211  \f]
212 
213  In the spirit of finite element method the constrain equation is multiplied
214  by shape functions and enforce using Lagrange multiplier method
215  \f[
216  \int_\Gamma \mathbf{N}^\mathsf{T}_\lambda
217  \left(
218  \frac{\mathbf{N}}{\|\mathbf{N}\|}\mathbf{N}_\mathbf{X}\cdot
219  (\overline{\mathbf{X}}-\overline{\mathbf{X}}_0)
220  \right)
221  \|\mathbf{N}\|
222  \textrm{d}\Gamma
223  =
224  \mathbf{0}.
225  \f]
226  Above equation is nonlinear, applying to it Taylor expansion, we can get form
227  which can be used with Newton interactive method \f[ \begin{split}
228  &\int_\Gamma \mathbf{N}^\mathsf{T}_\lambda
229  \left\{
230  \mathbf{N}\mathbf{N}_\mathbf{X}
231  +
232  \left(\mathbf{X}-\mathbf{X}_0\right) \cdot
233  \left(
234  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
235  -
236  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
237  \right)
238  \right\}
239  \textrm{d}\Gamma
240  \cdot
241  \delta
242  \overline{\mathbf{X}}\\
243  =
244  &\int_\Gamma \mathbf{N}^\mathsf{T}_\lambda
245  \mathbf{N}\cdot(\mathbf{X}-\mathbf{X}_0)
246  \textrm{d}\Gamma
247  \end{split}.
248  \f]
249  Equation expressing forces on shape as result of constrains, as result
250  Lagrange multiplier method have following form \f[ \begin{split}
251  &\int_\Gamma
252  \mathbf{N}^\mathsf{T}_\mathbf{X} \cdot \mathbf{N}
253  \mathbf{N}_\lambda
254  \textrm{d}\Gamma
255  \cdot
256  \delta\overline{\lambda}\\
257  +
258  &\int_\Gamma
259  \lambda
260  \mathbf{N}^\mathsf{T}_\mathbf{X}
261  \left(
262  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
263  -
264  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
265  \right)
266  \textrm{d}\Gamma
267  \delta\overline{\mathbf{X}}\\
268  =
269  &\int_\Gamma
270  \lambda
271  \mathbf{N}^\mathsf{T}_\mathbf{X} \cdot \mathbf{N}
272  \textrm{d}\Gamma
273  \end{split}
274  \f]
275 
276  Above equations are assembled into global system of equations as following
277  \f[
278  \left[
279  \begin{array}{cc}
280  \mathbf{K} + \mathbf{B} & \mathbf{C}^\mathsf{T} \\
281  \mathbf{C} + \mathbf{A} & 0
282  \end{array}
283  \right]
284  \left\{
285  \begin{array}{c}
286  \delta \overline{\mathbf{X}} \\
287  \delta \overline{\lambda}
288  \end{array}
289  \right\}=
290  \left[
291  \begin{array}{c}
292  \mathbf{f} - \mathbf{C}^\mathsf{T}\overline{\lambda} \\
293  \overline{\mathbf{r}}
294  \end{array}
295  \right]
296  \f]
297  where
298  \f[
299  \mathbf{C}=
300  \int_\Gamma
301  \mathbf{N}_\lambda^\mathsf{T}
302  \mathbf{N} \cdot
303  \mathbf{N}_\mathbf{X}
304  \textrm{d}\Gamma,
305  \f]
306  \f[
307  \mathbf{B}=
308  \int_\Gamma
309  \lambda
310  (\mathbf{X}-\mathbf{X}_0)\cdot
311  \left(
312  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
313  -
314  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
315  \right)
316  \textrm{d}\Gamma
317  \f]
318  and
319  \f[
320  \mathbf{A}=
321  \int_\Gamma
322  \mathbf{N}^\mathsf{T}_\lambda
323  \left(\mathbf{X}-\mathbf{X}_0\right) \cdot
324  \left(
325  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\xi}\right]\cdot\mathbf{B}_\eta
326  -
327  \textrm{Spin}\left[\frac{\partial\mathbf{X}}{\partial\eta}\right]\cdot\mathbf{B}_\xi
328  \right).
329  \f]
330 
331 */
333 
334  /** \brief Class implemented by user to detect face orientation
335 
336  If mesh generated is with surface mesher, usually you don't have to do
337  nothing, all elements on the surface have consistent orientation. In case of
338  internal faces or if you do something with mesh connectivity which breaks
339  orientation on the face, you have to implement method which will set
340  orientation to face.
341 
342  */
344 
346 
347  virtual MoFEMErrorCode
349  const FEMethod *fe_method_ptr) {
351  elementOrientation = 1;
353  }
354  };
355 
357 
359 
360  Mat B;
361  Vec F;
362 
364  : MoFEM::FaceElementForcesAndSourcesCore(m_field), B(PETSC_NULL),
365  F(PETSC_NULL) {}
366  int getRule(int order) { return 2 * order; };
367 
370 
371  CHKERR MoFEM::FaceElementForcesAndSourcesCore::preProcess();
372 
373  if (B != PETSC_NULL) {
374  snes_B = B;
375  }
376 
377  if (F != PETSC_NULL) {
378  snes_f = F;
379  }
380 
382  }
383  };
384 
385  boost::shared_ptr<MyTriangleFE> feRhsPtr, feLhsPtr;
386 
391 
393 
394  double aLpha;
397  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
398  "Get surface sliding constrains element scaling",
399  "none");
400  CHKERRQ(ierr);
401  CHKERR PetscOptionsScalar("-surface_sliding_alpha", "scaling parameter", "",
402  aLpha, &aLpha, PETSC_NULL);
403  ierr = PetscOptionsEnd();
404  CHKERRQ(ierr);
406  }
407 
409  DriverElementOrientation &orientation)
410  : mField(m_field), feRhsPtr(new MyTriangleFE(m_field)),
411  feLhsPtr(new MyTriangleFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr),
412  crackFrontOrientation(orientation), aLpha(1) {
413 
414  ierr = getOptions();
415  CHKERRABORT(PETSC_COMM_WORLD, ierr);
416  }
417 
418  struct OpJacobian
420 
421  const int tAg;
422  boost::shared_ptr<VectorDouble> activeVariablesPtr;
423  boost::shared_ptr<VectorDouble> resultsPtr;
424  boost::shared_ptr<MatrixDouble> jacobianPtr;
427 
428  const double &aLpha;
429 
430  OpJacobian(int tag, const std::string field_name,
431  boost::shared_ptr<VectorDouble> &active_variables_ptr,
432  boost::shared_ptr<VectorDouble> &results_ptr,
433  boost::shared_ptr<MatrixDouble> &jacobian_ptr,
434  DriverElementOrientation &orientation, bool evaluate_jacobian,
435  const double &alpha)
437  field_name, UserDataOperator::OPCOL),
438  tAg(tag), activeVariablesPtr(active_variables_ptr),
439  resultsPtr(results_ptr), jacobianPtr(jacobian_ptr),
440  oRientation(orientation), evaluateJacobian(evaluate_jacobian),
441  aLpha(alpha) {}
442 
443  MoFEMErrorCode doWork(int side, EntityType type,
446  if (type != MBVERTEX)
448 
449  VectorInt &indices = data.getIndices();
450 
451  trace_on(tAg);
452 
453  ublas::vector<adouble> lambda_dofs(3);
454  for (int dd = 0; dd != 3; ++dd) {
455  lambda_dofs[dd] <<= (*activeVariablesPtr)[dd];
456  }
457  ublas::vector<adouble> position_dofs(9);
458  for (int dd = 0; dd != 9; ++dd) {
459  position_dofs[dd] <<= (*activeVariablesPtr)[3 + dd];
460  }
461 
468 
470  getFEMethod());
472 
473  int nb_gauss_pts = data.getN().size1();
474  int nb_base_functions = data.getN().size2();
475  auto t_base1 = data.getFTensor0N();
476  auto t_base2 = data.getFTensor0N();
477  auto t_diff_base = data.getFTensor1DiffN<2>();
478  FTensor::Tensor1<adouble, 3> t_position;
479  FTensor::Tensor1<adouble, 3> t_position_ksi;
480  FTensor::Tensor1<adouble, 3> t_position_eta;
481  adouble lambda;
484 
485  ublas::vector<adouble> c_vec(3);
486  ublas::vector<adouble> f_vec(9);
487  c_vec.clear();
488  f_vec.clear();
489 
490  auto t_coord_ref = getFTensor1CoordsAtGaussPts();
491 
492  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
493 
494  FTensor::Tensor1<adouble *, 3> t_position_dofs(
495  &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
496  FTensor::Tensor0<adouble *> t_lambda_dof(&lambda_dofs[0]);
497 
498  t_position(i) = 0;
499  t_position_ksi(i) = 0;
500  t_position_eta(i) = 0;
501  lambda = 0;
502 
503  for (int bb = 0; bb != nb_base_functions; ++bb) {
504  t_position(i) += t_base1 * t_position_dofs(i);
505  t_position_ksi(i) += t_diff_base(N0) * t_position_dofs(i);
506  t_position_eta(i) += t_diff_base(N1) * t_position_dofs(i);
507  lambda += t_base1 * t_lambda_dof;
508  ++t_base1;
509  ++t_position_dofs;
510  ++t_lambda_dof;
511  ++t_diff_base;
512  }
513 
514  t_delta(i) = t_position(i) - t_coord_ref(i);
515  t_normal(k) = FTensor::cross(t_position_ksi(i), t_position_eta(j), k);
516 
517  double w = getGaussPts()(2, gg) * 0.5 * aLpha;
518  adouble val;
519  FTensor::Tensor0<adouble *> t_c(&c_vec[0]);
520  FTensor::Tensor1<adouble *, 3> t_f(&f_vec[0], &f_vec[1], &f_vec[2], 3);
521 
522  adouble area = sqrt(t_normal(i) * t_normal(i));
523 
524  for (int bb = 0; bb != nb_base_functions; ++bb) {
525  if (indices[bb] != -1) {
526  val = w * eo * t_base2;
527  t_c += val * t_normal(i) * t_delta(i);
528  val *= lambda;
529  t_f(i) += val * t_normal(i);
530  }
531  ++t_c;
532  ++t_f;
533  ++t_base2;
534  }
535 
536  ++t_coord_ref;
537  }
538 
539  for (int rr = 0; rr != 3; ++rr)
540  c_vec[rr] >>= (*resultsPtr)[rr];
541 
542  for (int rr = 0; rr != 9; ++rr)
543  f_vec(rr) >>= (*resultsPtr)[3 + rr];
544 
545 
546  trace_off();
547 
548  if (evaluateJacobian) {
549  double *jac_ptr[3 + 9];
550  for (int rr = 0; rr != 3 + 9; ++rr)
551  jac_ptr[rr] = &(*jacobianPtr)(rr, 0);
552 
553  // play recorder for jacobians
554  int r =
555  ::jacobian(tAg, 3 + 9, 3 + 9, &(*activeVariablesPtr)[0], jac_ptr);
556  if (r < 0)
557  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
558  "ADOL-C function evaluation with error");
559 
560  }
561 
563  }
564  };
565 
567  const std::string lagrange_multipliers_field_name,
568  const std::string material_field_name,
569  const double *alpha = nullptr) {
571 
572  if (alpha != nullptr) {
573  aLpha = *alpha;
574  }
575 
576  boost::shared_ptr<VectorDouble> active_variables_ptr(
577  new VectorDouble(3 + 9));
578  boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
579  boost::shared_ptr<MatrixDouble> jacobian_ptr(
580  new MatrixDouble(3 + 9, 3 + 9));
581 
582  feRhs.getOpPtrVector().clear();
584  lagrange_multipliers_field_name, active_variables_ptr));
586  material_field_name, active_variables_ptr));
587  feRhs.getOpPtrVector().push_back(new OpJacobian(
588  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
589  jacobian_ptr, crackFrontOrientation, false, aLpha));
590  feRhs.getOpPtrVector().push_back(
591  new OpAssembleRhs<3, 9>(lagrange_multipliers_field_name, results_ptr));
592  feRhs.getOpPtrVector().push_back(
593  new OpAssembleRhs<3, 9>(material_field_name, results_ptr));
594 
595  // Adding operators to calculate the left hand side
596  feLhs.getOpPtrVector().clear();
598  lagrange_multipliers_field_name, active_variables_ptr));
600  material_field_name, active_variables_ptr));
601  feLhs.getOpPtrVector().push_back(new OpJacobian(
602  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
603  jacobian_ptr, crackFrontOrientation, true, aLpha));
604  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
605  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
606  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
607  material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
608  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
609  material_field_name, material_field_name, jacobian_ptr));
610 
612  }
613 
616  const std::string lagrange_multipliers_field_name,
617  const std::string material_field_name) {
619 
620  boost::shared_ptr<VectorDouble> active_variables_ptr(
621  new VectorDouble(3 + 9));
622  boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(3 + 9));
623  boost::shared_ptr<MatrixDouble> jacobian_ptr(
624  new MatrixDouble(3 + 9, 3 + 9));
625 
626  // Adding operators to calculate the left hand side
627  feLhs.getOpPtrVector().clear();
629  lagrange_multipliers_field_name, active_variables_ptr));
631  material_field_name, active_variables_ptr));
632  feLhs.getOpPtrVector().push_back(new OpJacobian(
633  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
634  jacobian_ptr, crackFrontOrientation, true, aLpha));
635  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<3, 9>(
636  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
637 
639  }
640 };
641 
643 
645 
646  static MoFEMErrorCode createTag(moab::Interface &moab, Tag &th0, Tag &th1,
647  Tag &th2, Tag &th3) {
649  double def_val[] = {0, 0, 0};
650  CHKERR moab.tag_get_handle("EDGE_BASE0", 3, MB_TYPE_DOUBLE, th0,
651  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
652  CHKERR moab.tag_get_handle("EDGE_BASE1", 3, MB_TYPE_DOUBLE, th1,
653  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
654  int def_numb_val[] = {-1};
655  CHKERR moab.tag_get_handle("PATCH_NUMBER", 1, MB_TYPE_INTEGER, th2,
656  MB_TAG_CREAT | MB_TAG_SPARSE, def_numb_val);
657  int def_orientation_val[] = {1};
658  CHKERR moab.tag_get_handle("PATCH_ORIENTATION", 1, MB_TYPE_INTEGER, th3,
659  MB_TAG_CREAT | MB_TAG_SPARSE,
660  def_orientation_val);
661 
663  }
664 
665  static MoFEMErrorCode numberSurfaces(moab::Interface &moab, Range edges,
666  Range tris) {
668 
669  auto get_edges = [&](const Range &ents) {
670  Range edges;
671  CHKERR moab.get_adjacencies(ents, 1, false, edges,
672  moab::Interface::UNION);
673  return edges;
674  };
675 
676  auto get_face_adj = [&, get_edges](const Range &faces) {
677  Range adj_faces;
678  CHKERR moab.get_adjacencies(subtract(get_edges(faces), edges), 2, false,
679  adj_faces, moab::Interface::UNION);
680  return intersect(adj_faces, tris);
681  };
682 
683  auto get_patch = [&, get_face_adj](const EntityHandle face) {
684  Range patch_ents;
685  patch_ents.insert(face);
686  unsigned int nb0;
687  do {
688  nb0 = patch_ents.size();
689  patch_ents.merge(get_face_adj(patch_ents));
690  } while (nb0 != patch_ents.size());
691  return patch_ents;
692  };
693 
694  auto get_patches = [&]() {
695  std::vector<Range> patches;
696  while (!tris.empty()) {
697  patches.push_back(get_patch(tris[0]));
698  tris = subtract(tris, patches.back());
699  }
700  return patches;
701  };
702 
703  Tag th0, th1, th2, th3;
704  CHKERR createTag(moab, th0, th1, th2, th3);
705 
706  auto patches = get_patches();
707  int pp = 0;
708  for (auto patch : patches) {
709  // cerr << "pp: " << pp << endl;
710  // cerr << patch << endl;
711  std::vector<int> tags_vals(patch.size(), pp);
712  CHKERR moab.tag_set_data(th2, patch, &*tags_vals.begin());
713  ++pp;
714  }
715 
717  }
718 
719  static MoFEMErrorCode setTags(moab::Interface &moab, Range edges,
720  Range tris, bool number_pathes = true) {
722  Tag th0, th1, th2, th3;
723  CHKERR createTag(moab, th0, th1, th2, th3);
724  if (number_pathes) {
725  CHKERR numberSurfaces(moab, edges, tris);
726  }
727  for (auto edge : edges) {
728  Range adj_faces;
729  CHKERR moab.get_adjacencies(&edge, 1, 2, false, adj_faces);
730  adj_faces = intersect(adj_faces, tris);
731  if (adj_faces.size() != 2) {
732  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
733  "Should be 2 faces adjacent to edge but is %d",
734  adj_faces.size());
735  }
736  VectorDouble3 v[2] = {VectorDouble3(3), VectorDouble3(3)};
737  auto get_tensor_from_vec = [](VectorDouble3 &v) {
739  &v[0], &v[1], &v[2]);
740  };
741 
745 
746  std::array<double, 3> areas;
747  auto calculate_normals = [&, get_tensor_from_vec]() {
748  int ff = 0;
749  for (auto face : adj_faces) {
750  double &x = (v[ff])[0];
751  double &y = (v[ff])[1];
752  double &z = (v[ff])[2];
753  moab::Util::normal(&moab, face, x, y, z);
754  auto t_n = get_tensor_from_vec(v[ff]);
755  areas[ff] = sqrt(t_n(i) * t_n(i));
756  t_n(i) /= areas[ff];
757  int orientation;
758  CHKERR moab.tag_get_data(th3, &face, 1, &orientation);
759  if (orientation == -1) {
760  t_n(i) *= -1;
761  }
762  ++ff;
763  }
764  };
765  calculate_normals();
766 
767  auto get_patch_number = [&]() {
768  std::vector<int> p = {0, 0};
769  CHKERR moab.tag_get_data(th2, adj_faces, &*p.begin());
770  return p;
771  };
772 
773  std::array<EntityHandle, 2> faces_handles = {adj_faces[0],
774  adj_faces[1]};
775  auto order_normals = [&, get_patch_number]() {
776  auto p = get_patch_number();
777  if (p[0] < p[1]) {
778  v[0].swap(v[1]);
779  faces_handles[0] = adj_faces[1];
780  faces_handles[1] = adj_faces[0];
781  }
782  };
783  order_normals();
784 
786  auto t_n0 = get_tensor_from_vec(v[0]);
787  auto t_n1 = get_tensor_from_vec(v[1]);
788  t_cross(k) = FTensor::cross(t_n0(i), t_n1(j), k);
789 
790  auto get_base_for_coplanar_faces = [&]() {
792 
793  std::vector<EntityHandle> face_conn;
794  CHKERR moab.get_connectivity(&faces_handles[1], 1, face_conn, true);
795  std::vector<EntityHandle> edge_conn;
796  CHKERR moab.get_connectivity(&edge, 1, edge_conn, true);
797  std::sort(face_conn.begin(), face_conn.end());
798  std::sort(edge_conn.begin(), edge_conn.end());
799  int n = 0;
800  for (; n != 2; ++n)
801  if (face_conn[n] != edge_conn[n])
802  break;
803  VectorDouble3 coords_face(3);
804  CHKERR moab.get_coords(&face_conn[n], 1,
805  &*coords_face.data().begin());
806  VectorDouble6 coords_edge(6);
807  CHKERR moab.get_coords(&edge_conn[0], 2,
808  &*coords_edge.data().begin());
809  auto t_edge_n0 = FTensor::Tensor1<double, 3>(
810  coords_edge[0], coords_edge[1], coords_edge[2]);
811  auto t_edge_n1 = FTensor::Tensor1<double, 3>(
812  coords_edge[3], coords_edge[4], coords_edge[5]);
813  auto t_face_n = get_tensor_from_vec(coords_face);
814  t_face_n(i) -= t_edge_n0(i);
816  t_delta(i) = t_edge_n1(i) - t_edge_n0(i);
817  t_delta(i) /= sqrt(t_delta(i) * t_delta(i));
818  t_n1(k) = FTensor::cross(t_n0(i), t_delta(j), k);
819  if (t_n1(i) * t_face_n(i) < 0)
820  t_n1(i) *= -1;
821 
823  };
824 
825  auto get_base_for_not_planar_faces = [&]() {
827  t_n1(k) = FTensor::cross(t_n0(i), t_cross(j), k);
828  t_n1(i) /= sqrt(t_n1(i) * t_n1(i));
830  };
831 
832  // This a case when faces adjacent to edge are coplanar !!!
833  constexpr double tol = 1e-6;
834  if ((t_cross(i) * t_cross(i)) < tol)
835  CHKERR get_base_for_coplanar_faces();
836  else
837  CHKERR get_base_for_not_planar_faces();
838 
839  VectorDouble3 &v0 = v[0];
840  VectorDouble3 &v1 = v[1];
841  CHKERR moab.tag_set_data(th0, &edge, 1, &v0[0]);
842  CHKERR moab.tag_set_data(th1, &edge, 1, &v1[0]);
843  }
845  }
846 
847  static MoFEMErrorCode saveEdges(moab::Interface &moab, std::string name,
848  Range edges, Range *faces = nullptr) {
850  EntityHandle meshset;
851  Tag ths[4];
852  CHKERR createTag(moab, ths[0], ths[1], ths[2], ths[3]);
853  CHKERR moab.create_meshset(MESHSET_SET, meshset);
854  CHKERR moab.add_entities(meshset, edges);
855  if (faces != nullptr) {
856  CHKERR moab.add_entities(meshset, *faces);
857  }
858  CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1, ths, 4);
859  CHKERR moab.delete_entities(&meshset, 1);
861  }
862  };
863 
865 
867 
868  Mat B;
869  Vec F;
870 
872  : MoFEM::EdgeElementForcesAndSourcesCore(m_field), B(PETSC_NULL),
873  F(PETSC_NULL) {}
874  int getRule(int order) { return 2 * order; };
875 
878 
880 
881  if (B != PETSC_NULL) {
882  snes_B = B;
883  }
884 
885  if (F != PETSC_NULL) {
886  snes_f = F;
887  }
888 
890  }
891  };
892 
893  boost::shared_ptr<MyEdgeFE> feRhsPtr, feLhsPtr;
894 
896  MyEdgeFE &getLoopFeRhs() { return feRhs; }
898  MyEdgeFE &getLoopFeLhs() { return feLhs; }
899 
900  double aLpha;
901 
904  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "",
905  "Get edge sliding constrains element scaling",
906  "none");
907  CHKERRQ(ierr);
908  CHKERR PetscOptionsScalar("-edge_sliding_alpha", "scaling parameter", "",
909  aLpha, &aLpha, PETSC_NULL);
910  ierr = PetscOptionsEnd();
911  CHKERRQ(ierr);
913  }
914 
916  : mField(m_field), feRhsPtr(new MyEdgeFE(m_field)),
917  feLhsPtr(new MyEdgeFE(m_field)), feRhs(*feRhsPtr), feLhs(*feLhsPtr),
918  aLpha(1) {
919  ierr = getOptions();
920  CHKERRABORT(PETSC_COMM_WORLD, ierr);
921  }
922 
923  struct OpJacobian
925 
926  const int tAg;
927  boost::shared_ptr<VectorDouble> activeVariablesPtr;
928  boost::shared_ptr<VectorDouble> resultsPtr;
929  boost::shared_ptr<MatrixDouble> jacobianPtr;
931 
932  const double &aLpha;
933 
934  OpJacobian(int tag, const std::string field_name,
935  boost::shared_ptr<VectorDouble> &active_variables_ptr,
936  boost::shared_ptr<VectorDouble> &results_ptr,
937  boost::shared_ptr<MatrixDouble> &jacobian_ptr,
938  bool evaluate_jacobian, const double &alpha)
940  field_name, UserDataOperator::OPCOL),
941  tAg(tag), activeVariablesPtr(active_variables_ptr),
942  resultsPtr(results_ptr), jacobianPtr(jacobian_ptr),
943  evaluateJacobian(evaluate_jacobian), aLpha(alpha) {}
944 
945  MoFEMErrorCode doWork(int side, EntityType type,
948  if (type != MBVERTEX)
950 
955 
956  Tag th0, th1, th2, th3;
958  th1, th2, th3);
959  FTensor::Tensor1<double, 3> t_edge_base0, t_edge_base1;
960  EntityHandle fe_ent = getFEEntityHandle();
961  CHKERR getEdgeFE()->mField.get_moab().tag_get_data(th0, &fe_ent, 1,
962  &t_edge_base0(0));
963  CHKERR getEdgeFE()->mField.get_moab().tag_get_data(th1, &fe_ent, 1,
964  &t_edge_base1(0));
965 
966  VectorInt &indices = data.getIndices();
967 
968  trace_on(tAg);
969 
970  ublas::vector<adouble> lambda_dofs(4);
971  for (int dd = 0; dd != 4; ++dd) {
972  lambda_dofs[dd] <<= (*activeVariablesPtr)[dd];
973  }
974  ublas::vector<adouble> position_dofs(6);
975  for (int dd = 0; dd != 6; ++dd) {
976  position_dofs[dd] <<= (*activeVariablesPtr)[4 + dd];
977  }
978 
980  &position_dofs[0], &position_dofs[1], &position_dofs[2]);
982  &position_dofs[3], &position_dofs[4], &position_dofs[5]);
983 
985  t_tangent(i) = t_node1(i) - t_node0(i);
986  adouble l = sqrt(t_tangent(i) * t_tangent(i));
987  t_tangent(i) /= l;
988 
989  adouble t_dot0, t_dot1;
990  t_dot0 = t_edge_base0(i) * t_tangent(i);
991  t_dot1 = t_edge_base1(i) * t_tangent(i);
992 
993  FTensor::Tensor1<adouble, 3> t_base0, t_base1;
994  t_base0(i) = t_edge_base0(i) - t_dot0 * t_tangent(i);
995  t_base1(i) = t_edge_base1(i) - t_dot1 * t_tangent(i);
996  t_base0(i) /= sqrt(t_base0(i) * t_base0(i));
997  t_base1(i) /= sqrt(t_base1(i) * t_base1(i));
998 
999  auto t_base_fun1 = data.getFTensor0N();
1000  auto t_base_fun2 = data.getFTensor0N();
1001  FTensor::Tensor1<adouble, 3> t_position;
1004  auto t_coord_ref = getFTensor1CoordsAtGaussPts();
1005 
1006  ublas::vector<adouble> c_vec(4);
1007  ublas::vector<adouble> f_vec(6);
1008  c_vec.clear();
1009  f_vec.clear();
1010 
1011  int nb_gauss_pts = data.getN().size1();
1012  int nb_base_functions = data.getN().size2();
1013  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1014 
1015  FTensor::Tensor1<adouble *, 3> t_position_dofs(
1016  &position_dofs[0], &position_dofs[1], &position_dofs[2], 3);
1017  FTensor::Tensor1<adouble *, 2> t_lambda_dof(&lambda_dofs[0],
1018  &lambda_dofs[1], 2);
1019 
1020  t_position(i) = 0;
1021  t_lambda(j) = 0;
1022  for (int bb = 0; bb != nb_base_functions; ++bb) {
1023  t_position(i) += t_base_fun1 * t_position_dofs(i);
1024  t_lambda(j) += t_base_fun1 * t_lambda_dof(j);
1025  ++t_base_fun1;
1026  ++t_position_dofs;
1027  ++t_lambda_dof;
1028  }
1029 
1030  t_delta(i) = t_position(i) - t_coord_ref(i);
1031  adouble dot0 = t_base0(i) * t_delta(i);
1032  adouble dot1 = t_base1(i) * t_delta(i);
1033 
1034  adouble w = getGaussPts()(1, gg) * l * aLpha;
1035  adouble val, val1, val2;
1036  FTensor::Tensor1<adouble *, 2> t_c(&c_vec[0], &c_vec[1], 2);
1037  FTensor::Tensor1<adouble *, 3> t_f(&f_vec[0], &f_vec[1], &f_vec[2], 3);
1038  for (int bb = 0; bb != nb_base_functions; ++bb) {
1039  if (indices[2 * bb] != -1) {
1040  val = w * t_base_fun2;
1041  t_c(N0) += val * dot0;
1042  t_c(N1) += val * dot1;
1043  val1 = val * t_lambda(N0);
1044  val2 = val * t_lambda(N1);
1045  t_f(i) += val1 * t_base0(i) + val2 * t_base1(i);
1046  }
1047  ++t_c;
1048  ++t_f;
1049  ++t_base_fun2;
1050  }
1051 
1052  ++t_coord_ref;
1053  }
1054 
1055  for (int rr = 0; rr != 4; ++rr) {
1056  c_vec[rr] >>= (*resultsPtr)[rr];
1057  }
1058  for (int rr = 0; rr != 6; ++rr) {
1059  f_vec(rr) >>= (*resultsPtr)[4 + rr];
1060  }
1061 
1062  trace_off();
1063 
1064  if (evaluateJacobian) {
1065  double *jac_ptr[4 + 6];
1066  for (int rr = 0; rr != 4 + 6; ++rr) {
1067  jac_ptr[rr] = &(*jacobianPtr)(rr, 0);
1068  }
1069  // play recorder for jacobians
1070  int r =
1071  ::jacobian(tAg, 4 + 6, 4 + 6, &(*activeVariablesPtr)[0], jac_ptr);
1072  if (r < 0) {
1073  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1074  "ADOL-C function evaluation with error");
1075  }
1076  }
1077 
1079  }
1080  };
1081 
1082  MoFEMErrorCode setOperators(int tag, Range edges, Range faces,
1083  const std::string lagrange_multipliers_field_name,
1084  const std::string material_field_name) {
1087  edges, faces);
1088  CHKERR setOperators(tag, lagrange_multipliers_field_name,
1089  material_field_name);
1091  }
1092 
1094  const std::string lagrange_multipliers_field_name,
1095  const std::string material_field_name,
1096  const double *alpha = nullptr) {
1098 
1099  if (alpha) {
1100  aLpha = *alpha;
1101  }
1102 
1103  boost::shared_ptr<VectorDouble> active_variables_ptr(
1104  new VectorDouble(4 + 6));
1105  boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(4 + 6));
1106  boost::shared_ptr<MatrixDouble> jacobian_ptr(
1107  new MatrixDouble(4 + 6, 4 + 6));
1108 
1109  feRhs.getOpPtrVector().clear();
1110  feRhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
1111  lagrange_multipliers_field_name, active_variables_ptr));
1113  material_field_name, active_variables_ptr));
1114  feRhs.getOpPtrVector().push_back(new OpJacobian(
1115  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1116  jacobian_ptr, false, aLpha));
1117  feRhs.getOpPtrVector().push_back(
1118  new OpAssembleRhs<4, 6>(lagrange_multipliers_field_name, results_ptr));
1119  feRhs.getOpPtrVector().push_back(
1120  new OpAssembleRhs<4, 6>(material_field_name, results_ptr));
1121 
1122  // Adding operators to calculate the left hand side
1123  feLhs.getOpPtrVector().clear();
1124  feLhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
1125  lagrange_multipliers_field_name, active_variables_ptr));
1127  material_field_name, active_variables_ptr));
1128  feLhs.getOpPtrVector().push_back(new OpJacobian(
1129  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1130  jacobian_ptr, true, aLpha));
1131  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1132  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1133  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1134  material_field_name, lagrange_multipliers_field_name, jacobian_ptr));
1135  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1136  material_field_name, material_field_name, jacobian_ptr));
1137 
1139  }
1140 
1142  setOperatorsConstrainOnly(int tag, Range edges, Range faces,
1143  const std::string lagrange_multipliers_field_name,
1144  const std::string material_field_name) {
1146 
1148  edges, faces);
1149  CHKERR setOperatorsConstrainOnly(tag, lagrange_multipliers_field_name,
1150  material_field_name);
1152  }
1153 
1156  const std::string lagrange_multipliers_field_name,
1157  const std::string material_field_name) {
1159 
1160  boost::shared_ptr<VectorDouble> active_variables_ptr(
1161  new VectorDouble(4 + 6));
1162  boost::shared_ptr<VectorDouble> results_ptr(new VectorDouble(4 + 6));
1163  boost::shared_ptr<MatrixDouble> jacobian_ptr(
1164  new MatrixDouble(4 + 6, 4 + 6));
1165 
1166  // Adding operators to calculate the left hand side
1167  feLhs.getOpPtrVector().clear();
1168  feLhs.getOpPtrVector().push_back(new OpGetActiveDofsLambda(
1169  lagrange_multipliers_field_name, active_variables_ptr));
1171  material_field_name, active_variables_ptr));
1172  feLhs.getOpPtrVector().push_back(new OpJacobian(
1173  tag, lagrange_multipliers_field_name, active_variables_ptr, results_ptr,
1174  jacobian_ptr, true, aLpha));
1175  feLhs.getOpPtrVector().push_back(new OpAssembleLhs<4, 6>(
1176  lagrange_multipliers_field_name, material_field_name, jacobian_ptr));
1177 
1179  }
1180 };
1181 
1182 #endif // WITH_ADOL_C
1183 
1184 #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)
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
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:506
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
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:128
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:482
boost::shared_ptr< VectorDouble > activeVariablesPtr
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:513
UserDataOperator(const FieldSpace space, const char type=OPLAST, const bool symm=true)
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.
boost::shared_ptr< MyEdgeFE > feRhsPtr
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:90
Vec & snes_f
residual
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
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
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 MoFEMErrorCode saveEdges(moab::Interface &moab, std::string name, Range edges, Range *faces=nullptr)
boost::shared_ptr< VectorDouble > activeVariablesPtr
const FaceElementForcesAndSourcesCoreBase * getFaceFE()
return pointer to Generic Triangle Finite Element object
boost::shared_ptr< MyEdgeFE > feLhsPtr
CHKERRQ(ierr)
constexpr int order
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:66
MoFEMErrorCode setOperators(int tag, const std::string lagrange_multipliers_field_name, const std::string material_field_name, const double *alpha=nullptr)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
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)
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
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:601
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
virtual MoFEMErrorCode getElementOrientation(MoFEM::Interface &m_field, const FEMethod *fe_method_ptr)
boost::shared_ptr< MatrixDouble > jacobianPtr
MoFEMErrorCode preProcess()
function is run at the beginning of loop
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
ForcesAndSourcesCore::UserDataOperator UserDataOperator
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:71
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:87
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)
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:1791
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:412
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
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
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)
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
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:160
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
OpGetActiveDofsLambda(const std::string field_name, boost::shared_ptr< VectorDouble > &active_variables_ptr)