1498 {
1500
1501 auto &m_field = getPtrFE()->mField;
1502 auto fe_ent = getNumeredEntFiniteElementPtr()->getEnt();
1504
1507
1510
1511 const auto nb_gauss_pts = getGaussPts().size2();
1512
1513 auto t_disp_h1 = getFTensor1FromMat<3>(*
uH1Ptr);
1514 auto t_coords = getFTensor1CoordsAtGaussPts();
1515 auto t_traction = getFTensor1FromMat<3>(*
tractionPtr);
1516
1517 auto next = [&]() {
1518 ++t_disp_h1;
1519 ++t_traction;
1520 ++t_coords;
1521 };
1522
1523 auto get_ele_centre = [
i](
auto t_ele_coords) {
1525 t_ele_center(
i) = 0;
1526 for (int nn = 0; nn != 3; nn++) {
1527 t_ele_center(
i) += t_ele_coords(
i);
1528 ++t_ele_coords;
1529 }
1530 t_ele_center(
i) /= 3;
1531 return t_ele_center;
1532 };
1533
1534 auto get_ele_radius = [
i](
auto t_ele_center,
auto t_ele_coords) {
1536 t_n0(
i) = t_ele_center(
i) - t_ele_coords(
i);
1538 };
1539
1540 auto get_face_conn = [this](auto face) {
1541 const EntityHandle *conn;
1542 int num_nodes;
1544 face, conn, num_nodes, true),
1545 "get conn");
1546 if (num_nodes != 3) {
1548 }
1549 return conn;
1550 };
1551
1552 auto get_face_coords = [this](auto conn) {
1553 std::array<double, 9> coords;
1555 return coords;
1556 };
1557
1558 auto get_closet_face = [
this](
auto *point_ptr,
auto r) {
1560 std::vector<EntityHandle> faces_out;
1564 "get closest faces");
1565 return faces_out;
1566 };
1567
1568 auto get_faces_out = [this](auto *point_ptr, auto *unit_ray_ptr, auto radius,
1570 std::vector<double> distances_out;
1571 std::vector<EntityHandle> faces_out;
1573
1576 point_ptr, unit_ray_ptr, &radius),
1577
1578 "get closest faces");
1579 return std::make_pair(faces_out, distances_out);
1580 };
1581
1582 auto get_normal = [](auto &ele_coords) {
1585 return t_normal;
1586 };
1587
1588 auto make_map = [&](auto &face_out, auto &face_dist, auto &t_ray_point,
1589 auto &t_unit_ray, auto &t_master_coord) {
1592 std::map<double, EntityHandle>
m;
1593 for (auto ii = 0; ii != face_out.size(); ++ii) {
1594 auto face_conn = get_face_conn(face_out[ii]);
1597 t_face_normal.normalize();
1599 t_x(
i) = t_ray_point(
i) + t_unit_ray(
i) * face_dist[ii];
1602 t_x(
i) * t_face_normal(
j) - t_master_coord(
i) * t_unit_ray(
j);
1603 if (t_unit_ray(
i) * t_face_normal(
i) > std::cos(M_PI / 3)) {
1604 auto dot = std::sqrt(t_chi(
i,
j) * t_chi(
i,
j));
1605 m[dot] = face_out[ii];
1606 }
1607 }
1609 };
1610
1611 auto get_tag_data = [this](auto tag, auto face, auto &vec) {
1613 int tag_length;
1615 vec.resize(tag_length);
1617 &*vec.begin());
1619 };
1620
1621 auto create_tag = [this](const std::string tag_name, const int size) {
1622 double def_VAL[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
1626 tag_name.c_str(), size, MB_TYPE_DOUBLE,
th,
1627 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1630 }
1632 };
1633
1634 auto set_float_precision = [](const double x) {
1635 if (std::abs(x) < std::numeric_limits<float>::epsilon())
1636 return 0.;
1637 else
1638 return x;
1639 };
1640
1641
1642 auto save_scal_tag = [&](
auto &
th,
auto v,
const int gg) {
1645 v = set_float_precision(
v);
1647 }
1649 };
1650
1651
1652 auto get_fe_adjacencies = [this](auto fe_ent) {
1655 &fe_ent, 1, 2, false, adj_faces, moab::Interface::UNION),
1656 "get adj");
1657 std::set<int> adj_ids;
1658 for (auto f : adj_faces) {
1660 }
1661 return adj_ids;
1662 };
1663
1664 auto get_face_id = [this](auto face) {
1665 int id;
1668 return id;
1669 }
1670 return -1;
1671 };
1672
1673 auto get_body_id = [this](auto face) {
1674 int id;
1677 return id;
1678 }
1679 return -1;
1680 };
1681
1682 auto get_face_part = [this](auto face) {
1683 const moab::Core *core_mesh_ptr =
1684 dynamic_cast<const moab::Core *
>(&
contactTreePtr->getPostProcMesh());
1685 auto pcomm_post_proc_mesh =
1687 int part;
1689 pcomm_post_proc_mesh->part_tag(), &face, 1, &part) == MB_SUCCESS) {
1690 return part;
1691 }
1692 return -1;
1693 };
1694
1695 auto check_face = [&](auto face, auto fe_id, auto part) {
1696 auto face_id = get_face_id(face);
1697 auto face_part = get_face_part(face);
1698 if (face_id == fe_id && face_part == part)
1699 return true;
1700 return false;
1701 };
1702
1703
1706 auto save_vec_tag = [&](
auto &
th,
auto &t_d,
const int gg) {
1710 for (
auto &
a :
v.data())
1711 a = set_float_precision(
a);
1713 &*
v.data().begin());
1714 }
1716 };
1717
1718 Tag th_mark = create_tag(
"contact_mark", 1);
1719 Tag th_mark_slave = create_tag(
"contact_mark_slave", 1);
1720 Tag th_body_id = create_tag(
"contact_body_id", 1);
1721 Tag th_gap = create_tag(
"contact_gap", 1);
1722 Tag th_tn_master = create_tag(
"contact_tn_master", 1);
1723 Tag th_tn_slave = create_tag(
"contact_tn_slave", 1);
1724 Tag th_contact_traction = create_tag(
"contact_traction", 3);
1725 Tag th_contact_traction_master = create_tag(
"contact_traction_master", 3);
1726 Tag th_contact_traction_slave = create_tag(
"contact_traction_slave", 3);
1727 Tag th_c = create_tag(
"contact_c", 1);
1728 Tag th_normal = create_tag(
"contact_normal", 3);
1729 Tag th_dist = create_tag(
"contact_dip", 3);
1730
1731 auto t_ele_centre = get_ele_centre(getFTensor1Coords());
1732 auto ele_radius = get_ele_radius(t_ele_centre, getFTensor1Coords());
1733
1736 shadow_vec.clear();
1737
1738 auto adj_fe_ids = get_fe_adjacencies(fe_ent);
1739
1740 for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1741
1743 t_spatial_coords(
i) = t_coords(
i) + t_disp_h1(
i);
1744
1746 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1747 }
1748
1749 auto faces_close = get_closet_face(&t_spatial_coords(0), ele_radius);
1750 for (auto face_close : faces_close) {
1751 if (check_face(face_close, fe_id, m_field.get_comm_rank())) {
1752
1753 auto body_id = get_body_id(face_close);
1754
1755 auto master_face_conn = get_face_conn(face_close);
1756 std::array<double, 9> master_coords;
1759 master_coords.data());
1760 std::array<double, 9> master_traction;
1763 master_traction.data());
1764 auto t_normal_face_close = get_normal(master_coords);
1765 t_normal_face_close.normalize();
1766
1769 CHKERR save_scal_tag(th_mark,
m, gg);
1770 CHKERR save_scal_tag(th_body_id,
static_cast<double>(body_id), gg);
1771 CHKERR save_vec_tag(th_normal, t_normal_face_close, gg);
1772 CHKERR save_vec_tag(th_contact_traction, t_traction, gg);
1773 }
1774
1776 t_unit_ray(
i) = -t_normal_face_close(
i);
1779 t_spatial_coords(
i) -
1781
1782 constexpr double eps = 1e-3;
1783 auto [faces_out, faces_dist] =
1784 get_faces_out(&t_ray_point(0), &t_unit_ray(0),
1787
1788 auto m = make_map(faces_out, faces_dist, t_ray_point, t_unit_ray,
1789 t_spatial_coords);
1790 for (
auto m_it =
m.begin(); m_it !=
m.end(); ++m_it) {
1791 auto face = m_it->second;
1792 if (face != face_close) {
1793
1794 if (
1795
1796 (adj_fe_ids.find(get_face_id(face)) == adj_fe_ids.end() ||
1797 get_face_part(face) != m_field.get_comm_rank())
1798
1799 ) {
1800
1801 shadow_vec.push_back(ContactTree::FaceData());
1802 shadow_vec.back().gaussPtNb = gg;
1803
1804 auto slave_face_conn = get_face_conn(face);
1805 std::array<double, 9> slave_coords;
1808 slave_coords.data());
1809 auto t_normal_face = get_normal(slave_coords);
1810 std::array<double, 9> slave_tractions;
1813 slave_tractions.data());
1814
1815 auto t_master_point =
1816 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data());
1817 auto t_slave_point =
1818 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data());
1819 auto t_ray_point_data =
1820 getFTensor1FromPtr<3>(shadow_vec.back().rayPoint.data());
1821 auto t_unit_ray_data =
1822 getFTensor1FromPtr<3>(shadow_vec.back().unitRay.data());
1823
1824 t_slave_point(
i) = t_ray_point(
i) + m_it->first * t_unit_ray(
i);
1825
1826 auto eval_position = [&](auto &&t_elem_coords, auto &&t_point) {
1827 std::array<double, 2> loc_coords;
1830 &t_elem_coords(0, 0), &t_point(0), 1,
1831 loc_coords.data()),
1832 "get local coords");
1835 &loc_coords[0],
1836 &loc_coords[1], 1),
1837 "calc shape fun");
1841 t_point_out(
i) = t_shape_fun(
j) * t_elem_coords(
j,
i);
1842 return t_point_out;
1843 };
1844
1845 auto t_master_point_updated = eval_position(
1846 getFTensor2FromPtr<3, 3>(master_coords.data()),
1847 getFTensor1FromPtr<3>(shadow_vec.back().masterPoint.data()));
1848 t_master_point(
i) = t_master_point_updated(
i);
1849
1850 auto t_slave_point_updated = eval_position(
1851 getFTensor2FromPtr<3, 3>(slave_coords.data()),
1852 getFTensor1FromPtr<3>(shadow_vec.back().slavePoint.data()));
1853 t_slave_point(
i) = t_slave_point_updated(
i);
1854
1855 t_ray_point_data(
i) = t_ray_point(
i);
1856 t_unit_ray_data(
i) = t_unit_ray(
i);
1857
1858 std::copy(master_coords.begin(), master_coords.end(),
1859 shadow_vec.back().masterPointNodes.begin());
1860 std::copy(master_traction.begin(), master_traction.end(),
1861 shadow_vec.back().masterTractionNodes.begin());
1862 std::copy(slave_coords.begin(), slave_coords.end(),
1863 shadow_vec.back().slavePointNodes.begin());
1864 std::copy(slave_tractions.begin(), slave_tractions.end(),
1865 shadow_vec.back().slaveTractionNodes.begin());
1866
1867 shadow_vec.back().eleRadius = ele_radius;
1868
1869
1870
1871
1872
1873
1874
1875
1877 auto [gap, tn_master, tn_slave,
c, t_master_traction,
1878 t_slave_traction] =
1879 multiGetGap(&(shadow_vec.back()), t_spatial_coords);
1881 t_gap_vec(
i) = t_slave_point(
i) - t_spatial_coords(
i);
1882 CHKERR save_scal_tag(th_gap, gap, gg);
1883 CHKERR save_scal_tag(th_tn_master, tn_master, gg);
1884 CHKERR save_scal_tag(th_tn_slave, tn_slave, gg);
1885 CHKERR save_scal_tag(th_c,
c, gg);
1887 CHKERR save_scal_tag(th_mark_slave,
m, gg);
1888 CHKERR save_vec_tag(th_dist, t_gap_vec, gg);
1889 CHKERR save_vec_tag(th_contact_traction_master,
1890 t_master_traction, gg);
1891 CHKERR save_vec_tag(th_contact_traction_slave, t_slave_traction,
1892 gg);
1893 }
1894
1895 break;
1896 }
1897 }
1898 }
1899 break;
1900 }
1901 }
1902 next();
1903 }
1904
1906}
#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 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
#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
auto multiGetGap(ContactTree::FaceData *face_data_ptr, FTensor::Tensor1< T1, 3 > &t_spatial_coords)
VectorBoundedArray< double, 3 > VectorDouble3
auto id_from_handle(const EntityHandle h)
FTensor::Index< 'm', 3 > m
static auto getOrCreate(moab::Core *core_mesh_ptr)
