No Matches
/** \file mesh_cut.cpp
* \brief test for mesh cut interface
* \example mesh_cut.cpp
* \ingroup mesh_cut
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "mesh cutting\n\n";
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
PetscBool flg_myfile = PETSC_TRUE;
char mesh_file_name[255];
int surface_side_set = 200;
PetscBool flg_vol_block_set;
int vol_block_set = 1;
int edges_block_set = 1;
PetscBool flg_shift;
double shift[] = {0, 0, 0};
int nmax = 3;
int fraction_level = 2;
PetscBool squash_bits = PETSC_TRUE;
PetscBool set_coords = PETSC_TRUE;
PetscBool output_vtk = PETSC_TRUE;
int create_surface_side_set = 201;
PetscBool flg_create_surface_side_set;
int nb_ref_cut = 0;
int nb_ref_trim = 0;
PetscBool flg_tol;
double tol[] = {1e-2, 2e-1, 2e-1};
int nmax_tol = 3;
PetscBool debug = PETSC_FALSE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Mesh cut options", "none");
CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
CHKERR PetscOptionsInt("-surface_side_set", "surface side set", "",
surface_side_set, &surface_side_set, PETSC_NULL);
CHKERR PetscOptionsInt("-vol_block_set", "volume side set", "",
vol_block_set, &vol_block_set, &flg_vol_block_set);
CHKERR PetscOptionsInt("-edges_block_set", "edges block set", "",
CHKERR PetscOptionsInt("-vertex_block_set", "vertex block set", "",
CHKERR PetscOptionsRealArray("-shift", "shift surface by vector", "", shift,
&nmax, &flg_shift);
CHKERR PetscOptionsInt("-fraction_level", "fraction of merges merged", "",
fraction_level, &fraction_level, PETSC_NULL);
CHKERR PetscOptionsBool("-squash_bits", "true to squash bits at the end",
"", squash_bits, &squash_bits, PETSC_NULL);
CHKERR PetscOptionsBool("-set_coords", "true to set coords at the end", "",
set_coords, &set_coords, PETSC_NULL);
CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
output_vtk, &output_vtk, PETSC_NULL);
CHKERR PetscOptionsInt("-create_side_set", "crete side set", "",
create_surface_side_set, &create_surface_side_set,
CHKERR PetscOptionsInt("-nb_ref_cut", "nb refinements befor cut", "",
nb_ref_cut, &nb_ref_cut, PETSC_NULL);
CHKERR PetscOptionsInt("-nb_ref_trim", "nb refinements befor trim", "",
nb_ref_trim, &nb_ref_trim, PETSC_NULL);
CHKERR PetscOptionsRealArray(
"-tol", "tolerances tolCut, tolCutClose, tolTrim, tolTrimClose", "",
tol, &nmax_tol, &flg_tol);
CHKERR PetscOptionsBool("-debug", "debug (produces many VTK files)", "",
debug, &debug, PETSC_NULL);
ierr = PetscOptionsEnd();
if (flg_myfile != PETSC_TRUE) {
"*** ERROR -my_file (MESH FILE NEEDED)");
if (flg_shift && nmax != 3) {
if (flg_tol && nmax_tol != 3)
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);
const char *option;
option = "";
CHKERR moab.load_file(mesh_file_name, 0, option);
MoFEM::Core core(moab);
// get cut mesh interface
CutMeshInterface *cut_mesh;
CHKERR m_field.getInterface(cut_mesh);
// get meshset manager interface
MeshsetsManager *meshset_manager;
CHKERR m_field.getInterface(meshset_manager);
// get bit ref manager interface
BitRefManager *bit_ref_manager;
CHKERR m_field.getInterface(bit_ref_manager);
(core.getInterface<MeshsetsManager &, 0>()), SIDESET, it)) {
cout << *it << endl;
CHKERR bit_ref_manager->setBitRefLevelByType(0, MBTET,
// get surface entities form side set
if (meshset_manager->checkMeshset(surface_side_set, SIDESET))
CHKERR meshset_manager->getEntitiesByDimension(surface_side_set, SIDESET,
2, surface, true);
if (surface.empty())
// Set surface entities. If surface entities are from existing side set,
// copy those entities and do other geometrical transformations, like shift
// scale or streach, rotate.
if (meshset_manager->checkMeshset(surface_side_set, SIDESET))
CHKERR cut_mesh->copySurface(surface, NULL, shift);
"Side set not found %d", surface_side_set);
// Get geometric corner nodes and corner edges
Range fixed_edges, corner_nodes;
if (meshset_manager->checkMeshset(edges_block_set, BLOCKSET))
1, fixed_edges, true);
if (meshset_manager->checkMeshset(edges_block_set, SIDESET))
1, fixed_edges, true);
if (meshset_manager->checkMeshset(vertex_block_set, BLOCKSET))
CHKERR meshset_manager->getEntitiesByDimension(
vertex_block_set, BLOCKSET, 0, corner_nodes, true);
Range tets;
if (flg_vol_block_set) {
if (meshset_manager->checkMeshset(vol_block_set, BLOCKSET))
CHKERR meshset_manager->getEntitiesByDimension(vol_block_set, BLOCKSET,
3, tets, true);
"Block set %d not found", vol_block_set);
} else
CHKERR moab.get_entities_by_dimension(0, 3, tets, false);
CHKERR cut_mesh->setVolume(tets);
CHKERR cut_mesh->makeFront(true);
// Build tree, it is used to ask geometrical queries, i.e. to find edges to
// cut or trim.
CHKERR cut_mesh->buildTree();
// Refine mesh
CHKERR cut_mesh->refineMesh(0, nb_ref_cut, nb_ref_trim, &fixed_edges,
VERBOSE, false);
auto shift_after_ref = [&]() {
for (int ll = 1; ll != nb_ref_cut + nb_ref_trim + 1; ++ll)
CHKERR core.getInterface<BitRefManager>()->shiftRightBitRef(
nb_ref_cut + nb_ref_trim, mask, VERBOSE);
CHKERR shift_after_ref();
// Create tag storing nodal positions
double def_position[] = {0, 0, 0};
Tag th;
CHKERR moab.tag_get_handle("POSITION", 3, MB_TYPE_DOUBLE, th,
MB_TAG_CREAT | MB_TAG_SPARSE, def_position);
// Set tag values with coordinates of nodes
CHKERR cut_mesh->setTagData(th);
// Cut mesh, trim surface and merge bad edges
int first_bit = 1;
CHKERR cut_mesh->cutTrimAndMerge(first_bit, fraction_level, th, tol[0],
tol[1], tol[2], fixed_edges, corner_nodes,
true, false);
// Set coordinates for tag data
if (set_coords)
CHKERR cut_mesh->setCoords(th);
// Add tets from last level to block
if (flg_vol_block_set) {
EntityHandle meshset;
CHKERR meshset_manager->getMeshset(vol_block_set, BLOCKSET, meshset);
BitRefLevel().set(first_bit), BitRefLevel().set(), MBTET, meshset);
// Shift bits
if (squash_bits) {
BitRefLevel shift_mask;
for (int ll = 0; ll != first_bit; ++ll)
CHKERR core.getInterface<BitRefManager>()->shiftRightBitRef(
first_bit - 1, shift_mask, VERBOSE);
Range surface_verts;
CHKERR moab.get_connectivity(cut_mesh->getSurface(), surface_verts);
Range surface_edges;
CHKERR moab.get_adjacencies(cut_mesh->getSurface(), 1, false, surface_edges,
CHKERR moab.delete_entities(cut_mesh->getSurface());
CHKERR moab.delete_entities(surface_edges);
CHKERR moab.delete_entities(surface_verts);
if (flg_create_surface_side_set) {
// Check is meshset is there
create_surface_side_set, SIDESET))
CHKERR meshset_manager->addMeshset(SIDESET, create_surface_side_set);
MOFEM_LOG_C("WORLD", Sev::warning,
"Warring >>> sideset %d is on the mesh",
CHKERR meshset_manager->addEntitiesToMeshset(
SIDESET, create_surface_side_set, cut_mesh->getMergedSurfaces());
CHKERR moab.write_file("out.h5m");
if (output_vtk) {
EntityHandle meshset;
CHKERR moab.create_meshset(MESHSET_SET, meshset);
if (flg_vol_block_set) {
Range ents;
meshset_manager->getEntitiesByDimension(vol_block_set, BLOCKSET, 3,
ents, true);
CHKERR moab.add_entities(meshset, ents);
} else {
bit, BitRefLevel().set(), MBTET, meshset);
CHKERR moab.add_entities(meshset, cut_mesh->getMergedSurfaces());
CHKERR moab.write_file("out.vtk", "VTK", "", &meshset, 1);
CHKERR moab.delete_entities(&meshset, 1);
return 0;
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
static char help[]
int main()
Definition: adol-c_atom.cpp:46
static PetscErrorCode ierr
Definition: definitions.h:209
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:148
Definition: definitions.h:40
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
MoFEMErrorCode setBitRefLevelByType(const EntityHandle meshset, const EntityType type, const BitRefLevel bit, int verb=QUIET) const
Set the Bit Ref Level By Type object.
MoFEMErrorCode getEntitiesByTypeAndRefLevel(const BitRefLevel bit, const BitRefLevel mask, const EntityType type, const EntityHandle meshset, int verb=0) const
add all ents from ref level given by bit to meshset
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
MoFEMErrorCode addMeshset(const CubitBCType cubit_bc_type, const int ms_id, const std::string name="")
add cubit meshset
MoFEMErrorCode getEntitiesByDimension(const int ms_id, const unsigned int cubit_bc_type, const int dimension, Range &entities, const bool recursive=true) const
get entities from CUBIT/meshset of a particular entity dimension
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
auto bit
set bit
double tol
PetscBool output_vtk
PetscBool flg_myfile
int vertex_block_set
int edges_block_set
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.
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
Interface to cut meshes.
MoFEMErrorCode copySurface(const Range surface, Tag th=NULL, double *shift=NULL, double *origin=NULL, double *transform=NULL, const std::string save_mesh="")
copy surface entities
MoFEMErrorCode setCoords(Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
set coords from tag
const Range & getMergedSurfaces() const
MoFEMErrorCode buildTree()
build tree
const Range & getSurface() const
MoFEMErrorCode setVolume(const Range volume)
set volume entities
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
MoFEMErrorCode cutTrimAndMerge(int &first_bit, const int fraction_level, Tag th, const double tol_geometry, const double tol_cut_close, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
Cut, trim and merge.
MoFEMErrorCode makeFront(const bool debug=false)
Create front from the surface.
MoFEMErrorCode refineMesh(const int init_bit_level, const int surf_levels, const int front_levels, Range *fixed_edges=nullptr, int verb=QUIET, const bool debug=false)
Refine and set level sets.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode addEntitiesToMeshset(const CubitBCType cubit_bc_type, const int ms_id, const Range &ents)
add entities to cubit meshset
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.