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