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 "delaunay_info.h"
34 :
35 : namespace ug{
36 :
37 : template <class TAAPos>
38 0 : DelaunayInfo<TAAPos>::
39 : DelaunayInfo(Grid& g, TAAPos& aaPos,
40 : Grid::edge_traits::callback cbConstrainedEdge)
41 0 : : m_grid(g), m_aaPos(aaPos), m_minAngle(0),
42 0 : m_maxDot(1), m_cbConstrainedEdge(cbConstrainedEdge),
43 0 : m_candidateRecordingEnabled(false)
44 : {
45 0 : g.attach_to_vertices_dv(m_aMark, NONE);
46 0 : g.attach_to_edges_dv(m_aMark, NONE);
47 0 : g.attach_to_faces_dv(m_aMark, NONE);
48 0 : m_aaMark.access(g, m_aMark, true, true, true, false);
49 :
50 0 : g.attach_to_edges_dv(m_aCandidate, false);
51 : m_aaCandidateEDGE.access(g, m_aCandidate);
52 :
53 0 : if(!g.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES)){
54 : UG_LOG("WARNING in DelaunayInfo: enabling FACEOPT_AUTOGENERATE_EDGES.\n");
55 0 : g.enable_options(FACEOPT_AUTOGENERATE_EDGES);
56 : }
57 :
58 0 : g.register_observer(this, OT_VERTEX_OBSERVER | OT_EDGE_OBSERVER | OT_FACE_OBSERVER);
59 0 : }
60 :
61 : template <class TAAPos>
62 0 : DelaunayInfo<TAAPos>::
63 : ~DelaunayInfo()
64 : {
65 0 : m_grid.unregister_observer(this);
66 0 : enable_face_classification(0);
67 0 : m_grid.detach_from_vertices(m_aMark);
68 0 : m_grid.detach_from_edges(m_aMark);
69 0 : m_grid.detach_from_faces(m_aMark);
70 0 : m_grid.detach_from_edges(m_aCandidate);
71 0 : }
72 :
73 :
74 : template <class TAAPos>
75 0 : void DelaunayInfo<TAAPos>::
76 : set_mark(Face* f, typename DelaunayInfo<TAAPos>::Mark mark)
77 : {
78 0 : m_aaMark[f] = mark;
79 0 : if(face_classification_enabled()){
80 0 : FaceInfo* fi = m_aaFaceInfo[f];
81 0 : if((mark == INNER) || (mark == NEW_INNER)){
82 0 : if(!fi){
83 0 : m_aaFaceInfo[f] = fi = new FaceInfo();
84 0 : fi->f = f;
85 : }
86 :
87 0 : if(!fi->classified)
88 0 : classify_face(f);
89 : }
90 0 : else if(fi){
91 0 : if(fi->classified){
92 : // we can't delete the face-info now. instead we'll only set fi->f to NULL
93 : // so that it can be identified as invalid in m_faceQueue. Clean-up is performed later.
94 0 : fi->f = NULL;
95 : }
96 : else{
97 : // delete the associated face-info
98 0 : delete fi;
99 0 : m_aaFaceInfo[f] = NULL;
100 : }
101 : }
102 : }
103 0 : }
104 :
105 :
106 : template <class TAAPos>
107 0 : bool DelaunayInfo<TAAPos>::
108 : is_dart_segment(Edge* e)
109 : {
110 0 : Edge::ConstVertexArray vrts = e->vertices();
111 0 : return is_segment(e) && (mark(vrts[0]) == DART || mark(vrts[1]) == DART);
112 : }
113 :
114 : template <class TAAPos>
115 0 : bool DelaunayInfo<TAAPos>::
116 : is_new_dart_segment(Edge* e)
117 : {
118 0 : Edge::ConstVertexArray vrts = e->vertices();
119 0 : return mark(e) == NEW_SEGMENT && (mark(vrts[0]) == DART || mark(vrts[1]) == DART);
120 : }
121 :
122 :
123 : template <class TAAPos>
124 0 : bool DelaunayInfo<TAAPos>::
125 : is_dart_shell_segment(Edge* e)
126 : {
127 0 : Edge::ConstVertexArray vrts = e->vertices();
128 : return is_segment(e)
129 0 : && ((mark(vrts[0]) == DART && mark(vrts[1]) == SHELL)
130 0 : || (mark(vrts[0]) == SHELL && mark(vrts[1]) == DART));
131 : }
132 :
133 :
134 : template <class TAAPos>
135 0 : void DelaunayInfo<TAAPos>::
136 : push_candidate(Edge* e)
137 : {
138 0 : if(!is_candidate(e)){
139 : m_aaCandidateEDGE[e] = true;
140 : m_qEdgeCandidates.push(e);
141 : }
142 0 : }
143 :
144 :
145 : template <class TAAPos>
146 0 : Edge* DelaunayInfo<TAAPos>::
147 : pop_candidate()
148 : {
149 0 : Edge* e = m_qEdgeCandidates.front();
150 : m_qEdgeCandidates.pop();
151 : m_aaCandidateEDGE[e] = false;
152 0 : return e;
153 : }
154 :
155 : template <class TAAPos>
156 0 : void DelaunayInfo<TAAPos>::
157 : start_candidate_recording()
158 : {
159 0 : UG_COND_THROW(m_candidateRecordingEnabled,
160 : "Call stop_candidate_recording before recording new candidates!");
161 : m_recordedCandidates.clear();
162 0 : m_candidateRecordingEnabled = true;
163 0 : }
164 :
165 : template <class TAAPos>
166 0 : void DelaunayInfo<TAAPos>::
167 : stop_candidate_recording()
168 : {
169 0 : if(m_candidateRecordingEnabled){
170 0 : for(size_t i = 0; i < m_recordedCandidates.size(); ++i){
171 0 : if(m_recordedCandidates[i]){
172 0 : push_candidate(m_recordedCandidates[i]);
173 : }
174 : }
175 : m_recordedCandidates.clear();
176 0 : m_candidateRecordingEnabled = false;
177 : }
178 0 : }
179 :
180 :
181 : template <class TAAPos>
182 0 : void DelaunayInfo<TAAPos>::
183 : enable_face_classification(number minAngle)
184 : {
185 : bool faceClassificationWasEnabled = face_classification_enabled();
186 :
187 0 : m_maxDot = fabs(cos(deg_to_rad(minAngle)));
188 0 : if(minAngle > 0){
189 0 : if(faceClassificationWasEnabled){
190 : // careful - faceInfo->classified has to be set to false!
191 0 : m_faceQueue = FacePriorityQueue();
192 :
193 : // classification was already enabled.
194 0 : for(FaceIterator iter = m_grid.begin<Face>();
195 0 : iter != m_grid.end<Face>(); ++iter)
196 : {
197 0 : if(is_inner(*iter)){
198 0 : m_aaFaceInfo[*iter]->classified = false;
199 0 : classify_face(*iter);
200 : }
201 : }
202 : }
203 : else{
204 0 : m_grid.attach_to_faces_dv(m_aFaceInfo, NULL);
205 0 : m_aaFaceInfo.access(m_grid, m_aFaceInfo);
206 : // we have to create face-infos before classification
207 : for(FaceIterator iter = m_grid.begin<Face>();
208 0 : iter != m_grid.end<Face>(); ++iter)
209 : {
210 0 : if(is_inner(*iter)){
211 0 : m_aaFaceInfo[*iter] = new FaceInfo;
212 0 : m_aaFaceInfo[*iter]->f = *iter;
213 0 : classify_face(*iter);
214 : }
215 : }
216 : }
217 : }
218 0 : else if(faceClassificationWasEnabled){
219 : // unclassify faces.
220 0 : for(FaceIterator iter = m_grid.begin<Face>();
221 0 : iter != m_grid.end<Face>(); ++iter)
222 : {
223 0 : if(m_aaFaceInfo[*iter]){
224 0 : delete m_aaFaceInfo[*iter];
225 : }
226 : }
227 0 : m_faceQueue = FacePriorityQueue();
228 0 : m_grid.detach_from_faces(m_aFaceInfo);
229 : m_aaFaceInfo.invalidate();
230 : }
231 0 : m_minAngle = minAngle;
232 0 : }
233 :
234 : template <class TAAPos>
235 0 : bool DelaunayInfo<TAAPos>::
236 : classified_faces_left()
237 : {
238 : // pop illegal entries
239 0 : while(!m_faceQueue.empty()){
240 0 : FaceInfo* fi = m_faceQueue.top();
241 0 : if(fi->f == NULL){
242 : m_faceQueue.pop();
243 0 : delete fi;
244 : }
245 : else
246 : break;
247 : }
248 0 : return !m_faceQueue.empty();
249 : }
250 :
251 :
252 : template <class TAAPos>
253 0 : Face* DelaunayInfo<TAAPos>::
254 : pop_classified_face()
255 : {
256 : // pop illegal entries
257 0 : while(!m_faceQueue.empty()){
258 0 : FaceInfo* fi = m_faceQueue.top();
259 0 : if(fi->f == NULL){
260 : m_faceQueue.pop();
261 0 : delete fi;
262 : }
263 : else
264 : break;
265 : }
266 :
267 0 : if(m_faceQueue.empty())
268 : return NULL;
269 0 : FaceInfo* fi = m_faceQueue.top();
270 : m_faceQueue.pop();
271 0 : fi->classified = false;
272 0 : return fi->f;
273 : }
274 :
275 :
276 : template <class TAAPos>
277 0 : bool DelaunayInfo<TAAPos>::
278 : is_classified(Face* f)
279 : {
280 0 : if(face_classification_enabled()){
281 0 : return m_aaFaceInfo[f]->classified;
282 : }
283 : return false;
284 : }
285 :
286 :
287 : template <class TAAPos>
288 0 : bool DelaunayInfo<TAAPos>::
289 : is_classifiable(Face* f)
290 : {
291 : int subtended = -1;
292 : int numShell = 0;
293 0 : for(size_t i = 0; i < 3; ++i){
294 0 : if(mark(f->vertex(i)) == SHELL)
295 0 : ++numShell;
296 : else
297 0 : subtended = i;
298 : }
299 :
300 0 : vector_t v[3] = {m_aaPos[f->vertex(0)], m_aaPos[f->vertex(1)], m_aaPos[f->vertex(2)]};
301 :
302 0 : if(numShell == 3){
303 0 : number d1sq = VecDistanceSq(v[0], v[1]);
304 0 : number d2sq = VecDistanceSq(v[1], v[2]);
305 0 : number d3sq = VecDistanceSq(v[2], v[0]);
306 0 : if(d1sq < d2sq){
307 0 : if(d1sq < d3sq)
308 : subtended = 2;
309 : else
310 : subtended = 1;
311 : }
312 0 : else if(d2sq < d3sq)
313 : subtended = 0;
314 : else
315 : subtended = 1;
316 : }
317 :
318 0 : if(numShell >= 2){
319 : vector_t dir1, dir2;
320 0 : VecSubtract(dir1, v[(subtended+1)%3], v[subtended]);
321 0 : VecSubtract(dir2, v[(subtended+2)%3], v[subtended]);
322 0 : if(VecAngle(dir1, dir2) < PI / 3. + SMALL){
323 0 : return false;
324 : }
325 : }
326 : return true;
327 : }
328 :
329 :
330 : template <class TAAPos>
331 0 : bool DelaunayInfo<TAAPos>::
332 : classify_face(Face* f)
333 : {
334 0 : UG_COND_THROW(!face_classification_enabled(),
335 : "Face classification has to be enabled before calling 'classify_face'");
336 0 : FaceInfo* fi = m_aaFaceInfo[f];
337 0 : UG_COND_THROW(!fi, "Invalid face-info attached");
338 0 : UG_COND_THROW(fi->classified, "Only unclassified faces may be classified");
339 :
340 : // only triangles can be classified
341 0 : if(f->num_vertices() != 3){
342 0 : fi->classified = false;
343 0 : return false;
344 : }
345 :
346 0 : if(!is_classifiable(f))
347 : return false;
348 :
349 0 : vector_t& v1 = m_aaPos[f->vertex(0)];
350 0 : vector_t& v2 = m_aaPos[f->vertex(1)];
351 0 : vector_t& v3 = m_aaPos[f->vertex(2)];
352 :
353 : // if at least two of the vertices are SHELL vertices, we won't classify the triangle if
354 : // the subtended angle to the shortest edge between shell vertices is smaller than PI/3
355 : // {
356 : // int subtended = -1;
357 : // int numShell = 0;
358 : // for(size_t i = 0; i < 3; ++i){
359 : // if(mark(f->vertex(i)) == SHELL)
360 : // ++numShell;
361 : // else
362 : // subtended = i;
363 : // }
364 :
365 : // if(numShell == 3){
366 : // number d1sq = VecDistanceSq(v1, v2);
367 : // number d2sq = VecDistanceSq(v2, v3);
368 : // number d3sq = VecDistanceSq(v3, v1);
369 : // if(d1sq < d2sq){
370 : // if(d1sq < d3sq)
371 : // subtended = 2;
372 : // else
373 : // subtended = 1;
374 : // }
375 : // else if(d2sq < d3sq)
376 : // subtended = 0;
377 : // else
378 : // subtended = 1;
379 : // }
380 :
381 : // if(numShell >= 2){
382 : // vector_t dir1, dir2;
383 : // VecSubtract(dir1, m_aaPos[f->vertex((subtended+1)%3)], m_aaPos[f->vertex(subtended)]);
384 : // VecSubtract(dir2, m_aaPos[f->vertex((subtended+2)%3)], m_aaPos[f->vertex(subtended)]);
385 : // if(VecAngle(dir1, dir2) < PI / 3. + SMALL){
386 : // return false;
387 : // }
388 : // }
389 : // }
390 :
391 : // perform classification
392 : // calculate min angle
393 : vector_t v12, v13, v23;
394 0 : VecSubtract(v12, v2, v1); VecNormalize(v12, v12);
395 0 : VecSubtract(v13, v3, v1); VecNormalize(v13, v13);
396 0 : VecSubtract(v23, v3, v2); VecNormalize(v23, v23);
397 :
398 0 : number d1 = fabs(VecDot(v12, v13));
399 0 : number d2 = fabs(VecDot(v12, v23));
400 0 : number d3 = fabs(VecDot(v13, v23));
401 :
402 : number highestDot = 0;
403 0 : if(d1 > m_maxDot){
404 : // check edges
405 0 : Edge* e1 = m_grid.get_edge(f, 0);
406 0 : Edge* e2 = m_grid.get_edge(f, 2);
407 0 : if(!(is_segment(e1) && is_segment(e2)))
408 : highestDot = d1;
409 : }
410 :
411 0 : if((highestDot < m_maxDot) && (d2 > m_maxDot)){
412 : // check edges
413 0 : Edge* e1 = m_grid.get_edge(f, 0);
414 0 : Edge* e2 = m_grid.get_edge(f, 1);
415 0 : if(!(is_segment(e1) && is_segment(e2)))
416 : highestDot = d2;
417 : }
418 :
419 0 : if((highestDot < m_maxDot) && (d3 > m_maxDot)){
420 : // check edges
421 0 : Edge* e1 = m_grid.get_edge(f, 1);
422 0 : Edge* e2 = m_grid.get_edge(f, 2);
423 0 : if(!(is_segment(e1) && is_segment(e2)))
424 : highestDot = d3;
425 : }
426 :
427 0 : if(highestDot > m_maxDot){
428 : // calculate the radius of the circumcenter for the priority
429 : vector_t cc;
430 0 : if(TriangleCircumcenter(cc, v1, v2, v3)){
431 0 : fi->priority = VecDistanceSq(cc, v1);
432 0 : fi->classified = true;
433 0 : m_faceQueue.push(fi);
434 : }
435 : else{
436 : // UG_LOG("Couldn't calculate triangle-circumcenter. ");
437 : // UG_LOG("Ignoring triangle with center at: " << CalculateCenter(f, m_aaPos) << "\n");
438 0 : return false;
439 : }
440 : }
441 :
442 : //todo if the face-queue gets too large, we should do some clean up
443 : return true;
444 : }
445 :
446 :
447 : template <class TAAPos>
448 0 : void DelaunayInfo<TAAPos>::
449 : vertex_created(Grid* grid, Vertex* vrt, GridObject* pParent, bool replacesParent)
450 : {
451 :
452 : set_mark(vrt, NEW_INNER);
453 0 : if(pParent && (pParent->base_object_id() == EDGE)
454 0 : && is_segment(static_cast<Edge*>(pParent)))
455 : {
456 : set_mark(vrt, NEW_SEGMENT);
457 : }
458 0 : }
459 :
460 :
461 : template <class TAAPos>
462 0 : void DelaunayInfo<TAAPos>::
463 : edge_created(Grid* grid, Edge* e, GridObject* pParent, bool replacesParent)
464 : {
465 :
466 0 : set_mark(e, NEW_INNER);
467 0 : if(pParent && (pParent->base_object_id() == EDGE)
468 0 : && is_segment(static_cast<Edge*>(pParent)))
469 : {
470 : set_mark(e, NEW_SEGMENT);
471 : }
472 :
473 0 : if(m_candidateRecordingEnabled)
474 0 : m_recordedCandidates.push_back(e);
475 0 : }
476 :
477 :
478 : template <class TAAPos>
479 0 : void DelaunayInfo<TAAPos>::
480 : face_created(Grid* grid, Face* f, GridObject* pParent, bool replacesParent)
481 : {
482 0 : if(pParent && (pParent->base_object_id() == FACE)
483 0 : && is_inner(static_cast<Face*>(pParent)))
484 : {
485 0 : set_mark(f, NEW_INNER);
486 : }
487 0 : }
488 :
489 :
490 : template <class TAAPos>
491 0 : void DelaunayInfo<TAAPos>::
492 : edge_to_be_erased(Grid* grid, Edge* e, Edge* replacedBy)
493 : {
494 0 : if(m_candidateRecordingEnabled){
495 : // if e is a recorded candidate, we'll simply overwrite it
496 : // with the last valid candidate and shorten the candidates array.
497 0 : for(size_t i = 0; i < m_recordedCandidates.size(); ++i){
498 0 : if(m_recordedCandidates[i] == e){
499 0 : m_recordedCandidates[i] = m_recordedCandidates.back();
500 : m_recordedCandidates.pop_back();
501 : break;
502 : }
503 : }
504 : }
505 0 : }
506 :
507 : template <class TAAPos>
508 0 : void DelaunayInfo<TAAPos>::
509 : face_to_be_erased(Grid* grid, Face* f, Face* replacedBy)
510 : {
511 : // unmark the face.
512 0 : set_mark(f, NONE);
513 0 : }
514 :
515 :
516 : template class DelaunayInfo<Grid::VertexAttachmentAccessor<AVector3> >;
517 :
518 : }// end of namespace
|
