v0.14.0
ContactOps.hpp
Go to the documentation of this file.
1 
2 
3 /**
4  * \file ContactOps.hpp
5  * \example ContactOps.hpp
6  */
7 
8 #ifndef __CONTACTOPS_HPP__
9 #define __CONTACTOPS_HPP__
10 
11 namespace ContactOps {
12 
13 //! [Common data]
14 struct CommonData : public boost::enable_shared_from_this<CommonData> {
15  // MatrixDouble contactStress;
16  MatrixDouble contactTraction;
17  MatrixDouble contactDisp;
18 
19  VectorDouble sdfVals; ///< size is equal to number of gauss points on element
20  MatrixDouble gradsSdf; ///< nb of rows is equals to dimension, and nb of cols
21  ///< is equals to number of gauss points on element
22  MatrixDouble hessSdf; ///< nb of rows is equals to nb of element of symmetric
23  ///< matrix, and nb of cols is equals to number of gauss
24  ///< points on element
25  VectorDouble constraintVals;
26 
27  static SmartPetscObj<Vec>
28  totalTraction; // User have to release and create vector when appropiate.
29 
30  static auto createTotalTraction(MoFEM::Interface &m_field) {
31  constexpr int ghosts[] = {0, 1, 2};
32  totalTraction =
33  createGhostVector(m_field.get_comm(),
34 
35  (m_field.get_comm_rank() == 0) ? 3 : 0, 3,
36 
37  (m_field.get_comm_rank() == 0) ? 0 : 3, ghosts);
38  return totalTraction;
39  }
40 
41  static auto getFTensor1TotalTraction() {
42  if (CommonData::totalTraction) {
43  const double *t_ptr;
44  CHK_THROW_MESSAGE(VecGetArrayRead(CommonData::totalTraction, &t_ptr),
45  "get array");
46  FTensor::Tensor1<double, 3> t{t_ptr[0], t_ptr[1], t_ptr[2]};
47  CHK_THROW_MESSAGE(VecRestoreArrayRead(CommonData::totalTraction, &t_ptr),
48  "restore array");
49  return t;
50  } else {
51  return FTensor::Tensor1<double, 3>{0., 0., 0.};
52  }
53  }
54 
55  inline auto contactTractionPtr() {
56  return boost::shared_ptr<MatrixDouble>(shared_from_this(),
57  &contactTraction);
58  }
59 
60  inline auto contactDispPtr() {
61  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &contactDisp);
62  }
63 
64  inline auto sdfPtr() {
65  return boost::shared_ptr<VectorDouble>(shared_from_this(), &sdfVals);
66  }
67 
68  inline auto gradSdfPtr() {
69  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &gradsSdf);
70  }
71 
72  inline auto hessSdfPtr() {
73  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &hessSdf);
74  }
75 
76  inline auto constraintPtr() {
77  return boost::shared_ptr<VectorDouble>(shared_from_this(), &constraintVals);
78  }
79 };
80 
81 SmartPetscObj<Vec> CommonData::totalTraction;
82 
83 //! [Common data]
84 
85 //! [Surface distance function from python]
86 #ifdef PYTHON_SDF
87 struct SDFPython {
88  SDFPython() = default;
89  virtual ~SDFPython() = default;
90 
91  MoFEMErrorCode sdfInit(const std::string py_file) {
92  MoFEMFunctionBegin;
93  try {
94 
95  // create main module
96  auto main_module = bp::import("__main__");
97  mainNamespace = main_module.attr("__dict__");
98  bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
99  // create a reference to python function
100  sdfFun = mainNamespace["sdf"];
101  sdfGradFun = mainNamespace["grad_sdf"];
102  sdfHessFun = mainNamespace["hess_sdf"];
103 
104  } catch (bp::error_already_set const &) {
105  // print all other errors to stderr
106  PyErr_Print();
107  CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
108  }
109  MoFEMFunctionReturn(0);
110  };
111 
112  template <typename T>
113  inline std::vector<T>
114  py_list_to_std_vector(const boost::python::object &iterable) {
115  return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
116  boost::python::stl_input_iterator<T>());
117  }
118 
119  MoFEMErrorCode evalSdf(
120 
121  double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
122  np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
123  np::ndarray &sdf
124 
125  ) {
126  MoFEMFunctionBegin;
127  try {
128 
129  // call python function
130  sdf = bp::extract<np::ndarray>(
131  sdfFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
132 
133  } catch (bp::error_already_set const &) {
134  // print all other errors to stderr
135  PyErr_Print();
136  CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
137  }
138  MoFEMFunctionReturn(0);
139  }
140 
141  MoFEMErrorCode evalGradSdf(
142 
143  double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
144  np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
145  np::ndarray &grad_sdf
146 
147  ) {
148  MoFEMFunctionBegin;
149  try {
150 
151  // call python function
152  grad_sdf = bp::extract<np::ndarray>(
153  sdfGradFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
154 
155  } catch (bp::error_already_set const &) {
156  // print all other errors to stderr
157  PyErr_Print();
158  CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
159  }
160  MoFEMFunctionReturn(0);
161  }
162 
163  MoFEMErrorCode evalHessSdf(
164 
165  double delta_t, double t, np::ndarray x, np::ndarray y, np::ndarray z,
166  np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
167  np::ndarray &hess_sdf
168 
169  ) {
170  MoFEMFunctionBegin;
171  try {
172 
173  // call python function
174  hess_sdf = bp::extract<np::ndarray>(
175  sdfHessFun(delta_t, t, x, y, z, tx, ty, tz, block_id));
176 
177  } catch (bp::error_already_set const &) {
178  // print all other errors to stderr
179  PyErr_Print();
180  CHK_THROW_MESSAGE(MOFEM_OPERATION_UNSUCCESSFUL, "Python error");
181  }
182  MoFEMFunctionReturn(0);
183  }
184 
185 private:
186  bp::object mainNamespace;
187  bp::object sdfFun;
188  bp::object sdfGradFun;
189  bp::object sdfHessFun;
190 };
191 
192 static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
193 
194 inline np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts,
195  int id) {
196  auto dtype = np::dtype::get_builtin<double>();
197  auto size = bp::make_tuple(nb_gauss_pts);
198  auto stride = bp::make_tuple(3 * sizeof(double));
199  return (np::from_data(&data[id], dtype, size, stride, bp::object()));
200 };
201 #endif
202 //! [Surface distance function from python]
203 
204 using SurfaceDistanceFunction = boost::function<VectorDouble(
205  double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
206  MatrixDouble &normals_at_pts, int block_id)>;
207 
208 using GradSurfaceDistanceFunction = boost::function<MatrixDouble(
209  double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
210  MatrixDouble &normals_at_pts, int block_id)>;
211 
212 using HessSurfaceDistanceFunction = boost::function<MatrixDouble(
213  double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords,
214  MatrixDouble &normals_at_pts, int block_id)>;
215 
216 inline VectorDouble surface_distance_function(double delta_t, double t,
217  int nb_gauss_pts,
218  MatrixDouble &m_spatial_coords,
219  MatrixDouble &m_normals_at_pts,
220  int block_id) {
221 
222 #ifdef PYTHON_SDF
223  if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
224 
225  VectorDouble v_spatial_coords = m_spatial_coords.data();
226  VectorDouble v_normal_at_pts = m_normals_at_pts.data();
227 
228  bp::list python_coords;
229  bp::list python_normals;
230 
231  for (int idx = 0; idx < 3; ++idx) {
232  python_coords.append(
233  convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
234  python_normals.append(
235  convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
236  }
237 
238  np::ndarray np_sdf = np::empty(bp::make_tuple(nb_gauss_pts),
239  np::dtype::get_builtin<double>());
240  CHK_MOAB_THROW(sdf_ptr->evalSdf(delta_t, t,
241  bp::extract<np::ndarray>(python_coords[0]),
242  bp::extract<np::ndarray>(python_coords[1]),
243  bp::extract<np::ndarray>(python_coords[2]),
244  bp::extract<np::ndarray>(python_normals[0]),
245  bp::extract<np::ndarray>(python_normals[1]),
246  bp::extract<np::ndarray>(python_normals[2]),
247  block_id, np_sdf),
248  "Failed python call");
249 
250  double *sdf_val_ptr = reinterpret_cast<double *>(np_sdf.get_data());
251 
252  VectorDouble v_sdf;
253  v_sdf.resize(nb_gauss_pts, false);
254 
255  for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
256  v_sdf[gg] = *(sdf_val_ptr + gg);
257 
258  return v_sdf;
259  }
260 #endif
261  VectorDouble v_sdf;
262  v_sdf.resize(nb_gauss_pts, false);
263  m_spatial_coords.resize(3, nb_gauss_pts);
264  auto t_coords = getFTensor1FromMat<3>(m_spatial_coords);
265 
266  for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
267  v_sdf[gg] = t_coords(1) + 0.5;
268 
269  return v_sdf;
270 }
271 
272 inline MatrixDouble
273 grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
274  MatrixDouble &m_spatial_coords,
275  MatrixDouble &m_normals_at_pts, int block_id) {
276 #ifdef PYTHON_SDF
277  if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
278 
279  VectorDouble v_spatial_coords = m_spatial_coords.data();
280  VectorDouble v_normal_at_pts = m_normals_at_pts.data();
281 
282  bp::list python_coords;
283  bp::list python_normals;
284 
285  for (int idx = 0; idx < 3; ++idx) {
286  python_coords.append(
287  convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
288  python_normals.append(
289  convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
290  }
291 
292  np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
293  np::dtype::get_builtin<double>());
294  CHK_MOAB_THROW(sdf_ptr->evalGradSdf(
295  delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
296  bp::extract<np::ndarray>(python_coords[1]),
297  bp::extract<np::ndarray>(python_coords[2]),
298  bp::extract<np::ndarray>(python_normals[0]),
299  bp::extract<np::ndarray>(python_normals[1]),
300  bp::extract<np::ndarray>(python_normals[2]), block_id,
301  np_grad_sdf),
302  "Failed python call");
303 
304  double *grad_ptr = reinterpret_cast<double *>(np_grad_sdf.get_data());
305 
306  MatrixDouble m_grad_sdf;
307  m_grad_sdf.resize(3, nb_gauss_pts, false);
308  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
309  for (int idx = 0; idx < 3; ++idx)
310  m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
311  }
312  return m_grad_sdf;
313  }
314 #endif
315  MatrixDouble m_grad_sdf;
316  m_grad_sdf.resize(3, nb_gauss_pts, false);
317  FTensor::Index<'i', 3> i;
318  FTensor::Tensor1<double, 3> t_grad_sdf_set{0.0, 0.0, 1.0};
319  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
320 
321  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
322  t_grad_sdf(i) = t_grad_sdf_set(i);
323  ++t_grad_sdf;
324  }
325 
326  return m_grad_sdf;
327 }
328 
329 inline MatrixDouble
330 hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
331  MatrixDouble &m_spatial_coords,
332  MatrixDouble &m_normals_at_pts, int block_id) {
333 #ifdef PYTHON_SDF
334  if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
335 
336  VectorDouble v_spatial_coords = m_spatial_coords.data();
337  VectorDouble v_normal_at_pts = m_normals_at_pts.data();
338 
339  bp::list python_coords;
340  bp::list python_normals;
341 
342  for (int idx = 0; idx < 3; ++idx) {
343  python_coords.append(
344  convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
345  python_normals.append(
346  convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
347  };
348 
349  np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
350  np::dtype::get_builtin<double>());
351  CHK_MOAB_THROW(sdf_ptr->evalHessSdf(
352  delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
353  bp::extract<np::ndarray>(python_coords[1]),
354  bp::extract<np::ndarray>(python_coords[2]),
355  bp::extract<np::ndarray>(python_normals[0]),
356  bp::extract<np::ndarray>(python_normals[1]),
357  bp::extract<np::ndarray>(python_normals[2]), block_id,
358  np_hess_sdf),
359  "Failed python call");
360 
361  double *hess_ptr = reinterpret_cast<double *>(np_hess_sdf.get_data());
362 
363  MatrixDouble m_hess_sdf;
364  m_hess_sdf.resize(6, nb_gauss_pts, false);
365  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
366  for (int idx = 0; idx < 6; ++idx)
367  m_hess_sdf(idx, gg) =
368  *(hess_ptr + (6 * gg + idx));
369  }
370  return m_hess_sdf;
371  }
372 #endif
373  MatrixDouble m_hess_sdf;
374  m_hess_sdf.resize(6, nb_gauss_pts, false);
375  FTensor::Index<'i', 3> i;
376  FTensor::Index<'j', 3> j;
377  FTensor::Tensor2_symmetric<double, 3> t_hess_sdf_set{0., 0., 0., 0., 0., 0.};
378  auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
379 
380  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
381  t_hess_sdf(i, j) = t_hess_sdf_set(i, j);
382  ++t_hess_sdf;
383  }
384  return m_hess_sdf;
385 }
386 
387 
388 
389 template <int DIM, IntegrationType I, typename BoundaryEleOp>
390 struct OpAssembleTotalContactTractionImpl;
391 
392 template <int DIM, IntegrationType I, typename BoundaryEleOp>
393 struct OpEvaluateSDFImpl;
394 
395 template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
396 struct OpConstrainBoundaryRhsImpl;
397 
398 template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
399 struct OpConstrainBoundaryLhs_dUImpl;
400 
401 template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
402 struct OpConstrainBoundaryLhs_dTractionImpl;
403 
404 template <typename T1, typename T2, int DIM1, int DIM2>
405 inline auto get_spatial_coords(FTensor::Tensor1<T1, DIM1> &&t_coords,
406  FTensor::Tensor1<T2, DIM2> &&t_disp,
407  size_t nb_gauss_pts) {
408  MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
409  m_spatial_coords.clear();
410  auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
411  FTensor::Index<'i', DIM2> i;
412  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
413  t_spatial_coords(i) = t_coords(i) + t_disp(i);
414  ++t_spatial_coords;
415  ++t_coords;
416  ++t_disp;
417  }
418  return m_spatial_coords;
419 }
420 
421 template <typename T1, int DIM1>
422 inline auto get_normalize_normals(FTensor::Tensor1<T1, DIM1> &&t_normal_at_pts,
423  size_t nb_gauss_pts) {
424  MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
425  m_normals_at_pts.clear();
426  FTensor::Index<'i', DIM1> i;
427  auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
428  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
429  t_set_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
430  ++t_set_normal;
431  ++t_normal_at_pts;
432  }
433  return m_normals_at_pts;
434 }
435 
436 template <int DIM, typename BoundaryEleOp>
437 struct OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>
438  : public BoundaryEleOp {
439  OpAssembleTotalContactTractionImpl(
440  boost::shared_ptr<CommonData> common_data_ptr, double scale = 1,
441  bool is_axisymmetric = false);
442  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
443 
444 private:
445  boost::shared_ptr<CommonData> commonDataPtr;
446  const double scaleTraction;
447  bool isAxisymmetric;
448 };
449 
450 template <int DIM, typename BoundaryEleOp>
451 struct OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp> : public BoundaryEleOp {
452  OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
453  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
454 
455 private:
456  boost::shared_ptr<CommonData> commonDataPtr;
457 
458  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
459  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
460  grad_surface_distance_function;
461  HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
462  hess_surface_distance_function;
463 };
464 
465 template <int DIM, typename AssemblyBoundaryEleOp>
466 struct OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
467  : public AssemblyBoundaryEleOp {
468  OpConstrainBoundaryRhsImpl(const std::string field_name,
469  boost::shared_ptr<CommonData> common_data_ptr,
470  bool is_axisymmetric = false);
471  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
472 
473  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
474  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
475  grad_surface_distance_function;
476 
477 private:
478  boost::shared_ptr<CommonData> commonDataPtr;
479  bool isAxisymmetric;
480 };
481 
482 template <int DIM, typename AssemblyBoundaryEleOp>
483 struct OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
484  : public AssemblyBoundaryEleOp {
485  OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
486  const std::string col_field_name,
487  boost::shared_ptr<CommonData> common_data_ptr,
488  bool is_axisymmetric = false);
489  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
490  EntitiesFieldData::EntData &col_data);
491 
492  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
493  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
494  grad_surface_distance_function;
495  HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
496  hess_surface_distance_function;
497 
498  boost::shared_ptr<CommonData> commonDataPtr;
499  bool isAxisymmetric;
500 };
501 
502 template <int DIM, typename AssemblyBoundaryEleOp>
503 struct OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
504  : public AssemblyBoundaryEleOp {
505  OpConstrainBoundaryLhs_dTractionImpl(
506  const std::string row_field_name, const std::string col_field_name,
507  boost::shared_ptr<CommonData> common_data_ptr,
508  bool is_axisymmetric = false);
509  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
510  EntitiesFieldData::EntData &col_data);
511 
512  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
513  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
514  grad_surface_distance_function;
515 
516 private:
517  boost::shared_ptr<CommonData> commonDataPtr;
518  bool isAxisymmetric;
519 };
520 
521 template <typename BoundaryEleOp> struct ContactIntegrators {
522 
523  template <int DIM, IntegrationType I>
524  using OpAssembleTotalContactTraction =
525  OpAssembleTotalContactTractionImpl<DIM, I, BoundaryEleOp>;
526 
527  template <int DIM, IntegrationType I>
528  using OpEvaluateSDF = OpEvaluateSDFImpl<DIM, I, BoundaryEleOp>;
529 
530  template <AssemblyType A> struct Assembly {
531 
532  using AssemblyBoundaryEleOp =
533  typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
534 
535  template <int DIM, IntegrationType I>
536  using OpConstrainBoundaryRhs =
537  OpConstrainBoundaryRhsImpl<DIM, I, AssemblyBoundaryEleOp>;
538 
539  template <int DIM, IntegrationType I>
540  using OpConstrainBoundaryLhs_dU =
541  OpConstrainBoundaryLhs_dUImpl<DIM, I, AssemblyBoundaryEleOp>;
542 
543  template <int DIM, IntegrationType I>
544  using OpConstrainBoundaryLhs_dTraction =
545  OpConstrainBoundaryLhs_dTractionImpl<DIM, I, AssemblyBoundaryEleOp>;
546  };
547 };
548 
549 inline double sign(double x) {
550  constexpr auto eps = std::numeric_limits<float>::epsilon();
551  if (std::abs(x) < eps)
552  return 0;
553  else if (x > eps)
554  return 1;
555  else
556  return -1;
557 };
558 
559 inline double w(const double sdf, const double tn) {
560  return sdf - cn_contact * tn;
561 }
562 
563 /**
564  * @brief constrain function
565  *
566  * return 1 if negative sdf or positive tn
567  *
568  * @param sdf signed distance
569  * @param tn traction
570  * @return double
571  */
572 inline double constrain(double sdf, double tn) {
573  const auto s = sign(w(sdf, tn));
574  return (1 - s) / 2;
575 }
576 
577 template <int DIM, typename BoundaryEleOp>
578 OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::
579  OpAssembleTotalContactTractionImpl(
580  boost::shared_ptr<CommonData> common_data_ptr, double scale,
581  bool is_axisymmetric)
582  : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
583  commonDataPtr(common_data_ptr), scaleTraction(scale),
584  isAxisymmetric(is_axisymmetric) {}
585 
586 template <int DIM, typename BoundaryEleOp>
587 MoFEMErrorCode
588 OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
589  int side, EntityType type, EntData &data) {
590  MoFEMFunctionBegin;
591 
592  FTensor::Index<'i', DIM> i;
593  FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
594 
595  auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
596  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
597  auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
598 
599  const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
600  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
601  double jacobian = 1.;
602  if (isAxisymmetric) {
603  jacobian = 2. * M_PI * t_coords(0);
604  }
605  const double alpha = t_w * jacobian * BoundaryEleOp::getMeasure();
606  t_sum_t(i) += alpha * t_traction(i);
607  ++t_w;
608  ++t_traction;
609  ++t_coords;
610  }
611 
612  t_sum_t(i) *= scaleTraction;
613 
614  constexpr int ind[] = {0, 1, 2};
615  CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
616  ADD_VALUES);
617 
618  MoFEMFunctionReturn(0);
619 }
620 
621 template <int DIM, typename BoundaryEleOp>
622 OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::OpEvaluateSDFImpl(
623  boost::shared_ptr<CommonData> common_data_ptr)
624  : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
625  commonDataPtr(common_data_ptr) {}
626 
627 template <int DIM, typename BoundaryEleOp>
628 MoFEMErrorCode
629 OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::doWork(int side, EntityType type,
630  EntData &data) {
631  MoFEMFunctionBegin;
632 
633  const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
634  auto &sdf_vec = commonDataPtr->sdfVals;
635  auto &grad_mat = commonDataPtr->gradsSdf;
636  auto &hess_mat = commonDataPtr->hessSdf;
637  auto &constraint_vec = commonDataPtr->constraintVals;
638  auto &contactTraction_mat = commonDataPtr->contactTraction;
639 
640  sdf_vec.resize(nb_gauss_pts, false);
641  grad_mat.resize(DIM, nb_gauss_pts, false);
642  hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
643  constraint_vec.resize(nb_gauss_pts, false);
644 
645  auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
646 
647  auto t_sdf = getFTensor0FromVec(sdf_vec);
648  auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
649  auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
650  auto t_constraint = getFTensor0FromVec(constraint_vec);
651 
652  auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
653  auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
654  auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
655 
656  FTensor::Index<'i', DIM> i;
657  FTensor::Index<'j', DIM> j;
658 
659  auto ts_time = BoundaryEleOp::getTStime();
660  auto ts_time_step = BoundaryEleOp::getTStimeStep();
661 
662  auto m_spatial_coords = get_spatial_coords(
663  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
664  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
665  auto m_normals_at_pts = get_normalize_normals(
666  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
667 
668  // placeholder to pass boundary block id to python
669  int block_id = 0;
670 
671  auto v_sdf =
672  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
673  m_spatial_coords, m_normals_at_pts, block_id);
674 
675  auto m_grad_sdf =
676  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
677  m_spatial_coords, m_normals_at_pts, block_id);
678 
679  auto m_hess_sdf =
680  hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
681  m_spatial_coords, m_normals_at_pts, block_id);
682 
683  auto t_sdf_v = getFTensor0FromVec(v_sdf);
684  auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
685  auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
686 
687  auto next = [&]() {
688  ++t_sdf;
689  ++t_sdf_v;
690  ++t_grad_sdf;
691  ++t_grad_sdf_v;
692  ++t_hess_sdf;
693  ++t_hess_sdf_v;
694  ++t_disp;
695  ++t_traction;
696  ++t_constraint;
697  };
698 
699  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
700 
701  auto tn = -t_traction(i) * t_grad_sdf_v(i);
702  auto c = constrain(t_sdf_v, tn);
703 
704  t_sdf = t_sdf_v;
705  t_grad_sdf(i) = t_grad_sdf_v(i);
706  t_hess_sdf(i, j) = t_hess_sdf_v(i, j);
707  t_constraint = c;
708 
709  next();
710  }
711 
712  MoFEMFunctionReturn(0);
713 }
714 
715 template <int DIM, typename AssemblyBoundaryEleOp>
716 OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
717  OpConstrainBoundaryRhsImpl(const std::string field_name,
718  boost::shared_ptr<CommonData> common_data_ptr,
719  bool is_axisymmetric)
720  : AssemblyBoundaryEleOp(field_name, field_name,
721  AssemblyBoundaryEleOp::OPROW),
722  commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {}
723 
724 template <int DIM, typename AssemblyBoundaryEleOp>
725 MoFEMErrorCode
726 OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
727  EntitiesFieldData::EntData &data) {
728  MoFEMFunctionBegin;
729 
730  FTensor::Index<'i', DIM> i;
731  FTensor::Index<'j', DIM> j;
732  FTensor::Index<'k', DIM> k;
733  FTensor::Index<'l', DIM> l;
734 
735  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
736 
737  auto &nf = AssemblyBoundaryEleOp::locF;
738 
739  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
740 
741  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
742  auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
743  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
744  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
745 
746  size_t nb_base_functions = data.getN().size2() / 3;
747  auto t_base = data.getFTensor1N<3>();
748 
749  auto m_spatial_coords = get_spatial_coords(
750  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
751  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
752  auto m_normals_at_pts = get_normalize_normals(
753  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
754 
755  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
756 
757  auto ts_time = AssemblyBoundaryEleOp::getTStime();
758  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
759 
760  // placeholder to pass boundary block id to python
761  int block_id = 0;
762 
763  auto v_sdf =
764  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
765  m_spatial_coords, m_normals_at_pts, block_id);
766 
767  auto m_grad_sdf =
768  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
769  m_spatial_coords, m_normals_at_pts, block_id);
770 
771  auto t_sdf = getFTensor0FromVec(v_sdf);
772  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
773 
774  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
775 
776  auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
777 
778  double jacobian = 1.;
779  if (isAxisymmetric) {
780  jacobian = 2. * M_PI * t_coords(0);
781  }
782  const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
783 
784  auto tn = -t_traction(i) * t_grad_sdf(i);
785  auto c = constrain(t_sdf, tn);
786 
787  FTensor::Tensor2<double, DIM, DIM> t_cP;
788  t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
789  FTensor::Tensor2<double, DIM, DIM> t_cQ;
790  t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
791 
792  FTensor::Tensor1<double, DIM> t_rhs;
793  t_rhs(i) =
794 
795  t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
796 
797  +
798 
799  t_cP(i, j) * t_disp(j) +
800  c * (t_sdf * t_grad_sdf(i)); // add gap0 displacements
801 
802  size_t bb = 0;
803  for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
804  const double beta = alpha * (t_base(i) * t_normal(i));
805  t_nf(i) -= beta * t_rhs(i);
806 
807  ++t_nf;
808  ++t_base;
809  }
810  for (; bb < nb_base_functions; ++bb)
811  ++t_base;
812 
813  ++t_disp;
814  ++t_traction;
815  ++t_coords;
816  ++t_w;
817  ++t_normal;
818  ++t_sdf;
819  ++t_grad_sdf;
820  }
821 
822  MoFEMFunctionReturn(0);
823 }
824 
825 template <int DIM, typename AssemblyBoundaryEleOp>
826 OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
827  OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
828  const std::string col_field_name,
829  boost::shared_ptr<CommonData> common_data_ptr,
830  bool is_axisymmetric)
831  : AssemblyBoundaryEleOp(row_field_name, col_field_name,
832  AssemblyBoundaryEleOp::OPROWCOL),
833  commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
834  AssemblyBoundaryEleOp::sYmm = false;
835 }
836 
837 template <int DIM, typename AssemblyBoundaryEleOp>
838 MoFEMErrorCode
839 OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
840  EntitiesFieldData::EntData &row_data,
841  EntitiesFieldData::EntData &col_data) {
842  MoFEMFunctionBegin;
843 
844  FTensor::Index<'i', DIM> i;
845  FTensor::Index<'j', DIM> j;
846  FTensor::Index<'k', DIM> k;
847 
848  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
849  auto &locMat = AssemblyBoundaryEleOp::locMat;
850 
851  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
852  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
853  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
854 
855  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
856  auto t_row_base = row_data.getFTensor1N<3>();
857  size_t nb_face_functions = row_data.getN().size2() / 3;
858 
859  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
860 
861  auto m_spatial_coords = get_spatial_coords(
862  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
863  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
864  auto m_normals_at_pts = get_normalize_normals(
865  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
866 
867  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
868 
869  auto ts_time = AssemblyBoundaryEleOp::getTStime();
870  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
871 
872  // placeholder to pass boundary block id to python
873  int block_id = 0;
874 
875  auto v_sdf =
876  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
877  m_spatial_coords, m_normals_at_pts, block_id);
878 
879  auto m_grad_sdf =
880  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
881  m_spatial_coords, m_normals_at_pts, block_id);
882 
883  auto m_hess_sdf =
884  hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
885  m_spatial_coords, m_normals_at_pts, block_id);
886 
887  auto t_sdf = getFTensor0FromVec(v_sdf);
888  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
889  auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
890 
891  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
892 
893  double jacobian = 1.;
894  if (isAxisymmetric) {
895  jacobian = 2. * M_PI * t_coords(0);
896  }
897  const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
898 
899  auto tn = -t_traction(i) * t_grad_sdf(i);
900  auto c = constrain(t_sdf, tn);
901 
902  FTensor::Tensor2<double, DIM, DIM> t_cP;
903  t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
904  FTensor::Tensor2<double, DIM, DIM> t_cQ;
905  t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
906 
907  FTensor::Tensor2<double, DIM, DIM> t_res_dU;
908  t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
909 
910  if (c > 0) {
911  t_res_dU(i, j) +=
912  (c * cn_contact) *
913  (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
914  t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
915  c * t_sdf * t_hess_sdf(i, j);
916  }
917 
918  size_t rr = 0;
919  for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
920 
921  auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
922 
923  const double row_base = t_row_base(i) * t_normal(i);
924 
925  auto t_col_base = col_data.getFTensor0N(gg, 0);
926  for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
927  const double beta = alpha * row_base * t_col_base;
928 
929  t_mat(i, j) -= beta * t_res_dU(i, j);
930 
931  ++t_col_base;
932  ++t_mat;
933  }
934 
935  ++t_row_base;
936  }
937  for (; rr < nb_face_functions; ++rr)
938  ++t_row_base;
939 
940  ++t_traction;
941  ++t_coords;
942  ++t_w;
943  ++t_normal;
944  ++t_sdf;
945  ++t_grad_sdf;
946  ++t_hess_sdf;
947  }
948 
949  MoFEMFunctionReturn(0);
950 }
951 
952 template <int DIM, typename AssemblyBoundaryEleOp>
953 OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
954  OpConstrainBoundaryLhs_dTractionImpl(
955  const std::string row_field_name, const std::string col_field_name,
956  boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric)
957  : AssemblyBoundaryEleOp(row_field_name, col_field_name,
958  AssemblyBoundaryEleOp::OPROWCOL),
959  commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
960  AssemblyBoundaryEleOp::sYmm = false;
961 }
962 
963 template <int DIM, typename AssemblyBoundaryEleOp>
964 MoFEMErrorCode
965 OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
966  iNtegrate(EntitiesFieldData::EntData &row_data,
967  EntitiesFieldData::EntData &col_data) {
968  MoFEMFunctionBegin;
969 
970  FTensor::Index<'i', DIM> i;
971  FTensor::Index<'j', DIM> j;
972  FTensor::Index<'k', DIM> k;
973 
974  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
975  auto &locMat = AssemblyBoundaryEleOp::locMat;
976 
977  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
978  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
979  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
980 
981  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
982  auto t_row_base = row_data.getFTensor1N<3>();
983  size_t nb_face_functions = row_data.getN().size2() / 3;
984 
985  auto m_spatial_coords = get_spatial_coords(
986  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
987  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
988  auto m_normals_at_pts = get_normalize_normals(
989  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
990 
991  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
992 
993  auto ts_time = AssemblyBoundaryEleOp::getTStime();
994  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
995 
996  // placeholder to pass boundary block id to python
997  int block_id = 0;
998 
999  auto v_sdf =
1000  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1001  m_spatial_coords, m_normals_at_pts, block_id);
1002 
1003  auto m_grad_sdf =
1004  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1005  m_spatial_coords, m_normals_at_pts, block_id);
1006 
1007  auto t_sdf = getFTensor0FromVec(v_sdf);
1008  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
1009 
1010  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1011 
1012  double jacobian = 1.;
1013  if (isAxisymmetric) {
1014  jacobian = 2. * M_PI * t_coords(0);
1015  }
1016  const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
1017 
1018  auto tn = -t_traction(i) * t_grad_sdf(i);
1019  auto c = constrain(t_sdf, tn);
1020 
1021  FTensor::Tensor2<double, DIM, DIM> t_cP;
1022  t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
1023  FTensor::Tensor2<double, DIM, DIM> t_cQ;
1024  t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
1025 
1026  FTensor::Tensor2<double, DIM, DIM> t_res_dt;
1027  t_res_dt(i, j) = -cn_contact * t_cQ(i, j);
1028 
1029  size_t rr = 0;
1030  for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
1031 
1032  auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1033  const double row_base = t_row_base(i) * t_normal(i);
1034 
1035  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1036  for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
1037  const double col_base = t_col_base(i) * t_normal(i);
1038  const double beta = alpha * row_base * col_base;
1039 
1040  t_mat(i, j) -= beta * t_res_dt(i, j);
1041 
1042  ++t_col_base;
1043  ++t_mat;
1044  }
1045 
1046  ++t_row_base;
1047  }
1048  for (; rr < nb_face_functions; ++rr)
1049  ++t_row_base;
1050 
1051  ++t_traction;
1052  ++t_coords;
1053  ++t_w;
1054  ++t_normal;
1055  ++t_sdf;
1056  ++t_grad_sdf;
1057  }
1058 
1059  MoFEMFunctionReturn(0);
1060 }
1061 
1062 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
1063 MoFEMErrorCode opFactoryDomainRhs(
1064  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1065  std::string sigma, std::string u, bool is_axisymmetric = false) {
1066  MoFEMFunctionBegin;
1067 
1068  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1069  A>::template LinearForm<I>;
1070  using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
1071  using OpMixDivUCylRhs =
1072  typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
1073 
1074  using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
1075  using OpMixUTimesDivLambdaRhs =
1076  typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
1077  using OpMixUTimesLambdaRhs =
1078  typename B::template OpGradTimesTensor<1, DIM, DIM>;
1079 
1080  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1081  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1082  auto div_stress_ptr = boost::make_shared<MatrixDouble>();
1083  auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
1084 
1085  auto jacobian = [is_axisymmetric](const double r, const double,
1086  const double) {
1087  if (is_axisymmetric)
1088  return 2. * M_PI * r;
1089  else
1090  return 1.;
1091  };
1092 
1093  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1094  u, common_data_ptr->contactDispPtr()));
1095  pip.push_back(
1096  new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
1097 
1098  if (!is_axisymmetric) {
1099  pip.push_back(
1100  new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
1101  } else {
1102  pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
1103  sigma, div_stress_ptr));
1104  }
1105 
1106  pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
1107 
1108  if (!is_axisymmetric) {
1109  pip.push_back(
1110  new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
1111  } else {
1112  pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
1113  jacobian));
1114  }
1115 
1116  pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
1117  pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
1118  pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
1119 
1120  MoFEMFunctionReturn(0);
1121 }
1122 
1123 template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
1124  using OpMixLhs::OpMixLhs;
1125  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1126  EntityType col_type,
1127  EntitiesFieldData::EntData &row_data,
1128  EntitiesFieldData::EntData &col_data) {
1129  MoFEMFunctionBegin;
1130  auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
1131  auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
1132  // Only assemble side which correspond to edge entity on boundary
1133  if (side_fe_entity == side_fe_data) {
1134  CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
1135  col_data);
1136  }
1137  MoFEMFunctionReturn(0);
1138  }
1139 };
1140 
1141 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
1142 MoFEMErrorCode opFactoryBoundaryToDomainLhs(
1143  MoFEM::Interface &m_field,
1144  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1145  std::string fe_domain_name, std::string sigma, std::string u,
1146  std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
1147  bool is_axisymmetric = false) {
1148  MoFEMFunctionBegin;
1149 
1150  using DomainEleOp = typename DomainEle::UserDataOperator;
1151 
1152  auto op_loop_side = new OpLoopSide<DomainEle>(
1153  m_field, fe_domain_name, DIM, Sev::noisy,
1154  boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
1155  pip.push_back(op_loop_side);
1156 
1157  CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
1158  {H1, HDIV}, geom);
1159 
1160  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1161  A>::template BiLinearForm<I>;
1162 
1163  using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
1164  using OpMixDivUCylLhs =
1165  typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
1166  using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
1167 
1168  using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
1169  using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
1170  using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
1171 
1172  auto unity = []() { return 1; };
1173  auto jacobian = [is_axisymmetric](const double r, const double,
1174  const double) {
1175  if (is_axisymmetric)
1176  return 2. * M_PI * r;
1177  else
1178  return 1.;
1179  };
1180 
1181  if (!is_axisymmetric) {
1182  op_loop_side->getOpPtrVector().push_back(
1183  new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
1184  } else {
1185  op_loop_side->getOpPtrVector().push_back(
1186  new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
1187  }
1188  op_loop_side->getOpPtrVector().push_back(
1189  new OpLambdaGraULhsSide(sigma, u, unity, jacobian, true));
1190 
1191  op_loop_side->getSideFEPtr()->getRuleHook = rule;
1192  MoFEMFunctionReturn(0);
1193 }
1194 
1195 template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1196 MoFEMErrorCode opFactoryBoundaryLhs(
1197  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1198  std::string sigma, std::string u, bool is_axisymmetric = false) {
1199  MoFEMFunctionBegin;
1200 
1201  using C = ContactIntegrators<BoundaryEleOp>;
1202 
1203  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1204 
1205  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1206  u, common_data_ptr->contactDispPtr()));
1207  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1208  sigma, common_data_ptr->contactTractionPtr()));
1209  pip.push_back(
1210  new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
1211  DIM, GAUSS>(sigma, u, common_data_ptr, is_axisymmetric));
1212  pip.push_back(new typename C::template Assembly<A>::
1213  template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
1214  sigma, sigma, common_data_ptr, is_axisymmetric));
1215 
1216  MoFEMFunctionReturn(0);
1217 }
1218 
1219 template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1220 MoFEMErrorCode opFactoryBoundaryRhs(
1221  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1222  std::string sigma, std::string u, bool is_axisymmetric = false) {
1223  MoFEMFunctionBegin;
1224 
1225  using C = ContactIntegrators<BoundaryEleOp>;
1226 
1227  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1228 
1229  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1230  u, common_data_ptr->contactDispPtr()));
1231  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1232  sigma, common_data_ptr->contactTractionPtr()));
1233  pip.push_back(
1234  new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
1235  DIM, GAUSS>(sigma, common_data_ptr, is_axisymmetric));
1236 
1237  MoFEMFunctionReturn(0);
1238 }
1239 
1240 template <int DIM, IntegrationType I, typename BoundaryEleOp>
1241 MoFEMErrorCode opFactoryCalculateTraction(
1242  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1243  std::string sigma, bool is_axisymmetric = false) {
1244  MoFEMFunctionBegin;
1245 
1246  using C = ContactIntegrators<BoundaryEleOp>;
1247 
1248  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1249  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1250  sigma, common_data_ptr->contactTractionPtr()));
1251  pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
1252  common_data_ptr, 1. / scale, is_axisymmetric));
1253 
1254  MoFEMFunctionReturn(0);
1255 }
1256 
1257 }; // namespace ContactOps
1258 
1259 #endif // __CONTACTOPS_HPP__
ContactOps::get_spatial_coords
auto get_spatial_coords(FTensor::Tensor1< T1, DIM1 > &&t_coords, FTensor::Tensor1< T2, DIM2 > &&t_disp, size_t nb_gauss_pts)
Definition: ContactOps.hpp:405
NOSPACE
@ NOSPACE
Definition: definitions.h:83
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
ContactOps::CommonData::contactTractionPtr
auto contactTractionPtr()
Definition: ContactOps.hpp:55
OpMixLhs
ContactOps::OpAssembleTotalContactTractionImpl< DIM, GAUSS, BoundaryEleOp >::isAxisymmetric
bool isAxisymmetric
Definition: ContactOps.hpp:447
ContactOps::sign
double sign(double x)
Definition: ContactOps.hpp:549
ContactOps::CommonData::constraintVals
VectorDouble constraintVals
Definition: ContactOps.hpp:25
ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:498
ContactOps::CommonData::contactDisp
MatrixDouble contactDisp
Definition: ContactOps.hpp:17
sdf.hess_sdf
def hess_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:19
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:239
FTensor::Tensor1< double, 3 >
ContactOps::OpConstrainBoundaryRhsImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:478
ContactOps::hess_surface_distance_function
MatrixDouble hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:330
ContactOps::w
double w(const double sdf, const double tn)
Definition: ContactOps.hpp:559
ContactOps
Definition: EshelbianContact.hpp:10
ContactOps::constrain
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:572
is_axisymmetric
PetscBool is_axisymmetric
Definition: contact.cpp:137
ContactOps::EntData
EntitiesFieldData::EntData EntData
Definition: EshelbianContact.hpp:12
ContactOps::OpEvaluateSDFImpl
Definition: ContactOps.hpp:393
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:596
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
sdf
Definition: sdf.py:1
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
ContactOps::SurfaceDistanceFunction
boost::function< VectorDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> SurfaceDistanceFunction
[Common data]
Definition: ContactOps.hpp:206
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
ContactOps::HessSurfaceDistanceFunction
boost::function< MatrixDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> HessSurfaceDistanceFunction
Definition: ContactOps.hpp:214
ContactOps::OpEvaluateSDFImpl< DIM, GAUSS, BoundaryEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:456
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
MoFEM::ForcesAndSourcesCore::UserDataOperator::getTStimeStep
double getTStimeStep() const
Definition: ForcesAndSourcesCore.hpp:1207
ContactOps::opFactoryBoundaryToDomainLhs
MoFEMErrorCode opFactoryBoundaryToDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string fe_domain_name, std::string sigma, std::string u, std::string geom, ForcesAndSourcesCore::RuleHookFun rule, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1142
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1239
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::ForcesAndSourcesCore::UserDataOperator::getMeasure
double getMeasure() const
get measure of element
Definition: ForcesAndSourcesCore.hpp:1274
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1235
MoFEM::OpBaseImpl::locF
VectorDouble locF
local entity vector
Definition: FormsIntegrators.hpp:241
ContactOps::OpConstrainBoundaryLhs_dTractionImpl
Definition: ContactOps.hpp:402
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
ContactOps::OpConstrainBoundaryLhs_dTractionImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:517
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
MoFEM::createGhostVector
auto createGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
Definition: PetscSmartObj.hpp:175
ContactOps::CommonData::totalTraction
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:28
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
ContactOps::get_normalize_normals
auto get_normalize_normals(FTensor::Tensor1< T1, DIM1 > &&t_normal_at_pts, size_t nb_gauss_pts)
Definition: ContactOps.hpp:422
double
ContactOps::opFactoryBoundaryRhs
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1220
convert.type
type
Definition: convert.py:64
ContactOps::OpAssembleTotalContactTractionImpl< DIM, GAUSS, BoundaryEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:445
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
ContactOps::CommonData::hessSdf
MatrixDouble hessSdf
Definition: ContactOps.hpp:22
ContactOps::ContactIntegrators
Definition: ContactOps.hpp:521
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:128
ContactOps::CommonData::gradSdfPtr
auto gradSdfPtr()
Definition: ContactOps.hpp:68
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
ContactOps::CommonData::sdfVals
VectorDouble sdfVals
size is equal to number of gauss points on element
Definition: ContactOps.hpp:19
ContactOps::CommonData::gradsSdf
MatrixDouble gradsSdf
Definition: ContactOps.hpp:20
ContactOps::opFactoryDomainRhs
MoFEMErrorCode opFactoryDomainRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1063
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::isAxisymmetric
bool isAxisymmetric
Definition: ContactOps.hpp:499
ContactOps::CommonData::sdfPtr
auto sdfPtr()
Definition: ContactOps.hpp:64
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
BiLinearForm
ContactOps::CommonData
[Common data]
Definition: ContactOps.hpp:14
ContactOps::CommonData::contactDispPtr
auto contactDispPtr()
Definition: ContactOps.hpp:60
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
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 >
ContactOps::CommonData::createTotalTraction
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:30
ContactOps::cn_contact
double cn_contact
Definition: EshelbianContact.hpp:19
ContactOps::GradSurfaceDistanceFunction
boost::function< MatrixDouble(double delta_t, double t, int nb_gauss_pts, MatrixDouble &spatial_coords, MatrixDouble &normals_at_pts, int block_id)> GradSurfaceDistanceFunction
Definition: ContactOps.hpp:210
MoFEM::OpBaseImpl::nbCols
int nbCols
number if dof on column
Definition: FormsIntegrators.hpp:227
DomainEleOp
ContactOps::opFactoryBoundaryLhs
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1196
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition: FormsIntegrators.hpp:226
ContactOps::CommonData::constraintPtr
auto constraintPtr()
Definition: ContactOps.hpp:76
ContactOps::OpMixLhsSide
Definition: ContactOps.hpp:1123
ContactOps::OpConstrainBoundaryLhs_dTractionImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::isAxisymmetric
bool isAxisymmetric
Definition: ContactOps.hpp:518
ContactOps::OpConstrainBoundaryRhsImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::isAxisymmetric
bool isAxisymmetric
Definition: ContactOps.hpp:479
ContactOps::CommonData::contactTraction
MatrixDouble contactTraction
Definition: ContactOps.hpp:16
FTensor::kronecker_delta
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
Definition: Kronecker_Delta.hpp:81
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
ContactOps::ContactIntegrators::Assembly
Definition: ContactOps.hpp:530
ContactOps::CommonData::getFTensor1TotalTraction
static auto getFTensor1TotalTraction()
Definition: ContactOps.hpp:41
ContactOps::OpConstrainBoundaryRhsImpl
Definition: ContactOps.hpp:396
ContactOps::surface_distance_function
VectorDouble surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:216
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getFTensor1NormalsAtGaussPts
auto getFTensor1NormalsAtGaussPts()
get normal at integration points
Definition: FaceElementForcesAndSourcesCore.hpp:316
AssemblyBoundaryEleOp
ContactOps::OpAssembleTotalContactTractionImpl< DIM, GAUSS, BoundaryEleOp >::scaleTraction
const double scaleTraction
Definition: ContactOps.hpp:446
ContactOps::OpMixLhsSide::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: ContactOps.hpp:1125
ContactOps::OpAssembleTotalContactTractionImpl
Definition: ContactOps.hpp:390
ContactOps::grad_surface_distance_function
MatrixDouble grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts, MatrixDouble &m_spatial_coords, MatrixDouble &m_normals_at_pts, int block_id)
Definition: ContactOps.hpp:273
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
ContactOps::opFactoryCalculateTraction
MoFEMErrorCode opFactoryCalculateTraction(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1241
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
ContactOps::CommonData::hessSdfPtr
auto hessSdfPtr()
Definition: ContactOps.hpp:72
sdf.grad_sdf
def grad_sdf(t, x, y, z, tx, ty, tz)
Definition: sdf.py:15
OpCalculateVectorFieldGradient
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
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:1268
ContactOps::OpConstrainBoundaryLhs_dUImpl
Definition: ContactOps.hpp:399
MoFEM::ForcesAndSourcesCore::UserDataOperator::getTStime
double getTStime() const
Definition: ForcesAndSourcesCore.hpp:1199
Generated by Doxygen 1.8.17 and hosted at University of Glasgow