 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
#include "moab/verdict.h"#include "verdict_defines.hpp"#include "V_GaussIntegration.hpp"#include "VerdictVector.hpp"#include <memory.h>#include <cstddef>
Include dependency graph for V_TriMetric.cpp:Go to the source code of this file.
Defines | |
| #define | VERDICT_EXPORTS |
Functions | |
| static int | v_tri_get_weight (double &m11, double &m21, double &m12, double &m22) |
| C_FUNC_DEF void | v_set_tri_size (double size) |
| Sets average size (area) of tri, needed for v_tri_relative_size(...) | |
| C_FUNC_DEF void | v_set_tri_normal_func (ComputeNormal func) |
| Sets fuction pointer to calculate tri normal wrt surface. | |
| C_FUNC_DEF double | v_tri_edge_ratio (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_aspect_ratio (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_radius_ratio (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_aspect_frobenius (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_area (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_minimum_angle (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_maximum_angle (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_condition (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_scaled_jacobian (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_shape (int num_nodes, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_relative_size_squared (int, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_shape_and_size (int num_nodes, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF double | v_tri_distortion (int num_nodes, double coordinates[][3]) |
| Calculates tri metric. | |
| C_FUNC_DEF void | v_tri_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, TriMetricVals *metric_vals) |
| Calculates quality metrics for triangle elements. | |
Variables | |
| static double | verdict_tri_size = 0 |
| static ComputeNormal | compute_normal = NULL |
| #define VERDICT_EXPORTS |
Definition at line 23 of file V_TriMetric.cpp.
| C_FUNC_DEF void v_set_tri_normal_func | ( | ComputeNormal | func | ) |
Sets fuction pointer to calculate tri normal wrt surface.
Definition at line 63 of file V_TriMetric.cpp.
References compute_normal.
{
compute_normal = func;
}
| C_FUNC_DEF void v_set_tri_size | ( | double | size | ) |
Sets average size (area) of tri, needed for v_tri_relative_size(...)
sets the average area of a tri
Definition at line 58 of file V_TriMetric.cpp.
References size, and verdict_tri_size.
Referenced by moab::VerdictWrapper::set_size().
{
verdict_tri_size = size;
}
| C_FUNC_DEF double v_tri_area | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
The area of a tri
0.5 * jacobian at a node
Definition at line 278 of file V_TriMetric.cpp.
References VerdictVector::length(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), v_quad_jacobian(), and v_quad_quality().
{
// 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 );
}
| C_FUNC_DEF double v_tri_aspect_frobenius | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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
Definition at line 246 of file V_TriMetric.cpp.
References VerdictVector::length_squared(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure().
{
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 );
}
| C_FUNC_DEF double v_tri_aspect_ratio | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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"
Definition at line 160 of file V_TriMetric.cpp.
References VerdictVector::length(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure().
{
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 );
}
}
| C_FUNC_DEF double v_tri_condition | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
The condition of a tri
Condition number of the jacobian matrix at any corner
Definition at line 421 of file V_TriMetric.cpp.
References compute_normal, VerdictVector::length(), VERDICT_DBL_MAX, VERDICT_MIN, VerdictVector::x(), VerdictVector::y(), and VerdictVector::z().
Referenced by moab::VerdictWrapper::quality_measure(), v_quad_condition(), and v_tri_shape().
{
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 );
}
| C_FUNC_DEF double v_tri_distortion | ( | int | num_nodes, |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
The distortion of a tri
TODO: make a short definition of the distortion and comment below
Definition at line 588 of file V_TriMetric.cpp.
References GaussIntegration::calculate_derivative_at_nodes_2d_tri(), GaussIntegration::calculate_shape_function_2d_tri(), dot_product(), GaussIntegration::get_shape_func(), GaussIntegration::initialize(), VerdictVector::length(), maxNumberNodes, maxTotalNumberGaussPoints, VerdictVector::normalize(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tri_quality().
{
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 );
}
| C_FUNC_DEF double v_tri_edge_ratio | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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
Definition at line 76 of file V_TriMetric.cpp.
References VerdictVector::length_squared(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tri_quality().
{
// 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 );
}
}
| static int v_tri_get_weight | ( | double & | m11, |
| double & | m21, | ||
| double & | m12, | ||
| double & | m22 | ||
| ) | [static] |
get weights based on the average area of a set of tris
Definition at line 40 of file V_TriMetric.cpp.
References verdict_tri_size.
Referenced by v_tri_quality(), and v_tri_relative_size_squared().
{
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;
}
| C_FUNC_DEF double v_tri_maximum_angle | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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.
Definition at line 361 of file V_TriMetric.cpp.
References VerdictVector::interior_angle(), VerdictVector::length_squared(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), v_quad_maximum_angle(), and v_quad_quality().
{
// 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 );
}
| C_FUNC_DEF double v_tri_minimum_angle | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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.
Definition at line 303 of file V_TriMetric.cpp.
References VerdictVector::interior_angle(), VerdictVector::length_squared(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), v_quad_minimum_angle(), and v_quad_quality().
{
// 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 );
}
| C_FUNC_DEF void v_tri_quality | ( | int | num_nodes, |
| double | coordinates[][3], | ||
| unsigned int | metrics_request_flag, | ||
| TriMetricVals * | metric_vals | ||
| ) |
Calculates quality metrics for triangle elements.
tri_quality for calculating multiple tri metrics at once
using this method is generally faster than using the individual method multiple times.
Definition at line 777 of file V_TriMetric.cpp.
References TriMetricVals::area, TriMetricVals::aspect_frobenius, TriMetricVals::aspect_ratio, compute_normal, TriMetricVals::condition, determinant(), TriMetricVals::distortion, TriMetricVals::edge_ratio, interior_angle(), length(), VerdictVector::length(), TriMetricVals::maximum_angle, TriMetricVals::minimum_angle, TriMetricVals::radius_ratio, TriMetricVals::relative_size_squared, TriMetricVals::scaled_jacobian, VerdictVector::set(), TriMetricVals::shape, TriMetricVals::shape_and_size, size, V_TRI_AREA, V_TRI_ASPECT_FROBENIUS, V_TRI_CONDITION, V_TRI_DISTORTION, v_tri_distortion(), V_TRI_EDGE_RATIO, v_tri_edge_ratio(), v_tri_get_weight(), V_TRI_MAXIMUM_ANGLE, V_TRI_MINIMUM_ANGLE, V_TRI_RADIUS_RATIO, v_tri_radius_ratio(), V_TRI_RELATIVE_SIZE_SQUARED, V_TRI_SCALED_JACOBIAN, V_TRI_SHAPE, V_TRI_SHAPE_AND_SIZE, VERDICT_DBL_MAX, VERDICT_MAX, VERDICT_MIN, VerdictVector::x(), VerdictVector::y(), and VerdictVector::z().
Referenced by moab::VerdictWrapper::all_quality_measures().
{
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 );
}
}
| C_FUNC_DEF double v_tri_radius_ratio | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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"
Definition at line 206 of file V_TriMetric.cpp.
References VerdictVector::length_squared(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tri_quality().
{
// 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 );
}
| C_FUNC_DEF double v_tri_relative_size_squared | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
The relative size of a tri
Min(J,1/J) where J is the determinant of the weighted jacobian matrix.
Definition at line 534 of file V_TriMetric.cpp.
References determinant(), VerdictVector::length(), VerdictVector::set(), size, v_tri_get_weight(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tri_shape_and_size().
{
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 );
}
| C_FUNC_DEF double v_tri_scaled_jacobian | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
The scaled jacobian of a tri
minimum of the jacobian divided by the lengths of 2 edge vectors
Definition at line 461 of file V_TriMetric.cpp.
References compute_normal, moab::cross(), moab::GeomUtil::first(), length(), VerdictVector::length(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, VERDICT_MIN, VerdictVector::x(), VerdictVector::y(), and VerdictVector::z().
Referenced by moab::VerdictWrapper::quality_measure(), v_quad_quality(), and v_quad_scaled_jacobian().
{
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 );
}
| C_FUNC_DEF double v_tri_shape | ( | int | num_nodes, |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
The shape of a tri
2 / condition number of weighted jacobian matrix
Definition at line 515 of file V_TriMetric.cpp.
References v_tri_condition(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tri_shape_and_size().
{
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 );
}
| C_FUNC_DEF double v_tri_shape_and_size | ( | int | num_nodes, |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
The shape and size of a tri
Product of the Shape and Relative Size
Definition at line 570 of file V_TriMetric.cpp.
References size, v_tri_relative_size_squared(), v_tri_shape(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure().
{
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 );
}
ComputeNormal compute_normal = NULL [static] |
Definition at line 34 of file V_TriMetric.cpp.
Referenced by v_set_tri_normal_func(), v_tri_condition(), v_tri_quality(), and v_tri_scaled_jacobian().
double verdict_tri_size = 0 [static] |
Definition at line 33 of file V_TriMetric.cpp.
Referenced by v_set_tri_size(), and v_tri_get_weight().
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.