12                {
   13 
   15CutMeshInterface::query_interface(boost::typeindex::type_index type_index,
   16                                  UnknownInterface **iface) const {
   18  *iface = const_cast<CutMeshInterface *>(this);
   20}
   21 
   22CutMeshInterface::CutMeshInterface(const Core &core)
   23    : cOre(const_cast<
Core &>(core)) {
 
   24  lineSearchSteps = 10;
   25  nbMaxMergingCycles = 200;
   26  projectEntitiesQualityTrashold = 0.5;
   27}
   28 
   31  treeSurfPtr.reset();
   33}
   34 
   37  fRont = front;
   39}
   40 
   45}
   46 
   48                                             double *shift, double *origin,
   49                                             double *transform,
   50                                             const std::string save_mesh) {
   51  CoreInterface &m_field = cOre;
   52  moab::Interface &moab = m_field.get_moab();
   54  sUrface.clear();
   55  std::map<EntityHandle, EntityHandle> verts_map;
   56  for (Range::const_iterator tit = 
surface.begin(); tit != 
surface.end();
 
   57       tit++) {
   58    int num_nodes;
   60    CHKERR moab.get_connectivity(*tit, conn, num_nodes, 
true);
 
   62    if (th) {
   63      CHKERR moab.tag_get_data(th, conn, num_nodes, &coords(0, 0));
 
   64    } else {
   65      CHKERR moab.get_coords(conn, num_nodes, &coords(0, 0));
 
   66    }
   68    for (int nn = 0; nn != num_nodes; nn++) {
   69      if (verts_map.find(conn[nn]) != verts_map.end()) {
   70        new_verts[nn] = verts_map[conn[nn]];
   71      } else {
   72        if (transform) {
   73          ublas::matrix_row<MatrixDouble> mr(coords, nn);
   74          if (origin) {
   76                3, ublas::shallow_array_adaptor<double>(3, origin));
   77            mr = mr - vec_origin;
   78          }
   80              3, 3, ublas::shallow_array_adaptor<double>(9, transform));
   81          mr = prod(mat_transform, mr);
   82          if (origin) {
   84                3, ublas::shallow_array_adaptor<double>(3, origin));
   85            mr = mr + vec_origin;
   86          }
   87        }
   88        if (shift) {
   89          ublas::matrix_row<MatrixDouble> mr(coords, nn);
   91              3, ublas::shallow_array_adaptor<double>(3, shift));
   92          mr = mr + vec_shift;
   93        }
   94        CHKERR moab.create_vertex(&coords(nn, 0), new_verts[nn]);
 
   95        verts_map[conn[nn]] = new_verts[nn];
   96      }
   97    }
   99    CHKERR moab.create_element(MBTRI, new_verts, num_nodes, ele);
 
  100    sUrface.insert(ele);
  101  }
  102  if (!save_mesh.empty())
  103    CHKERR SaveData(m_field.get_moab())(save_mesh, sUrface);
 
  105}
  106 
  109  vOlume = volume;
  111}
  112 
  115  constrainSurface = surf;
  117}
  118 
  123}
  124 
  127  vOlume.merge(volume);
  129}
  130 
  132                                                        const double rel_tol,
  133                                                        const double abs_tol,
  136  CoreInterface &m_field = cOre;
  137  moab::Interface &moab = m_field.get_moab();
  139 
  140  
  141  Skinner skin(&moab);
  143  CHKERR skin.find_skin(0, sUrface, 
false, surface_skin);
 
  144 
  145  CHKERR snapSurfaceToEdges(surface_skin, fixed_edges, rel_tol, abs_tol, th,
 
  147 
  149}
  150 
  152                                                    const Range fixed_edges,
 
  153                                                    const double rel_tol,
  154                                                    const double abs_tol,
  156  CoreInterface &m_field = cOre;
  157  moab::Interface &moab = m_field.get_moab();
  160 
  161  map<EntityHandle, double> map_verts_length;
  162 
  163  for (auto f : surface_edges) {
  164    int num_nodes;
  166    CHKERR moab.get_connectivity(f, conn_skin, num_nodes, 
true);
 
  168    if (th)
  169      CHKERR moab.tag_get_data(th, conn_skin, num_nodes, &coords_skin[0]);
 
  170    else
  171      CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
 
  173        &coords_skin[0], &coords_skin[1], &coords_skin[2]);
  175        &coords_skin[3], &coords_skin[4], &coords_skin[5]);
  177    t_skin_delta(
i) = t_n1(
i) - t_n0(
i);
 
  178    const double skin_edge_length = sqrt(t_skin_delta(
i) * t_skin_delta(
i));
 
  179    for (int nn = 0; nn != 2; ++nn) {
  180      auto it = map_verts_length.find(conn_skin[nn]);
  181      if (it != map_verts_length.end())
  182        it->second = std::min(it->second, skin_edge_length);
  183      else
  184        map_verts_length[conn_skin[nn]] = skin_edge_length;
  185    }
  186  }
  187 
  188  for (
auto m : map_verts_length) {
 
  189 
  191    if (th)
  192      CHKERR moab.tag_get_data(th, &
m.first, 1, &t_n(0));
 
  193    else
  194      CHKERR moab.get_coords(&
m.first, 1, &t_n(0));
 
  195 
  196    double min_dist = rel_tol * 
m.second;
 
  198    CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
 
  199        &t_n(0), 1, fixed_edges, &min_dist, &t_min_coords(0));
  200 
  201    if (min_dist < rel_tol * 
m.second || min_dist < abs_tol) {
 
  203        MOFEM_LOG(
"WORLD", Sev::noisy) << 
"Snap " << min_dist;
 
  204      if (th)
  205        CHKERR moab.tag_set_data(th, &
m.first, 1, &t_min_coords(0));
 
  206      else
  207        CHKERR moab.set_coords(&
m.first, 1, &t_min_coords(0));
 
  208    }
  209  }
  210 
  212}
  213 
  215  CoreInterface &m_field = cOre;
  216  moab::Interface &moab = m_field.get_moab();
  218  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
  219      new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
  220  CHKERR treeSurfPtr->build(sUrface, rootSetSurf);
 
  222}
  223 
  225CutMeshInterface::cutOnly(
Range vol, 
const BitRefLevel cut_bit, 
Tag th,
 
  226                          const double tol_geometry, const double tol_cut_close,
  228                          const bool update_meshsets, 
const bool debug) {
 
  229  CoreInterface &m_field = cOre;
  230  moab::Interface &moab = m_field.get_moab();
  232 
  233  
  235  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_geometry,
 
  237  CHKERR cutEdgesInMiddle(cut_bit, cutNewVolumes, cutNewSurfaces,
 
  238                          cutNewVertices, 
debug);
 
  239 
  240  CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(constrainSurface,
 
  241                                                         constrainSurface);
  242  if (fixed_edges)
  243    CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(*fixed_edges,
 
  244                                                           *fixed_edges);
  245  if (corner_nodes)
  246    CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(*corner_nodes,
 
  247                                                           *corner_nodes);
  248  if (update_meshsets)
  249    CHKERR m_field.getInterface<MeshsetsManager>()
 
  250        ->updateAllMeshsetsByEntitiesChildren(cut_bit);
  251  CHKERR moveMidNodesOnCutEdges(th);
 
  252 
  255    if (fixed_edges)
  256      CHKERR SaveData(moab)(
"out_cut_new_fixed_edges.vtk", *fixed_edges);
 
  257  }
  258 
  260}
  261 
  263                                          const double tol_trim_close,
  266                                          const bool update_meshsets,
  268  CoreInterface &m_field = cOre;
  269  moab::Interface &moab = m_field.get_moab();
  271 
  272  
  273  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim_close, 
debug);
 
  275 
  276  CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(constrainSurface,
 
  277                                                         constrainSurface);
  278  if (fixed_edges)
  279    CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(*fixed_edges,
 
  280                                                           *fixed_edges);
  281 
  282  if (corner_nodes)
  283    CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(*corner_nodes,
 
  284                                                           *corner_nodes);
  285 
  286  if (update_meshsets)
  287    CHKERR m_field.getInterface<MeshsetsManager>()
 
  288        ->updateAllMeshsetsByEntitiesChildren(trim_bit);
  289 
  290  
  291  CHKERR moveMidNodesOnTrimmedEdges(th);
 
  292 
  293  
  294  CHKERR trimSurface(fixed_edges, corner_nodes, tol_trim_close, th, 
debug);
 
  295 
  299    CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
  300        trim_bit, 
BitRefLevel().set(), MBEDGE, bit2_edges);
 
  301    CHKERR SaveData(moab)(
"trim_fixed_edges.vtk",
 
  302                          intersect(*fixed_edges, bit2_edges));
  303  }
  304 
  306}
  307 
  309    int &first_bit, 
Tag th, 
const double tol_geometry,
 
  310    const double tol_cut_close, 
const double tol_trim_close, 
Range *fixed_edges,
 
  311    Range *corner_nodes, 
const bool update_meshsets, 
const bool debug) {
 
  312  CoreInterface &m_field = cOre;
  314 
  315  std::vector<BitRefLevel> bit_levels;
  316 
  317  auto get_back_bit_levels = [&]() {
  319    bit_levels.back().set(first_bit + bit_levels.size() - 1);
  320    return bit_levels.back();
  321  };
  322 
  324 
  325  CHKERR cutOnly(unite(cutSurfaceVolumes, cutFrontVolumes), cut_bit, th,
 
  326                 tol_geometry, tol_cut_close, fixed_edges, corner_nodes,
  327                 update_meshsets, 
debug);
 
  328 
  331    CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
  333    double min_q = 1;
  334    CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
 
  335    return min_q;
  336  };
  337 
  338  MOFEM_LOG_C(
"WORLD", Sev::inform, 
"Min quality cut %6.4g",
 
  339              get_min_quality(cut_bit, th));
  340 
  341  Range starting_volume = cutNewVolumes;
 
  342  Range starting_volume_nodes;
 
  343  CHKERR m_field.get_moab().get_connectivity(starting_volume,
 
  344                                             starting_volume_nodes, true);
  345  std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
  346  CHKERR m_field.get_moab().get_coords(starting_volume_nodes,
 
  347                                       &*staring_volume_coords.begin());
  348 
  349  if (th) {
  350    std::vector<double> staring_volume_th_coords(3 *
  351                                                 starting_volume_nodes.size());
  352    CHKERR m_field.get_moab().tag_get_data(th, starting_volume_nodes,
 
  353                                           &*staring_volume_th_coords.begin());
  354    CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
 
  355                                         &*staring_volume_th_coords.begin());
  356  }
  357 
  358  if (th)
  360 
  362 
  363  CHKERR trimOnly(trim_bit, th, tol_trim_close, fixed_edges, corner_nodes,
 
  364                  update_meshsets, 
debug);
 
  365 
  366  MOFEM_LOG_C(
"WORLD", Sev::inform, 
"Min quality trim %3.2g",
 
  367              get_min_quality(trim_bit, th));
  368 
  369  first_bit += bit_levels.size() - 1;
  370 
  371  if (th)
  372    CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
 
  373                                         &*staring_volume_coords.begin());
  374 
  376}
  377 
  379    int &first_bit, 
const int fraction_level, 
Tag th, 
const double tol_geometry,
 
  380    const double tol_cut_close, 
const double tol_trim_close, 
Range &fixed_edges,
 
  381    Range &corner_nodes, 
const bool update_meshsets, 
const bool debug) {
 
  382  CoreInterface &m_field = cOre;
  384 
  385  std::vector<BitRefLevel> bit_levels;
  386 
  387  auto get_back_bit_levels = [&]() {
  389    bit_levels.back().set(first_bit + bit_levels.size() - 1);
  390    return bit_levels.back();
  391  };
  392 
  394    CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
 
  395        "ents_not_in_database.vtk", "VTK", "");
  396  }
  397 
  398  CHKERR cutAndTrim(first_bit, th, tol_geometry, tol_cut_close, tol_trim_close,
 
  399                    &fixed_edges, &corner_nodes, update_meshsets, 
debug);
 
  401    CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
 
  402        "cut_trim_ents_not_in_database.vtk", "VTK", "");
  403 
  407 
  408  CHKERR mergeBadEdges(fraction_level, bit_level1, bit_level2, bit_level3,
 
  409                       getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
  410                       update_meshsets, 
debug);
 
  411 
  414    CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
  416    double min_q = 1;
  417    CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
 
  418    if (min_q < 0 && 
debug) {
 
  419      CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
 
  420          "negative_tets.vtk", "VTK", "", tets_level, th);
  421    }
  422    return min_q;
  423  };
  424 
  425  MOFEM_LOG_C(
"WORLD", Sev::inform, 
"Min quality node merge %6.4g",
 
  426              get_min_quality(bit_level3, th));
  427 
  428  CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(constrainSurface,
 
  429                                                         constrainSurface);
  430  CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(fixed_edges,
 
  431                                                         fixed_edges);
  432  CHKERR cOre.getInterface<BitRefManager>()->updateRangeByChildren(corner_nodes,
 
  433                                                         corner_nodes);
  434 
  435  first_bit += bit_levels.size() - 1;
  436 
  438    CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
 
  439        bit_level3, 
BitRefLevel().set(), MBTET, 
"out_tets_merged.vtk", 
"VTK",
 
  440        "");
  441    CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
 
  442        "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
  443    CHKERR SaveData(m_field.get_moab())(
"merged_surface.vtk", mergedSurfaces);
 
  444  }
  445 
  447}
  448 
  450  CoreInterface &m_field = cOre;
  451  moab::Interface &moab = m_field.get_moab();
  453  Skinner skin(&moab);
  455  CHKERR skin.find_skin(0, vOlume, 
false, tets_skin);
 
  456  Range tets_skin_edges;
 
  457  ErrorCode tmp_result;
  458  tmp_result = moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
  459                                    moab::Interface::UNION);
  460 
  461  if (MB_SUCCESS != tmp_result)
  463            "Duplicated edges: most likely the source of error is comming from "
  464            "adding the vertices of the cracking "
  465            "volume to a BLOCKSET rather than NODESET (corresponding to the "
  466            "input parameter-vertex_block_set)");
  467 
  469  CHKERR skin.find_skin(0, sUrface, 
false, surface_skin);
 
  470  fRont = subtract(surface_skin, tets_skin_edges);
  472    CHKERR SaveData(moab)(
"front.vtk", fRont);
 
  474}
  475 
  478  CoreInterface &m_field = cOre;
  479  moab::Interface &moab = m_field.get_moab();
  481  auto tools_interface = m_field.getInterface<Tools>();
  482 
  483  auto create_tag = [&](const std::string name, const int dim) {
  485    rval = moab.tag_get_handle(name.c_str(), 
th);
 
  486    if (rval == MB_SUCCESS)
  488    std::vector<double> def_val(dim, 0);
  489    CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
 
  490                               MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
  491 
  493  };
  494 
  495  auto set_vol = [&](
const Range &vol_verts, std::vector<double> &coords,
 
  496                     std::vector<double> &dist_surface_vec,
  497                     std::vector<double> &dist_surface_normal_vec,
  498                     std::vector<double> &dist_surface_signed_dist_vec) {
  500 
  501    coords.resize(3 * vol_verts.size());
  502    dist_surface_vec.resize(3 * vol_verts.size());
  503    dist_surface_normal_vec.resize(3 * vol_verts.size());
  504    dist_surface_signed_dist_vec.resize(vol_verts.size());
  505 
  506    CHKERR moab.get_coords(vol_verts, &*coords.begin());
 
  507    std::srand(0);
  508 
  509    for (
auto v : vol_verts) {
 
  510 
  511      const int index = vol_verts.index(
v);
 
  515      CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
 
  516                                              &point_out[0], facets_out);
  517 
  520      CHKERR tools_interface->getTriNormal(facets_out, &*n_first.begin());
 
  521 
  522      
  523      
  524      auto check_triangle_orientation = [&](
auto n) {
 
  525        int num_nodes;
  527        CHKERR moab.get_connectivity(facets_out, conn, num_nodes, 
true);
 
  529        CHKERR moab.get_coords(conn, 3, &coords(0, 0));
 
  531        center.clear();
  532        for (auto ii : {0, 1, 2})
  533          for (auto jj : {0, 1, 2})
  534            center(jj) += coords(ii, jj) / 3;
  535 
  536        std::vector<EntityHandle> triangles_out;
  537        std::vector<double> distance_out;
  538        auto a_max = norm_2(
n);
 
  541        const double ray_length = 2 * sqrt(a_max);
  542 
  543        CHKERR treeSurfPtr->ray_intersect_triangles(
 
  544            distance_out, triangles_out, rootSetSurf,
  545            std::numeric_limits<float>::epsilon(), &pt[0], &ray[0],
  546            &ray_length);
  547 
  548        if (triangles_out.size() > 1) {
  550          for (
auto t : triangles_out) {
 
  551            CHKERR tools_interface->getTriNormal(
t, &*nt.begin());
 
  552            double at = norm_2(nt);
  553            if (at > a_max) {
  555              a_max = at;
  556            }
  557          }
  558        }
  559 
  561      };
  562 
  563      auto n = check_triangle_orientation(n_first);
 
  565      const double dot = inner_prod(
delta, 
n);
 
  566 
  567 
  568      if (std::abs(dot) < std::numeric_limits<double>::epsilon()) {
  569        if (std::rand() % 2 == 0)
  570          delta += 
n * std::numeric_limits<double>::epsilon();
 
  571        else
  572          delta -= 
n * std::numeric_limits<double>::epsilon();
 
  573      }
  574 
  576      noalias(dist_vec) = 
delta;
 
  577 
  578      auto dist_normal_vec =
  580 
  581      dist_surface_signed_dist_vec[index] = dot;
  582      noalias(dist_normal_vec) = dot * 
n;
 
  583    }
  584 
  586  };
  587 
  588  std::vector<double> coords;
  589  std::vector<double> dist_surface_vec;
  590  std::vector<double> dist_surface_normal_vec;
  591  std::vector<double> dist_surface_signed_dist_vec;
  593  CHKERR moab.get_connectivity(vOlume, vol_verts, 
true);
 
  594 
  595  CHKERR set_vol(vol_verts, coords, dist_surface_vec, dist_surface_normal_vec,
 
  596                 dist_surface_signed_dist_vec);
  597 
  598  CHKERR moab.tag_set_data(create_tag(
"DIST_SURFACE_VECTOR", 3), vol_verts,
 
  599                           &*dist_surface_vec.begin());
  600  CHKERR moab.tag_set_data(create_tag(
"DIST_SURFACE_NORMAL_VECTOR", 3),
 
  601                           vol_verts, &*dist_surface_normal_vec.begin());
  602  CHKERR moab.tag_set_data(create_tag(
"DIST_SURFACE_NORMAL_SIGNED", 1),
 
  603                           vol_verts, &*dist_surface_signed_dist_vec.begin());
  604 
  606}
  607 
  609                                                      int verb,
  611  CoreInterface &m_field = cOre;
  612  moab::Interface &moab = m_field.get_moab();
  614 
  615  auto create_tag = [&](const std::string name, const int dim) {
  617    rval = moab.tag_get_handle(name.c_str(), 
th);
 
  618    if (rval == MB_SUCCESS)
  620    std::vector<double> def_val(dim, 0);
  621    CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
 
  622                               MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
  624  };
  625 
  627  CHKERR moab.get_connectivity(vol, vol_vertices, 
true);
 
  628 
  629  std::vector<double> min_distances_from_front(vol_vertices.size(), -1);
  630  std::vector<double> points_on_edges(3 * vol_vertices.size(), 0);
  631  std::vector<EntityHandle> closest_edges(vol_vertices.size(), 0);
  632 
  633  std::vector<double> coords(3 * vol_vertices.size());
  634  if (th)
  635    CHKERR moab.tag_get_data(th, vol_vertices, &*coords.begin());
 
  636  else
  637    CHKERR moab.get_coords(vol_vertices, &*coords.begin());
 
  638 
  639  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
 
  640      &*coords.begin(), vol_vertices.size(), fRont,
  641      &*min_distances_from_front.begin(), &*points_on_edges.begin(),
  642      &*closest_edges.begin());
  643 
  644  if (!points_on_edges.empty()) {
  645    for (
int i = 0; 
i != min_distances_from_front.size(); ++
i) {
 
  647      CHKERR moab.get_adjacencies(&closest_edges[0], 1, 2, 
false, faces);
 
  650      point_out -= point_in;
  651    }
  652  }
  653 
  654  auto th_dist_front_vec = create_tag("DIST_FRONT_VECTOR", 3);
  655  CHKERR moab.tag_set_data(th_dist_front_vec, vol_vertices,
 
  656                           &*points_on_edges.begin());
  657 
  659}
  660 
  662    Tag th, 
Range &vol_edges, 
const bool remove_adj_prims_edges, 
int verb,
 
  663    const bool debug, 
const std::string edges_file_name) {
 
  664  CoreInterface &m_field = cOre;
  665  moab::Interface &moab = m_field.get_moab();
  667 
  668  auto get_tag_dim = [&](
auto th) {
 
  669    int dim;
  670    CHKERR moab.tag_get_length(th, dim);
 
  671    return dim;
  672  };
  673  auto dim = get_tag_dim(th);
  674 
  675  auto get_tag_data = [&](
auto th, 
auto conn) {
 
  676    const void *ptr;
  677    CHKERR moab.tag_get_by_ptr(th, &conn, 1, &ptr);
 
  679        const_cast<double *>(static_cast<const double *>(ptr)), 3);
  680  };
  681 
  682  auto get_edge_ray = [&](auto conn) {
  684    CHKERR moab.get_coords(conn, 2, &coords[0]);
 
  688    return ray;
  689  };
  690 
  692  CHKERR moab.get_adjacencies(vOlume, 1, 
true, edges, moab::Interface::UNION);
 
  693 
  694  auto remove_prisms_edges = [&](
Range &edges) {
 
  697    CHKERR moab.get_adjacencies(edges, 3, 
false, prisms,
 
  698                                moab::Interface::UNION);
  699    prisms = prisms.subset_by_type(MBPRISM);
  701    CHKERR moab.get_connectivity(prisms, prisms_verts, 
true);
 
  703    CHKERR moab.get_adjacencies(prisms_verts, 1, 
false, prism_edges,
 
  704                                moab::Interface::UNION);
  705    edges = subtract(edges, prism_edges);
  707  };
  708  if (remove_adj_prims_edges)
  709    CHKERR remove_prisms_edges(edges);
 
  710 
  711  std::vector<EntityHandle> cut_edges;
  712  cut_edges.reserve(edges.size());
  713 
  714  auto get_cut_edges_vec = [&](
auto th, 
auto &cut_edges, 
auto e, 
auto &&conn) {
 
  716 
  717    auto ray = get_edge_ray(conn);
  718    const double length = norm_2(ray);
  719    ray /= length;
  720    const auto dist0 = get_tag_data(th, conn[0]);
  721    const auto dist1 = get_tag_data(th, conn[1]);
  722    const double max_dist = std::max(norm_2(dist0), norm_2(dist1));
  723    if (max_dist < length) {
  724      cut_edges.push_back(e);
  725    }
  726 
  728  };
  729 
  730  auto get_cut_edges_signed_dist = [&](
auto th, 
auto &cut_edges, 
auto e,
 
  731                                       auto &&conn) {
  733    auto get_tag_signed_dist = [&](auto conn) {
  734      double dist;
  735      CHKERR moab.tag_get_data(th, &conn, 1, &dist);
 
  736      return dist;
  737    };
  738    const auto dist0 = get_tag_signed_dist(conn[0]);
  739    const auto dist1 = get_tag_signed_dist(conn[1]);
  740    const double opposite = dist0 * dist1;
  741    if (opposite <= 0)
  742      cut_edges.push_back(e);
  744  };
  745 
  746  auto get_conn = [&](auto e) {
  747    int num_nodes;
  749    CHKERR moab.get_connectivity(e, conn, num_nodes, 
true);
 
  750    return conn;
  751  };
  752 
  753  if (dim == 3)
  754    for (auto e : edges)
  755      CHKERR get_cut_edges_vec(
th, cut_edges, e, get_conn(e));
 
  756  else
  757    for (auto e : edges)
  758      CHKERR get_cut_edges_signed_dist(
th, cut_edges, e, get_conn(e));
 
  759 
  760  CHKERR moab.get_adjacencies(&*cut_edges.begin(), cut_edges.size(), 3, 
false,
 
  761                              vol_edges, moab::Interface::UNION);
  762  vol_edges = intersect(vol_edges, vOlume);
  763 
  764  auto convert_to_ranger = [](
auto &
v) {
 
  766    r.insert_list(
v.begin(), 
v.end());
 
  768  };
  769 
  770  if (
debug && !edges_file_name.empty())
 
  771    CHKERR SaveData(m_field.get_moab())(edges_file_name,
 
  772                                        convert_to_ranger(cut_edges));
  773 
  775}
  776 
  779  CoreInterface &m_field = cOre;
  780  moab::Interface &moab = m_field.get_moab();
  782 
  783  CHKERR createFrontLevelSets(vOlume, 
nullptr, verb, 
debug);
 
  784  Tag th_dist_front_vec;
 
  785  CHKERR moab.tag_get_handle(
"DIST_FRONT_VECTOR", th_dist_front_vec);
 
  786  CHKERR findLevelSetVolumes(th_dist_front_vec, cutFrontVolumes, 
true, verb,
 
  787                             debug, 
"cutFrontEdges.vtk");
 
  788 
  790 
  791  Tag th_dist_surface_vec;
 
  792  CHKERR moab.tag_get_handle(
"DIST_SURFACE_VECTOR", th_dist_surface_vec);
 
  793  cutSurfaceVolumes.clear();
  794  CHKERR findLevelSetVolumes(th_dist_surface_vec, cutSurfaceVolumes, 
true, verb,
 
  795                             debug, 
"cutSurfaceEdges.vtk");
 
  796 
  798    CHKERR SaveData(m_field.get_moab())(
"level_sets.vtk", vOlume);
 
  800    CHKERR SaveData(m_field.get_moab())(
"cutSurfaceVolumes.vtk",
 
  801                                        cutSurfaceVolumes);
  803    CHKERR SaveData(m_field.get_moab())(
"cutFrontVolumes.vtk", cutFrontVolumes);
 
  804 
  806}
  807 
  808MoFEMErrorCode CutMeshInterface::refineMesh(
const int init_bit_level,
 
  809                                            const int surf_levels,
  810                                            const int front_levels,
  811                                            Range *fixed_edges, 
int verb,
 
  813  CoreInterface &m_field = cOre;
  814  moab::Interface &moab = m_field.get_moab();
  815  MeshRefinement *refiner;
  816  BitRefManager *bit_ref_manager;
  818  CHKERR m_field.getInterface(refiner);
 
  819  CHKERR m_field.getInterface(bit_ref_manager);
 
  820 
  821  auto add_bit = [&](
const int bit) {
 
  824                                           verb);
  826    for (auto d : {2, 1, 0}) {
  827      CHKERR moab.get_adjacencies(vOlume, d, 
true, adj_ents,
 
  828                                  moab::Interface::UNION);
  830                                             verb);
  831    }
  833  };
  834  CHKERR add_bit(init_bit_level);
 
  835 
  836  auto update_range = [&](
Range *r_ptr) {
 
  838    if (r_ptr) {
  840      CHKERR bit_ref_manager->updateRangeByChildren(*r_ptr, childs);
 
  841      r_ptr->merge(childs);
  842    }
  844  };
  845 
  849    CHKERR moab.get_connectivity(tets, verts, 
true);
 
  851    CHKERR moab.get_adjacencies(verts, 1, 
true, ref_edges,
 
  852                                moab::Interface::UNION);
  853 
  854    CHKERR refiner->addVerticesInTheMiddleOfEdges(ref_edges, 
bit);
 
  855    CHKERR refiner->refineTets(vOlume, 
bit, verb);
 
  856 
  857    CHKERR update_range(fixed_edges);
 
  858    CHKERR update_range(&vOlume);
 
  859    CHKERR m_field.getInterface<MeshsetsManager>()
 
  860        ->updateAllMeshsetsByEntitiesChildren(
bit);
 
  861 
  863    CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
  865    vOlume = intersect(vOlume, bit_tets);
  866 
  868    CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
  870    if (fixed_edges)
  871      *fixed_edges = intersect(*fixed_edges, bit_edges);
  872 
  874  };
  875 
  876  for (int ll = init_bit_level; ll != init_bit_level + surf_levels; ++ll) {
  879                  unite(cutSurfaceVolumes, cutFrontVolumes));
  880  }
  881 
  882  for (int ll = init_bit_level + surf_levels;
  883       ll != init_bit_level + surf_levels + front_levels; ++ll) {
  886  }
  887 
  889 
  891    CHKERR SaveData(m_field.get_moab())(
"refinedVolume.vtk", vOlume);
 
  892 
  894}
  895 
  898  CoreInterface &m_field = cOre;
  899  moab::Interface &moab = m_field.get_moab();
  901 
  902  edgesToCut.clear();
  903  cutEdges.clear();
  904 
  906  CHKERR moab.tag_get_handle(
"DIST_SURFACE_NORMAL_SIGNED", th_signed_dist);
 
  907 
  908  auto get_tag_edge_dist = [&](
auto th, 
auto conn) {
 
  909    std::array<double, 2> 
r;
 
  910    CHKERR moab.tag_get_data(th, conn, 2, 
r.data());
 
  912  };
  913 
  914  auto get_conn = [&](const auto e) {
  915    int num_nodes;
  917    CHKERR moab.get_connectivity(e, conn, num_nodes, 
true);
 
  918    return conn;
  919  };
  920 
  921  auto get_adj = [&moab](
const Range r, 
const int dim) {
 
  923    if (dim)
  924      CHKERR moab.get_adjacencies(r, dim, 
false, 
a, moab::Interface::UNION);
 
  925    else
  926      CHKERR moab.get_connectivity(r, 
a, 
true);
 
  928  };
  929 
  930  auto vol_edges = get_adj(vol, 1);
  931 
  932  aveLength = 0;
  933  maxLength = 0;
  934  int nb_ave_length = 0;
  935  for (auto e : vol_edges) {
  936 
  937    auto conn = get_conn(e);
  938    auto dist = get_tag_edge_dist(th_signed_dist, conn);
  939    const auto dist_max = std::max(dist[0], dist[1]);
  940    const auto dist_min = std::min(dist[0], dist[1]);
  941    const auto dot = dist[0] * dist[1];
  942 
  944    CHKERR moab.get_coords(conn, 2, &coords[0]);
 
  948    const double ray_length = norm_2(ray);
  949    ray /= ray_length;
  950 
  951    if (
  952 
  953        (dot < 0 && dist_max > std::numeric_limits<double>::epsilon()) ||
  954        (std::abs(dist_min) < std::numeric_limits<double>::epsilon() &&
  955         dist_max > std::numeric_limits<double>::epsilon())
  956 
  957    ) {
  958 
  959      
  960      const double s = dist[0] / (dist[0] - dist[1]);
  961      const double dist = s * ray_length;
  962 
  963      auto add_edge = [&](auto dist) {
  964        edgesToCut[e].dIst = dist;
  965        edgesToCut[e].lEngth = ray_length;
  966        edgesToCut[e].unitRayDir = ray;
  967        edgesToCut[e].rayPoint = n0;
  968        cutEdges.insert(e);
  969 
  970        aveLength += norm_2(ray);
  971        maxLength = fmax(maxLength, norm_2(ray));
  972        ++nb_ave_length;
  973      };
  974 
  975      add_edge(dist);
  976    }
  977  }
  978  aveLength /= nb_ave_length;
  979 
  981    CHKERR SaveData(m_field.get_moab())(
"cut_edges.vtk", cutEdges);
 
  982 
  984    CHKERR SaveData(m_field.get_moab())(
"cut_edges_vol.vtk",
 
  985                                        get_adj(cutEdges, 3));
  986 
  988}
  989 
  991    Range *fixed_edges, 
Range *corner_nodes, 
const double geometry_tol,
 
  992    const double close_tol, 
const int verb, 
const bool debug) {
 
  993  CoreInterface &m_field = cOre;
  994  moab::Interface &moab = m_field.get_moab();
  995  Skinner skin(&moab);
  997 
  998  auto get_ent_adj = [&moab](
const EntityHandle v, 
const int dim) {
 
 1000    if (dim)
 1001      CHKERR moab.get_adjacencies(&
v, 1, dim, 
true, 
a);
 
 1003  };
 1004 
 1005  auto get_adj = [&moab](
const Range r, 
const int dim) {
 
 1007    if (dim)
 1008      CHKERR moab.get_adjacencies(r, dim, 
false, 
a, moab::Interface::UNION);
 
 1009    else
 1010      CHKERR moab.get_connectivity(r, 
a, 
true);
 
 1012  };
 1013 
 1014  auto get_skin = [&skin](
const auto r) {
 
 1016    CHKERR skin.find_skin(0, r, 
false, s);
 
 1017    return s;
 1018  };
 1019 
 1020 
 1021  auto get_range = [](std::vector<EntityHandle> 
v) {
 
 1023    r.insert_list(
v.begin(), 
v.end());
 
 1025  };
 1026 
 1027  auto get_coords = [&](
const auto v) {
 
 1029    CHKERR moab.get_coords(&
v, 1, &coords[0]);
 
 1030    return coords;
 1031  };
 1032 
 1034  CHKERR moab.tag_get_handle(
"DIST_SURFACE_NORMAL_VECTOR", th_normal_dist);
 
 1035 
 1038    CHKERR moab.tag_get_data(th_normal_dist, &
v, 1, &*dist_vec.begin());
 
 1039    return dist_vec;
 1040  };
 1041 
 1043                           const double geometry_tol) {
 1044    auto get_tuple = [](
const EntityHandle e, 
const double dist,
 
 1045                        const double l) { 
return std::make_tuple(e, dist, 
l); };
 
 1046 
 1047    std::tuple<EntityHandle, double, double> min_tuple{0, 0, 0};
 1048 
 1049    for (auto e : edges) {
 1050      auto eit = edgesToCut.find(e);
 1051      if (eit != edgesToCut.end()) {
 1052 
 1053        auto get_dist = [&](auto eit) {
 1054          int num_nodes;
 1056          CHKERR moab.get_connectivity(eit->first, conn, num_nodes, 
true);
 
 1058            return eit->second.dIst;
 1059          else if (conn[1] == 
v)
 
 1060            return (eit->second.lEngth - eit->second.dIst);
 1061          else
 1063        };
 1064 
 1065        const auto d = get_dist(eit);
 
 1066        if (std::get<0>(min_tuple) == 0) {
 1067          min_tuple = get_tuple(e, d, eit->second.lEngth);
 1068        } else if (std::get<1>(min_tuple) > d) {
 1069          min_tuple = get_tuple(e, d, eit->second.lEngth);
 1070        }
 1071      }
 1072    }
 1073 
 1074    const auto geom_dist_vec = get_normal_dits(
v);
 
 1075    const auto geom_tol = norm_2(geom_dist_vec);
 1076    const auto l = std::get<2>(min_tuple);
 
 1077 
 1078    if (geom_tol < 
l * geometry_tol) {
 
 1079      return std::make_pair(get_coords(
v), 
l);
 
 1080 
 1081    } else {
 1082      const auto &
d = edgesToCut.at(std::get<0>(min_tuple));
 
 1085    }
 1086  };
 1087 
 1088  auto get_in_range = [](
auto v, 
auto &
r) { 
return (
r.find(
v) != 
r.end()); };
 
 1089 
 1090  auto project_nodes = [&](auto nodes_to_check) {
 1092      CHKERR SaveData(moab)(
"nodes_to_check.vtk", nodes_to_check);
 
 1093 
 1094    auto get_fix_e = [](auto fixed_edges) {
 1095      if (fixed_edges)
 1096        return *fixed_edges;
 1097      else
 1099    };
 1100 
 1101    const auto fix_e = get_fix_e(fixed_edges);
 1102    const auto skin_e = get_adj(unite(get_skin(vOlume), constrainSurface), 1);
 1103    const auto cut_fix = intersect(cutEdges, fix_e);
 1104    const auto cut_skin = intersect(cutEdges, skin_e);
 1105 
 1107    if (corner_nodes)
 1108      corner_n = intersect(*corner_nodes, nodes_to_check);
 1109    auto fixe_n = intersect(get_adj(fix_e, 0), nodes_to_check);
 1110    auto skin_n = intersect(get_adj(skin_e, 0), nodes_to_check);
 1111 
 1112    std::vector<std::pair<EntityHandle, TreeData>> vertices_on_cut_edges;
 1113    vertices_on_cut_edges.reserve(nodes_to_check.size());
 1114 
 1115    auto add_zero_vertex = [&](
const EntityHandle v, 
auto &&new_pos) {
 
 1116      auto coords = get_coords(
v);
 
 1117      auto ray = new_pos - coords;
 1118      auto l = norm_2(ray);
 
 1119      TreeData data;
 1121      if (
l > std::numeric_limits<double>::epsilon()) {
 
 1122        data.unitRayDir = ray / 
l;
 
 1123      } else {
 1124        data.dIst = 0;
 1125        data.unitRayDir = ray;
 1126      }
 1127      data.rayPoint = coords;
 1128      return std::make_pair(
v, data);
 
 1129    };
 1130 
 1131    auto intersect_v = [&](
const auto v, 
const auto r) {
 
 1132      return intersect(r, get_ent_adj(
v, 1));
 
 1133    };
 1134 
 1135    for (
auto v : nodes_to_check) {
 
 1136 
 1137      const auto e = intersect_v(
v, cutEdges);
 
 1138      if (!e.empty()) {
 1139 
 1140        if (get_in_range(
v, corner_n)) {
 
 1141          auto p = get_prj_point(
v, e, 0);
 
 1142          if (norm_2(get_coords(
v) - p.first) < close_tol * p.second) {
 
 1143            vertices_on_cut_edges.push_back(add_zero_vertex(
v, get_coords(
v)));
 
 1144            continue;
 1145          }
 1146 
 1147        } 
else if (get_in_range(
v, fixe_n)) {
 
 1148          const auto b = intersect_v(
v, cut_fix);
 
 1149          if (!b.empty()) {
 1150            auto p = get_prj_point(
v, b, geometry_tol);
 
 1151            if (norm_2(get_coords(
v) - p.first) < close_tol * p.second) {
 
 1152              vertices_on_cut_edges.push_back(add_zero_vertex(
v, p.first));
 
 1153              continue;
 1154            }
 1155 
 1156          } 
else if (get_in_range(
v, skin_n)) {
 
 1157            const auto b = intersect_v(
v, cut_skin);
 
 1158            if (!b.empty()) {
 1159              auto p = get_prj_point(
v, b, geometry_tol);
 
 1160              if (norm_2(get_coords(
v) - p.first) < close_tol * p.second) {
 
 1161                vertices_on_cut_edges.push_back(add_zero_vertex(
v, p.first));
 
 1162                continue;
 1163              }
 1164            }
 1165 
 1166          }
 1167        }
 1168 
 1169        auto p = get_prj_point(
v, e, geometry_tol);
 
 1170        if (norm_2(get_coords(
v) - p.first) < close_tol * p.second) {
 
 1171          if (get_in_range(
v, fixe_n) || get_in_range(
v, skin_n))
 
 1172            vertices_on_cut_edges.push_back(add_zero_vertex(
v, get_coords(
v)));
 
 1173          else
 1174            vertices_on_cut_edges.push_back(add_zero_vertex(
v, p.first));
 
 1175          }
 1176        }
 1177        
 1178    }
 1179 
 1180    auto get_distances = [&](auto &data) {
 1181      std::map<EntityHandle, double> dist_map;
 1182      if (!data.empty()) {
 1184        CHKERR moab.tag_get_handle(
"DIST_SURFACE_NORMAL_SIGNED", th);
 
 1185 
 1186        std::vector<EntityHandle> verts;
 1187        verts.reserve(data.size());
 1188        for (auto d : data)
 1189          verts.emplace_back(
d.first);
 
 1190        std::vector<double> distances(verts.size());
 1191        CHKERR moab.tag_get_data(th, &*verts.begin(), verts.size(),
 
 1192                                 &*distances.begin());
 1193        for (
size_t i = 0; 
i != verts.size(); ++
i)
 
 1194          dist_map[verts[
i]] = distances[
i];
 
 1195      }
 1196      return dist_map;
 1197    };
 1198 
 1199    auto dist_map = get_distances(vertices_on_cut_edges);
 1200    if (!dist_map.empty()) {
 1201      auto cmp = [&dist_map](
const auto &
a, 
const auto &b) {
 
 1202        return dist_map[
a.first] < dist_map[b.first];
 
 1203      };
 1204      std::sort(vertices_on_cut_edges.begin(), vertices_on_cut_edges.end(),
 1205                cmp);
 1206    }
 1207 
 1208    return vertices_on_cut_edges;
 1209  };
 1210 
 1211  auto get_min_quality =
 1212      [&](
const Range &adj_tets,
 
 1213          const map<EntityHandle, TreeData> &vertices_on_cut_edges,
 1214          const std::pair<EntityHandle, TreeData> *test_data_ptr) {
 1215        double min_q = 1;
 1216        for (
auto t : adj_tets) {
 
 1217          int num_nodes;
 1219          CHKERR m_field.get_moab().get_connectivity(
t, conn, num_nodes, 
true);
 
 1221          CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
 
 1222 
 1223          auto set_new_coord = [](
auto d) {
 
 1224            return d.rayPoint + 
d.dIst * 
d.unitRayDir;
 
 1225          };
 1226 
 1227          for (
auto n : {0, 1, 2, 3}) {
 
 1229            if (test_data_ptr) {
 1230              if (test_data_ptr->first == conn[
n]) {
 
 1231                noalias(n_coords) = set_new_coord(test_data_ptr->second);
 1232              }
 1233            }
 1234            auto mit = vertices_on_cut_edges.find(conn[
n]);
 
 1235            if (mit != vertices_on_cut_edges.end()) {
 1236              noalias(n_coords) = set_new_coord(mit->second);
 1237            }
 1238          }
 1239          min_q = std::min(min_q, Tools::volumeLengthQuality(&coords[0]));
 1240        }
 1241        return min_q;
 1242      };
 1243 
 1244  auto get_zero_distance_verts = [&](const auto &&vertices_on_cut_edges) {
 1245    std::vector<EntityHandle> zero_dist_vec;
 1246    zero_dist_vec.reserve(vertices_on_cut_edges.size());
 1247    for (
auto t : vertices_on_cut_edges) {
 
 1248 
 1249      auto adj_tet = intersect(vOlume, get_ent_adj(
t.first, 3));
 
 1250      const double q0 = get_min_quality(adj_tet, verticesOnCutEdges, nullptr);
 1251      const double q = get_min_quality(adj_tet, verticesOnCutEdges, &
t);
 
 1252      if (std::isnormal(q) && (q / q0) > projectEntitiesQualityTrashold) {
 1253        verticesOnCutEdges[
t.first] = 
t.second;
 
 1254        zero_dist_vec.push_back(
t.first);
 
 1255      }
 1256    }
 1257    return zero_dist_vec;
 1258  };
 1259 
 1260  auto vol_cut_ents = intersect(vOlume, get_adj(cutEdges, 3));
 1261 
 1262  auto get_zero_distant_ents = [&](auto bridge_ents, const int dim,
 1263                                   const int bridge_dim) {
 1264    auto ents =
 1265        intersect(get_adj(bridge_ents, dim), get_adj(vol_cut_ents, dim));
 1266    auto ents_to_remove =
 1267        get_adj(subtract(get_adj(ents, bridge_dim), bridge_ents), dim);
 1268    return subtract(ents, ents_to_remove);
 1269  };
 1270 
 1271  verticesOnCutEdges.clear();
 1272  zeroDistanceVerts.clear();
 1273  zeroDistanceVerts.merge(
 1274      get_range(get_zero_distance_verts(project_nodes(get_adj(cutEdges, 0)))));
 1275  zeroDistanceEnts = subtract(get_zero_distant_ents(zeroDistanceVerts, 2, 0),
 1276                              get_skin(vOlume));
 1277 
 1279    CHKERR SaveData(moab)(
"zero_distance_verts.vtk", zeroDistanceVerts);
 
 1281    CHKERR SaveData(moab)(
"zero_distance_ents.vtk", zeroDistanceEnts);
 
 1282 
 1283  for (auto f : zeroDistanceVerts) {
 1284    auto adj_edges = get_ent_adj(f, 1);
 1285    for (auto e : adj_edges) {
 1286      cutEdges.erase(e);
 1287      edgesToCut.erase(e);
 1288    }
 1289  }
 1290 
 1292    CHKERR SaveData(moab)(
"cut_edges_to_cut.vtk", cutEdges);
 
 1293 
 1295}
 1296 
 1302  CoreInterface &m_field = cOre;
 1303  moab::Interface &moab = m_field.get_moab();
 1304  MeshRefinement *refiner;
 1305  auto ref_ents_ptr = m_field.get_ref_ents();
 1307 
 1308  if (cutEdges.size() != edgesToCut.size())
 1310 
 1311  auto refine_mesh = [&]() {
 1313    CHKERR m_field.getInterface(refiner);
 
 1314    CHKERR refiner->addVerticesInTheMiddleOfEdges(cutEdges, 
bit);
 
 1317  };
 1318 
 1320 
 1322    if (cutEdges.size() != edgesToCut.size())
 1324 
 1325  cut_vols.clear();
 1326  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 1328  cut_surf.clear();
 1329  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 1330      bit, 
bit, MBTRI, cut_surf);
 
 1331 
 1333    CHKERR SaveData(moab)(
"cut_surf_from_bit.vtk", cut_surf);
 
 1334 
 1335  
 1336  cut_verts.clear();
 1337  cut_verts.merge(zeroDistanceVerts);
 1338 
 1339  for (
auto &
m : edgesToCut) {
 
 1340    auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
 1341        boost::make_tuple(MBVERTEX, 
m.first));
 
 1342    if (vit ==
 1343        ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
 1345              "No vertex on cut edges, that make no sense");
 1346    }
 1347    const boost::shared_ptr<RefEntity> &ref_ent = *vit;
 1348    if ((ref_ent->getBitRefLevel() & 
bit).any()) {
 
 1350      cut_verts.insert(vert);
 1351      verticesOnCutEdges[vert] = 
m.second;
 
 1352    } else {
 1353      SETERRQ(
 1355          "Vertex has wrong bit ref level %s",
 1356          boost::lexical_cast<std::string>(ref_ent->getBitRefLevel()).c_str());
 1357    }
 1358  }
 1359 
 1360  
 1362  Skinner skin(&moab);
 1363  CHKERR skin.find_skin(0, cut_vols, 
false, tets_skin);
 
 1364  tets_skin.merge(constrainSurface);
 1365 
 1366  
 1367  
 1368  
 1369  cut_surf.merge(zeroDistanceEnts.subset_by_type(MBTRI));
 1371  CHKERR moab.get_connectivity(cut_surf, diff_verts, 
true);
 
 1372  diff_verts = subtract(diff_verts, cut_verts);
 1373 
 1374  Range subtract_faces;
 
 1375  CHKERR moab.get_adjacencies(diff_verts, 2, 
false, subtract_faces,
 
 1376                              moab::Interface::UNION);
 1377  cut_surf = subtract(cut_surf, unite(subtract_faces, tets_skin));
 1378 
 1379  cut_verts.clear();
 1380  CHKERR moab.get_connectivity(cut_surf, cut_verts, 
true);
 
 1381 
 1382  
 1383  auto check_for_non_minfold = [&]() {
 1386    CHKERR moab.get_adjacencies(cut_surf, 1, 
false, surf_edges,
 
 1387                                moab::Interface::UNION);
 1388    for (auto e : surf_edges) {
 1389 
 1391      CHKERR moab.get_adjacencies(&e, 1, 2, 
false, faces);
 
 1392      faces = intersect(faces, cut_surf);
 1393      if (faces.size() > 2) {
 1394 
 1395        bool resolved = false;
 1396 
 1397        
 1399        CHKERR moab.get_connectivity(faces, nodes, 
true);
 
 1400        for (
auto n : nodes) {
 
 1402          CHKERR moab.get_adjacencies(&
n, 1, 2, 
false, adj_faces);
 
 1403          adj_faces = intersect(adj_faces, cut_surf);
 1404          if (adj_faces.size() == 1) {
 1405            cut_surf.erase(adj_faces[0]);
 1406            resolved = true;
 1407          }
 1408        }
 1409 
 1410        
 1412        CHKERR moab.get_adjacencies(faces, 1, 
false, adj_edges,
 
 1413                                    moab::Interface::UNION);
 1414        adj_edges = intersect(adj_edges, surf_edges);
 1415        adj_edges.erase(e);
 1416        for (auto other_e : adj_edges) {
 1418          CHKERR moab.get_adjacencies(&other_e, 1, 2, 
false, other_faces);
 
 1419          other_faces = intersect(other_faces, cut_surf);
 1420          if (other_faces.size() > 2) {
 1421            other_faces = intersect(other_faces, faces);
 1422            cut_surf = subtract(cut_surf, other_faces);
 1423            resolved = true;
 1424          }
 1425        }
 1426 
 1427        if (!resolved && !
debug)
 
 1429                  "Non-mainfold surface");
 1430 
 1431        cut_verts.clear();
 1432        CHKERR moab.get_connectivity(cut_surf, cut_verts, 
true);
 
 1433      }
 1434    }
 1436  };
 1437 
 1438  CHKERR check_for_non_minfold();
 
 1439 
 1441    CHKERR SaveData(moab)(
"cut_surf.vtk", cut_surf);
 
 1442 
 1444}
 1445 
 1447  CoreInterface &m_field = cOre;
 1448  moab::Interface &moab = m_field.get_moab();
 1450 
 1451  
 1452  for (
auto m : verticesOnCutEdges) {
 
 1453    double dist = 
m.second.dIst;
 
 1454    VectorDouble3 new_coors = 
m.second.rayPoint + dist * 
m.second.unitRayDir;
 
 1455    if (th)
 1456      CHKERR moab.tag_set_data(th, &
m.first, 1, &new_coors[0]);
 
 1457    else
 1458      CHKERR moab.set_coords(&
m.first, 1, &new_coors[0]);
 
 1459  }
 1460 
 1462}
 1463 
 1465  CoreInterface &m_field = cOre;
 1466  moab::Interface &moab = m_field.get_moab();
 1468  for (
auto &
v : verticesOnTrimEdges) {
 
 1469    double dist = 
v.second.dIst;
 
 1470    VectorDouble3 new_coors = 
v.second.rayPoint + dist * 
v.second.unitRayDir;
 
 1471    if (th)
 1472      CHKERR moab.tag_set_data(th, &
v.first, 1, &new_coors[0]);
 
 1473    else
 1474      CHKERR moab.set_coords(&
v.first, 1, &new_coors[0]);
 
 1475  }
 1477}
 1478 
 1481                                                 const double tol_trim_close,
 1483  CoreInterface &m_field = cOre;
 1484  moab::Interface &moab = m_field.get_moab();
 1486 
 1487  
 1488  Skinner skin(&moab);
 1490  CHKERR skin.find_skin(0, cutNewVolumes, 
false, tets_skin);
 
 1491  tets_skin.merge(constrainSurface);
 1492 
 1493  
 1494  Range tets_skin_verts;
 
 1495  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, 
true);
 
 1496  
 1497  Range tets_skin_edges;
 
 1498  CHKERR moab.get_adjacencies(tets_skin, 1, 
false, tets_skin_edges,
 
 1499                              moab::Interface::UNION);
 1500 
 1501  
 1502  Range cut_surface_edges;
 
 1503  CHKERR moab.get_adjacencies(cutNewSurfaces, 1, 
false, cut_surface_edges,
 
 1504                              moab::Interface::UNION);
 1505  Range cut_surface_verts;
 
 1506  CHKERR moab.get_connectivity(cutNewSurfaces, cut_surface_verts, 
true);
 
 1507 
 1509  if (corner_nodes)
 1510    corners = *corner_nodes;
 1511 
 1513  if (fixed_edges)
 1514    fix_edges = *fixed_edges;
 1515 
 1516  Range fixed_edges_verts;
 
 1517  if (fixed_edges)
 1518    CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts, 
true);
 
 1519 
 1521  if (fRont.empty())
 1522    CHKERR skin.find_skin(0, sUrface, 
false, surface_skin);
 
 1523  else
 1524    surface_skin = fRont;
 1525 
 1528    if (th)
 1529      CHKERR moab.tag_get_data(th, &
v, 1, &point[0]);
 
 1530    else
 1531      CHKERR moab.get_coords(&
v, 1, &point[0]);
 
 1532    return point;
 1533  };
 1534 
 1535  struct VertMap {
 1540        : 
d(
d), 
v(
v), e(e) {}
 
 1541  };
 1542 
 1543  typedef multi_index_container<
 1544      VertMap,
 1545      indexed_by<
 1546          ordered_non_unique<member<VertMap, const double, &VertMap::d>>,
 1547          ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::v>>,
 1548          ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::e>>
 1549 
 1550          >>
 1551      VertMapMultiIndex;
 1552 
 1553  VertMapMultiIndex verts_map;
 1554 
 1556                      const double dist) {
 1557    verts_map.insert(VertMap(dist, 
v, e));
 
 1558  };
 1559 
 1560  
 1561  trimEdges.clear();
 1562  edgesToTrim.clear();
 1563  verticesOnTrimEdges.clear();
 1564  trimNewVertices.clear();
 1565 
 1566  auto cut_this_edge = [&](
const EntityHandle e, 
const double length, 
auto &ray,
 
 1567                           auto &ray_point) {
 1568    trimEdges.insert(e);
 1569    edgesToTrim[e].dIst = 1;
 1570    edgesToTrim[e].lEngth = length;
 1571    edgesToTrim[e].unitRayDir.resize(3, false);
 1572    edgesToTrim[e].rayPoint.resize(3, false);
 1573    for (int ii = 0; ii != 3; ++ii)
 1574      edgesToTrim[e].unitRayDir[ii] = ray(ii);
 1575    for (int ii = 0; ii != 3; ++ii)
 1576      edgesToTrim[e].rayPoint[ii] = ray_point(ii);
 1577  };
 1578 
 1582 
 1583  int num_nodes;
 1584 
 1585  auto get_tag_data = [&](
auto th, 
auto conn) {
 
 1587    CHKERR moab.tag_get_data(th, &conn, 1, &
t(0));
 
 1589  };
 1590 
 1591  double max_edge_length = 0;
 1592 
 1593
 1594  if (!fRont.empty()) {
 1595    
 1596    treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
 1597        new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
 1598    CHKERR treeSurfPtr->build(cutNewSurfaces, rootSetSurf);
 
 1599 
 1600    for (auto s : surface_skin) {
 1601 
 1602      auto project_on_surface = [&](
auto &
t) {
 
 1604 
 1606        CHKERR treeSurfPtr->closest_to_location(&
t(0), rootSetSurf, &t_p(0),
 
 1607                                                facets_out);
 1608 
 1610        CHKERR m_field.getInterface<Tools>()->getTriNormal(facets_out,
 
 1611                                                           &t_normal(0));
 1612        t_normal(
i) /= sqrt(t_normal(
i) * t_normal(
i));
 
 1613        const double dot = t_normal(
i) * (t_p(
i) - 
t(
i));
 
 1614        t_p(
i) = 
t(
i) + dot * t_normal(
i);
 
 1615 
 1616        return t_p;
 1617      };
 1618 
 1620      CHKERR moab.get_connectivity(s, conn, num_nodes, 
true);
 
 1621 
 1623      CHKERR moab.get_coords(conn, num_nodes, &coords_front[0]);
 
 1624 
 1626                                         &coords_front[2]);
 1628                                         &coords_front[5]);
 1629 
 1630      auto t_p_f0 = project_on_surface(t_f0);
 1631      auto t_p_f1 = project_on_surface(t_f1);
 1632 
 1633      CHKERR moab.set_coords(&conn[0], 1, &t_p_f0(0));
 
 1634      CHKERR moab.set_coords(&conn[1], 1, &t_p_f1(0));
 
 1635    }
 1636  }
 1637 
 1639    CHKERR SaveData(moab)(
"surface_skin_to_trim.vtk", surface_skin);
 
 1640 
 1642  Tag th_dist_front_vec;
 
 1643  CHKERR moab.tag_get_handle(
"DIST_FRONT_VECTOR", th_dist_front_vec);
 
 1644 
 1645  
 1646  for (auto e : cut_surface_edges) {
 1647 
 1648    auto get_conn_edge = [&moab](auto e) {
 1649      
 1651      int num_nodes;
 1652      CHKERR moab.get_connectivity(e, conn_edge, num_nodes, 
true);
 
 1653      return conn_edge;
 1654    };
 1655 
 1656    auto get_coords_edge = [&moab](auto conn_edge) {
 1657      std::array<double, 6> coords_edge;
 1658      CHKERR moab.get_coords(conn_edge, 2, coords_edge.data());
 
 1659      return coords_edge;
 1660    };
 1661 
 1662    auto get_ftensor_coords = [](const double *ptr) {
 1664    };
 1665 
 1666    auto conn_edge = get_conn_edge(e);
 1667    auto coords_edge = get_coords_edge(conn_edge);
 1668 
 1669    auto t_dist0 = get_tag_data(th_dist_front_vec, conn_edge[0]);
 1670    auto t_dist1 = get_tag_data(th_dist_front_vec, conn_edge[1]);
 1671    auto t_e0 = get_ftensor_coords(&coords_edge[0]);
 1672    auto t_e1 = get_ftensor_coords(&coords_edge[3]);
 1673 
 1675    t_edge_delta(
i) = t_e1(
i) - t_e0(
i);
 
 1676    const double edge_length2 = t_edge_delta(
i) * t_edge_delta(
i);
 
 1677    const double edge_length = sqrt(edge_length2);
 1678    if (edge_length == 0)
 1680 
 1681    auto add_edge = [&](auto t_cut) {
 1683      if (t_cut < 0.5) {
 1684        t_ray(
i) = t_cut * t_edge_delta(
i);
 
 1685        add_vert(conn_edge[0], e, fabs(t_cut));
 1686        add_vert(conn_edge[1], e, fabs(t_cut - 1));
 1687        cut_this_edge(e, edge_length, t_ray, t_e0);
 1688      } else {
 1690        t_edge_point(
i) = t_e0(
i) + t_cut * t_edge_delta(
i);
 
 1691        t_ray(
i) = t_edge_point(
i) - t_e1(
i);
 
 1692        add_vert(conn_edge[0], e, fabs(t_cut));
 1693        add_vert(conn_edge[1], e, fabs(t_cut - 1));
 1694        cut_this_edge(e, edge_length, t_ray, t_e1);
 1695      }
 1696      max_edge_length = std::max(max_edge_length, edge_length);
 1697    };
 1698 
 1699    auto trim_cross_edges = [&]() {
 1700      if (std::min(t_dist0(
i) * t_dist0(
i), t_dist1(
i) * t_dist1(
i)) <
 
 1701          edge_length2) {
 1702 
 1703        auto make_pair = [&](
const double d, 
const double te) {
 
 1704          return std::make_pair(d, te);
 1705        };
 1706 
 1707        auto min_pair = make_pair(-1, -1);
 1708 
 1709        for (auto f : surface_skin) {
 1710 
 1711          auto conn_front = get_conn_edge(f);
 1712          auto coords_front = get_coords_edge(conn_front);
 1713          auto t_f0 = get_ftensor_coords(&coords_front[0]);
 1714          auto t_f1 = get_ftensor_coords(&coords_front[3]);
 1715 
 1716          double te;
 1717          double tf;
 1718 
 1719          auto res = Tools::minDistanceFromSegments(
 1720              &t_e0(0), &t_e1(0), &t_f0(0), &t_f1(0), &te, &tf);
 1721 
 1722          if (res != Tools::NO_SOLUTION) {
 1723 
 1724            auto check = [](
auto v) {
 
 1725              return (
v > -std::numeric_limits<double>::epsilon() &&
 
 1726                      (
v - 1) < std::numeric_limits<double>::epsilon());
 
 1727            };
 1728 
 1729            if (check(te) && check(tf)) {
 1730 
 1732              t_delta(
i) = (t_e0(
i) + te * t_edge_delta(
i)) -
 
 1733                           (t_f0(0) + tf * (t_f1(
i) - t_f0(
i)));
 
 1734 
 1735              t_delta(
i) /= edge_length;
 
 1736 
 1737              const double dot = t_delta(
i) * t_delta(
i);
 
 1738 
 1739              if (min_pair.first < 0 || min_pair.first > dot) {
 1740 
 1741                if (dot < tol_trim_close * tol_trim_close)
 1742                  min_pair = make_pair(dot, te);
 1743              }
 1744            }
 1745          }
 1746        }
 1747 
 1748        if (min_pair.first > -std::numeric_limits<double>::epsilon()) {
 1749          add_edge(min_pair.second);
 1750          return true;
 1751        }
 1752      }
 1753 
 1754      return false;
 1755    };
 1756 
 1757    if (!trim_cross_edges()) {
 1758 
 1759      const double dot = t_dist0(
i) * t_dist1(
i);
 
 1760      const double dot_direction0 = t_dist0(
i) * t_edge_delta(
i);
 
 1761      const double dot_direction1 = t_dist1(
i) * t_edge_delta(
i);
 
 1762 
 1763      if (dot < std::numeric_limits<double>::epsilon() &&
 1764          dot_direction0 > -std::numeric_limits<double>::epsilon() &&
 1765          dot_direction1 < std::numeric_limits<double>::epsilon()) {
 1766 
 1767        const double s0 = t_dist0(
i) * t_edge_delta(
i) / edge_length;
 
 1768        const double s1 = t_dist1(
i) * t_edge_delta(
i) / edge_length;
 
 1769        const double t_cut = s0 / (s0 - s1);
 1770 
 1771        add_edge(t_cut);
 1772      }
 1773    }
 1774  }
 1775 
 1777    CHKERR SaveData(moab)(
"edges_potentially_to_trim.vtk", cut_surface_edges);
 
 1778 
 1780    CHKERR SaveData(moab)(
"edges_selected_to_trim.vtk", trimEdges);
 
 1781 
 1784    double q0 = 1;
 1785    CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0, th);
 
 1786 
 1787    double min_q = 1;
 1788    for (
auto t : adj_tets) {
 
 1789      int num_nodes;
 1791      CHKERR m_field.get_moab().get_connectivity(
t, conn, num_nodes, 
true);
 
 1793      if (th)
 1794        CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
 
 1795      else
 1796        CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
 
 1797 
 1798      for (
int n = 0; 
n != 4; ++
n) {
 
 1801          noalias(n_coords) = new_pos;
 1802        } else {
 1803          auto m = verticesOnTrimEdges.find(conn[
n]);
 
 1804          if (
m != verticesOnTrimEdges.end())
 
 1805            noalias(n_coords) =
 1806                m->second.rayPoint + 
m->second.dIst * 
m->second.unitRayDir;
 
 1807        }
 1808      }
 1809 
 1810      const double q = Tools::volumeLengthQuality(&coords[0]);
 1811      if (!std::isnormal(q)) {
 1812        min_q = -2;
 1813        break;
 1814      }
 1815      min_q = std::min(min_q, q);
 1816    }
 1817    return min_q / q0;
 1818  };
 1819 
 1820  Range checked_verts;
 
 1822    if (!(trimNewVertices.find(
v) != trimNewVertices.end())) {
 
 1823      auto &
r = edgesToTrim.at(e);
 
 1828      CHKERR moab.get_adjacencies(&
v, 1, 3, 
false, adj_tets);
 
 1829      adj_tets = intersect(adj_tets, cutNewVolumes);
 1830      double q = get_quality_change(adj_tets, 
v, new_pos);
 
 1831      if (q > projectEntitiesQualityTrashold) {
 1833        double dist = norm_2(ray_dir);
 1834        verticesOnTrimEdges[
v].dIst = 1;
 
 1835        verticesOnTrimEdges[
v].lEngth = dist;
 
 1836        verticesOnTrimEdges[
v].unitRayDir = ray_dir;
 
 1837        verticesOnTrimEdges[
v].rayPoint = ray_point;
 
 1838        trimNewVertices.insert(
v);
 
 1839      }
 1840      checked_verts.insert(
v);
 
 1841    }
 1842  };
 1843 
 1845    if (!(trimNewVertices.find(
v) != trimNewVertices.end())) {
 
 1846      auto &
m = edgesToTrim.at(e);
 
 1847      verticesOnTrimEdges[
v] = 
m;
 
 1848      verticesOnTrimEdges[
v].rayPoint = get_point_coords(
v);
 
 1849      verticesOnTrimEdges[
v].lEngth = 0;
 
 1850      verticesOnTrimEdges[
v].dIst = 0;
 
 1851      trimNewVertices.insert(
v);
 
 1852      checked_verts.insert(
v);
 
 1853    }
 1854  };
 1855 
 1856  
 1857  
 1858 
 1859  
 1860  auto hi = verts_map.get<0>().upper_bound(tol_trim_close);
 1861  verts_map.get<0>().erase(hi, verts_map.end());
 1862 
 1863  auto remove_verts = [&](
Range nodes) {
 
 1864    for (
auto n : nodes) {
 
 1865      auto r = verts_map.get<1>().equal_range(
n);
 
 1866      verts_map.get<1>().erase(
r.first, 
r.second);
 
 1867    }
 1868  };
 1869 
 1870  auto trim_verts = [&](
const Range verts, 
const bool move) {
 
 1871    VertMapMultiIndex verts_map_tmp;
 1872    for (auto p = corners.pair_begin(); p != corners.pair_end(); ++p) {
 1873      auto lo = verts_map.get<1>().lower_bound(p->first);
 1874      auto hi = verts_map.get<1>().upper_bound(p->second);
 1875      verts_map_tmp.insert(lo, hi);
 1876    }
 1877    if (move) {
 1878      for (
auto &
m : verts_map_tmp.get<0>())
 
 1879        add_trim_vert(
m.
v, 
m.e);
 
 1880    } else {
 1881      for (
auto &
m : verts_map_tmp.get<0>())
 
 1882        add_no_move_trim(
m.
v, 
m.e);
 
 1883    }
 1884  };
 1885 
 1886  auto trim_edges = [&](
const Range ents, 
const bool move) {
 
 1887    VertMapMultiIndex verts_map_tmp;
 1888    for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
 1889      auto lo = verts_map.get<2>().lower_bound(p->first);
 1890      auto hi = verts_map.get<2>().upper_bound(p->second);
 1891      verts_map_tmp.insert(lo, hi);
 1892      verts_map.get<2>().erase(lo, hi);
 1893    }
 1894    if (move) {
 1895      for (
auto &
m : verts_map_tmp.get<0>())
 
 1896        add_trim_vert(
m.
v, 
m.e);
 
 1897    } else {
 1898      for (
auto &
m : verts_map_tmp.get<0>())
 
 1899        add_no_move_trim(
m.
v, 
m.e);
 
 1900    }
 1901  };
 1902 
 1903  auto intersect_self = [&](
Range &
a, 
const Range b) { 
a = intersect(
a, b); };
 
 1904 
 1905  map<std::string, Range> range_maps;
 1906 
 1907  
 1908  CHKERR skin.find_skin(0, cutNewSurfaces, 
false, range_maps[
"surface_skin"]);
 
 1909  
 1910  intersect_self(range_maps["surface_skin"], trimEdges);
 1911  
 1912  range_maps["fixed_edges_on_surface_skin"] =
 1913      intersect(range_maps["surface_skin"], fix_edges);
 1914 
 1915  
 1916  CHKERR moab.get_adjacencies(fixed_edges_verts, 1, 
false,
 
 1917                              range_maps["fixed_edges_verts_edges"],
 1918                              moab::Interface::UNION);
 1919  intersect_self(range_maps["fixed_edges_verts_edges"], trimEdges);
 1920 
 1921  CHKERR moab.get_connectivity(
 
 1922      range_maps["fixed_edges_verts_edges"],
 1923      range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
 1924  intersect_self(range_maps["fixed_edges_verts_edges_verts_on_the_skin"],
 1925                 tets_skin_verts);
 1926 
 1927  
 1928  trim_verts(corners, false);
 1929  remove_verts(corners);
 1930 
 1931  
 1932  trim_edges(range_maps["fixed_edges_on_surface_skin"], true);
 1933  remove_verts(range_maps["fixed_edges_on_surface_skin_verts"]);
 1934 
 1935  
 1936  trim_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
 1937  remove_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"]);
 1938 
 1939  
 1940  trim_edges(range_maps["surface_skin"], true);
 1941  trim_verts(tets_skin_verts, false);
 1942  remove_verts(tets_skin_verts);
 1943 
 1944  for (
auto &
m : verts_map.get<0>())
 
 1945    add_trim_vert(
m.
v, 
m.e);
 
 1946 
 1947  for (
auto m : verticesOnTrimEdges) {
 
 1950    CHKERR moab.get_adjacencies(&
v, 1, 1, 
false, adj);
 
 1951    for (auto e : adj) {
 1952      edgesToTrim.erase(e);
 1953      trimEdges.erase(e);
 1954    }
 1955  }
 1956 
 1958    if (!trimNewVertices.empty())
 1959      CHKERR SaveData(moab)(
"trim_close_vertices.vtk", trimNewVertices);
 
 1960 
 1962    if (!trimEdges.empty())
 1963      CHKERR SaveData(moab)(
"trim_edges.vtk", trimEdges);
 
 1964 
 1966}
 1967 
 1970  CoreInterface &m_field = cOre;
 1971  moab::Interface &moab = m_field.get_moab();
 1972  MeshRefinement *refiner;
 1973  auto ref_ents_ptr = m_field.get_ref_ents();
 1975 
 1976  CHKERR m_field.getInterface(refiner);
 
 1977  CHKERR refiner->addVerticesInTheMiddleOfEdges(trimEdges, 
bit);
 
 1979 
 1980  trimNewVolumes.clear();
 1981  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 1983 
 1984  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
 1985       mit != edgesToTrim.end(); mit++) {
 1986    auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
 1987        boost::make_tuple(MBVERTEX, mit->first));
 1988    if (vit ==
 1989        ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end())
 1991              "No vertex on trim edges, that make no sense");
 1992 
 1993    const boost::shared_ptr<RefEntity> &ref_ent = *vit;
 1994    if ((ref_ent->getBitRefLevel() & 
bit).any()) {
 
 1996      trimNewVertices.insert(vert);
 1997      verticesOnTrimEdges[vert] = mit->second;
 1998    }
 1999  }
 2000 
 2001  
 2002  trimNewSurfaces.clear();
 2003  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 2004      bit, 
bit, MBTRI, trimNewSurfaces);
 
 2005 
 2006  Range trim_new_surfaces_nodes;
 
 2007  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, 
true);
 
 2008  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
 2009  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
 2010  Range faces_not_on_surface;
 
 2011  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, 
false,
 
 2012                              faces_not_on_surface, moab::Interface::UNION);
 2013  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
 2014 
 2015  
 2016  Range all_surfaces_on_bit_level;
 
 2017  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 2019  all_surfaces_on_bit_level =
 2020      intersect(all_surfaces_on_bit_level, cutNewSurfaces);
 2021 
 2022  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
 2023 
 2025}
 2026 
 2028                                             Range *corner_nodes,
 
 2029                                             const double tol, 
Tag th,
 
 2031  CoreInterface &m_field = cOre;
 2032  moab::Interface &moab = m_field.get_moab();
 2033  Skinner skin(&moab);
 2034 
 2035  auto get_adj = [&moab](
const Range r, 
const int dim) {
 
 2037    if (dim)
 2038      CHKERR moab.get_adjacencies(r, dim, 
false, 
a, moab::Interface::UNION);
 
 2039    else
 2040      CHKERR moab.get_connectivity(r, 
a, 
true);
 
 2042  };
 2043 
 2044  auto get_skin = [&skin](
const auto r) {
 
 2046    CHKERR skin.find_skin(0, r, 
false, s);
 
 2047    return s;
 2048  };
 2049 
 2051 
 2052  auto trim_tets_skin = get_skin(trimNewVolumes);
 2053  auto trim_tets_skin_edges = get_adj(trim_tets_skin, 1);
 2054  auto trim_surface_edges = get_adj(trimNewSurfaces, 1);
 2055 
 2056  auto contarain_edges =
 2057      intersect(trim_surface_edges, get_adj(constrainSurface, 1));
 2058  if (fixed_edges)
 2059    contarain_edges.merge(
 2060        intersect(fixed_edges->subset_by_type(MBEDGE), trim_surface_edges));
 2061 
 2062  Range barrier_vertices(trimNewVertices);
 
 2063  if (corner_nodes) {
 2064    
 2065    barrier_vertices.merge(get_adj(
 2066        intersect(get_adj(intersect(*corner_nodes, barrier_vertices), 2),
 2067                  trimNewSurfaces),
 2068        0));
 2069  }
 2070 
 2071  auto get_nodes_with_one_node_on_fixed_edge_other_not = [&]() {
 2072    const auto fixed_edges_on_surface =
 2073        intersect(*fixed_edges, trim_surface_edges);
 2074    const auto skin_fixed_edges_on_surface = get_skin(fixed_edges_on_surface);
 2075    const auto barrier_nodes = subtract(skin_fixed_edges_on_surface,
 2076                                  get_adj(get_skin(trimNewSurfaces), 0));
 2077    return barrier_nodes;
 2078  };
 2079 
 2080  if (fixed_edges)
 2081    barrier_vertices.merge(get_nodes_with_one_node_on_fixed_edge_other_not());
 2082 
 2083  auto add_close_surface_barrier = [&]() {
 2085 
 2086    
 2087    
 2088    
 2089 
 2091 
 2092    auto test_edges =
 2093        subtract(trim_surface_edges, get_adj(constrainSurface, 1));
 2094    if (fixed_edges)
 2095      test_edges = subtract(test_edges, *fixed_edges);
 2096    auto trim_surface_nodes = get_adj(test_edges, 0);
 2097    trim_surface_nodes = subtract(trim_surface_nodes, barrier_vertices);
 2098    if (fixed_edges)
 2099      trim_surface_nodes =
 2100          subtract(trim_surface_nodes, get_adj(*fixed_edges, 0));
 2101 
 2103    trim_skin.merge(contarain_edges);
 2104    if (fRont.empty())
 2105      trim_skin.merge(get_skin(sUrface));
 2106    else
 2107      trim_skin.merge(fRont);
 2108    if (fixed_edges)
 2109      trim_skin.merge(intersect(*fixed_edges, trim_surface_edges));
 2110 
 2111    if (
debug && !trim_skin.empty())
 
 2112      CHKERR SaveData(m_field.get_moab())(
"trim_skin.vtk", trim_skin);
 
 2113 
 2115    if (th)
 2116      CHKERR moab.tag_get_data(th, trim_surface_nodes, &*coords.begin());
 
 2117    else
 2118      CHKERR moab.get_coords(trim_surface_nodes, &*coords.begin());
 
 2119 
 2120    if (!trim_skin.empty()) {
 2121 
 2122      std::vector<double> min_distances(trim_surface_nodes.size(), -1);
 2123      CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
 
 2124          &*coords.begin(), trim_surface_nodes.size(), trim_skin,
 2125          &*min_distances.begin());
 2126 
 2127      auto coords_ptr = &*coords.begin();
 2128      auto min_dist = &*min_distances.begin();
 2129 
 2130      std::vector<EntityHandle> add_nodes;
 2131      add_nodes.reserve(trim_surface_nodes.size());
 2132 
 2133      for (
auto n : trim_surface_nodes) {
 
 2134 
 2135        if ((*min_dist) > std::numeric_limits<double>::epsilon()) {
 2137 
 2139          CHKERR treeSurfPtr->closest_to_location(coords_ptr, rootSetSurf,
 
 2140                                                  &point_out[0], facets_out);
 2141 
 2144          CHKERR m_field.getInterface<Tools>()->getTriNormal(facets_out,
 
 2145                                                             &normal[0]);
 2146 
 2147          normal /= norm_2(normal);
 2149 
 2150          double dist = norm_2(
delta);
 
 2151          if (dist < 
tol * (*min_dist))
 
 2152            add_nodes.emplace_back(
n);
 
 2153        }
 2154 
 2155        coords_ptr += 3;
 2156        min_dist += 1;
 2157      }
 2158 
 2159      barrier_vertices.insert_list(add_nodes.begin(), add_nodes.end());
 2160    }
 2161 
 2163  };
 2164 
 2165  auto remove_faces_on_skin = [&]() {
 2167    Range skin_faces = intersect(trimNewSurfaces, trim_tets_skin);
 
 2168    if (!skin_faces.empty()) {
 2169      barrier_vertices.merge(get_adj(skin_faces, 0));
 2170      for (auto f : skin_faces)
 2171        trimNewSurfaces.erase(
f);
 
 2172    }
 2174  };
 2175 
 2176  auto get_trim_free_edges = [&]() {
 2177    
 2178    Range trim_surf_skin;
 
 2179    CHKERR skin.find_skin(0, trimNewSurfaces, 
false, trim_surf_skin);
 
 2180    trim_surf_skin = subtract(trim_surf_skin, trim_tets_skin_edges);
 2181 
 2182    Range trim_surf_skin_verts;
 
 2183    CHKERR moab.get_connectivity(trim_surf_skin, trim_surf_skin_verts, 
true);
 
 2184    for (auto e : barrier_vertices)
 2185      trim_surf_skin_verts.erase(e);
 2186 
 2188    CHKERR moab.get_adjacencies(trim_surf_skin_verts, 1, 
false, free_edges,
 
 2189                                moab::Interface::UNION);
 2190    free_edges =
 2191        subtract(intersect(free_edges, trim_surf_skin), contarain_edges);
 2192 
 2193    return free_edges;
 2194  };
 2195 
 2196  CHKERR remove_faces_on_skin();
 
 2197  CHKERR add_close_surface_barrier();
 
 2198 
 2199  if (
debug && !barrier_vertices.empty())
 
 2200    CHKERR SaveData(m_field.get_moab())(
"barrier_vertices.vtk",
 
 2201                                        barrier_vertices);
 2202 
 2203  int nn = 0;
 2205  while (!(out_edges = get_trim_free_edges()).empty()) {
 2206 
 2207    if (
debug && !trimNewSurfaces.empty())
 
 2208      CHKERR SaveData(m_field.get_moab())(
 
 2209          "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
 2210          trimNewSurfaces);
 2211 
 2212    if (
debug && !out_edges.empty())
 
 2213      CHKERR SaveData(m_field.get_moab())(
 
 2214          "trimNewSurfacesOutsideVerts_" +
 2215              boost::lexical_cast<std::string>(nn) + ".vtk",
 2216          out_edges);
 2217 
 2218    Range outside_faces;
 
 2219    CHKERR moab.get_adjacencies(out_edges, 2, 
false, outside_faces,
 
 2220                                moab::Interface::UNION);
 2221    trimNewSurfaces = subtract(trimNewSurfaces, outside_faces);
 2222    ++nn;
 2223  }
 2224 
 2225  if (
debug && !trimNewSurfaces.empty())
 
 2226    CHKERR SaveData(m_field.get_moab())(
 
 2227        "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
 2228        trimNewSurfaces);
 2229 
 2231}
 2232 
 2234CutMeshInterface::removePathologicalFrontTris(const BitRefLevel split_bit,
 2236  CoreInterface &m_field = cOre;
 2237  moab::Interface &moab = m_field.get_moab();
 2238  PrismInterface *interface;
 2240  CHKERR m_field.getInterface(interface);
 
 2241  
 2242  {
 2245    CHKERR moab.create_meshset(MESHSET_SET, meshset);
 
 2246    CHKERR moab.add_entities(meshset, ents);
 
 2247    CHKERR interface->findFacesWithThreeNodesOnInternalSurfaceSkin(
 
 2248        meshset, split_bit, true, front_tris);
 2249    CHKERR moab.delete_entities(&meshset, 1);
 
 2250    ents = subtract(ents, front_tris);
 2251  }
 2252  
 2254  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 2256  
 2257  Skinner skin(&moab);
 2259  CHKERR skin.find_skin(0, tets, 
false, tets_skin);
 
 2260  ents = subtract(ents, tets_skin);
 2261 
 2263}
 2264 
 2265MoFEMErrorCode CutMeshInterface::splitSides(
const BitRefLevel split_bit,
 
 2266                                            const BitRefLevel 
bit,
 
 2268  CoreInterface &m_field = cOre;
 2269  moab::Interface &moab = m_field.get_moab();
 2270  PrismInterface *interface;
 2272  CHKERR m_field.getInterface(interface);
 
 2274  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
 
 2275  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 2276      split_bit, 
BitRefLevel().set(), MBTET, meshset_volume);
 
 2278  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
 
 2279  CHKERR moab.add_entities(meshset_surface, ents);
 
 2280  CHKERR interface->getSides(meshset_surface, split_bit, 
true);
 
 2281  CHKERR interface->splitSides(meshset_volume, 
bit, meshset_surface, 
true,
 
 2282                               true);
 2283  CHKERR moab.delete_entities(&meshset_volume, 1);
 
 2284  CHKERR moab.delete_entities(&meshset_surface, 1);
 
 2285  if (th) {
 2287    ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 2289    for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
 2290      int num_nodes;
 2292      CHKERR moab.get_connectivity(*pit, conn, num_nodes, 
true);
 
 2294      CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
 
 2295      CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
 
 2296    }
 2297  }
 2299}
 2300 
 2303    const Range &fixed_edges, 
const Range &corner_nodes, 
Range &edges_to_merge,
 
 2304    Range &out_tets, 
Range &new_surf, 
Tag th, 
const bool update_meshsets,
 
 2305    const BitRefLevel *bit_ptr, 
const bool debug) {
 
 2306  CoreInterface &m_field = cOre;
 2307  moab::Interface &moab = m_field.get_moab();
 2309 
 2310
 2311
 2312
 2313  struct MergeNodes {
 2314    CoreInterface &mField;
 2315    const bool onlyIfImproveQuality;
 2317    bool updateMehsets;
 2318 
 2319    MergeNodes(CoreInterface &m_field, const bool only_if_improve_quality,
 2320               Tag th, 
bool update_mehsets)
 
 2321        : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
 2322          tH(
th), updateMehsets(update_mehsets) {
 
 2323      mField.getInterface(nodeMergerPtr);
 2324    }
 2325    NodeMergerInterface *nodeMergerPtr;
 2328                              bool add_child = true) {
 2329      moab::Interface &moab = mField.get_moab();
 2333      CHKERR moab.get_adjacencies(conn, 2, 3, 
false, vert_tets,
 
 2334                                  moab::Interface::UNION);
 2335      vert_tets = intersect(vert_tets, proc_tets);
 2337      CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
 
 2338                                       onlyIfImproveQuality, 0, line_search,
 2339                                       tH);
 2340 
 2341      if (add_child && nodeMergerPtr->getSuccessMerge()) {
 2342 
 2343        Range::iterator lo, hi = proc_tets.begin();
 2344        for (auto pt = vert_tets.pair_begin(); pt != vert_tets.pair_end();
 2345             ++pt) {
 2346          lo = proc_tets.lower_bound(hi, proc_tets.end(), pt->first);
 2347          if (lo != proc_tets.end()) {
 2348            hi = proc_tets.upper_bound(lo, proc_tets.end(), pt->second);
 2349            proc_tets.erase(lo, hi);
 2350          } else
 2351            break;
 2352        }
 2353        proc_tets.merge(out_tets);
 2354 
 2355        auto &parent_child_map = nodeMergerPtr->getParentChildMap();
 2356 
 2357        struct ChangeChild {
 2359          ChangeChild(
const EntityHandle child) : child(child) {}
 
 2360          void operator()(NodeMergerInterface::ParentChild &p) {
 2361            p.cHild = child;
 2362          }
 2363        };
 2364 
 2365        std::vector<decltype(parentsChildMap.get<0>().begin())> it_vec;
 2366        it_vec.reserve(parentsChildMap.size());
 2367 
 2368        for (auto &p : parent_child_map) {
 2369 
 2370          it_vec.clear();
 2371          for (auto it = parentsChildMap.get<0>().equal_range(p.pArent);
 2372               it.first != it.second; ++it.first)
 2373            it_vec.emplace_back(it.first);
 2374 
 2375          for (auto it = parentsChildMap.get<1>().equal_range(p.pArent);
 2376               it.first != it.second; ++it.first)
 2377            it_vec.emplace_back(parentsChildMap.project<0>(it.first));
 2378 
 2379          for (auto &it : it_vec)
 2380            parentsChildMap.modify(it, ChangeChild(p.cHild));
 2381        }
 2382 
 2383        parentsChildMap.insert(parent_child_map.begin(),
 2384                               parent_child_map.end());
 2385      }
 2387    }
 2388 
 2390                                       Range ¬_merged_edges,
 
 2391                                       bool add_child) {
 2392      moab::Interface &moab = mField.get_moab();
 2394      if (add_child) {
 2395 
 2396        std::vector<EntityHandle> parents_ents_vec(parentsChildMap.size());
 2397        for (auto &it : parentsChildMap)
 2398          parents_ents_vec.emplace_back(it.pArent);
 2400        parent_ents.insert_list(parents_ents_vec.begin(),
 2401                                parents_ents_vec.end());
 2402 
 2403        Range surf_parent_ents = intersect(new_surf, parent_ents);
 
 2404        new_surf = subtract(new_surf, surf_parent_ents);
 2405        Range child_surf_ents;
 
 2406        CHKERR updateRangeByChilds(parentsChildMap, surf_parent_ents,
 
 2407                                   child_surf_ents);
 2408        new_surf.merge(child_surf_ents);
 2409 
 2410        Range edges_to_merge_parent_ents =
 
 2411            intersect(edges_to_merge, parent_ents);
 2412        edges_to_merge = subtract(edges_to_merge, edges_to_merge_parent_ents);
 2413        Range merged_child_edge_ents;
 
 2414        CHKERR updateRangeByChilds(parentsChildMap, edges_to_merge_parent_ents,
 
 2415                                   merged_child_edge_ents);
 2416 
 2417        Range not_merged_edges_child_ents =
 
 2418            intersect(not_merged_edges, parent_ents);
 2419        not_merged_edges =
 2420            subtract(not_merged_edges, not_merged_edges_child_ents);
 2421        Range not_merged_child_edge_ents;
 
 2422        CHKERR updateRangeByChilds(parentsChildMap, not_merged_edges_child_ents,
 
 2423                                   not_merged_child_edge_ents);
 2424 
 2425        merged_child_edge_ents =
 2426            subtract(merged_child_edge_ents, not_merged_child_edge_ents);
 2427        edges_to_merge.merge(merged_child_edge_ents);
 2428        not_merged_edges.merge(not_merged_child_edge_ents);
 2429 
 2430        if (updateMehsets) {
 2432                   (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
 2435            CHKERR moab.get_entities_by_handle(cubit_meshset, meshset_ents,
 
 2436                                               true);
 2438            CHKERR updateRangeByChilds(parentsChildMap, meshset_ents,
 
 2439                                       child_ents);
 2440            CHKERR moab.add_entities(cubit_meshset, child_ents);
 
 2441          }
 2442        }
 2443      }
 2444 
 2446    };
 2447 
 2448  private:
 2449    NodeMergerInterface::ParentChildMap parentsChildMap;
 2450    std::vector<EntityHandle> childsVec;
 2451 
 2453        const NodeMergerInterface::ParentChildMap &parent_child_map,
 2456      childsVec.clear();
 2457      childsVec.reserve(parents.size());
 2458      for (auto pit = parents.pair_begin(); pit != parents.pair_end(); pit++) {
 2459        auto it = parent_child_map.lower_bound(pit->first);
 2460        if (it != parent_child_map.end()) {
 2461          for (auto hi_it = parent_child_map.upper_bound(pit->second);
 2462               it != hi_it; ++it)
 2463            childsVec.emplace_back(it->cHild);
 2464        }
 2465      }
 2466      childs.insert_list(childsVec.begin(), childsVec.end());
 2468    }
 2469  };
 2470 
 2471
 2472
 2473
 2474  struct LengthMap {
 2476    CoreInterface &mField;
 2477    moab::Interface &moab;
 2478    const double maxLength;
 2479    LengthMap(CoreInterface &m_field, 
Tag th, 
double max_length)
 
 2480        : tH(
th), mField(m_field), moab(m_field.get_moab()),
 
 2481          maxLength(max_length) {
 2482      gammaL = 1.;
 2483      gammaQ = 3.;
 2484      acceptedThrasholdMergeQuality = 0.5;
 2485    }
 2486 
 2487    double
 2488        gammaL; 
 2489    double
 2490        gammaQ; 
 2491    double acceptedThrasholdMergeQuality; 
 2492
 2493 
 2495                              LengthMapData_multi_index &length_map,
 2496                              double &ave) const {
 2497      int num_nodes;
 2499      std::array<double, 12> coords;
 2501      VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
 
 2502      VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
 
 2504 
 2505      struct NodeQuality {
 2507        double quality;
 2508        NodeQuality(
const EntityHandle node) : node(node), quality(1) {}
 
 2509      };
 2510 
 2511      typedef multi_index_container<
 2512          NodeQuality, indexed_by<
 2513 
 2514                           sequenced<>,
 2515 
 2516                           hashed_non_unique<tag<Ent_mi_tag>,
 2518                                                    &NodeQuality::node>>
 2519 
 2520                           >>
 2521          NodeQuality_sequence;
 2522 
 2523      NodeQuality_sequence node_quality_sequence;
 2524 
 2526      CHKERR moab.get_connectivity(tets, edges_nodes, 
false);
 
 2528      CHKERR moab.get_adjacencies(edges, 3, 
false, edges_tets,
 
 2529                                  moab::Interface::UNION);
 2530      edges_tets = intersect(edges_tets, tets);
 2531 
 2532      for (auto node : edges_nodes)
 2533        node_quality_sequence.emplace_back(node);
 2534 
 2535      for (auto tet : edges_tets) {
 2536 
 2537        CHKERR moab.get_connectivity(tet, conn, num_nodes, 
true);
 
 2538        if (tH)
 2539          CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
 
 2540        else
 2541          CHKERR moab.get_coords(conn, num_nodes, coords.data());
 
 2542 
 2543        const double q = Tools::volumeLengthQuality(coords.data());
 2544 
 2545        for (
auto n : {0, 1, 2, 3}) {
 
 2546          auto it = node_quality_sequence.get<1>().find(conn[
n]);
 
 2547          if (it != node_quality_sequence.get<1>().end())
 2548            const_cast<double &>(it->quality) = std::min(q, it->quality);
 2549        }
 2550      }
 2551 
 2552      for (auto edge : edges) {
 2553 
 2554        CHKERR moab.get_connectivity(edge, conn, num_nodes, 
true);
 
 2555 
 2556        if (tH)
 2557          CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
 
 2558        else
 2559          CHKERR moab.get_coords(conn, num_nodes, coords.data());
 
 2560 
 2561        double q = 1;
 2562        for (
auto n : {0, 1}) {
 
 2563          auto it = node_quality_sequence.get<1>().find(conn[
n]);
 
 2564          if (it != node_quality_sequence.get<1>().end())
 2565            q = std::min(q, it->quality);
 2566        }
 2567 
 2568        if (q < acceptedThrasholdMergeQuality) {
 2569          noalias(
delta) = (s0 - s1) / maxLength;
 
 2571          double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
 2572          length_map.insert(LengthMapData(val, q, edge));
 2573        }
 2574      }
 2575 
 2576      ave = 0;
 2577      for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
 2578               length_map.get<0>().begin();
 2579           mit != length_map.get<0>().end(); mit++) {
 2580        ave += mit->qUality;
 2581      }
 2582      ave /= length_map.size();
 2584    }
 2585  };
 2586 
 2587
 2588
 2589
 2590  struct Toplogy {
 2591 
 2592    CoreInterface &mField;
 2594    Toplogy(CoreInterface &m_field, 
Tag th) : mField(m_field), tH(
th) {}
 
 2595 
 2596    enum TYPE {
 2597      FREE = 0,
 2598      SKIN = 1 << 0,
 2599      SURFACE = 1 << 1,
 2600      SURFACE_SKIN = 1 << 2,
 2601      FRONT_ENDS = 1 << 3,
 2602      FIX_EDGES = 1 << 4,
 2603      FIX_CORNERS = 1 << 5
 2604    };
 2605 
 2606    typedef map<int, Range> SetsMap;
 2607 
 2609                                 const Range &fixed_edges,
 
 2610                                 const Range &corner_nodes,
 
 2611                                 const Range &constrain_surface,
 
 2612                                 SetsMap &sets_map) const {
 2613      moab::Interface &moab(mField.get_moab());
 2614      Skinner skin(&moab);
 2616 
 2617      sets_map[FIX_CORNERS].merge(corner_nodes);
 2619      CHKERR moab.get_connectivity(fixed_edges, fixed_verts, 
true);
 
 2620      sets_map[FIX_EDGES].swap(fixed_verts);
 2621 
 2623      CHKERR skin.find_skin(0, tets, 
false, tets_skin);
 
 2624      Range tets_skin_edges;
 
 2625      CHKERR moab.get_adjacencies(tets_skin, 1, 
false, tets_skin_edges,
 
 2626                                  moab::Interface::UNION);
 2627 
 2628      
 2629      Range constrain_surface_verts;
 
 2630      CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
 
 2631                                   true);
 2632      Range constrain_surface_edges;
 
 2633      CHKERR moab.get_adjacencies(constrain_surface, 1, 
false,
 
 2634                                  constrain_surface_edges,
 2635                                  moab::Interface::UNION);
 2636 
 2637      
 2640      Range front_in_the_body;
 
 2641      front_in_the_body = subtract(surface_skin, tets_skin_edges);
 2642      Range front_in_the_body_verts;
 
 2643      CHKERR moab.get_connectivity(front_in_the_body, front_in_the_body_verts,
 
 2644                                   true);
 2646      CHKERR skin.find_skin(0, front_in_the_body, 
false, front_ends);
 
 2647      front_ends.merge(
 2648          intersect(front_in_the_body_verts, constrain_surface_verts));
 2649      sets_map[FRONT_ENDS].swap(front_ends);
 2650 
 2651      Range surface_skin_verts;
 
 2652      CHKERR moab.get_connectivity(surface_skin, surface_skin_verts, 
true);
 
 2653      sets_map[SURFACE_SKIN].swap(surface_skin_verts);
 2654 
 2655      
 2656      Range surface_verts;
 
 2658      sets_map[SURFACE_SKIN].merge(
 2659          intersect(constrain_surface_verts, surface_verts));
 2660      sets_map[SURFACE].swap(surface_verts);
 2661 
 2662      
 2663      Range tets_skin_verts;
 
 2664      CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, 
true);
 
 2665      sets_map[SKIN].swap(tets_skin_verts);
 2666      sets_map[SKIN].merge(constrain_surface_verts);
 2667 
 2669      CHKERR moab.get_connectivity(tets, tets_verts, 
true);
 
 2670      sets_map[FREE].swap(tets_verts);
 2671 
 2673    }
 2674 
 2676                               Range &proc_tets)
 const {
 
 2677      moab::Interface &moab(mField.get_moab());
 2679      Range edges_to_merge_verts;
 
 2680      CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts, 
true);
 
 2681      Range edges_to_merge_verts_tets;
 
 2682      CHKERR moab.get_adjacencies(edges_to_merge_verts, 3, 
false,
 
 2683                                  edges_to_merge_verts_tets,
 2684                                  moab::Interface::UNION);
 2685      edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
 2686      proc_tets.swap(edges_to_merge_verts_tets);
 2688    }
 2689 
 2691                                  const Range &fixed_edges,
 
 2692                                  const Range &corner_nodes,
 
 2693                                  const Range &constrain_surface,
 
 2694                                  Range &edges_to_merge,
 
 2695                                  Range ¬_merged_edges) {
 
 2696      moab::Interface &moab(mField.get_moab());
 2698 
 2699      
 2700      Skinner skin(&moab);
 2702      CHKERR skin.find_skin(0, tets, 
false, tets_skin);
 
 2705 
 2706      
 2707      Range constrain_surface_verts;
 
 2708      CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
 
 2709                                   true);
 2710      Range constrain_surface_edges;
 
 2711      CHKERR moab.get_adjacencies(constrain_surface, 1, 
false,
 
 2712                                  constrain_surface_edges,
 2713                                  moab::Interface::UNION);
 2714 
 2715      
 2716      Range tets_skin_edges;
 
 2717      CHKERR moab.get_adjacencies(tets_skin, 1, 
false, tets_skin_edges,
 
 2718                                  moab::Interface::UNION);
 2719 
 2720      Range surface_front;
 
 2721      surface_front = subtract(surface_skin, tets_skin_edges);
 2722      Range surface_front_nodes;
 
 2723      CHKERR moab.get_connectivity(surface_front, surface_front_nodes, 
true);
 
 2724 
 2726      CHKERR skin.find_skin(0, surface_front, 
false, ends_nodes);
 
 2727      ends_nodes.merge(intersect(surface_front_nodes, constrain_surface_verts));
 2728 
 2729      
 2730      surface_skin.merge(constrain_surface);
 2731      tets_skin_edges.merge(constrain_surface_edges);
 2732 
 2733      
 2734      Range surface_edges;
 
 2735      CHKERR moab.get_adjacencies(
surface, 1, 
false, surface_edges,
 
 2736                                  moab::Interface::UNION);
 2737      
 2738      Range surface_edges_verts;
 
 2739      CHKERR moab.get_connectivity(surface_edges, surface_edges_verts, 
true);
 
 2740      
 2741      Range tets_skin_edges_verts;
 
 2742      CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
 
 2743                                   true);
 2744 
 2745      Range edges_to_remove;
 
 2746 
 2747      
 2748      {
 2749        Range ents_nodes_and_edges;
 
 2750        ents_nodes_and_edges.merge(tets_skin_edges_verts);
 2751        ents_nodes_and_edges.merge(tets_skin_edges);
 2752        CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
 
 2753                                        false);
 2754      }
 2755      edges_to_merge = subtract(edges_to_merge, edges_to_remove);
 2756      not_merged_edges.merge(edges_to_remove);
 2757 
 2758      
 2759      {
 2760        Range ents_nodes_and_edges;
 
 2761        ents_nodes_and_edges.merge(surface_edges_verts);
 2762        ents_nodes_and_edges.merge(surface_edges);
 2763        ents_nodes_and_edges.merge(tets_skin_edges_verts);
 2764        ents_nodes_and_edges.merge(tets_skin_edges);
 2765        CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
 
 2766                                        false);
 2767      }
 2768      edges_to_merge = subtract(edges_to_merge, edges_to_remove);
 2769      not_merged_edges.merge(edges_to_remove);
 2770 
 2771      
 2772      Range fixed_edges_nodes;
 
 2773      CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes, 
true);
 
 2774      {
 2775        Range ents_nodes_and_edges;
 
 2776        ents_nodes_and_edges.merge(fixed_edges_nodes);
 2777        ents_nodes_and_edges.merge(ends_nodes);
 2778        ents_nodes_and_edges.merge(corner_nodes);
 2779        ents_nodes_and_edges.merge(fixed_edges);
 2780        CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
 
 2781                                        false);
 2782      }
 2783      edges_to_merge = subtract(edges_to_merge, edges_to_remove);
 2784      not_merged_edges.merge(edges_to_remove);
 2785 
 2786      
 2787      CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove, 
false);
 
 2788      edges_to_merge = subtract(edges_to_merge, edges_to_remove);
 2789      not_merged_edges.merge(edges_to_remove);
 2790 
 2791      
 2792      {
 2793        Range ents_nodes_and_edges;
 
 2794        ents_nodes_and_edges.merge(surface_skin);
 2795        ents_nodes_and_edges.merge(fixed_edges_nodes);
 2796        CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
 
 2797                                        false);
 2798      }
 2799      edges_to_merge = subtract(edges_to_merge, edges_to_remove);
 2800      not_merged_edges.merge(edges_to_remove);
 2801 
 2802      
 2803      {
 2804        Range ents_nodes_and_edges;
 
 2805        ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
 2806        ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
 2807        CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
 
 2808                                        false);
 2809      }
 2810      edges_to_merge = subtract(edges_to_merge, edges_to_remove);
 2811      not_merged_edges.merge(edges_to_remove);
 2812 
 2813      
 2814      {
 2815        Range ents_nodes_and_edges;
 
 2816        ents_nodes_and_edges.merge(surface_front_nodes);
 2817        ents_nodes_and_edges.merge(surface_front);
 2818        ents_nodes_and_edges.merge(tets_skin_edges_verts);
 2819        ents_nodes_and_edges.merge(tets_skin_edges);
 2820        CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
 
 2821                                        false);
 2822      }
 2823      edges_to_merge = subtract(edges_to_merge, edges_to_remove);
 2824      not_merged_edges.merge(edges_to_remove);
 2825 
 2827    }
 2828 
 2829  private:
 2831                                            Range &edges_to_remove,
 
 2833      moab::Interface &moab(mField.get_moab());
 2835      
 2836      Range ents_nodes = ents.subset_by_type(MBVERTEX);
 
 2837      if (ents_nodes.empty()) {
 2838        CHKERR moab.get_connectivity(ents, ents_nodes, 
true);
 
 2839      }
 2840      
 2841      Range ents_nodes_edges;
 
 2842      CHKERR moab.get_adjacencies(ents_nodes, 1, 
false, ents_nodes_edges,
 
 2843                                  moab::Interface::UNION);
 2844      
 2845      Range ents_nodes_edges_nodes;
 
 2846      CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
 
 2847                                   true);
 2848      
 2849      ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
 2850      Range ents_nodes_edges_nodes_edges;
 
 2851      CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1, 
false,
 
 2852                                  ents_nodes_edges_nodes_edges,
 2853                                  moab::Interface::UNION);
 2854      
 2855      ents_nodes_edges =
 2856          subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
 2857      ents_nodes_edges =
 2858          subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
 2859 
 2860      edges_to_remove.swap(ents_nodes_edges);
 2862        CHKERR SaveData(moab)(
"edges_to_remove.vtk", edges_to_remove);
 
 2863      }
 2865    }
 2866  };
 2867 
 2868  Range not_merged_edges;
 
 2869  CHKERR Toplogy(m_field, th)
 
 2870      .removeBadEdges(
surface, tets, fixed_edges, corner_nodes,
 
 2871                      constrainSurface, edges_to_merge, not_merged_edges);
 2872  Toplogy::SetsMap sets_map;
 2873  CHKERR Toplogy(m_field, th)
 
 2874      .classifyVerts(
surface, tets, fixed_edges, corner_nodes, constrainSurface,
 
 2875                     sets_map);
 2877    for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
 2878         sit != sets_map.rend(); sit++) {
 2879      std::string name = "classification_verts_" +
 2880                         boost::lexical_cast<std::string>(sit->first) + ".vtk";
 2881      if (!sit->second.empty())
 2882        CHKERR SaveData(moab)(name, sit->second);
 
 2883    }
 2884  }
 2886  CHKERR Toplogy(m_field, th).getProcTets(tets, edges_to_merge, proc_tets);
 
 2887  out_tets = subtract(tets, proc_tets);
 2888 
 2889  if (bit_ptr) {
 2890    Range all_out_ents = out_tets;
 
 2891    for (
int dd = 2; 
dd >= 0; 
dd--) {
 
 2892      CHKERR moab.get_adjacencies(out_tets, dd, 
false, all_out_ents,
 
 2893                                  moab::Interface::UNION);
 2894    }
 2895    CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(all_out_ents,
 
 2896                                                                 *bit_ptr);
 2897  }
 2898 
 2899  int nb_nodes_merged = 0;
 2900  LengthMapData_multi_index length_map;
 2902 
 2903  auto save_merge_step = [&](
const int pp, 
const Range collapsed_edges) {
 
 2906    CHKERR moab.get_adjacencies(proc_tets, 2, 
false, adj_faces,
 
 2907                                moab::Interface::UNION);
 2908    std::string name;
 2909    name = "node_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
 2911        name, unite(intersect(new_surf, adj_faces), collapsed_edges));
 2912    name =
 2913        "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
 2915        name, unite(intersect(new_surf, adj_faces), edges_to_merge));
 2917  };
 2918 
 2921 
 2922  double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
 2923  for (int pp = 0; pp != nbMaxMergingCycles; ++pp) {
 2924 
 2925    int nb_nodes_merged_p = nb_nodes_merged;
 2926    length_map.clear();
 2927    min_pp = min_p;
 2928    min_p = min;
 2929    CHKERR LengthMap(m_field, th, aveLength)(proc_tets, edges_to_merge,
 
 2930                                             length_map, ave);
 2931 
 2932    if(!length_map.empty())
 2933      min = length_map.get<2>().begin()->qUality;
 2934    if (pp == 0) {
 2935      ave0 = ave;
 2936    }
 2937 
 2938    int nn = 0;
 2939    Range collapsed_edges;
 
 2940    MergeNodes merge_nodes(m_field, true, th, update_meshsets);
 2941 
 2942    for (auto mit = length_map.get<0>().begin();
 2943         mit != length_map.get<0>().end(); mit++, nn++) {
 2944 
 2945      if (!mit->skip) {
 2946 
 2947        auto get_father_and_mother =
 2950              int num_nodes;
 2952              CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes, 
true);
 
 2953              std::array<int, 2> conn_type = {0, 0};
 2954              for (int nn = 0; nn != 2; nn++) {
 2955                for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
 2956                     sit != sets_map.rend(); sit++) {
 2957                  if (sit->second.find(conn[nn]) != sit->second.end()) {
 2958                    conn_type[nn] |= sit->first;
 2959                  }
 2960                }
 2961              }
 2962              int type_father, type_mother;
 2963              if (conn_type[0] > conn_type[1]) {
 2964                father = conn[0];
 2965                mother = conn[1];
 2966                type_father = conn_type[0];
 2967                type_mother = conn_type[1];
 2968              } else {
 2969                father = conn[1];
 2970                mother = conn[0];
 2971                type_father = conn_type[1];
 2972                type_mother = conn_type[0];
 2973              }
 2974              if (type_father == type_mother) {
 2975                line_search = lineSearchSteps;
 2976              }
 2978            };
 2979 
 2980        int line_search = 0;
 2982        CHKERR get_father_and_mother(father, mother, line_search);
 
 2983        CHKERR merge_nodes.mergeNodes(line_search, father, mother, proc_tets);
 
 2984        if (m_field.getInterface<NodeMergerInterface>()->getSuccessMerge()) {
 2985          const EntityHandle father_and_mother[] = {father, mother};
 
 2987          CHKERR moab.get_adjacencies(father_and_mother, 1, 3, 
false, adj_tets);
 
 2988          Range adj_tets_nodes;
 
 2989          CHKERR moab.get_connectivity(adj_tets, adj_tets_nodes, 
true);
 
 2991          CHKERR moab.get_adjacencies(adj_tets_nodes, 1, 
false, adj_edges,
 
 2992                                      moab::Interface::UNION);
 2993          for (auto ait : adj_edges) {
 2994            auto miit = length_map.get<1>().find(ait);
 2995            if (miit != length_map.get<1>().end())
 2996              (const_cast<LengthMapData &>(*miit)).skip = true;
 2997          }
 2998          nb_nodes_merged++;
 2999          collapsed_edges.insert(mit->eDge);
 3000        }
 3001 
 3002        if (nn > static_cast<int>(length_map.size() / fraction_level))
 3003          break;
 3004        if (mit->qUality > ave)
 3005          break;
 3006      }
 3007    }
 3008 
 3009    CHKERR merge_nodes.updateRangeByChilds(new_surf, edges_to_merge,
 
 3010                                           not_merged_edges, true);
 3011 
 3013                "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e", pp,
 3014                nb_nodes_merged, ave, min);
 3015 
 3017      CHKERR save_merge_step(pp + 1, collapsed_edges);
 
 3018 
 3019    if (nb_nodes_merged == nb_nodes_merged_p)
 3020      break;
 3021    if (min > 1e-2 && min == min_pp)
 3022      break;
 3023    if (min > ave0)
 3024      break;
 3025 
 3027    CHKERR moab.get_adjacencies(proc_tets, 1, 
false, adj_edges,
 
 3028                                moab::Interface::UNION);
 3029    edges_to_merge = intersect(edges_to_merge, adj_edges);
 3030    CHKERR Toplogy(m_field, th)
 
 3031        .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
 3032                        constrainSurface, edges_to_merge, not_merged_edges);
 3033  }
 3034 
 3035  auto reconstruct_refined_ents = [&]() {
 3039  };
 3040 
 3041  
 3042  
 3043  
 3044  
 3045  
 3046  CHKERR reconstruct_refined_ents();
 
 3047 
 3048  if (bit_ptr)
 3049    CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
 
 3050                                                                 *bit_ptr);
 3051 
 3052  out_tets.merge(proc_tets);
 3054  CHKERR moab.get_adjacencies(out_tets, 2, 
false, adj_faces,
 
 3055                              moab::Interface::UNION);
 3056  new_surf = intersect(new_surf, adj_faces);
 3057 
 3059}
 3060 
 3062    const int fraction_level, const BitRefLevel cut_bit,
 3063    const BitRefLevel trim_bit, 
const BitRefLevel 
bit, 
const Range &
surface,
 
 3064    const Range &fixed_edges, 
const Range &corner_nodes, 
Tag th,
 
 3065    const bool update_meshsets, 
const bool debug) {
 
 3066  CoreInterface &m_field = cOre;
 3069  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 3070      trim_bit, 
BitRefLevel().set(), MBTET, tets_level);
 
 3071 
 3072  Range edges_to_merge;
 
 3073  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
 
 3074      trim_bit, cut_bit | trim_bit, edges_to_merge);
 3075 
 3076  
 3077  Range all_ents_not_in_database_before;
 
 3078  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
 
 3079      all_ents_not_in_database_before);
 3080 
 3081  edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
 3083    CHKERR SaveData(m_field.get_moab())(
"edges_to_merge.vtk", edges_to_merge);
 
 3084 
 3085  Range out_new_tets, out_new_surf;
 
 3086  CHKERR mergeBadEdges(fraction_level, tets_level, 
surface, fixed_edges,
 
 3087                       corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
 3089 
 3090  
 3091  Range all_ents_not_in_database_after;
 
 3092  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
 
 3093      all_ents_not_in_database_after);
 3094 
 3095  
 3096  all_ents_not_in_database_after =
 3097      subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
 3099  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
 
 3100                                                 false);
 3101  for (
auto m : meshsets)
 
 3102    CHKERR m_field.get_moab().remove_entities(
m,
 
 3103                                              all_ents_not_in_database_after);
 3104 
 3105  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
 3106 
 3107  mergedVolumes.swap(out_new_tets);
 3108  mergedSurfaces.swap(out_new_surf);
 3109 
 3111}
 3112 
 3114  CoreInterface &m_field = cOre;
 3115  moab::Interface &moab = m_field.get_moab();
 3119    CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
 
 3120  else
 3121    CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 3123  std::vector<double> coords(3 * nodes.size());
 3124  CHKERR moab.get_coords(nodes, &coords[0]);
 
 3125  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
 
 3127}
 3128 
 3130                                           const BitRefLevel mask) {
 3131  CoreInterface &m_field = cOre;
 3132  moab::Interface &moab = m_field.get_moab();
 3136    CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
 
 3137  else
 3138    CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
 
 3139        bit, mask, MBVERTEX, nodes);
 
 3140  std::vector<double> coords(3 * nodes.size());
 3141  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
 
 3142  CHKERR moab.set_coords(nodes, &coords[0]);
 
 3144}
 3145 
 3146MoFEMErrorCode CutMeshInterface::saveCutEdges(
const std::string prefix) {
 
 3147  CoreInterface &m_field = cOre;
 3148  moab::Interface &moab = m_field.get_moab();
 3150  CHKERR SaveData(moab)(prefix + 
"out_vol.vtk", vOlume);
 
 3151  CHKERR SaveData(moab)(prefix + 
"out_surface.vtk", sUrface);
 
 3152  CHKERR SaveData(moab)(prefix + 
"out_cut_edges.vtk", cutEdges);
 
 3153  CHKERR SaveData(moab)(prefix + 
"out_cut_new_volumes.vtk", cutNewVolumes);
 
 3154  CHKERR SaveData(moab)(prefix + 
"out_cut_new_surfaces.vtk", cutNewSurfaces);
 
 3155  CHKERR SaveData(moab)(prefix + 
"out_cut_zero_distance_ents.vtk",
 
 3156                        zeroDistanceEnts);
 3158}
 3159 
 3161  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
 3163  CHKERR SaveData(moab)(
"out_trim_surface.vtk", sUrface);
 
 3164  CHKERR SaveData(moab)(
"out_trim_new_volumes.vtk", trimNewVolumes);
 
 3165  CHKERR SaveData(moab)(
"out_trim_new_surfaces.vtk", trimNewSurfaces);
 
 3166  CHKERR SaveData(moab)(
"out_trim_edges.vtk", trimEdges);
 
 3168}
 3169 
 3170}
#define CutMeshFunctionBegin
#define MOFEM_LOG_C(channel, severity, format,...)
static PetscErrorCode ierr
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
#define MOFEM_LOG(channel, severity)
Log.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const double n
refractive index of diffusive medium
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
VectorBoundedArray< double, 6 > VectorDouble6
VectorShallowArrayAdaptor< double > VectorAdaptor
VectorBoundedArray< double, 12 > VectorDouble12
UBlasMatrix< double > MatrixDouble
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
UBlasVector< double > VectorDouble
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
MatrixBoundedArray< double, 9 > MatrixDouble3by3
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
constexpr double t
plate stiffness
static constexpr double delta
FTensor::Index< 'm', 3 > m