13                                 {
   15 
   16  auto core_log = logging::core::get();
   17  core_log->add_sink(
   21 
   22  try {
   23 
   24    
   26    PetscBool flg_file = PETSC_FALSE;
   27    PetscBool squash_bit_levels = PETSC_TRUE;
   29    PetscBool add_prisms = PETSC_FALSE;
   30    PetscBool split_corner_edges = PETSC_FALSE;
   31    char block_name[255] = "CRACK";
   32    PetscBool flg_block = PETSC_FALSE;
   33    int nb_sidesets = 10;
   34    int sidesets[nb_sidesets];
   35    PetscBool flg_list_of_sidesets = PETSC_FALSE;
   36 
   37    PetscOptionsBegin(PETSC_COMM_WORLD, "", "Split sides options",
   38                             "none");
   39    CHKERR PetscOptionsString(
"-my_file", 
"mesh file name", 
"", 
"mesh.h5m",
 
   41    CHKERR PetscOptionsString(
"-file_name", 
"mesh file name", 
"", 
"mesh.h5m",
 
   43    CHKERR PetscOptionsBool(
"-squash_bit_levels", 
"squash bit levels", 
"",
 
   44                            squash_bit_levels, &squash_bit_levels, NULL);
   45    CHKERR PetscOptionsIntArray(
"-side_sets", 
"get list of sidesets", 
"",
 
   46                                sidesets, &nb_sidesets, &flg_list_of_sidesets);
   47    CHKERR PetscOptionsString(
"-block_name", 
"block_name", 
"", 
"mesh.h5m",
 
   48                              block_name, 255, &flg_block);
   49    CHKERR PetscOptionsBool(
"-add_prisms", 
"add prisms", 
"", add_prisms,
 
   50                            &add_prisms, PETSC_NULLPTR);
   51    CHKERR PetscOptionsBool(
"-split_corner_edges", 
"if true output vtk file",
 
   52                            "", split_corner_edges, &split_corner_edges,
   53                            PETSC_NULLPTR);
   54    CHKERR PetscOptionsBool(
"-output_vtk", 
"if true output vtk file", 
"",
 
   56    PetscOptionsEnd();
   57 
   58    moab::Core mb_instance;
   59    moab::Interface &moab = mb_instance;
   60 
   61    const char *option;
   62    option = "";
   64 
   65    
   68 
   69    if (flg_file != PETSC_TRUE) {
   71              "*** ERROR -my_file (mesh file needed)");
   72    }
   73    if (flg_list_of_sidesets != PETSC_TRUE && flg_block != PETSC_TRUE) {
   75              "Block name or kist of sidesets not given ...");
   76    }
   77 
   78    
   80    
   82    
   84    
   86 
   87    for (
auto m : m_mng->getCubitMeshsetPtr(
SIDESET)) {
 
   89          << 
"Sideset on the mesh id = " << 
m->getMeshsetId();
 
   90    }
   91 
   93    std::vector<BitRefLevel> bit_levels;
   95 
   96    auto update_meshsets = [&]() {
  100        CHKERR bit_mng->updateMeshsetByEntitiesChildren(
 
  101            cubit_meshset, bit_levels.back(), cubit_meshset, MBMAXTYPE, true);
  102      }
  104    };
  105 
  106    
  107    if (split_corner_edges) {
  110 
  111      auto refine = [&](auto mit, auto meshset_of_edges_to_refine_ptr) {
  114        CHKERR moab.get_entities_by_type(mit->getMeshset(), MBTRI, faces, 
true);
 
  116        CHKERR moab.get_adjacencies(faces, 1, 
true, faces_edges,
 
  117                                    moab::Interface::UNION);
  118 
  120        CHKERR skin.find_skin(0, faces, 
false, skin_edges);
 
  122        CHKERR moab.get_connectivity(skin_edges, skin_verts, 
true);
 
  124        CHKERR moab.get_adjacencies(skin_verts, 1, 
true, vertex_edges,
 
  125                                    moab::Interface::UNION);
  126        Range vertex_edges_verts;
 
  127        CHKERR moab.get_connectivity(vertex_edges, vertex_edges_verts, 
true);
 
  128        vertex_edges_verts = subtract(vertex_edges_verts, skin_verts);
  129        Range vertex_edges_verts_edges;
 
  130        CHKERR moab.get_adjacencies(vertex_edges_verts, 1, 
true,
 
  131                                    vertex_edges_verts_edges,
  132                                    moab::Interface::UNION);
  133        vertex_edges = subtract(vertex_edges, vertex_edges_verts_edges);
  134        vertex_edges = subtract(vertex_edges, skin_edges);
  135        vertex_edges = intersect(vertex_edges, faces_edges);
  136        CHKERR moab.add_entities(*meshset_of_edges_to_refine_ptr, vertex_edges);
 
  138      };
  139 
  140      for (int mm = 0; mm != nb_sidesets; mm++) {
  141        CHKERR refine(m_mng->getCubitMeshsetPtr(sidesets[mm], 
SIDESET),
 
  142                      meshset_of_edges_to_refine_ptr);
  143      }
  144 
  145      for (
auto m : m_mng->getCubitMeshsetPtr(
std::regex(
 
  146 
  147               (boost::format("%s") % block_name).str()
  148 
  149                   ))
  150 
  151      ) {
  152        CHKERR refine(
m, meshset_of_edges_to_refine_ptr);
 
  153      }
  154 
  155      int nb_tris;
  156      CHKERR moab.get_number_entities_by_type(*meshset_of_edges_to_refine_ptr,
 
  157                                              MBEDGE, nb_tris, true);
  158      MOFEM_LOG(
"SPLIT", Sev::inform) << 
"Refine corner edges " << nb_tris;
 
  159 
  161        CHKERR moab.write_file(
"out_edges_to_refine.vtk", 
"VTK", 
"",
 
  162                               meshset_of_edges_to_refine_ptr->get_ptr(), 1);
  163 
  165      CHKERR refine_mng->addVerticesInTheMiddleOfEdges(
 
  166          *meshset_of_edges_to_refine_ptr, bit_levels.back());
  167      CHKERR refine_mng->refineTets(0, bit_levels.back());
 
  169    }
  170 
  171    auto split = [&](auto mit) {
  173      const auto cubit_meshset = mit->getMeshset();
  174      {
  175        
  177        CHKERR bit_mng->getEntitiesByTypeAndRefLevel(bit_levels.back(),
 
  179                                                     *ref_level_meshset_ptr);
  180        CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
 
  182            *ref_level_meshset_ptr);
  183        Range ref_level_tets;
 
  184        CHKERR moab.get_entities_by_handle(*ref_level_meshset_ptr,
 
  185                                           ref_level_tets, true);
  186 
  187        
  188        CHKERR interface_ptr->getSides(cubit_meshset, bit_levels.back(), 
true,
 
  189                                       0);
  190        
  192            << "Add bit level " << bit_levels.size();
  193        bit_levels.push_back(
BitRefLevel().set(bit_levels.size()));
 
  194        
  195        CHKERR interface_ptr->splitSides(*ref_level_meshset_ptr,
 
  196                                         bit_levels.back(), cubit_meshset,
  197                                         add_prisms, true, 0);
  198      }
  199      
  202    };
  203 
  204    for (int mm = 0; mm != nb_sidesets; mm++) {
  205      CHKERR split(m_mng->getCubitMeshsetPtr(sidesets[mm], 
SIDESET));
 
  206    }
  207 
  208    for (
auto m : m_mng->getCubitMeshsetPtr(
std::regex(
 
  209 
  210             (boost::format("%s") % block_name).str()
  211 
  212                 ))
  213 
  214    ) {
  216    }
  217 
  218    if (squash_bit_levels == PETSC_TRUE) {
  219      for (unsigned int ll = 0; ll != bit_levels.size() - 1; ll++) {
  221                                              true);
  222      }
  223      CHKERR bit_mng->shiftRightBitRef(bit_levels.size() - 1);
 
  224    }
  225 
  226    
  228    auto &list = m_mng->getMeshsetsMultindex();
  229 
  230    for (
auto &
m : list) {
 
  232      CHKERR moab.get_entities_by_handle(
m.getMeshset(), ents, 
false);
 
  233      meshset_ents.merge(ents);
  234    }
  235 
  236    Range not_in_meshset_ents;
 
  237    for (EntityType type = MBEDGE; 
type != MBENTITYSET; 
type++) {
 
  239      CHKERR moab.get_entities_by_type(0, type, ents, 
false);
 
  240      not_in_meshset_ents.merge(ents);
  241    }
  242 
  243    not_in_meshset_ents = subtract(not_in_meshset_ents, meshset_ents);
  244 
  245    CHKERR moab.delete_entities(not_in_meshset_ents);
 
  246 
  250      if (squash_bit_levels)
  252      else
  253        bit = bit_levels.back();
 
  255                                                   MBTET, *meshset_ptr);
  256      CHKERR moab.write_file(
"out.vtk", 
"VTK", 
"", meshset_ptr->get_ptr(), 1);
 
  257    }
  258 
  259    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, 
MYPCOMM_INDEX);
 
  260    if (pcomm == NULL)
  262              "Communicator should be allocated");
  263 
  264    CHKERR pcomm->assign_global_ids(0, 3, 1, 
false);
 
  265    CHKERR moab.write_file(
"out.h5m");
 
  266  }
  268 
  270 
  271  return 0;
  272}
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
FTensor::Index< 'm', 3 > m
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode delete_ents_by_bit_ref(const BitRefLevel bit, const BitRefLevel mask, const bool remove_parent=false, int verb=DEFAULT_VERBOSITY, MoFEMTypes mf=MF_ZERO)=0
Delete entities from MoFEM and MOAB database.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Mesh refinement interface.
Interface for managing meshsets containing materials and boundary conditions.
Create interface from given surface and insert flat prisms in-between.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.