Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include <vector>
34 : #include <sstream>
35 : #include "tetrahedralization.h"
36 : #include "../geom_obj_util/geom_obj_util.h"
37 : #include "../remove_duplicates_util.h"
38 :
39 : #ifdef UG_TETGEN
40 : #define TETLIBRARY
41 : #include "tetgen.h"
42 : #endif
43 :
44 : static char const* sTetgenMissing = "Tetrahedralization failed since Tetgen is not available in the current build.\n"
45 : "In order to compile with Tetgen, please install Tetgen through `ughub install tetgen`\n"
46 : "and enable Tetgen in your build using `cmake -Dtetgen=ON -DLINK_TETGEN=ON .`";
47 : using namespace std;
48 :
49 : namespace ug
50 : {
51 :
52 : #ifdef UG_TETGEN
53 : const char* VerbosityToTetgenParam(int verbosity)
54 : {
55 : if (verbosity <= 0)
56 : return "";
57 : else if(verbosity == 1)
58 : return "V";
59 : else if(verbosity == 2)
60 : return "VV";
61 : else
62 : return "VVV";
63 : }
64 : #endif
65 :
66 0 : static bool PerformTetrahedralization(Grid& grid,
67 : ISubsetHandler* pSH,
68 : number quality,
69 : bool preserveBnds,
70 : bool preserveAll,
71 : APosition& aPos,
72 : int verbosity)
73 : {
74 : #ifdef UG_TETGEN
75 : if(!grid.has_vertex_attachment(aPos))
76 : return false;
77 :
78 : // access position data
79 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPos);
80 : /*
81 : // make sure that no double-points exist
82 : RemoveDoubles(grid, grid.vertices_begin(), grid.vertices_end(), aPos, 1e-12);
83 : */
84 : // convert quadrilaterals to triangles
85 : if(grid.num<Quadrilateral>() > 0)
86 : Triangulate(grid, grid.begin<Quadrilateral>(), grid.end<Quadrilateral>(), &aaPos);
87 :
88 : if(grid.num_edges() > 0 || grid.num_faces() > 0){
89 : // we have to make sure that all vertices in the grid are actually connected
90 : // to a triangle. If this wouldn't be the case, tetgen would likely have problems.
91 : size_t numVrtsRemoved = 0;
92 : for(VertexIterator iter = grid.begin<Vertex>();
93 : iter != grid.end<Vertex>();)
94 : {
95 : Vertex* v = *iter;
96 : ++iter;
97 : if((NumAssociatedEdges(grid, v) == 0) || (NumAssociatedFaces(grid, v) == 0)){
98 : grid.erase(v);
99 : ++numVrtsRemoved;
100 : }
101 : }
102 :
103 : if(numVrtsRemoved > 0){
104 : UG_LOG("WARNING in Tetrahedralize: Removed " << numVrtsRemoved
105 : << " vertices which were not connected to any edge or face\n");
106 : }
107 : }
108 :
109 : // we have to make sure, that no face-duplicates exist in the grid
110 : RemoveDuplicates(grid, grid.begin<Face>(), grid.end<Face>());
111 :
112 : // attach an index to the vertices
113 : AInt aInd;
114 : grid.attach_to_vertices(aInd);
115 : Grid::VertexAttachmentAccessor<AInt> aaInd(grid, aInd);
116 :
117 : // datastructures to communicate with tetgenio
118 : tetgenio in, out;
119 :
120 : // setup points
121 : {
122 : in.numberofpoints = grid.num_vertices();
123 : in.pointlist = new REAL[in.numberofpoints*3];
124 :
125 : // copy position data
126 : int counter = 0;
127 :
128 : for(VertexIterator iter = grid.vertices_begin();
129 : iter != grid.vertices_end(); ++iter, ++counter)
130 : {
131 : aaInd[*iter] = counter;
132 : vector3& v = aaPos[*iter];
133 : // position types are casted to float, since this circumvents an
134 : // error that occurs on some geometries. Somehow tetgen constructs
135 : // selfintersecting facets otherwise (sometimes). I didn't really understand
136 : // this behaviour yet.
137 : //TODO: Think about how the following code could be improved.
138 :
139 : in.pointlist[counter * 3] = (float)v.x();
140 : in.pointlist[counter * 3 + 1] = (float)v.y();
141 : in.pointlist[counter * 3 + 2] = (float)v.z();
142 : }
143 : }
144 :
145 : // set up facets
146 : {/*
147 : if(pSH){
148 : // collect all edges that are assigned to subsets
149 : vector<Edge*> edges;
150 : for(int si = 0; si < pSH->num_subsets(); ++si){
151 : for(EdgeIterator iter = pSH->begin<Edge>(si);
152 : iter != pSH->end<Edge>(si); ++iter)
153 : {
154 : edges.push_back(*iter);
155 : }
156 : }
157 :
158 : in.numberofedges = (int)edges.size();
159 : in.edgelist = new int[in.numberofedges*2];
160 : in.edgemarkerlist = new int[in.numberofedges];
161 :
162 : for(size_t i = 0; i < edges.size(); ++i){
163 : in.edgelist[i*2] = aaInd[edges[i]->vertex(0)];
164 : in.edgelist[i*2 + 1] = aaInd[edges[i]->vertex(1)];
165 : in.edgemarkerlist[i] = pSH->get_subset_index(edges[i]);
166 : }
167 : UG_LOG("number of edges in: " << in.numberofedges << endl);
168 : }*/
169 :
170 : if(grid.num_faces() > 0){
171 : in.numberoffacets = grid.num_faces();
172 : in.facetlist = new tetgenio::facet[in.numberoffacets];
173 : in.facetmarkerlist = new int[in.numberoffacets];
174 :
175 : int counter = 0;
176 : for(FaceIterator iter = grid.faces_begin();
177 : iter != grid.faces_end(); ++iter, ++counter)
178 : {
179 : Face* f = *iter;
180 :
181 : // copy the face to the facetlist
182 : tetgenio::facet* tf = &in.facetlist[counter];
183 : tf->numberofpolygons = 1;
184 : tf->polygonlist = new tetgenio::polygon[tf->numberofpolygons];
185 : tf->numberofholes = 0;
186 : tf->holelist = NULL;
187 : tetgenio::polygon* p = &tf->polygonlist[0];
188 : p->numberofvertices = f->num_vertices();
189 : p->vertexlist = new int[p->numberofvertices];
190 : for(int i = 0; i < p->numberofvertices; ++i)
191 : p->vertexlist[i] = aaInd[f->vertex(i)];
192 :
193 : // set the face mark
194 : if(pSH)
195 : in.facetmarkerlist[counter] = pSH->get_subset_index(f);
196 : else
197 : in.facetmarkerlist[counter] = 0;
198 : }
199 : }
200 : }
201 :
202 : // the aInd attachment is no longer required
203 : grid.detach_from_vertices(aInd);
204 : aaInd.invalidate();
205 :
206 : // call tetrahedralization
207 : try{
208 : stringstream ss;
209 : ss << VerbosityToTetgenParam(verbosity);
210 : if(grid.num_faces() > 0){
211 : ss << "p";
212 : if(quality > SMALL)
213 : ss << "qq" << quality;
214 : if(preserveBnds || preserveAll)
215 : ss << "Y";
216 : if(preserveAll)
217 : ss << "Y"; // if inner bnds shall be preserved "YY" has to be passed to tetgen
218 : }
219 :
220 : ss << "Q";//"Q";
221 : tetrahedralize(const_cast<char*>(ss.str().c_str()), &in, &out);
222 : }
223 : catch(int errCode){
224 : UG_LOG(" aborting tetrahedralization. Received error: " << errCode << endl);
225 : return false;
226 : }
227 : /*
228 : if(out.numberofpoints < (int)grid.num_vertices()){
229 : LOG(" aborting tetrahedralization - bad number of points\n");
230 : return false;
231 : }
232 : */
233 : // add new vertices to the grid. store all vertices in a vector.
234 : vector<Vertex*> vVrts(out.numberofpoints);
235 : {
236 : int counter = 0;
237 : // add the old ones to the vector
238 : for(VertexIterator iter = grid.vertices_begin();
239 : iter != grid.vertices_end(); ++iter, ++counter)
240 : {
241 : aaPos[*iter].x() = out.pointlist[counter*3];
242 : aaPos[*iter].y() = out.pointlist[counter*3+1];
243 : aaPos[*iter].z() = out.pointlist[counter*3+2];
244 : vVrts[counter] = *iter;
245 : if(counter == out.numberofpoints -1){
246 : UG_LOG("WARNING: Unused points may remain!\n");
247 : break;
248 : }
249 : }
250 : // create new ones and add them to the vector
251 : for(; counter < out.numberofpoints; ++counter)
252 : {
253 : RegularVertex* v = *grid.create<RegularVertex>();
254 : aaPos[v].x() = out.pointlist[counter*3];
255 : aaPos[v].y() = out.pointlist[counter*3+1];
256 : aaPos[v].z() = out.pointlist[counter*3+2];
257 : vVrts[counter] = v;
258 : }
259 : }
260 :
261 : // erase edges if boundary segments were not preserved
262 : if(!preserveAll){
263 : grid.erase(grid.edges_begin(), grid.edges_end());
264 : }
265 :
266 : // add new faces
267 : grid.erase(grid.faces_begin(), grid.faces_end());
268 : for(int i = 0; i < out.numberoftrifaces; ++i)
269 : {
270 : Triangle* tri = *grid.create<Triangle>(TriangleDescriptor(vVrts[out.trifacelist[i*3]],
271 : vVrts[out.trifacelist[i*3 + 1]],
272 : vVrts[out.trifacelist[i*3 + 2]]));
273 :
274 : if(pSH && out.trifacemarkerlist)
275 : pSH->assign_subset(tri, out.trifacemarkerlist[i]);
276 : }
277 :
278 : if(out.numberoftetrahedra < 1)
279 : return false;
280 :
281 : // add new volumes
282 : for(int i = 0; i < out.numberoftetrahedra; ++i)
283 : {
284 : Tetrahedron* tet = *grid.create<Tetrahedron>(
285 : TetrahedronDescriptor(vVrts[out.tetrahedronlist[i*4]],
286 : vVrts[out.tetrahedronlist[i*4 + 1]],
287 : vVrts[out.tetrahedronlist[i*4 + 2]],
288 : vVrts[out.tetrahedronlist[i*4 + 3]]));
289 : if(pSH)
290 : pSH->assign_subset(tet, 0);
291 : }
292 :
293 : return true;
294 :
295 : #else
296 0 : UG_THROW("\nPerformTetrahedralization: " << sTetgenMissing);
297 : #endif
298 :
299 : }
300 :
301 :
302 0 : static bool PerformRetetrahedralization(Grid& grid,
303 : SubsetHandler& sh,
304 : number quality,
305 : bool preserveBnds,
306 : bool preserveAll,
307 : APosition& aPos,
308 : ANumber& aVolCon,
309 : bool applyVolumeConstraint,
310 : int verbosity)
311 : {
312 : #ifdef UG_TETGEN
313 : if(!grid.has_vertex_attachment(aPos))
314 : return false;
315 :
316 : if(!grid.has_volume_attachment(aVolCon))
317 : return false;
318 :
319 : if(grid.num<Tetrahedron>() == 0)
320 : return false;
321 :
322 : // access data
323 : Grid::VertexAttachmentAccessor<APosition> aaPos(grid, aPos);
324 : Grid::VolumeAttachmentAccessor<ANumber> aaVolCon(grid, aVolCon);
325 :
326 : // attach an index to the vertices
327 : AInt aInd;
328 : grid.attach_to_vertices(aInd);
329 : Grid::VertexAttachmentAccessor<AInt> aaInd(grid, aInd);
330 :
331 : // datastructures to communicate with tetgenio
332 : tetgenio in, out;
333 :
334 : // setup points
335 : {
336 : in.numberofpoints = grid.num_vertices();
337 : in.pointlist = new REAL[in.numberofpoints*3];
338 :
339 : // copy position data
340 : int counter = 0;
341 :
342 : for(VertexIterator iter = grid.vertices_begin();
343 : iter != grid.vertices_end(); ++iter, ++counter)
344 : {
345 : aaInd[*iter] = counter;
346 : vector3& v = aaPos[*iter];
347 : // position types are casted to float, since this circumvents an
348 : // error that occurs on some geometries. Somehow tetgen constructs
349 : // selfintersecting facets otherwise (sometimes). I didn't really understand
350 : // this behaviour yet.
351 : //TODO: Think about how the following code could be improved.
352 :
353 : in.pointlist[counter * 3] = (float)v.x();
354 : in.pointlist[counter * 3 + 1] = (float)v.y();
355 : in.pointlist[counter * 3 + 2] = (float)v.z();
356 : //in.pointlist[counter * 3] = v.x();
357 : //in.pointlist[counter * 3 + 1] = v.y();
358 : //in.pointlist[counter * 3 + 2] = v.z();
359 : }
360 : }
361 :
362 : // set up triangles
363 : {
364 : // count the total number of triangles in subsets
365 : int numTris = 0;
366 : for(int i = 0; i < sh.num_subsets(); ++i)
367 : numTris += sh.num<Triangle>(i);
368 :
369 : in.numberoftrifaces = numTris;
370 : in.trifacelist = new int[in.numberoftrifaces * 3];
371 : in.trifacemarkerlist = new int[in.numberoftrifaces];
372 :
373 : int curTri = 0;
374 : for(TriangleIterator iter = grid.begin<Triangle>();
375 : iter != grid.end<Triangle>(); ++iter)
376 : {
377 : Triangle* t = *iter;
378 :
379 : // ignore faces in subset -1
380 : if(sh.get_subset_index(t) == -1)
381 : continue;
382 :
383 : // add the triangle
384 : in.trifacelist[curTri * 3] = aaInd[t->vertex(0)];
385 : in.trifacelist[curTri * 3 + 1] = aaInd[t->vertex(1)];
386 : in.trifacelist[curTri * 3 + 2] = aaInd[t->vertex(2)];
387 :
388 : // set the face mark
389 : in.trifacemarkerlist[curTri] = sh.get_subset_index(t);
390 :
391 : ++curTri;
392 : }
393 : }
394 :
395 : // now fill the tetrahedron lists
396 : {
397 : in.numberoftetrahedra = grid.num<Tetrahedron>();
398 : in.tetrahedronlist = new int[in.numberoftetrahedra * 4];
399 : in.tetrahedronvolumelist = new REAL[in.numberoftetrahedra];
400 : in.numberoftetrahedronattributes = 1;
401 : in.tetrahedronattributelist = new REAL[in.numberoftetrahedra];
402 :
403 : // fill the lists
404 : int curTet = 0;
405 : for(Grid::traits<Tetrahedron>::iterator iter = grid.begin<Tetrahedron>();
406 : iter != grid.end<Tetrahedron>(); ++iter, ++curTet)
407 : {
408 : Tetrahedron* t = *iter;
409 : in.tetrahedronlist[curTet * 4] = aaInd[t->vertex(0)];
410 : in.tetrahedronlist[curTet * 4 + 1] = aaInd[t->vertex(1)];
411 : in.tetrahedronlist[curTet * 4 + 2] = aaInd[t->vertex(2)];
412 : in.tetrahedronlist[curTet * 4 + 3] = aaInd[t->vertex(3)];
413 :
414 : in.tetrahedronvolumelist[curTet] = aaVolCon[t];
415 : in.tetrahedronattributelist[curTet] = (REAL)sh.get_subset_index(t);
416 : }
417 : }
418 :
419 : // the aInd attachment is no longer required
420 : grid.detach_from_vertices(aInd);
421 : aaInd.invalidate();
422 :
423 : // call tetrahedralization
424 : try{
425 : stringstream ss;
426 : ss << VerbosityToTetgenParam(verbosity);
427 : if(applyVolumeConstraint)
428 : ss << "a";
429 : ss << "r";
430 : if(quality > SMALL)
431 : ss << "qq" << quality;
432 : if(preserveBnds || preserveAll)
433 : ss << "Y";
434 : if(preserveAll)
435 : ss << "Y"; // if inner bnds shall be preserved "YY" has to be passed to tetgen
436 :
437 : ss << "Q";
438 : tetrahedralize(const_cast<char*>(ss.str().c_str()), &in, &out);
439 : }
440 : catch(int errCode){
441 : UG_LOG(" aborting retetrahedralization. Received error: " << errCode << endl);
442 : return false;
443 : }
444 :
445 : // first we'll erase all existing tetrahedrons
446 : grid.erase(grid.begin<Tetrahedron>(), grid.end<Tetrahedron>());
447 :
448 : // add new vertices to the grid. store all vertices in a vector.
449 : vector<Vertex*> vVrts(out.numberofpoints);
450 : {
451 : int counter = 0;
452 : // add the old ones to the vector
453 : for(VertexIterator iter = grid.vertices_begin();
454 : iter != grid.vertices_end(); ++iter, ++counter)
455 : {
456 : aaPos[*iter].x() = out.pointlist[counter*3];
457 : aaPos[*iter].y() = out.pointlist[counter*3+1];
458 : aaPos[*iter].z() = out.pointlist[counter*3+2];
459 : vVrts[counter] = *iter;
460 : if(counter == out.numberofpoints -1){
461 : UG_LOG("WARNING: Unused points may remain!\n");
462 : break;
463 : }
464 : }
465 : // create new ones and add them to the vector
466 : for(; counter < out.numberofpoints; ++counter)
467 : {
468 : RegularVertex* v = *grid.create<RegularVertex>();
469 : aaPos[v].x() = out.pointlist[counter*3];
470 : aaPos[v].y() = out.pointlist[counter*3+1];
471 : aaPos[v].z() = out.pointlist[counter*3+2];
472 : vVrts[counter] = v;
473 : }
474 : }
475 :
476 : // erase edges if boundary segments were not preserved
477 : if(!preserveAll){
478 : grid.erase(grid.edges_begin(), grid.edges_end());
479 : }
480 :
481 : // add new faces
482 : grid.erase(grid.faces_begin(), grid.faces_end());
483 : for(int i = 0; i < out.numberoftrifaces; ++i)
484 : {
485 : Triangle* tri = *grid.create<Triangle>(TriangleDescriptor(vVrts[out.trifacelist[i*3]],
486 : vVrts[out.trifacelist[i*3 + 1]],
487 : vVrts[out.trifacelist[i*3 + 2]]));
488 :
489 : sh.assign_subset(tri, out.trifacemarkerlist[i]);
490 : }
491 :
492 : if(out.numberoftetrahedra < 1)
493 : return false;
494 :
495 : // add new volumes
496 : for(int i = 0; i < out.numberoftetrahedra; ++i)
497 : {
498 : Tetrahedron* tet = *grid.create<Tetrahedron>(
499 : TetrahedronDescriptor(vVrts[out.tetrahedronlist[i*4]],
500 : vVrts[out.tetrahedronlist[i*4 + 1]],
501 : vVrts[out.tetrahedronlist[i*4 + 2]],
502 : vVrts[out.tetrahedronlist[i*4 + 3]]));
503 :
504 : sh.assign_subset(tet, (int)out.tetrahedronattributelist[i]);
505 : }
506 :
507 : return true;
508 :
509 : #else
510 0 : UG_THROW("\nPerformRetetrahedralization: " << sTetgenMissing);
511 : return false;
512 : #endif
513 :
514 : }
515 :
516 :
517 0 : bool Tetrahedralize(Grid& grid, number quality, bool preserveBnds,
518 : bool preserveAll, APosition& aPos, int verbosity)
519 : {
520 0 : return PerformTetrahedralization(grid, NULL, quality, preserveBnds,
521 : preserveAll, aPos, verbosity);
522 : }
523 :
524 0 : bool Tetrahedralize(Grid& grid, ISubsetHandler& sh, number quality,
525 : bool preserveBnds, bool preserveAll, APosition& aPos,
526 : int verbosity)
527 : {
528 0 : return PerformTetrahedralization(grid, &sh, quality, preserveBnds,
529 : preserveAll, aPos, verbosity);
530 : }
531 :
532 0 : bool Retetrahedralize(Grid& grid, SubsetHandler& sh,
533 : ANumber& aVolumeConstraint,
534 : number quality,
535 : bool preserveBnds,
536 : bool preserveAll,
537 : APosition& aPos,
538 : bool applyVolumeConstraint,
539 : int verbosity)
540 : {
541 0 : return PerformRetetrahedralization(grid, sh, quality, preserveBnds,
542 : preserveAll, aPos, aVolumeConstraint,
543 : applyVolumeConstraint,
544 : verbosity);
545 : }
546 :
547 : }// end of namespace
|
