Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
V_TetMetric.cpp
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Module: $RCSfile: V_TetMetric.cpp,v $
4 
5  Copyright (c) 2006 Sandia Corporation.
6  All rights reserved.
7  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
8 
9  This software is distributed WITHOUT ANY WARRANTY; without even
10  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11  PURPOSE. See the above copyright notice for more information.
12 
13 =========================================================================*/
14 
15 /*
16  *
17  * TetMetric.cpp contains quality calculations for Tets
18  *
19  * This file is part of VERDICT
20  *
21  */
22 
23 #define VERDICT_EXPORTS
24 
25 #include "moab/verdict.h"
26 #include "verdict_defines.hpp"
27 #include "VerdictVector.hpp"
28 #include "V_GaussIntegration.hpp"
29 #include <memory.h>
30 
31 //! the average volume of a tet
32 double verdict_tet_size = 0;
33 
34 /*!
35  set the average volume of a tet
36 */
37 C_FUNC_DEF void v_set_tet_size( double size )
38 {
39  verdict_tet_size = size;
40 }
41 
42 /*!
43  get the weights based on the average size
44  of a tet
45 */
46 int get_weight( VerdictVector& w1, VerdictVector& w2, VerdictVector& w3 )
47 {
48  static const double rt3 = sqrt( 3.0 );
49  static const double root_of_2 = sqrt( 2.0 );
50 
51  w1.set( 1, 0, 0 );
52  w2.set( 0.5, 0.5 * rt3, 0 );
53  w3.set( 0.5, rt3 / 6.0, root_of_2 / rt3 );
54 
55  double scale = pow( 6. * verdict_tet_size / determinant( w1, w2, w3 ), 0.3333333333333 );
56 
57  w1 *= scale;
58  w2 *= scale;
59  w3 *= scale;
60 
61  return 1;
62 }
63 
64 /*!
65  the edge ratio of a tet
66 
67  NB (P. Pebay 01/22/07):
68  Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
69  minimum edge lengths
70 */
71 C_FUNC_DEF double v_tet_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
72 {
73  VerdictVector a, b, c, d, e, f;
74 
75  a.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
76  coordinates[1][2] - coordinates[0][2] );
77 
78  b.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
79  coordinates[2][2] - coordinates[1][2] );
80 
81  c.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
82  coordinates[0][2] - coordinates[2][2] );
83 
84  d.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
85  coordinates[3][2] - coordinates[0][2] );
86 
87  e.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
88  coordinates[3][2] - coordinates[1][2] );
89 
90  f.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
91  coordinates[3][2] - coordinates[2][2] );
92 
93  double a2 = a.length_squared();
94  double b2 = b.length_squared();
95  double c2 = c.length_squared();
96  double d2 = d.length_squared();
97  double e2 = e.length_squared();
98  double f2 = f.length_squared();
99 
100  double m2, M2, mab, mcd, mef, Mab, Mcd, Mef;
101 
102  if( a2 < b2 )
103  {
104  mab = a2;
105  Mab = b2;
106  }
107  else // b2 <= a2
108  {
109  mab = b2;
110  Mab = a2;
111  }
112  if( c2 < d2 )
113  {
114  mcd = c2;
115  Mcd = d2;
116  }
117  else // d2 <= c2
118  {
119  mcd = d2;
120  Mcd = c2;
121  }
122  if( e2 < f2 )
123  {
124  mef = e2;
125  Mef = f2;
126  }
127  else // f2 <= e2
128  {
129  mef = f2;
130  Mef = e2;
131  }
132 
133  m2 = mab < mcd ? mab : mcd;
134  m2 = m2 < mef ? m2 : mef;
135 
136  if( m2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
137 
138  M2 = Mab > Mcd ? Mab : Mcd;
139  M2 = M2 > Mef ? M2 : Mef;
140 
141  double edge_ratio = sqrt( M2 / m2 );
142 
143  if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
144  return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
145 }
146 
147 /*!
148  the scaled jacobian of a tet
149 
150  minimum of the jacobian divided by the lengths of 3 edge vectors
151 
152 */
153 C_FUNC_DEF double v_tet_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
154 {
155 
156  VerdictVector side0, side1, side2, side3, side4, side5;
157 
158  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
159  coordinates[1][2] - coordinates[0][2] );
160 
161  side1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
162  coordinates[2][2] - coordinates[1][2] );
163 
164  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
165  coordinates[0][2] - coordinates[2][2] );
166 
167  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
168  coordinates[3][2] - coordinates[0][2] );
169 
170  side4.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
171  coordinates[3][2] - coordinates[1][2] );
172 
173  side5.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
174  coordinates[3][2] - coordinates[2][2] );
175 
176  double jacobi;
177 
178  jacobi = side3 % ( side2 * side0 );
179 
180  // products of lengths squared of each edge attached to a node.
181  double length_squared[4] = { side0.length_squared() * side2.length_squared() * side3.length_squared(),
182  side0.length_squared() * side1.length_squared() * side4.length_squared(),
183  side1.length_squared() * side2.length_squared() * side5.length_squared(),
184  side3.length_squared() * side4.length_squared() * side5.length_squared() };
185  int which_node = 0;
186  if( length_squared[1] > length_squared[which_node] ) which_node = 1;
187  if( length_squared[2] > length_squared[which_node] ) which_node = 2;
188  if( length_squared[3] > length_squared[which_node] ) which_node = 3;
189 
190  double length_product = sqrt( length_squared[which_node] );
191  if( length_product < fabs( jacobi ) ) length_product = fabs( jacobi );
192 
193  if( length_product < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
194 
195  static const double root_of_2 = sqrt( 2.0 );
196 
197  return (double)( root_of_2 * jacobi / length_product );
198 }
199 
200 /*!
201  The radius ratio of a tet
202 
203  NB (P. Pebay 04/16/07):
204  CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed
205  sphere radius.
206  Note that this metric is similar to the aspect beta of a tet, except that
207  it does not return VERDICT_DBL_MAX if the element has negative orientation.
208 */
209 C_FUNC_DEF double v_tet_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
210 {
211 
212  // Determine side vectors
213  VerdictVector side[6];
214 
215  side[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
216  coordinates[1][2] - coordinates[0][2] );
217 
218  side[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
219  coordinates[2][2] - coordinates[1][2] );
220 
221  side[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
222  coordinates[0][2] - coordinates[2][2] );
223 
224  side[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
225  coordinates[3][2] - coordinates[0][2] );
226 
227  side[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
228  coordinates[3][2] - coordinates[1][2] );
229 
230  side[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
231  coordinates[3][2] - coordinates[2][2] );
232 
233  VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
234  side[2].length_squared() * ( side[3] * side[0] ) +
235  side[0].length_squared() * ( side[3] * side[2] );
236 
237  double area_sum;
238  area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
239  ( side[3] * side[2] ).length() ) *
240  0.5;
241 
242  double volume = v_tet_volume( 4, coordinates );
243 
244  if( fabs( volume ) < VERDICT_DBL_MIN )
245  return (double)VERDICT_DBL_MAX;
246  else
247  {
248  double radius_ratio;
249  radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );
250 
251  return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
252  }
253 }
254 
255 /*!
256  The radius ratio of a positively-oriented tet, a.k.a. "aspect beta"
257 
258  NB (P. Pebay 04/16/07):
259  CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed
260  sphere radius if the element has positive orientation.
261  Note that this metric is similar to the radius ratio of a tet, except that
262  it returns VERDICT_DBL_MAX if the element has negative orientation.
263 
264 */
265 C_FUNC_DEF double v_tet_aspect_beta( int /*num_nodes*/, double coordinates[][3] )
266 {
267 
268  // Determine side vectors
269  VerdictVector side[6];
270 
271  side[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
272  coordinates[1][2] - coordinates[0][2] );
273 
274  side[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
275  coordinates[2][2] - coordinates[1][2] );
276 
277  side[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
278  coordinates[0][2] - coordinates[2][2] );
279 
280  side[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
281  coordinates[3][2] - coordinates[0][2] );
282 
283  side[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
284  coordinates[3][2] - coordinates[1][2] );
285 
286  side[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
287  coordinates[3][2] - coordinates[2][2] );
288 
289  VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
290  side[2].length_squared() * ( side[3] * side[0] ) +
291  side[0].length_squared() * ( side[3] * side[2] );
292 
293  double area_sum;
294  area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
295  ( side[3] * side[2] ).length() ) *
296  0.5;
297 
298  double volume = v_tet_volume( 4, coordinates );
299 
300  if( volume < VERDICT_DBL_MIN )
301  return (double)VERDICT_DBL_MAX;
302  else
303  {
304  double radius_ratio;
305  radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );
306 
307  return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
308  }
309 }
310 
311 /*!
312  The aspect ratio of a tet
313 
314  NB (P. Pebay 01/22/07):
315  Hmax / (2 sqrt(6) r) where Hmax and r respectively denote the greatest edge
316  length and the inradius of the tetrahedron
317 */
318 C_FUNC_DEF double v_tet_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
319 {
320  static const double normal_coeff = sqrt( 6. ) / 12.;
321 
322  // Determine side vectors
323  VerdictVector ab, bc, ac, ad, bd, cd;
324 
325  ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
326  coordinates[1][2] - coordinates[0][2] );
327 
328  ac.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
329  coordinates[2][2] - coordinates[0][2] );
330 
331  ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
332  coordinates[3][2] - coordinates[0][2] );
333 
334  double detTet = ab % ( ac * ad );
335 
336  if( detTet < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
337 
338  bc.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
339  coordinates[2][2] - coordinates[1][2] );
340 
341  bd.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
342  coordinates[3][2] - coordinates[1][2] );
343 
344  cd.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
345  coordinates[3][2] - coordinates[2][2] );
346 
347  double ab2 = ab.length_squared();
348  double bc2 = bc.length_squared();
349  double ac2 = ac.length_squared();
350  double ad2 = ad.length_squared();
351  double bd2 = bd.length_squared();
352  double cd2 = cd.length_squared();
353 
354  double A = ab2 > bc2 ? ab2 : bc2;
355  double B = ac2 > ad2 ? ac2 : ad2;
356  double C = bd2 > cd2 ? bd2 : cd2;
357  double D = A > B ? A : B;
358  double hm = D > C ? sqrt( D ) : sqrt( C );
359 
360  bd = ab * bc;
361  A = bd.length();
362  bd = ab * ad;
363  B = bd.length();
364  bd = ac * ad;
365  C = bd.length();
366  bd = bc * cd;
367  D = bd.length();
368 
369  double aspect_ratio;
370  aspect_ratio = normal_coeff * hm * ( A + B + C + D ) / fabs( detTet );
371 
372  if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
373  return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
374 }
375 
376 /*!
377  the aspect gamma of a tet
378 
379  srms^3 / (8.48528137423857*V) where srms = sqrt(sum(Si^2)/6), where Si is the edge length
380 */
381 C_FUNC_DEF double v_tet_aspect_gamma( int /*num_nodes*/, double coordinates[][3] )
382 {
383 
384  // Determine side vectors
385  VerdictVector side0, side1, side2, side3, side4, side5;
386 
387  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
388  coordinates[1][2] - coordinates[0][2] );
389 
390  side1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
391  coordinates[2][2] - coordinates[1][2] );
392 
393  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
394  coordinates[0][2] - coordinates[2][2] );
395 
396  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
397  coordinates[3][2] - coordinates[0][2] );
398 
399  side4.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
400  coordinates[3][2] - coordinates[1][2] );
401 
402  side5.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
403  coordinates[3][2] - coordinates[2][2] );
404 
405  double volume = fabs( v_tet_volume( 4, coordinates ) );
406 
407  if( volume < VERDICT_DBL_MIN )
408  return (double)VERDICT_DBL_MAX;
409  else
410  {
411  double srms = sqrt( ( side0.length_squared() + side1.length_squared() + side2.length_squared() +
412  side3.length_squared() + side4.length_squared() + side5.length_squared() ) /
413  6.0 );
414 
415  double aspect_ratio_gamma = pow( srms, 3 ) / ( 8.48528137423857 * volume );
416  return (double)aspect_ratio_gamma;
417  }
418 }
419 
420 /*!
421  The aspect frobenius of a tet
422 
423  NB (P. Pebay 01/22/07):
424  Frobenius condition number when the reference element is regular
425 */
426 C_FUNC_DEF double v_tet_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
427 {
428  static const double normal_exp = 1. / 3.;
429 
430  VerdictVector ab, ac, ad;
431 
432  ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
433  coordinates[1][2] - coordinates[0][2] );
434 
435  ac.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
436  coordinates[2][2] - coordinates[0][2] );
437 
438  ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
439  coordinates[3][2] - coordinates[0][2] );
440 
441  double denominator = ab % ( ac * ad );
442  denominator *= denominator;
443  denominator *= 2.;
444  denominator = 3. * pow( denominator, normal_exp );
445 
446  if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
447 
448  double u[3];
449  ab.get_xyz( u );
450  double v[3];
451  ac.get_xyz( v );
452  double w[3];
453  ad.get_xyz( w );
454 
455  double numerator = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
456  numerator += v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
457  numerator += w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
458  numerator *= 1.5;
459  numerator -= v[0] * u[0] + v[1] * u[1] + v[2] * u[2];
460  numerator -= w[0] * u[0] + w[1] * u[1] + w[2] * u[2];
461  numerator -= w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
462 
463  double aspect_frobenius = numerator / denominator;
464 
465  if( aspect_frobenius > 0 ) return (double)VERDICT_MIN( aspect_frobenius, VERDICT_DBL_MAX );
466  return (double)VERDICT_MAX( aspect_frobenius, -VERDICT_DBL_MAX );
467 }
468 
469 /*!
470  The minimum angle of a tet
471 
472  NB (P. Pebay 01/22/07):
473  minimum nonoriented dihedral angle
474 */
475 C_FUNC_DEF double v_tet_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
476 {
477  static const double normal_coeff = 180. * .3183098861837906715377675267450287;
478 
479  // Determine side vectors
480  VerdictVector ab, bc, ad, cd;
481 
482  ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
483  coordinates[1][2] - coordinates[0][2] );
484 
485  ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
486  coordinates[3][2] - coordinates[0][2] );
487 
488  bc.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
489  coordinates[2][2] - coordinates[1][2] );
490 
491  cd.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
492  coordinates[3][2] - coordinates[2][2] );
493 
494  VerdictVector abc = ab * bc;
495  double nabc = abc.length();
496  VerdictVector abd = ab * ad;
497  double nabd = abd.length();
498  VerdictVector acd = ad * cd;
499  double nacd = acd.length();
500  VerdictVector bcd = bc * cd;
501  double nbcd = bcd.length();
502 
503  double alpha = acos( ( abc % abd ) / ( nabc * nabd ) );
504  double beta = acos( ( abc % acd ) / ( nabc * nacd ) );
505  double gamma = acos( ( abc % bcd ) / ( nabc * nbcd ) );
506  double delta = acos( ( abd % acd ) / ( nabd * nacd ) );
507  double epsilon = acos( ( abd % bcd ) / ( nabd * nbcd ) );
508  double zeta = acos( ( acd % bcd ) / ( nacd * nbcd ) );
509 
510  alpha = alpha < beta ? alpha : beta;
511  alpha = alpha < gamma ? alpha : gamma;
512  alpha = alpha < delta ? alpha : delta;
513  alpha = alpha < epsilon ? alpha : epsilon;
514  alpha = alpha < zeta ? alpha : zeta;
515  alpha *= normal_coeff;
516 
517  if( alpha < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
518 
519  if( alpha > 0 ) return (double)VERDICT_MIN( alpha, VERDICT_DBL_MAX );
520  return (double)VERDICT_MAX( alpha, -VERDICT_DBL_MAX );
521 }
522 
523 /*!
524  The collapse ratio of a tet
525 */
526 C_FUNC_DEF double v_tet_collapse_ratio( int /*num_nodes*/, double coordinates[][3] )
527 {
528  // Determine side vectors
529  VerdictVector e01, e02, e03, e12, e13, e23;
530 
531  e01.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
532  coordinates[1][2] - coordinates[0][2] );
533 
534  e02.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
535  coordinates[2][2] - coordinates[0][2] );
536 
537  e03.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
538  coordinates[3][2] - coordinates[0][2] );
539 
540  e12.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
541  coordinates[2][2] - coordinates[1][2] );
542 
543  e13.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
544  coordinates[3][2] - coordinates[1][2] );
545 
546  e23.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
547  coordinates[3][2] - coordinates[2][2] );
548 
549  double l[6];
550  l[0] = e01.length();
551  l[1] = e02.length();
552  l[2] = e03.length();
553  l[3] = e12.length();
554  l[4] = e13.length();
555  l[5] = e23.length();
556 
557  // Find longest edge for each bounding triangle of tetrahedron
558  double l012 = l[4] > l[0] ? l[4] : l[0];
559  l012 = l[1] > l012 ? l[1] : l012;
560  double l031 = l[0] > l[2] ? l[0] : l[2];
561  l031 = l[3] > l031 ? l[3] : l031;
562  double l023 = l[2] > l[1] ? l[2] : l[1];
563  l023 = l[5] > l023 ? l[5] : l023;
564  double l132 = l[4] > l[3] ? l[4] : l[3];
565  l132 = l[5] > l132 ? l[5] : l132;
566 
567  // Compute collapse ratio for each vertex/triangle pair
568  VerdictVector N;
569  double h, magN;
570  double cr;
571  double crMin;
572 
573  N = e01 * e02;
574  magN = N.length();
575  h = ( e03 % N ) / magN; // height of vertex 3 above 0-1-2
576  crMin = h / l012; // ratio of height to longest edge of 0-1-2
577 
578  N = e03 * e01;
579  magN = N.length();
580  h = ( e02 % N ) / magN; // height of vertex 2 above 0-3-1
581  cr = h / l031; // ratio of height to longest edge of 0-3-1
582  if( cr < crMin ) crMin = cr;
583 
584  N = e02 * e03;
585  magN = N.length();
586  h = ( e01 % N ) / magN; // height of vertex 1 above 0-2-3
587  cr = h / l023; // ratio of height to longest edge of 0-2-3
588  if( cr < crMin ) crMin = cr;
589 
590  N = e12 * e13;
591  magN = N.length();
592  h = ( e01 % N ) / magN; // height of vertex 0 above 1-3-2
593  cr = h / l132; // ratio of height to longest edge of 1-3-2
594  if( cr < crMin ) crMin = cr;
595 
596  if( crMin < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
597  if( crMin > 0 ) return (double)VERDICT_MIN( crMin, VERDICT_DBL_MAX );
598  return (double)VERDICT_MAX( crMin, -VERDICT_DBL_MAX );
599 }
600 
601 /*!
602  the volume of a tet
603 
604  1/6 * jacobian at a corner node
605 */
606 C_FUNC_DEF double v_tet_volume( int /*num_nodes*/, double coordinates[][3] )
607 {
608 
609  // Determine side vectors
610  VerdictVector side0, side2, side3;
611 
612  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
613  coordinates[1][2] - coordinates[0][2] );
614 
615  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
616  coordinates[0][2] - coordinates[2][2] );
617 
618  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
619  coordinates[3][2] - coordinates[0][2] );
620 
621  return (double)( ( side3 % ( side2 * side0 ) ) / 6.0 );
622 }
623 
624 /*!
625  the condition of a tet
626 
627  condition number of the jacobian matrix at any corner
628 */
629 C_FUNC_DEF double v_tet_condition( int /*num_nodes*/, double coordinates[][3] )
630 {
631 
632  double condition, term1, term2, det;
633  double rt3 = sqrt( 3.0 );
634  double rt6 = sqrt( 6.0 );
635 
636  VerdictVector side0, side2, side3;
637 
638  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
639  coordinates[1][2] - coordinates[0][2] );
640 
641  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
642  coordinates[0][2] - coordinates[2][2] );
643 
644  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
645  coordinates[3][2] - coordinates[0][2] );
646 
647  VerdictVector c_1, c_2, c_3;
648 
649  c_1 = side0;
650  c_2 = ( -2 * side2 - side0 ) / rt3;
651  c_3 = ( 3 * side3 + side2 - side0 ) / rt6;
652 
653  term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
654  term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) + ( c_2 * c_3 ) % ( c_2 * c_3 ) + ( c_1 * c_3 ) % ( c_1 * c_3 );
655  det = c_1 % ( c_2 * c_3 );
656 
657  if( det <= VERDICT_DBL_MIN )
658  return VERDICT_DBL_MAX;
659  else
660  condition = sqrt( term1 * term2 ) / ( 3.0 * det );
661 
662  return (double)condition;
663 }
664 
665 /*!
666  the jacobian of a tet
667 
668  TODO
669 */
670 C_FUNC_DEF double v_tet_jacobian( int /*num_nodes*/, double coordinates[][3] )
671 {
672  VerdictVector side0, side2, side3;
673 
674  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
675  coordinates[1][2] - coordinates[0][2] );
676 
677  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
678  coordinates[0][2] - coordinates[2][2] );
679 
680  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
681  coordinates[3][2] - coordinates[0][2] );
682 
683  return (double)( side3 % ( side2 * side0 ) );
684 }
685 
686 /*!
687  the shape of a tet
688 
689  3/ condition number of weighted jacobian matrix
690 */
691 C_FUNC_DEF double v_tet_shape( int /*num_nodes*/, double coordinates[][3] )
692 {
693 
694  static const double two_thirds = 2.0 / 3.0;
695  static const double root_of_2 = sqrt( 2.0 );
696 
697  VerdictVector edge0, edge2, edge3;
698 
699  edge0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
700  coordinates[1][2] - coordinates[0][2] );
701 
702  edge2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
703  coordinates[0][2] - coordinates[2][2] );
704 
705  edge3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
706  coordinates[3][2] - coordinates[0][2] );
707 
708  double jacobian = edge3 % ( edge2 * edge0 );
709  if( jacobian < VERDICT_DBL_MIN )
710  {
711  return (double)0.0;
712  }
713  double num = 3 * pow( root_of_2 * jacobian, two_thirds );
714  double den =
715  1.5 * ( edge0 % edge0 + edge2 % edge2 + edge3 % edge3 ) - ( edge0 % -edge2 + -edge2 % edge3 + edge3 % edge0 );
716 
717  if( den < VERDICT_DBL_MIN ) return (double)0.0;
718 
719  return (double)VERDICT_MAX( num / den, 0 );
720 }
721 
722 /*!
723  the relative size of a tet
724 
725  Min(J,1/J), where J is the determinant of the weighted Jacobian matrix
726 */
727 C_FUNC_DEF double v_tet_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
728 {
729  double size;
730  VerdictVector w1, w2, w3;
731  get_weight( w1, w2, w3 );
732  double avg_volume = ( w1 % ( w2 * w3 ) ) / 6.0;
733 
734  double volume = v_tet_volume( 4, coordinates );
735 
736  if( avg_volume < VERDICT_DBL_MIN )
737  return 0.0;
738  else
739  {
740  size = volume / avg_volume;
741  if( size <= VERDICT_DBL_MIN ) return 0.0;
742  if( size > 1 ) size = (double)( 1 ) / size;
743  }
744  return (double)( size * size );
745 }
746 
747 /*!
748  the shape and size of a tet
749 
750  Product of the shape and relative size
751 */
752 C_FUNC_DEF double v_tet_shape_and_size( int num_nodes, double coordinates[][3] )
753 {
754 
755  double shape, size;
756  shape = v_tet_shape( num_nodes, coordinates );
757  size = v_tet_relative_size_squared( num_nodes, coordinates );
758 
759  return (double)( shape * size );
760 }
761 
762 /*!
763  the distortion of a tet
764 */
765 C_FUNC_DEF double v_tet_distortion( int num_nodes, double coordinates[][3] )
766 {
767 
768  double distortion = VERDICT_DBL_MAX;
769  int number_of_gauss_points = 0;
770  if( num_nodes == 4 )
771  // for linear tet, the distortion is always 1 because
772  // straight edge tets are the target shape for tet
773  return 1.0;
774 
775  else if( num_nodes == 10 )
776  // use four integration points for quadratic tet
777  number_of_gauss_points = 4;
778 
779  int number_dims = 3;
780  int total_number_of_gauss_points = number_of_gauss_points;
781  // use is_tri=1 to indicate this is for tet in 3D
782  int is_tri = 1;
783 
784  double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
785  double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
786  double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
787  double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
788  double weight[maxTotalNumberGaussPoints];
789 
790  // create an object of GaussIntegration for tet
791  GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
792  GaussIntegration::calculate_shape_function_3d_tet();
793  GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight );
794 
795  // vector xxi is the derivative vector of coordinates w.r.t local xi coordinate in the
796  // computation space
797  // vector xet is the derivative vector of coordinates w.r.t local et coordinate in the
798  // computation space
799  // vector xze is the derivative vector of coordinates w.r.t local ze coordinate in the
800  // computation space
801  VerdictVector xxi, xet, xze, xin;
802 
803  double jacobian, minimum_jacobian;
804  double element_volume = 0.0;
805  minimum_jacobian = VERDICT_DBL_MAX;
806 
807  // calculate element volume
808  int ife, ja;
809  for( ife = 0; ife < total_number_of_gauss_points; ife++ )
810  {
811  xxi.set( 0.0, 0.0, 0.0 );
812  xet.set( 0.0, 0.0, 0.0 );
813  xze.set( 0.0, 0.0, 0.0 );
814 
815  for( ja = 0; ja < num_nodes; ja++ )
816  {
817  xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
818  xxi += dndy1[ife][ja] * xin;
819  xet += dndy2[ife][ja] * xin;
820  xze += dndy3[ife][ja] * xin;
821  }
822 
823  // determinant
824  jacobian = xxi % ( xet * xze );
825  if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
826 
827  element_volume += weight[ife] * jacobian;
828  } // element_volume is 6 times the actual volume
829 
830  // loop through all nodes
831  double dndy1_at_node[maxNumberNodes][maxNumberNodes];
832  double dndy2_at_node[maxNumberNodes][maxNumberNodes];
833  double dndy3_at_node[maxNumberNodes][maxNumberNodes];
834 
835  GaussIntegration::calculate_derivative_at_nodes_3d_tet( dndy1_at_node, dndy2_at_node, dndy3_at_node );
836  int node_id;
837  for( node_id = 0; node_id < num_nodes; node_id++ )
838  {
839  xxi.set( 0.0, 0.0, 0.0 );
840  xet.set( 0.0, 0.0, 0.0 );
841  xze.set( 0.0, 0.0, 0.0 );
842 
843  for( ja = 0; ja < num_nodes; ja++ )
844  {
845  xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
846  xxi += dndy1_at_node[node_id][ja] * xin;
847  xet += dndy2_at_node[node_id][ja] * xin;
848  xze += dndy3_at_node[node_id][ja] * xin;
849  }
850 
851  jacobian = xxi % ( xet * xze );
852  if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
853  }
854  distortion = minimum_jacobian / element_volume;
855 
856  return (double)distortion;
857 }
858 
859 /*!
860  the quality metrics of a tet
861 */
862 C_FUNC_DEF void v_tet_quality( int num_nodes,
863  double coordinates[][3],
864  unsigned int metrics_request_flag,
865  TetMetricVals* metric_vals )
866 {
867 
868  memset( metric_vals, 0, sizeof( TetMetricVals ) );
869 
870  /*
871 
872  node numbers and edge numbers below
873 
874 
875 
876  3
877  + edge 0 is node 0 to 1
878  +|+ edge 1 is node 1 to 2
879  3/ | \5 edge 2 is node 0 to 2
880  / 4| \ edge 3 is node 0 to 3
881  0 - -|- + 2 edge 4 is node 1 to 3
882  \ | + edge 5 is node 2 to 3
883  0\ | /1
884  +|/ edge 2 is behind edge 4
885  1
886 
887 
888  */
889 
890  // lets start with making the vectors
891  VerdictVector edges[6];
892  edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
893  coordinates[1][2] - coordinates[0][2] );
894 
895  edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
896  coordinates[2][2] - coordinates[1][2] );
897 
898  edges[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
899  coordinates[0][2] - coordinates[2][2] );
900 
901  edges[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
902  coordinates[3][2] - coordinates[0][2] );
903 
904  edges[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
905  coordinates[3][2] - coordinates[1][2] );
906 
907  edges[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
908  coordinates[3][2] - coordinates[2][2] );
909 
910  // common numbers
911  static const double root_of_2 = sqrt( 2.0 );
912 
913  // calculate the jacobian
914  static const int do_jacobian = V_TET_JACOBIAN | V_TET_VOLUME | V_TET_ASPECT_BETA | V_TET_ASPECT_GAMMA |
915  V_TET_SHAPE | V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE |
916  V_TET_SCALED_JACOBIAN | V_TET_CONDITION;
917  if( metrics_request_flag & do_jacobian )
918  {
919  metric_vals->jacobian = (double)( edges[3] % ( edges[2] * edges[0] ) );
920  }
921 
922  // calculate the volume
923  if( metrics_request_flag & V_TET_VOLUME )
924  {
925  metric_vals->volume = (double)( metric_vals->jacobian / 6.0 );
926  }
927 
928  // calculate aspect ratio
929  if( metrics_request_flag & V_TET_ASPECT_BETA )
930  {
931  double surface_area = ( ( edges[2] * edges[0] ).length() + ( edges[3] * edges[0] ).length() +
932  ( edges[4] * edges[1] ).length() + ( edges[3] * edges[2] ).length() ) *
933  0.5;
934 
935  VerdictVector numerator = edges[3].length_squared() * ( edges[2] * edges[0] ) +
936  edges[2].length_squared() * ( edges[3] * edges[0] ) +
937  edges[0].length_squared() * ( edges[3] * edges[2] );
938 
939  double volume = metric_vals->jacobian / 6.0;
940 
941  if( volume < VERDICT_DBL_MIN )
942  metric_vals->aspect_beta = (double)( VERDICT_DBL_MAX );
943  else
944  metric_vals->aspect_beta = (double)( numerator.length() * surface_area / ( 108 * volume * volume ) );
945  }
946 
947  // calculate the aspect gamma
948  if( metrics_request_flag & V_TET_ASPECT_GAMMA )
949  {
950  double volume = fabs( metric_vals->jacobian / 6.0 );
951  if( fabs( volume ) < VERDICT_DBL_MIN )
952  metric_vals->aspect_gamma = VERDICT_DBL_MAX;
953  else
954  {
955  double srms = sqrt( ( edges[0].length_squared() + edges[1].length_squared() + edges[2].length_squared() +
956  edges[3].length_squared() + edges[4].length_squared() + edges[5].length_squared() ) /
957  6.0 );
958 
959  // cube the srms
960  srms *= ( srms * srms );
961  metric_vals->aspect_gamma = (double)( srms / ( 8.48528137423857 * volume ) );
962  }
963  }
964 
965  // calculate the shape of the tet
966  if( metrics_request_flag & ( V_TET_SHAPE | V_TET_SHAPE_AND_SIZE ) )
967  {
968  // if the jacobian is non-positive, the shape is 0
969  if( metric_vals->jacobian < VERDICT_DBL_MIN )
970  {
971  metric_vals->shape = (double)0.0;
972  }
973  else
974  {
975  static const double two_thirds = 2.0 / 3.0;
976  double num = 3.0 * pow( root_of_2 * metric_vals->jacobian, two_thirds );
977  double den = 1.5 * ( edges[0] % edges[0] + edges[2] % edges[2] + edges[3] % edges[3] ) -
978  ( edges[0] % -edges[2] + -edges[2] % edges[3] + edges[3] % edges[0] );
979 
980  if( den < VERDICT_DBL_MIN )
981  metric_vals->shape = (double)0.0;
982  else
983  metric_vals->shape = (double)VERDICT_MAX( num / den, 0 );
984  }
985  }
986 
987  // calculate the relative size of the tet
988  if( metrics_request_flag & ( V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE ) )
989  {
990  VerdictVector w1, w2, w3;
991  get_weight( w1, w2, w3 );
992  double avg_vol = ( w1 % ( w2 * w3 ) ) / 6;
993 
994  if( avg_vol < VERDICT_DBL_MIN )
995  metric_vals->relative_size_squared = 0.0;
996  else
997  {
998  double tmp = metric_vals->jacobian / ( 6 * avg_vol );
999  if( tmp < VERDICT_DBL_MIN )
1000  metric_vals->relative_size_squared = 0.0;
1001  else
1002  {
1003  tmp *= tmp;
1004  metric_vals->relative_size_squared = (double)VERDICT_MIN( tmp, 1 / tmp );
1005  }
1006  }
1007  }
1008 
1009  // calculate the shape and size
1010  if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
1011  {
1012  metric_vals->shape_and_size = (double)( metric_vals->shape * metric_vals->relative_size_squared );
1013  }
1014 
1015  // calculate the scaled jacobian
1016  if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
1017  {
1018  // find out which node the normalized jacobian can be calculated at
1019  // and it will be the smaller than at other nodes
1020  double length_squared[4] = { edges[0].length_squared() * edges[2].length_squared() * edges[3].length_squared(),
1021  edges[0].length_squared() * edges[1].length_squared() * edges[4].length_squared(),
1022  edges[1].length_squared() * edges[2].length_squared() * edges[5].length_squared(),
1023  edges[3].length_squared() * edges[4].length_squared() *
1024  edges[5].length_squared() };
1025 
1026  int which_node = 0;
1027  if( length_squared[1] > length_squared[which_node] ) which_node = 1;
1028  if( length_squared[2] > length_squared[which_node] ) which_node = 2;
1029  if( length_squared[3] > length_squared[which_node] ) which_node = 3;
1030 
1031  // find the scaled jacobian at this node
1032  double length_product = sqrt( length_squared[which_node] );
1033  if( length_product < fabs( metric_vals->jacobian ) ) length_product = fabs( metric_vals->jacobian );
1034 
1035  if( length_product < VERDICT_DBL_MIN )
1036  metric_vals->scaled_jacobian = (double)VERDICT_DBL_MAX;
1037  else
1038  metric_vals->scaled_jacobian = (double)( root_of_2 * metric_vals->jacobian / length_product );
1039  }
1040 
1041  // calculate the condition number
1042  if( metrics_request_flag & V_TET_CONDITION )
1043  {
1044  static const double root_of_3 = sqrt( 3.0 );
1045  static const double root_of_6 = sqrt( 6.0 );
1046 
1047  VerdictVector c_1, c_2, c_3;
1048  c_1 = edges[0];
1049  c_2 = ( -2 * edges[2] - edges[0] ) / root_of_3;
1050  c_3 = ( 3 * edges[3] + edges[2] - edges[0] ) / root_of_6;
1051 
1052  double term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
1053  double term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) + ( c_2 * c_3 ) % ( c_2 * c_3 ) + ( c_3 * c_1 ) % ( c_3 * c_1 );
1054 
1055  double det = c_1 % ( c_2 * c_3 );
1056 
1057  if( det <= VERDICT_DBL_MIN )
1058  metric_vals->condition = (double)VERDICT_DBL_MAX;
1059  else
1060  metric_vals->condition = (double)( sqrt( term1 * term2 ) / ( 3.0 * det ) );
1061  }
1062 
1063  // calculate the distortion
1064  if( metrics_request_flag & V_TET_DISTORTION )
1065  {
1066  metric_vals->distortion = v_tet_distortion( num_nodes, coordinates );
1067  }
1068 
1069  // check for overflow
1070  if( metrics_request_flag & V_TET_ASPECT_BETA )
1071  {
1072  if( metric_vals->aspect_beta > 0 )
1073  metric_vals->aspect_beta = (double)VERDICT_MIN( metric_vals->aspect_beta, VERDICT_DBL_MAX );
1074  metric_vals->aspect_beta = (double)VERDICT_MAX( metric_vals->aspect_beta, -VERDICT_DBL_MAX );
1075  }
1076 
1077  if( metrics_request_flag & V_TET_ASPECT_GAMMA )
1078  {
1079  if( metric_vals->aspect_gamma > 0 )
1080  metric_vals->aspect_gamma = (double)VERDICT_MIN( metric_vals->aspect_gamma, VERDICT_DBL_MAX );
1081  metric_vals->aspect_gamma = (double)VERDICT_MAX( metric_vals->aspect_gamma, -VERDICT_DBL_MAX );
1082  }
1083 
1084  if( metrics_request_flag & V_TET_VOLUME )
1085  {
1086  if( metric_vals->volume > 0 ) metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
1087  metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
1088  }
1089 
1090  if( metrics_request_flag & V_TET_CONDITION )
1091  {
1092  if( metric_vals->condition > 0 )
1093  metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
1094  metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
1095  }
1096 
1097  if( metrics_request_flag & V_TET_JACOBIAN )
1098  {
1099  if( metric_vals->jacobian > 0 )
1100  metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
1101  metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
1102  }
1103 
1104  if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
1105  {
1106  if( metric_vals->scaled_jacobian > 0 )
1107  metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
1108  metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
1109  }
1110 
1111  if( metrics_request_flag & V_TET_SHAPE )
1112  {
1113  if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
1114  metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
1115  }
1116 
1117  if( metrics_request_flag & V_TET_RELATIVE_SIZE_SQUARED )
1118  {
1119  if( metric_vals->relative_size_squared > 0 )
1120  metric_vals->relative_size_squared =
1121  (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
1122  metric_vals->relative_size_squared =
1123  (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
1124  }
1125 
1126  if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
1127  {
1128  if( metric_vals->shape_and_size > 0 )
1129  metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
1130  metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
1131  }
1132 
1133  if( metrics_request_flag & V_TET_DISTORTION )
1134  {
1135  if( metric_vals->distortion > 0 )
1136  metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
1137  metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
1138  }
1139 
1140  if( metrics_request_flag & V_TET_ASPECT_RATIO ) metric_vals->aspect_ratio = v_tet_aspect_ratio( 4, coordinates );
1141 
1142  if( metrics_request_flag & V_TET_ASPECT_FROBENIUS )
1143  metric_vals->aspect_frobenius = v_tet_aspect_frobenius( 4, coordinates );
1144 
1145  if( metrics_request_flag & V_TET_MINIMUM_ANGLE ) metric_vals->minimum_angle = v_tet_minimum_angle( 4, coordinates );
1146 
1147  if( metrics_request_flag & V_TET_COLLAPSE_RATIO )
1148  metric_vals->collapse_ratio = v_tet_collapse_ratio( 4, coordinates );
1149 
1150  if( metrics_request_flag & V_TET_RADIUS_RATIO ) metric_vals->radius_ratio = v_tet_radius_ratio( 4, coordinates );
1151 }