v0.14.0
EshelbianContact.cpp
Go to the documentation of this file.
1 /**
2  * @file EshelbianContact.cpp
3  * @brief
4  * @date 2023-05-13
5  *
6  * @copyright Copyright (c) 2023
7  *
8  */
9 
10 extern "C" {
11 void tricircumcenter3d_tp(double a[3], double b[3], double c[3],
12  double circumcenter[3], double *xi, double *eta);
13 }
14 
15 namespace EshelbianPlasticity {
16 
17 auto checkSdf(EntityHandle fe_ent, std::map<int, Range> &sdf_map_range) {
18  for (auto &m_sdf : sdf_map_range) {
19  if (m_sdf.second.find(fe_ent) != m_sdf.second.end())
20  return m_sdf.first;
21  }
22  return -1;
23 }
24 
25 template <typename OP_PTR>
26 auto getSdf(OP_PTR op_ptr, MatrixDouble &contact_disp, int block_id,
27  bool eval_hessian) {
28 
29  auto nb_gauss_pts = op_ptr->getGaussPts().size2();
30 
31  auto ts_time = op_ptr->getTStime();
32  auto ts_time_step = op_ptr->getTStimeStep();
33 
34  auto m_spatial_coords = ContactOps::get_spatial_coords(
35  op_ptr->getFTensor1CoordsAtGaussPts(),
36  getFTensor1FromMat<3>(contact_disp), nb_gauss_pts);
37  auto m_normals_at_pts = ContactOps::get_normalize_normals(
38  op_ptr->getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
40  ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
41  block_id);
43  ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
44  block_id);
45 
46 #ifndef NDEBUG
47  if (v_sdf.size() != nb_gauss_pts)
49  "Wrong number of integration pts");
50  if (m_grad_sdf.size2() != nb_gauss_pts)
52  "Wrong number of integration pts");
53  if (m_grad_sdf.size1() != 3)
54  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Should be size of 3");
55 #endif // NDEBUG
56 
57  if (eval_hessian) {
59  ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
60  block_id);
61  return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
62  m_hess_sdf);
63  } else {
64 
65  return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
66  MatrixDouble(6, nb_gauss_pts, 0.));
67  }
68 };
69 
70 /**
71  * @brief Calculate points data on contact surfaces
72  *
73  * @tparam T1
74  * @param unit_ray
75  * @param point
76  * @param elem_point_nodes
77  * @param elem_traction_nodes
78  * @param t_spatial_coords
79  * @return auto
80  */
81 template <typename T1>
83 
84  std::array<double, 3> &unit_ray, std::array<double, 3> &point,
85 
86  std::array<double, 9> &elem_point_nodes,
87  std::array<double, 9> &elem_traction_nodes,
88 
89  FTensor::Tensor1<T1, 3> &t_spatial_coords) {
91 
92  auto t_unit_ray = getFTensor1FromPtr<3>(unit_ray.data());
93  auto t_point = getFTensor1FromPtr<3>(point.data());
94 
95  auto get_normal = [](auto &ele_coords) {
97  Tools::getTriNormal(ele_coords.data(), &t_normal(0));
98  return t_normal;
99  };
100 
101  auto t_normal = get_normal(elem_point_nodes);
102  t_normal(i) /= t_normal.l2();
103 
104  auto sn = t_normal(i) * t_point(i);
105  auto nm = t_normal(i) * t_spatial_coords(i);
106  auto nr = t_normal(i) * t_unit_ray(i);
107 
108  auto gamma = (sn - nm) / nr;
109 
110  FTensor::Tensor1<T1, 3> t_point_current;
111  t_point_current(i) = t_spatial_coords(i) + gamma * t_unit_ray(i);
112 
113  auto get_local_point_shape_functions = [&](auto &&t_elem_coords,
114  auto &t_point) {
115  std::array<T1, 2> loc_coords;
117  Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
118  &t_elem_coords(0, 0), &t_point(0), 1, loc_coords.data()),
119  "get local coords");
120  return FTensor::Tensor1<T1, 3>{N_MBTRI0(loc_coords[0], loc_coords[1]),
121  N_MBTRI1(loc_coords[0], loc_coords[1]),
122  N_MBTRI2(loc_coords[0], loc_coords[1])};
123  };
124 
125  auto eval_position = [&](auto &&t_field, auto &t_shape_fun) {
128  FTensor::Tensor1<T1, 3> t_point_out;
129  t_point_out(i) = t_shape_fun(j) * t_field(j, i);
130  return t_point_out;
131  };
132 
133  auto t_shape_fun = get_local_point_shape_functions(
134  getFTensor2FromPtr<3, 3>(elem_point_nodes.data()), t_point_current);
135  auto t_slave_point_updated = eval_position(
136  getFTensor2FromPtr<3, 3>(elem_point_nodes.data()), t_shape_fun);
137  auto t_traction_updated = eval_position(
138  getFTensor2FromPtr<3, 3>(elem_traction_nodes.data()), t_shape_fun);
139 
140  return std::make_tuple(t_slave_point_updated, t_traction_updated, t_normal);
141 };
142 
143 template <typename T1>
145 
146  ContactTree::FaceData *face_data_ptr,
147  FTensor::Tensor1<T1, 3> &t_spatial_coords
148 
149 ) {
150 
151  return multiPoint(face_data_ptr->unitRay, face_data_ptr->masterPoint,
152  face_data_ptr->masterPointNodes,
153  face_data_ptr->masterTractionNodes, t_spatial_coords);
154 };
155 
156 template <typename T1>
158 
159  ContactTree::FaceData *face_data_ptr,
160  FTensor::Tensor1<T1, 3> &t_spatial_coords
161 
162 ) {
163 
164  return multiPoint(face_data_ptr->unitRay, face_data_ptr->slavePoint,
165  face_data_ptr->slavePointNodes,
166  face_data_ptr->slaveTractionNodes, t_spatial_coords);
167 };
168 
169 /**
170  * Evaluate gap and tractions between master and slave points
171 */
172 template <typename T1>
173 auto multiGetGap(ContactTree::FaceData *face_data_ptr,
174  FTensor::Tensor1<T1, 3> &t_spatial_coords) {
175 
176  auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
177  multiSlavePoint(face_data_ptr, t_spatial_coords);
178  auto [t_master_point_current, t_master_traction_current, t_master_normal] =
179  multiMasterPoint(face_data_ptr, t_spatial_coords);
180 
182 
184  t_normal(i) = t_master_normal(i) - t_slave_normal(i);
185  t_normal.normalize();
186 
187  auto gap = t_normal(i) * (t_slave_point_current(i) - t_spatial_coords(i));
188  auto tn_master = t_master_traction_current(i) * t_normal(i);
189  auto tn_slave = t_slave_traction_current(i) * t_normal(i);
190  auto tn = std::max(-tn_master, tn_slave);
191 
192  return std::make_tuple(gap, tn_master, tn_slave,
193  ContactOps::constrain(gap, tn),
194  t_master_traction_current, t_slave_traction_current);
195 };
196 
198 
199 /** Caluclate rhs term for contact
200  */
201 template <typename T1, typename T2, typename T3>
203 
204  ContactTree::FaceData *face_data_ptr, FTensor::Tensor1<T1, 3> &t_coords,
205  FTensor::Tensor1<T2, 3> &t_spatial_coords,
206  FTensor::Tensor1<T3, 3> &t_master_traction, MultiPointRhsType type,
207  bool debug = false
208 
209 ) {
212 
214  t_u(i) = t_spatial_coords(i) - t_coords(i);
215 
217 
218  switch (type) {
219  case MultiPointRhsType::U: {
220 
221  auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
222  multiSlavePoint(face_data_ptr, t_spatial_coords);
223  auto [t_master_point_current, t_master_traction_current, t_master_normal] =
224  multiMasterPoint(face_data_ptr, t_spatial_coords);
225 
226  // average normal surface
228  t_normal(i) = t_master_normal(i) - t_slave_normal(i);
229  t_normal.normalize();
230 
231  // get projection operators
233  t_P(i, j) = t_normal(i) * t_normal(j);
235  t_Q(i, j) = kronecker_delta(i, j) - t_P(i, j);
236 
237  constexpr double beta = 0.5;
238 
239  auto zeta = 1e-8 * face_data_ptr->eleRadius;
240  auto alpha1 = ContactOps::alpha_contact_const;
242 
243  // this is regularised std::min(d,s)
244  auto f_min_gap = [zeta](auto d, auto s) {
245  return 0.5 * (d + s - std::sqrt((d - s) * (d - s) + zeta));
246  };
247  // this derivative of regularised std::min(d,s)
248  auto f_diff_min_gap = [zeta](auto d, auto s) {
249  return 0.5 * (1 - (d - s) / std::sqrt((d - s) * (d - s) + zeta));
250  };
251 
252  // add penalty on penetration side
253  auto f_barrier = [alpha1, alpha2, f_min_gap, f_diff_min_gap](auto g,
254  auto tn) {
255  auto d = alpha1 * g;
256  auto b1 =
257  0.5 * (tn + f_min_gap(d, tn) + f_diff_min_gap(d, tn) * (d - tn));
258  auto b2 = alpha2 * f_min_gap(g, 0) * g;
259  return b1 - b2;
260  };
261 
262  FTensor::Tensor1<T2, 3> t_gap_vec;
263  t_gap_vec(i) = beta * t_spatial_coords(i) +
264  (1 - beta) * t_slave_point_current(i) - t_spatial_coords(i);
266  t_traction_vec(i) =
267  -beta * t_master_traction(i) + (beta - 1) * t_slave_traction_current(i);
268 
269  auto t_gap = t_normal(i) * t_gap_vec(i);
270  auto t_tn = t_normal(i) * t_traction_vec(i);
271  auto barrier = f_barrier(t_gap, t_tn);
272 
274  // add penalty on penetration side
275  t_traction_bar(i) = t_normal(i) * barrier;
276 
277  t_rhs(i) =
278 
279  t_Q(i, j) * t_master_traction(j) +
280 
281  t_P(i, j) * (t_master_traction(j) - t_traction_bar(j));
282 
283  if (debug) {
284  auto is_nan_or_inf = [](double value) -> bool {
285  return std::isnan(value) || std::isinf(value);
286  };
287 
288  double v = std::complex<double>(t_rhs(i) * t_rhs(i)).real();
289  if (is_nan_or_inf(v)) {
290  MOFEM_LOG_CHANNEL("SELF");
291  MOFEM_LOG("SELF", Sev::error) << "t_rhs " << t_rhs;
292  CHK_MOAB_THROW(MOFEM_DATA_INCONSISTENCY, "rhs is nan or inf");
293  }
294  }
295 
296  } break;
298  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "not implemented");
299  break;
300  }
301 
302  return t_rhs;
303 };
304 
306  const std::string row_field_name,
307  boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
308  boost::shared_ptr<ContactTree> contact_tree_ptr,
309  boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
310  : ContactOps::AssemblyBoundaryEleOp(row_field_name, row_field_name,
311  ContactOps::BoundaryEleOp::OPROW),
312  commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr),
313  sdfMapRangePtr(sdf_map_range_ptr) {
314  CHK_THROW_MESSAGE(PetscOptionsGetScalar(PETSC_NULL, "", "-cn",
315  &ContactOps::cn_contact, PETSC_NULL),
316  "get cn failed");
319  PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_contact_const",
320  &ContactOps::alpha_contact_const, PETSC_NULL),
321  "get alpha contact failed");
323  PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_contact_quadratic",
325  "get alpha contact failed");
327  PetscOptionsGetScalar(PETSC_NULL, "", "-airplane_ray_distance",
328  &ContactOps::airplane_ray_distance, PETSC_NULL),
329  "get alpha contact failed");
330 
331 
332  MOFEM_LOG("EP", Sev::inform) << "cn " << ContactOps::cn_contact;
333  MOFEM_LOG("EP", Sev::inform)
334  << "alpha_contact_const " << ContactOps::alpha_contact_const;
335  MOFEM_LOG("EP", Sev::inform)
336  << "alpha_contact_quadratic " << ContactOps::alpha_contact_quadratic;
337  MOFEM_LOG("EP", Sev::inform)
338  << "airplane_ray_distance " << ContactOps::airplane_ray_distance;
339 }
340 
344 
349 
350  const size_t nb_gauss_pts = getGaussPts().size2();
351 
352 #ifndef NDEBUG
353  if (commonDataPtr->contactDisp.size2() != nb_gauss_pts) {
354  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
355  "Wrong number of integration pts %d != %d",
356  commonDataPtr->contactDisp.size2(), nb_gauss_pts);
357  }
358 #endif // !NDEBUG
359 
360  auto &nf = locF;
361  locF.clear();
362 
363  auto t_w = getFTensor0IntegrationWeight();
364  auto t_coords = getFTensor1CoordsAtGaussPts();
365  auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
366  auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
367 
368  // placeholder to pass boundary block id to python. default SDF is set on
369  // block = -1, one can choose different block by making block "CONTACT_SDF",
370  // then specific SDF can be set to that block.
371  auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
372  getSdf(this, commonDataPtr->contactDisp,
373  checkSdf(getFEEntityHandle(), *sdfMapRangePtr), false);
374 
375  auto t_sdf_v = getFTensor0FromVec(v_sdf);
376  auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
377  auto t_normalize_normal = getFTensor1FromMat<3>(m_normals_at_pts);
378 
379  auto next = [&]() {
380  ++t_w;
381  ++t_coords;
382  ++t_disp;
383  ++t_traction;
384  ++t_normalize_normal;
385  ++t_sdf_v;
386  ++t_grad_sdf_v;
387  };
388 
389  auto face_data_vec_ptr =
390  contactTreePtr->findFaceDataVecPtr(getFEEntityHandle());
391  auto face_gauss_pts_it = face_data_vec_ptr->begin();
392 
393  auto nb_base_functions = data.getN().size2();
394  auto t_base = data.getFTensor0N();
395  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
396 
398  auto face_data_ptr = contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
399  face_data_vec_ptr);
400 
401  auto check_face_contact = [&]() {
402  if (face_data_ptr) {
403  return true;
404  }
405  return false;
406  };
407 
408 #ifdef PYTHON_SDF
409  auto tn = t_traction(i) * t_grad_sdf_v(i);
410  auto c = ContactOps::constrain(t_sdf_v, tn);
411 #else
412  constexpr double c = 0;
413 #endif
414 
415  if (!c && check_face_contact()) {
416  FTensor::Tensor1<double, 3> t_spatial_coords;
417  t_spatial_coords(i) = t_coords(i) + t_disp(i);
418  auto t_rhs_tmp = multiPointRhs(face_data_ptr, t_coords, t_spatial_coords,
419  t_traction, MultiPointRhsType::U, true);
420  t_rhs(i) = t_rhs_tmp(i);
421 
422  } else {
423 
424 #ifdef PYTHON_SDF
425  auto inv_cn = ContactOps::alpha_contact_const;
426 
428  t_cP(i, j) = (c * t_grad_sdf_v(i)) * t_grad_sdf_v(j);
429  t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
430  t_rhs(i) =
431  t_cQ(i, j) * t_traction(j) + (c * inv_cn * t_sdf_v) * t_grad_sdf_v(i);
432 #else
433  t_rhs(i) = t_traction(i);
434 #endif
435  }
436 
437  auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
438  const double alpha = t_w * getMeasure();
439 
440  size_t bb = 0;
441  for (; bb != nbRows / 3; ++bb) {
442  const double beta = alpha * t_base;
443  t_nf(i) += beta * t_rhs(i);
444  ++t_nf;
445  ++t_base;
446  }
447  for (; bb < nb_base_functions; ++bb)
448  ++t_base;
449 
450  next();
451  }
452 
454 }
455 
457  const std::string row_field_name,
458  boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
459  boost::shared_ptr<ContactTree> contact_tree_ptr)
460  : ContactOps::AssemblyBoundaryEleOp(row_field_name, row_field_name,
461  ContactOps::BoundaryEleOp::OPROW),
462  commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr) {}
463 
467 
472 
473  const size_t nb_gauss_pts = getGaussPts().size2();
474 
475 #ifndef NDEBUG
476  if (commonDataPtr->contactDisp.size2() != nb_gauss_pts) {
477  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
478  "Wrong number of integration pts %d != %d",
479  commonDataPtr->contactDisp.size2(), nb_gauss_pts);
480  }
481 #endif // !NDEBUG
482 
483  auto &nf = locF;
484  locF.clear();
485 
486  auto t_w = getFTensor0IntegrationWeight();
487  auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
488  auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
489  auto t_coords = getFTensor1CoordsAtGaussPts();
490  auto t_material_normal = getFTensor1NormalsAtGaussPts();
491 
492  auto next = [&]() {
493  ++t_w;
494  ++t_disp;
495  ++t_traction;
496  ++t_coords;
497  ++t_material_normal;
498  };
499 
500  auto face_data_vec_ptr =
501  contactTreePtr->findFaceDataVecPtr(getFEEntityHandle());
502  auto face_gauss_pts_it = face_data_vec_ptr->begin();
503 
504  auto nb_base_functions = data.getN().size2() / 3;
505  auto t_base = data.getFTensor1N<3>();
506  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
507 
508  auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
509  const double alpha = t_w / 2.;
510 
511  size_t bb = 0;
512  for (; bb != nbRows / 3; ++bb) {
513  const double beta = alpha * t_base(i) * t_material_normal(i);
514  t_nf(i) -= beta * t_disp(i);
515  ++t_nf;
516  ++t_base;
517  }
518  for (; bb < nb_base_functions; ++bb)
519  ++t_base;
520 
521  next();
522  }
523 
525 }
526 
528  const std::string row_field_name, const std::string col_field_name,
529  boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
530  boost::shared_ptr<ContactTree> contact_tree_ptr,
531  boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
532  : ContactOps::AssemblyBoundaryEleOp(row_field_name, col_field_name,
533  ContactOps::BoundaryEleOp::OPROWCOL),
534  commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr),
535  sdfMapRangePtr(sdf_map_range_ptr) {
536 
537  sYmm = false;
538 }
539 
542  EntitiesFieldData::EntData &col_data) {
544 
545  using namespace ContactOps;
546 
550 
551  auto nb_rows = row_data.getIndices().size();
552  auto nb_cols = col_data.getIndices().size();
553 
554  auto &locMat = AssemblyBoundaryEleOp::locMat;
555  locMat.resize(nb_rows, nb_cols, false);
556  locMat.clear();
557 
558  if (nb_cols && nb_rows) {
559 
560  auto nb_gauss_pts = getGaussPts().size2();
561  auto t_w = getFTensor0IntegrationWeight();
562  auto t_coords = getFTensor1CoordsAtGaussPts();
563  auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
564  auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
565 
566  // placeholder to pass boundary block id to python
567  auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
568  getSdf(this, commonDataPtr->contactDisp,
569  checkSdf(getFEEntityHandle(), *sdfMapRangePtr), true);
570 
571  auto t_sdf_v = getFTensor0FromVec(v_sdf);
572  auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
573  auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
574  auto t_normalized_normal = getFTensor1FromMat<3>(m_normals_at_pts);
575 
576  auto next = [&]() {
577  ++t_w;
578  ++t_coords;
579  ++t_disp;
580  ++t_traction;
581  ++t_sdf_v;
582  ++t_grad_sdf_v;
583  ++t_hess_sdf_v;
584  ++t_normalized_normal;
585  };
586 
587  auto face_data_vec_ptr =
588  contactTreePtr->findFaceDataVecPtr(getFEEntityHandle());
589  auto face_gauss_pts_it = face_data_vec_ptr->begin();
590 
591  auto t_row_base = row_data.getFTensor0N();
592  auto nb_face_functions = row_data.getN().size2() / 3;
593  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
594 
595  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
596 
597  auto face_data_ptr = contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
598  face_data_vec_ptr);
599 
600  auto check_face_contact = [&]() {
601  if (face_data_ptr) {
602  return true;
603  }
604  return false;
605  };
606 
608 
609 #ifdef PYTHON_SDF
610  auto tn = t_traction(i) * t_grad_sdf_v(i);
611  auto c = ContactOps::constrain(t_sdf_v, tn);
612 #else
613  constexpr double c = 0;
614 #endif
615 
616  if (!c && check_face_contact()) {
617  FTensor::Tensor1<double, 3> t_spatial_coords;
618  t_spatial_coords(i) = t_coords(i) + t_disp(i);
619  constexpr double eps = std::numeric_limits<float>::epsilon();
620  for (auto ii = 0; ii < 3; ++ii) {
621  FTensor::Tensor1<std::complex<double>, 3> t_spatial_coords_cx{
622  t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2)};
623  t_spatial_coords_cx(ii) += eps * 1i;
624  auto t_rhs_tmp =
625  multiPointRhs(face_data_ptr, t_coords, t_spatial_coords_cx,
626  t_traction, MultiPointRhsType::U);
627  for (int jj = 0; jj != 3; ++jj) {
628  auto v = t_rhs_tmp(jj).imag();
629  t_res_dU(jj, ii) = v / eps;
630  }
631  }
632 
633  } else {
634 
635 #ifdef PYTHON_SDF
636 
637  auto inv_cn = ContactOps::alpha_contact_const;
638 
639  t_res_dU(i, j) =
640 
641  (-c) * (t_hess_sdf_v(i, j) * t_grad_sdf_v(k) * t_traction(k) +
642  t_grad_sdf_v(i) * t_hess_sdf_v(k, j) * t_traction(k))
643 
644  + (c * inv_cn) * (t_sdf_v * t_hess_sdf_v(i, j) +
645 
646  t_grad_sdf_v(j) * t_grad_sdf_v(i));
647 
648 #else
649  t_res_dU(i, j) = 0;
650 #endif
651  }
652 
653  auto alpha = t_w * getMeasure();
654 
655  size_t rr = 0;
656  for (; rr != nb_rows / 3; ++rr) {
657 
658  auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
659  auto t_col_base = col_data.getFTensor0N(gg, 0);
660 
661  for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
662  auto beta = alpha * t_row_base * t_col_base;
663  t_mat(i, j) += beta * t_res_dU(i, j);
664  ++t_col_base;
665  ++t_mat;
666  }
667 
668  ++t_row_base;
669  }
670  for (; rr < nb_face_functions; ++rr)
671  ++t_row_base;
672 
673  next();
674  }
675  }
676 
678 }
679 
681  const std::string row_field_name, const std::string col_field_name,
682  boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
683  boost::shared_ptr<ContactTree> contact_tree_ptr,
684  boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
685  : ContactOps::AssemblyBoundaryEleOp(row_field_name, col_field_name,
686  ContactOps::BoundaryEleOp::OPROWCOL),
687  commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr),
688  sdfMapRangePtr(sdf_map_range_ptr) {
689 
690  sYmm = false;
691 }
692 
695  EntitiesFieldData::EntData &col_data) {
697 
698  using namespace ContactOps;
699 
703 
704  auto nb_rows = row_data.getIndices().size();
705  auto nb_cols = col_data.getIndices().size();
706 
707  auto &locMat = AssemblyBoundaryEleOp::locMat;
708  locMat.resize(nb_rows, nb_cols, false);
709  locMat.clear();
710 
711  if (nb_cols && nb_rows) {
712 
713  const size_t nb_gauss_pts = getGaussPts().size2();
714 
715  auto t_w = getFTensor0IntegrationWeight();
716  auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
717  auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
718  auto t_coords = getFTensor1CoordsAtGaussPts();
719  auto t_material_normal = getFTensor1NormalsAtGaussPts();
720 
721  // placeholder to pass boundary block id to python
722  auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
723  getSdf(this, commonDataPtr->contactDisp,
724  checkSdf(getFEEntityHandle(), *sdfMapRangePtr), false);
725 
726  auto t_sdf_v = getFTensor0FromVec(v_sdf);
727  auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
728 
729  auto next = [&]() {
730  ++t_w;
731  ++t_disp;
732  ++t_traction;
733  ++t_coords;
734  ++t_material_normal;
735  ++t_sdf_v;
736  ++t_grad_sdf_v;
737  };
738 
739  auto face_data_vec_ptr =
740  contactTreePtr->findFaceDataVecPtr(getFEEntityHandle());
741  auto face_gauss_pts_it = face_data_vec_ptr->begin();
742 
743  auto t_row_base = row_data.getFTensor0N();
744  auto nb_face_functions = row_data.getN().size2();
745 
746  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
747 
748  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
749 
750  auto face_data_ptr = contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
751  face_data_vec_ptr);
752 
753  auto check_face_contact = [&]() {
754  if (face_data_ptr) {
755  return true;
756  }
757  return false;
758  };
759 
761 
762 #ifdef PYTHON_SDF
763  auto tn = t_traction(i) * t_grad_sdf_v(i);
764  auto c = ContactOps::constrain(t_sdf_v, tn);
765 #else
766  constexpr double c = 0;
767 #endif
768 
769  if (!c && check_face_contact()) {
770  FTensor::Tensor1<double, 3> t_spatial_coords;
771  t_spatial_coords(i) = t_coords(i) + t_disp(i);
772  constexpr double eps = std::numeric_limits<float>::epsilon();
773  for (auto ii = 0; ii != 3; ++ii) {
774  FTensor::Tensor1<std::complex<double>, 3> t_traction_cx{
775  t_traction(0), t_traction(1), t_traction(2)};
776  t_traction_cx(ii) += eps * 1i;
777  auto t_rhs_tmp =
778  multiPointRhs(face_data_ptr, t_coords, t_spatial_coords,
779  t_traction_cx, MultiPointRhsType::U);
780  for (int jj = 0; jj != 3; ++jj) {
781  auto v = t_rhs_tmp(jj).imag();
782  t_res_dP(jj, ii) = v / eps;
783  }
784  }
785  } else {
786 
787 #ifdef PYTHON_SDF
789  t_cP(i, j) = (c * t_grad_sdf_v(i)) * t_grad_sdf_v(j);
790  t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
791  t_res_dP(i, j) = t_cQ(i, j);
792 #else
793  t_res_dP(i, j) = t_kd(i, j);
794 #endif
795  }
796 
797  const double alpha = t_w / 2.;
798  size_t rr = 0;
799  for (; rr != nb_rows / 3; ++rr) {
800 
801  auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
802  auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
803 
804  for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
805  auto col_base = t_col_base(i) * t_material_normal(i);
806  const double beta = alpha * t_row_base * col_base;
807  t_mat(i, j) += beta * t_res_dP(i, j);
808  ++t_col_base;
809  ++t_mat;
810  }
811 
812  ++t_row_base;
813  }
814  for (; rr < nb_face_functions; ++rr)
815  ++t_row_base;
816 
817  next();
818  }
819  }
820 
822 }
823 
825  const std::string row_field_name, const std::string col_field_name,
826  boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
827  boost::shared_ptr<ContactTree> contact_tree_ptr)
828  : ContactOps::AssemblyBoundaryEleOp(row_field_name, col_field_name,
829  ContactOps::BoundaryEleOp::OPROWCOL),
830  commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr) {
831  sYmm = false;
832 }
833 
836  EntitiesFieldData::EntData &col_data) {
838 
839  using namespace ContactOps;
840 
844 
845  auto nb_rows = row_data.getIndices().size();
846  auto nb_cols = col_data.getIndices().size();
847 
848  auto &locMat = AssemblyBoundaryEleOp::locMat;
849  locMat.resize(nb_rows, nb_cols, false);
850  locMat.clear();
851 
852  if (nb_cols && nb_rows) {
853 
854  const size_t nb_gauss_pts = getGaussPts().size2();
855 
856  auto t_w = getFTensor0IntegrationWeight();
857  auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
858  auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
859  auto t_coords = getFTensor1CoordsAtGaussPts();
860  auto t_material_normal = getFTensor1NormalsAtGaussPts();
861 
862  auto next = [&]() {
863  ++t_w;
864  ++t_disp;
865  ++t_traction;
866  ++t_coords;
867  ++t_material_normal;
868  };
869 
870  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
871 
872  auto face_data_vec_ptr =
873  contactTreePtr->findFaceDataVecPtr(getFEEntityHandle());
874  auto face_gauss_pts_it = face_data_vec_ptr->begin();
875 
876  auto t_row_base = row_data.getFTensor1N<3>();
877  auto nb_face_functions = row_data.getN().size2() / 3;
878  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
879 
880  const auto alpha = t_w / 2.;
881 
882  size_t rr = 0;
883  for (; rr != nb_rows / 3; ++rr) {
884 
885  auto row_base = alpha * (t_row_base(i) * t_material_normal(i));
886 
887  auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
888  auto t_col_base = col_data.getFTensor0N(gg, 0);
889 
890  for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
891  const auto beta = row_base * t_col_base;
892  t_mat(i, j) -= beta * t_kd(i, j);
893  ++t_col_base;
894  ++t_mat;
895  }
896 
897  ++t_row_base;
898  }
899  for (; rr < nb_face_functions; ++rr)
900  ++t_row_base;
901 
902  next();
903  }
904  }
905 
907 }
908 
910  boost::shared_ptr<moab::Core> core_mesh_ptr,
911  int max_order, std::map<int, Range> &&body_map)
912  : Base(m_field, core_mesh_ptr, "contact"), maxOrder(max_order),
913  bodyMap(body_map) {
914 
915  auto ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
916  ref_ele_ptr->hoNodes =
917  PETSC_FALSE; ///< So far only linear geometry is implemented for contact
918 
919  CHK_THROW_MESSAGE(ref_ele_ptr->getOptions(optionsPrefix), "getOptions");
920  CHK_THROW_MESSAGE(ref_ele_ptr->generateReferenceElementMesh(),
921  "Error when generating reference element");
922 
923  MOFEM_LOG("EP", Sev::inform) << "Contact hoNodes " << ref_ele_ptr->hoNodes
924  ? "true"
925  : "false";
926  MOFEM_LOG("EP", Sev::inform)
927  << "Contact maxOrder " << ref_ele_ptr->defMaxLevel;
928 
929  refElementsMap[MBTRI] = ref_ele_ptr;
930 
931  int def_ele_id = -1;
932  CHKERR getPostProcMesh().tag_get_handle("ELE_ID", 1, MB_TYPE_INTEGER, thEleId,
933  MB_TAG_CREAT | MB_TAG_DENSE,
934  &def_ele_id);
935  CHKERR getPostProcMesh().tag_get_handle("BODY_ID", 1, MB_TYPE_INTEGER,
936  thBodyId, MB_TAG_CREAT | MB_TAG_DENSE,
937  &def_ele_id);
938 
939  const auto nb_nodes = (refElementsMap.at(MBTRI)->hoNodes) ? 6 : 3;
940  const auto nb_bases = NBVOLUMETET_L2(maxOrder);
941 
942  // std::vector<int> def_ids(3 * nb_bases, 0);
943  // CHKERR getPostProcMesh().tag_get_handle("IDS", 3 * nb_bases,
944  // MB_TYPE_INTEGER,
945  // thIds, MB_TAG_CREAT | MB_TAG_DENSE,
946  // &*def_ids.begin());
947  // std::vector<double> def_coeffs(3 * nb_bases, 0);
948  // CHKERR getPostProcMesh().tag_get_handle("COEFF", 3 * nb_bases,
949  // MB_TYPE_DOUBLE,
950  // thCoeff, MB_TAG_CREAT |
951  // MB_TAG_DENSE,
952  // &*def_coeffs.begin());
953  // std::vector<double> def_basses(nb_bases, 0);
954  // CHKERR getPostProcMesh().tag_get_handle("BASES", nb_bases, MB_TYPE_DOUBLE,
955  // thBases, MB_TAG_CREAT |
956  // MB_TAG_DENSE,
957  // &*def_basses.begin());
958 
959  std::array<double, 3> def_small_x{0., 0., 0.};
960  CHKERR getPostProcMesh().tag_get_handle("x", 3, MB_TYPE_DOUBLE, thSmallX,
961  MB_TAG_CREAT | MB_TAG_DENSE,
962  def_small_x.data());
963  CHKERR getPostProcMesh().tag_get_handle("X", 3, MB_TYPE_DOUBLE, thLargeX,
964  MB_TAG_CREAT | MB_TAG_DENSE,
965  def_small_x.data());
966 
967  std::array<double, 3> def_tractions{0., 0., 0.};
968  CHKERR getPostProcMesh().tag_get_handle(
969  "TRACTION", 3, MB_TYPE_DOUBLE, thTraction, MB_TAG_CREAT | MB_TAG_DENSE,
970  &*def_tractions.begin());
971 }
972 
975  CHKERR Base::preProcess();
977 }
978 
981 
982  CHKERR Base::postProcess();
983  shadowDataMap.clear();
984 
985  PetscBarrier(nullptr);
986 
987  ParallelComm *pcomm_post_proc_mesh =
988  ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
989  if (pcomm_post_proc_mesh == nullptr)
990  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "PComm not allocated");
991 
992  // MOFEM_LOG("EP", Sev::verbose) << "ContactTree broadcast entities";
993 
994  auto brodacts = [&](auto &brodacts_ents) {
996  Range recived_ents;
997  for (auto rr = 0; rr != mField.get_comm_size(); ++rr) {
998  pcomm_post_proc_mesh->broadcast_entities(
999  rr, rr == mField.get_comm_rank() ? brodacts_ents : recived_ents,
1000  false, true);
1001  }
1003  };
1004 
1005  Range brodacts_ents;
1006  CHKERR getPostProcMesh().get_entities_by_dimension(0, 2, brodacts_ents, true);
1007  CHKERR brodacts(brodacts_ents);
1008 
1009  // MOFEM_LOG("EP", Sev::verbose) << "ContactTree build tree";
1010 
1011  Range ents;
1012  CHKERR getPostProcMesh().get_entities_by_dimension(0, 2, ents, true);
1013  CHKERR buildTree(ents);
1014 
1015  // #ifndef NDEBUG
1016  CHKERR writeFile("debug_tree.h5m");
1017  // #endif // NDEBUG
1018 
1020 }
1021 
1024  // MOFEM_LOG("SELF", Sev::inform) << ents << endl;
1025 
1026  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
1027  new OrientedBoxTreeTool(&getPostProcMesh(), "ROOTSETSURF", true));
1028  CHKERR treeSurfPtr->build(ents, rootSetSurf);
1029 
1030  // CHKERR treeSurfPtr->stats(rootSetSurf, std::cerr);
1031 
1033 }
1034 
1036  boost::shared_ptr<ContactTree> contact_tree_ptr,
1037  boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
1038  boost::shared_ptr<MatrixDouble> u_h1_ptr)
1039  : UOP(NOSPACE, UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
1040  uH1Ptr(u_h1_ptr), commonDataPtr(common_data_ptr) {}
1041 
1045 
1046  auto get_body_id = [&](auto fe_ent) {
1047  for (auto &m : contactTreePtr->bodyMap) {
1048  if (m.second.find(fe_ent) != m.second.end()) {
1049  return m.first;
1050  }
1051  }
1052  return -1;
1053  };
1054 
1055  auto &moab_post_proc_mesh = contactTreePtr->getPostProcMesh();
1056  auto &post_proc_ents = contactTreePtr->getPostProcElements();
1057 
1058  auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1059  auto fe_id = id_from_handle(fe_ent);
1060  auto body_id = get_body_id(fe_ent);
1061  auto &map_gauss_pts = contactTreePtr->getMapGaussPts();
1062 
1063  CHKERR moab_post_proc_mesh.tag_clear_data(contactTreePtr->thEleId,
1064  post_proc_ents, &fe_id);
1065  CHKERR moab_post_proc_mesh.tag_clear_data(contactTreePtr->thBodyId,
1066  post_proc_ents, &body_id);
1067 
1068  auto nb_gauss_pts = getGaussPts().size2();
1069  auto t_u_h1 = getFTensor1FromMat<3>(*uH1Ptr);
1070  auto t_u_l2 = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
1071  auto t_coords = getFTensor1CoordsAtGaussPts();
1072 
1073  MatrixDouble x_h1(nb_gauss_pts, 3);
1074  auto t_x_h1 = getFTensor1FromPtr<3>(&x_h1(0, 0));
1075  MatrixDouble x_l2(nb_gauss_pts, 3);
1076  auto t_x_l2 = getFTensor1FromPtr<3>(&x_l2(0, 0));
1077  MatrixDouble tractions = trans(commonDataPtr->contactTraction);
1078  MatrixDouble coords = getCoordsAtGaussPts();
1079 
1081 
1082  // VectorDouble bases(nb_bases, 0);
1083  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1084  t_x_h1(i) = t_coords(i) + t_u_h1(i);
1085  t_x_l2(i) = t_coords(i) + t_u_l2(i);
1086 
1087  ++t_coords;
1088  ++t_u_h1;
1089  ++t_u_l2;
1090  ++t_x_h1;
1091  ++t_x_l2;
1092  }
1093 
1094  CHKERR moab_post_proc_mesh.set_coords(
1095  &*map_gauss_pts.begin(), map_gauss_pts.size(), &*x_h1.data().begin());
1096  CHKERR moab_post_proc_mesh.tag_set_data(
1097  contactTreePtr->thSmallX, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1098  &*x_h1.data().begin());
1099  CHKERR moab_post_proc_mesh.tag_set_data(
1100  contactTreePtr->thLargeX, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1101  &*coords.data().begin());
1102  CHKERR moab_post_proc_mesh.tag_set_data(
1103  contactTreePtr->thTraction, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1104  &*tractions.data().begin());
1105 
1107 }
1108 
1110  boost::shared_ptr<ContactTree> contact_tree_ptr,
1111  boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
1112  boost::shared_ptr<MatrixDouble> u_h1_ptr, Range r,
1113 
1114  moab::Interface *post_proc_mesh_ptr,
1115  std::vector<EntityHandle> *map_gauss_pts_ptr
1116 
1117  )
1118  : UOP(NOSPACE, UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
1119  commonDataPtr(common_data_ptr), uH1Ptr(u_h1_ptr),
1120  postProcMeshPtr(post_proc_mesh_ptr), mapGaussPtsPtr(map_gauss_pts_ptr),
1121  contactRange(r) {}
1122 
1126 
1127  auto &m_field = getPtrFE()->mField;
1128  auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1129  auto fe_id = id_from_handle(fe_ent);
1130 
1131  if (contactRange.find(fe_ent) == contactRange.end())
1133 
1136 
1137  const auto nb_gauss_pts = getGaussPts().size2();
1138 
1139  auto t_disp_h1 = getFTensor1FromMat<3>(*uH1Ptr);
1140  auto t_coords = getFTensor1CoordsAtGaussPts();
1141  auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
1142 
1143  auto next = [&]() {
1144  ++t_disp_h1;
1145  ++t_traction;
1146  ++t_coords;
1147  };
1148 
1149  auto get_ele_centre = [i](auto t_ele_coords) {
1150  FTensor::Tensor1<double, 3> t_ele_center;
1151  t_ele_center(i) = 0;
1152  for (int nn = 0; nn != 3; nn++) {
1153  t_ele_center(i) += t_ele_coords(i);
1154  ++t_ele_coords;
1155  }
1156  t_ele_center(i) /= 3;
1157  return t_ele_center;
1158  };
1159 
1160  auto get_ele_radius = [i](auto t_ele_center, auto t_ele_coords) {
1162  t_n0(i) = t_ele_center(i) - t_ele_coords(i);
1163  return t_n0.l2();
1164  };
1165 
1166  auto get_face_conn = [this](auto face) {
1167  const EntityHandle *conn;
1168  int num_nodes;
1169  CHK_MOAB_THROW(contactTreePtr->getPostProcMesh().get_connectivity(
1170  face, conn, num_nodes, true),
1171  "get conn");
1172  if (num_nodes != 3) {
1173  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "face is not a triangle");
1174  }
1175  return conn;
1176  };
1177 
1178  auto get_face_coords = [this](auto conn) {
1179  std::array<double, 9> coords;
1180  CHKERR contactTreePtr->getPostProcMesh().get_coords(conn, 3, coords.data());
1181  return coords;
1182  };
1183 
1184  auto get_closet_face = [this](auto *point_ptr, auto r) {
1185  FTensor::Tensor1<double, 3> t_point_out;
1186  std::vector<EntityHandle> faces_out;
1188  contactTreePtr->getTreeSurfPtr()->sphere_intersect_triangles(
1189  point_ptr, r / 8, contactTreePtr->getRootSetSurf(), faces_out),
1190  "get closest faces");
1191  return faces_out;
1192  };
1193 
1194  auto get_faces_out = [this](auto *point_ptr, auto *unit_ray_ptr, auto radius,
1195  auto eps) {
1196  std::vector<double> distances_out;
1197  std::vector<EntityHandle> faces_out;
1199 
1200  contactTreePtr->getTreeSurfPtr()->ray_intersect_triangles(
1201  distances_out, faces_out, contactTreePtr->getRootSetSurf(), eps,
1202  point_ptr, unit_ray_ptr, &radius),
1203 
1204  "get closest faces");
1205  return std::make_pair(faces_out, distances_out);
1206  };
1207 
1208  auto get_normal = [](auto &ele_coords) {
1209  FTensor::Tensor1<double, 3> t_normal;
1210  Tools::getTriNormal(ele_coords.data(), &t_normal(0));
1211  return t_normal;
1212  };
1213 
1214  auto make_map = [&](auto &face_out, auto &face_dist, auto &t_ray_point,
1215  auto &t_unit_ray, auto &t_master_coord) {
1218  std::map<double, EntityHandle> m;
1219  for (auto ii = 0; ii != face_out.size(); ++ii) {
1220  auto face_conn = get_face_conn(face_out[ii]);
1221  auto face_coords = get_face_coords(face_conn);
1222  auto t_face_normal = get_normal(face_coords);
1223  t_face_normal.normalize();
1225  t_x(i) = t_ray_point(i) + t_unit_ray(i) * face_dist[ii];
1227  t_chi(i, j) =
1228  t_x(i) * t_face_normal(j) - t_master_coord(i) * t_unit_ray(j);
1229  if (t_unit_ray(i) * t_face_normal(i) > std::cos(M_PI / 3)) {
1230  auto dot = std::sqrt(t_chi(i, j) * t_chi(i, j));
1231  m[dot] = face_out[ii];
1232  }
1233  }
1234  return m;
1235  };
1236 
1237  auto get_tag_data = [this](auto tag, auto face, auto &vec) {
1239  int tag_length;
1240  CHKERR contactTreePtr->getPostProcMesh().tag_get_length(tag, tag_length);
1241  vec.resize(tag_length);
1242  CHKERR contactTreePtr->getPostProcMesh().tag_get_data(tag, &face, 1,
1243  &*vec.begin());
1245  };
1246 
1247  auto create_tag = [this](const std::string tag_name, const int size) {
1248  double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1249  Tag th;
1250  if (postProcMeshPtr) {
1251  CHKERR postProcMeshPtr->tag_get_handle(
1252  tag_name.c_str(), size, MB_TYPE_DOUBLE, th,
1253  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1254  CHKERR postProcMeshPtr->tag_clear_data(th, &*mapGaussPtsPtr->begin(),
1255  mapGaussPtsPtr->size(), def_VAL);
1256  }
1257  return th;
1258  };
1259 
1260  auto set_float_precision = [](const double x) {
1261  if (std::abs(x) < std::numeric_limits<float>::epsilon())
1262  return 0.;
1263  else
1264  return x;
1265  };
1266 
1267  // scalars
1268  auto save_scal_tag = [&](auto &th, auto v, const int gg) {
1270  if (postProcMeshPtr) {
1271  v = set_float_precision(v);
1272  CHKERR postProcMeshPtr->tag_set_data(th, &(*mapGaussPtsPtr)[gg], 1, &v);
1273  }
1275  };
1276 
1277  // adjacencies
1278  auto get_fe_adjacencies = [this](auto fe_ent) {
1279  Range adj_faces;
1280  CHK_MOAB_THROW(getFaceFE()->mField.get_moab().get_adjacencies(
1281  &fe_ent, 1, 2, false, adj_faces, moab::Interface::UNION),
1282  "get adj");
1283  std::set<int> adj_ids;
1284  for (auto f : adj_faces) {
1285  adj_ids.insert(id_from_handle(f));
1286  }
1287  return adj_ids;
1288  };
1289 
1290  auto get_face_id = [this](auto face) {
1291  int id;
1292  if (contactTreePtr->getPostProcMesh().tag_get_data(
1293  contactTreePtr->thEleId, &face, 1, &id) == MB_SUCCESS) {
1294  return id;
1295  }
1296  return -1;
1297  };
1298 
1299  auto get_body_id = [this](auto face) {
1300  int id;
1301  if (contactTreePtr->getPostProcMesh().tag_get_data(
1302  contactTreePtr->thBodyId, &face, 1, &id) == MB_SUCCESS) {
1303  return id;
1304  }
1305  return -1;
1306  };
1307 
1308  auto get_face_part = [this](auto face) {
1309  ParallelComm *pcomm_post_proc_mesh = ParallelComm::get_pcomm(
1310  &(contactTreePtr->getPostProcMesh()), MYPCOMM_INDEX);
1311  int part;
1312  if (contactTreePtr->getPostProcMesh().tag_get_data(
1313  pcomm_post_proc_mesh->part_tag(), &face, 1, &part) == MB_SUCCESS) {
1314  return part;
1315  }
1316  return -1;
1317  };
1318 
1319  auto check_face = [&](auto face, auto fe_id, auto part) {
1320  auto face_id = get_face_id(face);
1321  auto face_part = get_face_part(face);
1322  if (face_id == fe_id && face_part == part)
1323  return true;
1324  return false;
1325  };
1326 
1327  // vectors
1328  VectorDouble3 v(3);
1329  FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
1330  auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
1332  if (postProcMeshPtr) {
1333  t_v(i) = t_d(i);
1334  for (auto &a : v.data())
1335  a = set_float_precision(a);
1336  CHKERR postProcMeshPtr->tag_set_data(th, &(*mapGaussPtsPtr)[gg], 1,
1337  &*v.data().begin());
1338  }
1340  };
1341 
1342  Tag th_mark = create_tag("contact_mark", 1);
1343  Tag th_mark_slave = create_tag("contact_mark_slave", 1);
1344  Tag th_body_id = create_tag("contact_body_id", 1);
1345  Tag th_gap = create_tag("contact_gap", 1);
1346  Tag th_tn_master = create_tag("contact_tn_master", 1);
1347  Tag th_tn_slave = create_tag("contact_tn_slave", 1);
1348  Tag th_contact_traction = create_tag("contact_traction", 3);
1349  Tag th_contact_traction_master = create_tag("contact_traction_master", 3);
1350  Tag th_contact_traction_slave = create_tag("contact_traction_slave", 3);
1351  Tag th_c = create_tag("contact_c", 1);
1352  Tag th_normal = create_tag("contact_normal", 3);
1353  Tag th_dist = create_tag("contact_dip", 3);
1354 
1355  auto t_ele_centre = get_ele_centre(getFTensor1Coords());
1356  auto ele_radius = get_ele_radius(t_ele_centre, getFTensor1Coords());
1357 
1358  contactTreePtr->shadowDataMap.clear();
1359  auto &shadow_vec = contactTreePtr->shadowDataMap[fe_ent];
1360  shadow_vec.clear();
1361 
1362  auto adj_fe_ids = get_fe_adjacencies(fe_ent);
1363 
1364  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1365 
1366  FTensor::Tensor1<double, 3> t_spatial_coords;
1367  t_spatial_coords(i) = t_coords(i) + t_disp_h1(i);
1368 
1369  if (postProcMeshPtr) {
1370  CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1371  }
1372 
1373  auto faces_close = get_closet_face(&t_spatial_coords(0), ele_radius);
1374  for (auto face_close : faces_close) {
1375  if (check_face(face_close, fe_id, m_field.get_comm_rank())) {
1376 
1377  auto body_id = get_body_id(face_close);
1378 
1379  auto master_face_conn = get_face_conn(face_close);
1380  std::array<double, 9> master_coords;
1381  CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
1382  contactTreePtr->thSmallX, master_face_conn, 3,
1383  master_coords.data());
1384  std::array<double, 9> master_traction;
1385  CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
1386  contactTreePtr->thTraction, master_face_conn, 3,
1387  master_traction.data());
1388  auto t_normal_face_close = get_normal(master_coords);
1389  t_normal_face_close.normalize();
1390 
1391  if (postProcMeshPtr) {
1392  double m = 1;
1393  CHKERR save_scal_tag(th_mark, m, gg);
1394  CHKERR save_scal_tag(th_body_id, static_cast<double>(body_id), gg);
1395  CHKERR save_vec_tag(th_normal, t_normal_face_close, gg);
1396  CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1397  }
1398 
1399  FTensor::Tensor1<double, 3> t_unit_ray;
1400  t_unit_ray(i) = -t_normal_face_close(i);
1401  FTensor::Tensor1<double, 3> t_ray_point;
1402  t_ray_point(i) =
1403  t_spatial_coords(i) -
1404  t_unit_ray(i) * ContactOps::airplane_ray_distance * ele_radius;
1405 
1406  constexpr double eps = 1e-3;
1407  auto [faces_out, faces_dist] =
1408  get_faces_out(&t_ray_point(0), &t_unit_ray(0),
1409  2 * ContactOps::airplane_ray_distance * ele_radius,
1410  eps * ele_radius);
1411 
1412  auto m = make_map(faces_out, faces_dist, t_ray_point, t_unit_ray,
1413  t_spatial_coords);
1414  for (auto m_it = m.begin(); m_it != m.end(); ++m_it) {
1415  auto face = m_it->second;
1416  if (face != face_close) {
1417 
1418  if (
1419 
1420  body_id != get_body_id(face) &&
1421  (adj_fe_ids.find(get_face_id(face)) == adj_fe_ids.end() ||
1422  get_face_part(face) != m_field.get_comm_rank())
1423 
1424  ) {
1425 
1426  shadow_vec.push_back(ContactTree::FaceData());
1427  shadow_vec.back().gaussPtNb = gg;
1428 
1429  auto slave_face_conn = get_face_conn(face);
1430  std::array<double, 9> slave_coords;
1431  CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
1432  contactTreePtr->thSmallX, slave_face_conn, 3,
1433  slave_coords.data());
1434  auto t_normal_face = get_normal(slave_coords);
1435  std::array<double, 9> slave_tractions;
1436  CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
1437  contactTreePtr->thTraction, slave_face_conn, 3,
1438  slave_tractions.data());
1439 
1440  auto t_master_point =
1441  getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data());
1442  auto t_slave_point =
1443  getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data());
1444  auto t_ray_point_data =
1445  getFTensor1FromPtr<3>(shadow_vec.back().rayPoint.data());
1446  auto t_unit_ray_data =
1447  getFTensor1FromPtr<3>(shadow_vec.back().unitRay.data());
1448 
1449  t_slave_point(i) = t_ray_point(i) + m_it->first * t_unit_ray(i);
1450 
1451  auto eval_position = [&](auto &&t_elem_coords, auto &&t_point) {
1452  std::array<double, 2> loc_coords;
1454  Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
1455  &t_elem_coords(0, 0), &t_point(0), 1,
1456  loc_coords.data()),
1457  "get local coords");
1458  FTensor::Tensor1<double, 3> t_shape_fun;
1459  CHK_THROW_MESSAGE(Tools::shapeFunMBTRI<0>(&t_shape_fun(0),
1460  &loc_coords[0],
1461  &loc_coords[1], 1),
1462  "calc shape fun");
1465  FTensor::Tensor1<double, 3> t_point_out;
1466  t_point_out(i) = t_shape_fun(j) * t_elem_coords(j, i);
1467  return t_point_out;
1468  };
1469 
1470  auto t_master_point_updated = eval_position(
1471  getFTensor2FromPtr<3, 3>(master_coords.data()),
1472  getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data()));
1473  t_master_point(i) = t_master_point_updated(i);
1474 
1475  auto t_slave_point_updated = eval_position(
1476  getFTensor2FromPtr<3, 3>(slave_coords.data()),
1477  getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data()));
1478  t_slave_point(i) = t_slave_point_updated(i);
1479 
1480  t_ray_point_data(i) = t_ray_point(i);
1481  t_unit_ray_data(i) = t_unit_ray(i);
1482 
1483  std::copy(master_coords.begin(), master_coords.end(),
1484  shadow_vec.back().masterPointNodes.begin());
1485  std::copy(master_traction.begin(), master_traction.end(),
1486  shadow_vec.back().masterTractionNodes.begin());
1487  std::copy(slave_coords.begin(), slave_coords.end(),
1488  shadow_vec.back().slavePointNodes.begin());
1489  std::copy(slave_tractions.begin(), slave_tractions.end(),
1490  shadow_vec.back().slaveTractionNodes.begin());
1491 
1492  shadow_vec.back().eleRadius = ele_radius;
1493 
1494  // CHKERR get_tag_data(contactTreePtr->thIds, face,
1495  // shadow_vec.back().dofsSlaveIds);
1496  // CHKERR get_tag_data(contactTreePtr->thCoeff, face,
1497  // shadow_vec.back().dofsSlaveCoeff);
1498  // CHKERR get_tag_data(contactTreePtr->thBases, face,
1499  // shadow_vec.back().baseSlaveFuncs);
1500 
1501  if (postProcMeshPtr) {
1502  auto [gap, tn_master, tn_slave, c, t_master_traction,
1503  t_slave_traction] =
1504  multiGetGap(&(shadow_vec.back()), t_spatial_coords);
1505  FTensor::Tensor1<double, 3> t_gap_vec;
1506  t_gap_vec(i) = t_slave_point(i) - t_spatial_coords(i);
1507  CHKERR save_scal_tag(th_gap, gap, gg);
1508  CHKERR save_scal_tag(th_tn_master, tn_master, gg);
1509  CHKERR save_scal_tag(th_tn_slave, tn_slave, gg);
1510  CHKERR save_scal_tag(th_c, c, gg);
1511  double m = 1;
1512  CHKERR save_scal_tag(th_mark_slave, m, gg);
1513  CHKERR save_vec_tag(th_dist, t_gap_vec, gg);
1514  CHKERR save_vec_tag(th_contact_traction_master,
1515  t_master_traction, gg);
1516  CHKERR save_vec_tag(th_contact_traction_slave, t_slave_traction,
1517  gg);
1518  }
1519 
1520  break;
1521  }
1522  }
1523  }
1524  break;
1525  }
1526  }
1527  next();
1528  }
1529 
1531 }
1532 
1533 } // namespace EshelbianPlasticity
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
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
ContactOps::airplane_ray_distance
double airplane_ray_distance
Definition: EshelbianContact.hpp:25
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP::commonDataPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: EshelbianContact.hpp:176
MoFEM::id_from_handle
auto id_from_handle(const EntityHandle h)
Definition: Templates.hpp:1868
g
constexpr double g
Definition: shallow_wave.cpp:63
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
EshelbianPlasticity::ContactTree::buildTree
MoFEMErrorCode buildTree(Range &ents)
Definition: EshelbianContact.cpp:1022
EshelbianPlasticity::OpTreeSearch::UOP
FaceElementForcesAndSourcesCore::UserDataOperator UOP
Definition: EshelbianContact.hpp:215
MoFEM::Types::VectorDouble3
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition: FormsIntegrators.hpp:249
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
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
EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU::OpConstrainBoundaryHDivLhs_dU
OpConstrainBoundaryHDivLhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< ContactTree > contact_tree_ptr)
Definition: EshelbianContact.cpp:824
ContactOps
Definition: EshelbianContact.hpp:10
ContactOps::constrain
double constrain(double sdf, double tn)
constrain function
Definition: ContactOps.hpp:603
EshelbianPlasticity::OpMoveNode::UOP
FaceElementForcesAndSourcesCore::UserDataOperator UOP
Definition: EshelbianContact.hpp:199
EshelbianPlasticity::OpConstrainBoundaryL2Rhs::sdfMapRangePtr
boost::shared_ptr< std::map< int, Range > > sdfMapRangePtr
Definition: EshelbianContact.hpp:130
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
EshelbianPlasticity::OpConstrainBoundaryHDivRhs::commonDataPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: EshelbianContact.hpp:143
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
EshelbianPlasticity::ContactTree::FaceData::masterPoint
std::array< double, 3 > masterPoint
Definition: EshelbianContact.hpp:57
FTensor::Tensor1::l2
T l2() const
Definition: Tensor1_value.hpp:84
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
EshelbianPlasticity::OpTreeSearch::contactTreePtr
boost::shared_ptr< ContactTree > contactTreePtr
Definition: EshelbianContact.hpp:229
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
EshelbianPlasticity::OpTreeSearch::OpTreeSearch
OpTreeSearch(boost::shared_ptr< ContactTree > contact_tree_ptr, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > u_h1_ptr, Range r, moab::Interface *post_proc_mesh_ptr, std::vector< EntityHandle > *map_gauss_pts_ptr)
Definition: EshelbianContact.cpp:1109
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU::commonDataPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: EshelbianContact.hpp:193
EshelbianPlasticity::ContactTree::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: EshelbianContact.cpp:979
EshelbianPlasticity::OpTreeSearch::contactRange
Range contactRange
Definition: EshelbianContact.hpp:236
MoFEM::PostProcBrokenMeshInMoabBase< FaceElementForcesAndSourcesCore >::writeFile
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
Definition: PostProcBrokenMeshInMoabBase.hpp:655
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:27
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP::contactTreePtr
boost::shared_ptr< ContactTree > contactTreePtr
Definition: EshelbianContact.hpp:177
EshelbianPlasticity::OpMoveNode::contactTreePtr
boost::shared_ptr< ContactTree > contactTreePtr
Definition: EshelbianContact.hpp:208
EshelbianPlasticity::checkSdf
auto checkSdf(EntityHandle fe_ent, std::map< int, Range > &sdf_map_range)
Definition: EshelbianContact.cpp:17
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
EshelbianPlasticity::ContactTree::FaceData
Definition: EshelbianContact.hpp:54
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
EshelbianPlasticity::ContactTree::ContactTree
ContactTree(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr, int max_order, std::map< int, Range > &&body_map)
Definition: EshelbianContact.cpp:909
EshelbianPlasticity::ContactTree::thSmallX
Tag thSmallX
Definition: EshelbianContact.hpp:107
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
EshelbianPlasticity::OpMoveNode::OpMoveNode
OpMoveNode(boost::shared_ptr< ContactTree > contact_tree_ptr, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > u_h1_ptr)
Definition: EshelbianContact.cpp:1035
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU::sdfMapRangePtr
boost::shared_ptr< std::map< int, Range > > sdfMapRangePtr
Definition: EshelbianContact.hpp:161
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::PostProcBrokenMeshInMoabBase< FaceElementForcesAndSourcesCore >::refElementsMap
std::map< EntityType, PostProcGenerateRefMeshPtr > refElementsMap
Definition: PostProcBrokenMeshInMoabBase.hpp:171
EshelbianPlasticity::ContactTree::FaceData::unitRay
std::array< double, 3 > unitRay
Definition: EshelbianContact.hpp:59
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
EshelbianPlasticity::ContactTree::thEleId
Tag thEleId
Definition: EshelbianContact.hpp:102
EshelbianPlasticity::OpConstrainBoundaryHDivRhs::contactTreePtr
boost::shared_ptr< ContactTree > contactTreePtr
Definition: EshelbianContact.hpp:144
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
EshelbianPlasticity::ContactTree::rootSetSurf
EntityHandle rootSetSurf
Definition: EshelbianContact.hpp:100
EshelbianPlasticity::ContactTree::shadowDataMap
MapFaceData shadowDataMap
Definition: EshelbianContact.hpp:111
EshelbianPlasticity::OpTreeSearch::postProcMeshPtr
moab::Interface * postProcMeshPtr
Definition: EshelbianContact.hpp:233
a
constexpr double a
Definition: approx_sphere.cpp:30
BoundaryEleOp
EshelbianPlasticity::ContactTree::FaceData::slavePointNodes
std::array< double, 9 > slavePointNodes
Definition: EshelbianContact.hpp:68
EshelbianPlasticity::ContactTree::thBodyId
Tag thBodyId
Definition: EshelbianContact.hpp:106
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
convert.type
type
Definition: convert.py:64
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU::commonDataPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: EshelbianContact.hpp:159
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:197
EshelbianPlasticity::OpConstrainBoundaryHDivRhs::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Definition: EshelbianContact.cpp:465
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: EshelbianContact.cpp:541
EshelbianPlasticity::MultiPointRhsType
MultiPointRhsType
Definition: EshelbianContact.cpp:197
eta
double eta
Definition: free_surface.cpp:170
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
EshelbianPlasticity::multiPoint
auto multiPoint(std::array< double, 3 > &unit_ray, std::array< double, 3 > &point, std::array< double, 9 > &elem_point_nodes, std::array< double, 9 > &elem_traction_nodes, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
Calculate points data on contact surfaces.
Definition: EshelbianContact.cpp:82
EshelbianPlasticity::ContactTree::FaceData::masterPointNodes
std::array< double, 9 > masterPointNodes
Definition: EshelbianContact.hpp:66
EshelbianPlasticity::OpConstrainBoundaryHDivRhs::OpConstrainBoundaryHDivRhs
OpConstrainBoundaryHDivRhs(const std::string row_field_name, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< ContactTree > contact_tree_ptr)
Definition: EshelbianContact.cpp:456
EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: EshelbianContact.cpp:835
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
EshelbianPlasticity::multiMasterPoint
auto multiMasterPoint(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
Definition: EshelbianContact.cpp:144
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
EshelbianPlasticity::ContactTree::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: EshelbianContact.cpp:973
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
EshelbianPlasticity::ContactTree::treeSurfPtr
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
Definition: EshelbianContact.hpp:99
EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU::contactTreePtr
boost::shared_ptr< ContactTree > contactTreePtr
Definition: EshelbianContact.hpp:194
ContactOps::alpha_contact_const
double alpha_contact_const
Definition: EshelbianContact.hpp:21
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
FTensor::Tensor1::normalize
Tensor1< T, Tensor_Dim > normalize()
Definition: Tensor1_value.hpp:77
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP::OpConstrainBoundaryL2Lhs_dP
OpConstrainBoundaryL2Lhs_dP(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< ContactTree > contact_tree_ptr, boost::shared_ptr< std::map< int, Range >> sdf_map_range_ptr=nullptr)
Definition: EshelbianContact.cpp:680
EshelbianPlasticity::OpConstrainBoundaryL2Rhs::OpConstrainBoundaryL2Rhs
OpConstrainBoundaryL2Rhs(const std::string row_field_name, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< ContactTree > contact_tree_ptr, boost::shared_ptr< std::map< int, Range >> sdf_map_range_ptr=nullptr)
Definition: EshelbianContact.cpp:305
N_MBTRI0
#define N_MBTRI0(x, y)
triangle shape function
Definition: fem_tools.h:46
EshelbianPlasticity::ContactTree::Base
PostProcBrokenMeshInMoabBase< FaceElementForcesAndSourcesCore > Base
Definition: EshelbianContact.hpp:37
ContactOps::cn_contact
double cn_contact
Definition: EshelbianContact.hpp:20
EshelbianPlasticity::ContactTree::thTraction
Tag thTraction
Definition: EshelbianContact.hpp:109
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
EshelbianPlasticity::ContactTree::FaceData::slavePoint
std::array< double, 3 > slavePoint
Definition: EshelbianContact.hpp:56
Range
EshelbianPlasticity::ContactTree::FaceData::masterTractionNodes
std::array< double, 9 > masterTractionNodes
Definition: EshelbianContact.hpp:67
face_coords
static const double face_coords[4][9]
Definition: forces_and_sources_h1_continuity_check.cpp:13
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP::sdfMapRangePtr
boost::shared_ptr< std::map< int, Range > > sdfMapRangePtr
Definition: EshelbianContact.hpp:178
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
EshelbianPlasticity::OpTreeSearch::uH1Ptr
boost::shared_ptr< MatrixDouble > uH1Ptr
Definition: EshelbianContact.hpp:231
EshelbianPlasticity::OpMoveNode::commonDataPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: EshelbianContact.hpp:210
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
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU::contactTreePtr
boost::shared_ptr< ContactTree > contactTreePtr
Definition: EshelbianContact.hpp:160
EshelbianPlasticity::OpMoveNode::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: EshelbianContact.cpp:1042
N_MBTRI1
#define N_MBTRI1(x, y)
triangle shape function
Definition: fem_tools.h:47
EshelbianPlasticity::OpTreeSearch::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Definition: EshelbianContact.cpp:1123
EshelbianPlasticity::multiPointRhs
auto multiPointRhs(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, 3 > &t_spatial_coords, FTensor::Tensor1< T3, 3 > &t_master_traction, MultiPointRhsType type, bool debug=false)
Definition: EshelbianContact.cpp:202
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
EshelbianPlasticity::OpMoveNode::uH1Ptr
boost::shared_ptr< MatrixDouble > uH1Ptr
Definition: EshelbianContact.hpp:209
tricircumcenter3d_tp
void tricircumcenter3d_tp(double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
EshelbianPlasticity::ContactTree::FaceData::eleRadius
double eleRadius
Definition: EshelbianContact.hpp:60
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
EshelbianPlasticity::multiGetGap
auto multiGetGap(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
Definition: EshelbianContact.cpp:173
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
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
EshelbianPlasticity::OpConstrainBoundaryL2Rhs::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Definition: EshelbianContact.cpp:342
EshelbianPlasticity::multiSlavePoint
auto multiSlavePoint(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
Definition: EshelbianContact.cpp:157
EshelbianPlasticity::OpTreeSearch::mapGaussPtsPtr
std::vector< EntityHandle > * mapGaussPtsPtr
Definition: EshelbianContact.hpp:234
EshelbianPlasticity::ContactTree::maxOrder
const int maxOrder
Definition: EshelbianContact.hpp:115
MoFEM::PostProcBrokenMeshInMoabBase< FaceElementForcesAndSourcesCore >::getPostProcMesh
auto & getPostProcMesh()
Get postprocessing mesh.
Definition: PostProcBrokenMeshInMoabBase.hpp:641
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
EshelbianPlasticity::OpConstrainBoundaryL2Rhs::contactTreePtr
boost::shared_ptr< ContactTree > contactTreePtr
Definition: EshelbianContact.hpp:129
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Definition: EshelbianContact.cpp:694
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
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::alpha_contact_quadratic
double alpha_contact_quadratic
Definition: EshelbianContact.hpp:22
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
EshelbianPlasticity::OpTreeSearch::commonDataPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: EshelbianContact.hpp:230
EshelbianPlasticity::OpConstrainBoundaryL2Rhs::commonDataPtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: EshelbianContact.hpp:128
EshelbianPlasticity::U
@ U
Definition: EshelbianContact.cpp:197
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
MoFEM::PostProcBrokenMeshInMoabBase< FaceElementForcesAndSourcesCore >::optionsPrefix
std::string optionsPrefix
Prefix for options.
Definition: PostProcBrokenMeshInMoabBase.hpp:176
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU::OpConstrainBoundaryL2Lhs_dU
OpConstrainBoundaryL2Lhs_dU(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< ContactTree > contact_tree_ptr, boost::shared_ptr< std::map< int, Range >> sdf_map_range_ptr=nullptr)
Definition: EshelbianContact.cpp:527
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
N_MBTRI2
#define N_MBTRI2(x, y)
triangle shape function
Definition: fem_tools.h:48
EshelbianPlasticity::getSdf
auto getSdf(OP_PTR op_ptr, MatrixDouble &contact_disp, int block_id, bool eval_hessian)
Definition: EshelbianContact.cpp:26
EshelbianPlasticity::ContactTree::FaceData::slaveTractionNodes
std::array< double, 9 > slaveTractionNodes
Definition: EshelbianContact.hpp:69
EshelbianPlasticity::ContactTree::thLargeX
Tag thLargeX
Definition: EshelbianContact.hpp:108