insert_constraint_prisms.cpp

Go to the documentation of this file.

/** \file insert_constraint_prisms.cpp
 * \brief Insert prism elements between paired master/slave surface meshsets
 *
 * Meshsets are paired by name, using prefixes:
 *   TIE_MASTER_<suffix>
 *   TIE_SLAVE_<suffix>
 *
 * Current limitations:
 * - only tetrahedral volume meshes are supported
 * - only triangular surface meshsets are supported
 */
 
#include <MoFEM.hpp>
 
using namespace MoFEM;
 
static char help[] =
    "Insert constraint prisms between paired master/slave triangle meshsets.\n"
    "\n"
    "Pairs are detected by meshset names:\n"
    "  TIE_MASTER_<suffix>\n"
    "  TIE_SLAVE_<suffix>\n"
    "\n"
    "Prisms are always created from the SLAVE surface.\n"
    "\n"
    "If the MASTER surface has more triangles than the SLAVE surface, the tool\n"
    "only warns and suggests renaming the TIE blocks so the finer surface is\n"
    "named SLAVE.\n"
    "\n"
    "Virtual refinement:\n"
    "  - -tie_ref_level N recursively refines each effective slave triangle N times\n"
    "    in-memory using 1-to-4 midpoint subdivision before ray casting\n"
    "  - -ray_length value overrides the automatic ray length\n"
    "\n"
    "Current limitations:\n"
    "  - only tet meshes are supported\n"
    "  - only triangle surface meshsets are supported\n\n";
 
namespace {
 
constexpr int SPACE_DIM = 3;
FTENSOR_INDEX(SPACE_DIM, i);
FTENSOR_INDEX(SPACE_DIM, j);
 
struct SurfaceMeshsetData {
  const CubitMeshSets *meshsetPtr = nullptr;
  std::string name;
  Range tris;
};
 
constexpr std::array<std::array<int, 3>, 3> TRI_PERMUTATIONS = {{
    {{0, 1, 2}},
    {{1, 2, 0}},
    {{2, 0, 1}}
}};
 
MoFEMErrorCode set_surface_meshset_triangles(
    moab::Interface &moab, const SurfaceMeshsetData &surface,
    const Range &desired_tris) {
  MoFEMFunctionBegin;
 
  if (!surface.meshsetPtr) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Can not update triangles for null meshset pointer");
  }
 
  const auto meshset = surface.meshsetPtr->meshset;
 
  // Keep TIE meshsets as pure surface meshsets. Prism-update utilities add
  // adjacent prisms/quads to any meshset sharing the corresponding faces or
  // edges, but for TIE detection we need these sets to contain only their
  // designated triangles.
  Range vols;
  CHKERR moab.get_entities_by_dimension(meshset, 3, vols, true);
  if (!vols.empty())
    CHKERR moab.remove_entities(meshset, vols);
 
  Range faces;
  CHKERR moab.get_entities_by_dimension(meshset, 2, faces, true);
  Range current_tris = faces.subset_by_type(MBTRI);
  Range other_faces = subtract(faces, current_tris);
  if (!other_faces.empty())
    CHKERR moab.remove_entities(meshset, other_faces);
 
  Range tris_to_remove = subtract(current_tris, desired_tris);
  if (!tris_to_remove.empty())
    CHKERR moab.remove_entities(meshset, tris_to_remove);
 
  Range tris_to_add = subtract(desired_tris, current_tris);
  if (!tris_to_add.empty())
    CHKERR moab.add_entities(meshset, tris_to_add);
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode get_range_bbox_diag(moab::Interface &moab, const Range &ents,
                                   double &diag) {
  MoFEMFunctionBegin;
 
  Range verts;
  CHKERR moab.get_connectivity(ents, verts, true);
  if (verts.empty())
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Can not calculate bounding box for empty vertex range");
 
  std::vector<double> coords(3 * verts.size());
  CHKERR moab.get_coords(verts, coords.data());
 
  double min_c[SPACE_DIM] = {std::numeric_limits<double>::max(),
                             std::numeric_limits<double>::max(),
                             std::numeric_limits<double>::max()};
  double max_c[SPACE_DIM] = {-std::numeric_limits<double>::max(),
                             -std::numeric_limits<double>::max(),
                             -std::numeric_limits<double>::max()};
 
  for (size_t vv = 0; vv != verts.size(); ++vv) {
    for (int dd = 0; dd != 3; ++dd) {
      min_c[dd] = std::min(min_c[dd], coords[3 * vv + dd]);
      max_c[dd] = std::max(max_c[dd], coords[3 * vv + dd]);
    }
  }
 
  FTensor::Tensor1<double, 3> t_diag;
  auto t_max_c = getFTensor1FromPtr<SPACE_DIM>(max_c);
  auto t_min_c = getFTensor1FromPtr<SPACE_DIM>(min_c);
  t_diag(i) = t_max_c(i) - t_min_c(i);
  diag = t_diag.l2();
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode validate_surface_meshset(moab::Interface &moab,
                                        const CubitMeshSets &meshset_data,
                                        const Range &tet_skin,
                                        SurfaceMeshsetData &surface_data) {
  MoFEMFunctionBegin;
 
  surface_data.meshsetPtr = &meshset_data;
  surface_data.name = meshset_data.getName();
 
  Range vols;
  CHKERR moab.get_entities_by_dimension(meshset_data.meshset, 3, vols, true);
  if (!vols.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Meshset %s contains 3D entities. Only surface meshsets are "
            "supported",
            surface_data.name.c_str());
  }
 
  Range faces;
  CHKERR moab.get_entities_by_dimension(meshset_data.meshset, 2, faces, true);
  if (faces.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Meshset %s does not contain surface entities",
            surface_data.name.c_str());
  }
 
  surface_data.tris = faces.subset_by_type(MBTRI);
  const Range other_faces = subtract(faces, surface_data.tris);
  if (!other_faces.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Meshset %s contains non-triangle surface entities. Only "
            "triangles are supported",
            surface_data.name.c_str());
  }
 
  const Range skin_tris = intersect(surface_data.tris, tet_skin);
  if (skin_tris.size() != surface_data.tris.size()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Meshset %s contains triangles that are not on the tet skin",
            surface_data.name.c_str());
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode get_coords_centroid(const double *coords, const int num_nodes,
                                   std::array<double, 3> &centroid) {
  MoFEMFunctionBegin;
  centroid = {0., 0., 0.};
  for (int nn = 0; nn != num_nodes; ++nn)
    for (int dd = 0; dd != SPACE_DIM; ++dd)
      centroid[dd] += coords[SPACE_DIM * nn + dd] / num_nodes;
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode get_skin_triangle_data(moab::Interface &moab,
                                      const EntityHandle tri,
                                      std::array<EntityHandle, 3> *tri_nodes,
                                      std::array<double, SPACE_DIM> &outward_normal) {
  MoFEMFunctionBegin;
 
  const EntityHandle *conn = nullptr;
  int num_nodes = 0;
  CHKERR moab.get_connectivity(tri, conn, num_nodes, true);
  if (num_nodes != 3)
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Only 3-node triangles are supported");
 
  if (tri_nodes) {
    for (int ii = 0; ii != 3; ++ii)
      (*tri_nodes)[ii] = conn[ii];
  }
 
  double tri_coords[9];
  CHKERR moab.get_coords(conn, 3, tri_coords);
 
  std::array<double, SPACE_DIM> centroid;
  CHKERR get_coords_centroid(tri_coords, 3, centroid);
  auto t_centroid = getFTensor1FromPtr<SPACE_DIM>(centroid.data());
 
  auto t_outward_normal = getFTensor1FromPtr<SPACE_DIM>(outward_normal.data());
  CHKERR Tools::getTriNormal(tri_coords, &t_outward_normal(0));
  const double tri_norm = t_outward_normal.l2();
  if (tri_norm <= std::numeric_limits<double>::epsilon()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Triangle %lu has zero normal", tri);
  }
  t_outward_normal(i) /= tri_norm;
 
  Range adj_tets;
  CHKERR moab.get_adjacencies(&tri, 1, 3, false, adj_tets);
  adj_tets = adj_tets.subset_by_type(MBTET);
  if (adj_tets.size() != 1) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Triangle %lu is expected to have exactly one adjacent tet", tri);
  }
 
  const EntityHandle tet = adj_tets[0];
  const EntityHandle *tet_conn = nullptr;
  int tet_num_nodes = 0;
  CHKERR moab.get_connectivity(tet, tet_conn, tet_num_nodes, true);
  if (tet_num_nodes != 4)
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Only 4-node tets are supported");
 
  double tet_coords[12];
  CHKERR moab.get_coords(tet_conn, 4, tet_coords);
  std::array<double, SPACE_DIM> tet_centroid;
  CHKERR get_coords_centroid(tet_coords, 4, tet_centroid);
  auto t_tet_centroid = getFTensor1FromPtr<SPACE_DIM>(tet_centroid.data());
 
  FTensor::Tensor1<double, 3> t_delta;
  t_delta(i) = t_tet_centroid(i) - t_centroid(i);
 
  if (t_outward_normal(i) * t_delta(i) > 0)
    t_outward_normal(i) *= -1.;
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode get_triangle_coords(moab::Interface &moab, const EntityHandle tri,
                                   std::array<double, 9> &coords) {
  MoFEMFunctionBegin;
  const EntityHandle *conn = nullptr;
  int num_nodes = 0;
  CHKERR moab.get_connectivity(tri, conn, num_nodes, true);
  if (num_nodes != 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Only 3-node triangles are supported");
  }
  CHKERR moab.get_coords(conn, 3, coords.data());
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode collect_virtual_refined_triangle_centroids(
    const int ref_level, MatrixDouble &reference_child_centroids) {
  MoFEMFunctionBegin;
 
  moab::Core moab_ref;
  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
  EntityHandle nodes[3];
  for (int nn = 0; nn != 3; ++nn)
    CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
  EntityHandle tri;
  CHKERR moab_ref.create_element(MBTRI, nodes, 3, tri);
 
  MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
  MoFEM::Interface &m_field_ref = mofem_ref_core;
  auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
  auto *bit_ref_manager = m_field_ref.getInterface<BitRefManager>();
 
  {
    Range tris(tri, tri);
    Range edges;
    CHKERR m_field_ref.get_moab().get_adjacencies(tris, 1, true, edges,
                                                  moab::Interface::UNION);
    CHKERR bit_ref_manager->setBitRefLevel(tris, BitRefLevel().set(0), false,
                                           VERBOSE);
  }
#ifndef NDEBUG
  const bool debug = true;
#else
  const bool debug = false;
#endif
 
  for (int ll = 0; ll != ref_level; ++ll) {
    Range current_tris;
    CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
        BitRefLevel().set(ll), BitRefLevel().set(), MBTRI, current_tris);
    const BitRefLevel bit = BitRefLevel().set(ll + 1);
    Range edges;
    CHKERR moab_ref.get_adjacencies(current_tris, 1, true, edges,
                                    moab::Interface::UNION);
    CHKERR m_ref->addVerticesInTheMiddleOfEdges(edges, bit);
    CHKERR m_ref->refineTris(current_tris, bit, QUIET, debug);
  }
 
  Range refined_tris;
  CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
      BitRefLevel().set(ref_level), BitRefLevel().set(), MBTRI, refined_tris);
  reference_child_centroids.resize(refined_tris.size(), SPACE_DIM, false);
  int rr = 0;
  for (auto child_tri : refined_tris) {
    std::array<double, 9> child_coords;
    std::array<double, 3> child_centroid;
    CHKERR get_triangle_coords(moab_ref, child_tri, child_coords);
    CHKERR get_coords_centroid(child_coords.data(), 3, child_centroid);
    for (int dd = 0; dd != SPACE_DIM; ++dd)
      reference_child_centroids(rr, dd) = child_centroid[dd];
    ++rr;
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode map_reference_centroids_to_triangle(
    const MatrixDouble &reference_child_centroids,
    const std::array<double, 9> &tri_coords, MatrixDouble &child_centroids) {
  MoFEMFunctionBegin;
 
  child_centroids.resize(reference_child_centroids.size1(), SPACE_DIM, false);
  std::array<double, 9> mutable_coords = tri_coords;
  auto t_p0 = getFTensor1FromPtr<SPACE_DIM>(&mutable_coords[0]);
  auto t_p1 = getFTensor1FromPtr<SPACE_DIM>(&mutable_coords[3]);
  auto t_p2 = getFTensor1FromPtr<SPACE_DIM>(&mutable_coords[6]);
 
  for (size_t rr = 0; rr != reference_child_centroids.size1(); ++rr) {
    const double n1 = reference_child_centroids(rr, 0);
    const double n2 = reference_child_centroids(rr, 1);
    const double n0 = 1. - n1 - n2;
    std::array<double, 3> child_centroid = {0., 0., 0.};
    auto t_child_centroid =
        getFTensor1FromPtr<SPACE_DIM>(child_centroid.data());
    t_child_centroid(i) = n0 * t_p0(i) + n1 * t_p1(i) + n2 * t_p2(i);
    for (int dd = 0; dd != SPACE_DIM; ++dd)
      child_centroids(rr, dd) = child_centroid[dd];
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode collect_target_triangles(
    OrientedBoxTreeTool &tree, const EntityHandle root_set,
    const std::array<double, 3> &source_centroid,
    const std::array<double, 3> &source_normal, const double ray_length,
    const double ray_tol,
    std::set<EntityHandle> &target_tris) {
  MoFEMFunctionBegin;
 
  auto collect_closest_hit = [&](const double ray_sign, EntityHandle &closest_tri,
                                 double &closest_dist) {
    MoFEMFunctionBeginHot;
 
    std::vector<double> distances;
    std::vector<EntityHandle> facets;
    double ray_point[SPACE_DIM] = {source_centroid[0], source_centroid[1],
                                   source_centroid[2]};
    double ray_dir[SPACE_DIM] = {ray_sign * source_normal[0],
                                 ray_sign * source_normal[1],
                                 ray_sign * source_normal[2]};
    double max_ray_length = ray_length;
 
    CHKERR tree.ray_intersect_triangles(distances, facets, root_set, ray_tol,
                                        ray_point, ray_dir, &max_ray_length);
 
    closest_tri = 0;
    closest_dist = std::numeric_limits<double>::max();
    for (size_t ii = 0; ii != facets.size(); ++ii) {
      const double dist = distances[ii];
      if (dist < -ray_tol)
        continue;
      if (dist < closest_dist) {
        closest_dist = dist;
        closest_tri = facets[ii];
      }
    }
 
    MoFEMFunctionReturnHot(0);
  };
 
  EntityHandle forward_tri = 0;
  EntityHandle reverse_tri = 0;
  double forward_dist = std::numeric_limits<double>::max();
  double reverse_dist = std::numeric_limits<double>::max();
  CHKERR collect_closest_hit(1., forward_tri, forward_dist);
  CHKERR collect_closest_hit(-1., reverse_tri, reverse_dist);
 
  if (forward_dist <= reverse_dist) {
    if (forward_tri)
      target_tris.insert(forward_tri);
  } else if (reverse_tri) {
    target_tris.insert(reverse_tri);
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode match_triangle_vertices(moab::Interface &moab,
                                       const std::array<EntityHandle, 3> &master_nodes,
                                       const std::array<EntityHandle, 3> &slave_nodes,
                                       std::array<int, 3> &perm) {
  MoFEMFunctionBegin;
 
  double master_coords[9];
  double slave_coords[9];
  CHKERR moab.get_coords(master_nodes.data(), 3, master_coords);
  CHKERR moab.get_coords(slave_nodes.data(), 3, slave_coords);
 
  double best_score = std::numeric_limits<double>::max();
  bool found = false;
 
  for (const auto &candidate_perm : TRI_PERMUTATIONS) {
    double score = 0;
 
    for (int ii = 0; ii != 3; ++ii) {
      auto t_master_c = getFTensor1FromPtr<3>(&master_coords[3 * ii]);
      auto t_slave_c =
          getFTensor1FromPtr<3>(&slave_coords[3 * candidate_perm[ii]]);
      FTensor::Tensor1<double, SPACE_DIM> t_delta;
      t_delta(j) = t_master_c(j) - t_slave_c(j);
      score += t_delta(j) * t_delta(j);
    }
 
    if (score < best_score) {
      best_score = score;
      perm = candidate_perm;
      found = true;
    }
  }
 
  if (!found) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Could not find a vertex mapping between paired master and slave "
            "triangles");
  }
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode create_prism_from_triangle_pair(
    PrismsFromSurfaceInterface &prisms_from_surface, const EntityHandle source_tri,
    const EntityHandle target_tri, EntityHandle &prism) {
  MoFEMFunctionBegin;
  MoFEM::Interface &m_field = prisms_from_surface.cOre;
  auto &moab = m_field.get_moab();
 
  std::array<EntityHandle, 3> source_nodes;
  std::array<double, SPACE_DIM> source_normal;
  CHKERR get_skin_triangle_data(moab, source_tri, &source_nodes, source_normal);
  auto t_source_normal = getFTensor1FromPtr<SPACE_DIM>(source_normal.data());
 
  const EntityHandle *target_conn = nullptr;
  int target_num_nodes = 0;
  CHKERR moab.get_connectivity(target_tri, target_conn, target_num_nodes, true);
  if (target_num_nodes != 3)
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Only 3-node target triangles are supported");
 
  const EntityHandle *source_conn = nullptr;
  int source_num_nodes = 0;
  CHKERR moab.get_connectivity(source_tri, source_conn, source_num_nodes, true);
  if (source_num_nodes != 3)
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
            "Only 3-node source triangles are supported");
 
  double source_coords[9];
  CHKERR moab.get_coords(source_conn, 3, source_coords);
  std::array<double, SPACE_DIM> raw_source_normal;
  auto t_raw_source_normal =
      getFTensor1FromPtr<SPACE_DIM>(raw_source_normal.data());
  CHKERR Tools::getTriNormal(source_coords, &t_raw_source_normal(0));
 
  std::array<EntityHandle, 3> prism_source_nodes = source_nodes;
  if (t_raw_source_normal(i) * t_source_normal(i) < 0)
    std::swap(prism_source_nodes[1], prism_source_nodes[2]);
 
  std::array<EntityHandle, 3> target_nodes = {target_conn[0], target_conn[1],
                                              target_conn[2]};
  std::array<int, 3> target_perm = {0, 1, 2};
  CHKERR match_triangle_vertices(moab, prism_source_nodes, target_nodes,
                                 target_perm);
 
  std::array<EntityHandle, 3> prism_target_nodes = {
      target_nodes[target_perm[0]], target_nodes[target_perm[1]],
      target_nodes[target_perm[2]]};
 
  CHKERR prisms_from_surface.createPrism(
      prism_source_nodes.data(), prism_target_nodes.data(),
      PrismsFromSurfaceInterface::NO_SWAP, prism);
 
  int side_number = -1;
  int sense = 0;
  int offset = 0;
 
  CHKERR moab.side_number(prism, source_tri, side_number, sense, offset);
  if (side_number != 3) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Created prism %lu does not use source triangle %lu as side 3",
            prism, source_tri);
  }
 
  CHKERR moab.side_number(prism, target_tri, side_number, sense, offset);
  if (side_number != 4) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Created prism %lu does not use target triangle %lu as side 4",
            prism, target_tri);
  }
 
  MoFEMFunctionReturn(0);
}
 
} // namespace
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    char mesh_file_name[255] = "mesh.h5m";
    char mesh_out_file[255] = "out.h5m";
    char master_prefix[255] = "TIE_MASTER_";
    char slave_prefix[255] = "TIE_SLAVE_";
    PetscBool flg_my_file = PETSC_FALSE;
    PetscBool flg_ray_length = PETSC_FALSE;
    PetscInt ref_level = 1;
    PetscReal ray_length_option = 0;
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "", "Insert constraint prisms", "none");
    CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
                              mesh_file_name, 255, &flg_my_file);
    if (flg_my_file != PETSC_TRUE)
      CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
                                mesh_file_name, 255, &flg_my_file);
    CHKERR PetscOptionsString("-output_file", "output mesh file name", "",
                              "out.h5m", mesh_out_file, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsString("-master_prefix", "master meshset name prefix",
                              "", "TIE_MASTER_", master_prefix, 255,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsString("-slave_prefix", "slave meshset name prefix", "",
                              "TIE_SLAVE_", slave_prefix, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsInt("-tie_ref_level",
                           "virtual recursive refinement level on slave "
                           "triangles",
                           "", ref_level, &ref_level, PETSC_NULLPTR);
    CHKERR PetscOptionsReal("-ray_length", "absolute ray length", "",
                            ray_length_option, &ray_length_option,
                            &flg_ray_length);
    PetscOptionsEnd();
 
    if (flg_my_file != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "error -my_file (-file_name) (mesh file needed)");
    }
    if (ref_level < 1) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "-tie_ref_level must be greater than or equal to 1");
    }
    if (flg_ray_length == PETSC_TRUE && ray_length_option <= 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
              "-ray_length must be greater than 0");
    }
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    if (pcomm == NULL)
      pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
 
    CHKERR moab.load_file(mesh_file_name);
 
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    Range tets;
    CHKERR moab.get_entities_by_type(0, MBTET, tets, false);
    if (tets.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "No tetrahedra found. Only tet meshes are currently supported");
    }
 
    Range other_vols;
    CHKERR moab.get_entities_by_dimension(0, 3, other_vols, false);
    other_vols = subtract(other_vols, tets);
    if (!other_vols.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
              "Only tetrahedral volume meshes are currently supported");
    }
 
    Range tet_skin;
    {
      Skinner skinner(&moab);
      CHKERR skinner.find_skin(0, tets, false, tet_skin);
      tet_skin = tet_skin.subset_by_type(MBTRI);
    }
 
    const std::string master_prefix_str(master_prefix);
    const std::string slave_prefix_str(slave_prefix);
 
    std::map<std::string, const CubitMeshSets *> master_sets;
    std::map<std::string, const CubitMeshSets *> slave_sets;
 
    for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, it)) {
      const std::string name = it->getName();
 
      if (name.compare(0, master_prefix_str.size(), master_prefix_str) == 0) {
        const auto suffix = name.substr(master_prefix_str.size());
        if (!master_sets.emplace(suffix, &*it).second) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "More than one master meshset uses suffix %s", suffix.c_str());
        }
      }
 
      if (name.compare(0, slave_prefix_str.size(), slave_prefix_str) == 0) {
        const auto suffix = name.substr(slave_prefix_str.size());
        if (!slave_sets.emplace(suffix, &*it).second) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "More than one slave meshset uses suffix %s", suffix.c_str());
        }
      }
    }
 
    if (master_sets.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
              "No master meshsets found with prefix %s",
              master_prefix_str.c_str());
    }
 
    for (const auto &slave_pair : slave_sets) {
      const auto &suffix = slave_pair.first;
      if (master_sets.find(suffix) == master_sets.end()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Found slave meshset with suffix %s but no matching master "
                "meshset",
                suffix.c_str());
      }
    }
 
    PrismsFromSurfaceInterface *prisms_from_surface = nullptr;
    CHKERR m_field.getInterface(prisms_from_surface);
 
    MatrixDouble reference_child_centroids;
    CHKERR collect_virtual_refined_triangle_centroids(ref_level,
                                                      reference_child_centroids);
 
    Range all_prisms;
 
    for (const auto &[suffix, master_ptr] : master_sets) {
      auto slave_it = slave_sets.find(suffix);
      if (slave_it == slave_sets.end()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Found master meshset with suffix %s but no matching slave "
                "meshset",
                suffix.c_str());
      }
 
      const auto *slave_ptr = slave_it->second;
 
      if (master_ptr->getBcType() != slave_ptr->getBcType()) {
        SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                "Master/slave meshset types differ for suffix %s",
                suffix.c_str());
      }
 
      SurfaceMeshsetData master_surface;
      CHKERR validate_surface_meshset(moab, *master_ptr, tet_skin,
                                      master_surface);
 
      SurfaceMeshsetData slave_surface;
      CHKERR validate_surface_meshset(moab, *slave_ptr, tet_skin, slave_surface);
 
      const auto master_tri_count = master_surface.tris.size();
      const auto slave_tri_count = slave_surface.tris.size();
      double bbox_diag = 0;
      CHKERR get_range_bbox_diag(moab,
                                 unite(master_surface.tris, slave_surface.tris),
                                 bbox_diag);
 
      const double geom_eps =
          std::max(1e-12 * std::max(1., bbox_diag), 10. *
                                                        std::numeric_limits<double>::epsilon());
      const double ray_tol = geom_eps;
      const double ray_length =
          (flg_ray_length == PETSC_TRUE) ? ray_length_option
                                         : std::max(10. * bbox_diag, 1.);
 
      OrientedBoxTreeTool master_tree(&moab, "ROOTSET_CONSTRAINT_PRISMS", true);
      EntityHandle master_root = 0;
      CHKERR master_tree.build(master_surface.tris, master_root);
 
      Range pair_prisms;
      size_t unmatched_slave_tris = 0;
 
      for (auto slave_tri : slave_surface.tris) {
        std::array<double, SPACE_DIM> slave_normal;
        CHKERR get_skin_triangle_data(moab, slave_tri, nullptr, slave_normal);
 
        std::array<double, 9> slave_coords;
        CHKERR get_triangle_coords(moab, slave_tri, slave_coords);
        MatrixDouble child_centroids;
        CHKERR map_reference_centroids_to_triangle(reference_child_centroids,
                                                   slave_coords,
                                                   child_centroids);
 
        std::set<EntityHandle> hit_master_tris;
        for (size_t rr = 0; rr != child_centroids.size1(); ++rr) {
          std::array<double, 3> child_centroid = {
              child_centroids(rr, 0), child_centroids(rr, 1),
              child_centroids(rr, 2)};
          CHKERR collect_target_triangles(master_tree, master_root,
                                          child_centroid, slave_normal,
                                          ray_length, ray_tol,
                                          hit_master_tris);
        }
 
        if (hit_master_tris.empty()) {
          ++unmatched_slave_tris;
          continue;
        }
 
        for (auto master_tri : hit_master_tris) {
          EntityHandle prism = 0;
          CHKERR create_prism_from_triangle_pair(*prisms_from_surface, slave_tri,
                                                 master_tri, prism);
          pair_prisms.insert(prism);
        }
      }
 
      CHKERR prisms_from_surface->updateMeshestByEdgeBlock(pair_prisms);
      CHKERR prisms_from_surface->updateMeshestByTriBlock(pair_prisms);
 
      CHKERR set_surface_meshset_triangles(moab, master_surface,
                                           master_surface.tris);
      CHKERR set_surface_meshset_triangles(moab, slave_surface,
                                           slave_surface.tris);
 
      all_prisms.merge(pair_prisms);
 
      MOFEM_LOG("WORLD", Sev::inform)
          << "Inserted " << pair_prisms.size()
          << " prisms for master/slave suffix " << suffix;
      if (master_tri_count > slave_tri_count) {
        MOFEM_LOG("WORLD", Sev::warning)
            << "Master surface has more triangles than slave surface for "
               "suffix "
            << suffix << " (" << master_tri_count << " > "
            << slave_tri_count
            << "); consider renaming the TIE blocks so the finer surface is "
               "named TIE_SLAVE_"
            << suffix;
      }
      MOFEM_LOG("WORLD", Sev::warning)
          << "Skipped " << unmatched_slave_tris
          << " slave triangles without master hits for suffix "
          << suffix;
    }
 
    MOFEM_LOG("WORLD", Sev::inform)
        << "Inserted total " << all_prisms.size() << " constraint prisms";
 
    CHKERR moab.write_file(mesh_out_file);
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}