v0.14.0
mesh_cut.cpp
Go to the documentation of this file.
1 /** \file mesh_cut.cpp
2  * \brief test for mesh cut interface
3  * \example mesh_cut.cpp
4  *
5  * \ingroup mesh_cut
6  */
7 
8 
9 
10 #include <MoFEM.hpp>
11 
12 using namespace MoFEM;
13 
14 static char help[] = "mesh cutting\n\n";
15 
16 int main(int argc, char *argv[]) {
17 
18  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
19 
20  try {
21 
22  PetscBool flg_myfile = PETSC_TRUE;
23  char mesh_file_name[255];
24  int surface_side_set = 200;
25  PetscBool flg_vol_block_set;
26  int vol_block_set = 1;
27  int edges_block_set = 1;
28  int vertex_block_set = 2;
29  PetscBool flg_shift;
30  double shift[] = {0, 0, 0};
31  int nmax = 3;
32  int fraction_level = 2;
33  PetscBool squash_bits = PETSC_TRUE;
34  PetscBool set_coords = PETSC_TRUE;
35  PetscBool output_vtk = PETSC_TRUE;
36  int create_surface_side_set = 201;
37  PetscBool flg_create_surface_side_set;
38  int nb_ref_cut = 0;
39  int nb_ref_trim = 0;
40  PetscBool flg_tol;
41  double tol[] = {1e-2, 2e-1, 2e-1};
42  int nmax_tol = 3;
43  PetscBool debug = PETSC_FALSE;
44 
45  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Mesh cut options", "none");
46  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
47  mesh_file_name, 255, &flg_myfile);
48  CHKERR PetscOptionsInt("-surface_side_set", "surface side set", "",
49  surface_side_set, &surface_side_set, PETSC_NULL);
50  CHKERR PetscOptionsInt("-vol_block_set", "volume side set", "",
51  vol_block_set, &vol_block_set, &flg_vol_block_set);
52  CHKERR PetscOptionsInt("-edges_block_set", "edges block set", "",
53  edges_block_set, &edges_block_set, PETSC_NULL);
54  CHKERR PetscOptionsInt("-vertex_block_set", "vertex block set", "",
55  vertex_block_set, &vertex_block_set, PETSC_NULL);
56  CHKERR PetscOptionsRealArray("-shift", "shift surface by vector", "", shift,
57  &nmax, &flg_shift);
58  CHKERR PetscOptionsInt("-fraction_level", "fraction of merges merged", "",
59  fraction_level, &fraction_level, PETSC_NULL);
60  CHKERR PetscOptionsBool("-squash_bits", "true to squash bits at the end",
61  "", squash_bits, &squash_bits, PETSC_NULL);
62  CHKERR PetscOptionsBool("-set_coords", "true to set coords at the end", "",
63  set_coords, &set_coords, PETSC_NULL);
64  CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
65  output_vtk, &output_vtk, PETSC_NULL);
66  CHKERR PetscOptionsInt("-create_side_set", "crete side set", "",
67  create_surface_side_set, &create_surface_side_set,
68  &flg_create_surface_side_set);
69  CHKERR PetscOptionsInt("-nb_ref_cut", "nb refinements befor cut", "",
70  nb_ref_cut, &nb_ref_cut, PETSC_NULL);
71  CHKERR PetscOptionsInt("-nb_ref_trim", "nb refinements befor trim", "",
72  nb_ref_trim, &nb_ref_trim, PETSC_NULL);
73  CHKERR PetscOptionsRealArray(
74  "-tol", "tolerances tolCut, tolCutClose, tolTrim, tolTrimClose", "",
75  tol, &nmax_tol, &flg_tol);
76  CHKERR PetscOptionsBool("-debug", "debug (produces many VTK files)", "",
77  debug, &debug, PETSC_NULL);
78  ierr = PetscOptionsEnd();
79  CHKERRG(ierr);
80 
81  if (flg_myfile != PETSC_TRUE) {
82  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
83  "*** ERROR -my_file (MESH FILE NEEDED)");
84  }
85  if (flg_shift && nmax != 3) {
86  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "three values expected");
87  }
88 
89  if (flg_tol && nmax_tol != 3)
90  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "four values expected");
91 
92  moab::Core mb_instance;
93  moab::Interface &moab = mb_instance;
94  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
95  if (pcomm == NULL)
96  pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
97 
98  const char *option;
99  option = "";
100  CHKERR moab.load_file(mesh_file_name, 0, option);
101 
102  MoFEM::Core core(moab);
103  MoFEM::CoreInterface &m_field =
104  *(core.getInterface<MoFEM::CoreInterface>());
105 
106  // get cut mesh interface
107  CutMeshInterface *cut_mesh;
108  CHKERR m_field.getInterface(cut_mesh);
109  // get meshset manager interface
110  MeshsetsManager *meshset_manager;
111  CHKERR m_field.getInterface(meshset_manager);
112  // get bit ref manager interface
113  BitRefManager *bit_ref_manager;
114  CHKERR m_field.getInterface(bit_ref_manager);
115 
116  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(
117  (core.getInterface<MeshsetsManager &, 0>()), SIDESET, it)) {
118  cout << *it << endl;
119  }
120 
121  CHKERR bit_ref_manager->setBitRefLevelByType(0, MBTET,
122  BitRefLevel().set(0));
123 
124  // get surface entities form side set
125  Range surface;
126  if (meshset_manager->checkMeshset(surface_side_set, SIDESET))
127  CHKERR meshset_manager->getEntitiesByDimension(surface_side_set, SIDESET,
128  2, surface, true);
129 
130  if (surface.empty())
131  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "No surface to cut");
132 
133  // Set surface entities. If surface entities are from existing side set,
134  // copy those entities and do other geometrical transformations, like shift
135  // scale or Stretch, rotate.
136  if (meshset_manager->checkMeshset(surface_side_set, SIDESET))
137  CHKERR cut_mesh->copySurface(surface, NULL, shift);
138  else
139  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
140  "Side set not found %d", surface_side_set);
141 
142  // Get geometric corner nodes and corner edges
143  Range fixed_edges, corner_nodes;
144  if (meshset_manager->checkMeshset(edges_block_set, BLOCKSET))
145  CHKERR meshset_manager->getEntitiesByDimension(edges_block_set, BLOCKSET,
146  1, fixed_edges, true);
147  if (meshset_manager->checkMeshset(edges_block_set, SIDESET))
148  CHKERR meshset_manager->getEntitiesByDimension(edges_block_set, SIDESET,
149  1, fixed_edges, true);
150  if (meshset_manager->checkMeshset(vertex_block_set, BLOCKSET))
151  CHKERR meshset_manager->getEntitiesByDimension(
152  vertex_block_set, BLOCKSET, 0, corner_nodes, true);
153 
154 
155  Range tets;
156  if (flg_vol_block_set) {
157  if (meshset_manager->checkMeshset(vol_block_set, BLOCKSET))
158  CHKERR meshset_manager->getEntitiesByDimension(vol_block_set, BLOCKSET,
159  3, tets, true);
160  else
161  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
162  "Block set %d not found", vol_block_set);
163  } else
164  CHKERR moab.get_entities_by_dimension(0, 3, tets, false);
165 
166  CHKERR cut_mesh->setVolume(tets);
167  CHKERR cut_mesh->makeFront(true);
168 
169  // Build tree, it is used to ask geometrical queries, i.e. to find edges to
170  // cut or trim.
171  CHKERR cut_mesh->buildTree();
172 
173  // Refine mesh
174  CHKERR cut_mesh->refineMesh(0, nb_ref_cut, nb_ref_trim, &fixed_edges,
175  VERBOSE, false);
176  auto shift_after_ref = [&]() {
177  MoFEMFunctionBegin;
178  BitRefLevel mask;
179  mask.set(0);
180  for (int ll = 1; ll != nb_ref_cut + nb_ref_trim + 1; ++ll)
181  mask.set(ll);
182  CHKERR core.getInterface<BitRefManager>()->shiftRightBitRef(
183  nb_ref_cut + nb_ref_trim, mask, VERBOSE);
184  MoFEMFunctionReturn(0);
185  };
186  CHKERR shift_after_ref();
187 
188  // Create tag storing nodal positions
189  double def_position[] = {0, 0, 0};
190  Tag th;
191  CHKERR moab.tag_get_handle("POSITION", 3, MB_TYPE_DOUBLE, th,
192  MB_TAG_CREAT | MB_TAG_SPARSE, def_position);
193  // Set tag values with coordinates of nodes
194  CHKERR cut_mesh->setTagData(th);
195 
196  // Cut mesh, trim surface and merge bad edges
197  int first_bit = 1;
198  CHKERR cut_mesh->cutTrimAndMerge(first_bit, fraction_level, th, tol[0],
199  tol[1], tol[2], fixed_edges, corner_nodes,
200  true, false);
201 
202  // Set coordinates for tag data
203  if (set_coords)
204  CHKERR cut_mesh->setCoords(th);
205 
206  // Add tets from last level to block
207  if (flg_vol_block_set) {
208  EntityHandle meshset;
209  CHKERR meshset_manager->getMeshset(vol_block_set, BLOCKSET, meshset);
210  CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
211  BitRefLevel().set(first_bit), BitRefLevel().set(), MBTET, meshset);
212  }
213 
214  // Shift bits
215  if (squash_bits) {
216  BitRefLevel shift_mask;
217  for (int ll = 0; ll != first_bit; ++ll)
218  shift_mask.set(ll);
219  CHKERR core.getInterface<BitRefManager>()->shiftRightBitRef(
220  first_bit - 1, shift_mask, VERBOSE);
221  }
222 
223  Range surface_verts;
224  CHKERR moab.get_connectivity(cut_mesh->getSurface(), surface_verts);
225  Range surface_edges;
226  CHKERR moab.get_adjacencies(cut_mesh->getSurface(), 1, false, surface_edges,
227  moab::Interface::UNION);
228  CHKERR moab.delete_entities(cut_mesh->getSurface());
229  CHKERR moab.delete_entities(surface_edges);
230  CHKERR moab.delete_entities(surface_verts);
231 
232  if (flg_create_surface_side_set) {
233  // Check is meshset is there
234  if (!core.getInterface<MeshsetsManager>()->checkMeshset(
235  create_surface_side_set, SIDESET))
236  CHKERR meshset_manager->addMeshset(SIDESET, create_surface_side_set);
237  else
238  MOFEM_LOG_C("WORLD", Sev::warning,
239  "Warring >>> sideset %d is on the mesh",
240  create_surface_side_set);
241 
242  CHKERR meshset_manager->addEntitiesToMeshset(
243  SIDESET, create_surface_side_set, cut_mesh->getMergedSurfaces());
244  }
245 
246  CHKERR moab.write_file("out.h5m");
247 
248  if (output_vtk) {
249  EntityHandle meshset;
250  CHKERR moab.create_meshset(MESHSET_SET, meshset);
251  if (flg_vol_block_set) {
252  Range ents;
253  meshset_manager->getEntitiesByDimension(vol_block_set, BLOCKSET, 3,
254  ents, true);
255  CHKERR moab.add_entities(meshset, ents);
256  } else {
257  BitRefLevel bit = BitRefLevel().set(0);
258  CHKERR bit_ref_manager->getEntitiesByTypeAndRefLevel(
259  bit, BitRefLevel().set(), MBTET, meshset);
260  }
261  CHKERR moab.add_entities(meshset, cut_mesh->getMergedSurfaces());
262  CHKERR moab.write_file("out.vtk", "VTK", "", &meshset, 1);
263  CHKERR moab.delete_entities(&meshset, 1);
264  }
265  }
266  CATCH_ERRORS;
267 
268  MoFEM::Core::Finalize();
269 
270  return 0;
271 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
SIDESET
@ SIDESET
Definition: definitions.h:160
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::MeshsetsManager::addEntitiesToMeshset
MoFEMErrorCode addEntitiesToMeshset(const CubitBCType cubit_bc_type, const int ms_id, const Range &ents)
add entities to cubit meshset
Definition: MeshsetsManager.cpp:418
surface
Definition: surface.py:1
EntityHandle
MoFEM::CutMeshInterface::setVolume
MoFEMErrorCode setVolume(const Range volume)
set volume entities
Definition: CutMeshInterface.cpp:108
MoFEM::BitRefManager::setBitRefLevelByType
MoFEMErrorCode setBitRefLevelByType(const EntityHandle meshset, const EntityType type, const BitRefLevel bit, int verb=QUIET) const
Set the Bit Ref Level By Type object.
Definition: BitRefManager.cpp:512
MoFEM::CutMeshInterface::setCoords
MoFEMErrorCode setCoords(Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
set coords from tag
Definition: CutMeshInterface.cpp:3130
MoFEM::CutMeshInterface::makeFront
MoFEMErrorCode makeFront(const bool debug=false)
Create front from the surface.
Definition: CutMeshInterface.cpp:450
MoFEM::CutMeshInterface::setTagData
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
Definition: CutMeshInterface.cpp:3114
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::CutMeshInterface::copySurface
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
Definition: CutMeshInterface.cpp:48
flg_myfile
PetscBool flg_myfile
Definition: mesh_smoothing.cpp:21
MoFEM.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::MeshsetsManager::getMeshset
MoFEMErrorCode getMeshset(const int ms_id, const unsigned int cubit_bc_type, EntityHandle &meshset) const
get meshset from CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:708
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
VERBOSE
@ VERBOSE
Definition: definitions.h:222
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CutMeshInterface::cutTrimAndMerge
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.
Definition: CutMeshInterface.cpp:379
vertex_block_set
int vertex_block_set
Definition: mesh_smoothing.cpp:24
MoFEM::CutMeshInterface
Interface to cut meshes.
Definition: CutMeshInterface.hpp:19
output_vtk
PetscBool output_vtk
Definition: mesh_smoothing.cpp:25
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CutMeshInterface::getMergedSurfaces
const Range & getMergedSurfaces() const
Definition: CutMeshInterface.hpp:483
MoFEM::MeshsetsManager::getEntitiesByDimension
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
Definition: MeshsetsManager.cpp:669
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::MeshsetsManager::addMeshset
MoFEMErrorCode addMeshset(const CubitBCType cubit_bc_type, const int ms_id, const std::string name="")
add cubit meshset
Definition: MeshsetsManager.cpp:385
surface.surface
def surface(x, y, z, eta)
Definition: surface.py:3
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
MoFEM::CutMeshInterface::buildTree
MoFEMErrorCode buildTree()
build tree
Definition: CutMeshInterface.cpp:215
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
Range
MoFEM::CoreTmp< 0 >::Initialize
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
edges_block_set
int edges_block_set
Definition: mesh_smoothing.cpp:23
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition: MeshsetsManager.hpp:71
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
MoFEM::CutMeshInterface::refineMesh
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.
Definition: CutMeshInterface.cpp:809
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface
Interface.
Definition: Interface.hpp:27
MoFEM::BitRefManager::getEntitiesByTypeAndRefLevel
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
Definition: BitRefManager.cpp:734
MoFEM::MeshsetsManager::checkMeshset
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:360
help
static char help[]
Definition: mesh_cut.cpp:14
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
MoFEM::CutMeshInterface::getSurface
const Range & getSurface() const
Definition: CutMeshInterface.hpp:466
main
int main(int argc, char *argv[])
Definition: mesh_cut.cpp:16
tol
double tol
Definition: mesh_smoothing.cpp:26
Generated by Doxygen 1.8.17 and hosted at University of Glasgow