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

  Module:    V_TetMetric.cpp

  Copyright (c) 2006 Sandia Corporation.
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

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


/*
 *
 * TetMetric.cpp contains quality calculations for Tets
 *
 * This file is part of VERDICT
 *
 */


#define VERDICT_EXPORTS

#include "verdict.h"
#include "verdict_defines.hpp"
#include "VerdictVector.hpp"
#include "V_GaussIntegration.hpp"
#include <memory.h>

//! the average volume of a tet
double verdict_tet_size = 0;

/*! 
  set the average volume of a tet
*/
C_FUNC_DEF void v_set_tet_size( double size )
{
  verdict_tet_size = size;
}

/*!
  get the weights based on the average size
  of a tet
*/
int get_weight ( VerdictVector &w1,
                 VerdictVector &w2,
                 VerdictVector &w3 )
{
  static const double rt3 = sqrt(3.0);
  static const double root_of_2 = sqrt(2.0);
  
  w1.set(1,0,0);
  w2.set(0.5, 0.5*rt3, 0 );
  w3.set(0.5, rt3/6.0, root_of_2/rt3); 

  double scale = pow( 6.*verdict_tet_size/determinant(w1,w2,w3),0.3333333333333);   

  w1 *= scale;
  w2 *= scale;
  w3 *= scale;
  
  return 1;
}

/*!
   the edge ratio of a tet

   NB (P. Pebay 01/22/07): 
     Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
     minimum edge lengths
*/
C_FUNC_DEF double v_tet_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
{
  VerdictVector a, b, c, d, e, f;

  a.set( coordinates[1][0] - coordinates[0][0],
         coordinates[1][1] - coordinates[0][1],
         coordinates[1][2] - coordinates[0][2] );
  
  b.set( coordinates[2][0] - coordinates[1][0],
         coordinates[2][1] - coordinates[1][1],
         coordinates[2][2] - coordinates[1][2] );
  
  c.set( coordinates[0][0] - coordinates[2][0],
         coordinates[0][1] - coordinates[2][1],
         coordinates[0][2] - coordinates[2][2] );
  
  d.set( coordinates[3][0] - coordinates[0][0],
         coordinates[3][1] - coordinates[0][1],
         coordinates[3][2] - coordinates[0][2] );
  
  e.set( coordinates[3][0] - coordinates[1][0],
         coordinates[3][1] - coordinates[1][1],
         coordinates[3][2] - coordinates[1][2] );
  
  f.set( coordinates[3][0] - coordinates[2][0],
         coordinates[3][1] - coordinates[2][1],
         coordinates[3][2] - coordinates[2][2] );
  
  double a2 = a.length_squared();
  double b2 = b.length_squared();
  double c2 = c.length_squared();
  double d2 = d.length_squared();
  double e2 = e.length_squared();
  double f2 = f.length_squared();

  double m2,M2,mab,mcd,mef,Mab,Mcd,Mef;

  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;
    }

  m2 = mab < mcd ? mab : mcd;
  m2 = m2  < mef ? m2  : mef;

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

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

  double edge_ratio = sqrt( M2 / m2 );
  
  if( edge_ratio > 0 )
    return (double) VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
  return (double) VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
}

/*!
  the scaled jacobian of a tet

  minimum of the jacobian divided by the lengths of 3 edge vectors

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

  VerdictVector side0, side1, side2, side3, side4, side5;

  side0.set( coordinates[1][0] - coordinates[0][0],
             coordinates[1][1] - coordinates[0][1],
             coordinates[1][2] - coordinates[0][2] );
  
  side1.set( coordinates[2][0] - coordinates[1][0],
             coordinates[2][1] - coordinates[1][1],
             coordinates[2][2] - coordinates[1][2] );
  
  side2.set( coordinates[0][0] - coordinates[2][0],
             coordinates[0][1] - coordinates[2][1],
             coordinates[0][2] - coordinates[2][2] );

  side3.set( coordinates[3][0] - coordinates[0][0],
             coordinates[3][1] - coordinates[0][1],
             coordinates[3][2] - coordinates[0][2] );
  
  side4.set( coordinates[3][0] - coordinates[1][0],
             coordinates[3][1] - coordinates[1][1],
             coordinates[3][2] - coordinates[1][2] );
  
  side5.set( coordinates[3][0] - coordinates[2][0],
             coordinates[3][1] - coordinates[2][1],
             coordinates[3][2] - coordinates[2][2] );
  
  double jacobi;

  jacobi = side3 % ( side2 * side0 );

  // products of lengths squared of each edge attached to a node.
  double length_squared[4] = {
    side0.length_squared() * side2.length_squared() * side3.length_squared(),
    side0.length_squared() * side1.length_squared() * side4.length_squared(),
    side1.length_squared() * side2.length_squared() * side5.length_squared(),
    side3.length_squared() * side4.length_squared() * side5.length_squared()
  };
  int which_node = 0;
  if(length_squared[1] > length_squared[which_node])
    which_node = 1;
  if(length_squared[2] > length_squared[which_node])
    which_node = 2;
  if(length_squared[3] > length_squared[which_node])
    which_node = 3;
  
  double length_product = sqrt( length_squared[which_node] );
  if(length_product < fabs(jacobi))
    length_product = fabs(jacobi);

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

  static const double root_of_2 = sqrt(2.0);

  return (double)(root_of_2 * jacobi / length_product);

}

/*!
  The radius ratio of a tet

  NB (P. Pebay 04/16/07):
    CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed 
    sphere radius.
    Note that this function is similar to the aspect beta of a tet, except that
    it does not return VERDICT_DBL_MAX if the element has negative orientation.
*/
C_FUNC_DEF double v_tet_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
{

  //Determine side vectors
  VerdictVector side[6];

  side[0].set( coordinates[1][0] - coordinates[0][0],
               coordinates[1][1] - coordinates[0][1],
               coordinates[1][2] - coordinates[0][2] );
  
  side[1].set( coordinates[2][0] - coordinates[1][0],
               coordinates[2][1] - coordinates[1][1],
               coordinates[2][2] - coordinates[1][2] );
  
  side[2].set( coordinates[0][0] - coordinates[2][0],
               coordinates[0][1] - coordinates[2][1],
               coordinates[0][2] - coordinates[2][2] );

  side[3].set( coordinates[3][0] - coordinates[0][0],
               coordinates[3][1] - coordinates[0][1],
               coordinates[3][2] - coordinates[0][2] );
  
  side[4].set( coordinates[3][0] - coordinates[1][0],
               coordinates[3][1] - coordinates[1][1],
               coordinates[3][2] - coordinates[1][2] );
  
  side[5].set( coordinates[3][0] - coordinates[2][0],
               coordinates[3][1] - coordinates[2][1],
               coordinates[3][2] - coordinates[2][2] );

  VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0]) +
                            side[2].length_squared() * ( side[3] * side[0]) +
                            side[0].length_squared() * ( side[3] * side[2]);

  double area_sum;
  area_sum = ((side[2] * side[0]).length() + 
              (side[3] * side[0]).length() +
              (side[4] * side[1]).length() + 
              (side[3] * side[2]).length() ) * 0.5;
  
  double volume = v_tet_volume(4, coordinates);
  
  if( fabs( volume ) < VERDICT_DBL_MIN ) 
    return (double)VERDICT_DBL_MAX;
  else
  {
    double radius_ratio;
    radius_ratio = numerator.length() * area_sum / (108 * volume * volume); 
    
    return (double) VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
  }

}

/*!
  The radius ratio of a positively-oriented tet, a.k.a. "aspect beta"

  NB (P. Pebay 04/16/07):  
    CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed 
    sphere radius if the element has positive orientation.
    Note that this function is similar to the radius ratio of a tet, except that
    it returns VERDICT_DBL_MAX if the element has negative orientation.

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

  //Determine side vectors
  VerdictVector side[6];

  side[0].set( coordinates[1][0] - coordinates[0][0],
               coordinates[1][1] - coordinates[0][1],
               coordinates[1][2] - coordinates[0][2] );
  
  side[1].set( coordinates[2][0] - coordinates[1][0],
               coordinates[2][1] - coordinates[1][1],
               coordinates[2][2] - coordinates[1][2] );
  
  side[2].set( coordinates[0][0] - coordinates[2][0],
               coordinates[0][1] - coordinates[2][1],
               coordinates[0][2] - coordinates[2][2] );

  side[3].set( coordinates[3][0] - coordinates[0][0],
               coordinates[3][1] - coordinates[0][1],
               coordinates[3][2] - coordinates[0][2] );
  
  side[4].set( coordinates[3][0] - coordinates[1][0],
               coordinates[3][1] - coordinates[1][1],
               coordinates[3][2] - coordinates[1][2] );
  
  side[5].set( coordinates[3][0] - coordinates[2][0],
               coordinates[3][1] - coordinates[2][1],
               coordinates[3][2] - coordinates[2][2] );

  VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0]) +
                            side[2].length_squared() * ( side[3] * side[0]) +
                            side[0].length_squared() * ( side[3] * side[2]);

  double area_sum;
  area_sum = ((side[2] * side[0]).length() + 
              (side[3] * side[0]).length() +
              (side[4] * side[1]).length() + 
              (side[3] * side[2]).length() ) * 0.5;
  
  double volume = v_tet_volume(4, coordinates);
  
  if( volume < VERDICT_DBL_MIN ) 
    return (double)VERDICT_DBL_MAX;
  else
  {
    double radius_ratio;
    radius_ratio = numerator.length() * area_sum / (108 * volume * volume); 
    
    return (double) VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
  }

}

/*!
  The aspect ratio of a tet

  NB (P. Pebay 01/22/07):
    Hmax / (2 sqrt(6) r) where Hmax and r respectively denote the greatest edge 
    length and the inradius of the tetrahedron
*/
C_FUNC_DEF double v_tet_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
{
  static const double normal_coeff = sqrt(6.) / 12.;

  //Determine side vectors
  VerdictVector ab, bc, ac, ad, bd, cd;

  ab.set( coordinates[1][0] - coordinates[0][0],
          coordinates[1][1] - coordinates[0][1],
          coordinates[1][2] - coordinates[0][2] );
  
  ac.set( coordinates[2][0] - coordinates[0][0],
          coordinates[2][1] - coordinates[0][1],
          coordinates[2][2] - coordinates[0][2] );

  ad.set( coordinates[3][0] - coordinates[0][0],
          coordinates[3][1] - coordinates[0][1],
          coordinates[3][2] - coordinates[0][2] );
  
  double detTet = ab % ( ac * ad );
  
  if( detTet < VERDICT_DBL_MIN ) 
    return (double)VERDICT_DBL_MAX;

  bc.set( coordinates[2][0] - coordinates[1][0],
          coordinates[2][1] - coordinates[1][1],
          coordinates[2][2] - coordinates[1][2] );
  
  bd.set( coordinates[3][0] - coordinates[1][0],
          coordinates[3][1] - coordinates[1][1],
          coordinates[3][2] - coordinates[1][2] );
  
  cd.set( coordinates[3][0] - coordinates[2][0],
          coordinates[3][1] - coordinates[2][1],
          coordinates[3][2] - coordinates[2][2] );

  double ab2 = ab.length_squared();
  double bc2 = bc.length_squared();
  double ac2 = ac.length_squared();
  double ad2 = ad.length_squared();
  double bd2 = bd.length_squared();
  double cd2 = cd.length_squared();

  double A = ab2 > bc2 ? ab2 : bc2;
  double B = ac2 > ad2 ? ac2 : ad2;
  double C = bd2 > cd2 ? bd2 : cd2;
  double D = A > B ? A : B;
  double hm = D > C ? sqrt( D ) : sqrt( C );

  bd = ab * bc;
  A = bd.length();
  bd = ab * ad;
  B = bd.length();
  bd = ac * ad;
  C = bd.length();
  bd = bc * cd;
  D = bd.length();
  
  double aspect_ratio;
  aspect_ratio = normal_coeff * hm * ( A + B + C + D ) / fabs( detTet );
  
  if( aspect_ratio > 0 )
    return (double) VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
  return (double) VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
}

/*!
  the aspect gamma of a tet

  srms^3 / (8.48528137423857*V) where srms = sqrt(sum(Si^2)/6), where Si is the edge length
*/
C_FUNC_DEF double v_tet_aspect_gamma( int /*num_nodes*/, double coordinates[][3] )
{

  //Determine side vectors
  VerdictVector side0, side1, side2, side3, side4, side5;

  side0.set( coordinates[1][0] - coordinates[0][0],
             coordinates[1][1] - coordinates[0][1], 
             coordinates[1][2] - coordinates[0][2] ); 

  side1.set( coordinates[2][0] - coordinates[1][0],
             coordinates[2][1] - coordinates[1][1],
             coordinates[2][2] - coordinates[1][2] );
  
  side2.set( coordinates[0][0] - coordinates[2][0],
             coordinates[0][1] - coordinates[2][1],
             coordinates[0][2] - coordinates[2][2] );

  side3.set( coordinates[3][0] - coordinates[0][0],
             coordinates[3][1] - coordinates[0][1],
             coordinates[3][2] - coordinates[0][2] );
  
  side4.set( coordinates[3][0] - coordinates[1][0],
             coordinates[3][1] - coordinates[1][1],
             coordinates[3][2] - coordinates[1][2] );
  
  side5.set( coordinates[3][0] - coordinates[2][0],
             coordinates[3][1] - coordinates[2][1],
             coordinates[3][2] - coordinates[2][2] );
  

  double volume = fabs( v_tet_volume(4, coordinates) );

  if( volume  < VERDICT_DBL_MIN )
    return (double)VERDICT_DBL_MAX;
  else
  {
    double srms = sqrt((side0.length_squared() + side1.length_squared() +
                        side2.length_squared() + side3.length_squared() +
                        side4.length_squared() + side5.length_squared()) / 6.0 );

    double aspect_ratio_gamma = pow(srms, 3) / (8.48528137423857 * volume );  
    return (double)aspect_ratio_gamma;
  }
}

/*!
  The aspect frobenius of a tet

  NB (P. Pebay 01/22/07):
    Frobenius condition number when the reference element is regular
*/
C_FUNC_DEF double v_tet_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
{
  static const double normal_exp = 1. / 3.;

  VerdictVector ab, ac, ad;

  ab.set( coordinates[1][0] - coordinates[0][0],
          coordinates[1][1] - coordinates[0][1],
          coordinates[1][2] - coordinates[0][2] );
  
  ac.set( coordinates[2][0] - coordinates[0][0],
          coordinates[2][1] - coordinates[0][1],
          coordinates[2][2] - coordinates[0][2] );

  ad.set( coordinates[3][0] - coordinates[0][0],
          coordinates[3][1] - coordinates[0][1],
          coordinates[3][2] - coordinates[0][2] );
  
  double denominator = ab % ( ac * ad );
  denominator *= denominator;
  denominator *= 2.;
  denominator = 3. * pow( denominator, normal_exp );

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

  double u[3];
  ab.get_xyz( u );
  double v[3];
  ac.get_xyz( v );
  double w[3];
  ad.get_xyz( w );

  double numerator = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
  numerator += v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
  numerator += w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
  numerator *= 1.5;
  numerator -= v[0] * u[0] + v[1] * u[1] + v[2] * u[2];
  numerator -= w[0] * u[0] + w[1] * u[1] + w[2] * u[2];
  numerator -= w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
  
  double aspect_frobenius = numerator / denominator;
  
  if( aspect_frobenius > 0 )
    return (double) VERDICT_MIN( aspect_frobenius, VERDICT_DBL_MAX );
  return (double) VERDICT_MAX( aspect_frobenius, -VERDICT_DBL_MAX );
}

/*!
  The minimum angle of a tet

  NB (P. Pebay 01/22/07):
    minimum nonoriented dihedral angle
*/
C_FUNC_DEF double v_tet_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
{
  static const double normal_coeff = 180. * .3183098861837906715377675267450287;

  //Determine side vectors
  VerdictVector ab, bc, ad, cd;

  ab.set( coordinates[1][0] - coordinates[0][0],
          coordinates[1][1] - coordinates[0][1],
          coordinates[1][2] - coordinates[0][2] );
  
  ad.set( coordinates[3][0] - coordinates[0][0],
          coordinates[3][1] - coordinates[0][1],
          coordinates[3][2] - coordinates[0][2] );
  
  bc.set( coordinates[2][0] - coordinates[1][0],
          coordinates[2][1] - coordinates[1][1],
          coordinates[2][2] - coordinates[1][2] );
  
  cd.set( coordinates[3][0] - coordinates[2][0],
          coordinates[3][1] - coordinates[2][1],
          coordinates[3][2] - coordinates[2][2] );

  VerdictVector abc = ab * bc;
  double nabc = abc.length();
  VerdictVector abd = ab * ad;
  double nabd = abd.length();
  VerdictVector acd = ad * cd;
  double nacd = acd.length();
  VerdictVector bcd = bc * cd;
  double nbcd = bcd.length();

  double alpha   = acos( ( abc % abd ) / ( nabc * nabd ) );
  double beta    = acos( ( abc % acd ) / ( nabc * nacd ) );
  double gamma   = acos( ( abc % bcd ) / ( nabc * nbcd ) );
  double delta   = acos( ( abd % acd ) / ( nabd * nacd ) );
  double epsilon = acos( ( abd % bcd ) / ( nabd * nbcd ) );
  double zeta    = acos( ( acd % bcd ) / ( nacd * nbcd ) );

  alpha = alpha < beta    ? alpha : beta;
  alpha = alpha < gamma   ? alpha : gamma;
  alpha = alpha < delta   ? alpha : delta;
  alpha = alpha < epsilon ? alpha : epsilon;
  alpha = alpha < zeta    ? alpha : zeta;
  alpha *= normal_coeff;

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

  if( alpha > 0 )
    return (double) VERDICT_MIN( alpha, VERDICT_DBL_MAX );
  return (double) VERDICT_MAX( alpha, -VERDICT_DBL_MAX );
}

/*!
  The collapse ratio of a tet
*/
C_FUNC_DEF double v_tet_collapse_ratio( int /*num_nodes*/, double coordinates[][3] )
{
  //Determine side vectors
  VerdictVector e01, e02, e03, e12, e13, e23;

  e01.set( coordinates[1][0] - coordinates[0][0],
           coordinates[1][1] - coordinates[0][1],
           coordinates[1][2] - coordinates[0][2] );
  
  e02.set( coordinates[2][0] - coordinates[0][0],
           coordinates[2][1] - coordinates[0][1],
           coordinates[2][2] - coordinates[0][2] );
 
  e03.set( coordinates[3][0] - coordinates[0][0],
           coordinates[3][1] - coordinates[0][1],
           coordinates[3][2] - coordinates[0][2] );
  
  e12.set( coordinates[2][0] - coordinates[1][0],
           coordinates[2][1] - coordinates[1][1],
           coordinates[2][2] - coordinates[1][2] );
  
  e13.set( coordinates[3][0] - coordinates[1][0],
           coordinates[3][1] - coordinates[1][1],
           coordinates[3][2] - coordinates[1][2] );
  
  e23.set( coordinates[3][0] - coordinates[2][0],
           coordinates[3][1] - coordinates[2][1],
           coordinates[3][2] - coordinates[2][2] );

  double l[6];
  l[0] = e01.length();
  l[1] = e02.length();
  l[2] = e03.length();
  l[3] = e12.length();
  l[4] = e13.length();
  l[5] = e23.length();

  // Find longest edge for each bounding triangle of tetrahedron
  double l012 = l[4] > l[0] ? l[4] : l[0]; l012 = l[1] > l012 ? l[1] : l012;
  double l031 = l[0] > l[2] ? l[0] : l[2]; l031 = l[3] > l031 ? l[3] : l031;
  double l023 = l[2] > l[1] ? l[2] : l[1]; l023 = l[5] > l023 ? l[5] : l023;
  double l132 = l[4] > l[3] ? l[4] : l[3]; l132 = l[5] > l132 ? l[5] : l132;

  // Compute collapse ratio for each vertex/triangle pair
  VerdictVector N;
  double h, magN;
  double cr;
  double crMin;

  N = e01 * e02;
  magN = N.length();
  h = ( e03 % N) / magN; // height of vertex 3 above 0-1-2
  crMin = h / l012;      // ratio of height to longest edge of 0-1-2

  N = e03 * e01 ;
  magN = N.length();
  h = ( e02 % N) / magN; // height of vertex 2 above 0-3-1
  cr = h / l031;         // ratio of height to longest edge of 0-3-1
  if ( cr < crMin ) crMin = cr;

  N = e02 * e03 ;
  magN = N.length();
  h = ( e01 % N) / magN; // height of vertex 1 above 0-2-3
  cr = h / l023;         // ratio of height to longest edge of 0-2-3
  if ( cr < crMin ) crMin = cr;

  N = e12 * e13 ;
  magN = N.length();
  h = ( e01 % N ) / magN; // height of vertex 0 above 1-3-2
  cr = h / l132;          // ratio of height to longest edge of 1-3-2
  if ( cr < crMin ) crMin = cr;

  if( crMin < VERDICT_DBL_MIN ) 
    return (double)VERDICT_DBL_MAX;
  if( crMin > 0 )
    return (double) VERDICT_MIN( crMin, VERDICT_DBL_MAX );
  return (double) VERDICT_MAX( crMin, -VERDICT_DBL_MAX );
}

/*!
  the volume of a tet

  1/6 * jacobian at a corner node
*/
C_FUNC_DEF double v_tet_volume( int /*num_nodes*/, double coordinates[][3] )
{

  //Determine side vectors
  VerdictVector side0, side2, side3;

  side0.set( coordinates[1][0] - coordinates[0][0],
             coordinates[1][1] - coordinates[0][1],
             coordinates[1][2] - coordinates[0][2] );
  
  side2.set( coordinates[0][0] - coordinates[2][0],
             coordinates[0][1] - coordinates[2][1],
             coordinates[0][2] - coordinates[2][2] );
  
  side3.set( coordinates[3][0] - coordinates[0][0],
             coordinates[3][1] - coordinates[0][1],
             coordinates[3][2] - coordinates[0][2] );

  return  (double)((side3 % (side2 * side0)) / 6.0);

}

/*!
  the condition of a tet

  condition number of the jacobian matrix at any corner
*/
C_FUNC_DEF double v_tet_condition( int /*num_nodes*/, double coordinates[][3] )
{

  double condition, term1, term2, det;
  double rt3 = sqrt(3.0);
  double rt6 = sqrt(6.0);

  VerdictVector side0, side2, side3;
  
  side0.set(coordinates[1][0] - coordinates[0][0],
            coordinates[1][1] - coordinates[0][1],
            coordinates[1][2] - coordinates[0][2]);

  side2.set(coordinates[0][0] - coordinates[2][0],
            coordinates[0][1] - coordinates[2][1],
            coordinates[0][2] - coordinates[2][2]);

  side3.set(coordinates[3][0] - coordinates[0][0],
            coordinates[3][1] - coordinates[0][1],
            coordinates[3][2] - coordinates[0][2]);

  VerdictVector c_1, c_2, c_3; 

  c_1 = side0;
  c_2 = (-2*side2-side0)/rt3;
  c_3 = (3*side3+side2-side0)/rt6;

  term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
  term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) + 
          ( c_2 * c_3 ) % ( c_2 * c_3 ) + 
          ( c_1 * c_3 ) % ( c_1 * c_3 );
  det = c_1 % ( c_2 * c_3 );
  
  if ( det <= VERDICT_DBL_MIN )
    return VERDICT_DBL_MAX;
  else 
    condition = sqrt( term1 * term2 ) /(3.0* det);
  
  return (double)condition;
}


/*!
  the jacobian of a tet

  TODO
*/
C_FUNC_DEF double v_tet_jacobian( int /*num_nodes*/, double coordinates[][3] )
{
  VerdictVector side0, side2, side3;

  side0.set( coordinates[1][0] - coordinates[0][0],
             coordinates[1][1] - coordinates[0][1],
             coordinates[1][2] - coordinates[0][2] );
  
  side2.set( coordinates[0][0] - coordinates[2][0],
             coordinates[0][1] - coordinates[2][1],
             coordinates[0][2] - coordinates[2][2] );

  side3.set( coordinates[3][0] - coordinates[0][0],
             coordinates[3][1] - coordinates[0][1],
             coordinates[3][2] - coordinates[0][2] );
  

  return (double)(side3 % (side2 * side0));

}


/*!
  the shape of a tet

  3/ condition number of weighted jacobian matrix
*/
C_FUNC_DEF double v_tet_shape( int /*num_nodes*/, double coordinates[][3] )
{

   static const double two_thirds = 2.0/3.0;
   static const double root_of_2 = sqrt(2.0);

   VerdictVector edge0, edge2, edge3;

  edge0.set(coordinates[1][0] - coordinates[0][0],
            coordinates[1][1] - coordinates[0][1],
            coordinates[1][2] - coordinates[0][2]);

  edge2.set(coordinates[0][0] - coordinates[2][0],
            coordinates[0][1] - coordinates[2][1],
            coordinates[0][2] - coordinates[2][2]);

  edge3.set(coordinates[3][0] - coordinates[0][0],
            coordinates[3][1] - coordinates[0][1],
            coordinates[3][2] - coordinates[0][2]);

  double jacobian = edge3 % (edge2 * edge0);
  if(jacobian < VERDICT_DBL_MIN){
    return (double)0.0;
  }
  double num = 3 * pow( root_of_2 * jacobian, two_thirds );
  double den = 1.5*(edge0%edge0  + edge2%edge2  + edge3%edge3)-
                   (edge0%-edge2 + -edge2%edge3 + edge3%edge0);

  if ( den < VERDICT_DBL_MIN ) 
    return (double)0.0;
    
  return (double)VERDICT_MAX( num/den, 0 );
}



/*!
  the relative size of a tet

  Min(J,1/J), where J is the determinant of the weighted Jacobian matrix
*/
C_FUNC_DEF double v_tet_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
{
  double size;
  VerdictVector w1, w2, w3;
  get_weight(w1,w2,w3);
  double avg_volume = (w1 % (w2 *w3))/6.0;
  
  double volume = v_tet_volume(4, coordinates);

  if( avg_volume < VERDICT_DBL_MIN )
    return 0.0;
  else
  {
    size = volume/avg_volume;
    if( size <= VERDICT_DBL_MIN )
      return 0.0;
    if ( size > 1 ) 
      size = (double)(1)/size;
  }
  return (double)(size*size);
}


/*!
  the shape and size of a tet

  Product of the shape and relative size
*/
C_FUNC_DEF double v_tet_shape_and_size( int num_nodes, double coordinates[][3] )
{
  
  double shape, size;
  shape = v_tet_shape( num_nodes, coordinates );
  size = v_tet_relative_size_squared (num_nodes, coordinates );  
  
  return (double)(shape * size);

}



/*!
  the distortion of a tet
*/
C_FUNC_DEF double v_tet_distortion( int num_nodes, double coordinates[][3] )
{

   double distortion = VERDICT_DBL_MAX;
   int   number_of_gauss_points=0;
   if (num_nodes ==4)
      // for linear tet, the distortion is always 1 because
      // straight edge tets are the target shape for tet
      return 1.0;

   else if (num_nodes ==10)
      //use four integration points for quadratic tet
      number_of_gauss_points = 4;

   int number_dims = 3;
   int total_number_of_gauss_points = number_of_gauss_points;
   // use is_tri=1 to indicate this is for tet in 3D
   int is_tri =1;

   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 for tet
   GaussIntegration::initialize(number_of_gauss_points,num_nodes, number_dims, is_tri);
   GaussIntegration::calculate_shape_function_3d_tet();
   GaussIntegration::get_shape_func(shape_function[0], dndy1[0], dndy2[0], dndy3[0],weight);

   // vector xxi is the derivative vector of coordinates w.r.t local xi coordinate in the
   // computation space
   // vector xet is the derivative vector of coordinates w.r.t local et coordinate in the
   // computation space
   // vector xze is the derivative vector of coordinates w.r.t local ze coordinate in the
   // computation space
   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;
      }

      // determinant
      jacobian = xxi % ( xet * xze );
      if (minimum_jacobian > jacobian)
         minimum_jacobian = jacobian;

      element_volume += weight[ife]*jacobian;
      }//element_volume is 6 times the actual volume

   // loop through all nodes
   double dndy1_at_node[maxNumberNodes][maxNumberNodes];
   double dndy2_at_node[maxNumberNodes][maxNumberNodes];
   double dndy3_at_node[maxNumberNodes][maxNumberNodes];


   GaussIntegration::calculate_derivative_at_nodes_3d_tet( 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 = xxi % ( xet * xze );
      if (minimum_jacobian > jacobian)
         minimum_jacobian = jacobian;

      }
   distortion = minimum_jacobian/element_volume;

   return (double)distortion;
}


/*!
  the quality functions of a tet
*/
C_FUNC_DEF void v_tet_quality( int num_nodes, double coordinates[][3], 
    unsigned int metrics_request_flag, TetMetricVals *metric_vals )
{

  memset( metric_vals, 0, sizeof(TetMetricVals) );

  /*
  
    node numbers and edge numbers below


    
             3 
             +            edge 0 is node 0 to 1
            +|+           edge 1 is node 1 to 2
          3/ | \5         edge 2 is node 0 to 2
          / 4|  \         edge 3 is node 0 to 3
        0 - -|- + 2       edge 4 is node 1 to 3
          \  |  +         edge 5 is node 2 to 3
          0\ | /1
            +|/           edge 2 is behind edge 4
             1 

             
  */

  // lets start with making the vectors
  VerdictVector edges[6];
  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[0][0] - coordinates[2][0],
                coordinates[0][1] - coordinates[2][1],
                coordinates[0][2] - coordinates[2][2] );

  edges[3].set( coordinates[3][0] - coordinates[0][0],
                coordinates[3][1] - coordinates[0][1],
                coordinates[3][2] - coordinates[0][2] );

  edges[4].set( coordinates[3][0] - coordinates[1][0],
                coordinates[3][1] - coordinates[1][1],
                coordinates[3][2] - coordinates[1][2] );

  edges[5].set( coordinates[3][0] - coordinates[2][0],
                coordinates[3][1] - coordinates[2][1],
                coordinates[3][2] - coordinates[2][2] );

  // common numbers
  static const double root_of_2 = sqrt(2.0);
 
  // calculate the jacobian 
  static const int do_jacobian = V_TET_JACOBIAN | V_TET_VOLUME | 
    V_TET_ASPECT_BETA | V_TET_ASPECT_GAMMA | V_TET_SHAPE | 
    V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE | 
    V_TET_SCALED_JACOBIAN | V_TET_CONDITION;
  if(metrics_request_flag & do_jacobian )
  {
    metric_vals->jacobian = (double)(edges[3] % (edges[2] * edges[0]));
  }
 
  // calculate the volume 
  if(metrics_request_flag & V_TET_VOLUME)
  {
    metric_vals->volume = (double)(metric_vals->jacobian / 6.0);
  }
  
  // calculate aspect ratio
  if(metrics_request_flag & V_TET_ASPECT_BETA)
  {
    double surface_area = ((edges[2] * edges[0]).length() + 
                           (edges[3] * edges[0]).length() +
                           (edges[4] * edges[1]).length() + 
                           (edges[3] * edges[2]).length() ) * 0.5;

    VerdictVector numerator = edges[3].length_squared() * ( edges[2] * edges[0] ) +
                              edges[2].length_squared() * ( edges[3] * edges[0] ) +
                              edges[0].length_squared() * ( edges[3] * edges[2] );

    double volume = metric_vals->jacobian / 6.0;

    if(volume < VERDICT_DBL_MIN )
      metric_vals->aspect_beta = (double)(VERDICT_DBL_MAX);
    else
      metric_vals->aspect_beta = 
        (double)( numerator.length() * surface_area/ (108*volume*volume) );
  }

  // calculate the aspect gamma 
  if(metrics_request_flag & V_TET_ASPECT_GAMMA)
  {
    double volume = fabs( metric_vals->jacobian / 6.0 );
    if( fabs( volume ) < VERDICT_DBL_MIN ) 
      metric_vals->aspect_gamma = VERDICT_DBL_MAX;
    else
    {
      double srms = sqrt((
            edges[0].length_squared() + edges[1].length_squared() +
            edges[2].length_squared() + edges[3].length_squared() +
            edges[4].length_squared() + edges[5].length_squared()
            ) / 6.0 );

      // cube the srms
      srms *= (srms * srms);
      metric_vals->aspect_gamma = (double)( srms / (8.48528137423857 * volume ));
    }
  }

  // calculate the shape of the tet
  if(metrics_request_flag & ( V_TET_SHAPE | V_TET_SHAPE_AND_SIZE ) )
  {
      //if the jacobian is non-positive, the shape is 0
    if(metric_vals->jacobian < VERDICT_DBL_MIN){
      metric_vals->shape = (double)0.0;
    }
    else{
      static const double two_thirds = 2.0/3.0;
      double num = 3.0 * pow(root_of_2 * metric_vals->jacobian, two_thirds);
      double den = 1.5 *
        (edges[0] % edges[0]  + edges[2] % edges[2]  + edges[3] % edges[3]) -
        (edges[0] % -edges[2] + -edges[2] % edges[3] + edges[3] % edges[0]);

      if( den < VERDICT_DBL_MIN )
        metric_vals->shape = (double)0.0;
      else
        metric_vals->shape = (double)VERDICT_MAX( num/den, 0 );
    }
    
  }
  
  // calculate the relative size of the tet
  if(metrics_request_flag & (V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE ))
  {
    VerdictVector w1, w2, w3;
    get_weight(w1,w2,w3);
    double avg_vol = (w1 % (w2 *w3))/6;
    
    if( avg_vol < VERDICT_DBL_MIN )
      metric_vals->relative_size_squared = 0.0; 
    else
    {
      double tmp = metric_vals->jacobian / (6*avg_vol);
      if( tmp < VERDICT_DBL_MIN )
        metric_vals->relative_size_squared = 0.0; 
      else
      {
        tmp *= tmp;
        metric_vals->relative_size_squared = (double)VERDICT_MIN(tmp, 1/tmp);
      }
    }
  }
  
  // calculate the shape and size
  if(metrics_request_flag & V_TET_SHAPE_AND_SIZE)
  {
    metric_vals->shape_and_size = (double)(metric_vals->shape * metric_vals->relative_size_squared);
  }
  
  // calculate the scaled jacobian
  if(metrics_request_flag & V_TET_SCALED_JACOBIAN)
  {
    //find out which node the normalized jacobian can be calculated at
    //and it will be the smaller than at other nodes
    double length_squared[4] = {
      edges[0].length_squared() * edges[2].length_squared() * edges[3].length_squared(),
      edges[0].length_squared() * edges[1].length_squared() * edges[4].length_squared(),
      edges[1].length_squared() * edges[2].length_squared() * edges[5].length_squared(),
      edges[3].length_squared() * edges[4].length_squared() * edges[5].length_squared()
    };
    
    int which_node = 0;
    if(length_squared[1] > length_squared[which_node])
      which_node = 1;
    if(length_squared[2] > length_squared[which_node])
      which_node = 2;
    if(length_squared[3] > length_squared[which_node])
      which_node = 3;

    // find the scaled jacobian at this node
    double length_product = sqrt( length_squared[which_node] );
    if(length_product < fabs(metric_vals->jacobian))
      length_product = fabs(metric_vals->jacobian);

    if( length_product < VERDICT_DBL_MIN )
      metric_vals->scaled_jacobian = (double) VERDICT_DBL_MAX; 
    else
      metric_vals->scaled_jacobian = 
        (double)(root_of_2 * metric_vals->jacobian / length_product);
  }
  
  // calculate the condition number
  if(metrics_request_flag & V_TET_CONDITION)
  {
    static const double root_of_3 = sqrt(3.0);
    static const double root_of_6 = sqrt(6.0);

    VerdictVector c_1, c_2, c_3; 
    c_1 = edges[0];
    c_2 = (-2*edges[2] - edges[0])/root_of_3;
    c_3 = (3*edges[3] + edges[2] - edges[0])/root_of_6;

    double term1 =  c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
    double term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) + 
                   ( c_2 * c_3 ) % ( c_2 * c_3 ) + 
                   ( c_3 * c_1 ) % ( c_3 * c_1 );

    double det = c_1 % ( c_2 * c_3 );

    if(det <= VERDICT_DBL_MIN)
      metric_vals->condition = (double)VERDICT_DBL_MAX; 
    else
      metric_vals->condition = (double)(sqrt(term1 * term2) / (3.0*det)); 
  }
    
  // calculate the distortion
  if(metrics_request_flag & V_TET_DISTORTION)
  {
    metric_vals->distortion = v_tet_distortion(num_nodes, coordinates);
  }

  //check for overflow
  if(metrics_request_flag & V_TET_ASPECT_BETA )
  {
    if( metric_vals->aspect_beta > 0 ) 
      metric_vals->aspect_beta = (double) VERDICT_MIN( metric_vals->aspect_beta, VERDICT_DBL_MAX );
    metric_vals->aspect_beta = (double) VERDICT_MAX( metric_vals->aspect_beta, -VERDICT_DBL_MAX );
  }

  if(metrics_request_flag & V_TET_ASPECT_GAMMA)
  {
    if( metric_vals->aspect_gamma > 0 ) 
      metric_vals->aspect_gamma = (double) VERDICT_MIN( metric_vals->aspect_gamma, VERDICT_DBL_MAX );
    metric_vals->aspect_gamma = (double) VERDICT_MAX( metric_vals->aspect_gamma, -VERDICT_DBL_MAX );
  }

  if(metrics_request_flag & V_TET_VOLUME)
  {
    if( metric_vals->volume > 0 ) 
      metric_vals->volume = (double) VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
    metric_vals->volume = (double) VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
  }

  if(metrics_request_flag & V_TET_CONDITION)
  {
    if( metric_vals->condition > 0 ) 
      metric_vals->condition = (double) VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
    metric_vals->condition = (double) VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
  }

  if(metrics_request_flag & V_TET_JACOBIAN)
  {
    if( metric_vals->jacobian > 0 ) 
      metric_vals->jacobian = (double) VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
    metric_vals->jacobian = (double) VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
  }

  if(metrics_request_flag & V_TET_SCALED_JACOBIAN)
  {
    if( metric_vals->scaled_jacobian > 0 ) 
      metric_vals->scaled_jacobian = (double) VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
    metric_vals->scaled_jacobian = (double) VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
  }

  if(metrics_request_flag & V_TET_SHAPE)
  {
    if( metric_vals->shape > 0 ) 
      metric_vals->shape = (double) VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
    metric_vals->shape = (double) VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
  }

  if(metrics_request_flag & V_TET_RELATIVE_SIZE_SQUARED)
  {
    if( metric_vals->relative_size_squared > 0 ) 
      metric_vals->relative_size_squared = (double) VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
    metric_vals->relative_size_squared = (double) VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
  }

  if(metrics_request_flag & V_TET_SHAPE_AND_SIZE)
  {
    if( metric_vals->shape_and_size > 0 ) 
      metric_vals->shape_and_size = (double) VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
    metric_vals->shape_and_size = (double) VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
  }

  if(metrics_request_flag & V_TET_DISTORTION)
  {
    if( metric_vals->distortion > 0 ) 
      metric_vals->distortion = (double) VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
    metric_vals->distortion = (double) VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
  }


}
