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);
354 const size_t nb_gauss_pts = getGaussPts().size2();
359 "Wrong number of integration pts %ld != %ld",
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] =
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);
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 = [&]() {
416#ifdef ENABLE_PYTHON_BINDING
418 if (ContactOps::sdfPythonWeakPtr.lock()) {
419 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
423 constexpr double c = 0;
426 if (!
c && check_face_contact()) {
428 t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
430 multiPointRhs(face_data_ptr, t_coords, t_spatial_coords, t_traction,
432 t_rhs(
i) = t_rhs_tmp(
i);
436#ifdef ENABLE_PYTHON_BINDING
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);
446 t_rhs(
i) = t_traction(
i);
449 t_rhs(
i) = t_traction(
i);
453 auto t_nf = getFTensor1FromPtr<3>(&nf[0]);
454 const double alpha = t_w * getMeasure();
457 for (; bb != nbRows / 3; ++bb) {
458 const double beta = alpha * t_base;
459 t_nf(
i) -= beta * t_rhs(
i);
463 for (; bb < nb_base_functions; ++bb)
562 EntitiesFieldData::EntData &col_data) {
571 auto nb_rows = row_data.getIndices().size();
572 auto nb_cols = col_data.getIndices().size();
574 auto &locMat = AssemblyBoundaryEleOp::locMat;
575 locMat.resize(nb_rows, nb_cols,
false);
578 if (nb_cols && nb_rows) {
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);
587 auto [block_id, m_normals_at_pts, v_sdf, m_grad_sdf, m_hess_sdf] =
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);
604 ++t_normalized_normal;
607 auto face_data_vec_ptr =
609 auto face_gauss_pts_it = face_data_vec_ptr->begin();
611 auto t_row_base = row_data.getFTensor0N();
612 auto nb_face_functions = row_data.getN().size2() / 3;
615 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
617 auto face_data_ptr =
contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
620 auto check_face_contact = [&]() {
632#ifdef ENABLE_PYTHON_BINDING
634 if (ContactOps::sdfPythonWeakPtr.lock()) {
635 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
639 constexpr double c = 0;
642 if (!
c && check_face_contact()) {
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) {
648 t_spatial_coords(0), t_spatial_coords(1), t_spatial_coords(2)};
649 t_spatial_coords_cx(ii) +=
eps * 1
i;
653 for (
int jj = 0; jj != 3; ++jj) {
654 auto v = t_rhs_tmp(jj).imag();
655 t_res_dU(jj, ii) =
v /
eps;
661#ifdef ENABLE_PYTHON_BINDING
663 if (ContactOps::sdfPythonWeakPtr.lock()) {
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))
670 + (
c * inv_cn) * (t_sdf_v * t_hess_sdf_v(
i,
j) +
672 t_grad_sdf_v(
j) * t_grad_sdf_v(
i));
681 auto alpha = t_w * getMeasure();
684 for (; rr != nb_rows / 3; ++rr) {
686 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
687 auto t_col_base = col_data.getFTensor0N(gg, 0);
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);
698 for (; rr < nb_face_functions; ++rr)
726 EntitiesFieldData::EntData &row_data,
727 EntitiesFieldData::EntData &col_data) {
736 auto nb_rows = row_data.getIndices().size();
737 auto nb_cols = col_data.getIndices().size();
739 auto &locMat = AssemblyBoundaryEleOp::locMat;
740 locMat.resize(nb_rows, nb_cols,
false);
743 if (nb_cols && nb_rows) {
745 const size_t nb_gauss_pts = OP::getGaussPts().size2();
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();
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);
758 auto t_sdf_v = getFTensor0FromVec(v_sdf);
759 auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
771 auto face_data_vec_ptr =
772 contactTreePtr->findFaceDataVecPtr(OP::getFEEntityHandle());
773 auto face_gauss_pts_it = face_data_vec_ptr->begin();
775 auto t_row_base = row_data.getFTensor0N();
776 auto nb_face_functions = row_data.getN().size2();
780 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
782 auto face_data_ptr = contactTreePtr->getFaceDataPtr(face_gauss_pts_it, gg,
785 auto check_face_contact = [&]() {
786 if (
checkSdf(OP::getFEEntityHandle(), *sdfMapRangePtr) != -1)
797#ifdef ENABLE_PYTHON_BINDING
799 if (ContactOps::sdfPythonWeakPtr.lock()) {
800 auto tn = t_traction(
i) * t_grad_sdf_v(
i);
804 constexpr double c = 0;
807 if (!
c && check_face_contact()) {
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 * 1
i;
818 for (
int jj = 0; jj != 3; ++jj) {
819 auto v = t_rhs_tmp(jj).imag();
820 t_res_dP(jj, ii) =
v /
eps;
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);
839 const double alpha = t_w / 2.;
841 for (; rr != nb_rows / 3; ++rr) {
843 auto t_mat = getFTensor2FromArray<3, 3, 3>(locMat, 3 * rr);
844 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
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);
856 for (; rr < nb_face_functions; ++rr)
1092 EntitiesFieldData::EntData &data) {
1095 auto get_body_id = [&](
auto fe_ent) {
1097 if (
m.second.find(fe_ent) !=
m.second.end()) {
1107 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1108 auto fe_id = id_from_handle(fe_ent);
1109 auto body_id = get_body_id(fe_ent);
1113 post_proc_ents, &fe_id);
1115 post_proc_ents, &body_id);
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();
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();
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);
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());
1173 EntitiesFieldData::EntData &data) {
1176 auto &m_field = getPtrFE()->mField;
1177 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1178 auto fe_id = id_from_handle(fe_ent);
1186 const auto nb_gauss_pts = getGaussPts().size2();
1188 auto t_disp_h1 = getFTensor1FromMat<3>(*
uH1Ptr);
1189 auto t_coords = getFTensor1CoordsAtGaussPts();
1190 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
1198 auto get_ele_centre = [
i](
auto t_ele_coords) {
1200 t_ele_center(
i) = 0;
1201 for (
int nn = 0; nn != 3; nn++) {
1202 t_ele_center(
i) += t_ele_coords(
i);
1205 t_ele_center(
i) /= 3;
1206 return t_ele_center;
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);
1215 auto get_face_conn = [
this](
auto face) {
1219 face, conn, num_nodes,
true),
1221 if (num_nodes != 3) {
1227 auto get_face_coords = [
this](
auto conn) {
1228 std::array<double, 9> coords;
1233 auto get_closet_face = [
this](
auto *point_ptr,
auto r) {
1235 std::vector<EntityHandle> faces_out;
1239 "get closest faces");
1243 auto get_faces_out = [
this](
auto *point_ptr,
auto *unit_ray_ptr,
auto radius,
1245 std::vector<double> distances_out;
1246 std::vector<EntityHandle> faces_out;
1251 point_ptr, unit_ray_ptr, &radius),
1253 "get closest faces");
1254 return std::make_pair(faces_out, distances_out);
1257 auto get_normal = [](
auto &ele_coords) {
1259 Tools::getTriNormal(ele_coords.data(), &t_normal(0));
1263 auto make_map = [&](
auto &face_out,
auto &face_dist,
auto &t_ray_point,
1264 auto &t_unit_ray,
auto &t_master_coord) {
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]);
1272 t_face_normal.normalize();
1274 t_x(
i) = t_ray_point(
i) + t_unit_ray(
i) * face_dist[ii];
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];
1286 auto get_tag_data = [
this](
auto tag,
auto face,
auto &vec) {
1290 vec.resize(tag_length);
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};
1301 tag_name.c_str(), size, MB_TYPE_DOUBLE, th,
1302 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1309 auto set_float_precision = [](
const double x) {
1310 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1317 auto save_scal_tag = [&](
auto &th,
auto v,
const int gg) {
1320 v = set_float_precision(
v);
1327 auto get_fe_adjacencies = [
this](
auto fe_ent) {
1330 &fe_ent, 1, 2,
false, adj_faces, moab::Interface::UNION),
1332 std::set<int> adj_ids;
1333 for (
auto f : adj_faces) {
1334 adj_ids.insert(id_from_handle(f));
1339 auto get_face_id = [
this](
auto face) {
1348 auto get_body_id = [
this](
auto face) {
1357 auto get_face_part = [
this](
auto face) {
1358 ParallelComm *pcomm_post_proc_mesh = ParallelComm::get_pcomm(
1362 pcomm_post_proc_mesh->part_tag(), &face, 1, &part) == MB_SUCCESS) {
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)
1379 auto save_vec_tag = [&](
auto &th,
auto &t_d,
const int gg) {
1383 for (
auto &
a :
v.data())
1384 a = set_float_precision(
a);
1386 &*
v.data().begin());
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);
1404 auto t_ele_centre = get_ele_centre(getFTensor1Coords());
1405 auto ele_radius = get_ele_radius(t_ele_centre, getFTensor1Coords());
1411 auto adj_fe_ids = get_fe_adjacencies(fe_ent);
1413 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1416 t_spatial_coords(
i) = t_coords(
i) + t_disp_h1(
i);
1419 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
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())) {
1426 auto body_id = get_body_id(face_close);
1428 auto master_face_conn = get_face_conn(face_close);
1429 std::array<double, 9> master_coords;
1432 master_coords.data());
1433 std::array<double, 9> master_traction;
1436 master_traction.data());
1437 auto t_normal_face_close = get_normal(master_coords);
1438 t_normal_face_close.normalize();
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);
1449 t_unit_ray(
i) = -t_normal_face_close(
i);
1452 t_spatial_coords(
i) -
1455 constexpr double eps = 1e-3;
1456 auto [faces_out, faces_dist] =
1457 get_faces_out(&t_ray_point(0), &t_unit_ray(0),
1461 auto m = make_map(faces_out, faces_dist, t_ray_point, t_unit_ray,
1463 for (
auto m_it =
m.begin(); m_it !=
m.end(); ++m_it) {
1464 auto face = m_it->second;
1465 if (face != face_close) {
1469 (adj_fe_ids.find(get_face_id(face)) == adj_fe_ids.end() ||
1470 get_face_part(face) != m_field.get_comm_rank())
1475 shadow_vec.back().gaussPtNb = gg;
1477 auto slave_face_conn = get_face_conn(face);
1478 std::array<double, 9> slave_coords;
1481 slave_coords.data());
1482 auto t_normal_face = get_normal(slave_coords);
1483 std::array<double, 9> slave_tractions;
1486 slave_tractions.data());
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());
1497 t_slave_point(
i) = t_ray_point(
i) + m_it->first * t_unit_ray(
i);
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,
1505 "get local coords");
1514 t_point_out(
i) = t_shape_fun(
j) * t_elem_coords(
j,
i);
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);
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);
1528 t_ray_point_data(
i) = t_ray_point(
i);
1529 t_unit_ray_data(
i) = t_unit_ray(
i);
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());
1540 shadow_vec.back().eleRadius = ele_radius;
1550 auto [gap, tn_master, tn_slave,
c, t_master_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);
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,