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

Module:    V_WedgeMetric.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.

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


/*
 *
 * WedgeMetric.cpp contains quality calculations for wedges
 *
 * This file is part of VERDICT
 *
 */


#include "verdict.h"
#include "VerdictVector.hpp"
#include <memory.h> 
#include "verdict_defines.hpp"
#include "V_GaussIntegration.hpp"
/*
   the wedge element


   5
   ^
   / \
   / | \
   / /2\ \
   6/_______\4
   | /   \ |
   |/_____\|
   3         1

 */



/*!

  calculate the volume of a wedge

  this is done by dividing the wedge into 3 tets
  and summing the volume of each tet

 */

C_FUNC_DEF double v_wedge_volume( int num_nodes, double coordinates[][3] )
{

  double volume = 0;
  VerdictVector side1, side2, side3;

  if ( num_nodes == 6 )
  {

    // divide the wedge into 3 tets and calculate each volume

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

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


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

    volume = side3 % (side1 * side2)  / 6;

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

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


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

    volume += side3 % (side1 * side2)  / 6;

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

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


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

    volume += side3 % (side1 * side2)  / 6;

  }

  return (double)volume;

}



C_FUNC_DEF void v_wedge_quality( int num_nodes, double coordinates[][3], 
    unsigned int metrics_request_flag, WedgeMetricVals *metric_vals )
{
  memset( metric_vals, 0, sizeof(WedgeMetricVals) );

  if(metrics_request_flag & V_WEDGE_VOLUME)
    metric_vals->volume = v_wedge_volume(num_nodes, coordinates);
  if(metrics_request_flag & V_WEDGE_EDGE_RATIO)
    metric_vals->edge_ratio = v_wedge_edge_ratio(num_nodes, coordinates);
  if(metrics_request_flag & V_WEDGE_MAX_ASPECT_FROBENIUS)
    metric_vals->max_aspect_frobenius = v_wedge_max_aspect_frobenius(num_nodes, coordinates);
  if(metrics_request_flag & V_WEDGE_MEAN_ASPECT_FROBENIUS)
    metric_vals->mean_aspect_frobenius = v_wedge_mean_aspect_frobenius(num_nodes, coordinates);
  if(metrics_request_flag & V_WEDGE_JACOBIAN)
    metric_vals->jacobian = v_wedge_jacobian(num_nodes, coordinates);
  if(metrics_request_flag & V_WEDGE_SCALED_JACOBIAN)
    metric_vals->scaled_jacobian = v_wedge_scaled_jacobian(num_nodes, coordinates);
  if(metrics_request_flag & V_WEDGE_DISTORTION)
    metric_vals->distortion = v_wedge_distortion(num_nodes, coordinates);
  if(metrics_request_flag & V_WEDGE_MAX_STRETCH)
    metric_vals->max_stretch = v_wedge_max_stretch(num_nodes, coordinates);
  if(metrics_request_flag & V_WEDGE_SHAPE)
    metric_vals->shape = v_wedge_shape(num_nodes, coordinates);
  if(metrics_request_flag & V_WEDGE_CONDITION)
    metric_vals->condition = v_wedge_condition(num_nodes, coordinates);
}


/* Edge ratio
   The edge ratio quality metric is the ratio of the longest to shortest edge of
   a wedge.
   q = L_max / L_min

   Dimension : 1
   Acceptable range : --
   Normal range : [1,DBL_MAX]
   Full range : [1,DBL_MAX]
   q for right, unit wedge : 1
   Reference : -
   */
C_FUNC_DEF double v_wedge_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
{
  VerdictVector a,b,c,d,e,f,g,h,i;

  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[4][0] - coordinates[3][0],
      coordinates[4][1] - coordinates[3][1],
      coordinates[4][2] - coordinates[3][2] );

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

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

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

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

  i.set( coordinates[5][0] - coordinates[2][0],
      coordinates[5][1] - coordinates[2][1],
      coordinates[5][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 g2 = g.length_squared();
  double h2 = h.length_squared();
  double i2 = i.length_squared();

  double max = a2, min = a2;

  if (max <= b2){max = b2;}
  if (b2 <= min){min = b2;}

  if (max <= c2){max = c2;}
  if (c2 <= min){min = c2;}

  if (max <= d2){max = d2;}
  if (d2 <= min){min = d2;}

  if (max <= e2){max = e2;}
  if (e2 <= min){min = e2;}

  if (max <= f2){max = f2;}
  if (f2 <= min){min = f2;}

  if (max <= g2){max = g2;}
  if (g2 <= min){min = g2;}

  if (max <= h2){max = h2;}
  if (h2 <= min){min = h2;}

  if (max <= i2){max = i2;}
  if (i2 <= min){min = i2;}

  double edge_ratio = sqrt( max / min );

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

}

/* For wedges, there is not a unique definition of the aspect Frobenius. Rather,
 * this metric uses the aspect Frobenius defined for tetrahedral (see section
 * 6.4) and is comparable in methodology to the maximum aspect Frobenius defined
 * for hexahedra (see section 7.7). This value is normalized for a unit wedge.

 q = max(F_0123, F_1204, F_2015, F_3540, F_4351, F_5432)

 This is also known as the wedge condition number.

 Dimension : 1
 Acceptable Range :
 Normal Range :
 Full Range :
 q for right, unit wedge : 1
 Reference : Adapted from section 7.7
 Verdict Function : v_wedge_max_aspect_frobenius or v_wedge_condition

 */

C_FUNC_DEF double v_wedge_max_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
{
  double mini_tris[4][3];
  double aspect1 = 0, aspect2 = 0, aspect3 = 0, aspect4 = 0, aspect6 = 0;
  int i = 0;
  // Take first tetrahedron
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[0][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[1][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[2][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[3][i];}

  aspect1 = v_tet_aspect_frobenius(4,mini_tris);

  //Take second tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[1][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[2][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[0][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[4][i];}

  aspect2 = v_tet_aspect_frobenius(4,mini_tris);

  //3rd tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[2][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[0][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[1][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[5][i];}

  aspect3 = v_tet_aspect_frobenius(4,mini_tris);

  //4th tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[3][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[5][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[4][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[0][i];}

  aspect4 = v_tet_aspect_frobenius(4,mini_tris);

  //5th tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[4][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[3][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[5][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[1][i];}

  v_tet_aspect_frobenius(4,mini_tris);

  //6th tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[5][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[4][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[3][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[2][i];}

  aspect6 = v_tet_aspect_frobenius(4,mini_tris);

  double max_aspect = VERDICT_MAX(aspect1,aspect2);
  max_aspect = VERDICT_MAX(max_aspect,aspect3);
  max_aspect = VERDICT_MAX(max_aspect,aspect4);
  max_aspect = VERDICT_MAX(max_aspect,aspect6);
  max_aspect = max_aspect/1.16477;

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

}
/*
   For wedges, there is not a unique definition of the aspect Frobenius. Rather,
   this metric uses the aspect Frobenius defined for tetrahedra (see section
   6.4) and is comparable in methodology to the mean aspect Frobenius defined
   for hexahedra (see section 7.8). This value is normalized for a unit wedge.

   q = 1/6 * (F_0123 + F_1204 + F+2015 + F_3540 + F_4351 + F_5432)

   Dimension : 1
   Acceptable Range :
   Normal Range :
   Full Range :
   q for right, unit wedge : 1
   Reference : Adapted from section 7.8
   Verdict Function : v_wedge_mean_aspect_frobenius

   */

C_FUNC_DEF double v_wedge_mean_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
{
  double mini_tris[4][3];
  double aspect1 = 0, aspect2 = 0, aspect3 = 0, aspect4 = 0, aspect5 = 0, aspect6 = 0;
  int i = 0;
  // Take first tetrahedron
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[0][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[1][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[2][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[3][i];}

  aspect1 = v_tet_aspect_frobenius(4,mini_tris);

  //Take second tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[1][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[2][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[0][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[4][i];}

  aspect2 = v_tet_aspect_frobenius(4,mini_tris);

  //3rd tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[2][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[0][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[1][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[5][i];}

  aspect3 = v_tet_aspect_frobenius(4,mini_tris);

  //4th tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[3][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[5][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[4][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[0][i];}

  aspect4 = v_tet_aspect_frobenius(4,mini_tris);

  //5th tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[4][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[3][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[5][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[1][i];}

  aspect5 = v_tet_aspect_frobenius(4,mini_tris);

  //6th tet
  for (i = 0; i < 3; i++){mini_tris[0][i] = coordinates[5][i];}
  for (i = 0; i < 3; i++){mini_tris[1][i] = coordinates[4][i];}
  for (i = 0; i < 3; i++){mini_tris[2][i] = coordinates[3][i];}
  for (i = 0; i < 3; i++){mini_tris[3][i] = coordinates[2][i];}

  aspect6 = v_tet_aspect_frobenius(4,mini_tris);

  double mean_aspect = (aspect1 + aspect2 + aspect3 + aspect4 + aspect5 + aspect6)/6;
  mean_aspect = mean_aspect/1.16477;

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

}

/* This is the minimum determinant of the Jacobian matrix evaluated at each
 * corner of the element.

 q = min[((L_2 X L_0) * L_3)_k]
 where ((L_2 X L_0) * L_3)_k is the determinant of the Jacobian of the
 tetrahedron defined at the kth corner node, and L_2, L_0 and L_3 are the edges
 defined according to the standard for tetrahedral elements.

 Dimension : L^3
 Acceptable Range : [0,DBL_MAX]
 Normal Range : [0,DBL_MAX]
 Full Range : [-DBL_MAX,DBL_MAX]
 q for right, unit wedge : sqrt(3)/2
 Reference : Adapted from section 6.10
 Verdict Function : v_wedge_jacobian
 */

C_FUNC_DEF double v_wedge_jacobian( int /*num_nodes*/, double coordinates[][3] )
{
  double min_jacobian = 0, current_jacobian = 0;
  VerdictVector vec1,vec2,vec3;

  // Node 0
  vec1.set( coordinates[1][0] - coordinates[0][0],
      coordinates[1][1] - coordinates[0][1],
      coordinates[1][2] - coordinates[0][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = current_jacobian;

  //node 1
  vec1.set( coordinates[2][0] - coordinates[1][0],
      coordinates[2][1] - coordinates[1][1],
      coordinates[2][2] - coordinates[1][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian,min_jacobian);

  //node 2
  vec1.set( coordinates[0][0] - coordinates[2][0],
      coordinates[0][1] - coordinates[2][1],
      coordinates[0][2] - coordinates[2][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian,min_jacobian);

  //node 3
  vec1.set( coordinates[0][0] - coordinates[3][0],
      coordinates[0][1] - coordinates[3][1],
      coordinates[0][2] - coordinates[3][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian,min_jacobian);

  //node 4
  vec1.set( coordinates[1][0] - coordinates[4][0],
      coordinates[1][1] - coordinates[4][1],
      coordinates[1][2] - coordinates[4][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian,min_jacobian);

  //node 5
  vec1.set( coordinates[3][0] - coordinates[5][0],
      coordinates[3][1] - coordinates[5][1],
      coordinates[3][2] - coordinates[5][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian,min_jacobian);

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

}

/* distortion is a measure of how well a particular wedge element maps to a
 * 'master' wedge with vertices:
 P0 - (0, 0, 0)
 P1 - (1, 0, 0)
 P2 - (1/2, sqrt(3)/2, 0)
 P3 - (0, 0, 1)
 P4 - (1, 0, 1)
 P2 - (1/2, sqrt(3)/2, 1)
 and volume (V_m).
 The behavior of the map is measured by sampling the determinant of the Jacobian
 at the vertices (k).  Thus the distortion is given by:
 q = ( min_k { det(J_k)} * V_m ) / V

 Dimension : 1
 Acceptable Range : [0.5,1]
 Normal Range : [0,1]
 Full Range : [-DBL_MAX,DBL_MAX]
 q for right, unit wedge : 1
 Reference : Adapted from section 7.3
 Verdict Function : v_wedge_distortion
 */

C_FUNC_DEF double v_wedge_distortion( int num_nodes, double coordinates[][3] )
{

  double jacobian = 0, distortion = 45, current_volume = 0, master_volume = 0.433013;

  jacobian = v_wedge_jacobian( num_nodes, coordinates );
  current_volume = v_wedge_volume( num_nodes, coordinates);
  distortion = jacobian*master_volume/current_volume/0.866025;

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

}

/*
   The stretch of a wedge element is here defined to be the maximum value of the
   stretch (S) of the three quadrilateral faces (see section 5.21):
   q = max[S_1043, S_1254, S_2035]

   Dimension : 1
   Acceptable Range :
   Normal Range :
   Full Range : [0,DBL_MAX]
   q for right, unit wedge : 1
   Reference : Adapted from section 5.21
   Verdict Function : v_wedge_max_stretch
   */

C_FUNC_DEF double v_wedge_max_stretch( int /*num_nodes*/, double coordinates[][3] )
{
  //This function finds the stretch of the 3 quadrilateral faces and returns the maximum value

  double stretch = 42, quad_face[4][3], stretch1 = 42, stretch2 = 42, stretch3 = 42;
  int i = 0;

  //first face
  for (i = 0; i < 3; i++){quad_face[0][i] = coordinates[0][i];}
  for (i = 0; i < 3; i++){quad_face[1][i] = coordinates[1][i];}
  for (i = 0; i < 3; i++){quad_face[2][i] = coordinates[4][i];}
  for (i = 0; i < 3; i++){quad_face[3][i] = coordinates[3][i];}
  stretch1 = v_quad_stretch(4,quad_face);

  //second face
  for (i = 0; i < 3; i++){quad_face[0][i] = coordinates[1][i];}
  for (i = 0; i < 3; i++){quad_face[1][i] = coordinates[2][i];}
  for (i = 0; i < 3; i++){quad_face[2][i] = coordinates[5][i];}
  for (i = 0; i < 3; i++){quad_face[3][i] = coordinates[4][i];}
  stretch2 = v_quad_stretch(4,quad_face);

  //third face
  for (i = 0; i < 3; i++){quad_face[0][i] = coordinates[2][i];}
  for (i = 0; i < 3; i++){quad_face[1][i] = coordinates[0][i];}
  for (i = 0; i < 3; i++){quad_face[2][i] = coordinates[3][i];}
  for (i = 0; i < 3; i++){quad_face[3][i] = coordinates[5][i];}
  stretch3 = v_quad_stretch(4,quad_face);

  stretch = VERDICT_MAX(stretch1,stretch2);
  stretch = VERDICT_MAX(stretch,stretch3);

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

}

/*
   This is the minimum determinant of the Jacobian matrix evaluated at each
   corner of the element, divided by the corresponding edge lengths and
   normalized to the unit wedge:
   q = min(  2 / sqrt(3) * ((L_2 X L_0) * L_3)_k / sqrt(mag(L_2) * mag(L_0) * mag(L_3)))
   where ((L_2 X L_0) * L_3)_k is the determinant of the Jacobian of the
   tetrahedron defined at the kth corner node, and L_2, L_0 and L_3 are the
   egdes defined according to the standard for tetrahedral elements.

   Dimension : 1
   Acceptable Range :
   Normal Range :
   Full Range : [?,DBL_MAX]
   q for right, unit wedge : 1
   Reference : Adapted from section 6.14 and 7.11
   Verdict Function : v_wedge_scaled_jacobian
   */

C_FUNC_DEF double v_wedge_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
{
  double min_jacobian = 0, current_jacobian = 0,lengths = 42;
  VerdictVector vec1,vec2,vec3;

  // Node 0
  vec1.set( coordinates[1][0] - coordinates[0][0],
      coordinates[1][1] - coordinates[0][1],
      coordinates[1][2] - coordinates[0][2] );

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

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

  lengths = sqrt(vec1.length_squared() * vec2.length_squared() * vec3.length_squared());

  current_jacobian = (vec2 % (vec1 * vec3));
  min_jacobian = current_jacobian/lengths;

  //node 1
  vec1.set( coordinates[2][0] - coordinates[1][0],
      coordinates[2][1] - coordinates[1][1],
      coordinates[2][2] - coordinates[1][2] );

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

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

  lengths = sqrt(vec1.length_squared() * vec2.length_squared() * vec3.length_squared());

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian/lengths,min_jacobian);

  //node 2
  vec1.set( coordinates[0][0] - coordinates[2][0],
      coordinates[0][1] - coordinates[2][1],
      coordinates[0][2] - coordinates[2][2] );

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

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

  lengths = sqrt(vec1.length_squared() * vec2.length_squared() * vec3.length_squared());

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian/lengths,min_jacobian);

  //node 3
  vec1.set( coordinates[0][0] - coordinates[3][0],
      coordinates[0][1] - coordinates[3][1],
      coordinates[0][2] - coordinates[3][2] );

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

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

  lengths = sqrt(vec1.length_squared() * vec2.length_squared() * vec3.length_squared());

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian/lengths,min_jacobian);

  //node 4
  vec1.set( coordinates[1][0] - coordinates[4][0],
      coordinates[1][1] - coordinates[4][1],
      coordinates[1][2] - coordinates[4][2] );

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

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

  lengths = sqrt(vec1.length_squared() * vec2.length_squared() * vec3.length_squared());

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian/lengths,min_jacobian);

  //node 5
  vec1.set( coordinates[3][0] - coordinates[5][0],
      coordinates[3][1] - coordinates[5][1],
      coordinates[3][2] - coordinates[5][2] );

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

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

  lengths = sqrt(vec1.length_squared() * vec2.length_squared() * vec3.length_squared());

  current_jacobian = vec2 % (vec1 * vec3);
  min_jacobian = VERDICT_MIN(current_jacobian/lengths,min_jacobian);

  min_jacobian *= 2/sqrt(3.0);

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

}

/*
   The shape metric is defined to be 3 divided by the minimum mean ratio of the
   Jacobian matrix evaluated at the element corners:
   q = 3 / min(i=0,1,...,6){ J_i ^ 2/3 / (mag(L_0) + mag(L_1) + mag(L_2) ) }
   where J_i is the Jacobian and L_0, L_1, L_2 are the sides of the tetrahedral
   formed at the ith corner.

   Dimension : 1
   Acceptable Range : [0.3,1]
   Normal Range : [0,1]
   Full Range : [0,1]
   q for right, unit wedge : 1
   Reference : Adapted from section 7.12
   Verdict Function : v_wedge_shape
   */

C_FUNC_DEF double v_wedge_shape( int /*num_nodes*/, double coordinates[][3] )
{
  double current_jacobian = 0, current_shape, norm_jacobi = 0;
  double min_shape = 1.0;
  static const double two_thirds = 2.0/3.0;
  VerdictVector vec1,vec2,vec3;

  // Node 0
  vec1.set( coordinates[1][0] - coordinates[0][0],
      coordinates[1][1] - coordinates[0][1],
      coordinates[1][2] - coordinates[0][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  norm_jacobi = current_jacobian*2.0/sqrt(3.0);
  current_shape = 3*pow(norm_jacobi,two_thirds)/(vec1.length_squared() + vec2.length_squared() + vec3.length_squared());
  min_shape = VERDICT_MIN(current_shape,min_shape);

  //node 1
  vec1.set( coordinates[2][0] - coordinates[1][0],
      coordinates[2][1] - coordinates[1][1],
      coordinates[2][2] - coordinates[1][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  norm_jacobi = current_jacobian*2.0/sqrt(3.0);
  current_shape = 3*pow(norm_jacobi,two_thirds)/(vec1.length_squared() + vec2.length_squared() + vec3.length_squared());
  min_shape = VERDICT_MIN(current_shape,min_shape);

  //node 2
  vec1.set( coordinates[0][0] - coordinates[2][0],
      coordinates[0][1] - coordinates[2][1],
      coordinates[0][2] - coordinates[2][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  norm_jacobi = current_jacobian*2.0/sqrt(3.0);
  current_shape = 3*pow(norm_jacobi,two_thirds)/(vec1.length_squared() + vec2.length_squared() + vec3.length_squared());
  min_shape = VERDICT_MIN(current_shape,min_shape);

  //node 3
  vec1.set( coordinates[0][0] - coordinates[3][0],
      coordinates[0][1] - coordinates[3][1],
      coordinates[0][2] - coordinates[3][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  norm_jacobi = current_jacobian*2.0/sqrt(3.0);
  current_shape = 3*pow(norm_jacobi,two_thirds)/(vec1.length_squared() + vec2.length_squared() + vec3.length_squared());
  min_shape = VERDICT_MIN(current_shape,min_shape);

  //node 4
  vec1.set( coordinates[1][0] - coordinates[4][0],
      coordinates[1][1] - coordinates[4][1],
      coordinates[1][2] - coordinates[4][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  norm_jacobi = current_jacobian*2.0/sqrt(3.0);
  current_shape = 3*pow(norm_jacobi,two_thirds)/(vec1.length_squared() + vec2.length_squared() + vec3.length_squared());
  min_shape = VERDICT_MIN(current_shape,min_shape);

  //node 5
  vec1.set( coordinates[3][0] - coordinates[5][0],
      coordinates[3][1] - coordinates[5][1],
      coordinates[3][2] - coordinates[5][2] );

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

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

  current_jacobian = vec2 % (vec1 * vec3);
  norm_jacobi = current_jacobian*2.0/sqrt(3.0);
  current_shape = 3*pow(norm_jacobi,two_thirds)/(vec1.length_squared() + vec2.length_squared() + vec3.length_squared());
  min_shape = VERDICT_MIN(current_shape,min_shape);

  if (min_shape < VERDICT_DBL_MIN)
    return 0;
  return min_shape;

}

/* For wedges, there is not a unique definition of the aspect Frobenius. Rather,
 * this metric uses the aspect Frobenius defined for tetrahedral (see section
 * 6.4) and is comparable in methodology to the maximum aspect Frobenius defined
 * for hexahedra (see section 7.7). This value is normalized for a unit wedge.

 q = max(F_0123, F_1204, F_2015, F_3540, F_4351, F_5432)

 This is also known as the wedge condition number.

 Dimension : 1
 Acceptable Range :
 Normal Range :
 Full Range :
 q for right, unit wedge : 1
 Reference : Adapted from section 7.7
 Verdict Function : v_wedge_max_aspect_frobenius or v_wedge_condition

 */
C_FUNC_DEF double v_wedge_condition( int /*num_nodes*/, double coordinates[][3] )
{
  return v_wedge_max_aspect_frobenius(6,coordinates);
}