Line data Source code
1 : /*
2 : * Copyright (c) 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 <algorithm>
34 : #include "delaunay_triangulation.h"
35 : #include "lib_grid/algorithms/grid_generation/triangle_fill_sweep_line.h"
36 : #include "lib_grid/algorithms/subset_util.h"
37 : #include "lib_grid/algorithms/polychain_util.h"
38 : //temporary
39 : #include "lib_grid/file_io/file_io.h"
40 :
41 : using namespace std;
42 :
43 : namespace ug{
44 :
45 : ////////////////////////////////////////////////////////////////////////////////
46 : class DelaunayDebugSaver
47 : {
48 : public:
49 0 : static DelaunayDebugSaver& inst()
50 : {
51 0 : static DelaunayDebugSaver dds;
52 0 : return dds;
53 : }
54 :
55 : void save(Grid& g, const char* msg)
56 : {
57 : if(m_saveEnabled){
58 : std::stringstream ss;
59 : ss << "dbg" << m_counter++ << ".ugx";
60 :
61 : UG_LOG(msg << ": debug-save to " << ss.str() << std::endl);
62 : SubsetHandler sh(g);
63 : AssignGridToSubset(g, sh, 0);
64 : SaveGridToFile(g, sh, ss.str().c_str());
65 : }
66 : }
67 :
68 : template <class TAAPos>
69 0 : void save(Grid& g, const char* msg, DelaunayInfo<TAAPos>& dinfo)
70 : {
71 0 : if(m_saveEnabled){
72 0 : std::stringstream ss;
73 0 : ss << "dbg" << m_counter++ << ".ugx";
74 :
75 0 : UG_LOG(msg << ": debug-save to " << ss.str() << std::endl);
76 0 : SubsetHandler sh(g);
77 :
78 : for(VertexIterator iter = g.begin<Vertex>();
79 0 : iter != g.end<Vertex>(); ++iter)
80 : {
81 0 : sh.assign_subset(*iter, dinfo.mark(*iter));
82 : }
83 :
84 : for(EdgeIterator iter = g.begin<Edge>();
85 0 : iter != g.end<Edge>(); ++iter)
86 : {
87 0 : sh.assign_subset(*iter, dinfo.mark(*iter));
88 : }
89 :
90 : for(FaceIterator iter = g.begin<Face>();
91 0 : iter != g.end<Face>(); ++iter)
92 : {
93 0 : sh.assign_subset(*iter, dinfo.mark(*iter));
94 : }
95 :
96 0 : sh.subset_info(0).name = "none";
97 0 : sh.subset_info(1).name = "inner";
98 0 : sh.subset_info(2).name = "new-inner";
99 0 : sh.subset_info(3).name = "segment";
100 0 : sh.subset_info(4).name = "new-segment";
101 0 : sh.subset_info(5).name = "dart";
102 0 : sh.subset_info(6).name = "shell";
103 :
104 0 : AssignSubsetColors(sh);
105 0 : SaveGridToFile(g, sh, ss.str().c_str());
106 0 : }
107 0 : }
108 :
109 0 : void enable_save(bool enable = true) {m_saveEnabled = enable;}
110 :
111 : private:
112 0 : DelaunayDebugSaver() : m_saveEnabled(false), m_counter(0) {}
113 :
114 : bool m_saveEnabled;
115 : int m_counter;
116 : };
117 :
118 : static void EnableDelaunayDebugSave(bool enable = true)
119 : {
120 0 : DelaunayDebugSaver::inst().enable_save(enable);
121 : }
122 :
123 : // static void DelaunayDebugSave(Grid& g, const char* msg)
124 : // {
125 : // DelaunayDebugSaver::inst().save(g, msg);
126 : // }
127 :
128 : template <class TAAPos>
129 : static void DelaunayDebugSave(Grid& g, const char* msg, DelaunayInfo<TAAPos>& dinfo)
130 : {
131 0 : DelaunayDebugSaver::inst().save(g, msg, dinfo);
132 0 : }
133 :
134 :
135 :
136 : ////////////////////////////////////////////////////////////////////////////////
137 : template <class TAAPos>
138 0 : bool MakeDelaunay(DelaunayInfo<TAAPos>& info)
139 : {
140 : using namespace std;
141 : typedef typename TAAPos::ValueType vector_t;
142 :
143 : Grid& grid = info.grid();
144 : Face* nbrFaces[2];
145 : vector<Edge*> edges;
146 :
147 : TAAPos& aaPos = info.position_accessor();
148 :
149 0 : while(info.candidates_left()){
150 0 : Edge* e = info.pop_candidate();
151 :
152 : // we only perform swaps on regular manifolds.
153 0 : if(GetAssociatedFaces(nbrFaces, grid, e, 2) == 2){
154 : // make sure that both neighbors are triangles
155 0 : if(nbrFaces[0]->num_vertices() != 3 || nbrFaces[1]->num_vertices() != 3)
156 0 : continue;
157 :
158 0 : Vertex* conVrt0 = GetConnectedVertex(e, nbrFaces[0]);
159 0 : Vertex* conVrt1 = GetConnectedVertex(e, nbrFaces[1]);
160 :
161 0 : vector_t& v0 = aaPos[e->vertex(0)];
162 0 : vector_t& v1 = aaPos[e->vertex(1)];
163 : vector_t& v2 = aaPos[conVrt0];
164 : vector_t& v3 = aaPos[conVrt1];
165 :
166 : vector_t cc1, cc2;
167 :
168 0 : bool cc1_ok = TriangleCircumcenter(cc1, v0, v1, v2);
169 0 : bool cc2_ok = TriangleCircumcenter(cc2, v1, v0, v3);
170 :
171 0 : number r1Sq = VecDistanceSq(cc1, v0);
172 0 : number r2Sq = VecDistanceSq(cc2, v0);
173 :
174 : // for stability reasons, we're checking against the smaller circle
175 0 : if(cc1_ok){
176 0 : if(cc2_ok){
177 0 : if(r1Sq <= r2Sq){
178 0 : if(r1Sq <= VecDistanceSq(cc1, v3))
179 0 : continue; //the edge is fine
180 : }
181 : else{
182 0 : if(r2Sq <= VecDistanceSq(cc2, v2))
183 0 : continue; //the edge is fine
184 : }
185 : }
186 : else{
187 0 : if(r1Sq <= VecDistanceSq(cc1, v3))
188 0 : continue; //the edge is fine
189 : }
190 : }
191 0 : else if(cc2_ok){
192 0 : if(r2Sq <= VecDistanceSq(cc2, v2))
193 0 : continue; //the edge is fine
194 : }
195 : else{
196 : UG_LOG("TriangleCircumcenter failed! Excpect non-delaunay output!\n");
197 : UG_LOG(" This is most likely caused by two degenerated triangles which "
198 : "share an edge.\n");
199 0 : UG_LOG("edge center: " << CalculateCenter(e, aaPos) << endl);
200 0 : return false;
201 : //continue;
202 : }
203 :
204 : // UG_LOG(" swap-it!\n");
205 :
206 : // before swapping, we have to make sure, that the generated triangle
207 : // won't be degenerated.
208 : // this is a costly operation... we check whether both circumcenters
209 : // of the new triangles can be calculated...
210 : /*
211 : if(!(TriangleCircumcenter(cc1, v0, v2, v3)
212 : && TriangleCircumcenter(cc2, v2, v1, v3)))
213 : {
214 : // we have to abort the swap
215 : continue;
216 : }
217 : */
218 : // ok - everything is fine. Now swap the edge
219 0 : Edge* eNew = SwapEdge(grid, e);
220 :
221 0 : if(!eNew){
222 : UG_LOG("An edge-swap failed. Expect degenerated or flipped triangles "
223 : "and a non-delaunay output!\n");
224 0 : UG_LOG("edge center: " << CalculateCenter(e, aaPos) << endl);
225 0 : return false;
226 : //continue;
227 : }
228 :
229 : e = eNew;
230 :
231 : //DelaunayDebugSave(grid, "Edge Swapped", info);
232 :
233 : // all edges of associated triangles are candidates again (except e)
234 0 : GetAssociatedFaces(nbrFaces, grid, e, 2);
235 0 : for(size_t i = 0; i < 2; ++i){
236 0 : CollectAssociated(edges, grid, nbrFaces[i]);
237 0 : for(size_t j = 0; j < edges.size(); ++j){
238 0 : if((edges[j] != e) && !(info.is_segment(edges[j])
239 0 : || info.is_candidate(edges[j])))
240 : {
241 0 : info.push_candidate(edges[j]);
242 : }
243 : }
244 : }
245 : }
246 : }
247 : return true;
248 0 : }
249 :
250 :
251 : static bool
252 0 : DelaunayLineLineIntersection(vector3& vOut,
253 : const vector3& lineFrom, const vector3& lineTo,
254 : const vector3& edgeVrt1, const vector3& edgeVrt2,
255 : vector3 areaNormal,
256 : number smallsq = SMALL_SQ)
257 : {
258 : number t;
259 :
260 : vector3 lineDir, edgeDir, planeNormal;
261 : VecSubtract(lineDir, lineTo, lineFrom);
262 0 : VecNormalize(lineDir, lineDir);
263 : VecSubtract(edgeDir, edgeVrt2, edgeVrt1);
264 0 : number threshold = smallsq * VecLengthSq(edgeDir);
265 0 : VecNormalize(edgeDir, edgeDir);
266 0 : VecCross(planeNormal, edgeDir, areaNormal);
267 0 : if(RayPlaneIntersection(vOut, t, lineFrom, lineDir, edgeVrt1, planeNormal)){
268 0 : if(t > 0 && (VecDistanceSq(lineFrom, vOut) < VecDistanceSq(lineFrom, lineTo) + threshold)){
269 0 : return true;
270 : }
271 : }
272 : return false;
273 : }
274 :
275 :
276 :
277 : template <class TAAPos>
278 0 : bool QualityGridGeneration(Grid& grid, DelaunayInfo<TAAPos>& info,
279 : number minAngle, int maxSteps)
280 : {
281 : EnableDelaunayDebugSave(false);
282 0 : minAngle = max<number>(minAngle, 0);
283 0 : minAngle = min<number>(minAngle, 60);
284 :
285 : // maxSteps = 20;
286 : // UG_LOG("DEBUG: Setting maxSteps to " << maxSteps << "\n");
287 :
288 : using namespace std;
289 :
290 : typedef typename TAAPos::ValueType vector_t;
291 : typedef DelaunayInfo<TAAPos> DI;
292 :
293 : TAAPos aaPos = info.position_accessor();
294 :
295 : int stepCount = 0;
296 :
297 : // helper to collect neighbors
298 : Face* nbrFaces[2];
299 : queue<Vertex*> qvrts; // used during splits
300 : vector<Vertex*> vrts;
301 : vector<Edge*> edges;
302 : vector<Edge*> closeEdges; // used during splits
303 : vector<Face*> faces;
304 :
305 0 : MakeDelaunay(info);
306 :
307 0 : if(minAngle > 0 && maxSteps != 0){
308 0 : const number offCenterThresholdAngle = 1.1 * minAngle;
309 0 : const number offCenterATan = atan(deg_to_rad(0.5 * offCenterThresholdAngle));
310 : // UG_LOG("off-center atan: " << offCenterATan << endl);
311 :
312 0 : info.enable_face_classification(minAngle);
313 :
314 : // while there are faces left which have to be improved
315 0 : while(info.classified_faces_left()){
316 : /*if(stepCount == 488){
317 : EnableDelaunayDebugSave();
318 : }*/
319 :
320 0 : Face* f = info.pop_classified_face();
321 0 : if(f->num_vertices() != 3)
322 0 : continue;
323 :
324 : //todo: the normal is only required for 3d-types. indeed this will crash
325 : // for 2d. This should be moved to a Ray-Line_Intersection3d test.
326 : vector3 faceNormal;
327 0 : CalculateNormal(faceNormal, f, aaPos);
328 : // we can't operate on degenerated faces. Let's hope, that this face
329 : // will improve during improvement of some non-degenerated face.
330 0 : if(VecLengthSq(faceNormal) < SMALL)
331 0 : return false;
332 :
333 0 : vector_t faceCenter = CalculateCenter(f, aaPos);
334 :
335 : // calculate triangle-circumcenter
336 0 : vector_t vpos[3] = {aaPos[f->vertex(0)], aaPos[f->vertex(1)], aaPos[f->vertex(2)]};
337 : vector_t cc;
338 :
339 0 : if(!TriangleCircumcenter(cc, vpos[0], vpos[1], vpos[2])){
340 : UG_LOG("Couldn't calculate triangle-circumcenter. Expect unexpected results!\n");
341 0 : UG_LOG("triangle: " << faceCenter << "\n");
342 : //SaveGridToFile(grid, "delaunay_debug.ugx");
343 : return false;
344 : //continue;
345 : }
346 :
347 : // Test whether an off-center would be more appropriate ('Alper Üngörs Off-Centers')
348 0 : if(offCenterATan > 0)
349 : {
350 0 : number edgeLenSq[3] = {VecDistanceSq(vpos[0], vpos[1]),
351 0 : VecDistanceSq(vpos[1], vpos[2]),
352 0 : VecDistanceSq(vpos[2], vpos[0])};
353 : size_t shortestEdge = 0;
354 0 : for(size_t iedge = 1; iedge < 3; ++iedge){
355 0 : if(edgeLenSq[iedge] < edgeLenSq[shortestEdge])
356 : shortestEdge = iedge;
357 : }
358 :
359 0 : size_t vrtInd[2] = {shortestEdge, (shortestEdge + 1) % 3};
360 0 : vector_t dir[2];
361 0 : for(size_t ivrt = 0; ivrt < 2 ; ++ivrt){
362 0 : VecSubtract(dir[ivrt], vpos[vrtInd[ivrt]], cc);
363 : }
364 :
365 0 : number angle = rad_to_deg(VecAngle(dir[0], dir[1]));
366 :
367 0 : if(angle < offCenterThresholdAngle){
368 : vector_t base;
369 : VecScaleAdd(base, 0.5, vpos[vrtInd[0]], 0.5, vpos[vrtInd[1]]);
370 : vector_t ndir;
371 : VecSubtract(ndir, cc, base);
372 0 : VecNormalize(ndir, ndir);
373 0 : VecScale(ndir, ndir, 0.5 * sqrt(edgeLenSq[shortestEdge]) / offCenterATan);
374 : // UG_LOG("Adjusting cc from: " << cc << " to: ");
375 : VecAdd(cc, base, ndir);
376 : // UG_LOG(cc << endl);
377 : }
378 : }
379 :
380 : // locate the triangle which contains cc. Do this by traversing edges
381 : // as required. Note that since the delaunay property holds, we're only
382 : // traversing edges in a circle, which does not contain any vertices.
383 : Edge* lastTraversedEdge = NULL;
384 : Face* curFace = f;
385 : //vector_t startPos = CalculateCenter(f, aaPos);
386 : vector_t rayDir;
387 : VecSubtract(rayDir, cc, faceCenter);
388 0 : Vertex* pointInserted = NULL;
389 :
390 : // UG_LOG("curFace: " << faceCenter << "\n");
391 0 : while(pointInserted == NULL){
392 : //UG_LOG("curTri: " << CalculateCenter(curFace, aaPos) << "\n");
393 : // to make things as robust as possible, we'll always intersect the
394 : // same line with each edge
395 : Edge* nextEdge = NULL;
396 : bool split = false;
397 :
398 : CollectAssociated(edges, grid, curFace);
399 0 : for(size_t i = 0; i < edges.size(); ++i){
400 0 : Edge* e = edges[i];
401 0 : if(e == lastTraversedEdge)
402 0 : continue;
403 :
404 : vector_t a;
405 0 : if(DelaunayLineLineIntersection(a, faceCenter, cc,
406 0 : aaPos[e->vertex(0)], aaPos[e->vertex(1)], faceNormal, SMALL_SQ))
407 : {
408 : // this edge will be traversed next.
409 : nextEdge = e;
410 0 : number threshold = VecDistanceSq(faceCenter, cc)
411 0 : - SMALL_SQ * EdgeLengthSq(e, aaPos);
412 :
413 : // check whether e has to be split, to avoid bad aspect ratios
414 0 : split = (VecDistanceSq(faceCenter, a) > threshold);
415 0 : break;
416 : }/*
417 : // This else-condition is not necessary - if they are parallel, that's not a problem...
418 : else{
419 : UG_LOG("Ray-Plane intersection failed in step " << stepCount << "... aborting.\n");
420 : UG_LOG(" curTri: " << CalculateCenter(curFace, aaPos)
421 : << ", curEdge: " << CalculateCenter(e, aaPos) << endl);
422 : UG_LOG(" face-normal: " << faceNormal << ", plane-normal: " << planeNormal
423 : << " ray-dir: " << rayDir << endl);
424 : return false;
425 : }*/
426 : }
427 :
428 0 : if(nextEdge){
429 : // check whether the edge has to be splitted
430 : const bool isSegment = info.is_segment(nextEdge);
431 0 : split |= isSegment;
432 0 : if(!split){
433 0 : int numNbrs = GetAssociatedFaces(nbrFaces, grid, nextEdge, 2);
434 0 : if(numNbrs != 2)
435 : split = true;
436 : else{
437 : // get the next face
438 0 : if(nbrFaces[0] == curFace)
439 0 : curFace = nbrFaces[1];
440 : else
441 : curFace = nbrFaces[0];
442 : // and set the last traversed edge
443 : lastTraversedEdge = nextEdge;
444 : }
445 : }
446 :
447 0 : if(split){
448 0 : vector_t center = CalculateCenter(nextEdge, aaPos);
449 :
450 0 : Vertex* vrt0 = nextEdge->vertex(0);
451 0 : Vertex* vrt1 = nextEdge->vertex(1);
452 0 : Vertex* edgeVrts[2] = {vrt0, vrt1};
453 0 : number radiusSq = VecDistanceSq(center, aaPos[vrt0]);
454 0 : number radius = sqrt(radiusSq);
455 0 : pointInserted = SplitEdge<RegularVertex>(grid, nextEdge, false);
456 :
457 0 : if(isSegment){
458 : // depending on the marks of the corners of nextEdge,
459 : // we may have to place the new point at a circular shell
460 : int dartInd = -1;
461 0 : for(int ivrt = 0; ivrt < 2; ++ivrt){
462 0 : if(info.mark(edgeVrts[ivrt]) == DI::DART){
463 : dartInd = ivrt;
464 : break;
465 : }
466 : }
467 0 : if(dartInd != -1){
468 0 : Vertex* dartVrt = edgeVrts[dartInd];
469 0 : Vertex* otherVrt = edgeVrts[(dartInd + 1) % 2];
470 : typename DI::Mark m = info.mark(otherVrt);
471 0 : if((m == DI::NEW_SEGMENT) || (m == DI::SHELL)){
472 : // the new vertex has to be placed on a circular shell!
473 : info.set_mark(pointInserted, DI::SHELL);
474 : number dist = VecDistance(center, aaPos[dartVrt]);
475 : number csCur = 1;
476 : number csNext;
477 :
478 0 : if(dist >= 1){
479 : csNext = 2;
480 0 : while(csNext < dist){
481 : csCur = csNext;
482 0 : csNext *= 2.;
483 : }
484 : }
485 : else{
486 : csNext = 0.5;
487 0 : while(csNext > dist){
488 : csCur = csNext;
489 0 : csNext *= 0.5;
490 : }
491 : }
492 :
493 : vector_t dir;
494 : VecSubtract(dir, center, aaPos[dartVrt]);
495 0 : VecNormalize(dir, dir);
496 :
497 0 : if(fabs(dist - csCur) < fabs(dist - csNext))
498 : VecScale(dir, dir, csCur);
499 : else
500 : VecScale(dir, dir, csNext);
501 :
502 : VecAdd(center, aaPos[dartVrt], dir);
503 :
504 : // we have to adjust the critical radius
505 0 : radiusSq = min(VecDistanceSq(center, aaPos[vrt0]),
506 0 : VecDistanceSq(center, aaPos[vrt1]));
507 0 : radius = sqrt(radiusSq);
508 : }
509 : }
510 :
511 : aaPos[pointInserted] = center;
512 :
513 : // we have to erase all vertices, which are marked and
514 : // which are closer to the inserted point than the
515 : // radius calculated above.
516 : //todo: And which are visible (not hidden by constrained edges)
517 : assert(qvrts.empty());
518 : qvrts.push(pointInserted);
519 :
520 : vrts.clear();
521 : closeEdges.clear();
522 :
523 0 : grid.begin_marking();
524 0 : while(!qvrts.empty()){
525 0 : Vertex* vrt = qvrts.front();
526 : qvrts.pop();
527 :
528 : // collect all edges connected to vrt
529 : CollectAssociated(edges, grid, vrt);
530 : // if the edge is not yet marked, we'll add it to the
531 : // closeEdges list.
532 0 : for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
533 0 : Edge* e = edges[i_edge];
534 0 : if(!grid.is_marked(e)){
535 : // check whether the edge intersects the critical circle
536 0 : if(DistancePointToLine(center, aaPos[e->vertex(0)], aaPos[e->vertex(1)])
537 : < radius)
538 : {
539 : // if the edge is a constrained edge, we'll push it to
540 : // closeEdges
541 0 : if(info.is_segment(e) && !EdgeContains(e, pointInserted))
542 0 : closeEdges.push_back(e);
543 :
544 : grid.mark(e);
545 :
546 : // the connected vertex has to be pushed to
547 : // the queue, regardless whether it lies in the circle
548 : // or not, since an edge from this vertex could
549 : // reenter the circle
550 0 : Vertex* vcon = GetConnectedVertex(e, vrt);
551 0 : if(grid.is_marked(vcon))
552 0 : continue;
553 : grid.mark(vcon);
554 : qvrts.push(vcon);
555 :
556 : // vrt0 and vrt1 won't be removed anyways
557 0 : if(vcon == vrt0 || vcon == vrt1)
558 0 : continue;
559 :
560 : // if the vertrex was created during remeshing
561 : // and lies in the circle, then we'll push it to vrts
562 : if((info.mark(vcon) == DI::NEW_INNER)
563 0 : && (VecDistanceSq(aaPos[vcon], center) < radiusSq))
564 : {
565 0 : vrts.push_back(vcon);
566 : }
567 : }
568 : }
569 : }
570 : }
571 0 : grid.end_marking();
572 :
573 : //UG_LOG(vrts.size() << " vertices have to be deleted!\n");
574 :
575 : // vrts now contains all vertices which lie in a circle of radius
576 : // 'radiusSq' around the center. We now have to check for each,
577 : // whether it is visible from the center, by checking whether
578 : // the line of sight intersects a constrained edge in
579 : // closeEdges.
580 : // if it is visible, then we'll delete it and locally remesh
581 : // the grid.
582 : // We want to record all new edges that are created during remeshing,
583 : // since they have to be considered as candidates for MakeDelaunay.
584 : // However, since some of them will possibly be deleted during
585 : // following erasures, we have to be careful here. This is
586 : // all handled by info.start/stop_candidate_recording
587 0 : info.start_candidate_recording();
588 0 : for(size_t i_vrts = 0; i_vrts < vrts.size(); ++i_vrts){
589 0 : Vertex* vrt = vrts[i_vrts];
590 : vector_t& v = aaPos[vrt];
591 : //UG_LOG("ERASING VRT at " << v << "...\n");
592 : bool intersects = false;
593 : //UG_LOG(" checking for intersections\n");
594 0 : for(size_t i_edge = 0; i_edge < closeEdges.size(); ++i_edge){
595 0 : vector_t& ev0 = aaPos[closeEdges[i_edge]->vertex(0)];
596 0 : vector_t& ev1 = aaPos[closeEdges[i_edge]->vertex(1)];
597 : vector_t a;
598 0 : if(DelaunayLineLineIntersection(a, faceCenter, v, ev0, ev1, faceNormal)){
599 : //if(LineLineIntersection3d(a, b, center, v, ev0, ev1)){
600 : //UG_LOG(" intersection!\n");
601 : intersects = true;
602 0 : break;
603 : }
604 : }
605 :
606 0 : if(intersects)
607 0 : continue;
608 :
609 : //UG_LOG(" gathering surrounding edges\n");
610 : // no intersection detected. erase the vertex and perform
611 : // local retriangulation. Add surrounding edges to
612 : // the 'edges' array
613 : // Store one face, which will be the parent for all new faces
614 : Face* parent = NULL;
615 : edges.clear();
616 : CollectAssociated(faces, grid, vrt);
617 0 : for(size_t i_face = 0; i_face < faces.size(); ++i_face){
618 0 : Face* f = faces[i_face];
619 : parent = f;
620 0 : for(size_t i = 0; i < f->num_edges(); ++i){
621 0 : Edge* e = grid.get_edge(f, i);
622 0 : if(!EdgeContains(e, vrt)){
623 0 : edges.push_back(e);
624 : //UG_LOG("surrounding edge: " << CalculateCenter(e, aaPos) << endl);
625 : }
626 : }
627 : }
628 : assert(parent);
629 :
630 : //UG_LOG(" retriangulation\n");
631 : // perform retriangulation
632 : // get the vertices of the poly-chain
633 : //UG_LOG("1");
634 : std::vector<Vertex*> vrtPolyChain;
635 0 : CreatePolyChain(vrtPolyChain, grid, edges.begin(), edges.end());
636 : //UG_LOG("2");
637 0 : std::vector<vector_t> posPolyChain(vrtPolyChain.size());
638 : std::vector<int> edgeChain, faceIndices;
639 0 : edgeChain.reserve(vrtPolyChain.size() * 2);
640 : //UG_LOG("3");
641 : size_t numVrts = vrtPolyChain.size();
642 0 : for(size_t i = 0; i < numVrts; ++i){
643 : //UG_LOG(aaPos[vrtPolyChain[i]]);
644 0 : posPolyChain[i] = aaPos[vrtPolyChain[i]];
645 0 : edgeChain.push_back(i);
646 0 : edgeChain.push_back((i+1) % numVrts);
647 : }
648 : /*
649 : ("edges.size: " << edges.size() << ", posPolyChain.size: " << posPolyChain.size());
650 : for(size_t i = 0; i < edgeChain.size(); i+=2){
651 : UG_LOG(" (" << edgeChain[i] << ", " << edgeChain[i+1] << ")");
652 : }
653 : UG_LOG(endl);
654 : */
655 :
656 : //UG_LOG("4");
657 :
658 0 : if(!TriangleFill_SweepLine(faceIndices, posPolyChain, edgeChain)){
659 0 : UG_LOG("Sweepline failed in step " << stepCount << "... aborting\n");
660 0 : UG_LOG(" While examining face " << faceCenter << endl);
661 : return false;
662 : }
663 : //UG_LOG("5");
664 : // create the faces
665 0 : for(size_t i = 0; i < faceIndices.size(); i+=3){
666 0 : grid.create<Triangle>(
667 0 : TriangleDescriptor(vrtPolyChain[faceIndices[i]],
668 : vrtPolyChain[faceIndices[i+1]],
669 : vrtPolyChain[faceIndices[i+2]]),
670 : parent);
671 : }
672 : //UG_LOG("6\n");
673 :
674 : //UG_LOG(" erase\n");
675 : // now erase vrt
676 0 : grid.erase(vrt);
677 :
678 : //DelaunayDebugSave(grid, "After TriangleFill", info);
679 : }
680 :
681 0 : info.stop_candidate_recording();
682 : }
683 : else{
684 : aaPos[pointInserted] = cc;
685 : }
686 : }
687 : }
688 : else{
689 : // we found the triangle, which contains cc. Insert the point
690 : // and perform local delaunay (not necessarily local...).
691 : //...
692 : // But first we have to check whether an edge of the triangle
693 : // would be encroached.
694 : // If this is the case, we'll insert the new vertex at the triangle's
695 : // center instead of the circum-center, to avoid skinny triangles
696 : CollectAssociated(edges, grid, curFace);
697 0 : for(size_t i = 0; i < edges.size(); ++i){
698 0 : Edge* e = edges[i];
699 0 : if(info.is_segment(e)){
700 0 : vector_t eCenter = CalculateCenter(e, aaPos);
701 0 : number eLenSq = EdgeLengthSq(e, aaPos);
702 0 : if(VecDistanceSq(eCenter, cc) < eLenSq / 4.){
703 0 : cc = CalculateCenter(curFace, aaPos);
704 0 : break;
705 : }
706 : }
707 : }
708 :
709 0 : Vertex* vrt0 = curFace->vertex(0);
710 0 : Vertex* vrt1 = curFace->vertex(1);
711 0 : Vertex* vrt2 = curFace->vertex(2);
712 :
713 : //UG_LOG("creating new elements...\n");
714 0 : Vertex* vrt = *grid.create<RegularVertex>(curFace);
715 : aaPos[vrt] = cc;
716 0 : grid.create<Triangle>(TriangleDescriptor(vrt0, vrt1, vrt), curFace);
717 0 : grid.create<Triangle>(TriangleDescriptor(vrt1, vrt2, vrt), curFace);
718 0 : grid.create<Triangle>(TriangleDescriptor(vrt2, vrt0, vrt), curFace);
719 : //UG_LOG("erasing old face\n");
720 0 : grid.erase(curFace);
721 :
722 : //DelaunayDebugSave(grid, "After InsertPoint", info);
723 :
724 0 : pointInserted = vrt;
725 : }
726 :
727 0 : if(pointInserted){
728 : //UG_LOG("temp-save to delaunay_debug.ugx\n");
729 : //SaveGridToFile(grid, "delaunay_debug.ugx");
730 :
731 : //UG_LOG("inserted point\n");
732 : // if a vertex was inserted, we'll have to perform a delaunay step.
733 : // Find new candidates by examining edges of associated triangles of vrt.
734 : // If an edge of such a triangle is connected to exactly 2
735 : // triangles, both marked, then it is a new candidate.
736 : CollectAssociated(faces, grid, pointInserted);
737 0 : for(size_t i_face = 0; i_face < faces.size(); ++i_face){
738 0 : if(!info.is_inner(faces[i_face]))
739 0 : continue;
740 :
741 : CollectAssociated(edges, grid, faces[i_face]);
742 0 : for(size_t i_edge = 0; i_edge < edges.size(); ++i_edge){
743 0 : Edge* e = edges[i_edge];
744 0 : if(info.is_candidate(e) || info.is_segment(e))
745 0 : continue;
746 :
747 : Face* nbrs[2];
748 0 : if(GetAssociatedFaces(nbrs, grid, e, 2) == 2){
749 0 : if(info.is_inner(nbrs[0])
750 0 : && info.is_inner(nbrs[1]))
751 : {
752 0 : info.push_candidate(e);
753 : }
754 : }
755 : }
756 : }
757 :
758 : // DelaunayDebugSave(grid, "Candidates Adjusted After Insert Point", info);
759 :
760 : //UG_LOG("redelaunaylizing\n");
761 : // UG_LOG("MakeDelaunay after InsertPoint\n");
762 0 : if(!MakeDelaunay(info)){
763 0 : UG_LOG("Make Delaunay failed in step " << stepCount << ".\n");
764 0 : UG_LOG(" While examining face " << faceCenter << endl);
765 : return false;
766 : }
767 : // UG_LOG("MakeDelaunay done\n");
768 : // UG_LOG("adaption step completed!\n");
769 : }
770 : }
771 :
772 : /*
773 : // check whether an illegal triangle has been inserted
774 : for(TriangleIterator iter = grid.begin<Triangle>();
775 : iter != grid.end<Triangle>(); ++iter)
776 : {
777 : vector3 n;
778 : CalculateNormal(n, *iter, aaPos);
779 : if(n.z() <= 0){
780 : UG_LOG("ATTENTION: Illegal triangle created in step " << stepCount << "\n");
781 : return false;
782 : }
783 : }
784 : */
785 0 : ++stepCount;
786 : DelaunayDebugSave(grid, "", info);
787 0 : if(stepCount == maxSteps)
788 : break;
789 : /*
790 : if((stepCount % 1000) == 0){
791 : UG_LOG("temp-save to delaunay_debug.ugx in step " << stepCount << endl);
792 : //DelaunayDebugSave(grid, "", info);
793 : SaveGridToFile(grid, "delaunay_debug.ugx");
794 : }
795 : */
796 : }
797 : }
798 : return true;
799 0 : }
800 :
801 :
802 :
803 : template bool MakeDelaunay(DelaunayInfo<Grid::VertexAttachmentAccessor<AVector3> >&);
804 :
805 : template bool QualityGridGeneration(Grid&,
806 : DelaunayInfo<Grid::VertexAttachmentAccessor<AVector3> >&,
807 : number, int);
808 :
809 : }// end of namespace
|
