v0.13.2
Loading...
Searching...
No Matches
Classes | Functions | Variables
mesh_cut_test.cpp File Reference
MoFEM » CutMeshInterface

test for mesh cut interface More...

#include <MoFEM.hpp>

Go to the source code of this file.

Classes

struct  TestBitLevel
 

Functions

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

Variables

static char help [] = "testing mesh cut test\n\n"
 

Detailed Description

test for mesh cut interface

Definition in file mesh_cut_test.cpp.

Function Documentation

◆ main()

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

Definition at line 32 of file mesh_cut_test.cpp.

32 {
33
34 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
35
36 try {
37
38 MOFEM_LOG_CHANNEL("WORLD");
39
40 PetscBool flg = PETSC_TRUE;
41 char mesh_file_name[255];
42 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
43 255, &flg);
44
45 if (flg != PETSC_TRUE)
46 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
47 "*** ERROR -my_file (MESH FILE NEEDED)");
48
49 int side_set = 200;
50 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-side_set", &side_set,
51 PETSC_NULL);
52
53 int restricted_side_set = 205;
54 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-restricted_side_set",
55 &restricted_side_set, PETSC_NULL);
56
57 double shift[] = {0, 0, 0};
58 int nmax = 3;
59 CHKERR PetscOptionsGetRealArray("", "-shift", shift, &nmax, &flg);
60 if (flg && nmax != 3)
61 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
62 "three values expected");
63
64 int fixed_edges_blockset = 100;
65 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-fixed_edges_blockset",
66 &fixed_edges_blockset, PETSC_NULL);
67
68 int corner_nodes_blockset = 1;
69 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-corner_nodes_blockset",
70 &corner_nodes_blockset, PETSC_NULL);
71
72 double tol[] = {0, 0, 0};
73 int nmax_tol = 3;
74 CHKERR PetscOptionsGetRealArray("", "-tol", tol, &nmax_tol, &flg);
75 if (flg && nmax_tol != 3)
76 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "four values expected");
77
78 PetscBool test = PETSC_FALSE;
79 CHKERR
80 PetscOptionsGetBool(PETSC_NULL, "", "-test", &test, PETSC_NULL);
81
82 moab::Core mb_instance;
83 moab::Interface &moab = mb_instance;
84 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
85 auto moab_comm_wrap =
86 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
87 if (pcomm == NULL)
88 pcomm =
89 new ParallelComm(&moab, moab_comm_wrap->get_comm());
90
91 const char *option;
92 option = "";
93 CHKERR moab.load_file(mesh_file_name, 0, option);
94
95 MoFEM::Core core(moab);
96 MoFEM::CoreInterface &m_field =
97 *(core.getInterface<MoFEM::CoreInterface>());
98
99 MOFEM_LOG_CHANNEL("WORLD");
100 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(
101 (*core.getInterface<MeshsetsManager>()), SIDESET, it)) {
102 MOFEM_LOG("WORLD", Sev::inform) << *it;
103 }
104
105 BitRefLevel bit_level0;
106 bit_level0.set(0);
107 BitRefLevel bit_last;
108 bit_last.set(BITREFLEVEL_SIZE - 1);
109
110 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
111 0, 3, BitRefLevel().set(0));
112 CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevelByDim(0, 3,
113 bit_last);
114
115 CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
116 bit_level0, BitRefLevel().set(), MBTET, "out_tets_init_level0.vtk",
117 "VTK", "");
118 CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
119 bit_last, BitRefLevel().set(), MBTET, "out_tets_bit_last.vtk", "VTK",
120 "");
121
122 int no_of_ents_not_in_database = -1;
123 Range ents_not_in_database;
124 if (test) {
125 core.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
126 ents_not_in_database);
127 no_of_ents_not_in_database = ents_not_in_database.size();
128 }
129
130 // Get BitRefManager interface,,
131 BitRefManager *bit_ref_manager;
132 CHKERR m_field.getInterface(bit_ref_manager);
133 // get meshset manager interface
134 MeshsetsManager *meshset_manager;
135 CHKERR m_field.getInterface(meshset_manager);
136 // get cut mesh interface
137 CutMeshInterface *cut_mesh;
138 CHKERR m_field.getInterface(cut_mesh);
139
140 // Get geometric corner nodes and corner edges
141 Range fixed_edges, corner_nodes;
142 if (meshset_manager->checkMeshset(fixed_edges_blockset, SIDESET)) {
143 CHKERR meshset_manager->getEntitiesByDimension(
144 fixed_edges_blockset, SIDESET, 1, fixed_edges, true);
145 }
146 if (meshset_manager->checkMeshset(fixed_edges_blockset, BLOCKSET)) {
147 CHKERR meshset_manager->getEntitiesByDimension(
148 fixed_edges_blockset, BLOCKSET, 1, fixed_edges, true);
149 }
150 if (meshset_manager->checkMeshset(corner_nodes_blockset, BLOCKSET)) {
151 CHKERR meshset_manager->getEntitiesByDimension(
152 corner_nodes_blockset, BLOCKSET, 0, corner_nodes, true);
153 }
154
155 // get surface entities form side set
156 Range restriced_surface;
157 if (meshset_manager->checkMeshset(restricted_side_set, SIDESET))
158 CHKERR meshset_manager->getEntitiesByDimension(
159 restricted_side_set, SIDESET, 2, restriced_surface, true);
160
161 Range surface;
162 if (meshset_manager->checkMeshset(side_set, SIDESET))
163 CHKERR meshset_manager->getEntitiesByDimension(side_set, SIDESET, 2,
164 surface, true);
165 if (surface.empty())
166 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID, "No surface to cut");
167 // Set surface entities. If surface entities are from existing side set,
168 // copy those entities and do other geometrical transformations, like shift
169 // scale or streach, rotate.
170 if (meshset_manager->checkMeshset(side_set, SIDESET))
171 CHKERR cut_mesh->copySurface(surface, NULL, shift, NULL, NULL,
172 "surface.vtk");
173 else
174 CHKERR cut_mesh->setSurface(surface);
175
176 cut_mesh->setConstrainSurface(restriced_surface);
177
178 Range tets;
179 CHKERR moab.get_entities_by_dimension(0, 3, tets, false);
180
181 CHKERR cut_mesh->setVolume(tets);
182 CHKERR cut_mesh->buildTree();
183 CHKERR cut_mesh->makeFront(true);
184 const int nb_ref_cut = 1;
185 const int nb_ref_trim = 1;
186 CHKERR cut_mesh->refineMesh(0, nb_ref_cut, nb_ref_trim, &fixed_edges,
187 VERBOSE, false);
188 auto shift_after_ref = [&]() {
189 MoFEMFunctionBegin;
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);
196 MoFEMFunctionReturn(0);
197 };
198 CHKERR shift_after_ref();
199
200 CHKERR m_field.getInterface<BitRefManager>()
201 ->writeEntitiesAllBitLevelsByType(BitRefLevel().set(), MBTET,
202 "all_bits.vtk", "VTK", "");
203
204 // Create tag storing nodal positions
205 double def_position[] = {0, 0, 0};
206 Tag th;
207 CHKERR moab.tag_get_handle("POSITION", 3, MB_TYPE_DOUBLE, th,
208 MB_TAG_CREAT | MB_TAG_SPARSE, def_position);
209 // Set tag values with coordinates of nodes
210 CHKERR cut_mesh->setTagData(th);
211
212 // Cut mesh, trim surface and merge bad edges
213 int first_bit = 1;
214 CHKERR cut_mesh->cutTrimAndMerge(first_bit, 5, th, tol[0], tol[1], tol[2],
215 fixed_edges, corner_nodes, true, false);
216
217 // Split faces
218 CHKERR cut_mesh->splitSides(BitRefLevel().set(first_bit),
219 BitRefLevel().set(first_bit + 1),
220 cut_mesh->getMergedSurfaces(), th);
221 CHKERR core.getInterface<MeshsetsManager>()
222 ->updateAllMeshsetsByEntitiesChildren(BitRefLevel().set(first_bit + 1));
223
224 CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
225 BitRefLevel().set(first_bit + 1), BitRefLevel().set(), MBTET,
226 "out_split_tets.vtk", "VTK", "");
227
228 CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
229 BitRefLevel().set(first_bit + 1), BitRefLevel().set(), MBPRISM,
230 "out_split_prism.vtk", "VTK", "");
231
232 if (test) {
233 for (int ll = 0; ll != first_bit + 2; ++ll)
234 CHKERR TestBitLevel(core.getInterface<BitRefManager>())(
235 BitRefLevel().set(ll), -1);
236 }
237
238 // Finally shift bits
239 BitRefLevel shift_mask;
240 for (int ll = 0; ll != first_bit; ++ll)
241 shift_mask.set(ll);
242 CHKERR core.getInterface<BitRefManager>()->shiftRightBitRef(
243 first_bit, shift_mask, VERBOSE);
244
245 // Set coordinates for tag data
246 CHKERR cut_mesh->setCoords(th);
247 CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
248 BitRefLevel().set(first_bit), BitRefLevel().set(), MBTET,
249 "out_tets_shift_level0.vtk", "VTK", "");
250 CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
251 BitRefLevel().set(first_bit + 1), BitRefLevel().set(), MBTET,
252 "out_tets_shift_level1.vtk", "VTK", "");
253 CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
254 BitRefLevel().set(first_bit + 1), BitRefLevel().set(), MBPRISM,
255 "out_tets_shift_level1_prism.vtk", "VTK", "");
256
257 Range surface_verts;
258 CHKERR moab.get_connectivity(cut_mesh->getSurface(), surface_verts);
259 Range adj_surface_edges;
260 CHKERR moab.get_adjacencies(cut_mesh->getSurface(), 1, false,
261 adj_surface_edges, moab::Interface::UNION);
262 CHKERR moab.delete_entities(cut_mesh->getSurface());
263 CHKERR moab.delete_entities(adj_surface_edges);
264 CHKERR moab.delete_entities(surface_verts);
265 CHKERR m_field.delete_ents_by_bit_ref(
266 BitRefLevel().set(0) | BitRefLevel().set(1),
267 BitRefLevel().set(0) | BitRefLevel().set(1), true, VERBOSE);
268
269 {
270 EntityHandle meshset;
271 CHKERR moab.create_meshset(MESHSET_SET, meshset);
272 Range tets;
273 CHKERR moab.get_entities_by_dimension(0, 3, tets, true);
274 CHKERR moab.add_entities(meshset, tets);
275 CHKERR moab.write_file("out.vtk", "VTK", "", &meshset, 1);
276 CHKERR moab.delete_entities(&meshset, 1);
277 }
278 CHKERR core.getInterface<BitRefManager>()->writeBitLevelByType(
279 bit_last, BitRefLevel().set(), MBTET, "out_tets_bit_last.vtk", "VTK",
280 "");
281 CHKERR core.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
282 "left_entities.vtk", "VTK", "");
283
284 MOFEM_LOG_CHANNEL("WORLD");
285 if (test) {
286 Range ents;
287 core.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(ents);
288 if (no_of_ents_not_in_database != static_cast<int>(ents.size())) {
289 MOFEM_LOG("WORLD", Sev::inform) << subtract(ents, ents_not_in_database);
290 EntityHandle meshset;
291 CHKERR moab.create_meshset(MESHSET_SET, meshset);
292 Range tets;
293 CHKERR moab.get_entities_by_dimension(0, 3, tets, true);
294 CHKERR moab.add_entities(meshset, subtract(ents, ents_not_in_database));
295 CHKERR moab.write_file("not_cleanded.vtk", "VTK", "", &meshset, 1);
296 CHKERR moab.delete_entities(&meshset, 1);
297 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
298 "Inconsistent number of ents %d!=%d",
299 no_of_ents_not_in_database, ents.size());
300 }
301 }
302 }
303 CATCH_ERRORS;
304
305 MoFEM::Core::Finalize();
306
307 return 0;
308}
VERBOSE
@ VERBOSE
Definition: definitions.h:209
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
BITREFLEVEL_SIZE
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:219
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
SIDESET
@ SIDESET
Definition: definitions.h:147
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
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:272
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:581
_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
help
static char help[]
Definition: mesh_cut_test.cpp:13
tol
double tol
Definition: mesh_smoothing.cpp:27
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
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::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::CoreInterface::delete_ents_by_bit_ref
virtual MoFEMErrorCode delete_ents_by_bit_ref(const BitRefLevel bit, const BitRefLevel mask, const bool remove_parent=false, int verb=DEFAULT_VERBOSITY, MoFEMTypes mf=MF_ZERO)=0
delete entities form mofem and moab database
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::CutMeshInterface
Interface to cut meshes.
Definition: CutMeshInterface.hpp:19
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
MoFEM::CutMeshInterface::setCoords
MoFEMErrorCode setCoords(Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
set coords from tag
Definition: CutMeshInterface.cpp:3142
MoFEM::CutMeshInterface::getMergedSurfaces
const Range & getMergedSurfaces() const
Definition: CutMeshInterface.hpp:483
MoFEM::CutMeshInterface::setSurface
MoFEMErrorCode setSurface(const Range surface)
set surface entities
Definition: CutMeshInterface.cpp:42
MoFEM::CutMeshInterface::buildTree
MoFEMErrorCode buildTree()
build tree
Definition: CutMeshInterface.cpp:215
MoFEM::CutMeshInterface::splitSides
MoFEMErrorCode splitSides(const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
split sides
Definition: CutMeshInterface.cpp:2278
MoFEM::CutMeshInterface::getSurface
const Range & getSurface() const
Definition: CutMeshInterface.hpp:466
MoFEM::CutMeshInterface::setConstrainSurface
MoFEMErrorCode setConstrainSurface(const Range surf)
Set the constrain surface object.
Definition: CutMeshInterface.cpp:114
MoFEM::CutMeshInterface::setVolume
MoFEMErrorCode setVolume(const Range volume)
set volume entities
Definition: CutMeshInterface.cpp:108
MoFEM::CutMeshInterface::setTagData
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
Definition: CutMeshInterface.cpp:3126
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
MoFEM::CutMeshInterface::makeFront
MoFEMErrorCode makeFront(const bool debug=false)
Create front from the surface.
Definition: CutMeshInterface.cpp:450
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
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93

Variable Documentation

◆ help

char help[] = "testing mesh cut test\n\n"
static

Definition at line 13 of file mesh_cut_test.cpp.

Generated by Doxygen 1.9.5 and hosted at University of Glasgow