Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : // C++ includes
19 : #include <iostream>
20 : #include <sstream>
21 :
22 : // Local includes
23 : #include "libmesh/libmesh_common.h"
24 : #include "libmesh/elem_quality.h"
25 : #include "libmesh/enum_elem_type.h"
26 : #include "libmesh/enum_elem_quality.h"
27 :
28 :
29 : namespace libMesh
30 : {
31 :
32 : // ------------------------------------------------------------
33 : // Quality function definitions
34 :
35 : /**
36 : * This function returns a string containing some name
37 : * for q. Useful for asking the enum what its name is.
38 : * I added this since you may want a simple way to attach
39 : * a name or description to the ElemQuality enums.
40 : * It can be removed if it is found to be useless.
41 : */
42 12 : std::string Quality::name (const ElemQuality q)
43 : {
44 1 : std::string its_name;
45 :
46 12 : switch (q)
47 : {
48 :
49 1 : case EDGE_LENGTH_RATIO:
50 1 : its_name = "Edge Length Ratio";
51 1 : break;
52 :
53 0 : case ASPECT_RATIO:
54 0 : its_name = "Aspect Ratio";
55 0 : break;
56 :
57 0 : case SKEW:
58 0 : its_name = "Skew";
59 0 : break;
60 :
61 0 : case SKEW_ANGLE:
62 0 : its_name = "Skew Angle";
63 0 : break;
64 :
65 0 : case SHEAR:
66 0 : its_name = "Shear";
67 0 : break;
68 :
69 0 : case SHAPE:
70 0 : its_name = "Shape";
71 0 : break;
72 :
73 0 : case MAX_ANGLE:
74 0 : its_name = "Maximum Angle";
75 0 : break;
76 :
77 0 : case MIN_ANGLE:
78 0 : its_name = "Minimum Angle";
79 0 : break;
80 :
81 0 : case MAX_DIHEDRAL_ANGLE:
82 0 : its_name = "Maximum Dihedral Angle";
83 0 : break;
84 :
85 0 : case MIN_DIHEDRAL_ANGLE:
86 0 : its_name = "Minimum Dihedral Angle";
87 0 : break;
88 :
89 0 : case CONDITION:
90 0 : its_name = "Condition Number";
91 0 : break;
92 :
93 0 : case DISTORTION:
94 0 : its_name = "Distortion";
95 0 : break;
96 :
97 0 : case TAPER:
98 0 : its_name = "Taper";
99 0 : break;
100 :
101 0 : case WARP:
102 0 : its_name = "Warp";
103 0 : break;
104 :
105 0 : case STRETCH:
106 0 : its_name = "Stretch";
107 0 : break;
108 :
109 0 : case DIAGONAL:
110 0 : its_name = "Diagonal";
111 0 : break;
112 :
113 0 : case ASPECT_RATIO_BETA:
114 0 : its_name = "AR Beta";
115 0 : break;
116 :
117 0 : case ASPECT_RATIO_GAMMA:
118 0 : its_name = "AR Gamma";
119 0 : break;
120 :
121 0 : case SIZE:
122 0 : its_name = "Size";
123 0 : break;
124 :
125 0 : case JACOBIAN:
126 0 : its_name = "Jacobian";
127 0 : break;
128 :
129 0 : case SCALED_JACOBIAN:
130 0 : its_name = "Scaled Jacobian";
131 0 : break;
132 :
133 0 : default:
134 0 : its_name = "Unknown";
135 0 : break;
136 : }
137 :
138 12 : return its_name;
139 : }
140 :
141 :
142 :
143 :
144 :
145 : /**
146 : * This function returns a string containing a short
147 : * description of q. Useful for asking the enum what
148 : * it computes.
149 : */
150 0 : std::string Quality::describe (const ElemQuality q)
151 : {
152 :
153 0 : std::ostringstream desc;
154 :
155 0 : switch (q)
156 : {
157 :
158 0 : case EDGE_LENGTH_RATIO:
159 : case ASPECT_RATIO:
160 : desc << "Max edge length ratio\n"
161 : << "at element center.\n"
162 : << '\n'
163 : << "Suggested ranges:\n"
164 : << "Hexes: (1 -> 4)\n"
165 0 : << "Quads: (1 -> 4)";
166 0 : break;
167 :
168 0 : case SKEW:
169 : desc << "Knupp's algebraic skew metric,\n"
170 : << "based on the nodal Jacobian\n"
171 : << "skew matrices. 1 is ideal,\n"
172 : << "smaller values are worse.\n"
173 : << '\n'
174 : << "Suggested ranges:\n"
175 : << "Hexes: (0.3 -> 1)\n"
176 0 : << "Quads: (0.3 -> 1)";
177 0 : break;
178 :
179 0 : case SKEW_ANGLE:
180 : desc << "Maximum |cos A|, where A\n"
181 : << "is the angle between edges\n"
182 : << "at element center.\n"
183 : << "NOTE: some degenerate elements\n"
184 : << "score 0 if zero-length along principal axis.\n"
185 : << "Suggested ranges:\n"
186 : << "Hexes: (0 -> 0.5)\n"
187 0 : << "Quads: (0 -> 0.5)";
188 0 : break;
189 :
190 0 : case SHEAR:
191 : desc << "LIBMESH_DIM / K(Js)\n"
192 : << '\n'
193 : << "LIBMESH_DIM = element dimension.\n"
194 : << "K(Js) = Condition number of \n"
195 : << " Jacobian skew matrix.\n"
196 : << '\n'
197 : << "Suggested ranges:\n"
198 : << "Hexes(LIBMESH_DIM=3): (0.3 -> 1)\n"
199 0 : << "Quads(LIBMESH_DIM=2): (0.3 -> 1)";
200 0 : break;
201 :
202 0 : case SHAPE:
203 : desc << "LIBMESH_DIM / K(Jw)\n"
204 : << '\n'
205 : << "LIBMESH_DIM = element dimension.\n"
206 : << "K(Jw) = Condition number of \n"
207 : << " weighted Jacobian\n"
208 : << " matrix.\n"
209 : << '\n'
210 : << "Suggested ranges:\n"
211 : << "Hexes(LIBMESH_DIM=3): (0.3 -> 1)\n"
212 : << "Tets(LIBMESH_DIM=3): (0.2 -> 1)\n"
213 0 : << "Quads(LIBMESH_DIM=2): (0.3 -> 1).";
214 0 : break;
215 :
216 0 : case MAX_ANGLE:
217 : desc << "Largest angle between all adjacent pairs of edges (in 2D, sides).\n"
218 : << '\n'
219 : << "Suggested ranges:\n"
220 : << "Quads: (90 -> 135)\n"
221 0 : << "Triangles: (60 -> 90)";
222 0 : break;
223 :
224 0 : case MIN_ANGLE:
225 : desc << "Smallest angle between all adjacent pairs of edges (in 2D, sides).\n"
226 : << '\n'
227 : << "Suggested ranges:\n"
228 : << "Quads: (45 -> 90)\n"
229 0 : << "Triangles: (30 -> 60)";
230 0 : break;
231 :
232 0 : case MAX_DIHEDRAL_ANGLE:
233 : desc << "Largest angle between all adjacent pairs of sides (in 2D, equivalent to MAX_ANGLE).\n"
234 : << "In 3D, this is the largest unoriented angle between adjacent side planes, in the range [0, 90].\n"
235 : << '\n'
236 : << "Suggested ranges:\n"
237 : << "Quads: (90 -> 135)\n"
238 : << "Triangles: (60 -> 90)\n"
239 0 : << "C0Polyhedra: (60 -> 90)";
240 0 : break;
241 :
242 0 : case MIN_DIHEDRAL_ANGLE:
243 : desc << "Smallest angle between all adjacent pairs of sides (in 2D, equivalent to MIN_ANGLE).\n"
244 : << "In 3D, this is the smallest unoriented angle between adjacent side planes, in the range [0, 90].\n"
245 : << '\n'
246 : << "Suggested ranges:\n"
247 : << "Quads: (45 -> 90)\n"
248 : << "Triangles: (30 -> 60)\n"
249 0 : << "C0Polyhedra: (30 -> 90)";
250 0 : break;
251 :
252 0 : case CONDITION:
253 : desc << "Maximum condition number of\n"
254 : << "the Jacobian matrix at each\n"
255 : << "corner, relative to an ideal\n"
256 : << "(regular) element. 1 is ideal,\n"
257 : << "larger values are worse.\n"
258 : << '\n'
259 : << "Suggested ranges:\n"
260 : << "Quads: (1 -> 4)\n"
261 : << "Hexes: (1 -> 8)\n"
262 : << "Tris: (1 -> 1.3)\n"
263 0 : << "Tets: (1 -> 3)";
264 0 : break;
265 :
266 0 : case DISTORTION:
267 : desc << "min |J| * A / <A>\n"
268 : << '\n'
269 : << "|J| = norm of Jacobian matrix\n"
270 : << " A = actual area\n"
271 : << "<A> = reference area\n"
272 : << '\n'
273 : << "Suggested ranges:\n"
274 : << "Quads: (0.6 -> 1), <A>=4\n"
275 : << "Hexes: (0.6 -> 1), <A>=8\n"
276 : << "Tris: (0.6 -> 1), <A>=1/2\n"
277 0 : << "Tets: (0.6 -> 1), <A>=1/6";
278 0 : break;
279 :
280 0 : case TAPER:
281 : desc << "Maximum ratio of lengths\n"
282 : << "derived from opposite edges.\n"
283 : << '\n'
284 : << "Suggested ranges:\n"
285 : << "Quads: (0.7 -> 1)\n"
286 0 : << "Hexes: (0.4 -> 1)";
287 0 : break;
288 :
289 0 : case WARP:
290 : desc << "cos D\n"
291 : << '\n'
292 : << "D = minimum dihedral angle\n"
293 : << " formed by diagonals.\n"
294 : << '\n'
295 : << "Suggested ranges:\n"
296 0 : << "Quads: (0.9 -> 1)";
297 0 : break;
298 :
299 0 : case STRETCH:
300 : desc << "Sqrt(3) * L_min / L_max\n"
301 : << '\n'
302 : << "L_min = minimum edge length.\n"
303 : << "L_max = maximum edge length.\n"
304 : << '\n'
305 : << "Suggested ranges:\n"
306 : << "Quads: (0.25 -> 1)\n"
307 0 : << "Hexes: (0.25 -> 1)";
308 0 : break;
309 :
310 0 : case DIAGONAL:
311 : desc << "D_min / D_max\n"
312 : << '\n'
313 : << "D_min = minimum diagonal.\n"
314 : << "D_max = maximum diagonal.\n"
315 : << '\n'
316 : << "Suggested ranges:\n"
317 0 : << "Hexes: (0.65 -> 1)";
318 0 : break;
319 :
320 0 : case ASPECT_RATIO_BETA:
321 : desc << "CR / (3 * IR)\n"
322 : << '\n'
323 : << "CR = circumsphere radius\n"
324 : << "IR = inscribed sphere radius\n"
325 : << '\n'
326 : << "Suggested ranges:\n"
327 0 : << "Tets: (1 -> 3)";
328 0 : break;
329 :
330 0 : case ASPECT_RATIO_GAMMA:
331 : desc << "S^(3/2) / 8.479670 * V\n"
332 : << '\n'
333 : << "S = sum(si*si/6)\n"
334 : << "si = edge length\n"
335 : << "V = volume\n"
336 : << '\n'
337 : << "Suggested ranges:\n"
338 0 : << "Tets: (1 -> 3)";
339 0 : break;
340 :
341 0 : case SIZE:
342 : desc << "Relative size: min(J, 1/J),\n"
343 : << "where J is the determinant of\n"
344 : << "the nodal Jacobian relative to\n"
345 : << "an ideal element of the same\n"
346 : << "volume. 1 for a uniform\n"
347 : << "(affine) element.\n"
348 : << '\n'
349 : << "Suggested ranges:\n"
350 : << "Quads: (0.3 -> 1)\n"
351 : << "Hexes: (0.5 -> 1)\n"
352 : << "Tris: (0.25 -> 1)\n"
353 0 : << "Tets: (0.2 -> 1)";
354 0 : break;
355 :
356 0 : case JACOBIAN:
357 : case SCALED_JACOBIAN:
358 : desc << "Minimum nodal Jacobian.\n"
359 : << "The nodal Jacobians are computed by taking the cross product (2D) or scalar product (3D) of the adjacent edges that meet at that node.\n"
360 : << "In the SCALED_JACOBIAN case, we also then divide by the lengths of each of the associated edges.\n"
361 : << "For Pyramid elements where four edges meet at the apex node, special handling is required.\n"
362 : << '\n'
363 : << "Suggested acceptable ranges (from Cubit documentation) for SCALED_JACOBIAN metric:\n"
364 : << "Quads/Hexes: (0.5 -> 1)\n"
365 0 : << "Tris/Tets: (0.2 -> 1.0)";
366 0 : break;
367 :
368 0 : default:
369 0 : desc << "Unknown";
370 0 : break;
371 : }
372 :
373 0 : return desc.str();
374 0 : }
375 :
376 :
377 22512 : std::vector<ElemQuality> Quality::valid(const ElemType t)
378 : {
379 1876 : std::vector<ElemQuality> v;
380 :
381 22512 : switch (t)
382 : {
383 14 : case NODEELEM:
384 : case EDGE2:
385 : case EDGE3:
386 : case EDGE4:
387 : {
388 : // None yet
389 14 : break;
390 : }
391 :
392 768 : case TRI3:
393 : case TRISHELL3:
394 : case TRI6:
395 : case TRI7:
396 : {
397 1344 : v = {
398 : CONDITION,
399 : DISTORTION,
400 : EDGE_LENGTH_RATIO,
401 : JACOBIAN,
402 : SCALED_JACOBIAN,
403 : MAX_ANGLE,
404 : MIN_ANGLE,
405 : MAX_DIHEDRAL_ANGLE,
406 : MIN_DIHEDRAL_ANGLE,
407 : SHAPE,
408 : SIZE
409 704 : };
410 :
411 768 : break;
412 : }
413 :
414 576 : case QUAD4:
415 : case QUADSHELL4:
416 : case QUAD8:
417 : case QUADSHELL8:
418 : case QUAD9:
419 : case QUADSHELL9:
420 : {
421 1008 : v = {
422 : ASPECT_RATIO,
423 : CONDITION,
424 : DISTORTION,
425 : EDGE_LENGTH_RATIO,
426 : JACOBIAN,
427 : SCALED_JACOBIAN,
428 : MAX_ANGLE,
429 : MIN_ANGLE,
430 : MAX_DIHEDRAL_ANGLE,
431 : MIN_DIHEDRAL_ANGLE,
432 : SHAPE,
433 : SHEAR,
434 : SIZE,
435 : SKEW,
436 : SKEW_ANGLE,
437 : STRETCH,
438 : TAPER,
439 : WARP
440 528 : };
441 :
442 576 : break;
443 : }
444 :
445 13824 : case TET4:
446 : case TET10:
447 : case TET14:
448 : {
449 24192 : v = {
450 : ASPECT_RATIO_BETA,
451 : ASPECT_RATIO_GAMMA,
452 : CONDITION,
453 : DISTORTION,
454 : JACOBIAN,
455 : SCALED_JACOBIAN,
456 : MAX_ANGLE,
457 : MIN_ANGLE,
458 : MAX_DIHEDRAL_ANGLE,
459 : MIN_DIHEDRAL_ANGLE,
460 : SHAPE,
461 : SIZE
462 12672 : };
463 :
464 13824 : break;
465 : }
466 :
467 576 : case HEX8:
468 : case HEX20:
469 : case HEX27:
470 : {
471 1008 : v = {
472 : ASPECT_RATIO,
473 : CONDITION,
474 : DIAGONAL,
475 : DISTORTION,
476 : JACOBIAN,
477 : SCALED_JACOBIAN,
478 : MAX_ANGLE,
479 : MIN_ANGLE,
480 : MAX_DIHEDRAL_ANGLE,
481 : MIN_DIHEDRAL_ANGLE,
482 : SHAPE,
483 : SHEAR,
484 : SIZE,
485 : SKEW,
486 : SKEW_ANGLE,
487 : STRETCH,
488 : TAPER
489 528 : };
490 :
491 576 : break;
492 : }
493 :
494 1920 : case PRISM6:
495 : case PRISM15:
496 : case PRISM18:
497 : case PRISM20:
498 : case PRISM21:
499 : {
500 3360 : v = {
501 : EDGE_LENGTH_RATIO,
502 : MAX_ANGLE,
503 : MIN_ANGLE,
504 : MAX_DIHEDRAL_ANGLE,
505 : MIN_DIHEDRAL_ANGLE,
506 1760 : };
507 :
508 1920 : break;
509 : }
510 :
511 4608 : case PYRAMID5:
512 : case PYRAMID13:
513 : case PYRAMID14:
514 : case PYRAMID18:
515 : {
516 8064 : v = {
517 : EDGE_LENGTH_RATIO,
518 : MAX_ANGLE,
519 : MIN_ANGLE,
520 : MAX_DIHEDRAL_ANGLE,
521 : MIN_DIHEDRAL_ANGLE,
522 4224 : };
523 :
524 4608 : break;
525 : }
526 :
527 36 : case C0POLYGON:
528 : {
529 63 : v = {
530 : EDGE_LENGTH_RATIO,
531 : JACOBIAN,
532 : SCALED_JACOBIAN,
533 : MAX_ANGLE,
534 : MIN_ANGLE
535 33 : };
536 :
537 36 : break;
538 : }
539 :
540 36 : case C0POLYHEDRON:
541 : {
542 : // The generic Jacobian metrics only inspect vertices with
543 : // exactly three adjacent edges, but arbitrary polyhedra may
544 : // have vertices of higher valence.
545 9 : v = {
546 : EDGE_LENGTH_RATIO,
547 : MAX_ANGLE,
548 : MIN_ANGLE,
549 : MAX_DIHEDRAL_ANGLE,
550 : MIN_DIHEDRAL_ANGLE
551 33 : };
552 :
553 36 : break;
554 : }
555 :
556 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
557 :
558 0 : case INFEDGE2:
559 : {
560 : // None yet
561 0 : break;
562 : }
563 :
564 0 : case INFQUAD4:
565 : case INFQUAD6:
566 : case INFHEX8:
567 : case INFHEX16:
568 : case INFHEX18:
569 : case INFPRISM6:
570 : case INFPRISM12:
571 : {
572 0 : v = {
573 : MAX_ANGLE,
574 : MIN_ANGLE,
575 : MAX_DIHEDRAL_ANGLE,
576 : MIN_DIHEDRAL_ANGLE,
577 0 : };
578 :
579 0 : break;
580 : }
581 :
582 : #endif
583 :
584 :
585 0 : default:
586 0 : libmesh_error_msg("Undefined element type!");
587 : }
588 :
589 22512 : return v;
590 : }
591 :
592 : } // namespace libMesh
|
