992 Range *fixed_edges,
Range *corner_nodes,
const double geometry_tol,
993 const double close_tol,
const int verb,
const bool debug) {
995 moab::Interface &moab = m_field.
get_moab();
999 auto get_ent_adj = [&moab](
const EntityHandle v,
const int dim) {
1002 CHKERR moab.get_adjacencies(&
v, 1, dim,
true,
a);
1006 auto get_adj = [&moab](
const Range r,
const int dim) {
1009 CHKERR moab.get_adjacencies(r, dim,
false,
a, moab::Interface::UNION);
1011 CHKERR moab.get_connectivity(r,
a,
true);
1015 auto get_skin = [&skin](
const auto r) {
1017 CHKERR skin.find_skin(0, r,
false, s);
1022 auto get_range = [](std::vector<EntityHandle>
v) {
1024 r.insert_list(
v.begin(),
v.end());
1028 auto get_coords = [&](
const auto v) {
1030 CHKERR moab.get_coords(&
v, 1, &coords[0]);
1035 CHKERR moab.tag_get_handle(
"DIST_SURFACE_NORMAL_VECTOR", th_normal_dist);
1039 CHKERR moab.tag_get_data(th_normal_dist, &
v, 1, &*dist_vec.begin());
1044 const double geometry_tol) {
1045 auto get_tuple = [](
const EntityHandle e,
const double dist,
1046 const double l) {
return std::make_tuple(e, dist,
l); };
1048 std::tuple<EntityHandle, double, double> min_tuple{0, 0, 0};
1050 for (
auto e : edges) {
1054 auto get_dist = [&](
auto eit) {
1057 CHKERR moab.get_connectivity(eit->first, conn, num_nodes,
true);
1059 return eit->second.dIst;
1060 else if (conn[1] ==
v)
1061 return (eit->second.lEngth - eit->second.dIst);
1066 const auto d = get_dist(eit);
1067 if (std::get<0>(min_tuple) == 0) {
1068 min_tuple = get_tuple(e, d, eit->second.lEngth);
1069 }
else if (std::get<1>(min_tuple) > d) {
1070 min_tuple = get_tuple(e, d, eit->second.lEngth);
1075 const auto geom_dist_vec = get_normal_dits(
v);
1076 const auto geom_tol = norm_2(geom_dist_vec);
1077 const auto l = std::get<2>(min_tuple);
1079 if (geom_tol <
l * geometry_tol) {
1080 return std::make_pair(get_coords(
v),
l);
1083 const auto &d =
edgesToCut.at(std::get<0>(min_tuple));
1084 return std::make_pair(
VectorDouble3(d.rayPoint + d.dIst * d.unitRayDir),
1089 auto get_in_range = [](
auto v,
auto &r) {
return (r.find(
v) != r.end()); };
1091 auto project_nodes = [&](
auto nodes_to_check) {
1095 auto get_fix_e = [](
auto fixed_edges) {
1097 return *fixed_edges;
1102 const auto fix_e = get_fix_e(fixed_edges);
1104 const auto cut_fix = intersect(
cutEdges, fix_e);
1105 const auto cut_skin = intersect(
cutEdges, skin_e);
1109 corner_n = intersect(*corner_nodes, nodes_to_check);
1110 auto fixe_n = intersect(get_adj(fix_e, 0), nodes_to_check);
1111 auto skin_n = intersect(get_adj(skin_e, 0), nodes_to_check);
1113 std::vector<std::pair<EntityHandle, TreeData>> vertices_on_cut_edges;
1114 vertices_on_cut_edges.reserve(nodes_to_check.size());
1116 auto add_zero_vertex = [&](
const EntityHandle v,
auto &&new_pos) {
1117 auto coords = get_coords(
v);
1118 auto ray = new_pos - coords;
1119 auto l = norm_2(ray);
1122 if (
l > std::numeric_limits<double>::epsilon()) {
1129 return std::make_pair(
v, data);
1132 auto intersect_v = [&](
const auto v,
const auto r) {
1133 return intersect(r, get_ent_adj(
v, 1));
1136 for (
auto v : nodes_to_check) {
1138 const auto e = intersect_v(
v,
cutEdges);
1141 if (get_in_range(
v, corner_n)) {
1142 auto p = get_prj_point(
v, e, 0);
1143 if (norm_2(get_coords(
v) - p.first) < close_tol * p.second) {
1144 vertices_on_cut_edges.push_back(add_zero_vertex(
v, get_coords(
v)));
1148 }
else if (get_in_range(
v, fixe_n)) {
1149 const auto b = intersect_v(
v, cut_fix);
1151 auto p = get_prj_point(
v, b, geometry_tol);
1152 if (norm_2(get_coords(
v) - p.first) < close_tol * p.second) {
1153 vertices_on_cut_edges.push_back(add_zero_vertex(
v, p.first));
1157 }
else if (get_in_range(
v, skin_n)) {
1158 const auto b = intersect_v(
v, cut_skin);
1160 auto p = get_prj_point(
v, b, geometry_tol);
1161 if (norm_2(get_coords(
v) - p.first) < close_tol * p.second) {
1162 vertices_on_cut_edges.push_back(add_zero_vertex(
v, p.first));
1170 auto p = get_prj_point(
v, e, geometry_tol);
1171 if (norm_2(get_coords(
v) - p.first) < close_tol * p.second) {
1172 if (get_in_range(
v, fixe_n) || get_in_range(
v, skin_n))
1173 vertices_on_cut_edges.push_back(add_zero_vertex(
v, get_coords(
v)));
1175 vertices_on_cut_edges.push_back(add_zero_vertex(
v, p.first));
1181 auto get_distances = [&](
auto &data) {
1182 std::map<EntityHandle, double> dist_map;
1183 if (!data.empty()) {
1185 CHKERR moab.tag_get_handle(
"DIST_SURFACE_NORMAL_SIGNED",
th);
1187 std::vector<EntityHandle> verts;
1188 verts.reserve(data.size());
1190 verts.emplace_back(d.first);
1191 std::vector<double> distances(verts.size());
1192 CHKERR moab.tag_get_data(
th, &*verts.begin(), verts.size(),
1193 &*distances.begin());
1194 for (
size_t i = 0;
i != verts.size(); ++
i)
1195 dist_map[verts[
i]] = distances[
i];
1200 auto dist_map = get_distances(vertices_on_cut_edges);
1201 if (!dist_map.empty()) {
1202 auto cmp = [&dist_map](
const auto &
a,
const auto &b) {
1203 return dist_map[
a.first] < dist_map[b.first];
1205 std::sort(vertices_on_cut_edges.begin(), vertices_on_cut_edges.end(),
1209 return vertices_on_cut_edges;
1212 auto get_min_quality =
1213 [&](
const Range &adj_tets,
1214 const map<EntityHandle, TreeData> &vertices_on_cut_edges,
1215 const std::pair<EntityHandle, TreeData> *test_data_ptr) {
1217 for (
auto t : adj_tets) {
1222 CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1224 auto set_new_coord = [](
auto d) {
1225 return d.rayPoint + d.dIst * d.unitRayDir;
1228 for (
auto n : {0, 1, 2, 3}) {
1230 if (test_data_ptr) {
1231 if (test_data_ptr->first == conn[
n]) {
1232 noalias(n_coords) = set_new_coord(test_data_ptr->second);
1235 auto mit = vertices_on_cut_edges.find(conn[
n]);
1236 if (mit != vertices_on_cut_edges.end()) {
1237 noalias(n_coords) = set_new_coord(mit->second);
1245 auto get_zero_distance_verts = [&](
const auto &&vertices_on_cut_edges) {
1246 std::vector<EntityHandle> zero_dist_vec;
1247 zero_dist_vec.reserve(vertices_on_cut_edges.size());
1248 for (
auto t : vertices_on_cut_edges) {
1250 auto adj_tet = intersect(
vOlume, get_ent_adj(
t.first, 3));
1255 zero_dist_vec.push_back(
t.first);
1258 return zero_dist_vec;
1263 auto get_zero_distant_ents = [&](
auto bridge_ents,
const int dim,
1264 const int bridge_dim) {
1266 intersect(get_adj(bridge_ents, dim), get_adj(vol_cut_ents, dim));
1267 auto ents_to_remove =
1268 get_adj(subtract(get_adj(ents, bridge_dim), bridge_ents), dim);
1269 return subtract(ents, ents_to_remove);
1275 get_range(get_zero_distance_verts(project_nodes(get_adj(
cutEdges, 0)))));
1285 auto adj_edges = get_ent_adj(f, 1);
1286 for (
auto e : adj_edges) {
1482 const double tol_trim_close,
1485 moab::Interface &moab = m_field.
get_moab();
1489 Skinner skin(&moab);
1495 Range tets_skin_verts;
1496 CHKERR moab.get_connectivity(tets_skin, tets_skin_verts,
true);
1498 Range tets_skin_edges;
1499 CHKERR moab.get_adjacencies(tets_skin, 1,
false, tets_skin_edges,
1500 moab::Interface::UNION);
1503 Range cut_surface_edges;
1505 moab::Interface::UNION);
1506 Range cut_surface_verts;
1511 corners = *corner_nodes;
1515 fix_edges = *fixed_edges;
1517 Range fixed_edges_verts;
1519 CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts,
true);
1525 surface_skin =
fRont;
1530 CHKERR moab.tag_get_data(
th, &
v, 1, &point[0]);
1532 CHKERR moab.get_coords(&
v, 1, &point[0]);
1541 : d(d),
v(
v), e(e) {}
1544 typedef multi_index_container<
1547 ordered_non_unique<member<VertMap, const double, &VertMap::d>>,
1548 ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::v>>,
1549 ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::e>>
1554 VertMapMultiIndex verts_map;
1557 const double dist) {
1558 verts_map.insert(VertMap(dist,
v, e));
1567 auto cut_this_edge = [&](
const EntityHandle e,
const double length,
auto &ray,
1574 for (
int ii = 0; ii != 3; ++ii)
1576 for (
int ii = 0; ii != 3; ++ii)
1586 auto get_tag_data = [&](
auto th,
auto conn) {
1588 CHKERR moab.tag_get_data(
th, &conn, 1, &
t(0));
1592 double max_edge_length = 0;
1595 if (!
fRont.empty()) {
1597 treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
1598 new OrientedBoxTreeTool(&moab,
"ROOTSETSURF",
true));
1601 for (
auto s : surface_skin) {
1603 auto project_on_surface = [&](
auto &
t) {
1613 t_normal(
i) /= sqrt(t_normal(
i) * t_normal(
i));
1614 const double dot = t_normal(
i) * (t_p(
i) -
t(
i));
1615 t_p(
i) =
t(
i) + dot * t_normal(
i);
1621 CHKERR moab.get_connectivity(s, conn, num_nodes,
true);
1624 CHKERR moab.get_coords(conn, num_nodes, &coords_front[0]);
1631 auto t_p_f0 = project_on_surface(t_f0);
1632 auto t_p_f1 = project_on_surface(t_f1);
1634 CHKERR moab.set_coords(&conn[0], 1, &t_p_f0(0));
1635 CHKERR moab.set_coords(&conn[1], 1, &t_p_f1(0));
1643 Tag th_dist_front_vec;
1644 CHKERR moab.tag_get_handle(
"DIST_FRONT_VECTOR", th_dist_front_vec);
1647 for (
auto e : cut_surface_edges) {
1649 auto get_conn_edge = [&moab](
auto e) {
1653 CHKERR moab.get_connectivity(e, conn_edge, num_nodes,
true);
1657 auto get_coords_edge = [&moab](
auto conn_edge) {
1658 std::array<double, 6> coords_edge;
1659 CHKERR moab.get_coords(conn_edge, 2, coords_edge.data());
1663 auto get_ftensor_coords = [](
const double *ptr) {
1667 auto conn_edge = get_conn_edge(e);
1668 auto coords_edge = get_coords_edge(conn_edge);
1670 auto t_dist0 = get_tag_data(th_dist_front_vec, conn_edge[0]);
1671 auto t_dist1 = get_tag_data(th_dist_front_vec, conn_edge[1]);
1672 auto t_e0 = get_ftensor_coords(&coords_edge[0]);
1673 auto t_e1 = get_ftensor_coords(&coords_edge[3]);
1676 t_edge_delta(
i) = t_e1(
i) - t_e0(
i);
1677 const double edge_length2 = t_edge_delta(
i) * t_edge_delta(
i);
1678 const double edge_length = sqrt(edge_length2);
1679 if (edge_length == 0)
1682 auto add_edge = [&](
auto t_cut) {
1685 t_ray(
i) = t_cut * t_edge_delta(
i);
1686 add_vert(conn_edge[0], e, fabs(t_cut));
1687 add_vert(conn_edge[1], e, fabs(t_cut - 1));
1688 cut_this_edge(e, edge_length, t_ray, t_e0);
1691 t_edge_point(
i) = t_e0(
i) + t_cut * t_edge_delta(
i);
1692 t_ray(
i) = t_edge_point(
i) - t_e1(
i);
1693 add_vert(conn_edge[0], e, fabs(t_cut));
1694 add_vert(conn_edge[1], e, fabs(t_cut - 1));
1695 cut_this_edge(e, edge_length, t_ray, t_e1);
1697 max_edge_length = std::max(max_edge_length, edge_length);
1700 auto trim_cross_edges = [&]() {
1701 if (std::min(t_dist0(
i) * t_dist0(
i), t_dist1(
i) * t_dist1(
i)) <
1704 auto make_pair = [&](
const double d,
const double te) {
1705 return std::make_pair(d, te);
1708 auto min_pair = make_pair(-1, -1);
1710 for (
auto f : surface_skin) {
1712 auto conn_front = get_conn_edge(f);
1713 auto coords_front = get_coords_edge(conn_front);
1714 auto t_f0 = get_ftensor_coords(&coords_front[0]);
1715 auto t_f1 = get_ftensor_coords(&coords_front[3]);
1721 &t_e0(0), &t_e1(0), &t_f0(0), &t_f1(0), &te, &tf);
1725 auto check = [](
auto v) {
1726 return (
v > -std::numeric_limits<double>::epsilon() &&
1727 (
v - 1) < std::numeric_limits<double>::epsilon());
1730 if (check(te) && check(tf)) {
1733 t_delta(
i) = (t_e0(
i) + te * t_edge_delta(
i)) -
1734 (t_f0(0) + tf * (t_f1(
i) - t_f0(
i)));
1736 t_delta(
i) /= edge_length;
1738 const double dot = t_delta(
i) * t_delta(
i);
1740 if (min_pair.first < 0 || min_pair.first > dot) {
1742 if (dot < tol_trim_close * tol_trim_close)
1743 min_pair = make_pair(dot, te);
1749 if (min_pair.first > -std::numeric_limits<double>::epsilon()) {
1750 add_edge(min_pair.second);
1758 if (!trim_cross_edges()) {
1760 const double dot = t_dist0(
i) * t_dist1(
i);
1761 const double dot_direction0 = t_dist0(
i) * t_edge_delta(
i);
1762 const double dot_direction1 = t_dist1(
i) * t_edge_delta(
i);
1764 if (dot < std::numeric_limits<double>::epsilon() &&
1765 dot_direction0 > -std::numeric_limits<double>::epsilon() &&
1766 dot_direction1 < std::numeric_limits<double>::epsilon()) {
1768 const double s0 = t_dist0(
i) * t_edge_delta(
i) / edge_length;
1769 const double s1 = t_dist1(
i) * t_edge_delta(
i) / edge_length;
1770 const double t_cut = s0 / (s0 - s1);
1778 CHKERR SaveData(moab)(
"edges_potentially_to_trim.vtk", cut_surface_edges);
1789 for (
auto t : adj_tets) {
1795 CHKERR moab.tag_get_data(
th, conn, num_nodes, &coords[0]);
1797 CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1799 for (
int n = 0;
n != 4; ++
n) {
1802 noalias(n_coords) = new_pos;
1807 m->second.rayPoint +
m->second.dIst *
m->second.unitRayDir;
1812 if (!std::isnormal(q)) {
1816 min_q = std::min(min_q, q);
1821 Range checked_verts;
1829 CHKERR moab.get_adjacencies(&
v, 1, 3,
false, adj_tets);
1831 double q = get_quality_change(adj_tets,
v, new_pos);
1834 double dist = norm_2(ray_dir);
1841 checked_verts.insert(
v);
1853 checked_verts.insert(
v);
1861 auto hi = verts_map.get<0>().upper_bound(tol_trim_close);
1862 verts_map.get<0>().erase(hi, verts_map.end());
1864 auto remove_verts = [&](
Range nodes) {
1865 for (
auto n : nodes) {
1866 auto r = verts_map.get<1>().equal_range(
n);
1867 verts_map.get<1>().erase(r.first, r.second);
1871 auto trim_verts = [&](
const Range verts,
const bool move) {
1872 VertMapMultiIndex verts_map_tmp;
1873 for (
auto p = corners.pair_begin(); p != corners.pair_end(); ++p) {
1874 auto lo = verts_map.get<1>().lower_bound(p->first);
1875 auto hi = verts_map.get<1>().upper_bound(p->second);
1876 verts_map_tmp.insert(lo, hi);
1879 for (
auto &
m : verts_map_tmp.get<0>())
1880 add_trim_vert(
m.v,
m.e);
1882 for (
auto &
m : verts_map_tmp.get<0>())
1883 add_no_move_trim(
m.v,
m.e);
1887 auto trim_edges = [&](
const Range ents,
const bool move) {
1888 VertMapMultiIndex verts_map_tmp;
1889 for (
auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
1890 auto lo = verts_map.get<2>().lower_bound(p->first);
1891 auto hi = verts_map.get<2>().upper_bound(p->second);
1892 verts_map_tmp.insert(lo, hi);
1893 verts_map.get<2>().erase(lo, hi);
1896 for (
auto &
m : verts_map_tmp.get<0>())
1897 add_trim_vert(
m.v,
m.e);
1899 for (
auto &
m : verts_map_tmp.get<0>())
1900 add_no_move_trim(
m.v,
m.e);
1904 auto intersect_self = [&](
Range &
a,
const Range b) {
a = intersect(
a, b); };
1906 map<std::string, Range> range_maps;
1911 intersect_self(range_maps[
"surface_skin"],
trimEdges);
1913 range_maps[
"fixed_edges_on_surface_skin"] =
1914 intersect(range_maps[
"surface_skin"], fix_edges);
1917 CHKERR moab.get_adjacencies(fixed_edges_verts, 1,
false,
1918 range_maps[
"fixed_edges_verts_edges"],
1919 moab::Interface::UNION);
1920 intersect_self(range_maps[
"fixed_edges_verts_edges"],
trimEdges);
1922 CHKERR moab.get_connectivity(
1923 range_maps[
"fixed_edges_verts_edges"],
1924 range_maps[
"fixed_edges_verts_edges_verts_on_the_skin"],
false);
1925 intersect_self(range_maps[
"fixed_edges_verts_edges_verts_on_the_skin"],
1929 trim_verts(corners,
false);
1930 remove_verts(corners);
1933 trim_edges(range_maps[
"fixed_edges_on_surface_skin"],
true);
1934 remove_verts(range_maps[
"fixed_edges_on_surface_skin_verts"]);
1937 trim_verts(range_maps[
"fixed_edges_verts_edges_verts_on_the_skin"],
false);
1938 remove_verts(range_maps[
"fixed_edges_verts_edges_verts_on_the_skin"]);
1941 trim_edges(range_maps[
"surface_skin"],
true);
1942 trim_verts(tets_skin_verts,
false);
1943 remove_verts(tets_skin_verts);
1945 for (
auto &
m : verts_map.get<0>())
1946 add_trim_vert(
m.v,
m.e);
1951 CHKERR moab.get_adjacencies(&
v, 1, 1,
false, adj);
1952 for (
auto e : adj) {
2029 Range *corner_nodes,
2033 moab::Interface &moab = m_field.
get_moab();
2034 Skinner skin(&moab);
2036 auto get_adj = [&moab](
const Range r,
const int dim) {
2039 CHKERR moab.get_adjacencies(r, dim,
false,
a, moab::Interface::UNION);
2041 CHKERR moab.get_connectivity(r,
a,
true);
2045 auto get_skin = [&skin](
const auto r) {
2047 CHKERR skin.find_skin(0, r,
false, s);
2054 auto trim_tets_skin_edges = get_adj(trim_tets_skin, 1);
2057 auto contarain_edges =
2060 contarain_edges.merge(
2061 intersect(fixed_edges->subset_by_type(MBEDGE), trim_surface_edges));
2066 barrier_vertices.merge(get_adj(
2067 intersect(get_adj(intersect(*corner_nodes, barrier_vertices), 2),
2072 auto get_nodes_with_one_node_on_fixed_edge_other_not = [&]() {
2073 const auto fixed_edges_on_surface =
2074 intersect(*fixed_edges, trim_surface_edges);
2075 const auto skin_fixed_edges_on_surface = get_skin(fixed_edges_on_surface);
2076 const auto barrier_nodes = subtract(skin_fixed_edges_on_surface,
2078 return barrier_nodes;
2082 barrier_vertices.merge(get_nodes_with_one_node_on_fixed_edge_other_not());
2084 auto add_close_surface_barrier = [&]() {
2096 test_edges = subtract(test_edges, *fixed_edges);
2097 auto trim_surface_nodes = get_adj(test_edges, 0);
2098 trim_surface_nodes = subtract(trim_surface_nodes, barrier_vertices);
2100 trim_surface_nodes =
2101 subtract(trim_surface_nodes, get_adj(*fixed_edges, 0));
2104 trim_skin.merge(contarain_edges);
2106 trim_skin.merge(get_skin(
sUrface));
2108 trim_skin.merge(
fRont);
2110 trim_skin.merge(intersect(*fixed_edges, trim_surface_edges));
2112 if (
debug && !trim_skin.empty())
2117 CHKERR moab.tag_get_data(
th, trim_surface_nodes, &*coords.begin());
2119 CHKERR moab.get_coords(trim_surface_nodes, &*coords.begin());
2121 if (!trim_skin.empty()) {
2123 std::vector<double> min_distances(trim_surface_nodes.size(), -1);
2125 &*coords.begin(), trim_surface_nodes.size(), trim_skin,
2126 &*min_distances.begin());
2128 auto coords_ptr = &*coords.begin();
2129 auto min_dist = &*min_distances.begin();
2131 std::vector<EntityHandle> add_nodes;
2132 add_nodes.reserve(trim_surface_nodes.size());
2134 for (
auto n : trim_surface_nodes) {
2136 if ((*min_dist) > std::numeric_limits<double>::epsilon()) {
2141 &point_out[0], facets_out);
2148 normal /= norm_2(normal);
2151 double dist = norm_2(
delta);
2152 if (dist <
tol * (*min_dist))
2153 add_nodes.emplace_back(
n);
2160 barrier_vertices.insert_list(add_nodes.begin(), add_nodes.end());
2166 auto remove_faces_on_skin = [&]() {
2169 if (!skin_faces.empty()) {
2170 barrier_vertices.merge(get_adj(skin_faces, 0));
2171 for (
auto f : skin_faces)
2177 auto get_trim_free_edges = [&]() {
2179 Range trim_surf_skin;
2181 trim_surf_skin = subtract(trim_surf_skin, trim_tets_skin_edges);
2183 Range trim_surf_skin_verts;
2184 CHKERR moab.get_connectivity(trim_surf_skin, trim_surf_skin_verts,
true);
2185 for (
auto e : barrier_vertices)
2186 trim_surf_skin_verts.erase(e);
2189 CHKERR moab.get_adjacencies(trim_surf_skin_verts, 1,
false, free_edges,
2190 moab::Interface::UNION);
2192 subtract(intersect(free_edges, trim_surf_skin), contarain_edges);
2197 CHKERR remove_faces_on_skin();
2198 CHKERR add_close_surface_barrier();
2200 if (
debug && !barrier_vertices.empty())
2206 while (!(out_edges = get_trim_free_edges()).empty()) {
2210 "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) +
".vtk",
2213 if (
debug && !out_edges.empty())
2215 "trimNewSurfacesOutsideVerts_" +
2216 boost::lexical_cast<std::string>(nn) +
".vtk",
2219 Range outside_faces;
2220 CHKERR moab.get_adjacencies(out_edges, 2,
false, outside_faces,
2221 moab::Interface::UNION);
2228 "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) +
".vtk",
2304 const Range &fixed_edges,
const Range &corner_nodes,
Range &edges_to_merge,
2308 moab::Interface &moab = m_field.
get_moab();
2316 const bool onlyIfImproveQuality;
2320 MergeNodes(
CoreInterface &m_field,
const bool only_if_improve_quality,
2321 Tag th,
bool update_mehsets)
2322 : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
2323 tH(
th), updateMehsets(update_mehsets) {
2329 bool add_child =
true) {
2330 moab::Interface &moab = mField.
get_moab();
2334 CHKERR moab.get_adjacencies(conn, 2, 3,
false, vert_tets,
2335 moab::Interface::UNION);
2336 vert_tets = intersect(vert_tets, proc_tets);
2339 onlyIfImproveQuality, 0, line_search,
2344 Range::iterator lo, hi = proc_tets.begin();
2345 for (
auto pt = vert_tets.pair_begin(); pt != vert_tets.pair_end();
2347 lo = proc_tets.lower_bound(hi, proc_tets.end(), pt->first);
2348 if (lo != proc_tets.end()) {
2349 hi = proc_tets.upper_bound(lo, proc_tets.end(), pt->second);
2350 proc_tets.erase(lo, hi);
2354 proc_tets.merge(out_tets);
2358 struct ChangeChild {
2360 ChangeChild(
const EntityHandle child) : child(child) {}
2366 std::vector<
decltype(parentsChildMap.get<0>().begin())> it_vec;
2367 it_vec.reserve(parentsChildMap.size());
2369 for (
auto &p : parent_child_map) {
2372 for (
auto it = parentsChildMap.get<0>().equal_range(p.pArent);
2373 it.first != it.second; ++it.first)
2374 it_vec.emplace_back(it.first);
2376 for (
auto it = parentsChildMap.get<1>().equal_range(p.pArent);
2377 it.first != it.second; ++it.first)
2378 it_vec.emplace_back(parentsChildMap.project<0>(it.first));
2380 for (
auto &it : it_vec)
2381 parentsChildMap.modify(it, ChangeChild(p.cHild));
2384 parentsChildMap.insert(parent_child_map.begin(),
2385 parent_child_map.end());
2391 Range ¬_merged_edges,
2393 moab::Interface &moab = mField.
get_moab();
2397 std::vector<EntityHandle> parents_ents_vec(parentsChildMap.size());
2398 for (
auto &it : parentsChildMap)
2399 parents_ents_vec.emplace_back(it.pArent);
2401 parent_ents.insert_list(parents_ents_vec.begin(),
2402 parents_ents_vec.end());
2404 Range surf_parent_ents = intersect(new_surf, parent_ents);
2405 new_surf = subtract(new_surf, surf_parent_ents);
2406 Range child_surf_ents;
2407 CHKERR updateRangeByChilds(parentsChildMap, surf_parent_ents,
2409 new_surf.merge(child_surf_ents);
2411 Range edges_to_merge_parent_ents =
2412 intersect(edges_to_merge, parent_ents);
2413 edges_to_merge = subtract(edges_to_merge, edges_to_merge_parent_ents);
2414 Range merged_child_edge_ents;
2415 CHKERR updateRangeByChilds(parentsChildMap, edges_to_merge_parent_ents,
2416 merged_child_edge_ents);
2418 Range not_merged_edges_child_ents =
2419 intersect(not_merged_edges, parent_ents);
2421 subtract(not_merged_edges, not_merged_edges_child_ents);
2422 Range not_merged_child_edge_ents;
2423 CHKERR updateRangeByChilds(parentsChildMap, not_merged_edges_child_ents,
2424 not_merged_child_edge_ents);
2426 merged_child_edge_ents =
2427 subtract(merged_child_edge_ents, not_merged_child_edge_ents);
2428 edges_to_merge.merge(merged_child_edge_ents);
2429 not_merged_edges.merge(not_merged_child_edge_ents);
2431 if (updateMehsets) {
2436 CHKERR moab.get_entities_by_handle(cubit_meshset, meshset_ents,
2439 CHKERR updateRangeByChilds(parentsChildMap, meshset_ents,
2441 CHKERR moab.add_entities(cubit_meshset, child_ents);
2451 std::vector<EntityHandle> childsVec;
2458 childsVec.reserve(parents.size());
2459 for (
auto pit = parents.pair_begin(); pit != parents.pair_end(); pit++) {
2460 auto it = parent_child_map.lower_bound(pit->first);
2461 if (it != parent_child_map.end()) {
2462 for (
auto hi_it = parent_child_map.upper_bound(pit->second);
2464 childsVec.emplace_back(it->cHild);
2467 childs.insert_list(childsVec.begin(), childsVec.end());
2478 moab::Interface &moab;
2479 const double maxLength;
2481 : tH(
th), mField(m_field), moab(m_field.
get_moab()),
2482 maxLength(max_length) {
2485 acceptedThrasholdMergeQuality = 0.5;
2492 double acceptedThrasholdMergeQuality;
2497 double &ave)
const {
2500 std::array<double, 12> coords;
2502 VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
2503 VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
2506 struct NodeQuality {
2509 NodeQuality(
const EntityHandle node) : node(node), quality(1) {}
2512 typedef multi_index_container<
2513 NodeQuality, indexed_by<
2517 hashed_non_unique<tag<Ent_mi_tag>,
2519 &NodeQuality::node>>
2522 NodeQuality_sequence;
2524 NodeQuality_sequence node_quality_sequence;
2527 CHKERR moab.get_connectivity(tets, edges_nodes,
false);
2529 CHKERR moab.get_adjacencies(edges, 3,
false, edges_tets,
2530 moab::Interface::UNION);
2531 edges_tets = intersect(edges_tets, tets);
2533 for (
auto node : edges_nodes)
2534 node_quality_sequence.emplace_back(node);
2536 for (
auto tet : edges_tets) {
2538 CHKERR moab.get_connectivity(tet, conn, num_nodes,
true);
2540 CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2542 CHKERR moab.get_coords(conn, num_nodes, coords.data());
2544 const double q = Tools::volumeLengthQuality(coords.data());
2546 for (
auto n : {0, 1, 2, 3}) {
2547 auto it = node_quality_sequence.get<1>().find(conn[
n]);
2548 if (it != node_quality_sequence.get<1>().end())
2549 const_cast<double &
>(it->quality) = std::min(q, it->quality);
2553 for (
auto edge : edges) {
2555 CHKERR moab.get_connectivity(edge, conn, num_nodes,
true);
2558 CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2560 CHKERR moab.get_coords(conn, num_nodes, coords.data());
2563 for (
auto n : {0, 1}) {
2564 auto it = node_quality_sequence.get<1>().find(conn[
n]);
2565 if (it != node_quality_sequence.get<1>().end())
2566 q = std::min(q, it->quality);
2569 if (q < acceptedThrasholdMergeQuality) {
2570 noalias(
delta) = (s0 - s1) / maxLength;
2572 double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
2578 for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
2579 length_map.get<0>().begin();
2580 mit != length_map.get<0>().end(); mit++) {
2581 ave += mit->qUality;
2583 ave /= length_map.size();
2601 SURFACE_SKIN = 1 << 2,
2602 FRONT_ENDS = 1 << 3,
2604 FIX_CORNERS = 1 << 5
2607 typedef map<int, Range> SetsMap;
2610 const Range &fixed_edges,
2611 const Range &corner_nodes,
2612 const Range &constrain_surface,
2613 SetsMap &sets_map)
const {
2614 moab::Interface &moab(mField.
get_moab());
2615 Skinner skin(&moab);
2618 sets_map[FIX_CORNERS].merge(corner_nodes);
2620 CHKERR moab.get_connectivity(fixed_edges, fixed_verts,
true);
2621 sets_map[FIX_EDGES].swap(fixed_verts);
2624 CHKERR skin.find_skin(0, tets,
false, tets_skin);
2625 Range tets_skin_edges;
2626 CHKERR moab.get_adjacencies(tets_skin, 1,
false, tets_skin_edges,
2627 moab::Interface::UNION);
2630 Range constrain_surface_verts;
2631 CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
2633 Range constrain_surface_edges;
2634 CHKERR moab.get_adjacencies(constrain_surface, 1,
false,
2635 constrain_surface_edges,
2636 moab::Interface::UNION);
2641 Range front_in_the_body;
2642 front_in_the_body = subtract(surface_skin, tets_skin_edges);
2643 Range front_in_the_body_verts;
2644 CHKERR moab.get_connectivity(front_in_the_body, front_in_the_body_verts,
2647 CHKERR skin.find_skin(0, front_in_the_body,
false, front_ends);
2649 intersect(front_in_the_body_verts, constrain_surface_verts));
2650 sets_map[FRONT_ENDS].swap(front_ends);
2652 Range surface_skin_verts;
2653 CHKERR moab.get_connectivity(surface_skin, surface_skin_verts,
true);
2654 sets_map[SURFACE_SKIN].swap(surface_skin_verts);
2657 Range surface_verts;
2659 sets_map[SURFACE_SKIN].merge(
2660 intersect(constrain_surface_verts, surface_verts));
2661 sets_map[SURFACE].swap(surface_verts);
2664 Range tets_skin_verts;
2665 CHKERR moab.get_connectivity(tets_skin, tets_skin_verts,
true);
2666 sets_map[SKIN].swap(tets_skin_verts);
2667 sets_map[SKIN].merge(constrain_surface_verts);
2670 CHKERR moab.get_connectivity(tets, tets_verts,
true);
2671 sets_map[FREE].swap(tets_verts);
2677 Range &proc_tets)
const {
2678 moab::Interface &moab(mField.
get_moab());
2680 Range edges_to_merge_verts;
2681 CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts,
true);
2682 Range edges_to_merge_verts_tets;
2683 CHKERR moab.get_adjacencies(edges_to_merge_verts, 3,
false,
2684 edges_to_merge_verts_tets,
2685 moab::Interface::UNION);
2686 edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
2687 proc_tets.swap(edges_to_merge_verts_tets);
2692 const Range &fixed_edges,
2693 const Range &corner_nodes,
2694 const Range &constrain_surface,
2695 Range &edges_to_merge,
2696 Range ¬_merged_edges) {
2697 moab::Interface &moab(mField.
get_moab());
2701 Skinner skin(&moab);
2703 CHKERR skin.find_skin(0, tets,
false, tets_skin);
2708 Range constrain_surface_verts;
2709 CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
2711 Range constrain_surface_edges;
2712 CHKERR moab.get_adjacencies(constrain_surface, 1,
false,
2713 constrain_surface_edges,
2714 moab::Interface::UNION);
2717 Range tets_skin_edges;
2718 CHKERR moab.get_adjacencies(tets_skin, 1,
false, tets_skin_edges,
2719 moab::Interface::UNION);
2721 Range surface_front;
2722 surface_front = subtract(surface_skin, tets_skin_edges);
2723 Range surface_front_nodes;
2724 CHKERR moab.get_connectivity(surface_front, surface_front_nodes,
true);
2727 CHKERR skin.find_skin(0, surface_front,
false, ends_nodes);
2728 ends_nodes.merge(intersect(surface_front_nodes, constrain_surface_verts));
2731 surface_skin.merge(constrain_surface);
2732 tets_skin_edges.merge(constrain_surface_edges);
2735 Range surface_edges;
2736 CHKERR moab.get_adjacencies(
surface, 1,
false, surface_edges,
2737 moab::Interface::UNION);
2739 Range surface_edges_verts;
2740 CHKERR moab.get_connectivity(surface_edges, surface_edges_verts,
true);
2742 Range tets_skin_edges_verts;
2743 CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
2746 Range edges_to_remove;
2750 Range ents_nodes_and_edges;
2751 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2752 ents_nodes_and_edges.merge(tets_skin_edges);
2753 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2756 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2757 not_merged_edges.merge(edges_to_remove);
2761 Range ents_nodes_and_edges;
2762 ents_nodes_and_edges.merge(surface_edges_verts);
2763 ents_nodes_and_edges.merge(surface_edges);
2764 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2765 ents_nodes_and_edges.merge(tets_skin_edges);
2766 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2769 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2770 not_merged_edges.merge(edges_to_remove);
2773 Range fixed_edges_nodes;
2774 CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes,
true);
2776 Range ents_nodes_and_edges;
2777 ents_nodes_and_edges.merge(fixed_edges_nodes);
2778 ents_nodes_and_edges.merge(ends_nodes);
2779 ents_nodes_and_edges.merge(corner_nodes);
2780 ents_nodes_and_edges.merge(fixed_edges);
2781 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2784 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2785 not_merged_edges.merge(edges_to_remove);
2788 CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove,
false);
2789 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2790 not_merged_edges.merge(edges_to_remove);
2794 Range ents_nodes_and_edges;
2795 ents_nodes_and_edges.merge(surface_skin);
2796 ents_nodes_and_edges.merge(fixed_edges_nodes);
2797 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2800 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2801 not_merged_edges.merge(edges_to_remove);
2805 Range ents_nodes_and_edges;
2806 ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2807 ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2808 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2811 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2812 not_merged_edges.merge(edges_to_remove);
2816 Range ents_nodes_and_edges;
2817 ents_nodes_and_edges.merge(surface_front_nodes);
2818 ents_nodes_and_edges.merge(surface_front);
2819 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2820 ents_nodes_and_edges.merge(tets_skin_edges);
2821 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2824 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2825 not_merged_edges.merge(edges_to_remove);
2832 Range &edges_to_remove,
2834 moab::Interface &moab(mField.
get_moab());
2837 Range ents_nodes = ents.subset_by_type(MBVERTEX);
2838 if (ents_nodes.empty()) {
2839 CHKERR moab.get_connectivity(ents, ents_nodes,
true);
2842 Range ents_nodes_edges;
2843 CHKERR moab.get_adjacencies(ents_nodes, 1,
false, ents_nodes_edges,
2844 moab::Interface::UNION);
2846 Range ents_nodes_edges_nodes;
2847 CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2850 ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2851 Range ents_nodes_edges_nodes_edges;
2852 CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1,
false,
2853 ents_nodes_edges_nodes_edges,
2854 moab::Interface::UNION);
2857 subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2859 subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2861 edges_to_remove.swap(ents_nodes_edges);
2869 Range not_merged_edges;
2871 .removeBadEdges(
surface, tets, fixed_edges, corner_nodes,
2873 Toplogy::SetsMap sets_map;
2878 for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2879 sit != sets_map.rend(); sit++) {
2880 std::string name =
"classification_verts_" +
2881 boost::lexical_cast<std::string>(sit->first) +
".vtk";
2882 if (!sit->second.empty())
2887 CHKERR Toplogy(m_field,
th).getProcTets(tets, edges_to_merge, proc_tets);
2888 out_tets = subtract(tets, proc_tets);
2891 Range all_out_ents = out_tets;
2892 for (
int dd = 2; dd >= 0; dd--) {
2893 CHKERR moab.get_adjacencies(out_tets, dd,
false, all_out_ents,
2894 moab::Interface::UNION);
2900 int nb_nodes_merged = 0;
2904 auto save_merge_step = [&](
const int pp,
const Range collapsed_edges) {
2907 CHKERR moab.get_adjacencies(proc_tets, 2,
false, adj_faces,
2908 moab::Interface::UNION);
2910 name =
"node_merge_step_" + boost::lexical_cast<std::string>(pp) +
".vtk";
2912 name, unite(intersect(new_surf, adj_faces), collapsed_edges));
2914 "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) +
".vtk";
2916 name, unite(intersect(new_surf, adj_faces), edges_to_merge));
2923 double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2926 int nb_nodes_merged_p = nb_nodes_merged;
2933 if(!length_map.empty())
2934 min = length_map.get<2>().begin()->qUality;
2940 Range collapsed_edges;
2941 MergeNodes merge_nodes(m_field,
true,
th, update_meshsets);
2943 for (
auto mit = length_map.get<0>().begin();
2944 mit != length_map.get<0>().end(); mit++, nn++) {
2948 auto get_father_and_mother =
2953 CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes,
true);
2954 std::array<int, 2> conn_type = {0, 0};
2955 for (
int nn = 0; nn != 2; nn++) {
2956 for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2957 sit != sets_map.rend(); sit++) {
2958 if (sit->second.find(conn[nn]) != sit->second.end()) {
2959 conn_type[nn] |= sit->first;
2963 int type_father, type_mother;
2964 if (conn_type[0] > conn_type[1]) {
2967 type_father = conn_type[0];
2968 type_mother = conn_type[1];
2972 type_father = conn_type[1];
2973 type_mother = conn_type[0];
2975 if (type_father == type_mother) {
2981 int line_search = 0;
2983 CHKERR get_father_and_mother(father, mother, line_search);
2984 CHKERR merge_nodes.mergeNodes(line_search, father, mother, proc_tets);
2986 const EntityHandle father_and_mother[] = {father, mother};
2988 CHKERR moab.get_adjacencies(father_and_mother, 1, 3,
false, adj_tets);
2989 Range adj_tets_nodes;
2990 CHKERR moab.get_connectivity(adj_tets, adj_tets_nodes,
true);
2992 CHKERR moab.get_adjacencies(adj_tets_nodes, 1,
false, adj_edges,
2993 moab::Interface::UNION);
2994 for (
auto ait : adj_edges) {
2995 auto miit = length_map.get<1>().find(ait);
2996 if (miit != length_map.get<1>().end())
3000 collapsed_edges.insert(mit->eDge);
3003 if (nn >
static_cast<int>(length_map.size() / fraction_level))
3005 if (mit->qUality > ave)
3010 CHKERR merge_nodes.updateRangeByChilds(new_surf, edges_to_merge,
3011 not_merged_edges,
true);
3014 "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e", pp,
3015 nb_nodes_merged, ave, min);
3018 CHKERR save_merge_step(pp + 1, collapsed_edges);
3020 if (nb_nodes_merged == nb_nodes_merged_p)
3022 if (min > 1e-2 && min == min_pp)
3028 CHKERR moab.get_adjacencies(proc_tets, 1,
false, adj_edges,
3029 moab::Interface::UNION);
3030 edges_to_merge = intersect(edges_to_merge, adj_edges);
3032 .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
3036 auto reconstruct_refined_ents = [&]() {
3047 CHKERR reconstruct_refined_ents();
3053 out_tets.merge(proc_tets);
3055 CHKERR moab.get_adjacencies(out_tets, 2,
false, adj_faces,
3056 moab::Interface::UNION);
3057 new_surf = intersect(new_surf, adj_faces);