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();
39 op_ptr->getFTensor1CoordsAtGaussPts(),
40 getFTensor1FromMat<3>(contact_disp), nb_gauss_pts);
42 op_ptr->getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
44 ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
47 ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
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)
63 ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
65 return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
69 return std::make_tuple(block_id, m_normals_at_pts, v_sdf, m_grad_sdf,
85 template <
typename T1>
88 std::array<double, 3> &unit_ray, std::array<double, 3> &point,
90 std::array<double, 9> &elem_point_nodes,
91 std::array<double, 9> &elem_traction_nodes,
96 auto t_unit_ray = getFTensor1FromPtr<3>(unit_ray.data());
97 auto t_point = getFTensor1FromPtr<3>(point.data());
99 auto get_normal = [](
auto &ele_coords) {
101 Tools::getTriNormal(ele_coords.data(), &t_normal(0));
105 auto t_normal = get_normal(elem_point_nodes);
106 t_normal(
i) /= t_normal.l2();
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);
112 auto gamma = (sn - nm) / nr;
115 t_point_current(
i) = t_spatial_coords(
i) + gamma * t_unit_ray(
i);
117 auto get_local_point_shape_functions = [&](
auto &&t_elem_coords,
119 std::array<T1, 2> loc_coords;
121 Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
122 &t_elem_coords(0, 0), &t_point(0), 1, loc_coords.data()),
125 N_MBTRI1(loc_coords[0], loc_coords[1]),
126 N_MBTRI2(loc_coords[0], loc_coords[1])};
129 auto eval_position = [&](
auto &&t_field,
auto &t_shape_fun) {
133 t_point_out(
i) = t_shape_fun(
j) * t_field(
j,
i);
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);
144 return std::make_tuple(t_slave_point_updated, t_traction_updated, t_normal);
147 template <
typename T1>
160 template <
typename T1>
176 template <
typename T1>
180 auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
182 auto [t_master_point_current, t_master_traction_current, t_master_normal] =
188 t_normal(
i) = t_master_normal(
i) - t_slave_normal(
i);
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);
196 return std::make_tuple(gap, tn_master, tn_slave,
198 t_master_traction_current, t_slave_traction_current);
205 template <
typename T1,
typename T2,
typename T3>
218 t_u(
i) = t_spatial_coords(
i) - t_coords(
i);
225 auto [t_slave_point_current, t_slave_traction_current, t_slave_normal] =
227 auto [t_master_point_current, t_master_traction_current, t_master_normal] =
232 t_normal(
i) = t_master_normal(
i) - t_slave_normal(
i);
237 t_P(
i,
j) = t_normal(
i) * t_normal(
j);
241 constexpr
double beta = 0.5;
248 auto f_min_gap = [
zeta](
auto d,
auto s) {
249 return 0.5 * (
d + s - std::sqrt((
d - s) * (
d - s) +
zeta));
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));
257 auto f_barrier = [alpha1, alpha2, f_min_gap, f_diff_min_gap](
auto g,
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;
267 t_gap_vec(
i) = beta * t_spatial_coords(
i) +
268 (1 - beta) * t_slave_point_current(
i) - t_spatial_coords(
i);
271 -beta * t_master_traction(
i) + (beta - 1) * t_slave_traction_current(
i);
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);
279 t_traction_bar(
i) = t_normal(
i) * barrier;
283 t_Q(
i,
j) * t_master_traction(
j) +
285 t_P(
i,
j) * (t_master_traction(
j) - t_traction_bar(
j));
288 auto is_nan_or_inf = [](
double value) ->
bool {
289 return std::isnan(value) || std::isinf(value);
292 double v = std::complex<double>(t_rhs(
i) * t_rhs(
i)).real();
293 if (is_nan_or_inf(
v)) {
295 MOFEM_LOG(
"SELF", Sev::error) <<
"t_rhs " << t_rhs;
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)
316 commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr),
317 sdfMapRangePtr(sdf_map_range_ptr) {
325 "get alpha contact failed");
329 "get alpha contact failed");
333 "get alpha contact failed");
354 const size_t nb_gauss_pts = getGaussPts().size2();
359 "Wrong number of integration pts %d != %d",
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);
375 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
380 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
381 auto t_normalize_normal = getFTensor1FromMat<3>(m_normals_at_pts);
388 ++t_normalize_normal;
393 auto face_data_vec_ptr =
395 auto face_gauss_pts_it = face_data_vec_ptr->begin();
397 auto nb_base_functions = data.getN().size2();
398 auto t_base = data.getFTensor0N();
399 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
402 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
405 auto check_face_contact = [&]() {
413 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
416 constexpr
double c = 0;
419 if (!
c && check_face_contact()) {
421 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
422 auto t_rhs_tmp =
multiPointRhs(face_data_ptr, t_coords, t_spatial_coords,
424 t_rhs(
i) = t_rhs_tmp(
i);
431 if (ContactOps::sdfPythonWeakPtr.lock()) {
434 t_cP(
i,
j) = (
c * t_grad_sdf_v(
i)) * t_grad_sdf_v(
j);
436 t_rhs(
i) = t_cQ(
i,
j) * t_traction(
j) +
437 (
c * inv_cn * t_sdf_v) * t_grad_sdf_v(
i);
439 t_rhs(
i) = t_traction(
i);
442 t_rhs(
i) = t_traction(
i);
446 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
447 const double alpha = t_w * getMeasure();
450 for (; bb != nbRows / 3; ++bb) {
451 const double beta = alpha * t_base;
452 t_nf(
i) -= beta * t_rhs(
i);
456 for (; bb < nb_base_functions; ++bb)
465 template <AssemblyType A>
468 boost::shared_ptr<std::vector<BrokenBaseSideData>>
469 broken_base_side_data,
470 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
471 boost::shared_ptr<ContactTree> contact_tree_ptr)
472 :
OP(broken_base_side_data), commonDataPtr(common_data_ptr),
473 contactTreePtr(contact_tree_ptr) {}
475 template <AssemblyType A>
485 const size_t nb_gauss_pts = OP::getGaussPts().size2();
488 if (commonDataPtr->contactDisp.size2() != nb_gauss_pts) {
490 "Wrong number of integration pts %d != %d",
491 commonDataPtr->contactDisp.size2(), nb_gauss_pts);
498 auto t_w = OP::getFTensor0IntegrationWeight();
499 auto t_material_normal = OP::getFTensor1NormalsAtGaussPts();
500 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
502 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
503 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
513 auto face_data_vec_ptr =
514 contactTreePtr->findFaceDataVecPtr(OP::getFEEntityHandle());
515 auto face_gauss_pts_it = face_data_vec_ptr->begin();
517 auto nb_base_functions = data.getN().size2() / 3;
518 auto t_base = data.getFTensor1N<3>();
519 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
521 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
522 const double alpha = t_w / 2.;
525 for (; bb != OP::nbRows / 3; ++bb) {
526 const double beta = alpha * t_base(
i) * t_material_normal(
i);
527 t_nf(
i) += beta * t_disp(
i);
531 for (; bb < nb_base_functions; ++bb)
541 const std::string row_field_name,
const std::string col_field_name,
542 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
543 boost::shared_ptr<ContactTree> contact_tree_ptr,
544 boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
547 commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr),
548 sdfMapRangePtr(sdf_map_range_ptr) {
564 auto nb_rows = row_data.getIndices().size();
565 auto nb_cols = col_data.getIndices().size();
568 locMat.resize(nb_rows, nb_cols,
false);
571 if (nb_cols && nb_rows) {
573 auto nb_gauss_pts = getGaussPts().size2();
574 auto t_w = getFTensor0IntegrationWeight();
575 auto t_coords = getFTensor1CoordsAtGaussPts();
576 auto t_disp = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
577 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
580 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
585 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
586 auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
587 auto t_normalized_normal = getFTensor1FromMat<3>(m_normals_at_pts);
597 ++t_normalized_normal;
600 auto face_data_vec_ptr =
602 auto face_gauss_pts_it = face_data_vec_ptr->begin();
604 auto t_row_base = row_data.getFTensor0N();
605 auto nb_face_functions = row_data.getN().size2() / 3;
608 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
610 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
613 auto check_face_contact = [&]() {
623 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
626 constexpr
double c = 0;
629 if (!
c && check_face_contact()) {
631 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
632 constexpr
double eps = std::numeric_limits<float>::epsilon();
633 for (
auto ii = 0; ii < 3; ++ii) {
635 t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2)};
636 t_spatial_coords_cx(ii) +=
eps * 1
i;
640 for (
int jj = 0; jj != 3; ++jj) {
641 auto v = t_rhs_tmp(jj).imag();
642 t_res_dU(jj, ii) =
v /
eps;
650 if (ContactOps::sdfPythonWeakPtr.lock()) {
654 (-
c) * (t_hess_sdf_v(
i,
j) * t_grad_sdf_v(
k) * t_traction(
k) +
655 t_grad_sdf_v(
i) * t_hess_sdf_v(
k,
j) * t_traction(
k))
657 + (
c * inv_cn) * (t_sdf_v * t_hess_sdf_v(
i,
j) +
659 t_grad_sdf_v(
j) * t_grad_sdf_v(
i));
668 auto alpha = t_w * getMeasure();
671 for (; rr != nb_rows / 3; ++rr) {
673 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
674 auto t_col_base = col_data.getFTensor0N(gg, 0);
676 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
677 auto beta = alpha * t_row_base * t_col_base;
678 t_mat(
i,
j) -= beta * t_res_dU(
i,
j);
685 for (; rr < nb_face_functions; ++rr)
695 template <AssemblyType A>
698 std::string row_field_name,
699 boost::shared_ptr<std::vector<BrokenBaseSideData>>
700 broken_base_side_data,
701 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
702 boost::shared_ptr<ContactTree> contact_tree_ptr,
703 boost::shared_ptr<std::map<int, Range>> sdf_map_range_ptr)
704 :
OP(row_field_name, broken_base_side_data, false, false),
710 template <AssemblyType A>
723 auto nb_rows = row_data.getIndices().size();
724 auto nb_cols = col_data.getIndices().size();
727 locMat.resize(nb_rows, nb_cols,
false);
730 if (nb_cols && nb_rows) {
732 const size_t nb_gauss_pts = OP::getGaussPts().size2();
734 auto t_w = OP::getFTensor0IntegrationWeight();
735 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
736 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
737 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
738 auto t_material_normal = OP::getFTensor1NormalsAtGaussPts();
741 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
742 getSdf(
this, commonDataPtr->contactDisp,
743 checkSdf(OP::getFEEntityHandle(), *sdfMapRangePtr),
false);
746 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
758 auto face_data_vec_ptr =
759 contactTreePtr->findFaceDataVecPtr(OP::getFEEntityHandle());
760 auto face_gauss_pts_it = face_data_vec_ptr->begin();
762 auto t_row_base = row_data.getFTensor0N();
763 auto nb_face_functions = row_data.getN().size2();
767 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
769 auto face_data_ptr = contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
772 auto check_face_contact = [&]() {
782 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
785 constexpr
double c = 0;
788 if (!
c && check_face_contact()) {
790 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
791 constexpr
double eps = std::numeric_limits<float>::epsilon();
792 for (
auto ii = 0; ii != 3; ++ii) {
794 t_traction(0), t_traction(1), t_traction(2)};
795 t_traction_cx(ii) +=
eps * 1
i;
799 for (
int jj = 0; jj != 3; ++jj) {
800 auto v = t_rhs_tmp(jj).imag();
801 t_res_dP(jj, ii) =
v /
eps;
807 if (ContactOps::sdfPythonWeakPtr.lock()) {
809 t_cP(
i,
j) = (
c * t_grad_sdf_v(
i)) * t_grad_sdf_v(
j);
811 t_res_dP(
i,
j) = t_cQ(
i,
j);
820 const double alpha = t_w / 2.;
822 for (; rr != nb_rows / 3; ++rr) {
824 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
825 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
827 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
828 auto col_base = t_col_base(
i) * t_material_normal(
i);
829 const double beta = alpha * t_row_base * col_base;
830 t_mat(
i,
j) -= beta * t_res_dP(
i,
j);
837 for (; rr < nb_face_functions; ++rr)
847 template <AssemblyType A>
850 boost::shared_ptr<std::vector<BrokenBaseSideData>>
851 broken_base_side_data,
852 std::string col_field_name,
853 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
854 boost::shared_ptr<ContactTree> contact_tree_ptr)
855 :
OP(col_field_name, broken_base_side_data, true, true),
856 commonDataPtr(common_data_ptr), contactTreePtr(contact_tree_ptr) {
860 template <AssemblyType A>
876 auto nb_rows = row_data.getIndices().size();
877 auto nb_cols = col_data.getIndices().size();
880 locMat.resize(nb_rows, nb_cols,
false);
883 if (nb_cols && nb_rows) {
885 const size_t nb_gauss_pts = OP::getGaussPts().size2();
887 auto t_w = OP::getFTensor0IntegrationWeight();
888 auto t_disp = getFTensor1FromMat<3>(commonDataPtr->contactDisp);
889 auto t_traction = getFTensor1FromMat<3>(commonDataPtr->contactTraction);
890 auto t_coords = OP::getFTensor1CoordsAtGaussPts();
891 auto t_material_normal = OP::getFTensor1NormalsAtGaussPts();
903 auto face_data_vec_ptr =
904 contactTreePtr->findFaceDataVecPtr(OP::getFEEntityHandle());
905 auto face_gauss_pts_it = face_data_vec_ptr->begin();
907 auto t_row_base = row_data.getFTensor1N<3>();
908 auto nb_face_functions = row_data.getN().size2() / 3;
909 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
911 const auto alpha = t_w / 2.;
914 for (; rr != nb_rows / 3; ++rr) {
916 auto row_base = alpha * (t_row_base(
i) * t_material_normal(
i));
918 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
919 auto t_col_base = col_data.getFTensor0N(gg, 0);
921 for (
size_t cc = 0; cc != nb_cols / 3; ++cc) {
922 const auto beta = row_base * t_col_base;
930 for (; rr < nb_face_functions; ++rr)
937 locMat = trans(locMat);
943 boost::shared_ptr<moab::Core> core_mesh_ptr,
944 int max_order, std::map<int, Range> &&body_map)
945 :
Base(m_field, core_mesh_ptr,
"contact"), maxOrder(max_order),
948 auto ref_ele_ptr = boost::make_shared<PostProcGenerateRefMesh<MBTRI>>();
949 ref_ele_ptr->hoNodes =
954 "Error when generating reference element");
956 MOFEM_LOG(
"EP", Sev::inform) <<
"Contact hoNodes " << ref_ele_ptr->hoNodes
960 <<
"Contact maxOrder " << ref_ele_ptr->defMaxLevel;
966 MB_TAG_CREAT | MB_TAG_DENSE,
969 thBodyId, MB_TAG_CREAT | MB_TAG_DENSE,
972 const auto nb_nodes = (
refElementsMap.at(MBTRI)->hoNodes) ? 6 : 3;
992 std::array<double, 3> def_small_x{0., 0., 0.};
994 MB_TAG_CREAT | MB_TAG_DENSE,
997 MB_TAG_CREAT | MB_TAG_DENSE,
1000 std::array<double, 3> def_tractions{0., 0., 0.};
1002 "TRACTION", 3, MB_TYPE_DOUBLE,
thTraction, MB_TAG_CREAT | MB_TAG_DENSE,
1003 &*def_tractions.begin());
1008 CHKERR Base::preProcess();
1015 CHKERR Base::postProcess();
1018 PetscBarrier(
nullptr);
1020 ParallelComm *pcomm_post_proc_mesh =
1022 if (pcomm_post_proc_mesh ==
nullptr)
1027 auto brodacts = [&](
auto &brodacts_ents) {
1031 pcomm_post_proc_mesh->broadcast_entities(
1038 Range brodacts_ents;
1040 CHKERR brodacts(brodacts_ents);
1059 treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
1069 boost::shared_ptr<ContactTree> contact_tree_ptr,
1070 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
1071 boost::shared_ptr<MatrixDouble> u_h1_ptr)
1072 :
UOP(
NOSPACE,
UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
1073 uH1Ptr(u_h1_ptr), commonDataPtr(common_data_ptr) {}
1079 auto get_body_id = [&](
auto fe_ent) {
1081 if (
m.second.find(fe_ent) !=
m.second.end()) {
1091 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1093 auto body_id = get_body_id(fe_ent);
1097 post_proc_ents, &fe_id);
1099 post_proc_ents, &body_id);
1101 auto nb_gauss_pts = getGaussPts().size2();
1102 auto t_u_h1 = getFTensor1FromMat<3>(*
uH1Ptr);
1103 auto t_u_l2 = getFTensor1FromMat<3>(
commonDataPtr->contactDisp);
1104 auto t_coords = getFTensor1CoordsAtGaussPts();
1107 auto t_x_h1 = getFTensor1FromPtr<3>(&x_h1(0, 0));
1109 auto t_x_l2 = getFTensor1FromPtr<3>(&x_l2(0, 0));
1116 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1117 t_x_h1(
i) = t_coords(
i) + t_u_h1(
i);
1118 t_x_l2(
i) = t_coords(
i) + t_u_l2(
i);
1127 CHKERR moab_post_proc_mesh.set_coords(
1128 &*map_gauss_pts.begin(), map_gauss_pts.size(), &*x_h1.data().begin());
1129 CHKERR moab_post_proc_mesh.tag_set_data(
1130 contactTreePtr->thSmallX, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1131 &*x_h1.data().begin());
1132 CHKERR moab_post_proc_mesh.tag_set_data(
1133 contactTreePtr->thLargeX, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1134 &*coords.data().begin());
1135 CHKERR moab_post_proc_mesh.tag_set_data(
1136 contactTreePtr->thTraction, &*map_gauss_pts.begin(), map_gauss_pts.size(),
1137 &*tractions.data().begin());
1143 boost::shared_ptr<ContactTree> contact_tree_ptr,
1144 boost::shared_ptr<ContactOps::CommonData> common_data_ptr,
1145 boost::shared_ptr<MatrixDouble> u_h1_ptr,
Range r,
1148 std::vector<EntityHandle> *map_gauss_pts_ptr
1151 :
UOP(
NOSPACE,
UOP::OPSPACE), contactTreePtr(contact_tree_ptr),
1152 commonDataPtr(common_data_ptr), uH1Ptr(u_h1_ptr),
1153 postProcMeshPtr(post_proc_mesh_ptr), mapGaussPtsPtr(map_gauss_pts_ptr),
1160 auto &m_field = getPtrFE()->mField;
1161 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1170 const auto nb_gauss_pts = getGaussPts().size2();
1172 auto t_disp_h1 = getFTensor1FromMat<3>(*
uH1Ptr);
1173 auto t_coords = getFTensor1CoordsAtGaussPts();
1174 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
1182 auto get_ele_centre = [
i](
auto t_ele_coords) {
1184 t_ele_center(
i) = 0;
1185 for (
int nn = 0; nn != 3; nn++) {
1186 t_ele_center(
i) += t_ele_coords(
i);
1189 t_ele_center(
i) /= 3;
1190 return t_ele_center;
1193 auto get_ele_radius = [
i](
auto t_ele_center,
auto t_ele_coords) {
1195 t_n0(
i) = t_ele_center(
i) - t_ele_coords(
i);
1199 auto get_face_conn = [
this](
auto face) {
1203 face, conn, num_nodes,
true),
1205 if (num_nodes != 3) {
1211 auto get_face_coords = [
this](
auto conn) {
1212 std::array<double, 9> coords;
1217 auto get_closet_face = [
this](
auto *point_ptr,
auto r) {
1219 std::vector<EntityHandle> faces_out;
1223 "get closest faces");
1227 auto get_faces_out = [
this](
auto *point_ptr,
auto *unit_ray_ptr,
auto radius,
1229 std::vector<double> distances_out;
1230 std::vector<EntityHandle> faces_out;
1235 point_ptr, unit_ray_ptr, &radius),
1237 "get closest faces");
1238 return std::make_pair(faces_out, distances_out);
1241 auto get_normal = [](
auto &ele_coords) {
1243 Tools::getTriNormal(ele_coords.data(), &t_normal(0));
1247 auto make_map = [&](
auto &face_out,
auto &face_dist,
auto &t_ray_point,
1248 auto &t_unit_ray,
auto &t_master_coord) {
1251 std::map<double, EntityHandle>
m;
1252 for (
auto ii = 0; ii != face_out.size(); ++ii) {
1253 auto face_conn = get_face_conn(face_out[ii]);
1256 t_face_normal.normalize();
1258 t_x(
i) = t_ray_point(
i) + t_unit_ray(
i) * face_dist[ii];
1261 t_x(
i) * t_face_normal(
j) - t_master_coord(
i) * t_unit_ray(
j);
1262 if (t_unit_ray(
i) * t_face_normal(
i) > std::cos(M_PI / 3)) {
1263 auto dot = std::sqrt(t_chi(
i,
j) * t_chi(
i,
j));
1264 m[dot] = face_out[ii];
1270 auto get_tag_data = [
this](
auto tag,
auto face,
auto &vec) {
1274 vec.resize(tag_length);
1280 auto create_tag = [
this](
const std::string tag_name,
const int size) {
1281 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1285 tag_name.c_str(), size, MB_TYPE_DOUBLE,
th,
1286 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1293 auto set_float_precision = [](
const double x) {
1294 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1301 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
1304 v = set_float_precision(
v);
1311 auto get_fe_adjacencies = [
this](
auto fe_ent) {
1314 &fe_ent, 1, 2,
false, adj_faces, moab::Interface::UNION),
1316 std::set<int> adj_ids;
1317 for (
auto f : adj_faces) {
1323 auto get_face_id = [
this](
auto face) {
1332 auto get_body_id = [
this](
auto face) {
1341 auto get_face_part = [
this](
auto face) {
1342 ParallelComm *pcomm_post_proc_mesh = ParallelComm::get_pcomm(
1346 pcomm_post_proc_mesh->part_tag(), &face, 1, &part) == MB_SUCCESS) {
1352 auto check_face = [&](
auto face,
auto fe_id,
auto part) {
1353 auto face_id = get_face_id(face);
1354 auto face_part = get_face_part(face);
1355 if (face_id == fe_id && face_part == part)
1363 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1367 for (
auto &
a :
v.data())
1368 a = set_float_precision(
a);
1370 &*
v.data().begin());
1375 Tag th_mark = create_tag(
"contact_mark", 1);
1376 Tag th_mark_slave = create_tag(
"contact_mark_slave", 1);
1377 Tag th_body_id = create_tag(
"contact_body_id", 1);
1378 Tag th_gap = create_tag(
"contact_gap", 1);
1379 Tag th_tn_master = create_tag(
"contact_tn_master", 1);
1380 Tag th_tn_slave = create_tag(
"contact_tn_slave", 1);
1381 Tag th_contact_traction = create_tag(
"contact_traction", 3);
1382 Tag th_contact_traction_master = create_tag(
"contact_traction_master", 3);
1383 Tag th_contact_traction_slave = create_tag(
"contact_traction_slave", 3);
1384 Tag th_c = create_tag(
"contact_c", 1);
1385 Tag th_normal = create_tag(
"contact_normal", 3);
1386 Tag th_dist = create_tag(
"contact_dip", 3);
1388 auto t_ele_centre = get_ele_centre(getFTensor1Coords());
1389 auto ele_radius = get_ele_radius(t_ele_centre, getFTensor1Coords());
1395 auto adj_fe_ids = get_fe_adjacencies(fe_ent);
1397 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1400 t_spatial_coords(
i) = t_coords(
i) + t_disp_h1(
i);
1403 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1406 auto faces_close = get_closet_face(&t_spatial_coords(0), ele_radius);
1407 for (
auto face_close : faces_close) {
1408 if (check_face(face_close, fe_id, m_field.get_comm_rank())) {
1410 auto body_id = get_body_id(face_close);
1412 auto master_face_conn = get_face_conn(face_close);
1413 std::array<double, 9> master_coords;
1416 master_coords.data());
1417 std::array<double, 9> master_traction;
1420 master_traction.data());
1421 auto t_normal_face_close = get_normal(master_coords);
1422 t_normal_face_close.normalize();
1426 CHKERR save_scal_tag(th_mark,
m, gg);
1427 CHKERR save_scal_tag(th_body_id,
static_cast<double>(body_id), gg);
1428 CHKERR save_vec_tag(th_normal, t_normal_face_close, gg);
1429 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1433 t_unit_ray(
i) = -t_normal_face_close(
i);
1436 t_spatial_coords(
i) -
1439 constexpr
double eps = 1e-3;
1440 auto [faces_out, faces_dist] =
1441 get_faces_out(&t_ray_point(0), &t_unit_ray(0),
1445 auto m = make_map(faces_out, faces_dist, t_ray_point, t_unit_ray,
1447 for (
auto m_it =
m.begin(); m_it !=
m.end(); ++m_it) {
1448 auto face = m_it->second;
1449 if (face != face_close) {
1453 body_id != get_body_id(face) &&
1454 (adj_fe_ids.find(get_face_id(face)) == adj_fe_ids.end() ||
1455 get_face_part(face) != m_field.get_comm_rank())
1460 shadow_vec.back().gaussPtNb = gg;
1462 auto slave_face_conn = get_face_conn(face);
1463 std::array<double, 9> slave_coords;
1466 slave_coords.data());
1467 auto t_normal_face = get_normal(slave_coords);
1468 std::array<double, 9> slave_tractions;
1471 slave_tractions.data());
1473 auto t_master_point =
1474 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data());
1475 auto t_slave_point =
1476 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data());
1477 auto t_ray_point_data =
1478 getFTensor1FromPtr<3>(shadow_vec.back().rayPoint.data());
1479 auto t_unit_ray_data =
1480 getFTensor1FromPtr<3>(shadow_vec.back().unitRay.data());
1482 t_slave_point(
i) = t_ray_point(
i) + m_it->first * t_unit_ray(
i);
1484 auto eval_position = [&](
auto &&t_elem_coords,
auto &&t_point) {
1485 std::array<double, 2> loc_coords;
1487 Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
1488 &t_elem_coords(0, 0), &t_point(0), 1,
1490 "get local coords");
1499 t_point_out(
i) = t_shape_fun(
j) * t_elem_coords(
j,
i);
1503 auto t_master_point_updated = eval_position(
1504 getFTensor2FromPtr<3, 3>(master_coords.data()),
1505 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data()));
1506 t_master_point(
i) = t_master_point_updated(
i);
1508 auto t_slave_point_updated = eval_position(
1509 getFTensor2FromPtr<3, 3>(slave_coords.data()),
1510 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data()));
1511 t_slave_point(
i) = t_slave_point_updated(
i);
1513 t_ray_point_data(
i) = t_ray_point(
i);
1514 t_unit_ray_data(
i) = t_unit_ray(
i);
1516 std::copy(master_coords.begin(), master_coords.end(),
1517 shadow_vec.back().masterPointNodes.begin());
1518 std::copy(master_traction.begin(), master_traction.end(),
1519 shadow_vec.back().masterTractionNodes.begin());
1520 std::copy(slave_coords.begin(), slave_coords.end(),
1521 shadow_vec.back().slavePointNodes.begin());
1522 std::copy(slave_tractions.begin(), slave_tractions.end(),
1523 shadow_vec.back().slaveTractionNodes.begin());
1525 shadow_vec.back().eleRadius = ele_radius;
1535 auto [gap, tn_master, tn_slave,
c, t_master_traction,
1537 multiGetGap(&(shadow_vec.back()), t_spatial_coords);
1539 t_gap_vec(
i) = t_slave_point(
i) - t_spatial_coords(
i);
1540 CHKERR save_scal_tag(th_gap, gap, gg);
1541 CHKERR save_scal_tag(th_tn_master, tn_master, gg);
1542 CHKERR save_scal_tag(th_tn_slave, tn_slave, gg);
1543 CHKERR save_scal_tag(th_c,
c, gg);
1545 CHKERR save_scal_tag(th_mark_slave,
m, gg);
1546 CHKERR save_vec_tag(th_dist, t_gap_vec, gg);
1547 CHKERR save_vec_tag(th_contact_traction_master,
1548 t_master_traction, gg);
1549 CHKERR save_vec_tag(th_contact_traction_slave, t_slave_traction,