insert_constraint_prisms.cpp File Reference

Insert prism elements between paired master/slave surface meshsets. More...

#include <MoFEM.hpp>

Go to the source code of this file.

Functions
int	main (int argc, char *argv[])

Variables
static char	help []

Detailed Description

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

Definition in file insert_constraint_prisms.cpp.

Function Documentation

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 524 of file insert_constraint_prisms.cpp.

                                 {
 
  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;
}

Variable Documentation

help

char help[]

static

Initial value:

=
    "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"

Definition at line 17 of file insert_constraint_prisms.cpp.