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 : #ifndef __H__LIB_GRID__EDGE_UTIL_IMPL__
34 : #define __H__LIB_GRID__EDGE_UTIL_IMPL__
35 :
36 : //#include "edge_util.h"
37 : #include <stack>
38 : #include <queue>
39 : #include <map>
40 : #include "lib_grid/grid/grid_util.h"
41 : #include "../orientation_util.h"
42 : #include "vertex_util.h"
43 : namespace ug
44 : {
45 :
46 : template <class face_iter_t>
47 : void GetInnerEdgesOfFaceSoup(
48 : std::vector<Edge*>& edgesOut,
49 : Grid& g,
50 : face_iter_t facesBegin,
51 : face_iter_t facesEnd)
52 : {
53 : using namespace std;
54 : map<Edge*, int> m;
55 :
56 : Grid::edge_traits::secure_container edges;
57 : for(face_iter_t fiter = facesBegin; fiter != facesEnd; ++fiter)
58 : {
59 : g.associated_elements(edges, *fiter);
60 : for(size_t i = 0; i < edges.size(); ++i)
61 : ++m[edges[i]];
62 : }
63 :
64 : for(map<Edge*, int>::iterator iter = m.begin(); iter != m.end(); ++iter)
65 : {
66 : if(iter->second > 1)
67 : edgesOut.push_back(iter->first);
68 : }
69 : }
70 :
71 : ////////////////////////////////////////////////////////////////////////
72 : template <class TAAPosVRT>
73 0 : inline number EdgeLengthSq(const EdgeVertices* e, TAAPosVRT& aaPos)
74 : {
75 0 : return VecDistanceSq(aaPos[e->vertex(0)], aaPos[e->vertex(1)]);
76 : }
77 :
78 : ////////////////////////////////////////////////////////////////////////
79 : template <class TAAPosVRT>
80 0 : inline number EdgeLength(const EdgeVertices* e, TAAPosVRT& aaPos)
81 : {
82 0 : return VecDistance(aaPos[e->vertex(0)], aaPos[e->vertex(1)]);
83 : }
84 :
85 : ////////////////////////////////////////////////////////////////////////
86 : // SplitEdge
87 : // see edge_operations.h for detailed description
88 : template<class TVertex>
89 : TVertex* SplitEdge(Grid& grid, Edge* e, bool bConservative)
90 : {
91 0 : return SplitEdge<TVertex>(grid, grid, e, NULL, bConservative);
92 : }
93 :
94 : ////////////////////////////////////////////////////////////////////////
95 : // SplitEdge
96 : // see edge_operations.h for detailed description
97 : template<class TVertex>
98 0 : TVertex* SplitEdge(Grid& destGrid, Grid& srcGrid, Edge* e,
99 : AVertex* paAssociatedVertices,
100 : bool bConservative)
101 : {
102 : TVertex* newVertex;
103 0 : if(&destGrid == &srcGrid)
104 0 : newVertex = *destGrid.create<TVertex>(e);
105 : else
106 0 : newVertex = *destGrid.create<TVertex>();
107 :
108 0 : if(CreateEdgeSplitGeometry(destGrid, srcGrid, e, newVertex, paAssociatedVertices))
109 : {
110 0 : if(!bConservative)
111 : {
112 : // erase unused elements.
113 0 : if(!srcGrid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES))
114 : {
115 : // we have to erase the faces manually
116 : // collect them
117 : std::vector<Face*> vFaces;
118 0 : CollectFaces(vFaces, srcGrid, e, false);
119 :
120 : // erase them
121 0 : for(std::vector<Face*>::iterator iter = vFaces.begin();
122 0 : iter != vFaces.end(); ++iter)
123 : {
124 0 : srcGrid.erase(*iter);
125 : }
126 0 : }
127 :
128 0 : if((!srcGrid.option_is_enabled(VOLOPT_AUTOGENERATE_EDGES)) &&
129 0 : (!srcGrid.option_is_enabled(VOLOPT_AUTOGENERATE_FACES)))
130 : {
131 : // we have to erase them manually
132 : // collect them
133 : std::vector<Volume*> vVolumes;
134 0 : CollectVolumes(vVolumes, srcGrid, e, false);
135 :
136 : // erase them
137 0 : for(std::vector<Volume*>::iterator iter = vVolumes.begin();
138 0 : iter != vVolumes.end(); ++iter)
139 : {
140 0 : srcGrid.erase(*iter);
141 : }
142 0 : }
143 : // erase the edge
144 0 : srcGrid.erase(e);
145 : }
146 :
147 : // return the new vertex
148 0 : return newVertex;
149 : }
150 :
151 : // something went wrong in CreateEdgeSplitGeometry.
152 : // erase the new vertex and return NULL
153 0 : destGrid.erase(newVertex);
154 0 : return NULL;
155 : }
156 :
157 : template<class TVertexPositionAttachmentAccessor>
158 : typename TVertexPositionAttachmentAccessor::ValueType
159 0 : CalculateCenter(const Edge* e, TVertexPositionAttachmentAccessor& aaPosVRT)
160 : {
161 : typename TVertexPositionAttachmentAccessor::ValueType v;
162 : // init v with 0.
163 : VecSet(v, 0);
164 :
165 : // sum up
166 0 : VecAdd(v, aaPosVRT[e->vertex(0)], aaPosVRT[e->vertex(1)]);
167 :
168 : // average
169 : VecScale(v, v, 0.5);
170 :
171 0 : return v;
172 : }
173 :
174 : ////////////////////////////////////////////////////////////////////////
175 : template<class TAAPosVRT, class TAAWeightVRT>
176 : UG_API
177 : typename TAAPosVRT::ValueType
178 : CalculateCenter(const EdgeVertices* e, TAAPosVRT& aaPos, TAAWeightVRT& aaWeight)
179 : {
180 : typename TAAPosVRT::ValueType v;
181 : typedef typename TAAWeightVRT::ValueType weight_t;
182 :
183 : // init v with 0.
184 : VecSet(v, 0);
185 :
186 : // sum up
187 : weight_t w0 = aaWeight[e->vertex(0)];
188 : weight_t w1 = aaWeight[e->vertex(1)];
189 :
190 : VecScaleAdd(v, w0, aaPos[e->vertex(0)], w1, aaPos[e->vertex(1)]);
191 :
192 : if((w0 + w1) != 0)
193 : VecScale(v, v, 1. / (number)(w0 + w1));
194 :
195 : return v;
196 : }
197 :
198 : template <class TEdgeIterator>
199 0 : void FixEdgeOrientation(Grid& grid, TEdgeIterator edgesBegin,
200 : TEdgeIterator edgesEnd)
201 : {
202 : using namespace std;
203 0 : grid.begin_marking();
204 :
205 : // we'll mark all edges.
206 0 : for(TEdgeIterator iter = edgesBegin; iter != edgesEnd; ++iter)
207 : grid.mark(*iter);
208 :
209 : stack<Edge*> candidates;
210 :
211 : // iterate through the edges
212 0 : for(TEdgeIterator iter = edgesBegin; iter != edgesEnd; ++iter)
213 : {
214 : // only marked edges have to be considered
215 0 : if(grid.is_marked(*iter)){
216 : grid.unmark(*iter);
217 0 : candidates.push(*iter);
218 :
219 0 : while(!candidates.empty()){
220 0 : Edge* e = candidates.top();
221 : candidates.pop();
222 :
223 : // iterate through all associated edges
224 0 : for(size_t i = 0; i < 2; ++i){
225 : for(Grid::AssociatedEdgeIterator assIter =
226 0 : grid.associated_edges_begin(e->vertex(i));
227 0 : assIter != grid.associated_edges_end(e->vertex(i)); ++assIter)
228 : {
229 0 : Edge* ae = *assIter;
230 : // fix orientation of marked associated edges.
231 : // those edges are new candidates afterwards
232 0 : if(grid.is_marked(ae)){
233 : // the local index of the vertex over which ae is connected to e has
234 : // to be different from the local index of the vertex in e.
235 0 : if(GetVertexIndex(ae, e->vertex(i)) == (int)i){
236 : // we have to flip the orientation
237 0 : grid.flip_orientation(ae);
238 : }
239 :
240 : // prepare the new candidate
241 : grid.unmark(ae);
242 : candidates.push(ae);
243 : }
244 : }
245 : }
246 : }
247 : }
248 : }
249 :
250 0 : grid.end_marking();
251 0 : }
252 :
253 :
254 : template <class TEdgeIterator>
255 : UG_API
256 : void AdjustEdgeOrientationToFaceOrientation(Grid& grid, TEdgeIterator edgesBegin,
257 : TEdgeIterator edgesEnd)
258 : {
259 : using namespace std;
260 :
261 : for(TEdgeIterator iter = edgesBegin; iter != edgesEnd; ++iter){
262 : Face* nbr;
263 : int numNbrs = GetAssociatedFaces(&nbr, grid, *iter, 1);
264 : if(numNbrs != 1)
265 : continue;
266 :
267 : if(!EdgeOrientationMatches(*iter, nbr)){
268 : grid.flip_orientation(*iter);
269 : }
270 : }
271 : }
272 :
273 : template <class TEdgeIterator>
274 : UG_API
275 : void AdjustEdgeOrientationToFaceOrientation(Grid& grid, TEdgeIterator edgesBegin,
276 : TEdgeIterator edgesEnd,
277 : Grid::face_traits::callback considerFace)
278 : {
279 : using namespace std;
280 :
281 : Grid::face_traits::secure_container faces;
282 :
283 : for(TEdgeIterator iter = edgesBegin; iter != edgesEnd; ++iter){
284 : Edge* e = *iter;
285 : grid.associated_elements(faces, e);
286 :
287 : size_t numConsidered = 0;
288 : Face* nbr = NULL;
289 : for(size_t i = 0; i < faces.size(); ++i){
290 : if(considerFace(faces[i])){
291 : ++numConsidered;
292 : nbr = faces[i];
293 : }
294 : }
295 :
296 : if((numConsidered == 1) && (!EdgeOrientationMatches(e, nbr))){
297 : grid.flip_orientation(e);
298 : }
299 : }
300 : }
301 :
302 :
303 : template <class TEdgeIterator, class TAAPosVRT>
304 : Edge* FindShortestEdge(TEdgeIterator edgesBegin, TEdgeIterator edgesEnd,
305 : TAAPosVRT& aaPos)
306 : {
307 : // if edgesBegin equals edgesEnd, then the list is empty and we can
308 : // immediately return NULL
309 : if(edgesBegin == edgesEnd)
310 : return NULL;
311 :
312 : // the first edge is the first candidate for the shortest edge.
313 : // We compare squares to avoid computation of the square root.
314 : Edge* shortestEdge = *edgesBegin;
315 : number shortestLen = EdgeLengthSq(shortestEdge, aaPos);
316 : ++edgesBegin;
317 :
318 : for(; edgesBegin != edgesEnd; ++edgesBegin){
319 : Edge* curEdge = *edgesBegin;
320 : number curLen = EdgeLengthSq(curEdge, aaPos);
321 : if(curLen < shortestLen){
322 : shortestEdge = curEdge;
323 : shortestLen = curLen;
324 : }
325 : }
326 :
327 : return shortestEdge;
328 : }
329 :
330 :
331 : template <class TEdgeIterator, class TAAPosVRT>
332 0 : Edge* FindLongestEdge(TEdgeIterator edgesBegin, TEdgeIterator edgesEnd,
333 : TAAPosVRT& aaPos)
334 : {
335 : // if edgesBegin equals edgesEnd, then the list is empty and we can
336 : // immediately return NULL
337 0 : if(edgesBegin == edgesEnd)
338 : return NULL;
339 :
340 : // the first edge is the first candidate for the shortest edge.
341 : // We compare squares to avoid computation of the square root.
342 0 : Edge* longestEdge = *edgesBegin;
343 0 : number longestLen = EdgeLengthSq(longestEdge, aaPos);
344 : ++edgesBegin;
345 :
346 0 : for(; edgesBegin != edgesEnd; ++edgesBegin){
347 0 : Edge* curEdge = *edgesBegin;
348 0 : number curLen = EdgeLengthSq(curEdge, aaPos);
349 0 : if(curLen > longestLen){
350 : longestEdge = curEdge;
351 : longestLen = curLen;
352 : }
353 : }
354 :
355 : return longestEdge;
356 : }
357 : // ////////////////////////////////////////////////////////////////////////////////
358 : // template <class TEdgeIterator>
359 : // void RemoveDoubleEdges(Grid& grid, TEdgeIterator edgesBegin, TEdgeIterator edgesEnd)
360 : // {
361 : // // iterate over all edges and check whether both associated vertices are already
362 : // // marked. If this is the case and a marked edge exists between those vertices,
363 : // // then the current edge is a double.
364 :
365 : // grid.begin_marking();
366 :
367 : // // the first is the double, the second the original
368 : // std::vector<std::pair<Edge*, Edge*> > doubles;
369 : // std::vector<Edge*> edges;
370 :
371 : // for(TEdgeIterator iter = edgesBegin; iter != edgesEnd; ++iter){
372 : // Edge* e = *iter;
373 : // if(!grid.is_marked(e)){
374 : // // check whether both vertices are marked
375 : // Vertex* v0 = e->vertex(0);
376 : // Vertex* v1 = e->vertex(1);
377 :
378 : // bool isDouble = false;
379 :
380 : // if(grid.is_marked(v0) && grid.is_marked(v1)){
381 : // // a necessary condition is met. However not yet sufficient.
382 : // // find marked edge between v0 and v1.
383 : // CollectAssociated(edges, grid, v0);
384 : // for(size_t i = 0; i < edges.size(); ++i){
385 : // Edge* te = edges[i];
386 : // if((te->vertex(0) == v1 || te->vertex(1) == v1)
387 : // && grid.is_marked(te))
388 : // {
389 : // // e is a double
390 : // isDouble = true;
391 : // doubles.push_back(std::make_pair(e, te));
392 : // break;
393 : // }
394 : // }
395 : // }
396 :
397 : // // finally mark e and its vertices (every processed edge is marked).
398 : // if(!isDouble){
399 : // grid.mark(e);
400 : // grid.mark(v0);
401 : // grid.mark(v1);
402 : // }
403 : // }
404 : // }
405 :
406 : // grid.end_marking();
407 :
408 : // // now erase all doubles
409 : // for(size_t i = 0; i < doubles.size(); ++i){
410 : // // this allows listeners to take actions
411 : // grid.objects_will_be_merged(doubles[i].second, doubles[i].second,
412 : // doubles[i].first);
413 : // grid.erase(doubles[i].first);
414 : // }
415 : // }
416 :
417 : ////////////////////////////////////////////////////////////////////////////////
418 : template <class EdgeIterator, class TAAPos>
419 : void MinimizeEdgeLength_SwapsOnly(Grid& grid, EdgeIterator edgesBegin,
420 : EdgeIterator edgesEnd, TAAPos& aaPos)
421 : {
422 : using namespace std;
423 :
424 : // helper to collect neighbors
425 : Face* nbrFaces[2];
426 : vector<Edge*> edges;
427 :
428 : // flipCandidates
429 : queue<Edge*> candidates;
430 :
431 : // sadly we can't use marking. Thats why we attach a simple byte to the edges,
432 : // which will tell whether an edge is already a candidate.
433 : AByte aIsCandidate;
434 : grid.attach_to_edges_dv(aIsCandidate, 0, false);
435 : Grid::AttachmentAccessor<Edge, AByte> aaIsCandidate(grid, aIsCandidate);
436 :
437 : // set up candidate array
438 : for(EdgeIterator iter = edgesBegin; iter != edgesEnd; ++iter){
439 : aaIsCandidate[*iter] = 1;
440 : candidates.push(*iter);
441 : }
442 :
443 :
444 : while(!candidates.empty()){
445 : Edge* e = candidates.front();
446 : candidates.pop();
447 : aaIsCandidate[e] = 0;
448 :
449 : // we only perform swaps on regular manifolds.
450 : if(GetAssociatedFaces(nbrFaces, grid, e, 2) == 2){
451 : // make sure that both neighbors are triangles
452 : if(nbrFaces[0]->num_vertices() != 3 || nbrFaces[1]->num_vertices() != 3)
453 : continue;
454 :
455 : // check whether a swap would make the edge shorter.
456 : Vertex* conVrt0 = GetConnectedVertex(e, nbrFaces[0]);
457 : Vertex* conVrt1 = GetConnectedVertex(e, nbrFaces[1]);
458 : if(VertexDistanceSq(conVrt0, conVrt1, aaPos) < EdgeLengthSq(e, aaPos))
459 : {
460 : // it'll be shorter
461 : // now make sure that associated triangles won't flip
462 : //todo: add support for 2d position attachments
463 : vector3 n0, n1;
464 : CalculateNormal(n0, nbrFaces[0], aaPos);
465 : CalculateNormal(n1, nbrFaces[1], aaPos);
466 : number oldDot = VecDot(n0, n1);
467 :
468 : FaceDescriptor ntri;
469 : ntri.set_num_vertices(3);
470 : ntri.set_vertex(0, e->vertex(0));
471 : ntri.set_vertex(1, conVrt1);
472 : ntri.set_vertex(2, conVrt0);
473 : CalculateNormal(n0, &ntri, aaPos);
474 :
475 : ntri.set_vertex(0, e->vertex(1));
476 : ntri.set_vertex(1, conVrt0);
477 : ntri.set_vertex(2, conVrt1);
478 : CalculateNormal(n1, &ntri, aaPos);
479 :
480 : number newDot = VecDot(n0, n1);
481 :
482 : // if both have the same sign, we're fine!
483 : if(oldDot * newDot < 0){
484 : continue;// not fine!
485 : }
486 :
487 : // ok - everything is fine. Now swap the edge
488 : e = SwapEdge(grid, e);
489 :
490 : UG_ASSERT(e, "SwapEdge did not produce a new edge.");
491 :
492 : // all edges of associated triangles are candidates again (except e)
493 : GetAssociatedFaces(nbrFaces, grid, e, 2);
494 : for(size_t i = 0; i < 2; ++i){
495 : CollectAssociated(edges, grid, nbrFaces[i]);
496 : for(size_t j = 0; j < edges.size(); ++j){
497 : if(edges[j] != e && (!aaIsCandidate[edges[j]])){
498 : candidates.push(edges[j]);
499 : aaIsCandidate[edges[j]] = 1;
500 : }
501 : }
502 : }
503 : }
504 : }
505 : }
506 :
507 : grid.detach_from_edges(aIsCandidate);
508 : }
509 :
510 : ////////////////////////////////////////////////////////////////////////
511 : template <class vector_t, class TAAPos>
512 : UG_API bool
513 0 : ContainsPoint(const EdgeVertices* e, const vector_t& p, TAAPos aaPos)
514 : {
515 0 : number center = (aaPos[e->vertex(0)].x() + aaPos[e->vertex(1)].x()) / 2.;
516 0 : number rad = fabs(aaPos[e->vertex(1)].x() - aaPos[e->vertex(0)].x()) / 2.;
517 :
518 0 : return (fabs(p.x() - center) < rad + SMALL);
519 : }
520 :
521 : ////////////////////////////////////////////////////////////////////////
522 : template <class TAAPosVRT>
523 : number CalculateAverageEdgeLength(Grid& grid, TAAPosVRT& aaPos)
524 : {
525 : ConstEdgeIterator eit = grid.begin<Edge>();
526 : ConstEdgeIterator eit_end = grid.end<Edge>();
527 : number avg_length = 0.;
528 : for (; eit != eit_end; ++eit) {
529 : avg_length += EdgeLength(*eit, aaPos);
530 : }
531 : return avg_length / grid.num_edges();
532 : }
533 :
534 :
535 : }// end of namespace
536 :
537 : #endif
|
