Line data Source code
1 : /*
2 : * Copyright (c) 2010-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__FACE_UTIL_IMPL__
34 : #define __H__LIB_GRID__FACE_UTIL_IMPL__
35 :
36 : #include "face_util.h"
37 :
38 : namespace ug
39 : {
40 :
41 : ////////////////////////////////////////////////////////////////////////
42 : template <class TIterator>
43 : number AreaFaceQuality(TIterator facesBegin, TIterator facesEnd,
44 : Grid::VertexAttachmentAccessor<APosition>& aaPos)
45 : {
46 : // if the area is empty return 0 (bad)
47 : if(facesBegin == facesEnd)
48 : return 0;
49 :
50 : // get the first
51 : number q = FaceQuality(*facesBegin, aaPos);
52 : ++facesBegin;
53 :
54 : // iterate over the others and find a worse one
55 : for(; facesBegin != facesEnd; ++facesBegin){
56 : number tq = FaceQuality(*facesBegin, aaPos);
57 : if(tq < q)
58 : q = tq;
59 : }
60 :
61 : // return the quality
62 : return q;
63 : }
64 :
65 : ////////////////////////////////////////////////////////////////////////
66 : inline void Triangulate(Grid& grid,
67 : Grid::VertexAttachmentAccessor<APosition>* paaPos)
68 : {
69 : Triangulate(grid, grid.begin<Quadrilateral>(),
70 : grid.end<Quadrilateral>(), paaPos);
71 : }
72 :
73 :
74 : ////////////////////////////////////////////////////////////////////////
75 : template<class TVertexPositionAttachmentAccessor>
76 : typename TVertexPositionAttachmentAccessor::ValueType
77 0 : CalculateCenter(const FaceVertices* f, TVertexPositionAttachmentAccessor& aaPosVRT)
78 : {
79 0 : const size_t numVrts = f->num_vertices();
80 : typename TVertexPositionAttachmentAccessor::ValueType v;
81 : // init v with 0.
82 : VecSet(v, 0);
83 :
84 0 : FaceVertices::ConstVertexArray vrts = f->vertices();
85 :
86 : // sum up
87 0 : for(size_t i = 0; i < numVrts; ++i)
88 : {
89 0 : VecAdd(v, v, aaPosVRT[vrts[i]]);
90 : }
91 :
92 : // average
93 0 : if(numVrts > 0)
94 0 : VecScale(v, v, 1./(number)numVrts);
95 :
96 0 : return v;
97 : }
98 :
99 :
100 : ////////////////////////////////////////////////////////////////////////
101 : template<class TAAPosVRT, class TAAWeightVRT>
102 : UG_API
103 : typename TAAPosVRT::ValueType
104 : CalculateCenter(const FaceVertices* f, TAAPosVRT& aaPos, TAAWeightVRT& aaWeight)
105 : {
106 : uint numVrts = f->num_vertices();
107 : typename TAAPosVRT::ValueType v;
108 : typedef typename TAAWeightVRT::ValueType weight_t;
109 : // init v with 0.
110 : VecSet(v, 0);
111 :
112 : FaceVertices::ConstVertexArray vrts = f->vertices();
113 :
114 : // sum up
115 : weight_t totalWeight = 0;
116 : for(uint i = 0; i < numVrts; ++i)
117 : {
118 : weight_t w = aaWeight[vrts[i]];
119 : VecScaleAppend(v, w, aaPos[vrts[i]]);
120 : totalWeight += w;
121 : }
122 :
123 : // average
124 : if(totalWeight != 0)
125 : VecScale(v, v, 1./(number)totalWeight);
126 :
127 : return v;
128 : }
129 :
130 : ////////////////////////////////////////////////////////////////////////
131 : template <class vector_t, class TAAPos>
132 : bool
133 0 : ContainsPoint(const FaceVertices* f, const vector_t& p, TAAPos aaPos)
134 : {
135 0 : switch(f->num_vertices()){
136 0 : case 3: return PointIsInsideTriangle(p, aaPos[f->vertex(0)],
137 0 : aaPos[f->vertex(1)],
138 0 : aaPos[f->vertex(2)]);
139 0 : case 4: return PointIsInsideQuadrilateral(p, aaPos[f->vertex(0)],
140 0 : aaPos[f->vertex(1)],
141 0 : aaPos[f->vertex(2)],
142 0 : aaPos[f->vertex(3)]);
143 0 : default:
144 0 : UG_THROW("Unknown face type with " << f->num_vertices()
145 : << " vertices encountered in ContainsPoint(...).");
146 : }
147 : return false;
148 : }
149 :
150 : ////////////////////////////////////////////////////////////////////////
151 : // project points to surface
152 : template <class TTriangleIterator, class TAAPosVRT>
153 : bool ProjectPointToSurface(vector3& vOut, const vector3& v, const vector3& n,
154 : TTriangleIterator trisBegin, TTriangleIterator trisEnd,
155 : TAAPosVRT& aaPos, bool compareNormals)
156 : {
157 : bool gotOne = false;
158 : number bestDist = 0; // value doesn't matter - will be overwritten later on.
159 :
160 : // iterate through all triangles and find the closest intersection
161 : for(TTriangleIterator iter = trisBegin; iter != trisEnd; ++iter)
162 : {
163 : Triangle* tri = *iter;
164 :
165 : const vector3& p1 = aaPos[tri->vertex(0)];
166 : const vector3& p2 = aaPos[tri->vertex(1)];
167 : const vector3& p3 = aaPos[tri->vertex(2)];
168 :
169 : vector3 tn;
170 : CalculateTriangleNormalNoNormalize(tn, p1, p2, p3);
171 :
172 : // if normal-check is enabled, we have to make sure, that the points
173 : // normal points into the same direction as the triangles normal.
174 : if(compareNormals){
175 : if(VecDot(tn, n) <= 0)
176 : continue;
177 : }
178 :
179 : number bc1, bc2;
180 : vector3 vTmp;
181 : number distance = DistancePointToTriangle(vTmp, bc1, bc2,
182 : v, p1, p2, p3, tn);
183 :
184 : if(gotOne){
185 : if(distance < bestDist){
186 : bestDist = distance;
187 : vOut = vTmp;
188 : }
189 : }
190 : else{
191 : gotOne = true;
192 : vOut = vTmp;
193 : bestDist = distance;
194 : }
195 : }
196 :
197 : return gotOne;
198 : }
199 :
200 : template <class TAAPosVRT>
201 0 : int PointFaceTest(vector3& v, Face* f, TAAPosVRT& aaPos)
202 : {
203 : vector3 n;
204 : vector3 dir;
205 :
206 0 : CalculateNormal(n, f, aaPos);
207 0 : VecSubtract(dir, v, aaPos[f->vertex(0)]);
208 :
209 : number d = VecDot(dir, n);
210 0 : if(d > 0)
211 : return 1;
212 0 : if(d < 0)
213 0 : return -1;
214 : return 0;
215 : }
216 :
217 : template <class TAAPosVRT>
218 0 : bool IsDegenerated(Face* f, TAAPosVRT& aaPos, number threshold)
219 : {
220 0 : number threshSQ = threshold * threshold;
221 0 : size_t numVrts = f->num_vertices();
222 :
223 0 : for(size_t i = 0; i < numVrts; ++i){
224 0 : if(VecDistanceSq(aaPos[f->vertex(i)], aaPos[f->vertex((i+1) % numVrts)])
225 : < threshSQ)
226 : {
227 : return true;
228 : }
229 : }
230 : return false;
231 : }
232 :
233 : ////////////////////////////////////////////////////////////////////////
234 : template <class TAAPosVRT>
235 0 : number FaceArea(FaceVertices* f, TAAPosVRT& aaPos)
236 : {
237 : number area = 0;
238 0 : for(size_t i = 2; i < f->num_vertices(); ++i){
239 0 : area += TriangleArea(aaPos[f->vertex(0)],
240 0 : aaPos[f->vertex(i - 1)],
241 0 : aaPos[f->vertex(i)]);
242 : }
243 0 : return area;
244 : }
245 :
246 : ////////////////////////////////////////////////////////////////////////
247 : template <class TIterator, class TAAPosVRT>
248 : number FaceArea(TIterator facesBegin, TIterator facesEnd, TAAPosVRT& aaPos)
249 : {
250 : number sum = 0.;
251 :
252 : for (; facesBegin != facesEnd; ++facesBegin)
253 : sum += FaceArea(*facesBegin, aaPos);
254 :
255 : return sum;
256 : }
257 :
258 : ////////////////////////////////////////////////////////////////////////
259 : template <class TIterator, class TAAPosVRT>
260 : Face* FindSmallestFace(TIterator facesBegin, TIterator facesEnd, TAAPosVRT& aaPos)
261 : {
262 : // if facesBegin equals facesEnd, then the list is empty and we can
263 : // immediately return NULL
264 : if(facesBegin == facesEnd)
265 : return NULL;
266 :
267 : // the first face is the first candidate for the smallest face.
268 : Face* smallestFace = *facesBegin;
269 : number smallestArea = FaceArea(smallestFace, aaPos);
270 : ++facesBegin;
271 :
272 : for(; facesBegin != facesEnd; ++facesBegin){
273 : Face* curFace = *facesBegin;
274 : number curArea = FaceArea(curFace, aaPos);
275 : if(curArea < smallestArea){
276 : smallestFace = curFace;
277 : smallestArea = curArea;
278 : }
279 : }
280 :
281 : return smallestFace;
282 : }
283 :
284 :
285 : }// end of namespace
286 :
287 : #endif
|
