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 : #include <vector>
34 : #include <stack>
35 : #include <cassert>
36 : #include "vertex_util.h"
37 : #include "face_util.h"
38 : #include "../attachment_util.h"
39 :
40 : using namespace std;
41 :
42 : namespace ug
43 : {
44 :
45 : ////////////////////////////////////////////////////////////////////////
46 0 : int GetFaceIndex(Volume* vol, Face* f)
47 : {
48 0 : size_t numFaces = vol->num_faces();
49 0 : FaceDescriptor fd;
50 0 : for(uint i = 0; i < numFaces; ++i)
51 : {
52 0 : vol->face_desc(i, fd);
53 0 : if(CompareVertices(f, &fd))
54 0 : return (int)i;
55 : }
56 : return -1;
57 : }
58 :
59 :
60 : ////////////////////////////////////////////////////////////////////////
61 : // CalculateNormal
62 0 : void CalculateNormal(vector3& vNormOut, const FaceVertices* face,
63 : Grid::AttachmentAccessor<Vertex, APosition>& aaPos)
64 : {
65 0 : if(face->num_vertices() == 3)
66 : {
67 0 : CalculateTriangleNormal(vNormOut, aaPos[face->vertex(0)],
68 0 : aaPos[face->vertex(1)], aaPos[face->vertex(2)]);
69 0 : return;
70 : }
71 0 : else if(face->num_vertices() == 4)
72 : {
73 : vector3 n1, n2;
74 0 : CalculateTriangleNormalNoNormalize(n1, aaPos[face->vertex(0)],
75 0 : aaPos[face->vertex(1)], aaPos[face->vertex(2)]);
76 0 : CalculateTriangleNormalNoNormalize(n2, aaPos[face->vertex(2)],
77 0 : aaPos[face->vertex(3)], aaPos[face->vertex(0)]);
78 : VecAdd(vNormOut, n1, n2);
79 0 : VecNormalize(vNormOut, vNormOut);
80 :
81 : return;
82 : }
83 0 : else if(face->num_vertices() > 4)
84 : {
85 0 : CalculateTriangleNormal(vNormOut, aaPos[face->vertex(0)],
86 0 : aaPos[face->vertex(1)], aaPos[face->vertex(2)]);
87 0 : return;
88 : }
89 :
90 : vNormOut = vector3(0, 0, 0);
91 0 : return;
92 : }
93 :
94 0 : void CalculateNormalNoNormalize(vector3& vNormOut, FaceVertices* face,
95 : Grid::AttachmentAccessor<Vertex, APosition>& aaPos)
96 : {
97 0 : if(face->num_vertices() == 3)
98 : {
99 0 : CalculateTriangleNormalNoNormalize(vNormOut, aaPos[face->vertex(0)],
100 0 : aaPos[face->vertex(1)], aaPos[face->vertex(2)]);
101 0 : return;
102 : }
103 0 : else if(face->num_vertices() == 4)
104 : {
105 : vector3 n1, n2;
106 0 : CalculateTriangleNormalNoNormalize(n1, aaPos[face->vertex(0)],
107 0 : aaPos[face->vertex(1)], aaPos[face->vertex(2)]);
108 0 : CalculateTriangleNormalNoNormalize(n2, aaPos[face->vertex(2)],
109 0 : aaPos[face->vertex(3)], aaPos[face->vertex(0)]);
110 : VecAdd(vNormOut, n1, n2);
111 : VecScale(vNormOut, vNormOut, 0.5);
112 :
113 : return;
114 : }
115 0 : else if(face->num_vertices() > 4)
116 : {
117 0 : CalculateTriangleNormalNoNormalize(vNormOut, aaPos[face->vertex(0)],
118 0 : aaPos[face->vertex(1)], aaPos[face->vertex(2)]);
119 0 : return;
120 : }
121 :
122 : vNormOut = vector3(0, 0, 0);
123 0 : return;
124 : }
125 :
126 : ////////////////////////////////////////////////////////////////////////
127 : // CalculateFaceNormals
128 0 : void CalculateFaceNormals(Grid& grid, const FaceIterator& facesBegin,
129 : const FaceIterator& facesEnd,
130 : AVector3& aPos, AVector3& aNorm)
131 : {
132 0 : if(!grid.has_vertex_attachment(aPos))
133 0 : return;
134 :
135 0 : if(!grid.has_face_attachment(aNorm))
136 0 : grid.attach_to_faces(aNorm);
137 :
138 : Grid::VertexAttachmentAccessor<AVector3> aaPos(grid, aPos);
139 : Grid::FaceAttachmentAccessor<AVector3> aaNorm(grid, aNorm);
140 :
141 : // iterate through the specified faces and calculate normals.
142 0 : for(FaceIterator iter = facesBegin; iter != facesEnd; ++iter)
143 0 : CalculateNormal(aaNorm[*iter], *iter, aaPos);
144 : }
145 :
146 : ////////////////////////////////////////////////////////////////////////
147 0 : int NumAssociatedVolumes(Grid& grid, Face* f)
148 : {
149 : // check if FACEOPT_STORE_ASSOCIATED_VOLUMES is enabled.
150 : // if so, use it to count the number of adjacent volumes.
151 : int counter = 0;
152 0 : if(grid.option_is_enabled(FACEOPT_STORE_ASSOCIATED_VOLUMES))
153 : {
154 0 : for(Grid::AssociatedVolumeIterator iter = grid.associated_volumes_begin(f);
155 0 : iter != grid.associated_volumes_end(f); iter++)
156 : {
157 0 : ++counter;
158 : }
159 : }
160 : else
161 : {
162 : // iterate over all volumes which are connected to the first vertex
163 : // and check if they contain the face...
164 0 : Grid::AssociatedVolumeIterator iterEnd = grid.associated_volumes_end(f->vertex(0));
165 0 : for(Grid::AssociatedVolumeIterator iter = grid.associated_volumes_begin(f->vertex(0));
166 0 : iter != iterEnd; iter++)
167 : {
168 0 : if(VolumeContains(*iter, f))
169 0 : ++counter;
170 : }
171 : }
172 :
173 0 : return counter;
174 : }
175 :
176 : ////////////////////////////////////////////////////////////////////////
177 0 : bool IsVolumeBoundaryFace(Grid& grid, Face* f)
178 : {
179 0 : return NumAssociatedVolumes(grid, f) == 1;
180 : }
181 :
182 0 : bool IsVolumeBoundaryFace(Grid& grid, ConstrainedFace* f)
183 : {
184 0 : static vector<Volume*> vols;//avoid repeated reallocation
185 0 : CollectVolumes(vols, grid, f);
186 0 : if(vols.empty())
187 0 : return true;
188 : return false;
189 : }
190 :
191 0 : bool IsVolumeBoundaryFace(Grid& grid, ConstrainingFace* f)
192 : {
193 0 : static vector<Volume*> vols;//avoid repeated reallocation
194 0 : CollectVolumes(vols, grid, f);
195 0 : if(vols.empty())
196 0 : return true;
197 : return false;
198 : }
199 :
200 :
201 : ////////////////////////////////////////////////////////////////////////
202 : // FaceQuality
203 0 : number FaceQuality(FaceVertices* f,
204 : Grid::VertexAttachmentAccessor<APosition> aaPos)
205 : {
206 : // take dot-products at the corners.
207 : number quality = 1;
208 0 : uint numVrts = f->num_vertices();
209 :
210 : vector3 v1, v2;
211 : number len1, len2;
212 :
213 : // calculate the direction from first to last vertex.
214 0 : VecSubtract(v1, aaPos[f->vertex(numVrts-1)], aaPos[f->vertex(0)]);
215 : len1 = VecLength(v1);
216 :
217 : //VecNormalize(v1, v1);
218 :
219 0 : for(uint i = 0; i < numVrts; ++i)
220 : {
221 0 : VecSubtract(v2, aaPos[f->vertex((i+1)%numVrts)], aaPos[f->vertex(i)]);
222 : len2 = VecLength(v2);
223 : //VecNormalize(v2, v2);
224 :
225 0 : number nQual = 1. - fabs(VecDot(v1, v2) / (len1 * len2));
226 0 : if(nQual < quality)
227 : quality = nQual;
228 : // v1 of the next iteration equals -v2 of this iteration
229 : VecScale(v1, v2, -1);
230 : len1 = len2;
231 : }
232 :
233 0 : if(numVrts == 3){
234 : // since at least one angle is <= 60, we have to normalize the return value
235 0 : return quality * 2.;
236 : }
237 : return quality;
238 : }
239 :
240 : ////////////////////////////////////////////////////////////////////////
241 : // TriangleQuality
242 0 : number TriangleQuality(vector3& v1, vector3& v2, vector3& v3)
243 : {
244 : // take dot-products at the corners.
245 : number quality = 1;
246 :
247 : // required for the iteration.
248 : vector3* pv[3];
249 0 : pv[0] = &v1; pv[1] = &v2; pv[2] = &v3;
250 :
251 : vector3 d1, d2;
252 : // calculate the direction from first to last vertex.
253 : VecSubtract(d1, v3, v1);
254 0 : VecNormalize(d1, d1);
255 :
256 0 : for(uint i = 0; i < 3; ++i)
257 : {
258 0 : VecSubtract(d2, *pv[i], *pv[(i+1)%3]);
259 0 : VecNormalize(d2, d2);
260 0 : number nQual = 1.f - fabs(VecDot(d1, d2));
261 0 : if(nQual < quality)
262 : quality = nQual;
263 : // v1 of the next iteration equals -v2 of this iteration
264 : VecScale(d1, d2, -1);
265 : }
266 :
267 : // since at least one angle is <= 60, we have to normalize the return value
268 0 : return quality * 2.;
269 : }
270 :
271 : ////////////////////////////////////////////////////////////////////////
272 : // Triangulate
273 0 : void Triangulate(Grid& grid, Quadrilateral* q,
274 : Grid::VertexAttachmentAccessor<APosition>* paaPos)
275 : {
276 : // there are two ways to triangulate a quadrilateral.
277 : // if quality is set to false, we simply choose the first way.
278 0 : if(!paaPos)
279 : {
280 0 : grid.create<Triangle>(TriangleDescriptor(q->vertex(0), q->vertex(1), q->vertex(2)), q);
281 0 : grid.create<Triangle>(TriangleDescriptor(q->vertex(2), q->vertex(3), q->vertex(0)), q);
282 : }
283 : else
284 : {
285 : Grid::VertexAttachmentAccessor<APosition>& aaPos = *paaPos;
286 :
287 : number q1 = 0;
288 : number q2 = 0;
289 :
290 : // first check whether normals in the corners should affect the direction of the edge
291 0 : vector3 n[4];
292 0 : for(int i = 0; i < 4; ++i)
293 0 : CalculateVertexNormal(n[i], grid, q->vertex(i), aaPos);
294 :
295 : // vector3 dir[2];
296 : // VecSubtract(dir[0], aaPos[q->vertex(2)], aaPos[q->vertex(0)]);
297 : // VecNormalize(dir[0], dir[0]);
298 : // VecSubtract(dir[1], aaPos[q->vertex(3)], aaPos[q->vertex(1)]);
299 : // VecNormalize(dir[1], dir[1]);
300 :
301 : // number maxDot[2] = {0, 0};
302 : // for(int i = 0; i < 4; ++i){
303 : // for(int j = 0; j < 2; ++j){
304 : // maxDot[j] = std::max(maxDot[j], fabs(VecDot(dir[j], n[i])));
305 : // }
306 : // }
307 :
308 : // int imin = maxDot[0] < maxDot[1] ? 0 : 1;
309 : // int imax = maxDot[0] < maxDot[1] ? 1 : 0;
310 :
311 : // if(maxDot[imax] > 0 && maxDot[imax] - maxDot[imin] > SMALL){
312 : // q1 = 1. - maxDot[0];
313 : // q2 = 1. - maxDot[1];
314 : // }
315 :
316 : q1 = VecDot(n[0], n[2]);
317 : q2 = VecDot(n[1], n[3]);
318 :
319 0 : if(fabs(q1-q2) < SMALL)
320 : {
321 : // if normals are equal we check which way produces the better triangles.
322 0 : q1 = std::min(
323 0 : TriangleQuality(aaPos[q->vertex(0)], aaPos[q->vertex(1)], aaPos[q->vertex(2)]),
324 0 : TriangleQuality(aaPos[q->vertex(2)], aaPos[q->vertex(3)], aaPos[q->vertex(0)]));
325 :
326 0 : q2 = std::min(
327 0 : TriangleQuality(aaPos[q->vertex(1)], aaPos[q->vertex(2)], aaPos[q->vertex(3)]),
328 0 : TriangleQuality(aaPos[q->vertex(3)], aaPos[q->vertex(0)], aaPos[q->vertex(1)]));
329 : }
330 :
331 0 : if(q1 >= q2)
332 : {
333 0 : grid.create<Triangle>(TriangleDescriptor(q->vertex(0), q->vertex(1), q->vertex(2)), q);
334 0 : grid.create<Triangle>(TriangleDescriptor(q->vertex(2), q->vertex(3), q->vertex(0)), q);
335 : }
336 : else
337 : {
338 0 : grid.create<Triangle>(TriangleDescriptor(q->vertex(1), q->vertex(2), q->vertex(3)), q);
339 0 : grid.create<Triangle>(TriangleDescriptor(q->vertex(3), q->vertex(0), q->vertex(1)), q);
340 : }
341 : }
342 :
343 0 : grid.erase(q);
344 0 : }
345 :
346 : ////////////////////////////////////////////////////////////////////////
347 : // Triangulate
348 0 : void Triangulate(Grid& grid,
349 : QuadrilateralIterator iterBegin,
350 : QuadrilateralIterator iterEnd,
351 : Grid::VertexAttachmentAccessor<APosition>* paaPos)
352 : {
353 0 : while(iterBegin != iterEnd)
354 : {
355 : Quadrilateral* q = *iterBegin;
356 : iterBegin++;
357 0 : Triangulate(grid, q, paaPos);
358 : }
359 0 : }
360 :
361 : ////////////////////////////////////////////////////////////////////////
362 : // GetNeighbours
363 0 : void GetNeighbours(std::vector<Face*>& vFacesOut, Grid& grid, Face* f,
364 : int side, bool clearContainer)
365 : {
366 : // in the current implementation this method requires, that all edges
367 : // are created for all faces.
368 : //TODO: improve this!
369 0 : if(!grid.option_is_enabled(FACEOPT_AUTOGENERATE_EDGES))
370 : {
371 : LOG("WARNING: autoenabling FACEOPT_AUTOGENERATE_EDGES in GetNeighbours(Face).\n");
372 0 : grid.enable_options(FACEOPT_AUTOGENERATE_EDGES);
373 : }
374 :
375 0 : if(clearContainer)
376 : vFacesOut.clear();
377 :
378 : // check if we can find an edge for the specified side
379 0 : Edge* e = grid.get_edge(f, side);
380 : assert(e && "edge not found though it should be there!");
381 :
382 : // get the connected faces
383 : vector<Face*> vFaces;
384 0 : CollectFaces(vFaces, grid, e);
385 :
386 : // push them to vFacesOut - except f
387 0 : for(size_t i = 0; i < vFaces.size(); ++i)
388 : {
389 0 : if(vFaces[i] != f)
390 0 : vFacesOut.push_back(vFaces[i]);
391 : }
392 0 : }
393 :
394 0 : void InsertCenterVertex(Grid& g, Face* f, Vertex* vrt, bool eraseOldFace)
395 : {
396 : // get the sides of the face and create new elements
397 0 : EdgeDescriptor ed;
398 0 : for(size_t i = 0; i < f->num_edges(); ++i){
399 0 : f->edge_desc(i, ed);
400 0 : g.create<Triangle>(TriangleDescriptor(ed.vertex(0), ed.vertex(1), vrt), f);
401 : }
402 :
403 0 : if(eraseOldFace)
404 0 : g.erase(f);
405 0 : }
406 :
407 : }// end of namespace
|
