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