V_HexMetric.cpp

/*=========================================================================

  Module:    V_HexMetric.cpp

  Copyright 2006 National Technology & Engineering Solutions of Sandia,
  LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS,
  the U.S. Government retains certain rights in this software.

  See LICENSE for details.

=========================================================================*/

/*
 *
 * HexMetric.cpp contains quality calculations for hexes
 *
 * This file is part of VERDICT
 *
 */

#include "verdict.h"
#include "VerdictVector.hpp"
#include "V_GaussIntegration.hpp"
#include "verdict_defines.hpp"
#include <V_HexMetric.hpp>
#include <memory.h>
#include <vector>
#include <algorithm>
#include <cmath>   // for std::isnan

#if defined(__BORLANDC__)
#pragma warn -8004 /* "assigned a value that is never used" */
#endif

namespace VERDICT_NAMESPACE
{

extern void quad_minimum_maximum_angle( double min_max_angles[2], double coordinates[][3] );

static const double two_thirds = 2.0/3.0;
static const double third=1.0/3.0;
static const double root_of_3 = sqrt(3.0);


//! weights based on the average size of a hex
static int hex_get_weight( VerdictVector &v1, VerdictVector &v2, VerdictVector &v3, double average_size)
{
  if ( average_size  == 0)
    return 0;

  v1.set(1,0,0);
  v2.set(0,1,0);
  v3.set(0,0,1);

  double scale = pow(average_size / (VerdictVector::Dot(v1, (v2 * v3))), 0.33333333333);
  v1 *= scale;
  v2 *= scale;
  v3 *= scale;

  return 1;
}


static const double HEX27_node_local_coord[27][3] =
{
  {-1,-1,-1},
  {1,-1,-1},
  {1,1,-1},
  {-1,1,-1},
  {-1,-1,1},
  {1,-1,1},
  {1,1,1},
  {-1,1,1},
  {0,-1,-1},
  {1,0,-1},
  {0,1,-1},
  {-1,0,-1},
  {-1,-1,0},
  {1,-1,0},
  {1,1,0},
  {-1,1,0},
  {0,-1,1},
  {1,0,1},
  {0,1,1},
  {-1,0,1},
  {0,0,0},
  {0,0,-1},
  {0,0,1},
  {-1,0,0},
  {1,0,0},
  {0,-1,0},
  {0,1,0}
};

static void HEX27_gradients_of_the_shape_functions_for_RST(const double rst[3], double dhdr[27],double dhds[27],double dhdt[27])
{
  double g1r = -0.5 * rst[0] * (1 - rst[0]);
  double g1s = -0.5 * rst[1] * (1 - rst[1]);
  double g1t = -0.5 * rst[2] * (1 - rst[2]);

  double g2r = (1 + rst[0])*(1 - rst[0]);
  double g2s = (1 + rst[1])*(1 - rst[1]);
  double g2t = (1 + rst[2])*(1 - rst[2]);

  double g3r = 0.5 * rst[0] * (1 + rst[0]);
  double g3s = 0.5 * rst[1] * (1 + rst[1]);
  double g3t = 0.5 * rst[2] * (1 + rst[2]);

  double g1r_r = rst[0] - 0.5;
  double g1s_s = rst[1] - 0.5;
  double g1t_t = rst[2] - 0.5;

  double g2r_r = -2*rst[0];
  double g2s_s = -2*rst[1];
  double g2t_t = -2*rst[2];

  double g3r_r = rst[0] + 0.5;
  double g3s_s = rst[1] + 0.5;
  double g3t_t = rst[2] + 0.5;

  //r-derivatives
  dhdr[0] = g1r_r * g1s * g1t;
  dhdr[1] = g3r_r * g1s * g1t;
  dhdr[2] = g3r_r * g3s * g1t;
  dhdr[3] = g1r_r * g3s * g1t;
  dhdr[4] = g1r_r * g1s * g3t;
  dhdr[5] = g3r_r * g1s * g3t;
  dhdr[6] = g3r_r * g3s * g3t;
  dhdr[7] = g1r_r * g3s * g3t;
  dhdr[8] = g2r_r * g1s * g1t;
  dhdr[9] = g3r_r * g2s * g1t;
  dhdr[10]= g2r_r * g3s * g1t;
  dhdr[11]= g1r_r * g2s * g1t;
  dhdr[16]= g2r_r * g1s * g3t;
  dhdr[17]= g3r_r * g2s * g3t;
  dhdr[18]= g2r_r * g3s * g3t;
  dhdr[19]= g1r_r * g2s * g3t;
  dhdr[12]= g1r_r * g1s * g2t;
  dhdr[13]= g3r_r * g1s * g2t;
  dhdr[14]= g3r_r * g3s * g2t;
  dhdr[15]= g1r_r * g3s * g2t;
  dhdr[23]= g1r_r * g2s * g2t;
  dhdr[24]= g3r_r * g2s * g2t;
  dhdr[25]= g2r_r * g1s * g2t;
  dhdr[26]= g2r_r * g3s * g2t;
  dhdr[21]= g2r_r * g2s * g1t;
  dhdr[22]= g2r_r * g2s * g3t;
  dhdr[20]= g2r_r * g2s * g2t;

  //s-derivatives
  dhds[0] = g1r * g1s_s * g1t;
  dhds[1] = g3r * g1s_s * g1t;
  dhds[2] = g3r * g3s_s * g1t;
  dhds[3] = g1r * g3s_s * g1t;
  dhds[4] = g1r * g1s_s * g3t;
  dhds[5] = g3r * g1s_s * g3t;
  dhds[6] = g3r * g3s_s * g3t;
  dhds[7] = g1r * g3s_s * g3t;
  dhds[8] = g2r * g1s_s * g1t;
  dhds[9] = g3r * g2s_s * g1t;
  dhds[10] = g2r * g3s_s * g1t;
  dhds[11] = g1r * g2s_s * g1t;
  dhds[16] = g2r * g1s_s * g3t;
  dhds[17] = g3r * g2s_s * g3t;
  dhds[18] = g2r * g3s_s * g3t;
  dhds[19] = g1r * g2s_s * g3t;
  dhds[12] = g1r * g1s_s * g2t;
  dhds[13] = g3r * g1s_s * g2t;
  dhds[14] = g3r * g3s_s * g2t;
  dhds[15] = g1r * g3s_s * g2t;
  dhds[23] = g1r * g2s_s * g2t;
  dhds[24] = g3r * g2s_s * g2t;
  dhds[25] = g2r * g1s_s * g2t;
  dhds[26] = g2r * g3s_s * g2t;
  dhds[21] = g2r * g2s_s * g1t;
  dhds[22] = g2r * g2s_s * g3t;
  dhds[20] = g2r * g2s_s * g2t;

  //t-derivatives
  dhdt[0] = g1r * g1s * g1t_t;
  dhdt[1] = g3r * g1s * g1t_t;
  dhdt[2] = g3r * g3s * g1t_t;
  dhdt[3] = g1r * g3s * g1t_t;
  dhdt[4] = g1r * g1s * g3t_t;
  dhdt[5] = g3r * g1s * g3t_t;
  dhdt[6] = g3r * g3s * g3t_t;
  dhdt[7] = g1r * g3s * g3t_t;
  dhdt[8] = g2r * g1s * g1t_t;
  dhdt[9] = g3r * g2s * g1t_t;
  dhdt[10] = g2r * g3s * g1t_t;
  dhdt[11] = g1r * g2s * g1t_t;
  dhdt[16] = g2r * g1s * g3t_t;
  dhdt[17] = g3r * g2s * g3t_t;
  dhdt[18] = g2r * g3s * g3t_t;
  dhdt[19] = g1r * g2s * g3t_t;
  dhdt[12] = g1r * g1s * g2t_t;
  dhdt[13] = g3r * g1s * g2t_t;
  dhdt[14] = g3r * g3s * g2t_t;
  dhdt[15] = g1r * g3s * g2t_t;
  dhdt[23] = g1r * g2s * g2t_t;
  dhdt[24] = g3r * g2s * g2t_t;
  dhdt[25] = g2r * g1s * g2t_t;
  dhdt[26] = g2r * g3s * g2t_t;
  dhdt[21] = g2r * g2s * g1t_t;
  dhdt[22] = g2r * g2s * g3t_t;
  dhdt[20] = g2r * g2s * g2t_t;

  for(int i=0; i<27; i++)
  {
    dhdr[i] *= 2;
    dhds[i] *= 2;
    dhdt[i] *= 2;
  }
}

#define make_hex_nodes(coord, pos)                      \
  for (int mhcii = 0; mhcii < 8; mhcii++ )              \
  {                                                     \
    pos[mhcii].set( coord[mhcii][0],                    \
        coord[mhcii][1],                                \
        coord[mhcii][2] );                              \
  }


#define make_edge_length_squares(edges, lengths)        \
{                                                       \
  for(int melii=0; melii<12; melii++)                   \
    lengths[melii] = edges[melii].length_squared();     \
}
  

//! make VerdictVectors from coordinates
static void make_hex_edges( double coordinates[][3], VerdictVector edges[12] )
{
  edges[0].set(
      coordinates[1][0] - coordinates[0][0],
      coordinates[1][1] - coordinates[0][1],
      coordinates[1][2] - coordinates[0][2]
      );
  edges[1].set(
      coordinates[2][0] - coordinates[1][0],
      coordinates[2][1] - coordinates[1][1],
      coordinates[2][2] - coordinates[1][2]
      );
  edges[2].set(
      coordinates[3][0] - coordinates[2][0],
      coordinates[3][1] - coordinates[2][1],
      coordinates[3][2] - coordinates[2][2]
      );
  edges[3].set(
      coordinates[0][0] - coordinates[3][0],
      coordinates[0][1] - coordinates[3][1],
      coordinates[0][2] - coordinates[3][2]
      );
  edges[4].set(
      coordinates[5][0] - coordinates[4][0],
      coordinates[5][1] - coordinates[4][1],
      coordinates[5][2] - coordinates[4][2]
      );
  edges[5].set(
      coordinates[6][0] - coordinates[5][0],
      coordinates[6][1] - coordinates[5][1],
      coordinates[6][2] - coordinates[5][2]
      );
  edges[6].set(
      coordinates[7][0] - coordinates[6][0],
      coordinates[7][1] - coordinates[6][1],
      coordinates[7][2] - coordinates[6][2]
      );
  edges[7].set(
      coordinates[4][0] - coordinates[7][0],
      coordinates[4][1] - coordinates[7][1],
      coordinates[4][2] - coordinates[7][2]
      );
  edges[8].set(
      coordinates[4][0] - coordinates[0][0],
      coordinates[4][1] - coordinates[0][1],
      coordinates[4][2] - coordinates[0][2]
      );
  edges[9].set(
      coordinates[5][0] - coordinates[1][0],
      coordinates[5][1] - coordinates[1][1],
      coordinates[5][2] - coordinates[1][2]
      );
  edges[10].set(
      coordinates[6][0] - coordinates[2][0],
      coordinates[6][1] - coordinates[2][1],
      coordinates[6][2] - coordinates[2][2]
      );
  edges[11].set(
      coordinates[7][0] - coordinates[3][0],
      coordinates[7][1] - coordinates[3][1],
      coordinates[7][2] - coordinates[3][2]
      );
}

#if 0 /* Not currently used and not exposed in verdict.h */
/*!
  localizes hex coordinates
*/
static void localize_hex_coordinates(double coordinates[][3], VerdictVector position[8] )
{

  int ii;
  for ( ii = 0; ii < 8; ii++ )
  {
    position[ii].set( coordinates[ii][0],
        coordinates[ii][1],
        coordinates[ii][2] );
  }
  
  // ... Make centroid of element the center of coordinate system
  VerdictVector point_1256 = position[1];
  point_1256 += position[2];
  point_1256 += position[5];
  point_1256 += position[6];

  VerdictVector point_0374 = position[0];
  point_0374 += position[3];
  point_0374 += position[7];
  point_0374 += position[4];

  VerdictVector centroid = point_1256;
  centroid += point_0374;
  centroid /= 8.0;

  int i;
  for ( i = 0; i < 8; i++)
    position[i] -= centroid;

  // ... Rotate element such that center of side 1-2 and 0-3 define X axis
  double DX = point_1256.x() - point_0374.x();
  double DY = point_1256.y() - point_0374.y();
  double DZ = point_1256.z() - point_0374.z();

  double AMAGX = sqrt(DX*DX + DZ*DZ);
  double AMAGY = sqrt(DX*DX + DY*DY + DZ*DZ);
  double FMAGX = AMAGX == 0.0 ? 1.0 : 0.0;
  double FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;

  double CZ = DX / (AMAGX + FMAGX) + FMAGX;
  double SZ = DZ / (AMAGX + FMAGX);
  double CY = AMAGX / (AMAGY + FMAGY) + FMAGY;
  double SY = DY / (AMAGY + FMAGY);

  double temp;
 
  for (i = 0; i < 8; i++) 
  {
    temp =  CY * CZ * position[i].x() + CY * SZ * position[i].z() +
      SY * position[i].y();
    position[i].y( -SY * CZ * position[i].x() - SY * SZ * position[i].z() +
        CY * position[i].y());
    position[i].z( -SZ * position[i].x() + CZ * position[i].z());
    position[i].x(temp);
  }


  // ... Now, rotate about Y
  VerdictVector delta = -position[0];
  delta -= position[1];
  delta += position[2];
  delta += position[3];
  delta -= position[4];
  delta -= position[5];
  delta += position[6];
  delta += position[7];

  DY = delta.y();
  DZ = delta.z();

  AMAGY = sqrt(DY*DY + DZ*DZ);
  FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;

  double CX = DY / (AMAGY + FMAGY) + FMAGY;
  double SX = DZ / (AMAGY + FMAGY);
  
  for (i = 0; i < 8; i++) 
    {
    temp =  CX * position[i].y() + SX * position[i].z();
    position[i].z(-SX * position[i].y() + CX * position[i].z());
    position[i].y(temp);
    }
}

static double safe_ratio3( const double numerator, 
    const double denominator,
    const double max_ratio )
{
    // this filter is essential for good running time in practice
  double return_value;

  const double filter_n = max_ratio * 1.0e-16;
  const double filter_d = 1.0e-16;
  if ( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d )
    {
    return_value = numerator / denominator;
    }
  else
    {
    return_value = fabs(numerator) / max_ratio >= fabs(denominator) ?
      ( (numerator >= 0.0 && denominator >= 0.0) ||
        (numerator < 0.0 && denominator < 0.0) ?
        max_ratio : -max_ratio )
      : numerator / denominator;
    }

  return return_value;
}
#endif /* Not currently used and not exposed in verdict.h */

static double safe_ratio( const double numerator, 
    const double denominator )
{

  double return_value;
  const double filter_n = VERDICT_DBL_MAX; 
  const double filter_d = VERDICT_DBL_MIN;
  if ( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d )
    {
    return_value = numerator / denominator;
    }
  else
    {
    return_value = VERDICT_DBL_MAX;
    }

  return return_value;
  
}



static double condition_comp(
  const VerdictVector &xxi, const VerdictVector &xet, const VerdictVector &xze )
{
  double det = VerdictVector::Dot(xxi, (xet * xze));
  
  if ( det <= VERDICT_DBL_MIN) { return VERDICT_DBL_MAX; }
  
  
  double term1 = VerdictVector::Dot(xxi, xxi) + VerdictVector::Dot(xet, xet) + VerdictVector::Dot(xze, xze);
  double term2 = VerdictVector::Dot((xxi*xet), (xxi*xet)) +
                 VerdictVector::Dot((xet*xze), (xet*xze)) +
                 VerdictVector::Dot((xze*xxi), (xze*xxi));
  
  return sqrt( term1 * term2 ) / det;
}



static double oddy_comp( const VerdictVector &xxi, 
    const VerdictVector &xet, 
    const VerdictVector &xze )
{
  
  double g11, g12, g13, g22, g23, g33, rt_g;
  
  g11 = VerdictVector::Dot(xxi, xxi);
  g12 = VerdictVector::Dot(xxi, xet);
  g13 = VerdictVector::Dot(xxi, xze);
  g22 = VerdictVector::Dot(xet, xet);
  g23 = VerdictVector::Dot(xet, xze);
  g33 = VerdictVector::Dot(xze, xze);
  rt_g = VerdictVector::Dot(xxi, (xet * xze));
  
  double oddy_metric;
  if ( rt_g > VERDICT_DBL_MIN ) 
  {
    double norm_G_squared = g11*g11 + 2.0*g12*g12 + 2.0*g13*g13 + g22*g22 + 2.0*g23*g23 +g33*g33;
    
    double norm_J_squared = g11 + g22 + g33;
    
    oddy_metric = ( norm_G_squared - third*norm_J_squared*norm_J_squared ) / pow( rt_g, 4.*third );
    
  }
  
  else
    oddy_metric = VERDICT_DBL_MAX;
  
  return oddy_metric;
  
}


//! calcualates edge lengths of a hex
static double hex_edge_length(int max_min, double coordinates[][3])
{
  double temp[3], edge[12];
  int i;
  
  //lengths^2 of edges
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[1][i] - coordinates[0][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[0] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[2][i] - coordinates[1][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[1] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[3][i] - coordinates[2][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[2] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[0][i] - coordinates[3][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[3] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[5][i] - coordinates[4][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[4] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[6][i] - coordinates[5][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[5] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[7][i] - coordinates[6][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[6] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[4][i] - coordinates[7][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[7] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[4][i] - coordinates[0][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[8] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[5][i] - coordinates[1][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[9] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[6][i] - coordinates[2][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[10] = sqrt( temp[0] + temp[1] + temp[2] );
  
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[7][i] - coordinates[3][i];
    temp[i] = temp[i] * temp[i];
  }
  edge[11] = sqrt( temp[0] + temp[1] + temp[2] );
  
  double _edge = edge[0];
  
  if ( max_min == 0)
  {
    for( i = 1; i<12; i++) 
      _edge = std::min( _edge, edge[i] ); 
    return _edge;
  }  
  else
  {
    for( i = 1; i<12; i++) 
      _edge = std::max( _edge, edge[i] );
    return _edge;
  }
  
}

static double diag_length(int max_min, double coordinates[][3]) 
{
  double temp[3], diag[4];
  int i;
  
  //lengths^2  f diag nals
  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[6][i] - coordinates[0][i];
    temp[i] = temp[i] * temp[i];
  }
  diag[0] = sqrt( temp[0] + temp[1] + temp[2] );

  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[4][i] - coordinates[2][i];
    temp[i] = temp[i] * temp[i];
  }
  diag[1] = sqrt( temp[0] + temp[1] + temp[2] );

  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[7][i] - coordinates[1][i];
    temp[i] = temp[i] * temp[i];
  }
  diag[2] = sqrt( temp[0] + temp[1] + temp[2] );

  for (i = 0; i < 3; i++ )
  {
    temp[i] = coordinates[5][i] - coordinates[3][i];
    temp[i] = temp[i] * temp[i];
  }
  diag[3] = sqrt( temp[0] + temp[1] + temp[2] );

  double diagonal = diag[0];
  if ( max_min == 0 )  //Return min diagonal
  { 
    for( i = 1; i<4; i++)
      diagonal = std::min( diagonal, diag[i] );
    return diagonal;
  }
  else          //Return max diagonal
  {
    for( i = 1; i<4; i++)
      diagonal = std::max( diagonal, diag[i] );
    return diagonal;  
  }

}

//! calculates efg values
static VerdictVector calc_hex_efg( int efg_index, VerdictVector coordinates[8])
{
  
  VerdictVector efg;
  
  switch(efg_index) {
    
    case 1:
      efg =  coordinates[1];
      efg += coordinates[2];
      efg += coordinates[5];
      efg += coordinates[6];
      efg -= coordinates[0];
      efg -= coordinates[3];
      efg -= coordinates[4];
      efg -= coordinates[7];
      break;
      
    case 2:
      efg =  coordinates[2];
      efg += coordinates[3];
      efg += coordinates[6];
      efg += coordinates[7];
      efg -= coordinates[0];
      efg -= coordinates[1];
      efg -= coordinates[4];
      efg -= coordinates[5];
      break;
      
    case 3:
      efg =  coordinates[4];
      efg += coordinates[5];
      efg += coordinates[6];
      efg += coordinates[7];
      efg -= coordinates[0];
      efg -= coordinates[1];
      efg -= coordinates[2];
      efg -= coordinates[3];
      break;

    case 12:
      efg =  coordinates[0];
      efg += coordinates[2];
      efg += coordinates[4];
      efg += coordinates[6];
      efg -= coordinates[1];
      efg -= coordinates[3];
      efg -= coordinates[5];
      efg -= coordinates[7];
      break;
      
    case 13:
      efg =  coordinates[0];
      efg += coordinates[3];
      efg += coordinates[5];
      efg += coordinates[6];
      efg -= coordinates[1];
      efg -= coordinates[2];
      efg -= coordinates[4];
      efg -= coordinates[7];
      break;

   case 23:
      efg =  coordinates[0];
      efg += coordinates[1];
      efg += coordinates[6];
      efg += coordinates[7];
      efg -= coordinates[2];
      efg -= coordinates[3];
      efg -= coordinates[4];
      efg -= coordinates[5];
      break;

   case 123:
      efg =  coordinates[0];
      efg += coordinates[2];
      efg += coordinates[5];
      efg += coordinates[7];
      efg -= coordinates[1];
      efg -= coordinates[5];
      efg -= coordinates[6];
      efg -= coordinates[2];
      break;
      
   default:
      efg.set(0,0,0);
      
  }
  
  return efg;
}

/*!
   the edge ratio of a hex

   NB (P. Pebay 01/23/07): 
     Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
     minimum edge lengths
*/
double hex_edge_ratio (int /*num_nodes*/, double coordinates[][3])
{

  VerdictVector edges[12];
  make_hex_edges(coordinates, edges);
  
  double a2 = edges[0].length_squared();
  double b2 = edges[1].length_squared();
  double c2 = edges[2].length_squared();
  double d2 = edges[3].length_squared();
  double e2 = edges[4].length_squared();
  double f2 = edges[5].length_squared();
  double g2 = edges[6].length_squared();
  double h2 = edges[7].length_squared();
  double i2 = edges[8].length_squared();
  double j2 = edges[9].length_squared();
  double k2 = edges[10].length_squared();
  double l2 = edges[11].length_squared();

  double mab,mcd,mef,Mab,Mcd,Mef;
  double mgh,mij,mkl,Mgh,Mij,Mkl;

  if ( a2 < b2 )
    {
      mab = a2;
      Mab = b2;
    }
  else // b2 <= a2
    {
      mab = b2;
      Mab = a2;
    }
  if ( c2 < d2 )
    {
      mcd = c2;
      Mcd = d2;
    }
  else // d2 <= c2
    {
      mcd = d2;
      Mcd = c2;
    }
  if ( e2 < f2 )
    {
      mef = e2;
      Mef = f2;
    }
  else // f2 <= e2
    {
      mef = f2;
      Mef = e2;
    }
  if ( g2 < h2 )
    {
      mgh = g2;
      Mgh = h2;
    }
  else // h2 <= g2
    {
      mgh = h2;
      Mgh = g2;
    }
  if ( i2 < j2 )
    {
      mij = i2;
      Mij = j2;
    }
  else // j2 <= i2
    {
      mij = j2;
      Mij = i2;
    }
  if ( k2 < l2 )
    {
      mkl = k2;
      Mkl = l2;
    }
  else // l2 <= k2
    {
      mkl = l2;
      Mkl = k2;
    }

  double m2;
  m2 = mab < mcd ? mab : mcd;
  m2 = m2  < mef ? m2  : mef;
  m2 = m2  < mgh ? m2  : mgh;
  m2 = m2  < mij ? m2  : mij;
  m2 = m2  < mkl ? m2  : mkl;

  if ( m2 < VERDICT_DBL_MIN ) 
    return (double)VERDICT_DBL_MAX;

  double M2;
  M2 = Mab > Mcd ? Mab : Mcd;
  M2 = M2  > Mef ? M2  : Mef;
  M2 = M2  > Mgh ? M2  : Mgh;
  M2 = M2  > Mij ? M2  : Mij;
  M2 = M2  > Mkl ? M2  : Mkl;
  m2 = m2  < mef ? m2  : mef;

  M2 = Mab > Mcd ? Mab : Mcd;
  M2 = M2  > Mef ? M2  : Mef;

  double edge_ratio = sqrt( M2 / m2 );
  
  if ( edge_ratio > 0 )
    return (double) std::min( edge_ratio, VERDICT_DBL_MAX );
  return (double) std::max( edge_ratio, -VERDICT_DBL_MAX );
}

/*!
  max edge ratio of a hex

  Maximum edge length ratio at hex center
*/
double hex_max_edge_ratio (int /*num_nodes*/, double coordinates[][3])
{
  double aspect;
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );
  
  double aspect_12, aspect_13, aspect_23;
  
  VerdictVector efg1 = calc_hex_efg( 1, node_pos);
  VerdictVector efg2 = calc_hex_efg( 2, node_pos);
  VerdictVector efg3 = calc_hex_efg( 3, node_pos);

  double mag_efg1 = efg1.length();
  double mag_efg2 = efg2.length();
  double mag_efg3 = efg3.length();
  
  aspect_12 = safe_ratio( std::max( mag_efg1, mag_efg2 ) , std::min( mag_efg1, mag_efg2 ) );
  aspect_13 = safe_ratio( std::max( mag_efg1, mag_efg3 ) , std::min( mag_efg1, mag_efg3 ) );
  aspect_23 = safe_ratio( std::max( mag_efg2, mag_efg3 ) , std::min( mag_efg2, mag_efg3 ) );

  aspect = std::max( {aspect_12, aspect_13, aspect_23 } );

  if ( aspect > 0 )
    return (double) std::min( aspect, VERDICT_DBL_MAX );
  return (double) std::max( aspect, -VERDICT_DBL_MAX );
 
}

double hex_equiangle_skew( int /*num_nodes*/, double coordinates[][3] )
{
  double quad[4][3];
  double min_angle=360.0;
  double max_angle=0.0;
  double min_max_angle[2];

  quad[0][0]=coordinates[0][0];
  quad[0][1]=coordinates[0][1];
  quad[0][2]=coordinates[0][2];
  quad[1][0]=coordinates[1][0];
  quad[1][1]=coordinates[1][1];
  quad[1][2]=coordinates[1][2];
  quad[2][0]=coordinates[5][0];
  quad[2][1]=coordinates[5][1];
  quad[2][2]=coordinates[5][2];
  quad[3][0]=coordinates[4][0];
  quad[3][1]=coordinates[4][1];
  quad[3][2]=coordinates[4][2];

  quad_minimum_maximum_angle(min_max_angle, quad);
  if(min_max_angle[0]<min_angle)
    min_angle=min_max_angle[0];
  if(min_max_angle[1]>max_angle)
    max_angle=min_max_angle[1];


  quad[0][0]=coordinates[1][0];
  quad[0][1]=coordinates[1][1];
  quad[0][2]=coordinates[1][2];
  quad[1][0]=coordinates[2][0];
  quad[1][1]=coordinates[2][1];
  quad[1][2]=coordinates[2][2];
  quad[2][0]=coordinates[6][0];
  quad[2][1]=coordinates[6][1];
  quad[2][2]=coordinates[6][2];
  quad[3][0]=coordinates[5][0];
  quad[3][1]=coordinates[5][1];
  quad[3][2]=coordinates[5][2];

  quad_minimum_maximum_angle(min_max_angle, quad);
  if(min_max_angle[0]<min_angle)
    min_angle=min_max_angle[0];
  if(min_max_angle[1]>max_angle)
    max_angle=min_max_angle[1];


  quad[0][0]=coordinates[2][0];
  quad[0][1]=coordinates[2][1];
  quad[0][2]=coordinates[2][2];
  quad[1][0]=coordinates[3][0];
  quad[1][1]=coordinates[3][1];
  quad[1][2]=coordinates[3][2];
  quad[2][0]=coordinates[7][0];
  quad[2][1]=coordinates[7][1];
  quad[2][2]=coordinates[7][2];
  quad[3][0]=coordinates[6][0];
  quad[3][1]=coordinates[6][1];
  quad[3][2]=coordinates[6][2];

  quad_minimum_maximum_angle(min_max_angle, quad);
  if(min_max_angle[0]<min_angle)
    min_angle=min_max_angle[0];
  if(min_max_angle[1]>max_angle)
    max_angle=min_max_angle[1];


  quad[0][0]=coordinates[3][0];
  quad[0][1]=coordinates[3][1];
  quad[0][2]=coordinates[3][2];
  quad[1][0]=coordinates[0][0];
  quad[1][1]=coordinates[0][1];
  quad[1][2]=coordinates[0][2];
  quad[2][0]=coordinates[4][0];
  quad[2][1]=coordinates[4][1];
  quad[2][2]=coordinates[4][2];
  quad[3][0]=coordinates[7][0];
  quad[3][1]=coordinates[7][1];
  quad[3][2]=coordinates[7][2];

  quad_minimum_maximum_angle(min_max_angle, quad);
  if(min_max_angle[0]<min_angle)
    min_angle=min_max_angle[0];
  if(min_max_angle[1]>max_angle)
    max_angle=min_max_angle[1];


  quad[0][0]=coordinates[4][0];
  quad[0][1]=coordinates[4][1];
  quad[0][2]=coordinates[4][2];
  quad[1][0]=coordinates[5][0];
  quad[1][1]=coordinates[5][1];
  quad[1][2]=coordinates[5][2];
  quad[2][0]=coordinates[6][0];
  quad[2][1]=coordinates[6][1];
  quad[2][2]=coordinates[6][2];
  quad[3][0]=coordinates[7][0];
  quad[3][1]=coordinates[7][1];
  quad[3][2]=coordinates[7][2];

  quad_minimum_maximum_angle(min_max_angle, quad);
  if(min_max_angle[0]<min_angle)
    min_angle=min_max_angle[0];
  if(min_max_angle[1]>max_angle)
    max_angle=min_max_angle[1];


  quad[0][0]=coordinates[3][0];
  quad[0][1]=coordinates[3][1];
  quad[0][2]=coordinates[3][2];
  quad[1][0]=coordinates[2][0];
  quad[1][1]=coordinates[2][1];
  quad[1][2]=coordinates[2][2];
  quad[2][0]=coordinates[1][0];
  quad[2][1]=coordinates[1][1];
  quad[2][2]=coordinates[1][2];
  quad[3][0]=coordinates[0][0];
  quad[3][1]=coordinates[0][1];
  quad[3][2]=coordinates[0][2];

  quad_minimum_maximum_angle(min_max_angle, quad);
  if(min_max_angle[0]<min_angle)
    min_angle=min_max_angle[0];
  if(min_max_angle[1]>max_angle)
    max_angle=min_max_angle[1];



  double skew_max= (max_angle-90.0)/90.0;
  double skew_min= (90.0-min_angle)/90.0;

  if(skew_max>skew_min)
    return skew_max;
  return skew_min;


}

/*!
  skew of a hex

  Maximum ||cosA|| where A is the angle between edges at hex center.
*/
double hex_skew( int /*num_nodes*/, double coordinates[][3] )
{
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );
  
  double skew_1, skew_2, skew_3;
  
  VerdictVector efg1 = calc_hex_efg( 1, node_pos);
  VerdictVector efg2 = calc_hex_efg( 2, node_pos);
  VerdictVector efg3 = calc_hex_efg( 3, node_pos);
 
  if ( efg1.normalize() <= VERDICT_DBL_MIN )
    return VERDICT_DBL_MAX;
  if ( efg2.normalize() <= VERDICT_DBL_MIN )
    return VERDICT_DBL_MAX;
  if ( efg3.normalize() <= VERDICT_DBL_MIN )
    return VERDICT_DBL_MAX;

  skew_1 = fabs(VerdictVector::Dot(efg1, efg2));
  skew_2 = fabs(VerdictVector::Dot(efg1, efg3));
  skew_3 = fabs(VerdictVector::Dot(efg2, efg3));

  double skew = std::max( {skew_1, skew_2, skew_3} );

  if ( skew > 0 )
    return (double) std::min( skew, VERDICT_DBL_MAX );
  return (double) std::max( skew, -VERDICT_DBL_MAX );
}

/*!
  taper of a hex

  Maximum ratio of lengths derived from opposite edges.
*/
double hex_taper( int /*num_nodes*/, double coordinates[][3] )
{
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );
  
  VerdictVector efg1 = calc_hex_efg( 1, node_pos);
  VerdictVector efg2 = calc_hex_efg( 2, node_pos);
  VerdictVector efg3 = calc_hex_efg( 3, node_pos);

  VerdictVector efg12 = calc_hex_efg( 12, node_pos);
  VerdictVector efg13 = calc_hex_efg( 13, node_pos);
  VerdictVector efg23 = calc_hex_efg( 23, node_pos);

  double taper_1 = fabs( safe_ratio( efg12.length(), std::min( efg1.length(), efg2.length())));
  double taper_2 = fabs( safe_ratio( efg13.length(), std::min( efg1.length(), efg3.length())));
  double taper_3 = fabs( safe_ratio( efg23.length(), std::min( efg2.length(), efg3.length())));

  double taper = std::max( {taper_1, taper_2, taper_3} );

  if ( taper > 0 )
    return (double) std::min( taper, VERDICT_DBL_MAX );
  return (double) std::max( taper, -VERDICT_DBL_MAX );
}

/*!
  volume of a hex
  Split the hex into 24 tets.
  sum the volume of each tet.
*/
double hex_volume( int /*num_nodes*/, double coordinates[][3] )
{

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  //define the nodes of each face of the hex
  int faces[6][4]=
  {
    {0,1,5,4},
    {1,2,6,5},
    {2,3,7,6},
    {3,0,4,7},
    {3,2,1,0},
    {4,5,6,7},
  };

  //calculate the center of each face
  VerdictVector fcenter[6];
  for(int f=0;f<6;f++)
  {
   fcenter[f]=(node_pos[faces[f][0]] + node_pos[faces[f][1]] + node_pos[faces[f][2]] + node_pos[faces[f][3]]) * 0.25;
  }

  //calculate the center of the hex
  VerdictVector hcenter = (node_pos[0] + node_pos[1] + node_pos[2] + node_pos[3] +
                         node_pos[4] + node_pos[5] + node_pos[6] + node_pos[7]) * 0.125;


  double volume=0.0;
  for(int i=0;i<6;i++)
  {
    //for each face calculate the vectors from the nodes and center of the face to the center of the hex.
    //These vectors define three of the sides of the tets.
    VerdictVector side[5];
    side[4]=hcenter-fcenter[i];//vector from center of face to center of hex.
    for(int s=0;s<4;s++)
    {
      side[s]=hcenter-node_pos[faces[i][s]];//vector from face node to center of hex.
    }

    //for each of the four tets that originate from this face.
    //calculate the volume of the tet.
    //This is done by calculating the triple product of three vectors that originate from a corner node of the tet.
    //This is also the jacobain at the corner node of the tet.
    //The volume is 1/6 of jacobian at a corner node.
    for(int j=0;j<3;j++)//first three tets
    {
      volume += (double)(VerdictVector::Dot(side[4], (side[j + 1] * side[j])) / 6.0);
    }
    volume += (double)(VerdictVector::Dot(side[4], (side[0] * side[3])) / 6.0);//fourth tet.
  }

  if ( volume > 0 )
    return (double) std::min( volume, VERDICT_DBL_MAX );
  return (double) std::max( volume, -VERDICT_DBL_MAX );
}

/*!
  stretch of a hex

  sqrt(3) * minimum edge length / maximum diagonal length
*/
double hex_stretch( int /*num_nodes*/, double coordinates[][3] )
{
  double min_edge = hex_edge_length( 0, coordinates );
  double max_diag = diag_length( 1, coordinates );  
  
  double stretch = root_of_3 * safe_ratio(min_edge, max_diag);

  if ( stretch > 0 )
    return (double) std::min( stretch, VERDICT_DBL_MAX );
  return (double) std::max( stretch, -VERDICT_DBL_MAX );
}

/*!
  diagonal ratio of a hex
  
  Minimum diagonal length / maximum diagonal length
*/
double hex_diagonal( int /*num_nodes*/, double coordinates[][3] )
{
  
  double min_diag = diag_length( 0, coordinates ); 
  double max_diag = diag_length( 1, coordinates );
  
  double diagonal = safe_ratio( min_diag, max_diag);

  if ( diagonal > 0 )
    return (double) std::min( diagonal, VERDICT_DBL_MAX );
  return (double) std::max( diagonal, -VERDICT_DBL_MAX );
}

#define SQR(x) ((x) * (x))

/*!
  dimension of a hex

  Pronto-specific characteristic length for stable time step calculation. 
  Char_length = Volume / 2 grad Volume
*/
double hex_dimension( int /*num_nodes*/, double coordinates[][3] )
{
  
  double gradop[9][4];

  double x1 = coordinates[0][0];
  double x2 = coordinates[1][0];
  double x3 = coordinates[2][0];
  double x4 = coordinates[3][0];
  double x5 = coordinates[4][0];
  double x6 = coordinates[5][0];
  double x7 = coordinates[6][0];
  double x8 = coordinates[7][0];

  double y1 = coordinates[0][1];
  double y2 = coordinates[1][1];
  double y3 = coordinates[2][1];
  double y4 = coordinates[3][1];
  double y5 = coordinates[4][1];
  double y6 = coordinates[5][1];
  double y7 = coordinates[6][1];
  double y8 = coordinates[7][1];

  double z1 = coordinates[0][2];
  double z2 = coordinates[1][2];
  double z3 = coordinates[2][2];
  double z4 = coordinates[3][2];
  double z5 = coordinates[4][2];
  double z6 = coordinates[5][2];
  double z7 = coordinates[6][2];
  double z8 = coordinates[7][2];

  double z24 = z2 - z4;
  double z52 = z5 - z2;
  double z45 = z4 - z5;
  gradop[1][1] = ( y2*(z6-z3-z45) + y3*z24 + y4*(z3-z8-z52)
      + y5*(z8-z6-z24) + y6*z52 + y8*z45 ) / 12.0;
         
  double z31 = z3 - z1;
  double z63 = z6 - z3;
  double z16 = z1 - z6;
  gradop[2][1] = ( y3*(z7-z4-z16) + y4*z31 + y1*(z4-z5-z63)
      + y6*(z5-z7-z31) + y7*z63 + y5*z16 ) / 12.0;

  double z42 = z4 - z2;
  double z74 = z7 - z4;
  double z27 = z2 - z7;
  gradop[3][1] = ( y4*(z8-z1-z27) + y1*z42 + y2*(z1-z6-z74)
      + y7*(z6-z8-z42) + y8*z74 + y6*z27 ) / 12.0;
             
  double z13 = z1 - z3;
  double z81 = z8 - z1;
  double z38 = z3 - z8;
  gradop[4][1] = ( y1*(z5-z2-z38) + y2*z13 + y3*(z2-z7-z81)
      + y8*(z7-z5-z13) + y5*z81 + y7*z38 ) / 12.0;
   
  double z86 = z8 - z6;
  double z18 = z1 - z8;
  double z61 = z6 - z1;
  gradop[5][1] = ( y8*(z4-z7-z61) + y7*z86 + y6*(z7-z2-z18)
      + y1*(z2-z4-z86) + y4*z18 + y2*z61 ) / 12.0;
     
  double z57 = z5 - z7;
  double z25 = z2 - z5;
  double z72 = z7 - z2;
  gradop[6][1] = ( y5*(z1-z8-z72) + y8*z57 + y7*(z8-z3-z25)
      + y2*(z3-z1-z57) + y1*z25 + y3*z72 ) / 12.0;
       
  double z68 = z6 - z8;
  double z36 = z3 - z6;
  double z83 = z8 - z3;
  gradop[7][1] = ( y6*(z2-z5-z83) + y5*z68 + y8*(z5-z4-z36)
      + y3*(z4-z2-z68) + y2*z36 + y4*z83 ) / 12.0;
         
  double z75 = z7 - z5;
  double z47 = z4 - z7;
  double z54 = z5 - z4;
  gradop[8][1] = ( y7*(z3-z6-z54) + y6*z75 + y5*(z6-z1-z47)
      + y4*(z1-z3-z75) + y3*z47 + y1*z54 ) / 12.0;
           
  double x24 = x2 - x4;
  double x52 = x5 - x2;
  double x45 = x4 - x5;
  gradop[1][2] = ( z2*(x6-x3-x45) + z3*x24 + z4*(x3-x8-x52)
      + z5*(x8-x6-x24) + z6*x52 + z8*x45 ) / 12.0;

  double x31 = x3 - x1;
  double x63 = x6 - x3;
  double x16 = x1 - x6;
  gradop[2][2] = ( z3*(x7-x4-x16) + z4*x31 + z1*(x4-x5-x63)
      + z6*(x5-x7-x31) + z7*x63 + z5*x16 ) / 12.0;

  double x42 = x4 - x2;
  double x74 = x7 - x4;
  double x27 = x2 - x7;
  gradop[3][2] = ( z4*(x8-x1-x27) + z1*x42 + z2*(x1-x6-x74)
      + z7*(x6-x8-x42) + z8*x74 + z6*x27 ) / 12.0;

  double x13 = x1 - x3;
  double x81 = x8 - x1;
  double x38 = x3 - x8;
  gradop[4][2] = ( z1*(x5-x2-x38) + z2*x13 + z3*(x2-x7-x81)
      + z8*(x7-x5-x13) + z5*x81 + z7*x38 ) / 12.0;

  double x86 = x8 - x6;
  double x18 = x1 - x8;
  double x61 = x6 - x1;
  gradop[5][2] = ( z8*(x4-x7-x61) + z7*x86 + z6*(x7-x2-x18)
      + z1*(x2-x4-x86) + z4*x18 + z2*x61 ) / 12.0;
                     
  double x57 = x5 - x7;
  double x25 = x2 - x5;
  double x72 = x7 - x2;
  gradop[6][2] = ( z5*(x1-x8-x72) + z8*x57 + z7*(x8-x3-x25)
      + z2*(x3-x1-x57) + z1*x25 + z3*x72 ) / 12.0;

  double x68 = x6 - x8;
  double x36 = x3 - x6;
  double x83 = x8 - x3;
  gradop[7][2] = ( z6*(x2-x5-x83) + z5*x68 + z8*(x5-x4-x36)
      + z3*(x4-x2-x68) + z2*x36 + z4*x83 ) / 12.0;

  double x75 = x7 - x5;
  double x47 = x4 - x7;
  double x54 = x5 - x4;
  gradop[8][2] = ( z7*(x3-x6-x54) + z6*x75 + z5*(x6-x1-x47)
      + z4*(x1-x3-x75) + z3*x47 + z1*x54 ) / 12.0;
           
  double y24 = y2 - y4;
  double y52 = y5 - y2;
  double y45 = y4 - y5;
  gradop[1][3] = ( x2*(y6-y3-y45) + x3*y24 + x4*(y3-y8-y52)
      + x5*(y8-y6-y24) + x6*y52 + x8*y45 ) / 12.0;
             
  double y31 = y3 - y1;
  double y63 = y6 - y3;
  double y16 = y1 - y6;
  gradop[2][3] = ( x3*(y7-y4-y16) + x4*y31 + x1*(y4-y5-y63)
      + x6*(y5-y7-y31) + x7*y63 + x5*y16 ) / 12.0;
               
  double y42 = y4 - y2;
  double y74 = y7 - y4;
  double y27 = y2 - y7;
  gradop[3][3] = ( x4*(y8-y1-y27) + x1*y42 + x2*(y1-y6-y74)
      + x7*(y6-y8-y42) + x8*y74 + x6*y27 ) / 12.0;
                 
  double y13 = y1 - y3;
  double y81 = y8 - y1;
  double y38 = y3 - y8;
  gradop[4][3] = ( x1*(y5-y2-y38) + x2*y13 + x3*(y2-y7-y81)
      + x8*(y7-y5-y13) + x5*y81 + x7*y38 ) / 12.0;
                   
  double y86 = y8 - y6;
  double y18 = y1 - y8;
  double y61 = y6 - y1;
  gradop[5][3] = ( x8*(y4-y7-y61) + x7*y86 + x6*(y7-y2-y18)
      + x1*(y2-y4-y86) + x4*y18 + x2*y61 ) / 12.0;
                     
  double y57 = y5 - y7;
  double y25 = y2 - y5;
  double y72 = y7 - y2;
  gradop[6][3] = ( x5*(y1-y8-y72) + x8*y57 + x7*(y8-y3-y25)
      + x2*(y3-y1-y57) + x1*y25 + x3*y72 ) / 12.0;
                       
  double y68 = y6 - y8;
  double y36 = y3 - y6;
  double y83 = y8 - y3;
  gradop[7][3] = ( x6*(y2-y5-y83) + x5*y68 + x8*(y5-y4-y36)
      + x3*(y4-y2-y68) + x2*y36 + x4*y83 ) / 12.0;
                         
  double y75 = y7 - y5;
  double y47 = y4 - y7;
  double y54 = y5 - y4;
  gradop[8][3] = ( x7*(y3-y6-y54) + x6*y75 + x5*(y6-y1-y47)
      + x4*(y1-y3-y75) + x3*y47 + x1*y54 ) / 12.0;

  //     calculate element volume and characteristic element aspect ratio
  //     (used in time step and hourglass control) - 



  double volume =  coordinates[0][0] * gradop[1][1]
    + coordinates[1][0] * gradop[2][1]
    + coordinates[2][0] * gradop[3][1]
    + coordinates[3][0] * gradop[4][1]
    + coordinates[4][0] * gradop[5][1]
    + coordinates[5][0] * gradop[6][1]
    + coordinates[6][0] * gradop[7][1]
    + coordinates[7][0] * gradop[8][1];
  double aspect = .5*SQR(volume) /
    ( SQR(gradop[1][1]) + SQR(gradop[2][1])
      + SQR(gradop[3][1]) + SQR(gradop[4][1])
      + SQR(gradop[5][1]) + SQR(gradop[6][1])
      + SQR(gradop[7][1]) + SQR(gradop[8][1])
      + SQR(gradop[1][2]) + SQR(gradop[2][2])
      + SQR(gradop[3][2]) + SQR(gradop[4][2])
      + SQR(gradop[5][2]) + SQR(gradop[6][2])
      + SQR(gradop[7][2]) + SQR(gradop[8][2])
      + SQR(gradop[1][3]) + SQR(gradop[2][3])
      + SQR(gradop[3][3]) + SQR(gradop[4][3])
      + SQR(gradop[5][3]) + SQR(gradop[6][3])
      + SQR(gradop[7][3]) + SQR(gradop[8][3]) );
     
  return (double)sqrt(aspect);
  
}

/*!
  oddy of a hex

  General distortion measure based on left Cauchy-Green Tensor
*/
double hex_oddy( int /*num_nodes*/, double coordinates[][3] )
{
  
  double oddy = 0.0, current_oddy;
  VerdictVector xxi, xet, xze;
  
  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );
  
  xxi = calc_hex_efg(1, node_pos);
  xet = calc_hex_efg(2, node_pos);
  xze = calc_hex_efg(3, node_pos);
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[1][0] - coordinates[0][0],
          coordinates[1][1] - coordinates[0][1],
          coordinates[1][2] - coordinates[0][2] );
  
  xet.set(coordinates[3][0] - coordinates[0][0],
          coordinates[3][1] - coordinates[0][1],
          coordinates[3][2] - coordinates[0][2] );
  
  xze.set(coordinates[4][0] - coordinates[0][0],
          coordinates[4][1] - coordinates[0][1],
          coordinates[4][2] - coordinates[0][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[2][0] - coordinates[1][0],
          coordinates[2][1] - coordinates[1][1],
          coordinates[2][2] - coordinates[1][2] );
  
  xet.set(coordinates[0][0] - coordinates[1][0],
          coordinates[0][1] - coordinates[1][1],
          coordinates[0][2] - coordinates[1][2] );
  
  xze.set(coordinates[5][0] - coordinates[1][0],
          coordinates[5][1] - coordinates[1][1],
          coordinates[5][2] - coordinates[1][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[3][0] - coordinates[2][0],
          coordinates[3][1] - coordinates[2][1],
          coordinates[3][2] - coordinates[2][2] );
  
  xet.set(coordinates[1][0] - coordinates[2][0],
          coordinates[1][1] - coordinates[2][1],
          coordinates[1][2] - coordinates[2][2] );
  
  xze.set(coordinates[6][0] - coordinates[2][0],
          coordinates[6][1] - coordinates[2][1],
          coordinates[6][2] - coordinates[2][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[0][0] - coordinates[3][0],
          coordinates[0][1] - coordinates[3][1],
          coordinates[0][2] - coordinates[3][2] );
  
  xet.set(coordinates[2][0] - coordinates[3][0],
          coordinates[2][1] - coordinates[3][1],
          coordinates[2][2] - coordinates[3][2] );
  
  xze.set(coordinates[7][0] - coordinates[3][0],
          coordinates[7][1] - coordinates[3][1],
          coordinates[7][2] - coordinates[3][2] );
 
 
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[7][0] - coordinates[4][0],
          coordinates[7][1] - coordinates[4][1],
          coordinates[7][2] - coordinates[4][2] );
  
  xet.set(coordinates[5][0] - coordinates[4][0],
          coordinates[5][1] - coordinates[4][1],
          coordinates[5][2] - coordinates[4][2] );
  
  xze.set(coordinates[0][0] - coordinates[4][0],
          coordinates[0][1] - coordinates[4][1],
          coordinates[0][2] - coordinates[4][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[4][0] - coordinates[5][0],
          coordinates[4][1] - coordinates[5][1],
          coordinates[4][2] - coordinates[5][2] );
  
  xet.set(coordinates[6][0] - coordinates[5][0],
          coordinates[6][1] - coordinates[5][1],
          coordinates[6][2] - coordinates[5][2] );
  
  xze.set(coordinates[1][0] - coordinates[5][0],
          coordinates[1][1] - coordinates[5][1],
          coordinates[1][2] - coordinates[5][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[5][0] - coordinates[6][0],
          coordinates[5][1] - coordinates[6][1],
          coordinates[5][2] - coordinates[6][2] );
  
  xet.set(coordinates[7][0] - coordinates[6][0],
          coordinates[7][1] - coordinates[6][1],
          coordinates[7][2] - coordinates[6][2] );
  
  xze.set(coordinates[2][0] - coordinates[6][0],
          coordinates[2][1] - coordinates[6][1],
          coordinates[2][2] - coordinates[6][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  xxi.set(coordinates[6][0] - coordinates[7][0],
          coordinates[6][1] - coordinates[7][1],
          coordinates[6][2] - coordinates[7][2] );
  
  xet.set(coordinates[4][0] - coordinates[7][0],
          coordinates[4][1] - coordinates[7][1],
          coordinates[4][2] - coordinates[7][2] );
  
  xze.set(coordinates[3][0] - coordinates[7][0],
          coordinates[3][1] - coordinates[7][1],
          coordinates[3][2] - coordinates[7][2] );
  
  
  current_oddy = oddy_comp( xxi, xet, xze);
  if ( current_oddy > oddy ) { oddy = current_oddy; }
  
  if ( oddy > 0 )
    return (double) std::min( oddy, VERDICT_DBL_MAX );
  return (double) std::max( oddy, -VERDICT_DBL_MAX );

}

/*!
   the average Frobenius aspect of a hex

   NB (P. Pebay 01/20/07): 
     this function is calculated by averaging the 8 Frobenius aspects at
     each corner of the hex, when the reference corner is right isosceles.
*/
double hex_med_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
{

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  VerdictVector xxi, xet, xze;

  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  double med_aspect_frobenius = condition_comp( xxi, xet, xze );

  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  med_aspect_frobenius += condition_comp( xxi, xet, xze );

  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  med_aspect_frobenius += condition_comp( xxi, xet, xze );

  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  med_aspect_frobenius += condition_comp( xxi, xet, xze );

  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  med_aspect_frobenius += condition_comp( xxi, xet, xze );

  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  med_aspect_frobenius += condition_comp( xxi, xet, xze );

  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  med_aspect_frobenius += condition_comp( xxi, xet, xze );

  // J(1,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  med_aspect_frobenius += condition_comp( xxi, xet, xze );
  
  if (med_aspect_frobenius >=  VERDICT_DBL_MAX)  return  VERDICT_DBL_MAX;
  if (med_aspect_frobenius <= -VERDICT_DBL_MAX)  return -VERDICT_DBL_MAX;

  return med_aspect_frobenius / 24.;
}

/*!
  maximum Frobenius condition number of a hex

  Maximum Frobenius condition number of the Jacobian matrix at 8 corners
   NB (P. Pebay 01/25/07): 
     this function is calculated by taking the maximum of the 8 Frobenius aspects at
     each corner of the hex, when the reference corner is right isosceles.
*/
double hex_max_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
{

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  VerdictVector xxi, xet, xze;

  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  double condition = condition_comp( xxi, xet, xze );

  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  double current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  // J(1,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  current_condition = condition_comp( xxi, xet, xze );
  if ( current_condition > condition ) { condition = current_condition; }

  if (condition >=  VERDICT_DBL_MAX)  return  VERDICT_DBL_MAX;
  if (condition <= -VERDICT_DBL_MAX)  return -VERDICT_DBL_MAX;
  
  return condition / 3.;
}

/*!
  The maximum Frobenius condition of a hex, a.k.a. condition
  NB (P. Pebay 01/25/07): 
     this method is maintained for backwards compatibility only.
     It will become deprecated at some point.

*/
double hex_condition( int /*num_nodes*/, double coordinates[][3] )
{

  return hex_max_aspect_frobenius(8, coordinates);
}

/*!
  jacobian of a hex

  Minimum pointwise volume of local map at 8 corners & center of hex
*/
double hex_jacobian( int num_nodes, double coordinates[][3] )
{
  if(num_nodes == 27)
  {
    double dhdr[27];
    double dhds[27];
    double dhdt[27];
    double min_determinant = VERDICT_DBL_MAX;

    for(int i=0; i<27; i++)
    {
      HEX27_gradients_of_the_shape_functions_for_RST(HEX27_node_local_coord[i], dhdr, dhds, dhdt);
      double jacobian[3][3] = {{0,0,0},{0,0,0},{0,0,0}};

      for(int j=0; j<27; j++)
      {
        jacobian[0][0]+=coordinates[j][0]*dhdr[j];
        jacobian[0][1]+=coordinates[j][0]*dhds[j];
        jacobian[0][2]+=coordinates[j][0]*dhdt[j];
        jacobian[1][0]+=coordinates[j][1]*dhdr[j];
        jacobian[1][1]+=coordinates[j][1]*dhds[j];
        jacobian[1][2]+=coordinates[j][1]*dhdt[j];
        jacobian[2][0]+=coordinates[j][2]*dhdr[j];
        jacobian[2][1]+=coordinates[j][2]*dhds[j];
        jacobian[2][2]+=coordinates[j][2]*dhdt[j];
      }
      double det = VerdictVector::Dot((VerdictVector(jacobian[0]) * VerdictVector(jacobian[1])), VerdictVector(jacobian[2]));
      min_determinant = std::min(det, min_determinant);
    }
    return min_determinant;
  }
  else
  {
    VerdictVector node_pos[8];
    make_hex_nodes ( coordinates, node_pos );

    double jacobian = VERDICT_DBL_MAX;
    double current_jacobian;
    VerdictVector xxi, xet, xze;

    xxi = calc_hex_efg(1, node_pos );
    xet = calc_hex_efg(2, node_pos );
    xze = calc_hex_efg(3, node_pos );


    current_jacobian = VerdictVector::Dot(xxi, (xet * xze)) / 64.0;
    if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(0,0,0):

    xxi = node_pos[1] - node_pos[0];
    xet = node_pos[3] - node_pos[0];
    xze = node_pos[4] - node_pos[0];

    current_jacobian = VerdictVector::Dot(xxi, (xet * xze));
    if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(1,0,0):

    xxi = node_pos[2] - node_pos[1];
    xet = node_pos[0] - node_pos[1];
    xze = node_pos[5] - node_pos[1];

    current_jacobian = VerdictVector::Dot(xxi, (xet * xze));
    if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(1,1,0):

    xxi = node_pos[3] - node_pos[2];
    xet = node_pos[1] - node_pos[2];
    xze = node_pos[6] - node_pos[2];

    current_jacobian = VerdictVector::Dot(xxi, (xet * xze));
    if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(0,1,0):

    xxi = node_pos[0] - node_pos[3];
    xet = node_pos[2] - node_pos[3];
    xze = node_pos[7] - node_pos[3];

    current_jacobian = VerdictVector::Dot(xxi, (xet * xze));
    if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(0,0,1):

    xxi = node_pos[7] - node_pos[4];
    xet = node_pos[5] - node_pos[4];
    xze = node_pos[0] - node_pos[4];

    current_jacobian = xxi % (xet * xze);
    if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(1,0,1):

    xxi = node_pos[4] - node_pos[5];
    xet = node_pos[6] - node_pos[5];
    xze = node_pos[1] - node_pos[5];

    current_jacobian = xxi % (xet * xze);
    if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(1,1,1):

    xxi = node_pos[5] - node_pos[6];
    xet = node_pos[7] - node_pos[6];
    xze = node_pos[2] - node_pos[6];

    current_jacobian = xxi % (xet * xze);
    if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(0,1,1):

    xxi = node_pos[6] - node_pos[7];
    xet = node_pos[4] - node_pos[7];
    xze = node_pos[3] - node_pos[7];

    current_jacobian = xxi % (xet * xze);
    if ( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    if ( jacobian > 0 )
      return (double) std::min( jacobian, VERDICT_DBL_MAX );
    return (double) std::max( jacobian, -VERDICT_DBL_MAX );
  }
}

/*!
  scaled jacobian of a hex

  Minimum Jacobian divided by the lengths of the 3 edge vectors
*/
double hex_scaled_jacobian( int num_nodes, double coordinates[][3] )
{

  double jacobi, min_norm_jac = VERDICT_DBL_MAX;

#if 0
  if(num_nodes == 27)
  {
    double dhdr[27];
    double dhds[27];
    double dhdt[27];

    for(int i=0; i<27; i++)
    {
      HEX27_gradients_of_the_shape_functions_for_RST(HEX27_node_local_coord[i], dhdr, dhds, dhdt);
      double jacobian[3][3] = {{0,0,0},{0,0,0},{0,0,0}};

      for(int j=0; j<27; j++)
      {
        jacobian[0][0]+=coordinates[j][0]*dhdr[j];
        jacobian[1][0]+=coordinates[j][0]*dhds[j];
        jacobian[2][0]+=coordinates[j][0]*dhdt[j];
        jacobian[0][1]+=coordinates[j][1]*dhdr[j];
        jacobian[1][1]+=coordinates[j][1]*dhds[j];
        jacobian[2][1]+=coordinates[j][1]*dhdt[j];
        jacobian[0][2]+=coordinates[j][2]*dhdr[j];
        jacobian[1][2]+=coordinates[j][2]*dhds[j];
        jacobian[2][2]+=coordinates[j][2]*dhdt[j];
      }
      VerdictVector xxi(jacobian[0]);
      VerdictVector xet(jacobian[1]);
      VerdictVector xze(jacobian[2]);
      jacobi = xxi % ( xet * xze );
      double scale = sqrt(xxi.length_squared() * xet.length_squared() * xze.length_squared());
      jacobi /= scale;
      if(jacobi < min_norm_jac)
      {
        min_norm_jac = jacobi;
      }
    }
  }
  else
#endif
  {

  double min_jacobi = VERDICT_DBL_MAX;
  double temp_norm_jac, lengths;
  double len1_sq, len2_sq, len3_sq; 
  VerdictVector xxi, xet, xze;

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  xxi = calc_hex_efg(1, node_pos );
  xet = calc_hex_efg(2, node_pos );
  xze = calc_hex_efg(3, node_pos );

  jacobi = VerdictVector::Dot(xxi, (xet * xze));
  if ( jacobi < min_jacobi) { min_jacobi = jacobi; }

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (double) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  
  if ( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;

  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  jacobi = VerdictVector::Dot(xxi, (xet * xze));
  if ( jacobi < min_jacobi ) { min_jacobi = jacobi; }

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (double) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if ( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;


  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  jacobi = VerdictVector::Dot(xxi, (xet * xze));
  if ( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (double) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if ( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;

  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  jacobi = VerdictVector::Dot(xxi, (xet * xze));
  if ( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (double) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if ( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;

  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  jacobi = VerdictVector::Dot(xxi, (xet * xze));
  if ( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (double) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if ( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else  
    temp_norm_jac = jacobi;

  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  jacobi = VerdictVector::Dot(xxi, (xet * xze));
  if ( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (double) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if ( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else 
    temp_norm_jac = jacobi;

  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  jacobi = VerdictVector::Dot(xxi, (xet * xze));
  if ( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (double) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if ( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else
    temp_norm_jac = jacobi;

  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  jacobi = VerdictVector::Dot(xxi, (xet * xze));
  if ( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (double) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if ( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else
    temp_norm_jac = jacobi;

  // J(0,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  jacobi = VerdictVector::Dot(xxi, (xet * xze));
  if ( jacobi < min_jacobi ) { min_jacobi = jacobi; }
  
  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return (double) VERDICT_DBL_MAX;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  temp_norm_jac = jacobi / lengths;
  if ( temp_norm_jac < min_norm_jac)
    min_norm_jac = temp_norm_jac;  
  else 
    temp_norm_jac = jacobi;
  }

  if ( min_norm_jac> 0 )
    return (double) std::min( min_norm_jac, VERDICT_DBL_MAX );
  return (double) std::max( min_norm_jac, -VERDICT_DBL_MAX );
}

/*!
  Nodal jacobian ratio of a hex
  Minimum nodal jacobian divided by the maximum.  Detects element skewness.
*/

inline std::pair<double /*min*/, double /*max*/> minMaxVal(const double a1, const double a2)
{
  return (a1 < a2) ? std::pair<double, double>(a1, a2) : std::pair<double, double>(a2, a1);
}

inline std::pair<double, double> minMaxValPair(const std::pair<double, double>& a1, const std::pair<double, double>& a2)
{
  return std::pair<double, double>(std::min(a1.first, a2.first), std::max(a1.second, a2.second));
}

double hex_nodal_jacobian_ratio( int num_nodes, double coordinates[][3])
{
  return verdict::hex_nodal_jacobian_ratio2( num_nodes, (double*)coordinates);
}

double hex_nodal_jacobian_ratio2( int /*num_nodes*/, double *coordinates)
{
  double Jdet8x[8];
  verdict::hex_nodal_jacobians(coordinates, Jdet8x);
  //
  //  Compute the minimum and maximum nodal determinates, use an optimal algorithm
  //
  // std::pair<double, double> minMaxResult = std::minmax_element(Jdet8x, Jdet8x+8);
 
  std::pair<double, double> m01(minMaxVal(Jdet8x[0], Jdet8x[1]));
  std::pair<double, double> m23(minMaxVal(Jdet8x[2], Jdet8x[3]));
  std::pair<double, double> m45(minMaxVal(Jdet8x[4], Jdet8x[5]));
  std::pair<double, double> m67(minMaxVal(Jdet8x[6], Jdet8x[7]));

  std::pair<double, double> m0123(minMaxValPair(m01, m23));
  std::pair<double, double> m4567(minMaxValPair(m45, m67));

  std::pair<double, double> m01234567(minMaxValPair(m0123, m4567));


  //
  //  Turn the determinates into a normalized quality ratio.
  //  If the maximum nodal jacobian is negative the element is fully inverted, return huge negative number
  //  Otherwise return ratio of the minimal nodal determinate to the maximum
  //

  if (m01234567.second <= VERDICT_DBL_MIN) {
    return -VERDICT_DBL_MAX;
  } else {
    return m01234567.first / m01234567.second;
  }

}

/*!
  shear of a hex
  
  3/Condition number of Jacobian Skew matrix
*/
double hex_shear( int /*num_nodes*/, double coordinates[][3] )
{

  double shear;
  double min_shear = 1.0; 
  VerdictVector xxi, xet, xze;
  double det, len1_sq, len2_sq, len3_sq, lengths;


  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths;
  min_shear = std::min( shear, min_shear );


  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );


  // J(0,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  len1_sq = xxi.length_squared();
  len2_sq = xet.length_squared();
  len3_sq = xze.length_squared();

  if ( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN ||
      len3_sq <= VERDICT_DBL_MIN)
    return 0;

  lengths = sqrt( len1_sq * len2_sq * len3_sq );
  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det < VERDICT_DBL_MIN ) { return 0; }

  shear = det / lengths; 
  min_shear = std::min( shear, min_shear );

  if ( min_shear <= VERDICT_DBL_MIN )
    min_shear = 0;
  
  if ( min_shear > 0 )
    return (double) std::min( min_shear, VERDICT_DBL_MAX );
  return (double) std::max( min_shear, -VERDICT_DBL_MAX );
}

/*!
  shape of a hex

  3/Condition number of weighted Jacobian matrix
*/
double hex_shape( int /*num_nodes*/, double coordinates[][3] )
{


  double det, shape;
  double min_shape = 1.0; 

  VerdictVector xxi, xet, xze;

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );


  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];

  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow(det, two_thirds) / (VerdictVector::Dot(xxi,xxi) + VerdictVector::Dot(xet,xet) + VerdictVector::Dot(xze,xze));
  else
    return 0;
  
  if ( shape < min_shape ) { min_shape = shape; }
  

  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow(det, two_thirds) / (VerdictVector::Dot(xxi,xxi) + VerdictVector::Dot(xet,xet) + VerdictVector::Dot(xze,xze));
  else
    return 0;
  
  if ( shape < min_shape ) { min_shape = shape; }
  

  // J(1,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow(det, two_thirds) / (VerdictVector::Dot(xxi,xxi) + VerdictVector::Dot(xet,xet) + VerdictVector::Dot(xze,xze));
  else
    return 0;
  
  if ( shape < min_shape ) { min_shape = shape; }
  

  // J(0,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow(det, two_thirds) / (VerdictVector::Dot(xxi, xxi) + VerdictVector::Dot(xet, xet) + VerdictVector::Dot(xze,xze));
  else
    return 0;
  
  if ( shape < min_shape ) { min_shape = shape; }
  

  // J(0,0,1):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow(det, two_thirds) / (VerdictVector::Dot(xxi,xxi) + VerdictVector::Dot(xet,xet) + VerdictVector::Dot(xze,xze));
  else
    return 0;
  
  if ( shape < min_shape ) { min_shape = shape; }
  

  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow(det, two_thirds) / (VerdictVector::Dot(xxi,xxi) + VerdictVector::Dot(xet,xet) + VerdictVector::Dot(xze,xze));
  else
    return 0;
  
  if ( shape < min_shape ) { min_shape = shape; }
  

  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow(det, two_thirds) / (VerdictVector::Dot(xxi, xxi) + VerdictVector::Dot(xet, xet) + VerdictVector::Dot(xze,xze));
  else
    return 0;
  
  if ( shape < min_shape ) { min_shape = shape; }
  

  // J(1,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  det = VerdictVector::Dot(xxi, (xet * xze));
  if ( det > VERDICT_DBL_MIN ) 
    shape = 3 * pow(det, two_thirds) / (VerdictVector::Dot(xxi, xxi) + VerdictVector::Dot(xet, xet) + VerdictVector::Dot(xze,xze));
  else
    return 0;
  
  if ( shape < min_shape ) { min_shape = shape; }

  if ( min_shape <= VERDICT_DBL_MIN )
    min_shape = 0;

  if ( min_shape > 0 )
    return (double) std::min( min_shape, VERDICT_DBL_MAX );
  return (double) std::max( min_shape, -VERDICT_DBL_MAX );
}

/*!
  relative size of a hex

  Min( J, 1/J ), where J is determinant of weighted Jacobian matrix
*/
double hex_relative_size_squared( int /*num_nodes*/, double coordinates[][3], double average_hex_volume )
{
  double size = 0;
  double tau; 

  VerdictVector xxi, xet, xze;
  double det, det_sum = 0;

  hex_get_weight( xxi, xet, xze, average_hex_volume );
 
  //This is the average relative size 
  double detw = VerdictVector::Dot(xxi, (xet * xze));

  if ( detw < VERDICT_DBL_MIN ) 
    return 0; 

  VerdictVector node_pos[8];
  make_hex_nodes ( coordinates, node_pos );

  // J(0,0,0):

  xxi = node_pos[1] - node_pos[0];
  xet = node_pos[3] - node_pos[0];
  xze = node_pos[4] - node_pos[0];
  
  det = VerdictVector::Dot(xxi, (xet * xze));
  det_sum += det;  


  // J(1,0,0):
  
  xxi = node_pos[2] - node_pos[1];
  xet = node_pos[0] - node_pos[1];
  xze = node_pos[5] - node_pos[1];

  det = VerdictVector::Dot(xxi, (xet * xze));
  det_sum += det;  


  // J(0,1,0):
  
  xxi = node_pos[3] - node_pos[2];
  xet = node_pos[1] - node_pos[2];
  xze = node_pos[6] - node_pos[2];

  det = VerdictVector::Dot(xxi, (xet * xze));
  det_sum += det;  

  // J(1,1,0):
  
  xxi = node_pos[0] - node_pos[3];
  xet = node_pos[2] - node_pos[3];
  xze = node_pos[7] - node_pos[3];

  det = VerdictVector::Dot(xxi, (xet * xze));
  det_sum += det;  


  // J(0,1,0):

  xxi = node_pos[7] - node_pos[4];
  xet = node_pos[5] - node_pos[4];
  xze = node_pos[0] - node_pos[4];

  det = VerdictVector::Dot(xxi, (xet * xze));
  det_sum += det;  


  // J(1,0,1):

  xxi = node_pos[4] - node_pos[5];
  xet = node_pos[6] - node_pos[5];
  xze = node_pos[1] - node_pos[5];

  det = VerdictVector::Dot(xxi, (xet * xze));
  det_sum += det;  


  // J(1,1,1):

  xxi = node_pos[5] - node_pos[6];
  xet = node_pos[7] - node_pos[6];
  xze = node_pos[2] - node_pos[6];

  det = VerdictVector::Dot(xxi, (xet * xze));
  det_sum += det;  


  // J(1,1,1):

  xxi = node_pos[6] - node_pos[7];
  xet = node_pos[4] - node_pos[7];
  xze = node_pos[3] - node_pos[7];

  det = VerdictVector::Dot(xxi, (xet * xze));
  det_sum += det;  


  if ( det_sum > VERDICT_DBL_MIN )
  {
    tau = det_sum / ( 8*detw);

    tau = std::min( tau, 1.0/tau); 

    size = tau*tau; 
  }

  if ( size > 0 )
    return (double) std::min( size, VERDICT_DBL_MAX );
  return (double) std::max( size, -VERDICT_DBL_MAX );
}

/*!
  shape and size of a hex

  Product of Shape and Relative Size
*/
double hex_shape_and_size( int num_nodes, double coordinates[][3], double average_hex_volume )
{
  double size = hex_relative_size_squared( num_nodes, coordinates, average_hex_volume );
  double shape = hex_shape( num_nodes, coordinates );

  double shape_size = size * shape;

  if ( shape_size > 0 )
    return (double) std::min( shape_size, VERDICT_DBL_MAX );
  return (double) std::max( shape_size, -VERDICT_DBL_MAX );

}

/*!
  shear and size of a hex

  Product of Shear and Relative Size
*/
double hex_shear_and_size( int num_nodes, double coordinates[][3], double average_hex_volume )
{
  double size = hex_relative_size_squared( num_nodes, coordinates, average_hex_volume );
  double shear = hex_shear( num_nodes, coordinates );

  double shear_size = shear * size; 

  if ( shear_size > 0 )
    return (double) std::min( shear_size, VERDICT_DBL_MAX );
  return (double) std::max( shear_size, -VERDICT_DBL_MAX );

}

/*!
  distortion of a hex
*/
double hex_distortion( int num_nodes, double coordinates[][3] )
{
  
  //use 2x2 gauss points for linear hex and 3x3 for 2nd order hex
  int number_of_gauss_points=0;
  if (num_nodes < 20)
  {
    //2x2 quadrature rule
    number_of_gauss_points = 2;
    num_nodes = 8;
  }
  else if (num_nodes >= 20)
  {
    //3x3 quadrature rule
    number_of_gauss_points = 3;
    num_nodes = 20;
  }
  
  int number_dimension = 3;
  int total_number_of_gauss_points = number_of_gauss_points
  *number_of_gauss_points*number_of_gauss_points;
  double distortion = VERDICT_DBL_MAX;
  
  // maxTotalNumberGaussPoints =27, maxNumberNodes = 20
  // they are defined in GaussIntegration.hpp
  // This is used to make these arrays static.
  // I tried dynamically allocated arrays but the new and delete
  // was very expensive
  
  double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
  double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
  double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
  double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
  double weight[maxTotalNumberGaussPoints];
  
  
  //create an object of GaussIntegration
  GaussIntegration gint;
  gint.initialize(number_of_gauss_points,num_nodes,number_dimension );
  gint.calculate_shape_function_3d_hex();
  gint.get_shape_func(shape_function[0], dndy1[0], dndy2[0], dndy3[0],weight);
  
  
  VerdictVector xxi, xet, xze, xin;
  
  double jacobian, minimum_jacobian;
  double element_volume =0.0;
  minimum_jacobian = VERDICT_DBL_MAX;
  // calculate element volume
  int ife, ja;
  for (ife=0;ife<total_number_of_gauss_points; ife++)
  {
    
    xxi.set(0.0,0.0,0.0);
    xet.set(0.0,0.0,0.0);
    xze.set(0.0,0.0,0.0);
    
    for (ja=0;ja<num_nodes;ja++)
    {
      xin.set(coordinates[ja][0], coordinates[ja][1], coordinates[ja][2]);
      xxi += dndy1[ife][ja]*xin;
      xet += dndy2[ife][ja]*xin;
      xze += dndy3[ife][ja]*xin;
    }
    
    jacobian = VerdictVector::Dot(xxi, (xet * xze));
    if (minimum_jacobian > jacobian)
      minimum_jacobian = jacobian;
    
    element_volume += weight[ife]*jacobian;
  }
  
  // loop through all nodes
  double dndy1_at_node[maxNumberNodes][maxNumberNodes];
  double dndy2_at_node[maxNumberNodes][maxNumberNodes];
  double dndy3_at_node[maxNumberNodes][maxNumberNodes];
  
  gint.calculate_derivative_at_nodes_3d( dndy1_at_node, dndy2_at_node, dndy3_at_node);
  int node_id;
  for (node_id=0;node_id<num_nodes; node_id++)
  {
    
    xxi.set(0.0,0.0,0.0);
    xet.set(0.0,0.0,0.0);
    xze.set(0.0,0.0,0.0);
    
    for (ja=0;ja<num_nodes;ja++)
    {
      xin.set(coordinates[ja][0], coordinates[ja][1], coordinates[ja][2]);
      xxi += dndy1_at_node[node_id][ja]*xin;
      xet += dndy2_at_node[node_id][ja]*xin;
      xze += dndy3_at_node[node_id][ja]*xin;
    }
    
    jacobian = VerdictVector::Dot(xxi, (xet * xze));
    if (minimum_jacobian > jacobian)
      minimum_jacobian = jacobian;
    
  }
  if (fabs(element_volume) > 0.0)
    distortion = minimum_jacobian/element_volume*8.;
  if      (distortion> VERDICT_DBL_MAX) distortion = VERDICT_DBL_MAX;
  else if (distortion<-VERDICT_DBL_MAX) distortion = -VERDICT_DBL_MAX;
  else if (std::isnan(distortion))           distortion = VERDICT_DBL_MAX; // 0/0, or should we return some other value?
  
  return (double)distortion;
}
  
double hex_timestep( int num_nodes, double coordinates[][3], 
                     double density,
                     double poissons_ratio,
                     double youngs_modulus )
{
  double char_length = hex_dimension( num_nodes, coordinates );
  double M = youngs_modulus*(1 - poissons_ratio) / ((1 - 2 * poissons_ratio)*(1 + poissons_ratio));
  double denominator = sqrt(M / density);

  return char_length / denominator;
}



/*
double hex_jac_normjac_oddy_cond( int choices[],
                      double coordinates[][3],
                      double answers[4]  )
{

  //Define variables
  int i;
 
  double xxi[3], xet[3], xze[3];
  double norm_jacobian = 0.0, current_norm_jac = 0.0; 
        double jacobian = 0.0, current_jacobian = 0.0;
  double oddy = 0.0, current_oddy = 0.0;  
  double condition = 0.0, current_condition = 0.0;


        for( i=0; i<3; i++)
          xxi[i] = calc_hex_efg( 2, i, coordinates );
        for( i=0; i<3; i++)
          xet[i] = calc_hex_efg( 3, i, coordinates );
        for( i=0; i<3; i++)
          xze[i] = calc_hex_efg( 6, i, coordinates );

  // norm jacobian and jacobian
  if ( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze  );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 

    current_jacobian /= 64.0;
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if ( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if ( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }


  for( i=0; i<3; i++)
  {
    xxi[i] = coordinates[1][i] - coordinates[0][i];
    xet[i] = coordinates[3][i] - coordinates[0][i];
    xze[i] = coordinates[4][i] - coordinates[0][i];
  } 

  // norm jacobian and jacobian
  if ( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze  );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if ( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if ( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

  
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[1][i] - coordinates[0][i];
          xet[i] = coordinates[2][i] - coordinates[1][i];
          xze[i] = coordinates[5][i] - coordinates[1][i];
  }

  // norm jacobian and jacobian
  if ( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi,  xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if ( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if ( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }


  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[2][i] - coordinates[3][i];
          xet[i] = coordinates[3][i] - coordinates[0][i];
          xze[i] = coordinates[7][i] - coordinates[3][i];
  }

  // norm jacobian and jacobian
  if ( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if ( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if ( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

  
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[5][i] - coordinates[4][i];
          xet[i] = coordinates[7][i] - coordinates[4][i];
          xze[i] = coordinates[4][i] - coordinates[0][i];
  }  


  // norm jacobian and jacobian
  if ( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if ( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if ( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

  
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[2][i] - coordinates[3][i];
          xet[i] = coordinates[2][i] - coordinates[1][i];
          xze[i] = coordinates[6][i] - coordinates[2][i];
  }  

  // norm jacobian and jacobian
  if ( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if ( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if ( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }


  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[5][i] - coordinates[4][i];
          xet[i] = coordinates[6][i] - coordinates[5][i];
          xze[i] = coordinates[5][i] - coordinates[1][i];
  }


  // norm jacobian and jacobian
  if ( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if ( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if ( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

    
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[6][i] - coordinates[7][i];
          xet[i] = coordinates[7][i] - coordinates[4][i];
          xze[i] = coordinates[7][i] - coordinates[3][i];
  }


  // norm jacobian and jacobian
  if ( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if ( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if ( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }
  }

  
  for( i=0; i<3; i++)
  {
          xxi[i] = coordinates[6][i] - coordinates[7][i];
          xet[i] = coordinates[6][i] - coordinates[5][i];
          xze[i] = coordinates[6][i] - coordinates[2][i];
  }

  // norm jacobian and jacobian
  if ( choices[0] || choices[1] )
  {
    current_jacobian = determinant( xxi, xet, xze );
    current_norm_jac = normalize_jacobian( current_jacobian,
            xxi, xet, xze );
    
    if (current_norm_jac < norm_jacobian) { norm_jacobian = current_norm_jac; } 
    if (current_jacobian < jacobian) { jacobian = current_jacobian; }
  }

  // oddy 
  if ( choices[2] )
  {
    current_oddy = oddy_comp( xxi, xet, xze );
      if ( current_oddy > oddy ) { oddy = current_oddy; }
  }

  // condition
  if ( choices[3] )
  {
    current_condition = condition_comp( xxi, xet, xze );
    if ( current_condition > condition ) { condition = current_condition; }

    condition /= 3.0;
  }

  
  answers[0] = jacobian;
  answers[1] = norm_jacobian;
  answers[2] = oddy;
  answers[3] = condition;

  return 1.0;

}
*/

} // namespace verdict