v0.14.0
Loading...
Searching...
No Matches
tetgen_interface.cpp
Go to the documentation of this file.
1/** \file tetgen_interface.cpp
2 * \brief test for TetGen interface
3 *
4 * \ingroup mesh_tetgen
5 */
6
7#include <MoFEM.hpp>
8
9using namespace MoFEM;
10
11static char help[] = "...\n\n";
12static int debug = 1;
13
14int main(int argc, char *argv[]) {
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;
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);
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
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
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 }
232
234}
int main()
#define CATCH_ERRORS
Catch errors.
@ SIDESET
@ BLOCKSET
#define CHKERR
Inline error check.
#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
char mesh_file_name[255]
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
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
Deprecated interface functions.
MoFEMErrorCode makePolygonFacet(Range &ents, Range &polygons, bool reduce_edges=false, Range *not_reducable_nodes=NULL, const double eps=1e-9, Tag th=NULL)
MoFEMErrorCode groupRegion_Triangle(Range &tris, std::vector< std::vector< Range > > &sorted, const double eps=1e-9)
Group surface triangles in planar regions.
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
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
MoFEMErrorCode tetRahedralize(char switches[], tetgenio &in, tetgenio &out)
run tetgen
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
MoFEMErrorCode getRegionData(tetGenMoab_Map &tetgen_moab_map, tetgenio &out, Range *ents=NULL, idxRange_Map *ents_map=NULL)
get region data to tetrahedral
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
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
MoFEMErrorCode setRegionData(std::vector< std::pair< EntityHandle, int > > &regions, tetgenio &in, Tag th=NULL)
set region data to tetrahedral
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
static char help[]
static int debug