 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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 }
1.9.1.