MOAB: Mesh Oriented datABase  (version 5.4.0)
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 <cstddef>
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 )
00339     {
00340         min_angle = sides[2].interior_angle( sides[1] );
00341     }
00342     else if( short_side == 1 )
00343     {
00344         min_angle = sides[0].interior_angle( sides[2] );
00345     }
00346     else
00347     {
00348         min_angle = sides[0].interior_angle( sides[3] );
00349     }
00350 
00351     if( min_angle > 0 ) return (double)VERDICT_MIN( min_angle, VERDICT_DBL_MAX );
00352     return (double)VERDICT_MAX( min_angle, -VERDICT_DBL_MAX );
00353 }
00354 
00355 /*!
00356   The maximum angle of a tri
00357 
00358   The maximum angle of a tri is the maximum angle between
00359   two adjacents sides out of all three corners of the triangle.
00360 */
00361 C_FUNC_DEF double v_tri_maximum_angle( int /*num_nodes*/, double coordinates[][3] )
00362 {
00363 
00364     // vectors for all the sides
00365     VerdictVector sides[4];
00366     sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00367                   coordinates[1][2] - coordinates[0][2] );
00368     sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00369                   coordinates[2][2] - coordinates[1][2] );
00370     sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00371                   coordinates[2][2] - coordinates[0][2] );
00372 
00373     // in case we need to find the interior angle
00374     // between sides 0 and 1
00375     sides[3] = -sides[1];
00376 
00377     // calculate the lengths squared of the sides
00378     double sides_lengths[3];
00379     sides_lengths[0] = sides[0].length_squared();
00380     sides_lengths[1] = sides[1].length_squared();
00381     sides_lengths[2] = sides[2].length_squared();
00382 
00383     if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 )
00384     {
00385         return 0.0;
00386     }
00387 
00388     // using the law of sines, we know that the maximum
00389     // angle is opposite of the longest side
00390 
00391     // find the longest side
00392     int short_side = 0;
00393     if( sides_lengths[1] > sides_lengths[0] ) short_side = 1;
00394     if( sides_lengths[2] > sides_lengths[short_side] ) short_side = 2;
00395 
00396     // from the longest side, calculate the angle of the
00397     // opposite angle
00398     double max_angle;
00399     if( short_side == 0 )
00400     {
00401         max_angle = sides[2].interior_angle( sides[1] );
00402     }
00403     else if( short_side == 1 )
00404     {
00405         max_angle = sides[0].interior_angle( sides[2] );
00406     }
00407     else
00408     {
00409         max_angle = sides[0].interior_angle( sides[3] );
00410     }
00411 
00412     if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
00413     return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
00414 }
00415 
00416 /*!
00417   The condition of a tri
00418 
00419   Condition number of the jacobian matrix at any corner
00420 */
00421 C_FUNC_DEF double v_tri_condition( int /*num_nodes*/, double coordinates[][3] )
00422 {
00423     static const double rt3 = sqrt( 3.0 );
00424 
00425     VerdictVector v1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00426                       coordinates[1][2] - coordinates[0][2] );
00427 
00428     VerdictVector v2( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00429                       coordinates[2][2] - coordinates[0][2] );
00430 
00431     VerdictVector tri_normal = v1 * v2;
00432     double areax2            = tri_normal.length();
00433 
00434     if( areax2 == 0.0 ) return (double)VERDICT_DBL_MAX;
00435 
00436     double condition = (double)( ( ( v1 % v1 ) + ( v2 % v2 ) - ( v1 % v2 ) ) / ( areax2 * rt3 ) );
00437 
00438     // check for inverted if we have access to the normal
00439     if( compute_normal )
00440     {
00441         // center of tri
00442         double point[3], surf_normal[3];
00443         point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
00444         point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
00445         point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
00446 
00447         // dot product
00448         compute_normal( point, surf_normal );
00449         if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
00450             0 )
00451             return (double)VERDICT_DBL_MAX;
00452     }
00453     return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
00454 }
00455 
00456 /*!
00457   The scaled jacobian of a tri
00458 
00459   minimum of the jacobian divided by the lengths of 2 edge vectors
00460 */
00461 C_FUNC_DEF double v_tri_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
00462 {
00463     static const double detw = 2. / sqrt( 3.0 );
00464     VerdictVector first, second;
00465     double jacobian;
00466 
00467     VerdictVector edge[3];
00468     edge[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00469                  coordinates[1][2] - coordinates[0][2] );
00470 
00471     edge[1].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00472                  coordinates[2][2] - coordinates[0][2] );
00473 
00474     edge[2].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00475                  coordinates[2][2] - coordinates[1][2] );
00476     first  = edge[1] - edge[0];
00477     second = edge[2] - edge[0];
00478 
00479     VerdictVector cross = first * second;
00480     jacobian            = cross.length();
00481 
00482     double max_edge_length_product;
00483     max_edge_length_product =
00484         VERDICT_MAX( edge[0].length() * edge[1].length(),
00485                      VERDICT_MAX( edge[1].length() * edge[2].length(), edge[0].length() * edge[2].length() ) );
00486 
00487     if( max_edge_length_product < VERDICT_DBL_MIN ) return (double)0.0;
00488 
00489     jacobian *= detw;
00490     jacobian /= max_edge_length_product;
00491 
00492     if( compute_normal )
00493     {
00494         // center of tri
00495         double point[3], surf_normal[3];
00496         point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
00497         point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
00498         point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
00499 
00500         // dot product
00501         compute_normal( point, surf_normal );
00502         if( ( cross.x() * surf_normal[0] + cross.y() * surf_normal[1] + cross.z() * surf_normal[2] ) < 0 )
00503             jacobian *= -1;
00504     }
00505 
00506     if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
00507     return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
00508 }
00509 
00510 /*!
00511   The shape of a tri
00512 
00513   2 / condition number of weighted jacobian matrix
00514 */
00515 C_FUNC_DEF double v_tri_shape( int num_nodes, double coordinates[][3] )
00516 {
00517     double condition = v_tri_condition( num_nodes, coordinates );
00518 
00519     double shape;
00520     if( condition <= VERDICT_DBL_MIN )
00521         shape = VERDICT_DBL_MAX;
00522     else
00523         shape = ( 1 / condition );
00524 
00525     if( shape > 0 ) return (double)VERDICT_MIN( shape, VERDICT_DBL_MAX );
00526     return (double)VERDICT_MAX( shape, -VERDICT_DBL_MAX );
00527 }
00528 
00529 /*!
00530   The relative size of a tri
00531 
00532   Min(J,1/J) where J is the determinant of the weighted jacobian matrix.
00533 */
00534 C_FUNC_DEF double v_tri_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
00535 {
00536     double w11, w21, w12, w22;
00537 
00538     VerdictVector xxi, xet, tri_normal;
00539 
00540     v_tri_get_weight( w11, w21, w12, w22 );
00541 
00542     double detw = determinant( w11, w21, w12, w22 );
00543 
00544     if( detw == 0.0 ) return 0.0;
00545 
00546     xxi.set( coordinates[0][0] - coordinates[1][0], coordinates[0][1] - coordinates[1][1],
00547              coordinates[0][2] - coordinates[1][2] );
00548 
00549     xet.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
00550              coordinates[0][2] - coordinates[2][2] );
00551 
00552     tri_normal = xxi * xet;
00553 
00554     double deta = tri_normal.length();
00555     if( deta == 0.0 || detw == 0.0 ) return 0.0;
00556 
00557     double size = pow( deta / detw, 2 );
00558 
00559     double rel_size = VERDICT_MIN( size, 1.0 / size );
00560 
00561     if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
00562     return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
00563 }
00564 
00565 /*!
00566   The shape and size of a tri
00567 
00568   Product of the Shape and Relative Size
00569 */
00570 C_FUNC_DEF double v_tri_shape_and_size( int num_nodes, double coordinates[][3] )
00571 {
00572     double size, shape;
00573 
00574     size  = v_tri_relative_size_squared( num_nodes, coordinates );
00575     shape = v_tri_shape( num_nodes, coordinates );
00576 
00577     double shape_and_size = size * shape;
00578 
00579     if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
00580     return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
00581 }
00582 
00583 /*!
00584   The distortion of a tri
00585 
00586 TODO:  make a short definition of the distortion and comment below
00587 */
00588 C_FUNC_DEF double v_tri_distortion( int num_nodes, double coordinates[][3] )
00589 {
00590 
00591     double distortion;
00592     int total_number_of_gauss_points = 0;
00593     VerdictVector aa, bb, cc, normal_at_point, xin;
00594     double element_area = 0.;
00595 
00596     aa.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00597             coordinates[1][2] - coordinates[0][2] );
00598 
00599     bb.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00600             coordinates[2][2] - coordinates[0][2] );
00601 
00602     VerdictVector tri_normal = aa * bb;
00603 
00604     int number_of_gauss_points = 0;
00605     if( num_nodes == 3 )
00606     {
00607         distortion = 1.0;
00608         return (double)distortion;
00609     }
00610 
00611     else if( num_nodes == 6 )
00612     {
00613         total_number_of_gauss_points = 6;
00614         number_of_gauss_points       = 6;
00615     }
00616 
00617     distortion = VERDICT_DBL_MAX;
00618     double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
00619     double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
00620     double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
00621     double weight[maxTotalNumberGaussPoints];
00622 
00623     // create an object of GaussIntegration
00624     int number_dims = 2;
00625     int is_tri      = 1;
00626     GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
00627     GaussIntegration::calculate_shape_function_2d_tri();
00628     GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
00629 
00630     // calculate element area
00631     int ife, ja;
00632     for( ife = 0; ife < total_number_of_gauss_points; ife++ )
00633     {
00634         aa.set( 0.0, 0.0, 0.0 );
00635         bb.set( 0.0, 0.0, 0.0 );
00636 
00637         for( ja = 0; ja < num_nodes; ja++ )
00638         {
00639             xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
00640             aa += dndy1[ife][ja] * xin;
00641             bb += dndy2[ife][ja] * xin;
00642         }
00643         normal_at_point = aa * bb;
00644         double jacobian = normal_at_point.length();
00645         element_area += weight[ife] * jacobian;
00646     }
00647 
00648     element_area *= 0.8660254;
00649     double dndy1_at_node[maxNumberNodes][maxNumberNodes];
00650     double dndy2_at_node[maxNumberNodes][maxNumberNodes];
00651 
00652     GaussIntegration::calculate_derivative_at_nodes_2d_tri( dndy1_at_node, dndy2_at_node );
00653 
00654     VerdictVector normal_at_nodes[7];
00655 
00656     // evaluate normal at nodes and distortion values at nodes
00657     int jai = 0;
00658     for( ja = 0; ja < num_nodes; ja++ )
00659     {
00660         aa.set( 0.0, 0.0, 0.0 );
00661         bb.set( 0.0, 0.0, 0.0 );
00662         for( jai = 0; jai < num_nodes; jai++ )
00663         {
00664             xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
00665             aa += dndy1_at_node[ja][jai] * xin;
00666             bb += dndy2_at_node[ja][jai] * xin;
00667         }
00668         normal_at_nodes[ja] = aa * bb;
00669         normal_at_nodes[ja].normalize();
00670     }
00671 
00672     // determine if element is flat
00673     bool flat_element = true;
00674     double dot_product;
00675 
00676     for( ja = 0; ja < num_nodes; ja++ )
00677     {
00678         dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
00679         if( fabs( dot_product ) < 0.99 )
00680         {
00681             flat_element = false;
00682             break;
00683         }
00684     }
00685 
00686     // take into consideration of the thickness of the element
00687     double thickness, thickness_gauss;
00688     double distrt;
00689     // get_tri_thickness(tri, element_area, thickness );
00690     thickness = 0.001 * sqrt( element_area );
00691 
00692     // set thickness gauss point location
00693     double zl = 0.5773502691896;
00694     if( flat_element ) zl = 0.0;
00695 
00696     int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
00697     double thickness_z;
00698 
00699     // loop on integration points
00700     int igz;
00701     for( ife = 0; ife < total_number_of_gauss_points; ife++ )
00702     {
00703         // loop on the thickness direction gauss points
00704         for( igz = 0; igz < no_gauss_pts_z; igz++ )
00705         {
00706             zl          = -zl;
00707             thickness_z = zl * thickness / 2.0;
00708 
00709             aa.set( 0.0, 0.0, 0.0 );
00710             bb.set( 0.0, 0.0, 0.0 );
00711             cc.set( 0.0, 0.0, 0.0 );
00712 
00713             for( ja = 0; ja < num_nodes; ja++ )
00714             {
00715                 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
00716                 xin += thickness_z * normal_at_nodes[ja];
00717                 aa += dndy1[ife][ja] * xin;
00718                 bb += dndy2[ife][ja] * xin;
00719                 thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
00720                 cc += thickness_gauss * normal_at_nodes[ja];
00721             }
00722 
00723             normal_at_point = aa * bb;
00724             distrt          = cc % normal_at_point;
00725             if( distrt < distortion ) distortion = distrt;
00726         }
00727     }
00728 
00729     // loop through nodal points
00730     for( ja = 0; ja < num_nodes; ja++ )
00731     {
00732         for( igz = 0; igz < no_gauss_pts_z; igz++ )
00733         {
00734             zl          = -zl;
00735             thickness_z = zl * thickness / 2.0;
00736 
00737             aa.set( 0.0, 0.0, 0.0 );
00738             bb.set( 0.0, 0.0, 0.0 );
00739             cc.set( 0.0, 0.0, 0.0 );
00740 
00741             for( jai = 0; jai < num_nodes; jai++ )
00742             {
00743                 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
00744                 xin += thickness_z * normal_at_nodes[ja];
00745                 aa += dndy1_at_node[ja][jai] * xin;
00746                 bb += dndy2_at_node[ja][jai] * xin;
00747                 if( jai == ja )
00748                     thickness_gauss = thickness / 2.0;
00749                 else
00750                     thickness_gauss = 0.;
00751                 cc += thickness_gauss * normal_at_nodes[jai];
00752             }
00753         }
00754 
00755         normal_at_point      = aa * bb;
00756         double sign_jacobian = ( tri_normal % normal_at_point ) > 0 ? 1. : -1.;
00757         distrt               = sign_jacobian * ( cc % normal_at_point );
00758 
00759         if( distrt < distortion ) distortion = distrt;
00760     }
00761     if( element_area * thickness != 0 )
00762         distortion *= 1. / ( element_area * thickness );
00763     else
00764         distortion *= 1.;
00765 
00766     if( distortion > 0 ) return (double)VERDICT_MIN( distortion, VERDICT_DBL_MAX );
00767     return (double)VERDICT_MAX( distortion, -VERDICT_DBL_MAX );
00768 }
00769 
00770 /*!
00771   tri_quality for calculating multiple tri metrics at once
00772 
00773   using this method is generally faster than using the individual
00774   method multiple times.
00775 
00776 */
00777 C_FUNC_DEF void v_tri_quality( int num_nodes,
00778                                double coordinates[][3],
00779                                unsigned int metrics_request_flag,
00780                                TriMetricVals* metric_vals )
00781 {
00782 
00783     memset( metric_vals, 0, sizeof( TriMetricVals ) );
00784 
00785     // for starts, lets set up some basic and common information
00786 
00787     /*  node numbers and side numbers used below
00788 
00789                2
00790                ++
00791               /  \
00792            2 /    \ 1
00793             /      \
00794            /        \
00795          0 ---------+ 1
00796                0
00797     */
00798 
00799     // vectors for each side
00800     VerdictVector sides[3];
00801     sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00802                   coordinates[1][2] - coordinates[0][2] );
00803     sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00804                   coordinates[2][2] - coordinates[1][2] );
00805     sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00806                   coordinates[2][2] - coordinates[0][2] );
00807     VerdictVector tri_normal = sides[0] * sides[2];
00808     // if we have access to normal information, check to see if the
00809     // element is inverted.  If we don't have the normal information
00810     // that we need for this, assume the element is not inverted.
00811     // This flag will be used for condition number, jacobian, shape,
00812     // and size and shape.
00813     bool is_inverted = false;
00814     if( compute_normal )
00815     {
00816         // center of tri
00817         double point[3], surf_normal[3];
00818         point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
00819         point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
00820         point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
00821         // dot product
00822         compute_normal( point, surf_normal );
00823         if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
00824             0 )
00825             is_inverted = true;
00826     }
00827 
00828     // lengths squared of each side
00829     double sides_lengths_squared[3];
00830     sides_lengths_squared[0] = sides[0].length_squared();
00831     sides_lengths_squared[1] = sides[1].length_squared();
00832     sides_lengths_squared[2] = sides[2].length_squared();
00833 
00834     // if we are doing angle calcuations
00835     if( metrics_request_flag & ( V_TRI_MINIMUM_ANGLE | V_TRI_MAXIMUM_ANGLE ) )
00836     {
00837         // which is short and long side
00838         int short_side = 0, long_side = 0;
00839 
00840         if( sides_lengths_squared[1] < sides_lengths_squared[0] ) short_side = 1;
00841         if( sides_lengths_squared[2] < sides_lengths_squared[short_side] ) short_side = 2;
00842 
00843         if( sides_lengths_squared[1] > sides_lengths_squared[0] ) long_side = 1;
00844         if( sides_lengths_squared[2] > sides_lengths_squared[long_side] ) long_side = 2;
00845 
00846         // calculate the minimum angle of the tri
00847         if( metrics_request_flag & V_TRI_MINIMUM_ANGLE )
00848         {
00849             if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
00850             {
00851                 metric_vals->minimum_angle = 0.0;
00852             }
00853             else if( short_side == 0 )
00854                 metric_vals->minimum_angle = (double)sides[2].interior_angle( sides[1] );
00855             else if( short_side == 1 )
00856                 metric_vals->minimum_angle = (double)sides[0].interior_angle( sides[2] );
00857             else
00858                 metric_vals->minimum_angle = (double)sides[0].interior_angle( -sides[1] );
00859         }
00860 
00861         // calculate the maximum angle of the tri
00862         if( metrics_request_flag & V_TRI_MAXIMUM_ANGLE )
00863         {
00864             if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
00865             {
00866                 metric_vals->maximum_angle = 0.0;
00867             }
00868             else if( long_side == 0 )
00869                 metric_vals->maximum_angle = (double)sides[2].interior_angle( sides[1] );
00870             else if( long_side == 1 )
00871                 metric_vals->maximum_angle = (double)sides[0].interior_angle( sides[2] );
00872             else
00873                 metric_vals->maximum_angle = (double)sides[0].interior_angle( -sides[1] );
00874         }
00875     }
00876 
00877     // calculate the area of the tri
00878     // the following metrics depend on area
00879     if( metrics_request_flag &
00880         ( V_TRI_AREA | V_TRI_SCALED_JACOBIAN | V_TRI_SHAPE | V_TRI_RELATIVE_SIZE_SQUARED | V_TRI_SHAPE_AND_SIZE ) )
00881     {
00882         metric_vals->area = (double)( ( sides[0] * sides[2] ).length() * 0.5 );
00883     }
00884 
00885     // calculate the aspect ratio
00886     if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )
00887     {
00888         // sum the lengths squared
00889         double srms = sides_lengths_squared[0] + sides_lengths_squared[1] + sides_lengths_squared[2];
00890 
00891         // calculate once and reuse
00892         static const double twoTimesRootOf3 = 2 * sqrt( 3.0 );
00893 
00894         double div = ( twoTimesRootOf3 * ( ( sides[0] * sides[2] ).length() ) );
00895 
00896         if( div == 0.0 )
00897             metric_vals->aspect_frobenius = (double)VERDICT_DBL_MAX;
00898         else
00899             metric_vals->aspect_frobenius = (double)( srms / div );
00900     }
00901 
00902     // calculate the scaled jacobian
00903     if( metrics_request_flag & V_TRI_SCALED_JACOBIAN )
00904     {
00905         // calculate once and reuse
00906         static const double twoOverRootOf3 = 2 / sqrt( 3.0 );
00907         // use the area from above
00908 
00909         double tmp = tri_normal.length() * twoOverRootOf3;
00910 
00911         // now scale it by the lengths of the sides
00912         double min_scaled_jac = VERDICT_DBL_MAX;
00913         double temp_scaled_jac;
00914         for( int i = 0; i < 3; i++ )
00915         {
00916             if( sides_lengths_squared[i % 3] == 0.0 || sides_lengths_squared[( i + 2 ) % 3] == 0.0 )
00917                 temp_scaled_jac = 0.0;
00918             else
00919                 temp_scaled_jac =
00920                     tmp / sqrt( sides_lengths_squared[i % 3] ) / sqrt( sides_lengths_squared[( i + 2 ) % 3] );
00921             if( temp_scaled_jac < min_scaled_jac ) min_scaled_jac = temp_scaled_jac;
00922         }
00923         // multiply by -1 if the normals are in opposite directions
00924         if( is_inverted )
00925         {
00926             min_scaled_jac *= -1;
00927         }
00928         metric_vals->scaled_jacobian = (double)min_scaled_jac;
00929     }
00930 
00931     // calculate the condition number
00932     if( metrics_request_flag & V_TRI_CONDITION )
00933     {
00934         // calculate once and reuse
00935         static const double rootOf3 = sqrt( 3.0 );
00936         // if it is inverted, the condition number is considered to be infinity.
00937         if( is_inverted )
00938         {
00939             metric_vals->condition = VERDICT_DBL_MAX;
00940         }
00941         else
00942         {
00943             double area2x = ( sides[0] * sides[2] ).length();
00944             if( area2x == 0.0 )
00945                 metric_vals->condition = (double)( VERDICT_DBL_MAX );
00946             else
00947                 metric_vals->condition = (double)( ( sides[0] % sides[0] + sides[2] % sides[2] - sides[0] % sides[2] ) /
00948                                                    ( area2x * rootOf3 ) );
00949         }
00950     }
00951 
00952     // calculate the shape
00953     if( metrics_request_flag & V_TRI_SHAPE || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
00954     {
00955         // if element is inverted, shape is zero.  We don't need to
00956         // calculate anything.
00957         if( is_inverted )
00958         {
00959             metric_vals->shape = 0.0;
00960         }
00961         else
00962         {  // otherwise, we calculate the shape
00963             // calculate once and reuse
00964             static const double rootOf3 = sqrt( 3.0 );
00965             // reuse area from before
00966             double area2x = metric_vals->area * 2;
00967             // dot products
00968             double dots[3] = { sides[0] % sides[0], sides[2] % sides[2], sides[0] % sides[2] };
00969 
00970             // add the dots
00971             double sum_dots = dots[0] + dots[1] - dots[2];
00972 
00973             // then the finale
00974             if( sum_dots == 0.0 )
00975                 metric_vals->shape = 0.0;
00976             else
00977                 metric_vals->shape = (double)( rootOf3 * area2x / sum_dots );
00978         }
00979     }
00980 
00981     // calculate relative size squared
00982     if( metrics_request_flag & V_TRI_RELATIVE_SIZE_SQUARED || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
00983     {
00984         // get weights
00985         double w11, w21, w12, w22;
00986         v_tri_get_weight( w11, w21, w12, w22 );
00987         // get the determinant
00988         double detw = determinant( w11, w21, w12, w22 );
00989         // use the area from above and divide with the determinant
00990         if( metric_vals->area == 0.0 || detw == 0.0 )
00991             metric_vals->relative_size_squared = 0.0;
00992         else
00993         {
00994             double size = metric_vals->area * 2.0 / detw;
00995             // square the size
00996             size *= size;
00997             // value ranges between 0 to 1
00998             metric_vals->relative_size_squared = (double)VERDICT_MIN( size, 1.0 / size );
00999         }
01000     }
01001 
01002     // calculate shape and size
01003     if( metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
01004     {
01005         metric_vals->shape_and_size = metric_vals->relative_size_squared * metric_vals->shape;
01006     }
01007 
01008     // calculate distortion
01009     if( metrics_request_flag & V_TRI_DISTORTION ) metric_vals->distortion = v_tri_distortion( num_nodes, coordinates );
01010 
01011     // take care of any over-flow problems
01012     if( metric_vals->aspect_frobenius > 0 )
01013         metric_vals->aspect_frobenius = (double)VERDICT_MIN( metric_vals->aspect_frobenius, VERDICT_DBL_MAX );
01014     else
01015         metric_vals->aspect_frobenius = (double)VERDICT_MAX( metric_vals->aspect_frobenius, -VERDICT_DBL_MAX );
01016 
01017     if( metric_vals->area > 0 )
01018         metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
01019     else
01020         metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
01021 
01022     if( metric_vals->minimum_angle > 0 )
01023         metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
01024     else
01025         metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
01026 
01027     if( metric_vals->maximum_angle > 0 )
01028         metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
01029     else
01030         metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
01031 
01032     if( metric_vals->condition > 0 )
01033         metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
01034     else
01035         metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
01036 
01037     if( metric_vals->shape > 0 )
01038         metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
01039     else
01040         metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
01041 
01042     if( metric_vals->scaled_jacobian > 0 )
01043         metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
01044     else
01045         metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
01046 
01047     if( metric_vals->relative_size_squared > 0 )
01048         metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
01049     else
01050         metric_vals->relative_size_squared =
01051             (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
01052 
01053     if( metric_vals->shape_and_size > 0 )
01054         metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
01055     else
01056         metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
01057 
01058     if( metric_vals->distortion > 0 )
01059         metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
01060     else
01061         metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
01062 
01063     if( metrics_request_flag & V_TRI_EDGE_RATIO )
01064     {
01065         metric_vals->edge_ratio = v_tri_edge_ratio( 3, coordinates );
01066     }
01067     if( metrics_request_flag & V_TRI_RADIUS_RATIO )
01068     {
01069         metric_vals->radius_ratio = v_tri_radius_ratio( 3, coordinates );
01070     }
01071     if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )  // there is no V_TRI_ASPECT_RATIO !
01072     {
01073         metric_vals->aspect_ratio = v_tri_radius_ratio( 3, coordinates );
01074     }
01075 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines