Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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 )
00721     {
00722         max_angle = 360 - max_angle;
00723     }
00724 
00725     if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
00726     return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
00727 }
00728 
00729 /*!
00730   the smallest angle of a quad
00731 
00732   smallest included quad angle (degrees)
00733 */
00734 C_FUNC_DEF double v_quad_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
00735 {
00736     // if this quad is a collapsed quad, then just
00737     // send it to the tri_smallest_angle routine
00738     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_minimum_angle( 3, coordinates );
00739 
00740     double angle;
00741     double min_angle = 360.0;
00742 
00743     VerdictVector edges[4];
00744     edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00745                   coordinates[1][2] - coordinates[0][2] );
00746     edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00747                   coordinates[2][2] - coordinates[1][2] );
00748     edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00749                   coordinates[3][2] - coordinates[2][2] );
00750     edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
00751                   coordinates[0][2] - coordinates[3][2] );
00752 
00753     // go around each node and calculate the angle
00754     // at each node
00755     double length[4];
00756     length[0] = edges[0].length();
00757     length[1] = edges[1].length();
00758     length[2] = edges[2].length();
00759     length[3] = edges[3].length();
00760 
00761     if( length[0] <= VERDICT_DBL_MIN || length[1] <= VERDICT_DBL_MIN || length[2] <= VERDICT_DBL_MIN ||
00762         length[3] <= VERDICT_DBL_MIN )
00763         return 360.0;
00764 
00765     angle     = acos( -( edges[0] % edges[1] ) / ( length[0] * length[1] ) );
00766     min_angle = VERDICT_MIN( angle, min_angle );
00767 
00768     angle     = acos( -( edges[1] % edges[2] ) / ( length[1] * length[2] ) );
00769     min_angle = VERDICT_MIN( angle, min_angle );
00770 
00771     angle     = acos( -( edges[2] % edges[3] ) / ( length[2] * length[3] ) );
00772     min_angle = VERDICT_MIN( angle, min_angle );
00773 
00774     angle     = acos( -( edges[3] % edges[0] ) / ( length[3] * length[0] ) );
00775     min_angle = VERDICT_MIN( angle, min_angle );
00776 
00777     min_angle = min_angle * 180.0 / VERDICT_PI;
00778 
00779     if( min_angle > 0 ) return (double)VERDICT_MIN( min_angle, VERDICT_DBL_MAX );
00780     return (double)VERDICT_MAX( min_angle, -VERDICT_DBL_MAX );
00781 }
00782 
00783 /*!
00784   the oddy of a quad
00785 
00786   general distortion measure based on left Cauchy-Green Tensor
00787 */
00788 C_FUNC_DEF double v_quad_oddy( int /*num_nodes*/, double coordinates[][3] )
00789 {
00790 
00791     double max_oddy = 0.;
00792 
00793     VerdictVector first, second, node_pos[4];
00794 
00795     double g, g11, g12, g22, cur_oddy;
00796     int i;
00797 
00798     for( i = 0; i < 4; i++ )
00799         node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );
00800 
00801     for( i = 0; i < 4; i++ )
00802     {
00803         first  = node_pos[i] - node_pos[( i + 1 ) % 4];
00804         second = node_pos[i] - node_pos[( i + 3 ) % 4];
00805 
00806         g11 = first % first;
00807         g12 = first % second;
00808         g22 = second % second;
00809         g   = g11 * g22 - g12 * g12;
00810 
00811         if( g < VERDICT_DBL_MIN )
00812             cur_oddy = VERDICT_DBL_MAX;
00813         else
00814             cur_oddy = ( ( g11 - g22 ) * ( g11 - g22 ) + 4. * g12 * g12 ) / 2. / g;
00815 
00816         max_oddy = VERDICT_MAX( max_oddy, cur_oddy );
00817     }
00818 
00819     if( max_oddy > 0 ) return (double)VERDICT_MIN( max_oddy, VERDICT_DBL_MAX );
00820     return (double)VERDICT_MAX( max_oddy, -VERDICT_DBL_MAX );
00821 }
00822 
00823 /*!
00824   the condition of a quad
00825 
00826   maximum condition number of the Jacobian matrix at 4 corners
00827 */
00828 C_FUNC_DEF double v_quad_condition( int /*num_nodes*/, double coordinates[][3] )
00829 {
00830 
00831     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_condition( 3, coordinates );
00832 
00833     double areas[4];
00834     signed_corner_areas( areas, coordinates );
00835 
00836     double max_condition = 0.;
00837 
00838     VerdictVector xxi, xet;
00839 
00840     double condition;
00841 
00842     for( int i = 0; i < 4; i++ )
00843     {
00844 
00845         xxi.set( coordinates[i][0] - coordinates[( i + 1 ) % 4][0], coordinates[i][1] - coordinates[( i + 1 ) % 4][1],
00846                  coordinates[i][2] - coordinates[( i + 1 ) % 4][2] );
00847 
00848         xet.set( coordinates[i][0] - coordinates[( i + 3 ) % 4][0], coordinates[i][1] - coordinates[( i + 3 ) % 4][1],
00849                  coordinates[i][2] - coordinates[( i + 3 ) % 4][2] );
00850 
00851         if( areas[i] < VERDICT_DBL_MIN )
00852             condition = VERDICT_DBL_MAX;
00853         else
00854             condition = ( xxi % xxi + xet % xet ) / areas[i];
00855 
00856         max_condition = VERDICT_MAX( max_condition, condition );
00857     }
00858 
00859     max_condition /= 2;
00860 
00861     if( max_condition > 0 ) return (double)VERDICT_MIN( max_condition, VERDICT_DBL_MAX );
00862     return (double)VERDICT_MAX( max_condition, -VERDICT_DBL_MAX );
00863 }
00864 
00865 /*!
00866   the jacobian of a quad
00867 
00868   minimum pointwise volume of local map at 4 corners and center of quad
00869 */
00870 C_FUNC_DEF double v_quad_jacobian( int /*num_nodes*/, double coordinates[][3] )
00871 {
00872 
00873     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return (double)( v_tri_area( 3, coordinates ) * 2.0 );
00874 
00875     double areas[4];
00876     signed_corner_areas( areas, coordinates );
00877 
00878     double jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) );
00879     if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
00880     return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
00881 }
00882 
00883 /*!
00884   scaled jacobian of a quad
00885 
00886   Minimum Jacobian divided by the lengths of the 2 edge vector
00887 */
00888 C_FUNC_DEF double v_quad_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
00889 {
00890     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_scaled_jacobian( 3, coordinates );
00891 
00892     double corner_areas[4], min_scaled_jac = VERDICT_DBL_MAX, scaled_jac;
00893     signed_corner_areas( corner_areas, coordinates );
00894 
00895     VerdictVector edges[4];
00896     make_quad_edges( edges, coordinates );
00897 
00898     double length[4];
00899     length[0] = edges[0].length();
00900     length[1] = edges[1].length();
00901     length[2] = edges[2].length();
00902     length[3] = edges[3].length();
00903 
00904     if( length[0] < VERDICT_DBL_MIN || length[1] < VERDICT_DBL_MIN || length[2] < VERDICT_DBL_MIN ||
00905         length[3] < VERDICT_DBL_MIN )
00906         return 0.0;
00907 
00908     scaled_jac     = corner_areas[0] / ( length[0] * length[3] );
00909     min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
00910 
00911     scaled_jac     = corner_areas[1] / ( length[1] * length[0] );
00912     min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
00913 
00914     scaled_jac     = corner_areas[2] / ( length[2] * length[1] );
00915     min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
00916 
00917     scaled_jac     = corner_areas[3] / ( length[3] * length[2] );
00918     min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
00919 
00920     if( min_scaled_jac > 0 ) return (double)VERDICT_MIN( min_scaled_jac, VERDICT_DBL_MAX );
00921     return (double)VERDICT_MAX( min_scaled_jac, -VERDICT_DBL_MAX );
00922 }
00923 
00924 /*!
00925   the shear of a quad
00926 
00927   2/Condition number of Jacobian Skew matrix
00928 */
00929 C_FUNC_DEF double v_quad_shear( int /*num_nodes*/, double coordinates[][3] )
00930 {
00931     double scaled_jacobian = v_quad_scaled_jacobian( 4, coordinates );
00932 
00933     if( scaled_jacobian <= VERDICT_DBL_MIN )
00934         return 0.0;
00935     else
00936         return (double)VERDICT_MIN( scaled_jacobian, VERDICT_DBL_MAX );
00937 }
00938 
00939 /*!
00940   the shape of a quad
00941 
00942    2/Condition number of weighted Jacobian matrix
00943 */
00944 C_FUNC_DEF double v_quad_shape( int /*num_nodes*/, double coordinates[][3] )
00945 {
00946 
00947     double corner_areas[4], min_shape = VERDICT_DBL_MAX, shape;
00948     signed_corner_areas( corner_areas, coordinates );
00949 
00950     VerdictVector edges[4];
00951     make_quad_edges( edges, coordinates );
00952 
00953     double length_squared[4];
00954     length_squared[0] = edges[0].length_squared();
00955     length_squared[1] = edges[1].length_squared();
00956     length_squared[2] = edges[2].length_squared();
00957     length_squared[3] = edges[3].length_squared();
00958 
00959     if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN ||
00960         length_squared[2] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN )
00961         return 0.0;
00962 
00963     shape     = corner_areas[0] / ( length_squared[0] + length_squared[3] );
00964     min_shape = VERDICT_MIN( shape, min_shape );
00965 
00966     shape     = corner_areas[1] / ( length_squared[1] + length_squared[0] );
00967     min_shape = VERDICT_MIN( shape, min_shape );
00968 
00969     shape     = corner_areas[2] / ( length_squared[2] + length_squared[1] );
00970     min_shape = VERDICT_MIN( shape, min_shape );
00971 
00972     shape     = corner_areas[3] / ( length_squared[3] + length_squared[2] );
00973     min_shape = VERDICT_MIN( shape, min_shape );
00974 
00975     min_shape *= 2;
00976 
00977     if( min_shape < VERDICT_DBL_MIN ) min_shape = 0;
00978 
00979     if( min_shape > 0 ) return (double)VERDICT_MIN( min_shape, VERDICT_DBL_MAX );
00980     return (double)VERDICT_MAX( min_shape, -VERDICT_DBL_MAX );
00981 }
00982 
00983 /*!
00984   the relative size of a quad
00985 
00986   Min( J, 1/J ), where J is determinant of weighted Jacobian matrix
00987 */
00988 C_FUNC_DEF double v_quad_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
00989 {
00990 
00991     double quad_area = v_quad_area( 4, coordinates );
00992     double rel_size  = 0;
00993 
00994     v_set_quad_size( quad_area );
00995     double w11, w21, w12, w22;
00996     get_weight( w11, w21, w12, w22 );
00997     double avg_area = determinant( w11, w21, w12, w22 );
00998 
00999     if( avg_area > VERDICT_DBL_MIN )
01000     {
01001 
01002         w11 = quad_area / avg_area;
01003 
01004         if( w11 > VERDICT_DBL_MIN )
01005         {
01006             rel_size = VERDICT_MIN( w11, 1 / w11 );
01007             rel_size *= rel_size;
01008         }
01009     }
01010 
01011     if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
01012     return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
01013 }
01014 
01015 /*!
01016   the relative shape and size of a quad
01017 
01018   Product of Shape and Relative Size
01019 */
01020 C_FUNC_DEF double v_quad_shape_and_size( int num_nodes, double coordinates[][3] )
01021 {
01022     double shape, size;
01023     size  = v_quad_relative_size_squared( num_nodes, coordinates );
01024     shape = v_quad_shape( num_nodes, coordinates );
01025 
01026     double shape_and_size = shape * size;
01027 
01028     if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
01029     return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
01030 }
01031 
01032 /*!
01033   the shear and size of a quad
01034 
01035   product of shear and relative size
01036 */
01037 C_FUNC_DEF double v_quad_shear_and_size( int num_nodes, double coordinates[][3] )
01038 {
01039     double shear, size;
01040     shear = v_quad_shear( num_nodes, coordinates );
01041     size  = v_quad_relative_size_squared( num_nodes, coordinates );
01042 
01043     double shear_and_size = shear * size;
01044 
01045     if( shear_and_size > 0 ) return (double)VERDICT_MIN( shear_and_size, VERDICT_DBL_MAX );
01046     return (double)VERDICT_MAX( shear_and_size, -VERDICT_DBL_MAX );
01047 }
01048 
01049 /*!
01050   the distortion of a quad
01051 */
01052 C_FUNC_DEF double v_quad_distortion( int num_nodes, double coordinates[][3] )
01053 {
01054     // To calculate distortion for linear and 2nd order quads
01055     // distortion = {min(|J|)/actual area}*{parent area}
01056     // parent area = 4 for a quad.
01057     // min |J| is the minimum over nodes and gaussian integration points
01058     // created by Ling Pan, CAT on 4/30/01
01059 
01060     double element_area = 0.0, distrt, thickness_gauss;
01061     double cur_jacobian = 0., sign_jacobian, jacobian;
01062     VerdictVector aa, bb, cc, normal_at_point, xin;
01063 
01064     // use 2x2 gauss points for linear quads and 3x3 for 2nd order quads
01065     int number_of_gauss_points = 0;
01066     if( num_nodes == 4 )
01067     {  // 2x2 quadrature rule
01068         number_of_gauss_points = 2;
01069     }
01070     else if( num_nodes == 8 )
01071     {  // 3x3 quadrature rule
01072         number_of_gauss_points = 3;
01073     }
01074 
01075     int total_number_of_gauss_points = number_of_gauss_points * number_of_gauss_points;
01076 
01077     VerdictVector face_normal = quad_normal( coordinates );
01078 
01079     double distortion = VERDICT_DBL_MAX;
01080 
01081     VerdictVector first, second;
01082 
01083     int i;
01084     // Will work out the case for collapsed quad later
01085     if( is_collapsed_quad( coordinates ) == VERDICT_TRUE )
01086     {
01087         for( i = 0; i < 3; i++ )
01088         {
01089 
01090             first.set( coordinates[i][0] - coordinates[( i + 1 ) % 3][0],
01091                        coordinates[i][1] - coordinates[( i + 1 ) % 3][1],
01092                        coordinates[i][2] - coordinates[( i + 1 ) % 3][2] );
01093 
01094             second.set( coordinates[i][0] - coordinates[( i + 2 ) % 3][0],
01095                         coordinates[i][1] - coordinates[( i + 2 ) % 3][1],
01096                         coordinates[i][2] - coordinates[( i + 2 ) % 3][2] );
01097 
01098             sign_jacobian = ( face_normal % ( first * second ) ) > 0 ? 1. : -1.;
01099             cur_jacobian  = sign_jacobian * ( first * second ).length();
01100             distortion    = VERDICT_MIN( distortion, cur_jacobian );
01101         }
01102         element_area = ( first * second ).length() / 2.0;
01103         distortion /= element_area;
01104     }
01105     else
01106     {
01107         double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
01108         double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
01109         double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
01110         double weight[maxTotalNumberGaussPoints];
01111 
01112         // create an object of GaussIntegration
01113         GaussIntegration::initialize( number_of_gauss_points, num_nodes );
01114         GaussIntegration::calculate_shape_function_2d_quad();
01115         GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
01116 
01117         // calculate element area
01118         int ife, ja;
01119         for( ife = 0; ife < total_number_of_gauss_points; ife++ )
01120         {
01121             aa.set( 0.0, 0.0, 0.0 );
01122             bb.set( 0.0, 0.0, 0.0 );
01123 
01124             for( ja = 0; ja < num_nodes; ja++ )
01125             {
01126                 xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
01127                 aa += dndy1[ife][ja] * xin;
01128                 bb += dndy2[ife][ja] * xin;
01129             }
01130             normal_at_point = aa * bb;
01131             jacobian        = normal_at_point.length();
01132             element_area += weight[ife] * jacobian;
01133         }
01134 
01135         double dndy1_at_node[maxNumberNodes][maxNumberNodes];
01136         double dndy2_at_node[maxNumberNodes][maxNumberNodes];
01137 
01138         GaussIntegration::calculate_derivative_at_nodes( dndy1_at_node, dndy2_at_node );
01139 
01140         VerdictVector normal_at_nodes[9];
01141 
01142         // evaluate normal at nodes and distortion values at nodes
01143         int jai;
01144         for( ja = 0; ja < num_nodes; ja++ )
01145         {
01146             aa.set( 0.0, 0.0, 0.0 );
01147             bb.set( 0.0, 0.0, 0.0 );
01148             for( jai = 0; jai < num_nodes; jai++ )
01149             {
01150                 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
01151                 aa += dndy1_at_node[ja][jai] * xin;
01152                 bb += dndy2_at_node[ja][jai] * xin;
01153             }
01154             normal_at_nodes[ja] = aa * bb;
01155             normal_at_nodes[ja].normalize();
01156         }
01157 
01158         // determine if element is flat
01159         bool flat_element = true;
01160         double dot_product;
01161 
01162         for( ja = 0; ja < num_nodes; ja++ )
01163         {
01164             dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
01165             if( fabs( dot_product ) < 0.99 )
01166             {
01167                 flat_element = false;
01168                 break;
01169             }
01170         }
01171 
01172         // take into consideration of the thickness of the element
01173         double thickness;
01174         // get_quad_thickness(face, element_area, thickness );
01175         thickness = 0.001 * sqrt( element_area );
01176 
01177         // set thickness gauss point location
01178         double zl = 0.5773502691896;
01179         if( flat_element ) zl = 0.0;
01180 
01181         int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
01182         double thickness_z;
01183         int igz;
01184         // loop on Gauss points
01185         for( ife = 0; ife < total_number_of_gauss_points; ife++ )
01186         {
01187             // loop on the thickness direction gauss points
01188             for( igz = 0; igz < no_gauss_pts_z; igz++ )
01189             {
01190                 zl          = -zl;
01191                 thickness_z = zl * thickness / 2.0;
01192 
01193                 aa.set( 0.0, 0.0, 0.0 );
01194                 bb.set( 0.0, 0.0, 0.0 );
01195                 cc.set( 0.0, 0.0, 0.0 );
01196 
01197                 for( ja = 0; ja < num_nodes; ja++ )
01198                 {
01199                     xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
01200                     xin += thickness_z * normal_at_nodes[ja];
01201                     aa += dndy1[ife][ja] * xin;
01202                     bb += dndy2[ife][ja] * xin;
01203                     thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
01204                     cc += thickness_gauss * normal_at_nodes[ja];
01205                 }
01206 
01207                 normal_at_point = aa * bb;
01208                 // jacobian = normal_at_point.length();
01209                 distrt = cc % normal_at_point;
01210                 if( distrt < distortion ) distortion = distrt;
01211             }
01212         }
01213 
01214         // loop through nodal points
01215         for( ja = 0; ja < num_nodes; ja++ )
01216         {
01217             for( igz = 0; igz < no_gauss_pts_z; igz++ )
01218             {
01219                 zl          = -zl;
01220                 thickness_z = zl * thickness / 2.0;
01221 
01222                 aa.set( 0.0, 0.0, 0.0 );
01223                 bb.set( 0.0, 0.0, 0.0 );
01224                 cc.set( 0.0, 0.0, 0.0 );
01225 
01226                 for( jai = 0; jai < num_nodes; jai++ )
01227                 {
01228                     xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
01229                     xin += thickness_z * normal_at_nodes[ja];
01230                     aa += dndy1_at_node[ja][jai] * xin;
01231                     bb += dndy2_at_node[ja][jai] * xin;
01232                     if( jai == ja )
01233                         thickness_gauss = thickness / 2.0;
01234                     else
01235                         thickness_gauss = 0.;
01236                     cc += thickness_gauss * normal_at_nodes[jai];
01237                 }
01238             }
01239             normal_at_point = aa * bb;
01240             sign_jacobian   = ( face_normal % normal_at_point ) > 0 ? 1. : -1.;
01241             distrt          = sign_jacobian * ( cc % normal_at_point );
01242 
01243             if( distrt < distortion ) distortion = distrt;
01244         }
01245 
01246         if( element_area * thickness != 0 )
01247             distortion *= 8. / ( element_area * thickness );
01248         else
01249             distortion *= 8.;
01250     }
01251 
01252     return (double)distortion;
01253 }
01254 
01255 /*!
01256   multiple quality measures of a quad
01257 */
01258 C_FUNC_DEF void v_quad_quality( int num_nodes,
01259                                 double coordinates[][3],
01260                                 unsigned int metrics_request_flag,
01261                                 QuadMetricVals* metric_vals )
01262 {
01263 
01264     memset( metric_vals, 0, sizeof( QuadMetricVals ) );
01265 
01266     // for starts, lets set up some basic and common information
01267 
01268     /*  node numbers and side numbers used below
01269 
01270                     2
01271               3 +--------- 2
01272                /         +
01273               /          |
01274            3 /           | 1
01275             /            |
01276            +             |
01277          0 -------------+ 1
01278                0
01279     */
01280 
01281     // vectors for each side
01282     VerdictVector edges[4];
01283     make_quad_edges( edges, coordinates );
01284 
01285     double areas[4];
01286     signed_corner_areas( areas, coordinates );
01287 
01288     double lengths[4];
01289     lengths[0] = edges[0].length();
01290     lengths[1] = edges[1].length();
01291     lengths[2] = edges[2].length();
01292     lengths[3] = edges[3].length();
01293 
01294     VerdictBoolean is_collapsed = is_collapsed_quad( coordinates );
01295 
01296     // handle collapsed quads metrics here
01297     if( is_collapsed == VERDICT_TRUE && metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE |
01298                                                                  V_QUAD_JACOBIAN | V_QUAD_SCALED_JACOBIAN ) )
01299     {
01300         if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
01301             metric_vals->minimum_angle = v_tri_minimum_angle( 3, coordinates );
01302         if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
01303             metric_vals->maximum_angle = v_tri_maximum_angle( 3, coordinates );
01304         if( metrics_request_flag & V_QUAD_JACOBIAN )
01305             metric_vals->jacobian = (double)( v_tri_area( 3, coordinates ) * 2.0 );
01306         if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
01307             metric_vals->jacobian = (double)( v_tri_scaled_jacobian( 3, coordinates ) * 2.0 );
01308     }
01309 
01310     // calculate both largest and smallest angles
01311     if( metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE ) && is_collapsed == VERDICT_FALSE )
01312     {
01313         // gather the angles
01314         double angles[4];
01315         angles[0] = acos( -( edges[0] % edges[1] ) / ( lengths[0] * lengths[1] ) );
01316         angles[1] = acos( -( edges[1] % edges[2] ) / ( lengths[1] * lengths[2] ) );
01317         angles[2] = acos( -( edges[2] % edges[3] ) / ( lengths[2] * lengths[3] ) );
01318         angles[3] = acos( -( edges[3] % edges[0] ) / ( lengths[3] * lengths[0] ) );
01319 
01320         if( lengths[0] <= VERDICT_DBL_MIN || lengths[1] <= VERDICT_DBL_MIN || lengths[2] <= VERDICT_DBL_MIN ||
01321             lengths[3] <= VERDICT_DBL_MIN )
01322         {
01323             metric_vals->minimum_angle = 360.0;
01324             metric_vals->maximum_angle = 0.0;
01325         }
01326         else
01327         {
01328             // if smallest angle, find the smallest angle
01329             if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
01330             {
01331                 metric_vals->minimum_angle = VERDICT_DBL_MAX;
01332                 for( int i = 0; i < 4; i++ )
01333                     metric_vals->minimum_angle = VERDICT_MIN( angles[i], metric_vals->minimum_angle );
01334                 metric_vals->minimum_angle *= 180.0 / VERDICT_PI;
01335             }
01336             // if largest angle, find the largest angle
01337             if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
01338             {
01339                 metric_vals->maximum_angle = 0.0;
01340                 for( int i = 0; i < 4; i++ )
01341                     metric_vals->maximum_angle = VERDICT_MAX( angles[i], metric_vals->maximum_angle );
01342                 metric_vals->maximum_angle *= 180.0 / VERDICT_PI;
01343 
01344                 if( areas[0] < 0 || areas[1] < 0 || areas[2] < 0 || areas[3] < 0 )
01345                     metric_vals->maximum_angle = 360 - metric_vals->maximum_angle;
01346             }
01347         }
01348     }
01349 
01350     // handle max_edge_ratio, skew, taper, and area together
01351     if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
01352     {
01353         // get principle axes
01354         VerdictVector principal_axes[2];
01355         principal_axes[0] = edges[0] - edges[2];
01356         principal_axes[1] = edges[1] - edges[3];
01357 
01358         if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
01359         {
01360             double len1 = principal_axes[0].length();
01361             double len2 = principal_axes[1].length();
01362 
01363             // calculate the max_edge_ratio ratio
01364             if( metrics_request_flag & V_QUAD_MAX_EDGE_RATIO )
01365             {
01366                 if( len1 < VERDICT_DBL_MIN || len2 < VERDICT_DBL_MIN )
01367                     metric_vals->max_edge_ratio = VERDICT_DBL_MAX;
01368                 else
01369                     metric_vals->max_edge_ratio = VERDICT_MAX( len1 / len2, len2 / len1 );
01370             }
01371 
01372             // calculate the taper
01373             if( metrics_request_flag & V_QUAD_TAPER )
01374             {
01375                 double min_length = VERDICT_MIN( len1, len2 );
01376 
01377                 VerdictVector cross_derivative = edges[1] + edges[3];
01378 
01379                 if( min_length < VERDICT_DBL_MIN )
01380                     metric_vals->taper = VERDICT_DBL_MAX;
01381                 else
01382                     metric_vals->taper = cross_derivative.length() / min_length;
01383             }
01384 
01385             // calculate the skew
01386             if( metrics_request_flag & V_QUAD_SKEW )
01387             {
01388                 if( principal_axes[0].normalize() < VERDICT_DBL_MIN || principal_axes[1].normalize() < VERDICT_DBL_MIN )
01389                     metric_vals->skew = 0.0;
01390                 else
01391                     metric_vals->skew = fabs( principal_axes[0] % principal_axes[1] );
01392             }
01393         }
01394     }
01395 
01396     // calculate the area
01397     if( metrics_request_flag & ( V_QUAD_AREA | V_QUAD_RELATIVE_SIZE_SQUARED ) )
01398     {
01399         metric_vals->area = 0.25 * ( areas[0] + areas[1] + areas[2] + areas[3] );
01400     }
01401 
01402     // calculate the relative size
01403     if( metrics_request_flag & ( V_QUAD_RELATIVE_SIZE_SQUARED | V_QUAD_SHAPE_AND_SIZE | V_QUAD_SHEAR_AND_SIZE ) )
01404     {
01405         double quad_area = v_quad_area( 4, coordinates );
01406         v_set_quad_size( quad_area );
01407         double w11, w21, w12, w22;
01408         get_weight( w11, w21, w12, w22 );
01409         double avg_area = determinant( w11, w21, w12, w22 );
01410 
01411         if( avg_area < VERDICT_DBL_MIN )
01412             metric_vals->relative_size_squared = 0.0;
01413         else
01414             metric_vals->relative_size_squared =
01415                 pow( VERDICT_MIN( metric_vals->area / avg_area, avg_area / metric_vals->area ), 2 );
01416     }
01417 
01418     // calculate the jacobian
01419     if( metrics_request_flag & V_QUAD_JACOBIAN )
01420     {
01421         metric_vals->jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) );
01422     }
01423 
01424     if( metrics_request_flag & ( V_QUAD_SCALED_JACOBIAN | V_QUAD_SHEAR | V_QUAD_SHEAR_AND_SIZE ) )
01425     {
01426         double scaled_jac, min_scaled_jac = VERDICT_DBL_MAX;
01427 
01428         if( lengths[0] < VERDICT_DBL_MIN || lengths[1] < VERDICT_DBL_MIN || lengths[2] < VERDICT_DBL_MIN ||
01429             lengths[3] < VERDICT_DBL_MIN )
01430         {
01431             metric_vals->scaled_jacobian = 0.0;
01432             metric_vals->shear           = 0.0;
01433         }
01434         else
01435         {
01436             scaled_jac     = areas[0] / ( lengths[0] * lengths[3] );
01437             min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
01438 
01439             scaled_jac     = areas[1] / ( lengths[1] * lengths[0] );
01440             min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
01441 
01442             scaled_jac     = areas[2] / ( lengths[2] * lengths[1] );
01443             min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
01444 
01445             scaled_jac     = areas[3] / ( lengths[3] * lengths[2] );
01446             min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );
01447 
01448             metric_vals->scaled_jacobian = min_scaled_jac;
01449 
01450             // what the heck...set shear as well
01451             if( min_scaled_jac <= VERDICT_DBL_MIN )
01452                 metric_vals->shear = 0.0;
01453             else
01454                 metric_vals->shear = min_scaled_jac;
01455         }
01456     }
01457 
01458     if( metrics_request_flag & ( V_QUAD_WARPAGE | V_QUAD_ODDY ) )
01459     {
01460         VerdictVector corner_normals[4];
01461         corner_normals[0] = edges[3] * edges[0];
01462         corner_normals[1] = edges[0] * edges[1];
01463         corner_normals[2] = edges[1] * edges[2];
01464         corner_normals[3] = edges[2] * edges[3];
01465 
01466         if( metrics_request_flag & V_QUAD_ODDY )
01467         {
01468             double oddy, max_oddy = 0.0;
01469 
01470             double diff, dot_prod;
01471 
01472             double length_squared[4];
01473             length_squared[0] = corner_normals[0].length_squared();
01474             length_squared[1] = corner_normals[1].length_squared();
01475             length_squared[2] = corner_normals[2].length_squared();
01476             length_squared[3] = corner_normals[3].length_squared();
01477 
01478             if( length_squared[0] < VERDICT_DBL_MIN || length_squared[1] < VERDICT_DBL_MIN ||
01479                 length_squared[2] < VERDICT_DBL_MIN || length_squared[3] < VERDICT_DBL_MIN )
01480                 metric_vals->oddy = VERDICT_DBL_MAX;
01481             else
01482             {
01483                 diff     = ( lengths[0] * lengths[0] ) - ( lengths[1] * lengths[1] );
01484                 dot_prod = edges[0] % edges[1];
01485                 oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[1] );
01486                 max_oddy = VERDICT_MAX( oddy, max_oddy );
01487 
01488                 diff     = ( lengths[1] * lengths[1] ) - ( lengths[2] * lengths[2] );
01489                 dot_prod = edges[1] % edges[2];
01490                 oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[2] );
01491                 max_oddy = VERDICT_MAX( oddy, max_oddy );
01492 
01493                 diff     = ( lengths[2] * lengths[2] ) - ( lengths[3] * lengths[3] );
01494                 dot_prod = edges[2] % edges[3];
01495                 oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[3] );
01496                 max_oddy = VERDICT_MAX( oddy, max_oddy );
01497 
01498                 diff     = ( lengths[3] * lengths[3] ) - ( lengths[0] * lengths[0] );
01499                 dot_prod = edges[3] % edges[0];
01500                 oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[0] );
01501                 max_oddy = VERDICT_MAX( oddy, max_oddy );
01502 
01503                 metric_vals->oddy = max_oddy;
01504             }
01505         }
01506 
01507         if( metrics_request_flag & V_QUAD_WARPAGE )
01508         {
01509             if( corner_normals[0].normalize() < VERDICT_DBL_MIN || corner_normals[1].normalize() < VERDICT_DBL_MIN ||
01510                 corner_normals[2].normalize() < VERDICT_DBL_MIN || corner_normals[3].normalize() < VERDICT_DBL_MIN )
01511                 metric_vals->warpage = VERDICT_DBL_MAX;
01512             else
01513             {
01514                 metric_vals->warpage =
01515                     pow( VERDICT_MIN( corner_normals[0] % corner_normals[2], corner_normals[1] % corner_normals[3] ),
01516                          3 );
01517             }
01518         }
01519     }
01520 
01521     if( metrics_request_flag & V_QUAD_STRETCH )
01522     {
01523         VerdictVector temp;
01524 
01525         temp.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
01526                   coordinates[2][2] - coordinates[0][2] );
01527         double diag02 = temp.length_squared();
01528 
01529         temp.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
01530                   coordinates[3][2] - coordinates[1][2] );
01531         double diag13 = temp.length_squared();
01532 
01533         static const double QUAD_STRETCH_FACTOR = sqrt( 2.0 );
01534 
01535         // 'diag02' is now the max diagonal of the quad
01536         diag02 = VERDICT_MAX( diag02, diag13 );
01537 
01538         if( diag02 < VERDICT_DBL_MIN )
01539             metric_vals->stretch = VERDICT_DBL_MAX;
01540         else
01541             metric_vals->stretch =
01542                 QUAD_STRETCH_FACTOR *
01543                 VERDICT_MIN( VERDICT_MIN( lengths[0], lengths[1] ), VERDICT_MIN( lengths[2], lengths[3] ) ) /
01544                 sqrt( diag02 );
01545     }
01546 
01547     if( metrics_request_flag & ( V_QUAD_CONDITION | V_QUAD_SHAPE | V_QUAD_SHAPE_AND_SIZE ) )
01548     {
01549         double lengths_squared[4];
01550         lengths_squared[0] = edges[0].length_squared();
01551         lengths_squared[1] = edges[1].length_squared();
01552         lengths_squared[2] = edges[2].length_squared();
01553         lengths_squared[3] = edges[3].length_squared();
01554 
01555         if( areas[0] < VERDICT_DBL_MIN || areas[1] < VERDICT_DBL_MIN || areas[2] < VERDICT_DBL_MIN ||
01556             areas[3] < VERDICT_DBL_MIN )
01557         {
01558             metric_vals->condition = VERDICT_DBL_MAX;
01559             metric_vals->shape     = 0.0;
01560         }
01561         else
01562         {
01563             double max_condition = 0.0, condition;
01564 
01565             condition     = ( lengths_squared[0] + lengths_squared[3] ) / areas[0];
01566             max_condition = VERDICT_MAX( max_condition, condition );
01567 
01568             condition     = ( lengths_squared[1] + lengths_squared[0] ) / areas[1];
01569             max_condition = VERDICT_MAX( max_condition, condition );
01570 
01571             condition     = ( lengths_squared[2] + lengths_squared[1] ) / areas[2];
01572             max_condition = VERDICT_MAX( max_condition, condition );
01573 
01574             condition     = ( lengths_squared[3] + lengths_squared[2] ) / areas[3];
01575             max_condition = VERDICT_MAX( max_condition, condition );
01576 
01577             metric_vals->condition = 0.5 * max_condition;
01578             metric_vals->shape     = 2 / max_condition;
01579         }
01580     }
01581 
01582     if( metrics_request_flag & V_QUAD_AREA )
01583     {
01584         if( metric_vals->area > 0 ) metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
01585         metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
01586     }
01587 
01588     if( metrics_request_flag & V_QUAD_MAX_EDGE_RATIO )
01589     {
01590         if( metric_vals->max_edge_ratio > 0 )
01591             metric_vals->max_edge_ratio = (double)VERDICT_MIN( metric_vals->max_edge_ratio, VERDICT_DBL_MAX );
01592         metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX );
01593     }
01594 
01595     if( metrics_request_flag & V_QUAD_CONDITION )
01596     {
01597         if( metric_vals->condition > 0 )
01598             metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
01599         metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
01600     }
01601 
01602     // calculate distortion
01603     if( metrics_request_flag & V_QUAD_DISTORTION )
01604     {
01605         metric_vals->distortion = v_quad_distortion( num_nodes, coordinates );
01606 
01607         if( metric_vals->distortion > 0 )
01608             metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
01609         metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
01610     }
01611 
01612     if( metrics_request_flag & V_QUAD_JACOBIAN )
01613     {
01614         if( metric_vals->jacobian > 0 )
01615             metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
01616         metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
01617     }
01618 
01619     if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
01620     {
01621         if( metric_vals->maximum_angle > 0 )
01622             metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
01623         metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
01624     }
01625 
01626     if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
01627     {
01628         if( metric_vals->minimum_angle > 0 )
01629             metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
01630         metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
01631     }
01632 
01633     if( metrics_request_flag & V_QUAD_ODDY )
01634     {
01635         if( metric_vals->oddy > 0 ) metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX );
01636         metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX );
01637     }
01638 
01639     if( metrics_request_flag & V_QUAD_RELATIVE_SIZE_SQUARED )
01640     {
01641         if( metric_vals->relative_size_squared > 0 )
01642             metric_vals->relative_size_squared =
01643                 (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
01644         metric_vals->relative_size_squared =
01645             (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
01646     }
01647 
01648     if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
01649     {
01650         if( metric_vals->scaled_jacobian > 0 )
01651             metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
01652         metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
01653     }
01654 
01655     if( metrics_request_flag & V_QUAD_SHEAR )
01656     {
01657         if( metric_vals->shear > 0 ) metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX );
01658         metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX );
01659     }
01660 
01661     // calculate shear and size
01662     // reuse values from above
01663     if( metrics_request_flag & V_QUAD_SHEAR_AND_SIZE )
01664     {
01665         metric_vals->shear_and_size = metric_vals->shear * metric_vals->relative_size_squared;
01666 
01667         if( metric_vals->shear_and_size > 0 )
01668             metric_vals->shear_and_size = (double)VERDICT_MIN( metric_vals->shear_and_size, VERDICT_DBL_MAX );
01669         metric_vals->shear_and_size = (double)VERDICT_MAX( metric_vals->shear_and_size, -VERDICT_DBL_MAX );
01670     }
01671 
01672     if( metrics_request_flag & V_QUAD_SHAPE )
01673     {
01674         if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
01675         metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
01676     }
01677 
01678     // calculate shape and size
01679     // reuse values from above
01680     if( metrics_request_flag & V_QUAD_SHAPE_AND_SIZE )
01681     {
01682         metric_vals->shape_and_size = metric_vals->shape * metric_vals->relative_size_squared;
01683 
01684         if( metric_vals->shape_and_size > 0 )
01685             metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
01686         metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
01687     }
01688 
01689     if( metrics_request_flag & V_QUAD_SKEW )
01690     {
01691         if( metric_vals->skew > 0 ) metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX );
01692         metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX );
01693     }
01694 
01695     if( metrics_request_flag & V_QUAD_STRETCH )
01696     {
01697         if( metric_vals->stretch > 0 )
01698             metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX );
01699         metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX );
01700     }
01701 
01702     if( metrics_request_flag & V_QUAD_TAPER )
01703     {
01704         if( metric_vals->taper > 0 ) metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX );
01705         metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX );
01706     }
01707 
01708     if( metrics_request_flag & V_QUAD_WARPAGE )
01709     {
01710         if( metric_vals->warpage > 0 )
01711             metric_vals->warpage = (double)VERDICT_MIN( metric_vals->warpage, VERDICT_DBL_MAX );
01712         metric_vals->warpage = (double)VERDICT_MAX( metric_vals->warpage, -VERDICT_DBL_MAX );
01713     }
01714 
01715     if( metrics_request_flag & V_QUAD_MED_ASPECT_FROBENIUS )
01716         metric_vals->med_aspect_frobenius = v_quad_med_aspect_frobenius( 4, coordinates );
01717     if( metrics_request_flag & V_QUAD_MAX_ASPECT_FROBENIUS )
01718         metric_vals->max_aspect_frobenius = v_quad_max_aspect_frobenius( 4, coordinates );
01719     if( metrics_request_flag & V_QUAD_EDGE_RATIO ) metric_vals->edge_ratio = v_quad_edge_ratio( 4, coordinates );
01720     if( metrics_request_flag & V_QUAD_ASPECT_RATIO ) metric_vals->aspect_ratio = v_quad_aspect_ratio( 4, coordinates );
01721     if( metrics_request_flag & V_QUAD_RADIUS_RATIO ) metric_vals->radius_ratio = v_quad_radius_ratio( 4, coordinates );
01722 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines