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);