MOAB: Mesh Oriented datABase  (version 5.4.1)
Go to the documentation of this file.
00001 /*=========================================================================
00002
00003   Module:    $RCSfile: V_QuadMetric.cpp,v$
00004
00005   Copyright (c) 2006 Sandia Corporation.
00008
00009      This software is distributed WITHOUT ANY WARRANTY; without even
00010      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00012
00013 =========================================================================*/
00014
00015 /*
00016  *
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
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 {
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];
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];
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 {
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];
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];
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 /*!
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 */
00387 {
00388     static const double normal_coeff = 1. / ( 2. * sqrt( 2. ) );
00389
00390     VerdictVector edges[4];
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
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];
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];
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 /*!
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 /*!
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 /*!
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];
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
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;
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 {
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];
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];
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
00992     double rel_size  = 0;
00993
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 );
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,
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];
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 |
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
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
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
01398     {
01399         metric_vals->area = 0.25 * ( areas[0] + areas[1] + areas[2] + areas[3] );
01400     }
01401
01402     // calculate the relative size
01404     {
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
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
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 =
01543                 VERDICT_MIN( VERDICT_MIN( lengths[0], lengths[1] ), VERDICT_MIN( lengths[2], lengths[3] ) ) /
01544                 sqrt( diag02 );
01545     }
01546
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 );