v0.14.0
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 stretch, 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 }

Variable Documentation

◆ help

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

Definition at line 13 of file mesh_cut_test.cpp.

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
surface
Definition: surface.py:1
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
MoFEM::CutMeshInterface::setVolume
MoFEMErrorCode setVolume(const Range volume)
set volume entities
Definition: CutMeshInterface.cpp:108
MoFEM::CutMeshInterface::setCoords
MoFEMErrorCode setCoords(Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
set coords from tag
Definition: CutMeshInterface.cpp:3130
help
static char help[]
Definition: mesh_cut_test.cpp:13
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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
MoFEM::CutMeshInterface
Interface to cut meshes.
Definition: CutMeshInterface.hpp:19
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
surface.surface
def surface(x, y, z, eta)
Definition: surface.py:3
MoFEM::CutMeshInterface::setConstrainSurface
MoFEMErrorCode setConstrainSurface(const Range surf)
Set the constrain surface object.
Definition: CutMeshInterface.cpp:114
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::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
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
Range
MoFEM::PetscOptionsGetRealArray
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
Definition: DeprecatedPetsc.hpp:192
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_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
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
BITREFLEVEL_SIZE
#define BITREFLEVEL_SIZE
max number of refinements
Definition: definitions.h:232
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::CutMeshInterface::splitSides
MoFEMErrorCode splitSides(const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
split sides
Definition: CutMeshInterface.cpp:2266
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
TestBitLevel
Definition: mesh_cut_test.cpp:15
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::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
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::CutMeshInterface::getSurface
const Range & getSurface() const
Definition: CutMeshInterface.hpp:466
tol
double tol
Definition: mesh_smoothing.cpp:26
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow