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