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