V_TriMetric.cpp

Go to the documentation of this file.

/*=========================================================================
 
 Module: $RCSfile: V_TriMetric.cpp,v $
 
 Copyright (c) 2007 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.
 
=========================================================================*/
 
/*
 *
 * TriMetric.cpp contains quality calculations for Tris
 *
 * This file is part of VERDICT
 *
 */
 
#define VERDICT_EXPORTS
 
#include "moab/verdict.h"
#include "verdict_defines.hpp"
#include "V_GaussIntegration.hpp"
#include "VerdictVector.hpp"
#include <memory.h>
#include <cstddef>
 
// the average area of a tri
static double verdict_tri_size = 0;
static ComputeNormal compute_normal = NULL;
 
/*!
 get weights based on the average area of a set of
 tris
*/
static int v_tri_get_weight( double& m11, double& m21, double& m12, double& m22 )
{
 static const double rootOf3 = sqrt( 3.0 );
 m11 = 1;
 m21 = 0;
 m12 = 0.5;
 m22 = 0.5 * rootOf3;
 double scale = sqrt( 2.0 * verdict_tri_size / ( m11 * m22 - m21 * m12 ) );
 
 m11 *= scale;
 m21 *= scale;
 m12 *= scale;
 m22 *= scale;
 
 return 1;
}
 
/*! sets the average area of a tri */
C_FUNC_DEF void v_set_tri_size( double size )
{
 verdict_tri_size = size;
}
 
C_FUNC_DEF void v_set_tri_normal_func( ComputeNormal func )
{
 compute_normal = func;
}
 
/*!
 the edge ratio of a triangle
 
 NB (P. Pebay 01/14/07):
 Hmax / Hmin where Hmax and Hmin are respectively the maximum and the
 minimum edge lengths
 
*/
C_FUNC_DEF double v_tri_edge_ratio( int /*num_nodes*/, double coordinates[][3] )
{
 
 // three vectors for each side
 VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 
 VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
 coordinates[2][2] - coordinates[1][2] );
 
 VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
 coordinates[0][2] - coordinates[2][2] );
 
 double a2 = a.length_squared();
 double b2 = b.length_squared();
 double c2 = c.length_squared();
 
 double m2, M2;
 if( a2 < b2 )
 {
 if( b2 < c2 )
 {
 m2 = a2;
 M2 = c2;
 }
 else // b2 <= a2
 {
 if( a2 < c2 )
 {
 m2 = a2;
 M2 = b2;
 }
 else // c2 <= a2
 {
 m2 = c2;
 M2 = b2;
 }
 }
 }
 else // b2 <= a2
 {
 if( a2 < c2 )
 {
 m2 = b2;
 M2 = c2;
 }
 else // c2 <= a2
 {
 if( b2 < c2 )
 {
 m2 = b2;
 M2 = a2;
 }
 else // c2 <= b2
 {
 m2 = c2;
 M2 = a2;
 }
 }
 }
 
 if( m2 < VERDICT_DBL_MIN )
 return (double)VERDICT_DBL_MAX;
 else
 {
 double edge_ratio;
 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 aspect ratio of a triangle
 
 NB (P. Pebay 01/14/07):
 Hmax / ( 2.0 * sqrt(3.0) * IR) where Hmax is the maximum edge length
 and IR is the inradius
 
 note that previous incarnations of verdict used "v_tri_aspect_ratio" to denote
 what is now called "v_tri_aspect_frobenius"
 
*/
C_FUNC_DEF double v_tri_aspect_ratio( int /*num_nodes*/, double coordinates[][3] )
{
 static const double normal_coeff = sqrt( 3. ) / 6.;
 
 // three vectors for each side
 VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 
 VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
 coordinates[2][2] - coordinates[1][2] );
 
 VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
 coordinates[0][2] - coordinates[2][2] );
 
 double a1 = a.length();
 double b1 = b.length();
 double c1 = c.length();
 
 double hm = a1 > b1 ? a1 : b1;
 hm = hm > c1 ? hm : c1;
 
 VerdictVector ab = a * b;
 double denominator = ab.length();
 
 if( denominator < VERDICT_DBL_MIN )
 return (double)VERDICT_DBL_MAX;
 else
 {
 double aspect_ratio;
 aspect_ratio = normal_coeff * hm * ( a1 + b1 + c1 ) / denominator;
 
 if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
 }
}
 
/*!
 the radius ratio of a triangle
 
 NB (P. Pebay 01/13/07):
 CR / (3.0*IR) where CR is the circumradius and IR is the inradius
 
 this quality metric is also known to VERDICT, for tetrahedral elements only,
 a the "aspect beta"
 
*/
C_FUNC_DEF double v_tri_radius_ratio( int /*num_nodes*/, double coordinates[][3] )
{
 
 // three vectors for each side
 VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 
 VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
 coordinates[2][2] - coordinates[1][2] );
 
 VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
 coordinates[0][2] - coordinates[2][2] );
 
 double a2 = a.length_squared();
 double b2 = b.length_squared();
 double c2 = c.length_squared();
 
 VerdictVector ab = a * b;
 double denominator = ab.length_squared();
 
 if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
 
 double radius_ratio;
 radius_ratio = .25 * a2 * b2 * c2 * ( a2 + b2 + c2 ) / denominator;
 
 if( radius_ratio > 0 ) return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( radius_ratio, -VERDICT_DBL_MAX );
}
 
/*!
 the Frobenius aspect of a tri
 
 srms^2/(2 * sqrt(3.0) * area)
 where srms^2 is sum of the lengths squared
 
 NB (P. Pebay 01/14/07):
 this method was called "aspect ratio" in earlier incarnations of VERDICT
 
*/
 
C_FUNC_DEF double v_tri_aspect_frobenius( int /*num_nodes*/, double coordinates[][3] )
{
 static const double two_times_root_of_3 = 2 * sqrt( 3.0 );
 
 // three vectors for each side
 VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 
 VerdictVector side2( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
 coordinates[2][2] - coordinates[1][2] );
 
 VerdictVector side3( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
 coordinates[0][2] - coordinates[2][2] );
 
 // sum the lengths squared of each side
 double srms = ( side1.length_squared() + side2.length_squared() + side3.length_squared() );
 
 // find two times the area of the triangle by cross product
 double areaX2 = ( ( side1 * ( -side3 ) ).length() );
 
 if( areaX2 == 0.0 ) return (double)VERDICT_DBL_MAX;
 
 double aspect = (double)( srms / ( two_times_root_of_3 * ( areaX2 ) ) );
 if( aspect > 0 ) return (double)VERDICT_MIN( aspect, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( aspect, -VERDICT_DBL_MAX );
}
 
/*!
 The area of a tri
 
 0.5 * jacobian at a node
*/
C_FUNC_DEF double v_tri_area( int /*num_nodes*/, double coordinates[][3] )
{
 // two vectors for two sides
 VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 
 VerdictVector side3( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
 coordinates[2][2] - coordinates[0][2] );
 
 // the cross product of the two vectors representing two sides of the
 // triangle
 VerdictVector tmp = side1 * side3;
 
 // return the magnitude of the vector divided by two
 double area = 0.5 * tmp.length();
 if( area > 0 ) return (double)VERDICT_MIN( area, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( area, -VERDICT_DBL_MAX );
}
 
/*!
 The minimum angle of a tri
 
 The minimum angle of a tri is the minimum angle between
 two adjacents sides out of all three corners of the triangle.
*/
C_FUNC_DEF double v_tri_minimum_angle( int /*num_nodes*/, double coordinates[][3] )
{
 
 // vectors for all the sides
 VerdictVector sides[4];
 sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
 coordinates[2][2] - coordinates[1][2] );
 sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
 coordinates[2][2] - coordinates[0][2] );
 
 // in case we need to find the interior angle
 // between sides 0 and 1
 sides[3] = -sides[1];
 
 // calculate the lengths squared of the sides
 double sides_lengths[3];
 sides_lengths[0] = sides[0].length_squared();
 sides_lengths[1] = sides[1].length_squared();
 sides_lengths[2] = sides[2].length_squared();
 
 if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 ) return 0.0;
 
 // using the law of sines, we know that the minimum
 // angle is opposite of the shortest side
 
 // find the shortest side
 int short_side = 0;
 if( sides_lengths[1] < sides_lengths[0] ) short_side = 1;
 if( sides_lengths[2] < sides_lengths[short_side] ) short_side = 2;
 
 // from the shortest side, calculate the angle of the
 // opposite angle
 double min_angle;
 if( short_side == 0 )
 {
 min_angle = sides[2].interior_angle( sides[1] );
 }
 else if( short_side == 1 )
 {
 min_angle = sides[0].interior_angle( sides[2] );
 }
 else
 {
 min_angle = sides[0].interior_angle( sides[3] );
 }
 
 if( min_angle > 0 ) return (double)VERDICT_MIN( min_angle, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( min_angle, -VERDICT_DBL_MAX );
}
 
/*!
 The maximum angle of a tri
 
 The maximum angle of a tri is the maximum angle between
 two adjacents sides out of all three corners of the triangle.
*/
C_FUNC_DEF double v_tri_maximum_angle( int /*num_nodes*/, double coordinates[][3] )
{
 
 // vectors for all the sides
 VerdictVector sides[4];
 sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
 coordinates[2][2] - coordinates[1][2] );
 sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
 coordinates[2][2] - coordinates[0][2] );
 
 // in case we need to find the interior angle
 // between sides 0 and 1
 sides[3] = -sides[1];
 
 // calculate the lengths squared of the sides
 double sides_lengths[3];
 sides_lengths[0] = sides[0].length_squared();
 sides_lengths[1] = sides[1].length_squared();
 sides_lengths[2] = sides[2].length_squared();
 
 if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 )
 {
 return 0.0;
 }
 
 // using the law of sines, we know that the maximum
 // angle is opposite of the longest side
 
 // find the longest side
 int short_side = 0;
 if( sides_lengths[1] > sides_lengths[0] ) short_side = 1;
 if( sides_lengths[2] > sides_lengths[short_side] ) short_side = 2;
 
 // from the longest side, calculate the angle of the
 // opposite angle
 double max_angle;
 if( short_side == 0 )
 {
 max_angle = sides[2].interior_angle( sides[1] );
 }
 else if( short_side == 1 )
 {
 max_angle = sides[0].interior_angle( sides[2] );
 }
 else
 {
 max_angle = sides[0].interior_angle( sides[3] );
 }
 
 if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
}
 
/*!
 The condition of a tri
 
 Condition number of the jacobian matrix at any corner
*/
C_FUNC_DEF double v_tri_condition( int /*num_nodes*/, double coordinates[][3] )
{
 static const double rt3 = sqrt( 3.0 );
 
 VerdictVector v1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 
 VerdictVector v2( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
 coordinates[2][2] - coordinates[0][2] );
 
 VerdictVector tri_normal = v1 * v2;
 double areax2 = tri_normal.length();
 
 if( areax2 == 0.0 ) return (double)VERDICT_DBL_MAX;
 
 double condition = (double)( ( ( v1 % v1 ) + ( v2 % v2 ) - ( v1 % v2 ) ) / ( areax2 * rt3 ) );
 
 // check for inverted if we have access to the normal
 if( compute_normal )
 {
 // center of tri
 double point[3], surf_normal[3];
 point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
 point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
 point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
 
 // dot product
 compute_normal( point, surf_normal );
 if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
 0 )
 return (double)VERDICT_DBL_MAX;
 }
 return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
}
 
/*!
 The scaled jacobian of a tri
 
 minimum of the jacobian divided by the lengths of 2 edge vectors
*/
C_FUNC_DEF double v_tri_scaled_jacobian( int /*num_nodes*/, double coordinates[][3] )
{
 static const double detw = 2. / sqrt( 3.0 );
 VerdictVector first, second;
 double jacobian;
 
 VerdictVector edge[3];
 edge[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 
 edge[1].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
 coordinates[2][2] - coordinates[0][2] );
 
 edge[2].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
 coordinates[2][2] - coordinates[1][2] );
 first = edge[1] - edge[0];
 second = edge[2] - edge[0];
 
 VerdictVector cross = first * second;
 jacobian = cross.length();
 
 double max_edge_length_product;
 max_edge_length_product =
 VERDICT_MAX( edge[0].length() * edge[1].length(),
 VERDICT_MAX( edge[1].length() * edge[2].length(), edge[0].length() * edge[2].length() ) );
 
 if( max_edge_length_product < VERDICT_DBL_MIN ) return (double)0.0;
 
 jacobian *= detw;
 jacobian /= max_edge_length_product;
 
 if( compute_normal )
 {
 // center of tri
 double point[3], surf_normal[3];
 point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
 point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
 point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
 
 // dot product
 compute_normal( point, surf_normal );
 if( ( cross.x() * surf_normal[0] + cross.y() * surf_normal[1] + cross.z() * surf_normal[2] ) < 0 )
 jacobian *= -1;
 }
 
 if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
}
 
/*!
 The shape of a tri
 
 2 / condition number of weighted jacobian matrix
*/
C_FUNC_DEF double v_tri_shape( int num_nodes, double coordinates[][3] )
{
 double condition = v_tri_condition( num_nodes, coordinates );
 
 double shape;
 if( condition <= VERDICT_DBL_MIN )
 shape = VERDICT_DBL_MAX;
 else
 shape = ( 1 / condition );
 
 if( shape > 0 ) return (double)VERDICT_MIN( shape, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( shape, -VERDICT_DBL_MAX );
}
 
/*!
 The relative size of a tri
 
 Min(J,1/J) where J is the determinant of the weighted jacobian matrix.
*/
C_FUNC_DEF double v_tri_relative_size_squared( int /*num_nodes*/, double coordinates[][3] )
{
 double w11, w21, w12, w22;
 
 VerdictVector xxi, xet, tri_normal;
 
 v_tri_get_weight( w11, w21, w12, w22 );
 
 double detw = determinant( w11, w21, w12, w22 );
 
 if( detw == 0.0 ) return 0.0;
 
 xxi.set( coordinates[0][0] - coordinates[1][0], coordinates[0][1] - coordinates[1][1],
 coordinates[0][2] - coordinates[1][2] );
 
 xet.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
 coordinates[0][2] - coordinates[2][2] );
 
 tri_normal = xxi * xet;
 
 double deta = tri_normal.length();
 if( deta == 0.0 || detw == 0.0 ) return 0.0;
 
 double size = pow( deta / detw, 2 );
 
 double rel_size = VERDICT_MIN( size, 1.0 / size );
 
 if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
}
 
/*!
 The shape and size of a tri
 
 Product of the Shape and Relative Size
*/
C_FUNC_DEF double v_tri_shape_and_size( int num_nodes, double coordinates[][3] )
{
 double size, shape;
 
 size = v_tri_relative_size_squared( num_nodes, coordinates );
 shape = v_tri_shape( num_nodes, coordinates );
 
 double shape_and_size = size * shape;
 
 if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
}
 
/*!
 The distortion of a tri
 
TODO: make a short definition of the distortion and comment below
*/
C_FUNC_DEF double v_tri_distortion( int num_nodes, double coordinates[][3] )
{
 
 double distortion;
 int total_number_of_gauss_points = 0;
 VerdictVector aa, bb, cc, normal_at_point, xin;
 double element_area = 0.;
 
 aa.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 
 bb.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
 coordinates[2][2] - coordinates[0][2] );
 
 VerdictVector tri_normal = aa * bb;
 
 int number_of_gauss_points = 0;
 if( num_nodes == 3 )
 {
 distortion = 1.0;
 return (double)distortion;
 }
 
 else if( num_nodes == 6 )
 {
 total_number_of_gauss_points = 6;
 number_of_gauss_points = 6;
 }
 
 distortion = VERDICT_DBL_MAX;
 double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
 double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
 double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
 double weight[maxTotalNumberGaussPoints];
 
 // create an object of GaussIntegration
 int number_dims = 2;
 int is_tri = 1;
 GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
 GaussIntegration::calculate_shape_function_2d_tri();
 GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
 
 // calculate element area
 int ife, ja;
 for( ife = 0; ife < total_number_of_gauss_points; ife++ )
 {
 aa.set( 0.0, 0.0, 0.0 );
 bb.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] );
 aa += dndy1[ife][ja] * xin;
 bb += dndy2[ife][ja] * xin;
 }
 normal_at_point = aa * bb;
 double jacobian = normal_at_point.length();
 element_area += weight[ife] * jacobian;
 }
 
 element_area *= 0.8660254;
 double dndy1_at_node[maxNumberNodes][maxNumberNodes];
 double dndy2_at_node[maxNumberNodes][maxNumberNodes];
 
 GaussIntegration::calculate_derivative_at_nodes_2d_tri( dndy1_at_node, dndy2_at_node );
 
 VerdictVector normal_at_nodes[7];
 
 // evaluate normal at nodes and distortion values at nodes
 int jai = 0;
 for( ja = 0; ja < num_nodes; ja++ )
 {
 aa.set( 0.0, 0.0, 0.0 );
 bb.set( 0.0, 0.0, 0.0 );
 for( jai = 0; jai < num_nodes; jai++ )
 {
 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
 aa += dndy1_at_node[ja][jai] * xin;
 bb += dndy2_at_node[ja][jai] * xin;
 }
 normal_at_nodes[ja] = aa * bb;
 normal_at_nodes[ja].normalize();
 }
 
 // determine if element is flat
 bool flat_element = true;
 double dot_product;
 
 for( ja = 0; ja < num_nodes; ja++ )
 {
 dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
 if( fabs( dot_product ) < 0.99 )
 {
 flat_element = false;
 break;
 }
 }
 
 // take into consideration of the thickness of the element
 double thickness, thickness_gauss;
 double distrt;
 // get_tri_thickness(tri, element_area, thickness );
 thickness = 0.001 * sqrt( element_area );
 
 // set thickness gauss point location
 double zl = 0.5773502691896;
 if( flat_element ) zl = 0.0;
 
 int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
 double thickness_z;
 
 // loop on integration points
 int igz;
 for( ife = 0; ife < total_number_of_gauss_points; ife++ )
 {
 // loop on the thickness direction gauss points
 for( igz = 0; igz < no_gauss_pts_z; igz++ )
 {
 zl = -zl;
 thickness_z = zl * thickness / 2.0;
 
 aa.set( 0.0, 0.0, 0.0 );
 bb.set( 0.0, 0.0, 0.0 );
 cc.set( 0.0, 0.0, 0.0 );
 
 for( ja = 0; ja < num_nodes; ja++ )
 {
 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
 xin += thickness_z * normal_at_nodes[ja];
 aa += dndy1[ife][ja] * xin;
 bb += dndy2[ife][ja] * xin;
 thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
 cc += thickness_gauss * normal_at_nodes[ja];
 }
 
 normal_at_point = aa * bb;
 distrt = cc % normal_at_point;
 if( distrt < distortion ) distortion = distrt;
 }
 }
 
 // loop through nodal points
 for( ja = 0; ja < num_nodes; ja++ )
 {
 for( igz = 0; igz < no_gauss_pts_z; igz++ )
 {
 zl = -zl;
 thickness_z = zl * thickness / 2.0;
 
 aa.set( 0.0, 0.0, 0.0 );
 bb.set( 0.0, 0.0, 0.0 );
 cc.set( 0.0, 0.0, 0.0 );
 
 for( jai = 0; jai < num_nodes; jai++ )
 {
 xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
 xin += thickness_z * normal_at_nodes[ja];
 aa += dndy1_at_node[ja][jai] * xin;
 bb += dndy2_at_node[ja][jai] * xin;
 if( jai == ja )
 thickness_gauss = thickness / 2.0;
 else
 thickness_gauss = 0.;
 cc += thickness_gauss * normal_at_nodes[jai];
 }
 }
 
 normal_at_point = aa * bb;
 double sign_jacobian = ( tri_normal % normal_at_point ) > 0 ? 1. : -1.;
 distrt = sign_jacobian * ( cc % normal_at_point );
 
 if( distrt < distortion ) distortion = distrt;
 }
 if( element_area * thickness != 0 )
 distortion *= 1. / ( element_area * thickness );
 else
 distortion *= 1.;
 
 if( distortion > 0 ) return (double)VERDICT_MIN( distortion, VERDICT_DBL_MAX );
 return (double)VERDICT_MAX( distortion, -VERDICT_DBL_MAX );
}
 
/*!
 tri_quality for calculating multiple tri metrics at once
 
 using this method is generally faster than using the individual
 method multiple times.
 
*/
C_FUNC_DEF void v_tri_quality( int num_nodes,
 double coordinates[][3],
 unsigned int metrics_request_flag,
 TriMetricVals* metric_vals )
{
 
 memset( metric_vals, 0, sizeof( TriMetricVals ) );
 
 // for starts, lets set up some basic and common information
 
 /* node numbers and side numbers used below
 
 2
 ++
 / \
 2 / \ 1
 / \
 / \
 0 ---------+ 1
 0
 */
 
 // vectors for each side
 VerdictVector sides[3];
 sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
 coordinates[1][2] - coordinates[0][2] );
 sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
 coordinates[2][2] - coordinates[1][2] );
 sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
 coordinates[2][2] - coordinates[0][2] );
 VerdictVector tri_normal = sides[0] * sides[2];
 // if we have access to normal information, check to see if the
 // element is inverted. If we don't have the normal information
 // that we need for this, assume the element is not inverted.
 // This flag will be used for condition number, jacobian, shape,
 // and size and shape.
 bool is_inverted = false;
 if( compute_normal )
 {
 // center of tri
 double point[3], surf_normal[3];
 point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
 point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
 point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
 // dot product
 compute_normal( point, surf_normal );
 if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
 0 )
 is_inverted = true;
 }
 
 // lengths squared of each side
 double sides_lengths_squared[3];
 sides_lengths_squared[0] = sides[0].length_squared();
 sides_lengths_squared[1] = sides[1].length_squared();
 sides_lengths_squared[2] = sides[2].length_squared();
 
 // if we are doing angle calcuations
 if( metrics_request_flag & ( V_TRI_MINIMUM_ANGLE | V_TRI_MAXIMUM_ANGLE ) )
 {
 // which is short and long side
 int short_side = 0, long_side = 0;
 
 if( sides_lengths_squared[1] < sides_lengths_squared[0] ) short_side = 1;
 if( sides_lengths_squared[2] < sides_lengths_squared[short_side] ) short_side = 2;
 
 if( sides_lengths_squared[1] > sides_lengths_squared[0] ) long_side = 1;
 if( sides_lengths_squared[2] > sides_lengths_squared[long_side] ) long_side = 2;
 
 // calculate the minimum angle of the tri
 if( metrics_request_flag & V_TRI_MINIMUM_ANGLE )
 {
 if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
 {
 metric_vals->minimum_angle = 0.0;
 }
 else if( short_side == 0 )
 metric_vals->minimum_angle = (double)sides[2].interior_angle( sides[1] );
 else if( short_side == 1 )
 metric_vals->minimum_angle = (double)sides[0].interior_angle( sides[2] );
 else
 metric_vals->minimum_angle = (double)sides[0].interior_angle( -sides[1] );
 }
 
 // calculate the maximum angle of the tri
 if( metrics_request_flag & V_TRI_MAXIMUM_ANGLE )
 {
 if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
 {
 metric_vals->maximum_angle = 0.0;
 }
 else if( long_side == 0 )
 metric_vals->maximum_angle = (double)sides[2].interior_angle( sides[1] );
 else if( long_side == 1 )
 metric_vals->maximum_angle = (double)sides[0].interior_angle( sides[2] );
 else
 metric_vals->maximum_angle = (double)sides[0].interior_angle( -sides[1] );
 }
 }
 
 // calculate the area of the tri
 // the following metrics depend on area
 if( metrics_request_flag &
 ( V_TRI_AREA | V_TRI_SCALED_JACOBIAN | V_TRI_SHAPE | V_TRI_RELATIVE_SIZE_SQUARED | V_TRI_SHAPE_AND_SIZE ) )
 {
 metric_vals->area = (double)( ( sides[0] * sides[2] ).length() * 0.5 );
 }
 
 // calculate the aspect ratio
 if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )
 {
 // sum the lengths squared
 double srms = sides_lengths_squared[0] + sides_lengths_squared[1] + sides_lengths_squared[2];
 
 // calculate once and reuse
 static const double twoTimesRootOf3 = 2 * sqrt( 3.0 );
 
 double div = ( twoTimesRootOf3 * ( ( sides[0] * sides[2] ).length() ) );
 
 if( div == 0.0 )
 metric_vals->aspect_frobenius = (double)VERDICT_DBL_MAX;
 else
 metric_vals->aspect_frobenius = (double)( srms / div );
 }
 
 // calculate the scaled jacobian
 if( metrics_request_flag & V_TRI_SCALED_JACOBIAN )
 {
 // calculate once and reuse
 static const double twoOverRootOf3 = 2 / sqrt( 3.0 );
 // use the area from above
 
 double tmp = tri_normal.length() * twoOverRootOf3;
 
 // now scale it by the lengths of the sides
 double min_scaled_jac = VERDICT_DBL_MAX;
 double temp_scaled_jac;
 for( int i = 0; i < 3; i++ )
 {
 if( sides_lengths_squared[i % 3] == 0.0 || sides_lengths_squared[( i + 2 ) % 3] == 0.0 )
 temp_scaled_jac = 0.0;
 else
 temp_scaled_jac =
 tmp / sqrt( sides_lengths_squared[i % 3] ) / sqrt( sides_lengths_squared[( i + 2 ) % 3] );
 if( temp_scaled_jac < min_scaled_jac ) min_scaled_jac = temp_scaled_jac;
 }
 // multiply by -1 if the normals are in opposite directions
 if( is_inverted )
 {
 min_scaled_jac *= -1;
 }
 metric_vals->scaled_jacobian = (double)min_scaled_jac;
 }
 
 // calculate the condition number
 if( metrics_request_flag & V_TRI_CONDITION )
 {
 // calculate once and reuse
 static const double rootOf3 = sqrt( 3.0 );
 // if it is inverted, the condition number is considered to be infinity.
 if( is_inverted )
 {
 metric_vals->condition = VERDICT_DBL_MAX;
 }
 else
 {
 double area2x = ( sides[0] * sides[2] ).length();
 if( area2x == 0.0 )
 metric_vals->condition = (double)( VERDICT_DBL_MAX );
 else
 metric_vals->condition = (double)( ( sides[0] % sides[0] + sides[2] % sides[2] - sides[0] % sides[2] ) /
 ( area2x * rootOf3 ) );
 }
 }
 
 // calculate the shape
 if( metrics_request_flag & V_TRI_SHAPE || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
 {
 // if element is inverted, shape is zero. We don't need to
 // calculate anything.
 if( is_inverted )
 {
 metric_vals->shape = 0.0;
 }
 else
 { // otherwise, we calculate the shape
 // calculate once and reuse
 static const double rootOf3 = sqrt( 3.0 );
 // reuse area from before
 double area2x = metric_vals->area * 2;
 // dot products
 double dots[3] = { sides[0] % sides[0], sides[2] % sides[2], sides[0] % sides[2] };
 
 // add the dots
 double sum_dots = dots[0] + dots[1] - dots[2];
 
 // then the finale
 if( sum_dots == 0.0 )
 metric_vals->shape = 0.0;
 else
 metric_vals->shape = (double)( rootOf3 * area2x / sum_dots );
 }
 }
 
 // calculate relative size squared
 if( metrics_request_flag & V_TRI_RELATIVE_SIZE_SQUARED || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
 {
 // get weights
 double w11, w21, w12, w22;
 v_tri_get_weight( w11, w21, w12, w22 );
 // get the determinant
 double detw = determinant( w11, w21, w12, w22 );
 // use the area from above and divide with the determinant
 if( metric_vals->area == 0.0 || detw == 0.0 )
 metric_vals->relative_size_squared = 0.0;
 else
 {
 double size = metric_vals->area * 2.0 / detw;
 // square the size
 size *= size;
 // value ranges between 0 to 1
 metric_vals->relative_size_squared = (double)VERDICT_MIN( size, 1.0 / size );
 }
 }
 
 // calculate shape and size
 if( metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
 {
 metric_vals->shape_and_size = metric_vals->relative_size_squared * metric_vals->shape;
 }
 
 // calculate distortion
 if( metrics_request_flag & V_TRI_DISTORTION ) metric_vals->distortion = v_tri_distortion( num_nodes, coordinates );
 
 // take care of any over-flow problems
 if( metric_vals->aspect_frobenius > 0 )
 metric_vals->aspect_frobenius = (double)VERDICT_MIN( metric_vals->aspect_frobenius, VERDICT_DBL_MAX );
 else
 metric_vals->aspect_frobenius = (double)VERDICT_MAX( metric_vals->aspect_frobenius, -VERDICT_DBL_MAX );
 
 if( metric_vals->area > 0 )
 metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
 else
 metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
 
 if( metric_vals->minimum_angle > 0 )
 metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
 else
 metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
 
 if( metric_vals->maximum_angle > 0 )
 metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
 else
 metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
 
 if( metric_vals->condition > 0 )
 metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
 else
 metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
 
 if( metric_vals->shape > 0 )
 metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
 else
 metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
 
 if( metric_vals->scaled_jacobian > 0 )
 metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
 else
 metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
 
 if( metric_vals->relative_size_squared > 0 )
 metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
 else
 metric_vals->relative_size_squared =
 (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
 
 if( metric_vals->shape_and_size > 0 )
 metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
 else
 metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
 
 if( metric_vals->distortion > 0 )
 metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
 else
 metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
 
 if( metrics_request_flag & V_TRI_EDGE_RATIO )
 {
 metric_vals->edge_ratio = v_tri_edge_ratio( 3, coordinates );
 }
 if( metrics_request_flag & V_TRI_RADIUS_RATIO )
 {
 metric_vals->radius_ratio = v_tri_radius_ratio( 3, coordinates );
 }
 if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS ) // there is no V_TRI_ASPECT_RATIO !
 {
 metric_vals->aspect_ratio = v_tri_radius_ratio( 3, coordinates );
 }
}