v0.9.0
Functions | Variables
mesh_cut.cpp File Reference

test for mesh cut interface More...

#include <MoFEM.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "mesh cutting\n\n"
 

Detailed Description

test for mesh cut interface

Definition in file mesh_cut.cpp.

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
mesh_cut.cpp.

Definition at line 28 of file mesh_cut.cpp.

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

Variable Documentation

◆ help

char help[] = "mesh cutting\n\n"
static
Examples
mesh_cut.cpp.

Definition at line 26 of file mesh_cut.cpp.