Split sidesets. More...
#include <MoFEM.hpp>
Functions
int	main (int argc, char *argv[])
Variables
static char	help [] = "...\n\n"

constexpr bool	debug = false
Detailed Description

Split sidesets.
Definition in file split_sideset.cpp.
Function Documentation

◆ main()

int main	(	int	argc,
		char *	argv[]
	)
Examples: split_sideset.cpp.
Definition at line 15 of file split_sideset.cpp.
                                  {
   MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
  
   auto core_log = logging::core::get();
   core_log->add_sink(
       LogManager::createSink(LogManager::getStrmWorld(), "SPLIT"));
   LogManager::setLog("SPLIT");
   MOFEM_LOG_TAG("SPLIT", "split");
  
   try {
  
     // global variables
     char mesh_file_name[255];
     PetscBool flg_file = PETSC_FALSE;
     PetscBool squash_bit_levels = PETSC_TRUE;
     PetscBool flg_list_of_sidesets = PETSC_FALSE;
     PetscBool output_vtk = PETSC_TRUE;
     PetscBool add_prisms = PETSC_FALSE;
     PetscBool split_corner_edges = PETSC_FALSE;
     int nb_sidesets = 10;
     int sidesets[nb_sidesets];
  
     CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Split sides options",
                              "none");
     CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
                               mesh_file_name, 255, &flg_file);
     CHKERR PetscOptionsBool("-squash_bit_levels", "squash bit levels", "",
                             squash_bit_levels, &squash_bit_levels, NULL);
     CHKERR PetscOptionsIntArray("-side_sets", "get list of sidesets", "",
                                 sidesets, &nb_sidesets, &flg_list_of_sidesets);
     CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
                             output_vtk, &output_vtk, PETSC_NULL);
     CHKERR PetscOptionsBool("-add_prisms", "if true outout vtk file", "",
                             add_prisms, &add_prisms, PETSC_NULL);
     CHKERR PetscOptionsBool("-split_corner_edges", "if true outout vtk file",
                             "", split_corner_edges, &split_corner_edges,
                             PETSC_NULL);
     CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
                             output_vtk, &output_vtk, PETSC_NULL);
     ierr = PetscOptionsEnd();
     CHKERRG(ierr);
  
     moab::Core mb_instance;
     moab::Interface &moab = mb_instance;
  
     const char *option;
     option = "";
     CHKERR moab.load_file(mesh_file_name, 0, option);
  
     // Create MoFEM  database
     MoFEM::Core core(moab);
     MoFEM::Interface &m_field = core;
  
     if (flg_file != PETSC_TRUE) {
       SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
               "*** ERROR -my_file (mesh file needed)");
     }
     if (flg_list_of_sidesets != PETSC_TRUE) {
       SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
               "List of sidesets not given -my_side_sets ...");
     }
  
     // Get interface to meshsets manager
     auto m_mng = m_field.getInterface<MeshsetsManager>();
     // Get interface for splitting manager
     auto interface_ptr = m_field.getInterface<PrismInterface>();
     // Managing bits
     auto bit_mng = m_field.getInterface<BitRefManager>();
     // Refine mesh
     auto refine_mng = m_field.getInterface<MeshRefinement>();
  
     auto &meshsets_index = m_mng->getMeshsetsMultindex();
     auto &m_by_type = meshsets_index.get<CubitMeshsetMaskedType_mi_tag>();
     auto &m_by_id_and_type =
         meshsets_index.get<Composite_Cubit_msId_And_MeshsetType_mi_tag>();
  
     for (auto mit = m_by_type.lower_bound(SIDESET);
          mit != m_by_type.upper_bound(SIDESET); mit++) {
       MOFEM_LOG("SPLIT", Sev::verbose)
           << "Sideset on the mesh id = " << mit->getMeshsetId();
     }
  
     CHKERR bit_mng->setBitRefLevelByDim(0, 3, BitRefLevel().set(0));
     std::vector<BitRefLevel> bit_levels;
     bit_levels.push_back(BitRefLevel().set(0));
  
     auto update_meshsets = [&]() {
       MoFEMFunctionBegin;
       for (_IT_CUBITMESHSETS_FOR_LOOP_(m_field, ciit)) {
         EntityHandle cubit_meshset = ciit->meshset;
         CHKERR bit_mng->updateMeshsetByEntitiesChildren(
             cubit_meshset, bit_levels.back(), cubit_meshset, MBMAXTYPE, true);
       }
       MoFEMFunctionReturn(0);
     };
  
     // refine corner mesh
     if (split_corner_edges) {
       Skinner skin(&m_field.get_moab());
       auto meshset_of_edges_to_refine_ptr = get_temp_meshset_ptr(moab);
  
       for (int mm = 0; mm != nb_sidesets; mm++) {
  
         // find side set
         auto mit =
             m_by_id_and_type.find(boost::make_tuple(sidesets[mm], SIDESET));
         if (mit == m_by_id_and_type.end())
           SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                    "No sideset in database id = %d", sidesets[mm]);
  
         Range faces;
         CHKERR moab.get_entities_by_type(mit->getMeshset(), MBTRI, faces, true);
         Range faces_edges;
         CHKERR moab.get_adjacencies(faces, 1, true, faces_edges,
                                     moab::Interface::UNION);
  
         Range skin_edges;
         CHKERR skin.find_skin(0, faces, false, skin_edges);
         Range skin_verts;
         CHKERR moab.get_connectivity(skin_edges, skin_verts, true);
         Range vertex_edges;
         CHKERR moab.get_adjacencies(skin_verts, 1, true, vertex_edges,
                                     moab::Interface::UNION);
         Range vertex_edges_verts;
         CHKERR moab.get_connectivity(vertex_edges, vertex_edges_verts, true);
         vertex_edges_verts = subtract(vertex_edges_verts, skin_verts);
         Range vertex_edges_verts_edges;
         CHKERR moab.get_adjacencies(vertex_edges_verts, 1, true,
                                     vertex_edges_verts_edges,
                                     moab::Interface::UNION);
         vertex_edges = subtract(vertex_edges, vertex_edges_verts_edges);
         vertex_edges = subtract(vertex_edges, skin_edges);
         vertex_edges = intersect(vertex_edges, faces_edges);
         CHKERR moab.add_entities(*meshset_of_edges_to_refine_ptr, vertex_edges);
       }
  
       int nb_tris;
       CHKERR moab.get_number_entities_by_type(*meshset_of_edges_to_refine_ptr,
                                               MBEDGE, nb_tris, true);
       MOFEM_LOG("SPLIT", Sev::inform) << "Refine corner edges " << nb_tris;
  
       if (debug)
         CHKERR moab.write_file("out_edges_to_refine.vtk", "VTK", "",
                                meshset_of_edges_to_refine_ptr->get_ptr(), 1);
  
       bit_levels.push_back(BitRefLevel().set(1));
       CHKERR refine_mng->addVerticesInTheMiddleOfEdges(
           *meshset_of_edges_to_refine_ptr, bit_levels.back());
       CHKERR refine_mng->refineTets(0, bit_levels.back());
       CHKERR update_meshsets();
     }
  
     // iterate sideset and split
     for (int mm = 0; mm != nb_sidesets; mm++) {
  
       // find side set
       auto mit =
           m_by_id_and_type.find(boost::make_tuple(sidesets[mm], SIDESET));
       if (mit == m_by_id_and_type.end())
         SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
                  "No sideset in database id = %d", sidesets[mm]);
  
       MOFEM_LOG("SPLIT", Sev::inform)
           << "Split sideset " << mit->getMeshsetId();
  
       const auto cubit_meshset = mit->getMeshset();
       {
         // get tet entities form back bit_level
         auto ref_level_meshset_ptr = get_temp_meshset_ptr(moab);
         CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit_levels.back(),
                                                      BitRefLevel().set(), MBTET,
                                                      *ref_level_meshset_ptr);
         CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
             bit_levels.back(), BitRefLevel().set(), MBPRISM,
             *ref_level_meshset_ptr);
         Range ref_level_tets;
         CHKERR moab.get_entities_by_handle(*ref_level_meshset_ptr,
                                            ref_level_tets, true);
  
         // get faces and test to split
         CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), true,
                                        0);
         // set new bit level
         MOFEM_LOG("SPLIT", Sev::verbose)
             << "Add bit level " << bit_levels.size();
         bit_levels.push_back(BitRefLevel().set(bit_levels.size()));
         // split faces and
         CHKERR interface_ptr->splitSides(*ref_level_meshset_ptr,
                                          bit_levels.back(), cubit_meshset,
                                          add_prisms, true, 0);
       }
       // Update cubit meshsets
       CHKERR update_meshsets();
     }
  
     if (squash_bit_levels == PETSC_TRUE) {
       for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
         CHKERR m_field.delete_ents_by_bit_ref(bit_levels[ll], bit_levels[ll],
                                               true);
       }
       CHKERR bit_mng->shiftRightBitRef(bit_levels.size() - 1);
     }
  
     if (output_vtk) {
       auto meshset_ptr = get_temp_meshset_ptr(moab);
       BitRefLevel bit;
       if (squash_bit_levels)
         bit = bit_levels[0];
       else
         bit = bit_levels.back();
       CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit, BitRefLevel().set(),
                                                    MBTET, *meshset_ptr);
       CHKERR moab.write_file("out.vtk", "VTK", "", meshset_ptr->get_ptr(), 1);
     }
  
     ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
     if (pcomm == NULL)
       SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
               "Communicator should be allocated");
  
     CHKERR pcomm->assign_global_ids(0, 3, 1, false);
     CHKERR moab.write_file("out.h5m");
   }
   CATCH_ERRORS;
  
   MoFEM::Core::Finalize();
  
   return 0;
 }