12 double circumcenter[3],
double *xi,
double *
eta);
18 for (
auto &m_sdf : sdf_map_range) {
19 if (m_sdf.second.find(fe_ent) != m_sdf.second.end())
25 template <
typename OP_PTR>
29 auto nb_gauss_pts = op_ptr->getGaussPts().size2();
31 auto ts_time = op_ptr->getTStime();
32 auto ts_time_step = op_ptr->getTStimeStep();
35 op_ptr->getFTensor1CoordsAtGaussPts(),
36 getFTensor1FromMat<3>(contact_disp), nb_gauss_pts);
38 op_ptr->getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
40 ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
43 ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
47 if (v_sdf.size() != nb_gauss_pts)
49 "Wrong number of integration pts");
50 if (m_grad_sdf.size2() != nb_gauss_pts)
52 "Wrong number of integration pts");
53 if (m_grad_sdf.size1() != 3)
59 ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
61 return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
65 return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
81 template <
typename T1>
84 std::array<double, 3> &unit_ray, std::array<double, 3> &point,
86 std::array<double, 9> &elem_point_nodes,
87 std::array<double, 9> &elem_traction_nodes,
92 auto t_unit_ray = getFTensor1FromPtr<3>(unit_ray.data());
93 auto t_point = getFTensor1FromPtr<3>(point.data());
95 auto get_normal = [](
auto &ele_coords) {
97 Tools::getTriNormal(ele_coords.data(), &t_normal(0));
101 auto t_normal = get_normal(elem_point_nodes);
102 t_normal(
i) /= t_normal.l2();
104 auto sn = t_normal(
i) * t_point(
i);
105 auto nm = t_normal(
i) * t_spatial_coords(
i);
106 auto nr = t_normal(
i) * t_unit_ray(
i);
108 auto gamma = (sn - nm) / nr;
111 t_point_current(
i) = t_spatial_coords(
i) + gamma * t_unit_ray(
i);
113 auto get_local_point_shape_functions = [&](
auto &&t_elem_coords,
115 std::array<T1, 2> loc_coords;
117 Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
118 &t_elem_coords(0, 0), &t_point(0), 1, loc_coords.data()),
121 N_MBTRI1(loc_coords[0], loc_coords[1]),
122 N_MBTRI2(loc_coords[0], loc_coords[1])};
125 auto eval_position = [&](
auto &&t_field,
auto &t_shape_fun) {
129 t_point_out(
i) = t_shape_fun(
j) * t_field(
j,
i);
133 auto t_shape_fun = get_local_point_shape_functions(
134 getFTensor2FromPtr<3, 3>(elem_point_nodes.data()), t_point_current);
135 auto t_slave_point_updated = eval_position(
136 getFTensor2FromPtr<3, 3>(elem_point_nodes.data()), t_shape_fun);
137 auto t_traction_updated = eval_position(
138 getFTensor2FromPtr<3, 3>(elem_traction_nodes.data()), t_shape_fun);
140 return std::make_tuple(t_slave_point_updated, t_traction_updated, t_normal);
143 template <
typename T1>
156 template <
typename T1>
172 template <
typename T1>
176 auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
178 auto [t_master_point_current, t_master_traction_current, t_master_normal] =
184 t_normal(
i) = t_master_normal(
i) - t_slave_normal(
i);
187 auto gap = t_normal(
i) * (t_slave_point_current(
i) - t_spatial_coords(
i));
188 auto tn_master = t_master_traction_current(
i) * t_normal(
i);
189 auto tn_slave = t_slave_traction_current(
i) * t_normal(
i);
190 auto tn = std::max(-tn_master, tn_slave);
192 return std::make_tuple(gap, tn_master, tn_slave,
194 t_master_traction_current, t_slave_traction_current);
201 template <
typename T1,
typename T2,
typename T3>
214 t_u(
i) = t_spatial_coords(
i) - t_coords(
i);
221 auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
223 auto [t_master_point_current, t_master_traction_current, t_master_normal] =
228 t_normal(
i) = t_master_normal(
i) - t_slave_normal(
i);
233 t_P(
i,
j) = t_normal(
i) * t_normal(
j);
237 constexpr
double beta = 0.5;
244 auto f_min_gap = [
zeta](
auto d,
auto s) {
245 return 0.5 * (
d + s - std::sqrt((
d - s) * (
d - s) +
zeta));
248 auto f_diff_min_gap = [
zeta](
auto d,
auto s) {
249 return 0.5 * (1 - (
d - s) / std::sqrt((
d - s) * (
d - s) +
zeta));
253 auto f_barrier = [alpha1, alpha2, f_min_gap, f_diff_min_gap](
auto g,
257 0.5 * (tn + f_min_gap(
d, tn) + f_diff_min_gap(
d, tn) * (
d - tn));
258 auto b2 = alpha2 * f_min_gap(
g, 0) *
g;
263 t_gap_vec(
i) = beta * t_spatial_coords(
i) +
264 (1 - beta) * t_slave_point_current(
i) - t_spatial_coords(
i);
267 -beta * t_master_traction(
i) + (beta - 1) * t_slave_traction_current(
i);
269 auto t_gap = t_normal(
i) * t_gap_vec(
i);
270 auto t_tn = t_normal(
i) * t_traction_vec(
i);
271 auto barrier = f_barrier(t_gap, t_tn);
275 t_traction_bar(
i) = t_normal(
i) * barrier;
279 t_Q(
i,
j) * t_master_traction(
j) +
281 t_P(
i,
j) * (t_master_traction(
j) - t_traction_bar(
j));
284 auto is_nan_or_inf = [](
double value) ->
bool {
285 return std::isnan(value) || std::isinf(value);
288 double v = std::complex<double>(t_rhs(
i) * t_rhs(
i)).real();
289 if (is_nan_or_inf(
v)) {
291 MOFEM_LOG(
"SELF", Sev::error) <<
"t_rhs " << t_rhs;
306 const std::string row_field_name,
307 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
308 boost::shared_ptr<ContactTree> contact_tree_ptr,
309 boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
312 commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr),
313 sdfMapRangePtr(sdf_map_range_ptr) {
321 "get alpha contact failed");
325 "get alpha contact failed");
329 "get alpha contact failed");
350 const size_t nb_gauss_pts = getGaussPts().size2();
355 "Wrong number of integration pts %d != %d",
363 auto t_w = getFTensor0IntegrationWeight();
364 auto t_coords = getFTensor1CoordsAtGaussPts();
365 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
366 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
371 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
376 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
377 auto t_normalize_normal = getFTensor1FromMat<3>(m_normals_at_pts);
384 ++t_normalize_normal;
389 auto face_data_vec_ptr =
391 auto face_gauss_pts_it = face_data_vec_ptr->begin();
393 auto nb_base_functions = data.getN().size2();
394 auto t_base = data.getFTensor0N();
395 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
398 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
401 auto check_face_contact = [&]() {
409 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
412 constexpr
double c = 0;
415 if (!
c && check_face_contact()) {
417 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
418 auto t_rhs_tmp =
multiPointRhs(face_data_ptr, t_coords, t_spatial_coords,
420 t_rhs(
i) = t_rhs_tmp(
i);
428 t_cP(
i,
j) = (
c * t_grad_sdf_v(
i)) * t_grad_sdf_v(
j);
431 t_cQ(
i,
j) * t_traction(
j) + (
c * inv_cn * t_sdf_v) * t_grad_sdf_v(
i);
433 t_rhs(
i) = t_traction(
i);
437 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
438 const double alpha = t_w * getMeasure();
441 for (; bb != nbRows / 3; ++bb) {
442 const double beta = alpha * t_base;
443 t_nf(
i) += beta * t_rhs(
i);
447 for (; bb < nb_base_functions; ++bb)
457 const std::string row_field_name,
458 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
459 boost::shared_ptr<ContactTree> contact_tree_ptr)
462 commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr) {}
473 const size_t nb_gauss_pts = getGaussPts().size2();
478 "Wrong number of integration pts %d != %d",
486 auto t_w = getFTensor0IntegrationWeight();
487 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
488 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
489 auto t_coords = getFTensor1CoordsAtGaussPts();
490 auto t_material_normal = getFTensor1NormalsAtGaussPts();
500 auto face_data_vec_ptr =
502 auto face_gauss_pts_it = face_data_vec_ptr->begin();
504 auto nb_base_functions = data.getN().size2() / 3;
505 auto t_base = data.getFTensor1N<3>();
506 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
508 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
509 const double alpha = t_w / 2.;
512 for (; bb != nbRows / 3; ++bb) {
513 const double beta = alpha * t_base(
i) * t_material_normal(
i);
514 t_nf(
i) -= beta * t_disp(
i);
518 for (; bb < nb_base_functions; ++bb)
528 const std::string row_field_name,
const std::string col_field_name,
529 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
530 boost::shared_ptr<ContactTree> contact_tree_ptr,
531 boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
534 commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr),
535 sdfMapRangePtr(sdf_map_range_ptr) {
551 auto nb_rows = row_data.getIndices().size();
552 auto nb_cols = col_data.getIndices().size();
555 locMat.resize(nb_rows, nb_cols,
false);
558 if (nb_cols && nb_rows) {
560 auto nb_gauss_pts = getGaussPts().size2();
561 auto t_w = getFTensor0IntegrationWeight();
562 auto t_coords = getFTensor1CoordsAtGaussPts();
563 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
564 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
567 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
572 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
573 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
574 auto t_normalized_normal = getFTensor1FromMat<3>(m_normals_at_pts);
584 ++t_normalized_normal;
587 auto face_data_vec_ptr =
589 auto face_gauss_pts_it = face_data_vec_ptr->begin();
591 auto t_row_base = row_data.getFTensor0N();
592 auto nb_face_functions = row_data.getN().size2() / 3;
595 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
597 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
600 auto check_face_contact = [&]() {
610 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
613 constexpr
double c = 0;
616 if (!
c && check_face_contact()) {
618 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
619 constexpr
double eps = std::numeric_limits<float>::epsilon();
620 for (
auto ii = 0; ii < 3; ++ii) {
622 t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2)};
623 t_spatial_coords_cx(ii) +=
eps * 1
i;
627 for (
int jj = 0; jj != 3; ++jj) {
628 auto v = t_rhs_tmp(jj).imag();
629 t_res_dU(jj, ii) =
v /
eps;
641 (-
c) * (t_hess_sdf_v(
i,
j) * t_grad_sdf_v(
k) * t_traction(
k) +
642 t_grad_sdf_v(
i) * t_hess_sdf_v(
k,
j) * t_traction(
k))
644 + (
c * inv_cn) * (t_sdf_v * t_hess_sdf_v(
i,
j) +
646 t_grad_sdf_v(
j) * t_grad_sdf_v(
i));
653 auto alpha = t_w * getMeasure();
656 for (; rr != nb_rows / 3; ++rr) {
658 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
659 auto t_col_base = col_data.getFTensor0N(gg, 0);
661 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
662 auto beta = alpha * t_row_base * t_col_base;
663 t_mat(
i,
j) += beta * t_res_dU(
i,
j);
670 for (; rr < nb_face_functions; ++rr)
681 const std::string row_field_name,
const std::string col_field_name,
682 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
683 boost::shared_ptr<ContactTree> contact_tree_ptr,
684 boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
704 auto nb_rows = row_data.getIndices().size();
705 auto nb_cols = col_data.getIndices().size();
708 locMat.resize(nb_rows, nb_cols,
false);
711 if (nb_cols && nb_rows) {
713 const size_t nb_gauss_pts = getGaussPts().size2();
715 auto t_w = getFTensor0IntegrationWeight();
716 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
717 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
718 auto t_coords = getFTensor1CoordsAtGaussPts();
719 auto t_material_normal = getFTensor1NormalsAtGaussPts();
722 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
727 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
739 auto face_data_vec_ptr =
741 auto face_gauss_pts_it = face_data_vec_ptr->begin();
743 auto t_row_base = row_data.getFTensor0N();
744 auto nb_face_functions = row_data.getN().size2();
748 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
750 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
753 auto check_face_contact = [&]() {
763 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
766 constexpr
double c = 0;
769 if (!
c && check_face_contact()) {
771 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
772 constexpr
double eps = std::numeric_limits<float>::epsilon();
773 for (
auto ii = 0; ii != 3; ++ii) {
775 t_traction(0), t_traction(1), t_traction(2)};
776 t_traction_cx(ii) +=
eps * 1
i;
780 for (
int jj = 0; jj != 3; ++jj) {
781 auto v = t_rhs_tmp(jj).imag();
782 t_res_dP(jj, ii) =
v /
eps;
789 t_cP(
i,
j) = (
c * t_grad_sdf_v(
i)) * t_grad_sdf_v(
j);
791 t_res_dP(
i,
j) = t_cQ(
i,
j);
797 const double alpha = t_w / 2.;
799 for (; rr != nb_rows / 3; ++rr) {
801 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
802 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
804 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
805 auto col_base = t_col_base(
i) * t_material_normal(
i);
806 const double beta = alpha * t_row_base * col_base;
807 t_mat(
i,
j) += beta * t_res_dP(
i,
j);
814 for (; rr < nb_face_functions; ++rr)
825 const std::string row_field_name,
const std::string col_field_name,
826 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
827 boost::shared_ptr<ContactTree> contact_tree_ptr)
845 auto nb_rows = row_data.getIndices().size();
846 auto nb_cols = col_data.getIndices().size();
849 locMat.resize(nb_rows, nb_cols,
false);
852 if (nb_cols && nb_rows) {
854 const size_t nb_gauss_pts = getGaussPts().size2();
856 auto t_w = getFTensor0IntegrationWeight();
857 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
858 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
859 auto t_coords = getFTensor1CoordsAtGaussPts();
860 auto t_material_normal = getFTensor1NormalsAtGaussPts();
872 auto face_data_vec_ptr =
874 auto face_gauss_pts_it = face_data_vec_ptr->begin();
876 auto t_row_base = row_data.getFTensor1N<3>();
877 auto nb_face_functions = row_data.getN().size2() / 3;
878 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
880 const auto alpha = t_w / 2.;
883 for (; rr != nb_rows / 3; ++rr) {
885 auto row_base = alpha * (t_row_base(
i) * t_material_normal(
i));
887 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
888 auto t_col_base = col_data.getFTensor0N(gg, 0);
890 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
891 const auto beta = row_base * t_col_base;
899 for (; rr < nb_face_functions; ++rr)
910 boost::shared_ptr<moab::Core> core_mesh_ptr,
911 int max_order, std::map<int, Range> &&body_map)
912 :
Base(m_field, core_mesh_ptr,
"contact"), maxOrder(max_order),
915 auto ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
916 ref_ele_ptr->hoNodes =
921 "Error when generating reference element");
923 MOFEM_LOG(
"EP", Sev::inform) <<
"Contact hoNodes " << ref_ele_ptr->hoNodes
927 <<
"Contact maxOrder " << ref_ele_ptr->defMaxLevel;
933 MB_TAG_CREAT | MB_TAG_DENSE,
936 thBodyId, MB_TAG_CREAT | MB_TAG_DENSE,
939 const auto nb_nodes = (
refElementsMap.at(MBTRI)->hoNodes) ? 6 : 3;
959 std::array<double, 3> def_small_x{0., 0., 0.};
961 MB_TAG_CREAT | MB_TAG_DENSE,
964 MB_TAG_CREAT | MB_TAG_DENSE,
967 std::array<double, 3> def_tractions{0., 0., 0.};
969 "TRACTION", 3, MB_TYPE_DOUBLE,
thTraction, MB_TAG_CREAT | MB_TAG_DENSE,
970 &*def_tractions.begin());
975 CHKERR Base::preProcess();
982 CHKERR Base::postProcess();
985 PetscBarrier(
nullptr);
987 ParallelComm *pcomm_post_proc_mesh =
989 if (pcomm_post_proc_mesh ==
nullptr)
994 auto brodacts = [&](
auto &brodacts_ents) {
998 pcomm_post_proc_mesh->broadcast_entities(
1005 Range brodacts_ents;
1007 CHKERR brodacts(brodacts_ents);
1026 treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
1036 boost::shared_ptr<ContactTree> contact_tree_ptr,
1037 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
1038 boost::shared_ptr<MatrixDouble> u_h1_ptr)
1039 :
UOP(
NOSPACE,
UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
1040 uH1Ptr(u_h1_ptr), commonDataPtr(common_data_ptr) {}
1046 auto get_body_id = [&](
auto fe_ent) {
1048 if (
m.second.find(fe_ent) !=
m.second.end()) {
1058 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1060 auto body_id = get_body_id(fe_ent);
1064 post_proc_ents, &fe_id);
1066 post_proc_ents, &body_id);
1068 auto nb_gauss_pts = getGaussPts().size2();
1069 auto t_u_h1 = getFTensor1FromMat<3>(*
uH1Ptr);
1070 auto t_u_l2 = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
1071 auto t_coords = getFTensor1CoordsAtGaussPts();
1074 auto t_x_h1 = getFTensor1FromPtr<3>(&x_h1(0, 0));
1076 auto t_x_l2 = getFTensor1FromPtr<3>(&x_l2(0, 0));
1083 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1084 t_x_h1(
i) = t_coords(
i) + t_u_h1(
i);
1085 t_x_l2(
i) = t_coords(
i) + t_u_l2(
i);
1094 CHKERR moab_post_proc_mesh.set_coords(
1095 &*map_gauss_pts.begin(), map_gauss_pts.size(), &*x_h1.data().begin());
1096 CHKERR moab_post_proc_mesh.tag_set_data(
1097 contactTreePtr->thSmallX, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1098 &*x_h1.data().begin());
1099 CHKERR moab_post_proc_mesh.tag_set_data(
1100 contactTreePtr->thLargeX, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1101 &*coords.data().begin());
1102 CHKERR moab_post_proc_mesh.tag_set_data(
1103 contactTreePtr->thTraction, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1104 &*tractions.data().begin());
1110 boost::shared_ptr<ContactTree> contact_tree_ptr,
1111 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
1112 boost::shared_ptr<MatrixDouble> u_h1_ptr,
Range r,
1115 std::vector<EntityHandle> *map_gauss_pts_ptr
1118 :
UOP(
NOSPACE,
UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
1119 commonDataPtr(common_data_ptr), uH1Ptr(u_h1_ptr),
1120 postProcMeshPtr(post_proc_mesh_ptr), mapGaussPtsPtr(map_gauss_pts_ptr),
1127 auto &m_field = getPtrFE()->mField;
1128 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1137 const auto nb_gauss_pts = getGaussPts().size2();
1139 auto t_disp_h1 = getFTensor1FromMat<3>(*
uH1Ptr);
1140 auto t_coords = getFTensor1CoordsAtGaussPts();
1141 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
1149 auto get_ele_centre = [
i](
auto t_ele_coords) {
1151 t_ele_center(
i) = 0;
1152 for (
int nn = 0; nn != 3; nn++) {
1153 t_ele_center(
i) += t_ele_coords(
i);
1156 t_ele_center(
i) /= 3;
1157 return t_ele_center;
1160 auto get_ele_radius = [
i](
auto t_ele_center,
auto t_ele_coords) {
1162 t_n0(
i) = t_ele_center(
i) - t_ele_coords(
i);
1166 auto get_face_conn = [
this](
auto face) {
1170 face, conn, num_nodes,
true),
1172 if (num_nodes != 3) {
1178 auto get_face_coords = [
this](
auto conn) {
1179 std::array<double, 9> coords;
1184 auto get_closet_face = [
this](
auto *point_ptr,
auto r) {
1186 std::vector<EntityHandle> faces_out;
1190 "get closest faces");
1194 auto get_faces_out = [
this](
auto *point_ptr,
auto *unit_ray_ptr,
auto radius,
1196 std::vector<double> distances_out;
1197 std::vector<EntityHandle> faces_out;
1202 point_ptr, unit_ray_ptr, &radius),
1204 "get closest faces");
1205 return std::make_pair(faces_out, distances_out);
1208 auto get_normal = [](
auto &ele_coords) {
1210 Tools::getTriNormal(ele_coords.data(), &t_normal(0));
1214 auto make_map = [&](
auto &face_out,
auto &face_dist,
auto &t_ray_point,
1215 auto &t_unit_ray,
auto &t_master_coord) {
1218 std::map<double, EntityHandle>
m;
1219 for (
auto ii = 0; ii != face_out.size(); ++ii) {
1220 auto face_conn = get_face_conn(face_out[ii]);
1223 t_face_normal.normalize();
1225 t_x(
i) = t_ray_point(
i) + t_unit_ray(
i) * face_dist[ii];
1228 t_x(
i) * t_face_normal(
j) - t_master_coord(
i) * t_unit_ray(
j);
1229 if (t_unit_ray(
i) * t_face_normal(
i) > std::cos(M_PI / 3)) {
1230 auto dot = std::sqrt(t_chi(
i,
j) * t_chi(
i,
j));
1231 m[dot] = face_out[ii];
1237 auto get_tag_data = [
this](
auto tag,
auto face,
auto &vec) {
1241 vec.resize(tag_length);
1247 auto create_tag = [
this](
const std::string tag_name,
const int size) {
1248 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1252 tag_name.c_str(), size, MB_TYPE_DOUBLE,
th,
1253 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1260 auto set_float_precision = [](
const double x) {
1261 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1268 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
1271 v = set_float_precision(
v);
1278 auto get_fe_adjacencies = [
this](
auto fe_ent) {
1281 &fe_ent, 1, 2,
false, adj_faces, moab::Interface::UNION),
1283 std::set<int> adj_ids;
1284 for (
auto f : adj_faces) {
1290 auto get_face_id = [
this](
auto face) {
1299 auto get_body_id = [
this](
auto face) {
1308 auto get_face_part = [
this](
auto face) {
1309 ParallelComm *pcomm_post_proc_mesh = ParallelComm::get_pcomm(
1313 pcomm_post_proc_mesh->part_tag(), &face, 1, &part) == MB_SUCCESS) {
1319 auto check_face = [&](
auto face,
auto fe_id,
auto part) {
1320 auto face_id = get_face_id(face);
1321 auto face_part = get_face_part(face);
1322 if (face_id == fe_id && face_part == part)
1330 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1334 for (
auto &
a :
v.data())
1335 a = set_float_precision(
a);
1337 &*
v.data().begin());
1342 Tag th_mark = create_tag(
"contact_mark", 1);
1343 Tag th_mark_slave = create_tag(
"contact_mark_slave", 1);
1344 Tag th_body_id = create_tag(
"contact_body_id", 1);
1345 Tag th_gap = create_tag(
"contact_gap", 1);
1346 Tag th_tn_master = create_tag(
"contact_tn_master", 1);
1347 Tag th_tn_slave = create_tag(
"contact_tn_slave", 1);
1348 Tag th_contact_traction = create_tag(
"contact_traction", 3);
1349 Tag th_contact_traction_master = create_tag(
"contact_traction_master", 3);
1350 Tag th_contact_traction_slave = create_tag(
"contact_traction_slave", 3);
1351 Tag th_c = create_tag(
"contact_c", 1);
1352 Tag th_normal = create_tag(
"contact_normal", 3);
1353 Tag th_dist = create_tag(
"contact_dip", 3);
1355 auto t_ele_centre = get_ele_centre(getFTensor1Coords());
1356 auto ele_radius = get_ele_radius(t_ele_centre, getFTensor1Coords());
1362 auto adj_fe_ids = get_fe_adjacencies(fe_ent);
1364 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1367 t_spatial_coords(
i) = t_coords(
i) + t_disp_h1(
i);
1370 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1373 auto faces_close = get_closet_face(&t_spatial_coords(0), ele_radius);
1374 for (
auto face_close : faces_close) {
1375 if (check_face(face_close, fe_id, m_field.get_comm_rank())) {
1377 auto body_id = get_body_id(face_close);
1379 auto master_face_conn = get_face_conn(face_close);
1380 std::array<double, 9> master_coords;
1383 master_coords.data());
1384 std::array<double, 9> master_traction;
1387 master_traction.data());
1388 auto t_normal_face_close = get_normal(master_coords);
1389 t_normal_face_close.normalize();
1393 CHKERR save_scal_tag(th_mark,
m, gg);
1394 CHKERR save_scal_tag(th_body_id,
static_cast<double>(body_id), gg);
1395 CHKERR save_vec_tag(th_normal, t_normal_face_close, gg);
1396 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1400 t_unit_ray(
i) = -t_normal_face_close(
i);
1403 t_spatial_coords(
i) -
1406 constexpr
double eps = 1e-3;
1407 auto [faces_out, faces_dist] =
1408 get_faces_out(&t_ray_point(0), &t_unit_ray(0),
1412 auto m = make_map(faces_out, faces_dist, t_ray_point, t_unit_ray,
1414 for (
auto m_it =
m.begin(); m_it !=
m.end(); ++m_it) {
1415 auto face = m_it->second;
1416 if (face != face_close) {
1420 body_id != get_body_id(face) &&
1421 (adj_fe_ids.find(get_face_id(face)) == adj_fe_ids.end() ||
1422 get_face_part(face) != m_field.get_comm_rank())
1427 shadow_vec.back().gaussPtNb = gg;
1429 auto slave_face_conn = get_face_conn(face);
1430 std::array<double, 9> slave_coords;
1433 slave_coords.data());
1434 auto t_normal_face = get_normal(slave_coords);
1435 std::array<double, 9> slave_tractions;
1438 slave_tractions.data());
1440 auto t_master_point =
1441 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data());
1442 auto t_slave_point =
1443 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data());
1444 auto t_ray_point_data =
1445 getFTensor1FromPtr<3>(shadow_vec.back().rayPoint.data());
1446 auto t_unit_ray_data =
1447 getFTensor1FromPtr<3>(shadow_vec.back().unitRay.data());
1449 t_slave_point(
i) = t_ray_point(
i) + m_it->first * t_unit_ray(
i);
1451 auto eval_position = [&](
auto &&t_elem_coords,
auto &&t_point) {
1452 std::array<double, 2> loc_coords;
1454 Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
1455 &t_elem_coords(0, 0), &t_point(0), 1,
1457 "get local coords");
1466 t_point_out(
i) = t_shape_fun(
j) * t_elem_coords(
j,
i);
1470 auto t_master_point_updated = eval_position(
1471 getFTensor2FromPtr<3, 3>(master_coords.data()),
1472 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data()));
1473 t_master_point(
i) = t_master_point_updated(
i);
1475 auto t_slave_point_updated = eval_position(
1476 getFTensor2FromPtr<3, 3>(slave_coords.data()),
1477 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data()));
1478 t_slave_point(
i) = t_slave_point_updated(
i);
1480 t_ray_point_data(
i) = t_ray_point(
i);
1481 t_unit_ray_data(
i) = t_unit_ray(
i);
1483 std::copy(master_coords.begin(), master_coords.end(),
1484 shadow_vec.back().masterPointNodes.begin());
1485 std::copy(master_traction.begin(), master_traction.end(),
1486 shadow_vec.back().masterTractionNodes.begin());
1487 std::copy(slave_coords.begin(), slave_coords.end(),
1488 shadow_vec.back().slavePointNodes.begin());
1489 std::copy(slave_tractions.begin(), slave_tractions.end(),
1490 shadow_vec.back().slaveTractionNodes.begin());
1492 shadow_vec.back().eleRadius = ele_radius;
1502 auto [gap, tn_master, tn_slave,
c, t_master_traction,
1504 multiGetGap(&(shadow_vec.back()), t_spatial_coords);
1506 t_gap_vec(
i) = t_slave_point(
i) - t_spatial_coords(
i);
1507 CHKERR save_scal_tag(th_gap, gap, gg);
1508 CHKERR save_scal_tag(th_tn_master, tn_master, gg);
1509 CHKERR save_scal_tag(th_tn_slave, tn_slave, gg);
1510 CHKERR save_scal_tag(th_c,
c, gg);
1512 CHKERR save_scal_tag(th_mark_slave,
m, gg);
1513 CHKERR save_vec_tag(th_dist, t_gap_vec, gg);
1514 CHKERR save_vec_tag(th_contact_traction_master,
1515 t_master_traction, gg);
1516 CHKERR save_vec_tag(th_contact_traction_slave, t_slave_traction,