v0.14.0
Loading...
Searching...
No Matches
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[] )

Definition at line 16 of file mesh_cut.cpp.

16 {
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",
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", "",
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();
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
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
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))
146 1, fixed_edges, true);
147 if (meshset_manager->checkMeshset(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 = [&]() {
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);
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 }
267
269
270 return 0;
271}
#define MOFEM_LOG_C(channel, severity, format,...)
@ VERBOSE
#define CATCH_ERRORS
Catch errors.
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ SIDESET
@ BLOCKSET
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_INVALID_DATA
Definition definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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
#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.
auto bit
set bit
static char help[]
Definition mesh_cut.cpp:14
double tol
PetscBool output_vtk
PetscBool flg_myfile
int vertex_block_set
int edges_block_set
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
static const bool debug
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 reference to pointer of interface.

Variable Documentation

◆ help

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

Definition at line 14 of file mesh_cut.cpp.