1 /*=========================================================================
2
3 Module: V_TetMetric.cpp
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 *
18 * TetMetric.cpp contains quality calculations for Tets
19 *
20 * This file is part of VERDICT
21 *
22 */
23
24
25 #include "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 static double v_tet_size = 0;
33
34 /*!
35 set the average volume of a tet
36 */
v_set_tet_size(double size)37 C_FUNC_DEF void v_set_tet_size( double size )
38 {
39 v_tet_size = size;
40 }
41
42 /*!
43 get the weights based on the average size
44 of a tet
45 */
v_tet_get_weight(VerdictVector & w1,VerdictVector & w2,VerdictVector & w3)46 static int v_tet_get_weight (
47 VerdictVector &w1, VerdictVector &w2, VerdictVector &w3 )
48 {
49 static const double rt3 = sqrt(3.0);
50 static const double root_of_2 = sqrt(2.0);
51
52 w1.set(1,0,0);
53 w2.set(0.5, 0.5*rt3, 0 );
54 w3.set(0.5, rt3/6.0, root_of_2/rt3);
55
56 double scale = pow( 6.*v_tet_size/v_determinant(w1,w2,w3),0.3333333333333);
57
58 w1 *= scale;
59 w2 *= scale;
60 w3 *= scale;
61
62 return 1;
63 }
64
65 /*!
66 the edge ratio of a tet
67
68 NB (P. Pebay 01/22/07):
69 Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
70 minimum edge lengths
71 */
v_tet_edge_ratio(int,double coordinates[][3])72 C_FUNC_DEF double v_tet_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
73 {
74 VerdictVector a, b, c, d, e, f;
75
76 a.set( coordinates[1][0] - coordinates[0][0],
77 coordinates[1][1] - coordinates[0][1],
78 coordinates[1][2] - coordinates[0][2] );
79
80 b.set( coordinates[2][0] - coordinates[1][0],
81 coordinates[2][1] - coordinates[1][1],
82 coordinates[2][2] - coordinates[1][2] );
83
84 c.set( coordinates[0][0] - coordinates[2][0],
85 coordinates[0][1] - coordinates[2][1],
86 coordinates[0][2] - coordinates[2][2] );
87
88 d.set( coordinates[3][0] - coordinates[0][0],
89 coordinates[3][1] - coordinates[0][1],
90 coordinates[3][2] - coordinates[0][2] );
91
92 e.set( coordinates[3][0] - coordinates[1][0],
93 coordinates[3][1] - coordinates[1][1],
94 coordinates[3][2] - coordinates[1][2] );
95
96 f.set( coordinates[3][0] - coordinates[2][0],
97 coordinates[3][1] - coordinates[2][1],
98 coordinates[3][2] - coordinates[2][2] );
99
100 double a2 = a.length_squared();
101 double b2 = b.length_squared();
102 double c2 = c.length_squared();
103 double d2 = d.length_squared();
104 double e2 = e.length_squared();
105 double f2 = f.length_squared();
106
107 double m2,M2,mab,mcd,mef,Mab,Mcd,Mef;
108
109 if ( a2 < b2 )
110 {
111 mab = a2;
112 Mab = b2;
113 }
114 else // b2 <= a2
115 {
116 mab = b2;
117 Mab = a2;
118 }
119 if ( c2 < d2 )
120 {
121 mcd = c2;
122 Mcd = d2;
123 }
124 else // d2 <= c2
125 {
126 mcd = d2;
127 Mcd = c2;
128 }
129 if ( e2 < f2 )
130 {
131 mef = e2;
132 Mef = f2;
133 }
134 else // f2 <= e2
135 {
136 mef = f2;
137 Mef = e2;
138 }
139
140 m2 = mab < mcd ? mab : mcd;
141 m2 = m2 < mef ? m2 : mef;
142
143 if( m2 < VERDICT_DBL_MIN )
144 return (double)VERDICT_DBL_MAX;
145
146 M2 = Mab > Mcd ? Mab : Mcd;
147 M2 = M2 > Mef ? M2 : Mef;
148
149 double edge_ratio = sqrt( M2 / m2 );
150
151 if( edge_ratio > 0 )
152 return (double) VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
153 return (double) VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
154 }
155
156 /*!
157 the scaled jacobian of a tet
158
159 minimum of the jacobian divided by the lengths of 3 edge vectors
160
161 */
v_tet_scaled_jacobian(int,double coordinates[][3])162 C_FUNC_DEF double v_tet_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
163 {
164
165 VerdictVector side0, side1, side2, side3, side4, side5;
166
167 side0.set( coordinates[1][0] - coordinates[0][0],
168 coordinates[1][1] - coordinates[0][1],
169 coordinates[1][2] - coordinates[0][2] );
170
171 side1.set( coordinates[2][0] - coordinates[1][0],
172 coordinates[2][1] - coordinates[1][1],
173 coordinates[2][2] - coordinates[1][2] );
174
175 side2.set( coordinates[0][0] - coordinates[2][0],
176 coordinates[0][1] - coordinates[2][1],
177 coordinates[0][2] - coordinates[2][2] );
178
179 side3.set( coordinates[3][0] - coordinates[0][0],
180 coordinates[3][1] - coordinates[0][1],
181 coordinates[3][2] - coordinates[0][2] );
182
183 side4.set( coordinates[3][0] - coordinates[1][0],
184 coordinates[3][1] - coordinates[1][1],
185 coordinates[3][2] - coordinates[1][2] );
186
187 side5.set( coordinates[3][0] - coordinates[2][0],
188 coordinates[3][1] - coordinates[2][1],
189 coordinates[3][2] - coordinates[2][2] );
190
191 double jacobi;
192
193 jacobi = side3 % ( side2 * side0 );
194
195 // products of lengths squared of each edge attached to a node.
196 double length_squared[4] = {
197 side0.length_squared() * side2.length_squared() * side3.length_squared(),
198 side0.length_squared() * side1.length_squared() * side4.length_squared(),
199 side1.length_squared() * side2.length_squared() * side5.length_squared(),
200 side3.length_squared() * side4.length_squared() * side5.length_squared()
201 };
202 int which_node = 0;
203 if(length_squared[1] > length_squared[which_node])
204 which_node = 1;
205 if(length_squared[2] > length_squared[which_node])
206 which_node = 2;
207 if(length_squared[3] > length_squared[which_node])
208 which_node = 3;
209
210 double length_product = sqrt( length_squared[which_node] );
211 if(length_product < fabs(jacobi))
212 length_product = fabs(jacobi);
213
214 if( length_product < VERDICT_DBL_MIN )
215 return (double) VERDICT_DBL_MAX;
216
217 static const double root_of_2 = sqrt(2.0);
218
219 return (double)(root_of_2 * jacobi / length_product);
220
221 }
222
223 /*!
224 The radius ratio of a tet
225
226 NB (P. Pebay 04/16/07):
227 CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed
228 sphere radius.
229 Note that this function is similar to the aspect beta of a tet, except that
230 it does not return VERDICT_DBL_MAX if the element has negative orientation.
231 */
v_tet_radius_ratio(int,double coordinates[][3])232 C_FUNC_DEF double v_tet_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
233 {
234
235 //Determine side vectors
236 VerdictVector side[6];
237
238 side[0].set( coordinates[1][0] - coordinates[0][0],
239 coordinates[1][1] - coordinates[0][1],
240 coordinates[1][2] - coordinates[0][2] );
241
242 side[1].set( coordinates[2][0] - coordinates[1][0],
243 coordinates[2][1] - coordinates[1][1],
244 coordinates[2][2] - coordinates[1][2] );
245
246 side[2].set( coordinates[0][0] - coordinates[2][0],
247 coordinates[0][1] - coordinates[2][1],
248 coordinates[0][2] - coordinates[2][2] );
249
250 side[3].set( coordinates[3][0] - coordinates[0][0],
251 coordinates[3][1] - coordinates[0][1],
252 coordinates[3][2] - coordinates[0][2] );
253
254 side[4].set( coordinates[3][0] - coordinates[1][0],
255 coordinates[3][1] - coordinates[1][1],
256 coordinates[3][2] - coordinates[1][2] );
257
258 side[5].set( coordinates[3][0] - coordinates[2][0],
259 coordinates[3][1] - coordinates[2][1],
260 coordinates[3][2] - coordinates[2][2] );
261
262 VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0]) +
263 side[2].length_squared() * ( side[3] * side[0]) +
264 side[0].length_squared() * ( side[3] * side[2]);
265
266 double area_sum;
267 area_sum = ((side[2] * side[0]).length() +
268 (side[3] * side[0]).length() +
269 (side[4] * side[1]).length() +
270 (side[3] * side[2]).length() ) * 0.5;
271
272 double volume = v_tet_volume(4, coordinates);
273
274 if( fabs( volume ) < VERDICT_DBL_MIN )
275 return (double)VERDICT_DBL_MAX;
276 else
277 {
278 double radius_ratio;
279 radius_ratio = numerator.length() * area_sum / (108 * volume * volume);
280
281 return (double) VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
282 }
283
284 }
285
286 /*!
287 The radius ratio of a positively-oriented tet, a.k.a. "aspect beta"
288
289 NB (P. Pebay 04/16/07):
290 CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed
291 sphere radius if the element has positive orientation.
292 Note that this function is similar to the radius ratio of a tet, except that
293 it returns VERDICT_DBL_MAX if the element has negative orientation.
294 NB (J. Pouderoux 01/27/15)
295 This will return VERDICT_DBL_MAX when the volume of the tetrahedron is ill-
296 conditioned. Previously, this would only happen when the volume was small
297 and positive, but now ill-conditioned inverted tetrahedra are also included.
298
299 */
v_tet_aspect_beta(int,double coordinates[][3])300 C_FUNC_DEF double v_tet_aspect_beta( int /*num_nodes*/, double coordinates[][3] )
301 {
302
303 //Determine side vectors
304 VerdictVector side[6];
305
306 side[0].set( coordinates[1][0] - coordinates[0][0],
307 coordinates[1][1] - coordinates[0][1],
308 coordinates[1][2] - coordinates[0][2] );
309
310 side[1].set( coordinates[2][0] - coordinates[1][0],
311 coordinates[2][1] - coordinates[1][1],
312 coordinates[2][2] - coordinates[1][2] );
313
314 side[2].set( coordinates[0][0] - coordinates[2][0],
315 coordinates[0][1] - coordinates[2][1],
316 coordinates[0][2] - coordinates[2][2] );
317
318 side[3].set( coordinates[3][0] - coordinates[0][0],
319 coordinates[3][1] - coordinates[0][1],
320 coordinates[3][2] - coordinates[0][2] );
321
322 side[4].set( coordinates[3][0] - coordinates[1][0],
323 coordinates[3][1] - coordinates[1][1],
324 coordinates[3][2] - coordinates[1][2] );
325
326 side[5].set( coordinates[3][0] - coordinates[2][0],
327 coordinates[3][1] - coordinates[2][1],
328 coordinates[3][2] - coordinates[2][2] );
329
330 VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0]) +
331 side[2].length_squared() * ( side[3] * side[0]) +
332 side[0].length_squared() * ( side[3] * side[2]);
333
334 double area_sum;
335 area_sum = ((side[2] * side[0]).length() +
336 (side[3] * side[0]).length() +
337 (side[4] * side[1]).length() +
338 (side[3] * side[2]).length() ) * 0.5;
339
340 double volume = v_tet_volume(4, coordinates);
341
342 if( fabs( volume ) < VERDICT_DBL_MIN )
343 return (double)VERDICT_DBL_MAX;
344 else
345 {
346 double radius_ratio;
347 radius_ratio = numerator.length() * area_sum / (108 * volume * volume);
348
349 return (double) VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
350 }
351
352 }
353
354 /*!
355 The aspect ratio of a tet
356
357 NB (P. Pebay 01/22/07):
358 Hmax / (2 sqrt(6) r) where Hmax and r respectively denote the greatest edge
359 length and the inradius of the tetrahedron
360 NB (J. Pouderoux 01/27/15)
361 This will return VERDICT_DBL_MAX when the volume of the tetrahedron is ill-
362 conditioned. Previously, this would only happen when the volume was small
363 and positive, but now ill-conditioned inverted tetrahedra are also included.
364 */
v_tet_aspect_ratio(int,double coordinates[][3])365 C_FUNC_DEF double v_tet_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
366 {
367 static const double normal_coeff = sqrt(6.) / 12.;
368
369 //Determine side vectors
370 VerdictVector ab, bc, ac, ad, bd, cd;
371
372 ab.set( coordinates[1][0] - coordinates[0][0],
373 coordinates[1][1] - coordinates[0][1],
374 coordinates[1][2] - coordinates[0][2] );
375
376 ac.set( coordinates[2][0] - coordinates[0][0],
377 coordinates[2][1] - coordinates[0][1],
378 coordinates[2][2] - coordinates[0][2] );
379
380 ad.set( coordinates[3][0] - coordinates[0][0],
381 coordinates[3][1] - coordinates[0][1],
382 coordinates[3][2] - coordinates[0][2] );
383
384 double detTet = ab % ( ac * ad );
385
386 if( fabs( detTet ) < VERDICT_DBL_MIN )
387 return (double)VERDICT_DBL_MAX;
388
389 bc.set( coordinates[2][0] - coordinates[1][0],
390 coordinates[2][1] - coordinates[1][1],
391 coordinates[2][2] - coordinates[1][2] );
392
393 bd.set( coordinates[3][0] - coordinates[1][0],
394 coordinates[3][1] - coordinates[1][1],
395 coordinates[3][2] - coordinates[1][2] );
396
397 cd.set( coordinates[3][0] - coordinates[2][0],
398 coordinates[3][1] - coordinates[2][1],
399 coordinates[3][2] - coordinates[2][2] );
400
401 double ab2 = ab.length_squared();
402 double bc2 = bc.length_squared();
403 double ac2 = ac.length_squared();
404 double ad2 = ad.length_squared();
405 double bd2 = bd.length_squared();
406 double cd2 = cd.length_squared();
407
408 double A = ab2 > bc2 ? ab2 : bc2;
409 double B = ac2 > ad2 ? ac2 : ad2;
410 double C = bd2 > cd2 ? bd2 : cd2;
411 double D = A > B ? A : B;
412 double hm = D > C ? sqrt( D ) : sqrt( C );
413
414 bd = ab * bc;
415 A = bd.length();
416 bd = ab * ad;
417 B = bd.length();
418 bd = ac * ad;
419 C = bd.length();
420 bd = bc * cd;
421 D = bd.length();
422
423 double aspect_ratio;
424 aspect_ratio = normal_coeff * hm * ( A + B + C + D ) / fabs( detTet );
425
426 if( aspect_ratio > 0 )
427 return (double) VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
428 return (double) VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
429 }
430
431 /*!
432 the aspect gamma of a tet
433
434 srms^3 / (8.48528137423857*V) where srms = sqrt(sum(Si^2)/6), where Si is the edge length
435 */
v_tet_aspect_gamma(int,double coordinates[][3])436 C_FUNC_DEF double v_tet_aspect_gamma( int /*num_nodes*/, double coordinates[][3] )
437 {
438
439 //Determine side vectors
440 VerdictVector side0, side1, side2, side3, side4, side5;
441
442 side0.set( coordinates[1][0] - coordinates[0][0],
443 coordinates[1][1] - coordinates[0][1],
444 coordinates[1][2] - coordinates[0][2] );
445
446 side1.set( coordinates[2][0] - coordinates[1][0],
447 coordinates[2][1] - coordinates[1][1],
448 coordinates[2][2] - coordinates[1][2] );
449
450 side2.set( coordinates[0][0] - coordinates[2][0],
451 coordinates[0][1] - coordinates[2][1],
452 coordinates[0][2] - coordinates[2][2] );
453
454 side3.set( coordinates[3][0] - coordinates[0][0],
455 coordinates[3][1] - coordinates[0][1],
456 coordinates[3][2] - coordinates[0][2] );
457
458 side4.set( coordinates[3][0] - coordinates[1][0],
459 coordinates[3][1] - coordinates[1][1],
460 coordinates[3][2] - coordinates[1][2] );
461
462 side5.set( coordinates[3][0] - coordinates[2][0],
463 coordinates[3][1] - coordinates[2][1],
464 coordinates[3][2] - coordinates[2][2] );
465
466
467 double volume = fabs( v_tet_volume(4, coordinates) );
468
469 if( volume < VERDICT_DBL_MIN )
470 return (double)VERDICT_DBL_MAX;
471 else
472 {
473 double srms = sqrt((side0.length_squared() + side1.length_squared() +
474 side2.length_squared() + side3.length_squared() +
475 side4.length_squared() + side5.length_squared()) / 6.0 );
476
477 double aspect_ratio_gamma = pow(srms, 3) / (8.48528137423857 * volume );
478 return (double)aspect_ratio_gamma;
479 }
480 }
481
482 /*!
483 The aspect frobenius of a tet
484
485 NB (P. Pebay 01/22/07):
486 Frobenius condition number when the reference element is regular
487 */
v_tet_aspect_frobenius(int,double coordinates[][3])488 C_FUNC_DEF double v_tet_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
489 {
490 static const double normal_exp = 1. / 3.;
491
492 VerdictVector ab, ac, ad;
493
494 ab.set( coordinates[1][0] - coordinates[0][0],
495 coordinates[1][1] - coordinates[0][1],
496 coordinates[1][2] - coordinates[0][2] );
497
498 ac.set( coordinates[2][0] - coordinates[0][0],
499 coordinates[2][1] - coordinates[0][1],
500 coordinates[2][2] - coordinates[0][2] );
501
502 ad.set( coordinates[3][0] - coordinates[0][0],
503 coordinates[3][1] - coordinates[0][1],
504 coordinates[3][2] - coordinates[0][2] );
505
506 double denominator = ab % ( ac * ad );
507 denominator *= denominator;
508 denominator *= 2.;
509 denominator = 3. * pow( denominator, normal_exp );
510
511 if( denominator < VERDICT_DBL_MIN )
512 return (double)VERDICT_DBL_MAX;
513
514 double u[3];
515 ab.get_xyz( u );
516 double v[3];
517 ac.get_xyz( v );
518 double w[3];
519 ad.get_xyz( w );
520
521 double numerator = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
522 numerator += v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
523 numerator += w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
524 numerator *= 1.5;
525 numerator -= v[0] * u[0] + v[1] * u[1] + v[2] * u[2];
526 numerator -= w[0] * u[0] + w[1] * u[1] + w[2] * u[2];
527 numerator -= w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
528
529 double aspect_frobenius = numerator / denominator;
530
531 if( aspect_frobenius > 0 )
532 return (double) VERDICT_MIN( aspect_frobenius, VERDICT_DBL_MAX );
533 return (double) VERDICT_MAX( aspect_frobenius, -VERDICT_DBL_MAX );
534 }
535
536 /*!
537 The minimum angle of a tet
538
539 NB (P. Pebay 01/22/07):
540 minimum nonoriented dihedral angle
541 */
v_tet_minimum_angle(int,double coordinates[][3])542 C_FUNC_DEF double v_tet_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
543 {
544 static const double normal_coeff = 180. * .3183098861837906715377675267450287;
545
546 //Determine side vectors
547 VerdictVector ab, bc, ad, cd;
548
549 ab.set( coordinates[1][0] - coordinates[0][0],
550 coordinates[1][1] - coordinates[0][1],
551 coordinates[1][2] - coordinates[0][2] );
552
553 ad.set( coordinates[3][0] - coordinates[0][0],
554 coordinates[3][1] - coordinates[0][1],
555 coordinates[3][2] - coordinates[0][2] );
556
557 bc.set( coordinates[2][0] - coordinates[1][0],
558 coordinates[2][1] - coordinates[1][1],
559 coordinates[2][2] - coordinates[1][2] );
560
561 cd.set( coordinates[3][0] - coordinates[2][0],
562 coordinates[3][1] - coordinates[2][1],
563 coordinates[3][2] - coordinates[2][2] );
564
565 VerdictVector abc = ab * bc;
566 double nabc = abc.length();
567 VerdictVector abd = ab * ad;
568 double nabd = abd.length();
569 VerdictVector acd = ad * cd;
570 double nacd = acd.length();
571 VerdictVector bcd = bc * cd;
572 double nbcd = bcd.length();
573
574 double alpha = acos( ( abc % abd ) / ( nabc * nabd ) );
575 double beta = acos( ( abc % acd ) / ( nabc * nacd ) );
576 double gamma = acos( ( abc % bcd ) / ( nabc * nbcd ) );
577 double delta = acos( ( abd % acd ) / ( nabd * nacd ) );
578 double epsilon = acos( ( abd % bcd ) / ( nabd * nbcd ) );
579 double zeta = acos( ( acd % bcd ) / ( nacd * nbcd ) );
580
581 alpha = alpha < beta ? alpha : beta;
582 alpha = alpha < gamma ? alpha : gamma;
583 alpha = alpha < delta ? alpha : delta;
584 alpha = alpha < epsilon ? alpha : epsilon;
585 alpha = alpha < zeta ? alpha : zeta;
586 alpha *= normal_coeff;
587
588 if( alpha < VERDICT_DBL_MIN )
589 return (double)VERDICT_DBL_MAX;
590
591 if( alpha > 0 )
592 return (double) VERDICT_MIN( alpha, VERDICT_DBL_MAX );
593 return (double) VERDICT_MAX( alpha, -VERDICT_DBL_MAX );
594 }
595
596 /*!
597 The collapse ratio of a tet
598
599 NB (J. Pouderoux 01/27/15)
600 This will return VERDICT_DBL_MAX when the volume of the tetrahedron is ill-
601 conditioned. Previously, this would only happen when the volume was small
602 and positive, but now ill-conditioned inverted tetrahedra are also included.
603 */
v_tet_collapse_ratio(int,double coordinates[][3])604 C_FUNC_DEF double v_tet_collapse_ratio( int /*num_nodes*/, double coordinates[][3] )
605 {
606 //Determine side vectors
607 VerdictVector e01, e02, e03, e12, e13, e23;
608
609 e01.set( coordinates[1][0] - coordinates[0][0],
610 coordinates[1][1] - coordinates[0][1],
611 coordinates[1][2] - coordinates[0][2] );
612
613 e02.set( coordinates[2][0] - coordinates[0][0],
614 coordinates[2][1] - coordinates[0][1],
615 coordinates[2][2] - coordinates[0][2] );
616
617 e03.set( coordinates[3][0] - coordinates[0][0],
618 coordinates[3][1] - coordinates[0][1],
619 coordinates[3][2] - coordinates[0][2] );
620
621 e12.set( coordinates[2][0] - coordinates[1][0],
622 coordinates[2][1] - coordinates[1][1],
623 coordinates[2][2] - coordinates[1][2] );
624
625 e13.set( coordinates[3][0] - coordinates[1][0],
626 coordinates[3][1] - coordinates[1][1],
627 coordinates[3][2] - coordinates[1][2] );
628
629 e23.set( coordinates[3][0] - coordinates[2][0],
630 coordinates[3][1] - coordinates[2][1],
631 coordinates[3][2] - coordinates[2][2] );
632
633 double l[6];
634 l[0] = e01.length();
635 l[1] = e02.length();
636 l[2] = e03.length();
637 l[3] = e12.length();
638 l[4] = e13.length();
639 l[5] = e23.length();
640
641 // Find longest edge for each bounding triangle of tetrahedron
642 double l012 = l[4] > l[0] ? l[4] : l[0]; l012 = l[1] > l012 ? l[1] : l012;
643 double l031 = l[0] > l[2] ? l[0] : l[2]; l031 = l[3] > l031 ? l[3] : l031;
644 double l023 = l[2] > l[1] ? l[2] : l[1]; l023 = l[5] > l023 ? l[5] : l023;
645 double l132 = l[4] > l[3] ? l[4] : l[3]; l132 = l[5] > l132 ? l[5] : l132;
646
647 // Compute collapse ratio for each vertex/triangle pair
648 VerdictVector N;
649 double h, magN;
650 double cr;
651 double crMin;
652
653 N = e01 * e02;
654 magN = N.length();
655 h = ( e03 % N) / magN; // height of vertex 3 above 0-1-2
656 crMin = h / l012; // ratio of height to longest edge of 0-1-2
657
658 N = e03 * e01 ;
659 magN = N.length();
660 h = ( e02 % N) / magN; // height of vertex 2 above 0-3-1
661 cr = h / l031; // ratio of height to longest edge of 0-3-1
662 if ( cr < crMin ) crMin = cr;
663
664 N = e02 * e03 ;
665 magN = N.length();
666 h = ( e01 % N) / magN; // height of vertex 1 above 0-2-3
667 cr = h / l023; // ratio of height to longest edge of 0-2-3
668 if ( cr < crMin ) crMin = cr;
669
670 N = e12 * e13 ;
671 magN = N.length();
672 h = ( e01 % N ) / magN; // height of vertex 0 above 1-3-2
673 cr = h / l132; // ratio of height to longest edge of 1-3-2
674 if ( cr < crMin ) crMin = cr;
675
676 if( fabs( crMin ) < VERDICT_DBL_MIN )
677 return (double)VERDICT_DBL_MAX;
678 if( crMin > 0 )
679 return (double) VERDICT_MIN( crMin, VERDICT_DBL_MAX );
680 return (double) VERDICT_MAX( crMin, -VERDICT_DBL_MAX );
681 }
682
v_tet_equivolume_skew(int num_nodes,double coordinates[][3])683 C_FUNC_DEF double v_tet_equivolume_skew( int num_nodes, double coordinates[][3] )
684 {
685
686 //- Find the vectors from the origin to each of the nodes on the tet.
687 VerdictVector vectA(coordinates[0][0],coordinates[0][1],coordinates[0][2]);
688 VerdictVector vectB(coordinates[1][0],coordinates[1][1],coordinates[1][2]);
689 VerdictVector vectC(coordinates[2][0],coordinates[2][1],coordinates[2][2]);
690 VerdictVector vectD(coordinates[3][0],coordinates[3][1],coordinates[3][2]);
691
692 VerdictVector vectAB = vectB - vectA;
693 VerdictVector vectAC = vectC - vectA;
694 VerdictVector vectAD = vectD - vectA;
695
696
697 double sq_lengthAB = vectAB.length_squared();
698 double sq_lengthAC = vectAC.length_squared();
699 double sq_lengthAD = vectAD.length_squared();
700
701 VerdictVector cpBC=vectAB * vectAC;
702 VerdictVector cpDB=vectAD * vectAB ;
703 VerdictVector cpCD=vectAC * vectAD;
704
705
706
707 VerdictVector num=sq_lengthAD*cpBC+sq_lengthAC*cpDB+sq_lengthAB*cpCD;
708 double den=2*vectAB%cpCD;
709
710 double circumradius=num.length()/den;
711
712
713 double volume=v_tet_volume(num_nodes,coordinates);
714 double optimal_length=circumradius/sqrt(double(3.0)/8.0);
715 double optimal_volume=(1.0/12.0)*sqrt(double(2.0))*pow(optimal_length,3);
716
717 return (optimal_volume-volume)/optimal_volume;
718 }
719
v_tet_squish_index(int,double coordinates[][3])720 C_FUNC_DEF double v_tet_squish_index( int /*num_nodes*/, double coordinates[][3] )
721 {
722 VerdictVector vectA(coordinates[0][0],coordinates[0][1],coordinates[0][2]);
723 VerdictVector vectB(coordinates[1][0],coordinates[1][1],coordinates[1][2]);
724 VerdictVector vectC(coordinates[2][0],coordinates[2][1],coordinates[2][2]);
725 VerdictVector vectD(coordinates[3][0],coordinates[3][1],coordinates[3][2]);
726
727
728 VerdictVector tetCenter=vectA+vectB+vectC+vectD;
729 tetCenter/=4.0;
730
731 /* top view
732
733 C
734 /|\
735 / 5 \
736 2 / D \ 1
737 / 3/ \4 \
738 /_/ \_\
739 A-----------B
740 0
741 */
742
743
744
745 VerdictVector side[6];
746
747 side[0].set( vectA,vectB);
748 side[1].set( vectB,vectC);
749 side[2].set( vectC,vectA);
750 side[3].set( vectA,vectD);
751 side[4].set( vectB,vectD);
752 side[5].set( vectC,vectD);
753
754
755 double maxSquishIndex=0;
756 double squishIndex=0;
757 VerdictVector faceCenter;
758 VerdictVector centerCenterVector;
759 VerdictVector faceAreaVector;
760
761 //face 1
762 faceCenter=(vectA+vectB+vectD)/3.0;
763 centerCenterVector=faceCenter-tetCenter;
764 faceAreaVector=0.5*(side[0]*side[4]);
765
766 squishIndex=1-(faceAreaVector%centerCenterVector)/(faceAreaVector.length()*centerCenterVector.length());
767 if(squishIndex>maxSquishIndex)
768 maxSquishIndex=squishIndex;
769
770 //face 2
771 faceCenter=(vectB+vectC+vectD)/3.0;
772 centerCenterVector=faceCenter-tetCenter;
773 faceAreaVector=0.5*(side[1]*side[5]);
774
775 squishIndex=1-(faceAreaVector%centerCenterVector)/(faceAreaVector.length()*centerCenterVector.length());
776 if(squishIndex>maxSquishIndex)
777 maxSquishIndex=squishIndex;
778
779 //face 3
780 faceCenter=(vectA+vectC+vectD)/3.0;
781 centerCenterVector=faceCenter-tetCenter;
782 faceAreaVector=0.5*(side[2]*side[3]);
783
784 squishIndex=1-(faceAreaVector%centerCenterVector)/(faceAreaVector.length()*centerCenterVector.length());
785 if(squishIndex>maxSquishIndex)
786 maxSquishIndex=squishIndex;
787
788 //face 4
789 faceCenter=(vectA+vectB+vectC)/3.0;
790 centerCenterVector=faceCenter-tetCenter;
791 faceAreaVector=0.5*(side[1]*side[0]);
792
793 squishIndex=1-(faceAreaVector%centerCenterVector)/(faceAreaVector.length()*centerCenterVector.length());
794 if(squishIndex>maxSquishIndex)
795 maxSquishIndex=squishIndex;
796
797
798
799 return maxSquishIndex;
800 }
801
802
803
804 /*!
805 the volume of a tet
806
807 1/6 * jacobian at a corner node
808 */
v_tet_volume(int,double coordinates[][3])809 C_FUNC_DEF double v_tet_volume( int /*num_nodes*/, double coordinates[][3] )
810 {
811
812 //Determine side vectors
813 VerdictVector side0, side2, side3;
814
815 side0.set( coordinates[1][0] - coordinates[0][0],
816 coordinates[1][1] - coordinates[0][1],
817 coordinates[1][2] - coordinates[0][2] );
818
819 side2.set( coordinates[0][0] - coordinates[2][0],
820 coordinates[0][1] - coordinates[2][1],
821 coordinates[0][2] - coordinates[2][2] );
822
823 side3.set( coordinates[3][0] - coordinates[0][0],
824 coordinates[3][1] - coordinates[0][1],
825 coordinates[3][2] - coordinates[0][2] );
826
827 return (double)((side3 % (side2 * side0)) / 6.0);
828
829 }
830
831 /*!
832 the condition of a tet
833
834 condition number of the jacobian matrix at any corner
835
836 NB (J. Pouderoux 01/27/15)
837 This will return VERDICT_DBL_MAX when the volume of the tetrahedron is ill-
838 conditioned. Previously, this would only happen when the volume was small
839 and positive, but now ill-conditioned inverted tetrahedra are also included.
840 */
v_tet_condition(int,double coordinates[][3])841 C_FUNC_DEF double v_tet_condition( int /*num_nodes*/, double coordinates[][3] )
842 {
843
844 double condition, term1, term2, det;
845 double rt3 = sqrt(3.0);
846 double rt6 = sqrt(6.0);
847
848 VerdictVector side0, side2, side3;
849
850 side0.set(coordinates[1][0] - coordinates[0][0],
851 coordinates[1][1] - coordinates[0][1],
852 coordinates[1][2] - coordinates[0][2]);
853
854 side2.set(coordinates[0][0] - coordinates[2][0],
855 coordinates[0][1] - coordinates[2][1],
856 coordinates[0][2] - coordinates[2][2]);
857
858 side3.set(coordinates[3][0] - coordinates[0][0],
859 coordinates[3][1] - coordinates[0][1],
860 coordinates[3][2] - coordinates[0][2]);
861
862 VerdictVector c_1, c_2, c_3;
863
864 c_1 = side0;
865 c_2 = (-2*side2-side0)/rt3;
866 c_3 = (3*side3+side2-side0)/rt6;
867
868 term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
869 term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) +
870 ( c_2 * c_3 ) % ( c_2 * c_3 ) +
871 ( c_1 * c_3 ) % ( c_1 * c_3 );
872 det = c_1 % ( c_2 * c_3 );
873
874 if ( fabs( det ) <= VERDICT_DBL_MIN )
875 return VERDICT_DBL_MAX;
876 else
877 condition = sqrt( term1 * term2 ) /(3.0* det);
878
879 return (double)condition;
880 }
881
882
883 /*!
884 the jacobian of a tet
885
886 TODO
887 */
v_tet_jacobian(int,double coordinates[][3])888 C_FUNC_DEF double v_tet_jacobian( int /*num_nodes*/, double coordinates[][3] )
889 {
890 VerdictVector side0, side2, side3;
891
892 side0.set( coordinates[1][0] - coordinates[0][0],
893 coordinates[1][1] - coordinates[0][1],
894 coordinates[1][2] - coordinates[0][2] );
895
896 side2.set( coordinates[0][0] - coordinates[2][0],
897 coordinates[0][1] - coordinates[2][1],
898 coordinates[0][2] - coordinates[2][2] );
899
900 side3.set( coordinates[3][0] - coordinates[0][0],
901 coordinates[3][1] - coordinates[0][1],
902 coordinates[3][2] - coordinates[0][2] );
903
904
905 return (double)(side3 % (side2 * side0));
906
907 }
908
909
910 /*!
911 the shape of a tet
912
913 3/ condition number of weighted jacobian matrix
914 */
v_tet_shape(int,double coordinates[][3])915 C_FUNC_DEF double v_tet_shape( int /*num_nodes*/, double coordinates[][3] )
916 {
917
918 static const double two_thirds = 2.0/3.0;
919 static const double root_of_2 = sqrt(2.0);
920
921 VerdictVector edge0, edge2, edge3;
922
923 edge0.set(coordinates[1][0] - coordinates[0][0],
924 coordinates[1][1] - coordinates[0][1],
925 coordinates[1][2] - coordinates[0][2]);
926
927 edge2.set(coordinates[0][0] - coordinates[2][0],
928 coordinates[0][1] - coordinates[2][1],
929 coordinates[0][2] - coordinates[2][2]);
930
931 edge3.set(coordinates[3][0] - coordinates[0][0],
932 coordinates[3][1] - coordinates[0][1],
933 coordinates[3][2] - coordinates[0][2]);
934
935 double jacobian = edge3 % (edge2 * edge0);
936 if(jacobian < VERDICT_DBL_MIN){
937 return (double)0.0;
938 }
939 double num = 3 * pow( root_of_2 * jacobian, two_thirds );
940 double den = 1.5*(edge0%edge0 + edge2%edge2 + edge3%edge3)-
941 (edge0%-edge2 + -edge2%edge3 + edge3%edge0);
942
943 if ( den < VERDICT_DBL_MIN )
944 return (double)0.0;
945
946 return (double)VERDICT_MAX( num/den, 0 );
947 }
948
949
950
951 /*!
952 the relative size of a tet
953
954 Min(J,1/J), where J is the determinant of the weighted Jacobian matrix
955 */
v_tet_relative_size_squared(int,double coordinates[][3])956 C_FUNC_DEF double v_tet_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
957 {
958 double size;
959 VerdictVector w1, w2, w3;
960 v_tet_get_weight(w1,w2,w3);
961 double avg_volume = (w1 % (w2 *w3))/6.0;
962
963 double volume = v_tet_volume(4, coordinates);
964
965 if( avg_volume < VERDICT_DBL_MIN )
966 return 0.0;
967 else
968 {
969 size = volume/avg_volume;
970 if( size <= VERDICT_DBL_MIN )
971 return 0.0;
972 if ( size > 1 )
973 size = (double)(1)/size;
974 }
975 return (double)(size*size);
976 }
977
978
979 /*!
980 the shape and size of a tet
981
982 Product of the shape and relative size
983 */
v_tet_shape_and_size(int num_nodes,double coordinates[][3])984 C_FUNC_DEF double v_tet_shape_and_size( int num_nodes, double coordinates[][3] )
985 {
986
987 double shape, size;
988 shape = v_tet_shape( num_nodes, coordinates );
989 size = v_tet_relative_size_squared (num_nodes, coordinates );
990
991 return (double)(shape * size);
992
993 }
994
995
996
997 /*!
998 the distortion of a tet
999 */
v_tet_distortion(int num_nodes,double coordinates[][3])1000 C_FUNC_DEF double v_tet_distortion( int num_nodes, double coordinates[][3] )
1001 {
1002
1003 double distortion = VERDICT_DBL_MAX;
1004 int number_of_gauss_points=0;
1005 if (num_nodes ==4)
1006 // for linear tet, the distortion is always 1 because
1007 // straight edge tets are the target shape for tet
1008 return 1.0;
1009
1010 else if (num_nodes ==10)
1011 //use four integration points for quadratic tet
1012 number_of_gauss_points = 4;
1013
1014 int number_dims = 3;
1015 int total_number_of_gauss_points = number_of_gauss_points;
1016 // use is_tri=1 to indicate this is for tet in 3D
1017 int is_tri =1;
1018
1019 double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
1020 double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
1021 double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
1022 double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
1023 double weight[maxTotalNumberGaussPoints];
1024
1025 // create an object of GaussIntegration for tet
1026 GaussIntegration::initialize(number_of_gauss_points,num_nodes, number_dims, is_tri);
1027 GaussIntegration::calculate_shape_function_3d_tet();
1028 GaussIntegration::get_shape_func(shape_function[0], dndy1[0], dndy2[0], dndy3[0],weight);
1029
1030 // vector xxi is the derivative vector of coordinates w.r.t local xi coordinate in the
1031 // computation space
1032 // vector xet is the derivative vector of coordinates w.r.t local et coordinate in the
1033 // computation space
1034 // vector xze is the derivative vector of coordinates w.r.t local ze coordinate in the
1035 // computation space
1036 VerdictVector xxi, xet, xze, xin;
1037
1038 double jacobian, minimum_jacobian;
1039 double element_volume =0.0;
1040 minimum_jacobian = VERDICT_DBL_MAX;
1041
1042 // calculate element volume
1043 int ife, ja;
1044 for (ife=0;ife<total_number_of_gauss_points; ife++)
1045 {
1046 xxi.set(0.0,0.0,0.0);
1047 xet.set(0.0,0.0,0.0);
1048 xze.set(0.0,0.0,0.0);
1049
1050 for (ja=0;ja<num_nodes;ja++)
1051 {
1052 xin.set(coordinates[ja][0], coordinates[ja][1], coordinates[ja][2]);
1053 xxi += dndy1[ife][ja]*xin;
1054 xet += dndy2[ife][ja]*xin;
1055 xze += dndy3[ife][ja]*xin;
1056 }
1057
1058 // determinant
1059 jacobian = xxi % ( xet * xze );
1060 if (minimum_jacobian > jacobian)
1061 minimum_jacobian = jacobian;
1062
1063 element_volume += weight[ife]*jacobian;
1064 }//element_volume is 6 times the actual volume
1065
1066 // loop through all nodes
1067 double dndy1_at_node[maxNumberNodes][maxNumberNodes];
1068 double dndy2_at_node[maxNumberNodes][maxNumberNodes];
1069 double dndy3_at_node[maxNumberNodes][maxNumberNodes];
1070
1071
1072 GaussIntegration::calculate_derivative_at_nodes_3d_tet( dndy1_at_node,
1073 dndy2_at_node,
1074 dndy3_at_node);
1075 int node_id;
1076 for (node_id=0;node_id<num_nodes; node_id++)
1077 {
1078 xxi.set(0.0,0.0,0.0);
1079 xet.set(0.0,0.0,0.0);
1080 xze.set(0.0,0.0,0.0);
1081
1082 for (ja=0;ja<num_nodes;ja++)
1083 {
1084 xin.set(coordinates[ja][0], coordinates[ja][1], coordinates[ja][2]);
1085 xxi += dndy1_at_node[node_id][ja]*xin;
1086 xet += dndy2_at_node[node_id][ja]*xin;
1087 xze += dndy3_at_node[node_id][ja]*xin;
1088 }
1089
1090 jacobian = xxi % ( xet * xze );
1091 if (minimum_jacobian > jacobian)
1092 minimum_jacobian = jacobian;
1093
1094 }
1095 distortion = minimum_jacobian/element_volume;
1096
1097 return (double)distortion;
1098 }
1099
1100
1101 /*!
1102 the quality functions of a tet
1103 */
v_tet_quality(int num_nodes,double coordinates[][3],unsigned int metrics_request_flag,TetMetricVals * metric_vals)1104 C_FUNC_DEF void v_tet_quality( int num_nodes, double coordinates[][3],
1105 unsigned int metrics_request_flag, TetMetricVals *metric_vals )
1106 {
1107
1108 memset( metric_vals, 0, sizeof(TetMetricVals) );
1109
1110 /*
1111
1112 node numbers and edge numbers below
1113
1114
1115
1116 3
1117 + edge 0 is node 0 to 1
1118 +|+ edge 1 is node 1 to 2
1119 3/ | \5 edge 2 is node 0 to 2
1120 / 4| \ edge 3 is node 0 to 3
1121 0 - -|- + 2 edge 4 is node 1 to 3
1122 \ | + edge 5 is node 2 to 3
1123 0\ | /1
1124 +|/ edge 2 is behind edge 4
1125 1
1126
1127
1128 */
1129
1130 // lets start with making the vectors
1131 VerdictVector edges[6];
1132 edges[0].set( coordinates[1][0] - coordinates[0][0],
1133 coordinates[1][1] - coordinates[0][1],
1134 coordinates[1][2] - coordinates[0][2] );
1135
1136 edges[1].set( coordinates[2][0] - coordinates[1][0],
1137 coordinates[2][1] - coordinates[1][1],
1138 coordinates[2][2] - coordinates[1][2] );
1139
1140 edges[2].set( coordinates[0][0] - coordinates[2][0],
1141 coordinates[0][1] - coordinates[2][1],
1142 coordinates[0][2] - coordinates[2][2] );
1143
1144 edges[3].set( coordinates[3][0] - coordinates[0][0],
1145 coordinates[3][1] - coordinates[0][1],
1146 coordinates[3][2] - coordinates[0][2] );
1147
1148 edges[4].set( coordinates[3][0] - coordinates[1][0],
1149 coordinates[3][1] - coordinates[1][1],
1150 coordinates[3][2] - coordinates[1][2] );
1151
1152 edges[5].set( coordinates[3][0] - coordinates[2][0],
1153 coordinates[3][1] - coordinates[2][1],
1154 coordinates[3][2] - coordinates[2][2] );
1155
1156 // common numbers
1157 static const double root_of_2 = sqrt(2.0);
1158
1159 // calculate the jacobian
1160 static const int do_jacobian = V_TET_JACOBIAN | V_TET_VOLUME |
1161 V_TET_ASPECT_BETA | V_TET_ASPECT_GAMMA | V_TET_SHAPE |
1162 V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE |
1163 V_TET_SCALED_JACOBIAN | V_TET_CONDITION;
1164 if(metrics_request_flag & do_jacobian )
1165 {
1166 metric_vals->jacobian = (double)(edges[3] % (edges[2] * edges[0]));
1167 }
1168
1169 // calculate the volume
1170 if(metrics_request_flag & V_TET_VOLUME)
1171 {
1172 metric_vals->volume = (double)(metric_vals->jacobian / 6.0);
1173 }
1174
1175 // calculate aspect ratio
1176 if(metrics_request_flag & V_TET_ASPECT_BETA)
1177 {
1178 double surface_area = ((edges[2] * edges[0]).length() +
1179 (edges[3] * edges[0]).length() +
1180 (edges[4] * edges[1]).length() +
1181 (edges[3] * edges[2]).length() ) * 0.5;
1182
1183 VerdictVector numerator = edges[3].length_squared() * ( edges[2] * edges[0] ) +
1184 edges[2].length_squared() * ( edges[3] * edges[0] ) +
1185 edges[0].length_squared() * ( edges[3] * edges[2] );
1186
1187 double volume = metric_vals->jacobian / 6.0;
1188
1189 if(volume < VERDICT_DBL_MIN )
1190 metric_vals->aspect_beta = (double)(VERDICT_DBL_MAX);
1191 else
1192 metric_vals->aspect_beta =
1193 (double)( numerator.length() * surface_area/ (108*volume*volume) );
1194 }
1195
1196 // calculate the aspect gamma
1197 if(metrics_request_flag & V_TET_ASPECT_GAMMA)
1198 {
1199 double volume = fabs( metric_vals->jacobian / 6.0 );
1200 if( fabs( volume ) < VERDICT_DBL_MIN )
1201 metric_vals->aspect_gamma = VERDICT_DBL_MAX;
1202 else
1203 {
1204 double srms = sqrt((
1205 edges[0].length_squared() + edges[1].length_squared() +
1206 edges[2].length_squared() + edges[3].length_squared() +
1207 edges[4].length_squared() + edges[5].length_squared()
1208 ) / 6.0 );
1209
1210 // cube the srms
1211 srms *= (srms * srms);
1212 metric_vals->aspect_gamma = (double)( srms / (8.48528137423857 * volume ));
1213 }
1214 }
1215
1216 // calculate the shape of the tet
1217 if(metrics_request_flag & ( V_TET_SHAPE | V_TET_SHAPE_AND_SIZE ) )
1218 {
1219 //if the jacobian is non-positive, the shape is 0
1220 if(metric_vals->jacobian < VERDICT_DBL_MIN){
1221 metric_vals->shape = (double)0.0;
1222 }
1223 else{
1224 static const double two_thirds = 2.0/3.0;
1225 double num = 3.0 * pow(root_of_2 * metric_vals->jacobian, two_thirds);
1226 double den = 1.5 *
1227 (edges[0] % edges[0] + edges[2] % edges[2] + edges[3] % edges[3]) -
1228 (edges[0] % -edges[2] + -edges[2] % edges[3] + edges[3] % edges[0]);
1229
1230 if( den < VERDICT_DBL_MIN )
1231 metric_vals->shape = (double)0.0;
1232 else
1233 metric_vals->shape = (double)VERDICT_MAX( num/den, 0 );
1234 }
1235
1236 }
1237
1238 // calculate the relative size of the tet
1239 if(metrics_request_flag & (V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE ))
1240 {
1241 VerdictVector w1, w2, w3;
1242 v_tet_get_weight(w1,w2,w3);
1243 double avg_vol = (w1 % (w2 *w3))/6;
1244
1245 if( avg_vol < VERDICT_DBL_MIN )
1246 metric_vals->relative_size_squared = 0.0;
1247 else
1248 {
1249 double tmp = metric_vals->jacobian / (6*avg_vol);
1250 if( tmp < VERDICT_DBL_MIN )
1251 metric_vals->relative_size_squared = 0.0;
1252 else
1253 {
1254 tmp *= tmp;
1255 metric_vals->relative_size_squared = (double)VERDICT_MIN(tmp, 1/tmp);
1256 }
1257 }
1258 }
1259
1260 // calculate the shape and size
1261 if(metrics_request_flag & V_TET_SHAPE_AND_SIZE)
1262 {
1263 metric_vals->shape_and_size = (double)(metric_vals->shape * metric_vals->relative_size_squared);
1264 }
1265
1266 // calculate the scaled jacobian
1267 if(metrics_request_flag & V_TET_SCALED_JACOBIAN)
1268 {
1269 //find out which node the normalized jacobian can be calculated at
1270 //and it will be the smaller than at other nodes
1271 double length_squared[4] = {
1272 edges[0].length_squared() * edges[2].length_squared() * edges[3].length_squared(),
1273 edges[0].length_squared() * edges[1].length_squared() * edges[4].length_squared(),
1274 edges[1].length_squared() * edges[2].length_squared() * edges[5].length_squared(),
1275 edges[3].length_squared() * edges[4].length_squared() * edges[5].length_squared()
1276 };
1277
1278 int which_node = 0;
1279 if(length_squared[1] > length_squared[which_node])
1280 which_node = 1;
1281 if(length_squared[2] > length_squared[which_node])
1282 which_node = 2;
1283 if(length_squared[3] > length_squared[which_node])
1284 which_node = 3;
1285
1286 // find the scaled jacobian at this node
1287 double length_product = sqrt( length_squared[which_node] );
1288 if(length_product < fabs(metric_vals->jacobian))
1289 length_product = fabs(metric_vals->jacobian);
1290
1291 if( length_product < VERDICT_DBL_MIN )
1292 metric_vals->scaled_jacobian = (double) VERDICT_DBL_MAX;
1293 else
1294 metric_vals->scaled_jacobian =
1295 (double)(root_of_2 * metric_vals->jacobian / length_product);
1296 }
1297
1298 // calculate the condition number
1299 if(metrics_request_flag & V_TET_CONDITION)
1300 {
1301 static const double root_of_3 = sqrt(3.0);
1302 static const double root_of_6 = sqrt(6.0);
1303
1304 VerdictVector c_1, c_2, c_3;
1305 c_1 = edges[0];
1306 c_2 = (-2*edges[2] - edges[0])/root_of_3;
1307 c_3 = (3*edges[3] + edges[2] - edges[0])/root_of_6;
1308
1309 double term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
1310 double term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) +
1311 ( c_2 * c_3 ) % ( c_2 * c_3 ) +
1312 ( c_3 * c_1 ) % ( c_3 * c_1 );
1313
1314 double det = c_1 % ( c_2 * c_3 );
1315
1316 if(det <= VERDICT_DBL_MIN)
1317 metric_vals->condition = (double)VERDICT_DBL_MAX;
1318 else
1319 metric_vals->condition = (double)(sqrt(term1 * term2) / (3.0*det));
1320 }
1321
1322 // calculate the distortion
1323 if(metrics_request_flag & V_TET_DISTORTION)
1324 {
1325 metric_vals->distortion = v_tet_distortion(num_nodes, coordinates);
1326 }
1327
1328
1329 // calculate the equivolume skew
1330 if(metrics_request_flag & V_TET_EQUIVOLUME_SKEW)
1331 {
1332 metric_vals->equivolume_skew = v_tet_equivolume_skew(num_nodes, coordinates);
1333 }
1334
1335
1336 // calculate the squish index
1337 if(metrics_request_flag & V_TET_SQUISH_INDEX)
1338 {
1339 metric_vals->squish_index = v_tet_squish_index(num_nodes, coordinates);
1340 }
1341
1342
1343 //check for overflow
1344 if(metrics_request_flag & V_TET_ASPECT_BETA )
1345 {
1346 if( metric_vals->aspect_beta > 0 )
1347 metric_vals->aspect_beta = (double) VERDICT_MIN( metric_vals->aspect_beta, VERDICT_DBL_MAX );
1348 metric_vals->aspect_beta = (double) VERDICT_MAX( metric_vals->aspect_beta, -VERDICT_DBL_MAX );
1349 }
1350
1351 if(metrics_request_flag & V_TET_ASPECT_GAMMA)
1352 {
1353 if( metric_vals->aspect_gamma > 0 )
1354 metric_vals->aspect_gamma = (double) VERDICT_MIN( metric_vals->aspect_gamma, VERDICT_DBL_MAX );
1355 metric_vals->aspect_gamma = (double) VERDICT_MAX( metric_vals->aspect_gamma, -VERDICT_DBL_MAX );
1356 }
1357
1358 if(metrics_request_flag & V_TET_VOLUME)
1359 {
1360 if( metric_vals->volume > 0 )
1361 metric_vals->volume = (double) VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
1362 metric_vals->volume = (double) VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
1363 }
1364
1365 if(metrics_request_flag & V_TET_CONDITION)
1366 {
1367 if( metric_vals->condition > 0 )
1368 metric_vals->condition = (double) VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
1369 metric_vals->condition = (double) VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
1370 }
1371
1372 if(metrics_request_flag & V_TET_JACOBIAN)
1373 {
1374 if( metric_vals->jacobian > 0 )
1375 metric_vals->jacobian = (double) VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
1376 metric_vals->jacobian = (double) VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
1377 }
1378
1379 if(metrics_request_flag & V_TET_SCALED_JACOBIAN)
1380 {
1381 if( metric_vals->scaled_jacobian > 0 )
1382 metric_vals->scaled_jacobian = (double) VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
1383 metric_vals->scaled_jacobian = (double) VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
1384 }
1385
1386 if(metrics_request_flag & V_TET_SHAPE)
1387 {
1388 if( metric_vals->shape > 0 )
1389 metric_vals->shape = (double) VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
1390 metric_vals->shape = (double) VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
1391 }
1392
1393 if(metrics_request_flag & V_TET_RELATIVE_SIZE_SQUARED)
1394 {
1395 if( metric_vals->relative_size_squared > 0 )
1396 metric_vals->relative_size_squared = (double) VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
1397 metric_vals->relative_size_squared = (double) VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
1398 }
1399
1400 if(metrics_request_flag & V_TET_SHAPE_AND_SIZE)
1401 {
1402 if( metric_vals->shape_and_size > 0 )
1403 metric_vals->shape_and_size = (double) VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
1404 metric_vals->shape_and_size = (double) VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
1405 }
1406
1407 if(metrics_request_flag & V_TET_DISTORTION)
1408 {
1409 if( metric_vals->distortion > 0 )
1410 metric_vals->distortion = (double) VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
1411 metric_vals->distortion = (double) VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
1412 }
1413
1414 if(metrics_request_flag & V_TET_EQUIVOLUME_SKEW)
1415 {
1416 if( metric_vals->equivolume_skew > 0 )
1417 metric_vals->equivolume_skew = (double) VERDICT_MIN( metric_vals->equivolume_skew, VERDICT_DBL_MAX );
1418 metric_vals->equivolume_skew = (double) VERDICT_MAX( metric_vals->equivolume_skew, -VERDICT_DBL_MAX );
1419 }
1420
1421 if(metrics_request_flag & V_TET_SQUISH_INDEX)
1422 {
1423 if( metric_vals->squish_index > 0 )
1424 metric_vals->squish_index = (double) VERDICT_MIN( metric_vals->squish_index, VERDICT_DBL_MAX );
1425 metric_vals->squish_index = (double) VERDICT_MAX( metric_vals->squish_index, -VERDICT_DBL_MAX );
1426 }
1427 }
1428