1145 auto &m_field = getPtrFE()->mField;
1146 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1155 const auto nb_gauss_pts = getGaussPts().size2();
1157 auto t_disp_h1 = getFTensor1FromMat<3>(*
uH1Ptr);
1158 auto t_coords = getFTensor1CoordsAtGaussPts();
1159 auto t_traction = getFTensor1FromMat<3>(
commonDataPtr->contactTraction);
1167 auto get_ele_centre = [
i](
auto t_ele_coords) {
1169 t_ele_center(
i) = 0;
1170 for (
int nn = 0; nn != 3; nn++) {
1171 t_ele_center(
i) += t_ele_coords(
i);
1174 t_ele_center(
i) /= 3;
1175 return t_ele_center;
1178 auto get_ele_radius = [
i](
auto t_ele_center,
auto t_ele_coords) {
1180 t_n0(
i) = t_ele_center(
i) - t_ele_coords(
i);
1184 auto get_face_conn = [
this](
auto face) {
1188 face, conn, num_nodes,
true),
1190 if (num_nodes != 3) {
1196 auto get_face_coords = [
this](
auto conn) {
1197 std::array<double, 9> coords;
1202 auto get_closet_face = [
this](
auto *point_ptr,
auto r) {
1204 std::vector<EntityHandle> faces_out;
1208 "get closest faces");
1212 auto get_faces_out = [
this](
auto *point_ptr,
auto *unit_ray_ptr,
auto radius,
1214 std::vector<double> distances_out;
1215 std::vector<EntityHandle> faces_out;
1220 point_ptr, unit_ray_ptr, &radius),
1222 "get closest faces");
1223 return std::make_pair(faces_out, distances_out);
1226 auto get_normal = [](
auto &ele_coords) {
1228 Tools::getTriNormal(ele_coords.data(), &t_normal(0));
1232 auto make_map = [&](
auto &face_out,
auto &face_dist,
auto &t_ray_point,
1233 auto &t_unit_ray,
auto &t_master_coord) {
1236 std::map<double, EntityHandle>
m;
1237 for (
auto ii = 0; ii != face_out.size(); ++ii) {
1238 auto face_conn = get_face_conn(face_out[ii]);
1241 t_face_normal.normalize();
1243 t_x(
i) = t_ray_point(
i) + t_unit_ray(
i) * face_dist[ii];
1246 t_x(
i) * t_face_normal(
j) - t_master_coord(
i) * t_unit_ray(
j);
1247 if (t_unit_ray(
i) * t_face_normal(
i) > std::cos(M_PI / 3)) {
1248 auto dot = std::sqrt(t_chi(
i,
j) * t_chi(
i,
j));
1249 m[dot] = face_out[ii];
1255 auto get_tag_data = [
this](
auto tag,
auto face,
auto &vec) {
1259 vec.resize(tag_length);
1265 auto create_tag = [
this](
const std::string tag_name,
const int size) {
1266 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1270 tag_name.c_str(), size, MB_TYPE_DOUBLE,
th,
1271 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1278 auto set_float_precision = [](
const double x) {
1279 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1286 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
1289 v = set_float_precision(
v);
1296 auto get_fe_adjacencies = [
this](
auto fe_ent) {
1299 &fe_ent, 1, 2,
false, adj_faces, moab::Interface::UNION),
1301 std::set<int> adj_ids;
1302 for (
auto f : adj_faces) {
1308 auto get_face_id = [
this](
auto face) {
1317 auto get_body_id = [
this](
auto face) {
1326 auto get_face_part = [
this](
auto face) {
1327 ParallelComm *pcomm_post_proc_mesh = ParallelComm::get_pcomm(
1331 pcomm_post_proc_mesh->part_tag(), &face, 1, &part) == MB_SUCCESS) {
1337 auto check_face = [&](
auto face,
auto fe_id,
auto part) {
1338 auto face_id = get_face_id(face);
1339 auto face_part = get_face_part(face);
1340 if (face_id == fe_id && face_part == part)
1348 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1352 for (
auto &
a :
v.data())
1353 a = set_float_precision(
a);
1355 &*
v.data().begin());
1360 Tag th_mark = create_tag(
"contact_mark", 1);
1361 Tag th_mark_slave = create_tag(
"contact_mark_slave", 1);
1362 Tag th_body_id = create_tag(
"contact_body_id", 1);
1363 Tag th_gap = create_tag(
"contact_gap", 1);
1364 Tag th_tn_master = create_tag(
"contact_tn_master", 1);
1365 Tag th_tn_slave = create_tag(
"contact_tn_slave", 1);
1366 Tag th_contact_traction = create_tag(
"contact_traction", 3);
1367 Tag th_contact_traction_master = create_tag(
"contact_traction_master", 3);
1368 Tag th_contact_traction_slave = create_tag(
"contact_traction_slave", 3);
1369 Tag th_c = create_tag(
"contact_c", 1);
1370 Tag th_normal = create_tag(
"contact_normal", 3);
1371 Tag th_dist = create_tag(
"contact_dip", 3);
1373 auto t_ele_centre = get_ele_centre(getFTensor1Coords());
1374 auto ele_radius = get_ele_radius(t_ele_centre, getFTensor1Coords());
1380 auto adj_fe_ids = get_fe_adjacencies(fe_ent);
1382 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
1385 t_spatial_coords(
i) = t_coords(
i) + t_disp_h1(
i);
1388 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1391 auto faces_close = get_closet_face(&t_spatial_coords(0), ele_radius);
1392 for (
auto face_close : faces_close) {
1393 if (check_face(face_close, fe_id, m_field.get_comm_rank())) {
1395 auto body_id = get_body_id(face_close);
1397 auto master_face_conn = get_face_conn(face_close);
1398 std::array<double, 9> master_coords;
1401 master_coords.data());
1402 std::array<double, 9> master_traction;
1405 master_traction.data());
1406 auto t_normal_face_close = get_normal(master_coords);
1407 t_normal_face_close.normalize();
1411 CHKERR save_scal_tag(th_mark,
m, gg);
1412 CHKERR save_scal_tag(th_body_id,
static_cast<double>(body_id), gg);
1413 CHKERR save_vec_tag(th_normal, t_normal_face_close, gg);
1414 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1418 t_unit_ray(
i) = -t_normal_face_close(
i);
1421 t_spatial_coords(
i) -
1424 constexpr
double eps = 1e-3;
1425 auto [faces_out, faces_dist] =
1426 get_faces_out(&t_ray_point(0), &t_unit_ray(0),
1430 auto m = make_map(faces_out, faces_dist, t_ray_point, t_unit_ray,
1432 for (
auto m_it =
m.begin(); m_it !=
m.end(); ++m_it) {
1433 auto face = m_it->second;
1434 if (face != face_close) {
1438 body_id != get_body_id(face) &&
1439 (adj_fe_ids.find(get_face_id(face)) == adj_fe_ids.end() ||
1440 get_face_part(face) != m_field.get_comm_rank())
1444 shadow_vec.push_back(ContactTree::FaceData());
1445 shadow_vec.back().gaussPtNb = gg;
1447 auto slave_face_conn = get_face_conn(face);
1448 std::array<double, 9> slave_coords;
1451 slave_coords.data());
1452 auto t_normal_face = get_normal(slave_coords);
1453 std::array<double, 9> slave_tractions;
1456 slave_tractions.data());
1458 auto t_master_point =
1459 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data());
1460 auto t_slave_point =
1461 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data());
1462 auto t_ray_point_data =
1463 getFTensor1FromPtr<3>(shadow_vec.back().rayPoint.data());
1464 auto t_unit_ray_data =
1465 getFTensor1FromPtr<3>(shadow_vec.back().unitRay.data());
1467 t_slave_point(
i) = t_ray_point(
i) + m_it->first * t_unit_ray(
i);
1469 auto eval_position = [&](
auto &&t_elem_coords,
auto &&t_point) {
1470 std::array<double, 2> loc_coords;
1472 Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
1473 &t_elem_coords(0, 0), &t_point(0), 1,
1475 "get local coords");
1484 t_point_out(
i) = t_shape_fun(
j) * t_elem_coords(
j,
i);
1488 auto t_master_point_updated = eval_position(
1489 getFTensor2FromPtr<3, 3>(master_coords.data()),
1490 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data()));
1491 t_master_point(
i) = t_master_point_updated(
i);
1493 auto t_slave_point_updated = eval_position(
1494 getFTensor2FromPtr<3, 3>(slave_coords.data()),
1495 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data()));
1496 t_slave_point(
i) = t_slave_point_updated(
i);
1498 t_ray_point_data(
i) = t_ray_point(
i);
1499 t_unit_ray_data(
i) = t_unit_ray(
i);
1501 std::copy(master_coords.begin(), master_coords.end(),
1502 shadow_vec.back().masterPointNodes.begin());
1503 std::copy(master_traction.begin(), master_traction.end(),
1504 shadow_vec.back().masterTractionNodes.begin());
1505 std::copy(slave_coords.begin(), slave_coords.end(),
1506 shadow_vec.back().slavePointNodes.begin());
1507 std::copy(slave_tractions.begin(), slave_tractions.end(),
1508 shadow_vec.back().slaveTractionNodes.begin());
1510 shadow_vec.back().eleRadius = ele_radius;
1520 auto [gap, tn_master, tn_slave,
c, t_master_traction,
1522 multiGetGap(&(shadow_vec.back()), t_spatial_coords);
1524 t_gap_vec(
i) = t_slave_point(
i) - t_spatial_coords(
i);
1525 CHKERR save_scal_tag(th_gap, gap, gg);
1526 CHKERR save_scal_tag(th_tn_master, tn_master, gg);
1527 CHKERR save_scal_tag(th_tn_slave, tn_slave, gg);
1528 CHKERR save_scal_tag(th_c,
c, gg);
1530 CHKERR save_scal_tag(th_mark_slave,
m, gg);
1531 CHKERR save_vec_tag(th_dist, t_gap_vec, gg);
1532 CHKERR save_vec_tag(th_contact_traction_master,
1533 t_master_traction, gg);
1534 CHKERR save_vec_tag(th_contact_traction_slave, t_slave_traction,