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  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
983 
984  auto m_spatial_coords = get_spatial_coords(
985  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
986  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
987  auto m_normals_at_pts = get_normalize_normals(
988  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
989 
990  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
991 
992  auto ts_time = AssemblyBoundaryEleOp::getTStime();
993  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
994 
995  // placeholder to pass boundary block id to python
996  int block_id = 0;
997 
998  auto v_sdf =
999  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1000  m_spatial_coords, m_normals_at_pts, block_id);
1001 
1002  auto m_grad_sdf =
1003  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1004  m_spatial_coords, m_normals_at_pts, block_id);
1005 
1006  auto m_hess_sdf =
1007  hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1008  m_spatial_coords, m_normals_at_pts, block_id);
1009 
1010  auto t_sdf = getFTensor0FromVec(v_sdf);
1011  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
1012  auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
1013 
1014  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1015 
1016  double jacobian = 1.;
1017  if (isAxisymmetric) {
1018  jacobian = 2. * M_PI * t_coords(0);
1019  }
1020  const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
1021 
1022  auto tn = -t_traction(i) * t_grad_sdf(i);
1023  auto c = constrain(t_sdf, tn);
1024 
1025  FTensor::Tensor2<double, DIM, DIM> t_cP;
1026  t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
1027  FTensor::Tensor2<double, DIM, DIM> t_cQ;
1028  t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
1029 
1030  FTensor::Tensor2<double, DIM, DIM> t_res_dU;
1031  t_res_dU(i, j) = kronecker_delta(i, j) + t_cP(i, j);
1032 
1033  if (c > 0) {
1034  t_res_dU(i, j) +=
1035  (c * cn_contact) *
1036  (t_hess_sdf(i, j) * (t_grad_sdf(k) * t_traction(k)) +
1037  t_grad_sdf(i) * t_hess_sdf(k, j) * t_traction(k)) +
1038  c * t_sdf * t_hess_sdf(i, j);
1039  }
1040 
1041  size_t rr = 0;
1042  for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
1043 
1044  auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1045 
1046  const double row_base = t_row_base(i) * t_normal(i);
1047 
1048  auto t_col_base = col_data.getFTensor0N(gg, 0);
1049  for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
1050  const double beta = alpha * row_base * t_col_base;
1051 
1052  t_mat(i, j) -= beta * t_res_dU(i, j);
1053 
1054  ++t_col_base;
1055  ++t_mat;
1056  }
1057 
1058  ++t_row_base;
1059  }
1060  for (; rr < nb_face_functions; ++rr)
1061  ++t_row_base;
1062 
1063  ++t_traction;
1064  ++t_coords;
1065  ++t_w;
1066  ++t_normal;
1067  ++t_sdf;
1068  ++t_grad_sdf;
1069  ++t_hess_sdf;
1070  }
1071 
1072  MoFEMFunctionReturn(0);
1073 }
1074 
1075 template <int DIM, typename AssemblyBoundaryEleOp>
1076 OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
1077  OpConstrainBoundaryLhs_dTractionImpl(
1078  const std::string row_field_name, const std::string col_field_name,
1079  boost::shared_ptr<CommonData> common_data_ptr, bool is_axisymmetric)
1080  : AssemblyBoundaryEleOp(row_field_name, col_field_name,
1081  AssemblyBoundaryEleOp::OPROWCOL),
1082  commonDataPtr(common_data_ptr), isAxisymmetric(is_axisymmetric) {
1083  AssemblyBoundaryEleOp::sYmm = false;
1084 }
1085 
1086 template <int DIM, typename AssemblyBoundaryEleOp>
1087 MoFEMErrorCode
1088 OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
1089  iNtegrate(EntitiesFieldData::EntData &row_data,
1090  EntitiesFieldData::EntData &col_data) {
1091  MoFEMFunctionBegin;
1092 
1093  FTensor::Index<'i', DIM> i;
1094  FTensor::Index<'j', DIM> j;
1095  FTensor::Index<'k', DIM> k;
1096 
1097  const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
1098  auto &locMat = AssemblyBoundaryEleOp::locMat;
1099 
1100  auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
1101  auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
1102  auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
1103 
1104  auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
1105  auto t_row_base = row_data.getFTensor1N<3>();
1106  size_t nb_face_functions = row_data.getN().size2() / 3;
1107 
1108  auto m_spatial_coords = get_spatial_coords(
1109  BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
1110  getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
1111  auto m_normals_at_pts = get_normalize_normals(
1112  BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
1113 
1114  auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
1115 
1116  auto ts_time = AssemblyBoundaryEleOp::getTStime();
1117  auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
1118 
1119  // placeholder to pass boundary block id to python
1120  int block_id = 0;
1121 
1122  auto v_sdf =
1123  surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1124  m_spatial_coords, m_normals_at_pts, block_id);
1125 
1126  auto m_grad_sdf =
1127  gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
1128  m_spatial_coords, m_normals_at_pts, block_id);
1129 
1130  auto t_sdf = getFTensor0FromVec(v_sdf);
1131  auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
1132 
1133  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
1134 
1135  double jacobian = 1.;
1136  if (isAxisymmetric) {
1137  jacobian = 2. * M_PI * t_coords(0);
1138  }
1139  const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
1140 
1141  auto tn = -t_traction(i) * t_grad_sdf(i);
1142  auto c = constrain(t_sdf, tn);
1143 
1144  FTensor::Tensor2<double, DIM, DIM> t_cP;
1145  t_cP(i, j) = (c * t_grad_sdf(i)) * t_grad_sdf(j);
1146  FTensor::Tensor2<double, DIM, DIM> t_cQ;
1147  t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
1148 
1149  FTensor::Tensor2<double, DIM, DIM> t_res_dt;
1150  t_res_dt(i, j) = -cn_contact * t_cQ(i, j);
1151 
1152  size_t rr = 0;
1153  for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
1154 
1155  auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
1156  const double row_base = t_row_base(i) * t_normal(i);
1157 
1158  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
1159  for (size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols / DIM; ++cc) {
1160  const double col_base = t_col_base(i) * t_normal(i);
1161  const double beta = alpha * row_base * col_base;
1162 
1163  t_mat(i, j) -= beta * t_res_dt(i, j);
1164 
1165  ++t_col_base;
1166  ++t_mat;
1167  }
1168 
1169  ++t_row_base;
1170  }
1171  for (; rr < nb_face_functions; ++rr)
1172  ++t_row_base;
1173 
1174  ++t_traction;
1175  ++t_coords;
1176  ++t_w;
1177  ++t_normal;
1178  ++t_sdf;
1179  ++t_grad_sdf;
1180  }
1181 
1182  MoFEMFunctionReturn(0);
1183 }
1184 
1185 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
1186 MoFEMErrorCode opFactoryDomainRhs(
1187  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1188  std::string sigma, std::string u, bool is_axisymmetric = false) {
1189  MoFEMFunctionBegin;
1190 
1191  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1192  A>::template LinearForm<I>;
1193  using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
1194  using OpMixDivUCylRhs =
1195  typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
1196 
1197  using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
1198  using OpMixUTimesDivLambdaRhs =
1199  typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
1200  using OpMixUTimesLambdaRhs =
1201  typename B::template OpGradTimesTensor<1, DIM, DIM>;
1202 
1203  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1204  auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
1205  auto div_stress_ptr = boost::make_shared<MatrixDouble>();
1206  auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
1207 
1208  auto jacobian = [is_axisymmetric](const double r, const double,
1209  const double) {
1210  if (is_axisymmetric)
1211  return 2. * M_PI * r;
1212  else
1213  return 1.;
1214  };
1215 
1216  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1217  u, common_data_ptr->contactDispPtr()));
1218  pip.push_back(
1219  new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
1220 
1221  if (!is_axisymmetric) {
1222  pip.push_back(
1223  new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
1224  } else {
1225  pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
1226  sigma, div_stress_ptr));
1227  }
1228 
1229  pip.push_back(new OpCalculateVectorFieldGradient<DIM, DIM>(u, mat_grad_ptr));
1230 
1231  if (!is_axisymmetric) {
1232  pip.push_back(
1233  new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
1234  } else {
1235  pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
1236  jacobian));
1237  }
1238 
1239  pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
1240  pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
1241  pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
1242 
1243  MoFEMFunctionReturn(0);
1244 }
1245 
1246 template <typename OpMixLhs> struct OpMixLhsSide : public OpMixLhs {
1247  using OpMixLhs::OpMixLhs;
1248  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
1249  EntityType col_type,
1250  EntitiesFieldData::EntData &row_data,
1251  EntitiesFieldData::EntData &col_data) {
1252  MoFEMFunctionBegin;
1253  auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
1254  auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
1255  // Only assemble side which correspond to edge entity on boundary
1256  if (side_fe_entity == side_fe_data) {
1257  CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
1258  col_data);
1259  }
1260  MoFEMFunctionReturn(0);
1261  }
1262 };
1263 
1264 template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
1265 MoFEMErrorCode opFactoryBoundaryToDomainLhs(
1266  MoFEM::Interface &m_field,
1267  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1268  std::string fe_domain_name, std::string sigma, std::string u,
1269  std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
1270  bool is_axisymmetric = false) {
1271  MoFEMFunctionBegin;
1272 
1273  using DomainEleOp = typename DomainEle::UserDataOperator;
1274 
1275  auto op_loop_side = new OpLoopSide<DomainEle>(
1276  m_field, fe_domain_name, DIM, Sev::noisy,
1277  boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
1278  pip.push_back(op_loop_side);
1279 
1280  CHKERR AddHOOps<DIM, DIM, DIM>::add(op_loop_side->getOpPtrVector(),
1281  {H1, HDIV}, geom);
1282 
1283  using B = typename FormsIntegrators<DomainEleOp>::template Assembly<
1284  A>::template BiLinearForm<I>;
1285 
1286  using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
1287  using OpMixDivUCylLhs =
1288  typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
1289  using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
1290 
1291  using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
1292  using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
1293  using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
1294 
1295  auto unity = []() { return 1; };
1296  auto jacobian = [is_axisymmetric](const double r, const double,
1297  const double) {
1298  if (is_axisymmetric)
1299  return 2. * M_PI * r;
1300  else
1301  return 1.;
1302  };
1303 
1304  if (!is_axisymmetric) {
1305  op_loop_side->getOpPtrVector().push_back(
1306  new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
1307  } else {
1308  op_loop_side->getOpPtrVector().push_back(
1309  new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
1310  }
1311  op_loop_side->getOpPtrVector().push_back(
1312  new OpLambdaGraULhsSide(sigma, u, unity, jacobian, true));
1313 
1314  op_loop_side->getSideFEPtr()->getRuleHook = rule;
1315  MoFEMFunctionReturn(0);
1316 }
1317 
1318 template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1319 MoFEMErrorCode opFactoryBoundaryLhs(
1320  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1321  std::string sigma, std::string u, bool is_axisymmetric = false) {
1322  MoFEMFunctionBegin;
1323 
1324  using C = ContactIntegrators<BoundaryEleOp>;
1325 
1326  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1327 
1328  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1329  u, common_data_ptr->contactDispPtr()));
1330  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1331  sigma, common_data_ptr->contactTractionPtr()));
1332  pip.push_back(
1333  new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
1334  DIM, GAUSS>(sigma, u, common_data_ptr, is_axisymmetric));
1335  pip.push_back(new typename C::template Assembly<A>::
1336  template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
1337  sigma, sigma, common_data_ptr, is_axisymmetric));
1338 
1339  MoFEMFunctionReturn(0);
1340 }
1341 
1342 template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
1343 MoFEMErrorCode opFactoryBoundaryRhs(
1344  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1345  std::string sigma, std::string u, bool is_axisymmetric = false) {
1346  MoFEMFunctionBegin;
1347 
1348  using C = ContactIntegrators<BoundaryEleOp>;
1349 
1350  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1351 
1352  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1353  u, common_data_ptr->contactDispPtr()));
1354  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1355  sigma, common_data_ptr->contactTractionPtr()));
1356  pip.push_back(
1357  new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
1358  DIM, GAUSS>(sigma, common_data_ptr, is_axisymmetric));
1359 
1360  MoFEMFunctionReturn(0);
1361 }
1362 
1363 template <int DIM, IntegrationType I, typename BoundaryEleOp>
1364 MoFEMErrorCode opFactoryCalculateTraction(
1365  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1366  std::string sigma, bool is_axisymmetric = false) {
1367  MoFEMFunctionBegin;
1368 
1369  using C = ContactIntegrators<BoundaryEleOp>;
1370 
1371  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1372  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1373  sigma, common_data_ptr->contactTractionPtr()));
1374  pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
1375  common_data_ptr, 1. / scale, is_axisymmetric));
1376 
1377  MoFEMFunctionReturn(0);
1378 }
1379 
1380 template <int DIM, IntegrationType I, typename BoundaryEleOp>
1381 MoFEMErrorCode opFactoryCalculateArea(
1382  boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
1383  OpLoopSide<SideEle> *op_loop_side, std::string sigma, std::string u,
1384  bool is_axisymmetric = false,
1385  boost::shared_ptr<Range> contact_range_ptr = nullptr) {
1386  MoFEMFunctionBegin;
1387  using C = ContactIntegrators<BoundaryEleOp>;
1388 
1389  auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1390 
1391  op_loop_side->getOpPtrVector().push_back(
1392  new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
1393  "U", common_data_ptr->contactDispGradPtr()));
1394 
1395  if (contact_range_ptr) {
1396  pip.push_back(new OpCalculateVectorFieldValues<DIM>(
1397  u, common_data_ptr->contactDispPtr()));
1398  pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
1399  sigma, common_data_ptr->contactTractionPtr()));
1400  pip.push_back(op_loop_side);
1401  pip.push_back(new typename C::template OpAssembleTotalContactArea<DIM, I>(
1402  common_data_ptr, is_axisymmetric, contact_range_ptr));
1403  }
1404  MoFEMFunctionReturn(0);
1405 }
1406 
1407 }; // namespace ContactOps
1408 
1409 #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
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
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
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
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:1265
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
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:1343
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:1381
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:1186
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
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: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
DomainEleOp
ContactOps::opFactoryBoundaryLhs
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
Definition: ContactOps.hpp:1319
ContactOps::CommonData::constraintPtr
auto constraintPtr()
Definition: ContactOps.hpp:83
ContactOps::OpMixLhsSide
Definition: ContactOps.hpp:1246
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:1248
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:1364
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:394
Generated by Doxygen 1.8.17 and hosted at University of Glasgow