MOAB: Mesh Oriented datABase  (version 5.2.1)
V_QuadMetric.cpp
Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Module:    $RCSfile: V_QuadMetric.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  * QuadMetric.cpp contains quality calculations for Quads
00018  *
00019  * This file is part of VERDICT
00020  *
00021  */
00022 
00023 #define VERDICT_EXPORTS
00024 
00025 #include "moab/verdict.h"
00026 #include "VerdictVector.hpp"
00027 #include "verdict_defines.hpp"
00028 #include "V_GaussIntegration.hpp"
00029 #include <memory.h>
00030 
00031 //! the average area of a quad
00032 double verdict_quad_size = 0;
00033 
00034 /*!
00035   weights based on the average size of a quad
00036 */
00037 int get_weight( double& m11, double& m21, double& m12, double& m22 )
00038 {
00039 
00040     m11 = 1;
00041     m21 = 0;
00042     m12 = 0;
00043     m22 = 1;
00044 
00045     double scale = sqrt( verdict_quad_size / ( m11 * m22 - m21 * m12 ) );
00046 
00047     m11 *= scale;
00048     m21 *= scale;
00049     m12 *= scale;
00050     m22 *= scale;
00051 
00052     return 1;
00053 }
00054 
00055 //! return the average area of a quad
00056 C_FUNC_DEF void v_set_quad_size( double size )
00057 {
00058     verdict_quad_size = size;
00059 }
00060 
00061 //! returns whether the quad is collapsed or not
00062 VerdictBoolean is_collapsed_quad( double coordinates[][3] )
00063 {
00064     if( coordinates[3][0] == coordinates[2][0] && coordinates[3][1] == coordinates[2][1] &&
00065         coordinates[3][2] == coordinates[2][2] )
00066         return VERDICT_TRUE;
00067 
00068     else
00069         return VERDICT_FALSE;
00070 }
00071 
00072 void make_quad_edges( VerdictVector edges[4], double coordinates[][3] )
00073 {
00074 
00075     edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00076                   coordinates[1][2] - coordinates[0][2] );
00077     edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00078                   coordinates[2][2] - coordinates[1][2] );
00079     edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00080                   coordinates[3][2] - coordinates[2][2] );
00081     edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
00082                   coordinates[0][2] - coordinates[3][2] );
00083 }
00084 
00085 void signed_corner_areas( double areas[4], double coordinates[][3] )
00086 {
00087     VerdictVector edges[4];
00088     make_quad_edges( edges, coordinates );
00089 
00090     VerdictVector corner_normals[4];
00091     corner_normals[0] = edges[3] * edges[0];
00092     corner_normals[1] = edges[0] * edges[1];
00093     corner_normals[2] = edges[1] * edges[2];
00094     corner_normals[3] = edges[2] * edges[3];
00095 
00096     // principal axes
00097     VerdictVector principal_axes[2];
00098     principal_axes[0] = edges[0] - edges[2];
00099     principal_axes[1] = edges[1] - edges[3];
00100 
00101     // quad center unit normal
00102     VerdictVector unit_center_normal;
00103     unit_center_normal = principal_axes[0] * principal_axes[1];
00104     unit_center_normal.normalize();
00105 
00106     areas[0] = unit_center_normal % corner_normals[0];
00107     areas[1] = unit_center_normal % corner_normals[1];
00108     areas[2] = unit_center_normal % corner_normals[2];
00109     areas[3] = unit_center_normal % corner_normals[3];
00110 }
00111 
00112 /*!
00113   localize the coordinates of a quad
00114 
00115   localizing puts the centriod of the quad
00116   at the orgin and also rotates the quad
00117   such that edge (0,1) is aligned with the x axis
00118   and the quad normal lines up with the y axis.
00119 
00120 */
00121 void localize_quad_coordinates( VerdictVector nodes[4] )
00122 {
00123     int i;
00124     VerdictVector global[4] = { nodes[0], nodes[1], nodes[2], nodes[3] };
00125 
00126     VerdictVector center = ( global[0] + global[1] + global[2] + global[3] ) / 4.0;
00127     for( i = 0; i < 4; i++ )
00128         global[i] -= center;
00129 
00130     VerdictVector vector_diff;
00131     VerdictVector vector_sum;
00132     VerdictVector ref_point( 0.0, 0.0, 0.0 );
00133     VerdictVector tmp_vector, normal( 0.0, 0.0, 0.0 );
00134     VerdictVector vector1, vector2;
00135 
00136     for( i = 0; i < 4; i++ )
00137     {
00138         vector1     = global[i];
00139         vector2     = global[( i + 1 ) % 4];
00140         vector_diff = vector2 - vector1;
00141         ref_point += vector1;
00142         vector_sum = vector1 + vector2;
00143 
00144         tmp_vector.set( vector_diff.y() * vector_sum.z(), vector_diff.z() * vector_sum.x(),
00145                         vector_diff.x() * vector_sum.y() );
00146         normal += tmp_vector;
00147     }
00148 
00149     normal.normalize();
00150     normal *= -1.0;
00151 
00152     VerdictVector local_x_axis = global[1] - global[0];
00153     local_x_axis.normalize();
00154 
00155     VerdictVector local_y_axis = normal * local_x_axis;
00156     local_y_axis.normalize();
00157 
00158     for( i = 0; i < 4; i++ )
00159     {
00160         nodes[i].x( global[i] % local_x_axis );
00161         nodes[i].y( global[i] % local_y_axis );
00162         nodes[i].z( global[i] % normal );
00163     }
00164 }
00165 
00166 /*!
00167   moves and rotates the quad such that it enables us to
00168   use components of ef's
00169 */
00170 void localize_quad_for_ef( VerdictVector node_pos[4] )
00171 {
00172 
00173     VerdictVector centroid( node_pos[0] );
00174     centroid += node_pos[1];
00175     centroid += node_pos[2];
00176     centroid += node_pos[3];
00177 
00178     centroid /= 4.0;
00179 
00180     node_pos[0] -= centroid;
00181     node_pos[1] -= centroid;
00182     node_pos[2] -= centroid;
00183     node_pos[3] -= centroid;
00184 
00185     VerdictVector rotate = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
00186     rotate.normalize();
00187 
00188     double cosine = rotate.x();
00189     double sine   = rotate.y();
00190 
00191     double xnew;
00192 
00193     for( int i = 0; i < 4; i++ )
00194     {
00195         xnew = cosine * node_pos[i].x() + sine * node_pos[i].y();
00196         node_pos[i].y( -sine * node_pos[i].x() + cosine * node_pos[i].y() );
00197         node_pos[i].x( xnew );
00198     }
00199 }
00200 
00201 /*!
00202   returns the normal vector of a quad
00203 */
00204 VerdictVector quad_normal( double coordinates[][3] )
00205 {
00206     // get normal at node 0
00207     VerdictVector edge0, edge1;
00208 
00209     edge0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00210                coordinates[1][2] - coordinates[0][2] );
00211 
00212     edge1.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
00213                coordinates[3][2] - coordinates[0][2] );
00214 
00215     VerdictVector norm0 = edge0 * edge1;
00216     norm0.normalize();
00217 
00218     // because some faces may have obtuse angles, check another normal at
00219     // node 2 for consistent sense
00220 
00221     edge0.set( coordinates[2][0] - coordinates[3][0], coordinates[2][1] - coordinates[3][1],
00222                coordinates[2][2] - coordinates[3][2] );
00223 
00224     edge1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00225                coordinates[2][2] - coordinates[1][2] );
00226 
00227     VerdictVector norm2 = edge0 * edge1;
00228     norm2.normalize();
00229 
00230     // if these two agree, we are done, else test a third to decide
00231 
00232     if( ( norm0 % norm2 ) > 0.0 )
00233     {
00234         norm0 += norm2;
00235         norm0 *= 0.5;
00236         return norm0;
00237     }
00238 
00239     // test normal at node1
00240 
00241     edge0.set( coordinates[1][0] - coordinates[2][0], coordinates[1][1] - coordinates[2][1],
00242                coordinates[1][2] - coordinates[2][2] );
00243 
00244     edge1.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00245                coordinates[1][2] - coordinates[0][2] );
00246 
00247     VerdictVector norm1 = edge0 * edge1;
00248     norm1.normalize();
00249 
00250     if( ( norm0 % norm1 ) > 0.0 )
00251     {
00252         norm0 += norm1;
00253         norm0 *= 0.5;
00254         return norm0;
00255     }
00256     else
00257     {
00258         norm2 += norm1;
00259         norm2 *= 0.5;
00260         return norm2;
00261     }
00262 }
00263 
00264 /*!
00265    the edge ratio of a quad
00266 
00267    NB (P. Pebay 01/19/07):
00268      Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
00269      minimum edge lengths
00270 */
00271 C_FUNC_DEF double v_quad_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
00272 {
00273     VerdictVector edges[4];
00274     make_quad_edges( edges, coordinates );
00275 
00276     double a2 = edges[0].length_squared();
00277     double b2 = edges[1].length_squared();
00278     double c2 = edges[2].length_squared();
00279     double d2 = edges[3].length_squared();
00280 
00281     double mab, Mab, mcd, Mcd, m2, M2;
00282     if( a2 < b2 )
00283     {
00284         mab = a2;
00285         Mab = b2;
00286     }
00287     else  // b2 <= a2
00288     {
00289         mab = b2;
00290         Mab = a2;
00291     }
00292     if( c2 < d2 )
00293     {
00294         mcd = c2;
00295         Mcd = d2;
00296     }
00297     else  // d2 <= c2
00298     {
00299         mcd = d2;
00300         Mcd = c2;
00301     }
00302     m2 = mab < mcd ? mab : mcd;
00303     M2 = Mab > Mcd ? Mab : Mcd;
00304 
00305     if( m2 < VERDICT_DBL_MIN )
00306         return (double)VERDICT_DBL_MAX;
00307     else
00308     {
00309         double edge_ratio = sqrt( M2 / m2 );
00310 
00311         if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
00312         return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
00313     }
00314 }
00315 
00316 /*!
00317   maximum of edge ratio of a quad
00318 
00319   maximum edge length ratio at quad center
00320 */
00321 C_FUNC_DEF double v_quad_max_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
00322 {
00323     VerdictVector quad_nodes[4];
00324     quad_nodes[0].set( coordinates[0][0], coordinates[0][1], coordinates[0][2] );
00325     quad_nodes[1].set( coordinates[1][0], coordinates[1][1], coordinates[1][2] );
00326     quad_nodes[2].set( coordinates[2][0], coordinates[2][1], coordinates[2][2] );
00327     quad_nodes[3].set( coordinates[3][0], coordinates[3][1], coordinates[3][2] );
00328 
00329     VerdictVector principal_axes[2];
00330     principal_axes[0] = quad_nodes[1] + quad_nodes[2] - quad_nodes[0] - quad_nodes[3];
00331     principal_axes[1] = quad_nodes[2] + quad_nodes[3] - quad_nodes[0] - quad_nodes[1];
00332 
00333     double len1 = principal_axes[0].length();
00334     double len2 = principal_axes[1].length();
00335 
00336     if( len1 < VERDICT_DBL_MIN || len2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00337 
00338     double max_edge_ratio = VERDICT_MAX( len1 / len2, len2 / len1 );
00339 
00340     if( max_edge_ratio > 0 ) return (double)VERDICT_MIN( max_edge_ratio, VERDICT_DBL_MAX );
00341     return (double)VERDICT_MAX( max_edge_ratio, -VERDICT_DBL_MAX );
00342 }
00343 
00344 /*!
00345    the aspect ratio of a quad
00346 
00347    NB (P. Pebay 01/20/07):
00348      this is a generalization of the triangle aspect ratio
00349      using Heron's formula.
00350 */
00351 C_FUNC_DEF double v_quad_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
00352 {
00353 
00354     VerdictVector edges[4];
00355     make_quad_edges( edges, coordinates );
00356 
00357     double a1 = edges[0].length();
00358     double b1 = edges[1].length();
00359     double c1 = edges[2].length();
00360     double d1 = edges[3].length();
00361 
00362     double ma = a1 > b1 ? a1 : b1;
00363     double mb = c1 > d1 ? c1 : d1;
00364     double hm = ma > mb ? ma : mb;
00365 
00366     VerdictVector ab   = edges[0] * edges[1];
00367     VerdictVector cd   = edges[2] * edges[3];
00368     double denominator = ab.length() + cd.length();
00369 
00370     if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00371 
00372     double aspect_ratio = .5 * hm * ( a1 + b1 + c1 + d1 ) / denominator;
00373 
00374     if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
00375     return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
00376 }
00377 
00378 /*!
00379    the radius ratio of a quad
00380 
00381    NB (P. Pebay 01/19/07):
00382      this metric is called "radius ratio" by extension of a concept that does
00383      not exist in general with quads -- although a different name should probably
00384      be used in the future.
00385 */
00386 C_FUNC_DEF double v_quad_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
00387 {
00388     static const double normal_coeff = 1. / ( 2. * sqrt( 2. ) );
00389 
00390     VerdictVector edges[4];
00391     make_quad_edges( edges, coordinates );
00392 
00393     double a2 = edges[0].length_squared();
00394     double b2 = edges[1].length_squared();
00395     double c2 = edges[2].length_squared();
00396     double d2 = edges[3].length_squared();
00397 
00398     VerdictVector diag;
00399     diag.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00400               coordinates[2][2] - coordinates[0][2] );
00401     double m2 = diag.length_squared();
00402 
00403     diag.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00404               coordinates[3][2] - coordinates[1][2] );
00405     double n2 = diag.length_squared();
00406 
00407     double t0 = a2 > b2 ? a2 : b2;
00408     double t1 = c2 > d2 ? c2 : d2;
00409     double t2 = m2 > n2 ? m2 : n2;
00410     double h2 = t0 > t1 ? t0 : t1;
00411     h2        = h2 > t2 ? h2 : t2;
00412 
00413     VerdictVector ab = edges[0] * edges[1];
00414     VerdictVector bc = edges[1] * edges[2];
00415     VerdictVector cd = edges[2] * edges[3];
00416     VerdictVector da = edges[3] * edges[0];
00417 
00418     t0        = da.length();
00419     t1        = ab.length();
00420     t2        = bc.length();
00421     double t3 = cd.length();
00422 
00423     t0 = t0 < t1 ? t0 : t1;
00424     t2 = t2 < t3 ? t2 : t3;
00425     t0 = t0 < t2 ? t0 : t2;
00426 
00427     if( t0 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00428 
00429     double radius_ratio = normal_coeff * sqrt( ( a2 + b2 + c2 + d2 ) * h2 ) / t0;
00430 
00431     if( radius_ratio > 0 ) return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
00432     return (double)VERDICT_MAX( radius_ratio, -VERDICT_DBL_MAX );
00433 }
00434 
00435 /*!
00436    the average Frobenius aspect of a quad
00437 
00438    NB (P. Pebay 01/20/07):
00439      this metric is calculated by averaging the 4 Frobenius aspects at
00440      each corner of the quad, when the reference triangle is right isosceles.
00441 */
00442 C_FUNC_DEF double v_quad_med_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
00443 {
00444 
00445     VerdictVector edges[4];
00446     make_quad_edges( edges, coordinates );
00447 
00448     double a2 = edges[0].length_squared();
00449     double b2 = edges[1].length_squared();
00450     double c2 = edges[2].length_squared();
00451     double d2 = edges[3].length_squared();
00452 
00453     VerdictVector ab = edges[0] * edges[1];
00454     VerdictVector bc = edges[1] * edges[2];
00455     VerdictVector cd = edges[2] * edges[3];
00456     VerdictVector da = edges[3] * edges[0];
00457 
00458     double ab1 = ab.length();
00459     double bc1 = bc.length();
00460     double cd1 = cd.length();
00461     double da1 = da.length();
00462 
00463     if( ab1 < VERDICT_DBL_MIN || bc1 < VERDICT_DBL_MIN || cd1 < VERDICT_DBL_MIN || da1 < VERDICT_DBL_MIN )
00464         return (double)VERDICT_DBL_MAX;
00465 
00466     double qsum = ( a2 + b2 ) / ab1;
00467     qsum += ( b2 + c2 ) / bc1;
00468     qsum += ( c2 + d2 ) / cd1;
00469     qsum += ( d2 + a2 ) / da1;
00470 
00471     double med_aspect_frobenius = .125 * qsum;
00472 
00473     if( med_aspect_frobenius > 0 ) return (double)VERDICT_MIN( med_aspect_frobenius, VERDICT_DBL_MAX );
00474     return (double)VERDICT_MAX( med_aspect_frobenius, -VERDICT_DBL_MAX );
00475 }
00476 
00477 /*!
00478    the maximum Frobenius aspect of a quad
00479 
00480    NB (P. Pebay 01/20/07):
00481      this metric is calculated by taking the maximum of the 4 Frobenius aspects at
00482      each corner of the quad, when the reference triangle is right isosceles.
00483 */
00484 C_FUNC_DEF double v_quad_max_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
00485 {
00486 
00487     VerdictVector edges[4];
00488     make_quad_edges( edges, coordinates );
00489 
00490     double a2 = edges[0].length_squared();
00491     double b2 = edges[1].length_squared();
00492     double c2 = edges[2].length_squared();
00493     double d2 = edges[3].length_squared();
00494 
00495     VerdictVector ab = edges[0] * edges[1];
00496     VerdictVector bc = edges[1] * edges[2];
00497     VerdictVector cd = edges[2] * edges[3];
00498     VerdictVector da = edges[3] * edges[0];
00499 
00500     double ab1 = ab.length();
00501     double bc1 = bc.length();
00502     double cd1 = cd.length();
00503     double da1 = da.length();
00504 
00505     if( ab1 < VERDICT_DBL_MIN || bc1 < VERDICT_DBL_MIN || cd1 < VERDICT_DBL_MIN || da1 < VERDICT_DBL_MIN )
00506         return (double)VERDICT_DBL_MAX;
00507 
00508     double qmax = ( a2 + b2 ) / ab1;
00509 
00510     double qcur = ( b2 + c2 ) / bc1;
00511     qmax        = qmax > qcur ? qmax : qcur;
00512 
00513     qcur = ( c2 + d2 ) / cd1;
00514     qmax = qmax > qcur ? qmax : qcur;
00515 
00516     qcur = ( d2 + a2 ) / da1;
00517     qmax = qmax > qcur ? qmax : qcur;
00518 
00519     double max_aspect_frobenius = .5 * qmax;
00520 
00521     if( max_aspect_frobenius > 0 ) return (double)VERDICT_MIN( max_aspect_frobenius, VERDICT_DBL_MAX );
00522     return (double)VERDICT_MAX( max_aspect_frobenius, -VERDICT_DBL_MAX );
00523 }
00524 
00525 /*!
00526   skew of a quad
00527 
00528   maximum ||cos A|| where A is the angle between edges at quad center
00529 */
00530 C_FUNC_DEF double v_quad_skew( int /*num_nodes*/, double coordinates[][3] )
00531 {
00532     VerdictVector node_pos[4];
00533     for( int i = 0; i < 4; i++ )
00534         node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
00535 
00536     VerdictVector principle_axes[2];
00537     principle_axes[0] = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
00538     principle_axes[1] = node_pos[2] + node_pos[3] - node_pos[0] - node_pos[1];
00539 
00540     if( principle_axes[0].normalize() < VERDICT_DBL_MIN ) return 0.0;
00541     if( principle_axes[1].normalize() < VERDICT_DBL_MIN ) return 0.0;
00542 
00543     double skew = fabs( principle_axes[0] % principle_axes[1] );
00544 
00545     return (double)VERDICT_MIN( skew, VERDICT_DBL_MAX );
00546 }
00547 
00548 /*!
00549   taper of a quad
00550 
00551   maximum ratio of lengths derived from opposite edges
00552 */
00553 C_FUNC_DEF double v_quad_taper( int /*num_nodes*/, double coordinates[][3] )
00554 {
00555     VerdictVector node_pos[4];
00556     for( int i = 0; i < 4; i++ )
00557         node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
00558 
00559     VerdictVector principle_axes[2];
00560     principle_axes[0] = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
00561     principle_axes[1] = node_pos[2] + node_pos[3] - node_pos[0] - node_pos[1];
00562 
00563     VerdictVector cross_derivative = node_pos[0] + node_pos[2] - node_pos[1] - node_pos[3];
00564 
00565     double lengths[2];
00566     lengths[0] = principle_axes[0].length();
00567     lengths[1] = principle_axes[1].length();
00568 
00569     // get min length
00570     lengths[0] = VERDICT_MIN( lengths[0], lengths[1] );
00571 
00572     if( lengths[0] < VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
00573 
00574     double taper = cross_derivative.length() / lengths[0];
00575     return (double)VERDICT_MIN( taper, VERDICT_DBL_MAX );
00576 }
00577 
00578 /*!
00579   warpage of a quad
00580 
00581   deviation of element from planarity
00582 */
00583 C_FUNC_DEF double v_quad_warpage( int /*num_nodes*/, double coordinates[][3] )
00584 {
00585 
00586     VerdictVector edges[4];
00587     make_quad_edges( edges, coordinates );
00588 
00589     VerdictVector corner_normals[4];
00590     corner_normals[0] = edges[3] * edges[0];
00591     corner_normals[1] = edges[0] * edges[1];
00592     corner_normals[2] = edges[1] * edges[2];
00593     corner_normals[3] = edges[2] * edges[3];
00594 
00595     if( corner_normals[0].normalize() < VERDICT_DBL_MIN || corner_normals[1].normalize() < VERDICT_DBL_MIN ||
00596         corner_normals[2].normalize() < VERDICT_DBL_MIN || corner_normals[3].normalize() < VERDICT_DBL_MIN )
00597         return (double)VERDICT_DBL_MIN;
00598 
00599     double warpage =
00600         pow( VERDICT_MIN( corner_normals[0] % corner_normals[2], corner_normals[1] % corner_normals[3] ), 3 );
00601 
00602     if( warpage > 0 ) return (double)VERDICT_MIN( warpage, VERDICT_DBL_MAX );
00603     return (double)VERDICT_MAX( warpage, -VERDICT_DBL_MAX );
00604 }
00605 
00606 /*!
00607   the area of a quad
00608 
00609   jacobian at quad center
00610 */
00611 C_FUNC_DEF double v_quad_area( int /*num_nodes*/, double coordinates[][3] )
00612 {
00613 
00614     double corner_areas[4];
00615     signed_corner_areas( corner_areas, coordinates );
00616 
00617     double area = 0.25 * ( corner_areas[0] + corner_areas[1] + corner_areas[2] + corner_areas[3] );
00618 
00619     if( area > 0 ) return (double)VERDICT_MIN( area, VERDICT_DBL_MAX );
00620     return (double)VERDICT_MAX( area, -VERDICT_DBL_MAX );
00621 }
00622 
00623 /*!
00624   the stretch of a quad
00625 
00626   sqrt(2) * minimum edge length / maximum diagonal length
00627 */
00628 C_FUNC_DEF double v_quad_stretch( int /*num_nodes*/, double coordinates[][3] )
00629 {
00630     VerdictVector edges[4], temp;
00631     make_quad_edges( edges, coordinates );
00632 
00633     double lengths_squared[4];
00634     lengths_squared[0] = edges[0].length_squared();
00635     lengths_squared[1] = edges[1].length_squared();
00636     lengths_squared[2] = edges[2].length_squared();
00637     lengths_squared[3] = edges[3].length_squared();
00638 
00639     temp.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
00640               coordinates[2][2] - coordinates[0][2] );
00641     double diag02 = temp.length_squared();
00642 
00643     temp.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
00644               coordinates[3][2] - coordinates[1][2] );
00645     double diag13 = temp.length_squared();
00646 
00647     static const double QUAD_STRETCH_FACTOR = sqrt( 2.0 );
00648 
00649     // 'diag02' is now the max diagonal of the quad
00650     diag02 = VERDICT_MAX( diag02, diag13 );
00651 
00652     if( diag02 < VERDICT_DBL_MIN )
00653         return (double)VERDICT_DBL_MAX;
00654     else
00655     {
00656         double stretch =
00657             (double)( QUAD_STRETCH_FACTOR * sqrt( VERDICT_MIN( VERDICT_MIN( lengths_squared[0], lengths_squared[1] ),
00658                                                                VERDICT_MIN( lengths_squared[2], lengths_squared[3] ) ) /
00659                                                   diag02 ) );
00660 
00661         return (double)VERDICT_MIN( stretch, VERDICT_DBL_MAX );
00662     }
00663 }
00664 
00665 /*!
00666   the largest angle of a quad
00667 
00668   largest included quad area (degrees)
00669 */
00670 C_FUNC_DEF double v_quad_maximum_angle( int /*num_nodes*/, double coordinates[][3] )
00671 {
00672 
00673     // if this is a collapsed quad, just pass it on to
00674     // the tri_largest_angle routine
00675     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_maximum_angle( 3, coordinates );
00676 
00677     double angle;
00678     double max_angle = 0.0;
00679 
00680     VerdictVector edges[4];
00681     edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00682                   coordinates[1][2] - coordinates[0][2] );
00683     edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00684                   coordinates[2][2] - coordinates[1][2] );
00685     edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00686                   coordinates[3][2] - coordinates[2][2] );
00687     edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
00688                   coordinates[0][2] - coordinates[3][2] );
00689 
00690     // go around each node and calculate the angle
00691     // at each node
00692     double length[4];
00693     length[0] = edges[0].length();
00694     length[1] = edges[1].length();
00695     length[2] = edges[2].length();
00696     length[3] = edges[3].length();
00697 
00698     if( length[0] <= VERDICT_DBL_MIN || length[1] <= VERDICT_DBL_MIN || length[2] <= VERDICT_DBL_MIN ||
00699         length[3] <= VERDICT_DBL_MIN )
00700         return 0.0;
00701 
00702     angle     = acos( -( edges[0] % edges[1] ) / ( length[0] * length[1] ) );
00703     max_angle = VERDICT_MAX( angle, max_angle );
00704 
00705     angle     = acos( -( edges[1] % edges[2] ) / ( length[1] * length[2] ) );
00706     max_angle = VERDICT_MAX( angle, max_angle );
00707 
00708     angle     = acos( -( edges[2] % edges[3] ) / ( length[2] * length[3] ) );
00709     max_angle = VERDICT_MAX( angle, max_angle );
00710 
00711     angle     = acos( -( edges[3] % edges[0] ) / ( length[3] * length[0] ) );
00712     max_angle = VERDICT_MAX( angle, max_angle );
00713 
00714     max_angle = max_angle * 180.0 / VERDICT_PI;
00715 
00716     // if any signed areas are < 0, then you are getting the wrong angle
00717     double areas[4];
00718     signed_corner_areas( areas, coordinates );
00719 
00720     if( areas[0] < 0 || areas[1] < 0 || areas[2] < 0 || areas[3] < 0 ) { max_angle = 360 - max_angle; }
00721 
00722     if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
00723     return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
00724 }
00725 
00726 /*!
00727   the smallest angle of a quad
00728 
00729   smallest included quad angle (degrees)
00730 */
00731 C_FUNC_DEF double v_quad_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
00732 {
00733     // if this quad is a collapsed quad, then just
00734     // send it to the tri_smallest_angle routine
00735     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_minimum_angle( 3, coordinates );
00736 
00737     double angle;
00738     double min_angle = 360.0;
00739 
00740     VerdictVector edges[4];
00741     edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00742                   coordinates[1][2] - coordinates[0][2] );
00743     edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00744                   coordinates[2][2] - coordinates[1][2] );
00745     edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00746                   coordinates[3][2] - coordinates[2][2] );
00747     edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
00748                   coordinates[0][2] - coordinates[3][2] );
00749 
00750     // go around each node and calculate the angle
00751     // at each node
00752     double length[4];
00753     length[0] = edges[0].length();
00754     length[1] = edges[1].length();
00755     length[2] = edges[2].length();
00756     length[3] = edges[3].length();
00757 
00758     if( length[0] <= VERDICT_DBL_MIN || length[1] <= VERDICT_DBL_MIN || length[2] <= VERDICT_DBL_MIN ||
00759         length[3] <= VERDICT_DBL_MIN )
00760         return 360.0;
00761 
00762     angle     = acos( -( edges[0] % edges[1] ) / ( length[0] * length[1] ) );
00763     min_angle = VERDICT_MIN( angle, min_angle );
00764 
00765     angle     = acos( -( edges[1] % edges[2] ) / ( length[1] * length[2] ) );
00766     min_angle = VERDICT_MIN( angle, min_angle );
00767 
00768     angle     = acos( -( edges[2] % edges[3] ) / ( length[2] * length[3] ) );
00769     min_angle = VERDICT_MIN( angle, min_angle );
00770 
00771     angle     = acos( -( edges[3] % edges[0] ) / ( length[3] * length[0] ) );
00772     min_angle = VERDICT_MIN( angle, min_angle );
00773 
00774     min_angle = min_angle * 180.0 / VERDICT_PI;
00775 
00776     if( min_angle > 0 ) return (double)VERDICT_MIN( min_angle, VERDICT_DBL_MAX );
00777     return (double)VERDICT_MAX( min_angle, -VERDICT_DBL_MAX );
00778 }
00779 
00780 /*!
00781   the oddy of a quad
00782 
00783   general distortion measure based on left Cauchy-Green Tensor
00784 */
00785 C_FUNC_DEF double v_quad_oddy( int /*num_nodes*/, double coordinates[][3] )
00786 {
00787 
00788     double max_oddy = 0.;
00789 
00790     VerdictVector first, second, node_pos[4];
00791 
00792     double g, g11, g12, g22, cur_oddy;
00793     int i;
00794 
00795     for( i = 0; i < 4; i++ )
00796         node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
00797 
00798     for( i = 0; i < 4; i++ )
00799     {
00800         first  = node_pos[i] - node_pos[( i + 1 ) % 4];
00801         second = node_pos[i] - node_pos[( i + 3 ) % 4];
00802 
00803         g11 = first % first;
00804         g12 = first % second;
00805         g22 = second % second;
00806         g   = g11 * g22 - g12 * g12;
00807 
00808         if( g < VERDICT_DBL_MIN )
00809             cur_oddy = VERDICT_DBL_MAX;
00810         else
00811             cur_oddy = ( ( g11 - g22 ) * ( g11 - g22 ) + 4. * g12 * g12 ) / 2. / g;
00812 
00813         max_oddy = VERDICT_MAX( max_oddy, cur_oddy );
00814     }
00815 
00816     if( max_oddy > 0 ) return (double)VERDICT_MIN( max_oddy, VERDICT_DBL_MAX );
00817     return (double)VERDICT_MAX( max_oddy, -VERDICT_DBL_MAX );
00818 }
00819 
00820 /*!
00821   the condition of a quad
00822 
00823   maximum condition number of the Jacobian matrix at 4 corners
00824 */
00825 C_FUNC_DEF double v_quad_condition( int /*num_nodes*/, double coordinates[][3] )
00826 {
00827 
00828     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_condition( 3, coordinates );
00829 
00830     double areas[4];
00831     signed_corner_areas( areas, coordinates );
00832 
00833     double max_condition = 0.;
00834 
00835     VerdictVector xxi, xet;
00836 
00837     double condition;
00838 
00839     for( int i = 0; i < 4; i++ )
00840     {
00841 
00842         xxi.set( coordinates[i][0] - coordinates[( i + 1 ) % 4][0], coordinates[i][1] - coordinates[( i + 1 ) % 4][1],
00843                  coordinates[i][2] - coordinates[( i + 1 ) % 4][2] );
00844 
00845         xet.set( coordinates[i][0] - coordinates[( i + 3 ) % 4][0], coordinates[i][1] - coordinates[( i + 3 ) % 4][1],
00846                  coordinates[i][2] - coordinates[( i + 3 ) % 4][2] );
00847 
00848         if( areas[i] < VERDICT_DBL_MIN )
00849             condition = VERDICT_DBL_MAX;
00850         else
00851             condition = ( xxi % xxi + xet % xet ) / areas[i];
00852 
00853         max_condition = VERDICT_MAX( max_condition, condition );
00854     }
00855 
00856     max_condition /= 2;
00857 
00858     if( max_condition > 0 ) return (double)VERDICT_MIN( max_condition, VERDICT_DBL_MAX );
00859     return (double)VERDICT_MAX( max_condition, -VERDICT_DBL_MAX );
00860 }
00861 
00862 /*!
00863   the jacobian of a quad
00864 
00865   minimum pointwise volume of local map at 4 corners and center of quad
00866 */
00867 C_FUNC_DEF double v_quad_jacobian( int /*num_nodes*/, double coordinates[][3] )
00868 {
00869 
00870     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return (double)( v_tri_area( 3, coordinates ) * 2.0 );
00871 
00872     double areas[4];
00873     signed_corner_areas( areas, coordinates );
00874 
00875     double jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) );
00876     if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
00877     return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
00878 }
00879 
00880 /*!
00881   scaled jacobian of a quad
00882 
00883   Minimum Jacobian divided by the lengths of the 2 edge vector
00884 */
00885 C_FUNC_DEF double v_quad_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
00886 {
00887     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_scaled_jacobian( 3, coordinates );
00888 
00889     double corner_areas[4], min_scaled_jac = VERDICT_DBL_MAX, scaled_jac;
00890     signed_corner_areas( corner_areas, coordinates );
00891 
00892     VerdictVector edges[4];
00893     make_quad_edges( edges, coordinates );
00894 
00895     double length[4];
00896     length[0] = edges[0].length();
00897     length[1] = edges[1].length();
00898     length[2] = edges[2].length();
00899     length[3] = edges[3].length();
00900 
00901     if( length[0] < VERDICT_DBL_MIN || length[1] < VERDICT_DBL_MIN || length[2] < VERDICT_DBL_MIN ||
00902         length[3] < VERDICT_DBL_MIN )
00903         return 0.0;
00904 
00905     scaled_jac     = corner_areas[0] / ( length[0] * length[3] );
00906     min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
00907 
00908     scaled_jac     = corner_areas[1] / ( length[1] * length[0] );
00909     min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
00910 
00911     scaled_jac     = corner_areas[2] / ( length[2] * length[1] );
00912     min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
00913 
00914     scaled_jac     = corner_areas[3] / ( length[3] * length[2] );
00915     min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
00916 
00917     if( min_scaled_jac > 0 ) return (double)VERDICT_MIN( min_scaled_jac, VERDICT_DBL_MAX );
00918     return (double)VERDICT_MAX( min_scaled_jac, -VERDICT_DBL_MAX );
00919 }
00920 
00921 /*!
00922   the shear of a quad
00923 
00924   2/Condition number of Jacobian Skew matrix
00925 */
00926 C_FUNC_DEF double v_quad_shear( int /*num_nodes*/, double coordinates[][3] )
00927 {
00928     double scaled_jacobian = v_quad_scaled_jacobian( 4, coordinates );
00929 
00930     if( scaled_jacobian <= VERDICT_DBL_MIN )
00931         return 0.0;
00932     else
00933         return (double)VERDICT_MIN( scaled_jacobian, VERDICT_DBL_MAX );
00934 }
00935 
00936 /*!
00937   the shape of a quad
00938 
00939    2/Condition number of weighted Jacobian matrix
00940 */
00941 C_FUNC_DEF double v_quad_shape( int /*num_nodes*/, double coordinates[][3] )
00942 {
00943 
00944     double corner_areas[4], min_shape = VERDICT_DBL_MAX, shape;
00945     signed_corner_areas( corner_areas, coordinates );
00946 
00947     VerdictVector edges[4];
00948     make_quad_edges( edges, coordinates );
00949 
00950     double length_squared[4];
00951     length_squared[0] = edges[0].length_squared();
00952     length_squared[1] = edges[1].length_squared();
00953     length_squared[2] = edges[2].length_squared();
00954     length_squared[3] = edges[3].length_squared();
00955 
00956     if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN ||
00957         length_squared[2] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN )
00958         return 0.0;
00959 
00960     shape     = corner_areas[0] / ( length_squared[0] + length_squared[3] );
00961     min_shape = VERDICT_MIN( shape, min_shape );
00962 
00963     shape     = corner_areas[1] / ( length_squared[1] + length_squared[0] );
00964     min_shape = VERDICT_MIN( shape, min_shape );
00965 
00966     shape     = corner_areas[2] / ( length_squared[2] + length_squared[1] );
00967     min_shape = VERDICT_MIN( shape, min_shape );
00968 
00969     shape     = corner_areas[3] / ( length_squared[3] + length_squared[2] );
00970     min_shape = VERDICT_MIN( shape, min_shape );
00971 
00972     min_shape *= 2;
00973 
00974     if( min_shape < VERDICT_DBL_MIN ) min_shape = 0;
00975 
00976     if( min_shape > 0 ) return (double)VERDICT_MIN( min_shape, VERDICT_DBL_MAX );
00977     return (double)VERDICT_MAX( min_shape, -VERDICT_DBL_MAX );
00978 }
00979 
00980 /*!
00981   the relative size of a quad
00982 
00983   Min( J, 1/J ), where J is determinant of weighted Jacobian matrix
00984 */
00985 C_FUNC_DEF double v_quad_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
00986 {
00987 
00988     double quad_area = v_quad_area( 4, coordinates );
00989     double rel_size  = 0;
00990 
00991     v_set_quad_size( quad_area );
00992     double w11, w21, w12, w22;
00993     get_weight( w11, w21, w12, w22 );
00994     double avg_area = determinant( w11, w21, w12, w22 );
00995 
00996     if( avg_area > VERDICT_DBL_MIN )
00997     {
00998 
00999         w11 = quad_area / avg_area;
01000 
01001         if( w11 > VERDICT_DBL_MIN )
01002         {
01003             rel_size = VERDICT_MIN( w11, 1 / w11 );
01004             rel_size *= rel_size;
01005         }
01006     }
01007 
01008     if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
01009     return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
01010 }
01011 
01012 /*!
01013   the relative shape and size of a quad
01014 
01015   Product of Shape and Relative Size
01016 */
01017 C_FUNC_DEF double v_quad_shape_and_size( int num_nodes, double coordinates[][3] )
01018 {
01019     double shape, size;
01020     size  = v_quad_relative_size_squared( num_nodes, coordinates );
01021     shape = v_quad_shape( num_nodes, coordinates );
01022 
01023     double shape_and_size = shape * size;
01024 
01025     if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
01026     return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
01027 }
01028 
01029 /*!
01030   the shear and size of a quad
01031 
01032   product of shear and relative size
01033 */
01034 C_FUNC_DEF double v_quad_shear_and_size( int num_nodes, double coordinates[][3] )
01035 {
01036     double shear, size;
01037     shear = v_quad_shear( num_nodes, coordinates );
01038     size  = v_quad_relative_size_squared( num_nodes, coordinates );
01039 
01040     double shear_and_size = shear * size;
01041 
01042     if( shear_and_size > 0 ) return (double)VERDICT_MIN( shear_and_size, VERDICT_DBL_MAX );
01043     return (double)VERDICT_MAX( shear_and_size, -VERDICT_DBL_MAX );
01044 }
01045 
01046 /*!
01047   the distortion of a quad
01048 */
01049 C_FUNC_DEF double v_quad_distortion( int num_nodes, double coordinates[][3] )
01050 {
01051     // To calculate distortion for linear and 2nd order quads
01052     // distortion = {min(|J|)/actual area}*{parent area}
01053     // parent area = 4 for a quad.
01054     // min |J| is the minimum over nodes and gaussian integration points
01055     // created by Ling Pan, CAT on 4/30/01
01056 
01057     double element_area = 0.0, distrt, thickness_gauss;
01058     double cur_jacobian = 0., sign_jacobian, jacobian;
01059     VerdictVector aa, bb, cc, normal_at_point, xin;
01060 
01061     // use 2x2 gauss points for linear quads and 3x3 for 2nd order quads
01062     int number_of_gauss_points = 0;
01063     if( num_nodes == 4 )
01064     {  // 2x2 quadrature rule
01065         number_of_gauss_points = 2;
01066     }
01067     else if( num_nodes == 8 )
01068     {  // 3x3 quadrature rule
01069         number_of_gauss_points = 3;
01070     }
01071 
01072     int total_number_of_gauss_points = number_of_gauss_points * number_of_gauss_points;
01073 
01074     VerdictVector face_normal = quad_normal( coordinates );
01075 
01076     double distortion = VERDICT_DBL_MAX;
01077 
01078     VerdictVector first, second;
01079 
01080     int i;
01081     // Will work out the case for collapsed quad later
01082     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE )
01083     {
01084         for( i = 0; i < 3; i++ )
01085         {
01086 
01087             first.set( coordinates[i][0] - coordinates[( i + 1 ) % 3][0],
01088                        coordinates[i][1] - coordinates[( i + 1 ) % 3][1],
01089                        coordinates[i][2] - coordinates[( i + 1 ) % 3][2] );
01090 
01091             second.set( coordinates[i][0] - coordinates[( i + 2 ) % 3][0],
01092                         coordinates[i][1] - coordinates[( i + 2 ) % 3][1],
01093                         coordinates[i][2] - coordinates[( i + 2 ) % 3][2] );
01094 
01095             sign_jacobian = ( face_normal % ( first * second ) ) > 0 ? 1. : -1.;
01096             cur_jacobian  = sign_jacobian * ( first * second ).length();
01097             distortion    = VERDICT_MIN( distortion, cur_jacobian );
01098         }
01099         element_area = ( first * second ).length() / 2.0;
01100         distortion /= element_area;
01101     }
01102     else
01103     {
01104         double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
01105         double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
01106         double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
01107         double weight[maxTotalNumberGaussPoints];
01108 
01109         // create an object of GaussIntegration
01110         GaussIntegration::initialize( number_of_gauss_points, num_nodes );
01111         GaussIntegration::calculate_shape_function_2d_quad();
01112         GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
01113 
01114         // calculate element area
01115         int ife, ja;
01116         for( ife = 0; ife < total_number_of_gauss_points; ife++ )
01117         {
01118             aa.set( 0.0, 0.0, 0.0 );
01119             bb.set( 0.0, 0.0, 0.0 );
01120 
01121             for( ja = 0; ja < num_nodes; ja++ )
01122             {
01123                 xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
01124                 aa += dndy1[ife][ja] * xin;
01125                 bb += dndy2[ife][ja] * xin;
01126             }
01127             normal_at_point = aa * bb;
01128             jacobian        = normal_at_point.length();
01129             element_area += weight[ife] * jacobian;
01130         }
01131 
01132         double dndy1_at_node[maxNumberNodes][maxNumberNodes];
01133         double dndy2_at_node[maxNumberNodes][maxNumberNodes];
01134 
01135         GaussIntegration::calculate_derivative_at_nodes( dndy1_at_node, dndy2_at_node );
01136 
01137         VerdictVector normal_at_nodes[9];
01138 
01139         // evaluate normal at nodes and distortion values at nodes
01140         int jai;
01141         for( ja = 0; ja < num_nodes; ja++ )
01142         {
01143             aa.set( 0.0, 0.0, 0.0 );
01144             bb.set( 0.0, 0.0, 0.0 );
01145             for( jai = 0; jai < num_nodes; jai++ )
01146             {
01147                 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
01148                 aa += dndy1_at_node[ja][jai] * xin;
01149                 bb += dndy2_at_node[ja][jai] * xin;
01150             }
01151             normal_at_nodes[ja] = aa * bb;
01152             normal_at_nodes[ja].normalize();
01153         }
01154 
01155         // determine if element is flat
01156         bool flat_element = true;
01157         double dot_product;
01158 
01159         for( ja = 0; ja < num_nodes; ja++ )
01160         {
01161             dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
01162             if( fabs( dot_product ) < 0.99 )
01163             {
01164                 flat_element = false;
01165                 break;
01166             }
01167         }
01168 
01169         // take into consideration of the thickness of the element
01170         double thickness;
01171         // get_quad_thickness(face, element_area, thickness );
01172         thickness = 0.001 * sqrt( element_area );
01173 
01174         // set thickness gauss point location
01175         double zl = 0.5773502691896;
01176         if( flat_element ) zl = 0.0;
01177 
01178         int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
01179         double thickness_z;
01180         int igz;
01181         // loop on Gauss points
01182         for( ife = 0; ife < total_number_of_gauss_points; ife++ )
01183         {
01184             // loop on the thickness direction gauss points
01185             for( igz = 0; igz < no_gauss_pts_z; igz++ )
01186             {
01187                 zl          = -zl;
01188                 thickness_z = zl * thickness / 2.0;
01189 
01190                 aa.set( 0.0, 0.0, 0.0 );
01191                 bb.set( 0.0, 0.0, 0.0 );
01192                 cc.set( 0.0, 0.0, 0.0 );
01193 
01194                 for( ja = 0; ja < num_nodes; ja++ )
01195                 {
01196                     xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
01197                     xin += thickness_z * normal_at_nodes[ja];
01198                     aa += dndy1[ife][ja] * xin;
01199                     bb += dndy2[ife][ja] * xin;
01200                     thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
01201                     cc += thickness_gauss * normal_at_nodes[ja];
01202                 }
01203 
01204                 normal_at_point = aa * bb;
01205                 // jacobian = normal_at_point.length();
01206                 distrt = cc % normal_at_point;
01207                 if( distrt < distortion ) distortion = distrt;
01208             }
01209         }
01210 
01211         // loop through nodal points
01212         for( ja = 0; ja < num_nodes; ja++ )
01213         {
01214             for( igz = 0; igz < no_gauss_pts_z; igz++ )
01215             {
01216                 zl          = -zl;
01217                 thickness_z = zl * thickness / 2.0;
01218 
01219                 aa.set( 0.0, 0.0, 0.0 );
01220                 bb.set( 0.0, 0.0, 0.0 );
01221                 cc.set( 0.0, 0.0, 0.0 );
01222 
01223                 for( jai = 0; jai < num_nodes; jai++ )
01224                 {
01225                     xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
01226                     xin += thickness_z * normal_at_nodes[ja];
01227                     aa += dndy1_at_node[ja][jai] * xin;
01228                     bb += dndy2_at_node[ja][jai] * xin;
01229                     if( jai == ja )
01230                         thickness_gauss = thickness / 2.0;
01231                     else
01232                         thickness_gauss = 0.;
01233                     cc += thickness_gauss * normal_at_nodes[jai];
01234                 }
01235             }
01236             normal_at_point = aa * bb;
01237             sign_jacobian   = ( face_normal % normal_at_point ) > 0 ? 1. : -1.;
01238             distrt          = sign_jacobian * ( cc % normal_at_point );
01239 
01240             if( distrt < distortion ) distortion = distrt;
01241         }
01242 
01243         if( element_area * thickness != 0 )
01244             distortion *= 8. / ( element_area * thickness );
01245         else
01246             distortion *= 8.;
01247     }
01248 
01249     return (double)distortion;
01250 }
01251 
01252 /*!
01253   multiple quality measures of a quad
01254 */
01255 C_FUNC_DEF void v_quad_quality( int num_nodes, double coordinates[][3], unsigned int metrics_request_flag,
01256                                 QuadMetricVals* metric_vals )
01257 {
01258 
01259     memset( metric_vals, 0, sizeof( QuadMetricVals ) );
01260 
01261     // for starts, lets set up some basic and common information
01262 
01263     /*  node numbers and side numbers used below
01264 
01265                     2
01266               3 +--------- 2
01267                /         +
01268               /          |
01269            3 /           | 1
01270             /            |
01271            +             |
01272          0 -------------+ 1
01273                0
01274     */
01275 
01276     // vectors for each side
01277     VerdictVector edges[4];
01278     make_quad_edges( edges, coordinates );
01279 
01280     double areas[4];
01281     signed_corner_areas( areas, coordinates );
01282 
01283     double lengths[4];
01284     lengths[0] = edges[0].length();
01285     lengths[1] = edges[1].length();
01286     lengths[2] = edges[2].length();
01287     lengths[3] = edges[3].length();
01288 
01289     VerdictBoolean is_collapsed = is_collapsed_quad( coordinates );
01290 
01291     // handle collapsed quads metrics here
01292     if( is_collapsed == VERDICT_TRUE && metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE |
01293                                                                  V_QUAD_JACOBIAN | V_QUAD_SCALED_JACOBIAN ) )
01294     {
01295         if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
01296             metric_vals->minimum_angle = v_tri_minimum_angle( 3, coordinates );
01297         if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
01298             metric_vals->maximum_angle = v_tri_maximum_angle( 3, coordinates );
01299         if( metrics_request_flag & V_QUAD_JACOBIAN )
01300             metric_vals->jacobian = (double)( v_tri_area( 3, coordinates ) * 2.0 );
01301         if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
01302             metric_vals->jacobian = (double)( v_tri_scaled_jacobian( 3, coordinates ) * 2.0 );
01303     }
01304 
01305     // calculate both largest and smallest angles
01306     if( metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE ) && is_collapsed == VERDICT_FALSE )
01307     {
01308         // gather the angles
01309         double angles[4];
01310         angles[0] = acos( -( edges[0] % edges[1] ) / ( lengths[0] * lengths[1] ) );
01311         angles[1] = acos( -( edges[1] % edges[2] ) / ( lengths[1] * lengths[2] ) );
01312         angles[2] = acos( -( edges[2] % edges[3] ) / ( lengths[2] * lengths[3] ) );
01313         angles[3] = acos( -( edges[3] % edges[0] ) / ( lengths[3] * lengths[0] ) );
01314 
01315         if( lengths[0] <= VERDICT_DBL_MIN || lengths[1] <= VERDICT_DBL_MIN || lengths[2] <= VERDICT_DBL_MIN ||
01316             lengths[3] <= VERDICT_DBL_MIN )
01317         {
01318             metric_vals->minimum_angle = 360.0;
01319             metric_vals->maximum_angle = 0.0;
01320         }
01321         else
01322         {
01323             // if smallest angle, find the smallest angle
01324             if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
01325             {
01326                 metric_vals->minimum_angle = VERDICT_DBL_MAX;
01327                 for( int i = 0; i < 4; i++ )
01328                     metric_vals->minimum_angle = VERDICT_MIN( angles[i], metric_vals->minimum_angle );
01329                 metric_vals->minimum_angle *= 180.0 / VERDICT_PI;
01330             }
01331             // if largest angle, find the largest angle
01332             if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
01333             {
01334                 metric_vals->maximum_angle = 0.0;
01335                 for( int i = 0; i < 4; i++ )
01336                     metric_vals->maximum_angle = VERDICT_MAX( angles[i], metric_vals->maximum_angle );
01337                 metric_vals->maximum_angle *= 180.0 / VERDICT_PI;
01338 
01339                 if( areas[0] < 0 || areas[1] < 0 || areas[2] < 0 || areas[3] < 0 )
01340                     metric_vals->maximum_angle = 360 - metric_vals->maximum_angle;
01341             }
01342         }
01343     }
01344 
01345     // handle max_edge_ratio, skew, taper, and area together
01346     if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
01347     {
01348         // get principle axes
01349         VerdictVector principal_axes[2];
01350         principal_axes[0] = edges[0] - edges[2];
01351         principal_axes[1] = edges[1] - edges[3];
01352 
01353         if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
01354         {
01355             double len1 = principal_axes[0].length();
01356             double len2 = principal_axes[1].length();
01357 
01358             // calculate the max_edge_ratio ratio
01359             if( metrics_request_flag & V_QUAD_MAX_EDGE_RATIO )
01360             {
01361                 if( len1 < VERDICT_DBL_MIN || len2 < VERDICT_DBL_MIN )
01362                     metric_vals->max_edge_ratio = VERDICT_DBL_MAX;
01363                 else
01364                     metric_vals->max_edge_ratio = VERDICT_MAX( len1 / len2, len2 / len1 );
01365             }
01366 
01367             // calculate the taper
01368             if( metrics_request_flag & V_QUAD_TAPER )
01369             {
01370                 double min_length = VERDICT_MIN( len1, len2 );
01371 
01372                 VerdictVector cross_derivative = edges[1] + edges[3];
01373 
01374                 if( min_length < VERDICT_DBL_MIN )
01375                     metric_vals->taper = VERDICT_DBL_MAX;
01376                 else
01377                     metric_vals->taper = cross_derivative.length() / min_length;
01378             }
01379 
01380             // calculate the skew
01381             if( metrics_request_flag & V_QUAD_SKEW )
01382             {
01383                 if( principal_axes[0].normalize() < VERDICT_DBL_MIN || principal_axes[1].normalize() < VERDICT_DBL_MIN )
01384                     metric_vals->skew = 0.0;
01385                 else
01386                     metric_vals->skew = fabs( principal_axes[0] % principal_axes[1] );
01387             }
01388         }
01389     }
01390 
01391     // calculate the area
01392     if( metrics_request_flag & ( V_QUAD_AREA | V_QUAD_RELATIVE_SIZE_SQUARED ) )
01393     { metric_vals->area = 0.25 * ( areas[0] + areas[1] + areas[2] + areas[3] ); }
01394 
01395     // calculate the relative size
01396     if( metrics_request_flag & ( V_QUAD_RELATIVE_SIZE_SQUARED | V_QUAD_SHAPE_AND_SIZE | V_QUAD_SHEAR_AND_SIZE ) )
01397     {
01398         double quad_area = v_quad_area( 4, coordinates );
01399         v_set_quad_size( quad_area );
01400         double w11, w21, w12, w22;
01401         get_weight( w11, w21, w12, w22 );
01402         double avg_area = determinant( w11, w21, w12, w22 );
01403 
01404         if( avg_area < VERDICT_DBL_MIN )
01405             metric_vals->relative_size_squared = 0.0;
01406         else
01407             metric_vals->relative_size_squared =
01408                 pow( VERDICT_MIN( metric_vals->area / avg_area, avg_area / metric_vals->area ), 2 );
01409     }
01410 
01411     // calculate the jacobian
01412     if( metrics_request_flag & V_QUAD_JACOBIAN )
01413     { metric_vals->jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) ); }
01414 
01415     if( metrics_request_flag & ( V_QUAD_SCALED_JACOBIAN | V_QUAD_SHEAR | V_QUAD_SHEAR_AND_SIZE ) )
01416     {
01417         double scaled_jac, min_scaled_jac = VERDICT_DBL_MAX;
01418 
01419         if( lengths[0] < VERDICT_DBL_MIN || lengths[1] < VERDICT_DBL_MIN || lengths[2] < VERDICT_DBL_MIN ||
01420             lengths[3] < VERDICT_DBL_MIN )
01421         {
01422             metric_vals->scaled_jacobian = 0.0;
01423             metric_vals->shear           = 0.0;
01424         }
01425         else
01426         {
01427             scaled_jac     = areas[0] / ( lengths[0] * lengths[3] );
01428             min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
01429 
01430             scaled_jac     = areas[1] / ( lengths[1] * lengths[0] );
01431             min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
01432 
01433             scaled_jac     = areas[2] / ( lengths[2] * lengths[1] );
01434             min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
01435 
01436             scaled_jac     = areas[3] / ( lengths[3] * lengths[2] );
01437             min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
01438 
01439             metric_vals->scaled_jacobian = min_scaled_jac;
01440 
01441             // what the heck...set shear as well
01442             if( min_scaled_jac <= VERDICT_DBL_MIN )
01443                 metric_vals->shear = 0.0;
01444             else
01445                 metric_vals->shear = min_scaled_jac;
01446         }
01447     }
01448 
01449     if( metrics_request_flag & ( V_QUAD_WARPAGE | V_QUAD_ODDY ) )
01450     {
01451         VerdictVector corner_normals[4];
01452         corner_normals[0] = edges[3] * edges[0];
01453         corner_normals[1] = edges[0] * edges[1];
01454         corner_normals[2] = edges[1] * edges[2];
01455         corner_normals[3] = edges[2] * edges[3];
01456 
01457         if( metrics_request_flag & V_QUAD_ODDY )
01458         {
01459             double oddy, max_oddy = 0.0;
01460 
01461             double diff, dot_prod;
01462 
01463             double length_squared[4];
01464             length_squared[0] = corner_normals[0].length_squared();
01465             length_squared[1] = corner_normals[1].length_squared();
01466             length_squared[2] = corner_normals[2].length_squared();
01467             length_squared[3] = corner_normals[3].length_squared();
01468 
01469             if( length_squared[0] < VERDICT_DBL_MIN || length_squared[1] < VERDICT_DBL_MIN ||
01470                 length_squared[2] < VERDICT_DBL_MIN || length_squared[3] < VERDICT_DBL_MIN )
01471                 metric_vals->oddy = VERDICT_DBL_MAX;
01472             else
01473             {
01474                 diff     = ( lengths[0] * lengths[0] ) - ( lengths[1] * lengths[1] );
01475                 dot_prod = edges[0] % edges[1];
01476                 oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[1] );
01477                 max_oddy = VERDICT_MAX( oddy, max_oddy );
01478 
01479                 diff     = ( lengths[1] * lengths[1] ) - ( lengths[2] * lengths[2] );
01480                 dot_prod = edges[1] % edges[2];
01481                 oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[2] );
01482                 max_oddy = VERDICT_MAX( oddy, max_oddy );
01483 
01484                 diff     = ( lengths[2] * lengths[2] ) - ( lengths[3] * lengths[3] );
01485                 dot_prod = edges[2] % edges[3];
01486                 oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[3] );
01487                 max_oddy = VERDICT_MAX( oddy, max_oddy );
01488 
01489                 diff     = ( lengths[3] * lengths[3] ) - ( lengths[0] * lengths[0] );
01490                 dot_prod = edges[3] % edges[0];
01491                 oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[0] );
01492                 max_oddy = VERDICT_MAX( oddy, max_oddy );
01493 
01494                 metric_vals->oddy = max_oddy;
01495             }
01496         }
01497 
01498         if( metrics_request_flag & V_QUAD_WARPAGE )
01499         {
01500             if( corner_normals[0].normalize() < VERDICT_DBL_MIN || corner_normals[1].normalize() < VERDICT_DBL_MIN ||
01501                 corner_normals[2].normalize() < VERDICT_DBL_MIN || corner_normals[3].normalize() < VERDICT_DBL_MIN )
01502                 metric_vals->warpage = VERDICT_DBL_MAX;
01503             else
01504             {
01505                 metric_vals->warpage =
01506                     pow( VERDICT_MIN( corner_normals[0] % corner_normals[2], corner_normals[1] % corner_normals[3] ),
01507                          3 );
01508             }
01509         }
01510     }
01511 
01512     if( metrics_request_flag & V_QUAD_STRETCH )
01513     {
01514         VerdictVector temp;
01515 
01516         temp.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
01517                   coordinates[2][2] - coordinates[0][2] );
01518         double diag02 = temp.length_squared();
01519 
01520         temp.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
01521                   coordinates[3][2] - coordinates[1][2] );
01522         double diag13 = temp.length_squared();
01523 
01524         static const double QUAD_STRETCH_FACTOR = sqrt( 2.0 );
01525 
01526         // 'diag02' is now the max diagonal of the quad
01527         diag02 = VERDICT_MAX( diag02, diag13 );
01528 
01529         if( diag02 < VERDICT_DBL_MIN )
01530             metric_vals->stretch = VERDICT_DBL_MAX;
01531         else
01532             metric_vals->stretch =
01533                 QUAD_STRETCH_FACTOR *
01534                 VERDICT_MIN( VERDICT_MIN( lengths[0], lengths[1] ), VERDICT_MIN( lengths[2], lengths[3] ) ) /
01535                 sqrt( diag02 );
01536     }
01537 
01538     if( metrics_request_flag & ( V_QUAD_CONDITION | V_QUAD_SHAPE | V_QUAD_SHAPE_AND_SIZE ) )
01539     {
01540         double lengths_squared[4];
01541         lengths_squared[0] = edges[0].length_squared();
01542         lengths_squared[1] = edges[1].length_squared();
01543         lengths_squared[2] = edges[2].length_squared();
01544         lengths_squared[3] = edges[3].length_squared();
01545 
01546         if( areas[0] < VERDICT_DBL_MIN || areas[1] < VERDICT_DBL_MIN || areas[2] < VERDICT_DBL_MIN ||
01547             areas[3] < VERDICT_DBL_MIN )
01548         {
01549             metric_vals->condition = VERDICT_DBL_MAX;
01550             metric_vals->shape     = 0.0;
01551         }
01552         else
01553         {
01554             double max_condition = 0.0, condition;
01555 
01556             condition     = ( lengths_squared[0] + lengths_squared[3] ) / areas[0];
01557             max_condition = VERDICT_MAX( max_condition, condition );
01558 
01559             condition     = ( lengths_squared[1] + lengths_squared[0] ) / areas[1];
01560             max_condition = VERDICT_MAX( max_condition, condition );
01561 
01562             condition     = ( lengths_squared[2] + lengths_squared[1] ) / areas[2];
01563             max_condition = VERDICT_MAX( max_condition, condition );
01564 
01565             condition     = ( lengths_squared[3] + lengths_squared[2] ) / areas[3];
01566             max_condition = VERDICT_MAX( max_condition, condition );
01567 
01568             metric_vals->condition = 0.5 * max_condition;
01569             metric_vals->shape     = 2 / max_condition;
01570         }
01571     }
01572 
01573     if( metrics_request_flag & V_QUAD_AREA )
01574     {
01575         if( metric_vals->area > 0 ) metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
01576         metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
01577     }
01578 
01579     if( metrics_request_flag & V_QUAD_MAX_EDGE_RATIO )
01580     {
01581         if( metric_vals->max_edge_ratio > 0 )
01582             metric_vals->max_edge_ratio = (double)VERDICT_MIN( metric_vals->max_edge_ratio, VERDICT_DBL_MAX );
01583         metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX );
01584     }
01585 
01586     if( metrics_request_flag & V_QUAD_CONDITION )
01587     {
01588         if( metric_vals->condition > 0 )
01589             metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
01590         metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
01591     }
01592 
01593     // calculate distortion
01594     if( metrics_request_flag & V_QUAD_DISTORTION )
01595     {
01596         metric_vals->distortion = v_quad_distortion( num_nodes, coordinates );
01597 
01598         if( metric_vals->distortion > 0 )
01599             metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
01600         metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
01601     }
01602 
01603     if( metrics_request_flag & V_QUAD_JACOBIAN )
01604     {
01605         if( metric_vals->jacobian > 0 )
01606             metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
01607         metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
01608     }
01609 
01610     if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
01611     {
01612         if( metric_vals->maximum_angle > 0 )
01613             metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
01614         metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
01615     }
01616 
01617     if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
01618     {
01619         if( metric_vals->minimum_angle > 0 )
01620             metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
01621         metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
01622     }
01623 
01624     if( metrics_request_flag & V_QUAD_ODDY )
01625     {
01626         if( metric_vals->oddy > 0 ) metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX );
01627         metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX );
01628     }
01629 
01630     if( metrics_request_flag & V_QUAD_RELATIVE_SIZE_SQUARED )
01631     {
01632         if( metric_vals->relative_size_squared > 0 )
01633             metric_vals->relative_size_squared =
01634                 (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
01635         metric_vals->relative_size_squared =
01636             (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
01637     }
01638 
01639     if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
01640     {
01641         if( metric_vals->scaled_jacobian > 0 )
01642             metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
01643         metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
01644     }
01645 
01646     if( metrics_request_flag & V_QUAD_SHEAR )
01647     {
01648         if( metric_vals->shear > 0 ) metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX );
01649         metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX );
01650     }
01651 
01652     // calculate shear and size
01653     // reuse values from above
01654     if( metrics_request_flag & V_QUAD_SHEAR_AND_SIZE )
01655     {
01656         metric_vals->shear_and_size = metric_vals->shear * metric_vals->relative_size_squared;
01657 
01658         if( metric_vals->shear_and_size > 0 )
01659             metric_vals->shear_and_size = (double)VERDICT_MIN( metric_vals->shear_and_size, VERDICT_DBL_MAX );
01660         metric_vals->shear_and_size = (double)VERDICT_MAX( metric_vals->shear_and_size, -VERDICT_DBL_MAX );
01661     }
01662 
01663     if( metrics_request_flag & V_QUAD_SHAPE )
01664     {
01665         if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
01666         metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
01667     }
01668 
01669     // calculate shape and size
01670     // reuse values from above
01671     if( metrics_request_flag & V_QUAD_SHAPE_AND_SIZE )
01672     {
01673         metric_vals->shape_and_size = metric_vals->shape * metric_vals->relative_size_squared;
01674 
01675         if( metric_vals->shape_and_size > 0 )
01676             metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
01677         metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
01678     }
01679 
01680     if( metrics_request_flag & V_QUAD_SKEW )
01681     {
01682         if( metric_vals->skew > 0 ) metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX );
01683         metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX );
01684     }
01685 
01686     if( metrics_request_flag & V_QUAD_STRETCH )
01687     {
01688         if( metric_vals->stretch > 0 )
01689             metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX );
01690         metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX );
01691     }
01692 
01693     if( metrics_request_flag & V_QUAD_TAPER )
01694     {
01695         if( metric_vals->taper > 0 ) metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX );
01696         metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX );
01697     }
01698 
01699     if( metrics_request_flag & V_QUAD_WARPAGE )
01700     {
01701         if( metric_vals->warpage > 0 )
01702             metric_vals->warpage = (double)VERDICT_MIN( metric_vals->warpage, VERDICT_DBL_MAX );
01703         metric_vals->warpage = (double)VERDICT_MAX( metric_vals->warpage, -VERDICT_DBL_MAX );
01704     }
01705 
01706     if( metrics_request_flag & V_QUAD_MED_ASPECT_FROBENIUS )
01707         metric_vals->med_aspect_frobenius = v_quad_med_aspect_frobenius( 4, coordinates );
01708     if( metrics_request_flag & V_QUAD_MAX_ASPECT_FROBENIUS )
01709         metric_vals->max_aspect_frobenius = v_quad_max_aspect_frobenius( 4, coordinates );
01710     if( metrics_request_flag & V_QUAD_EDGE_RATIO ) metric_vals->edge_ratio = v_quad_edge_ratio( 4, coordinates );
01711     if( metrics_request_flag & V_QUAD_ASPECT_RATIO ) metric_vals->aspect_ratio = v_quad_aspect_ratio( 4, coordinates );
01712     if( metrics_request_flag & V_QUAD_RADIUS_RATIO ) metric_vals->radius_ratio = v_quad_radius_ratio( 4, coordinates );
01713 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines