Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
V_HexMetric.cpp
Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Module:    $RCSfile: V_HexMetric.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  * HexMetric.cpp contains quality calculations for hexes
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 "V_GaussIntegration.hpp"
00028 #include "verdict_defines.hpp"
00029 #include <memory.h>
00030 
00031 #if defined( __BORLANDC__ )
00032 #pragma warn - 8004 /* "assigned a value that is never used" */
00033 #endif
00034 
00035 //! the average volume of a hex
00036 double verdict_hex_size = 0;
00037 
00038 //! weights based on the average size of a hex
00039 int v_hex_get_weight( VerdictVector& v1, VerdictVector& v2, VerdictVector& v3 )
00040 {
00041 
00042     if( verdict_hex_size == 0 ) return 0;
00043 
00044     v1.set( 1, 0, 0 );
00045     v2.set( 0, 1, 0 );
00046     v3.set( 0, 0, 1 );
00047 
00048     double scale = pow( verdict_hex_size / ( v1 % ( v2 * v3 ) ), 0.33333333333 );
00049     v1 *= scale;
00050     v2 *= scale;
00051     v3 *= scale;
00052 
00053     return 1;
00054 }
00055 
00056 //! returns the average volume of a hex
00057 C_FUNC_DEF void v_set_hex_size( double size )
00058 {
00059     verdict_hex_size = size;
00060 }
00061 
00062 #define make_hex_nodes( coord, pos )                                                         \
00063     for( int mhcii = 0; mhcii < 8; mhcii++ )                                                 \
00064     {                                                                                        \
00065         ( pos )[mhcii].set( ( coord )[mhcii][0], ( coord )[mhcii][1], ( coord )[mhcii][2] ); \
00066     }
00067 
00068 #define make_edge_length_squares( edges, lengths )                  \
00069     {                                                               \
00070         for( int melii = 0; melii < 12; melii++ )                   \
00071             ( lengths )[melii] = ( edges )[melii].length_squared(); \
00072     }
00073 
00074 //! make VerdictVectors from coordinates
00075 void make_hex_edges( double coordinates[][3], VerdictVector edges[12] )
00076 {
00077     edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
00078                   coordinates[1][2] - coordinates[0][2] );
00079     edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
00080                   coordinates[2][2] - coordinates[1][2] );
00081     edges[2].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
00082                   coordinates[3][2] - coordinates[2][2] );
00083     edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
00084                   coordinates[0][2] - coordinates[3][2] );
00085     edges[4].set( coordinates[5][0] - coordinates[4][0], coordinates[5][1] - coordinates[4][1],
00086                   coordinates[5][2] - coordinates[4][2] );
00087     edges[5].set( coordinates[6][0] - coordinates[5][0], coordinates[6][1] - coordinates[5][1],
00088                   coordinates[6][2] - coordinates[5][2] );
00089     edges[6].set( coordinates[7][0] - coordinates[6][0], coordinates[7][1] - coordinates[6][1],
00090                   coordinates[7][2] - coordinates[6][2] );
00091     edges[7].set( coordinates[4][0] - coordinates[7][0], coordinates[4][1] - coordinates[7][1],
00092                   coordinates[4][2] - coordinates[7][2] );
00093     edges[8].set( coordinates[4][0] - coordinates[0][0], coordinates[4][1] - coordinates[0][1],
00094                   coordinates[4][2] - coordinates[0][2] );
00095     edges[9].set( coordinates[5][0] - coordinates[1][0], coordinates[5][1] - coordinates[1][1],
00096                   coordinates[5][2] - coordinates[1][2] );
00097     edges[10].set( coordinates[6][0] - coordinates[2][0], coordinates[6][1] - coordinates[2][1],
00098                    coordinates[6][2] - coordinates[2][2] );
00099     edges[11].set( coordinates[7][0] - coordinates[3][0], coordinates[7][1] - coordinates[3][1],
00100                    coordinates[7][2] - coordinates[3][2] );
00101 }
00102 
00103 /*!
00104   localizes hex coordinates
00105 */
00106 void localize_hex_coordinates( double coordinates[][3], VerdictVector position[8] )
00107 {
00108 
00109     int ii;
00110     for( ii = 0; ii < 8; ii++ )
00111     {
00112         position[ii].set( coordinates[ii][0], coordinates[ii][1], coordinates[ii][2] );
00113     }
00114 
00115     // ... Make centroid of element the center of coordinate system
00116     VerdictVector point_1256 = position[1];
00117     point_1256 += position[2];
00118     point_1256 += position[5];
00119     point_1256 += position[6];
00120 
00121     VerdictVector point_0374 = position[0];
00122     point_0374 += position[3];
00123     point_0374 += position[7];
00124     point_0374 += position[4];
00125 
00126     VerdictVector centroid = point_1256;
00127     centroid += point_0374;
00128     centroid /= 8.0;
00129 
00130     int i;
00131     for( i = 0; i < 8; i++ )
00132         position[i] -= centroid;
00133 
00134     // ... Rotate element such that center of side 1-2 and 0-3 define X axis
00135     double DX = point_1256.x() - point_0374.x();
00136     double DY = point_1256.y() - point_0374.y();
00137     double DZ = point_1256.z() - point_0374.z();
00138 
00139     double AMAGX = sqrt( DX * DX + DZ * DZ );
00140     double AMAGY = sqrt( DX * DX + DY * DY + DZ * DZ );
00141     double FMAGX = AMAGX == 0.0 ? 1.0 : 0.0;
00142     double FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;
00143 
00144     double CZ = DX / ( AMAGX + FMAGX ) + FMAGX;
00145     double SZ = DZ / ( AMAGX + FMAGX );
00146     double CY = AMAGX / ( AMAGY + FMAGY ) + FMAGY;
00147     double SY = DY / ( AMAGY + FMAGY );
00148 
00149     double temp;
00150 
00151     for( i = 0; i < 8; i++ )
00152     {
00153         temp = CY * CZ * position[i].x() + CY * SZ * position[i].z() + SY * position[i].y();
00154         position[i].y( -SY * CZ * position[i].x() - SY * SZ * position[i].z() + CY * position[i].y() );
00155         position[i].z( -SZ * position[i].x() + CZ * position[i].z() );
00156         position[i].x( temp );
00157     }
00158 
00159     // ... Now, rotate about Y
00160     VerdictVector delta = -position[0];
00161     delta -= position[1];
00162     delta += position[2];
00163     delta += position[3];
00164     delta -= position[4];
00165     delta -= position[5];
00166     delta += position[6];
00167     delta += position[7];
00168 
00169     DY = delta.y();
00170     DZ = delta.z();
00171 
00172     AMAGY = sqrt( DY * DY + DZ * DZ );
00173     FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;
00174 
00175     double CX = DY / ( AMAGY + FMAGY ) + FMAGY;
00176     double SX = DZ / ( AMAGY + FMAGY );
00177 
00178     for( i = 0; i < 8; i++ )
00179     {
00180         temp = CX * position[i].y() + SX * position[i].z();
00181         position[i].z( -SX * position[i].y() + CX * position[i].z() );
00182         position[i].y( temp );
00183     }
00184 }
00185 
00186 double safe_ratio3( const double numerator, const double denominator, const double max_ratio )
00187 {
00188     // this filter is essential for good running time in practice
00189     double return_value;
00190 
00191     const double filter_n = max_ratio * 1.0e-16;
00192     const double filter_d = 1.0e-16;
00193     if( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d )
00194     {
00195         return_value = numerator / denominator;
00196     }
00197     else
00198     {
00199         return_value = fabs( numerator ) / max_ratio >= fabs( denominator )
00200                            ? ( ( numerator >= 0.0 && denominator >= 0.0 ) || ( numerator < 0.0 && denominator < 0.0 )
00201                                    ? max_ratio
00202                                    : -max_ratio )
00203                            : numerator / denominator;
00204     }
00205 
00206     return return_value;
00207 }
00208 
00209 double safe_ratio( const double numerator, const double denominator )
00210 {
00211 
00212     double return_value;
00213     const double filter_n = VERDICT_DBL_MAX;
00214     const double filter_d = VERDICT_DBL_MIN;
00215     if( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d )
00216     {
00217         return_value = numerator / denominator;
00218     }
00219     else
00220     {
00221         return_value = VERDICT_DBL_MAX;
00222     }
00223 
00224     return return_value;
00225 }
00226 
00227 double condition_comp( const VerdictVector& xxi, const VerdictVector& xet, const VerdictVector& xze )
00228 {
00229     double det = xxi % ( xet * xze );
00230 
00231     if( det <= VERDICT_DBL_MIN )
00232     {
00233         return VERDICT_DBL_MAX;
00234     }
00235 
00236     double term1 = xxi % xxi + xet % xet + xze % xze;
00237     double term2 =
00238         ( ( xxi * xet ) % ( xxi * xet ) ) + ( ( xet * xze ) % ( xet * xze ) ) + ( ( xze * xxi ) % ( xze * xxi ) );
00239 
00240     return sqrt( term1 * term2 ) / det;
00241 }
00242 
00243 double oddy_comp( const VerdictVector& xxi, const VerdictVector& xet, const VerdictVector& xze )
00244 {
00245     static const double third = 1.0 / 3.0;
00246 
00247     double g11, g12, g13, g22, g23, g33, rt_g;
00248 
00249     g11  = xxi % xxi;
00250     g12  = xxi % xet;
00251     g13  = xxi % xze;
00252     g22  = xet % xet;
00253     g23  = xet % xze;
00254     g33  = xze % xze;
00255     rt_g = xxi % ( xet * xze );
00256 
00257     double oddy_metric;
00258     if( rt_g > VERDICT_DBL_MIN )
00259     {
00260         double norm_G_squared = g11 * g11 + 2.0 * g12 * g12 + 2.0 * g13 * g13 + g22 * g22 + 2.0 * g23 * g23 + g33 * g33;
00261 
00262         double norm_J_squared = g11 + g22 + g33;
00263 
00264         oddy_metric = ( norm_G_squared - third * norm_J_squared * norm_J_squared ) / pow( rt_g, 4. * third );
00265     }
00266 
00267     else
00268         oddy_metric = VERDICT_DBL_MAX;
00269 
00270     return oddy_metric;
00271 }
00272 
00273 //! calcualates edge lengths of a hex
00274 double hex_edge_length( int max_min, double coordinates[][3] )
00275 {
00276     double temp[3], edge[12];
00277     int i;
00278 
00279     // lengths^2 of edges
00280     for( i = 0; i < 3; i++ )
00281     {
00282         temp[i] = coordinates[1][i] - coordinates[0][i];
00283         temp[i] = temp[i] * temp[i];
00284     }
00285     edge[0] = sqrt( temp[0] + temp[1] + temp[2] );
00286 
00287     for( i = 0; i < 3; i++ )
00288     {
00289         temp[i] = coordinates[2][i] - coordinates[1][i];
00290         temp[i] = temp[i] * temp[i];
00291     }
00292     edge[1] = sqrt( temp[0] + temp[1] + temp[2] );
00293 
00294     for( i = 0; i < 3; i++ )
00295     {
00296         temp[i] = coordinates[3][i] - coordinates[2][i];
00297         temp[i] = temp[i] * temp[i];
00298     }
00299     edge[2] = sqrt( temp[0] + temp[1] + temp[2] );
00300 
00301     for( i = 0; i < 3; i++ )
00302     {
00303         temp[i] = coordinates[0][i] - coordinates[3][i];
00304         temp[i] = temp[i] * temp[i];
00305     }
00306     edge[3] = sqrt( temp[0] + temp[1] + temp[2] );
00307 
00308     for( i = 0; i < 3; i++ )
00309     {
00310         temp[i] = coordinates[5][i] - coordinates[4][i];
00311         temp[i] = temp[i] * temp[i];
00312     }
00313     edge[4] = sqrt( temp[0] + temp[1] + temp[2] );
00314 
00315     for( i = 0; i < 3; i++ )
00316     {
00317         temp[i] = coordinates[6][i] - coordinates[5][i];
00318         temp[i] = temp[i] * temp[i];
00319     }
00320     edge[5] = sqrt( temp[0] + temp[1] + temp[2] );
00321 
00322     for( i = 0; i < 3; i++ )
00323     {
00324         temp[i] = coordinates[7][i] - coordinates[6][i];
00325         temp[i] = temp[i] * temp[i];
00326     }
00327     edge[6] = sqrt( temp[0] + temp[1] + temp[2] );
00328 
00329     for( i = 0; i < 3; i++ )
00330     {
00331         temp[i] = coordinates[4][i] - coordinates[7][i];
00332         temp[i] = temp[i] * temp[i];
00333     }
00334     edge[7] = sqrt( temp[0] + temp[1] + temp[2] );
00335 
00336     for( i = 0; i < 3; i++ )
00337     {
00338         temp[i] = coordinates[4][i] - coordinates[0][i];
00339         temp[i] = temp[i] * temp[i];
00340     }
00341     edge[8] = sqrt( temp[0] + temp[1] + temp[2] );
00342 
00343     for( i = 0; i < 3; i++ )
00344     {
00345         temp[i] = coordinates[5][i] - coordinates[1][i];
00346         temp[i] = temp[i] * temp[i];
00347     }
00348     edge[9] = sqrt( temp[0] + temp[1] + temp[2] );
00349 
00350     for( i = 0; i < 3; i++ )
00351     {
00352         temp[i] = coordinates[6][i] - coordinates[2][i];
00353         temp[i] = temp[i] * temp[i];
00354     }
00355     edge[10] = sqrt( temp[0] + temp[1] + temp[2] );
00356 
00357     for( i = 0; i < 3; i++ )
00358     {
00359         temp[i] = coordinates[7][i] - coordinates[3][i];
00360         temp[i] = temp[i] * temp[i];
00361     }
00362     edge[11] = sqrt( temp[0] + temp[1] + temp[2] );
00363 
00364     double _edge = edge[0];
00365 
00366     if( max_min == 0 )
00367     {
00368         for( i = 1; i < 12; i++ )
00369             _edge = VERDICT_MIN( _edge, edge[i] );
00370         return (double)_edge;
00371     }
00372     else
00373     {
00374         for( i = 1; i < 12; i++ )
00375             _edge = VERDICT_MAX( _edge, edge[i] );
00376         return (double)_edge;
00377     }
00378 }
00379 
00380 double diag_length( int max_min, double coordinates[][3] )
00381 {
00382     double temp[3], diag[4];
00383     int i;
00384 
00385     // lengths^2  f diag nals
00386     for( i = 0; i < 3; i++ )
00387     {
00388         temp[i] = coordinates[6][i] - coordinates[0][i];
00389         temp[i] = temp[i] * temp[i];
00390     }
00391     diag[0] = sqrt( temp[0] + temp[1] + temp[2] );
00392 
00393     for( i = 0; i < 3; i++ )
00394     {
00395         temp[i] = coordinates[4][i] - coordinates[2][i];
00396         temp[i] = temp[i] * temp[i];
00397     }
00398     diag[1] = sqrt( temp[0] + temp[1] + temp[2] );
00399 
00400     for( i = 0; i < 3; i++ )
00401     {
00402         temp[i] = coordinates[7][i] - coordinates[1][i];
00403         temp[i] = temp[i] * temp[i];
00404     }
00405     diag[2] = sqrt( temp[0] + temp[1] + temp[2] );
00406 
00407     for( i = 0; i < 3; i++ )
00408     {
00409         temp[i] = coordinates[5][i] - coordinates[3][i];
00410         temp[i] = temp[i] * temp[i];
00411     }
00412     diag[3] = sqrt( temp[0] + temp[1] + temp[2] );
00413 
00414     double diagonal = diag[0];
00415     if( max_min == 0 )  // Return min diagonal
00416     {
00417         for( i = 1; i < 4; i++ )
00418             diagonal = VERDICT_MIN( diagonal, diag[i] );
00419         return (double)diagonal;
00420     }
00421     else  // Return max diagonal
00422     {
00423         for( i = 1; i < 4; i++ )
00424             diagonal = VERDICT_MAX( diagonal, diag[i] );
00425         return (double)diagonal;
00426     }
00427 }
00428 
00429 //! calculates efg values
00430 VerdictVector calc_hex_efg( int efg_index, VerdictVector coordinates[8] )
00431 {
00432 
00433     VerdictVector efg;
00434 
00435     switch( efg_index )
00436     {
00437 
00438         case 1:
00439             efg = coordinates[1];
00440             efg += coordinates[2];
00441             efg += coordinates[5];
00442             efg += coordinates[6];
00443             efg -= coordinates[0];
00444             efg -= coordinates[3];
00445             efg -= coordinates[4];
00446             efg -= coordinates[7];
00447             break;
00448 
00449         case 2:
00450             efg = coordinates[2];
00451             efg += coordinates[3];
00452             efg += coordinates[6];
00453             efg += coordinates[7];
00454             efg -= coordinates[0];
00455             efg -= coordinates[1];
00456             efg -= coordinates[4];
00457             efg -= coordinates[5];
00458             break;
00459 
00460         case 3:
00461             efg = coordinates[4];
00462             efg += coordinates[5];
00463             efg += coordinates[6];
00464             efg += coordinates[7];
00465             efg -= coordinates[0];
00466             efg -= coordinates[1];
00467             efg -= coordinates[2];
00468             efg -= coordinates[3];
00469             break;
00470 
00471         case 12:
00472             efg = coordinates[0];
00473             efg += coordinates[2];
00474             efg += coordinates[4];
00475             efg += coordinates[6];
00476             efg -= coordinates[1];
00477             efg -= coordinates[3];
00478             efg -= coordinates[5];
00479             efg -= coordinates[7];
00480             break;
00481 
00482         case 13:
00483             efg = coordinates[0];
00484             efg += coordinates[3];
00485             efg += coordinates[5];
00486             efg += coordinates[6];
00487             efg -= coordinates[1];
00488             efg -= coordinates[2];
00489             efg -= coordinates[4];
00490             efg -= coordinates[7];
00491             break;
00492 
00493         case 23:
00494             efg = coordinates[0];
00495             efg += coordinates[1];
00496             efg += coordinates[6];
00497             efg += coordinates[7];
00498             efg -= coordinates[2];
00499             efg -= coordinates[3];
00500             efg -= coordinates[4];
00501             efg -= coordinates[5];
00502             break;
00503 
00504         case 123:
00505             efg = coordinates[0];
00506             efg += coordinates[2];
00507             efg += coordinates[5];
00508             efg += coordinates[7];
00509             efg -= coordinates[1];
00510             efg -= coordinates[5];
00511             efg -= coordinates[6];
00512             efg -= coordinates[2];
00513             break;
00514 
00515         default:
00516             efg.set( 0, 0, 0 );
00517     }
00518 
00519     return efg;
00520 }
00521 
00522 /*!
00523    the edge ratio of a hex
00524 
00525    NB (P. Pebay 01/23/07):
00526      Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
00527      minimum edge lengths
00528 */
00529 C_FUNC_DEF double v_hex_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
00530 {
00531 
00532     VerdictVector edges[12];
00533     make_hex_edges( coordinates, edges );
00534 
00535     double a2 = edges[0].length_squared();
00536     double b2 = edges[1].length_squared();
00537     double c2 = edges[2].length_squared();
00538     double d2 = edges[3].length_squared();
00539     double e2 = edges[4].length_squared();
00540     double f2 = edges[5].length_squared();
00541     double g2 = edges[6].length_squared();
00542     double h2 = edges[7].length_squared();
00543     double i2 = edges[8].length_squared();
00544     double j2 = edges[9].length_squared();
00545     double k2 = edges[10].length_squared();
00546     double l2 = edges[11].length_squared();
00547 
00548     double mab, mcd, mef, Mab, Mcd, Mef;
00549     double mgh, mij, mkl, Mgh, Mij, Mkl;
00550 
00551     if( a2 < b2 )
00552     {
00553         mab = a2;
00554         Mab = b2;
00555     }
00556     else  // b2 <= a2
00557     {
00558         mab = b2;
00559         Mab = a2;
00560     }
00561     if( c2 < d2 )
00562     {
00563         mcd = c2;
00564         Mcd = d2;
00565     }
00566     else  // d2 <= c2
00567     {
00568         mcd = d2;
00569         Mcd = c2;
00570     }
00571     if( e2 < f2 )
00572     {
00573         mef = e2;
00574         Mef = f2;
00575     }
00576     else  // f2 <= e2
00577     {
00578         mef = f2;
00579         Mef = e2;
00580     }
00581     if( g2 < h2 )
00582     {
00583         mgh = g2;
00584         Mgh = h2;
00585     }
00586     else  // h2 <= g2
00587     {
00588         mgh = h2;
00589         Mgh = g2;
00590     }
00591     if( i2 < j2 )
00592     {
00593         mij = i2;
00594         Mij = j2;
00595     }
00596     else  // j2 <= i2
00597     {
00598         mij = j2;
00599         Mij = i2;
00600     }
00601     if( k2 < l2 )
00602     {
00603         mkl = k2;
00604         Mkl = l2;
00605     }
00606     else  // l2 <= k2
00607     {
00608         mkl = l2;
00609         Mkl = k2;
00610     }
00611 
00612     double m2;
00613     m2 = mab < mcd ? mab : mcd;
00614     m2 = m2 < mef ? m2 : mef;
00615     m2 = m2 < mgh ? m2 : mgh;
00616     m2 = m2 < mij ? m2 : mij;
00617     m2 = m2 < mkl ? m2 : mkl;
00618 
00619     if( m2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
00620 
00621     double M2;
00622     M2 = Mab > Mcd ? Mab : Mcd;
00623     M2 = M2 > Mef ? M2 : Mef;
00624     M2 = M2 > Mgh ? M2 : Mgh;
00625     M2 = M2 > Mij ? M2 : Mij;
00626     M2 = M2 > Mkl ? M2 : Mkl;
00627     m2 = m2 < mef ? m2 : mef;
00628 
00629     M2 = Mab > Mcd ? Mab : Mcd;
00630     M2 = M2 > Mef ? M2 : Mef;
00631 
00632     double edge_ratio = sqrt( M2 / m2 );
00633 
00634     if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
00635     return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
00636 }
00637 
00638 /*!
00639   max edge ratio of a hex
00640 
00641   Maximum edge length ratio at hex center
00642 */
00643 C_FUNC_DEF double v_hex_max_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
00644 {
00645     double aspect;
00646     VerdictVector node_pos[8];
00647     make_hex_nodes( coordinates, node_pos );
00648 
00649     double aspect_12, aspect_13, aspect_23;
00650 
00651     VerdictVector efg1 = calc_hex_efg( 1, node_pos );
00652     VerdictVector efg2 = calc_hex_efg( 2, node_pos );
00653     VerdictVector efg3 = calc_hex_efg( 3, node_pos );
00654 
00655     double mag_efg1 = efg1.length();
00656     double mag_efg2 = efg2.length();
00657     double mag_efg3 = efg3.length();
00658 
00659     aspect_12 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg2 ), VERDICT_MIN( mag_efg1, mag_efg2 ) );
00660     aspect_13 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg3 ), VERDICT_MIN( mag_efg1, mag_efg3 ) );
00661     aspect_23 = safe_ratio( VERDICT_MAX( mag_efg2, mag_efg3 ), VERDICT_MIN( mag_efg2, mag_efg3 ) );
00662 
00663     aspect = VERDICT_MAX( aspect_12, VERDICT_MAX( aspect_13, aspect_23 ) );
00664 
00665     if( aspect > 0 ) return (double)VERDICT_MIN( aspect, VERDICT_DBL_MAX );
00666     return (double)VERDICT_MAX( aspect, -VERDICT_DBL_MAX );
00667 }
00668 
00669 /*!
00670   skew of a hex
00671 
00672   Maximum ||cosA|| where A is the angle between edges at hex center.
00673 */
00674 C_FUNC_DEF double v_hex_skew( int /*num_nodes*/, double coordinates[][3] )
00675 {
00676     VerdictVector node_pos[8];
00677     make_hex_nodes( coordinates, node_pos );
00678 
00679     double skew_1, skew_2, skew_3;
00680 
00681     VerdictVector efg1 = calc_hex_efg( 1, node_pos );
00682     VerdictVector efg2 = calc_hex_efg( 2, node_pos );
00683     VerdictVector efg3 = calc_hex_efg( 3, node_pos );
00684 
00685     if( efg1.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
00686     if( efg2.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
00687     if( efg3.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
00688 
00689     skew_1 = fabs( efg1 % efg2 );
00690     skew_2 = fabs( efg1 % efg3 );
00691     skew_3 = fabs( efg2 % efg3 );
00692 
00693     double skew = ( VERDICT_MAX( skew_1, VERDICT_MAX( skew_2, skew_3 ) ) );
00694 
00695     if( skew > 0 ) return (double)VERDICT_MIN( skew, VERDICT_DBL_MAX );
00696     return (double)VERDICT_MAX( skew, -VERDICT_DBL_MAX );
00697 }
00698 
00699 /*!
00700   taper of a hex
00701 
00702   Maximum ratio of lengths derived from opposite edges.
00703 */
00704 C_FUNC_DEF double v_hex_taper( int /*num_nodes*/, double coordinates[][3] )
00705 {
00706     VerdictVector node_pos[8];
00707     make_hex_nodes( coordinates, node_pos );
00708 
00709     VerdictVector efg1 = calc_hex_efg( 1, node_pos );
00710     VerdictVector efg2 = calc_hex_efg( 2, node_pos );
00711     VerdictVector efg3 = calc_hex_efg( 3, node_pos );
00712 
00713     VerdictVector efg12 = calc_hex_efg( 12, node_pos );
00714     VerdictVector efg13 = calc_hex_efg( 13, node_pos );
00715     VerdictVector efg23 = calc_hex_efg( 23, node_pos );
00716 
00717     double taper_1 = fabs( safe_ratio( efg12.length(), VERDICT_MIN( efg1.length(), efg2.length() ) ) );
00718     double taper_2 = fabs( safe_ratio( efg13.length(), VERDICT_MIN( efg1.length(), efg3.length() ) ) );
00719     double taper_3 = fabs( safe_ratio( efg23.length(), VERDICT_MIN( efg2.length(), efg3.length() ) ) );
00720 
00721     double taper = (double)VERDICT_MAX( taper_1, VERDICT_MAX( taper_2, taper_3 ) );
00722 
00723     if( taper > 0 ) return (double)VERDICT_MIN( taper, VERDICT_DBL_MAX );
00724     return (double)VERDICT_MAX( taper, -VERDICT_DBL_MAX );
00725 }
00726 
00727 /*!
00728   volume of a hex
00729 
00730   Jacobian at hex center
00731 */
00732 C_FUNC_DEF double v_hex_volume( int /*num_nodes*/, double coordinates[][3] )
00733 {
00734     VerdictVector node_pos[8];
00735     make_hex_nodes( coordinates, node_pos );
00736 
00737     VerdictVector efg1 = calc_hex_efg( 1, node_pos );
00738     VerdictVector efg2 = calc_hex_efg( 2, node_pos );
00739     VerdictVector efg3 = calc_hex_efg( 3, node_pos );
00740 
00741     double volume;
00742     volume = (double)( efg1 % ( efg2 * efg3 ) ) / 64.0;
00743 
00744     if( volume > 0 ) return (double)VERDICT_MIN( volume, VERDICT_DBL_MAX );
00745     return (double)VERDICT_MAX( volume, -VERDICT_DBL_MAX );
00746 }
00747 
00748 /*!
00749   stretch of a hex
00750 
00751   sqrt(3) * minimum edge length / maximum diagonal length
00752 */
00753 C_FUNC_DEF double v_hex_stretch( int /*num_nodes*/, double coordinates[][3] )
00754 {
00755     static const double HEX_STRETCH_SCALE_FACTOR = sqrt( 3.0 );
00756 
00757     double min_edge = hex_edge_length( 0, coordinates );
00758     double max_diag = diag_length( 1, coordinates );
00759 
00760     double stretch = HEX_STRETCH_SCALE_FACTOR * safe_ratio( min_edge, max_diag );
00761 
00762     if( stretch > 0 ) return (double)VERDICT_MIN( stretch, VERDICT_DBL_MAX );
00763     return (double)VERDICT_MAX( stretch, -VERDICT_DBL_MAX );
00764 }
00765 
00766 /*!
00767   diagonal ratio of a hex
00768 
00769   Minimum diagonal length / maximum diagonal length
00770 */
00771 C_FUNC_DEF double v_hex_diagonal( int /*num_nodes*/, double coordinates[][3] )
00772 {
00773 
00774     double min_diag = diag_length( 0, coordinates );
00775     double max_diag = diag_length( 1, coordinates );
00776 
00777     double diagonal = safe_ratio( min_diag, max_diag );
00778 
00779     if( diagonal > 0 ) return (double)VERDICT_MIN( diagonal, VERDICT_DBL_MAX );
00780     return (double)VERDICT_MAX( diagonal, -VERDICT_DBL_MAX );
00781 }
00782 
00783 #define SQR( x ) ( ( x ) * ( x ) )
00784 
00785 /*!
00786   dimension of a hex
00787 
00788   Pronto-specific characteristic length for stable time step calculation.
00789   Char_length = Volume / 2 grad Volume
00790 */
00791 C_FUNC_DEF double v_hex_dimension( int /*num_nodes*/, double coordinates[][3] )
00792 {
00793 
00794     double gradop[9][4];
00795 
00796     double x1 = coordinates[0][0];
00797     double x2 = coordinates[1][0];
00798     double x3 = coordinates[2][0];
00799     double x4 = coordinates[3][0];
00800     double x5 = coordinates[4][0];
00801     double x6 = coordinates[5][0];
00802     double x7 = coordinates[6][0];
00803     double x8 = coordinates[7][0];
00804 
00805     double y1 = coordinates[0][1];
00806     double y2 = coordinates[1][1];
00807     double y3 = coordinates[2][1];
00808     double y4 = coordinates[3][1];
00809     double y5 = coordinates[4][1];
00810     double y6 = coordinates[5][1];
00811     double y7 = coordinates[6][1];
00812     double y8 = coordinates[7][1];
00813 
00814     double z1 = coordinates[0][2];
00815     double z2 = coordinates[1][2];
00816     double z3 = coordinates[2][2];
00817     double z4 = coordinates[3][2];
00818     double z5 = coordinates[4][2];
00819     double z6 = coordinates[5][2];
00820     double z7 = coordinates[6][2];
00821     double z8 = coordinates[7][2];
00822 
00823     double z24 = z2 - z4;
00824     double z52 = z5 - z2;
00825     double z45 = z4 - z5;
00826     gradop[1][1] =
00827         ( y2 * ( z6 - z3 - z45 ) + y3 * z24 + y4 * ( z3 - z8 - z52 ) + y5 * ( z8 - z6 - z24 ) + y6 * z52 + y8 * z45 ) /
00828         12.0;
00829 
00830     double z31 = z3 - z1;
00831     double z63 = z6 - z3;
00832     double z16 = z1 - z6;
00833     gradop[2][1] =
00834         ( y3 * ( z7 - z4 - z16 ) + y4 * z31 + y1 * ( z4 - z5 - z63 ) + y6 * ( z5 - z7 - z31 ) + y7 * z63 + y5 * z16 ) /
00835         12.0;
00836 
00837     double z42 = z4 - z2;
00838     double z74 = z7 - z4;
00839     double z27 = z2 - z7;
00840     gradop[3][1] =
00841         ( y4 * ( z8 - z1 - z27 ) + y1 * z42 + y2 * ( z1 - z6 - z74 ) + y7 * ( z6 - z8 - z42 ) + y8 * z74 + y6 * z27 ) /
00842         12.0;
00843 
00844     double z13 = z1 - z3;
00845     double z81 = z8 - z1;
00846     double z38 = z3 - z8;
00847     gradop[4][1] =
00848         ( y1 * ( z5 - z2 - z38 ) + y2 * z13 + y3 * ( z2 - z7 - z81 ) + y8 * ( z7 - z5 - z13 ) + y5 * z81 + y7 * z38 ) /
00849         12.0;
00850 
00851     double z86 = z8 - z6;
00852     double z18 = z1 - z8;
00853     double z61 = z6 - z1;
00854     gradop[5][1] =
00855         ( y8 * ( z4 - z7 - z61 ) + y7 * z86 + y6 * ( z7 - z2 - z18 ) + y1 * ( z2 - z4 - z86 ) + y4 * z18 + y2 * z61 ) /
00856         12.0;
00857 
00858     double z57 = z5 - z7;
00859     double z25 = z2 - z5;
00860     double z72 = z7 - z2;
00861     gradop[6][1] =
00862         ( y5 * ( z1 - z8 - z72 ) + y8 * z57 + y7 * ( z8 - z3 - z25 ) + y2 * ( z3 - z1 - z57 ) + y1 * z25 + y3 * z72 ) /
00863         12.0;
00864 
00865     double z68 = z6 - z8;
00866     double z36 = z3 - z6;
00867     double z83 = z8 - z3;
00868     gradop[7][1] =
00869         ( y6 * ( z2 - z5 - z83 ) + y5 * z68 + y8 * ( z5 - z4 - z36 ) + y3 * ( z4 - z2 - z68 ) + y2 * z36 + y4 * z83 ) /
00870         12.0;
00871 
00872     double z75 = z7 - z5;
00873     double z47 = z4 - z7;
00874     double z54 = z5 - z4;
00875     gradop[8][1] =
00876         ( y7 * ( z3 - z6 - z54 ) + y6 * z75 + y5 * ( z6 - z1 - z47 ) + y4 * ( z1 - z3 - z75 ) + y3 * z47 + y1 * z54 ) /
00877         12.0;
00878 
00879     double x24 = x2 - x4;
00880     double x52 = x5 - x2;
00881     double x45 = x4 - x5;
00882     gradop[1][2] =
00883         ( z2 * ( x6 - x3 - x45 ) + z3 * x24 + z4 * ( x3 - x8 - x52 ) + z5 * ( x8 - x6 - x24 ) + z6 * x52 + z8 * x45 ) /
00884         12.0;
00885 
00886     double x31 = x3 - x1;
00887     double x63 = x6 - x3;
00888     double x16 = x1 - x6;
00889     gradop[2][2] =
00890         ( z3 * ( x7 - x4 - x16 ) + z4 * x31 + z1 * ( x4 - x5 - x63 ) + z6 * ( x5 - x7 - x31 ) + z7 * x63 + z5 * x16 ) /
00891         12.0;
00892 
00893     double x42 = x4 - x2;
00894     double x74 = x7 - x4;
00895     double x27 = x2 - x7;
00896     gradop[3][2] =
00897         ( z4 * ( x8 - x1 - x27 ) + z1 * x42 + z2 * ( x1 - x6 - x74 ) + z7 * ( x6 - x8 - x42 ) + z8 * x74 + z6 * x27 ) /
00898         12.0;
00899 
00900     double x13 = x1 - x3;
00901     double x81 = x8 - x1;
00902     double x38 = x3 - x8;
00903     gradop[4][2] =
00904         ( z1 * ( x5 - x2 - x38 ) + z2 * x13 + z3 * ( x2 - x7 - x81 ) + z8 * ( x7 - x5 - x13 ) + z5 * x81 + z7 * x38 ) /
00905         12.0;
00906 
00907     double x86 = x8 - x6;
00908     double x18 = x1 - x8;
00909     double x61 = x6 - x1;
00910     gradop[5][2] =
00911         ( z8 * ( x4 - x7 - x61 ) + z7 * x86 + z6 * ( x7 - x2 - x18 ) + z1 * ( x2 - x4 - x86 ) + z4 * x18 + z2 * x61 ) /
00912         12.0;
00913 
00914     double x57 = x5 - x7;
00915     double x25 = x2 - x5;
00916     double x72 = x7 - x2;
00917     gradop[6][2] =
00918         ( z5 * ( x1 - x8 - x72 ) + z8 * x57 + z7 * ( x8 - x3 - x25 ) + z2 * ( x3 - x1 - x57 ) + z1 * x25 + z3 * x72 ) /
00919         12.0;
00920 
00921     double x68 = x6 - x8;
00922     double x36 = x3 - x6;
00923     double x83 = x8 - x3;
00924     gradop[7][2] =
00925         ( z6 * ( x2 - x5 - x83 ) + z5 * x68 + z8 * ( x5 - x4 - x36 ) + z3 * ( x4 - x2 - x68 ) + z2 * x36 + z4 * x83 ) /
00926         12.0;
00927 
00928     double x75 = x7 - x5;
00929     double x47 = x4 - x7;
00930     double x54 = x5 - x4;
00931     gradop[8][2] =
00932         ( z7 * ( x3 - x6 - x54 ) + z6 * x75 + z5 * ( x6 - x1 - x47 ) + z4 * ( x1 - x3 - x75 ) + z3 * x47 + z1 * x54 ) /
00933         12.0;
00934 
00935     double y24 = y2 - y4;
00936     double y52 = y5 - y2;
00937     double y45 = y4 - y5;
00938     gradop[1][3] =
00939         ( x2 * ( y6 - y3 - y45 ) + x3 * y24 + x4 * ( y3 - y8 - y52 ) + x5 * ( y8 - y6 - y24 ) + x6 * y52 + x8 * y45 ) /
00940         12.0;
00941 
00942     double y31 = y3 - y1;
00943     double y63 = y6 - y3;
00944     double y16 = y1 - y6;
00945     gradop[2][3] =
00946         ( x3 * ( y7 - y4 - y16 ) + x4 * y31 + x1 * ( y4 - y5 - y63 ) + x6 * ( y5 - y7 - y31 ) + x7 * y63 + x5 * y16 ) /
00947         12.0;
00948 
00949     double y42 = y4 - y2;
00950     double y74 = y7 - y4;
00951     double y27 = y2 - y7;
00952     gradop[3][3] =
00953         ( x4 * ( y8 - y1 - y27 ) + x1 * y42 + x2 * ( y1 - y6 - y74 ) + x7 * ( y6 - y8 - y42 ) + x8 * y74 + x6 * y27 ) /
00954         12.0;
00955 
00956     double y13 = y1 - y3;
00957     double y81 = y8 - y1;
00958     double y38 = y3 - y8;
00959     gradop[4][3] =
00960         ( x1 * ( y5 - y2 - y38 ) + x2 * y13 + x3 * ( y2 - y7 - y81 ) + x8 * ( y7 - y5 - y13 ) + x5 * y81 + x7 * y38 ) /
00961         12.0;
00962 
00963     double y86 = y8 - y6;
00964     double y18 = y1 - y8;
00965     double y61 = y6 - y1;
00966     gradop[5][3] =
00967         ( x8 * ( y4 - y7 - y61 ) + x7 * y86 + x6 * ( y7 - y2 - y18 ) + x1 * ( y2 - y4 - y86 ) + x4 * y18 + x2 * y61 ) /
00968         12.0;
00969 
00970     double y57 = y5 - y7;
00971     double y25 = y2 - y5;
00972     double y72 = y7 - y2;
00973     gradop[6][3] =
00974         ( x5 * ( y1 - y8 - y72 ) + x8 * y57 + x7 * ( y8 - y3 - y25 ) + x2 * ( y3 - y1 - y57 ) + x1 * y25 + x3 * y72 ) /
00975         12.0;
00976 
00977     double y68 = y6 - y8;
00978     double y36 = y3 - y6;
00979     double y83 = y8 - y3;
00980     gradop[7][3] =
00981         ( x6 * ( y2 - y5 - y83 ) + x5 * y68 + x8 * ( y5 - y4 - y36 ) + x3 * ( y4 - y2 - y68 ) + x2 * y36 + x4 * y83 ) /
00982         12.0;
00983 
00984     double y75 = y7 - y5;
00985     double y47 = y4 - y7;
00986     double y54 = y5 - y4;
00987     gradop[8][3] =
00988         ( x7 * ( y3 - y6 - y54 ) + x6 * y75 + x5 * ( y6 - y1 - y47 ) + x4 * ( y1 - y3 - y75 ) + x3 * y47 + x1 * y54 ) /
00989         12.0;
00990 
00991     //     calculate element volume and characteristic element aspect ratio
00992     //     (used in time step and hourglass control) -
00993 
00994     double volume = coordinates[0][0] * gradop[1][1] + coordinates[1][0] * gradop[2][1] +
00995                     coordinates[2][0] * gradop[3][1] + coordinates[3][0] * gradop[4][1] +
00996                     coordinates[4][0] * gradop[5][1] + coordinates[5][0] * gradop[6][1] +
00997                     coordinates[6][0] * gradop[7][1] + coordinates[7][0] * gradop[8][1];
00998     double aspect =
00999         .5 * SQR( volume ) /
01000         ( SQR( gradop[1][1] ) + SQR( gradop[2][1] ) + SQR( gradop[3][1] ) + SQR( gradop[4][1] ) + SQR( gradop[5][1] ) +
01001           SQR( gradop[6][1] ) + SQR( gradop[7][1] ) + SQR( gradop[8][1] ) + SQR( gradop[1][2] ) + SQR( gradop[2][2] ) +
01002           SQR( gradop[3][2] ) + SQR( gradop[4][2] ) + SQR( gradop[5][2] ) + SQR( gradop[6][2] ) + SQR( gradop[7][2] ) +
01003           SQR( gradop[8][2] ) + SQR( gradop[1][3] ) + SQR( gradop[2][3] ) + SQR( gradop[3][3] ) + SQR( gradop[4][3] ) +
01004           SQR( gradop[5][3] ) + SQR( gradop[6][3] ) + SQR( gradop[7][3] ) + SQR( gradop[8][3] ) );
01005 
01006     return (double)sqrt( aspect );
01007 }
01008 
01009 /*!
01010   oddy of a hex
01011 
01012   General distortion measure based on left Cauchy-Green Tensor
01013 */
01014 C_FUNC_DEF double v_hex_oddy( int /*num_nodes*/, double coordinates[][3] )
01015 {
01016 
01017     double oddy = 0.0, current_oddy;
01018     VerdictVector xxi, xet, xze;
01019 
01020     VerdictVector node_pos[8];
01021     make_hex_nodes( coordinates, node_pos );
01022 
01023     xxi = calc_hex_efg( 1, node_pos );
01024     xet = calc_hex_efg( 2, node_pos );
01025     xze = calc_hex_efg( 3, node_pos );
01026 
01027     current_oddy = oddy_comp( xxi, xet, xze );
01028     if( current_oddy > oddy )
01029     {
01030         oddy = current_oddy;
01031     }
01032 
01033     xxi.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
01034              coordinates[1][2] - coordinates[0][2] );
01035 
01036     xet.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
01037              coordinates[3][2] - coordinates[0][2] );
01038 
01039     xze.set( coordinates[4][0] - coordinates[0][0], coordinates[4][1] - coordinates[0][1],
01040              coordinates[4][2] - coordinates[0][2] );
01041 
01042     current_oddy = oddy_comp( xxi, xet, xze );
01043     if( current_oddy > oddy )
01044     {
01045         oddy = current_oddy;
01046     }
01047 
01048     xxi.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
01049              coordinates[2][2] - coordinates[1][2] );
01050 
01051     xet.set( coordinates[0][0] - coordinates[1][0], coordinates[0][1] - coordinates[1][1],
01052              coordinates[0][2] - coordinates[1][2] );
01053 
01054     xze.set( coordinates[5][0] - coordinates[1][0], coordinates[5][1] - coordinates[1][1],
01055              coordinates[5][2] - coordinates[1][2] );
01056 
01057     current_oddy = oddy_comp( xxi, xet, xze );
01058     if( current_oddy > oddy )
01059     {
01060         oddy = current_oddy;
01061     }
01062 
01063     xxi.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
01064              coordinates[3][2] - coordinates[2][2] );
01065 
01066     xet.set( coordinates[1][0] - coordinates[2][0], coordinates[1][1] - coordinates[2][1],
01067              coordinates[1][2] - coordinates[2][2] );
01068 
01069     xze.set( coordinates[6][0] - coordinates[2][0], coordinates[6][1] - coordinates[2][1],
01070              coordinates[6][2] - coordinates[2][2] );
01071 
01072     current_oddy = oddy_comp( xxi, xet, xze );
01073     if( current_oddy > oddy )
01074     {
01075         oddy = current_oddy;
01076     }
01077 
01078     xxi.set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
01079              coordinates[0][2] - coordinates[3][2] );
01080 
01081     xet.set( coordinates[2][0] - coordinates[3][0], coordinates[2][1] - coordinates[3][1],
01082              coordinates[2][2] - coordinates[3][2] );
01083 
01084     xze.set( coordinates[7][0] - coordinates[3][0], coordinates[7][1] - coordinates[3][1],
01085              coordinates[7][2] - coordinates[3][2] );
01086 
01087     current_oddy = oddy_comp( xxi, xet, xze );
01088     if( current_oddy > oddy )
01089     {
01090         oddy = current_oddy;
01091     }
01092 
01093     xxi.set( coordinates[7][0] - coordinates[4][0], coordinates[7][1] - coordinates[4][1],
01094              coordinates[7][2] - coordinates[4][2] );
01095 
01096     xet.set( coordinates[5][0] - coordinates[4][0], coordinates[5][1] - coordinates[4][1],
01097              coordinates[5][2] - coordinates[4][2] );
01098 
01099     xze.set( coordinates[0][0] - coordinates[4][0], coordinates[0][1] - coordinates[4][1],
01100              coordinates[0][2] - coordinates[4][2] );
01101 
01102     current_oddy = oddy_comp( xxi, xet, xze );
01103     if( current_oddy > oddy )
01104     {
01105         oddy = current_oddy;
01106     }
01107 
01108     xxi.set( coordinates[4][0] - coordinates[5][0], coordinates[4][1] - coordinates[5][1],
01109              coordinates[4][2] - coordinates[5][2] );
01110 
01111     xet.set( coordinates[6][0] - coordinates[5][0], coordinates[6][1] - coordinates[5][1],
01112              coordinates[6][2] - coordinates[5][2] );
01113 
01114     xze.set( coordinates[1][0] - coordinates[5][0], coordinates[1][1] - coordinates[5][1],
01115              coordinates[1][2] - coordinates[5][2] );
01116 
01117     current_oddy = oddy_comp( xxi, xet, xze );
01118     if( current_oddy > oddy )
01119     {
01120         oddy = current_oddy;
01121     }
01122 
01123     xxi.set( coordinates[5][0] - coordinates[6][0], coordinates[5][1] - coordinates[6][1],
01124              coordinates[5][2] - coordinates[6][2] );
01125 
01126     xet.set( coordinates[7][0] - coordinates[6][0], coordinates[7][1] - coordinates[6][1],
01127              coordinates[7][2] - coordinates[6][2] );
01128 
01129     xze.set( coordinates[2][0] - coordinates[6][0], coordinates[2][1] - coordinates[6][1],
01130              coordinates[2][2] - coordinates[6][2] );
01131 
01132     current_oddy = oddy_comp( xxi, xet, xze );
01133     if( current_oddy > oddy )
01134     {
01135         oddy = current_oddy;
01136     }
01137 
01138     xxi.set( coordinates[6][0] - coordinates[7][0], coordinates[6][1] - coordinates[7][1],
01139              coordinates[6][2] - coordinates[7][2] );
01140 
01141     xet.set( coordinates[4][0] - coordinates[7][0], coordinates[4][1] - coordinates[7][1],
01142              coordinates[4][2] - coordinates[7][2] );
01143 
01144     xze.set( coordinates[3][0] - coordinates[7][0], coordinates[3][1] - coordinates[7][1],
01145              coordinates[3][2] - coordinates[7][2] );
01146 
01147     current_oddy = oddy_comp( xxi, xet, xze );
01148     if( current_oddy > oddy )
01149     {
01150         oddy = current_oddy;
01151     }
01152 
01153     if( oddy > 0 ) return (double)VERDICT_MIN( oddy, VERDICT_DBL_MAX );
01154     return (double)VERDICT_MAX( oddy, -VERDICT_DBL_MAX );
01155 }
01156 
01157 /*!
01158    the average Frobenius aspect of a hex
01159 
01160    NB (P. Pebay 01/20/07):
01161      this metric is calculated by averaging the 8 Frobenius aspects at
01162      each corner of the hex, when the reference corner is right isosceles.
01163 */
01164 C_FUNC_DEF double v_hex_med_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
01165 {
01166 
01167     VerdictVector node_pos[8];
01168     make_hex_nodes( coordinates, node_pos );
01169 
01170     double med_aspect_frobenius = 0.;
01171     VerdictVector xxi, xet, xze;
01172 
01173     // J(0,0,0):
01174 
01175     xxi = node_pos[1] - node_pos[0];
01176     xet = node_pos[3] - node_pos[0];
01177     xze = node_pos[4] - node_pos[0];
01178 
01179     med_aspect_frobenius += condition_comp( xxi, xet, xze );
01180     // J(1,0,0):
01181 
01182     xxi = node_pos[2] - node_pos[1];
01183     xet = node_pos[0] - node_pos[1];
01184     xze = node_pos[5] - node_pos[1];
01185 
01186     med_aspect_frobenius += condition_comp( xxi, xet, xze );
01187     // J(1,1,0):
01188 
01189     xxi = node_pos[3] - node_pos[2];
01190     xet = node_pos[1] - node_pos[2];
01191     xze = node_pos[6] - node_pos[2];
01192 
01193     med_aspect_frobenius += condition_comp( xxi, xet, xze );
01194     // J(0,1,0):
01195 
01196     xxi = node_pos[0] - node_pos[3];
01197     xet = node_pos[2] - node_pos[3];
01198     xze = node_pos[7] - node_pos[3];
01199 
01200     med_aspect_frobenius += condition_comp( xxi, xet, xze );
01201     // J(0,0,1):
01202 
01203     xxi = node_pos[7] - node_pos[4];
01204     xet = node_pos[5] - node_pos[4];
01205     xze = node_pos[0] - node_pos[4];
01206 
01207     med_aspect_frobenius += condition_comp( xxi, xet, xze );
01208     // J(1,0,1):
01209 
01210     xxi = node_pos[4] - node_pos[5];
01211     xet = node_pos[6] - node_pos[5];
01212     xze = node_pos[1] - node_pos[5];
01213 
01214     med_aspect_frobenius += condition_comp( xxi, xet, xze );
01215     // J(1,1,1):
01216 
01217     xxi = node_pos[5] - node_pos[6];
01218     xet = node_pos[7] - node_pos[6];
01219     xze = node_pos[2] - node_pos[6];
01220 
01221     med_aspect_frobenius += condition_comp( xxi, xet, xze );
01222     // J(1,1,1):
01223 
01224     xxi = node_pos[6] - node_pos[7];
01225     xet = node_pos[4] - node_pos[7];
01226     xze = node_pos[3] - node_pos[7];
01227 
01228     med_aspect_frobenius += condition_comp( xxi, xet, xze );
01229     med_aspect_frobenius /= 24.;
01230 
01231     if( med_aspect_frobenius > 0 ) return (double)VERDICT_MIN( med_aspect_frobenius, VERDICT_DBL_MAX );
01232     return (double)VERDICT_MAX( med_aspect_frobenius, -VERDICT_DBL_MAX );
01233 }
01234 
01235 /*!
01236   maximum Frobenius condition number of a hex
01237 
01238   Maximum Frobenius condition number of the Jacobian matrix at 8 corners
01239    NB (P. Pebay 01/25/07):
01240      this metric is calculated by taking the maximum of the 8 Frobenius aspects at
01241      each corner of the hex, when the reference corner is right isosceles.
01242 */
01243 C_FUNC_DEF double v_hex_max_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
01244 {
01245 
01246     VerdictVector node_pos[8];
01247     make_hex_nodes( coordinates, node_pos );
01248 
01249     double condition = 0.0, current_condition;
01250     VerdictVector xxi, xet, xze;
01251 
01252     xxi = calc_hex_efg( 1, node_pos );
01253     xet = calc_hex_efg( 2, node_pos );
01254     xze = calc_hex_efg( 3, node_pos );
01255 
01256     current_condition = condition_comp( xxi, xet, xze );
01257     if( current_condition > condition )
01258     {
01259         condition = current_condition;
01260     }
01261 
01262     // J(0,0,0):
01263 
01264     xxi = node_pos[1] - node_pos[0];
01265     xet = node_pos[3] - node_pos[0];
01266     xze = node_pos[4] - node_pos[0];
01267 
01268     current_condition = condition_comp( xxi, xet, xze );
01269     if( current_condition > condition )
01270     {
01271         condition = current_condition;
01272     }
01273 
01274     // J(1,0,0):
01275 
01276     xxi = node_pos[2] - node_pos[1];
01277     xet = node_pos[0] - node_pos[1];
01278     xze = node_pos[5] - node_pos[1];
01279 
01280     current_condition = condition_comp( xxi, xet, xze );
01281     if( current_condition > condition )
01282     {
01283         condition = current_condition;
01284     }
01285 
01286     // J(1,1,0):
01287 
01288     xxi = node_pos[3] - node_pos[2];
01289     xet = node_pos[1] - node_pos[2];
01290     xze = node_pos[6] - node_pos[2];
01291 
01292     current_condition = condition_comp( xxi, xet, xze );
01293     if( current_condition > condition )
01294     {
01295         condition = current_condition;
01296     }
01297 
01298     // J(0,1,0):
01299 
01300     xxi = node_pos[0] - node_pos[3];
01301     xet = node_pos[2] - node_pos[3];
01302     xze = node_pos[7] - node_pos[3];
01303 
01304     current_condition = condition_comp( xxi, xet, xze );
01305     if( current_condition > condition )
01306     {
01307         condition = current_condition;
01308     }
01309 
01310     // J(0,0,1):
01311 
01312     xxi = node_pos[7] - node_pos[4];
01313     xet = node_pos[5] - node_pos[4];
01314     xze = node_pos[0] - node_pos[4];
01315 
01316     current_condition = condition_comp( xxi, xet, xze );
01317     if( current_condition > condition )
01318     {
01319         condition = current_condition;
01320     }
01321 
01322     // J(1,0,1):
01323 
01324     xxi = node_pos[4] - node_pos[5];
01325     xet = node_pos[6] - node_pos[5];
01326     xze = node_pos[1] - node_pos[5];
01327 
01328     current_condition = condition_comp( xxi, xet, xze );
01329     if( current_condition > condition )
01330     {
01331         condition = current_condition;
01332     }
01333 
01334     // J(1,1,1):
01335 
01336     xxi = node_pos[5] - node_pos[6];
01337     xet = node_pos[7] - node_pos[6];
01338     xze = node_pos[2] - node_pos[6];
01339 
01340     current_condition = condition_comp( xxi, xet, xze );
01341     if( current_condition > condition )
01342     {
01343         condition = current_condition;
01344     }
01345 
01346     // J(1,1,1):
01347 
01348     xxi = node_pos[6] - node_pos[7];
01349     xet = node_pos[4] - node_pos[7];
01350     xze = node_pos[3] - node_pos[7];
01351 
01352     current_condition = condition_comp( xxi, xet, xze );
01353     if( current_condition > condition )
01354     {
01355         condition = current_condition;
01356     }
01357 
01358     condition /= 3.0;
01359 
01360     if( condition > 0 ) return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
01361     return (double)VERDICT_MAX( condition, -VERDICT_DBL_MAX );
01362 }
01363 
01364 /*!
01365   The maximum Frobenius condition of a hex, a.k.a. condition
01366   NB (P. Pebay 01/25/07):
01367      this method is maintained for backwards compatibility only.
01368      It will become deprecated at some point.
01369 
01370 */
01371 C_FUNC_DEF double v_hex_condition( int /*num_nodes*/, double coordinates[][3] )
01372 {
01373 
01374     return v_hex_max_aspect_frobenius( 8, coordinates );
01375 }
01376 
01377 /*!
01378   jacobian of a hex
01379 
01380   Minimum pointwise volume of local map at 8 corners & center of hex
01381 */
01382 C_FUNC_DEF double v_hex_jacobian( int /*num_nodes*/, double coordinates[][3] )
01383 {
01384 
01385     VerdictVector node_pos[8];
01386     make_hex_nodes( coordinates, node_pos );
01387 
01388     double jacobian = VERDICT_DBL_MAX;
01389     double current_jacobian;
01390     VerdictVector xxi, xet, xze;
01391 
01392     xxi = calc_hex_efg( 1, node_pos );
01393     xet = calc_hex_efg( 2, node_pos );
01394     xze = calc_hex_efg( 3, node_pos );
01395 
01396     current_jacobian = xxi % ( xet * xze ) / 64.0;
01397     if( current_jacobian < jacobian )
01398     {
01399         jacobian = current_jacobian;
01400     }
01401 
01402     // J(0,0,0):
01403 
01404     xxi = node_pos[1] - node_pos[0];
01405     xet = node_pos[3] - node_pos[0];
01406     xze = node_pos[4] - node_pos[0];
01407 
01408     current_jacobian = xxi % ( xet * xze );
01409     if( current_jacobian < jacobian )
01410     {
01411         jacobian = current_jacobian;
01412     }
01413 
01414     // J(1,0,0):
01415 
01416     xxi = node_pos[2] - node_pos[1];
01417     xet = node_pos[0] - node_pos[1];
01418     xze = node_pos[5] - node_pos[1];
01419 
01420     current_jacobian = xxi % ( xet * xze );
01421     if( current_jacobian < jacobian )
01422     {
01423         jacobian = current_jacobian;
01424     }
01425 
01426     // J(1,1,0):
01427 
01428     xxi = node_pos[3] - node_pos[2];
01429     xet = node_pos[1] - node_pos[2];
01430     xze = node_pos[6] - node_pos[2];
01431 
01432     current_jacobian = xxi % ( xet * xze );
01433     if( current_jacobian < jacobian )
01434     {
01435         jacobian = current_jacobian;
01436     }
01437 
01438     // J(0,1,0):
01439 
01440     xxi = node_pos[0] - node_pos[3];
01441     xet = node_pos[2] - node_pos[3];
01442     xze = node_pos[7] - node_pos[3];
01443 
01444     current_jacobian = xxi % ( xet * xze );
01445     if( current_jacobian < jacobian )
01446     {
01447         jacobian = current_jacobian;
01448     }
01449 
01450     // J(0,0,1):
01451 
01452     xxi = node_pos[7] - node_pos[4];
01453     xet = node_pos[5] - node_pos[4];
01454     xze = node_pos[0] - node_pos[4];
01455 
01456     current_jacobian = xxi % ( xet * xze );
01457     if( current_jacobian < jacobian )
01458     {
01459         jacobian = current_jacobian;
01460     }
01461 
01462     // J(1,0,1):
01463 
01464     xxi = node_pos[4] - node_pos[5];
01465     xet = node_pos[6] - node_pos[5];
01466     xze = node_pos[1] - node_pos[5];
01467 
01468     current_jacobian = xxi % ( xet * xze );
01469     if( current_jacobian < jacobian )
01470     {
01471         jacobian = current_jacobian;
01472     }
01473 
01474     // J(1,1,1):
01475 
01476     xxi = node_pos[5] - node_pos[6];
01477     xet = node_pos[7] - node_pos[6];
01478     xze = node_pos[2] - node_pos[6];
01479 
01480     current_jacobian = xxi % ( xet * xze );
01481     if( current_jacobian < jacobian )
01482     {
01483         jacobian = current_jacobian;
01484     }
01485 
01486     // J(0,1,1):
01487 
01488     xxi = node_pos[6] - node_pos[7];
01489     xet = node_pos[4] - node_pos[7];
01490     xze = node_pos[3] - node_pos[7];
01491 
01492     current_jacobian = xxi % ( xet * xze );
01493     if( current_jacobian < jacobian )
01494     {
01495         jacobian = current_jacobian;
01496     }
01497 
01498     if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
01499     return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
01500 }
01501 
01502 /*!
01503   scaled jacobian of a hex
01504 
01505   Minimum Jacobian divided by the lengths of the 3 edge vectors
01506 */
01507 C_FUNC_DEF double v_hex_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
01508 {
01509 
01510     double jacobi, min_norm_jac = VERDICT_DBL_MAX;
01511     // double min_jacobi = VERDICT_DBL_MAX;
01512     double temp_norm_jac, lengths;
01513     double len1_sq, len2_sq, len3_sq;
01514     VerdictVector xxi, xet, xze;
01515 
01516     VerdictVector node_pos[8];
01517     make_hex_nodes( coordinates, node_pos );
01518 
01519     xxi = calc_hex_efg( 1, node_pos );
01520     xet = calc_hex_efg( 2, node_pos );
01521     xze = calc_hex_efg( 3, node_pos );
01522 
01523     jacobi = xxi % ( xet * xze );
01524     // if( jacobi < min_jacobi) { min_jacobi = jacobi; }
01525 
01526     len1_sq = xxi.length_squared();
01527     len2_sq = xet.length_squared();
01528     len3_sq = xze.length_squared();
01529 
01530     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
01531         return (double)VERDICT_DBL_MAX;
01532 
01533     lengths       = sqrt( len1_sq * len2_sq * len3_sq );
01534     temp_norm_jac = jacobi / lengths;
01535 
01536     if( temp_norm_jac < min_norm_jac )
01537         min_norm_jac = temp_norm_jac;
01538     else
01539         temp_norm_jac = jacobi;
01540 
01541     // J(0,0,0):
01542 
01543     xxi = node_pos[1] - node_pos[0];
01544     xet = node_pos[3] - node_pos[0];
01545     xze = node_pos[4] - node_pos[0];
01546 
01547     jacobi = xxi % ( xet * xze );
01548     // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
01549 
01550     len1_sq = xxi.length_squared();
01551     len2_sq = xet.length_squared();
01552     len3_sq = xze.length_squared();
01553 
01554     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
01555         return (double)VERDICT_DBL_MAX;
01556 
01557     lengths       = sqrt( len1_sq * len2_sq * len3_sq );
01558     temp_norm_jac = jacobi / lengths;
01559     if( temp_norm_jac < min_norm_jac )
01560         min_norm_jac = temp_norm_jac;
01561     else
01562         temp_norm_jac = jacobi;
01563 
01564     // J(1,0,0):
01565 
01566     xxi = node_pos[2] - node_pos[1];
01567     xet = node_pos[0] - node_pos[1];
01568     xze = node_pos[5] - node_pos[1];
01569 
01570     jacobi = xxi % ( xet * xze );
01571     // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
01572 
01573     len1_sq = xxi.length_squared();
01574     len2_sq = xet.length_squared();
01575     len3_sq = xze.length_squared();
01576 
01577     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
01578         return (double)VERDICT_DBL_MAX;
01579 
01580     lengths       = sqrt( len1_sq * len2_sq * len3_sq );
01581     temp_norm_jac = jacobi / lengths;
01582     if( temp_norm_jac < min_norm_jac )
01583         min_norm_jac = temp_norm_jac;
01584     else
01585         temp_norm_jac = jacobi;
01586 
01587     // J(1,1,0):
01588 
01589     xxi = node_pos[3] - node_pos[2];
01590     xet = node_pos[1] - node_pos[2];
01591     xze = node_pos[6] - node_pos[2];
01592 
01593     jacobi = xxi % ( xet * xze );
01594     // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
01595 
01596     len1_sq = xxi.length_squared();
01597     len2_sq = xet.length_squared();
01598     len3_sq = xze.length_squared();
01599 
01600     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
01601         return (double)VERDICT_DBL_MAX;
01602 
01603     lengths       = sqrt( len1_sq * len2_sq * len3_sq );
01604     temp_norm_jac = jacobi / lengths;
01605     if( temp_norm_jac < min_norm_jac )
01606         min_norm_jac = temp_norm_jac;
01607     else
01608         temp_norm_jac = jacobi;
01609 
01610     // J(0,1,0):
01611 
01612     xxi = node_pos[0] - node_pos[3];
01613     xet = node_pos[2] - node_pos[3];
01614     xze = node_pos[7] - node_pos[3];
01615 
01616     jacobi = xxi % ( xet * xze );
01617     // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
01618 
01619     len1_sq = xxi.length_squared();
01620     len2_sq = xet.length_squared();
01621     len3_sq = xze.length_squared();
01622 
01623     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
01624         return (double)VERDICT_DBL_MAX;
01625 
01626     lengths       = sqrt( len1_sq * len2_sq * len3_sq );
01627     temp_norm_jac = jacobi / lengths;
01628     if( temp_norm_jac < min_norm_jac )
01629         min_norm_jac = temp_norm_jac;
01630     else
01631         temp_norm_jac = jacobi;
01632 
01633     // J(0,0,1):
01634 
01635     xxi = node_pos[7] - node_pos[4];
01636     xet = node_pos[5] - node_pos[4];
01637     xze = node_pos[0] - node_pos[4];
01638 
01639     jacobi = xxi % ( xet * xze );
01640     // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
01641 
01642     len1_sq = xxi.length_squared();
01643     len2_sq = xet.length_squared();
01644     len3_sq = xze.length_squared();
01645 
01646     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
01647         return (double)VERDICT_DBL_MAX;
01648 
01649     lengths       = sqrt( len1_sq * len2_sq * len3_sq );
01650     temp_norm_jac = jacobi / lengths;
01651     if( temp_norm_jac < min_norm_jac )
01652         min_norm_jac = temp_norm_jac;
01653     else
01654         temp_norm_jac = jacobi;
01655 
01656     // J(1,0,1):
01657 
01658     xxi = node_pos[4] - node_pos[5];
01659     xet = node_pos[6] - node_pos[5];
01660     xze = node_pos[1] - node_pos[5];
01661 
01662     jacobi = xxi % ( xet * xze );
01663     // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
01664 
01665     len1_sq = xxi.length_squared();
01666     len2_sq = xet.length_squared();
01667     len3_sq = xze.length_squared();
01668 
01669     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
01670         return (double)VERDICT_DBL_MAX;
01671 
01672     lengths       = sqrt( len1_sq * len2_sq * len3_sq );
01673     temp_norm_jac = jacobi / lengths;
01674     if( temp_norm_jac < min_norm_jac )
01675         min_norm_jac = temp_norm_jac;
01676     else
01677         temp_norm_jac = jacobi;
01678 
01679     // J(1,1,1):
01680 
01681     xxi = node_pos[5] - node_pos[6];
01682     xet = node_pos[7] - node_pos[6];
01683     xze = node_pos[2] - node_pos[6];
01684 
01685     jacobi = xxi % ( xet * xze );
01686     // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
01687 
01688     len1_sq = xxi.length_squared();
01689     len2_sq = xet.length_squared();
01690     len3_sq = xze.length_squared();
01691 
01692     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
01693         return (double)VERDICT_DBL_MAX;
01694 
01695     lengths       = sqrt( len1_sq * len2_sq * len3_sq );
01696     temp_norm_jac = jacobi / lengths;
01697     if( temp_norm_jac < min_norm_jac )
01698         min_norm_jac = temp_norm_jac;
01699     else
01700         temp_norm_jac = jacobi;
01701 
01702     // J(0,1,1):
01703 
01704     xxi = node_pos[6] - node_pos[7];
01705     xet = node_pos[4] - node_pos[7];
01706     xze = node_pos[3] - node_pos[7];
01707 
01708     jacobi = xxi % ( xet * xze );
01709     // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }
01710 
01711     len1_sq = xxi.length_squared();
01712     len2_sq = xet.length_squared();
01713     len3_sq = xze.length_squared();
01714 
01715     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
01716         return (double)VERDICT_DBL_MAX;
01717 
01718     lengths       = sqrt( len1_sq * len2_sq * len3_sq );
01719     temp_norm_jac = jacobi / lengths;
01720     if( temp_norm_jac < min_norm_jac ) min_norm_jac = temp_norm_jac;
01721     // else
01722     // temp_norm_jac = jacobi;
01723 
01724     if( min_norm_jac > 0 ) return (double)VERDICT_MIN( min_norm_jac, VERDICT_DBL_MAX );
01725     return (double)VERDICT_MAX( min_norm_jac, -VERDICT_DBL_MAX );
01726 }
01727 
01728 /*!
01729   shear of a hex
01730 
01731   3/Condition number of Jacobian Skew matrix
01732 */
01733 C_FUNC_DEF double v_hex_shear( int /*num_nodes*/, double coordinates[][3] )
01734 {
01735 
01736     double shear;
01737     double min_shear = 1.0;
01738     VerdictVector xxi, xet, xze;
01739     double det, len1_sq, len2_sq, len3_sq, lengths;
01740 
01741     VerdictVector node_pos[8];
01742     make_hex_nodes( coordinates, node_pos );
01743 
01744     // J(0,0,0):
01745 
01746     xxi = node_pos[1] - node_pos[0];
01747     xet = node_pos[3] - node_pos[0];
01748     xze = node_pos[4] - node_pos[0];
01749 
01750     len1_sq = xxi.length_squared();
01751     len2_sq = xet.length_squared();
01752     len3_sq = xze.length_squared();
01753 
01754     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
01755 
01756     lengths = sqrt( len1_sq * len2_sq * len3_sq );
01757     det     = xxi % ( xet * xze );
01758     if( det < VERDICT_DBL_MIN )
01759     {
01760         return 0;
01761     }
01762 
01763     shear     = det / lengths;
01764     min_shear = VERDICT_MIN( shear, min_shear );
01765 
01766     // J(1,0,0):
01767 
01768     xxi = node_pos[2] - node_pos[1];
01769     xet = node_pos[0] - node_pos[1];
01770     xze = node_pos[5] - node_pos[1];
01771 
01772     len1_sq = xxi.length_squared();
01773     len2_sq = xet.length_squared();
01774     len3_sq = xze.length_squared();
01775 
01776     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
01777 
01778     lengths = sqrt( len1_sq * len2_sq * len3_sq );
01779     det     = xxi % ( xet * xze );
01780     if( det < VERDICT_DBL_MIN )
01781     {
01782         return 0;
01783     }
01784 
01785     shear     = det / lengths;
01786     min_shear = VERDICT_MIN( shear, min_shear );
01787 
01788     // J(1,1,0):
01789 
01790     xxi = node_pos[3] - node_pos[2];
01791     xet = node_pos[1] - node_pos[2];
01792     xze = node_pos[6] - node_pos[2];
01793 
01794     len1_sq = xxi.length_squared();
01795     len2_sq = xet.length_squared();
01796     len3_sq = xze.length_squared();
01797 
01798     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
01799 
01800     lengths = sqrt( len1_sq * len2_sq * len3_sq );
01801     det     = xxi % ( xet * xze );
01802     if( det < VERDICT_DBL_MIN )
01803     {
01804         return 0;
01805     }
01806 
01807     shear     = det / lengths;
01808     min_shear = VERDICT_MIN( shear, min_shear );
01809 
01810     // J(0,1,0):
01811 
01812     xxi = node_pos[0] - node_pos[3];
01813     xet = node_pos[2] - node_pos[3];
01814     xze = node_pos[7] - node_pos[3];
01815 
01816     len1_sq = xxi.length_squared();
01817     len2_sq = xet.length_squared();
01818     len3_sq = xze.length_squared();
01819 
01820     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
01821 
01822     lengths = sqrt( len1_sq * len2_sq * len3_sq );
01823     det     = xxi % ( xet * xze );
01824     if( det < VERDICT_DBL_MIN )
01825     {
01826         return 0;
01827     }
01828 
01829     shear     = det / lengths;
01830     min_shear = VERDICT_MIN( shear, min_shear );
01831 
01832     // J(0,0,1):
01833 
01834     xxi = node_pos[7] - node_pos[4];
01835     xet = node_pos[5] - node_pos[4];
01836     xze = node_pos[0] - node_pos[4];
01837 
01838     len1_sq = xxi.length_squared();
01839     len2_sq = xet.length_squared();
01840     len3_sq = xze.length_squared();
01841 
01842     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
01843 
01844     lengths = sqrt( len1_sq * len2_sq * len3_sq );
01845     det     = xxi % ( xet * xze );
01846     if( det < VERDICT_DBL_MIN )
01847     {
01848         return 0;
01849     }
01850 
01851     shear     = det / lengths;
01852     min_shear = VERDICT_MIN( shear, min_shear );
01853 
01854     // J(1,0,1):
01855 
01856     xxi = node_pos[4] - node_pos[5];
01857     xet = node_pos[6] - node_pos[5];
01858     xze = node_pos[1] - node_pos[5];
01859 
01860     len1_sq = xxi.length_squared();
01861     len2_sq = xet.length_squared();
01862     len3_sq = xze.length_squared();
01863 
01864     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
01865 
01866     lengths = sqrt( len1_sq * len2_sq * len3_sq );
01867     det     = xxi % ( xet * xze );
01868     if( det < VERDICT_DBL_MIN )
01869     {
01870         return 0;
01871     }
01872 
01873     shear     = det / lengths;
01874     min_shear = VERDICT_MIN( shear, min_shear );
01875 
01876     // J(1,1,1):
01877 
01878     xxi = node_pos[5] - node_pos[6];
01879     xet = node_pos[7] - node_pos[6];
01880     xze = node_pos[2] - node_pos[6];
01881 
01882     len1_sq = xxi.length_squared();
01883     len2_sq = xet.length_squared();
01884     len3_sq = xze.length_squared();
01885 
01886     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
01887 
01888     lengths = sqrt( len1_sq * len2_sq * len3_sq );
01889     det     = xxi % ( xet * xze );
01890     if( det < VERDICT_DBL_MIN )
01891     {
01892         return 0;
01893     }
01894 
01895     shear     = det / lengths;
01896     min_shear = VERDICT_MIN( shear, min_shear );
01897 
01898     // J(0,1,1):
01899 
01900     xxi = node_pos[6] - node_pos[7];
01901     xet = node_pos[4] - node_pos[7];
01902     xze = node_pos[3] - node_pos[7];
01903 
01904     len1_sq = xxi.length_squared();
01905     len2_sq = xet.length_squared();
01906     len3_sq = xze.length_squared();
01907 
01908     if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;
01909 
01910     lengths = sqrt( len1_sq * len2_sq * len3_sq );
01911     det     = xxi % ( xet * xze );
01912     if( det < VERDICT_DBL_MIN )
01913     {
01914         return 0;
01915     }
01916 
01917     shear     = det / lengths;
01918     min_shear = VERDICT_MIN( shear, min_shear );
01919 
01920     if( min_shear <= VERDICT_DBL_MIN ) min_shear = 0;
01921 
01922     if( min_shear > 0 ) return (double)VERDICT_MIN( min_shear, VERDICT_DBL_MAX );
01923     return (double)VERDICT_MAX( min_shear, -VERDICT_DBL_MAX );
01924 }
01925 
01926 /*!
01927   shape of a hex
01928 
01929   3/Condition number of weighted Jacobian matrix
01930 */
01931 C_FUNC_DEF double v_hex_shape( int /*num_nodes*/, double coordinates[][3] )
01932 {
01933 
01934     double det, shape;
01935     double min_shape               = 1.0;
01936     static const double two_thirds = 2.0 / 3.0;
01937 
01938     VerdictVector xxi, xet, xze;
01939 
01940     VerdictVector node_pos[8];
01941     make_hex_nodes( coordinates, node_pos );
01942 
01943     // J(0,0,0):
01944 
01945     xxi = node_pos[1] - node_pos[0];
01946     xet = node_pos[3] - node_pos[0];
01947     xze = node_pos[4] - node_pos[0];
01948 
01949     det = xxi % ( xet * xze );
01950     if( det > VERDICT_DBL_MIN )
01951         shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
01952     else
01953         return 0;
01954 
01955     if( shape < min_shape )
01956     {
01957         min_shape = shape;
01958     }
01959 
01960     // J(1,0,0):
01961 
01962     xxi = node_pos[2] - node_pos[1];
01963     xet = node_pos[0] - node_pos[1];
01964     xze = node_pos[5] - node_pos[1];
01965 
01966     det = xxi % ( xet * xze );
01967     if( det > VERDICT_DBL_MIN )
01968         shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
01969     else
01970         return 0;
01971 
01972     if( shape < min_shape )
01973     {
01974         min_shape = shape;
01975     }
01976 
01977     // J(1,1,0):
01978 
01979     xxi = node_pos[3] - node_pos[2];
01980     xet = node_pos[1] - node_pos[2];
01981     xze = node_pos[6] - node_pos[2];
01982 
01983     det = xxi % ( xet * xze );
01984     if( det > VERDICT_DBL_MIN )
01985         shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
01986     else
01987         return 0;
01988 
01989     if( shape < min_shape )
01990     {
01991         min_shape = shape;
01992     }
01993 
01994     // J(0,1,0):
01995 
01996     xxi = node_pos[0] - node_pos[3];
01997     xet = node_pos[2] - node_pos[3];
01998     xze = node_pos[7] - node_pos[3];
01999 
02000     det = xxi % ( xet * xze );
02001     if( det > VERDICT_DBL_MIN )
02002         shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
02003     else
02004         return 0;
02005 
02006     if( shape < min_shape )
02007     {
02008         min_shape = shape;
02009     }
02010 
02011     // J(0,0,1):
02012 
02013     xxi = node_pos[7] - node_pos[4];
02014     xet = node_pos[5] - node_pos[4];
02015     xze = node_pos[0] - node_pos[4];
02016 
02017     det = xxi % ( xet * xze );
02018     if( det > VERDICT_DBL_MIN )
02019         shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
02020     else
02021         return 0;
02022 
02023     if( shape < min_shape )
02024     {
02025         min_shape = shape;
02026     }
02027 
02028     // J(1,0,1):
02029 
02030     xxi = node_pos[4] - node_pos[5];
02031     xet = node_pos[6] - node_pos[5];
02032     xze = node_pos[1] - node_pos[5];
02033 
02034     det = xxi % ( xet * xze );
02035     if( det > VERDICT_DBL_MIN )
02036         shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
02037     else
02038         return 0;
02039 
02040     if( shape < min_shape )
02041     {
02042         min_shape = shape;
02043     }
02044 
02045     // J(1,1,1):
02046 
02047     xxi = node_pos[5] - node_pos[6];
02048     xet = node_pos[7] - node_pos[6];
02049     xze = node_pos[2] - node_pos[6];
02050 
02051     det = xxi % ( xet * xze );
02052     if( det > VERDICT_DBL_MIN )
02053         shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
02054     else
02055         return 0;
02056 
02057     if( shape < min_shape )
02058     {
02059         min_shape = shape;
02060     }
02061 
02062     // J(1,1,1):
02063 
02064     xxi = node_pos[6] - node_pos[7];
02065     xet = node_pos[4] - node_pos[7];
02066     xze = node_pos[3] - node_pos[7];
02067 
02068     det = xxi % ( xet * xze );
02069     if( det > VERDICT_DBL_MIN )
02070         shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
02071     else
02072         return 0;
02073 
02074     if( shape < min_shape )
02075     {
02076         min_shape = shape;
02077     }
02078 
02079     if( min_shape <= VERDICT_DBL_MIN ) min_shape = 0;
02080 
02081     if( min_shape > 0 ) return (double)VERDICT_MIN( min_shape, VERDICT_DBL_MAX );
02082     return (double)VERDICT_MAX( min_shape, -VERDICT_DBL_MAX );
02083 }
02084 
02085 /*!
02086   relative size of a hex
02087 
02088   Min( J, 1/J ), where J is determinant of weighted Jacobian matrix
02089 */
02090 C_FUNC_DEF double v_hex_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
02091 {
02092     double size = 0;
02093     double tau;
02094 
02095     VerdictVector xxi, xet, xze;
02096     double det, det_sum = 0;
02097 
02098     v_hex_get_weight( xxi, xet, xze );
02099 
02100     // This is the average relative size
02101     double detw = xxi % ( xet * xze );
02102 
02103     if( detw < VERDICT_DBL_MIN ) return 0;
02104 
02105     VerdictVector node_pos[8];
02106     make_hex_nodes( coordinates, node_pos );
02107 
02108     // J(0,0,0):
02109 
02110     xxi = node_pos[1] - node_pos[0];
02111     xet = node_pos[3] - node_pos[0];
02112     xze = node_pos[4] - node_pos[0];
02113 
02114     det = xxi % ( xet * xze );
02115     det_sum += det;
02116 
02117     // J(1,0,0):
02118 
02119     xxi = node_pos[2] - node_pos[1];
02120     xet = node_pos[0] - node_pos[1];
02121     xze = node_pos[5] - node_pos[1];
02122 
02123     det = xxi % ( xet * xze );
02124     det_sum += det;
02125 
02126     // J(0,1,0):
02127 
02128     xxi = node_pos[3] - node_pos[2];
02129     xet = node_pos[1] - node_pos[2];
02130     xze = node_pos[6] - node_pos[2];
02131 
02132     det = xxi % ( xet * xze );
02133     det_sum += det;
02134 
02135     // J(1,1,0):
02136 
02137     xxi = node_pos[0] - node_pos[3];
02138     xet = node_pos[2] - node_pos[3];
02139     xze = node_pos[7] - node_pos[3];
02140 
02141     det = xxi % ( xet * xze );
02142     det_sum += det;
02143 
02144     // J(0,1,0):
02145 
02146     xxi = node_pos[7] - node_pos[4];
02147     xet = node_pos[5] - node_pos[4];
02148     xze = node_pos[0] - node_pos[4];
02149 
02150     det = xxi % ( xet * xze );
02151     det_sum += det;
02152 
02153     // J(1,0,1):
02154 
02155     xxi = node_pos[4] - node_pos[5];
02156     xet = node_pos[6] - node_pos[5];
02157     xze = node_pos[1] - node_pos[5];
02158 
02159     det = xxi % ( xet * xze );
02160     det_sum += det;
02161 
02162     // J(1,1,1):
02163 
02164     xxi = node_pos[5] - node_pos[6];
02165     xet = node_pos[7] - node_pos[6];
02166     xze = node_pos[2] - node_pos[6];
02167 
02168     det = xxi % ( xet * xze );
02169     det_sum += det;
02170 
02171     // J(1,1,1):
02172 
02173     xxi = node_pos[6] - node_pos[7];
02174     xet = node_pos[4] - node_pos[7];
02175     xze = node_pos[3] - node_pos[7];
02176 
02177     det = xxi % ( xet * xze );
02178     det_sum += det;
02179 
02180     if( det_sum > VERDICT_DBL_MIN )
02181     {
02182         tau = det_sum / ( 8 * detw );
02183 
02184         tau = VERDICT_MIN( tau, 1.0 / tau );
02185 
02186         size = tau * tau;
02187     }
02188 
02189     if( size > 0 ) return (double)VERDICT_MIN( size, VERDICT_DBL_MAX );
02190     return (double)VERDICT_MAX( size, -VERDICT_DBL_MAX );
02191 }
02192 
02193 /*!
02194   shape and size of a hex
02195 
02196   Product of Shape and Relative Size
02197 */
02198 C_FUNC_DEF double v_hex_shape_and_size( int num_nodes, double coordinates[][3] )
02199 {
02200     double size  = v_hex_relative_size_squared( num_nodes, coordinates );
02201     double shape = v_hex_shape( num_nodes, coordinates );
02202 
02203     double shape_size = size * shape;
02204 
02205     if( shape_size > 0 ) return (double)VERDICT_MIN( shape_size, VERDICT_DBL_MAX );
02206     return (double)VERDICT_MAX( shape_size, -VERDICT_DBL_MAX );
02207 }
02208 
02209 /*!
02210   shear and size of a hex
02211 
02212   Product of Shear and Relative Size
02213 */
02214 C_FUNC_DEF double v_hex_shear_and_size( int num_nodes, double coordinates[][3] )
02215 {
02216     double size  = v_hex_relative_size_squared( num_nodes, coordinates );
02217     double shear = v_hex_shear( num_nodes, coordinates );
02218 
02219     double shear_size = shear * size;
02220 
02221     if( shear_size > 0 ) return (double)VERDICT_MIN( shear_size, VERDICT_DBL_MAX );
02222     return (double)VERDICT_MAX( shear_size, -VERDICT_DBL_MAX );
02223 }
02224 
02225 /*!
02226   distortion of a hex
02227 */
02228 C_FUNC_DEF double v_hex_distortion( int num_nodes, double coordinates[][3] )
02229 {
02230 
02231     // use 2x2 gauss points for linear hex and 3x3 for 2nd order hex
02232     int number_of_gauss_points = 0;
02233     if( num_nodes == 8 )
02234         // 2x2 quadrature rule
02235         number_of_gauss_points = 2;
02236     else if( num_nodes == 20 )
02237         // 3x3 quadrature rule
02238         number_of_gauss_points = 3;
02239 
02240     int number_dimension             = 3;
02241     int total_number_of_gauss_points = number_of_gauss_points * number_of_gauss_points * number_of_gauss_points;
02242     double distortion                = VERDICT_DBL_MAX;
02243 
02244     // maxTotalNumberGaussPoints =27, maxNumberNodes = 20
02245     // they are defined in GaussIntegration.hpp
02246     // This is used to make these arrays static.
02247     // I tried dynamically allocated arrays but the new and delete
02248     // was very expensive
02249 
02250     double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
02251     double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
02252     double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
02253     double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
02254     double weight[maxTotalNumberGaussPoints];
02255 
02256     // create an object of GaussIntegration
02257     GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dimension );
02258     GaussIntegration::calculate_shape_function_3d_hex();
02259     GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight );
02260 
02261     VerdictVector xxi, xet, xze, xin;
02262 
02263     double jacobian, minimum_jacobian;
02264     double element_volume = 0.0;
02265     minimum_jacobian      = VERDICT_DBL_MAX;
02266     // calculate element volume
02267     int ife, ja;
02268     for( ife = 0; ife < total_number_of_gauss_points; ife++ )
02269     {
02270 
02271         xxi.set( 0.0, 0.0, 0.0 );
02272         xet.set( 0.0, 0.0, 0.0 );
02273         xze.set( 0.0, 0.0, 0.0 );
02274 
02275         for( ja = 0; ja < num_nodes; ja++ )
02276         {
02277             xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
02278             xxi += dndy1[ife][ja] * xin;
02279             xet += dndy2[ife][ja] * xin;
02280             xze += dndy3[ife][ja] * xin;
02281         }
02282 
02283         jacobian = xxi % ( xet * xze );
02284         if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
02285 
02286         element_volume += weight[ife] * jacobian;
02287     }
02288 
02289     // loop through all nodes
02290     double dndy1_at_node[maxNumberNodes][maxNumberNodes];
02291     double dndy2_at_node[maxNumberNodes][maxNumberNodes];
02292     double dndy3_at_node[maxNumberNodes][maxNumberNodes];
02293 
02294     GaussIntegration::calculate_derivative_at_nodes_3d( dndy1_at_node, dndy2_at_node, dndy3_at_node );
02295     int node_id;
02296     for( node_id = 0; node_id < num_nodes; node_id++ )
02297     {
02298 
02299         xxi.set( 0.0, 0.0, 0.0 );
02300         xet.set( 0.0, 0.0, 0.0 );
02301         xze.set( 0.0, 0.0, 0.0 );
02302 
02303         for( ja = 0; ja < num_nodes; ja++ )
02304         {
02305             xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
02306             xxi += dndy1_at_node[node_id][ja] * xin;
02307             xet += dndy2_at_node[node_id][ja] * xin;
02308             xze += dndy3_at_node[node_id][ja] * xin;
02309         }
02310 
02311         jacobian = xxi % ( xet * xze );
02312         if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
02313     }
02314     distortion = minimum_jacobian / element_volume * 8.;
02315     return (double)distortion;
02316 }
02317 
02318 /*
02319 C_FUNC_DEF double hex_jac_normjac_oddy_cond( int choices[],
02320                       double coordinates[][3],
02321                       double answers[4]  )
02322 {
02323 
02324   //Define variables
02325   int i;
02326 
02327   double xxi[3], xet[3], xze[3];
02328   double norm_jacobian = 0.0, current_norm_jac = 0.0;
02329         double jacobian = 0.0, current_jacobian = 0.0;
02330   double oddy = 0.0, current_oddy = 0.0;
02331   double condition = 0.0, current_condition = 0.0;
02332 
02333 
02334         for( i=0; i<3; i++)
02335           xxi[i] = calc_hex_efg( 2, i, coordinates );
02336         for( i=0; i<3; i++)
02337           xet[i] = calc_hex_efg( 3, i, coordinates );
02338         for( i=0; i<3; i++)
02339           xze[i] = calc_hex_efg( 6, i, coordinates );
02340 
02341   // norm jacobian and jacobian
02342   if( choices[0] || choices[1] )
02343   {
02344     current_jacobian = determinant( xxi, xet, xze  );
02345     current_norm_jac = normalize_jacobian( current_jacobian,
02346             xxi, xet, xze );
02347 
02348     if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
02349 
02350     current_jacobian /= 64.0;
02351     if (current_jacobian < jacobian) { jacobian = current_jacobian; }
02352   }
02353 
02354   // oddy
02355   if( choices[2] )
02356   {
02357     current_oddy = oddy_comp( xxi, xet, xze );
02358       if ( current_oddy > oddy ) { oddy = current_oddy; }
02359   }
02360 
02361   // condition
02362   if( choices[3] )
02363   {
02364     current_condition = condition_comp( xxi, xet, xze );
02365     if ( current_condition > condition ) { condition = current_condition; }
02366   }
02367 
02368 
02369   for( i=0; i<3; i++)
02370   {
02371     xxi[i] = coordinates[1][i] - coordinates[0][i];
02372     xet[i] = coordinates[3][i] - coordinates[0][i];
02373     xze[i] = coordinates[4][i] - coordinates[0][i];
02374   }
02375 
02376   // norm jacobian and jacobian
02377   if( choices[0] || choices[1] )
02378   {
02379     current_jacobian = determinant( xxi, xet, xze  );
02380     current_norm_jac = normalize_jacobian( current_jacobian,
02381             xxi, xet, xze );
02382 
02383     if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
02384     if (current_jacobian < jacobian) { jacobian = current_jacobian; }
02385   }
02386 
02387   // oddy
02388   if( choices[2] )
02389   {
02390     current_oddy = oddy_comp( xxi, xet, xze );
02391       if ( current_oddy > oddy ) { oddy = current_oddy; }
02392   }
02393 
02394   // condition
02395   if( choices[3] )
02396   {
02397     current_condition = condition_comp( xxi, xet, xze );
02398     if ( current_condition > condition ) { condition = current_condition; }
02399   }
02400 
02401 
02402   for( i=0; i<3; i++)
02403   {
02404           xxi[i] = coordinates[1][i] - coordinates[0][i];
02405           xet[i] = coordinates[2][i] - coordinates[1][i];
02406           xze[i] = coordinates[5][i] - coordinates[1][i];
02407   }
02408 
02409   // norm jacobian and jacobian
02410   if( choices[0] || choices[1] )
02411   {
02412     current_jacobian = determinant( xxi,  xet, xze );
02413     current_norm_jac = normalize_jacobian( current_jacobian,
02414             xxi, xet, xze );
02415 
02416     if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
02417     if (current_jacobian < jacobian) { jacobian = current_jacobian; }
02418   }
02419 
02420   // oddy
02421   if( choices[2] )
02422   {
02423     current_oddy = oddy_comp( xxi, xet, xze );
02424       if ( current_oddy > oddy ) { oddy = current_oddy; }
02425   }
02426 
02427   // condition
02428   if( choices[3] )
02429   {
02430     current_condition = condition_comp( xxi, xet, xze );
02431     if ( current_condition > condition ) { condition = current_condition; }
02432   }
02433 
02434 
02435   for( i=0; i<3; i++)
02436   {
02437           xxi[i] = coordinates[2][i] - coordinates[3][i];
02438           xet[i] = coordinates[3][i] - coordinates[0][i];
02439           xze[i] = coordinates[7][i] - coordinates[3][i];
02440   }
02441 
02442   // norm jacobian and jacobian
02443   if( choices[0] || choices[1] )
02444   {
02445     current_jacobian = determinant( xxi, xet, xze );
02446     current_norm_jac = normalize_jacobian( current_jacobian,
02447             xxi, xet, xze );
02448 
02449     if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
02450     if (current_jacobian < jacobian) { jacobian = current_jacobian; }
02451   }
02452 
02453   // oddy
02454   if( choices[2] )
02455   {
02456     current_oddy = oddy_comp( xxi, xet, xze );
02457       if ( current_oddy > oddy ) { oddy = current_oddy; }
02458   }
02459 
02460   // condition
02461   if( choices[3] )
02462   {
02463     current_condition = condition_comp( xxi, xet, xze );
02464     if ( current_condition > condition ) { condition = current_condition; }
02465   }
02466 
02467 
02468   for( i=0; i<3; i++)
02469   {
02470           xxi[i] = coordinates[5][i] - coordinates[4][i];
02471           xet[i] = coordinates[7][i] - coordinates[4][i];
02472           xze[i] = coordinates[4][i] - coordinates[0][i];
02473   }
02474 
02475 
02476   // norm jacobian and jacobian
02477   if( choices[0] || choices[1] )
02478   {
02479     current_jacobian = determinant( xxi, xet, xze );
02480     current_norm_jac = normalize_jacobian( current_jacobian,
02481             xxi, xet, xze );
02482 
02483     if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
02484     if (current_jacobian < jacobian) { jacobian = current_jacobian; }
02485   }
02486 
02487   // oddy
02488   if( choices[2] )
02489   {
02490     current_oddy = oddy_comp( xxi, xet, xze );
02491       if ( current_oddy > oddy ) { oddy = current_oddy; }
02492   }
02493 
02494   // condition
02495   if( choices[3] )
02496   {
02497     current_condition = condition_comp( xxi, xet, xze );
02498     if ( current_condition > condition ) { condition = current_condition; }
02499   }
02500 
02501 
02502   for( i=0; i<3; i++)
02503   {
02504           xxi[i] = coordinates[2][i] - coordinates[3][i];
02505           xet[i] = coordinates[2][i] - coordinates[1][i];
02506           xze[i] = coordinates[6][i] - coordinates[2][i];
02507   }
02508 
02509   // norm jacobian and jacobian
02510   if( choices[0] || choices[1] )
02511   {
02512     current_jacobian = determinant( xxi, xet, xze );
02513     current_norm_jac = normalize_jacobian( current_jacobian,
02514             xxi, xet, xze );
02515 
02516     if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
02517     if (current_jacobian < jacobian) { jacobian = current_jacobian; }
02518   }
02519 
02520   // oddy
02521   if( choices[2] )
02522   {
02523     current_oddy = oddy_comp( xxi, xet, xze );
02524       if ( current_oddy > oddy ) { oddy = current_oddy; }
02525   }
02526 
02527   // condition
02528   if( choices[3] )
02529   {
02530     current_condition = condition_comp( xxi, xet, xze );
02531     if ( current_condition > condition ) { condition = current_condition; }
02532   }
02533 
02534 
02535   for( i=0; i<3; i++)
02536   {
02537           xxi[i] = coordinates[5][i] - coordinates[4][i];
02538           xet[i] = coordinates[6][i] - coordinates[5][i];
02539           xze[i] = coordinates[5][i] - coordinates[1][i];
02540   }
02541 
02542 
02543   // norm jacobian and jacobian
02544   if( choices[0] || choices[1] )
02545   {
02546     current_jacobian = determinant( xxi, xet, xze );
02547     current_norm_jac = normalize_jacobian( current_jacobian,
02548             xxi, xet, xze );
02549 
02550     if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
02551     if (current_jacobian < jacobian) { jacobian = current_jacobian; }
02552   }
02553 
02554   // oddy
02555   if( choices[2] )
02556   {
02557     current_oddy = oddy_comp( xxi, xet, xze );
02558       if ( current_oddy > oddy ) { oddy = current_oddy; }
02559   }
02560 
02561   // condition
02562   if( choices[3] )
02563   {
02564     current_condition = condition_comp( xxi, xet, xze );
02565     if ( current_condition > condition ) { condition = current_condition; }
02566   }
02567 
02568 
02569   for( i=0; i<3; i++)
02570   {
02571           xxi[i] = coordinates[6][i] - coordinates[7][i];
02572           xet[i] = coordinates[7][i] - coordinates[4][i];
02573           xze[i] = coordinates[7][i] - coordinates[3][i];
02574   }
02575 
02576 
02577   // norm jacobian and jacobian
02578   if( choices[0] || choices[1] )
02579   {
02580     current_jacobian = determinant( xxi, xet, xze );
02581     current_norm_jac = normalize_jacobian( current_jacobian,
02582             xxi, xet, xze );
02583 
02584     if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
02585     if (current_jacobian < jacobian) { jacobian = current_jacobian; }
02586   }
02587 
02588   // oddy
02589   if( choices[2] )
02590   {
02591     current_oddy = oddy_comp( xxi, xet, xze );
02592       if ( current_oddy > oddy ) { oddy = current_oddy; }
02593   }
02594 
02595   // condition
02596   if( choices[3] )
02597   {
02598     current_condition = condition_comp( xxi, xet, xze );
02599     if ( current_condition > condition ) { condition = current_condition; }
02600   }
02601 
02602 
02603   for( i=0; i<3; i++)
02604   {
02605           xxi[i] = coordinates[6][i] - coordinates[7][i];
02606           xet[i] = coordinates[6][i] - coordinates[5][i];
02607           xze[i] = coordinates[6][i] - coordinates[2][i];
02608   }
02609 
02610   // norm jacobian and jacobian
02611   if( choices[0] || choices[1] )
02612   {
02613     current_jacobian = determinant( xxi, xet, xze );
02614     current_norm_jac = normalize_jacobian( current_jacobian,
02615             xxi, xet, xze );
02616 
02617     if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; }
02618     if (current_jacobian < jacobian) { jacobian = current_jacobian; }
02619   }
02620 
02621   // oddy
02622   if( choices[2] )
02623   {
02624     current_oddy = oddy_comp( xxi, xet, xze );
02625       if ( current_oddy > oddy ) { oddy = current_oddy; }
02626   }
02627 
02628   // condition
02629   if( choices[3] )
02630   {
02631     current_condition = condition_comp( xxi, xet, xze );
02632     if ( current_condition > condition ) { condition = current_condition; }
02633 
02634     condition /= 3.0;
02635   }
02636 
02637 
02638   answers[0] = jacobian;
02639   answers[1] = norm_jacobian;
02640   answers[2] = oddy;
02641   answers[3] = condition;
02642 
02643   return 1.0;
02644 
02645 }
02646 */
02647 
02648 /*!
02649   multiple quality metrics of a hex
02650 */
02651 C_FUNC_DEF void v_hex_quality( int num_nodes,
02652                                double coordinates[][3],
02653                                unsigned int metrics_request_flag,
02654                                HexMetricVals* metric_vals )
02655 {
02656     memset( metric_vals, 0, sizeof( HexMetricVals ) );
02657 
02658     // max edge ratio, skew, taper, and volume
02659     if( metrics_request_flag & ( V_HEX_MAX_EDGE_RATIO | V_HEX_SKEW | V_HEX_TAPER ) )
02660     {
02661         VerdictVector node_pos[8];
02662         make_hex_nodes( coordinates, node_pos );
02663 
02664         VerdictVector efg1, efg2, efg3;
02665         efg1 = calc_hex_efg( 1, node_pos );
02666         efg2 = calc_hex_efg( 2, node_pos );
02667         efg3 = calc_hex_efg( 3, node_pos );
02668 
02669         if( metrics_request_flag & V_HEX_MAX_EDGE_RATIO )
02670         {
02671             double max_edge_ratio_12, max_edge_ratio_13, max_edge_ratio_23;
02672 
02673             double mag_efg1 = efg1.length();
02674             double mag_efg2 = efg2.length();
02675             double mag_efg3 = efg3.length();
02676 
02677             max_edge_ratio_12 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg2 ), VERDICT_MIN( mag_efg1, mag_efg2 ) );
02678             max_edge_ratio_13 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg3 ), VERDICT_MIN( mag_efg1, mag_efg3 ) );
02679             max_edge_ratio_23 = safe_ratio( VERDICT_MAX( mag_efg2, mag_efg3 ), VERDICT_MIN( mag_efg2, mag_efg3 ) );
02680 
02681             metric_vals->max_edge_ratio =
02682                 (double)VERDICT_MAX( max_edge_ratio_12, VERDICT_MAX( max_edge_ratio_13, max_edge_ratio_23 ) );
02683         }
02684 
02685         if( metrics_request_flag & V_HEX_SKEW )
02686         {
02687 
02688             VerdictVector vec1 = efg1;
02689             VerdictVector vec2 = efg2;
02690             VerdictVector vec3 = efg3;
02691 
02692             if( vec1.normalize() <= VERDICT_DBL_MIN || vec2.normalize() <= VERDICT_DBL_MIN ||
02693                 vec3.normalize() <= VERDICT_DBL_MIN )
02694             {
02695                 metric_vals->skew = (double)VERDICT_DBL_MAX;
02696             }
02697             else
02698             {
02699                 double skewx = fabs( vec1 % vec2 );
02700                 double skewy = fabs( vec1 % vec3 );
02701                 double skewz = fabs( vec2 % vec3 );
02702 
02703                 metric_vals->skew = (double)( VERDICT_MAX( skewx, VERDICT_MAX( skewy, skewz ) ) );
02704             }
02705         }
02706 
02707         if( metrics_request_flag & V_HEX_TAPER )
02708         {
02709             VerdictVector efg12 = calc_hex_efg( 12, node_pos );
02710             VerdictVector efg13 = calc_hex_efg( 13, node_pos );
02711             VerdictVector efg23 = calc_hex_efg( 23, node_pos );
02712 
02713             double taperx = fabs( safe_ratio( efg12.length(), VERDICT_MIN( efg1.length(), efg2.length() ) ) );
02714             double tapery = fabs( safe_ratio( efg13.length(), VERDICT_MIN( efg1.length(), efg3.length() ) ) );
02715             double taperz = fabs( safe_ratio( efg23.length(), VERDICT_MIN( efg2.length(), efg3.length() ) ) );
02716 
02717             metric_vals->taper = (double)VERDICT_MAX( taperx, VERDICT_MAX( tapery, taperz ) );
02718         }
02719     }
02720 
02721     if( metrics_request_flag & V_HEX_VOLUME )
02722     {
02723         metric_vals->volume = v_hex_volume( 8, coordinates );
02724     }
02725 
02726     if( metrics_request_flag &
02727         ( V_HEX_JACOBIAN | V_HEX_SCALED_JACOBIAN | V_HEX_CONDITION | V_HEX_SHEAR | V_HEX_SHAPE |
02728           V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) )
02729     {
02730 
02731         static const double two_thirds = 2.0 / 3.0;
02732         VerdictVector edges[12];
02733         // the length squares
02734         double length_squared[12];
02735         // make vectors from coordinates
02736         make_hex_edges( coordinates, edges );
02737 
02738         // calculate the length squares if we need to
02739         if( metrics_request_flag &
02740             ( V_HEX_JACOBIAN | V_HEX_SHEAR | V_HEX_SCALED_JACOBIAN | V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE |
02741               V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) )
02742             make_edge_length_squares( edges, length_squared );
02743 
02744         double jacobian = VERDICT_DBL_MAX, scaled_jacobian = VERDICT_DBL_MAX, condition = 0.0, shear = 1.0, shape = 1.0,
02745                oddy = 0.0;
02746         double current_jacobian, current_scaled_jacobian, current_condition, current_shape, detw = 0, det_sum = 0,
02747                                                                                             current_oddy;
02748         VerdictBoolean rel_size_error = VERDICT_FALSE;
02749 
02750         VerdictVector xxi, xet, xze;
02751 
02752         // get weights if we need based on average size of a hex
02753         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
02754         {
02755             v_hex_get_weight( xxi, xet, xze );
02756             detw = xxi % ( xet * xze );
02757             if( detw < VERDICT_DBL_MIN ) rel_size_error = VERDICT_TRUE;
02758         }
02759 
02760         xxi = edges[0] - edges[2] + edges[4] - edges[6];
02761         xet = edges[1] - edges[3] + edges[5] - edges[7];
02762         xze = edges[8] + edges[9] + edges[10] + edges[11];
02763 
02764         current_jacobian = xxi % ( xet * xze ) / 64.0;
02765         if( current_jacobian < jacobian ) jacobian = current_jacobian;
02766 
02767         if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
02768         {
02769             current_jacobian *= 64.0;
02770             current_scaled_jacobian =
02771                 current_jacobian / sqrt( xxi.length_squared() * xet.length_squared() * xze.length_squared() );
02772             if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
02773         }
02774 
02775         if( metrics_request_flag & V_HEX_CONDITION )
02776         {
02777             current_condition = condition_comp( xxi, xet, xze );
02778             if( current_condition > condition )
02779             {
02780                 condition = current_condition;
02781             }
02782         }
02783 
02784         if( metrics_request_flag & V_HEX_ODDY )
02785         {
02786             current_oddy = oddy_comp( xxi, xet, xze );
02787             if( current_oddy > oddy )
02788             {
02789                 oddy = current_oddy;
02790             }
02791         }
02792 
02793         // J(0,0,0)
02794         current_jacobian = edges[0] % ( -edges[3] * edges[8] );
02795         if( current_jacobian < jacobian ) jacobian = current_jacobian;
02796 
02797         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
02798         {
02799             det_sum += current_jacobian;
02800         }
02801 
02802         if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
02803         {
02804             if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN ||
02805                 length_squared[8] <= VERDICT_DBL_MIN )
02806             {
02807                 current_scaled_jacobian = VERDICT_DBL_MAX;
02808             }
02809             else
02810             {
02811                 current_scaled_jacobian =
02812                     current_jacobian / sqrt( length_squared[0] * length_squared[3] * length_squared[8] );
02813             }
02814             if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
02815         }
02816 
02817         if( metrics_request_flag & V_HEX_CONDITION )
02818         {
02819             current_condition = condition_comp( edges[0], -edges[3], edges[8] );
02820             if( current_condition > condition )
02821             {
02822                 condition = current_condition;
02823             }
02824         }
02825 
02826         if( metrics_request_flag & V_HEX_ODDY )
02827         {
02828             current_oddy = oddy_comp( edges[0], -edges[3], edges[8] );
02829             if( current_oddy > oddy )
02830             {
02831                 oddy = current_oddy;
02832             }
02833         }
02834 
02835         if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
02836         {
02837             if( current_jacobian > VERDICT_DBL_MIN )
02838                 current_shape = 3 * pow( current_jacobian, two_thirds ) /
02839                                 ( length_squared[0] + length_squared[3] + length_squared[8] );
02840             else
02841                 current_shape = 0;
02842 
02843             if( current_shape < shape )
02844             {
02845                 shape = current_shape;
02846             }
02847         }
02848 
02849         // J(1,0,0)
02850         current_jacobian = edges[1] % ( -edges[0] * edges[9] );
02851         if( current_jacobian < jacobian ) jacobian = current_jacobian;
02852 
02853         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
02854         {
02855             det_sum += current_jacobian;
02856         }
02857 
02858         if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
02859         {
02860             if( length_squared[1] <= VERDICT_DBL_MIN || length_squared[0] <= VERDICT_DBL_MIN ||
02861                 length_squared[9] <= VERDICT_DBL_MIN )
02862             {
02863                 current_scaled_jacobian = VERDICT_DBL_MAX;
02864             }
02865             else
02866             {
02867                 current_scaled_jacobian =
02868                     current_jacobian / sqrt( length_squared[1] * length_squared[0] * length_squared[9] );
02869             }
02870             if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
02871         }
02872 
02873         if( metrics_request_flag & V_HEX_CONDITION )
02874         {
02875             current_condition = condition_comp( edges[1], -edges[0], edges[9] );
02876             if( current_condition > condition )
02877             {
02878                 condition = current_condition;
02879             }
02880         }
02881 
02882         if( metrics_request_flag & V_HEX_ODDY )
02883         {
02884             current_oddy = oddy_comp( edges[1], -edges[0], edges[9] );
02885             if( current_oddy > oddy )
02886             {
02887                 oddy = current_oddy;
02888             }
02889         }
02890 
02891         if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
02892         {
02893             if( current_jacobian > VERDICT_DBL_MIN )
02894                 current_shape = 3 * pow( current_jacobian, two_thirds ) /
02895                                 ( length_squared[1] + length_squared[0] + length_squared[9] );
02896             else
02897                 current_shape = 0;
02898 
02899             if( current_shape < shape )
02900             {
02901                 shape = current_shape;
02902             }
02903         }
02904 
02905         // J(1,1,0)
02906         current_jacobian = edges[2] % ( -edges[1] * edges[10] );
02907         if( current_jacobian < jacobian ) jacobian = current_jacobian;
02908 
02909         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
02910         {
02911             det_sum += current_jacobian;
02912         }
02913 
02914         if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
02915         {
02916             if( length_squared[2] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN ||
02917                 length_squared[10] <= VERDICT_DBL_MIN )
02918             {
02919                 current_scaled_jacobian = VERDICT_DBL_MAX;
02920             }
02921             else
02922             {
02923                 current_scaled_jacobian =
02924                     current_jacobian / sqrt( length_squared[2] * length_squared[1] * length_squared[10] );
02925             }
02926             if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
02927         }
02928 
02929         if( metrics_request_flag & V_HEX_CONDITION )
02930         {
02931             current_condition = condition_comp( edges[2], -edges[1], edges[10] );
02932             if( current_condition > condition )
02933             {
02934                 condition = current_condition;
02935             }
02936         }
02937 
02938         if( metrics_request_flag & V_HEX_ODDY )
02939         {
02940             current_oddy = oddy_comp( edges[2], -edges[1], edges[10] );
02941             if( current_oddy > oddy )
02942             {
02943                 oddy = current_oddy;
02944             }
02945         }
02946 
02947         if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
02948         {
02949             if( current_jacobian > VERDICT_DBL_MIN )
02950                 current_shape = 3 * pow( current_jacobian, two_thirds ) /
02951                                 ( length_squared[2] + length_squared[1] + length_squared[10] );
02952             else
02953                 current_shape = 0;
02954 
02955             if( current_shape < shape )
02956             {
02957                 shape = current_shape;
02958             }
02959         }
02960 
02961         // J(0,1,0)
02962         current_jacobian = edges[3] % ( -edges[2] * edges[11] );
02963         if( current_jacobian < jacobian ) jacobian = current_jacobian;
02964 
02965         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
02966         {
02967             det_sum += current_jacobian;
02968         }
02969 
02970         if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
02971         {
02972             if( length_squared[3] <= VERDICT_DBL_MIN || length_squared[2] <= VERDICT_DBL_MIN ||
02973                 length_squared[11] <= VERDICT_DBL_MIN )
02974             {
02975                 current_scaled_jacobian = VERDICT_DBL_MAX;
02976             }
02977             else
02978             {
02979                 current_scaled_jacobian =
02980                     current_jacobian / sqrt( length_squared[3] * length_squared[2] * length_squared[11] );
02981             }
02982             if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
02983         }
02984 
02985         if( metrics_request_flag & V_HEX_CONDITION )
02986         {
02987             current_condition = condition_comp( edges[3], -edges[2], edges[11] );
02988             if( current_condition > condition )
02989             {
02990                 condition = current_condition;
02991             }
02992         }
02993 
02994         if( metrics_request_flag & V_HEX_ODDY )
02995         {
02996             current_oddy = oddy_comp( edges[3], -edges[2], edges[11] );
02997             if( current_oddy > oddy )
02998             {
02999                 oddy = current_oddy;
03000             }
03001         }
03002 
03003         if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
03004         {
03005             if( current_jacobian > VERDICT_DBL_MIN )
03006                 current_shape = 3 * pow( current_jacobian, two_thirds ) /
03007                                 ( length_squared[3] + length_squared[2] + length_squared[11] );
03008             else
03009                 current_shape = 0;
03010 
03011             if( current_shape < shape )
03012             {
03013                 shape = current_shape;
03014             }
03015         }
03016 
03017         // J(0,0,1)
03018         current_jacobian = edges[4] % ( -edges[8] * -edges[7] );
03019         if( current_jacobian < jacobian ) jacobian = current_jacobian;
03020 
03021         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
03022         {
03023             det_sum += current_jacobian;
03024         }
03025 
03026         if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
03027         {
03028             if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[8] <= VERDICT_DBL_MIN ||
03029                 length_squared[7] <= VERDICT_DBL_MIN )
03030             {
03031                 current_scaled_jacobian = VERDICT_DBL_MAX;
03032             }
03033             else
03034             {
03035                 current_scaled_jacobian =
03036                     current_jacobian / sqrt( length_squared[4] * length_squared[8] * length_squared[7] );
03037             }
03038             if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
03039         }
03040 
03041         if( metrics_request_flag & V_HEX_CONDITION )
03042         {
03043             current_condition = condition_comp( edges[4], -edges[8], -edges[7] );
03044             if( current_condition > condition )
03045             {
03046                 condition = current_condition;
03047             }
03048         }
03049 
03050         if( metrics_request_flag & V_HEX_ODDY )
03051         {
03052             current_oddy = oddy_comp( edges[4], -edges[8], -edges[7] );
03053             if( current_oddy > oddy )
03054             {
03055                 oddy = current_oddy;
03056             }
03057         }
03058 
03059         if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
03060         {
03061             if( current_jacobian > VERDICT_DBL_MIN )
03062                 current_shape = 3 * pow( current_jacobian, two_thirds ) /
03063                                 ( length_squared[4] + length_squared[8] + length_squared[7] );
03064             else
03065                 current_shape = 0;
03066 
03067             if( current_shape < shape )
03068             {
03069                 shape = current_shape;
03070             }
03071         }
03072 
03073         // J(1,0,1)
03074         current_jacobian = -edges[4] % ( edges[5] * -edges[9] );
03075         if( current_jacobian < jacobian ) jacobian = current_jacobian;
03076 
03077         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
03078         {
03079             det_sum += current_jacobian;
03080         }
03081 
03082         if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
03083         {
03084             if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[5] <= VERDICT_DBL_MIN ||
03085                 length_squared[9] <= VERDICT_DBL_MIN )
03086             {
03087                 current_scaled_jacobian = VERDICT_DBL_MAX;
03088             }
03089             else
03090             {
03091                 current_scaled_jacobian =
03092                     current_jacobian / sqrt( length_squared[4] * length_squared[5] * length_squared[9] );
03093             }
03094             if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
03095         }
03096 
03097         if( metrics_request_flag & V_HEX_CONDITION )
03098         {
03099             current_condition = condition_comp( -edges[4], edges[5], -edges[9] );
03100             if( current_condition > condition )
03101             {
03102                 condition = current_condition;
03103             }
03104         }
03105 
03106         if( metrics_request_flag & V_HEX_ODDY )
03107         {
03108             current_oddy = oddy_comp( -edges[4], edges[5], -edges[9] );
03109             if( current_oddy > oddy )
03110             {
03111                 oddy = current_oddy;
03112             }
03113         }
03114 
03115         if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
03116         {
03117             if( current_jacobian > VERDICT_DBL_MIN )
03118                 current_shape = 3 * pow( current_jacobian, two_thirds ) /
03119                                 ( length_squared[4] + length_squared[5] + length_squared[9] );
03120             else
03121                 current_shape = 0;
03122 
03123             if( current_shape < shape )
03124             {
03125                 shape = current_shape;
03126             }
03127         }
03128 
03129         // J(1,1,1)
03130         current_jacobian = -edges[5] % ( edges[6] * -edges[10] );
03131         if( current_jacobian < jacobian ) jacobian = current_jacobian;
03132 
03133         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
03134         {
03135             det_sum += current_jacobian;
03136         }
03137 
03138         if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
03139         {
03140             if( length_squared[5] <= VERDICT_DBL_MIN || length_squared[6] <= VERDICT_DBL_MIN ||
03141                 length_squared[10] <= VERDICT_DBL_MIN )
03142             {
03143                 current_scaled_jacobian = VERDICT_DBL_MAX;
03144             }
03145             else
03146             {
03147                 current_scaled_jacobian =
03148                     current_jacobian / sqrt( length_squared[5] * length_squared[6] * length_squared[10] );
03149             }
03150             if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
03151         }
03152 
03153         if( metrics_request_flag & V_HEX_CONDITION )
03154         {
03155             current_condition = condition_comp( -edges[5], edges[6], -edges[10] );
03156             if( current_condition > condition )
03157             {
03158                 condition = current_condition;
03159             }
03160         }
03161 
03162         if( metrics_request_flag & V_HEX_ODDY )
03163         {
03164             current_oddy = oddy_comp( -edges[5], edges[6], -edges[10] );
03165             if( current_oddy > oddy )
03166             {
03167                 oddy = current_oddy;
03168             }
03169         }
03170 
03171         if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
03172         {
03173             if( current_jacobian > VERDICT_DBL_MIN )
03174                 current_shape = 3 * pow( current_jacobian, two_thirds ) /
03175                                 ( length_squared[5] + length_squared[6] + length_squared[10] );
03176             else
03177                 current_shape = 0;
03178 
03179             if( current_shape < shape )
03180             {
03181                 shape = current_shape;
03182             }
03183         }
03184 
03185         // J(0,1,1)
03186         current_jacobian = -edges[6] % ( edges[7] * -edges[11] );
03187         if( current_jacobian < jacobian ) jacobian = current_jacobian;
03188 
03189         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
03190         {
03191             det_sum += current_jacobian;
03192         }
03193 
03194         if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
03195         {
03196             if( length_squared[6] <= VERDICT_DBL_MIN || length_squared[7] <= VERDICT_DBL_MIN ||
03197                 length_squared[11] <= VERDICT_DBL_MIN )
03198             {
03199                 current_scaled_jacobian = VERDICT_DBL_MAX;
03200             }
03201             else
03202             {
03203                 current_scaled_jacobian =
03204                     current_jacobian / sqrt( length_squared[6] * length_squared[7] * length_squared[11] );
03205             }
03206             if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
03207         }
03208 
03209         if( metrics_request_flag & V_HEX_CONDITION )
03210         {
03211             current_condition = condition_comp( -edges[6], edges[7], -edges[11] );
03212             if( current_condition > condition )
03213             {
03214                 condition = current_condition;
03215             }
03216         }
03217 
03218         if( metrics_request_flag & V_HEX_ODDY )
03219         {
03220             current_oddy = oddy_comp( -edges[6], edges[7], -edges[11] );
03221             if( current_oddy > oddy )
03222             {
03223                 oddy = current_oddy;
03224             }
03225         }
03226 
03227         if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
03228         {
03229             if( current_jacobian > VERDICT_DBL_MIN )
03230                 current_shape = 3 * pow( current_jacobian, two_thirds ) /
03231                                 ( length_squared[6] + length_squared[7] + length_squared[11] );
03232             else
03233                 current_shape = 0;
03234 
03235             if( current_shape < shape )
03236             {
03237                 shape = current_shape;
03238             }
03239         }
03240 
03241         if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
03242         {
03243             if( det_sum > VERDICT_DBL_MIN && rel_size_error != VERDICT_TRUE )
03244             {
03245                 double tau                         = det_sum / ( 8 * detw );
03246                 metric_vals->relative_size_squared = (double)VERDICT_MIN( tau * tau, 1.0 / tau / tau );
03247             }
03248             else
03249                 metric_vals->relative_size_squared = 0.0;
03250         }
03251 
03252         // set values from above calculations
03253         if( metrics_request_flag & V_HEX_JACOBIAN ) metric_vals->jacobian = (double)jacobian;
03254 
03255         if( metrics_request_flag & V_HEX_SCALED_JACOBIAN ) metric_vals->scaled_jacobian = (double)scaled_jacobian;
03256 
03257         if( metrics_request_flag & V_HEX_CONDITION ) metric_vals->condition = (double)( condition / 3.0 );
03258 
03259         if( metrics_request_flag & V_HEX_SHEAR )
03260         {
03261             if( shear < VERDICT_DBL_MIN )  // shear has range 0 to +1
03262                 shear = 0;
03263             metric_vals->shear = (double)shear;
03264         }
03265 
03266         if( metrics_request_flag & V_HEX_SHAPE ) metric_vals->shape = (double)shape;
03267 
03268         if( metrics_request_flag & V_HEX_SHAPE_AND_SIZE )
03269             metric_vals->shape_and_size = (double)( shape * metric_vals->relative_size_squared );
03270 
03271         if( metrics_request_flag & V_HEX_SHEAR_AND_SIZE )
03272             metric_vals->shear_and_size = (double)( shear * metric_vals->relative_size_squared );
03273 
03274         if( metrics_request_flag & V_HEX_ODDY ) metric_vals->oddy = (double)oddy;
03275 
03276         if( metrics_request_flag & V_HEX_STRETCH )
03277         {
03278             static const double HEX_STRETCH_SCALE_FACTOR = sqrt( 3.0 );
03279             double min_edge                              = length_squared[0];
03280             for( int j = 1; j < 12; j++ )
03281                 min_edge = VERDICT_MIN( min_edge, length_squared[j] );
03282 
03283             double max_diag = diag_length( 1, coordinates );
03284 
03285             metric_vals->stretch = (double)( HEX_STRETCH_SCALE_FACTOR * ( safe_ratio( sqrt( min_edge ), max_diag ) ) );
03286         }
03287     }
03288 
03289     if( metrics_request_flag & V_HEX_DIAGONAL ) metric_vals->diagonal = v_hex_diagonal( num_nodes, coordinates );
03290     if( metrics_request_flag & V_HEX_DIMENSION ) metric_vals->dimension = v_hex_dimension( num_nodes, coordinates );
03291     if( metrics_request_flag & V_HEX_DISTORTION ) metric_vals->distortion = v_hex_distortion( num_nodes, coordinates );
03292 
03293     // take care of any overflow problems
03294     // max_edge_ratio
03295     if( metric_vals->max_edge_ratio > 0 )
03296         metric_vals->max_edge_ratio = (double)VERDICT_MIN( metric_vals->max_edge_ratio, VERDICT_DBL_MAX );
03297     else
03298         metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX );
03299 
03300     // skew
03301     if( metric_vals->skew > 0 )
03302         metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX );
03303     else
03304         metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX );
03305 
03306     // taper
03307     if( metric_vals->taper > 0 )
03308         metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX );
03309     else
03310         metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX );
03311 
03312     // volume
03313     if( metric_vals->volume > 0 )
03314         metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
03315     else
03316         metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
03317 
03318     // stretch
03319     if( metric_vals->stretch > 0 )
03320         metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX );
03321     else
03322         metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX );
03323 
03324     // diagonal
03325     if( metric_vals->diagonal > 0 )
03326         metric_vals->diagonal = (double)VERDICT_MIN( metric_vals->diagonal, VERDICT_DBL_MAX );
03327     else
03328         metric_vals->diagonal = (double)VERDICT_MAX( metric_vals->diagonal, -VERDICT_DBL_MAX );
03329 
03330     // dimension
03331     if( metric_vals->dimension > 0 )
03332         metric_vals->dimension = (double)VERDICT_MIN( metric_vals->dimension, VERDICT_DBL_MAX );
03333     else
03334         metric_vals->dimension = (double)VERDICT_MAX( metric_vals->dimension, -VERDICT_DBL_MAX );
03335 
03336     // oddy
03337     if( metric_vals->oddy > 0 )
03338         metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX );
03339     else
03340         metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX );
03341 
03342     // condition
03343     if( metric_vals->condition > 0 )
03344         metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
03345     else
03346         metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
03347 
03348     // jacobian
03349     if( metric_vals->jacobian > 0 )
03350         metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
03351     else
03352         metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
03353 
03354     // scaled_jacobian
03355     if( metric_vals->scaled_jacobian > 0 )
03356         metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
03357     else
03358         metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
03359 
03360     // shear
03361     if( metric_vals->shear > 0 )
03362         metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX );
03363     else
03364         metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX );
03365 
03366     // shape
03367     if( metric_vals->shape > 0 )
03368         metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
03369     else
03370         metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
03371 
03372     // relative_size_squared
03373     if( metric_vals->relative_size_squared > 0 )
03374         metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
03375     else
03376         metric_vals->relative_size_squared =
03377             (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
03378 
03379     // shape_and_size
03380     if( metric_vals->shape_and_size > 0 )
03381         metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
03382     else
03383         metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
03384 
03385     // shear_and_size
03386     if( metric_vals->shear_and_size > 0 )
03387         metric_vals->shear_and_size = (double)VERDICT_MIN( metric_vals->shear_and_size, VERDICT_DBL_MAX );
03388     else
03389         metric_vals->shear_and_size = (double)VERDICT_MAX( metric_vals->shear_and_size, -VERDICT_DBL_MAX );
03390 
03391     // distortion
03392     if( metric_vals->distortion > 0 )
03393         metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
03394     else
03395         metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
03396 
03397     if( metrics_request_flag & V_HEX_MED_ASPECT_FROBENIUS )
03398         metric_vals->med_aspect_frobenius = v_hex_med_aspect_frobenius( 8, coordinates );
03399 
03400     if( metrics_request_flag & V_HEX_EDGE_RATIO ) metric_vals->edge_ratio = v_hex_edge_ratio( 8, coordinates );
03401 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines