Line data Source code
1 : /*
2 : * Copyright (c) 2009-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Markus Breit
4 : * Date: 2018-05-25
5 : *
6 : * This file is part of UG4.
7 : *
8 : * UG4 is free software: you can redistribute it and/or modify it under the
9 : * terms of the GNU Lesser General Public License version 3 (as published by the
10 : * Free Software Foundation) with the following additional attribution
11 : * requirements (according to LGPL/GPL v3 §7):
12 : *
13 : * (1) The following notice must be displayed in the Appropriate Legal Notices
14 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
15 : *
16 : * (2) The following notice must be displayed at a prominent place in the
17 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
18 : *
19 : * (3) The following bibliography is recommended for citation and must be
20 : * preserved in all covered files:
21 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
22 : * parallel geometric multigrid solver on hierarchically distributed grids.
23 : * Computing and visualization in science 16, 4 (2013), 151-164"
24 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
25 : * flexible software system for simulating pde based models on high performance
26 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
27 : *
28 : * This program is distributed in the hope that it will be useful,
29 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
30 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 : * GNU Lesser General Public License for more details.
32 : */
33 :
34 : #include "lib_grid/algorithms/element_side_util.h" // for GetOpposingSide
35 : #include "lib_grid/grid/grid_util.h" // for CompareVertices
36 :
37 :
38 : namespace ug {
39 :
40 :
41 : template <typename TAAPos>
42 : AnisotropyState is_anisotropic
43 : (
44 : Edge* elem,
45 : const TAAPos& aaPos,
46 : number thresholdRatio
47 : )
48 : {
49 : return ISOTROPIC;
50 : }
51 :
52 :
53 :
54 : template <typename TAAPos>
55 0 : static AnisotropyState is_anisotropic
56 : (
57 : Quadrilateral* q,
58 : const TAAPos& aaPos,
59 : number thresholdRatio
60 : )
61 : {
62 : // check whether elem is anisotropic
63 0 : number sideLength02 = VertexDistance(q->vertex(0), q->vertex(1), aaPos)
64 0 : + VertexDistance(q->vertex(2), q->vertex(3), aaPos);
65 :
66 0 : number sideLength13 = VertexDistance(q->vertex(1), q->vertex(2), aaPos)
67 0 : + VertexDistance(q->vertex(3), q->vertex(0), aaPos);
68 :
69 0 : if (sideLength02 < thresholdRatio * sideLength13)
70 : return QUAD_SHORTX;
71 0 : if (sideLength13 < thresholdRatio * sideLength02)
72 0 : return QUAD_SHORTY;
73 :
74 : return ISOTROPIC;
75 : }
76 :
77 :
78 :
79 : template <typename TAAPos>
80 : AnisotropyState is_anisotropic
81 : (
82 : Face* elem,
83 : const TAAPos& aaPos,
84 : number thresholdRatio
85 : )
86 : {
87 : Quadrilateral* q = dynamic_cast<Quadrilateral*>(elem);
88 : if (!q)
89 : return ISOTROPIC;
90 :
91 : return is_anisotropic(q, aaPos, thresholdRatio);
92 : }
93 :
94 :
95 :
96 :
97 : template <typename TAAPos>
98 0 : static AnisotropyState is_anisotropic
99 : (
100 : Prism* p,
101 : const TAAPos& aaPos,
102 : number thresholdRatio
103 : )
104 : {
105 : // check whether elem is anisotropic
106 0 : number length = VertexDistance(p->vertex(0), p->vertex(3), aaPos)
107 0 : + VertexDistance(p->vertex(1), p->vertex(4), aaPos)
108 0 : + VertexDistance(p->vertex(2), p->vertex(5), aaPos);
109 0 : number width = VertexDistance(p->vertex(0), p->vertex(1), aaPos)
110 0 : + VertexDistance(p->vertex(1), p->vertex(2), aaPos)
111 0 : + VertexDistance(p->vertex(2), p->vertex(0), aaPos)
112 0 : + VertexDistance(p->vertex(3), p->vertex(4), aaPos)
113 0 : + VertexDistance(p->vertex(4), p->vertex(5), aaPos)
114 0 : + VertexDistance(p->vertex(5), p->vertex(3), aaPos);
115 :
116 0 : number ratio = width ? 2.0*length/width : std::numeric_limits<number>::max();
117 :
118 : // flat case
119 0 : if (ratio < thresholdRatio)
120 : return PRISM_FLAT;
121 :
122 : // long case
123 0 : if (ratio*thresholdRatio > 1)
124 0 : return PRISM_LONG;
125 :
126 : return ISOTROPIC;
127 : }
128 :
129 :
130 :
131 : template <typename TAAPos>
132 0 : static AnisotropyState is_anisotropic
133 : (
134 : Hexahedron* hex,
135 : const TAAPos& aaPos,
136 : number thresholdRatio
137 : )
138 : {
139 0 : number length1 = VertexDistance(hex->vertex(0), hex->vertex(1), aaPos)
140 0 : + VertexDistance(hex->vertex(2), hex->vertex(3), aaPos)
141 0 : + VertexDistance(hex->vertex(4), hex->vertex(5), aaPos)
142 0 : + VertexDistance(hex->vertex(6), hex->vertex(7), aaPos);
143 0 : number length2 = VertexDistance(hex->vertex(0), hex->vertex(3), aaPos)
144 0 : + VertexDistance(hex->vertex(1), hex->vertex(2), aaPos)
145 0 : + VertexDistance(hex->vertex(4), hex->vertex(7), aaPos)
146 0 : + VertexDistance(hex->vertex(5), hex->vertex(6), aaPos);
147 0 : number length3 = VertexDistance(hex->vertex(0), hex->vertex(4), aaPos)
148 0 : + VertexDistance(hex->vertex(1), hex->vertex(5), aaPos)
149 0 : + VertexDistance(hex->vertex(2), hex->vertex(6), aaPos)
150 0 : + VertexDistance(hex->vertex(3), hex->vertex(7), aaPos);
151 :
152 : bool shortx = false;
153 : bool shorty = false;
154 : bool shortz = false;
155 :
156 0 : if (length1 < thresholdRatio * length2 || length1 < thresholdRatio * length3)
157 : shortx = true;
158 :
159 0 : if (length2 < thresholdRatio * length1 || length2 < thresholdRatio * length3)
160 : shorty = true;
161 :
162 0 : if (length3 < thresholdRatio * length1 || length3 < thresholdRatio * length2)
163 : shortz = true;
164 :
165 0 : if (shortx)
166 : {
167 0 : if (shorty)
168 : return HEX_SHORTXY;
169 0 : if (shortz)
170 : return HEX_SHORTXZ;
171 0 : return HEX_SHORTX;
172 : }
173 :
174 0 : if (shorty)
175 : {
176 0 : if (shortz)
177 : return HEX_SHORTYZ;
178 0 : return HEX_SHORTY;
179 : }
180 :
181 0 : if (shortz)
182 0 : return HEX_SHORTZ;
183 :
184 : return ISOTROPIC;
185 : }
186 :
187 :
188 :
189 : template <typename TAAPos>
190 : AnisotropyState is_anisotropic
191 : (
192 : Volume* elem,
193 : const TAAPos& aaPos,
194 : number thresholdRatio
195 : )
196 : {
197 : // treat prism case
198 : Prism* prism = dynamic_cast<Prism*>(elem);
199 : if (prism)
200 : return is_anisotropic(prism, aaPos, thresholdRatio);
201 :
202 : // treat hexahedron case
203 : Hexahedron* hex = dynamic_cast<Hexahedron*>(elem);
204 : if (hex)
205 : return is_anisotropic(hex, aaPos, thresholdRatio);
206 :
207 : // other cases do not exist
208 : return ISOTROPIC;
209 : }
210 :
211 :
212 :
213 : template <typename TAAPos>
214 : AnisotropyState close_sides_of_anisotropic_elem
215 : (
216 : Edge* elem,
217 : Grid& grid,
218 : const TAAPos& aaPos,
219 : number thresholdRatio,
220 : std::vector<Vertex*>& sidesOut
221 : )
222 : {
223 : return ISOTROPIC;
224 : }
225 :
226 :
227 : template <typename TAAPos>
228 0 : AnisotropyState close_sides_of_anisotropic_elem
229 : (
230 : Face* elem,
231 : Grid& grid,
232 : const TAAPos& aaPos,
233 : number thresholdRatio,
234 : std::vector<Edge*>& sidesOut
235 : )
236 : {
237 : // check whether this is a quadrilateral
238 0 : Quadrilateral* q = dynamic_cast<Quadrilateral*>(elem);
239 0 : if (!q)
240 : return ISOTROPIC;
241 :
242 : // check whether element is anisotropic (and which case)
243 0 : AnisotropyState state = is_anisotropic(q, aaPos, thresholdRatio);
244 0 : if (state == ISOTROPIC)
245 : return state;
246 :
247 0 : if (state == QUAD_SHORTX)
248 : {
249 : Grid::SecureEdgeContainer assEd;
250 : grid.associated_elements_sorted(assEd, q);
251 0 : sidesOut.push_back(assEd[1]);
252 0 : sidesOut.push_back(assEd[3]);
253 :
254 : return state;
255 : }
256 :
257 0 : if (state == QUAD_SHORTY)
258 : {
259 : Grid::SecureEdgeContainer assEd;
260 : grid.associated_elements_sorted(assEd, q);
261 0 : sidesOut.push_back(assEd[0]);
262 0 : sidesOut.push_back(assEd[2]);
263 :
264 : return state;
265 : }
266 :
267 : return state;
268 : }
269 :
270 :
271 :
272 : template <typename TAAPos>
273 : AnisotropyState close_sides_of_anisotropic_elem
274 : (
275 : Volume* elem,
276 : Grid& grid,
277 : const TAAPos& aaPos,
278 : number thresholdRatio,
279 : std::vector<Face*>& sidesOut
280 : )
281 : {
282 : // treat prism case
283 : Prism* prism = dynamic_cast<Prism*>(elem);
284 : if (prism)
285 : {
286 : AnisotropyState state = is_anisotropic(prism, aaPos, thresholdRatio);
287 : if (state == ISOTROPIC)
288 : return state;
289 :
290 : // flat case
291 : if (state == PRISM_FLAT)
292 : {
293 : Face* side = grid.get_side(prism, 0);
294 : if (side)
295 : sidesOut.push_back(side);
296 :
297 : side = grid.get_side(prism, 4);
298 : if (side)
299 : sidesOut.push_back(side);
300 : }
301 :
302 : // long case
303 : if (state == PRISM_LONG)
304 : {
305 : Face* side = grid.get_side(prism, 1);
306 : if (side)
307 : sidesOut.push_back(side);
308 :
309 : side = grid.get_side(prism, 2);
310 : if (side)
311 : sidesOut.push_back(side);
312 :
313 : side = grid.get_side(prism, 3);
314 : if (side)
315 : sidesOut.push_back(side);
316 : }
317 :
318 : return state;
319 : }
320 :
321 :
322 : // treat heaxahedron case
323 : Hexahedron* hex = dynamic_cast<Hexahedron*>(elem);
324 : if (hex)
325 : {
326 : AnisotropyState state = is_anisotropic(hex, aaPos, thresholdRatio);
327 : if (state == ISOTROPIC)
328 : return state;
329 :
330 :
331 : // short in x direction
332 : if (state == HEX_SHORTX || state == HEX_SHORTXY || state == HEX_SHORTXZ)
333 : {
334 : Face* side = grid.get_side(hex, 2);
335 : if (side)
336 : sidesOut.push_back(side);
337 :
338 : side = grid.get_side(hex, 4);
339 : if (side)
340 : sidesOut.push_back(side);
341 : }
342 :
343 : // short in y direction
344 : if (state == HEX_SHORTY || state == HEX_SHORTXY || state == HEX_SHORTYZ)
345 : {
346 : Face* side = grid.get_side(hex, 1);
347 : if (side)
348 : sidesOut.push_back(side);
349 :
350 : side = grid.get_side(hex, 3);
351 : if (side)
352 : sidesOut.push_back(side);
353 : }
354 :
355 : // short in z direction
356 : if (state == HEX_SHORTZ || state == HEX_SHORTXZ || state == HEX_SHORTYZ)
357 : {
358 : Face* side = grid.get_side(hex, 0);
359 : if (side)
360 : sidesOut.push_back(side);
361 :
362 : side = grid.get_side(hex, 5);
363 : if (side)
364 : sidesOut.push_back(side);
365 : }
366 :
367 : return state;
368 : }
369 :
370 :
371 : // other elements cannot be anisotropic
372 : return ISOTROPIC;
373 : }
374 :
375 :
376 :
377 :
378 :
379 : template <typename TAAPos>
380 : AnisotropyState long_edges_of_anisotropic_elem
381 : (
382 : Edge* elem,
383 : Grid& grid,
384 : const TAAPos& aaPos,
385 : number thresholdRatio,
386 : std::vector<Edge*>& longEdges
387 : )
388 : {
389 : return ISOTROPIC;
390 : }
391 :
392 :
393 : template <typename TAAPos>
394 : AnisotropyState long_edges_of_anisotropic_elem
395 : (
396 : Face* elem,
397 : Grid& grid,
398 : const TAAPos& aaPos,
399 : number thresholdRatio,
400 : std::vector<Edge*>& longEdges
401 : )
402 : {
403 0 : return close_sides_of_anisotropic_elem(elem, grid, aaPos, thresholdRatio, longEdges);
404 : }
405 :
406 :
407 : template <typename TAAPos>
408 0 : AnisotropyState long_edges_of_anisotropic_elem
409 : (
410 : Volume* elem,
411 : Grid& grid,
412 : const TAAPos& aaPos,
413 : number thresholdRatio,
414 : std::vector<Edge*>& longEdges
415 : )
416 : {
417 : // treat prism case
418 0 : Prism* prism = dynamic_cast<Prism*>(elem);
419 0 : if (prism)
420 : {
421 0 : AnisotropyState state = is_anisotropic(prism, aaPos, thresholdRatio);
422 0 : if (state == ISOTROPIC)
423 : return state;
424 :
425 : // flat case
426 0 : if (state == PRISM_FLAT)
427 : {
428 : Grid::SecureEdgeContainer assEd;
429 : grid.associated_elements_sorted(assEd, prism);
430 :
431 0 : UG_COND_THROW(assEd.size() < 9, "Prism needs to have 9 edges, but only "
432 : << assEd.size() << " found.");
433 :
434 0 : longEdges.reserve(longEdges.size() + 6);
435 0 : longEdges.push_back(assEd[0]);
436 0 : longEdges.push_back(assEd[1]);
437 0 : longEdges.push_back(assEd[2]);
438 0 : longEdges.push_back(assEd[6]);
439 0 : longEdges.push_back(assEd[7]);
440 0 : longEdges.push_back(assEd[8]);
441 :
442 : return state;
443 : }
444 :
445 : // long case
446 0 : if (state == PRISM_LONG)
447 : {
448 : Grid::SecureEdgeContainer assEd;
449 : grid.associated_elements_sorted(assEd, prism);
450 :
451 0 : UG_COND_THROW(assEd.size() < 9, "Prism needs to have 9 edges, but only "
452 : << assEd.size() << " found.");
453 :
454 0 : longEdges.reserve(longEdges.size() + 3);
455 0 : longEdges.push_back(assEd[3]);
456 0 : longEdges.push_back(assEd[4]);
457 0 : longEdges.push_back(assEd[5]);
458 :
459 : return state;
460 : }
461 :
462 : return state;
463 : }
464 :
465 :
466 : // treat heaxahedron case
467 0 : Hexahedron* hex = dynamic_cast<Hexahedron*>(elem);
468 0 : if (hex)
469 : {
470 0 : AnisotropyState state = is_anisotropic(hex, aaPos, thresholdRatio);
471 : if (state == ISOTROPIC)
472 : return state;
473 :
474 :
475 : // short in x direction
476 : if (state == HEX_SHORTX)
477 : {
478 : Grid::SecureEdgeContainer assEd;
479 : grid.associated_elements_sorted(assEd, hex);
480 :
481 0 : UG_COND_THROW(assEd.size() < 12, "Hexahedron needs to have 12 edges, but only "
482 : << assEd.size() << " found.");
483 :
484 0 : longEdges.reserve(longEdges.size() + 8);
485 0 : longEdges.push_back(assEd[1]);
486 0 : longEdges.push_back(assEd[3]);
487 0 : longEdges.push_back(assEd[4]);
488 0 : longEdges.push_back(assEd[5]);
489 0 : longEdges.push_back(assEd[6]);
490 0 : longEdges.push_back(assEd[7]);
491 0 : longEdges.push_back(assEd[8]);
492 0 : longEdges.push_back(assEd[11]);
493 :
494 : return state;
495 : }
496 :
497 : // short in y direction
498 : if (state == HEX_SHORTY)
499 : {
500 : Grid::SecureEdgeContainer assEd;
501 : grid.associated_elements_sorted(assEd, hex);
502 :
503 0 : UG_COND_THROW(assEd.size() < 12, "Hexahedron needs to have 12 edges, but only "
504 : << assEd.size() << " found.");
505 :
506 0 : longEdges.reserve(longEdges.size() + 8);
507 0 : longEdges.push_back(assEd[0]);
508 0 : longEdges.push_back(assEd[2]);
509 0 : longEdges.push_back(assEd[4]);
510 0 : longEdges.push_back(assEd[5]);
511 0 : longEdges.push_back(assEd[6]);
512 0 : longEdges.push_back(assEd[7]);
513 0 : longEdges.push_back(assEd[8]);
514 0 : longEdges.push_back(assEd[10]);
515 :
516 : return state;
517 : }
518 :
519 : // short in z direction
520 : if (state == HEX_SHORTZ)
521 : {
522 : Grid::SecureEdgeContainer assEd;
523 : grid.associated_elements_sorted(assEd, hex);
524 :
525 0 : UG_COND_THROW(assEd.size() < 12, "Hexahedron needs to have 12 edges, but only "
526 : << assEd.size() << " found.");
527 :
528 0 : longEdges.reserve(longEdges.size() + 8);
529 0 : longEdges.push_back(assEd[0]);
530 0 : longEdges.push_back(assEd[1]);
531 0 : longEdges.push_back(assEd[2]);
532 0 : longEdges.push_back(assEd[3]);
533 0 : longEdges.push_back(assEd[8]);
534 0 : longEdges.push_back(assEd[9]);
535 0 : longEdges.push_back(assEd[10]);
536 0 : longEdges.push_back(assEd[11]);
537 :
538 : return state;
539 : }
540 :
541 : // short in x and y direction
542 : if (state == HEX_SHORTXY)
543 : {
544 : Grid::SecureEdgeContainer assEd;
545 : grid.associated_elements_sorted(assEd, hex);
546 :
547 0 : UG_COND_THROW(assEd.size() < 12, "Hexahedron needs to have 12 edges, but only "
548 : << assEd.size() << " found.");
549 :
550 0 : longEdges.reserve(longEdges.size() + 4);
551 0 : longEdges.push_back(assEd[4]);
552 0 : longEdges.push_back(assEd[5]);
553 0 : longEdges.push_back(assEd[6]);
554 0 : longEdges.push_back(assEd[7]);
555 :
556 : return state;
557 : }
558 :
559 : // short in x and z direction
560 : if (state == HEX_SHORTXZ)
561 : {
562 : Grid::SecureEdgeContainer assEd;
563 : grid.associated_elements_sorted(assEd, hex);
564 :
565 0 : UG_COND_THROW(assEd.size() < 12, "Hexahedron needs to have 12 edges, but only "
566 : << assEd.size() << " found.");
567 :
568 0 : longEdges.reserve(longEdges.size() + 4);
569 0 : longEdges.push_back(assEd[1]);
570 0 : longEdges.push_back(assEd[3]);
571 0 : longEdges.push_back(assEd[9]);
572 0 : longEdges.push_back(assEd[11]);
573 :
574 : return state;
575 : }
576 :
577 : // short in y and z direction
578 : if (state == HEX_SHORTYZ)
579 : {
580 : Grid::SecureEdgeContainer assEd;
581 : grid.associated_elements_sorted(assEd, hex);
582 :
583 0 : UG_COND_THROW(assEd.size() < 12, "Hexahedron needs to have 12 edges, but only "
584 : << assEd.size() << " found.");
585 :
586 0 : longEdges.reserve(longEdges.size() + 4);
587 0 : longEdges.push_back(assEd[0]);
588 0 : longEdges.push_back(assEd[2]);
589 0 : longEdges.push_back(assEd[8]);
590 0 : longEdges.push_back(assEd[10]);
591 :
592 : return state;
593 : }
594 :
595 : return state;
596 : }
597 :
598 :
599 : // other elements cannot be anisotropic
600 : return ISOTROPIC;
601 : }
602 :
603 : } // namespace ug
604 :
|
