elem_quality.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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"
27 
28 
29 namespace libMesh
30 {
31 
32 // ------------------------------------------------------------
33 // Quality function definitions
34 
42 std::string Quality::name (const ElemQuality q)
43 {
44  std::string its_name;
45 
46  switch (q)
47  {
48 
49  case ASPECT_RATIO:
50  its_name = "Aspect Ratio";
51  break;
52 
53  case SKEW:
54  its_name = "Skew";
55  break;
56 
57  case SHEAR:
58  its_name = "Shear";
59  break;
60 
61  case SHAPE:
62  its_name = "Shape";
63  break;
64 
65  case MAX_ANGLE:
66  its_name = "Maximum Angle";
67  break;
68 
69  case MIN_ANGLE:
70  its_name = "Minimum Angle";
71  break;
72 
73  case CONDITION:
74  its_name = "Condition Number";
75  break;
76 
77  case DISTORTION:
78  its_name = "Distortion";
79  break;
80 
81  case TAPER:
82  its_name = "Taper";
83  break;
84 
85  case WARP:
86  its_name = "Warp";
87  break;
88 
89  case STRETCH:
90  its_name = "Stretch";
91  break;
92 
93  case DIAGONAL:
94  its_name = "Diagonal";
95  break;
96 
97  case ASPECT_RATIO_BETA:
98  its_name = "AR Beta";
99  break;
100 
101  case ASPECT_RATIO_GAMMA:
102  its_name = "AR Gamma";
103  break;
104 
105  case SIZE:
106  its_name = "Size";
107  break;
108 
109  case JACOBIAN:
110  its_name = "Jacobian";
111  break;
112 
113  default:
114  its_name = "Unknown";
115  break;
116  }
117 
118  return its_name;
119 }
120 
121 
122 
123 
124 
130 std::string Quality::describe (const ElemQuality q)
131 {
132 
133  std::ostringstream desc;
134 
135  switch (q)
136  {
137 
138  case ASPECT_RATIO:
139  desc << "Max edge length ratio\n"
140  << "at element center.\n"
141  << '\n'
142  << "Suggested ranges:\n"
143  << "Hexes: (1 -> 4)\n"
144  << "Quads: (1 -> 4)";
145  break;
146 
147  case SKEW:
148  desc << "Maximum |cos A|, where A\n"
149  << "is the angle between edges\n"
150  << "at element center.\n"
151  << '\n'
152  << "Suggested ranges:\n"
153  << "Hexes: (0 -> 0.5)\n"
154  << "Quads: (0 -> 0.5)";
155  break;
156 
157  case SHEAR:
158  desc << "LIBMESH_DIM / K(Js)\n"
159  << '\n'
160  << "LIBMESH_DIM = element dimension.\n"
161  << "K(Js) = Condition number of \n"
162  << " Jacobian skew matrix.\n"
163  << '\n'
164  << "Suggested ranges:\n"
165  << "Hexes(LIBMESH_DIM=3): (0.3 -> 1)\n"
166  << "Quads(LIBMESH_DIM=2): (0.3 -> 1)";
167  break;
168 
169  case SHAPE:
170  desc << "LIBMESH_DIM / K(Jw)\n"
171  << '\n'
172  << "LIBMESH_DIM = element dimension.\n"
173  << "K(Jw) = Condition number of \n"
174  << " weighted Jacobian\n"
175  << " matrix.\n"
176  << '\n'
177  << "Suggested ranges:\n"
178  << "Hexes(LIBMESH_DIM=3): (0.3 -> 1)\n"
179  << "Tets(LIBMESH_DIM=3): (0.2 -> 1)\n"
180  << "Quads(LIBMESH_DIM=2): (0.3 -> 1).";
181  break;
182 
183  case MAX_ANGLE:
184  desc << "Largest included angle.\n"
185  << '\n'
186  << "Suggested ranges:\n"
187  << "Quads: (90 -> 135)\n"
188  << "Triangles: (60 -> 90)";
189  break;
190 
191  case MIN_ANGLE:
192  desc << "Smallest included angle.\n"
193  << '\n'
194  << "Suggested ranges:\n"
195  << "Quads: (45 -> 90)\n"
196  << "Triangles: (30 -> 60)";
197  break;
198 
199  case CONDITION:
200  desc << "Condition number of the\n"
201  << "Jacobian matrix.\n"
202  << '\n'
203  << "Suggested ranges:\n"
204  << "Quads: (1 -> 4)\n"
205  << "Hexes: (1 -> 8)\n"
206  << "Tris: (1 -> 1.3)\n"
207  << "Tets: (1 -> 3)";
208  break;
209 
210  case DISTORTION:
211  desc << "min |J| * A / <A>\n"
212  << '\n'
213  << "|J| = norm of Jacobian matrix\n"
214  << " A = actual area\n"
215  << "<A> = reference area\n"
216  << '\n'
217  << "Suggested ranges:\n"
218  << "Quads: (0.6 -> 1), <A>=4\n"
219  << "Hexes: (0.6 -> 1), <A>=8\n"
220  << "Tris: (0.6 -> 1), <A>=1/2\n"
221  << "Tets: (0.6 -> 1), <A>=1/6";
222  break;
223 
224  case TAPER:
225  desc << "Maximum ratio of lengths\n"
226  << "derived from opposite edges.\n"
227  << '\n'
228  << "Suggested ranges:\n"
229  << "Quads: (0.7 -> 1)\n"
230  << "Hexes: (0.4 -> 1)";
231  break;
232 
233  case WARP:
234  desc << "cos D\n"
235  << '\n'
236  << "D = minimum dihedral angle\n"
237  << " formed by diagonals.\n"
238  << '\n'
239  << "Suggested ranges:\n"
240  << "Quads: (0.9 -> 1)";
241  break;
242 
243  case STRETCH:
244  desc << "Sqrt(3) * L_min / L_max\n"
245  << '\n'
246  << "L_min = minimum edge length.\n"
247  << "L_max = maximum edge length.\n"
248  << '\n'
249  << "Suggested ranges:\n"
250  << "Quads: (0.25 -> 1)\n"
251  << "Hexes: (0.25 -> 1)";
252  break;
253 
254  case DIAGONAL:
255  desc << "D_min / D_max\n"
256  << '\n'
257  << "D_min = minimum diagonal.\n"
258  << "D_max = maximum diagonal.\n"
259  << '\n'
260  << "Suggested ranges:\n"
261  << "Hexes: (0.65 -> 1)";
262  break;
263 
264  case ASPECT_RATIO_BETA:
265  desc << "CR / (3 * IR)\n"
266  << '\n'
267  << "CR = circumsphere radius\n"
268  << "IR = inscribed sphere radius\n"
269  << '\n'
270  << "Suggested ranges:\n"
271  << "Tets: (1 -> 3)";
272  break;
273 
274  case ASPECT_RATIO_GAMMA:
275  desc << "S^(3/2) / 8.479670 * V\n"
276  << '\n'
277  << "S = sum(si*si/6)\n"
278  << "si = edge length\n"
279  << "V = volume\n"
280  << '\n'
281  << "Suggested ranges:\n"
282  << "Tets: (1 -> 3)";
283  break;
284 
285  case SIZE:
286  desc << "min (|J|, |1/J|)\n"
287  << '\n'
288  << "|J| = norm of Jacobian matrix.\n"
289  << '\n'
290  << "Suggested ranges:\n"
291  << "Quads: (0.3 -> 1)\n"
292  << "Hexes: (0.5 -> 1)\n"
293  << "Tris: (0.25 -> 1)\n"
294  << "Tets: (0.2 -> 1)";
295  break;
296 
297  case JACOBIAN:
298  desc << "Minimum Jacobian divided by\n"
299  << "the lengths of the LIBMESH_DIM\n"
300  << "largest edge vectors.\n"
301  << '\n'
302  << "LIBMESH_DIM = element dimension.\n"
303  << '\n'
304  << "Suggested ranges:\n"
305  << "Quads: (0.5 -> 1)\n"
306  << "Hexes: (0.5 -> 1)\n"
307  << "Tris: (0.5 -> 1.155)\n"
308  << "Tets: (0.5 -> 1.414)";
309  break;
310 
311  default:
312  desc << "Unknown";
313  break;
314  }
315 
316  return desc.str();
317 }
318 
319 
324 std::vector<ElemQuality> Quality::valid(const ElemType t)
325 {
326  std::vector<ElemQuality> v;
327 
328  switch (t)
329  {
330  case EDGE2:
331  case EDGE3:
332  case EDGE4:
333  {
334  // None yet
335  break;
336  }
337 
338  case TRI3:
339  case TRISHELL3:
340  case TRI6:
341  {
342  v.resize(7);
343  v[0] = MAX_ANGLE;
344  v[1] = MIN_ANGLE;
345  v[2] = CONDITION;
346  v[3] = JACOBIAN;
347  v[4] = SIZE;
348  v[5] = SHAPE;
349  v[6] = DISTORTION;
350  break;
351  }
352 
353  case QUAD4:
354  case QUADSHELL4:
355  case QUAD8:
356  case QUADSHELL8:
357  case QUAD9:
358  {
359  v.resize(13);
360  v[0] = ASPECT_RATIO;
361  v[1] = SKEW;
362  v[2] = TAPER;
363  v[3] = WARP;
364  v[4] = STRETCH;
365  v[5] = MIN_ANGLE;
366  v[6] = MAX_ANGLE;
367  v[7] = CONDITION;
368  v[8] = JACOBIAN;
369  v[9] = SHEAR;
370  v[10] = SHAPE;
371  v[11] = SIZE;
372  v[12] = DISTORTION;
373  break;
374  }
375 
376  case TET4:
377  case TET10:
378  {
379  v.resize(7);
380  v[0] = ASPECT_RATIO_BETA;
381  v[1] = ASPECT_RATIO_GAMMA;
382  v[2] = CONDITION;
383  v[3] = JACOBIAN;
384  v[4] = SHAPE;
385  v[5] = SIZE;
386  v[6] = DISTORTION;
387  break;
388  }
389 
390  case HEX8:
391  case HEX20:
392  case HEX27:
393  {
394  v.resize(11);
395  v[0] = ASPECT_RATIO;
396  v[1] = SKEW;
397  v[2] = SHEAR;
398  v[3] = SHAPE;
399  v[4] = CONDITION;
400  v[5] = JACOBIAN;
401  v[6] = DISTORTION;
402  v[7] = TAPER;
403  v[8] = STRETCH;
404  v[9] = DIAGONAL;
405  v[10] = SIZE;
406  break;
407  }
408 
409  case PRISM6:
410  case PRISM18:
411  {
412  // None yet
413  break;
414  }
415 
416  case PYRAMID5:
417  case PYRAMID13:
418  case PYRAMID14:
419  {
420  // None yet
421  break;
422  }
423 
424 
425 
426 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
427 
428  case INFEDGE2:
429  {
430  // None yet
431  break;
432  }
433 
434  case INFQUAD4:
435  case INFQUAD6:
436  {
437  // None yet
438  break;
439  }
440 
441  case INFHEX8:
442  case INFHEX16:
443  case INFHEX18:
444  {
445  // None yet
446  break;
447  }
448 
449  case INFPRISM6:
450  case INFPRISM12:
451  {
452  // None yet
453  break;
454  }
455 
456 #endif
457 
458 
459  default:
460  libmesh_error_msg("Undefined element type!");
461  }
462 
463  return v;
464 }
465 
466 } // namespace libMesh
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
std::string describe(const ElemQuality q)
Definition: elem_quality.C:130
std::vector< ElemQuality > valid(const ElemType t)
Definition: elem_quality.C:324