Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Authors: Sebastian Reiter, Martin Stepniewski
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 "volume_util.h"
34 : #include "lib_grid/lg_base.h"
35 : #include "edge_util.h"
36 : #include "lib_grid/grid_objects/tetrahedron_rules.h"
37 :
38 : using namespace std;
39 :
40 : namespace ug
41 : {
42 :
43 : ////////////////////////////////////////////////////////////////////////
44 : // GetNeighbours - sreiter
45 0 : void GetNeighbours(std::vector<Volume*>& vVolsOut, Grid& grid, Volume* v,
46 : int side, bool clearContainer)
47 : {
48 0 : if(clearContainer)
49 : vVolsOut.clear();
50 :
51 : // if VOLOPT_AUTOGENERATE_FACES and FACEOPT_STORE_ASSOCIATED_VOLUMES are
52 : // activated, we may use them to find the connected volume quite fast.
53 0 : if(grid.option_is_enabled(VOLOPT_AUTOGENERATE_FACES
54 : | FACEOPT_STORE_ASSOCIATED_VOLUMES))
55 : {
56 0 : Face* f = grid.get_face(v, side);
57 0 : Grid::AssociatedVolumeIterator iterEnd = grid.associated_volumes_end(f);
58 0 : for(Grid::AssociatedVolumeIterator iter = grid.associated_volumes_begin(f);
59 0 : iter != iterEnd; ++iter)
60 : {
61 0 : if(*iter != v)
62 0 : vVolsOut.push_back(*iter);
63 : }
64 :
65 : return;
66 : }
67 :
68 : // we can't assume that associated faces exist.
69 : // we have to find the neighbour by hand.
70 : // mark all vertices of the side
71 0 : grid.begin_marking();
72 :
73 0 : FaceDescriptor fd;
74 0 : v->face_desc(side, fd);
75 : uint numFaceVrts = fd.num_vertices();
76 0 : for(uint i = 0; i < numFaceVrts; ++ i)
77 0 : grid.mark(fd.vertex(i));
78 :
79 : // iterate over associated volumes of the first vertex and count
80 : // the number of marked vertices it contains.
81 : Vertex* vrt = fd.vertex(0);
82 0 : Grid::AssociatedVolumeIterator iterEnd = grid.associated_volumes_end(vrt);
83 0 : for(Grid::AssociatedVolumeIterator iter = grid.associated_volumes_begin(vrt);
84 0 : iter != iterEnd; ++iter)
85 : {
86 0 : Volume* vol = *iter;
87 0 : if(vol != v){
88 : size_t count = 0;
89 0 : uint numVrts = vol->num_vertices();
90 0 : for(uint i = 0; i < numVrts; ++i){
91 0 : if(grid.is_marked(vol->vertex(i)))
92 0 : ++count;
93 : }
94 :
95 : // if the number of marked vertices in vol matches the
96 : // number of vertices of the specified side, we consider
97 : // the volume to be a neighbout of that side.
98 0 : if(count == numFaceVrts)
99 0 : vVolsOut.push_back(vol);
100 : }
101 : }
102 :
103 0 : grid.end_marking();
104 : }
105 :
106 :
107 : ////////////////////////////////////////////////////////////////////////////////////////////
108 : // CalculateVolumeMinHeight - mstepnie
109 0 : number CalculateMinVolumeHeight(Tetrahedron* tet,
110 : Grid::VertexAttachmentAccessor<APosition>& aaPos)
111 : {
112 : vector3 vfaceNormal;
113 : number minHeight = std::numeric_limits<double>::max();
114 : number tmpMinHeight;
115 :
116 :
117 : // Iterate over all tetrahedron vertices and calculate height to opposing face
118 0 : for(size_t i = 0; i < 4; ++i)
119 : {
120 0 : Vertex* vrt = tet->vertex(i);
121 0 : FaceDescriptor fd = tet->face_desc(tet->get_opposing_object(vrt).second);
122 :
123 0 : CalculateNormal(vfaceNormal, &fd, aaPos);
124 0 : tmpMinHeight = DistancePointToPlane(aaPos[vrt], aaPos[fd.vertex(0)], vfaceNormal);
125 :
126 0 : if(tmpMinHeight < minHeight)
127 : {
128 : minHeight = tmpMinHeight;
129 : }
130 : }
131 :
132 0 : return minHeight;
133 : }
134 :
135 :
136 : ////////////////////////////////////////////////////////////////////////////////////////////
137 : // CalculateTetrahedronAspectRatio - mstepnie
138 0 : number CalculateTetrahedronAspectRatio(Grid& grid, Tetrahedron* tet,
139 : Grid::VertexAttachmentAccessor<AVector3>& aaPos)
140 : {
141 : /*
142 : * optimal Aspect Ratio of a regular tetrahedron with edge lengths a:
143 : * Q = hmin/lmax = sqrt(2/3)*a / a = 0.8164...
144 : *
145 : * Info: return value is normalized by factor sqrt(3/2)
146 : * (s. Shewchuk 2002)
147 : */
148 :
149 : number aspectRatio;
150 : number maxEdgeLength;
151 : number minTetrahedronHeight;
152 :
153 : // Collect tetrahedron edges, find longest edge and calculate its length
154 : vector<Edge*> edges;
155 : CollectAssociated(edges, grid, tet);
156 0 : Edge* longestEdge = FindLongestEdge(edges.begin(), edges.end(), aaPos);
157 0 : maxEdgeLength = EdgeLength(longestEdge, aaPos);
158 :
159 : // Calculate the minimal tetrahedron height
160 0 : minTetrahedronHeight = CalculateMinVolumeHeight(tet, aaPos);
161 :
162 : // Calculate the aspect ratio
163 0 : aspectRatio = std::sqrt(3/2.0) * minTetrahedronHeight / maxEdgeLength;
164 :
165 0 : return aspectRatio;
166 0 : }
167 :
168 : ////////////////////////////////////////////////////////////////////////////////////////////
169 : // CalculateHexahedronAspectRatio
170 : ////////////////////////////////////////////////////////////////////////////////////////////
171 0 : number CalculateHexahedronAspectRatio(Grid& grid, Hexahedron* hex,
172 : Grid::VertexAttachmentAccessor<AVector3>& aaPos)
173 : {
174 : /*
175 : * Another important indicator of the mesh quality is the aspect ratio.
176 : * The aspect ratio is a measure of the stretching of a cell. It is computed
177 : * as the ratio of the maximum value to the minimum value of any of the
178 : * following distances: the normal distances between the cell centroid and
179 : * face centroids (computed as a dot product of the distance vector and the
180 : * face normal), and the distances between the cell centroid and nodes.
181 : * For a unit cube (see Figure 5.23: Calculating the Aspect Ratio for a Unit
182 : * Cube), the maximum distance is 0.866, and the minimum distance is 0.5,
183 : * so the aspect ratio is 1.732. This type of definition can be applied on
184 : * any type of mesh, including polyhedral.
185 : *
186 : * \see https://www.sharcnet.ca/Software/Ansys/17.0/en-us/help/flu_ug/flu_ug_mesh_quality.html
187 : *
188 : */
189 : vector<Face*> faces;
190 : CollectAssociated(faces, grid, hex);
191 0 : ug::vector3 center = CalculateCenter(hex, aaPos);
192 :
193 : /// Max distance to face center
194 : std::vector<number> As;
195 0 : for (size_t i = 0; i < faces.size(); i++)
196 : {
197 : vector3 dir;
198 0 : vector3 faceCenter = CalculateCenter(faces[i], aaPos);
199 : VecSubtract(dir, faceCenter, center);
200 : vector3 n;
201 0 : CalculateNormal(n, faces[i], aaPos);
202 0 : number dist = VecDot(dir, n);
203 0 : As.push_back(dist);
204 : }
205 0 : number max = *std::max_element(As.begin(), As.end());
206 :
207 : /// Min distance to node of hexaeder
208 : std::vector<number> Bs;
209 0 : for (size_t i = 0; i < hex->num_vertices(); i++)
210 : {
211 0 : number dist = VecDistance(aaPos[hex->vertex(i)], center);
212 0 : Bs.push_back(dist);
213 : }
214 0 : number min = *std::min_element(Bs.begin(), Bs.end());
215 :
216 : UG_COND_WARNING(std::abs(min) < SMALL, "Near 0-length minimum distance detected.");
217 :
218 0 : return max / min;
219 0 : }
220 :
221 : ////////////////////////////////////////////////////////////////////////////////////////////
222 : // CalculateHexahedronAspectRatio
223 : // order of vertices should be the same as described in \sa PyramidDescriptor
224 : // v1, v2, v3, v4: bottom-vertices in counterclockwise order (if viewed from the top).
225 : // v5: top-vertex.
226 : ////////////////////////////////////////////////////////////////////////////////////////////
227 0 : number CalculatePyramidAspectRatio
228 : (
229 : Grid& grid,
230 : Pyramid* pyr,
231 : Grid::VertexAttachmentAccessor<AVector3>& aaPos
232 : )
233 : {
234 : // average edge length of base of pyramid
235 : number avg_edge_length = 0;
236 0 : for (size_t i = 0; i < 4; i++)
237 : {
238 0 : avg_edge_length += VecDistance(aaPos[pyr->vertex(i%4)], aaPos[pyr->vertex((i+1)%4)]);
239 : }
240 0 : avg_edge_length /= 4;
241 :
242 : // distance from base to top of pyramid
243 0 : const Face* const face = grid.get_element(FaceDescriptor(pyr->vertex(0), pyr->vertex(1), pyr->vertex(2), pyr->vertex(3)));
244 : vector3 normal;
245 0 : CalculateNormal(normal, face, aaPos);
246 0 : const number distance = VecDot(normal, aaPos[pyr->vertex(4)]);
247 :
248 : // AR of pyramid is ratio between distance and average edge length of base
249 0 : return distance / avg_edge_length;
250 : }
251 :
252 :
253 : ////////////////////////////////////////////////////////////////////////////////////////////
254 : // CalculateTetrahedronRootMeanSquareFaceArea - mstepnie
255 : ////////////////////////////////////////////////////////////////////////////////////////////
256 : UG_API
257 0 : number CalculateTetrahedronRootMeanSquareFaceArea(Grid& grid,
258 : Tetrahedron* tet,
259 : Grid::VertexAttachmentAccessor<APosition>& aaPos)
260 : {
261 : // Collect tetrahedron faces
262 : vector<Face*> faces;
263 : CollectAssociated(faces, grid, tet);
264 :
265 : number A;
266 : number A_rms = 0.0;
267 :
268 0 : for(size_t i = 0; i < faces.size(); ++i)
269 : {
270 0 : A = FaceArea(faces[i], aaPos);
271 0 : A_rms += A*A;
272 : }
273 :
274 0 : A_rms *= 0.25;
275 0 : A_rms = std::sqrt(A_rms);
276 :
277 0 : return A_rms;
278 0 : }
279 :
280 :
281 : ////////////////////////////////////////////////////////////////////////////////////////////
282 : // CalculateTetrahedronVolToRMSFaceAreaRatio - mstepnie
283 : ////////////////////////////////////////////////////////////////////////////////////////////
284 : UG_API
285 0 : number CalculateTetrahedronVolToRMSFaceAreaRatio(Grid& grid,
286 : Tetrahedron* tet,
287 : Grid::VertexAttachmentAccessor<APosition>& aaPos)
288 : {
289 : /*
290 : * optimal volume to root-mean-square face area ratio of a
291 : * regular tetrahedron with edge lengths a:
292 : * Q = V/A_rms^(3/2)
293 : *
294 : * Info: return value is normalized by factor pow(3, 7/4.0) / 2.0 / sqrt(2);
295 : * (s. Shewchuk 2002)
296 : */
297 :
298 0 : number vol = CalculateTetrahedronVolume(aaPos[tet->vertex(0)],
299 0 : aaPos[tet->vertex(1)],
300 0 : aaPos[tet->vertex(2)],
301 0 : aaPos[tet->vertex(3)]);
302 :
303 0 : number A_rms = CalculateTetrahedronRootMeanSquareFaceArea(grid, tet, aaPos);
304 : number normalization = std::pow(3, 7/4.0) / 2.0 / std::sqrt(2);
305 :
306 0 : return normalization * vol / std::pow(A_rms, 3/2.0);
307 : }
308 :
309 : ////////////////////////////////////////////////////////////////////////////////////////////
310 : // CalculateHexahedronVolToRMSFaceAreaRatio - sgrein
311 : ////////////////////////////////////////////////////////////////////////////////////////////
312 : UG_API
313 0 : number CalculateHexahedronVolToRMSFaceAreaRatio(Grid& grid,
314 : Hexahedron* hex,
315 : Grid::VertexAttachmentAccessor<APosition>& aaPos)
316 : {
317 0 : UG_THROW("ERROR in CalculateHexahedronVolToRMSFaceAreaRatio: Not yet implemented");
318 : return NAN;
319 : }
320 :
321 :
322 : ////////////////////////////////////////////////////////////////////////
323 : // IntersectPlaneWithTetrahedron
324 0 : size_t IntersectPlaneWithTetrahedron
325 : (
326 : vector3 intsOut[4],
327 : const vector3& planePoint,
328 : const vector3& planeNormal,
329 : const vector3 t[4]
330 : )
331 : {
332 : using namespace tet_rules;
333 : size_t numIntersections = 0;
334 : int intersectingEdges[4];
335 0 : for(size_t ie = 0; ie < (size_t) NUM_EDGES; ++ie)
336 : {
337 0 : vector3 v0 = t[EDGE_VRT_INDS[ie][0]];
338 0 : vector3 v1 = t[EDGE_VRT_INDS[ie][1]];
339 : vector3 dir;
340 : VecSubtract(dir, v1, v0);
341 :
342 : vector3 p;
343 : number s;
344 0 : if(RayPlaneIntersection(p, s, v0, dir, planePoint, planeNormal))
345 : {
346 0 : if(s > 0 && s < 1)
347 : {
348 0 : if (numIntersections >= 4) // to avoid numerical artefacts
349 0 : UG_THROW ("IntersectPlaneWithTetrahedron:"
350 : " Illegal number of intersections of a plane with a tetrahedron.");
351 0 : intsOut[numIntersections] = p;
352 0 : intersectingEdges[numIntersections] = ie;
353 0 : ++numIntersections;
354 : }
355 : }
356 : }
357 :
358 0 : if(numIntersections == 4)
359 : {
360 0 : if(intersectingEdges[1] == OPPOSED_EDGE[intersectingEdges[0]])
361 : swap(intsOut[1], intsOut[2]);
362 0 : else if(intersectingEdges[2] == OPPOSED_EDGE[intersectingEdges[1]])
363 : swap(intsOut[0], intsOut[1]);
364 : }
365 :
366 0 : return numIntersections;
367 : }
368 :
369 0 : void InsertCenterVertex(Grid& g, Volume* vol, Vertex* vrt, bool eraseOldVol)
370 : {
371 : // get the sides of the volume and create new elements
372 0 : FaceDescriptor fd;
373 0 : for(size_t i = 0; i < vol->num_faces(); ++i){
374 0 : vol->face_desc(i, fd);
375 0 : if(fd.num_vertices() == 3){
376 : // create a tetrahedron
377 0 : g.create<Tetrahedron>(TetrahedronDescriptor(fd.vertex(2), fd.vertex(1),
378 : fd.vertex(0), vrt), vol);
379 : }
380 0 : else if(fd.num_vertices() == 4){
381 : // create a pyramid
382 0 : g.create<Pyramid>(PyramidDescriptor(fd.vertex(3), fd.vertex(2),
383 : fd.vertex(1), fd.vertex(0), vrt), vol);
384 : }
385 : else{
386 0 : UG_THROW("Unsupported face type in InsertCenterVertex (#Corners "
387 : << fd.num_vertices() << ")");
388 : }
389 : }
390 :
391 0 : if(eraseOldVol)
392 0 : g.erase(vol);
393 0 : }
394 :
395 : /*
396 : double CMesh::calculate_volume_gauss() {
397 :
398 : double volume = 0.0;
399 :
400 : for(uint i = 0; i < number_of_triangles; i++) {
401 :
402 : uint a = triangles[i].a;
403 : uint b = triangles[i].b;
404 : uint c = triangles[i].c;
405 :
406 : //printf("%f %f %f\n", vertices[b].x(), vertices[b].y(), vertices[b].z());
407 :
408 : double x = (vertices[b].y() - vertices[a].y()) * (vertices[c].z() - vertices[a].z()) - (vertices[b].z() - vertices[a].z()) * (vertices[c].y() - vertices[a].y());
409 : double y = (vertices[b].z() - vertices[a].z()) * (vertices[c].x() - vertices[a].x()) - (vertices[b].x() - vertices[a].x()) * (vertices[c].z() - vertices[a].z());
410 : double z = (vertices[b].x() - vertices[a].x()) * (vertices[c].y() - vertices[a].y()) - (vertices[b].y() - vertices[a].y()) * (vertices[c].x() - vertices[a].x());
411 :
412 : double length = sqrt(x * x + y * y + z * z);
413 :
414 : double surface = 0.5 * length;
415 :
416 : if(length > 0.0) {
417 :
418 : x /= length;
419 : y /= length;
420 : z /= length;
421 :
422 : double sx = (vertices[a].x() + vertices[b].x() + vertices[c].x()) / 3.0;
423 : double sy = (vertices[a].y() + vertices[b].y() + vertices[c].y()) / 3.0;
424 : double sz = (vertices[a].z() + vertices[b].z() + vertices[c].z()) / 3.0;
425 :
426 : volume += 1.0 / 3.0 * surface * (sx * x + sy * y + sz * z);
427 :
428 : }
429 :
430 :
431 : }
432 :
433 : return volume;
434 :
435 :
436 :
437 : }
438 : */
439 : }// end of namespace
440 :
441 :
|
