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 SHEAR:
62 0 : its_name = "Shear";
63 0 : break;
64 :
65 0 : case SHAPE:
66 0 : its_name = "Shape";
67 0 : break;
68 :
69 0 : case MAX_ANGLE:
70 0 : its_name = "Maximum Angle";
71 0 : break;
72 :
73 0 : case MIN_ANGLE:
74 0 : its_name = "Minimum Angle";
75 0 : break;
76 :
77 0 : case MAX_DIHEDRAL_ANGLE:
78 0 : its_name = "Maximum Dihedral Angle";
79 0 : break;
80 :
81 0 : case MIN_DIHEDRAL_ANGLE:
82 0 : its_name = "Minimum Dihedral Angle";
83 0 : break;
84 :
85 0 : case CONDITION:
86 0 : its_name = "Condition Number";
87 0 : break;
88 :
89 0 : case DISTORTION:
90 0 : its_name = "Distortion";
91 0 : break;
92 :
93 0 : case TAPER:
94 0 : its_name = "Taper";
95 0 : break;
96 :
97 0 : case WARP:
98 0 : its_name = "Warp";
99 0 : break;
100 :
101 0 : case STRETCH:
102 0 : its_name = "Stretch";
103 0 : break;
104 :
105 0 : case DIAGONAL:
106 0 : its_name = "Diagonal";
107 0 : break;
108 :
109 0 : case ASPECT_RATIO_BETA:
110 0 : its_name = "AR Beta";
111 0 : break;
112 :
113 0 : case ASPECT_RATIO_GAMMA:
114 0 : its_name = "AR Gamma";
115 0 : break;
116 :
117 0 : case SIZE:
118 0 : its_name = "Size";
119 0 : break;
120 :
121 0 : case JACOBIAN:
122 0 : its_name = "Jacobian";
123 0 : break;
124 :
125 0 : case SCALED_JACOBIAN:
126 0 : its_name = "Scaled Jacobian";
127 0 : break;
128 :
129 0 : default:
130 0 : its_name = "Unknown";
131 0 : break;
132 : }
133 :
134 12 : return its_name;
135 : }
136 :
137 :
138 :
139 :
140 :
141 : /**
142 : * This function returns a string containing a short
143 : * description of q. Useful for asking the enum what
144 : * it computes.
145 : */
146 0 : std::string Quality::describe (const ElemQuality q)
147 : {
148 :
149 0 : std::ostringstream desc;
150 :
151 0 : switch (q)
152 : {
153 :
154 0 : case EDGE_LENGTH_RATIO:
155 : case ASPECT_RATIO:
156 : desc << "Max edge length ratio\n"
157 : << "at element center.\n"
158 : << '\n'
159 : << "Suggested ranges:\n"
160 : << "Hexes: (1 -> 4)\n"
161 0 : << "Quads: (1 -> 4)";
162 0 : break;
163 :
164 0 : case SKEW:
165 : desc << "Maximum |cos A|, where A\n"
166 : << "is the angle between edges\n"
167 : << "at element center.\n"
168 : << '\n'
169 : << "Suggested ranges:\n"
170 : << "Hexes: (0 -> 0.5)\n"
171 0 : << "Quads: (0 -> 0.5)";
172 0 : break;
173 :
174 0 : case SHEAR:
175 : desc << "LIBMESH_DIM / K(Js)\n"
176 : << '\n'
177 : << "LIBMESH_DIM = element dimension.\n"
178 : << "K(Js) = Condition number of \n"
179 : << " Jacobian skew matrix.\n"
180 : << '\n'
181 : << "Suggested ranges:\n"
182 : << "Hexes(LIBMESH_DIM=3): (0.3 -> 1)\n"
183 0 : << "Quads(LIBMESH_DIM=2): (0.3 -> 1)";
184 0 : break;
185 :
186 0 : case SHAPE:
187 : desc << "LIBMESH_DIM / K(Jw)\n"
188 : << '\n'
189 : << "LIBMESH_DIM = element dimension.\n"
190 : << "K(Jw) = Condition number of \n"
191 : << " weighted Jacobian\n"
192 : << " matrix.\n"
193 : << '\n'
194 : << "Suggested ranges:\n"
195 : << "Hexes(LIBMESH_DIM=3): (0.3 -> 1)\n"
196 : << "Tets(LIBMESH_DIM=3): (0.2 -> 1)\n"
197 0 : << "Quads(LIBMESH_DIM=2): (0.3 -> 1).";
198 0 : break;
199 :
200 0 : case MAX_ANGLE:
201 : desc << "Largest angle between all adjacent pairs of edges (in 2D, sides).\n"
202 : << '\n'
203 : << "Suggested ranges:\n"
204 : << "Quads: (90 -> 135)\n"
205 0 : << "Triangles: (60 -> 90)";
206 0 : break;
207 :
208 0 : case MIN_ANGLE:
209 : desc << "Smallest angle between all adjacent pairs of edges (in 2D, sides).\n"
210 : << '\n'
211 : << "Suggested ranges:\n"
212 : << "Quads: (45 -> 90)\n"
213 0 : << "Triangles: (30 -> 60)";
214 0 : break;
215 :
216 0 : case MAX_DIHEDRAL_ANGLE:
217 : desc << "Largest angle between all adjacent pairs of sides (in 2D, equivalent to MAX_ANGLE).\n"
218 : << "In 3D, this is the largest unoriented angle between adjacent side planes, in the range [0, 90].\n"
219 : << '\n'
220 : << "Suggested ranges:\n"
221 : << "Quads: (90 -> 135)\n"
222 : << "Triangles: (60 -> 90)\n"
223 0 : << "C0Polyhedra: (60 -> 90)";
224 0 : break;
225 :
226 0 : case MIN_DIHEDRAL_ANGLE:
227 : desc << "Smallest angle between all adjacent pairs of sides (in 2D, equivalent to MIN_ANGLE).\n"
228 : << "In 3D, this is the smallest unoriented angle between adjacent side planes, in the range [0, 90].\n"
229 : << '\n'
230 : << "Suggested ranges:\n"
231 : << "Quads: (45 -> 90)\n"
232 : << "Triangles: (30 -> 60)\n"
233 0 : << "C0Polyhedra: (30 -> 90)";
234 0 : break;
235 :
236 0 : case CONDITION:
237 : desc << "Condition number of the\n"
238 : << "Jacobian matrix.\n"
239 : << '\n'
240 : << "Suggested ranges:\n"
241 : << "Quads: (1 -> 4)\n"
242 : << "Hexes: (1 -> 8)\n"
243 : << "Tris: (1 -> 1.3)\n"
244 0 : << "Tets: (1 -> 3)";
245 0 : break;
246 :
247 0 : case DISTORTION:
248 : desc << "min |J| * A / <A>\n"
249 : << '\n'
250 : << "|J| = norm of Jacobian matrix\n"
251 : << " A = actual area\n"
252 : << "<A> = reference area\n"
253 : << '\n'
254 : << "Suggested ranges:\n"
255 : << "Quads: (0.6 -> 1), <A>=4\n"
256 : << "Hexes: (0.6 -> 1), <A>=8\n"
257 : << "Tris: (0.6 -> 1), <A>=1/2\n"
258 0 : << "Tets: (0.6 -> 1), <A>=1/6";
259 0 : break;
260 :
261 0 : case TAPER:
262 : desc << "Maximum ratio of lengths\n"
263 : << "derived from opposite edges.\n"
264 : << '\n'
265 : << "Suggested ranges:\n"
266 : << "Quads: (0.7 -> 1)\n"
267 0 : << "Hexes: (0.4 -> 1)";
268 0 : break;
269 :
270 0 : case WARP:
271 : desc << "cos D\n"
272 : << '\n'
273 : << "D = minimum dihedral angle\n"
274 : << " formed by diagonals.\n"
275 : << '\n'
276 : << "Suggested ranges:\n"
277 0 : << "Quads: (0.9 -> 1)";
278 0 : break;
279 :
280 0 : case STRETCH:
281 : desc << "Sqrt(3) * L_min / L_max\n"
282 : << '\n'
283 : << "L_min = minimum edge length.\n"
284 : << "L_max = maximum edge length.\n"
285 : << '\n'
286 : << "Suggested ranges:\n"
287 : << "Quads: (0.25 -> 1)\n"
288 0 : << "Hexes: (0.25 -> 1)";
289 0 : break;
290 :
291 0 : case DIAGONAL:
292 : desc << "D_min / D_max\n"
293 : << '\n'
294 : << "D_min = minimum diagonal.\n"
295 : << "D_max = maximum diagonal.\n"
296 : << '\n'
297 : << "Suggested ranges:\n"
298 0 : << "Hexes: (0.65 -> 1)";
299 0 : break;
300 :
301 0 : case ASPECT_RATIO_BETA:
302 : desc << "CR / (3 * IR)\n"
303 : << '\n'
304 : << "CR = circumsphere radius\n"
305 : << "IR = inscribed sphere radius\n"
306 : << '\n'
307 : << "Suggested ranges:\n"
308 0 : << "Tets: (1 -> 3)";
309 0 : break;
310 :
311 0 : case ASPECT_RATIO_GAMMA:
312 : desc << "S^(3/2) / 8.479670 * V\n"
313 : << '\n'
314 : << "S = sum(si*si/6)\n"
315 : << "si = edge length\n"
316 : << "V = volume\n"
317 : << '\n'
318 : << "Suggested ranges:\n"
319 0 : << "Tets: (1 -> 3)";
320 0 : break;
321 :
322 0 : case SIZE:
323 : desc << "min (|J|, |1/J|)\n"
324 : << '\n'
325 : << "|J| = norm of Jacobian matrix.\n"
326 : << '\n'
327 : << "Suggested ranges:\n"
328 : << "Quads: (0.3 -> 1)\n"
329 : << "Hexes: (0.5 -> 1)\n"
330 : << "Tris: (0.25 -> 1)\n"
331 0 : << "Tets: (0.2 -> 1)";
332 0 : break;
333 :
334 0 : case JACOBIAN:
335 : case SCALED_JACOBIAN:
336 : desc << "Minimum nodal Jacobian.\n"
337 : << "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"
338 : << "In the SCALED_JACOBIAN case, we also then divide by the lengths of each of the associated edges.\n"
339 : << "For Pyramid elements where four edges meet at the apex node, special handling is required.\n"
340 : << '\n'
341 : << "Suggested acceptable ranges (from Cubit documentation) for SCALED_JACOBIAN metric:\n"
342 : << "Quads/Hexes: (0.5 -> 1)\n"
343 0 : << "Tris/Tets: (0.2 -> 1.0)";
344 0 : break;
345 :
346 0 : default:
347 0 : desc << "Unknown";
348 0 : break;
349 : }
350 :
351 0 : return desc.str();
352 0 : }
353 :
354 :
355 24 : std::vector<ElemQuality> Quality::valid(const ElemType t)
356 : {
357 2 : std::vector<ElemQuality> v;
358 :
359 24 : switch (t)
360 : {
361 0 : case EDGE2:
362 : case EDGE3:
363 : case EDGE4:
364 : {
365 : // None yet
366 0 : break;
367 : }
368 :
369 0 : case TRI3:
370 : case TRISHELL3:
371 : case TRI6:
372 : case TRI7:
373 : {
374 0 : v = {
375 : CONDITION,
376 : DISTORTION,
377 : EDGE_LENGTH_RATIO,
378 : JACOBIAN,
379 : SCALED_JACOBIAN,
380 : MAX_ANGLE,
381 : MIN_ANGLE,
382 : MAX_DIHEDRAL_ANGLE,
383 : MIN_DIHEDRAL_ANGLE,
384 : SHAPE,
385 : SIZE
386 0 : };
387 :
388 0 : break;
389 : }
390 :
391 0 : case QUAD4:
392 : case QUADSHELL4:
393 : case QUAD8:
394 : case QUADSHELL8:
395 : case QUAD9:
396 : case QUADSHELL9:
397 : {
398 0 : v = {
399 : ASPECT_RATIO,
400 : CONDITION,
401 : DISTORTION,
402 : EDGE_LENGTH_RATIO,
403 : JACOBIAN,
404 : SCALED_JACOBIAN,
405 : MAX_ANGLE,
406 : MIN_ANGLE,
407 : MAX_DIHEDRAL_ANGLE,
408 : MIN_DIHEDRAL_ANGLE,
409 : SHAPE,
410 : SHEAR,
411 : SIZE,
412 : SKEW,
413 : STRETCH,
414 : TAPER,
415 : WARP
416 0 : };
417 :
418 0 : break;
419 : }
420 :
421 0 : case TET4:
422 : case TET10:
423 : case TET14:
424 : {
425 0 : v = {
426 : ASPECT_RATIO_BETA,
427 : ASPECT_RATIO_GAMMA,
428 : CONDITION,
429 : DISTORTION,
430 : JACOBIAN,
431 : SCALED_JACOBIAN,
432 : MAX_ANGLE,
433 : MIN_ANGLE,
434 : MAX_DIHEDRAL_ANGLE,
435 : MIN_DIHEDRAL_ANGLE,
436 : SHAPE,
437 : SIZE
438 0 : };
439 :
440 0 : break;
441 : }
442 :
443 0 : case HEX8:
444 : case HEX20:
445 : case HEX27:
446 : {
447 0 : v = {
448 : ASPECT_RATIO,
449 : CONDITION,
450 : DIAGONAL,
451 : DISTORTION,
452 : JACOBIAN,
453 : SCALED_JACOBIAN,
454 : MAX_ANGLE,
455 : MIN_ANGLE,
456 : MAX_DIHEDRAL_ANGLE,
457 : MIN_DIHEDRAL_ANGLE,
458 : SHAPE,
459 : SHEAR,
460 : SIZE,
461 : SKEW,
462 : STRETCH,
463 : TAPER
464 0 : };
465 :
466 0 : break;
467 : }
468 :
469 0 : case PRISM6:
470 : case PRISM18:
471 : case PRISM20:
472 : case PRISM21:
473 : {
474 0 : v = {
475 : EDGE_LENGTH_RATIO,
476 : MAX_ANGLE,
477 : MIN_ANGLE,
478 : MAX_DIHEDRAL_ANGLE,
479 : MIN_DIHEDRAL_ANGLE,
480 0 : };
481 :
482 0 : break;
483 : }
484 :
485 0 : case PYRAMID5:
486 : case PYRAMID13:
487 : case PYRAMID14:
488 : case PYRAMID18:
489 : {
490 0 : v = {
491 : EDGE_LENGTH_RATIO,
492 : MAX_ANGLE,
493 : MIN_ANGLE,
494 : MAX_DIHEDRAL_ANGLE,
495 : MIN_DIHEDRAL_ANGLE,
496 0 : };
497 :
498 0 : break;
499 : }
500 :
501 12 : case C0POLYGON:
502 : {
503 21 : v = {
504 : EDGE_LENGTH_RATIO,
505 : JACOBIAN,
506 : SCALED_JACOBIAN,
507 : MAX_ANGLE,
508 : MIN_ANGLE
509 11 : };
510 :
511 12 : break;
512 : }
513 :
514 12 : case C0POLYHEDRON:
515 : {
516 : // The generic Jacobian metrics only inspect vertices with
517 : // exactly three adjacent edges, but arbitrary polyhedra may
518 : // have vertices of higher valence.
519 3 : v = {
520 : EDGE_LENGTH_RATIO,
521 : MAX_ANGLE,
522 : MIN_ANGLE,
523 : MAX_DIHEDRAL_ANGLE,
524 : MIN_DIHEDRAL_ANGLE
525 11 : };
526 :
527 12 : break;
528 : }
529 :
530 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
531 :
532 0 : case INFEDGE2:
533 : {
534 : // None yet
535 0 : break;
536 : }
537 :
538 0 : case INFQUAD4:
539 : case INFQUAD6:
540 : case INFHEX8:
541 : case INFHEX16:
542 : case INFHEX18:
543 : case INFPRISM6:
544 : case INFPRISM12:
545 : {
546 0 : v = {
547 : MAX_ANGLE,
548 : MIN_ANGLE,
549 : MAX_DIHEDRAL_ANGLE,
550 : MIN_DIHEDRAL_ANGLE,
551 0 : };
552 :
553 0 : break;
554 : }
555 :
556 : #endif
557 :
558 :
559 0 : default:
560 0 : libmesh_error_msg("Undefined element type!");
561 : }
562 :
563 24 : return v;
564 : }
565 :
566 : } // namespace libMesh
|
