Line data Source code
1 : /*
2 : * Copyright (c) 2017: 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__UG_orientation_util_impl
34 : #define __H__UG_orientation_util_impl
35 :
36 : #include "common/math/ugmath.h"
37 : #include "geom_obj_util/face_util.h"
38 : #include "geom_obj_util/volume_util.h"
39 :
40 : #include <stack>
41 : #include <vector>
42 :
43 :
44 : namespace ug{
45 :
46 :
47 : ////////////////////////////////////////////////////////////////////////
48 : // InvertOrientation
49 : template <class iter_t>
50 : void InvertOrientation(Grid& grid, iter_t elemsBegin,
51 : iter_t elemsEnd)
52 : {
53 : for(iter_t iter = elemsBegin; iter != elemsEnd; ++iter)
54 : grid.flip_orientation(*iter);
55 : }
56 :
57 :
58 : ////////////////////////////////////////////////////////////////////////
59 : // FixFaceOrientation
60 : template <class TFaceIterator>
61 0 : void FixFaceOrientation(Grid& grid, TFaceIterator facesBegin,
62 : TFaceIterator facesEnd)
63 : {
64 : using namespace std;
65 :
66 : {
67 : // make sure that boundary faces are oriented outwards
68 : Grid::volume_traits::secure_container vols;
69 0 : for(TFaceIterator iter = facesBegin; iter != facesEnd; ++iter){
70 : grid.associated_elements(vols, *iter);
71 0 : if(vols.size() == 1 && !OrientationMatches(*iter, vols[0])){
72 0 : grid.flip_orientation(*iter);
73 : }
74 : }
75 : }
76 :
77 : // the rest only considers manifold faces
78 : // we use marks to differentiate between processed and unprocessed faces
79 0 : grid.begin_marking();
80 :
81 : // we have to mark all faces between facesBegin and facesEnd initially,
82 : // to differentiate them from the other faces in the grid.
83 0 : for(TFaceIterator iter = facesBegin; iter != facesEnd; ++iter){
84 0 : if(NumAssociatedVolumes(grid, *iter) == 0)
85 : grid.mark(*iter);
86 : }
87 :
88 : // this edge descriptor will be used multiple times
89 0 : EdgeDescriptor ed;
90 :
91 : // containers to store neighbour elements.
92 : vector<Face*> vNeighbours;
93 :
94 : // stack that stores candidates
95 : stack<Face*> stkFaces;
96 :
97 : // we'll iterate through all faces
98 0 : while(facesBegin != facesEnd)
99 : {
100 : // candidates are empty at this point.
101 : // if the face is unprocessed it is a new candidate
102 0 : if(grid.is_marked(*facesBegin))
103 : {
104 : // mark it as candidate (by removing the mark)
105 : grid.unmark(*facesBegin);
106 0 : stkFaces.push(*facesBegin);
107 :
108 : // while the stack is not empty
109 0 : while(!stkFaces.empty())
110 : {
111 : // get the candidate
112 0 : Face* f = stkFaces.top();
113 : stkFaces.pop();
114 :
115 : // get the neighbours for each side
116 0 : for(size_t i = 0; i < f->num_edges(); ++i)
117 : {
118 0 : f->edge_desc(i, ed);
119 0 : GetNeighbours(vNeighbours, grid, f, i);
120 :
121 : // fix orientation of unprocessed neighbours.
122 0 : for(size_t j = 0; j < vNeighbours.size(); ++j)
123 : {
124 0 : Face* fn = vNeighbours[j];
125 0 : if(grid.is_marked(fn))
126 : {
127 : // check whether the orientation of f and fn differs.
128 0 : if(EdgeOrientationMatches(&ed, fn))
129 : {
130 : // the orientation of ed is the same as the orientation
131 : // of an edge in fn.
132 : // the faces thus have different orientation.
133 0 : grid.flip_orientation(fn);
134 : }
135 :
136 : // mark the face as processed and add it to the stack
137 : grid.unmark(fn);
138 : stkFaces.push(fn);
139 : }
140 : }
141 : }
142 : }
143 : }
144 :
145 : // check the next face
146 : ++facesBegin;
147 : }
148 :
149 0 : grid.end_marking();
150 0 : }
151 :
152 :
153 : template<class TAAPosVRT>
154 : bool
155 0 : CheckOrientation(Volume* vol, TAAPosVRT& aaPosVRT)
156 : {
157 : // some typedefs
158 : typedef typename TAAPosVRT::ValueType vector_t;
159 :
160 : // First calculate the center of the volume
161 0 : vector_t volCenter = CalculateCenter(vol, aaPosVRT);
162 :
163 : // now check for each side whether it points away from the center.
164 0 : size_t numFaces = vol->num_faces();
165 0 : FaceDescriptor fd;
166 : vector_t normal;
167 0 : for(size_t i = 0; i < numFaces; ++i){
168 0 : vol->face_desc(i, fd);
169 0 : CalculateNormal(normal, &fd, aaPosVRT);
170 :
171 : // in order to best approximate quadrilateral faces, we'll calculate the
172 : // center of the face and compare that to the volCenter.
173 0 : vector_t faceCenter = CalculateCenter(&fd, aaPosVRT);
174 :
175 : // now compare normal and center
176 : vector_t dir;
177 : VecSubtract(dir, faceCenter, volCenter);
178 0 : if(VecDot(dir, normal) < 0)
179 0 : return false;
180 : }
181 :
182 : // all center / normal checks succeeded. Orientation is fine.
183 : return true;
184 : }
185 :
186 : template<class TVolIterator, class TAAPosVRT>
187 : int
188 : FixOrientation(Grid& grid, TVolIterator volsBegin, TVolIterator volsEnd,
189 : TAAPosVRT& aaPosVRT)
190 : {
191 : int numFlips = 0;
192 : // iterate through all volumes
193 : for(TVolIterator iter = volsBegin; iter != volsEnd; ++iter){
194 : // check whether the orientation is fine
195 : if(!CheckOrientation(*iter, aaPosVRT)){
196 : grid.flip_orientation(*iter);
197 : ++numFlips;
198 : }
199 : }
200 :
201 : return numFlips;
202 : }
203 :
204 :
205 : }// end of namespace
206 :
207 : #endif //__H__UG_orientation_util_impl
|
