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