v0.15.0
Loading...
Searching...
No Matches
EshelbianPlasticity::OpTreeSearch Struct Reference

#include "users_modules/eshelbian_plasticity/src/EshelbianContact.hpp"

Inheritance diagram for EshelbianPlasticity::OpTreeSearch:
[legend]
Collaboration diagram for EshelbianPlasticity::OpTreeSearch:
[legend]

Public Types

using UOP = FaceElementForcesAndSourcesCore::UserDataOperator
 

Public Member Functions

 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)
 
MoFEMErrorCode doWork (int side, EntityType type, EntitiesFieldData::EntData &data)
 

Protected Attributes

boost::shared_ptr< ContactTreecontactTreePtr
 
boost::shared_ptr< ContactOps::CommonDatacommonDataPtr
 
boost::shared_ptr< MatrixDouble > uH1Ptr
 
moab::Interface * postProcMeshPtr = nullptr
 
std::vector< EntityHandle > * mapGaussPtsPtr = nullptr
 
Range contactRange
 

Detailed Description

Definition at line 237 of file EshelbianContact.hpp.

Member Typedef Documentation

◆ UOP

using EshelbianPlasticity::OpTreeSearch::UOP = FaceElementForcesAndSourcesCore::UserDataOperator

Definition at line 239 of file EshelbianContact.hpp.

Constructor & Destructor Documentation

◆ OpTreeSearch()

EshelbianPlasticity::OpTreeSearch::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 )

Definition at line 1158 of file EshelbianContact.cpp.

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) {}
@ NOSPACE
Definition definitions.h:83
boost::shared_ptr< ContactTree > contactTreePtr
std::vector< EntityHandle > * mapGaussPtsPtr
boost::shared_ptr< MatrixDouble > uH1Ptr
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
FaceElementForcesAndSourcesCore::UserDataOperator UOP

Member Function Documentation

◆ doWork()

MoFEMErrorCode EshelbianPlasticity::OpTreeSearch::doWork ( int side,
EntityType type,
EntitiesFieldData::EntData & data )

Definition at line 1172 of file EshelbianContact.cpp.

1173 {
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}
constexpr double a
static const double eps
#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()
#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 double face_coords[4][9]
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< 'j', 3 > j
double airplane_ray_distance
auto multiGetGap(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
auto id_from_handle(const EntityHandle h)
int r
Definition sdf.py:8
FTensor::Index< 'm', 3 > m

Member Data Documentation

◆ commonDataPtr

boost::shared_ptr<ContactOps::CommonData> EshelbianPlasticity::OpTreeSearch::commonDataPtr
protected

Definition at line 254 of file EshelbianContact.hpp.

◆ contactRange

Range EshelbianPlasticity::OpTreeSearch::contactRange
protected

Definition at line 260 of file EshelbianContact.hpp.

◆ contactTreePtr

boost::shared_ptr<ContactTree> EshelbianPlasticity::OpTreeSearch::contactTreePtr
protected

Definition at line 253 of file EshelbianContact.hpp.

◆ mapGaussPtsPtr

std::vector<EntityHandle>* EshelbianPlasticity::OpTreeSearch::mapGaussPtsPtr = nullptr
protected

Definition at line 258 of file EshelbianContact.hpp.

◆ postProcMeshPtr

moab::Interface* EshelbianPlasticity::OpTreeSearch::postProcMeshPtr = nullptr
protected

Definition at line 257 of file EshelbianContact.hpp.

◆ uH1Ptr

boost::shared_ptr<MatrixDouble> EshelbianPlasticity::OpTreeSearch::uH1Ptr
protected

Definition at line 255 of file EshelbianContact.hpp.


The documentation for this struct was generated from the following files: