MOAB: Mesh Oriented datABase  (version 5.2.1)
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 ) { return (double)0.0; }
00710     double num = 3 * pow( root_of_2 * jacobian, two_thirds );
00711     double den =
00712         1.5 * ( edge0 % edge0 + edge2 % edge2 + edge3 % edge3 ) - ( edge0 % -edge2 + -edge2 % edge3 + edge3 % edge0 );
00713 
00714     if( den < VERDICT_DBL_MIN ) return (double)0.0;
00715 
00716     return (double)VERDICT_MAX( num / den, 0 );
00717 }
00718 
00719 /*!
00720   the relative size of a tet
00721 
00722   Min(J,1/J), where J is the determinant of the weighted Jacobian matrix
00723 */
00724 C_FUNC_DEF double v_tet_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
00725 {
00726     double size;
00727     VerdictVector w1, w2, w3;
00728     get_weight( w1, w2, w3 );
00729     double avg_volume = ( w1 % ( w2 * w3 ) ) / 6.0;
00730 
00731     double volume = v_tet_volume( 4, coordinates );
00732 
00733     if( avg_volume < VERDICT_DBL_MIN )
00734         return 0.0;
00735     else
00736     {
00737         size = volume / avg_volume;
00738         if( size <= VERDICT_DBL_MIN ) return 0.0;
00739         if( size > 1 ) size = (double)( 1 ) / size;
00740     }
00741     return (double)( size * size );
00742 }
00743 
00744 /*!
00745   the shape and size of a tet
00746 
00747   Product of the shape and relative size
00748 */
00749 C_FUNC_DEF double v_tet_shape_and_size( int num_nodes, double coordinates[][3] )
00750 {
00751 
00752     double shape, size;
00753     shape = v_tet_shape( num_nodes, coordinates );
00754     size  = v_tet_relative_size_squared( num_nodes, coordinates );
00755 
00756     return (double)( shape * size );
00757 }
00758 
00759 /*!
00760   the distortion of a tet
00761 */
00762 C_FUNC_DEF double v_tet_distortion( int num_nodes, double coordinates[][3] )
00763 {
00764 
00765     double distortion          = VERDICT_DBL_MAX;
00766     int number_of_gauss_points = 0;
00767     if( num_nodes == 4 )
00768         // for linear tet, the distortion is always 1 because
00769         // straight edge tets are the target shape for tet
00770         return 1.0;
00771 
00772     else if( num_nodes == 10 )
00773         // use four integration points for quadratic tet
00774         number_of_gauss_points = 4;
00775 
00776     int number_dims                  = 3;
00777     int total_number_of_gauss_points = number_of_gauss_points;
00778     // use is_tri=1 to indicate this is for tet in 3D
00779     int is_tri = 1;
00780 
00781     double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
00782     double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
00783     double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
00784     double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
00785     double weight[maxTotalNumberGaussPoints];
00786 
00787     // create an object of GaussIntegration for tet
00788     GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
00789     GaussIntegration::calculate_shape_function_3d_tet();
00790     GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight );
00791 
00792     // vector xxi is the derivative vector of coordinates w.r.t local xi coordinate in the
00793     // computation space
00794     // vector xet is the derivative vector of coordinates w.r.t local et coordinate in the
00795     // computation space
00796     // vector xze is the derivative vector of coordinates w.r.t local ze coordinate in the
00797     // computation space
00798     VerdictVector xxi, xet, xze, xin;
00799 
00800     double jacobian, minimum_jacobian;
00801     double element_volume = 0.0;
00802     minimum_jacobian      = VERDICT_DBL_MAX;
00803 
00804     // calculate element volume
00805     int ife, ja;
00806     for( ife = 0; ife < total_number_of_gauss_points; ife++ )
00807     {
00808         xxi.set( 0.0, 0.0, 0.0 );
00809         xet.set( 0.0, 0.0, 0.0 );
00810         xze.set( 0.0, 0.0, 0.0 );
00811 
00812         for( ja = 0; ja < num_nodes; ja++ )
00813         {
00814             xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
00815             xxi += dndy1[ife][ja] * xin;
00816             xet += dndy2[ife][ja] * xin;
00817             xze += dndy3[ife][ja] * xin;
00818         }
00819 
00820         // determinant
00821         jacobian = xxi % ( xet * xze );
00822         if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
00823 
00824         element_volume += weight[ife] * jacobian;
00825     }  // element_volume is 6 times the actual volume
00826 
00827     // loop through all nodes
00828     double dndy1_at_node[maxNumberNodes][maxNumberNodes];
00829     double dndy2_at_node[maxNumberNodes][maxNumberNodes];
00830     double dndy3_at_node[maxNumberNodes][maxNumberNodes];
00831 
00832     GaussIntegration::calculate_derivative_at_nodes_3d_tet( dndy1_at_node, dndy2_at_node, dndy3_at_node );
00833     int node_id;
00834     for( node_id = 0; node_id < num_nodes; node_id++ )
00835     {
00836         xxi.set( 0.0, 0.0, 0.0 );
00837         xet.set( 0.0, 0.0, 0.0 );
00838         xze.set( 0.0, 0.0, 0.0 );
00839 
00840         for( ja = 0; ja < num_nodes; ja++ )
00841         {
00842             xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
00843             xxi += dndy1_at_node[node_id][ja] * xin;
00844             xet += dndy2_at_node[node_id][ja] * xin;
00845             xze += dndy3_at_node[node_id][ja] * xin;
00846         }
00847 
00848         jacobian = xxi % ( xet * xze );
00849         if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
00850     }
00851     distortion = minimum_jacobian / element_volume;
00852 
00853     return (double)distortion;
00854 }
00855 
00856 /*!
00857   the quality metrics of a tet
00858 */
00859 C_FUNC_DEF void v_tet_quality( int num_nodes, double coordinates[][3], unsigned int metrics_request_flag,
00860                                TetMetricVals* metric_vals )
00861 {
00862 
00863     memset( metric_vals, 0, sizeof( TetMetricVals ) );
00864 
00865     /*
00866 
00867       node numbers and edge numbers below
00868 
00869 
00870 
00871                3
00872                +            edge 0 is node 0 to 1
00873               +|+           edge 1 is node 1 to 2
00874             3/ | \5         edge 2 is node 0 to 2
00875             / 4|  \         edge 3 is node 0 to 3
00876           0 - -|- + 2       edge 4 is node 1 to 3
00877             \  |  +         edge 5 is node 2 to 3
00878             0\ | /1
00879               +|/           edge 2 is behind edge 4
00880                1
00881 
00882 
00883     */
00884 
00885     // lets start with making the vectors
00886     VerdictVector edges[6];
00887     edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00888                   coordinates[1][2] - coordinates[0][2] );
00889 
00890     edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00891                   coordinates[2][2] - coordinates[1][2] );
00892 
00893     edges[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00894                   coordinates[0][2] - coordinates[2][2] );
00895 
00896     edges[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00897                   coordinates[3][2] - coordinates[0][2] );
00898 
00899     edges[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00900                   coordinates[3][2] - coordinates[1][2] );
00901 
00902     edges[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00903                   coordinates[3][2] - coordinates[2][2] );
00904 
00905     // common numbers
00906     static const double root_of_2 = sqrt( 2.0 );
00907 
00908     // calculate the jacobian
00909     static const int do_jacobian = V_TET_JACOBIAN | V_TET_VOLUME | V_TET_ASPECT_BETA | V_TET_ASPECT_GAMMA |
00910                                    V_TET_SHAPE | V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE |
00911                                    V_TET_SCALED_JACOBIAN | V_TET_CONDITION;
00912     if( metrics_request_flag & do_jacobian ) { metric_vals->jacobian = (double)( edges[3] % ( edges[2] * edges[0] ) ); }
00913 
00914     // calculate the volume
00915     if( metrics_request_flag & V_TET_VOLUME ) { metric_vals->volume = (double)( metric_vals->jacobian / 6.0 ); }
00916 
00917     // calculate aspect ratio
00918     if( metrics_request_flag & V_TET_ASPECT_BETA )
00919     {
00920         double surface_area = ( ( edges[2] * edges[0] ).length() + ( edges[3] * edges[0] ).length() +
00921                                 ( edges[4] * edges[1] ).length() + ( edges[3] * edges[2] ).length() ) *
00922                               0.5;
00923 
00924         VerdictVector numerator = edges[3].length_squared() * ( edges[2] * edges[0] ) +
00925                                   edges[2].length_squared() * ( edges[3] * edges[0] ) +
00926                                   edges[0].length_squared() * ( edges[3] * edges[2] );
00927 
00928         double volume = metric_vals->jacobian / 6.0;
00929 
00930         if( volume < VERDICT_DBL_MIN )
00931             metric_vals->aspect_beta = (double)( VERDICT_DBL_MAX );
00932         else
00933             metric_vals->aspect_beta = (double)( numerator.length() * surface_area / ( 108 * volume * volume ) );
00934     }
00935 
00936     // calculate the aspect gamma
00937     if( metrics_request_flag & V_TET_ASPECT_GAMMA )
00938     {
00939         double volume = fabs( metric_vals->jacobian / 6.0 );
00940         if( fabs( volume ) < VERDICT_DBL_MIN )
00941             metric_vals->aspect_gamma = VERDICT_DBL_MAX;
00942         else
00943         {
00944             double srms = sqrt( ( edges[0].length_squared() + edges[1].length_squared() + edges[2].length_squared() +
00945                                   edges[3].length_squared() + edges[4].length_squared() + edges[5].length_squared() ) /
00946                                 6.0 );
00947 
00948             // cube the srms
00949             srms *= ( srms * srms );
00950             metric_vals->aspect_gamma = (double)( srms / ( 8.48528137423857 * volume ) );
00951         }
00952     }
00953 
00954     // calculate the shape of the tet
00955     if( metrics_request_flag & ( V_TET_SHAPE | V_TET_SHAPE_AND_SIZE ) )
00956     {
00957         // if the jacobian is non-positive, the shape is 0
00958         if( metric_vals->jacobian < VERDICT_DBL_MIN ) { metric_vals->shape = (double)0.0; }
00959         else
00960         {
00961             static const double two_thirds = 2.0 / 3.0;
00962             double num                     = 3.0 * pow( root_of_2 * metric_vals->jacobian, two_thirds );
00963             double den                     = 1.5 * ( edges[0] % edges[0] + edges[2] % edges[2] + edges[3] % edges[3] ) -
00964                          ( edges[0] % -edges[2] + -edges[2] % edges[3] + edges[3] % edges[0] );
00965 
00966             if( den < VERDICT_DBL_MIN )
00967                 metric_vals->shape = (double)0.0;
00968             else
00969                 metric_vals->shape = (double)VERDICT_MAX( num / den, 0 );
00970         }
00971     }
00972 
00973     // calculate the relative size of the tet
00974     if( metrics_request_flag & ( V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE ) )
00975     {
00976         VerdictVector w1, w2, w3;
00977         get_weight( w1, w2, w3 );
00978         double avg_vol = ( w1 % ( w2 * w3 ) ) / 6;
00979 
00980         if( avg_vol < VERDICT_DBL_MIN )
00981             metric_vals->relative_size_squared = 0.0;
00982         else
00983         {
00984             double tmp = metric_vals->jacobian / ( 6 * avg_vol );
00985             if( tmp < VERDICT_DBL_MIN )
00986                 metric_vals->relative_size_squared = 0.0;
00987             else
00988             {
00989                 tmp *= tmp;
00990                 metric_vals->relative_size_squared = (double)VERDICT_MIN( tmp, 1 / tmp );
00991             }
00992         }
00993     }
00994 
00995     // calculate the shape and size
00996     if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
00997     { metric_vals->shape_and_size = (double)( metric_vals->shape * metric_vals->relative_size_squared ); }
00998 
00999     // calculate the scaled jacobian
01000     if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
01001     {
01002         // find out which node the normalized jacobian can be calculated at
01003         // and it will be the smaller than at other nodes
01004         double length_squared[4] = { edges[0].length_squared() * edges[2].length_squared() * edges[3].length_squared(),
01005                                      edges[0].length_squared() * edges[1].length_squared() * edges[4].length_squared(),
01006                                      edges[1].length_squared() * edges[2].length_squared() * edges[5].length_squared(),
01007                                      edges[3].length_squared() * edges[4].length_squared() *
01008                                          edges[5].length_squared() };
01009 
01010         int which_node = 0;
01011         if( length_squared[1] > length_squared[which_node] ) which_node = 1;
01012         if( length_squared[2] > length_squared[which_node] ) which_node = 2;
01013         if( length_squared[3] > length_squared[which_node] ) which_node = 3;
01014 
01015         // find the scaled jacobian at this node
01016         double length_product = sqrt( length_squared[which_node] );
01017         if( length_product < fabs( metric_vals->jacobian ) ) length_product = fabs( metric_vals->jacobian );
01018 
01019         if( length_product < VERDICT_DBL_MIN )
01020             metric_vals->scaled_jacobian = (double)VERDICT_DBL_MAX;
01021         else
01022             metric_vals->scaled_jacobian = (double)( root_of_2 * metric_vals->jacobian / length_product );
01023     }
01024 
01025     // calculate the condition number
01026     if( metrics_request_flag & V_TET_CONDITION )
01027     {
01028         static const double root_of_3 = sqrt( 3.0 );
01029         static const double root_of_6 = sqrt( 6.0 );
01030 
01031         VerdictVector c_1, c_2, c_3;
01032         c_1 = edges[0];
01033         c_2 = ( -2 * edges[2] - edges[0] ) / root_of_3;
01034         c_3 = ( 3 * edges[3] + edges[2] - edges[0] ) / root_of_6;
01035 
01036         double term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
01037         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 );
01038 
01039         double det = c_1 % ( c_2 * c_3 );
01040 
01041         if( det <= VERDICT_DBL_MIN )
01042             metric_vals->condition = (double)VERDICT_DBL_MAX;
01043         else
01044             metric_vals->condition = (double)( sqrt( term1 * term2 ) / ( 3.0 * det ) );
01045     }
01046 
01047     // calculate the distortion
01048     if( metrics_request_flag & V_TET_DISTORTION )
01049     { metric_vals->distortion = v_tet_distortion( num_nodes, coordinates ); }
01050 
01051     // check for overflow
01052     if( metrics_request_flag & V_TET_ASPECT_BETA )
01053     {
01054         if( metric_vals->aspect_beta > 0 )
01055             metric_vals->aspect_beta = (double)VERDICT_MIN( metric_vals->aspect_beta, VERDICT_DBL_MAX );
01056         metric_vals->aspect_beta = (double)VERDICT_MAX( metric_vals->aspect_beta, -VERDICT_DBL_MAX );
01057     }
01058 
01059     if( metrics_request_flag & V_TET_ASPECT_GAMMA )
01060     {
01061         if( metric_vals->aspect_gamma > 0 )
01062             metric_vals->aspect_gamma = (double)VERDICT_MIN( metric_vals->aspect_gamma, VERDICT_DBL_MAX );
01063         metric_vals->aspect_gamma = (double)VERDICT_MAX( metric_vals->aspect_gamma, -VERDICT_DBL_MAX );
01064     }
01065 
01066     if( metrics_request_flag & V_TET_VOLUME )
01067     {
01068         if( metric_vals->volume > 0 ) metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
01069         metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
01070     }
01071 
01072     if( metrics_request_flag & V_TET_CONDITION )
01073     {
01074         if( metric_vals->condition > 0 )
01075             metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
01076         metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
01077     }
01078 
01079     if( metrics_request_flag & V_TET_JACOBIAN )
01080     {
01081         if( metric_vals->jacobian > 0 )
01082             metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
01083         metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
01084     }
01085 
01086     if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
01087     {
01088         if( metric_vals->scaled_jacobian > 0 )
01089             metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
01090         metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
01091     }
01092 
01093     if( metrics_request_flag & V_TET_SHAPE )
01094     {
01095         if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
01096         metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
01097     }
01098 
01099     if( metrics_request_flag & V_TET_RELATIVE_SIZE_SQUARED )
01100     {
01101         if( metric_vals->relative_size_squared > 0 )
01102             metric_vals->relative_size_squared =
01103                 (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
01104         metric_vals->relative_size_squared =
01105             (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
01106     }
01107 
01108     if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
01109     {
01110         if( metric_vals->shape_and_size > 0 )
01111             metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
01112         metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
01113     }
01114 
01115     if( metrics_request_flag & V_TET_DISTORTION )
01116     {
01117         if( metric_vals->distortion > 0 )
01118             metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
01119         metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
01120     }
01121 
01122     if( metrics_request_flag & V_TET_ASPECT_RATIO ) metric_vals->aspect_ratio = v_tet_aspect_ratio( 4, coordinates );
01123 
01124     if( metrics_request_flag & V_TET_ASPECT_FROBENIUS )
01125         metric_vals->aspect_frobenius = v_tet_aspect_frobenius( 4, coordinates );
01126 
01127     if( metrics_request_flag & V_TET_MINIMUM_ANGLE ) metric_vals->minimum_angle = v_tet_minimum_angle( 4, coordinates );
01128 
01129     if( metrics_request_flag & V_TET_COLLAPSE_RATIO )
01130         metric_vals->collapse_ratio = v_tet_collapse_ratio( 4, coordinates );
01131 
01132     if( metrics_request_flag & V_TET_RADIUS_RATIO ) metric_vals->radius_ratio = v_tet_radius_ratio( 4, coordinates );
01133 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines