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  ++t_coords;
269  }
270 
271  return v_sdf;
272 }
273 
274 inline MatrixDouble
275 grad_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
276  MatrixDouble &m_spatial_coords,
277  MatrixDouble &m_normals_at_pts, int block_id) {
278 #ifdef PYTHON_SDF
279  if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
280 
281  VectorDouble v_spatial_coords = m_spatial_coords.data();
282  VectorDouble v_normal_at_pts = m_normals_at_pts.data();
283 
284  bp::list python_coords;
285  bp::list python_normals;
286 
287  for (int idx = 0; idx < 3; ++idx) {
288  python_coords.append(
289  convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
290  python_normals.append(
291  convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
292  }
293 
294  np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
295  np::dtype::get_builtin<double>());
296  CHK_MOAB_THROW(sdf_ptr->evalGradSdf(
297  delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
298  bp::extract<np::ndarray>(python_coords[1]),
299  bp::extract<np::ndarray>(python_coords[2]),
300  bp::extract<np::ndarray>(python_normals[0]),
301  bp::extract<np::ndarray>(python_normals[1]),
302  bp::extract<np::ndarray>(python_normals[2]), block_id,
303  np_grad_sdf),
304  "Failed python call");
305 
306  double *grad_ptr = reinterpret_cast<double *>(np_grad_sdf.get_data());
307 
308  MatrixDouble m_grad_sdf;
309  m_grad_sdf.resize(3, nb_gauss_pts, false);
310  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
311  for (int idx = 0; idx < 3; ++idx)
312  m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
313  }
314  return m_grad_sdf;
315  }
316 #endif
317  MatrixDouble m_grad_sdf;
318  m_grad_sdf.resize(3, nb_gauss_pts, false);
319  FTensor::Index<'i', 3> i;
320  FTensor::Tensor1<double, 3> t_grad_sdf_set{0.0, 1.0, 0.0};
321  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
322 
323  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
324  t_grad_sdf(i) = t_grad_sdf_set(i);
325  ++t_grad_sdf;
326  }
327 
328  return m_grad_sdf;
329 }
330 
331 inline MatrixDouble
332 hess_surface_distance_function(double delta_t, double t, int nb_gauss_pts,
333  MatrixDouble &m_spatial_coords,
334  MatrixDouble &m_normals_at_pts, int block_id) {
335 #ifdef PYTHON_SDF
336  if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
337 
338  VectorDouble v_spatial_coords = m_spatial_coords.data();
339  VectorDouble v_normal_at_pts = m_normals_at_pts.data();
340 
341  bp::list python_coords;
342  bp::list python_normals;
343 
344  for (int idx = 0; idx < 3; ++idx) {
345  python_coords.append(
346  convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
347  python_normals.append(
348  convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
349  };
350 
351  np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
352  np::dtype::get_builtin<double>());
353  CHK_MOAB_THROW(sdf_ptr->evalHessSdf(
354  delta_t, t, bp::extract<np::ndarray>(python_coords[0]),
355  bp::extract<np::ndarray>(python_coords[1]),
356  bp::extract<np::ndarray>(python_coords[2]),
357  bp::extract<np::ndarray>(python_normals[0]),
358  bp::extract<np::ndarray>(python_normals[1]),
359  bp::extract<np::ndarray>(python_normals[2]), block_id,
360  np_hess_sdf),
361  "Failed python call");
362 
363  double *hess_ptr = reinterpret_cast<double *>(np_hess_sdf.get_data());
364 
365  MatrixDouble m_hess_sdf;
366  m_hess_sdf.resize(6, nb_gauss_pts, false);
367  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
368  for (int idx = 0; idx < 6; ++idx)
369  m_hess_sdf(idx, gg) =
370  *(hess_ptr + (6 * gg + idx));
371  }
372  return m_hess_sdf;
373  }
374 #endif
375  MatrixDouble m_hess_sdf;
376  m_hess_sdf.resize(6, nb_gauss_pts, false);
377  FTensor::Index<'i', 3> i;
378  FTensor::Index<'j', 3> j;
379  FTensor::Tensor2_symmetric<double, 3> t_hess_sdf_set{0., 0., 0., 0., 0., 0.};
380  auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
381 
382  for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
383  t_hess_sdf(i, j) = t_hess_sdf_set(i, j);
384  ++t_hess_sdf;
385  }
386  return m_hess_sdf;
387 }
388 
389 
390 
391 template <int DIM, IntegrationType I, typename BoundaryEleOp>
392 struct OpAssembleTotalContactTractionImpl;
393 
394 template <int DIM, IntegrationType I, typename BoundaryEleOp>
395 struct OpEvaluateSDFImpl;
396 
397 template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
398 struct OpConstrainBoundaryRhsImpl;
399 
400 template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
401 struct OpConstrainBoundaryLhs_dUImpl;
402 
403 template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
404 struct OpConstrainBoundaryLhs_dTractionImpl;
405 
406 template <typename T1, typename T2, int DIM1, int DIM2>
407 inline auto get_spatial_coords(FTensor::Tensor1<T1, DIM1> &&t_coords,
408  FTensor::Tensor1<T2, DIM2> &&t_disp,
409  size_t nb_gauss_pts) {
410  MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
411  m_spatial_coords.clear();
412  auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
413  FTensor::Index<'i', DIM2> i;
414  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
415  t_spatial_coords(i) = t_coords(i) + t_disp(i);
416  ++t_spatial_coords;
417  ++t_coords;
418  ++t_disp;
419  }
420  return m_spatial_coords;
421 }
422 
423 template <typename T1, int DIM1>
424 inline auto get_normalize_normals(FTensor::Tensor1<T1, DIM1> &&t_normal_at_pts,
425  size_t nb_gauss_pts) {
426  MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
427  m_normals_at_pts.clear();
428  FTensor::Index<'i', DIM1> i;
429  auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
430  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
431  t_set_normal(i) = t_normal_at_pts(i) / t_normal_at_pts.l2();
432  ++t_set_normal;
433  ++t_normal_at_pts;
434  }
435  return m_normals_at_pts;
436 }
437 
438 template <int DIM, typename BoundaryEleOp>
439 struct OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>
440  : public BoundaryEleOp {
441  OpAssembleTotalContactTractionImpl(
442  boost::shared_ptr<CommonData> common_data_ptr, double scale = 1,
443  bool is_axisymmetric = false);
444  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
445 
446 private:
447  boost::shared_ptr<CommonData> commonDataPtr;
448  const double scaleTraction;
449  bool isAxisymmetric;
450 };
451 
452 template <int DIM, typename BoundaryEleOp>
453 struct OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp> : public BoundaryEleOp {
454  OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
455  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
456 
457 private:
458  boost::shared_ptr<CommonData> commonDataPtr;
459 
460  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
461  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
462  grad_surface_distance_function;
463  HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
464  hess_surface_distance_function;
465 };
466 
467 template <int DIM, typename AssemblyBoundaryEleOp>
468 struct OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
469  : public AssemblyBoundaryEleOp {
470  OpConstrainBoundaryRhsImpl(const std::string field_name,
471  boost::shared_ptr<CommonData> common_data_ptr,
472  bool is_axisymmetric = false);
473  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data);
474 
475  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
476  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
477  grad_surface_distance_function;
478 
479 private:
480  boost::shared_ptr<CommonData> commonDataPtr;
481  bool isAxisymmetric;
482 };
483 
484 template <int DIM, typename AssemblyBoundaryEleOp>
485 struct OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
486  : public AssemblyBoundaryEleOp {
487  OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
488  const std::string col_field_name,
489  boost::shared_ptr<CommonData> common_data_ptr,
490  bool is_axisymmetric = false);
491  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
492  EntitiesFieldData::EntData &col_data);
493 
494  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
495  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
496  grad_surface_distance_function;
497  HessSurfaceDistanceFunction hessSurfaceDistanceFunction =
498  hess_surface_distance_function;
499 
500  boost::shared_ptr<CommonData> commonDataPtr;
501  bool isAxisymmetric;
502 };
503 
504 template <int DIM, typename AssemblyBoundaryEleOp>
505 struct OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>
506  : public AssemblyBoundaryEleOp {
507  OpConstrainBoundaryLhs_dTractionImpl(
508  const std::string row_field_name, const std::string col_field_name,
509  boost::shared_ptr<CommonData> common_data_ptr,
510  bool is_axisymmetric = false);
511  MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
512  EntitiesFieldData::EntData &col_data);
513 
514  SurfaceDistanceFunction surfaceDistanceFunction = surface_distance_function;
515  GradSurfaceDistanceFunction gradSurfaceDistanceFunction =
516  grad_surface_distance_function;
517 
518 private:
519  boost::shared_ptr<CommonData> commonDataPtr;
520  bool isAxisymmetric;
521 };
522 
523 template <typename BoundaryEleOp> struct ContactIntegrators {
524 
525  template <int DIM, IntegrationType I>
526  using OpAssembleTotalContactTraction =
527  OpAssembleTotalContactTractionImpl<DIM, I, BoundaryEleOp>;
528 
529  template <int DIM, IntegrationType I>
530  using OpEvaluateSDF = OpEvaluateSDFImpl<DIM, I, BoundaryEleOp>;
531 
532  template <AssemblyType A> struct Assembly {
533 
534  using AssemblyBoundaryEleOp =
535  typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
536 
537  template <int DIM, IntegrationType I>
538  using OpConstrainBoundaryRhs =
539  OpConstrainBoundaryRhsImpl<DIM, I, AssemblyBoundaryEleOp>;
540 
541  template <int DIM, IntegrationType I>
542  using OpConstrainBoundaryLhs_dU =
543  OpConstrainBoundaryLhs_dUImpl<DIM, I, AssemblyBoundaryEleOp>;
544 
545  template <int DIM, IntegrationType I>
546  using OpConstrainBoundaryLhs_dTraction =
547  OpConstrainBoundaryLhs_dTractionImpl<DIM, I, AssemblyBoundaryEleOp>;
548  };
549 };
550 
551 inline double sign(double x) {
552  constexpr auto eps = std::numeric_limits<float>::epsilon();
553  if (std::abs(x) < eps)
554  return 0;
555  else if (x > eps)
556  return 1;
557  else
558  return -1;
559 };
560 
561 inline double w(const double sdf, const double tn) {
562  return sdf - cn_contact * tn;
563 }
564 
565 /**
566  * @brief constrain function
567  *
568  * return 1 if negative sdf or positive tn
569  *
570  * @param sdf signed distance
571  * @param tn traction
572  * @return double
573  */
574 inline double constrain(double sdf, double tn) {
575  const auto s = sign(w(sdf, tn));
576  return (1 - s) / 2;
577 }
578 
579 template <int DIM, typename BoundaryEleOp>
580 OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::
581  OpAssembleTotalContactTractionImpl(
582  boost::shared_ptr<CommonData> common_data_ptr, double scale,
583  bool is_axisymmetric)
584  : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
585  commonDataPtr(common_data_ptr), scaleTraction(scale),
586  isAxisymmetric(is_axisymmetric) {}
587 
588 template <int DIM, typename BoundaryEleOp>
589 MoFEMErrorCode
590 OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
591  int side, EntityType type, EntData &data) {
592  MoFEMFunctionBegin;
593 
594  FTensor::Index<'i', DIM> i;
595  FTensor::Tensor1<double, 3> t_sum_t{0., 0., 0.};
596 
597  auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
598  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
599  auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
600 
601  const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
602  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
603  double jacobian = 1.;
604  if (isAxisymmetric) {
605  jacobian = 2. * M_PI * t_coords(0);
606  }
607  const double alpha = t_w * jacobian * BoundaryEleOp::getMeasure();
608  t_sum_t(i) += alpha * t_traction(i);
609  ++t_w;
610  ++t_traction;
611  ++t_coords;
612  }
613 
614  t_sum_t(i) *= scaleTraction;
615 
616  constexpr int ind[] = {0, 1, 2};
617  CHKERR VecSetValues(commonDataPtr->totalTraction, 3, ind, &t_sum_t(0),
618  ADD_VALUES);
619 
620  MoFEMFunctionReturn(0);
621 }
622 
623 template <int DIM, typename BoundaryEleOp>
624 OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::OpEvaluateSDFImpl(
625  boost::shared_ptr<CommonData> common_data_ptr)
626  : BoundaryEleOp(NOSPACE, BoundaryEleOp::OPSPACE),
627  commonDataPtr(common_data_ptr) {}
628 
629 template <int DIM, typename BoundaryEleOp>
630 MoFEMErrorCode
631 OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::doWork(int side, EntityType type,
632  EntData &data) {
633  MoFEMFunctionBegin;
634 
635  const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
636  auto &sdf_vec = commonDataPtr->sdfVals;
637  auto &grad_mat = commonDataPtr->gradsSdf;
638  auto &hess_mat = commonDataPtr->hessSdf;
639  auto &constraint_vec = commonDataPtr->constraintVals;
640  auto &contactTraction_mat = commonDataPtr->contactTraction;
641 
642  sdf_vec.resize(nb_gauss_pts, false);
643  grad_mat.resize(DIM, nb_gauss_pts, false);
644  hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
645  constraint_vec.resize(nb_gauss_pts, false);
646 
647  auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
648 
649  auto t_sdf = getFTensor0FromVec(sdf_vec);
650  auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
651  auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
652  auto t_constraint = getFTensor0FromVec(constraint_vec);
653 
654  auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
655  auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
656  auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
657 
658  FTensor::Index<'i', DIM> i;
659  FTensor::Index<'j', DIM> j;
660 
661  auto ts_time = BoundaryEleOp::getTStime();
662  auto ts_time_step = BoundaryEleOp::getTStimeStep();
663 
664  auto m_spatial_coords = get_spatial_coords(
665  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
666  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
667  auto m_normals_at_pts = get_normalize_normals(
668  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
669 
670  // placeholder to pass boundary block id to python
671  int block_id = 0;
672 
673  auto v_sdf =
674  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
675  m_spatial_coords, m_normals_at_pts, block_id);
676 
677  auto m_grad_sdf =
678  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
679  m_spatial_coords, m_normals_at_pts, block_id);
680 
681  auto m_hess_sdf =
682  hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
683  m_spatial_coords, m_normals_at_pts, block_id);
684 
685  auto t_sdf_v = getFTensor0FromVec(v_sdf);
686  auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
687  auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
688 
689  auto next = [&]() {
690  ++t_sdf;
691  ++t_sdf_v;
692  ++t_grad_sdf;
693  ++t_grad_sdf_v;
694  ++t_hess_sdf;
695  ++t_hess_sdf_v;
696  ++t_disp;
697  ++t_traction;
698  ++t_constraint;
699  };
700 
701  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
702 
703  auto tn = -t_traction(i) * t_grad_sdf_v(i);
704  auto c = constrain(t_sdf_v, tn);
705 
706  t_sdf = t_sdf_v;
707  t_grad_sdf(i) = t_grad_sdf_v(i);
708  t_hess_sdf(i, j) = t_hess_sdf_v(i, j);
709  t_constraint = c;
710 
711  next();
712  }
713 
714  MoFEMFunctionReturn(0);
715 }
716 
717 template <int DIM, typename AssemblyBoundaryEleOp>
718 OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
719  OpConstrainBoundaryRhsImpl(const std::string field_name,
720  boost::shared_ptr<CommonData> common_data_ptr,
721  bool is_axisymmetric)
722  : AssemblyBoundaryEleOp(field_name, field_name,
723  AssemblyBoundaryEleOp::OPROW),
724  commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {}
725 
726 template <int DIM, typename AssemblyBoundaryEleOp>
727 MoFEMErrorCode
728 OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
729  EntitiesFieldData::EntData &data) {
730  MoFEMFunctionBegin;
731 
732  FTensor::Index<'i', DIM> i;
733  FTensor::Index<'j', DIM> j;
734  FTensor::Index<'k', DIM> k;
735  FTensor::Index<'l', DIM> l;
736 
737  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
738 
739  auto &nf = AssemblyBoundaryEleOp::locF;
740 
741  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
742 
743  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
744  auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
745  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
746  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
747 
748  size_t nb_base_functions = data.getN().size2() / 3;
749  auto t_base = data.getFTensor1N<3>();
750 
751  auto m_spatial_coords = get_spatial_coords(
752  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
753  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
754  auto m_normals_at_pts = get_normalize_normals(
755  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
756 
757  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
758 
759  auto ts_time = AssemblyBoundaryEleOp::getTStime();
760  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
761 
762  // placeholder to pass boundary block id to python
763  int block_id = 0;
764 
765  auto v_sdf =
766  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
767  m_spatial_coords, m_normals_at_pts, block_id);
768 
769  auto m_grad_sdf =
770  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
771  m_spatial_coords, m_normals_at_pts, block_id);
772 
773  auto t_sdf = getFTensor0FromVec(v_sdf);
774  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
775 
776  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
777 
778  auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
779 
780  double jacobian = 1.;
781  if (isAxisymmetric) {
782  jacobian = 2. * M_PI * t_coords(0);
783  }
784  const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
785 
786  auto tn = -t_traction(i) * t_grad_sdf(i);
787  auto c = constrain(t_sdf, tn);
788 
789  FTensor::Tensor2<double, DIM, DIM> t_cP;
790  t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
791  FTensor::Tensor2<double, DIM, DIM> t_cQ;
792  t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
793 
794  FTensor::Tensor1<double, DIM> t_rhs;
795  t_rhs(i) =
796 
797  t_cQ(i, j) * (t_disp(j) - cn_contact * t_traction(j))
798 
799  +
800 
801  t_cP(i, j) * t_disp(j) +
802  c * (t_sdf * t_grad_sdf(i)); // add gap0 displacements
803 
804  size_t bb = 0;
805  for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
806  const double beta = alpha * (t_base(i) * t_normal(i));
807  t_nf(i) -= beta * t_rhs(i);
808 
809  ++t_nf;
810  ++t_base;
811  }
812  for (; bb < nb_base_functions; ++bb)
813  ++t_base;
814 
815  ++t_disp;
816  ++t_traction;
817  ++t_coords;
818  ++t_w;
819  ++t_normal;
820  ++t_sdf;
821  ++t_grad_sdf;
822  }
823 
824  MoFEMFunctionReturn(0);
825 }
826 
827 template <int DIM, typename AssemblyBoundaryEleOp>
828 OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
829  OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
830  const std::string col_field_name,
831  boost::shared_ptr<CommonData> common_data_ptr,
832  bool is_axisymmetric)
833  : AssemblyBoundaryEleOp(row_field_name, col_field_name,
834  AssemblyBoundaryEleOp::OPROWCOL),
835  commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
836  AssemblyBoundaryEleOp::sYmm = false;
837 }
838 
839 template <int DIM, typename AssemblyBoundaryEleOp>
840 MoFEMErrorCode
841 OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
842  EntitiesFieldData::EntData &row_data,
843  EntitiesFieldData::EntData &col_data) {
844  MoFEMFunctionBegin;
845 
846  FTensor::Index<'i', DIM> i;
847  FTensor::Index<'j', DIM> j;
848  FTensor::Index<'k', DIM> k;
849 
850  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
851  auto &locMat = AssemblyBoundaryEleOp::locMat;
852 
853  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
854  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
855  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
856 
857  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
858  auto t_row_base = row_data.getFTensor1N<3>();
859  size_t nb_face_functions = row_data.getN().size2() / 3;
860 
861  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
862 
863  auto m_spatial_coords = get_spatial_coords(
864  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
865  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
866  auto m_normals_at_pts = get_normalize_normals(
867  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
868 
869  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
870 
871  auto ts_time = AssemblyBoundaryEleOp::getTStime();
872  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
873 
874  // placeholder to pass boundary block id to python
875  int block_id = 0;
876 
877  auto v_sdf =
878  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
879  m_spatial_coords, m_normals_at_pts, block_id);
880 
881  auto m_grad_sdf =
882  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
883  m_spatial_coords, m_normals_at_pts, block_id);
884 
885  auto m_hess_sdf =
886  hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
887  m_spatial_coords, m_normals_at_pts, block_id);
888 
889  auto t_sdf = getFTensor0FromVec(v_sdf);
890  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
891  auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
892 
893  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
894 
895  double jacobian = 1.;
896  if (isAxisymmetric) {
897  jacobian = 2. * M_PI * t_coords(0);
898  }
899  const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
900 
901  auto tn = -t_traction(i) * t_grad_sdf(i);
902  auto c = constrain(t_sdf, tn);
903 
904  FTensor::Tensor2<double, DIM, DIM> t_cP;
905  t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
906  FTensor::Tensor2<double, DIM, DIM> t_cQ;
907  t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
908 
909  FTensor::Tensor2<double, DIM, DIM> t_res_dU;
910  t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
911 
912  if (c > 0) {
913  t_res_dU(i, j) +=
914  (c * cn_contact) *
915  (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
916  t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
917  c * t_sdf * t_hess_sdf(i, j);
918  }
919 
920  size_t rr = 0;
921  for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
922 
923  auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
924 
925  const double row_base = t_row_base(i) * t_normal(i);
926 
927  auto t_col_base = col_data.getFTensor0N(gg, 0);
928  for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
929  const double beta = alpha * row_base * t_col_base;
930 
931  t_mat(i, j) -= beta * t_res_dU(i, j);
932 
933  ++t_col_base;
934  ++t_mat;
935  }
936 
937  ++t_row_base;
938  }
939  for (; rr < nb_face_functions; ++rr)
940  ++t_row_base;
941 
942  ++t_traction;
943  ++t_coords;
944  ++t_w;
945  ++t_normal;
946  ++t_sdf;
947  ++t_grad_sdf;
948  ++t_hess_sdf;
949  }
950 
951  MoFEMFunctionReturn(0);
952 }
953 
954 template <int DIM, typename AssemblyBoundaryEleOp>
955 OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
956  OpConstrainBoundaryLhs_dTractionImpl(
957  const std::string row_field_name, const std::string col_field_name,
958  boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric)
959  : AssemblyBoundaryEleOp(row_field_name, col_field_name,
960  AssemblyBoundaryEleOp::OPROWCOL),
961  commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
962  AssemblyBoundaryEleOp::sYmm = false;
963 }
964 
965 template <int DIM, typename AssemblyBoundaryEleOp>
966 MoFEMErrorCode
967 OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
968  iNtegrate(EntitiesFieldData::EntData &row_data,
969  EntitiesFieldData::EntData &col_data) {
970  MoFEMFunctionBegin;
971 
972  FTensor::Index<'i', DIM> i;
973  FTensor::Index<'j', DIM> j;
974  FTensor::Index<'k', DIM> k;
975 
976  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
977  auto &locMat = AssemblyBoundaryEleOp::locMat;
978 
979  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
980  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
981  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
982 
983  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
984  auto t_row_base = row_data.getFTensor1N<3>();
985  size_t nb_face_functions = row_data.getN().size2() / 3;
986 
987  auto m_spatial_coords = get_spatial_coords(
988  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
989  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
990  auto m_normals_at_pts = get_normalize_normals(
991  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
992 
993  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
994 
995  auto ts_time = AssemblyBoundaryEleOp::getTStime();
996  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
997 
998  // placeholder to pass boundary block id to python
999  int block_id = 0;
1000 
1001  auto v_sdf =
1002  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1003  m_spatial_coords, m_normals_at_pts, block_id);
1004 
1005  auto m_grad_sdf =
1006  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1007  m_spatial_coords, m_normals_at_pts, block_id);
1008 
1009  auto t_sdf = getFTensor0FromVec(v_sdf);
1010  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_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_dt;
1029  t_res_dt(i, j) = -cn_contact * t_cQ(i, j);
1030 
1031  size_t rr = 0;
1032  for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
1033 
1034  auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1035  const double row_base = t_row_base(i) * t_normal(i);
1036 
1037  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1038  for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
1039  const double col_base = t_col_base(i) * t_normal(i);
1040  const double beta = alpha * row_base * col_base;
1041 
1042  t_mat(i, j) -= beta * t_res_dt(i, j);
1043 
1044  ++t_col_base;
1045  ++t_mat;
1046  }
1047 
1048  ++t_row_base;
1049  }
1050  for (; rr < nb_face_functions; ++rr)
1051  ++t_row_base;
1052 
1053  ++t_traction;
1054  ++t_coords;
1055  ++t_w;
1056  ++t_normal;
1057  ++t_sdf;
1058  ++t_grad_sdf;
1059  }
1060 
1061  MoFEMFunctionReturn(0);
1062 }
1063 
1064 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
1065 MoFEMErrorCode opFactoryDomainRhs(
1066  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1067  std::string sigma, std::string u, bool is_axisymmetric = false) {
1068  MoFEMFunctionBegin;
1069 
1070  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1071  A>::template LinearForm<I>;
1072  using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
1073  using OpMixDivUCylRhs =
1074  typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
1075 
1076  using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
1077  using OpMixUTimesDivLambdaRhs =
1078  typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
1079  using OpMixUTimesLambdaRhs =
1080  typename B::template OpGradTimesTensor<1, DIM, DIM>;
1081 
1082  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1083  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1084  auto div_stress_ptr = boost::make_shared<MatrixDouble>();
1085  auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
1086 
1087  auto jacobian = [is_axisymmetric](const double r, const double,
1088  const double) {
1089  if (is_axisymmetric)
1090  return 2. * M_PI * r;
1091  else
1092  return 1.;
1093  };
1094 
1095  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1096  u, common_data_ptr->contactDispPtr()));
1097  pip.push_back(
1098  new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
1099 
1100  if (!is_axisymmetric) {
1101  pip.push_back(
1102  new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
1103  } else {
1104  pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
1105  sigma, div_stress_ptr));
1106  }
1107 
1108  pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
1109 
1110  if (!is_axisymmetric) {
1111  pip.push_back(
1112  new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
1113  } else {
1114  pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
1115  jacobian));
1116  }
1117 
1118  pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
1119  pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
1120  pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
1121 
1122  MoFEMFunctionReturn(0);
1123 }
1124 
1125 template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
1126  using OpMixLhs::OpMixLhs;
1127  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1128  EntityType col_type,
1129  EntitiesFieldData::EntData &row_data,
1130  EntitiesFieldData::EntData &col_data) {
1131  MoFEMFunctionBegin;
1132  auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
1133  auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
1134  // Only assemble side which correspond to edge entity on boundary
1135  if (side_fe_entity == side_fe_data) {
1136  CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
1137  col_data);
1138  }
1139  MoFEMFunctionReturn(0);
1140  }
1141 };
1142 
1143 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
1144 MoFEMErrorCode opFactoryBoundaryToDomainLhs(
1145  MoFEM::Interface &m_field,
1146  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1147  std::string fe_domain_name, std::string sigma, std::string u,
1148  std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
1149  bool is_axisymmetric = false) {
1150  MoFEMFunctionBegin;
1151 
1152  using DomainEleOp = typename DomainEle::UserDataOperator;
1153 
1154  auto op_loop_side = new OpLoopSide<DomainEle>(
1155  m_field, fe_domain_name, DIM, Sev::noisy,
1156  boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
1157  pip.push_back(op_loop_side);
1158 
1159  CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
1160  {H1, HDIV}, geom);
1161 
1162  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1163  A>::template BiLinearForm<I>;
1164 
1165  using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
1166  using OpMixDivUCylLhs =
1167  typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
1168  using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
1169 
1170  using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
1171  using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
1172  using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
1173 
1174  auto unity = []() { return 1; };
1175  auto jacobian = [is_axisymmetric](const double r, const double,
1176  const double) {
1177  if (is_axisymmetric)
1178  return 2. * M_PI * r;
1179  else
1180  return 1.;
1181  };
1182 
1183  if (!is_axisymmetric) {
1184  op_loop_side->getOpPtrVector().push_back(
1185  new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
1186  } else {
1187  op_loop_side->getOpPtrVector().push_back(
1188  new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
1189  }
1190  op_loop_side->getOpPtrVector().push_back(
1191  new OpLambdaGraULhsSide(sigma, u, unity, jacobian, true));
1192 
1193  op_loop_side->getSideFEPtr()->getRuleHook = rule;
1194  MoFEMFunctionReturn(0);
1195 }
1196 
1197 template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1198 MoFEMErrorCode opFactoryBoundaryLhs(
1199  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1200  std::string sigma, std::string u, bool is_axisymmetric = false) {
1201  MoFEMFunctionBegin;
1202 
1203  using C = ContactIntegrators<BoundaryEleOp>;
1204 
1205  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1206 
1207  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1208  u, common_data_ptr->contactDispPtr()));
1209  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1210  sigma, common_data_ptr->contactTractionPtr()));
1211  pip.push_back(
1212  new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
1213  DIM, GAUSS>(sigma, u, common_data_ptr, is_axisymmetric));
1214  pip.push_back(new typename C::template Assembly<A>::
1215  template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
1216  sigma, sigma, common_data_ptr, is_axisymmetric));
1217 
1218  MoFEMFunctionReturn(0);
1219 }
1220 
1221 template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1222 MoFEMErrorCode opFactoryBoundaryRhs(
1223  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1224  std::string sigma, std::string u, bool is_axisymmetric = false) {
1225  MoFEMFunctionBegin;
1226 
1227  using C = ContactIntegrators<BoundaryEleOp>;
1228 
1229  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1230 
1231  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1232  u, common_data_ptr->contactDispPtr()));
1233  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1234  sigma, common_data_ptr->contactTractionPtr()));
1235  pip.push_back(
1236  new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
1237  DIM, GAUSS>(sigma, common_data_ptr, is_axisymmetric));
1238 
1239  MoFEMFunctionReturn(0);
1240 }
1241 
1242 template <int DIM, IntegrationType I, typename BoundaryEleOp>
1243 MoFEMErrorCode opFactoryCalculateTraction(
1244  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1245  std::string sigma, bool is_axisymmetric = false) {
1246  MoFEMFunctionBegin;
1247 
1248  using C = ContactIntegrators<BoundaryEleOp>;
1249 
1250  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1251  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1252  sigma, common_data_ptr->contactTractionPtr()));
1253  pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
1254  common_data_ptr, 1. / scale, is_axisymmetric));
1255 
1256  MoFEMFunctionReturn(0);
1257 }
1258 
1259 }; // namespace ContactOps
1260 
1261 #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:407
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:55
OpMixLhs
ContactOps::OpAssembleTotalContactTractionImpl< DIM, GAUSS, BoundaryEleOp >::isAxisymmetric
bool isAxisymmetric
Definition: ContactOps.hpp:449
ContactOps::sign
double sign(double x)
Definition: ContactOps.hpp:551
ContactOps::CommonData::constraintVals
VectorDouble constraintVals
Definition: ContactOps.hpp:25
ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:500
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< double, 3 >
ContactOps::OpConstrainBoundaryRhsImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:480
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:332
ContactOps::w
double w(const double sdf, const double tn)
Definition: ContactOps.hpp:561
ContactOps
Definition: EshelbianContact.hpp:10
ContactOps::constrain
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:574
is_axisymmetric
PetscBool is_axisymmetric
Definition: contact.cpp:90
ContactOps::EntData
EntitiesFieldData::EntData EntData
Definition: EshelbianContact.hpp:12
ContactOps::OpEvaluateSDFImpl
Definition: ContactOps.hpp:395
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
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:458
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
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:1144
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1579
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpBaseImpl::locF
VectorDouble locF
local entity vector
Definition: FormsIntegrators.hpp:251
ContactOps::OpConstrainBoundaryLhs_dTractionImpl
Definition: ContactOps.hpp:404
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:519
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
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:179
BoundaryEleOp
ContactOps::CommonData::totalTraction
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:28
ContactOps::get_normalize_normals
auto get_normalize_normals(FTensor::Tensor1< T1, DIM1 > &&t_normal_at_pts, size_t nb_gauss_pts)
Definition: ContactOps.hpp:424
double
ContactOps::opFactoryBoundaryRhs
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1222
convert.type
type
Definition: convert.py:64
ContactOps::OpAssembleTotalContactTractionImpl< DIM, GAUSS, BoundaryEleOp >::commonDataPtr
boost::shared_ptr< CommonData > commonDataPtr
Definition: ContactOps.hpp:447
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:523
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
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:1065
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
ContactOps::OpConstrainBoundaryLhs_dUImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::isAxisymmetric
bool isAxisymmetric
Definition: ContactOps.hpp:501
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: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:1198
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition: FormsIntegrators.hpp:236
ContactOps::CommonData::constraintPtr
auto constraintPtr()
Definition: ContactOps.hpp:76
ContactOps::OpMixLhsSide
Definition: ContactOps.hpp:1125
ContactOps::OpConstrainBoundaryLhs_dTractionImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::isAxisymmetric
bool isAxisymmetric
Definition: ContactOps.hpp:520
ContactOps::OpConstrainBoundaryRhsImpl< DIM, GAUSS, AssemblyBoundaryEleOp >::isAxisymmetric
bool isAxisymmetric
Definition: ContactOps.hpp:481
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:532
ContactOps::CommonData::getFTensor1TotalTraction
static auto getFTensor1TotalTraction()
Definition: ContactOps.hpp:41
ContactOps::OpConstrainBoundaryRhsImpl
Definition: ContactOps.hpp:398
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
AssemblyBoundaryEleOp
ContactOps::OpAssembleTotalContactTractionImpl< DIM, GAUSS, BoundaryEleOp >::scaleTraction
const double scaleTraction
Definition: ContactOps.hpp:448
sdf_wavy_2d.ind
float ind
Definition: sdf_wavy_2d.py:8
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:1127
ContactOps::OpAssembleTotalContactTractionImpl
Definition: ContactOps.hpp:392
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:275
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:1243
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: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:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
ContactOps::OpConstrainBoundaryLhs_dUImpl
Definition: ContactOps.hpp:401
Generated by Doxygen 1.8.17 and hosted at University of Glasgow