MOAB: Mesh Oriented datABase  (version 5.3.0)
V_HexMetric.cpp File Reference
#include "moab/verdict.h"
#include "VerdictVector.hpp"
#include "V_GaussIntegration.hpp"
#include "verdict_defines.hpp"
#include <memory.h>
+ Include dependency graph for V_HexMetric.cpp:

Go to the source code of this file.

Defines

#define VERDICT_EXPORTS
#define make_hex_nodes(coord, pos)
#define make_edge_length_squares(edges, lengths)
#define SQR(x)   ( ( x ) * ( x ) )

Functions

int v_hex_get_weight (VerdictVector &v1, VerdictVector &v2, VerdictVector &v3)
 weights based on the average size of a hex
C_FUNC_DEF void v_set_hex_size (double size)
 returns the average volume of a hex
void make_hex_edges (double coordinates[][3], VerdictVector edges[12])
 make VerdictVectors from coordinates
void localize_hex_coordinates (double coordinates[][3], VerdictVector position[8])
double safe_ratio3 (const double numerator, const double denominator, const double max_ratio)
double safe_ratio (const double numerator, const double denominator)
double condition_comp (const VerdictVector &xxi, const VerdictVector &xet, const VerdictVector &xze)
double oddy_comp (const VerdictVector &xxi, const VerdictVector &xet, const VerdictVector &xze)
double hex_edge_length (int max_min, double coordinates[][3])
 calcualates edge lengths of a hex
double diag_length (int max_min, double coordinates[][3])
VerdictVector calc_hex_efg (int efg_index, VerdictVector coordinates[8])
 calculates efg values
C_FUNC_DEF double v_hex_edge_ratio (int, double coordinates[][3])
 Calculates hex edge ratio metric.
C_FUNC_DEF double v_hex_max_edge_ratio (int, double coordinates[][3])
 Calculates hex maximum of edge ratio.
C_FUNC_DEF double v_hex_skew (int, double coordinates[][3])
 Calculates hex skew metric.
C_FUNC_DEF double v_hex_taper (int, double coordinates[][3])
 Calculates hex taper metric.
C_FUNC_DEF double v_hex_volume (int, double coordinates[][3])
 Calculates hex volume.
C_FUNC_DEF double v_hex_stretch (int, double coordinates[][3])
 Calculates hex stretch metric.
C_FUNC_DEF double v_hex_diagonal (int, double coordinates[][3])
 Calculates hex diagonal metric.
C_FUNC_DEF double v_hex_dimension (int, double coordinates[][3])
 Calculates hex dimension metric.
C_FUNC_DEF double v_hex_oddy (int, double coordinates[][3])
 Calculates hex oddy metric.
C_FUNC_DEF double v_hex_med_aspect_frobenius (int, double coordinates[][3])
 Calculates hex condition metric.
C_FUNC_DEF double v_hex_max_aspect_frobenius (int, double coordinates[][3])
 Calculates hex condition metric.
C_FUNC_DEF double v_hex_condition (int, double coordinates[][3])
C_FUNC_DEF double v_hex_jacobian (int, double coordinates[][3])
 Calculates hex jacobian metric.
C_FUNC_DEF double v_hex_scaled_jacobian (int, double coordinates[][3])
 Calculates hex scaled jacobian metric.
C_FUNC_DEF double v_hex_shear (int, double coordinates[][3])
 Calculates hex shear metric.
C_FUNC_DEF double v_hex_shape (int, double coordinates[][3])
 Calculates hex shape metric.
C_FUNC_DEF double v_hex_relative_size_squared (int, double coordinates[][3])
 Calculates hex relative size metric.
C_FUNC_DEF double v_hex_shape_and_size (int num_nodes, double coordinates[][3])
 Calculates hex shape-size metric.
C_FUNC_DEF double v_hex_shear_and_size (int num_nodes, double coordinates[][3])
 Calculates hex shear-size metric.
C_FUNC_DEF double v_hex_distortion (int num_nodes, double coordinates[][3])
 Calculates hex distortion metric.
C_FUNC_DEF void v_hex_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, HexMetricVals *metric_vals)
 Calculates quality metrics for hexahedral elements.

Variables

double verdict_hex_size = 0
 the average volume of a hex

Define Documentation

#define make_edge_length_squares (   edges,
  lengths 
)
Value:
{                                                               \
        for( int melii = 0; melii < 12; melii++ )                   \
            ( lengths )[melii] = ( edges )[melii].length_squared(); \
    }

Definition at line 68 of file V_HexMetric.cpp.

Referenced by v_hex_quality().

#define make_hex_nodes (   coord,
  pos 
)
Value:
for( int mhcii = 0; mhcii < 8; mhcii++ )                                                 \
    {                                                                                        \
        ( pos )[mhcii].set( ( coord )[mhcii][0], ( coord )[mhcii][1], ( coord )[mhcii][2] ); \
    }

Definition at line 62 of file V_HexMetric.cpp.

Referenced by v_hex_jacobian(), v_hex_max_aspect_frobenius(), v_hex_max_edge_ratio(), v_hex_med_aspect_frobenius(), v_hex_oddy(), v_hex_quality(), v_hex_relative_size_squared(), v_hex_scaled_jacobian(), v_hex_shape(), v_hex_shear(), v_hex_skew(), v_hex_taper(), and v_hex_volume().

#define SQR (   x)    ( ( x ) * ( x ) )

Definition at line 774 of file V_HexMetric.cpp.

Referenced by v_hex_dimension().

#define VERDICT_EXPORTS

Definition at line 23 of file V_HexMetric.cpp.


Function Documentation

VerdictVector calc_hex_efg ( int  efg_index,
VerdictVector  coordinates[8] 
)

calculates efg values

Definition at line 421 of file V_HexMetric.cpp.

References VerdictVector::set().

Referenced by v_hex_jacobian(), v_hex_max_aspect_frobenius(), v_hex_max_edge_ratio(), v_hex_oddy(), v_hex_quality(), v_hex_scaled_jacobian(), v_hex_skew(), v_hex_taper(), and v_hex_volume().

{

    VerdictVector efg;

    switch( efg_index )
    {

        case 1:
            efg = coordinates[1];
            efg += coordinates[2];
            efg += coordinates[5];
            efg += coordinates[6];
            efg -= coordinates[0];
            efg -= coordinates[3];
            efg -= coordinates[4];
            efg -= coordinates[7];
            break;

        case 2:
            efg = coordinates[2];
            efg += coordinates[3];
            efg += coordinates[6];
            efg += coordinates[7];
            efg -= coordinates[0];
            efg -= coordinates[1];
            efg -= coordinates[4];
            efg -= coordinates[5];
            break;

        case 3:
            efg = coordinates[4];
            efg += coordinates[5];
            efg += coordinates[6];
            efg += coordinates[7];
            efg -= coordinates[0];
            efg -= coordinates[1];
            efg -= coordinates[2];
            efg -= coordinates[3];
            break;

        case 12:
            efg = coordinates[0];
            efg += coordinates[2];
            efg += coordinates[4];
            efg += coordinates[6];
            efg -= coordinates[1];
            efg -= coordinates[3];
            efg -= coordinates[5];
            efg -= coordinates[7];
            break;

        case 13:
            efg = coordinates[0];
            efg += coordinates[3];
            efg += coordinates[5];
            efg += coordinates[6];
            efg -= coordinates[1];
            efg -= coordinates[2];
            efg -= coordinates[4];
            efg -= coordinates[7];
            break;

        case 23:
            efg = coordinates[0];
            efg += coordinates[1];
            efg += coordinates[6];
            efg += coordinates[7];
            efg -= coordinates[2];
            efg -= coordinates[3];
            efg -= coordinates[4];
            efg -= coordinates[5];
            break;

        case 123:
            efg = coordinates[0];
            efg += coordinates[2];
            efg += coordinates[5];
            efg += coordinates[7];
            efg -= coordinates[1];
            efg -= coordinates[5];
            efg -= coordinates[6];
            efg -= coordinates[2];
            break;

        default:
            efg.set( 0, 0, 0 );
    }

    return efg;
}
double condition_comp ( const VerdictVector xxi,
const VerdictVector xet,
const VerdictVector xze 
)

Definition at line 221 of file V_HexMetric.cpp.

References MBMesquite::det(), VERDICT_DBL_MAX, and VERDICT_DBL_MIN.

Referenced by v_hex_max_aspect_frobenius(), v_hex_med_aspect_frobenius(), and v_hex_quality().

{
    double det = xxi % ( xet * xze );

    if( det <= VERDICT_DBL_MIN ) { return VERDICT_DBL_MAX; }

    double term1 = xxi % xxi + xet % xet + xze % xze;
    double term2 =
        ( ( xxi * xet ) % ( xxi * xet ) ) + ( ( xet * xze ) % ( xet * xze ) ) + ( ( xze * xxi ) % ( xze * xxi ) );

    return sqrt( term1 * term2 ) / det;
}
double diag_length ( int  max_min,
double  coordinates[][3] 
)

Definition at line 371 of file V_HexMetric.cpp.

References VERDICT_MAX, and VERDICT_MIN.

Referenced by v_hex_diagonal(), v_hex_quality(), and v_hex_stretch().

{
    double temp[3], diag[4];
    int i;

    // lengths^2  f diag nals
    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[6][i] - coordinates[0][i];
        temp[i] = temp[i] * temp[i];
    }
    diag[0] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[4][i] - coordinates[2][i];
        temp[i] = temp[i] * temp[i];
    }
    diag[1] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[7][i] - coordinates[1][i];
        temp[i] = temp[i] * temp[i];
    }
    diag[2] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[5][i] - coordinates[3][i];
        temp[i] = temp[i] * temp[i];
    }
    diag[3] = sqrt( temp[0] + temp[1] + temp[2] );

    double diagonal = diag[0];
    if( max_min == 0 )  // Return min diagonal
    {
        for( i = 1; i < 4; i++ )
            diagonal = VERDICT_MIN( diagonal, diag[i] );
        return (double)diagonal;
    }
    else  // Return max diagonal
    {
        for( i = 1; i < 4; i++ )
            diagonal = VERDICT_MAX( diagonal, diag[i] );
        return (double)diagonal;
    }
}
double hex_edge_length ( int  max_min,
double  coordinates[][3] 
)

calcualates edge lengths of a hex

Definition at line 265 of file V_HexMetric.cpp.

References VERDICT_MAX, and VERDICT_MIN.

Referenced by v_hex_stretch().

{
    double temp[3], edge[12];
    int i;

    // lengths^2 of edges
    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[1][i] - coordinates[0][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[0] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[2][i] - coordinates[1][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[1] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[3][i] - coordinates[2][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[2] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[0][i] - coordinates[3][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[3] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[5][i] - coordinates[4][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[4] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[6][i] - coordinates[5][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[5] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[7][i] - coordinates[6][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[6] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[4][i] - coordinates[7][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[7] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[4][i] - coordinates[0][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[8] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[5][i] - coordinates[1][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[9] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[6][i] - coordinates[2][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[10] = sqrt( temp[0] + temp[1] + temp[2] );

    for( i = 0; i < 3; i++ )
    {
        temp[i] = coordinates[7][i] - coordinates[3][i];
        temp[i] = temp[i] * temp[i];
    }
    edge[11] = sqrt( temp[0] + temp[1] + temp[2] );

    double _edge = edge[0];

    if( max_min == 0 )
    {
        for( i = 1; i < 12; i++ )
            _edge = VERDICT_MIN( _edge, edge[i] );
        return (double)_edge;
    }
    else
    {
        for( i = 1; i < 12; i++ )
            _edge = VERDICT_MAX( _edge, edge[i] );
        return (double)_edge;
    }
}
void localize_hex_coordinates ( double  coordinates[][3],
VerdictVector  position[8] 
)

localizes hex coordinates

Definition at line 106 of file V_HexMetric.cpp.

References VerdictVector::set(), VerdictVector::x(), VerdictVector::y(), z, and VerdictVector::z().

{

    int ii;
    for( ii = 0; ii < 8; ii++ )
    {
        position[ii].set( coordinates[ii][0], coordinates[ii][1], coordinates[ii][2] );
    }

    // ... Make centroid of element the center of coordinate system
    VerdictVector point_1256 = position[1];
    point_1256 += position[2];
    point_1256 += position[5];
    point_1256 += position[6];

    VerdictVector point_0374 = position[0];
    point_0374 += position[3];
    point_0374 += position[7];
    point_0374 += position[4];

    VerdictVector centroid = point_1256;
    centroid += point_0374;
    centroid /= 8.0;

    int i;
    for( i = 0; i < 8; i++ )
        position[i] -= centroid;

    // ... Rotate element such that center of side 1-2 and 0-3 define X axis
    double DX = point_1256.x() - point_0374.x();
    double DY = point_1256.y() - point_0374.y();
    double DZ = point_1256.z() - point_0374.z();

    double AMAGX = sqrt( DX * DX + DZ * DZ );
    double AMAGY = sqrt( DX * DX + DY * DY + DZ * DZ );
    double FMAGX = AMAGX == 0.0 ? 1.0 : 0.0;
    double FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;

    double CZ = DX / ( AMAGX + FMAGX ) + FMAGX;
    double SZ = DZ / ( AMAGX + FMAGX );
    double CY = AMAGX / ( AMAGY + FMAGY ) + FMAGY;
    double SY = DY / ( AMAGY + FMAGY );

    double temp;

    for( i = 0; i < 8; i++ )
    {
        temp = CY * CZ * position[i].x() + CY * SZ * position[i].z() + SY * position[i].y();
        position[i].y( -SY * CZ * position[i].x() - SY * SZ * position[i].z() + CY * position[i].y() );
        position[i].z( -SZ * position[i].x() + CZ * position[i].z() );
        position[i].x( temp );
    }

    // ... Now, rotate about Y
    VerdictVector delta = -position[0];
    delta -= position[1];
    delta += position[2];
    delta += position[3];
    delta -= position[4];
    delta -= position[5];
    delta += position[6];
    delta += position[7];

    DY = delta.y();
    DZ = delta.z();

    AMAGY = sqrt( DY * DY + DZ * DZ );
    FMAGY = AMAGY == 0.0 ? 1.0 : 0.0;

    double CX = DY / ( AMAGY + FMAGY ) + FMAGY;
    double SX = DZ / ( AMAGY + FMAGY );

    for( i = 0; i < 8; i++ )
    {
        temp = CX * position[i].y() + SX * position[i].z();
        position[i].z( -SX * position[i].y() + CX * position[i].z() );
        position[i].y( temp );
    }
}
void make_hex_edges ( double  coordinates[][3],
VerdictVector  edges[12] 
)

make VerdictVectors from coordinates

Definition at line 75 of file V_HexMetric.cpp.

References VerdictVector::set().

Referenced by v_hex_edge_ratio(), and v_hex_quality().

{
    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[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
                  coordinates[3][2] - coordinates[2][2] );
    edges[3].set( coordinates[0][0] - coordinates[3][0], coordinates[0][1] - coordinates[3][1],
                  coordinates[0][2] - coordinates[3][2] );
    edges[4].set( coordinates[5][0] - coordinates[4][0], coordinates[5][1] - coordinates[4][1],
                  coordinates[5][2] - coordinates[4][2] );
    edges[5].set( coordinates[6][0] - coordinates[5][0], coordinates[6][1] - coordinates[5][1],
                  coordinates[6][2] - coordinates[5][2] );
    edges[6].set( coordinates[7][0] - coordinates[6][0], coordinates[7][1] - coordinates[6][1],
                  coordinates[7][2] - coordinates[6][2] );
    edges[7].set( coordinates[4][0] - coordinates[7][0], coordinates[4][1] - coordinates[7][1],
                  coordinates[4][2] - coordinates[7][2] );
    edges[8].set( coordinates[4][0] - coordinates[0][0], coordinates[4][1] - coordinates[0][1],
                  coordinates[4][2] - coordinates[0][2] );
    edges[9].set( coordinates[5][0] - coordinates[1][0], coordinates[5][1] - coordinates[1][1],
                  coordinates[5][2] - coordinates[1][2] );
    edges[10].set( coordinates[6][0] - coordinates[2][0], coordinates[6][1] - coordinates[2][1],
                   coordinates[6][2] - coordinates[2][2] );
    edges[11].set( coordinates[7][0] - coordinates[3][0], coordinates[7][1] - coordinates[3][1],
                   coordinates[7][2] - coordinates[3][2] );
}
double oddy_comp ( const VerdictVector xxi,
const VerdictVector xet,
const VerdictVector xze 
)

Definition at line 234 of file V_HexMetric.cpp.

References VERDICT_DBL_MAX, and VERDICT_DBL_MIN.

Referenced by v_hex_oddy(), and v_hex_quality().

{
    static const double third = 1.0 / 3.0;

    double g11, g12, g13, g22, g23, g33, rt_g;

    g11  = xxi % xxi;
    g12  = xxi % xet;
    g13  = xxi % xze;
    g22  = xet % xet;
    g23  = xet % xze;
    g33  = xze % xze;
    rt_g = xxi % ( xet * xze );

    double oddy_metric;
    if( rt_g > VERDICT_DBL_MIN )
    {
        double norm_G_squared = g11 * g11 + 2.0 * g12 * g12 + 2.0 * g13 * g13 + g22 * g22 + 2.0 * g23 * g23 + g33 * g33;

        double norm_J_squared = g11 + g22 + g33;

        oddy_metric = ( norm_G_squared - third * norm_J_squared * norm_J_squared ) / pow( rt_g, 4. * third );
    }

    else
        oddy_metric = VERDICT_DBL_MAX;

    return oddy_metric;
}
double safe_ratio ( const double  numerator,
const double  denominator 
)

Definition at line 206 of file V_HexMetric.cpp.

References VERDICT_DBL_MAX, and VERDICT_DBL_MIN.

Referenced by v_hex_diagonal(), v_hex_max_edge_ratio(), v_hex_quality(), v_hex_stretch(), and v_hex_taper().

{

    double return_value;
    const double filter_n = VERDICT_DBL_MAX;
    const double filter_d = VERDICT_DBL_MIN;
    if( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d ) { return_value = numerator / denominator; }
    else
    {
        return_value = VERDICT_DBL_MAX;
    }

    return return_value;
}
double safe_ratio3 ( const double  numerator,
const double  denominator,
const double  max_ratio 
)

Definition at line 186 of file V_HexMetric.cpp.

{
    // this filter is essential for good running time in practice
    double return_value;

    const double filter_n = max_ratio * 1.0e-16;
    const double filter_d = 1.0e-16;
    if( fabs( numerator ) <= filter_n && fabs( denominator ) >= filter_d ) { return_value = numerator / denominator; }
    else
    {
        return_value = fabs( numerator ) / max_ratio >= fabs( denominator )
                           ? ( ( numerator >= 0.0 && denominator >= 0.0 ) || ( numerator < 0.0 && denominator < 0.0 )
                                   ? max_ratio
                                   : -max_ratio )
                           : numerator / denominator;
    }

    return return_value;
}
C_FUNC_DEF double v_hex_condition ( int  ,
double  coordinates[][3] 
)

The maximum Frobenius condition of a hex, a.k.a. condition NB (P. Pebay 01/25/07): this method is maintained for backwards compatibility only. It will become deprecated at some point.

Definition at line 1308 of file V_HexMetric.cpp.

References v_hex_max_aspect_frobenius().

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

{

    return v_hex_max_aspect_frobenius( 8, coordinates );
}
C_FUNC_DEF double v_hex_diagonal ( int  ,
double  coordinates[][3] 
)

Calculates hex diagonal metric.

diagonal ratio of a hex

Minimum diagonal length / maximum diagonal length

Definition at line 762 of file V_HexMetric.cpp.

References diag_length(), safe_ratio(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{

    double min_diag = diag_length( 0, coordinates );
    double max_diag = diag_length( 1, coordinates );

    double diagonal = safe_ratio( min_diag, max_diag );

    if( diagonal > 0 ) return (double)VERDICT_MIN( diagonal, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( diagonal, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_dimension ( int  ,
double  coordinates[][3] 
)

Calculates hex dimension metric.

dimension of a hex

Pronto-specific characteristic length for stable time step calculation. Char_length = Volume / 2 grad Volume

Definition at line 782 of file V_HexMetric.cpp.

References SQR.

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

{

    double gradop[9][4];

    double x1 = coordinates[0][0];
    double x2 = coordinates[1][0];
    double x3 = coordinates[2][0];
    double x4 = coordinates[3][0];
    double x5 = coordinates[4][0];
    double x6 = coordinates[5][0];
    double x7 = coordinates[6][0];
    double x8 = coordinates[7][0];

    double y1 = coordinates[0][1];
    double y2 = coordinates[1][1];
    double y3 = coordinates[2][1];
    double y4 = coordinates[3][1];
    double y5 = coordinates[4][1];
    double y6 = coordinates[5][1];
    double y7 = coordinates[6][1];
    double y8 = coordinates[7][1];

    double z1 = coordinates[0][2];
    double z2 = coordinates[1][2];
    double z3 = coordinates[2][2];
    double z4 = coordinates[3][2];
    double z5 = coordinates[4][2];
    double z6 = coordinates[5][2];
    double z7 = coordinates[6][2];
    double z8 = coordinates[7][2];

    double z24 = z2 - z4;
    double z52 = z5 - z2;
    double z45 = z4 - z5;
    gradop[1][1] =
        ( y2 * ( z6 - z3 - z45 ) + y3 * z24 + y4 * ( z3 - z8 - z52 ) + y5 * ( z8 - z6 - z24 ) + y6 * z52 + y8 * z45 ) /
        12.0;

    double z31 = z3 - z1;
    double z63 = z6 - z3;
    double z16 = z1 - z6;
    gradop[2][1] =
        ( y3 * ( z7 - z4 - z16 ) + y4 * z31 + y1 * ( z4 - z5 - z63 ) + y6 * ( z5 - z7 - z31 ) + y7 * z63 + y5 * z16 ) /
        12.0;

    double z42 = z4 - z2;
    double z74 = z7 - z4;
    double z27 = z2 - z7;
    gradop[3][1] =
        ( y4 * ( z8 - z1 - z27 ) + y1 * z42 + y2 * ( z1 - z6 - z74 ) + y7 * ( z6 - z8 - z42 ) + y8 * z74 + y6 * z27 ) /
        12.0;

    double z13 = z1 - z3;
    double z81 = z8 - z1;
    double z38 = z3 - z8;
    gradop[4][1] =
        ( y1 * ( z5 - z2 - z38 ) + y2 * z13 + y3 * ( z2 - z7 - z81 ) + y8 * ( z7 - z5 - z13 ) + y5 * z81 + y7 * z38 ) /
        12.0;

    double z86 = z8 - z6;
    double z18 = z1 - z8;
    double z61 = z6 - z1;
    gradop[5][1] =
        ( y8 * ( z4 - z7 - z61 ) + y7 * z86 + y6 * ( z7 - z2 - z18 ) + y1 * ( z2 - z4 - z86 ) + y4 * z18 + y2 * z61 ) /
        12.0;

    double z57 = z5 - z7;
    double z25 = z2 - z5;
    double z72 = z7 - z2;
    gradop[6][1] =
        ( y5 * ( z1 - z8 - z72 ) + y8 * z57 + y7 * ( z8 - z3 - z25 ) + y2 * ( z3 - z1 - z57 ) + y1 * z25 + y3 * z72 ) /
        12.0;

    double z68 = z6 - z8;
    double z36 = z3 - z6;
    double z83 = z8 - z3;
    gradop[7][1] =
        ( y6 * ( z2 - z5 - z83 ) + y5 * z68 + y8 * ( z5 - z4 - z36 ) + y3 * ( z4 - z2 - z68 ) + y2 * z36 + y4 * z83 ) /
        12.0;

    double z75 = z7 - z5;
    double z47 = z4 - z7;
    double z54 = z5 - z4;
    gradop[8][1] =
        ( y7 * ( z3 - z6 - z54 ) + y6 * z75 + y5 * ( z6 - z1 - z47 ) + y4 * ( z1 - z3 - z75 ) + y3 * z47 + y1 * z54 ) /
        12.0;

    double x24 = x2 - x4;
    double x52 = x5 - x2;
    double x45 = x4 - x5;
    gradop[1][2] =
        ( z2 * ( x6 - x3 - x45 ) + z3 * x24 + z4 * ( x3 - x8 - x52 ) + z5 * ( x8 - x6 - x24 ) + z6 * x52 + z8 * x45 ) /
        12.0;

    double x31 = x3 - x1;
    double x63 = x6 - x3;
    double x16 = x1 - x6;
    gradop[2][2] =
        ( z3 * ( x7 - x4 - x16 ) + z4 * x31 + z1 * ( x4 - x5 - x63 ) + z6 * ( x5 - x7 - x31 ) + z7 * x63 + z5 * x16 ) /
        12.0;

    double x42 = x4 - x2;
    double x74 = x7 - x4;
    double x27 = x2 - x7;
    gradop[3][2] =
        ( z4 * ( x8 - x1 - x27 ) + z1 * x42 + z2 * ( x1 - x6 - x74 ) + z7 * ( x6 - x8 - x42 ) + z8 * x74 + z6 * x27 ) /
        12.0;

    double x13 = x1 - x3;
    double x81 = x8 - x1;
    double x38 = x3 - x8;
    gradop[4][2] =
        ( z1 * ( x5 - x2 - x38 ) + z2 * x13 + z3 * ( x2 - x7 - x81 ) + z8 * ( x7 - x5 - x13 ) + z5 * x81 + z7 * x38 ) /
        12.0;

    double x86 = x8 - x6;
    double x18 = x1 - x8;
    double x61 = x6 - x1;
    gradop[5][2] =
        ( z8 * ( x4 - x7 - x61 ) + z7 * x86 + z6 * ( x7 - x2 - x18 ) + z1 * ( x2 - x4 - x86 ) + z4 * x18 + z2 * x61 ) /
        12.0;

    double x57 = x5 - x7;
    double x25 = x2 - x5;
    double x72 = x7 - x2;
    gradop[6][2] =
        ( z5 * ( x1 - x8 - x72 ) + z8 * x57 + z7 * ( x8 - x3 - x25 ) + z2 * ( x3 - x1 - x57 ) + z1 * x25 + z3 * x72 ) /
        12.0;

    double x68 = x6 - x8;
    double x36 = x3 - x6;
    double x83 = x8 - x3;
    gradop[7][2] =
        ( z6 * ( x2 - x5 - x83 ) + z5 * x68 + z8 * ( x5 - x4 - x36 ) + z3 * ( x4 - x2 - x68 ) + z2 * x36 + z4 * x83 ) /
        12.0;

    double x75 = x7 - x5;
    double x47 = x4 - x7;
    double x54 = x5 - x4;
    gradop[8][2] =
        ( z7 * ( x3 - x6 - x54 ) + z6 * x75 + z5 * ( x6 - x1 - x47 ) + z4 * ( x1 - x3 - x75 ) + z3 * x47 + z1 * x54 ) /
        12.0;

    double y24 = y2 - y4;
    double y52 = y5 - y2;
    double y45 = y4 - y5;
    gradop[1][3] =
        ( x2 * ( y6 - y3 - y45 ) + x3 * y24 + x4 * ( y3 - y8 - y52 ) + x5 * ( y8 - y6 - y24 ) + x6 * y52 + x8 * y45 ) /
        12.0;

    double y31 = y3 - y1;
    double y63 = y6 - y3;
    double y16 = y1 - y6;
    gradop[2][3] =
        ( x3 * ( y7 - y4 - y16 ) + x4 * y31 + x1 * ( y4 - y5 - y63 ) + x6 * ( y5 - y7 - y31 ) + x7 * y63 + x5 * y16 ) /
        12.0;

    double y42 = y4 - y2;
    double y74 = y7 - y4;
    double y27 = y2 - y7;
    gradop[3][3] =
        ( x4 * ( y8 - y1 - y27 ) + x1 * y42 + x2 * ( y1 - y6 - y74 ) + x7 * ( y6 - y8 - y42 ) + x8 * y74 + x6 * y27 ) /
        12.0;

    double y13 = y1 - y3;
    double y81 = y8 - y1;
    double y38 = y3 - y8;
    gradop[4][3] =
        ( x1 * ( y5 - y2 - y38 ) + x2 * y13 + x3 * ( y2 - y7 - y81 ) + x8 * ( y7 - y5 - y13 ) + x5 * y81 + x7 * y38 ) /
        12.0;

    double y86 = y8 - y6;
    double y18 = y1 - y8;
    double y61 = y6 - y1;
    gradop[5][3] =
        ( x8 * ( y4 - y7 - y61 ) + x7 * y86 + x6 * ( y7 - y2 - y18 ) + x1 * ( y2 - y4 - y86 ) + x4 * y18 + x2 * y61 ) /
        12.0;

    double y57 = y5 - y7;
    double y25 = y2 - y5;
    double y72 = y7 - y2;
    gradop[6][3] =
        ( x5 * ( y1 - y8 - y72 ) + x8 * y57 + x7 * ( y8 - y3 - y25 ) + x2 * ( y3 - y1 - y57 ) + x1 * y25 + x3 * y72 ) /
        12.0;

    double y68 = y6 - y8;
    double y36 = y3 - y6;
    double y83 = y8 - y3;
    gradop[7][3] =
        ( x6 * ( y2 - y5 - y83 ) + x5 * y68 + x8 * ( y5 - y4 - y36 ) + x3 * ( y4 - y2 - y68 ) + x2 * y36 + x4 * y83 ) /
        12.0;

    double y75 = y7 - y5;
    double y47 = y4 - y7;
    double y54 = y5 - y4;
    gradop[8][3] =
        ( x7 * ( y3 - y6 - y54 ) + x6 * y75 + x5 * ( y6 - y1 - y47 ) + x4 * ( y1 - y3 - y75 ) + x3 * y47 + x1 * y54 ) /
        12.0;

    //     calculate element volume and characteristic element aspect ratio
    //     (used in time step and hourglass control) -

    double volume = coordinates[0][0] * gradop[1][1] + coordinates[1][0] * gradop[2][1] +
                    coordinates[2][0] * gradop[3][1] + coordinates[3][0] * gradop[4][1] +
                    coordinates[4][0] * gradop[5][1] + coordinates[5][0] * gradop[6][1] +
                    coordinates[6][0] * gradop[7][1] + coordinates[7][0] * gradop[8][1];
    double aspect =
        .5 * SQR( volume ) /
        ( SQR( gradop[1][1] ) + SQR( gradop[2][1] ) + SQR( gradop[3][1] ) + SQR( gradop[4][1] ) + SQR( gradop[5][1] ) +
          SQR( gradop[6][1] ) + SQR( gradop[7][1] ) + SQR( gradop[8][1] ) + SQR( gradop[1][2] ) + SQR( gradop[2][2] ) +
          SQR( gradop[3][2] ) + SQR( gradop[4][2] ) + SQR( gradop[5][2] ) + SQR( gradop[6][2] ) + SQR( gradop[7][2] ) +
          SQR( gradop[8][2] ) + SQR( gradop[1][3] ) + SQR( gradop[2][3] ) + SQR( gradop[3][3] ) + SQR( gradop[4][3] ) +
          SQR( gradop[5][3] ) + SQR( gradop[6][3] ) + SQR( gradop[7][3] ) + SQR( gradop[8][3] ) );

    return (double)sqrt( aspect );
}
C_FUNC_DEF double v_hex_distortion ( int  num_nodes,
double  coordinates[][3] 
)

Calculates hex distortion metric.

distortion of a hex

Definition at line 2090 of file V_HexMetric.cpp.

References GaussIntegration::calculate_derivative_at_nodes_3d(), GaussIntegration::calculate_shape_function_3d_hex(), GaussIntegration::get_shape_func(), GaussIntegration::initialize(), maxNumberNodes, maxTotalNumberGaussPoints, VerdictVector::set(), and VERDICT_DBL_MAX.

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

{

    // use 2x2 gauss points for linear hex and 3x3 for 2nd order hex
    int number_of_gauss_points = 0;
    if( num_nodes == 8 )
        // 2x2 quadrature rule
        number_of_gauss_points = 2;
    else if( num_nodes == 20 )
        // 3x3 quadrature rule
        number_of_gauss_points = 3;

    int number_dimension             = 3;
    int total_number_of_gauss_points = number_of_gauss_points * number_of_gauss_points * number_of_gauss_points;
    double distortion                = VERDICT_DBL_MAX;

    // maxTotalNumberGaussPoints =27, maxNumberNodes = 20
    // they are defined in GaussIntegration.hpp
    // This is used to make these arrays static.
    // I tried dynamically allocated arrays but the new and delete
    // was very expensive

    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
    GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dimension );
    GaussIntegration::calculate_shape_function_3d_hex();
    GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight );

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

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

        element_volume += weight[ife] * jacobian;
    }

    // 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( 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 * 8.;
    return (double)distortion;
}
C_FUNC_DEF double v_hex_edge_ratio ( int  ,
double  coordinates[][3] 
)

Calculates hex edge ratio metric.

the edge ratio of a hex

NB (P. Pebay 01/23/07): Hmax / Hmin where Hmax and Hmin are respectively the maximum and the minimum edge lengths

Definition at line 520 of file V_HexMetric.cpp.

References b2, edges, VerdictVector::length_squared(), make_hex_edges(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{

    VerdictVector edges[12];
    make_hex_edges( coordinates, edges );

    double a2 = edges[0].length_squared();
    double b2 = edges[1].length_squared();
    double c2 = edges[2].length_squared();
    double d2 = edges[3].length_squared();
    double e2 = edges[4].length_squared();
    double f2 = edges[5].length_squared();
    double g2 = edges[6].length_squared();
    double h2 = edges[7].length_squared();
    double i2 = edges[8].length_squared();
    double j2 = edges[9].length_squared();
    double k2 = edges[10].length_squared();
    double l2 = edges[11].length_squared();

    double mab, mcd, mef, Mab, Mcd, Mef;
    double mgh, mij, mkl, Mgh, Mij, Mkl;

    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;
    }
    if( g2 < h2 )
    {
        mgh = g2;
        Mgh = h2;
    }
    else  // h2 <= g2
    {
        mgh = h2;
        Mgh = g2;
    }
    if( i2 < j2 )
    {
        mij = i2;
        Mij = j2;
    }
    else  // j2 <= i2
    {
        mij = j2;
        Mij = i2;
    }
    if( k2 < l2 )
    {
        mkl = k2;
        Mkl = l2;
    }
    else  // l2 <= k2
    {
        mkl = l2;
        Mkl = k2;
    }

    double m2;
    m2 = mab < mcd ? mab : mcd;
    m2 = m2 < mef ? m2 : mef;
    m2 = m2 < mgh ? m2 : mgh;
    m2 = m2 < mij ? m2 : mij;
    m2 = m2 < mkl ? m2 : mkl;

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

    double M2;
    M2 = Mab > Mcd ? Mab : Mcd;
    M2 = M2 > Mef ? M2 : Mef;
    M2 = M2 > Mgh ? M2 : Mgh;
    M2 = M2 > Mij ? M2 : Mij;
    M2 = M2 > Mkl ? M2 : Mkl;
    m2 = m2 < mef ? m2 : mef;

    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 );
}
int v_hex_get_weight ( VerdictVector v1,
VerdictVector v2,
VerdictVector v3 
)

weights based on the average size of a hex

Definition at line 39 of file V_HexMetric.cpp.

References scale(), VerdictVector::set(), and verdict_hex_size.

Referenced by v_hex_quality(), and v_hex_relative_size_squared().

{

    if( verdict_hex_size == 0 ) return 0;

    v1.set( 1, 0, 0 );
    v2.set( 0, 1, 0 );
    v3.set( 0, 0, 1 );

    double scale = pow( verdict_hex_size / ( v1 % ( v2 * v3 ) ), 0.33333333333 );
    v1 *= scale;
    v2 *= scale;
    v3 *= scale;

    return 1;
}
C_FUNC_DEF double v_hex_jacobian ( int  ,
double  coordinates[][3] 
)

Calculates hex jacobian metric.

jacobian of a hex

Minimum pointwise volume of local map at 8 corners & center of hex

Definition at line 1319 of file V_HexMetric.cpp.

References calc_hex_efg(), make_hex_nodes, VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{

    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    double jacobian = VERDICT_DBL_MAX;
    double current_jacobian;
    VerdictVector xxi, xet, xze;

    xxi = calc_hex_efg( 1, node_pos );
    xet = calc_hex_efg( 2, node_pos );
    xze = calc_hex_efg( 3, node_pos );

    current_jacobian = xxi % ( xet * xze ) / 64.0;
    if( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(0,0,0):

    xxi = node_pos[1] - node_pos[0];
    xet = node_pos[3] - node_pos[0];
    xze = node_pos[4] - node_pos[0];

    current_jacobian = xxi % ( xet * xze );
    if( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(1,0,0):

    xxi = node_pos[2] - node_pos[1];
    xet = node_pos[0] - node_pos[1];
    xze = node_pos[5] - node_pos[1];

    current_jacobian = xxi % ( xet * xze );
    if( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(1,1,0):

    xxi = node_pos[3] - node_pos[2];
    xet = node_pos[1] - node_pos[2];
    xze = node_pos[6] - node_pos[2];

    current_jacobian = xxi % ( xet * xze );
    if( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(0,1,0):

    xxi = node_pos[0] - node_pos[3];
    xet = node_pos[2] - node_pos[3];
    xze = node_pos[7] - node_pos[3];

    current_jacobian = xxi % ( xet * xze );
    if( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(0,0,1):

    xxi = node_pos[7] - node_pos[4];
    xet = node_pos[5] - node_pos[4];
    xze = node_pos[0] - node_pos[4];

    current_jacobian = xxi % ( xet * xze );
    if( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(1,0,1):

    xxi = node_pos[4] - node_pos[5];
    xet = node_pos[6] - node_pos[5];
    xze = node_pos[1] - node_pos[5];

    current_jacobian = xxi % ( xet * xze );
    if( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(1,1,1):

    xxi = node_pos[5] - node_pos[6];
    xet = node_pos[7] - node_pos[6];
    xze = node_pos[2] - node_pos[6];

    current_jacobian = xxi % ( xet * xze );
    if( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    // J(0,1,1):

    xxi = node_pos[6] - node_pos[7];
    xet = node_pos[4] - node_pos[7];
    xze = node_pos[3] - node_pos[7];

    current_jacobian = xxi % ( xet * xze );
    if( current_jacobian < jacobian ) { jacobian = current_jacobian; }

    if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_max_aspect_frobenius ( int  ,
double  coordinates[][3] 
)

Calculates hex condition metric.

maximum Frobenius condition number of a hex

Maximum Frobenius condition number of the Jacobian matrix at 8 corners NB (P. Pebay 01/25/07): this metric is calculated by taking the maximum of the 8 Frobenius aspects at each corner of the hex, when the reference corner is right isosceles.

Definition at line 1207 of file V_HexMetric.cpp.

References calc_hex_efg(), condition_comp(), make_hex_nodes, VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{

    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    double condition = 0.0, current_condition;
    VerdictVector xxi, xet, xze;

    xxi = calc_hex_efg( 1, node_pos );
    xet = calc_hex_efg( 2, node_pos );
    xze = calc_hex_efg( 3, node_pos );

    current_condition = condition_comp( xxi, xet, xze );
    if( current_condition > condition ) { condition = current_condition; }

    // J(0,0,0):

    xxi = node_pos[1] - node_pos[0];
    xet = node_pos[3] - node_pos[0];
    xze = node_pos[4] - node_pos[0];

    current_condition = condition_comp( xxi, xet, xze );
    if( current_condition > condition ) { condition = current_condition; }

    // J(1,0,0):

    xxi = node_pos[2] - node_pos[1];
    xet = node_pos[0] - node_pos[1];
    xze = node_pos[5] - node_pos[1];

    current_condition = condition_comp( xxi, xet, xze );
    if( current_condition > condition ) { condition = current_condition; }

    // J(1,1,0):

    xxi = node_pos[3] - node_pos[2];
    xet = node_pos[1] - node_pos[2];
    xze = node_pos[6] - node_pos[2];

    current_condition = condition_comp( xxi, xet, xze );
    if( current_condition > condition ) { condition = current_condition; }

    // J(0,1,0):

    xxi = node_pos[0] - node_pos[3];
    xet = node_pos[2] - node_pos[3];
    xze = node_pos[7] - node_pos[3];

    current_condition = condition_comp( xxi, xet, xze );
    if( current_condition > condition ) { condition = current_condition; }

    // J(0,0,1):

    xxi = node_pos[7] - node_pos[4];
    xet = node_pos[5] - node_pos[4];
    xze = node_pos[0] - node_pos[4];

    current_condition = condition_comp( xxi, xet, xze );
    if( current_condition > condition ) { condition = current_condition; }

    // J(1,0,1):

    xxi = node_pos[4] - node_pos[5];
    xet = node_pos[6] - node_pos[5];
    xze = node_pos[1] - node_pos[5];

    current_condition = condition_comp( xxi, xet, xze );
    if( current_condition > condition ) { condition = current_condition; }

    // J(1,1,1):

    xxi = node_pos[5] - node_pos[6];
    xet = node_pos[7] - node_pos[6];
    xze = node_pos[2] - node_pos[6];

    current_condition = condition_comp( xxi, xet, xze );
    if( current_condition > condition ) { condition = current_condition; }

    // J(1,1,1):

    xxi = node_pos[6] - node_pos[7];
    xet = node_pos[4] - node_pos[7];
    xze = node_pos[3] - node_pos[7];

    current_condition = condition_comp( xxi, xet, xze );
    if( current_condition > condition ) { condition = current_condition; }

    condition /= 3.0;

    if( condition > 0 ) return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( condition, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_max_edge_ratio ( int  ,
double  coordinates[][3] 
)

Calculates hex maximum of edge ratio.

max edge ratio of a hex

Maximum edge length ratio at hex center

Definition at line 634 of file V_HexMetric.cpp.

References calc_hex_efg(), VerdictVector::length(), make_hex_nodes, safe_ratio(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    double aspect;
    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    double aspect_12, aspect_13, aspect_23;

    VerdictVector efg1 = calc_hex_efg( 1, node_pos );
    VerdictVector efg2 = calc_hex_efg( 2, node_pos );
    VerdictVector efg3 = calc_hex_efg( 3, node_pos );

    double mag_efg1 = efg1.length();
    double mag_efg2 = efg2.length();
    double mag_efg3 = efg3.length();

    aspect_12 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg2 ), VERDICT_MIN( mag_efg1, mag_efg2 ) );
    aspect_13 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg3 ), VERDICT_MIN( mag_efg1, mag_efg3 ) );
    aspect_23 = safe_ratio( VERDICT_MAX( mag_efg2, mag_efg3 ), VERDICT_MIN( mag_efg2, mag_efg3 ) );

    aspect = VERDICT_MAX( aspect_12, VERDICT_MAX( aspect_13, aspect_23 ) );

    if( aspect > 0 ) return (double)VERDICT_MIN( aspect, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( aspect, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_med_aspect_frobenius ( int  ,
double  coordinates[][3] 
)

Calculates hex condition metric.

the average Frobenius aspect of a hex

NB (P. Pebay 01/20/07): this metric is calculated by averaging the 8 Frobenius aspects at each corner of the hex, when the reference corner is right isosceles.

Definition at line 1128 of file V_HexMetric.cpp.

References condition_comp(), make_hex_nodes, VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{

    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    double med_aspect_frobenius = 0.;
    VerdictVector xxi, xet, xze;

    // J(0,0,0):

    xxi = node_pos[1] - node_pos[0];
    xet = node_pos[3] - node_pos[0];
    xze = node_pos[4] - node_pos[0];

    med_aspect_frobenius += condition_comp( xxi, xet, xze );
    // J(1,0,0):

    xxi = node_pos[2] - node_pos[1];
    xet = node_pos[0] - node_pos[1];
    xze = node_pos[5] - node_pos[1];

    med_aspect_frobenius += condition_comp( xxi, xet, xze );
    // J(1,1,0):

    xxi = node_pos[3] - node_pos[2];
    xet = node_pos[1] - node_pos[2];
    xze = node_pos[6] - node_pos[2];

    med_aspect_frobenius += condition_comp( xxi, xet, xze );
    // J(0,1,0):

    xxi = node_pos[0] - node_pos[3];
    xet = node_pos[2] - node_pos[3];
    xze = node_pos[7] - node_pos[3];

    med_aspect_frobenius += condition_comp( xxi, xet, xze );
    // J(0,0,1):

    xxi = node_pos[7] - node_pos[4];
    xet = node_pos[5] - node_pos[4];
    xze = node_pos[0] - node_pos[4];

    med_aspect_frobenius += condition_comp( xxi, xet, xze );
    // J(1,0,1):

    xxi = node_pos[4] - node_pos[5];
    xet = node_pos[6] - node_pos[5];
    xze = node_pos[1] - node_pos[5];

    med_aspect_frobenius += condition_comp( xxi, xet, xze );
    // J(1,1,1):

    xxi = node_pos[5] - node_pos[6];
    xet = node_pos[7] - node_pos[6];
    xze = node_pos[2] - node_pos[6];

    med_aspect_frobenius += condition_comp( xxi, xet, xze );
    // J(1,1,1):

    xxi = node_pos[6] - node_pos[7];
    xet = node_pos[4] - node_pos[7];
    xze = node_pos[3] - node_pos[7];

    med_aspect_frobenius += condition_comp( xxi, xet, xze );
    med_aspect_frobenius /= 24.;

    if( med_aspect_frobenius > 0 ) return (double)VERDICT_MIN( med_aspect_frobenius, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( med_aspect_frobenius, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_oddy ( int  ,
double  coordinates[][3] 
)

Calculates hex oddy metric.

oddy of a hex

General distortion measure based on left Cauchy-Green Tensor

Definition at line 1005 of file V_HexMetric.cpp.

References calc_hex_efg(), make_hex_nodes, oddy_comp(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{

    double oddy = 0.0, current_oddy;
    VerdictVector xxi, xet, xze;

    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    xxi = calc_hex_efg( 1, node_pos );
    xet = calc_hex_efg( 2, node_pos );
    xze = calc_hex_efg( 3, node_pos );

    current_oddy = oddy_comp( xxi, xet, xze );
    if( current_oddy > oddy ) { oddy = current_oddy; }

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

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

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

    current_oddy = oddy_comp( xxi, xet, xze );
    if( current_oddy > oddy ) { oddy = current_oddy; }

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

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

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

    current_oddy = oddy_comp( xxi, xet, xze );
    if( current_oddy > oddy ) { oddy = current_oddy; }

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

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

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

    current_oddy = oddy_comp( xxi, xet, xze );
    if( current_oddy > oddy ) { oddy = current_oddy; }

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

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

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

    current_oddy = oddy_comp( xxi, xet, xze );
    if( current_oddy > oddy ) { oddy = current_oddy; }

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

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

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

    current_oddy = oddy_comp( xxi, xet, xze );
    if( current_oddy > oddy ) { oddy = current_oddy; }

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

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

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

    current_oddy = oddy_comp( xxi, xet, xze );
    if( current_oddy > oddy ) { oddy = current_oddy; }

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

    xet.set( coordinates[7][0] - coordinates[6][0], coordinates[7][1] - coordinates[6][1],
             coordinates[7][2] - coordinates[6][2] );

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

    current_oddy = oddy_comp( xxi, xet, xze );
    if( current_oddy > oddy ) { oddy = current_oddy; }

    xxi.set( coordinates[6][0] - coordinates[7][0], coordinates[6][1] - coordinates[7][1],
             coordinates[6][2] - coordinates[7][2] );

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

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

    current_oddy = oddy_comp( xxi, xet, xze );
    if( current_oddy > oddy ) { oddy = current_oddy; }

    if( oddy > 0 ) return (double)VERDICT_MIN( oddy, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( oddy, -VERDICT_DBL_MAX );
}
C_FUNC_DEF void v_hex_quality ( int  num_nodes,
double  coordinates[][3],
unsigned int  metrics_request_flag,
HexMetricVals metric_vals 
)

Calculates quality metrics for hexahedral elements.

multiple quality metrics of a hex

Definition at line 2513 of file V_HexMetric.cpp.

References calc_hex_efg(), HexMetricVals::condition, condition_comp(), diag_length(), HexMetricVals::diagonal, HexMetricVals::dimension, HexMetricVals::distortion, HexMetricVals::edge_ratio, edges, HexMetricVals::jacobian, VerdictVector::length(), length_squared(), VerdictVector::length_squared(), make_edge_length_squares, make_hex_edges(), make_hex_nodes, HexMetricVals::max_edge_ratio, HexMetricVals::med_aspect_frobenius, VerdictVector::normalize(), HexMetricVals::oddy, oddy_comp(), HexMetricVals::relative_size_squared, safe_ratio(), HexMetricVals::scaled_jacobian, HexMetricVals::shape, HexMetricVals::shape_and_size, HexMetricVals::shear, HexMetricVals::shear_and_size, HexMetricVals::skew, HexMetricVals::stretch, HexMetricVals::taper, V_HEX_CONDITION, V_HEX_DIAGONAL, v_hex_diagonal(), V_HEX_DIMENSION, v_hex_dimension(), V_HEX_DISTORTION, v_hex_distortion(), V_HEX_EDGE_RATIO, v_hex_edge_ratio(), v_hex_get_weight(), V_HEX_JACOBIAN, V_HEX_MAX_EDGE_RATIO, V_HEX_MED_ASPECT_FROBENIUS, v_hex_med_aspect_frobenius(), V_HEX_ODDY, V_HEX_RELATIVE_SIZE_SQUARED, V_HEX_SCALED_JACOBIAN, V_HEX_SHAPE, V_HEX_SHAPE_AND_SIZE, V_HEX_SHEAR, V_HEX_SHEAR_AND_SIZE, V_HEX_SKEW, V_HEX_STRETCH, V_HEX_TAPER, V_HEX_VOLUME, v_hex_volume(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_FALSE, VERDICT_MAX, VERDICT_MIN, VERDICT_TRUE, and HexMetricVals::volume.

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

{
    memset( metric_vals, 0, sizeof( HexMetricVals ) );

    // max edge ratio, skew, taper, and volume
    if( metrics_request_flag & ( V_HEX_MAX_EDGE_RATIO | V_HEX_SKEW | V_HEX_TAPER ) )
    {
        VerdictVector node_pos[8];
        make_hex_nodes( coordinates, node_pos );

        VerdictVector efg1, efg2, efg3;
        efg1 = calc_hex_efg( 1, node_pos );
        efg2 = calc_hex_efg( 2, node_pos );
        efg3 = calc_hex_efg( 3, node_pos );

        if( metrics_request_flag & V_HEX_MAX_EDGE_RATIO )
        {
            double max_edge_ratio_12, max_edge_ratio_13, max_edge_ratio_23;

            double mag_efg1 = efg1.length();
            double mag_efg2 = efg2.length();
            double mag_efg3 = efg3.length();

            max_edge_ratio_12 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg2 ), VERDICT_MIN( mag_efg1, mag_efg2 ) );
            max_edge_ratio_13 = safe_ratio( VERDICT_MAX( mag_efg1, mag_efg3 ), VERDICT_MIN( mag_efg1, mag_efg3 ) );
            max_edge_ratio_23 = safe_ratio( VERDICT_MAX( mag_efg2, mag_efg3 ), VERDICT_MIN( mag_efg2, mag_efg3 ) );

            metric_vals->max_edge_ratio =
                (double)VERDICT_MAX( max_edge_ratio_12, VERDICT_MAX( max_edge_ratio_13, max_edge_ratio_23 ) );
        }

        if( metrics_request_flag & V_HEX_SKEW )
        {

            VerdictVector vec1 = efg1;
            VerdictVector vec2 = efg2;
            VerdictVector vec3 = efg3;

            if( vec1.normalize() <= VERDICT_DBL_MIN || vec2.normalize() <= VERDICT_DBL_MIN ||
                vec3.normalize() <= VERDICT_DBL_MIN )
            {
                metric_vals->skew = (double)VERDICT_DBL_MAX;
            }
            else
            {
                double skewx = fabs( vec1 % vec2 );
                double skewy = fabs( vec1 % vec3 );
                double skewz = fabs( vec2 % vec3 );

                metric_vals->skew = (double)( VERDICT_MAX( skewx, VERDICT_MAX( skewy, skewz ) ) );
            }
        }

        if( metrics_request_flag & V_HEX_TAPER )
        {
            VerdictVector efg12 = calc_hex_efg( 12, node_pos );
            VerdictVector efg13 = calc_hex_efg( 13, node_pos );
            VerdictVector efg23 = calc_hex_efg( 23, node_pos );

            double taperx = fabs( safe_ratio( efg12.length(), VERDICT_MIN( efg1.length(), efg2.length() ) ) );
            double tapery = fabs( safe_ratio( efg13.length(), VERDICT_MIN( efg1.length(), efg3.length() ) ) );
            double taperz = fabs( safe_ratio( efg23.length(), VERDICT_MIN( efg2.length(), efg3.length() ) ) );

            metric_vals->taper = (double)VERDICT_MAX( taperx, VERDICT_MAX( tapery, taperz ) );
        }
    }

    if( metrics_request_flag & V_HEX_VOLUME ) { metric_vals->volume = v_hex_volume( 8, coordinates ); }

    if( metrics_request_flag &
        ( V_HEX_JACOBIAN | V_HEX_SCALED_JACOBIAN | V_HEX_CONDITION | V_HEX_SHEAR | V_HEX_SHAPE |
          V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) )
    {

        static const double two_thirds = 2.0 / 3.0;
        VerdictVector edges[12];
        // the length squares
        double length_squared[12];
        // make vectors from coordinates
        make_hex_edges( coordinates, edges );

        // calculate the length squares if we need to
        if( metrics_request_flag &
            ( V_HEX_JACOBIAN | V_HEX_SHEAR | V_HEX_SCALED_JACOBIAN | V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE |
              V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHEAR_AND_SIZE | V_HEX_STRETCH ) )
            make_edge_length_squares( edges, length_squared );

        double jacobian = VERDICT_DBL_MAX, scaled_jacobian = VERDICT_DBL_MAX, condition = 0.0, shear = 1.0, shape = 1.0,
               oddy = 0.0;
        double current_jacobian, current_scaled_jacobian, current_condition, current_shape, detw = 0, det_sum = 0,
                                                                                            current_oddy;
        VerdictBoolean rel_size_error = VERDICT_FALSE;

        VerdictVector xxi, xet, xze;

        // get weights if we need based on average size of a hex
        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            v_hex_get_weight( xxi, xet, xze );
            detw = xxi % ( xet * xze );
            if( detw < VERDICT_DBL_MIN ) rel_size_error = VERDICT_TRUE;
        }

        xxi = edges[0] - edges[2] + edges[4] - edges[6];
        xet = edges[1] - edges[3] + edges[5] - edges[7];
        xze = edges[8] + edges[9] + edges[10] + edges[11];

        current_jacobian = xxi % ( xet * xze ) / 64.0;
        if( current_jacobian < jacobian ) jacobian = current_jacobian;

        if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
        {
            current_jacobian *= 64.0;
            current_scaled_jacobian =
                current_jacobian / sqrt( xxi.length_squared() * xet.length_squared() * xze.length_squared() );
            if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
        }

        if( metrics_request_flag & V_HEX_CONDITION )
        {
            current_condition = condition_comp( xxi, xet, xze );
            if( current_condition > condition ) { condition = current_condition; }
        }

        if( metrics_request_flag & V_HEX_ODDY )
        {
            current_oddy = oddy_comp( xxi, xet, xze );
            if( current_oddy > oddy ) { oddy = current_oddy; }
        }

        // J(0,0,0)
        current_jacobian = edges[0] % ( -edges[3] * edges[8] );
        if( current_jacobian < jacobian ) jacobian = current_jacobian;

        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            det_sum += current_jacobian;
        }

        if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
        {
            if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN ||
                length_squared[8] <= VERDICT_DBL_MIN )
            {
                current_scaled_jacobian = VERDICT_DBL_MAX;
            }
            else
            {
                current_scaled_jacobian =
                    current_jacobian / sqrt( length_squared[0] * length_squared[3] * length_squared[8] );
            }
            if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
        }

        if( metrics_request_flag & V_HEX_CONDITION )
        {
            current_condition = condition_comp( edges[0], -edges[3], edges[8] );
            if( current_condition > condition ) { condition = current_condition; }
        }

        if( metrics_request_flag & V_HEX_ODDY )
        {
            current_oddy = oddy_comp( edges[0], -edges[3], edges[8] );
            if( current_oddy > oddy ) { oddy = current_oddy; }
        }

        if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
        {
            if( current_jacobian > VERDICT_DBL_MIN )
                current_shape = 3 * pow( current_jacobian, two_thirds ) /
                                ( length_squared[0] + length_squared[3] + length_squared[8] );
            else
                current_shape = 0;

            if( current_shape < shape ) { shape = current_shape; }
        }

        // J(1,0,0)
        current_jacobian = edges[1] % ( -edges[0] * edges[9] );
        if( current_jacobian < jacobian ) jacobian = current_jacobian;

        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            det_sum += current_jacobian;
        }

        if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
        {
            if( length_squared[1] <= VERDICT_DBL_MIN || length_squared[0] <= VERDICT_DBL_MIN ||
                length_squared[9] <= VERDICT_DBL_MIN )
            {
                current_scaled_jacobian = VERDICT_DBL_MAX;
            }
            else
            {
                current_scaled_jacobian =
                    current_jacobian / sqrt( length_squared[1] * length_squared[0] * length_squared[9] );
            }
            if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
        }

        if( metrics_request_flag & V_HEX_CONDITION )
        {
            current_condition = condition_comp( edges[1], -edges[0], edges[9] );
            if( current_condition > condition ) { condition = current_condition; }
        }

        if( metrics_request_flag & V_HEX_ODDY )
        {
            current_oddy = oddy_comp( edges[1], -edges[0], edges[9] );
            if( current_oddy > oddy ) { oddy = current_oddy; }
        }

        if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
        {
            if( current_jacobian > VERDICT_DBL_MIN )
                current_shape = 3 * pow( current_jacobian, two_thirds ) /
                                ( length_squared[1] + length_squared[0] + length_squared[9] );
            else
                current_shape = 0;

            if( current_shape < shape ) { shape = current_shape; }
        }

        // J(1,1,0)
        current_jacobian = edges[2] % ( -edges[1] * edges[10] );
        if( current_jacobian < jacobian ) jacobian = current_jacobian;

        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            det_sum += current_jacobian;
        }

        if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
        {
            if( length_squared[2] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN ||
                length_squared[10] <= VERDICT_DBL_MIN )
            {
                current_scaled_jacobian = VERDICT_DBL_MAX;
            }
            else
            {
                current_scaled_jacobian =
                    current_jacobian / sqrt( length_squared[2] * length_squared[1] * length_squared[10] );
            }
            if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
        }

        if( metrics_request_flag & V_HEX_CONDITION )
        {
            current_condition = condition_comp( edges[2], -edges[1], edges[10] );
            if( current_condition > condition ) { condition = current_condition; }
        }

        if( metrics_request_flag & V_HEX_ODDY )
        {
            current_oddy = oddy_comp( edges[2], -edges[1], edges[10] );
            if( current_oddy > oddy ) { oddy = current_oddy; }
        }

        if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
        {
            if( current_jacobian > VERDICT_DBL_MIN )
                current_shape = 3 * pow( current_jacobian, two_thirds ) /
                                ( length_squared[2] + length_squared[1] + length_squared[10] );
            else
                current_shape = 0;

            if( current_shape < shape ) { shape = current_shape; }
        }

        // J(0,1,0)
        current_jacobian = edges[3] % ( -edges[2] * edges[11] );
        if( current_jacobian < jacobian ) jacobian = current_jacobian;

        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            det_sum += current_jacobian;
        }

        if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
        {
            if( length_squared[3] <= VERDICT_DBL_MIN || length_squared[2] <= VERDICT_DBL_MIN ||
                length_squared[11] <= VERDICT_DBL_MIN )
            {
                current_scaled_jacobian = VERDICT_DBL_MAX;
            }
            else
            {
                current_scaled_jacobian =
                    current_jacobian / sqrt( length_squared[3] * length_squared[2] * length_squared[11] );
            }
            if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
        }

        if( metrics_request_flag & V_HEX_CONDITION )
        {
            current_condition = condition_comp( edges[3], -edges[2], edges[11] );
            if( current_condition > condition ) { condition = current_condition; }
        }

        if( metrics_request_flag & V_HEX_ODDY )
        {
            current_oddy = oddy_comp( edges[3], -edges[2], edges[11] );
            if( current_oddy > oddy ) { oddy = current_oddy; }
        }

        if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
        {
            if( current_jacobian > VERDICT_DBL_MIN )
                current_shape = 3 * pow( current_jacobian, two_thirds ) /
                                ( length_squared[3] + length_squared[2] + length_squared[11] );
            else
                current_shape = 0;

            if( current_shape < shape ) { shape = current_shape; }
        }

        // J(0,0,1)
        current_jacobian = edges[4] % ( -edges[8] * -edges[7] );
        if( current_jacobian < jacobian ) jacobian = current_jacobian;

        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            det_sum += current_jacobian;
        }

        if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
        {
            if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[8] <= VERDICT_DBL_MIN ||
                length_squared[7] <= VERDICT_DBL_MIN )
            {
                current_scaled_jacobian = VERDICT_DBL_MAX;
            }
            else
            {
                current_scaled_jacobian =
                    current_jacobian / sqrt( length_squared[4] * length_squared[8] * length_squared[7] );
            }
            if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
        }

        if( metrics_request_flag & V_HEX_CONDITION )
        {
            current_condition = condition_comp( edges[4], -edges[8], -edges[7] );
            if( current_condition > condition ) { condition = current_condition; }
        }

        if( metrics_request_flag & V_HEX_ODDY )
        {
            current_oddy = oddy_comp( edges[4], -edges[8], -edges[7] );
            if( current_oddy > oddy ) { oddy = current_oddy; }
        }

        if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
        {
            if( current_jacobian > VERDICT_DBL_MIN )
                current_shape = 3 * pow( current_jacobian, two_thirds ) /
                                ( length_squared[4] + length_squared[8] + length_squared[7] );
            else
                current_shape = 0;

            if( current_shape < shape ) { shape = current_shape; }
        }

        // J(1,0,1)
        current_jacobian = -edges[4] % ( edges[5] * -edges[9] );
        if( current_jacobian < jacobian ) jacobian = current_jacobian;

        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            det_sum += current_jacobian;
        }

        if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
        {
            if( length_squared[4] <= VERDICT_DBL_MIN || length_squared[5] <= VERDICT_DBL_MIN ||
                length_squared[9] <= VERDICT_DBL_MIN )
            {
                current_scaled_jacobian = VERDICT_DBL_MAX;
            }
            else
            {
                current_scaled_jacobian =
                    current_jacobian / sqrt( length_squared[4] * length_squared[5] * length_squared[9] );
            }
            if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
        }

        if( metrics_request_flag & V_HEX_CONDITION )
        {
            current_condition = condition_comp( -edges[4], edges[5], -edges[9] );
            if( current_condition > condition ) { condition = current_condition; }
        }

        if( metrics_request_flag & V_HEX_ODDY )
        {
            current_oddy = oddy_comp( -edges[4], edges[5], -edges[9] );
            if( current_oddy > oddy ) { oddy = current_oddy; }
        }

        if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
        {
            if( current_jacobian > VERDICT_DBL_MIN )
                current_shape = 3 * pow( current_jacobian, two_thirds ) /
                                ( length_squared[4] + length_squared[5] + length_squared[9] );
            else
                current_shape = 0;

            if( current_shape < shape ) { shape = current_shape; }
        }

        // J(1,1,1)
        current_jacobian = -edges[5] % ( edges[6] * -edges[10] );
        if( current_jacobian < jacobian ) jacobian = current_jacobian;

        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            det_sum += current_jacobian;
        }

        if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
        {
            if( length_squared[5] <= VERDICT_DBL_MIN || length_squared[6] <= VERDICT_DBL_MIN ||
                length_squared[10] <= VERDICT_DBL_MIN )
            {
                current_scaled_jacobian = VERDICT_DBL_MAX;
            }
            else
            {
                current_scaled_jacobian =
                    current_jacobian / sqrt( length_squared[5] * length_squared[6] * length_squared[10] );
            }
            if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
        }

        if( metrics_request_flag & V_HEX_CONDITION )
        {
            current_condition = condition_comp( -edges[5], edges[6], -edges[10] );
            if( current_condition > condition ) { condition = current_condition; }
        }

        if( metrics_request_flag & V_HEX_ODDY )
        {
            current_oddy = oddy_comp( -edges[5], edges[6], -edges[10] );
            if( current_oddy > oddy ) { oddy = current_oddy; }
        }

        if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
        {
            if( current_jacobian > VERDICT_DBL_MIN )
                current_shape = 3 * pow( current_jacobian, two_thirds ) /
                                ( length_squared[5] + length_squared[6] + length_squared[10] );
            else
                current_shape = 0;

            if( current_shape < shape ) { shape = current_shape; }
        }

        // J(0,1,1)
        current_jacobian = -edges[6] % ( edges[7] * -edges[11] );
        if( current_jacobian < jacobian ) jacobian = current_jacobian;

        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            det_sum += current_jacobian;
        }

        if( metrics_request_flag & ( V_HEX_SCALED_JACOBIAN | V_HEX_SHEAR | V_HEX_SHEAR_AND_SIZE ) )
        {
            if( length_squared[6] <= VERDICT_DBL_MIN || length_squared[7] <= VERDICT_DBL_MIN ||
                length_squared[11] <= VERDICT_DBL_MIN )
            {
                current_scaled_jacobian = VERDICT_DBL_MAX;
            }
            else
            {
                current_scaled_jacobian =
                    current_jacobian / sqrt( length_squared[6] * length_squared[7] * length_squared[11] );
            }
            if( current_scaled_jacobian < scaled_jacobian ) shear = scaled_jacobian = current_scaled_jacobian;
        }

        if( metrics_request_flag & V_HEX_CONDITION )
        {
            current_condition = condition_comp( -edges[6], edges[7], -edges[11] );
            if( current_condition > condition ) { condition = current_condition; }
        }

        if( metrics_request_flag & V_HEX_ODDY )
        {
            current_oddy = oddy_comp( -edges[6], edges[7], -edges[11] );
            if( current_oddy > oddy ) { oddy = current_oddy; }
        }

        if( metrics_request_flag & ( V_HEX_SHAPE | V_HEX_SHAPE_AND_SIZE ) )
        {
            if( current_jacobian > VERDICT_DBL_MIN )
                current_shape = 3 * pow( current_jacobian, two_thirds ) /
                                ( length_squared[6] + length_squared[7] + length_squared[11] );
            else
                current_shape = 0;

            if( current_shape < shape ) { shape = current_shape; }
        }

        if( metrics_request_flag & ( V_HEX_RELATIVE_SIZE_SQUARED | V_HEX_SHAPE_AND_SIZE | V_HEX_SHEAR_AND_SIZE ) )
        {
            if( det_sum > VERDICT_DBL_MIN && rel_size_error != VERDICT_TRUE )
            {
                double tau                         = det_sum / ( 8 * detw );
                metric_vals->relative_size_squared = (double)VERDICT_MIN( tau * tau, 1.0 / tau / tau );
            }
            else
                metric_vals->relative_size_squared = 0.0;
        }

        // set values from above calculations
        if( metrics_request_flag & V_HEX_JACOBIAN ) metric_vals->jacobian = (double)jacobian;

        if( metrics_request_flag & V_HEX_SCALED_JACOBIAN ) metric_vals->scaled_jacobian = (double)scaled_jacobian;

        if( metrics_request_flag & V_HEX_CONDITION ) metric_vals->condition = (double)( condition / 3.0 );

        if( metrics_request_flag & V_HEX_SHEAR )
        {
            if( shear < VERDICT_DBL_MIN )  // shear has range 0 to +1
                shear = 0;
            metric_vals->shear = (double)shear;
        }

        if( metrics_request_flag & V_HEX_SHAPE ) metric_vals->shape = (double)shape;

        if( metrics_request_flag & V_HEX_SHAPE_AND_SIZE )
            metric_vals->shape_and_size = (double)( shape * metric_vals->relative_size_squared );

        if( metrics_request_flag & V_HEX_SHEAR_AND_SIZE )
            metric_vals->shear_and_size = (double)( shear * metric_vals->relative_size_squared );

        if( metrics_request_flag & V_HEX_ODDY ) metric_vals->oddy = (double)oddy;

        if( metrics_request_flag & V_HEX_STRETCH )
        {
            static const double HEX_STRETCH_SCALE_FACTOR = sqrt( 3.0 );
            double min_edge                              = length_squared[0];
            for( int j = 1; j < 12; j++ )
                min_edge = VERDICT_MIN( min_edge, length_squared[j] );

            double max_diag = diag_length( 1, coordinates );

            metric_vals->stretch = (double)( HEX_STRETCH_SCALE_FACTOR * ( safe_ratio( sqrt( min_edge ), max_diag ) ) );
        }
    }

    if( metrics_request_flag & V_HEX_DIAGONAL ) metric_vals->diagonal = v_hex_diagonal( num_nodes, coordinates );
    if( metrics_request_flag & V_HEX_DIMENSION ) metric_vals->dimension = v_hex_dimension( num_nodes, coordinates );
    if( metrics_request_flag & V_HEX_DISTORTION ) metric_vals->distortion = v_hex_distortion( num_nodes, coordinates );

    // take care of any overflow problems
    // max_edge_ratio
    if( metric_vals->max_edge_ratio > 0 )
        metric_vals->max_edge_ratio = (double)VERDICT_MIN( metric_vals->max_edge_ratio, VERDICT_DBL_MAX );
    else
        metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX );

    // skew
    if( metric_vals->skew > 0 )
        metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX );
    else
        metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX );

    // taper
    if( metric_vals->taper > 0 )
        metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX );
    else
        metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX );

    // volume
    if( metric_vals->volume > 0 )
        metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
    else
        metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );

    // stretch
    if( metric_vals->stretch > 0 )
        metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX );
    else
        metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX );

    // diagonal
    if( metric_vals->diagonal > 0 )
        metric_vals->diagonal = (double)VERDICT_MIN( metric_vals->diagonal, VERDICT_DBL_MAX );
    else
        metric_vals->diagonal = (double)VERDICT_MAX( metric_vals->diagonal, -VERDICT_DBL_MAX );

    // dimension
    if( metric_vals->dimension > 0 )
        metric_vals->dimension = (double)VERDICT_MIN( metric_vals->dimension, VERDICT_DBL_MAX );
    else
        metric_vals->dimension = (double)VERDICT_MAX( metric_vals->dimension, -VERDICT_DBL_MAX );

    // oddy
    if( metric_vals->oddy > 0 )
        metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX );
    else
        metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX );

    // condition
    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 );

    // jacobian
    if( metric_vals->jacobian > 0 )
        metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
    else
        metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );

    // scaled_jacobian
    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 );

    // shear
    if( metric_vals->shear > 0 )
        metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX );
    else
        metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX );

    // shape
    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 );

    // 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 );
    else
        metric_vals->relative_size_squared =
            (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );

    // 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 );
    else
        metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );

    // shear_and_size
    if( metric_vals->shear_and_size > 0 )
        metric_vals->shear_and_size = (double)VERDICT_MIN( metric_vals->shear_and_size, VERDICT_DBL_MAX );
    else
        metric_vals->shear_and_size = (double)VERDICT_MAX( metric_vals->shear_and_size, -VERDICT_DBL_MAX );

    // distortion
    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_HEX_MED_ASPECT_FROBENIUS )
        metric_vals->med_aspect_frobenius = v_hex_med_aspect_frobenius( 8, coordinates );

    if( metrics_request_flag & V_HEX_EDGE_RATIO ) metric_vals->edge_ratio = v_hex_edge_ratio( 8, coordinates );
}
C_FUNC_DEF double v_hex_relative_size_squared ( int  ,
double  coordinates[][3] 
)

Calculates hex relative size metric.

relative size of a hex

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

Definition at line 1952 of file V_HexMetric.cpp.

References MBMesquite::det(), make_hex_nodes, size, v_hex_get_weight(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure(), v_hex_shape_and_size(), and v_hex_shear_and_size().

{
    double size = 0;
    double tau;

    VerdictVector xxi, xet, xze;
    double det, det_sum = 0;

    v_hex_get_weight( xxi, xet, xze );

    // This is the average relative size
    double detw = xxi % ( xet * xze );

    if( detw < VERDICT_DBL_MIN ) return 0;

    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    // J(0,0,0):

    xxi = node_pos[1] - node_pos[0];
    xet = node_pos[3] - node_pos[0];
    xze = node_pos[4] - node_pos[0];

    det = xxi % ( xet * xze );
    det_sum += det;

    // J(1,0,0):

    xxi = node_pos[2] - node_pos[1];
    xet = node_pos[0] - node_pos[1];
    xze = node_pos[5] - node_pos[1];

    det = xxi % ( xet * xze );
    det_sum += det;

    // J(0,1,0):

    xxi = node_pos[3] - node_pos[2];
    xet = node_pos[1] - node_pos[2];
    xze = node_pos[6] - node_pos[2];

    det = xxi % ( xet * xze );
    det_sum += det;

    // J(1,1,0):

    xxi = node_pos[0] - node_pos[3];
    xet = node_pos[2] - node_pos[3];
    xze = node_pos[7] - node_pos[3];

    det = xxi % ( xet * xze );
    det_sum += det;

    // J(0,1,0):

    xxi = node_pos[7] - node_pos[4];
    xet = node_pos[5] - node_pos[4];
    xze = node_pos[0] - node_pos[4];

    det = xxi % ( xet * xze );
    det_sum += det;

    // J(1,0,1):

    xxi = node_pos[4] - node_pos[5];
    xet = node_pos[6] - node_pos[5];
    xze = node_pos[1] - node_pos[5];

    det = xxi % ( xet * xze );
    det_sum += det;

    // J(1,1,1):

    xxi = node_pos[5] - node_pos[6];
    xet = node_pos[7] - node_pos[6];
    xze = node_pos[2] - node_pos[6];

    det = xxi % ( xet * xze );
    det_sum += det;

    // J(1,1,1):

    xxi = node_pos[6] - node_pos[7];
    xet = node_pos[4] - node_pos[7];
    xze = node_pos[3] - node_pos[7];

    det = xxi % ( xet * xze );
    det_sum += det;

    if( det_sum > VERDICT_DBL_MIN )
    {
        tau = det_sum / ( 8 * detw );

        tau = VERDICT_MIN( tau, 1.0 / tau );

        size = tau * tau;
    }

    if( size > 0 ) return (double)VERDICT_MIN( size, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( size, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_scaled_jacobian ( int  ,
double  coordinates[][3] 
)

Calculates hex scaled jacobian metric.

scaled jacobian of a hex

Minimum Jacobian divided by the lengths of the 3 edge vectors

Definition at line 1417 of file V_HexMetric.cpp.

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

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

{

    double jacobi, min_norm_jac = VERDICT_DBL_MAX;
    // double min_jacobi = VERDICT_DBL_MAX;
    double temp_norm_jac, lengths;
    double len1_sq, len2_sq, len3_sq;
    VerdictVector xxi, xet, xze;

    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    xxi = calc_hex_efg( 1, node_pos );
    xet = calc_hex_efg( 2, node_pos );
    xze = calc_hex_efg( 3, node_pos );

    jacobi = xxi % ( xet * xze );
    // if( jacobi < min_jacobi) { min_jacobi = jacobi; }

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    lengths       = sqrt( len1_sq * len2_sq * len3_sq );
    temp_norm_jac = jacobi / lengths;

    if( temp_norm_jac < min_norm_jac )
        min_norm_jac = temp_norm_jac;
    else
        temp_norm_jac = jacobi;

    // J(0,0,0):

    xxi = node_pos[1] - node_pos[0];
    xet = node_pos[3] - node_pos[0];
    xze = node_pos[4] - node_pos[0];

    jacobi = xxi % ( xet * xze );
    // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    lengths       = sqrt( len1_sq * len2_sq * len3_sq );
    temp_norm_jac = jacobi / lengths;
    if( temp_norm_jac < min_norm_jac )
        min_norm_jac = temp_norm_jac;
    else
        temp_norm_jac = jacobi;

    // J(1,0,0):

    xxi = node_pos[2] - node_pos[1];
    xet = node_pos[0] - node_pos[1];
    xze = node_pos[5] - node_pos[1];

    jacobi = xxi % ( xet * xze );
    // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    lengths       = sqrt( len1_sq * len2_sq * len3_sq );
    temp_norm_jac = jacobi / lengths;
    if( temp_norm_jac < min_norm_jac )
        min_norm_jac = temp_norm_jac;
    else
        temp_norm_jac = jacobi;

    // J(1,1,0):

    xxi = node_pos[3] - node_pos[2];
    xet = node_pos[1] - node_pos[2];
    xze = node_pos[6] - node_pos[2];

    jacobi = xxi % ( xet * xze );
    // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    lengths       = sqrt( len1_sq * len2_sq * len3_sq );
    temp_norm_jac = jacobi / lengths;
    if( temp_norm_jac < min_norm_jac )
        min_norm_jac = temp_norm_jac;
    else
        temp_norm_jac = jacobi;

    // J(0,1,0):

    xxi = node_pos[0] - node_pos[3];
    xet = node_pos[2] - node_pos[3];
    xze = node_pos[7] - node_pos[3];

    jacobi = xxi % ( xet * xze );
    // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    lengths       = sqrt( len1_sq * len2_sq * len3_sq );
    temp_norm_jac = jacobi / lengths;
    if( temp_norm_jac < min_norm_jac )
        min_norm_jac = temp_norm_jac;
    else
        temp_norm_jac = jacobi;

    // J(0,0,1):

    xxi = node_pos[7] - node_pos[4];
    xet = node_pos[5] - node_pos[4];
    xze = node_pos[0] - node_pos[4];

    jacobi = xxi % ( xet * xze );
    // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    lengths       = sqrt( len1_sq * len2_sq * len3_sq );
    temp_norm_jac = jacobi / lengths;
    if( temp_norm_jac < min_norm_jac )
        min_norm_jac = temp_norm_jac;
    else
        temp_norm_jac = jacobi;

    // J(1,0,1):

    xxi = node_pos[4] - node_pos[5];
    xet = node_pos[6] - node_pos[5];
    xze = node_pos[1] - node_pos[5];

    jacobi = xxi % ( xet * xze );
    // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    lengths       = sqrt( len1_sq * len2_sq * len3_sq );
    temp_norm_jac = jacobi / lengths;
    if( temp_norm_jac < min_norm_jac )
        min_norm_jac = temp_norm_jac;
    else
        temp_norm_jac = jacobi;

    // J(1,1,1):

    xxi = node_pos[5] - node_pos[6];
    xet = node_pos[7] - node_pos[6];
    xze = node_pos[2] - node_pos[6];

    jacobi = xxi % ( xet * xze );
    // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    lengths       = sqrt( len1_sq * len2_sq * len3_sq );
    temp_norm_jac = jacobi / lengths;
    if( temp_norm_jac < min_norm_jac )
        min_norm_jac = temp_norm_jac;
    else
        temp_norm_jac = jacobi;

    // J(0,1,1):

    xxi = node_pos[6] - node_pos[7];
    xet = node_pos[4] - node_pos[7];
    xze = node_pos[3] - node_pos[7];

    jacobi = xxi % ( xet * xze );
    // if( jacobi < min_jacobi ) { min_jacobi = jacobi; }

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    lengths       = sqrt( len1_sq * len2_sq * len3_sq );
    temp_norm_jac = jacobi / lengths;
    if( temp_norm_jac < min_norm_jac ) min_norm_jac = temp_norm_jac;
    // else
    // temp_norm_jac = jacobi;

    if( min_norm_jac > 0 ) return (double)VERDICT_MIN( min_norm_jac, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( min_norm_jac, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_shape ( int  ,
double  coordinates[][3] 
)

Calculates hex shape metric.

shape of a hex

3/Condition number of weighted Jacobian matrix

Definition at line 1817 of file V_HexMetric.cpp.

References MBMesquite::det(), make_hex_nodes, VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{

    double det, shape;
    double min_shape               = 1.0;
    static const double two_thirds = 2.0 / 3.0;

    VerdictVector xxi, xet, xze;

    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    // J(0,0,0):

    xxi = node_pos[1] - node_pos[0];
    xet = node_pos[3] - node_pos[0];
    xze = node_pos[4] - node_pos[0];

    det = xxi % ( xet * xze );
    if( det > VERDICT_DBL_MIN )
        shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
    else
        return 0;

    if( shape < min_shape ) { min_shape = shape; }

    // J(1,0,0):

    xxi = node_pos[2] - node_pos[1];
    xet = node_pos[0] - node_pos[1];
    xze = node_pos[5] - node_pos[1];

    det = xxi % ( xet * xze );
    if( det > VERDICT_DBL_MIN )
        shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
    else
        return 0;

    if( shape < min_shape ) { min_shape = shape; }

    // J(1,1,0):

    xxi = node_pos[3] - node_pos[2];
    xet = node_pos[1] - node_pos[2];
    xze = node_pos[6] - node_pos[2];

    det = xxi % ( xet * xze );
    if( det > VERDICT_DBL_MIN )
        shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
    else
        return 0;

    if( shape < min_shape ) { min_shape = shape; }

    // J(0,1,0):

    xxi = node_pos[0] - node_pos[3];
    xet = node_pos[2] - node_pos[3];
    xze = node_pos[7] - node_pos[3];

    det = xxi % ( xet * xze );
    if( det > VERDICT_DBL_MIN )
        shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
    else
        return 0;

    if( shape < min_shape ) { min_shape = shape; }

    // J(0,0,1):

    xxi = node_pos[7] - node_pos[4];
    xet = node_pos[5] - node_pos[4];
    xze = node_pos[0] - node_pos[4];

    det = xxi % ( xet * xze );
    if( det > VERDICT_DBL_MIN )
        shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
    else
        return 0;

    if( shape < min_shape ) { min_shape = shape; }

    // J(1,0,1):

    xxi = node_pos[4] - node_pos[5];
    xet = node_pos[6] - node_pos[5];
    xze = node_pos[1] - node_pos[5];

    det = xxi % ( xet * xze );
    if( det > VERDICT_DBL_MIN )
        shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
    else
        return 0;

    if( shape < min_shape ) { min_shape = shape; }

    // J(1,1,1):

    xxi = node_pos[5] - node_pos[6];
    xet = node_pos[7] - node_pos[6];
    xze = node_pos[2] - node_pos[6];

    det = xxi % ( xet * xze );
    if( det > VERDICT_DBL_MIN )
        shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
    else
        return 0;

    if( shape < min_shape ) { min_shape = shape; }

    // J(1,1,1):

    xxi = node_pos[6] - node_pos[7];
    xet = node_pos[4] - node_pos[7];
    xze = node_pos[3] - node_pos[7];

    det = xxi % ( xet * xze );
    if( det > VERDICT_DBL_MIN )
        shape = 3 * pow( det, two_thirds ) / ( xxi % xxi + xet % xet + xze % xze );
    else
        return 0;

    if( shape < min_shape ) { min_shape = shape; }

    if( min_shape <= VERDICT_DBL_MIN ) min_shape = 0;

    if( min_shape > 0 ) return (double)VERDICT_MIN( min_shape, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( min_shape, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_shape_and_size ( int  num_nodes,
double  coordinates[][3] 
)

Calculates hex shape-size metric.

shape and size of a hex

Product of Shape and Relative Size

Definition at line 2060 of file V_HexMetric.cpp.

References size, v_hex_relative_size_squared(), v_hex_shape(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    double size  = v_hex_relative_size_squared( num_nodes, coordinates );
    double shape = v_hex_shape( num_nodes, coordinates );

    double shape_size = size * shape;

    if( shape_size > 0 ) return (double)VERDICT_MIN( shape_size, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( shape_size, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_shear ( int  ,
double  coordinates[][3] 
)

Calculates hex shear metric.

shear of a hex

3/Condition number of Jacobian Skew matrix

Definition at line 1643 of file V_HexMetric.cpp.

References MBMesquite::det(), VerdictVector::length_squared(), make_hex_nodes, VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{

    double shear;
    double min_shear = 1.0;
    VerdictVector xxi, xet, xze;
    double det, len1_sq, len2_sq, len3_sq, lengths;

    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    // J(0,0,0):

    xxi = node_pos[1] - node_pos[0];
    xet = node_pos[3] - node_pos[0];
    xze = node_pos[4] - node_pos[0];

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;

    lengths = sqrt( len1_sq * len2_sq * len3_sq );
    det     = xxi % ( xet * xze );
    if( det < VERDICT_DBL_MIN ) { return 0; }

    shear     = det / lengths;
    min_shear = VERDICT_MIN( shear, min_shear );

    // J(1,0,0):

    xxi = node_pos[2] - node_pos[1];
    xet = node_pos[0] - node_pos[1];
    xze = node_pos[5] - node_pos[1];

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;

    lengths = sqrt( len1_sq * len2_sq * len3_sq );
    det     = xxi % ( xet * xze );
    if( det < VERDICT_DBL_MIN ) { return 0; }

    shear     = det / lengths;
    min_shear = VERDICT_MIN( shear, min_shear );

    // J(1,1,0):

    xxi = node_pos[3] - node_pos[2];
    xet = node_pos[1] - node_pos[2];
    xze = node_pos[6] - node_pos[2];

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;

    lengths = sqrt( len1_sq * len2_sq * len3_sq );
    det     = xxi % ( xet * xze );
    if( det < VERDICT_DBL_MIN ) { return 0; }

    shear     = det / lengths;
    min_shear = VERDICT_MIN( shear, min_shear );

    // J(0,1,0):

    xxi = node_pos[0] - node_pos[3];
    xet = node_pos[2] - node_pos[3];
    xze = node_pos[7] - node_pos[3];

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;

    lengths = sqrt( len1_sq * len2_sq * len3_sq );
    det     = xxi % ( xet * xze );
    if( det < VERDICT_DBL_MIN ) { return 0; }

    shear     = det / lengths;
    min_shear = VERDICT_MIN( shear, min_shear );

    // J(0,0,1):

    xxi = node_pos[7] - node_pos[4];
    xet = node_pos[5] - node_pos[4];
    xze = node_pos[0] - node_pos[4];

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;

    lengths = sqrt( len1_sq * len2_sq * len3_sq );
    det     = xxi % ( xet * xze );
    if( det < VERDICT_DBL_MIN ) { return 0; }

    shear     = det / lengths;
    min_shear = VERDICT_MIN( shear, min_shear );

    // J(1,0,1):

    xxi = node_pos[4] - node_pos[5];
    xet = node_pos[6] - node_pos[5];
    xze = node_pos[1] - node_pos[5];

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;

    lengths = sqrt( len1_sq * len2_sq * len3_sq );
    det     = xxi % ( xet * xze );
    if( det < VERDICT_DBL_MIN ) { return 0; }

    shear     = det / lengths;
    min_shear = VERDICT_MIN( shear, min_shear );

    // J(1,1,1):

    xxi = node_pos[5] - node_pos[6];
    xet = node_pos[7] - node_pos[6];
    xze = node_pos[2] - node_pos[6];

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;

    lengths = sqrt( len1_sq * len2_sq * len3_sq );
    det     = xxi % ( xet * xze );
    if( det < VERDICT_DBL_MIN ) { return 0; }

    shear     = det / lengths;
    min_shear = VERDICT_MIN( shear, min_shear );

    // J(0,1,1):

    xxi = node_pos[6] - node_pos[7];
    xet = node_pos[4] - node_pos[7];
    xze = node_pos[3] - node_pos[7];

    len1_sq = xxi.length_squared();
    len2_sq = xet.length_squared();
    len3_sq = xze.length_squared();

    if( len1_sq <= VERDICT_DBL_MIN || len2_sq <= VERDICT_DBL_MIN || len3_sq <= VERDICT_DBL_MIN ) return 0;

    lengths = sqrt( len1_sq * len2_sq * len3_sq );
    det     = xxi % ( xet * xze );
    if( det < VERDICT_DBL_MIN ) { return 0; }

    shear     = det / lengths;
    min_shear = VERDICT_MIN( shear, min_shear );

    if( min_shear <= VERDICT_DBL_MIN ) min_shear = 0;

    if( min_shear > 0 ) return (double)VERDICT_MIN( min_shear, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( min_shear, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_shear_and_size ( int  num_nodes,
double  coordinates[][3] 
)

Calculates hex shear-size metric.

shear and size of a hex

Product of Shear and Relative Size

Definition at line 2076 of file V_HexMetric.cpp.

References size, v_hex_relative_size_squared(), v_hex_shear(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    double size  = v_hex_relative_size_squared( num_nodes, coordinates );
    double shear = v_hex_shear( num_nodes, coordinates );

    double shear_size = shear * size;

    if( shear_size > 0 ) return (double)VERDICT_MIN( shear_size, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( shear_size, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_skew ( int  ,
double  coordinates[][3] 
)

Calculates hex skew metric.

skew of a hex

Maximum ||cosA|| where A is the angle between edges at hex center.

Definition at line 665 of file V_HexMetric.cpp.

References calc_hex_efg(), make_hex_nodes, VerdictVector::normalize(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{
    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    double skew_1, skew_2, skew_3;

    VerdictVector efg1 = calc_hex_efg( 1, node_pos );
    VerdictVector efg2 = calc_hex_efg( 2, node_pos );
    VerdictVector efg3 = calc_hex_efg( 3, node_pos );

    if( efg1.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
    if( efg2.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;
    if( efg3.normalize() <= VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;

    skew_1 = fabs( efg1 % efg2 );
    skew_2 = fabs( efg1 % efg3 );
    skew_3 = fabs( efg2 % efg3 );

    double skew = ( VERDICT_MAX( skew_1, VERDICT_MAX( skew_2, skew_3 ) ) );

    if( skew > 0 ) return (double)VERDICT_MIN( skew, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( skew, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_stretch ( int  ,
double  coordinates[][3] 
)

Calculates hex stretch metric.

stretch of a hex

sqrt(3) * minimum edge length / maximum diagonal length

Definition at line 744 of file V_HexMetric.cpp.

References diag_length(), hex_edge_length(), safe_ratio(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    static const double HEX_STRETCH_SCALE_FACTOR = sqrt( 3.0 );

    double min_edge = hex_edge_length( 0, coordinates );
    double max_diag = diag_length( 1, coordinates );

    double stretch = HEX_STRETCH_SCALE_FACTOR * safe_ratio( min_edge, max_diag );

    if( stretch > 0 ) return (double)VERDICT_MIN( stretch, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( stretch, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_taper ( int  ,
double  coordinates[][3] 
)

Calculates hex taper metric.

taper of a hex

Maximum ratio of lengths derived from opposite edges.

Definition at line 695 of file V_HexMetric.cpp.

References calc_hex_efg(), VerdictVector::length(), make_hex_nodes, safe_ratio(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    VerdictVector efg1 = calc_hex_efg( 1, node_pos );
    VerdictVector efg2 = calc_hex_efg( 2, node_pos );
    VerdictVector efg3 = calc_hex_efg( 3, node_pos );

    VerdictVector efg12 = calc_hex_efg( 12, node_pos );
    VerdictVector efg13 = calc_hex_efg( 13, node_pos );
    VerdictVector efg23 = calc_hex_efg( 23, node_pos );

    double taper_1 = fabs( safe_ratio( efg12.length(), VERDICT_MIN( efg1.length(), efg2.length() ) ) );
    double taper_2 = fabs( safe_ratio( efg13.length(), VERDICT_MIN( efg1.length(), efg3.length() ) ) );
    double taper_3 = fabs( safe_ratio( efg23.length(), VERDICT_MIN( efg2.length(), efg3.length() ) ) );

    double taper = (double)VERDICT_MAX( taper_1, VERDICT_MAX( taper_2, taper_3 ) );

    if( taper > 0 ) return (double)VERDICT_MIN( taper, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( taper, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_hex_volume ( int  ,
double  coordinates[][3] 
)

Calculates hex volume.

volume of a hex

Jacobian at hex center

Definition at line 723 of file V_HexMetric.cpp.

References calc_hex_efg(), make_hex_nodes, VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    VerdictVector node_pos[8];
    make_hex_nodes( coordinates, node_pos );

    VerdictVector efg1 = calc_hex_efg( 1, node_pos );
    VerdictVector efg2 = calc_hex_efg( 2, node_pos );
    VerdictVector efg3 = calc_hex_efg( 3, node_pos );

    double volume;
    volume = (double)( efg1 % ( efg2 * efg3 ) ) / 64.0;

    if( volume > 0 ) return (double)VERDICT_MIN( volume, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( volume, -VERDICT_DBL_MAX );
}
C_FUNC_DEF void v_set_hex_size ( double  size)

returns the average volume of a hex

Sets average size (volume) of hex, needed for v_hex_relative_size(...)

Definition at line 57 of file V_HexMetric.cpp.

References size, and verdict_hex_size.

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


Variable Documentation

double verdict_hex_size = 0

the average volume of a hex

Definition at line 36 of file V_HexMetric.cpp.

Referenced by v_hex_get_weight(), and v_set_hex_size().

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines