v0.14.0
Functions | Variables
tetgen_interface.cpp File Reference
MoFEM » TetGen interface

test for TetGen interface More...

#include <MoFEM.hpp>

Go to the source code of this file.

Functions

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

Variables

static char help [] = "...\n\n"
 
static int debug = 1
 

Detailed Description

test for TetGen interface

Definition in file tetgen_interface.cpp.

Function Documentation

◆ main()

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

Definition at line 14 of file tetgen_interface.cpp.

14  {
15 
16  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
17 
18  try {
19 
20  moab::Core mb_instance;
21  moab::Interface &moab = mb_instance;
22  int rank;
23  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
24 
25  // Read parameters from line command
26  PetscBool flg = PETSC_TRUE;
27  char mesh_file_name[255];
28 #if PETSC_VERSION_GE(3, 6, 4)
29  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
30  255, &flg);
31 #else
32  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
33  mesh_file_name, 255, &flg);
34 #endif
35  if (flg != PETSC_TRUE) {
36  SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
37  }
38 
39  // Read mesh to MOAB
40  const char *option;
41  option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
42  CHKERR moab.load_file(mesh_file_name, 0, option);
43 
44  // Create MoFEM (Joseph) databas
45  MoFEM::Core core(moab);
46  MoFEM::Interface &m_field = core;
47 
48  BitRefLevel bit_level0;
49  bit_level0.set(0);
50  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
51  0, 3, bit_level0);
52 
53  Range nodes, tets;
54  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes, false);
55  CHKERR moab.get_entities_by_type(0, MBTET, tets, false);
56 
57  // mapping between MoAB and TetGen
58  tetgenio in, out;
59  std::map<EntityHandle, unsigned long> moab_tetgen_map;
60  std::map<unsigned long, EntityHandle> tetgen_moab_map;
61 
62  // get TetGen interface
63  TetGenInterface *tetgen_iface;
64  CHKERR m_field.getInterface(tetgen_iface);
65 
66  // set MoAB nodes to TetGen data structure
67  CHKERR tetgen_iface->inData(nodes, in, moab_tetgen_map, tetgen_moab_map);
68  std::map<int, Range> types_ents;
69  types_ents[TetGenInterface::RIDGEVERTEX].merge(nodes);
70  CHKERR tetgen_iface->setGeomData(in, moab_tetgen_map, tetgen_moab_map,
71  types_ents);
72 
73  char switches[] = ""; // TetGen switches, look to TegGen user manual
74  CHKERR tetgen_iface->tetRahedralize(switches, in,
75  out); // make a TetGen mesh
76  BitRefLevel bit_level1;
77  bit_level1.set(1);
78  // get mesh from TetGen and set bit level to 1
79  CHKERR tetgen_iface->outData(in, out, moab_tetgen_map, tetgen_moab_map,
80  bit_level1, false, true);
81 
82  // save results
83  EntityHandle meshset_level1;
84  CHKERR moab.create_meshset(MESHSET_SET, meshset_level1);
85  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
86  bit_level1, BitRefLevel().set(), MBTET, meshset_level1);
87  if (debug)
88  CHKERR moab.write_file("level1.vtk", "VTK", "", &meshset_level1, 1);
89 
90  // clean data
91  in.deinitialize();
92  out.deinitialize();
93  // initialize TetGen data structures
94  in.initialize();
95  out.initialize();
96 
97  // clear data stratus used by MoFEM
98  moab_tetgen_map.clear();
99  tetgen_moab_map.clear();
100 
101  // get mesh skin
102  Skinner skin(&moab);
103  Range outer_surface_skin;
104  CHKERR skin.find_skin(0, tets, false, outer_surface_skin);
105 
106  // set nodes to TetGen data structures
107  CHKERR tetgen_iface->inData(nodes, in, moab_tetgen_map, tetgen_moab_map);
108 
109  Range side_set_faces;
110  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, SIDESET, sit)) {
111  Range faces;
112  CHKERR moab.get_entities_by_type(sit->meshset, MBTRI, faces, true);
113  side_set_faces.merge(faces);
114  }
115 
116  outer_surface_skin = subtract(outer_surface_skin, side_set_faces);
117  std::vector<std::vector<Range>> sorted_outer_surface_skin;
118  int nb_ents = outer_surface_skin.size();
119  // sort surface elements in planar groups
120  CHKERR tetgen_iface->groupRegion_Triangle(outer_surface_skin,
121  sorted_outer_surface_skin, 1e-10);
122  if (debug > 0) {
123  std::cout << " number of triangle entities " << nb_ents
124  << " number of faces disjoint regions "
125  << sorted_outer_surface_skin.size() << std::endl;
126  for (unsigned int vv = 0; vv < sorted_outer_surface_skin.size(); vv++) {
127  std::cout << "\tnb of disjoint region " << vv
128  << " nb of no-planar subregions "
129  << sorted_outer_surface_skin[vv].size() << std::endl;
130  for (unsigned int vvv = 0; vvv < sorted_outer_surface_skin[vv].size();
131  vvv++) {
132  std::cout << "\t\tnb. of subregion " << vvv
133  << " nb. elements in subregion "
134  << sorted_outer_surface_skin[vv][vvv].size() << std::endl;
135  }
136  }
137  }
138 
139  std::vector<std::pair<Range, int>>
140  markers; // set markers to surface elements
141  std::vector<std::vector<Range>>::iterator vit =
142  sorted_outer_surface_skin.begin();
143  for (; vit != sorted_outer_surface_skin.end(); vit++) {
144  std::vector<Range>::iterator viit = vit->begin();
145  for (; viit != vit->end(); viit++) {
146  Range polygons;
147  CHKERR tetgen_iface->makePolygonFacet(*viit, polygons);
148  Range aa;
149  CHKERR moab.get_connectivity(*viit, aa, true);
150  Range viit_edges;
151  CHKERR m_field.get_moab().get_adjacencies(*viit, 1, false, viit_edges,
152  moab::Interface::UNION);
153  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, SIDESET, sit)) {
154  Range faces;
155  CHKERR moab.get_entities_by_type(sit->meshset, MBTRI, faces, true);
156  Range faces_edges;
157  CHKERR m_field.get_moab().get_adjacencies(
158  faces, 1, false, faces_edges, moab::Interface::UNION);
159  aa.merge(intersect(faces_edges, viit_edges));
160  }
161  markers.push_back(std::pair<Range, int>(unite(polygons, aa), -1));
162  // markers.push_back(std::pair<Range,int>(polygons,-1));
163  }
164  }
165 
166  // set markers to side Cubit sets
167  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, SIDESET, sit)) {
168  int id = sit->getMeshsetId();
169  Range faces;
170  CHKERR moab.get_entities_by_type(sit->meshset, MBTRI, faces, true);
171  markers.push_back(std::pair<Range, int>(faces, id));
172  }
173 
174  // CHKERR tetgen_iface->inData(nodes,in,moab_tetgen_map,tetgen_moab_map);
175  // set face markers
176  CHKERR tetgen_iface->setFaceData(markers, in, moab_tetgen_map,
177  tetgen_moab_map);
178 
179  std::vector<std::pair<EntityHandle, int>> regions; // list of regions
180  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
181  int id = bit->getMeshsetId();
182  Range tets;
183  CHKERR moab.get_entities_by_type(bit->meshset, MBTET, tets, true);
184  regions.push_back(std::pair<EntityHandle, int>(*tets.begin(), -id));
185  }
186  // set volume regions
187  CHKERR tetgen_iface->setRegionData(regions, in);
188 
189  // print mesh in tetgeb format
190  if (debug > 0) {
191  char tetgen_in_file_name[] = "in";
192  in.save_nodes(tetgen_in_file_name);
193  in.save_elements(tetgen_in_file_name);
194  in.save_faces(tetgen_in_file_name);
195  in.save_edges(tetgen_in_file_name);
196  in.save_poly(tetgen_in_file_name);
197  }
198 
199  // make TetGen mesh
200  char switches2[] = "pYA"; // TetGen line command switches, look to TetGen
201  // manual for details
202  CHKERR tetgen_iface->tetRahedralize(switches2, in, out); // make a mesh
203  BitRefLevel bit_level2;
204  bit_level2.set(2);
205  // get mesh form TetGen and store it on second bit refined level
206  CHKERR tetgen_iface->outData(in, out, moab_tetgen_map, tetgen_moab_map,
207  bit_level2, false, true);
208  // set markers to triangles
209  CHKERR tetgen_iface->getTriangleMarkers(tetgen_moab_map, out);
210  // set refions to triangles
211  CHKERR tetgen_iface->getRegionData(tetgen_moab_map, out);
212 
213  // char tetgen_out_file_name[] = "out";
214  // out.save_elements(tetgen_out_file_name);
215 
216  // post-process results, save a mesh and skin
217  EntityHandle meshset_level2;
218  CHKERR moab.create_meshset(MESHSET_SET, meshset_level2);
219  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
220  bit_level2, BitRefLevel().set(), MBTET, meshset_level2);
221  if (debug)
222  CHKERR moab.write_file("level2.vtk", "VTK", "", &meshset_level2, 1);
223  EntityHandle meshset_skin_level2;
224  CHKERR moab.create_meshset(MESHSET_SET, meshset_skin_level2);
225  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
226  bit_level2, BitRefLevel().set(), MBTRI, meshset_skin_level2);
227  if (debug)
228  CHKERR moab.write_file("level_skin2.vtk", "VTK", "", &meshset_skin_level2,
229  1);
230  }
231  CATCH_ERRORS;
232 
233  MoFEM::Core::Finalize();
234 }

Variable Documentation

◆ debug

int debug = 1
static

Definition at line 12 of file tetgen_interface.cpp.

◆ help

char help[] = "...\n\n"
static

Definition at line 11 of file tetgen_interface.cpp.

MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
SIDESET
@ SIDESET
Definition: definitions.h:147
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
EntityHandle
MoFEM::TetGenInterface::getRegionData
MoFEMErrorCode getRegionData(tetGenMoab_Map &tetgen_moab_map, tetgenio &out, Range *ents=NULL, idxRange_Map *ents_map=NULL)
get region data to tetrahedral
Definition: TetGenInterface.cpp:606
MoFEM::TetGenInterface::setFaceData
MoFEMErrorCode setFaceData(std::vector< std::pair< Range, int > > &markers, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map)
set markers to faces
Definition: TetGenInterface.cpp:454
help
static char help[]
Definition: tetgen_interface.cpp:11
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::TetGenInterface::tetRahedralize
MoFEMErrorCode tetRahedralize(char switches[], tetgenio &in, tetgenio &out)
run tetgen
Definition: TetGenInterface.cpp:646
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::TetGenInterface::outData
MoFEMErrorCode outData(tetgenio &in, tetgenio &out, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Range *ents=NULL, bool id_in_tags=false, bool error_if_created=false, bool assume_first_nodes_the_same=false, Tag th=nullptr)
get entities for TetGen data structure
Definition: TetGenInterface.cpp:246
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MoFEM::TetGenInterface
TetGen interface.
Definition: TetGenInterface.hpp:23
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
debug
static int debug
Definition: tetgen_interface.cpp:12
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
_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:1094
MoFEM::TetGenInterface::groupRegion_Triangle
MoFEMErrorCode groupRegion_Triangle(Range &tris, std::vector< std::vector< Range > > &sorted, const double eps=1e-9)
Group surface triangles in planar regions.
Definition: TetGenInterface.cpp:812
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
MoFEM::TetGenInterface::inData
MoFEMErrorCode inData(Range &ents, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Tag th=NULL)
create TetGen data structure form range of moab entities
Definition: TetGenInterface.cpp:38
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::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::TetGenInterface::setRegionData
MoFEMErrorCode setRegionData(std::vector< std::pair< EntityHandle, int > > &regions, tetgenio &in, Tag th=NULL)
set region data to tetrahedral
Definition: TetGenInterface.cpp:556
MoFEM::TetGenInterface::setGeomData
MoFEMErrorCode setGeomData(tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, std::map< int, Range > &type_ents)
set point tags and type
Definition: TetGenInterface.cpp:197
MoFEM::TetGenInterface::getTriangleMarkers
MoFEMErrorCode getTriangleMarkers(tetGenMoab_Map &tetgen_moab_map, tetgenio &out, Range *ents=NULL, idxRange_Map *ents_map=NULL, bool only_non_zero=true)
get markers to faces
Definition: TetGenInterface.cpp:509
MoFEM::TetGenInterface::makePolygonFacet
MoFEMErrorCode makePolygonFacet(Range &ents, Range &polygons, bool reduce_edges=false, Range *not_reducable_nodes=NULL, const double eps=1e-9, Tag th=NULL)
Definition: TetGenInterface.cpp:842
Generated by Doxygen 1.8.17 and hosted at University of Glasgow