Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
V_TetMetric.cpp
Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Module:    $RCSfile: V_TetMetric.cpp,v $
00004 
00005   Copyright (c) 2006 Sandia Corporation.
00006   All rights reserved.
00007   See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
00008 
00009      This software is distributed WITHOUT ANY WARRANTY; without even
00010      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00011      PURPOSE.  See the above copyright notice for more information.
00012 
00013 =========================================================================*/
00014 
00015 /*
00016  *
00017  * TetMetric.cpp contains quality calculations for Tets
00018  *
00019  * This file is part of VERDICT
00020  *
00021  */
00022 
00023 #define VERDICT_EXPORTS
00024 
00025 #include "moab/verdict.h"
00026 #include "verdict_defines.hpp"
00027 #include "VerdictVector.hpp"
00028 #include "V_GaussIntegration.hpp"
00029 #include <memory.h>
00030 
00031 //! the average volume of a tet
00032 double verdict_tet_size = 0;
00033 
00034 /*!
00035   set the average volume of a tet
00036 */
00037 C_FUNC_DEF void v_set_tet_size( double size )
00038 {
00039     verdict_tet_size = size;
00040 }
00041 
00042 /*!
00043   get the weights based on the average size
00044   of a tet
00045 */
00046 int get_weight( VerdictVector& w1, VerdictVector& w2, VerdictVector& w3 )
00047 {
00048     static const double rt3       = sqrt( 3.0 );
00049     static const double root_of_2 = sqrt( 2.0 );
00050 
00051     w1.set( 1, 0, 0 );
00052     w2.set( 0.5, 0.5 * rt3, 0 );
00053     w3.set( 0.5, rt3 / 6.0, root_of_2 / rt3 );
00054 
00055     double scale = pow( 6. * verdict_tet_size / determinant( w1, w2, w3 ), 0.3333333333333 );
00056 
00057     w1 *= scale;
00058     w2 *= scale;
00059     w3 *= scale;
00060 
00061     return 1;
00062 }
00063 
00064 /*!
00065    the edge ratio of a tet
00066 
00067    NB (P. Pebay 01/22/07):
00068      Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
00069      minimum edge lengths
00070 */
00071 C_FUNC_DEF double v_tet_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
00072 {
00073     VerdictVector a, b, c, d, e, f;
00074 
00075     a.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00076            coordinates[1][2] - coordinates[0][2] );
00077 
00078     b.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00079            coordinates[2][2] - coordinates[1][2] );
00080 
00081     c.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00082            coordinates[0][2] - coordinates[2][2] );
00083 
00084     d.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00085            coordinates[3][2] - coordinates[0][2] );
00086 
00087     e.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00088            coordinates[3][2] - coordinates[1][2] );
00089 
00090     f.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00091            coordinates[3][2] - coordinates[2][2] );
00092 
00093     double a2 = a.length_squared();
00094     double b2 = b.length_squared();
00095     double c2 = c.length_squared();
00096     double d2 = d.length_squared();
00097     double e2 = e.length_squared();
00098     double f2 = f.length_squared();
00099 
00100     double m2, M2, mab, mcd, mef, Mab, Mcd, Mef;
00101 
00102     if( a2 < b2 )
00103     {
00104         mab = a2;
00105         Mab = b2;
00106     }
00107     else  // b2 <= a2
00108     {
00109         mab = b2;
00110         Mab = a2;
00111     }
00112     if( c2 < d2 )
00113     {
00114         mcd = c2;
00115         Mcd = d2;
00116     }
00117     else  // d2 <= c2
00118     {
00119         mcd = d2;
00120         Mcd = c2;
00121     }
00122     if( e2 < f2 )
00123     {
00124         mef = e2;
00125         Mef = f2;
00126     }
00127     else  // f2 <= e2
00128     {
00129         mef = f2;
00130         Mef = e2;
00131     }
00132 
00133     m2 = mab < mcd ? mab : mcd;
00134     m2 = m2 < mef ? m2 : mef;
00135 
00136     if( m2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00137 
00138     M2 = Mab > Mcd ? Mab : Mcd;
00139     M2 = M2 > Mef ? M2 : Mef;
00140 
00141     double edge_ratio = sqrt( M2 / m2 );
00142 
00143     if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
00144     return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
00145 }
00146 
00147 /*!
00148   the scaled jacobian of a tet
00149 
00150   minimum of the jacobian divided by the lengths of 3 edge vectors
00151 
00152 */
00153 C_FUNC_DEF double v_tet_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
00154 {
00155 
00156     VerdictVector side0, side1, side2, side3, side4, side5;
00157 
00158     side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00159                coordinates[1][2] - coordinates[0][2] );
00160 
00161     side1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00162                coordinates[2][2] - coordinates[1][2] );
00163 
00164     side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00165                coordinates[0][2] - coordinates[2][2] );
00166 
00167     side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00168                coordinates[3][2] - coordinates[0][2] );
00169 
00170     side4.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00171                coordinates[3][2] - coordinates[1][2] );
00172 
00173     side5.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00174                coordinates[3][2] - coordinates[2][2] );
00175 
00176     double jacobi;
00177 
00178     jacobi = side3 % ( side2 * side0 );
00179 
00180     // products of lengths squared of each edge attached to a node.
00181     double length_squared[4] = { side0.length_squared() * side2.length_squared() * side3.length_squared(),
00182                                  side0.length_squared() * side1.length_squared() * side4.length_squared(),
00183                                  side1.length_squared() * side2.length_squared() * side5.length_squared(),
00184                                  side3.length_squared() * side4.length_squared() * side5.length_squared() };
00185     int which_node           = 0;
00186     if( length_squared[1] > length_squared[which_node] ) which_node = 1;
00187     if( length_squared[2] > length_squared[which_node] ) which_node = 2;
00188     if( length_squared[3] > length_squared[which_node] ) which_node = 3;
00189 
00190     double length_product = sqrt( length_squared[which_node] );
00191     if( length_product < fabs( jacobi ) ) length_product = fabs( jacobi );
00192 
00193     if( length_product < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00194 
00195     static const double root_of_2 = sqrt( 2.0 );
00196 
00197     return (double)( root_of_2 * jacobi / length_product );
00198 }
00199 
00200 /*!
00201   The radius ratio of a tet
00202 
00203   NB (P. Pebay 04/16/07):
00204     CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed
00205     sphere radius.
00206     Note that this metric is similar to the aspect beta of a tet, except that
00207     it does not return VERDICT_DBL_MAX if the element has negative orientation.
00208 */
00209 C_FUNC_DEF double v_tet_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
00210 {
00211 
00212     // Determine side vectors
00213     VerdictVector side[6];
00214 
00215     side[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00216                  coordinates[1][2] - coordinates[0][2] );
00217 
00218     side[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00219                  coordinates[2][2] - coordinates[1][2] );
00220 
00221     side[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00222                  coordinates[0][2] - coordinates[2][2] );
00223 
00224     side[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00225                  coordinates[3][2] - coordinates[0][2] );
00226 
00227     side[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00228                  coordinates[3][2] - coordinates[1][2] );
00229 
00230     side[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00231                  coordinates[3][2] - coordinates[2][2] );
00232 
00233     VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
00234                               side[2].length_squared() * ( side[3] * side[0] ) +
00235                               side[0].length_squared() * ( side[3] * side[2] );
00236 
00237     double area_sum;
00238     area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
00239                  ( side[3] * side[2] ).length() ) *
00240                0.5;
00241 
00242     double volume = v_tet_volume( 4, coordinates );
00243 
00244     if( fabs( volume ) < VERDICT_DBL_MIN )
00245         return (double)VERDICT_DBL_MAX;
00246     else
00247     {
00248         double radius_ratio;
00249         radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );
00250 
00251         return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
00252     }
00253 }
00254 
00255 /*!
00256   The radius ratio of a positively-oriented tet, a.k.a. "aspect beta"
00257 
00258   NB (P. Pebay 04/16/07):
00259     CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed
00260     sphere radius if the element has positive orientation.
00261     Note that this metric is similar to the radius ratio of a tet, except that
00262     it returns VERDICT_DBL_MAX if the element has negative orientation.
00263 
00264 */
00265 C_FUNC_DEF double v_tet_aspect_beta( int /*num_nodes*/, double coordinates[][3] )
00266 {
00267 
00268     // Determine side vectors
00269     VerdictVector side[6];
00270 
00271     side[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00272                  coordinates[1][2] - coordinates[0][2] );
00273 
00274     side[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00275                  coordinates[2][2] - coordinates[1][2] );
00276 
00277     side[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00278                  coordinates[0][2] - coordinates[2][2] );
00279 
00280     side[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00281                  coordinates[3][2] - coordinates[0][2] );
00282 
00283     side[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00284                  coordinates[3][2] - coordinates[1][2] );
00285 
00286     side[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00287                  coordinates[3][2] - coordinates[2][2] );
00288 
00289     VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
00290                               side[2].length_squared() * ( side[3] * side[0] ) +
00291                               side[0].length_squared() * ( side[3] * side[2] );
00292 
00293     double area_sum;
00294     area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
00295                  ( side[3] * side[2] ).length() ) *
00296                0.5;
00297 
00298     double volume = v_tet_volume( 4, coordinates );
00299 
00300     if( volume < VERDICT_DBL_MIN )
00301         return (double)VERDICT_DBL_MAX;
00302     else
00303     {
00304         double radius_ratio;
00305         radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );
00306 
00307         return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
00308     }
00309 }
00310 
00311 /*!
00312   The aspect ratio of a tet
00313 
00314   NB (P. Pebay 01/22/07):
00315     Hmax / (2 sqrt(6) r) where Hmax and r respectively denote the greatest edge
00316     length and the inradius of the tetrahedron
00317 */
00318 C_FUNC_DEF double v_tet_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
00319 {
00320     static const double normal_coeff = sqrt( 6. ) / 12.;
00321 
00322     // Determine side vectors
00323     VerdictVector ab, bc, ac, ad, bd, cd;
00324 
00325     ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00326             coordinates[1][2] - coordinates[0][2] );
00327 
00328     ac.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00329             coordinates[2][2] - coordinates[0][2] );
00330 
00331     ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00332             coordinates[3][2] - coordinates[0][2] );
00333 
00334     double detTet = ab % ( ac * ad );
00335 
00336     if( detTet < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00337 
00338     bc.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00339             coordinates[2][2] - coordinates[1][2] );
00340 
00341     bd.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00342             coordinates[3][2] - coordinates[1][2] );
00343 
00344     cd.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00345             coordinates[3][2] - coordinates[2][2] );
00346 
00347     double ab2 = ab.length_squared();
00348     double bc2 = bc.length_squared();
00349     double ac2 = ac.length_squared();
00350     double ad2 = ad.length_squared();
00351     double bd2 = bd.length_squared();
00352     double cd2 = cd.length_squared();
00353 
00354     double A  = ab2 > bc2 ? ab2 : bc2;
00355     double B  = ac2 > ad2 ? ac2 : ad2;
00356     double C  = bd2 > cd2 ? bd2 : cd2;
00357     double D  = A > B ? A : B;
00358     double hm = D > C ? sqrt( D ) : sqrt( C );
00359 
00360     bd = ab * bc;
00361     A  = bd.length();
00362     bd = ab * ad;
00363     B  = bd.length();
00364     bd = ac * ad;
00365     C  = bd.length();
00366     bd = bc * cd;
00367     D  = bd.length();
00368 
00369     double aspect_ratio;
00370     aspect_ratio = normal_coeff * hm * ( A + B + C + D ) / fabs( detTet );
00371 
00372     if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
00373     return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
00374 }
00375 
00376 /*!
00377   the aspect gamma of a tet
00378 
00379   srms^3 / (8.48528137423857*V) where srms = sqrt(sum(Si^2)/6), where Si is the edge length
00380 */
00381 C_FUNC_DEF double v_tet_aspect_gamma( int /*num_nodes*/, double coordinates[][3] )
00382 {
00383 
00384     // Determine side vectors
00385     VerdictVector side0, side1, side2, side3, side4, side5;
00386 
00387     side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00388                coordinates[1][2] - coordinates[0][2] );
00389 
00390     side1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00391                coordinates[2][2] - coordinates[1][2] );
00392 
00393     side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00394                coordinates[0][2] - coordinates[2][2] );
00395 
00396     side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00397                coordinates[3][2] - coordinates[0][2] );
00398 
00399     side4.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00400                coordinates[3][2] - coordinates[1][2] );
00401 
00402     side5.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00403                coordinates[3][2] - coordinates[2][2] );
00404 
00405     double volume = fabs( v_tet_volume( 4, coordinates ) );
00406 
00407     if( volume < VERDICT_DBL_MIN )
00408         return (double)VERDICT_DBL_MAX;
00409     else
00410     {
00411         double srms = sqrt( ( side0.length_squared() + side1.length_squared() + side2.length_squared() +
00412                               side3.length_squared() + side4.length_squared() + side5.length_squared() ) /
00413                             6.0 );
00414 
00415         double aspect_ratio_gamma = pow( srms, 3 ) / ( 8.48528137423857 * volume );
00416         return (double)aspect_ratio_gamma;
00417     }
00418 }
00419 
00420 /*!
00421   The aspect frobenius of a tet
00422 
00423   NB (P. Pebay 01/22/07):
00424     Frobenius condition number when the reference element is regular
00425 */
00426 C_FUNC_DEF double v_tet_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
00427 {
00428     static const double normal_exp = 1. / 3.;
00429 
00430     VerdictVector ab, ac, ad;
00431 
00432     ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00433             coordinates[1][2] - coordinates[0][2] );
00434 
00435     ac.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00436             coordinates[2][2] - coordinates[0][2] );
00437 
00438     ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00439             coordinates[3][2] - coordinates[0][2] );
00440 
00441     double denominator = ab % ( ac * ad );
00442     denominator *= denominator;
00443     denominator *= 2.;
00444     denominator = 3. * pow( denominator, normal_exp );
00445 
00446     if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00447 
00448     double u[3];
00449     ab.get_xyz( u );
00450     double v[3];
00451     ac.get_xyz( v );
00452     double w[3];
00453     ad.get_xyz( w );
00454 
00455     double numerator = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
00456     numerator += v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
00457     numerator += w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
00458     numerator *= 1.5;
00459     numerator -= v[0] * u[0] + v[1] * u[1] + v[2] * u[2];
00460     numerator -= w[0] * u[0] + w[1] * u[1] + w[2] * u[2];
00461     numerator -= w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
00462 
00463     double aspect_frobenius = numerator / denominator;
00464 
00465     if( aspect_frobenius > 0 ) return (double)VERDICT_MIN( aspect_frobenius, VERDICT_DBL_MAX );
00466     return (double)VERDICT_MAX( aspect_frobenius, -VERDICT_DBL_MAX );
00467 }
00468 
00469 /*!
00470   The minimum angle of a tet
00471 
00472   NB (P. Pebay 01/22/07):
00473     minimum nonoriented dihedral angle
00474 */
00475 C_FUNC_DEF double v_tet_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
00476 {
00477     static const double normal_coeff = 180. * .3183098861837906715377675267450287;
00478 
00479     // Determine side vectors
00480     VerdictVector ab, bc, ad, cd;
00481 
00482     ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00483             coordinates[1][2] - coordinates[0][2] );
00484 
00485     ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00486             coordinates[3][2] - coordinates[0][2] );
00487 
00488     bc.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00489             coordinates[2][2] - coordinates[1][2] );
00490 
00491     cd.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00492             coordinates[3][2] - coordinates[2][2] );
00493 
00494     VerdictVector abc = ab * bc;
00495     double nabc       = abc.length();
00496     VerdictVector abd = ab * ad;
00497     double nabd       = abd.length();
00498     VerdictVector acd = ad * cd;
00499     double nacd       = acd.length();
00500     VerdictVector bcd = bc * cd;
00501     double nbcd       = bcd.length();
00502 
00503     double alpha   = acos( ( abc % abd ) / ( nabc * nabd ) );
00504     double beta    = acos( ( abc % acd ) / ( nabc * nacd ) );
00505     double gamma   = acos( ( abc % bcd ) / ( nabc * nbcd ) );
00506     double delta   = acos( ( abd % acd ) / ( nabd * nacd ) );
00507     double epsilon = acos( ( abd % bcd ) / ( nabd * nbcd ) );
00508     double zeta    = acos( ( acd % bcd ) / ( nacd * nbcd ) );
00509 
00510     alpha = alpha < beta ? alpha : beta;
00511     alpha = alpha < gamma ? alpha : gamma;
00512     alpha = alpha < delta ? alpha : delta;
00513     alpha = alpha < epsilon ? alpha : epsilon;
00514     alpha = alpha < zeta ? alpha : zeta;
00515     alpha *= normal_coeff;
00516 
00517     if( alpha < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00518 
00519     if( alpha > 0 ) return (double)VERDICT_MIN( alpha, VERDICT_DBL_MAX );
00520     return (double)VERDICT_MAX( alpha, -VERDICT_DBL_MAX );
00521 }
00522 
00523 /*!
00524   The collapse ratio of a tet
00525 */
00526 C_FUNC_DEF double v_tet_collapse_ratio( int /*num_nodes*/, double coordinates[][3] )
00527 {
00528     // Determine side vectors
00529     VerdictVector e01, e02, e03, e12, e13, e23;
00530 
00531     e01.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00532              coordinates[1][2] - coordinates[0][2] );
00533 
00534     e02.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00535              coordinates[2][2] - coordinates[0][2] );
00536 
00537     e03.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00538              coordinates[3][2] - coordinates[0][2] );
00539 
00540     e12.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00541              coordinates[2][2] - coordinates[1][2] );
00542 
00543     e13.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00544              coordinates[3][2] - coordinates[1][2] );
00545 
00546     e23.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00547              coordinates[3][2] - coordinates[2][2] );
00548 
00549     double l[6];
00550     l[0] = e01.length();
00551     l[1] = e02.length();
00552     l[2] = e03.length();
00553     l[3] = e12.length();
00554     l[4] = e13.length();
00555     l[5] = e23.length();
00556 
00557     // Find longest edge for each bounding triangle of tetrahedron
00558     double l012 = l[4] > l[0] ? l[4] : l[0];
00559     l012        = l[1] > l012 ? l[1] : l012;
00560     double l031 = l[0] > l[2] ? l[0] : l[2];
00561     l031        = l[3] > l031 ? l[3] : l031;
00562     double l023 = l[2] > l[1] ? l[2] : l[1];
00563     l023        = l[5] > l023 ? l[5] : l023;
00564     double l132 = l[4] > l[3] ? l[4] : l[3];
00565     l132        = l[5] > l132 ? l[5] : l132;
00566 
00567     // Compute collapse ratio for each vertex/triangle pair
00568     VerdictVector N;
00569     double h, magN;
00570     double cr;
00571     double crMin;
00572 
00573     N     = e01 * e02;
00574     magN  = N.length();
00575     h     = ( e03 % N ) / magN;  // height of vertex 3 above 0-1-2
00576     crMin = h / l012;            // ratio of height to longest edge of 0-1-2
00577 
00578     N    = e03 * e01;
00579     magN = N.length();
00580     h    = ( e02 % N ) / magN;  // height of vertex 2 above 0-3-1
00581     cr   = h / l031;            // ratio of height to longest edge of 0-3-1
00582     if( cr < crMin ) crMin = cr;
00583 
00584     N    = e02 * e03;
00585     magN = N.length();
00586     h    = ( e01 % N ) / magN;  // height of vertex 1 above 0-2-3
00587     cr   = h / l023;            // ratio of height to longest edge of 0-2-3
00588     if( cr < crMin ) crMin = cr;
00589 
00590     N    = e12 * e13;
00591     magN = N.length();
00592     h    = ( e01 % N ) / magN;  // height of vertex 0 above 1-3-2
00593     cr   = h / l132;            // ratio of height to longest edge of 1-3-2
00594     if( cr < crMin ) crMin = cr;
00595 
00596     if( crMin < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00597     if( crMin > 0 ) return (double)VERDICT_MIN( crMin, VERDICT_DBL_MAX );
00598     return (double)VERDICT_MAX( crMin, -VERDICT_DBL_MAX );
00599 }
00600 
00601 /*!
00602   the volume of a tet
00603 
00604   1/6 * jacobian at a corner node
00605 */
00606 C_FUNC_DEF double v_tet_volume( int /*num_nodes*/, double coordinates[][3] )
00607 {
00608 
00609     // Determine side vectors
00610     VerdictVector side0, side2, side3;
00611 
00612     side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00613                coordinates[1][2] - coordinates[0][2] );
00614 
00615     side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00616                coordinates[0][2] - coordinates[2][2] );
00617 
00618     side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00619                coordinates[3][2] - coordinates[0][2] );
00620 
00621     return (double)( ( side3 % ( side2 * side0 ) ) / 6.0 );
00622 }
00623 
00624 /*!
00625   the condition of a tet
00626 
00627   condition number of the jacobian matrix at any corner
00628 */
00629 C_FUNC_DEF double v_tet_condition( int /*num_nodes*/, double coordinates[][3] )
00630 {
00631 
00632     double condition, term1, term2, det;
00633     double rt3 = sqrt( 3.0 );
00634     double rt6 = sqrt( 6.0 );
00635 
00636     VerdictVector side0, side2, side3;
00637 
00638     side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00639                coordinates[1][2] - coordinates[0][2] );
00640 
00641     side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00642                coordinates[0][2] - coordinates[2][2] );
00643 
00644     side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00645                coordinates[3][2] - coordinates[0][2] );
00646 
00647     VerdictVector c_1, c_2, c_3;
00648 
00649     c_1 = side0;
00650     c_2 = ( -2 * side2 - side0 ) / rt3;
00651     c_3 = ( 3 * side3 + side2 - side0 ) / rt6;
00652 
00653     term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
00654     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 );
00655     det   = c_1 % ( c_2 * c_3 );
00656 
00657     if( det <= VERDICT_DBL_MIN )
00658         return VERDICT_DBL_MAX;
00659     else
00660         condition = sqrt( term1 * term2 ) / ( 3.0 * det );
00661 
00662     return (double)condition;
00663 }
00664 
00665 /*!
00666   the jacobian of a tet
00667 
00668   TODO
00669 */
00670 C_FUNC_DEF double v_tet_jacobian( int /*num_nodes*/, double coordinates[][3] )
00671 {
00672     VerdictVector side0, side2, side3;
00673 
00674     side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00675                coordinates[1][2] - coordinates[0][2] );
00676 
00677     side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00678                coordinates[0][2] - coordinates[2][2] );
00679 
00680     side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00681                coordinates[3][2] - coordinates[0][2] );
00682 
00683     return (double)( side3 % ( side2 * side0 ) );
00684 }
00685 
00686 /*!
00687   the shape of a tet
00688 
00689   3/ condition number of weighted jacobian matrix
00690 */
00691 C_FUNC_DEF double v_tet_shape( int /*num_nodes*/, double coordinates[][3] )
00692 {
00693 
00694     static const double two_thirds = 2.0 / 3.0;
00695     static const double root_of_2  = sqrt( 2.0 );
00696 
00697     VerdictVector edge0, edge2, edge3;
00698 
00699     edge0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00700                coordinates[1][2] - coordinates[0][2] );
00701 
00702     edge2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00703                coordinates[0][2] - coordinates[2][2] );
00704 
00705     edge3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00706                coordinates[3][2] - coordinates[0][2] );
00707 
00708     double jacobian = edge3 % ( edge2 * edge0 );
00709     if( jacobian < VERDICT_DBL_MIN )
00710     {
00711         return (double)0.0;
00712     }
00713     double num = 3 * pow( root_of_2 * jacobian, two_thirds );
00714     double den =
00715         1.5 * ( edge0 % edge0 + edge2 % edge2 + edge3 % edge3 ) - ( edge0 % -edge2 + -edge2 % edge3 + edge3 % edge0 );
00716 
00717     if( den < VERDICT_DBL_MIN ) return (double)0.0;
00718 
00719     return (double)VERDICT_MAX( num / den, 0 );
00720 }
00721 
00722 /*!
00723   the relative size of a tet
00724 
00725   Min(J,1/J), where J is the determinant of the weighted Jacobian matrix
00726 */
00727 C_FUNC_DEF double v_tet_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
00728 {
00729     double size;
00730     VerdictVector w1, w2, w3;
00731     get_weight( w1, w2, w3 );
00732     double avg_volume = ( w1 % ( w2 * w3 ) ) / 6.0;
00733 
00734     double volume = v_tet_volume( 4, coordinates );
00735 
00736     if( avg_volume < VERDICT_DBL_MIN )
00737         return 0.0;
00738     else
00739     {
00740         size = volume / avg_volume;
00741         if( size <= VERDICT_DBL_MIN ) return 0.0;
00742         if( size > 1 ) size = (double)( 1 ) / size;
00743     }
00744     return (double)( size * size );
00745 }
00746 
00747 /*!
00748   the shape and size of a tet
00749 
00750   Product of the shape and relative size
00751 */
00752 C_FUNC_DEF double v_tet_shape_and_size( int num_nodes, double coordinates[][3] )
00753 {
00754 
00755     double shape, size;
00756     shape = v_tet_shape( num_nodes, coordinates );
00757     size  = v_tet_relative_size_squared( num_nodes, coordinates );
00758 
00759     return (double)( shape * size );
00760 }
00761 
00762 /*!
00763   the distortion of a tet
00764 */
00765 C_FUNC_DEF double v_tet_distortion( int num_nodes, double coordinates[][3] )
00766 {
00767 
00768     double distortion          = VERDICT_DBL_MAX;
00769     int number_of_gauss_points = 0;
00770     if( num_nodes == 4 )
00771         // for linear tet, the distortion is always 1 because
00772         // straight edge tets are the target shape for tet
00773         return 1.0;
00774 
00775     else if( num_nodes == 10 )
00776         // use four integration points for quadratic tet
00777         number_of_gauss_points = 4;
00778 
00779     int number_dims                  = 3;
00780     int total_number_of_gauss_points = number_of_gauss_points;
00781     // use is_tri=1 to indicate this is for tet in 3D
00782     int is_tri = 1;
00783 
00784     double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
00785     double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
00786     double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
00787     double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
00788     double weight[maxTotalNumberGaussPoints];
00789 
00790     // create an object of GaussIntegration for tet
00791     GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
00792     GaussIntegration::calculate_shape_function_3d_tet();
00793     GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight );
00794 
00795     // vector xxi is the derivative vector of coordinates w.r.t local xi coordinate in the
00796     // computation space
00797     // vector xet is the derivative vector of coordinates w.r.t local et coordinate in the
00798     // computation space
00799     // vector xze is the derivative vector of coordinates w.r.t local ze coordinate in the
00800     // computation space
00801     VerdictVector xxi, xet, xze, xin;
00802 
00803     double jacobian, minimum_jacobian;
00804     double element_volume = 0.0;
00805     minimum_jacobian      = VERDICT_DBL_MAX;
00806 
00807     // calculate element volume
00808     int ife, ja;
00809     for( ife = 0; ife < total_number_of_gauss_points; ife++ )
00810     {
00811         xxi.set( 0.0, 0.0, 0.0 );
00812         xet.set( 0.0, 0.0, 0.0 );
00813         xze.set( 0.0, 0.0, 0.0 );
00814 
00815         for( ja = 0; ja < num_nodes; ja++ )
00816         {
00817             xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
00818             xxi += dndy1[ife][ja] * xin;
00819             xet += dndy2[ife][ja] * xin;
00820             xze += dndy3[ife][ja] * xin;
00821         }
00822 
00823         // determinant
00824         jacobian = xxi % ( xet * xze );
00825         if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
00826 
00827         element_volume += weight[ife] * jacobian;
00828     }  // element_volume is 6 times the actual volume
00829 
00830     // loop through all nodes
00831     double dndy1_at_node[maxNumberNodes][maxNumberNodes];
00832     double dndy2_at_node[maxNumberNodes][maxNumberNodes];
00833     double dndy3_at_node[maxNumberNodes][maxNumberNodes];
00834 
00835     GaussIntegration::calculate_derivative_at_nodes_3d_tet( dndy1_at_node, dndy2_at_node, dndy3_at_node );
00836     int node_id;
00837     for( node_id = 0; node_id < num_nodes; node_id++ )
00838     {
00839         xxi.set( 0.0, 0.0, 0.0 );
00840         xet.set( 0.0, 0.0, 0.0 );
00841         xze.set( 0.0, 0.0, 0.0 );
00842 
00843         for( ja = 0; ja < num_nodes; ja++ )
00844         {
00845             xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
00846             xxi += dndy1_at_node[node_id][ja] * xin;
00847             xet += dndy2_at_node[node_id][ja] * xin;
00848             xze += dndy3_at_node[node_id][ja] * xin;
00849         }
00850 
00851         jacobian = xxi % ( xet * xze );
00852         if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
00853     }
00854     distortion = minimum_jacobian / element_volume;
00855 
00856     return (double)distortion;
00857 }
00858 
00859 /*!
00860   the quality metrics of a tet
00861 */
00862 C_FUNC_DEF void v_tet_quality( int num_nodes,
00863                                double coordinates[][3],
00864                                unsigned int metrics_request_flag,
00865                                TetMetricVals* metric_vals )
00866 {
00867 
00868     memset( metric_vals, 0, sizeof( TetMetricVals ) );
00869 
00870     /*
00871 
00872       node numbers and edge numbers below
00873 
00874 
00875 
00876                3
00877                +            edge 0 is node 0 to 1
00878               +|+           edge 1 is node 1 to 2
00879             3/ | \5         edge 2 is node 0 to 2
00880             / 4|  \         edge 3 is node 0 to 3
00881           0 - -|- + 2       edge 4 is node 1 to 3
00882             \  |  +         edge 5 is node 2 to 3
00883             0\ | /1
00884               +|/           edge 2 is behind edge 4
00885                1
00886 
00887 
00888     */
00889 
00890     // lets start with making the vectors
00891     VerdictVector edges[6];
00892     edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00893                   coordinates[1][2] - coordinates[0][2] );
00894 
00895     edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00896                   coordinates[2][2] - coordinates[1][2] );
00897 
00898     edges[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00899                   coordinates[0][2] - coordinates[2][2] );
00900 
00901     edges[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00902                   coordinates[3][2] - coordinates[0][2] );
00903 
00904     edges[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00905                   coordinates[3][2] - coordinates[1][2] );
00906 
00907     edges[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00908                   coordinates[3][2] - coordinates[2][2] );
00909 
00910     // common numbers
00911     static const double root_of_2 = sqrt( 2.0 );
00912 
00913     // calculate the jacobian
00914     static const int do_jacobian = V_TET_JACOBIAN | V_TET_VOLUME | V_TET_ASPECT_BETA | V_TET_ASPECT_GAMMA |
00915                                    V_TET_SHAPE | V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE |
00916                                    V_TET_SCALED_JACOBIAN | V_TET_CONDITION;
00917     if( metrics_request_flag & do_jacobian )
00918     {
00919         metric_vals->jacobian = (double)( edges[3] % ( edges[2] * edges[0] ) );
00920     }
00921 
00922     // calculate the volume
00923     if( metrics_request_flag & V_TET_VOLUME )
00924     {
00925         metric_vals->volume = (double)( metric_vals->jacobian / 6.0 );
00926     }
00927 
00928     // calculate aspect ratio
00929     if( metrics_request_flag & V_TET_ASPECT_BETA )
00930     {
00931         double surface_area = ( ( edges[2] * edges[0] ).length() + ( edges[3] * edges[0] ).length() +
00932                                 ( edges[4] * edges[1] ).length() + ( edges[3] * edges[2] ).length() ) *
00933                               0.5;
00934 
00935         VerdictVector numerator = edges[3].length_squared() * ( edges[2] * edges[0] ) +
00936                                   edges[2].length_squared() * ( edges[3] * edges[0] ) +
00937                                   edges[0].length_squared() * ( edges[3] * edges[2] );
00938 
00939         double volume = metric_vals->jacobian / 6.0;
00940 
00941         if( volume < VERDICT_DBL_MIN )
00942             metric_vals->aspect_beta = (double)( VERDICT_DBL_MAX );
00943         else
00944             metric_vals->aspect_beta = (double)( numerator.length() * surface_area / ( 108 * volume * volume ) );
00945     }
00946 
00947     // calculate the aspect gamma
00948     if( metrics_request_flag & V_TET_ASPECT_GAMMA )
00949     {
00950         double volume = fabs( metric_vals->jacobian / 6.0 );
00951         if( fabs( volume ) < VERDICT_DBL_MIN )
00952             metric_vals->aspect_gamma = VERDICT_DBL_MAX;
00953         else
00954         {
00955             double srms = sqrt( ( edges[0].length_squared() + edges[1].length_squared() + edges[2].length_squared() +
00956                                   edges[3].length_squared() + edges[4].length_squared() + edges[5].length_squared() ) /
00957                                 6.0 );
00958 
00959             // cube the srms
00960             srms *= ( srms * srms );
00961             metric_vals->aspect_gamma = (double)( srms / ( 8.48528137423857 * volume ) );
00962         }
00963     }
00964 
00965     // calculate the shape of the tet
00966     if( metrics_request_flag & ( V_TET_SHAPE | V_TET_SHAPE_AND_SIZE ) )
00967     {
00968         // if the jacobian is non-positive, the shape is 0
00969         if( metric_vals->jacobian < VERDICT_DBL_MIN )
00970         {
00971             metric_vals->shape = (double)0.0;
00972         }
00973         else
00974         {
00975             static const double two_thirds = 2.0 / 3.0;
00976             double num                     = 3.0 * pow( root_of_2 * metric_vals->jacobian, two_thirds );
00977             double den                     = 1.5 * ( edges[0] % edges[0] + edges[2] % edges[2] + edges[3] % edges[3] ) -
00978                          ( edges[0] % -edges[2] + -edges[2] % edges[3] + edges[3] % edges[0] );
00979 
00980             if( den < VERDICT_DBL_MIN )
00981                 metric_vals->shape = (double)0.0;
00982             else
00983                 metric_vals->shape = (double)VERDICT_MAX( num / den, 0 );
00984         }
00985     }
00986 
00987     // calculate the relative size of the tet
00988     if( metrics_request_flag & ( V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE ) )
00989     {
00990         VerdictVector w1, w2, w3;
00991         get_weight( w1, w2, w3 );
00992         double avg_vol = ( w1 % ( w2 * w3 ) ) / 6;
00993 
00994         if( avg_vol < VERDICT_DBL_MIN )
00995             metric_vals->relative_size_squared = 0.0;
00996         else
00997         {
00998             double tmp = metric_vals->jacobian / ( 6 * avg_vol );
00999             if( tmp < VERDICT_DBL_MIN )
01000                 metric_vals->relative_size_squared = 0.0;
01001             else
01002             {
01003                 tmp *= tmp;
01004                 metric_vals->relative_size_squared = (double)VERDICT_MIN( tmp, 1 / tmp );
01005             }
01006         }
01007     }
01008 
01009     // calculate the shape and size
01010     if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
01011     {
01012         metric_vals->shape_and_size = (double)( metric_vals->shape * metric_vals->relative_size_squared );
01013     }
01014 
01015     // calculate the scaled jacobian
01016     if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
01017     {
01018         // find out which node the normalized jacobian can be calculated at
01019         // and it will be the smaller than at other nodes
01020         double length_squared[4] = { edges[0].length_squared() * edges[2].length_squared() * edges[3].length_squared(),
01021                                      edges[0].length_squared() * edges[1].length_squared() * edges[4].length_squared(),
01022                                      edges[1].length_squared() * edges[2].length_squared() * edges[5].length_squared(),
01023                                      edges[3].length_squared() * edges[4].length_squared() *
01024                                          edges[5].length_squared() };
01025 
01026         int which_node = 0;
01027         if( length_squared[1] > length_squared[which_node] ) which_node = 1;
01028         if( length_squared[2] > length_squared[which_node] ) which_node = 2;
01029         if( length_squared[3] > length_squared[which_node] ) which_node = 3;
01030 
01031         // find the scaled jacobian at this node
01032         double length_product = sqrt( length_squared[which_node] );
01033         if( length_product < fabs( metric_vals->jacobian ) ) length_product = fabs( metric_vals->jacobian );
01034 
01035         if( length_product < VERDICT_DBL_MIN )
01036             metric_vals->scaled_jacobian = (double)VERDICT_DBL_MAX;
01037         else
01038             metric_vals->scaled_jacobian = (double)( root_of_2 * metric_vals->jacobian / length_product );
01039     }
01040 
01041     // calculate the condition number
01042     if( metrics_request_flag & V_TET_CONDITION )
01043     {
01044         static const double root_of_3 = sqrt( 3.0 );
01045         static const double root_of_6 = sqrt( 6.0 );
01046 
01047         VerdictVector c_1, c_2, c_3;
01048         c_1 = edges[0];
01049         c_2 = ( -2 * edges[2] - edges[0] ) / root_of_3;
01050         c_3 = ( 3 * edges[3] + edges[2] - edges[0] ) / root_of_6;
01051 
01052         double term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
01053         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 );
01054 
01055         double det = c_1 % ( c_2 * c_3 );
01056 
01057         if( det <= VERDICT_DBL_MIN )
01058             metric_vals->condition = (double)VERDICT_DBL_MAX;
01059         else
01060             metric_vals->condition = (double)( sqrt( term1 * term2 ) / ( 3.0 * det ) );
01061     }
01062 
01063     // calculate the distortion
01064     if( metrics_request_flag & V_TET_DISTORTION )
01065     {
01066         metric_vals->distortion = v_tet_distortion( num_nodes, coordinates );
01067     }
01068 
01069     // check for overflow
01070     if( metrics_request_flag & V_TET_ASPECT_BETA )
01071     {
01072         if( metric_vals->aspect_beta > 0 )
01073             metric_vals->aspect_beta = (double)VERDICT_MIN( metric_vals->aspect_beta, VERDICT_DBL_MAX );
01074         metric_vals->aspect_beta = (double)VERDICT_MAX( metric_vals->aspect_beta, -VERDICT_DBL_MAX );
01075     }
01076 
01077     if( metrics_request_flag & V_TET_ASPECT_GAMMA )
01078     {
01079         if( metric_vals->aspect_gamma > 0 )
01080             metric_vals->aspect_gamma = (double)VERDICT_MIN( metric_vals->aspect_gamma, VERDICT_DBL_MAX );
01081         metric_vals->aspect_gamma = (double)VERDICT_MAX( metric_vals->aspect_gamma, -VERDICT_DBL_MAX );
01082     }
01083 
01084     if( metrics_request_flag & V_TET_VOLUME )
01085     {
01086         if( metric_vals->volume > 0 ) metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
01087         metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
01088     }
01089 
01090     if( metrics_request_flag & V_TET_CONDITION )
01091     {
01092         if( metric_vals->condition > 0 )
01093             metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
01094         metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
01095     }
01096 
01097     if( metrics_request_flag & V_TET_JACOBIAN )
01098     {
01099         if( metric_vals->jacobian > 0 )
01100             metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
01101         metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
01102     }
01103 
01104     if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
01105     {
01106         if( metric_vals->scaled_jacobian > 0 )
01107             metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
01108         metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
01109     }
01110 
01111     if( metrics_request_flag & V_TET_SHAPE )
01112     {
01113         if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
01114         metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
01115     }
01116 
01117     if( metrics_request_flag & V_TET_RELATIVE_SIZE_SQUARED )
01118     {
01119         if( metric_vals->relative_size_squared > 0 )
01120             metric_vals->relative_size_squared =
01121                 (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
01122         metric_vals->relative_size_squared =
01123             (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
01124     }
01125 
01126     if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
01127     {
01128         if( metric_vals->shape_and_size > 0 )
01129             metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
01130         metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
01131     }
01132 
01133     if( metrics_request_flag & V_TET_DISTORTION )
01134     {
01135         if( metric_vals->distortion > 0 )
01136             metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
01137         metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
01138     }
01139 
01140     if( metrics_request_flag & V_TET_ASPECT_RATIO ) metric_vals->aspect_ratio = v_tet_aspect_ratio( 4, coordinates );
01141 
01142     if( metrics_request_flag & V_TET_ASPECT_FROBENIUS )
01143         metric_vals->aspect_frobenius = v_tet_aspect_frobenius( 4, coordinates );
01144 
01145     if( metrics_request_flag & V_TET_MINIMUM_ANGLE ) metric_vals->minimum_angle = v_tet_minimum_angle( 4, coordinates );
01146 
01147     if( metrics_request_flag & V_TET_COLLAPSE_RATIO )
01148         metric_vals->collapse_ratio = v_tet_collapse_ratio( 4, coordinates );
01149 
01150     if( metrics_request_flag & V_TET_RADIUS_RATIO ) metric_vals->radius_ratio = v_tet_radius_ratio( 4, coordinates );
01151 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines