V_TetMetric.cpp File Reference

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

Include dependency graph for V_TetMetric.cpp:

Go to the source code of this file.

Macros
#define	VERDICT_EXPORTS

Functions
C_FUNC_DEF void	v_set_tet_size (double size)
	Sets average size (volume) of tet, needed for v_tet_relative_size(...) More...
int	get_weight (VerdictVector &w1, VerdictVector &w2, VerdictVector &w3)
C_FUNC_DEF double	v_tet_edge_ratio (int, double coordinates[][3])
	Calculates tet edge ratio metric. More...
C_FUNC_DEF double	v_tet_scaled_jacobian (int, double coordinates[][3])
	Calculates tet scaled jacobian. More...
C_FUNC_DEF double	v_tet_radius_ratio (int, double coordinates[][3])
	Calculates tet radius ratio metric. More...
C_FUNC_DEF double	v_tet_aspect_beta (int, double coordinates[][3])
	Calculates the radius ratio metric of a positively oriented tet. More...
C_FUNC_DEF double	v_tet_aspect_ratio (int, double coordinates[][3])
	Calculates tet aspect ratio metric. More...
C_FUNC_DEF double	v_tet_aspect_gamma (int, double coordinates[][3])
	Calculates tet aspect gamma metric. More...
C_FUNC_DEF double	v_tet_aspect_frobenius (int, double coordinates[][3])
	Calculates tet aspect frobenius metric. More...
C_FUNC_DEF double	v_tet_minimum_angle (int, double coordinates[][3])
	Calculates tet minimum dihedral angle. More...
C_FUNC_DEF double	v_tet_collapse_ratio (int, double coordinates[][3])
	Calculates tet collapse ratio metric. More...
C_FUNC_DEF double	v_tet_volume (int, double coordinates[][3])
	Calculates tet volume. More...
C_FUNC_DEF double	v_tet_condition (int, double coordinates[][3])
	Calculates tet condition metric. More...
C_FUNC_DEF double	v_tet_jacobian (int, double coordinates[][3])
	Calculates tet jacobian. More...
C_FUNC_DEF double	v_tet_shape (int, double coordinates[][3])
	Calculates tet shape metric. More...
C_FUNC_DEF double	v_tet_relative_size_squared (int, double coordinates[][3])
	Calculates tet relative size metric. More...
C_FUNC_DEF double	v_tet_shape_and_size (int num_nodes, double coordinates[][3])
	Calculates tet shape-size metric. More...
C_FUNC_DEF double	v_tet_distortion (int num_nodes, double coordinates[][3])
	Calculates tet distortion metric. More...
C_FUNC_DEF void	v_tet_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, TetMetricVals *metric_vals)
	Calculates quality metrics for tetrahedral elements. More...

Variables
double	verdict_tet_size = 0
	the average volume of a tet More...

Macro Definition Documentation

VERDICT_EXPORTS

#define VERDICT_EXPORTS

Definition at line 23 of file V_TetMetric.cpp.

Function Documentation

get_weight()

int get_weight	(	VerdictVector &	w1,
		VerdictVector &	w2,
		VerdictVector &	w3
	)

get the weights based on the average size of a tet

Definition at line 46 of file V_TetMetric.cpp.

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

References determinant(), VerdictVector::set(), and verdict_tet_size.

Referenced by v_tet_quality(), and v_tet_relative_size_squared().

v_set_tet_size()

C_FUNC_DEF void v_set_tet_size

(

double

size

)

Sets average size (volume) of tet, needed for v_tet_relative_size(...)

set the average volume of a tet

Definition at line 37 of file V_TetMetric.cpp.

{
 verdict_tet_size = size;
}

References size, and verdict_tet_size.

Referenced by moab::VerdictWrapper::set_size().

v_tet_aspect_beta()

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

Calculates the radius ratio metric of a positively oriented tet.

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 metric is similar to the radius ratio of a tet, except that it returns VERDICT_DBL_MAX if the element has negative orientation.

Definition at line 265 of file V_TetMetric.cpp.

{
 
 // 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 );
 }
}

References VerdictVector::length(), length(), VerdictVector::length_squared(), length_squared(), VerdictVector::set(), v_tet_volume(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure().

v_tet_aspect_frobenius()

C_FUNC_DEF double v_tet_aspect_frobenius	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet aspect frobenius metric.

The aspect frobenius of a tet

NB (P. Pebay 01/22/07): Frobenius condition number when the reference element is regular

Definition at line 426 of file V_TetMetric.cpp.

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

References VerdictVector::get_xyz(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_quality().

v_tet_aspect_gamma()

C_FUNC_DEF double v_tet_aspect_gamma	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet aspect gamma metric.

the aspect gamma of a tet

srms^3 / (8.48528137423857*V) where srms = sqrt(sum(Si^2)/6), where Si is the edge length

Definition at line 381 of file V_TetMetric.cpp.

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

References VerdictVector::length_squared(), VerdictVector::set(), v_tet_volume(), VERDICT_DBL_MAX, and VERDICT_DBL_MIN.

Referenced by moab::VerdictWrapper::quality_measure().

v_tet_aspect_ratio()

C_FUNC_DEF double v_tet_aspect_ratio	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet aspect ratio metric.

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

Definition at line 318 of file V_TetMetric.cpp.

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

References VerdictVector::length(), VerdictVector::length_squared(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_quality().

v_tet_collapse_ratio()

C_FUNC_DEF double v_tet_collapse_ratio	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet collapse ratio metric.

The collapse ratio of a tet

Definition at line 526 of file V_TetMetric.cpp.

{
 // 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 );
}

References VerdictVector::length(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_quality().

v_tet_condition()

C_FUNC_DEF double v_tet_condition	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet condition metric.

the condition of a tet

condition number of the jacobian matrix at any corner

Definition at line 629 of file V_TetMetric.cpp.

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

References VerdictVector::set(), VERDICT_DBL_MAX, and VERDICT_DBL_MIN.

Referenced by moab::VerdictWrapper::quality_measure().

v_tet_distortion()

C_FUNC_DEF double v_tet_distortion	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet distortion metric.

the distortion of a tet

Definition at line 765 of file V_TetMetric.cpp.

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

References GaussIntegration::calculate_derivative_at_nodes_3d_tet(), GaussIntegration::calculate_shape_function_3d_tet(), GaussIntegration::get_shape_func(), GaussIntegration::initialize(), maxNumberNodes, maxTotalNumberGaussPoints, VerdictVector::set(), and VERDICT_DBL_MAX.

Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_quality().

v_tet_edge_ratio()

C_FUNC_DEF double v_tet_edge_ratio	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet edge ratio metric.

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

Definition at line 71 of file V_TetMetric.cpp.

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

References VerdictVector::length_squared(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure().

v_tet_jacobian()

C_FUNC_DEF double v_tet_jacobian	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet jacobian.

the jacobian of a tet

TODO

Definition at line 670 of file V_TetMetric.cpp.

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

References VerdictVector::set().

Referenced by moab::VerdictWrapper::quality_measure().

v_tet_minimum_angle()

C_FUNC_DEF double v_tet_minimum_angle	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet minimum dihedral angle.

The minimum angle of a tet

NB (P. Pebay 01/22/07): minimum nonoriented dihedral angle

Definition at line 475 of file V_TetMetric.cpp.

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

References VerdictVector::length(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_quality().

v_tet_quality()

C_FUNC_DEF void v_tet_quality	(	int	num_nodes,
		double	coordinates[][3],
		unsigned int	metrics_request_flag,
		TetMetricVals *	metric_vals
	)

Calculates quality metrics for tetrahedral elements.

the quality metrics of a tet

Definition at line 862 of file V_TetMetric.cpp.

{
 
 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 );
 }
 
 if( metrics_request_flag & V_TET_ASPECT_RATIO ) metric_vals->aspect_ratio = v_tet_aspect_ratio( 4, coordinates );
 
 if( metrics_request_flag & V_TET_ASPECT_FROBENIUS )
 metric_vals->aspect_frobenius = v_tet_aspect_frobenius( 4, coordinates );
 
 if( metrics_request_flag & V_TET_MINIMUM_ANGLE ) metric_vals->minimum_angle = v_tet_minimum_angle( 4, coordinates );
 
 if( metrics_request_flag & V_TET_COLLAPSE_RATIO )
 metric_vals->collapse_ratio = v_tet_collapse_ratio( 4, coordinates );
 
 if( metrics_request_flag & V_TET_RADIUS_RATIO ) metric_vals->radius_ratio = v_tet_radius_ratio( 4, coordinates );
}

References TetMetricVals::aspect_beta, TetMetricVals::aspect_frobenius, TetMetricVals::aspect_gamma, TetMetricVals::aspect_ratio, TetMetricVals::collapse_ratio, TetMetricVals::condition, TetMetricVals::distortion, get_weight(), TetMetricVals::jacobian, VerdictVector::length(), length(), VerdictVector::length_squared(), length_squared(), TetMetricVals::minimum_angle, TetMetricVals::radius_ratio, TetMetricVals::relative_size_squared, TetMetricVals::scaled_jacobian, VerdictVector::set(), TetMetricVals::shape, TetMetricVals::shape_and_size, V_TET_ASPECT_BETA, V_TET_ASPECT_FROBENIUS, v_tet_aspect_frobenius(), V_TET_ASPECT_GAMMA, V_TET_ASPECT_RATIO, v_tet_aspect_ratio(), V_TET_COLLAPSE_RATIO, v_tet_collapse_ratio(), V_TET_CONDITION, V_TET_DISTORTION, v_tet_distortion(), V_TET_JACOBIAN, V_TET_MINIMUM_ANGLE, v_tet_minimum_angle(), V_TET_RADIUS_RATIO, v_tet_radius_ratio(), V_TET_RELATIVE_SIZE_SQUARED, V_TET_SCALED_JACOBIAN, V_TET_SHAPE, V_TET_SHAPE_AND_SIZE, V_TET_VOLUME, VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, VERDICT_MIN, and TetMetricVals::volume.

Referenced by moab::VerdictWrapper::all_quality_measures().

v_tet_radius_ratio()

C_FUNC_DEF double v_tet_radius_ratio	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet radius ratio metric.

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 metric is similar to the aspect beta of a tet, except that it does not return VERDICT_DBL_MAX if the element has negative orientation.

Definition at line 209 of file V_TetMetric.cpp.

{
 
 // 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 );
 }
}

References VerdictVector::length(), length(), VerdictVector::length_squared(), length_squared(), VerdictVector::set(), v_tet_volume(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_quality().

v_tet_relative_size_squared()

C_FUNC_DEF double v_tet_relative_size_squared	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet relative size metric.

the relative size of a tet

Min(J,1/J), where J is the determinant of the weighted Jacobian matrix

Definition at line 727 of file V_TetMetric.cpp.

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

References get_weight(), size, v_tet_volume(), and VERDICT_DBL_MIN.

Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_shape_and_size().

v_tet_scaled_jacobian()

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

Calculates tet scaled jacobian.

the scaled jacobian of a tet

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

Definition at line 153 of file V_TetMetric.cpp.

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

References VerdictVector::length_squared(), length_squared(), VerdictVector::set(), VERDICT_DBL_MAX, and VERDICT_DBL_MIN.

Referenced by moab::VerdictWrapper::quality_measure().

v_tet_shape()

C_FUNC_DEF double v_tet_shape	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet shape metric.

the shape of a tet

3/ condition number of weighted jacobian matrix

Definition at line 691 of file V_TetMetric.cpp.

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

References VerdictVector::set(), VERDICT_DBL_MIN, and VERDICT_MAX.

Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_shape_and_size().

v_tet_shape_and_size()

C_FUNC_DEF double v_tet_shape_and_size	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet shape-size metric.

the shape and size of a tet

Product of the shape and relative size

Definition at line 752 of file V_TetMetric.cpp.

{
 
 double shape, size;
 shape = v_tet_shape( num_nodes, coordinates );
 size = v_tet_relative_size_squared( num_nodes, coordinates );
 
 return (double)( shape * size );
}

References size, v_tet_relative_size_squared(), and v_tet_shape().

Referenced by moab::VerdictWrapper::quality_measure().

v_tet_volume()

C_FUNC_DEF double v_tet_volume	(	int	num_nodes,
		double	coordinates[][3]
	)

Calculates tet volume.

the volume of a tet

1/6 * jacobian at a corner node

Definition at line 606 of file V_TetMetric.cpp.

{
 
 // 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 );
}

References VerdictVector::set().

Referenced by moab::VerdictWrapper::quality_measure(), v_tet_aspect_beta(), v_tet_aspect_gamma(), v_tet_radius_ratio(), and v_tet_relative_size_squared().

Variable Documentation

verdict_tet_size

double verdict_tet_size = 0

the average volume of a tet

Definition at line 32 of file V_TetMetric.cpp.

Referenced by get_weight(), and v_set_tet_size().