No Matches
/** \file split_sideset.cpp
\brief Split sidesets
\example split_sideset.cpp
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
constexpr bool debug = false;
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
auto core_log = logging::core::get();
LogManager::createSink(LogManager::getStrmWorld(), "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",
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,
CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
output_vtk, &output_vtk, PETSC_NULL);
ierr = PetscOptionsEnd();
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) {
"*** ERROR -my_file (mesh file needed)");
if (flg_list_of_sidesets != PETSC_TRUE) {
"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 =
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;
auto update_meshsets = [&]() {
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);
// 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())
"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,
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,
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 = 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);
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())
"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,
CHKERR bit_mng->getEntitiesByTypeAndRefLevel(
bit_levels.back(), BitRefLevel().set(), MBPRISM,
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,
// set new bit level
MOFEM_LOG("SPLIT", Sev::verbose)
<< "Add bit level " << 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],
CHKERR bit_mng->shiftRightBitRef(bit_levels.size() - 1);
if (output_vtk) {
auto meshset_ptr = get_temp_meshset_ptr(moab);
if (squash_bit_levels)
bit = bit_levels[0];
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)
"Communicator should be allocated");
CHKERR pcomm->assign_global_ids(0, 3, 1, false);
CHKERR moab.write_file("out.h5m");
return 0;
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static PetscErrorCode ierr
Catch errors.
Definition: definitions.h:372
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
Definition: definitions.h:147
Definition: definitions.h:31
Definition: definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
static const bool debug
#define MOFEM_LOG(channel, severity)
Definition: LogManager.hpp:308
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
Iterator that loops over all the Cubit MeshSets in a moFEM field.
auto bit
set bit
PetscBool output_vtk
char mesh_file_name[255]
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
Managing BitRefLevels.
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 form mofem and moab database
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
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 refernce to pointer of interface.