v0.8.13
mesh_cut.cpp
/** \file mesh_cut.cpp
* \brief test for mesh cut interface
* \example mesh_cut.cpp
*
* \ingroup mesh_cut
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#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;
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 side set", "",
CHKERR PetscOptionsInt("-vertex_block_set", "vertex side 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,
&flg_create_surface_side_set);
ierr = PetscOptionsEnd(); CHKERRG(ierr);
if (flg_myfile != PETSC_TRUE) {
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"*** ERROR -my_file (MESH FILE NEEDED)");
}
if (flg_shift && nmax != 3) {
SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
"three values expected");
}
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 = ""; //"PARALLEL=BCAST";//;DEBUG_IO";
CHKERR moab.load_file(mesh_file_name, 0, option);
// 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;
}
BitRefLevel bit_level0;
bit_level0.set(0);
CHKERR bit_ref_manager->setBitRefLevelByType(0, MBTET, bit_level0);
// get surface entities form side set
Range surface;
if (meshset_manager->checkMeshset(surface_side_set, SIDESET)) {
CHKERR meshset_manager->getEntitiesByDimension(surface_side_set, SIDESET,
2, surface, true);
}
if (surface.empty()) {
SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "No surface to cut");
}
// 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);
} else {
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Side set not found %d", surface_side_set);
}
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);
} else {
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Block set %d not found", vol_block_set);
}
} else {
CHKERR moab.get_entities_by_dimension(0, 3, tets, false);
}
CHKERR cut_mesh->setVolume(tets);
// Build tree, it is used to ask geometrical queries, i.e. to find edges to
// cut or trim.
CHKERR cut_mesh->buildTree();
BitRefLevel bit_level1; // Cut level
bit_level1.set(1);
BitRefLevel bit_level2; // Trim level
bit_level2.set(2);
BitRefLevel bit_level3; // Merge level
bit_level3.set(3);
BitRefLevel bit_level4; // TetGen level
bit_level4.set(4);
// 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);
// 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(vertex_block_set, BLOCKSET)) {
CHKERR meshset_manager->getEntitiesByDimension(
vertex_block_set, BLOCKSET, 0, corner_nodes, true);
}
// Cut mesh, trim surface and merge bad edges
CHKERR cut_mesh->cutTrimAndMerge(fraction_level, bit_level1, bit_level2,
bit_level3, th, 1e-2, 1e-1, 1e-1, 1e-2,
fixed_edges, corner_nodes, true, true);
// Improve mesh with tetgen
#ifdef WITH_TETGEN
// Switches controling TetGen
vector<string> switches;
switches.push_back("rp180YsqORJS0VV");
CHKERR cut_mesh->rebuildMeshWithTetGen(switches, bit_level3, bit_level4,
cut_mesh->getMergedSurfaces(),
fixed_edges, corner_nodes, th);
#else
bit_level4 = bit_level3;
const_cast<Range &>(cut_mesh->getTetgenSurfaces()) =
cut_mesh->getMergedSurfaces();
#endif // WITH_TETGEN
// Add tets from last level to block
if(flg_vol_block_set) {
EntityHandle meshset;
CHKERR meshset_manager->getMeshset(vol_block_set, BLOCKSET, meshset);
bit_level4, BitRefLevel().set(), MBTET, meshset);
}
// Shift bits
if (squash_bits) {
BitRefLevel shift_mask;
shift_mask.set(0);
shift_mask.set(1);
shift_mask.set(2);
shift_mask.set(3);
shift_mask.set(4);
CHKERR core.getInterface<BitRefManager>()->shiftRightBitRef(4, shift_mask,
}
if (set_coords) {
// Set coordinates for tag data
CHKERR cut_mesh->setCoords(th);
}
Range surface_verts;
CHKERR moab.get_connectivity(cut_mesh->getSurface(), surface_verts);
CHKERR moab.delete_entities(cut_mesh->getSurface());
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);
} else {
PetscPrintf(PETSC_COMM_SELF,
"<<< Warring >>> sideset %d is on the mesh\n",
create_surface_side_set);
}
CHKERR meshset_manager->addEntitiesToMeshset(
SIDESET, create_surface_side_set, cut_mesh->getTetgenSurfaces());
}
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 {
if (squash_bits)
bit = bit_level0;
else
bit = bit_level4;
bit_level4, BitRefLevel().set(), MBTET, meshset);
}
CHKERR moab.add_entities(meshset, cut_mesh->getTetgenSurfaces());
CHKERR moab.write_file("out.vtk","VTK","",&meshset,1);
CHKERR moab.delete_entities(&meshset,1);
}
}
return 0;
}