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