MOAB: Mesh Oriented datABase  (version 5.2.1)
V_TriMetric.cpp
Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Module:    $RCSfile: V_TriMetric.cpp,v $
00004 
00005   Copyright (c) 2007 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  * TriMetric.cpp contains quality calculations for Tris
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 "V_GaussIntegration.hpp"
00028 #include "VerdictVector.hpp"
00029 #include <memory.h>
00030 #include <stddef.h>
00031 
00032 // the average area of a tri
00033 static double verdict_tri_size      = 0;
00034 static ComputeNormal compute_normal = NULL;
00035 
00036 /*!
00037   get weights based on the average area of a set of
00038   tris
00039 */
00040 static int v_tri_get_weight( double& m11, double& m21, double& m12, double& m22 )
00041 {
00042     static const double rootOf3 = sqrt( 3.0 );
00043     m11                         = 1;
00044     m21                         = 0;
00045     m12                         = 0.5;
00046     m22                         = 0.5 * rootOf3;
00047     double scale                = sqrt( 2.0 * verdict_tri_size / ( m11 * m22 - m21 * m12 ) );
00048 
00049     m11 *= scale;
00050     m21 *= scale;
00051     m12 *= scale;
00052     m22 *= scale;
00053 
00054     return 1;
00055 }
00056 
00057 /*! sets the average area of a tri */
00058 C_FUNC_DEF void v_set_tri_size( double size )
00059 {
00060     verdict_tri_size = size;
00061 }
00062 
00063 C_FUNC_DEF void v_set_tri_normal_func( ComputeNormal func )
00064 {
00065     compute_normal = func;
00066 }
00067 
00068 /*!
00069    the edge ratio of a triangle
00070 
00071    NB (P. Pebay 01/14/07):
00072      Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
00073      minimum edge lengths
00074 
00075 */
00076 C_FUNC_DEF double v_tri_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
00077 {
00078 
00079     // three vectors for each side
00080     VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00081                      coordinates[1][2] - coordinates[0][2] );
00082 
00083     VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00084                      coordinates[2][2] - coordinates[1][2] );
00085 
00086     VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00087                      coordinates[0][2] - coordinates[2][2] );
00088 
00089     double a2 = a.length_squared();
00090     double b2 = b.length_squared();
00091     double c2 = c.length_squared();
00092 
00093     double m2, M2;
00094     if( a2 < b2 )
00095     {
00096         if( b2 < c2 )
00097         {
00098             m2 = a2;
00099             M2 = c2;
00100         }
00101         else  // b2 <= a2
00102         {
00103             if( a2 < c2 )
00104             {
00105                 m2 = a2;
00106                 M2 = b2;
00107             }
00108             else  // c2 <= a2
00109             {
00110                 m2 = c2;
00111                 M2 = b2;
00112             }
00113         }
00114     }
00115     else  // b2 <= a2
00116     {
00117         if( a2 < c2 )
00118         {
00119             m2 = b2;
00120             M2 = c2;
00121         }
00122         else  // c2 <= a2
00123         {
00124             if( b2 < c2 )
00125             {
00126                 m2 = b2;
00127                 M2 = a2;
00128             }
00129             else  // c2 <= b2
00130             {
00131                 m2 = c2;
00132                 M2 = a2;
00133             }
00134         }
00135     }
00136 
00137     if( m2 < VERDICT_DBL_MIN )
00138         return (double)VERDICT_DBL_MAX;
00139     else
00140     {
00141         double edge_ratio;
00142         edge_ratio = sqrt( M2 / m2 );
00143 
00144         if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
00145         return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
00146     }
00147 }
00148 
00149 /*!
00150    the aspect ratio of a triangle
00151 
00152    NB (P. Pebay 01/14/07):
00153      Hmax / ( 2.0 * sqrt(3.0) * IR) where Hmax is the maximum edge length
00154      and IR is the inradius
00155 
00156      note that previous incarnations of verdict used "v_tri_aspect_ratio" to denote
00157      what is now called "v_tri_aspect_frobenius"
00158 
00159 */
00160 C_FUNC_DEF double v_tri_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
00161 {
00162     static const double normal_coeff = sqrt( 3. ) / 6.;
00163 
00164     // three vectors for each side
00165     VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00166                      coordinates[1][2] - coordinates[0][2] );
00167 
00168     VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00169                      coordinates[2][2] - coordinates[1][2] );
00170 
00171     VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00172                      coordinates[0][2] - coordinates[2][2] );
00173 
00174     double a1 = a.length();
00175     double b1 = b.length();
00176     double c1 = c.length();
00177 
00178     double hm = a1 > b1 ? a1 : b1;
00179     hm        = hm > c1 ? hm : c1;
00180 
00181     VerdictVector ab   = a * b;
00182     double denominator = ab.length();
00183 
00184     if( denominator < VERDICT_DBL_MIN )
00185         return (double)VERDICT_DBL_MAX;
00186     else
00187     {
00188         double aspect_ratio;
00189         aspect_ratio = normal_coeff * hm * ( a1 + b1 + c1 ) / denominator;
00190 
00191         if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
00192         return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
00193     }
00194 }
00195 
00196 /*!
00197    the radius ratio of a triangle
00198 
00199    NB (P. Pebay 01/13/07):
00200      CR / (3.0*IR) where CR is the circumradius and IR is the inradius
00201 
00202      this quality metric is also known to VERDICT, for tetrahedral elements only,
00203      a the "aspect beta"
00204 
00205 */
00206 C_FUNC_DEF double v_tri_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
00207 {
00208 
00209     // three vectors for each side
00210     VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00211                      coordinates[1][2] - coordinates[0][2] );
00212 
00213     VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00214                      coordinates[2][2] - coordinates[1][2] );
00215 
00216     VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00217                      coordinates[0][2] - coordinates[2][2] );
00218 
00219     double a2 = a.length_squared();
00220     double b2 = b.length_squared();
00221     double c2 = c.length_squared();
00222 
00223     VerdictVector ab   = a * b;
00224     double denominator = ab.length_squared();
00225 
00226     if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00227 
00228     double radius_ratio;
00229     radius_ratio = .25 * a2 * b2 * c2 * ( a2 + b2 + c2 ) / denominator;
00230 
00231     if( radius_ratio > 0 ) return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
00232     return (double)VERDICT_MAX( radius_ratio, -VERDICT_DBL_MAX );
00233 }
00234 
00235 /*!
00236    the Frobenius aspect of a tri
00237 
00238    srms^2/(2 * sqrt(3.0) * area)
00239    where srms^2 is sum of the lengths squared
00240 
00241    NB (P. Pebay 01/14/07):
00242      this method was called "aspect ratio" in earlier incarnations of VERDICT
00243 
00244 */
00245 
00246 C_FUNC_DEF double v_tri_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
00247 {
00248     static const double two_times_root_of_3 = 2 * sqrt( 3.0 );
00249 
00250     // three vectors for each side
00251     VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00252                          coordinates[1][2] - coordinates[0][2] );
00253 
00254     VerdictVector side2( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00255                          coordinates[2][2] - coordinates[1][2] );
00256 
00257     VerdictVector side3( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00258                          coordinates[0][2] - coordinates[2][2] );
00259 
00260     // sum the lengths squared of each side
00261     double srms = ( side1.length_squared() + side2.length_squared() + side3.length_squared() );
00262 
00263     // find two times the area of the triangle by cross product
00264     double areaX2 = ( ( side1 * ( -side3 ) ).length() );
00265 
00266     if( areaX2 == 0.0 ) return (double)VERDICT_DBL_MAX;
00267 
00268     double aspect = (double)( srms / ( two_times_root_of_3 * ( areaX2 ) ) );
00269     if( aspect > 0 ) return (double)VERDICT_MIN( aspect, VERDICT_DBL_MAX );
00270     return (double)VERDICT_MAX( aspect, -VERDICT_DBL_MAX );
00271 }
00272 
00273 /*!
00274   The area of a tri
00275 
00276   0.5 * jacobian at a node
00277 */
00278 C_FUNC_DEF double v_tri_area( int /*num_nodes*/, double coordinates[][3] )
00279 {
00280     // two vectors for two sides
00281     VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00282                          coordinates[1][2] - coordinates[0][2] );
00283 
00284     VerdictVector side3( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00285                          coordinates[2][2] - coordinates[0][2] );
00286 
00287     // the cross product of the two vectors representing two sides of the
00288     // triangle
00289     VerdictVector tmp = side1 * side3;
00290 
00291     // return the magnitude of the vector divided by two
00292     double area = 0.5 * tmp.length();
00293     if( area > 0 ) return (double)VERDICT_MIN( area, VERDICT_DBL_MAX );
00294     return (double)VERDICT_MAX( area, -VERDICT_DBL_MAX );
00295 }
00296 
00297 /*!
00298   The minimum angle of a tri
00299 
00300   The minimum angle of a tri is the minimum angle between
00301   two adjacents sides out of all three corners of the triangle.
00302 */
00303 C_FUNC_DEF double v_tri_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
00304 {
00305 
00306     // vectors for all the sides
00307     VerdictVector sides[4];
00308     sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00309                   coordinates[1][2] - coordinates[0][2] );
00310     sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00311                   coordinates[2][2] - coordinates[1][2] );
00312     sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00313                   coordinates[2][2] - coordinates[0][2] );
00314 
00315     // in case we need to find the interior angle
00316     // between sides 0 and 1
00317     sides[3] = -sides[1];
00318 
00319     // calculate the lengths squared of the sides
00320     double sides_lengths[3];
00321     sides_lengths[0] = sides[0].length_squared();
00322     sides_lengths[1] = sides[1].length_squared();
00323     sides_lengths[2] = sides[2].length_squared();
00324 
00325     if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 ) return 0.0;
00326 
00327     // using the law of sines, we know that the minimum
00328     // angle is opposite of the shortest side
00329 
00330     // find the shortest side
00331     int short_side = 0;
00332     if( sides_lengths[1] < sides_lengths[0] ) short_side = 1;
00333     if( sides_lengths[2] < sides_lengths[short_side] ) short_side = 2;
00334 
00335     // from the shortest side, calculate the angle of the
00336     // opposite angle
00337     double min_angle;
00338     if( short_side == 0 ) { min_angle = sides[2].interior_angle( sides[1] ); }
00339     else if( short_side == 1 )
00340     {
00341         min_angle = sides[0].interior_angle( sides[2] );
00342     }
00343     else
00344     {
00345         min_angle = sides[0].interior_angle( sides[3] );
00346     }
00347 
00348     if( min_angle > 0 ) return (double)VERDICT_MIN( min_angle, VERDICT_DBL_MAX );
00349     return (double)VERDICT_MAX( min_angle, -VERDICT_DBL_MAX );
00350 }
00351 
00352 /*!
00353   The maximum angle of a tri
00354 
00355   The maximum angle of a tri is the maximum angle between
00356   two adjacents sides out of all three corners of the triangle.
00357 */
00358 C_FUNC_DEF double v_tri_maximum_angle( int /*num_nodes*/, double coordinates[][3] )
00359 {
00360 
00361     // vectors for all the sides
00362     VerdictVector sides[4];
00363     sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00364                   coordinates[1][2] - coordinates[0][2] );
00365     sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00366                   coordinates[2][2] - coordinates[1][2] );
00367     sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00368                   coordinates[2][2] - coordinates[0][2] );
00369 
00370     // in case we need to find the interior angle
00371     // between sides 0 and 1
00372     sides[3] = -sides[1];
00373 
00374     // calculate the lengths squared of the sides
00375     double sides_lengths[3];
00376     sides_lengths[0] = sides[0].length_squared();
00377     sides_lengths[1] = sides[1].length_squared();
00378     sides_lengths[2] = sides[2].length_squared();
00379 
00380     if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 ) { return 0.0; }
00381 
00382     // using the law of sines, we know that the maximum
00383     // angle is opposite of the longest side
00384 
00385     // find the longest side
00386     int short_side = 0;
00387     if( sides_lengths[1] > sides_lengths[0] ) short_side = 1;
00388     if( sides_lengths[2] > sides_lengths[short_side] ) short_side = 2;
00389 
00390     // from the longest side, calculate the angle of the
00391     // opposite angle
00392     double max_angle;
00393     if( short_side == 0 ) { max_angle = sides[2].interior_angle( sides[1] ); }
00394     else if( short_side == 1 )
00395     {
00396         max_angle = sides[0].interior_angle( sides[2] );
00397     }
00398     else
00399     {
00400         max_angle = sides[0].interior_angle( sides[3] );
00401     }
00402 
00403     if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
00404     return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
00405 }
00406 
00407 /*!
00408   The condition of a tri
00409 
00410   Condition number of the jacobian matrix at any corner
00411 */
00412 C_FUNC_DEF double v_tri_condition( int /*num_nodes*/, double coordinates[][3] )
00413 {
00414     static const double rt3 = sqrt( 3.0 );
00415 
00416     VerdictVector v1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00417                       coordinates[1][2] - coordinates[0][2] );
00418 
00419     VerdictVector v2( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00420                       coordinates[2][2] - coordinates[0][2] );
00421 
00422     VerdictVector tri_normal = v1 * v2;
00423     double areax2            = tri_normal.length();
00424 
00425     if( areax2 == 0.0 ) return (double)VERDICT_DBL_MAX;
00426 
00427     double condition = (double)( ( ( v1 % v1 ) + ( v2 % v2 ) - ( v1 % v2 ) ) / ( areax2 * rt3 ) );
00428 
00429     // check for inverted if we have access to the normal
00430     if( compute_normal )
00431     {
00432         // center of tri
00433         double point[3], surf_normal[3];
00434         point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
00435         point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
00436         point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
00437 
00438         // dot product
00439         compute_normal( point, surf_normal );
00440         if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
00441             0 )
00442             return (double)VERDICT_DBL_MAX;
00443     }
00444     return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
00445 }
00446 
00447 /*!
00448   The scaled jacobian of a tri
00449 
00450   minimum of the jacobian divided by the lengths of 2 edge vectors
00451 */
00452 C_FUNC_DEF double v_tri_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
00453 {
00454     static const double detw = 2. / sqrt( 3.0 );
00455     VerdictVector first, second;
00456     double jacobian;
00457 
00458     VerdictVector edge[3];
00459     edge[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00460                  coordinates[1][2] - coordinates[0][2] );
00461 
00462     edge[1].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00463                  coordinates[2][2] - coordinates[0][2] );
00464 
00465     edge[2].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00466                  coordinates[2][2] - coordinates[1][2] );
00467     first  = edge[1] - edge[0];
00468     second = edge[2] - edge[0];
00469 
00470     VerdictVector cross = first * second;
00471     jacobian            = cross.length();
00472 
00473     double max_edge_length_product;
00474     max_edge_length_product =
00475         VERDICT_MAX( edge[0].length() * edge[1].length(),
00476                      VERDICT_MAX( edge[1].length() * edge[2].length(), edge[0].length() * edge[2].length() ) );
00477 
00478     if( max_edge_length_product < VERDICT_DBL_MIN ) return (double)0.0;
00479 
00480     jacobian *= detw;
00481     jacobian /= max_edge_length_product;
00482 
00483     if( compute_normal )
00484     {
00485         // center of tri
00486         double point[3], surf_normal[3];
00487         point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
00488         point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
00489         point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
00490 
00491         // dot product
00492         compute_normal( point, surf_normal );
00493         if( ( cross.x() * surf_normal[0] + cross.y() * surf_normal[1] + cross.z() * surf_normal[2] ) < 0 )
00494             jacobian *= -1;
00495     }
00496 
00497     if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
00498     return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
00499 }
00500 
00501 /*!
00502   The shape of a tri
00503 
00504   2 / condition number of weighted jacobian matrix
00505 */
00506 C_FUNC_DEF double v_tri_shape( int num_nodes, double coordinates[][3] )
00507 {
00508     double condition = v_tri_condition( num_nodes, coordinates );
00509 
00510     double shape;
00511     if( condition <= VERDICT_DBL_MIN )
00512         shape = VERDICT_DBL_MAX;
00513     else
00514         shape = ( 1 / condition );
00515 
00516     if( shape > 0 ) return (double)VERDICT_MIN( shape, VERDICT_DBL_MAX );
00517     return (double)VERDICT_MAX( shape, -VERDICT_DBL_MAX );
00518 }
00519 
00520 /*!
00521   The relative size of a tri
00522 
00523   Min(J,1/J) where J is the determinant of the weighted jacobian matrix.
00524 */
00525 C_FUNC_DEF double v_tri_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
00526 {
00527     double w11, w21, w12, w22;
00528 
00529     VerdictVector xxi, xet, tri_normal;
00530 
00531     v_tri_get_weight( w11, w21, w12, w22 );
00532 
00533     double detw = determinant( w11, w21, w12, w22 );
00534 
00535     if( detw == 0.0 ) return 0.0;
00536 
00537     xxi.set( coordinates[0][0] - coordinates[1][0], coordinates[0][1] - coordinates[1][1],
00538              coordinates[0][2] - coordinates[1][2] );
00539 
00540     xet.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00541              coordinates[0][2] - coordinates[2][2] );
00542 
00543     tri_normal = xxi * xet;
00544 
00545     double deta = tri_normal.length();
00546     if( deta == 0.0 || detw == 0.0 ) return 0.0;
00547 
00548     double size = pow( deta / detw, 2 );
00549 
00550     double rel_size = VERDICT_MIN( size, 1.0 / size );
00551 
00552     if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
00553     return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
00554 }
00555 
00556 /*!
00557   The shape and size of a tri
00558 
00559   Product of the Shape and Relative Size
00560 */
00561 C_FUNC_DEF double v_tri_shape_and_size( int num_nodes, double coordinates[][3] )
00562 {
00563     double size, shape;
00564 
00565     size  = v_tri_relative_size_squared( num_nodes, coordinates );
00566     shape = v_tri_shape( num_nodes, coordinates );
00567 
00568     double shape_and_size = size * shape;
00569 
00570     if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
00571     return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
00572 }
00573 
00574 /*!
00575   The distortion of a tri
00576 
00577 TODO:  make a short definition of the distortion and comment below
00578 */
00579 C_FUNC_DEF double v_tri_distortion( int num_nodes, double coordinates[][3] )
00580 {
00581 
00582     double distortion;
00583     int total_number_of_gauss_points = 0;
00584     VerdictVector aa, bb, cc, normal_at_point, xin;
00585     double element_area = 0.;
00586 
00587     aa.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00588             coordinates[1][2] - coordinates[0][2] );
00589 
00590     bb.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00591             coordinates[2][2] - coordinates[0][2] );
00592 
00593     VerdictVector tri_normal = aa * bb;
00594 
00595     int number_of_gauss_points = 0;
00596     if( num_nodes == 3 )
00597     {
00598         distortion = 1.0;
00599         return (double)distortion;
00600     }
00601 
00602     else if( num_nodes == 6 )
00603     {
00604         total_number_of_gauss_points = 6;
00605         number_of_gauss_points       = 6;
00606     }
00607 
00608     distortion = VERDICT_DBL_MAX;
00609     double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
00610     double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
00611     double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
00612     double weight[maxTotalNumberGaussPoints];
00613 
00614     // create an object of GaussIntegration
00615     int number_dims = 2;
00616     int is_tri      = 1;
00617     GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
00618     GaussIntegration::calculate_shape_function_2d_tri();
00619     GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
00620 
00621     // calculate element area
00622     int ife, ja;
00623     for( ife = 0; ife < total_number_of_gauss_points; ife++ )
00624     {
00625         aa.set( 0.0, 0.0, 0.0 );
00626         bb.set( 0.0, 0.0, 0.0 );
00627 
00628         for( ja = 0; ja < num_nodes; ja++ )
00629         {
00630             xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
00631             aa += dndy1[ife][ja] * xin;
00632             bb += dndy2[ife][ja] * xin;
00633         }
00634         normal_at_point = aa * bb;
00635         double jacobian = normal_at_point.length();
00636         element_area += weight[ife] * jacobian;
00637     }
00638 
00639     element_area *= 0.8660254;
00640     double dndy1_at_node[maxNumberNodes][maxNumberNodes];
00641     double dndy2_at_node[maxNumberNodes][maxNumberNodes];
00642 
00643     GaussIntegration::calculate_derivative_at_nodes_2d_tri( dndy1_at_node, dndy2_at_node );
00644 
00645     VerdictVector normal_at_nodes[7];
00646 
00647     // evaluate normal at nodes and distortion values at nodes
00648     int jai = 0;
00649     for( ja = 0; ja < num_nodes; ja++ )
00650     {
00651         aa.set( 0.0, 0.0, 0.0 );
00652         bb.set( 0.0, 0.0, 0.0 );
00653         for( jai = 0; jai < num_nodes; jai++ )
00654         {
00655             xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
00656             aa += dndy1_at_node[ja][jai] * xin;
00657             bb += dndy2_at_node[ja][jai] * xin;
00658         }
00659         normal_at_nodes[ja] = aa * bb;
00660         normal_at_nodes[ja].normalize();
00661     }
00662 
00663     // determine if element is flat
00664     bool flat_element = true;
00665     double dot_product;
00666 
00667     for( ja = 0; ja < num_nodes; ja++ )
00668     {
00669         dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
00670         if( fabs( dot_product ) < 0.99 )
00671         {
00672             flat_element = false;
00673             break;
00674         }
00675     }
00676 
00677     // take into consideration of the thickness of the element
00678     double thickness, thickness_gauss;
00679     double distrt;
00680     // get_tri_thickness(tri, element_area, thickness );
00681     thickness = 0.001 * sqrt( element_area );
00682 
00683     // set thickness gauss point location
00684     double zl = 0.5773502691896;
00685     if( flat_element ) zl = 0.0;
00686 
00687     int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
00688     double thickness_z;
00689 
00690     // loop on integration points
00691     int igz;
00692     for( ife = 0; ife < total_number_of_gauss_points; ife++ )
00693     {
00694         // loop on the thickness direction gauss points
00695         for( igz = 0; igz < no_gauss_pts_z; igz++ )
00696         {
00697             zl          = -zl;
00698             thickness_z = zl * thickness / 2.0;
00699 
00700             aa.set( 0.0, 0.0, 0.0 );
00701             bb.set( 0.0, 0.0, 0.0 );
00702             cc.set( 0.0, 0.0, 0.0 );
00703 
00704             for( ja = 0; ja < num_nodes; ja++ )
00705             {
00706                 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
00707                 xin += thickness_z * normal_at_nodes[ja];
00708                 aa += dndy1[ife][ja] * xin;
00709                 bb += dndy2[ife][ja] * xin;
00710                 thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
00711                 cc += thickness_gauss * normal_at_nodes[ja];
00712             }
00713 
00714             normal_at_point = aa * bb;
00715             distrt          = cc % normal_at_point;
00716             if( distrt < distortion ) distortion = distrt;
00717         }
00718     }
00719 
00720     // loop through nodal points
00721     for( ja = 0; ja < num_nodes; ja++ )
00722     {
00723         for( igz = 0; igz < no_gauss_pts_z; igz++ )
00724         {
00725             zl          = -zl;
00726             thickness_z = zl * thickness / 2.0;
00727 
00728             aa.set( 0.0, 0.0, 0.0 );
00729             bb.set( 0.0, 0.0, 0.0 );
00730             cc.set( 0.0, 0.0, 0.0 );
00731 
00732             for( jai = 0; jai < num_nodes; jai++ )
00733             {
00734                 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
00735                 xin += thickness_z * normal_at_nodes[ja];
00736                 aa += dndy1_at_node[ja][jai] * xin;
00737                 bb += dndy2_at_node[ja][jai] * xin;
00738                 if( jai == ja )
00739                     thickness_gauss = thickness / 2.0;
00740                 else
00741                     thickness_gauss = 0.;
00742                 cc += thickness_gauss * normal_at_nodes[jai];
00743             }
00744         }
00745 
00746         normal_at_point      = aa * bb;
00747         double sign_jacobian = ( tri_normal % normal_at_point ) > 0 ? 1. : -1.;
00748         distrt               = sign_jacobian * ( cc % normal_at_point );
00749 
00750         if( distrt < distortion ) distortion = distrt;
00751     }
00752     if( element_area * thickness != 0 )
00753         distortion *= 1. / ( element_area * thickness );
00754     else
00755         distortion *= 1.;
00756 
00757     if( distortion > 0 ) return (double)VERDICT_MIN( distortion, VERDICT_DBL_MAX );
00758     return (double)VERDICT_MAX( distortion, -VERDICT_DBL_MAX );
00759 }
00760 
00761 /*!
00762   tri_quality for calculating multiple tri metrics at once
00763 
00764   using this method is generally faster than using the individual
00765   method multiple times.
00766 
00767 */
00768 C_FUNC_DEF void v_tri_quality( int num_nodes, double coordinates[][3], unsigned int metrics_request_flag,
00769                                TriMetricVals* metric_vals )
00770 {
00771 
00772     memset( metric_vals, 0, sizeof( TriMetricVals ) );
00773 
00774     // for starts, lets set up some basic and common information
00775 
00776     /*  node numbers and side numbers used below
00777 
00778                2
00779                ++
00780               /  \
00781            2 /    \ 1
00782             /      \
00783            /        \
00784          0 ---------+ 1
00785                0
00786     */
00787 
00788     // vectors for each side
00789     VerdictVector sides[3];
00790     sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00791                   coordinates[1][2] - coordinates[0][2] );
00792     sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00793                   coordinates[2][2] - coordinates[1][2] );
00794     sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00795                   coordinates[2][2] - coordinates[0][2] );
00796     VerdictVector tri_normal = sides[0] * sides[2];
00797     // if we have access to normal information, check to see if the
00798     // element is inverted.  If we don't have the normal information
00799     // that we need for this, assume the element is not inverted.
00800     // This flag will be used for condition number, jacobian, shape,
00801     // and size and shape.
00802     bool is_inverted = false;
00803     if( compute_normal )
00804     {
00805         // center of tri
00806         double point[3], surf_normal[3];
00807         point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
00808         point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
00809         point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
00810         // dot product
00811         compute_normal( point, surf_normal );
00812         if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
00813             0 )
00814             is_inverted = true;
00815     }
00816 
00817     // lengths squared of each side
00818     double sides_lengths_squared[3];
00819     sides_lengths_squared[0] = sides[0].length_squared();
00820     sides_lengths_squared[1] = sides[1].length_squared();
00821     sides_lengths_squared[2] = sides[2].length_squared();
00822 
00823     // if we are doing angle calcuations
00824     if( metrics_request_flag & ( V_TRI_MINIMUM_ANGLE | V_TRI_MAXIMUM_ANGLE ) )
00825     {
00826         // which is short and long side
00827         int short_side = 0, long_side = 0;
00828 
00829         if( sides_lengths_squared[1] < sides_lengths_squared[0] ) short_side = 1;
00830         if( sides_lengths_squared[2] < sides_lengths_squared[short_side] ) short_side = 2;
00831 
00832         if( sides_lengths_squared[1] > sides_lengths_squared[0] ) long_side = 1;
00833         if( sides_lengths_squared[2] > sides_lengths_squared[long_side] ) long_side = 2;
00834 
00835         // calculate the minimum angle of the tri
00836         if( metrics_request_flag & V_TRI_MINIMUM_ANGLE )
00837         {
00838             if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
00839             { metric_vals->minimum_angle = 0.0; }
00840             else if( short_side == 0 )
00841                 metric_vals->minimum_angle = (double)sides[2].interior_angle( sides[1] );
00842             else if( short_side == 1 )
00843                 metric_vals->minimum_angle = (double)sides[0].interior_angle( sides[2] );
00844             else
00845                 metric_vals->minimum_angle = (double)sides[0].interior_angle( -sides[1] );
00846         }
00847 
00848         // calculate the maximum angle of the tri
00849         if( metrics_request_flag & V_TRI_MAXIMUM_ANGLE )
00850         {
00851             if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
00852             { metric_vals->maximum_angle = 0.0; }
00853             else if( long_side == 0 )
00854                 metric_vals->maximum_angle = (double)sides[2].interior_angle( sides[1] );
00855             else if( long_side == 1 )
00856                 metric_vals->maximum_angle = (double)sides[0].interior_angle( sides[2] );
00857             else
00858                 metric_vals->maximum_angle = (double)sides[0].interior_angle( -sides[1] );
00859         }
00860     }
00861 
00862     // calculate the area of the tri
00863     // the following metrics depend on area
00864     if( metrics_request_flag &
00865         ( V_TRI_AREA | V_TRI_SCALED_JACOBIAN | V_TRI_SHAPE | V_TRI_RELATIVE_SIZE_SQUARED | V_TRI_SHAPE_AND_SIZE ) )
00866     { metric_vals->area = (double)( ( sides[0] * sides[2] ).length() * 0.5 ); }
00867 
00868     // calculate the aspect ratio
00869     if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )
00870     {
00871         // sum the lengths squared
00872         double srms = sides_lengths_squared[0] + sides_lengths_squared[1] + sides_lengths_squared[2];
00873 
00874         // calculate once and reuse
00875         static const double twoTimesRootOf3 = 2 * sqrt( 3.0 );
00876 
00877         double div = ( twoTimesRootOf3 * ( ( sides[0] * sides[2] ).length() ) );
00878 
00879         if( div == 0.0 )
00880             metric_vals->aspect_frobenius = (double)VERDICT_DBL_MAX;
00881         else
00882             metric_vals->aspect_frobenius = (double)( srms / div );
00883     }
00884 
00885     // calculate the scaled jacobian
00886     if( metrics_request_flag & V_TRI_SCALED_JACOBIAN )
00887     {
00888         // calculate once and reuse
00889         static const double twoOverRootOf3 = 2 / sqrt( 3.0 );
00890         // use the area from above
00891 
00892         double tmp = tri_normal.length() * twoOverRootOf3;
00893 
00894         // now scale it by the lengths of the sides
00895         double min_scaled_jac = VERDICT_DBL_MAX;
00896         double temp_scaled_jac;
00897         for( int i = 0; i < 3; i++ )
00898         {
00899             if( sides_lengths_squared[i % 3] == 0.0 || sides_lengths_squared[( i + 2 ) % 3] == 0.0 )
00900                 temp_scaled_jac = 0.0;
00901             else
00902                 temp_scaled_jac =
00903                     tmp / sqrt( sides_lengths_squared[i % 3] ) / sqrt( sides_lengths_squared[( i + 2 ) % 3] );
00904             if( temp_scaled_jac < min_scaled_jac ) min_scaled_jac = temp_scaled_jac;
00905         }
00906         // multiply by -1 if the normals are in opposite directions
00907         if( is_inverted ) { min_scaled_jac *= -1; }
00908         metric_vals->scaled_jacobian = (double)min_scaled_jac;
00909     }
00910 
00911     // calculate the condition number
00912     if( metrics_request_flag & V_TRI_CONDITION )
00913     {
00914         // calculate once and reuse
00915         static const double rootOf3 = sqrt( 3.0 );
00916         // if it is inverted, the condition number is considered to be infinity.
00917         if( is_inverted ) { metric_vals->condition = VERDICT_DBL_MAX; }
00918         else
00919         {
00920             double area2x = ( sides[0] * sides[2] ).length();
00921             if( area2x == 0.0 )
00922                 metric_vals->condition = (double)( VERDICT_DBL_MAX );
00923             else
00924                 metric_vals->condition = (double)( ( sides[0] % sides[0] + sides[2] % sides[2] - sides[0] % sides[2] ) /
00925                                                    ( area2x * rootOf3 ) );
00926         }
00927     }
00928 
00929     // calculate the shape
00930     if( metrics_request_flag & V_TRI_SHAPE || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
00931     {
00932         // if element is inverted, shape is zero.  We don't need to
00933         // calculate anything.
00934         if( is_inverted ) { metric_vals->shape = 0.0; }
00935         else
00936         {  // otherwise, we calculate the shape
00937             // calculate once and reuse
00938             static const double rootOf3 = sqrt( 3.0 );
00939             // reuse area from before
00940             double area2x = metric_vals->area * 2;
00941             // dot products
00942             double dots[3] = { sides[0] % sides[0], sides[2] % sides[2], sides[0] % sides[2] };
00943 
00944             // add the dots
00945             double sum_dots = dots[0] + dots[1] - dots[2];
00946 
00947             // then the finale
00948             if( sum_dots == 0.0 )
00949                 metric_vals->shape = 0.0;
00950             else
00951                 metric_vals->shape = (double)( rootOf3 * area2x / sum_dots );
00952         }
00953     }
00954 
00955     // calculate relative size squared
00956     if( metrics_request_flag & V_TRI_RELATIVE_SIZE_SQUARED || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
00957     {
00958         // get weights
00959         double w11, w21, w12, w22;
00960         v_tri_get_weight( w11, w21, w12, w22 );
00961         // get the determinant
00962         double detw = determinant( w11, w21, w12, w22 );
00963         // use the area from above and divide with the determinant
00964         if( metric_vals->area == 0.0 || detw == 0.0 )
00965             metric_vals->relative_size_squared = 0.0;
00966         else
00967         {
00968             double size = metric_vals->area * 2.0 / detw;
00969             // square the size
00970             size *= size;
00971             // value ranges between 0 to 1
00972             metric_vals->relative_size_squared = (double)VERDICT_MIN( size, 1.0 / size );
00973         }
00974     }
00975 
00976     // calculate shape and size
00977     if( metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
00978     { metric_vals->shape_and_size = metric_vals->relative_size_squared * metric_vals->shape; }
00979 
00980     // calculate distortion
00981     if( metrics_request_flag & V_TRI_DISTORTION ) metric_vals->distortion = v_tri_distortion( num_nodes, coordinates );
00982 
00983     // take care of any over-flow problems
00984     if( metric_vals->aspect_frobenius > 0 )
00985         metric_vals->aspect_frobenius = (double)VERDICT_MIN( metric_vals->aspect_frobenius, VERDICT_DBL_MAX );
00986     else
00987         metric_vals->aspect_frobenius = (double)VERDICT_MAX( metric_vals->aspect_frobenius, -VERDICT_DBL_MAX );
00988 
00989     if( metric_vals->area > 0 )
00990         metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
00991     else
00992         metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
00993 
00994     if( metric_vals->minimum_angle > 0 )
00995         metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
00996     else
00997         metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
00998 
00999     if( metric_vals->maximum_angle > 0 )
01000         metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
01001     else
01002         metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
01003 
01004     if( metric_vals->condition > 0 )
01005         metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
01006     else
01007         metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
01008 
01009     if( metric_vals->shape > 0 )
01010         metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
01011     else
01012         metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
01013 
01014     if( metric_vals->scaled_jacobian > 0 )
01015         metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
01016     else
01017         metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
01018 
01019     if( metric_vals->relative_size_squared > 0 )
01020         metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
01021     else
01022         metric_vals->relative_size_squared =
01023             (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
01024 
01025     if( metric_vals->shape_and_size > 0 )
01026         metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
01027     else
01028         metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
01029 
01030     if( metric_vals->distortion > 0 )
01031         metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
01032     else
01033         metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
01034 
01035     if( metrics_request_flag & V_TRI_EDGE_RATIO ) { metric_vals->edge_ratio = v_tri_edge_ratio( 3, coordinates ); }
01036     if( metrics_request_flag & V_TRI_RADIUS_RATIO )
01037     { metric_vals->radius_ratio = v_tri_radius_ratio( 3, coordinates ); }
01038     if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )  // there is no V_TRI_ASPECT_RATIO !
01039     { metric_vals->aspect_ratio = v_tri_radius_ratio( 3, coordinates ); }
01040 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines