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 : << '\n'
184 : << "Suggested ranges:\n"
185 : << "Hexes: (0 -> 0.5)\n"
186 0 : << "Quads: (0 -> 0.5)";
187 0 : break;
188 :
189 0 : case SHEAR:
190 : desc << "LIBMESH_DIM / K(Js)\n"
191 : << '\n'
192 : << "LIBMESH_DIM = element dimension.\n"
193 : << "K(Js) = Condition number of \n"
194 : << " Jacobian skew matrix.\n"
195 : << '\n'
196 : << "Suggested ranges:\n"
197 : << "Hexes(LIBMESH_DIM=3): (0.3 -> 1)\n"
198 0 : << "Quads(LIBMESH_DIM=2): (0.3 -> 1)";
199 0 : break;
200 :
201 0 : case SHAPE:
202 : desc << "LIBMESH_DIM / K(Jw)\n"
203 : << '\n'
204 : << "LIBMESH_DIM = element dimension.\n"
205 : << "K(Jw) = Condition number of \n"
206 : << " weighted Jacobian\n"
207 : << " matrix.\n"
208 : << '\n'
209 : << "Suggested ranges:\n"
210 : << "Hexes(LIBMESH_DIM=3): (0.3 -> 1)\n"
211 : << "Tets(LIBMESH_DIM=3): (0.2 -> 1)\n"
212 0 : << "Quads(LIBMESH_DIM=2): (0.3 -> 1).";
213 0 : break;
214 :
215 0 : case MAX_ANGLE:
216 : desc << "Largest angle between all adjacent pairs of edges (in 2D, sides).\n"
217 : << '\n'
218 : << "Suggested ranges:\n"
219 : << "Quads: (90 -> 135)\n"
220 0 : << "Triangles: (60 -> 90)";
221 0 : break;
222 :
223 0 : case MIN_ANGLE:
224 : desc << "Smallest angle between all adjacent pairs of edges (in 2D, sides).\n"
225 : << '\n'
226 : << "Suggested ranges:\n"
227 : << "Quads: (45 -> 90)\n"
228 0 : << "Triangles: (30 -> 60)";
229 0 : break;
230 :
231 0 : case MAX_DIHEDRAL_ANGLE:
232 : desc << "Largest angle between all adjacent pairs of sides (in 2D, equivalent to MAX_ANGLE).\n"
233 : << "In 3D, this is the largest unoriented angle between adjacent side planes, in the range [0, 90].\n"
234 : << '\n'
235 : << "Suggested ranges:\n"
236 : << "Quads: (90 -> 135)\n"
237 : << "Triangles: (60 -> 90)\n"
238 0 : << "C0Polyhedra: (60 -> 90)";
239 0 : break;
240 :
241 0 : case MIN_DIHEDRAL_ANGLE:
242 : desc << "Smallest angle between all adjacent pairs of sides (in 2D, equivalent to MIN_ANGLE).\n"
243 : << "In 3D, this is the smallest unoriented angle between adjacent side planes, in the range [0, 90].\n"
244 : << '\n'
245 : << "Suggested ranges:\n"
246 : << "Quads: (45 -> 90)\n"
247 : << "Triangles: (30 -> 60)\n"
248 0 : << "C0Polyhedra: (30 -> 90)";
249 0 : break;
250 :
251 0 : case CONDITION:
252 : desc << "Maximum condition number of\n"
253 : << "the Jacobian matrix at each\n"
254 : << "corner. 1 is ideal, larger\n"
255 : << "values are worse.\n"
256 : << '\n'
257 : << "Suggested ranges:\n"
258 : << "Quads: (1 -> 4)\n"
259 : << "Hexes: (1 -> 8)\n"
260 : << "Tris: (1 -> 1.3)\n"
261 0 : << "Tets: (1 -> 3)";
262 0 : break;
263 :
264 0 : case DISTORTION:
265 : desc << "min |J| * A / <A>\n"
266 : << '\n'
267 : << "|J| = norm of Jacobian matrix\n"
268 : << " A = actual area\n"
269 : << "<A> = reference area\n"
270 : << '\n'
271 : << "Suggested ranges:\n"
272 : << "Quads: (0.6 -> 1), <A>=4\n"
273 : << "Hexes: (0.6 -> 1), <A>=8\n"
274 : << "Tris: (0.6 -> 1), <A>=1/2\n"
275 0 : << "Tets: (0.6 -> 1), <A>=1/6";
276 0 : break;
277 :
278 0 : case TAPER:
279 : desc << "Maximum ratio of lengths\n"
280 : << "derived from opposite edges.\n"
281 : << '\n'
282 : << "Suggested ranges:\n"
283 : << "Quads: (0.7 -> 1)\n"
284 0 : << "Hexes: (0.4 -> 1)";
285 0 : break;
286 :
287 0 : case WARP:
288 : desc << "cos D\n"
289 : << '\n'
290 : << "D = minimum dihedral angle\n"
291 : << " formed by diagonals.\n"
292 : << '\n'
293 : << "Suggested ranges:\n"
294 0 : << "Quads: (0.9 -> 1)";
295 0 : break;
296 :
297 0 : case STRETCH:
298 : desc << "Sqrt(3) * L_min / L_max\n"
299 : << '\n'
300 : << "L_min = minimum edge length.\n"
301 : << "L_max = maximum edge length.\n"
302 : << '\n'
303 : << "Suggested ranges:\n"
304 : << "Quads: (0.25 -> 1)\n"
305 0 : << "Hexes: (0.25 -> 1)";
306 0 : break;
307 :
308 0 : case DIAGONAL:
309 : desc << "D_min / D_max\n"
310 : << '\n'
311 : << "D_min = minimum diagonal.\n"
312 : << "D_max = maximum diagonal.\n"
313 : << '\n'
314 : << "Suggested ranges:\n"
315 0 : << "Hexes: (0.65 -> 1)";
316 0 : break;
317 :
318 0 : case ASPECT_RATIO_BETA:
319 : desc << "CR / (3 * IR)\n"
320 : << '\n'
321 : << "CR = circumsphere radius\n"
322 : << "IR = inscribed sphere radius\n"
323 : << '\n'
324 : << "Suggested ranges:\n"
325 0 : << "Tets: (1 -> 3)";
326 0 : break;
327 :
328 0 : case ASPECT_RATIO_GAMMA:
329 : desc << "S^(3/2) / 8.479670 * V\n"
330 : << '\n'
331 : << "S = sum(si*si/6)\n"
332 : << "si = edge length\n"
333 : << "V = volume\n"
334 : << '\n'
335 : << "Suggested ranges:\n"
336 0 : << "Tets: (1 -> 3)";
337 0 : break;
338 :
339 0 : case SIZE:
340 : desc << "Relative size: min(J, 1/J),\n"
341 : << "where J is the determinant\n"
342 : << "of the nodal Jacobian relative\n"
343 : << "to the unit reference element.\n"
344 : << '\n'
345 : << "Suggested ranges:\n"
346 : << "Quads: (0.3 -> 1)\n"
347 : << "Hexes: (0.5 -> 1)\n"
348 : << "Tris: (0.25 -> 1)\n"
349 0 : << "Tets: (0.2 -> 1)";
350 0 : break;
351 :
352 0 : case JACOBIAN:
353 : case SCALED_JACOBIAN:
354 : desc << "Minimum nodal Jacobian.\n"
355 : << "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"
356 : << "In the SCALED_JACOBIAN case, we also then divide by the lengths of each of the associated edges.\n"
357 : << "For Pyramid elements where four edges meet at the apex node, special handling is required.\n"
358 : << '\n'
359 : << "Suggested acceptable ranges (from Cubit documentation) for SCALED_JACOBIAN metric:\n"
360 : << "Quads/Hexes: (0.5 -> 1)\n"
361 0 : << "Tris/Tets: (0.2 -> 1.0)";
362 0 : break;
363 :
364 0 : default:
365 0 : desc << "Unknown";
366 0 : break;
367 : }
368 :
369 0 : return desc.str();
370 0 : }
371 :
372 :
373 24 : std::vector<ElemQuality> Quality::valid(const ElemType t)
374 : {
375 2 : std::vector<ElemQuality> v;
376 :
377 24 : switch (t)
378 : {
379 0 : case EDGE2:
380 : case EDGE3:
381 : case EDGE4:
382 : {
383 : // None yet
384 0 : break;
385 : }
386 :
387 0 : case TRI3:
388 : case TRISHELL3:
389 : case TRI6:
390 : case TRI7:
391 : {
392 0 : v = {
393 : CONDITION,
394 : DISTORTION,
395 : EDGE_LENGTH_RATIO,
396 : JACOBIAN,
397 : SCALED_JACOBIAN,
398 : MAX_ANGLE,
399 : MIN_ANGLE,
400 : MAX_DIHEDRAL_ANGLE,
401 : MIN_DIHEDRAL_ANGLE,
402 : SHAPE,
403 : SIZE
404 0 : };
405 :
406 0 : break;
407 : }
408 :
409 0 : case QUAD4:
410 : case QUADSHELL4:
411 : case QUAD8:
412 : case QUADSHELL8:
413 : case QUAD9:
414 : case QUADSHELL9:
415 : {
416 0 : v = {
417 : ASPECT_RATIO,
418 : CONDITION,
419 : DISTORTION,
420 : EDGE_LENGTH_RATIO,
421 : JACOBIAN,
422 : SCALED_JACOBIAN,
423 : MAX_ANGLE,
424 : MIN_ANGLE,
425 : MAX_DIHEDRAL_ANGLE,
426 : MIN_DIHEDRAL_ANGLE,
427 : SHAPE,
428 : SHEAR,
429 : SIZE,
430 : SKEW,
431 : SKEW_ANGLE,
432 : STRETCH,
433 : TAPER,
434 : WARP
435 0 : };
436 :
437 0 : break;
438 : }
439 :
440 0 : case TET4:
441 : case TET10:
442 : case TET14:
443 : {
444 0 : v = {
445 : ASPECT_RATIO_BETA,
446 : ASPECT_RATIO_GAMMA,
447 : CONDITION,
448 : DISTORTION,
449 : JACOBIAN,
450 : SCALED_JACOBIAN,
451 : MAX_ANGLE,
452 : MIN_ANGLE,
453 : MAX_DIHEDRAL_ANGLE,
454 : MIN_DIHEDRAL_ANGLE,
455 : SHAPE,
456 : SIZE
457 0 : };
458 :
459 0 : break;
460 : }
461 :
462 0 : case HEX8:
463 : case HEX20:
464 : case HEX27:
465 : {
466 0 : v = {
467 : ASPECT_RATIO,
468 : CONDITION,
469 : DIAGONAL,
470 : DISTORTION,
471 : JACOBIAN,
472 : SCALED_JACOBIAN,
473 : MAX_ANGLE,
474 : MIN_ANGLE,
475 : MAX_DIHEDRAL_ANGLE,
476 : MIN_DIHEDRAL_ANGLE,
477 : SHAPE,
478 : SHEAR,
479 : SIZE,
480 : SKEW,
481 : SKEW_ANGLE,
482 : STRETCH,
483 : TAPER
484 0 : };
485 :
486 0 : break;
487 : }
488 :
489 0 : case PRISM6:
490 : case PRISM18:
491 : case PRISM20:
492 : case PRISM21:
493 : {
494 0 : v = {
495 : EDGE_LENGTH_RATIO,
496 : MAX_ANGLE,
497 : MIN_ANGLE,
498 : MAX_DIHEDRAL_ANGLE,
499 : MIN_DIHEDRAL_ANGLE,
500 0 : };
501 :
502 0 : break;
503 : }
504 :
505 0 : case PYRAMID5:
506 : case PYRAMID13:
507 : case PYRAMID14:
508 : case PYRAMID18:
509 : {
510 0 : v = {
511 : EDGE_LENGTH_RATIO,
512 : MAX_ANGLE,
513 : MIN_ANGLE,
514 : MAX_DIHEDRAL_ANGLE,
515 : MIN_DIHEDRAL_ANGLE,
516 0 : };
517 :
518 0 : break;
519 : }
520 :
521 12 : case C0POLYGON:
522 : {
523 21 : v = {
524 : EDGE_LENGTH_RATIO,
525 : JACOBIAN,
526 : SCALED_JACOBIAN,
527 : MAX_ANGLE,
528 : MIN_ANGLE
529 11 : };
530 :
531 12 : break;
532 : }
533 :
534 12 : case C0POLYHEDRON:
535 : {
536 : // The generic Jacobian metrics only inspect vertices with
537 : // exactly three adjacent edges, but arbitrary polyhedra may
538 : // have vertices of higher valence.
539 3 : v = {
540 : EDGE_LENGTH_RATIO,
541 : MAX_ANGLE,
542 : MIN_ANGLE,
543 : MAX_DIHEDRAL_ANGLE,
544 : MIN_DIHEDRAL_ANGLE
545 11 : };
546 :
547 12 : break;
548 : }
549 :
550 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
551 :
552 0 : case INFEDGE2:
553 : {
554 : // None yet
555 0 : break;
556 : }
557 :
558 0 : case INFQUAD4:
559 : case INFQUAD6:
560 : case INFHEX8:
561 : case INFHEX16:
562 : case INFHEX18:
563 : case INFPRISM6:
564 : case INFPRISM12:
565 : {
566 0 : v = {
567 : MAX_ANGLE,
568 : MIN_ANGLE,
569 : MAX_DIHEDRAL_ANGLE,
570 : MIN_DIHEDRAL_ANGLE,
571 0 : };
572 :
573 0 : break;
574 : }
575 :
576 : #endif
577 :
578 :
579 0 : default:
580 0 : libmesh_error_msg("Undefined element type!");
581 : }
582 :
583 24 : return v;
584 : }
585 :
586 : } // namespace libMesh
|
