v0.15.0
Loading...
Searching...
No Matches
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
10extern "C" {
11void tricircumcenter3d_tp(double a[3], double b[3], double c[3],
12 double circumcenter[3], double *xi, double *eta);
13}
14
15namespace EshelbianPlasticity {
16
17auto 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
25template <typename OP_PTR>
26auto 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 */
85template <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) {
94 FTensor::Index<'i', 3> i;
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) {
130 FTensor::Index<'i', 3> i;
131 FTensor::Index<'j', 3> j;
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
147template <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
160template <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*/
176template <typename T1>
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
185 FTensor::Index<'i', 3> i;
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 */
205template <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) {
214 FTensor::Index<'i', 3> i;
215 FTensor::Index<'j', 3> j;
216
218 t_u(i) = t_spatial_coords(i) - t_coords(i);
219
221
222 switch (type) {
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;
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;
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_NULLPTR, "", "-cn",
319 &ContactOps::cn_contact, PETSC_NULLPTR),
320 "get cn failed");
323 PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_contact_const",
324 &ContactOps::alpha_contact_const, PETSC_NULLPTR),
325 "get alpha contact failed");
327 PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_contact_quadratic",
329 "get alpha contact failed");
331 PetscOptionsGetScalar(PETSC_NULLPTR, "", "-airplane_ray_distance",
332 &ContactOps::airplane_ray_distance, PETSC_NULLPTR),
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
345MoFEMErrorCode
346OpConstrainBoundaryL2Rhs::iNtegrate(EntitiesFieldData::EntData &data) {
348
349 FTensor::Index<'i', 3> i;
350 FTensor::Index<'j', 3> j;
351 FTensor::Index<'k', 3> k;
352 FTensor::Index<'l', 3> l;
353
354 const size_t nb_gauss_pts = getGaussPts().size2();
355
356#ifndef NDEBUG
357 if (commonDataPtr->contactDisp.size2() != nb_gauss_pts) {
358 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
359 "Wrong number of integration pts %ld != %ld",
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
407 if (checkSdf(getFEEntityHandle(), *sdfMapRangePtr) != -1)
408 return false;
409
410 if (face_data_ptr) {
411 return true;
412 }
413 return false;
414 };
415
416#ifdef ENABLE_PYTHON_BINDING
417 double c = 0.;
418 if (ContactOps::sdfPythonWeakPtr.lock()) {
419 auto tn = t_traction(i) * t_grad_sdf_v(i);
420 c = ContactOps::constrain(t_sdf_v, tn);
421 }
422#else
423 constexpr double c = 0;
424#endif
425
426 if (!c && check_face_contact()) {
427 FTensor::Tensor1<double, 3> t_spatial_coords;
428 t_spatial_coords(i) = t_coords(i) + t_disp(i);
429 auto t_rhs_tmp =
430 multiPointRhs(face_data_ptr, t_coords, t_spatial_coords, t_traction,
432 t_rhs(i) = t_rhs_tmp(i);
433
434 } else {
435
436#ifdef ENABLE_PYTHON_BINDING
438
439 if (ContactOps::sdfPythonWeakPtr.lock()) {
441 t_cP(i, j) = (c * t_grad_sdf_v(i)) * t_grad_sdf_v(j);
442 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
443 t_rhs(i) = t_cQ(i, j) * t_traction(j) +
444 (c * inv_cn * t_sdf_v) * t_grad_sdf_v(i);
445 } else {
446 t_rhs(i) = t_traction(i);
447 }
448#else
449 t_rhs(i) = t_traction(i);
450#endif
451 }
452
453 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
454 const double alpha = t_w * getMeasure();
455
456 size_t bb = 0;
457 for (; bb != nbRows / 3; ++bb) {
458 const double beta = alpha * t_base;
459 t_nf(i) -= beta * t_rhs(i);
460 ++t_nf;
461 ++t_base;
462 }
463 for (; bb < nb_base_functions; ++bb)
464 ++t_base;
465
466 next();
467 }
468
470}
471
472template <AssemblyType A>
475 boost::shared_ptr<std::vector<BrokenBaseSideData>>
476 broken_base_side_data,
477 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
478 boost::shared_ptr<ContactTree> contact_tree_ptr)
479 : OP(broken_base_side_data), commonDataPtr(common_data_ptr),
480 contactTreePtr(contact_tree_ptr) {}
481
482template <AssemblyType A>
484 EntitiesFieldData::EntData &data) {
486
491
492 const size_t nb_gauss_pts = OP::getGaussPts().size2();
493
494#ifndef NDEBUG
495 if (commonDataPtr->contactDisp.size2() != nb_gauss_pts) {
496 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
497 "Wrong number of integration pts %ld != %ld",
498 commonDataPtr->contactDisp.size2(), nb_gauss_pts);
499 }
500#endif // !NDEBUG
501
502 auto &nf = OP::locF;
503 OP::locF.clear();
504
505 auto t_w = OP::getFTensor0IntegrationWeight();
506 auto t_material_normal = OP::getFTensor1NormalsAtGaussPts();
507 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
508
509 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
510 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
511
512 auto next = [&]() {
513 ++t_w;
514 ++t_disp;
515 ++t_traction;
516 ++t_coords;
517 ++t_material_normal;
518 };
519
520 auto face_data_vec_ptr =
521 contactTreePtr->findFaceDataVecPtr(OP::getFEEntityHandle());
522 auto face_gauss_pts_it = face_data_vec_ptr->begin();
523
524 auto nb_base_functions = data.getN().size2() / 3;
525 auto t_base = data.getFTensor1N<3>();
526 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
527
528 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
529 const double alpha = t_w / 2.;
530
531 size_t bb = 0;
532 for (; bb != OP::nbRows / 3; ++bb) {
533 const double beta = alpha * t_base(i) * t_material_normal(i);
534 t_nf(i) += beta * t_disp(i);
535 ++t_nf;
536 ++t_base;
537 }
538 for (; bb < nb_base_functions; ++bb)
539 ++t_base;
540
541 next();
542 }
543
545}
546
548 const std::string row_field_name, const std::string col_field_name,
549 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
550 boost::shared_ptr<ContactTree> contact_tree_ptr,
551 boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
552 : ContactOps::AssemblyBoundaryEleOp(row_field_name, col_field_name,
553 ContactOps::BoundaryEleOp::OPROWCOL),
554 commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr),
555 sdfMapRangePtr(sdf_map_range_ptr) {
556
557 sYmm = false;
558}
559
560MoFEMErrorCode
561OpConstrainBoundaryL2Lhs_dU::iNtegrate(EntitiesFieldData::EntData &row_data,
562 EntitiesFieldData::EntData &col_data) {
564
565 using namespace ContactOps;
566
567 FTensor::Index<'i', 3> i;
568 FTensor::Index<'j', 3> j;
569 FTensor::Index<'k', 3> k;
570
571 auto nb_rows = row_data.getIndices().size();
572 auto nb_cols = col_data.getIndices().size();
573
574 auto &locMat = AssemblyBoundaryEleOp::locMat;
575 locMat.resize(nb_rows, nb_cols, false);
576 locMat.clear();
577
578 if (nb_cols && nb_rows) {
579
580 auto nb_gauss_pts = getGaussPts().size2();
581 auto t_w = getFTensor0IntegrationWeight();
582 auto t_coords = getFTensor1CoordsAtGaussPts();
583 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
584 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
585
586 // placeholder to pass boundary block id to python
587 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
588 getSdf(this, commonDataPtr->contactDisp,
589 checkSdf(getFEEntityHandle(), *sdfMapRangePtr), true);
590
591 auto t_sdf_v = getFTensor0FromVec(v_sdf);
592 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
593 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
594 auto t_normalized_normal = getFTensor1FromMat<3>(m_normals_at_pts);
595
596 auto next = [&]() {
597 ++t_w;
598 ++t_coords;
599 ++t_disp;
600 ++t_traction;
601 ++t_sdf_v;
602 ++t_grad_sdf_v;
603 ++t_hess_sdf_v;
604 ++t_normalized_normal;
605 };
606
607 auto face_data_vec_ptr =
608 contactTreePtr->findFaceDataVecPtr(getFEEntityHandle());
609 auto face_gauss_pts_it = face_data_vec_ptr->begin();
610
611 auto t_row_base = row_data.getFTensor0N();
612 auto nb_face_functions = row_data.getN().size2() / 3;
613 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
614
615 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
616
617 auto face_data_ptr = contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
618 face_data_vec_ptr);
619
620 auto check_face_contact = [&]() {
621 if (checkSdf(getFEEntityHandle(), *sdfMapRangePtr) != -1)
622 return false;
623
624 if (face_data_ptr) {
625 return true;
626 }
627 return false;
628 };
629
631
632#ifdef ENABLE_PYTHON_BINDING
633 double c = 0.;
634 if (ContactOps::sdfPythonWeakPtr.lock()) {
635 auto tn = t_traction(i) * t_grad_sdf_v(i);
636 c = ContactOps::constrain(t_sdf_v, tn);
637 }
638#else
639 constexpr double c = 0;
640#endif
641
642 if (!c && check_face_contact()) {
643 FTensor::Tensor1<double, 3> t_spatial_coords;
644 t_spatial_coords(i) = t_coords(i) + t_disp(i);
645 constexpr double eps = std::numeric_limits<float>::epsilon();
646 for (auto ii = 0; ii < 3; ++ii) {
647 FTensor::Tensor1<std::complex<double>, 3> t_spatial_coords_cx{
648 t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2)};
649 t_spatial_coords_cx(ii) += eps * 1i;
650 auto t_rhs_tmp =
651 multiPointRhs(face_data_ptr, t_coords, t_spatial_coords_cx,
652 t_traction, MultiPointRhsType::U);
653 for (int jj = 0; jj != 3; ++jj) {
654 auto v = t_rhs_tmp(jj).imag();
655 t_res_dU(jj, ii) = v / eps;
656 }
657 }
658
659 } else {
660
661#ifdef ENABLE_PYTHON_BINDING
662
663 if (ContactOps::sdfPythonWeakPtr.lock()) {
665 t_res_dU(i, j) =
666
667 (-c) * (t_hess_sdf_v(i, j) * t_grad_sdf_v(k) * t_traction(k) +
668 t_grad_sdf_v(i) * t_hess_sdf_v(k, j) * t_traction(k))
669
670 + (c * inv_cn) * (t_sdf_v * t_hess_sdf_v(i, j) +
671
672 t_grad_sdf_v(j) * t_grad_sdf_v(i));
673 } else {
674 t_res_dU(i, j) = 0;
675 }
676#else
677 t_res_dU(i, j) = 0;
678#endif
679 }
680
681 auto alpha = t_w * getMeasure();
682
683 size_t rr = 0;
684 for (; rr != nb_rows / 3; ++rr) {
685
686 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
687 auto t_col_base = col_data.getFTensor0N(gg, 0);
688
689 for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
690 auto beta = alpha * t_row_base * t_col_base;
691 t_mat(i, j) -= beta * t_res_dU(i, j);
692 ++t_col_base;
693 ++t_mat;
694 }
695
696 ++t_row_base;
697 }
698 for (; rr < nb_face_functions; ++rr)
699 ++t_row_base;
700
701 next();
702 }
703 }
704
706}
707
708template <AssemblyType A>
711 std::string row_field_name,
712 boost::shared_ptr<std::vector<BrokenBaseSideData>>
713 broken_base_side_data,
714 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
715 boost::shared_ptr<ContactTree> contact_tree_ptr,
716 boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
717 : OP(row_field_name, broken_base_side_data, false, false),
718 commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr),
719 sdfMapRangePtr(sdf_map_range_ptr) {
720 OP::sYmm = false;
721}
722
723template <AssemblyType A>
724MoFEMErrorCode
726 EntitiesFieldData::EntData &row_data,
727 EntitiesFieldData::EntData &col_data) {
729
730 using namespace ContactOps;
731
732 FTensor::Index<'i', 3> i;
733 FTensor::Index<'j', 3> j;
734 FTensor::Index<'k', 3> k;
735
736 auto nb_rows = row_data.getIndices().size();
737 auto nb_cols = col_data.getIndices().size();
738
739 auto &locMat = AssemblyBoundaryEleOp::locMat;
740 locMat.resize(nb_rows, nb_cols, false);
741 locMat.clear();
742
743 if (nb_cols && nb_rows) {
744
745 const size_t nb_gauss_pts = OP::getGaussPts().size2();
746
747 auto t_w = OP::getFTensor0IntegrationWeight();
748 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
749 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
750 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
751 auto t_material_normal = OP::getFTensor1NormalsAtGaussPts();
752
753 // placeholder to pass boundary block id to python
754 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
755 getSdf(this, commonDataPtr->contactDisp,
756 checkSdf(OP::getFEEntityHandle(), *sdfMapRangePtr), false);
757
758 auto t_sdf_v = getFTensor0FromVec(v_sdf);
759 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
760
761 auto next = [&]() {
762 ++t_w;
763 ++t_disp;
764 ++t_traction;
765 ++t_coords;
766 ++t_material_normal;
767 ++t_sdf_v;
768 ++t_grad_sdf_v;
769 };
770
771 auto face_data_vec_ptr =
772 contactTreePtr->findFaceDataVecPtr(OP::getFEEntityHandle());
773 auto face_gauss_pts_it = face_data_vec_ptr->begin();
774
775 auto t_row_base = row_data.getFTensor0N();
776 auto nb_face_functions = row_data.getN().size2();
777
778 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
779
780 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
781
782 auto face_data_ptr = contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
783 face_data_vec_ptr);
784
785 auto check_face_contact = [&]() {
786 if (checkSdf(OP::getFEEntityHandle(), *sdfMapRangePtr) != -1)
787 return false;
788
789 if (face_data_ptr) {
790 return true;
791 }
792 return false;
793 };
794
796
797#ifdef ENABLE_PYTHON_BINDING
798 double c = 0.;
799 if (ContactOps::sdfPythonWeakPtr.lock()) {
800 auto tn = t_traction(i) * t_grad_sdf_v(i);
801 c = ContactOps::constrain(t_sdf_v, tn);
802 }
803#else
804 constexpr double c = 0;
805#endif
806
807 if (!c && check_face_contact()) {
808 FTensor::Tensor1<double, 3> t_spatial_coords;
809 t_spatial_coords(i) = t_coords(i) + t_disp(i);
810 constexpr double eps = std::numeric_limits<float>::epsilon();
811 for (auto ii = 0; ii != 3; ++ii) {
813 t_traction(0), t_traction(1), t_traction(2)};
814 t_traction_cx(ii) += eps * 1i;
815 auto t_rhs_tmp =
816 multiPointRhs(face_data_ptr, t_coords, t_spatial_coords,
817 t_traction_cx, MultiPointRhsType::U);
818 for (int jj = 0; jj != 3; ++jj) {
819 auto v = t_rhs_tmp(jj).imag();
820 t_res_dP(jj, ii) = v / eps;
821 }
822 }
823 } else {
824
825#ifdef ENABLE_PYTHON_BINDING
826 if (ContactOps::sdfPythonWeakPtr.lock()) {
828 t_cP(i, j) = (c * t_grad_sdf_v(i)) * t_grad_sdf_v(j);
829 t_cQ(i, j) = kronecker_delta(i, j) - t_cP(i, j);
830 t_res_dP(i, j) = t_cQ(i, j);
831 } else {
832 t_res_dP(i, j) = t_kd(i, j);
833 }
834#else
835 t_res_dP(i, j) = t_kd(i, j);
836#endif
837 }
838
839 const double alpha = t_w / 2.;
840 size_t rr = 0;
841 for (; rr != nb_rows / 3; ++rr) {
842
843 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
844 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
845
846 for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
847 auto col_base = t_col_base(i) * t_material_normal(i);
848 const double beta = alpha * t_row_base * col_base;
849 t_mat(i, j) -= beta * t_res_dP(i, j);
850 ++t_col_base;
851 ++t_mat;
852 }
853
854 ++t_row_base;
855 }
856 for (; rr < nb_face_functions; ++rr)
857 ++t_row_base;
858
859 next();
860 }
861 }
862
864}
865
866template <AssemblyType A>
869 boost::shared_ptr<std::vector<BrokenBaseSideData>>
870 broken_base_side_data,
871 std::string col_field_name,
872 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
873 boost::shared_ptr<ContactTree> contact_tree_ptr)
874 : OP(col_field_name, broken_base_side_data, true, true),
875 commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr) {
876 OP::sYmm = false;
877}
878
879template <AssemblyType A>
880MoFEMErrorCode
882 EntitiesFieldData::EntData &col_data,
883 EntitiesFieldData::EntData &row_data) {
885
886 // Note: col_data and row_data are swapped in this function, we going to
887 // transpose locMat at the end
888
889 using namespace ContactOps;
890
894
895 auto nb_rows = row_data.getIndices().size();
896 auto nb_cols = col_data.getIndices().size();
897
898 auto &locMat = AssemblyBoundaryEleOp::locMat;
899 locMat.resize(nb_rows, nb_cols, false);
900 locMat.clear();
901
902 if (nb_cols && nb_rows) {
903
904 const size_t nb_gauss_pts = OP::getGaussPts().size2();
905
906 auto t_w = OP::getFTensor0IntegrationWeight();
907 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
908 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
909 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
910 auto t_material_normal = OP::getFTensor1NormalsAtGaussPts();
911
912 auto next = [&]() {
913 ++t_w;
914 ++t_disp;
915 ++t_traction;
916 ++t_coords;
917 ++t_material_normal;
918 };
919
920 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
921
922 auto face_data_vec_ptr =
923 contactTreePtr->findFaceDataVecPtr(OP::getFEEntityHandle());
924 auto face_gauss_pts_it = face_data_vec_ptr->begin();
925
926 auto t_row_base = row_data.getFTensor1N<3>();
927 auto nb_face_functions = row_data.getN().size2() / 3;
928 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
929
930 const auto alpha = t_w / 2.;
931
932 size_t rr = 0;
933 for (; rr != nb_rows / 3; ++rr) {
934
935 auto row_base = alpha * (t_row_base(i) * t_material_normal(i));
936
937 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
938 auto t_col_base = col_data.getFTensor0N(gg, 0);
939
940 for (size_t cc = 0; cc != nb_cols / 3; ++cc) {
941 const auto beta = row_base * t_col_base;
942 t_mat(i, j) += beta * t_kd(i, j);
943 ++t_col_base;
944 ++t_mat;
945 }
946
947 ++t_row_base;
948 }
949 for (; rr < nb_face_functions; ++rr)
950 ++t_row_base;
951
952 next();
953 }
954 }
955
956 locMat = trans(locMat);
957
959}
960
962 boost::shared_ptr<moab::Core> core_mesh_ptr,
963 int max_order, std::map<int, Range> &&body_map)
964 : Base(m_field, core_mesh_ptr, "contact"), maxOrder(max_order),
965 bodyMap(body_map) {
966
967 auto ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
968 ref_ele_ptr->hoNodes =
969 PETSC_FALSE; ///< So far only linear geometry is implemented for contact
970
971 CHK_THROW_MESSAGE(ref_ele_ptr->getOptions(optionsPrefix), "getOptions");
972 CHK_THROW_MESSAGE(ref_ele_ptr->generateReferenceElementMesh(),
973 "Error when generating reference element");
974
975 MOFEM_LOG("EP", Sev::inform) << "Contact hoNodes " << ref_ele_ptr->hoNodes
976 ? "true"
977 : "false";
978 MOFEM_LOG("EP", Sev::inform)
979 << "Contact maxOrder " << ref_ele_ptr->defMaxLevel;
980
981 refElementsMap[MBTRI] = ref_ele_ptr;
982
983 int def_ele_id = -1;
984 CHKERR getPostProcMesh().tag_get_handle("ELE_ID", 1, MB_TYPE_INTEGER, thEleId,
985 MB_TAG_CREAT | MB_TAG_DENSE,
986 &def_ele_id);
987 CHKERR getPostProcMesh().tag_get_handle("BODY_ID", 1, MB_TYPE_INTEGER,
988 thBodyId, MB_TAG_CREAT | MB_TAG_DENSE,
989 &def_ele_id);
990
991 // std::vector<int> def_ids(3 * nb_bases, 0);
992 // CHKERR getPostProcMesh().tag_get_handle("IDS", 3 * nb_bases,
993 // MB_TYPE_INTEGER,
994 // thIds, MB_TAG_CREAT | MB_TAG_DENSE,
995 // &*def_ids.begin());
996 // std::vector<double> def_coeffs(3 * nb_bases, 0);
997 // CHKERR getPostProcMesh().tag_get_handle("COEFF", 3 * nb_bases,
998 // MB_TYPE_DOUBLE,
999 // thCoeff, MB_TAG_CREAT |
1000 // MB_TAG_DENSE,
1001 // &*def_coeffs.begin());
1002 // std::vector<double> def_basses(nb_bases, 0);
1003 // CHKERR getPostProcMesh().tag_get_handle("BASES", nb_bases, MB_TYPE_DOUBLE,
1004 // thBases, MB_TAG_CREAT |
1005 // MB_TAG_DENSE,
1006 // &*def_basses.begin());
1007
1008 std::array<double, 3> def_small_x{0., 0., 0.};
1009 CHKERR getPostProcMesh().tag_get_handle("x", 3, MB_TYPE_DOUBLE, thSmallX,
1010 MB_TAG_CREAT | MB_TAG_DENSE,
1011 def_small_x.data());
1012 CHKERR getPostProcMesh().tag_get_handle("X", 3, MB_TYPE_DOUBLE, thLargeX,
1013 MB_TAG_CREAT | MB_TAG_DENSE,
1014 def_small_x.data());
1015
1016 std::array<double, 3> def_tractions{0., 0., 0.};
1017 CHKERR getPostProcMesh().tag_get_handle(
1018 "TRACTION", 3, MB_TYPE_DOUBLE, thTraction, MB_TAG_CREAT | MB_TAG_DENSE,
1019 &*def_tractions.begin());
1020}
1021
1022MoFEMErrorCode ContactTree::preProcess() {
1024 CHKERR Base::preProcess();
1026}
1027
1030
1031 CHKERR Base::postProcess();
1032 shadowDataMap.clear();
1033
1034 PetscBarrier(nullptr);
1035
1036 ParallelComm *pcomm_post_proc_mesh =
1037 ParallelComm::get_pcomm(&getPostProcMesh(), MYPCOMM_INDEX);
1038 if (pcomm_post_proc_mesh == nullptr)
1039 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "PComm not allocated");
1040
1041 // MOFEM_LOG("EP", Sev::verbose) << "ContactTree broadcast entities";
1042
1043 auto brodacts = [&](auto &brodacts_ents) {
1045 Range recived_ents;
1046 for (auto rr = 0; rr != mField.get_comm_size(); ++rr) {
1047 pcomm_post_proc_mesh->broadcast_entities(
1048 rr, rr == mField.get_comm_rank() ? brodacts_ents : recived_ents,
1049 false, true);
1050 }
1052 };
1053
1054 Range brodacts_ents;
1055 CHKERR getPostProcMesh().get_entities_by_dimension(0, 2, brodacts_ents, true);
1056 CHKERR brodacts(brodacts_ents);
1057
1058 // MOFEM_LOG("EP", Sev::verbose) << "ContactTree build tree";
1059
1060 Range ents;
1061 CHKERR getPostProcMesh().get_entities_by_dimension(0, 2, ents, true);
1062 CHKERR buildTree(ents);
1063
1064#ifndef NDEBUG
1065 CHKERR writeFile("debug_tree.h5m");
1066#endif // NDEBUG
1067
1069}
1070
1071MoFEMErrorCode ContactTree::buildTree(Range &ents) {
1073 // MOFEM_LOG("SELF", Sev::inform) << ents << endl;
1074
1075 treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
1076 new OrientedBoxTreeTool(&getPostProcMesh(), "ROOTSETSURF", true));
1077 CHKERR treeSurfPtr->build(ents, rootSetSurf);
1078
1079 // CHKERR treeSurfPtr->stats(rootSetSurf, std::cerr);
1080
1082}
1083
1085 boost::shared_ptr<ContactTree> contact_tree_ptr,
1086 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
1087 boost::shared_ptr<MatrixDouble> u_h1_ptr)
1088 : UOP(NOSPACE, UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
1089 uH1Ptr(u_h1_ptr), commonDataPtr(common_data_ptr) {}
1090
1091MoFEMErrorCode OpMoveNode::doWork(int side, EntityType type,
1092 EntitiesFieldData::EntData &data) {
1094
1095 auto get_body_id = [&](auto fe_ent) {
1096 for (auto &m : contactTreePtr->bodyMap) {
1097 if (m.second.find(fe_ent) != m.second.end()) {
1098 return m.first;
1099 }
1100 }
1101 return -1;
1102 };
1103
1104 auto &moab_post_proc_mesh = contactTreePtr->getPostProcMesh();
1105 auto &post_proc_ents = contactTreePtr->getPostProcElements();
1106
1107 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1108 auto fe_id = id_from_handle(fe_ent);
1109 auto body_id = get_body_id(fe_ent);
1110 auto &map_gauss_pts = contactTreePtr->getMapGaussPts();
1111
1112 CHKERR moab_post_proc_mesh.tag_clear_data(contactTreePtr->thEleId,
1113 post_proc_ents, &fe_id);
1114 CHKERR moab_post_proc_mesh.tag_clear_data(contactTreePtr->thBodyId,
1115 post_proc_ents, &body_id);
1116
1117 auto nb_gauss_pts = getGaussPts().size2();
1118 auto t_u_h1 = getFTensor1FromMat<3>(*uH1Ptr);
1119 auto t_u_l2 = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
1120 auto t_coords = getFTensor1CoordsAtGaussPts();
1121
1122 MatrixDouble x_h1(nb_gauss_pts, 3);
1123 auto t_x_h1 = getFTensor1FromPtr<3>(&x_h1(0, 0));
1124 MatrixDouble x_l2(nb_gauss_pts, 3);
1125 auto t_x_l2 = getFTensor1FromPtr<3>(&x_l2(0, 0));
1126 MatrixDouble tractions = trans(commonDataPtr->contactTraction);
1127 MatrixDouble coords = getCoordsAtGaussPts();
1128
1129 FTensor::Index<'i', 3> i;
1130
1131 // VectorDouble bases(nb_bases, 0);
1132 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1133 t_x_h1(i) = t_coords(i) + t_u_h1(i);
1134 t_x_l2(i) = t_coords(i) + t_u_l2(i);
1135
1136 ++t_coords;
1137 ++t_u_h1;
1138 ++t_u_l2;
1139 ++t_x_h1;
1140 ++t_x_l2;
1141 }
1142
1143 CHKERR moab_post_proc_mesh.set_coords(
1144 &*map_gauss_pts.begin(), map_gauss_pts.size(), &*x_h1.data().begin());
1145 CHKERR moab_post_proc_mesh.tag_set_data(
1146 contactTreePtr->thSmallX, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1147 &*x_h1.data().begin());
1148 CHKERR moab_post_proc_mesh.tag_set_data(
1149 contactTreePtr->thLargeX, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1150 &*coords.data().begin());
1151 CHKERR moab_post_proc_mesh.tag_set_data(
1152 contactTreePtr->thTraction, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1153 &*tractions.data().begin());
1154
1156}
1157
1159 boost::shared_ptr<ContactTree> contact_tree_ptr,
1160 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
1161 boost::shared_ptr<MatrixDouble> u_h1_ptr, Range r,
1162
1163 moab::Interface *post_proc_mesh_ptr,
1164 std::vector<EntityHandle> *map_gauss_pts_ptr
1165
1166 )
1167 : UOP(NOSPACE, UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
1168 commonDataPtr(common_data_ptr), uH1Ptr(u_h1_ptr),
1169 postProcMeshPtr(post_proc_mesh_ptr), mapGaussPtsPtr(map_gauss_pts_ptr),
1170 contactRange(r) {}
1171
1172MoFEMErrorCode OpTreeSearch::doWork(int side, EntityType type,
1173 EntitiesFieldData::EntData &data) {
1175
1176 auto &m_field = getPtrFE()->mField;
1177 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1178 auto fe_id = id_from_handle(fe_ent);
1179
1180 if (contactRange.find(fe_ent) == contactRange.end())
1182
1183 FTensor::Index<'i', 3> i;
1184 FTensor::Index<'j', 3> j;
1185
1186 const auto nb_gauss_pts = getGaussPts().size2();
1187
1188 auto t_disp_h1 = getFTensor1FromMat<3>(*uH1Ptr);
1189 auto t_coords = getFTensor1CoordsAtGaussPts();
1190 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
1191
1192 auto next = [&]() {
1193 ++t_disp_h1;
1194 ++t_traction;
1195 ++t_coords;
1196 };
1197
1198 auto get_ele_centre = [i](auto t_ele_coords) {
1199 FTensor::Tensor1<double, 3> t_ele_center;
1200 t_ele_center(i) = 0;
1201 for (int nn = 0; nn != 3; nn++) {
1202 t_ele_center(i) += t_ele_coords(i);
1203 ++t_ele_coords;
1204 }
1205 t_ele_center(i) /= 3;
1206 return t_ele_center;
1207 };
1208
1209 auto get_ele_radius = [i](auto t_ele_center, auto t_ele_coords) {
1211 t_n0(i) = t_ele_center(i) - t_ele_coords(i);
1212 return t_n0.l2();
1213 };
1214
1215 auto get_face_conn = [this](auto face) {
1216 const EntityHandle *conn;
1217 int num_nodes;
1218 CHK_MOAB_THROW(contactTreePtr->getPostProcMesh().get_connectivity(
1219 face, conn, num_nodes, true),
1220 "get conn");
1221 if (num_nodes != 3) {
1222 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "face is not a triangle");
1223 }
1224 return conn;
1225 };
1226
1227 auto get_face_coords = [this](auto conn) {
1228 std::array<double, 9> coords;
1229 CHKERR contactTreePtr->getPostProcMesh().get_coords(conn, 3, coords.data());
1230 return coords;
1231 };
1232
1233 auto get_closet_face = [this](auto *point_ptr, auto r) {
1234 FTensor::Tensor1<double, 3> t_point_out;
1235 std::vector<EntityHandle> faces_out;
1237 contactTreePtr->getTreeSurfPtr()->sphere_intersect_triangles(
1238 point_ptr, r / 8, contactTreePtr->getRootSetSurf(), faces_out),
1239 "get closest faces");
1240 return faces_out;
1241 };
1242
1243 auto get_faces_out = [this](auto *point_ptr, auto *unit_ray_ptr, auto radius,
1244 auto eps) {
1245 std::vector<double> distances_out;
1246 std::vector<EntityHandle> faces_out;
1248
1249 contactTreePtr->getTreeSurfPtr()->ray_intersect_triangles(
1250 distances_out, faces_out, contactTreePtr->getRootSetSurf(), eps,
1251 point_ptr, unit_ray_ptr, &radius),
1252
1253 "get closest faces");
1254 return std::make_pair(faces_out, distances_out);
1255 };
1256
1257 auto get_normal = [](auto &ele_coords) {
1259 Tools::getTriNormal(ele_coords.data(), &t_normal(0));
1260 return t_normal;
1261 };
1262
1263 auto make_map = [&](auto &face_out, auto &face_dist, auto &t_ray_point,
1264 auto &t_unit_ray, auto &t_master_coord) {
1265 FTensor::Index<'i', 3> i;
1266 FTensor::Index<'j', 3> j;
1267 std::map<double, EntityHandle> m;
1268 for (auto ii = 0; ii != face_out.size(); ++ii) {
1269 auto face_conn = get_face_conn(face_out[ii]);
1270 auto face_coords = get_face_coords(face_conn);
1271 auto t_face_normal = get_normal(face_coords);
1272 t_face_normal.normalize();
1274 t_x(i) = t_ray_point(i) + t_unit_ray(i) * face_dist[ii];
1276 t_chi(i, j) =
1277 t_x(i) * t_face_normal(j) - t_master_coord(i) * t_unit_ray(j);
1278 if (t_unit_ray(i) * t_face_normal(i) > std::cos(M_PI / 3)) {
1279 auto dot = std::sqrt(t_chi(i, j) * t_chi(i, j));
1280 m[dot] = face_out[ii];
1281 }
1282 }
1283 return m;
1284 };
1285
1286 auto get_tag_data = [this](auto tag, auto face, auto &vec) {
1288 int tag_length;
1289 CHKERR contactTreePtr->getPostProcMesh().tag_get_length(tag, tag_length);
1290 vec.resize(tag_length);
1291 CHKERR contactTreePtr->getPostProcMesh().tag_get_data(tag, &face, 1,
1292 &*vec.begin());
1294 };
1295
1296 auto create_tag = [this](const std::string tag_name, const int size) {
1297 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1298 Tag th;
1299 if (postProcMeshPtr) {
1300 CHKERR postProcMeshPtr->tag_get_handle(
1301 tag_name.c_str(), size, MB_TYPE_DOUBLE, th,
1302 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1303 CHKERR postProcMeshPtr->tag_clear_data(th, &*mapGaussPtsPtr->begin(),
1304 mapGaussPtsPtr->size(), def_VAL);
1305 }
1306 return th;
1307 };
1308
1309 auto set_float_precision = [](const double x) {
1310 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1311 return 0.;
1312 else
1313 return x;
1314 };
1315
1316 // scalars
1317 auto save_scal_tag = [&](auto &th, auto v, const int gg) {
1319 if (postProcMeshPtr) {
1320 v = set_float_precision(v);
1321 CHKERR postProcMeshPtr->tag_set_data(th, &(*mapGaussPtsPtr)[gg], 1, &v);
1322 }
1324 };
1325
1326 // adjacencies
1327 auto get_fe_adjacencies = [this](auto fe_ent) {
1328 Range adj_faces;
1329 CHK_MOAB_THROW(getFaceFE()->mField.get_moab().get_adjacencies(
1330 &fe_ent, 1, 2, false, adj_faces, moab::Interface::UNION),
1331 "get adj");
1332 std::set<int> adj_ids;
1333 for (auto f : adj_faces) {
1334 adj_ids.insert(id_from_handle(f));
1335 }
1336 return adj_ids;
1337 };
1338
1339 auto get_face_id = [this](auto face) {
1340 int id;
1341 if (contactTreePtr->getPostProcMesh().tag_get_data(
1342 contactTreePtr->thEleId, &face, 1, &id) == MB_SUCCESS) {
1343 return id;
1344 }
1345 return -1;
1346 };
1347
1348 auto get_body_id = [this](auto face) {
1349 int id;
1350 if (contactTreePtr->getPostProcMesh().tag_get_data(
1351 contactTreePtr->thBodyId, &face, 1, &id) == MB_SUCCESS) {
1352 return id;
1353 }
1354 return -1;
1355 };
1356
1357 auto get_face_part = [this](auto face) {
1358 ParallelComm *pcomm_post_proc_mesh = ParallelComm::get_pcomm(
1359 &(contactTreePtr->getPostProcMesh()), MYPCOMM_INDEX);
1360 int part;
1361 if (contactTreePtr->getPostProcMesh().tag_get_data(
1362 pcomm_post_proc_mesh->part_tag(), &face, 1, &part) == MB_SUCCESS) {
1363 return part;
1364 }
1365 return -1;
1366 };
1367
1368 auto check_face = [&](auto face, auto fe_id, auto part) {
1369 auto face_id = get_face_id(face);
1370 auto face_part = get_face_part(face);
1371 if (face_id == fe_id && face_part == part)
1372 return true;
1373 return false;
1374 };
1375
1376 // vectors
1377 VectorDouble3 v(3);
1378 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_v(&v[0], &v[1], &v[2]);
1379 auto save_vec_tag = [&](auto &th, auto &t_d, const int gg) {
1381 if (postProcMeshPtr) {
1382 t_v(i) = t_d(i);
1383 for (auto &a : v.data())
1384 a = set_float_precision(a);
1385 CHKERR postProcMeshPtr->tag_set_data(th, &(*mapGaussPtsPtr)[gg], 1,
1386 &*v.data().begin());
1387 }
1389 };
1390
1391 Tag th_mark = create_tag("contact_mark", 1);
1392 Tag th_mark_slave = create_tag("contact_mark_slave", 1);
1393 Tag th_body_id = create_tag("contact_body_id", 1);
1394 Tag th_gap = create_tag("contact_gap", 1);
1395 Tag th_tn_master = create_tag("contact_tn_master", 1);
1396 Tag th_tn_slave = create_tag("contact_tn_slave", 1);
1397 Tag th_contact_traction = create_tag("contact_traction", 3);
1398 Tag th_contact_traction_master = create_tag("contact_traction_master", 3);
1399 Tag th_contact_traction_slave = create_tag("contact_traction_slave", 3);
1400 Tag th_c = create_tag("contact_c", 1);
1401 Tag th_normal = create_tag("contact_normal", 3);
1402 Tag th_dist = create_tag("contact_dip", 3);
1403
1404 auto t_ele_centre = get_ele_centre(getFTensor1Coords());
1405 auto ele_radius = get_ele_radius(t_ele_centre, getFTensor1Coords());
1406
1407 contactTreePtr->shadowDataMap.clear();
1408 auto &shadow_vec = contactTreePtr->shadowDataMap[fe_ent];
1409 shadow_vec.clear();
1410
1411 auto adj_fe_ids = get_fe_adjacencies(fe_ent);
1412
1413 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1414
1415 FTensor::Tensor1<double, 3> t_spatial_coords;
1416 t_spatial_coords(i) = t_coords(i) + t_disp_h1(i);
1417
1418 if (postProcMeshPtr) {
1419 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1420 }
1421
1422 auto faces_close = get_closet_face(&t_spatial_coords(0), ele_radius);
1423 for (auto face_close : faces_close) {
1424 if (check_face(face_close, fe_id, m_field.get_comm_rank())) {
1425
1426 auto body_id = get_body_id(face_close);
1427
1428 auto master_face_conn = get_face_conn(face_close);
1429 std::array<double, 9> master_coords;
1430 CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
1431 contactTreePtr->thSmallX, master_face_conn, 3,
1432 master_coords.data());
1433 std::array<double, 9> master_traction;
1434 CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
1435 contactTreePtr->thTraction, master_face_conn, 3,
1436 master_traction.data());
1437 auto t_normal_face_close = get_normal(master_coords);
1438 t_normal_face_close.normalize();
1439
1440 if (postProcMeshPtr) {
1441 double m = 1;
1442 CHKERR save_scal_tag(th_mark, m, gg);
1443 CHKERR save_scal_tag(th_body_id, static_cast<double>(body_id), gg);
1444 CHKERR save_vec_tag(th_normal, t_normal_face_close, gg);
1445 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1446 }
1447
1448 FTensor::Tensor1<double, 3> t_unit_ray;
1449 t_unit_ray(i) = -t_normal_face_close(i);
1450 FTensor::Tensor1<double, 3> t_ray_point;
1451 t_ray_point(i) =
1452 t_spatial_coords(i) -
1453 t_unit_ray(i) * ContactOps::airplane_ray_distance * ele_radius;
1454
1455 constexpr double eps = 1e-3;
1456 auto [faces_out, faces_dist] =
1457 get_faces_out(&t_ray_point(0), &t_unit_ray(0),
1458 2 * ContactOps::airplane_ray_distance * ele_radius,
1459 eps * ele_radius);
1460
1461 auto m = make_map(faces_out, faces_dist, t_ray_point, t_unit_ray,
1462 t_spatial_coords);
1463 for (auto m_it = m.begin(); m_it != m.end(); ++m_it) {
1464 auto face = m_it->second;
1465 if (face != face_close) {
1466
1467 if (
1468
1469 (adj_fe_ids.find(get_face_id(face)) == adj_fe_ids.end() ||
1470 get_face_part(face) != m_field.get_comm_rank())
1471
1472 ) {
1473
1474 shadow_vec.push_back(ContactTree::FaceData());
1475 shadow_vec.back().gaussPtNb = gg;
1476
1477 auto slave_face_conn = get_face_conn(face);
1478 std::array<double, 9> slave_coords;
1479 CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
1480 contactTreePtr->thSmallX, slave_face_conn, 3,
1481 slave_coords.data());
1482 auto t_normal_face = get_normal(slave_coords);
1483 std::array<double, 9> slave_tractions;
1484 CHKERR contactTreePtr->getPostProcMesh().tag_get_data(
1485 contactTreePtr->thTraction, slave_face_conn, 3,
1486 slave_tractions.data());
1487
1488 auto t_master_point =
1489 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data());
1490 auto t_slave_point =
1491 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data());
1492 auto t_ray_point_data =
1493 getFTensor1FromPtr<3>(shadow_vec.back().rayPoint.data());
1494 auto t_unit_ray_data =
1495 getFTensor1FromPtr<3>(shadow_vec.back().unitRay.data());
1496
1497 t_slave_point(i) = t_ray_point(i) + m_it->first * t_unit_ray(i);
1498
1499 auto eval_position = [&](auto &&t_elem_coords, auto &&t_point) {
1500 std::array<double, 2> loc_coords;
1502 Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
1503 &t_elem_coords(0, 0), &t_point(0), 1,
1504 loc_coords.data()),
1505 "get local coords");
1506 FTensor::Tensor1<double, 3> t_shape_fun;
1507 CHK_THROW_MESSAGE(Tools::shapeFunMBTRI<0>(&t_shape_fun(0),
1508 &loc_coords[0],
1509 &loc_coords[1], 1),
1510 "calc shape fun");
1511 FTensor::Index<'i', 3> i;
1512 FTensor::Index<'j', 3> j;
1513 FTensor::Tensor1<double, 3> t_point_out;
1514 t_point_out(i) = t_shape_fun(j) * t_elem_coords(j, i);
1515 return t_point_out;
1516 };
1517
1518 auto t_master_point_updated = eval_position(
1519 getFTensor2FromPtr<3, 3>(master_coords.data()),
1520 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data()));
1521 t_master_point(i) = t_master_point_updated(i);
1522
1523 auto t_slave_point_updated = eval_position(
1524 getFTensor2FromPtr<3, 3>(slave_coords.data()),
1525 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data()));
1526 t_slave_point(i) = t_slave_point_updated(i);
1527
1528 t_ray_point_data(i) = t_ray_point(i);
1529 t_unit_ray_data(i) = t_unit_ray(i);
1530
1531 std::copy(master_coords.begin(), master_coords.end(),
1532 shadow_vec.back().masterPointNodes.begin());
1533 std::copy(master_traction.begin(), master_traction.end(),
1534 shadow_vec.back().masterTractionNodes.begin());
1535 std::copy(slave_coords.begin(), slave_coords.end(),
1536 shadow_vec.back().slavePointNodes.begin());
1537 std::copy(slave_tractions.begin(), slave_tractions.end(),
1538 shadow_vec.back().slaveTractionNodes.begin());
1539
1540 shadow_vec.back().eleRadius = ele_radius;
1541
1542 // CHKERR get_tag_data(contactTreePtr->thIds, face,
1543 // shadow_vec.back().dofsSlaveIds);
1544 // CHKERR get_tag_data(contactTreePtr->thCoeff, face,
1545 // shadow_vec.back().dofsSlaveCoeff);
1546 // CHKERR get_tag_data(contactTreePtr->thBases, face,
1547 // shadow_vec.back().baseSlaveFuncs);
1548
1549 if (postProcMeshPtr) {
1550 auto [gap, tn_master, tn_slave, c, t_master_traction,
1551 t_slave_traction] =
1552 multiGetGap(&(shadow_vec.back()), t_spatial_coords);
1554 t_gap_vec(i) = t_slave_point(i) - t_spatial_coords(i);
1555 CHKERR save_scal_tag(th_gap, gap, gg);
1556 CHKERR save_scal_tag(th_tn_master, tn_master, gg);
1557 CHKERR save_scal_tag(th_tn_slave, tn_slave, gg);
1558 CHKERR save_scal_tag(th_c, c, gg);
1559 double m = 1;
1560 CHKERR save_scal_tag(th_mark_slave, m, gg);
1561 CHKERR save_vec_tag(th_dist, t_gap_vec, gg);
1562 CHKERR save_vec_tag(th_contact_traction_master,
1563 t_master_traction, gg);
1564 CHKERR save_vec_tag(th_contact_traction_slave, t_slave_traction,
1565 gg);
1566 }
1567
1568 break;
1569 }
1570 }
1571 }
1572 break;
1573 }
1574 }
1575 next();
1576 }
1577
1579}
1580
1581} // namespace EshelbianPlasticity
void tricircumcenter3d_tp(double a[3], double b[3], double c[3], double circumcenter[3], double *xi, double *eta)
#define FTENSOR_INDEX(DIM, I)
constexpr double a
static const double eps
constexpr int SPACE_DIM
Kronecker Delta class.
Tensor1< T, Tensor_Dim > normalize()
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ NOSPACE
Definition definitions.h:83
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
static const bool debug
#define N_MBTRI1(x, y)
triangle shape function
Definition fem_tools.h:47
#define N_MBTRI0(x, y)
triangle shape function
Definition fem_tools.h:46
#define N_MBTRI2(x, y)
triangle shape function
Definition fem_tools.h:48
static const double face_coords[4][9]
constexpr auto t_kd
double eta
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
auto get_normalize_normals(FTensor::Tensor1< T1, DIM1 > &&t_normal_at_pts, size_t nb_gauss_pts)
double cn_contact
Definition contact.cpp:99
double alpha_contact_quadratic
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)
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)
double alpha_contact_const
double airplane_ray_distance
auto get_spatial_coords(FTensor::Tensor1< T1, DIM1 > &&t_coords, FTensor::Tensor1< T2, DIM2 > &&t_disp, size_t nb_gauss_pts)
double constrain(double sdf, double tn)
constrain 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)
auto checkSdf(EntityHandle fe_ent, std::map< int, Range > &sdf_map_range)
auto multiMasterPoint(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
auto multiGetGap(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
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)
auto multiSlavePoint(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
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.
auto getSdf(OP_PTR op_ptr, MatrixDouble &contact_disp, int block_id, bool eval_hessian)
double zeta
Viscous hardening.
Definition plastic.cpp:130
constexpr double g
FTensor::Index< 'm', 3 > m
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
PostProcBrokenMeshInMoabBase< FaceElementForcesAndSourcesCore > Base
MoFEMErrorCode buildTree(Range &ents)
ContactTree(MoFEM::Interface &m_field, boost::shared_ptr< moab::Core > core_mesh_ptr, int max_order, std::map< int, Range > &&body_map)
boost::shared_ptr< std::map< int, Range > > sdfMapRangePtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
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)
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
boost::shared_ptr< std::map< int, Range > > sdfMapRangePtr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
boost::shared_ptr< ContactTree > contactTreePtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
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)
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
boost::shared_ptr< ContactTree > contactTreePtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
FaceElementForcesAndSourcesCore::UserDataOperator UOP
boost::shared_ptr< MatrixDouble > uH1Ptr
OpMoveNode(boost::shared_ptr< ContactTree > contact_tree_ptr, boost::shared_ptr< ContactOps::CommonData > common_data_ptr, boost::shared_ptr< MatrixDouble > u_h1_ptr)
boost::shared_ptr< ContactTree > contactTreePtr
std::vector< EntityHandle > * mapGaussPtsPtr
boost::shared_ptr< MatrixDouble > uH1Ptr
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)
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
FaceElementForcesAndSourcesCore::UserDataOperator UOP
Deprecated interface functions.