MOAB: Mesh Oriented datABase  (version 5.2.1)
verdict.h File Reference

Header file for verdict library that calculates metrics for finite elements. Also see: Main Page. More...

+ This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes

struct  HexMetricVals
struct  EdgeMetricVals
struct  KnifeMetricVals
struct  QuadMetricVals
struct  PyramidMetricVals
struct  WedgeMetricVals
struct  TetMetricVals
struct  TriMetricVals

Defines

#define VERDICT_VERSION   120
#define VERDICT_DBL_MIN   1.0E-30
#define VERDICT_DBL_MAX   1.0E+30
#define VERDICT_PI   3.1415926535897932384626
#define C_FUNC_DEF
Hex bit fields
#define V_HEX_MAX_EDGE_RATIO
#define V_HEX_SKEW
#define V_HEX_TAPER
#define V_HEX_VOLUME
#define V_HEX_STRETCH
#define V_HEX_DIAGONAL
#define V_HEX_DIMENSION
#define V_HEX_ODDY
#define V_HEX_MAX_ASPECT_FROBENIUS
#define V_HEX_CONDITION
#define V_HEX_JACOBIAN
#define V_HEX_SCALED_JACOBIAN
#define V_HEX_SHEAR
#define V_HEX_SHAPE
#define V_HEX_RELATIVE_SIZE_SQUARED
#define V_HEX_SHAPE_AND_SIZE
#define V_HEX_SHEAR_AND_SIZE
#define V_HEX_DISTORTION
#define V_HEX_EDGE_RATIO
#define V_HEX_MED_ASPECT_FROBENIUS
#define V_HEX_ALL
#define V_HEX_TRADITIONAL
#define V_HEX_DIAGNOSTIC
#define V_HEX_ALGEBRAIC
#define V_HEX_ROBINSON
Tet bit fields
#define V_TET_RADIUS_RATIO
#define V_TET_ASPECT_BETA
#define V_TET_ASPECT_GAMMA
#define V_TET_VOLUME
#define V_TET_CONDITION
#define V_TET_JACOBIAN
#define V_TET_SCALED_JACOBIAN
#define V_TET_SHAPE
#define V_TET_RELATIVE_SIZE_SQUARED
#define V_TET_SHAPE_AND_SIZE
#define V_TET_DISTORTION
#define V_TET_EDGE_RATIO
#define V_TET_ASPECT_RATIO
#define V_TET_ASPECT_FROBENIUS
#define V_TET_MINIMUM_ANGLE
#define V_TET_COLLAPSE_RATIO
#define V_TET_ALL
#define V_TET_TRADITIONAL
#define V_TET_DIAGNOSTIC
#define V_TET_ALGEBRAIC
Pyramid bit fields
#define V_PYRAMID_VOLUME
Wedge bit fields
#define V_WEDGE_VOLUME
Knife bit fields
#define V_KNIFE_VOLUME
Quad bit fields
#define V_QUAD_MAX_EDGE_RATIO
#define V_QUAD_SKEW
#define V_QUAD_TAPER
#define V_QUAD_WARPAGE
#define V_QUAD_AREA
#define V_QUAD_STRETCH
#define V_QUAD_MINIMUM_ANGLE
#define V_QUAD_MAXIMUM_ANGLE
#define V_QUAD_ODDY
#define V_QUAD_CONDITION
#define V_QUAD_JACOBIAN
#define V_QUAD_SCALED_JACOBIAN
#define V_QUAD_SHEAR
#define V_QUAD_SHAPE
#define V_QUAD_RELATIVE_SIZE_SQUARED
#define V_QUAD_SHAPE_AND_SIZE
#define V_QUAD_SHEAR_AND_SIZE
#define V_QUAD_DISTORTION
#define V_QUAD_EDGE_RATIO
#define V_QUAD_ASPECT_RATIO
#define V_QUAD_RADIUS_RATIO
#define V_QUAD_MED_ASPECT_FROBENIUS
#define V_QUAD_MAX_ASPECT_FROBENIUS
#define V_QUAD_ALL
#define V_QUAD_TRADITIONAL
#define V_QUAD_DIAGNOSTIC
#define V_QUAD_ALGEBRAIC
#define V_QUAD_ROBINSON
Tri bit fields
#define V_TRI_ASPECT_FROBENIUS
#define V_TRI_AREA
#define V_TRI_MINIMUM_ANGLE
#define V_TRI_MAXIMUM_ANGLE
#define V_TRI_CONDITION
#define V_TRI_SCALED_JACOBIAN
#define V_TRI_SHAPE
#define V_TRI_RELATIVE_SIZE_SQUARED
#define V_TRI_SHAPE_AND_SIZE
#define V_TRI_DISTORTION
#define V_TRI_RADIUS_RATIO
#define V_TRI_EDGE_RATIO
#define V_TRI_ALL
#define V_TRI_TRADITIONAL
#define V_TRI_DIAGNOSTIC
#define V_TRI_ALGEBRAIC
#define V_EDGE_LENGTH

Typedefs

typedef double(* VerdictFunction )(int, double[][3])
typedef int(* ComputeNormal )(double point[3], double normal[3])

Functions

C_FUNC_DEF void v_hex_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, struct HexMetricVals *metric_vals)
 Calculates quality metrics for hexahedral elements.
C_FUNC_DEF void v_tet_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, struct TetMetricVals *metric_vals)
 Calculates quality metrics for tetrahedral elements.
C_FUNC_DEF void v_pyramid_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, struct PyramidMetricVals *metric_vals)
 Calculates quality metrics for pyramid elements.
C_FUNC_DEF void v_wedge_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, struct WedgeMetricVals *metric_vals)
 Calculates quality metrics for wedge elements.
C_FUNC_DEF void v_knife_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, struct KnifeMetricVals *metric_vals)
 Calculates quality metrics for knife elements.
C_FUNC_DEF void v_quad_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, struct QuadMetricVals *metric_vals)
 Calculates quality metrics for quadrilateral elements.
C_FUNC_DEF void v_tri_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, struct TriMetricVals *metric_vals)
 Calculates quality metrics for triangle elements.
C_FUNC_DEF void v_edge_quality (int num_nodes, double coordinates[][3], unsigned int metrics_request_flag, struct EdgeMetricVals *metric_vals)
 Calculates quality metrics for edge elements.
C_FUNC_DEF void v_set_hex_size (double size)
 Sets average size (volume) of hex, needed for v_hex_relative_size(...)
C_FUNC_DEF double v_hex_edge_ratio (int num_nodes, double coordinates[][3])
 Calculates hex edge ratio metric.
C_FUNC_DEF double v_hex_max_edge_ratio (int num_nodes, double coordinates[][3])
 Calculates hex maximum of edge ratio.
C_FUNC_DEF double v_hex_skew (int num_nodes, double coordinates[][3])
 Calculates hex skew metric.
C_FUNC_DEF double v_hex_taper (int num_nodes, double coordinates[][3])
 Calculates hex taper metric.
C_FUNC_DEF double v_hex_volume (int num_nodes, double coordinates[][3])
 Calculates hex volume.
C_FUNC_DEF double v_hex_stretch (int num_nodes, double coordinates[][3])
 Calculates hex stretch metric.
C_FUNC_DEF double v_hex_diagonal (int num_nodes, double coordinates[][3])
 Calculates hex diagonal metric.
C_FUNC_DEF double v_hex_dimension (int num_nodes, double coordinates[][3])
 Calculates hex dimension metric.
C_FUNC_DEF double v_hex_oddy (int num_nodes, double coordinates[][3])
 Calculates hex oddy metric.
C_FUNC_DEF double v_hex_med_aspect_frobenius (int num_nodes, double coordinates[][3])
 Calculates hex condition metric.
C_FUNC_DEF double v_hex_max_aspect_frobenius (int num_nodes, double coordinates[][3])
 Calculates hex condition metric.
C_FUNC_DEF double v_hex_condition (int num_nodes, double coordinates[][3])
C_FUNC_DEF double v_hex_jacobian (int num_nodes, double coordinates[][3])
 Calculates hex jacobian metric.
C_FUNC_DEF double v_hex_scaled_jacobian (int num_nodes, double coordinates[][3])
 Calculates hex scaled jacobian metric.
C_FUNC_DEF double v_hex_shear (int num_nodes, double coordinates[][3])
 Calculates hex shear metric.
C_FUNC_DEF double v_hex_shape (int num_nodes, double coordinates[][3])
 Calculates hex shape metric.
C_FUNC_DEF double v_hex_relative_size_squared (int num_nodes, 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_set_tet_size (double size)
 Sets average size (volume) of tet, needed for v_tet_relative_size(...)
C_FUNC_DEF double v_tet_edge_ratio (int num_nodes, double coordinates[][3])
 Calculates tet edge ratio metric.
C_FUNC_DEF double v_tet_radius_ratio (int num_nodes, double coordinates[][3])
 Calculates tet radius ratio metric.
C_FUNC_DEF double v_tet_aspect_beta (int num_nodes, double coordinates[][3])
 Calculates the radius ratio metric of a positively oriented tet.
C_FUNC_DEF double v_tet_aspect_ratio (int num_nodes, double coordinates[][3])
 Calculates tet aspect ratio metric.
C_FUNC_DEF double v_tet_aspect_gamma (int num_nodes, double coordinates[][3])
 Calculates tet aspect gamma metric.
C_FUNC_DEF double v_tet_aspect_frobenius (int num_nodes, double coordinates[][3])
 Calculates tet aspect frobenius metric.
C_FUNC_DEF double v_tet_minimum_angle (int num_nodes, double coordinates[][3])
 Calculates tet minimum dihedral angle.
C_FUNC_DEF double v_tet_collapse_ratio (int num_nodes, double coordinates[][3])
 Calculates tet collapse ratio metric.
C_FUNC_DEF double v_tet_volume (int num_nodes, double coordinates[][3])
 Calculates tet volume.
C_FUNC_DEF double v_tet_condition (int num_nodes, double coordinates[][3])
 Calculates tet condition metric.
C_FUNC_DEF double v_tet_jacobian (int num_nodes, double coordinates[][3])
 Calculates tet jacobian.
C_FUNC_DEF double v_tet_scaled_jacobian (int num_nodes, double coordinates[][3])
 Calculates tet scaled jacobian.
C_FUNC_DEF double v_tet_shape (int num_nodes, double coordinates[][3])
 Calculates tet shape metric.
C_FUNC_DEF double v_tet_relative_size_squared (int num_nodes, double coordinates[][3])
 Calculates tet relative size metric.
C_FUNC_DEF double v_tet_shape_and_size (int num_nodes, double coordinates[][3])
 Calculates tet shape-size metric.
C_FUNC_DEF double v_tet_distortion (int num_nodes, double coordinates[][3])
 Calculates tet distortion metric.
C_FUNC_DEF double v_pyramid_volume (int num_nodes, double coordinates[][3])
 Calculates pyramid volume.
C_FUNC_DEF double v_wedge_volume (int num_nodes, double coordinates[][3])
 Calculates wedge volume.
C_FUNC_DEF double v_knife_volume (int num_nodes, double coordinates[][3])
 Calculates knife volume.
C_FUNC_DEF double v_edge_length (int num_nodes, double coordinates[][3])
 Calculates edge length.
C_FUNC_DEF void v_set_quad_size (double size)
 Sets average size (area) of quad, needed for v_quad_relative_size(...)
C_FUNC_DEF double v_quad_edge_ratio (int num_nodes, double coordinates[][3])
 Calculates quad edge ratio.
C_FUNC_DEF double v_quad_max_edge_ratio (int num_nodes, double coordinates[][3])
 Calculates quad maximum of edge ratio.
C_FUNC_DEF double v_quad_aspect_ratio (int num_nodes, double coordinates[][3])
 Calculates quad aspect ratio.
C_FUNC_DEF double v_quad_radius_ratio (int num_nodes, double coordinates[][3])
 Calculates quad radius ratio.
C_FUNC_DEF double v_quad_med_aspect_frobenius (int num_nodes, double coordinates[][3])
 Calculates quad average Frobenius aspect.
C_FUNC_DEF double v_quad_max_aspect_frobenius (int num_nodes, double coordinates[][3])
 Calculates quad maximum Frobenius aspect.
C_FUNC_DEF double v_quad_skew (int num_nodes, double coordinates[][3])
 Calculates quad skew metric.
C_FUNC_DEF double v_quad_taper (int num_nodes, double coordinates[][3])
 Calculates quad taper metric.
C_FUNC_DEF double v_quad_warpage (int num_nodes, double coordinates[][3])
 Calculates quad warpage metric.
C_FUNC_DEF double v_quad_area (int num_nodes, double coordinates[][3])
 Calculates quad area.
C_FUNC_DEF double v_quad_stretch (int num_nodes, double coordinates[][3])
 Calculates quad strech metric.
C_FUNC_DEF double v_quad_minimum_angle (int num_nodes, double coordinates[][3])
 Calculates quad's smallest angle.
C_FUNC_DEF double v_quad_maximum_angle (int num_nodes, double coordinates[][3])
 Calculates quad's largest angle.
C_FUNC_DEF double v_quad_oddy (int num_nodes, double coordinates[][3])
 Calculates quad oddy metric.
C_FUNC_DEF double v_quad_condition (int num_nodes, double coordinates[][3])
 Calculates quad condition number metric.
C_FUNC_DEF double v_quad_jacobian (int num_nodes, double coordinates[][3])
 Calculates quad jacobian.
C_FUNC_DEF double v_quad_scaled_jacobian (int num_nodes, double coordinates[][3])
 Calculates quad scaled jacobian.
C_FUNC_DEF double v_quad_shear (int num_nodes, double coordinates[][3])
 Calculates quad shear metric.
C_FUNC_DEF double v_quad_shape (int num_nodes, double coordinates[][3])
 Calculates quad shape metric.
C_FUNC_DEF double v_quad_relative_size_squared (int num_nodes, double coordinates[][3])
 Calculates quad relative size metric.
C_FUNC_DEF double v_quad_shape_and_size (int num_nodes, double coordinates[][3])
 Calculates quad shape-size metric.
C_FUNC_DEF double v_quad_shear_and_size (int num_nodes, double coordinates[][3])
 Calculates quad shear-size metric.
C_FUNC_DEF double v_quad_distortion (int num_nodes, double coordinates[][3])
 Calculates quad distortion metric.
C_FUNC_DEF void v_set_tri_size (double size)
 Sets average size (area) of tri, needed for v_tri_relative_size(...)
C_FUNC_DEF void v_set_tri_normal_func (ComputeNormal func)
 Sets fuction pointer to calculate tri normal wrt surface.
C_FUNC_DEF double v_tri_edge_ratio (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_aspect_ratio (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_radius_ratio (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_aspect_frobenius (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_area (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_minimum_angle (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_maximum_angle (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_condition (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_scaled_jacobian (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_relative_size_squared (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_shape (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_shape_and_size (int num_nodes, double coordinates[][3])
 Calculates tri metric.
C_FUNC_DEF double v_tri_distortion (int num_nodes, double coordinates[][3])
 Calculates tri metric.

Detailed Description

Header file for verdict library that calculates metrics for finite elements. Also see: Main Page.

verdict.h is the header file for applications/libraries to include to compute quality metrics.

This file is part of VERDICT

Definition in file verdict.h.


Define Documentation

#define C_FUNC_DEF

Definition at line 59 of file verdict.h.

#define V_EDGE_LENGTH

Definition at line 508 of file verdict.h.

Referenced by edge_quality().

#define V_HEX_ALGEBRAIC

Definition at line 363 of file verdict.h.

#define V_HEX_ALL

Definition at line 347 of file verdict.h.

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

#define V_HEX_CONDITION

Definition at line 336 of file verdict.h.

Referenced by v_hex_quality().

Definition at line 361 of file verdict.h.

#define V_HEX_DIAGONAL

Definition at line 332 of file verdict.h.

Referenced by v_hex_quality().

#define V_HEX_DIMENSION

Definition at line 333 of file verdict.h.

Referenced by v_hex_quality().

Definition at line 344 of file verdict.h.

Referenced by v_hex_quality().

Definition at line 345 of file verdict.h.

Referenced by v_hex_quality().

#define V_HEX_JACOBIAN

Definition at line 337 of file verdict.h.

Referenced by v_hex_quality().

Definition at line 335 of file verdict.h.

Definition at line 327 of file verdict.h.

Referenced by v_hex_quality().

Definition at line 346 of file verdict.h.

Referenced by v_hex_quality().

#define V_HEX_ODDY

Definition at line 334 of file verdict.h.

Referenced by v_hex_quality().

Definition at line 341 of file verdict.h.

Referenced by v_hex_quality().

#define V_HEX_ROBINSON
Value:
(V_HEX_SKEW + \
                                     V_HEX_TAPER)

Definition at line 369 of file verdict.h.

Definition at line 338 of file verdict.h.

Referenced by v_hex_quality().

#define V_HEX_SHAPE

Definition at line 340 of file verdict.h.

Referenced by v_hex_quality().

Definition at line 342 of file verdict.h.

Referenced by v_hex_quality().

#define V_HEX_SHEAR

Definition at line 339 of file verdict.h.

Referenced by v_hex_quality().

Definition at line 343 of file verdict.h.

Referenced by v_hex_quality().

#define V_HEX_SKEW

Definition at line 328 of file verdict.h.

Referenced by v_hex_quality().

#define V_HEX_STRETCH

Definition at line 331 of file verdict.h.

Referenced by v_hex_quality().

#define V_HEX_TAPER

Definition at line 329 of file verdict.h.

Referenced by v_hex_quality().

Definition at line 349 of file verdict.h.

#define V_HEX_VOLUME

Definition at line 330 of file verdict.h.

Referenced by v_hex_quality().

#define V_KNIFE_VOLUME

Definition at line 423 of file verdict.h.

Referenced by v_knife_quality().

Definition at line 411 of file verdict.h.

Referenced by v_pyramid_quality().

Definition at line 468 of file verdict.h.

#define V_QUAD_ALL

Definition at line 452 of file verdict.h.

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

#define V_QUAD_AREA

Definition at line 433 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 448 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 438 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 466 of file verdict.h.

Definition at line 446 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 447 of file verdict.h.

Referenced by v_quad_quality().

#define V_QUAD_JACOBIAN

Definition at line 439 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 451 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 429 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 436 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 450 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 435 of file verdict.h.

Referenced by v_quad_quality().

#define V_QUAD_ODDY

Definition at line 437 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 449 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 443 of file verdict.h.

Referenced by v_quad_quality().

#define V_QUAD_ROBINSON
Value:
(V_QUAD_MAX_EDGE_RATIO + \
                                     V_QUAD_SKEW   + \
                                     V_QUAD_TAPER)

Definition at line 473 of file verdict.h.

Definition at line 440 of file verdict.h.

Referenced by v_quad_quality().

#define V_QUAD_SHAPE

Definition at line 442 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 444 of file verdict.h.

Referenced by v_quad_quality().

#define V_QUAD_SHEAR

Definition at line 441 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 445 of file verdict.h.

Referenced by v_quad_quality().

#define V_QUAD_SKEW

Definition at line 430 of file verdict.h.

Referenced by v_quad_quality().

#define V_QUAD_STRETCH

Definition at line 434 of file verdict.h.

Referenced by v_quad_quality().

#define V_QUAD_TAPER

Definition at line 431 of file verdict.h.

Referenced by v_quad_quality().

Definition at line 454 of file verdict.h.

#define V_QUAD_WARPAGE

Definition at line 432 of file verdict.h.

Referenced by v_quad_quality().

#define V_TET_ALGEBRAIC
Value:
(V_TET_SHAPE                  + \
                                     V_TET_RELATIVE_SIZE_SQUARED  + \
                                     V_TET_SHAPE_AND_SIZE)

Definition at line 402 of file verdict.h.

#define V_TET_ALL

Definition at line 392 of file verdict.h.

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

Definition at line 377 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 389 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 378 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 388 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 391 of file verdict.h.

Referenced by v_tet_quality().

#define V_TET_CONDITION

Definition at line 380 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 400 of file verdict.h.

Definition at line 386 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 387 of file verdict.h.

#define V_TET_JACOBIAN

Definition at line 381 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 390 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 376 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 384 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 382 of file verdict.h.

Referenced by v_tet_quality().

#define V_TET_SHAPE

Definition at line 383 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 385 of file verdict.h.

Referenced by v_tet_quality().

Definition at line 394 of file verdict.h.

#define V_TET_VOLUME

Definition at line 379 of file verdict.h.

Referenced by v_tet_quality().

#define V_TRI_ALGEBRAIC
Value:
(V_TRI_SHAPE + \
                                     V_TRI_SHAPE_AND_SIZE + \
                                     V_TRI_RELATIVE_SIZE_SQUARED)

Definition at line 504 of file verdict.h.

#define V_TRI_ALL

Definition at line 494 of file verdict.h.

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

#define V_TRI_AREA

Definition at line 483 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 482 of file verdict.h.

Referenced by v_tri_quality().

#define V_TRI_CONDITION

Definition at line 486 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 502 of file verdict.h.

Definition at line 491 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 493 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 485 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 484 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 492 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 489 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 487 of file verdict.h.

Referenced by v_tri_quality().

#define V_TRI_SHAPE

Definition at line 488 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 490 of file verdict.h.

Referenced by v_tri_quality().

Definition at line 496 of file verdict.h.

#define V_WEDGE_VOLUME

Definition at line 417 of file verdict.h.

Referenced by v_wedge_quality().

#define VERDICT_DBL_MAX   1.0E+30

Definition at line 38 of file verdict.h.

Referenced by condition_comp(), oddy_comp(), safe_ratio(), v_hex_diagonal(), v_hex_distortion(), v_hex_edge_ratio(), 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_shape_and_size(), v_hex_shear(), v_hex_shear_and_size(), v_hex_skew(), v_hex_stretch(), v_hex_taper(), v_hex_volume(), v_quad_area(), v_quad_aspect_ratio(), v_quad_condition(), v_quad_distortion(), v_quad_edge_ratio(), v_quad_jacobian(), v_quad_max_aspect_frobenius(), v_quad_max_edge_ratio(), v_quad_maximum_angle(), v_quad_med_aspect_frobenius(), v_quad_minimum_angle(), v_quad_oddy(), v_quad_quality(), v_quad_radius_ratio(), v_quad_relative_size_squared(), v_quad_scaled_jacobian(), v_quad_shape(), v_quad_shape_and_size(), v_quad_shear(), v_quad_shear_and_size(), v_quad_skew(), v_quad_stretch(), v_quad_taper(), v_quad_warpage(), v_tet_aspect_beta(), v_tet_aspect_frobenius(), v_tet_aspect_gamma(), v_tet_aspect_ratio(), v_tet_collapse_ratio(), v_tet_condition(), v_tet_distortion(), v_tet_edge_ratio(), v_tet_minimum_angle(), v_tet_quality(), v_tet_radius_ratio(), v_tet_scaled_jacobian(), v_tri_area(), v_tri_aspect_frobenius(), v_tri_aspect_ratio(), v_tri_condition(), v_tri_distortion(), v_tri_edge_ratio(), v_tri_maximum_angle(), v_tri_minimum_angle(), v_tri_quality(), v_tri_radius_ratio(), v_tri_relative_size_squared(), v_tri_scaled_jacobian(), v_tri_shape(), and v_tri_shape_and_size().

#define VERDICT_VERSION   120

Definition at line 34 of file verdict.h.


Typedef Documentation

typedef int(* ComputeNormal)(double point[3], double normal[3])

Definition at line 70 of file verdict.h.

typedef double(* VerdictFunction)(int, double[][3])

Definition at line 69 of file verdict.h.


Function Documentation

C_FUNC_DEF double v_edge_length ( int  ,
double  coordinates[][3] 
)

Calculates edge length.

length of and edge length is calculated by taking the distance between the end nodes

Definition at line 32 of file V_EdgeMetric.cpp.

References z.

Referenced by moab::VerdictWrapper::all_quality_measures(), edge_quality(), and moab::VerdictWrapper::quality_measure().

{

    double x = coordinates[1][0] - coordinates[0][0];
    double y = coordinates[1][1] - coordinates[0][1];
    double z = coordinates[1][2] - coordinates[0][2];
    return (double)( sqrt( x * x + y * y + z * z ) );
}
C_FUNC_DEF void v_edge_quality ( int  num_nodes,
double  coordinates[][3],
unsigned int  metrics_request_flag,
struct EdgeMetricVals metric_vals 
)

Calculates quality metrics for edge elements.

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.

Minimum diagonal length / maximum diagonal length. Reference --- Unknown

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.

Pronto-specific characteristic length for stable time step calculation. Char_length = Volume / 2 grad Volume. Reference --- L.M. Taylor, and D.P. Flanagan, Pronto3D - A Three Dimensional Transient Solid Dynamics Program, SAND87-1912, Sandia National Laboratories, 1989.

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.

{min(|J|)/actual volume}*parent volume, parent volume = 8 for hex. Reference --- SDRC/IDEAS Simulation: Finite Element Modeling--User's Guide

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.

Hmax / Hmin where Hmax and Hmin are respectively the maximum and the minimum edge lengths

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 );
}
C_FUNC_DEF double v_hex_jacobian ( int  ,
double  coordinates[][3] 
)

Calculates hex jacobian metric.

Minimum pointwise volume of local map at 8 corners & center of hex. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

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 the Jacobian matrix at 8 corners. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

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.

Maximum edge length ratio at hex center. Reference --- L.M. Taylor, and D.P. Flanagan, Pronto3D - A Three Dimensional Transient Solid Dynamics Program, SAND87-1912, Sandia National Laboratories, 1989.

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.

Average Frobenius condition number of the Jacobian matrix at 8 corners.

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.

3/Mean Ratio of weighted Jacobian matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

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.

Minimum Jacobian divided by the lengths of the 3 edge vectors. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

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.

3/Mean Ratio of weighted Jacobian matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

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.

Product of Shape and Relative Size. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

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.

3/Mean Ratio of Jacobian Skew matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

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.

Product of Shear and Relative Size. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

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.

Maximum |cos A| where A is the angle between edges at hex center. Reference --- L.M. Taylor, and D.P. Flanagan, Pronto3D - A Three Dimensional Transient Solid Dynamics Program, SAND87-1912, Sandia National Laboratories, 1989.

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.

Sqrt(3) * minimum edge length / maximum diagonal length. Reference --- FIMESH code

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.

Maximum ratio of lengths derived from opposite edges. Reference --- L.M. Taylor, and D.P. Flanagan, Pronto3D - A Three Dimensional Transient Solid Dynamics Program, SAND87-1912, Sandia National Laboratories, 1989.

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.

Jacobian at hex center. Reference --- L.M. Taylor, and D.P. Flanagan, Pronto3D - A Three Dimensional Transient Solid Dynamics Program, SAND87-1912, Sandia National Laboratories, 1989.

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_knife_quality ( int  num_nodes,
double  coordinates[][3],
unsigned int  metrics_request_flag,
KnifeMetricVals metric_vals 
)

Calculates quality metrics for knife elements.

calculate the quality metrics of a knife element.

There is only one, but we put this here to be consistent with functions for other element types. Who knows if we'll add more metrics.

Definition at line 117 of file V_KnifeMetric.cpp.

References V_KNIFE_VOLUME, v_knife_volume(), and KnifeMetricVals::volume.

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

    if( metrics_request_flag & V_KNIFE_VOLUME ) metric_vals->volume = v_knife_volume( num_nodes, coordinates );
}
C_FUNC_DEF double v_knife_volume ( int  num_nodes,
double  coordinates[][3] 
)

Calculates knife volume.

calculates the volume of a knife element

this is done by dividing the knife into 4 tets and summing the volumes of each.

Definition at line 58 of file V_KnifeMetric.cpp.

References VerdictVector::set().

Referenced by moab::VerdictWrapper::all_quality_measures(), moab::VerdictWrapper::quality_measure(), and v_knife_quality().

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

    if( num_nodes == 7 )
    {

        // divide the knife into 4 tets and calculate the volume

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

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

        side1.set( coordinates[5][0] - coordinates[1][0], coordinates[5][1] - coordinates[1][1],
                   coordinates[5][2] - coordinates[1][2] );
        side2.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
                   coordinates[3][2] - coordinates[1][2] );
        side3.set( coordinates[4][0] - coordinates[1][0], coordinates[4][1] - coordinates[1][1],
                   coordinates[4][2] - coordinates[1][2] );

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

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

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

        side1.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
                   coordinates[3][2] - coordinates[1][2] );
        side2.set( coordinates[5][0] - coordinates[1][0], coordinates[5][1] - coordinates[1][1],
                   coordinates[5][2] - coordinates[1][2] );
        side3.set( coordinates[6][0] - coordinates[1][0], coordinates[6][1] - coordinates[1][1],
                   coordinates[6][2] - coordinates[1][2] );

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

    return (double)volume;
}
C_FUNC_DEF void v_pyramid_quality ( int  num_nodes,
double  coordinates[][3],
unsigned int  metrics_request_flag,
struct PyramidMetricVals metric_vals 
)

Calculates quality metrics for pyramid elements.

Definition at line 96 of file V_PyramidMetric.cpp.

References V_PYRAMID_VOLUME, v_pyramid_volume(), and PyramidMetricVals::volume.

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

    if( metrics_request_flag & V_PYRAMID_VOLUME ) metric_vals->volume = v_pyramid_volume( num_nodes, coordinates );
}
C_FUNC_DEF double v_pyramid_volume ( int  num_nodes,
double  coordinates[][3] 
)

Calculates pyramid volume.

the volume of a pyramid

the volume is calculated by dividing the pyramid into 2 tets and summing the volumes of the 2 tets.

Definition at line 59 of file V_PyramidMetric.cpp.

References VerdictVector::set().

Referenced by v_pyramid_quality().

{

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

    if( num_nodes == 5 )
    {
        // divide the pyramid into 2 tets and calculate each

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

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

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

        // volume of the first tet
        volume = ( side3 % ( side1 * side2 ) ) / 6.0;

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

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

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

        // volume of the second tet
        volume += ( side3 % ( side1 * side2 ) ) / 6.0;
    }
    return (double)volume;
}
C_FUNC_DEF double v_quad_area ( int  ,
double  coordinates[][3] 
)

Calculates quad area.

Jacobian at quad center. Reference --- J. Robinson, CRE Method of element testing and the Jacobian shape parameters, Eng. Comput., Vol 4, 1987.

the area of a quad

jacobian at quad center

Definition at line 611 of file V_QuadMetric.cpp.

References signed_corner_areas(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure(), v_quad_quality(), and v_quad_relative_size_squared().

{

    double corner_areas[4];
    signed_corner_areas( corner_areas, coordinates );

    double area = 0.25 * ( corner_areas[0] + corner_areas[1] + corner_areas[2] + corner_areas[3] );

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

Calculates quad aspect ratio.

aspect ratio Reference --- P. P. Pebay, Planar Quadrangle Quality Measures, Eng. Comp., 2004, 20(2):157-173

the aspect ratio of a quad

NB (P. Pebay 01/20/07): this is a generalization of the triangle aspect ratio using Heron's formula.

Definition at line 351 of file V_QuadMetric.cpp.

References b1, edges, VerdictVector::length(), make_quad_edges(), mb, VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{

    VerdictVector edges[4];
    make_quad_edges( edges, coordinates );

    double a1 = edges[0].length();
    double b1 = edges[1].length();
    double c1 = edges[2].length();
    double d1 = edges[3].length();

    double ma = a1 > b1 ? a1 : b1;
    double mb = c1 > d1 ? c1 : d1;
    double hm = ma > mb ? ma : mb;

    VerdictVector ab   = edges[0] * edges[1];
    VerdictVector cd   = edges[2] * edges[3];
    double denominator = ab.length() + cd.length();

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

    double aspect_ratio = .5 * hm * ( a1 + b1 + c1 + d1 ) / denominator;

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

Calculates quad condition number metric.

Maximum condition number of the Jacobian matrix at 4 corners. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

the condition of a quad

maximum condition number of the Jacobian matrix at 4 corners

Definition at line 825 of file V_QuadMetric.cpp.

References is_collapsed_quad(), VerdictVector::set(), signed_corner_areas(), v_tri_condition(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, VERDICT_MIN, and VERDICT_TRUE.

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

{

    if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_condition( 3, coordinates );

    double areas[4];
    signed_corner_areas( areas, coordinates );

    double max_condition = 0.;

    VerdictVector xxi, xet;

    double condition;

    for( int i = 0; i < 4; i++ )
    {

        xxi.set( coordinates[i][0] - coordinates[( i + 1 ) % 4][0], coordinates[i][1] - coordinates[( i + 1 ) % 4][1],
                 coordinates[i][2] - coordinates[( i + 1 ) % 4][2] );

        xet.set( coordinates[i][0] - coordinates[( i + 3 ) % 4][0], coordinates[i][1] - coordinates[( i + 3 ) % 4][1],
                 coordinates[i][2] - coordinates[( i + 3 ) % 4][2] );

        if( areas[i] < VERDICT_DBL_MIN )
            condition = VERDICT_DBL_MAX;
        else
            condition = ( xxi % xxi + xet % xet ) / areas[i];

        max_condition = VERDICT_MAX( max_condition, condition );
    }

    max_condition /= 2;

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

Calculates quad distortion metric.

{min(|J|)/actual area}*parent area, parent area = 4 for quad. Reference --- SDRC/IDEAS Simulation: Finite Element Modeling--User's Guide

the distortion of a quad

Definition at line 1049 of file V_QuadMetric.cpp.

References GaussIntegration::calculate_derivative_at_nodes(), GaussIntegration::calculate_shape_function_2d_quad(), dot_product(), moab::GeomUtil::first(), GaussIntegration::get_shape_func(), GaussIntegration::initialize(), is_collapsed_quad(), length(), VerdictVector::length(), maxNumberNodes, maxTotalNumberGaussPoints, VerdictVector::normalize(), quad_normal(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_MIN, and VERDICT_TRUE.

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

{
    // To calculate distortion for linear and 2nd order quads
    // distortion = {min(|J|)/actual area}*{parent area}
    // parent area = 4 for a quad.
    // min |J| is the minimum over nodes and gaussian integration points
    // created by Ling Pan, CAT on 4/30/01

    double element_area = 0.0, distrt, thickness_gauss;
    double cur_jacobian = 0., sign_jacobian, jacobian;
    VerdictVector aa, bb, cc, normal_at_point, xin;

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

    int total_number_of_gauss_points = number_of_gauss_points * number_of_gauss_points;

    VerdictVector face_normal = quad_normal( coordinates );

    double distortion = VERDICT_DBL_MAX;

    VerdictVector first, second;

    int i;
    // Will work out the case for collapsed quad later
    if( is_collapsed_quad( coordinates ) == VERDICT_TRUE )
    {
        for( i = 0; i < 3; i++ )
        {

            first.set( coordinates[i][0] - coordinates[( i + 1 ) % 3][0],
                       coordinates[i][1] - coordinates[( i + 1 ) % 3][1],
                       coordinates[i][2] - coordinates[( i + 1 ) % 3][2] );

            second.set( coordinates[i][0] - coordinates[( i + 2 ) % 3][0],
                        coordinates[i][1] - coordinates[( i + 2 ) % 3][1],
                        coordinates[i][2] - coordinates[( i + 2 ) % 3][2] );

            sign_jacobian = ( face_normal % ( first * second ) ) > 0 ? 1. : -1.;
            cur_jacobian  = sign_jacobian * ( first * second ).length();
            distortion    = VERDICT_MIN( distortion, cur_jacobian );
        }
        element_area = ( first * second ).length() / 2.0;
        distortion /= element_area;
    }
    else
    {
        double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
        double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
        double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
        double weight[maxTotalNumberGaussPoints];

        // create an object of GaussIntegration
        GaussIntegration::initialize( number_of_gauss_points, num_nodes );
        GaussIntegration::calculate_shape_function_2d_quad();
        GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );

        // calculate element area
        int ife, ja;
        for( ife = 0; ife < total_number_of_gauss_points; ife++ )
        {
            aa.set( 0.0, 0.0, 0.0 );
            bb.set( 0.0, 0.0, 0.0 );

            for( ja = 0; ja < num_nodes; ja++ )
            {
                xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
                aa += dndy1[ife][ja] * xin;
                bb += dndy2[ife][ja] * xin;
            }
            normal_at_point = aa * bb;
            jacobian        = normal_at_point.length();
            element_area += weight[ife] * jacobian;
        }

        double dndy1_at_node[maxNumberNodes][maxNumberNodes];
        double dndy2_at_node[maxNumberNodes][maxNumberNodes];

        GaussIntegration::calculate_derivative_at_nodes( dndy1_at_node, dndy2_at_node );

        VerdictVector normal_at_nodes[9];

        // evaluate normal at nodes and distortion values at nodes
        int jai;
        for( ja = 0; ja < num_nodes; ja++ )
        {
            aa.set( 0.0, 0.0, 0.0 );
            bb.set( 0.0, 0.0, 0.0 );
            for( jai = 0; jai < num_nodes; jai++ )
            {
                xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
                aa += dndy1_at_node[ja][jai] * xin;
                bb += dndy2_at_node[ja][jai] * xin;
            }
            normal_at_nodes[ja] = aa * bb;
            normal_at_nodes[ja].normalize();
        }

        // determine if element is flat
        bool flat_element = true;
        double dot_product;

        for( ja = 0; ja < num_nodes; ja++ )
        {
            dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
            if( fabs( dot_product ) < 0.99 )
            {
                flat_element = false;
                break;
            }
        }

        // take into consideration of the thickness of the element
        double thickness;
        // get_quad_thickness(face, element_area, thickness );
        thickness = 0.001 * sqrt( element_area );

        // set thickness gauss point location
        double zl = 0.5773502691896;
        if( flat_element ) zl = 0.0;

        int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
        double thickness_z;
        int igz;
        // loop on Gauss points
        for( ife = 0; ife < total_number_of_gauss_points; ife++ )
        {
            // loop on the thickness direction gauss points
            for( igz = 0; igz < no_gauss_pts_z; igz++ )
            {
                zl          = -zl;
                thickness_z = zl * thickness / 2.0;

                aa.set( 0.0, 0.0, 0.0 );
                bb.set( 0.0, 0.0, 0.0 );
                cc.set( 0.0, 0.0, 0.0 );

                for( ja = 0; ja < num_nodes; ja++ )
                {
                    xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
                    xin += thickness_z * normal_at_nodes[ja];
                    aa += dndy1[ife][ja] * xin;
                    bb += dndy2[ife][ja] * xin;
                    thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
                    cc += thickness_gauss * normal_at_nodes[ja];
                }

                normal_at_point = aa * bb;
                // jacobian = normal_at_point.length();
                distrt = cc % normal_at_point;
                if( distrt < distortion ) distortion = distrt;
            }
        }

        // loop through nodal points
        for( ja = 0; ja < num_nodes; ja++ )
        {
            for( igz = 0; igz < no_gauss_pts_z; igz++ )
            {
                zl          = -zl;
                thickness_z = zl * thickness / 2.0;

                aa.set( 0.0, 0.0, 0.0 );
                bb.set( 0.0, 0.0, 0.0 );
                cc.set( 0.0, 0.0, 0.0 );

                for( jai = 0; jai < num_nodes; jai++ )
                {
                    xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
                    xin += thickness_z * normal_at_nodes[ja];
                    aa += dndy1_at_node[ja][jai] * xin;
                    bb += dndy2_at_node[ja][jai] * xin;
                    if( jai == ja )
                        thickness_gauss = thickness / 2.0;
                    else
                        thickness_gauss = 0.;
                    cc += thickness_gauss * normal_at_nodes[jai];
                }
            }
            normal_at_point = aa * bb;
            sign_jacobian   = ( face_normal % normal_at_point ) > 0 ? 1. : -1.;
            distrt          = sign_jacobian * ( cc % normal_at_point );

            if( distrt < distortion ) distortion = distrt;
        }

        if( element_area * thickness != 0 )
            distortion *= 8. / ( element_area * thickness );
        else
            distortion *= 8.;
    }

    return (double)distortion;
}
C_FUNC_DEF double v_quad_edge_ratio ( int  ,
double  coordinates[][3] 
)

Calculates quad edge ratio.

edge ratio Reference --- P. P. Pebay, Planar Quadrangle Quality Measures, Eng. Comp., 2004, 20(2):157-173

the edge ratio of a quad

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

Definition at line 271 of file V_QuadMetric.cpp.

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

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

{
    VerdictVector edges[4];
    make_quad_edges( edges, coordinates );

    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 mab, Mab, mcd, Mcd, m2, M2;
    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;
    }
    m2 = mab < mcd ? mab : mcd;
    M2 = Mab > Mcd ? Mab : Mcd;

    if( m2 < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;
    else
    {
        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 );
    }
}
C_FUNC_DEF double v_quad_jacobian ( int  ,
double  coordinates[][3] 
)

Calculates quad jacobian.

Minimum pointwise volume of local map at 4 corners & center of quad. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

the jacobian of a quad

minimum pointwise volume of local map at 4 corners and center of quad

Definition at line 867 of file V_QuadMetric.cpp.

References is_collapsed_quad(), signed_corner_areas(), v_tri_area(), VERDICT_DBL_MAX, VERDICT_MAX, VERDICT_MIN, and VERDICT_TRUE.

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

{

    if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return (double)( v_tri_area( 3, coordinates ) * 2.0 );

    double areas[4];
    signed_corner_areas( areas, coordinates );

    double jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) );
    if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_quad_max_aspect_frobenius ( int  ,
double  coordinates[][3] 
)

Calculates quad maximum Frobenius aspect.

average Frobenius aspect Reference --- P. P. Pebay, Planar Quadrangle Quality Measures, Eng. Comp., 2004, 20(2):157-173

the maximum Frobenius aspect of a quad

NB (P. Pebay 01/20/07): this metric is calculated by taking the maximum of the 4 Frobenius aspects at each corner of the quad, when the reference triangle is right isosceles.

Definition at line 484 of file V_QuadMetric.cpp.

References b2, da(), edges, VerdictVector::length(), VerdictVector::length_squared(), make_quad_edges(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{

    VerdictVector edges[4];
    make_quad_edges( edges, coordinates );

    double a2 = edges[0].length_squared();
    double b2 = edges[1].length_squared();
    double c2 = edges[2].length_squared();
    double d2 = edges[3].length_squared();

    VerdictVector ab = edges[0] * edges[1];
    VerdictVector bc = edges[1] * edges[2];
    VerdictVector cd = edges[2] * edges[3];
    VerdictVector da = edges[3] * edges[0];

    double ab1 = ab.length();
    double bc1 = bc.length();
    double cd1 = cd.length();
    double da1 = da.length();

    if( ab1 < VERDICT_DBL_MIN || bc1 < VERDICT_DBL_MIN || cd1 < VERDICT_DBL_MIN || da1 < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    double qmax = ( a2 + b2 ) / ab1;

    double qcur = ( b2 + c2 ) / bc1;
    qmax        = qmax > qcur ? qmax : qcur;

    qcur = ( c2 + d2 ) / cd1;
    qmax = qmax > qcur ? qmax : qcur;

    qcur = ( d2 + a2 ) / da1;
    qmax = qmax > qcur ? qmax : qcur;

    double max_aspect_frobenius = .5 * qmax;

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

Calculates quad maximum of edge ratio.

Maximum edge length ratio at quad center. Reference --- J. Robinson, CRE Method of element testing and the Jacobian shape parameters, Eng. Comput., Vol 4, 1987.

maximum of edge ratio of a quad

maximum edge length ratio at quad center

Definition at line 321 of file V_QuadMetric.cpp.

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

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

{
    VerdictVector quad_nodes[4];
    quad_nodes[0].set( coordinates[0][0], coordinates[0][1], coordinates[0][2] );
    quad_nodes[1].set( coordinates[1][0], coordinates[1][1], coordinates[1][2] );
    quad_nodes[2].set( coordinates[2][0], coordinates[2][1], coordinates[2][2] );
    quad_nodes[3].set( coordinates[3][0], coordinates[3][1], coordinates[3][2] );

    VerdictVector principal_axes[2];
    principal_axes[0] = quad_nodes[1] + quad_nodes[2] - quad_nodes[0] - quad_nodes[3];
    principal_axes[1] = quad_nodes[2] + quad_nodes[3] - quad_nodes[0] - quad_nodes[1];

    double len1 = principal_axes[0].length();
    double len2 = principal_axes[1].length();

    if( len1 < VERDICT_DBL_MIN || len2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;

    double max_edge_ratio = VERDICT_MAX( len1 / len2, len2 / len1 );

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

Calculates quad's largest angle.

Largest included quad angle (degrees). Reference --- Unknown.

the largest angle of a quad

largest included quad area (degrees)

Definition at line 670 of file V_QuadMetric.cpp.

References moab::angle(), edges, is_collapsed_quad(), length(), VerdictVector::length(), VerdictVector::set(), signed_corner_areas(), v_tri_maximum_angle(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, VERDICT_MIN, VERDICT_PI, and VERDICT_TRUE.

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

{

    // if this is a collapsed quad, just pass it on to
    // the tri_largest_angle routine
    if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_maximum_angle( 3, coordinates );

    double angle;
    double max_angle = 0.0;

    VerdictVector edges[4];
    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] );

    // go around each node and calculate the angle
    // at each node
    double length[4];
    length[0] = edges[0].length();
    length[1] = edges[1].length();
    length[2] = edges[2].length();
    length[3] = edges[3].length();

    if( length[0] <= VERDICT_DBL_MIN || length[1] <= VERDICT_DBL_MIN || length[2] <= VERDICT_DBL_MIN ||
        length[3] <= VERDICT_DBL_MIN )
        return 0.0;

    angle     = acos( -( edges[0] % edges[1] ) / ( length[0] * length[1] ) );
    max_angle = VERDICT_MAX( angle, max_angle );

    angle     = acos( -( edges[1] % edges[2] ) / ( length[1] * length[2] ) );
    max_angle = VERDICT_MAX( angle, max_angle );

    angle     = acos( -( edges[2] % edges[3] ) / ( length[2] * length[3] ) );
    max_angle = VERDICT_MAX( angle, max_angle );

    angle     = acos( -( edges[3] % edges[0] ) / ( length[3] * length[0] ) );
    max_angle = VERDICT_MAX( angle, max_angle );

    max_angle = max_angle * 180.0 / VERDICT_PI;

    // if any signed areas are < 0, then you are getting the wrong angle
    double areas[4];
    signed_corner_areas( areas, coordinates );

    if( areas[0] < 0 || areas[1] < 0 || areas[2] < 0 || areas[3] < 0 ) { max_angle = 360 - max_angle; }

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

Calculates quad average Frobenius aspect.

average Frobenius aspect Reference --- P. P. Pebay, Planar Quadrangle Quality Measures, Eng. Comp., 2004, 20(2):157-173

the average Frobenius aspect of a quad

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

Definition at line 442 of file V_QuadMetric.cpp.

References b2, da(), edges, VerdictVector::length(), VerdictVector::length_squared(), make_quad_edges(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{

    VerdictVector edges[4];
    make_quad_edges( edges, coordinates );

    double a2 = edges[0].length_squared();
    double b2 = edges[1].length_squared();
    double c2 = edges[2].length_squared();
    double d2 = edges[3].length_squared();

    VerdictVector ab = edges[0] * edges[1];
    VerdictVector bc = edges[1] * edges[2];
    VerdictVector cd = edges[2] * edges[3];
    VerdictVector da = edges[3] * edges[0];

    double ab1 = ab.length();
    double bc1 = bc.length();
    double cd1 = cd.length();
    double da1 = da.length();

    if( ab1 < VERDICT_DBL_MIN || bc1 < VERDICT_DBL_MIN || cd1 < VERDICT_DBL_MIN || da1 < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;

    double qsum = ( a2 + b2 ) / ab1;
    qsum += ( b2 + c2 ) / bc1;
    qsum += ( c2 + d2 ) / cd1;
    qsum += ( d2 + a2 ) / da1;

    double med_aspect_frobenius = .125 * qsum;

    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_quad_minimum_angle ( int  ,
double  coordinates[][3] 
)

Calculates quad's smallest angle.

Smallest included quad angle (degrees). Reference --- Unknown.

the smallest angle of a quad

smallest included quad angle (degrees)

Definition at line 731 of file V_QuadMetric.cpp.

References moab::angle(), edges, is_collapsed_quad(), length(), VerdictVector::length(), VerdictVector::set(), v_tri_minimum_angle(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, VERDICT_MIN, VERDICT_PI, and VERDICT_TRUE.

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

{
    // if this quad is a collapsed quad, then just
    // send it to the tri_smallest_angle routine
    if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_minimum_angle( 3, coordinates );

    double angle;
    double min_angle = 360.0;

    VerdictVector edges[4];
    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] );

    // go around each node and calculate the angle
    // at each node
    double length[4];
    length[0] = edges[0].length();
    length[1] = edges[1].length();
    length[2] = edges[2].length();
    length[3] = edges[3].length();

    if( length[0] <= VERDICT_DBL_MIN || length[1] <= VERDICT_DBL_MIN || length[2] <= VERDICT_DBL_MIN ||
        length[3] <= VERDICT_DBL_MIN )
        return 360.0;

    angle     = acos( -( edges[0] % edges[1] ) / ( length[0] * length[1] ) );
    min_angle = VERDICT_MIN( angle, min_angle );

    angle     = acos( -( edges[1] % edges[2] ) / ( length[1] * length[2] ) );
    min_angle = VERDICT_MIN( angle, min_angle );

    angle     = acos( -( edges[2] % edges[3] ) / ( length[2] * length[3] ) );
    min_angle = VERDICT_MIN( angle, min_angle );

    angle     = acos( -( edges[3] % edges[0] ) / ( length[3] * length[0] ) );
    min_angle = VERDICT_MIN( angle, min_angle );

    min_angle = min_angle * 180.0 / VERDICT_PI;

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

Calculates quad oddy metric.

the oddy of a quad

general distortion measure based on left Cauchy-Green Tensor

Definition at line 785 of file V_QuadMetric.cpp.

References moab::GeomUtil::first(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{

    double max_oddy = 0.;

    VerdictVector first, second, node_pos[4];

    double g, g11, g12, g22, cur_oddy;
    int i;

    for( i = 0; i < 4; i++ )
        node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );

    for( i = 0; i < 4; i++ )
    {
        first  = node_pos[i] - node_pos[( i + 1 ) % 4];
        second = node_pos[i] - node_pos[( i + 3 ) % 4];

        g11 = first % first;
        g12 = first % second;
        g22 = second % second;
        g   = g11 * g22 - g12 * g12;

        if( g < VERDICT_DBL_MIN )
            cur_oddy = VERDICT_DBL_MAX;
        else
            cur_oddy = ( ( g11 - g22 ) * ( g11 - g22 ) + 4. * g12 * g12 ) / 2. / g;

        max_oddy = VERDICT_MAX( max_oddy, cur_oddy );
    }

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

Calculates quality metrics for quadrilateral elements.

multiple quality measures of a quad

Definition at line 1255 of file V_QuadMetric.cpp.

References QuadMetricVals::area, QuadMetricVals::aspect_ratio, QuadMetricVals::condition, determinant(), QuadMetricVals::distortion, QuadMetricVals::edge_ratio, edges, get_weight(), is_collapsed_quad(), QuadMetricVals::jacobian, VerdictVector::length(), length_squared(), VerdictVector::length_squared(), make_quad_edges(), QuadMetricVals::max_aspect_frobenius, QuadMetricVals::max_edge_ratio, QuadMetricVals::maximum_angle, QuadMetricVals::med_aspect_frobenius, QuadMetricVals::minimum_angle, normalize(), QuadMetricVals::oddy, QuadMetricVals::radius_ratio, QuadMetricVals::relative_size_squared, QuadMetricVals::scaled_jacobian, VerdictVector::set(), QuadMetricVals::shape, QuadMetricVals::shape_and_size, QuadMetricVals::shear, QuadMetricVals::shear_and_size, signed_corner_areas(), QuadMetricVals::skew, QuadMetricVals::stretch, QuadMetricVals::taper, V_QUAD_AREA, v_quad_area(), V_QUAD_ASPECT_RATIO, v_quad_aspect_ratio(), V_QUAD_CONDITION, V_QUAD_DISTORTION, v_quad_distortion(), V_QUAD_EDGE_RATIO, v_quad_edge_ratio(), V_QUAD_JACOBIAN, V_QUAD_MAX_ASPECT_FROBENIUS, v_quad_max_aspect_frobenius(), V_QUAD_MAX_EDGE_RATIO, V_QUAD_MAXIMUM_ANGLE, V_QUAD_MED_ASPECT_FROBENIUS, v_quad_med_aspect_frobenius(), V_QUAD_MINIMUM_ANGLE, V_QUAD_ODDY, V_QUAD_RADIUS_RATIO, v_quad_radius_ratio(), V_QUAD_RELATIVE_SIZE_SQUARED, V_QUAD_SCALED_JACOBIAN, V_QUAD_SHAPE, V_QUAD_SHAPE_AND_SIZE, V_QUAD_SHEAR, V_QUAD_SHEAR_AND_SIZE, V_QUAD_SKEW, V_QUAD_STRETCH, V_QUAD_TAPER, V_QUAD_WARPAGE, v_set_quad_size(), v_tri_area(), v_tri_maximum_angle(), v_tri_minimum_angle(), v_tri_scaled_jacobian(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_FALSE, VERDICT_MAX, VERDICT_MIN, VERDICT_PI, VERDICT_TRUE, and QuadMetricVals::warpage.

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

{

    memset( metric_vals, 0, sizeof( QuadMetricVals ) );

    // for starts, lets set up some basic and common information

    /*  node numbers and side numbers used below

                    2
              3 +--------- 2
               /         +
              /          |
           3 /           | 1
            /            |
           +             |
         0 -------------+ 1
               0
    */

    // vectors for each side
    VerdictVector edges[4];
    make_quad_edges( edges, coordinates );

    double areas[4];
    signed_corner_areas( areas, coordinates );

    double lengths[4];
    lengths[0] = edges[0].length();
    lengths[1] = edges[1].length();
    lengths[2] = edges[2].length();
    lengths[3] = edges[3].length();

    VerdictBoolean is_collapsed = is_collapsed_quad( coordinates );

    // handle collapsed quads metrics here
    if( is_collapsed == VERDICT_TRUE && metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE |
                                                                 V_QUAD_JACOBIAN | V_QUAD_SCALED_JACOBIAN ) )
    {
        if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
            metric_vals->minimum_angle = v_tri_minimum_angle( 3, coordinates );
        if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
            metric_vals->maximum_angle = v_tri_maximum_angle( 3, coordinates );
        if( metrics_request_flag & V_QUAD_JACOBIAN )
            metric_vals->jacobian = (double)( v_tri_area( 3, coordinates ) * 2.0 );
        if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
            metric_vals->jacobian = (double)( v_tri_scaled_jacobian( 3, coordinates ) * 2.0 );
    }

    // calculate both largest and smallest angles
    if( metrics_request_flag & ( V_QUAD_MINIMUM_ANGLE | V_QUAD_MAXIMUM_ANGLE ) && is_collapsed == VERDICT_FALSE )
    {
        // gather the angles
        double angles[4];
        angles[0] = acos( -( edges[0] % edges[1] ) / ( lengths[0] * lengths[1] ) );
        angles[1] = acos( -( edges[1] % edges[2] ) / ( lengths[1] * lengths[2] ) );
        angles[2] = acos( -( edges[2] % edges[3] ) / ( lengths[2] * lengths[3] ) );
        angles[3] = acos( -( edges[3] % edges[0] ) / ( lengths[3] * lengths[0] ) );

        if( lengths[0] <= VERDICT_DBL_MIN || lengths[1] <= VERDICT_DBL_MIN || lengths[2] <= VERDICT_DBL_MIN ||
            lengths[3] <= VERDICT_DBL_MIN )
        {
            metric_vals->minimum_angle = 360.0;
            metric_vals->maximum_angle = 0.0;
        }
        else
        {
            // if smallest angle, find the smallest angle
            if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
            {
                metric_vals->minimum_angle = VERDICT_DBL_MAX;
                for( int i = 0; i < 4; i++ )
                    metric_vals->minimum_angle = VERDICT_MIN( angles[i], metric_vals->minimum_angle );
                metric_vals->minimum_angle *= 180.0 / VERDICT_PI;
            }
            // if largest angle, find the largest angle
            if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
            {
                metric_vals->maximum_angle = 0.0;
                for( int i = 0; i < 4; i++ )
                    metric_vals->maximum_angle = VERDICT_MAX( angles[i], metric_vals->maximum_angle );
                metric_vals->maximum_angle *= 180.0 / VERDICT_PI;

                if( areas[0] < 0 || areas[1] < 0 || areas[2] < 0 || areas[3] < 0 )
                    metric_vals->maximum_angle = 360 - metric_vals->maximum_angle;
            }
        }
    }

    // handle max_edge_ratio, skew, taper, and area together
    if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
    {
        // get principle axes
        VerdictVector principal_axes[2];
        principal_axes[0] = edges[0] - edges[2];
        principal_axes[1] = edges[1] - edges[3];

        if( metrics_request_flag & ( V_QUAD_MAX_EDGE_RATIO | V_QUAD_SKEW | V_QUAD_TAPER ) )
        {
            double len1 = principal_axes[0].length();
            double len2 = principal_axes[1].length();

            // calculate the max_edge_ratio ratio
            if( metrics_request_flag & V_QUAD_MAX_EDGE_RATIO )
            {
                if( len1 < VERDICT_DBL_MIN || len2 < VERDICT_DBL_MIN )
                    metric_vals->max_edge_ratio = VERDICT_DBL_MAX;
                else
                    metric_vals->max_edge_ratio = VERDICT_MAX( len1 / len2, len2 / len1 );
            }

            // calculate the taper
            if( metrics_request_flag & V_QUAD_TAPER )
            {
                double min_length = VERDICT_MIN( len1, len2 );

                VerdictVector cross_derivative = edges[1] + edges[3];

                if( min_length < VERDICT_DBL_MIN )
                    metric_vals->taper = VERDICT_DBL_MAX;
                else
                    metric_vals->taper = cross_derivative.length() / min_length;
            }

            // calculate the skew
            if( metrics_request_flag & V_QUAD_SKEW )
            {
                if( principal_axes[0].normalize() < VERDICT_DBL_MIN || principal_axes[1].normalize() < VERDICT_DBL_MIN )
                    metric_vals->skew = 0.0;
                else
                    metric_vals->skew = fabs( principal_axes[0] % principal_axes[1] );
            }
        }
    }

    // calculate the area
    if( metrics_request_flag & ( V_QUAD_AREA | V_QUAD_RELATIVE_SIZE_SQUARED ) )
    { metric_vals->area = 0.25 * ( areas[0] + areas[1] + areas[2] + areas[3] ); }

    // calculate the relative size
    if( metrics_request_flag & ( V_QUAD_RELATIVE_SIZE_SQUARED | V_QUAD_SHAPE_AND_SIZE | V_QUAD_SHEAR_AND_SIZE ) )
    {
        double quad_area = v_quad_area( 4, coordinates );
        v_set_quad_size( quad_area );
        double w11, w21, w12, w22;
        get_weight( w11, w21, w12, w22 );
        double avg_area = determinant( w11, w21, w12, w22 );

        if( avg_area < VERDICT_DBL_MIN )
            metric_vals->relative_size_squared = 0.0;
        else
            metric_vals->relative_size_squared =
                pow( VERDICT_MIN( metric_vals->area / avg_area, avg_area / metric_vals->area ), 2 );
    }

    // calculate the jacobian
    if( metrics_request_flag & V_QUAD_JACOBIAN )
    { metric_vals->jacobian = VERDICT_MIN( VERDICT_MIN( areas[0], areas[1] ), VERDICT_MIN( areas[2], areas[3] ) ); }

    if( metrics_request_flag & ( V_QUAD_SCALED_JACOBIAN | V_QUAD_SHEAR | V_QUAD_SHEAR_AND_SIZE ) )
    {
        double scaled_jac, min_scaled_jac = VERDICT_DBL_MAX;

        if( lengths[0] < VERDICT_DBL_MIN || lengths[1] < VERDICT_DBL_MIN || lengths[2] < VERDICT_DBL_MIN ||
            lengths[3] < VERDICT_DBL_MIN )
        {
            metric_vals->scaled_jacobian = 0.0;
            metric_vals->shear           = 0.0;
        }
        else
        {
            scaled_jac     = areas[0] / ( lengths[0] * lengths[3] );
            min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

            scaled_jac     = areas[1] / ( lengths[1] * lengths[0] );
            min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

            scaled_jac     = areas[2] / ( lengths[2] * lengths[1] );
            min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

            scaled_jac     = areas[3] / ( lengths[3] * lengths[2] );
            min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

            metric_vals->scaled_jacobian = min_scaled_jac;

            // what the heck...set shear as well
            if( min_scaled_jac <= VERDICT_DBL_MIN )
                metric_vals->shear = 0.0;
            else
                metric_vals->shear = min_scaled_jac;
        }
    }

    if( metrics_request_flag & ( V_QUAD_WARPAGE | V_QUAD_ODDY ) )
    {
        VerdictVector corner_normals[4];
        corner_normals[0] = edges[3] * edges[0];
        corner_normals[1] = edges[0] * edges[1];
        corner_normals[2] = edges[1] * edges[2];
        corner_normals[3] = edges[2] * edges[3];

        if( metrics_request_flag & V_QUAD_ODDY )
        {
            double oddy, max_oddy = 0.0;

            double diff, dot_prod;

            double length_squared[4];
            length_squared[0] = corner_normals[0].length_squared();
            length_squared[1] = corner_normals[1].length_squared();
            length_squared[2] = corner_normals[2].length_squared();
            length_squared[3] = corner_normals[3].length_squared();

            if( length_squared[0] < VERDICT_DBL_MIN || length_squared[1] < VERDICT_DBL_MIN ||
                length_squared[2] < VERDICT_DBL_MIN || length_squared[3] < VERDICT_DBL_MIN )
                metric_vals->oddy = VERDICT_DBL_MAX;
            else
            {
                diff     = ( lengths[0] * lengths[0] ) - ( lengths[1] * lengths[1] );
                dot_prod = edges[0] % edges[1];
                oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[1] );
                max_oddy = VERDICT_MAX( oddy, max_oddy );

                diff     = ( lengths[1] * lengths[1] ) - ( lengths[2] * lengths[2] );
                dot_prod = edges[1] % edges[2];
                oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[2] );
                max_oddy = VERDICT_MAX( oddy, max_oddy );

                diff     = ( lengths[2] * lengths[2] ) - ( lengths[3] * lengths[3] );
                dot_prod = edges[2] % edges[3];
                oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[3] );
                max_oddy = VERDICT_MAX( oddy, max_oddy );

                diff     = ( lengths[3] * lengths[3] ) - ( lengths[0] * lengths[0] );
                dot_prod = edges[3] % edges[0];
                oddy     = ( ( diff * diff ) + 4 * dot_prod * dot_prod ) / ( 2 * length_squared[0] );
                max_oddy = VERDICT_MAX( oddy, max_oddy );

                metric_vals->oddy = max_oddy;
            }
        }

        if( metrics_request_flag & V_QUAD_WARPAGE )
        {
            if( corner_normals[0].normalize() < VERDICT_DBL_MIN || corner_normals[1].normalize() < VERDICT_DBL_MIN ||
                corner_normals[2].normalize() < VERDICT_DBL_MIN || corner_normals[3].normalize() < VERDICT_DBL_MIN )
                metric_vals->warpage = VERDICT_DBL_MAX;
            else
            {
                metric_vals->warpage =
                    pow( VERDICT_MIN( corner_normals[0] % corner_normals[2], corner_normals[1] % corner_normals[3] ),
                         3 );
            }
        }
    }

    if( metrics_request_flag & V_QUAD_STRETCH )
    {
        VerdictVector temp;

        temp.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
                  coordinates[2][2] - coordinates[0][2] );
        double diag02 = temp.length_squared();

        temp.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
                  coordinates[3][2] - coordinates[1][2] );
        double diag13 = temp.length_squared();

        static const double QUAD_STRETCH_FACTOR = sqrt( 2.0 );

        // 'diag02' is now the max diagonal of the quad
        diag02 = VERDICT_MAX( diag02, diag13 );

        if( diag02 < VERDICT_DBL_MIN )
            metric_vals->stretch = VERDICT_DBL_MAX;
        else
            metric_vals->stretch =
                QUAD_STRETCH_FACTOR *
                VERDICT_MIN( VERDICT_MIN( lengths[0], lengths[1] ), VERDICT_MIN( lengths[2], lengths[3] ) ) /
                sqrt( diag02 );
    }

    if( metrics_request_flag & ( V_QUAD_CONDITION | V_QUAD_SHAPE | V_QUAD_SHAPE_AND_SIZE ) )
    {
        double lengths_squared[4];
        lengths_squared[0] = edges[0].length_squared();
        lengths_squared[1] = edges[1].length_squared();
        lengths_squared[2] = edges[2].length_squared();
        lengths_squared[3] = edges[3].length_squared();

        if( areas[0] < VERDICT_DBL_MIN || areas[1] < VERDICT_DBL_MIN || areas[2] < VERDICT_DBL_MIN ||
            areas[3] < VERDICT_DBL_MIN )
        {
            metric_vals->condition = VERDICT_DBL_MAX;
            metric_vals->shape     = 0.0;
        }
        else
        {
            double max_condition = 0.0, condition;

            condition     = ( lengths_squared[0] + lengths_squared[3] ) / areas[0];
            max_condition = VERDICT_MAX( max_condition, condition );

            condition     = ( lengths_squared[1] + lengths_squared[0] ) / areas[1];
            max_condition = VERDICT_MAX( max_condition, condition );

            condition     = ( lengths_squared[2] + lengths_squared[1] ) / areas[2];
            max_condition = VERDICT_MAX( max_condition, condition );

            condition     = ( lengths_squared[3] + lengths_squared[2] ) / areas[3];
            max_condition = VERDICT_MAX( max_condition, condition );

            metric_vals->condition = 0.5 * max_condition;
            metric_vals->shape     = 2 / max_condition;
        }
    }

    if( metrics_request_flag & V_QUAD_AREA )
    {
        if( metric_vals->area > 0 ) metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
        metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_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 );
        metric_vals->max_edge_ratio = (double)VERDICT_MAX( metric_vals->max_edge_ratio, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_CONDITION )
    {
        if( metric_vals->condition > 0 )
            metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
        metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
    }

    // calculate distortion
    if( metrics_request_flag & V_QUAD_DISTORTION )
    {
        metric_vals->distortion = v_quad_distortion( num_nodes, coordinates );

        if( metric_vals->distortion > 0 )
            metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
        metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_JACOBIAN )
    {
        if( metric_vals->jacobian > 0 )
            metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
        metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_MAXIMUM_ANGLE )
    {
        if( metric_vals->maximum_angle > 0 )
            metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
        metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_MINIMUM_ANGLE )
    {
        if( metric_vals->minimum_angle > 0 )
            metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
        metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_ODDY )
    {
        if( metric_vals->oddy > 0 ) metric_vals->oddy = (double)VERDICT_MIN( metric_vals->oddy, VERDICT_DBL_MAX );
        metric_vals->oddy = (double)VERDICT_MAX( metric_vals->oddy, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_RELATIVE_SIZE_SQUARED )
    {
        if( metric_vals->relative_size_squared > 0 )
            metric_vals->relative_size_squared =
                (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
        metric_vals->relative_size_squared =
            (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_SCALED_JACOBIAN )
    {
        if( metric_vals->scaled_jacobian > 0 )
            metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
        metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_SHEAR )
    {
        if( metric_vals->shear > 0 ) metric_vals->shear = (double)VERDICT_MIN( metric_vals->shear, VERDICT_DBL_MAX );
        metric_vals->shear = (double)VERDICT_MAX( metric_vals->shear, -VERDICT_DBL_MAX );
    }

    // calculate shear and size
    // reuse values from above
    if( metrics_request_flag & V_QUAD_SHEAR_AND_SIZE )
    {
        metric_vals->shear_and_size = metric_vals->shear * metric_vals->relative_size_squared;

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

    if( metrics_request_flag & V_QUAD_SHAPE )
    {
        if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
        metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
    }

    // calculate shape and size
    // reuse values from above
    if( metrics_request_flag & V_QUAD_SHAPE_AND_SIZE )
    {
        metric_vals->shape_and_size = metric_vals->shape * metric_vals->relative_size_squared;

        if( metric_vals->shape_and_size > 0 )
            metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
        metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_SKEW )
    {
        if( metric_vals->skew > 0 ) metric_vals->skew = (double)VERDICT_MIN( metric_vals->skew, VERDICT_DBL_MAX );
        metric_vals->skew = (double)VERDICT_MAX( metric_vals->skew, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_STRETCH )
    {
        if( metric_vals->stretch > 0 )
            metric_vals->stretch = (double)VERDICT_MIN( metric_vals->stretch, VERDICT_DBL_MAX );
        metric_vals->stretch = (double)VERDICT_MAX( metric_vals->stretch, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_TAPER )
    {
        if( metric_vals->taper > 0 ) metric_vals->taper = (double)VERDICT_MIN( metric_vals->taper, VERDICT_DBL_MAX );
        metric_vals->taper = (double)VERDICT_MAX( metric_vals->taper, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_WARPAGE )
    {
        if( metric_vals->warpage > 0 )
            metric_vals->warpage = (double)VERDICT_MIN( metric_vals->warpage, VERDICT_DBL_MAX );
        metric_vals->warpage = (double)VERDICT_MAX( metric_vals->warpage, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_QUAD_MED_ASPECT_FROBENIUS )
        metric_vals->med_aspect_frobenius = v_quad_med_aspect_frobenius( 4, coordinates );
    if( metrics_request_flag & V_QUAD_MAX_ASPECT_FROBENIUS )
        metric_vals->max_aspect_frobenius = v_quad_max_aspect_frobenius( 4, coordinates );
    if( metrics_request_flag & V_QUAD_EDGE_RATIO ) metric_vals->edge_ratio = v_quad_edge_ratio( 4, coordinates );
    if( metrics_request_flag & V_QUAD_ASPECT_RATIO ) metric_vals->aspect_ratio = v_quad_aspect_ratio( 4, coordinates );
    if( metrics_request_flag & V_QUAD_RADIUS_RATIO ) metric_vals->radius_ratio = v_quad_radius_ratio( 4, coordinates );
}
C_FUNC_DEF double v_quad_radius_ratio ( int  ,
double  coordinates[][3] 
)

Calculates quad radius ratio.

radius ratio Reference --- P. P. Pebay, Planar Quadrangle Quality Measures, Eng. Comp., 2004, 20(2):157-173

the radius ratio of a quad

NB (P. Pebay 01/19/07): this metric is called "radius ratio" by extension of a concept that does not exist in general with quads -- although a different name should probably be used in the future.

Definition at line 386 of file V_QuadMetric.cpp.

References b2, da(), edges, VerdictVector::length(), VerdictVector::length_squared(), make_quad_edges(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{
    static const double normal_coeff = 1. / ( 2. * sqrt( 2. ) );

    VerdictVector edges[4];
    make_quad_edges( edges, coordinates );

    double a2 = edges[0].length_squared();
    double b2 = edges[1].length_squared();
    double c2 = edges[2].length_squared();
    double d2 = edges[3].length_squared();

    VerdictVector diag;
    diag.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
              coordinates[2][2] - coordinates[0][2] );
    double m2 = diag.length_squared();

    diag.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
              coordinates[3][2] - coordinates[1][2] );
    double n2 = diag.length_squared();

    double t0 = a2 > b2 ? a2 : b2;
    double t1 = c2 > d2 ? c2 : d2;
    double t2 = m2 > n2 ? m2 : n2;
    double h2 = t0 > t1 ? t0 : t1;
    h2        = h2 > t2 ? h2 : t2;

    VerdictVector ab = edges[0] * edges[1];
    VerdictVector bc = edges[1] * edges[2];
    VerdictVector cd = edges[2] * edges[3];
    VerdictVector da = edges[3] * edges[0];

    t0        = da.length();
    t1        = ab.length();
    t2        = bc.length();
    double t3 = cd.length();

    t0 = t0 < t1 ? t0 : t1;
    t2 = t2 < t3 ? t2 : t3;
    t0 = t0 < t2 ? t0 : t2;

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

    double radius_ratio = normal_coeff * sqrt( ( a2 + b2 + c2 + d2 ) * h2 ) / t0;

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

Calculates quad relative size metric.

Min( J, 1/J ), where J is determinant of weighted Jacobian matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

the relative size of a quad

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

Definition at line 985 of file V_QuadMetric.cpp.

References determinant(), get_weight(), v_quad_area(), v_set_quad_size(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

Referenced by moab::VerdictWrapper::quality_measure(), v_quad_shape_and_size(), and v_quad_shear_and_size().

{

    double quad_area = v_quad_area( 4, coordinates );
    double rel_size  = 0;

    v_set_quad_size( quad_area );
    double w11, w21, w12, w22;
    get_weight( w11, w21, w12, w22 );
    double avg_area = determinant( w11, w21, w12, w22 );

    if( avg_area > VERDICT_DBL_MIN )
    {

        w11 = quad_area / avg_area;

        if( w11 > VERDICT_DBL_MIN )
        {
            rel_size = VERDICT_MIN( w11, 1 / w11 );
            rel_size *= rel_size;
        }
    }

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

Calculates quad scaled jacobian.

Minimum Jacobian divided by the lengths of the 2 edge vectors. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

scaled jacobian of a quad

Minimum Jacobian divided by the lengths of the 2 edge vector

Definition at line 885 of file V_QuadMetric.cpp.

References edges, is_collapsed_quad(), length(), VerdictVector::length(), make_quad_edges(), signed_corner_areas(), v_tri_scaled_jacobian(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, VERDICT_MIN, and VERDICT_TRUE.

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

{
    if( is_collapsed_quad( coordinates ) == VERDICT_TRUE ) return v_tri_scaled_jacobian( 3, coordinates );

    double corner_areas[4], min_scaled_jac = VERDICT_DBL_MAX, scaled_jac;
    signed_corner_areas( corner_areas, coordinates );

    VerdictVector edges[4];
    make_quad_edges( edges, coordinates );

    double length[4];
    length[0] = edges[0].length();
    length[1] = edges[1].length();
    length[2] = edges[2].length();
    length[3] = edges[3].length();

    if( length[0] < VERDICT_DBL_MIN || length[1] < VERDICT_DBL_MIN || length[2] < VERDICT_DBL_MIN ||
        length[3] < VERDICT_DBL_MIN )
        return 0.0;

    scaled_jac     = corner_areas[0] / ( length[0] * length[3] );
    min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

    scaled_jac     = corner_areas[1] / ( length[1] * length[0] );
    min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

    scaled_jac     = corner_areas[2] / ( length[2] * length[1] );
    min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

    scaled_jac     = corner_areas[3] / ( length[3] * length[2] );
    min_scaled_jac = VERDICT_MIN( scaled_jac, min_scaled_jac );

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

Calculates quad shape metric.

2/Condition number of weighted Jacobian matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

the shape of a quad

2/Condition number of weighted Jacobian matrix

Definition at line 941 of file V_QuadMetric.cpp.

References edges, length_squared(), VerdictVector::length_squared(), make_quad_edges(), signed_corner_areas(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{

    double corner_areas[4], min_shape = VERDICT_DBL_MAX, shape;
    signed_corner_areas( corner_areas, coordinates );

    VerdictVector edges[4];
    make_quad_edges( edges, coordinates );

    double length_squared[4];
    length_squared[0] = edges[0].length_squared();
    length_squared[1] = edges[1].length_squared();
    length_squared[2] = edges[2].length_squared();
    length_squared[3] = edges[3].length_squared();

    if( length_squared[0] <= VERDICT_DBL_MIN || length_squared[1] <= VERDICT_DBL_MIN ||
        length_squared[2] <= VERDICT_DBL_MIN || length_squared[3] <= VERDICT_DBL_MIN )
        return 0.0;

    shape     = corner_areas[0] / ( length_squared[0] + length_squared[3] );
    min_shape = VERDICT_MIN( shape, min_shape );

    shape     = corner_areas[1] / ( length_squared[1] + length_squared[0] );
    min_shape = VERDICT_MIN( shape, min_shape );

    shape     = corner_areas[2] / ( length_squared[2] + length_squared[1] );
    min_shape = VERDICT_MIN( shape, min_shape );

    shape     = corner_areas[3] / ( length_squared[3] + length_squared[2] );
    min_shape = VERDICT_MIN( shape, min_shape );

    min_shape *= 2;

    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_quad_shape_and_size ( int  num_nodes,
double  coordinates[][3] 
)

Calculates quad shape-size metric.

Product of Shape and Relative Size. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

the relative shape and size of a quad

Product of Shape and Relative Size

Definition at line 1017 of file V_QuadMetric.cpp.

References size, v_quad_relative_size_squared(), v_quad_shape(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    double shape, size;
    size  = v_quad_relative_size_squared( num_nodes, coordinates );
    shape = v_quad_shape( num_nodes, coordinates );

    double shape_and_size = shape * size;

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

Calculates quad shear metric.

2/Condition number of Jacobian Skew matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

the shear of a quad

2/Condition number of Jacobian Skew matrix

Definition at line 926 of file V_QuadMetric.cpp.

References v_quad_scaled_jacobian(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, and VERDICT_MIN.

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

{
    double scaled_jacobian = v_quad_scaled_jacobian( 4, coordinates );

    if( scaled_jacobian <= VERDICT_DBL_MIN )
        return 0.0;
    else
        return (double)VERDICT_MIN( scaled_jacobian, VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_quad_shear_and_size ( int  num_nodes,
double  coordinates[][3] 
)

Calculates quad shear-size metric.

Product of Shear and Relative Size. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

the shear and size of a quad

product of shear and relative size

Definition at line 1034 of file V_QuadMetric.cpp.

References size, v_quad_relative_size_squared(), v_quad_shear(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    double shear, size;
    shear = v_quad_shear( num_nodes, coordinates );
    size  = v_quad_relative_size_squared( num_nodes, coordinates );

    double shear_and_size = shear * size;

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

Calculates quad skew metric.

Maximum |cos A| where A is the angle between edges at quad center. Reference --- J. Robinson, CRE Method of element testing and the Jacobian shape parameters, Eng. Comput., Vol 4, 1987.

skew of a quad

maximum ||cos A|| where A is the angle between edges at quad center

Definition at line 530 of file V_QuadMetric.cpp.

References normalize(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, and VERDICT_MIN.

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

{
    VerdictVector node_pos[4];
    for( int i = 0; i < 4; i++ )
        node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );

    VerdictVector principle_axes[2];
    principle_axes[0] = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
    principle_axes[1] = node_pos[2] + node_pos[3] - node_pos[0] - node_pos[1];

    if( principle_axes[0].normalize() < VERDICT_DBL_MIN ) return 0.0;
    if( principle_axes[1].normalize() < VERDICT_DBL_MIN ) return 0.0;

    double skew = fabs( principle_axes[0] % principle_axes[1] );

    return (double)VERDICT_MIN( skew, VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_quad_stretch ( int  ,
double  coordinates[][3] 
)

Calculates quad strech metric.

Sqrt(2) * minimum edge length / maximum diagonal length. Reference --- FIMESH code.

the stretch of a quad

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

Definition at line 628 of file V_QuadMetric.cpp.

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

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

{
    VerdictVector edges[4], temp;
    make_quad_edges( edges, coordinates );

    double lengths_squared[4];
    lengths_squared[0] = edges[0].length_squared();
    lengths_squared[1] = edges[1].length_squared();
    lengths_squared[2] = edges[2].length_squared();
    lengths_squared[3] = edges[3].length_squared();

    temp.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
              coordinates[2][2] - coordinates[0][2] );
    double diag02 = temp.length_squared();

    temp.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
              coordinates[3][2] - coordinates[1][2] );
    double diag13 = temp.length_squared();

    static const double QUAD_STRETCH_FACTOR = sqrt( 2.0 );

    // 'diag02' is now the max diagonal of the quad
    diag02 = VERDICT_MAX( diag02, diag13 );

    if( diag02 < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;
    else
    {
        double stretch =
            (double)( QUAD_STRETCH_FACTOR * sqrt( VERDICT_MIN( VERDICT_MIN( lengths_squared[0], lengths_squared[1] ),
                                                               VERDICT_MIN( lengths_squared[2], lengths_squared[3] ) ) /
                                                  diag02 ) );

        return (double)VERDICT_MIN( stretch, VERDICT_DBL_MAX );
    }
}
C_FUNC_DEF double v_quad_taper ( int  ,
double  coordinates[][3] 
)

Calculates quad taper metric.

Maximum ratio of lengths derived from opposite edges. Reference --- J. Robinson, CRE Method of element testing and the Jacobian shape parameters, Eng. Comput., Vol 4, 1987.

taper of a quad

maximum ratio of lengths derived from opposite edges

Definition at line 553 of file V_QuadMetric.cpp.

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

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

{
    VerdictVector node_pos[4];
    for( int i = 0; i < 4; i++ )
        node_pos[i].set( coordinates[i][0], coordinates[i][1], coordinates[i][2] );

    VerdictVector principle_axes[2];
    principle_axes[0] = node_pos[1] + node_pos[2] - node_pos[3] - node_pos[0];
    principle_axes[1] = node_pos[2] + node_pos[3] - node_pos[0] - node_pos[1];

    VerdictVector cross_derivative = node_pos[0] + node_pos[2] - node_pos[1] - node_pos[3];

    double lengths[2];
    lengths[0] = principle_axes[0].length();
    lengths[1] = principle_axes[1].length();

    // get min length
    lengths[0] = VERDICT_MIN( lengths[0], lengths[1] );

    if( lengths[0] < VERDICT_DBL_MIN ) return VERDICT_DBL_MAX;

    double taper = cross_derivative.length() / lengths[0];
    return (double)VERDICT_MIN( taper, VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_quad_warpage ( int  ,
double  coordinates[][3] 
)

Calculates quad warpage metric.

Cosine of Minimum Dihedral Angle formed by Planes Intersecting in Diagonals. Reference --- J. Robinson, CRE Method of element testing and the Jacobian shape parameters, Eng. Comput., Vol 4, 1987.

warpage of a quad

deviation of element from planarity

Definition at line 583 of file V_QuadMetric.cpp.

References edges, make_quad_edges(), normalize(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{

    VerdictVector edges[4];
    make_quad_edges( edges, coordinates );

    VerdictVector corner_normals[4];
    corner_normals[0] = edges[3] * edges[0];
    corner_normals[1] = edges[0] * edges[1];
    corner_normals[2] = edges[1] * edges[2];
    corner_normals[3] = edges[2] * edges[3];

    if( corner_normals[0].normalize() < VERDICT_DBL_MIN || corner_normals[1].normalize() < VERDICT_DBL_MIN ||
        corner_normals[2].normalize() < VERDICT_DBL_MIN || corner_normals[3].normalize() < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MIN;

    double warpage =
        pow( VERDICT_MIN( corner_normals[0] % corner_normals[2], corner_normals[1] % corner_normals[3] ), 3 );

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

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().

C_FUNC_DEF void v_set_quad_size ( double  size)

Sets average size (area) of quad, needed for v_quad_relative_size(...)

Definition at line 56 of file V_QuadMetric.cpp.

References size, and verdict_quad_size.

Referenced by moab::VerdictWrapper::set_size(), v_quad_quality(), and v_quad_relative_size_squared().

C_FUNC_DEF void v_set_tet_size ( double  size)

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

set the average volume of a tet

Definition at line 37 of file V_TetMetric.cpp.

References size, and verdict_tet_size.

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

Sets fuction pointer to calculate tri normal wrt surface.

Definition at line 63 of file V_TriMetric.cpp.

References compute_normal.

{
    compute_normal = func;
}
C_FUNC_DEF void v_set_tri_size ( double  size)

Sets average size (area) of tri, needed for v_tri_relative_size(...)

sets the average area of a tri

Definition at line 58 of file V_TriMetric.cpp.

References size, and verdict_tri_size.

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

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

Calculates the radius ratio metric of a positively oriented tet.

CR / (3.0 * IR) where CR = circumsphere radius, IR = inscribed sphere radius if the element is positively-oriented. Reference --- V. N. Parthasarathy et al, A comparison of tetrahedron quality measures, Finite Elem. Anal. Des., Vol 15(1993), 255-261.

The radius ratio of a positively-oriented tet, a.k.a. "aspect beta"

NB (P. Pebay 04/16/07): CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed sphere radius if the element has positive orientation. Note that this metric is similar to the radius ratio of a tet, except that it returns VERDICT_DBL_MAX if the element has negative orientation.

Definition at line 265 of file V_TetMetric.cpp.

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

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

{

    // Determine side vectors
    VerdictVector side[6];

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

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

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

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

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

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

    VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
                              side[2].length_squared() * ( side[3] * side[0] ) +
                              side[0].length_squared() * ( side[3] * side[2] );

    double area_sum;
    area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
                 ( side[3] * side[2] ).length() ) *
               0.5;

    double volume = v_tet_volume( 4, coordinates );

    if( volume < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;
    else
    {
        double radius_ratio;
        radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );

        return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
    }
}
C_FUNC_DEF double v_tet_aspect_frobenius ( int  ,
double  coordinates[][3] 
)

Calculates tet aspect frobenius metric.

Frobenius condition number when the reference element is regular Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

The aspect frobenius of a tet

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

Definition at line 426 of file V_TetMetric.cpp.

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

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

{
    static const double normal_exp = 1. / 3.;

    VerdictVector ab, ac, ad;

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

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

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

    double denominator = ab % ( ac * ad );
    denominator *= denominator;
    denominator *= 2.;
    denominator = 3. * pow( denominator, normal_exp );

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

    double u[3];
    ab.get_xyz( u );
    double v[3];
    ac.get_xyz( v );
    double w[3];
    ad.get_xyz( w );

    double numerator = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
    numerator += v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
    numerator += w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
    numerator *= 1.5;
    numerator -= v[0] * u[0] + v[1] * u[1] + v[2] * u[2];
    numerator -= w[0] * u[0] + w[1] * u[1] + w[2] * u[2];
    numerator -= w[0] * v[0] + w[1] * v[1] + w[2] * v[2];

    double aspect_frobenius = numerator / denominator;

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

Calculates tet aspect gamma metric.

Srms**3 / (8.479670*V) where Srms = sqrt(Sum(Si**2)/6), Si = edge length. Reference --- V. N. Parthasarathy et al, A comparison of tetrahedron quality measures, Finite Elem. Anal. Des., Vol 15(1993), 255-261.

the aspect gamma of a tet

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

Definition at line 381 of file V_TetMetric.cpp.

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

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

{

    // Determine side vectors
    VerdictVector side0, side1, side2, side3, side4, side5;

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

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

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

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

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

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

    double volume = fabs( v_tet_volume( 4, coordinates ) );

    if( volume < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;
    else
    {
        double srms = sqrt( ( side0.length_squared() + side1.length_squared() + side2.length_squared() +
                              side3.length_squared() + side4.length_squared() + side5.length_squared() ) /
                            6.0 );

        double aspect_ratio_gamma = pow( srms, 3 ) / ( 8.48528137423857 * volume );
        return (double)aspect_ratio_gamma;
    }
}
C_FUNC_DEF double v_tet_aspect_ratio ( int  ,
double  coordinates[][3] 
)

Calculates tet aspect ratio metric.

Hmax / (2 sqrt(6) r) where Hmax and r respectively denote the greatest edge length and the inradius of the tetrahedron Reference --- P. Frey and P.-L. George, Meshing, Hermes (2000).

The aspect ratio of a tet

NB (P. Pebay 01/22/07): Hmax / (2 sqrt(6) r) where Hmax and r respectively denote the greatest edge length and the inradius of the tetrahedron

Definition at line 318 of file V_TetMetric.cpp.

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

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

{
    static const double normal_coeff = sqrt( 6. ) / 12.;

    // Determine side vectors
    VerdictVector ab, bc, ac, ad, bd, cd;

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

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

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

    double detTet = ab % ( ac * ad );

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

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

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

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

    double ab2 = ab.length_squared();
    double bc2 = bc.length_squared();
    double ac2 = ac.length_squared();
    double ad2 = ad.length_squared();
    double bd2 = bd.length_squared();
    double cd2 = cd.length_squared();

    double A  = ab2 > bc2 ? ab2 : bc2;
    double B  = ac2 > ad2 ? ac2 : ad2;
    double C  = bd2 > cd2 ? bd2 : cd2;
    double D  = A > B ? A : B;
    double hm = D > C ? sqrt( D ) : sqrt( C );

    bd = ab * bc;
    A  = bd.length();
    bd = ab * ad;
    B  = bd.length();
    bd = ac * ad;
    C  = bd.length();
    bd = bc * cd;
    D  = bd.length();

    double aspect_ratio;
    aspect_ratio = normal_coeff * hm * ( A + B + C + D ) / fabs( detTet );

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

Calculates tet collapse ratio metric.

Collapse ratio

The collapse ratio of a tet

Definition at line 526 of file V_TetMetric.cpp.

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

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

{
    // Determine side vectors
    VerdictVector e01, e02, e03, e12, e13, e23;

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

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

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

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

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

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

    double l[6];
    l[0] = e01.length();
    l[1] = e02.length();
    l[2] = e03.length();
    l[3] = e12.length();
    l[4] = e13.length();
    l[5] = e23.length();

    // Find longest edge for each bounding triangle of tetrahedron
    double l012 = l[4] > l[0] ? l[4] : l[0];
    l012        = l[1] > l012 ? l[1] : l012;
    double l031 = l[0] > l[2] ? l[0] : l[2];
    l031        = l[3] > l031 ? l[3] : l031;
    double l023 = l[2] > l[1] ? l[2] : l[1];
    l023        = l[5] > l023 ? l[5] : l023;
    double l132 = l[4] > l[3] ? l[4] : l[3];
    l132        = l[5] > l132 ? l[5] : l132;

    // Compute collapse ratio for each vertex/triangle pair
    VerdictVector N;
    double h, magN;
    double cr;
    double crMin;

    N     = e01 * e02;
    magN  = N.length();
    h     = ( e03 % N ) / magN;  // height of vertex 3 above 0-1-2
    crMin = h / l012;            // ratio of height to longest edge of 0-1-2

    N    = e03 * e01;
    magN = N.length();
    h    = ( e02 % N ) / magN;  // height of vertex 2 above 0-3-1
    cr   = h / l031;            // ratio of height to longest edge of 0-3-1
    if( cr < crMin ) crMin = cr;

    N    = e02 * e03;
    magN = N.length();
    h    = ( e01 % N ) / magN;  // height of vertex 1 above 0-2-3
    cr   = h / l023;            // ratio of height to longest edge of 0-2-3
    if( cr < crMin ) crMin = cr;

    N    = e12 * e13;
    magN = N.length();
    h    = ( e01 % N ) / magN;  // height of vertex 0 above 1-3-2
    cr   = h / l132;            // ratio of height to longest edge of 1-3-2
    if( cr < crMin ) crMin = cr;

    if( crMin < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
    if( crMin > 0 ) return (double)VERDICT_MIN( crMin, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( crMin, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_tet_condition ( int  ,
double  coordinates[][3] 
)

Calculates tet condition metric.

Condition number of the Jacobian matrix at any corner. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

the condition of a tet

condition number of the jacobian matrix at any corner

Definition at line 629 of file V_TetMetric.cpp.

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

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

{

    double condition, term1, term2, det;
    double rt3 = sqrt( 3.0 );
    double rt6 = sqrt( 6.0 );

    VerdictVector side0, side2, side3;

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

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

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

    VerdictVector c_1, c_2, c_3;

    c_1 = side0;
    c_2 = ( -2 * side2 - side0 ) / rt3;
    c_3 = ( 3 * side3 + side2 - side0 ) / rt6;

    term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
    term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) + ( c_2 * c_3 ) % ( c_2 * c_3 ) + ( c_1 * c_3 ) % ( c_1 * c_3 );
    det   = c_1 % ( c_2 * c_3 );

    if( det <= VERDICT_DBL_MIN )
        return VERDICT_DBL_MAX;
    else
        condition = sqrt( term1 * term2 ) / ( 3.0 * det );

    return (double)condition;
}
C_FUNC_DEF double v_tet_distortion ( int  num_nodes,
double  coordinates[][3] 
)

Calculates tet distortion metric.

{min(|J|)/actual volume}*parent volume, parent volume = 1/6 for tet. Reference --- SDRC/IDEAS Simulation: Finite Element Modeling--User's Guide

the distortion of a tet

Definition at line 762 of file V_TetMetric.cpp.

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

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

{

    double distortion          = VERDICT_DBL_MAX;
    int number_of_gauss_points = 0;
    if( num_nodes == 4 )
        // for linear tet, the distortion is always 1 because
        // straight edge tets are the target shape for tet
        return 1.0;

    else if( num_nodes == 10 )
        // use four integration points for quadratic tet
        number_of_gauss_points = 4;

    int number_dims                  = 3;
    int total_number_of_gauss_points = number_of_gauss_points;
    // use is_tri=1 to indicate this is for tet in 3D
    int is_tri = 1;

    double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
    double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
    double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
    double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
    double weight[maxTotalNumberGaussPoints];

    // create an object of GaussIntegration for tet
    GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
    GaussIntegration::calculate_shape_function_3d_tet();
    GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight );

    // vector xxi is the derivative vector of coordinates w.r.t local xi coordinate in the
    // computation space
    // vector xet is the derivative vector of coordinates w.r.t local et coordinate in the
    // computation space
    // vector xze is the derivative vector of coordinates w.r.t local ze coordinate in the
    // computation space
    VerdictVector xxi, xet, xze, xin;

    double jacobian, minimum_jacobian;
    double element_volume = 0.0;
    minimum_jacobian      = VERDICT_DBL_MAX;

    // calculate element volume
    int ife, ja;
    for( ife = 0; ife < total_number_of_gauss_points; ife++ )
    {
        xxi.set( 0.0, 0.0, 0.0 );
        xet.set( 0.0, 0.0, 0.0 );
        xze.set( 0.0, 0.0, 0.0 );

        for( ja = 0; ja < num_nodes; ja++ )
        {
            xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
            xxi += dndy1[ife][ja] * xin;
            xet += dndy2[ife][ja] * xin;
            xze += dndy3[ife][ja] * xin;
        }

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

        element_volume += weight[ife] * jacobian;
    }  // element_volume is 6 times the actual volume

    // loop through all nodes
    double dndy1_at_node[maxNumberNodes][maxNumberNodes];
    double dndy2_at_node[maxNumberNodes][maxNumberNodes];
    double dndy3_at_node[maxNumberNodes][maxNumberNodes];

    GaussIntegration::calculate_derivative_at_nodes_3d_tet( dndy1_at_node, dndy2_at_node, dndy3_at_node );
    int node_id;
    for( node_id = 0; node_id < num_nodes; node_id++ )
    {
        xxi.set( 0.0, 0.0, 0.0 );
        xet.set( 0.0, 0.0, 0.0 );
        xze.set( 0.0, 0.0, 0.0 );

        for( ja = 0; ja < num_nodes; ja++ )
        {
            xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
            xxi += dndy1_at_node[node_id][ja] * xin;
            xet += dndy2_at_node[node_id][ja] * xin;
            xze += dndy3_at_node[node_id][ja] * xin;
        }

        jacobian = xxi % ( xet * xze );
        if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
    }
    distortion = minimum_jacobian / element_volume;

    return (double)distortion;
}
C_FUNC_DEF double v_tet_edge_ratio ( int  ,
double  coordinates[][3] 
)

Calculates tet edge ratio metric.

Hmax / Hmin where Hmax and Hmin are respectively the maximum and the minimum edge lengths

the edge ratio of a tet

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

Definition at line 71 of file V_TetMetric.cpp.

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

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

{
    VerdictVector a, b, c, d, e, f;

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

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

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

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

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

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

    double a2 = a.length_squared();
    double b2 = b.length_squared();
    double c2 = c.length_squared();
    double d2 = d.length_squared();
    double e2 = e.length_squared();
    double f2 = f.length_squared();

    double m2, M2, mab, mcd, mef, Mab, Mcd, Mef;

    if( a2 < b2 )
    {
        mab = a2;
        Mab = b2;
    }
    else  // b2 <= a2
    {
        mab = b2;
        Mab = a2;
    }
    if( c2 < d2 )
    {
        mcd = c2;
        Mcd = d2;
    }
    else  // d2 <= c2
    {
        mcd = d2;
        Mcd = c2;
    }
    if( e2 < f2 )
    {
        mef = e2;
        Mef = f2;
    }
    else  // f2 <= e2
    {
        mef = f2;
        Mef = e2;
    }

    m2 = mab < mcd ? mab : mcd;
    m2 = m2 < mef ? m2 : mef;

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

    M2 = Mab > Mcd ? Mab : Mcd;
    M2 = M2 > Mef ? M2 : Mef;

    double edge_ratio = sqrt( M2 / m2 );

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

Calculates tet jacobian.

Minimum pointwise volume at any corner. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

the jacobian of a tet

TODO

Definition at line 670 of file V_TetMetric.cpp.

References VerdictVector::set().

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

{
    VerdictVector side0, side2, side3;

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

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

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

    return (double)( side3 % ( side2 * side0 ) );
}
C_FUNC_DEF double v_tet_minimum_angle ( int  ,
double  coordinates[][3] 
)

Calculates tet minimum dihedral angle.

Minimum (nonoriented) dihedral angle of a tetrahedron, expressed in degrees.

The minimum angle of a tet

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

Definition at line 475 of file V_TetMetric.cpp.

References beta, epsilon, VerdictVector::length(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, VERDICT_MIN, and MBMesquite::zeta.

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

{
    static const double normal_coeff = 180. * .3183098861837906715377675267450287;

    // Determine side vectors
    VerdictVector ab, bc, ad, cd;

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

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

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

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

    VerdictVector abc = ab * bc;
    double nabc       = abc.length();
    VerdictVector abd = ab * ad;
    double nabd       = abd.length();
    VerdictVector acd = ad * cd;
    double nacd       = acd.length();
    VerdictVector bcd = bc * cd;
    double nbcd       = bcd.length();

    double alpha   = acos( ( abc % abd ) / ( nabc * nabd ) );
    double beta    = acos( ( abc % acd ) / ( nabc * nacd ) );
    double gamma   = acos( ( abc % bcd ) / ( nabc * nbcd ) );
    double delta   = acos( ( abd % acd ) / ( nabd * nacd ) );
    double epsilon = acos( ( abd % bcd ) / ( nabd * nbcd ) );
    double zeta    = acos( ( acd % bcd ) / ( nacd * nbcd ) );

    alpha = alpha < beta ? alpha : beta;
    alpha = alpha < gamma ? alpha : gamma;
    alpha = alpha < delta ? alpha : delta;
    alpha = alpha < epsilon ? alpha : epsilon;
    alpha = alpha < zeta ? alpha : zeta;
    alpha *= normal_coeff;

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

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

Calculates quality metrics for tetrahedral elements.

the quality metrics of a tet

Definition at line 859 of file V_TetMetric.cpp.

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

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

{

    memset( metric_vals, 0, sizeof( TetMetricVals ) );

    /*

      node numbers and edge numbers below



               3
               +            edge 0 is node 0 to 1
              +|+           edge 1 is node 1 to 2
            3/ | \5         edge 2 is node 0 to 2
            / 4|  \         edge 3 is node 0 to 3
          0 - -|- + 2       edge 4 is node 1 to 3
            \  |  +         edge 5 is node 2 to 3
            0\ | /1
              +|/           edge 2 is behind edge 4
               1


    */

    // lets start with making the vectors
    VerdictVector edges[6];
    edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                  coordinates[1][2] - coordinates[0][2] );

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

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

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

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

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

    // common numbers
    static const double root_of_2 = sqrt( 2.0 );

    // calculate the jacobian
    static const int do_jacobian = V_TET_JACOBIAN | V_TET_VOLUME | V_TET_ASPECT_BETA | V_TET_ASPECT_GAMMA |
                                   V_TET_SHAPE | V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE |
                                   V_TET_SCALED_JACOBIAN | V_TET_CONDITION;
    if( metrics_request_flag & do_jacobian ) { metric_vals->jacobian = (double)( edges[3] % ( edges[2] * edges[0] ) ); }

    // calculate the volume
    if( metrics_request_flag & V_TET_VOLUME ) { metric_vals->volume = (double)( metric_vals->jacobian / 6.0 ); }

    // calculate aspect ratio
    if( metrics_request_flag & V_TET_ASPECT_BETA )
    {
        double surface_area = ( ( edges[2] * edges[0] ).length() + ( edges[3] * edges[0] ).length() +
                                ( edges[4] * edges[1] ).length() + ( edges[3] * edges[2] ).length() ) *
                              0.5;

        VerdictVector numerator = edges[3].length_squared() * ( edges[2] * edges[0] ) +
                                  edges[2].length_squared() * ( edges[3] * edges[0] ) +
                                  edges[0].length_squared() * ( edges[3] * edges[2] );

        double volume = metric_vals->jacobian / 6.0;

        if( volume < VERDICT_DBL_MIN )
            metric_vals->aspect_beta = (double)( VERDICT_DBL_MAX );
        else
            metric_vals->aspect_beta = (double)( numerator.length() * surface_area / ( 108 * volume * volume ) );
    }

    // calculate the aspect gamma
    if( metrics_request_flag & V_TET_ASPECT_GAMMA )
    {
        double volume = fabs( metric_vals->jacobian / 6.0 );
        if( fabs( volume ) < VERDICT_DBL_MIN )
            metric_vals->aspect_gamma = VERDICT_DBL_MAX;
        else
        {
            double srms = sqrt( ( edges[0].length_squared() + edges[1].length_squared() + edges[2].length_squared() +
                                  edges[3].length_squared() + edges[4].length_squared() + edges[5].length_squared() ) /
                                6.0 );

            // cube the srms
            srms *= ( srms * srms );
            metric_vals->aspect_gamma = (double)( srms / ( 8.48528137423857 * volume ) );
        }
    }

    // calculate the shape of the tet
    if( metrics_request_flag & ( V_TET_SHAPE | V_TET_SHAPE_AND_SIZE ) )
    {
        // if the jacobian is non-positive, the shape is 0
        if( metric_vals->jacobian < VERDICT_DBL_MIN ) { metric_vals->shape = (double)0.0; }
        else
        {
            static const double two_thirds = 2.0 / 3.0;
            double num                     = 3.0 * pow( root_of_2 * metric_vals->jacobian, two_thirds );
            double den                     = 1.5 * ( edges[0] % edges[0] + edges[2] % edges[2] + edges[3] % edges[3] ) -
                         ( edges[0] % -edges[2] + -edges[2] % edges[3] + edges[3] % edges[0] );

            if( den < VERDICT_DBL_MIN )
                metric_vals->shape = (double)0.0;
            else
                metric_vals->shape = (double)VERDICT_MAX( num / den, 0 );
        }
    }

    // calculate the relative size of the tet
    if( metrics_request_flag & ( V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE ) )
    {
        VerdictVector w1, w2, w3;
        get_weight( w1, w2, w3 );
        double avg_vol = ( w1 % ( w2 * w3 ) ) / 6;

        if( avg_vol < VERDICT_DBL_MIN )
            metric_vals->relative_size_squared = 0.0;
        else
        {
            double tmp = metric_vals->jacobian / ( 6 * avg_vol );
            if( tmp < VERDICT_DBL_MIN )
                metric_vals->relative_size_squared = 0.0;
            else
            {
                tmp *= tmp;
                metric_vals->relative_size_squared = (double)VERDICT_MIN( tmp, 1 / tmp );
            }
        }
    }

    // calculate the shape and size
    if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
    { metric_vals->shape_and_size = (double)( metric_vals->shape * metric_vals->relative_size_squared ); }

    // calculate the scaled jacobian
    if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
    {
        // find out which node the normalized jacobian can be calculated at
        // and it will be the smaller than at other nodes
        double length_squared[4] = { edges[0].length_squared() * edges[2].length_squared() * edges[3].length_squared(),
                                     edges[0].length_squared() * edges[1].length_squared() * edges[4].length_squared(),
                                     edges[1].length_squared() * edges[2].length_squared() * edges[5].length_squared(),
                                     edges[3].length_squared() * edges[4].length_squared() *
                                         edges[5].length_squared() };

        int which_node = 0;
        if( length_squared[1] > length_squared[which_node] ) which_node = 1;
        if( length_squared[2] > length_squared[which_node] ) which_node = 2;
        if( length_squared[3] > length_squared[which_node] ) which_node = 3;

        // find the scaled jacobian at this node
        double length_product = sqrt( length_squared[which_node] );
        if( length_product < fabs( metric_vals->jacobian ) ) length_product = fabs( metric_vals->jacobian );

        if( length_product < VERDICT_DBL_MIN )
            metric_vals->scaled_jacobian = (double)VERDICT_DBL_MAX;
        else
            metric_vals->scaled_jacobian = (double)( root_of_2 * metric_vals->jacobian / length_product );
    }

    // calculate the condition number
    if( metrics_request_flag & V_TET_CONDITION )
    {
        static const double root_of_3 = sqrt( 3.0 );
        static const double root_of_6 = sqrt( 6.0 );

        VerdictVector c_1, c_2, c_3;
        c_1 = edges[0];
        c_2 = ( -2 * edges[2] - edges[0] ) / root_of_3;
        c_3 = ( 3 * edges[3] + edges[2] - edges[0] ) / root_of_6;

        double term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
        double term2 = ( c_1 * c_2 ) % ( c_1 * c_2 ) + ( c_2 * c_3 ) % ( c_2 * c_3 ) + ( c_3 * c_1 ) % ( c_3 * c_1 );

        double det = c_1 % ( c_2 * c_3 );

        if( det <= VERDICT_DBL_MIN )
            metric_vals->condition = (double)VERDICT_DBL_MAX;
        else
            metric_vals->condition = (double)( sqrt( term1 * term2 ) / ( 3.0 * det ) );
    }

    // calculate the distortion
    if( metrics_request_flag & V_TET_DISTORTION )
    { metric_vals->distortion = v_tet_distortion( num_nodes, coordinates ); }

    // check for overflow
    if( metrics_request_flag & V_TET_ASPECT_BETA )
    {
        if( metric_vals->aspect_beta > 0 )
            metric_vals->aspect_beta = (double)VERDICT_MIN( metric_vals->aspect_beta, VERDICT_DBL_MAX );
        metric_vals->aspect_beta = (double)VERDICT_MAX( metric_vals->aspect_beta, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_ASPECT_GAMMA )
    {
        if( metric_vals->aspect_gamma > 0 )
            metric_vals->aspect_gamma = (double)VERDICT_MIN( metric_vals->aspect_gamma, VERDICT_DBL_MAX );
        metric_vals->aspect_gamma = (double)VERDICT_MAX( metric_vals->aspect_gamma, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_VOLUME )
    {
        if( metric_vals->volume > 0 ) metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
        metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_CONDITION )
    {
        if( metric_vals->condition > 0 )
            metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
        metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_JACOBIAN )
    {
        if( metric_vals->jacobian > 0 )
            metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
        metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
    {
        if( metric_vals->scaled_jacobian > 0 )
            metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
        metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_SHAPE )
    {
        if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
        metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_RELATIVE_SIZE_SQUARED )
    {
        if( metric_vals->relative_size_squared > 0 )
            metric_vals->relative_size_squared =
                (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
        metric_vals->relative_size_squared =
            (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
    {
        if( metric_vals->shape_and_size > 0 )
            metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
        metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_DISTORTION )
    {
        if( metric_vals->distortion > 0 )
            metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
        metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
    }

    if( metrics_request_flag & V_TET_ASPECT_RATIO ) metric_vals->aspect_ratio = v_tet_aspect_ratio( 4, coordinates );

    if( metrics_request_flag & V_TET_ASPECT_FROBENIUS )
        metric_vals->aspect_frobenius = v_tet_aspect_frobenius( 4, coordinates );

    if( metrics_request_flag & V_TET_MINIMUM_ANGLE ) metric_vals->minimum_angle = v_tet_minimum_angle( 4, coordinates );

    if( metrics_request_flag & V_TET_COLLAPSE_RATIO )
        metric_vals->collapse_ratio = v_tet_collapse_ratio( 4, coordinates );

    if( metrics_request_flag & V_TET_RADIUS_RATIO ) metric_vals->radius_ratio = v_tet_radius_ratio( 4, coordinates );
}
C_FUNC_DEF double v_tet_radius_ratio ( int  ,
double  coordinates[][3] 
)

Calculates tet radius ratio metric.

CR / (3.0 * IR) where CR = circumsphere radius, IR = inscribed sphere radius. Reference --- V. N. Parthasarathy et al, A comparison of tetrahedron quality measures, Finite Elem. Anal. Des., Vol 15(1993), 255-261.

The radius ratio of a tet

NB (P. Pebay 04/16/07): CR / (3.0 * IR) where CR is the circumsphere radius and IR is the inscribed sphere radius. Note that this metric is similar to the aspect beta of a tet, except that it does not return VERDICT_DBL_MAX if the element has negative orientation.

Definition at line 209 of file V_TetMetric.cpp.

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

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

{

    // Determine side vectors
    VerdictVector side[6];

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

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

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

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

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

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

    VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
                              side[2].length_squared() * ( side[3] * side[0] ) +
                              side[0].length_squared() * ( side[3] * side[2] );

    double area_sum;
    area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
                 ( side[3] * side[2] ).length() ) *
               0.5;

    double volume = v_tet_volume( 4, coordinates );

    if( fabs( volume ) < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;
    else
    {
        double radius_ratio;
        radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );

        return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
    }
}
C_FUNC_DEF double v_tet_relative_size_squared ( int  ,
double  coordinates[][3] 
)

Calculates tet relative size metric.

Min( J, 1/J ), where J is determinant of weighted Jacobian matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

the relative size of a tet

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

Definition at line 724 of file V_TetMetric.cpp.

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

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

{
    double size;
    VerdictVector w1, w2, w3;
    get_weight( w1, w2, w3 );
    double avg_volume = ( w1 % ( w2 * w3 ) ) / 6.0;

    double volume = v_tet_volume( 4, coordinates );

    if( avg_volume < VERDICT_DBL_MIN )
        return 0.0;
    else
    {
        size = volume / avg_volume;
        if( size <= VERDICT_DBL_MIN ) return 0.0;
        if( size > 1 ) size = (double)( 1 ) / size;
    }
    return (double)( size * size );
}
C_FUNC_DEF double v_tet_scaled_jacobian ( int  ,
double  coordinates[][3] 
)

Calculates tet scaled jacobian.

Minimum Jacobian divided by the lengths of 3 edge vectors Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

the scaled jacobian of a tet

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

Definition at line 153 of file V_TetMetric.cpp.

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

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

{

    VerdictVector side0, side1, side2, side3, side4, side5;

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

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

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

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

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

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

    double jacobi;

    jacobi = side3 % ( side2 * side0 );

    // products of lengths squared of each edge attached to a node.
    double length_squared[4] = { side0.length_squared() * side2.length_squared() * side3.length_squared(),
                                 side0.length_squared() * side1.length_squared() * side4.length_squared(),
                                 side1.length_squared() * side2.length_squared() * side5.length_squared(),
                                 side3.length_squared() * side4.length_squared() * side5.length_squared() };
    int which_node           = 0;
    if( length_squared[1] > length_squared[which_node] ) which_node = 1;
    if( length_squared[2] > length_squared[which_node] ) which_node = 2;
    if( length_squared[3] > length_squared[which_node] ) which_node = 3;

    double length_product = sqrt( length_squared[which_node] );
    if( length_product < fabs( jacobi ) ) length_product = fabs( jacobi );

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

    static const double root_of_2 = sqrt( 2.0 );

    return (double)( root_of_2 * jacobi / length_product );
}
C_FUNC_DEF double v_tet_shape ( int  ,
double  coordinates[][3] 
)

Calculates tet shape metric.

3/Mean Ratio of weighted Jacobian matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

the shape of a tet

3/ condition number of weighted jacobian matrix

Definition at line 691 of file V_TetMetric.cpp.

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

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

{

    static const double two_thirds = 2.0 / 3.0;
    static const double root_of_2  = sqrt( 2.0 );

    VerdictVector edge0, edge2, edge3;

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

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

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

    double jacobian = edge3 % ( edge2 * edge0 );
    if( jacobian < VERDICT_DBL_MIN ) { return (double)0.0; }
    double num = 3 * pow( root_of_2 * jacobian, two_thirds );
    double den =
        1.5 * ( edge0 % edge0 + edge2 % edge2 + edge3 % edge3 ) - ( edge0 % -edge2 + -edge2 % edge3 + edge3 % edge0 );

    if( den < VERDICT_DBL_MIN ) return (double)0.0;

    return (double)VERDICT_MAX( num / den, 0 );
}
C_FUNC_DEF double v_tet_shape_and_size ( int  num_nodes,
double  coordinates[][3] 
)

Calculates tet shape-size metric.

Product of Shape and Relative Size. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

the shape and size of a tet

Product of the shape and relative size

Definition at line 749 of file V_TetMetric.cpp.

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

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

{

    double shape, size;
    shape = v_tet_shape( num_nodes, coordinates );
    size  = v_tet_relative_size_squared( num_nodes, coordinates );

    return (double)( shape * size );
}
C_FUNC_DEF double v_tet_volume ( int  ,
double  coordinates[][3] 
)

Calculates tet volume.

(1/6) * Jacobian at corner node. Reference --- V. N. Parthasarathy et al, A comparison of tetrahedron quality measures, Finite Elem. Anal. Des., Vol 15(1993), 255-261.

the volume of a tet

1/6 * jacobian at a corner node

Definition at line 606 of file V_TetMetric.cpp.

References VerdictVector::set().

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

{

    // Determine side vectors
    VerdictVector side0, side2, side3;

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

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

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

    return (double)( ( side3 % ( side2 * side0 ) ) / 6.0 );
}
C_FUNC_DEF double v_tri_area ( int  ,
double  coordinates[][3] 
)

Calculates tri metric.

Maximum included angle in triangle

The area of a tri

0.5 * jacobian at a node

Definition at line 278 of file V_TriMetric.cpp.

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

Referenced by moab::VerdictWrapper::quality_measure(), v_quad_jacobian(), and v_quad_quality().

{
    // two vectors for two sides
    VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                         coordinates[1][2] - coordinates[0][2] );

    VerdictVector side3( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
                         coordinates[2][2] - coordinates[0][2] );

    // the cross product of the two vectors representing two sides of the
    // triangle
    VerdictVector tmp = side1 * side3;

    // return the magnitude of the vector divided by two
    double area = 0.5 * tmp.length();
    if( area > 0 ) return (double)VERDICT_MIN( area, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( area, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_tri_aspect_frobenius ( int  ,
double  coordinates[][3] 
)

Calculates tri metric.

Frobenius aspect

the Frobenius aspect of a tri

srms^2/(2 * sqrt(3.0) * area) where srms^2 is sum of the lengths squared

NB (P. Pebay 01/14/07): this method was called "aspect ratio" in earlier incarnations of VERDICT

Definition at line 246 of file V_TriMetric.cpp.

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

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

{
    static const double two_times_root_of_3 = 2 * sqrt( 3.0 );

    // three vectors for each side
    VerdictVector side1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                         coordinates[1][2] - coordinates[0][2] );

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

    VerdictVector side3( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
                         coordinates[0][2] - coordinates[2][2] );

    // sum the lengths squared of each side
    double srms = ( side1.length_squared() + side2.length_squared() + side3.length_squared() );

    // find two times the area of the triangle by cross product
    double areaX2 = ( ( side1 * ( -side3 ) ).length() );

    if( areaX2 == 0.0 ) return (double)VERDICT_DBL_MAX;

    double aspect = (double)( srms / ( two_times_root_of_3 * ( areaX2 ) ) );
    if( aspect > 0 ) return (double)VERDICT_MIN( aspect, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( aspect, -VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_tri_aspect_ratio ( int  ,
double  coordinates[][3] 
)

Calculates tri metric.

aspect ratio Reference --- P. P. Pebay & T. J. Baker, Analysis of Triangle Quality Measures, AMS Math. Comp., 2003, 72(244):1817-1839

the aspect ratio of a triangle

NB (P. Pebay 01/14/07): Hmax / ( 2.0 * sqrt(3.0) * IR) where Hmax is the maximum edge length and IR is the inradius

note that previous incarnations of verdict used "v_tri_aspect_ratio" to denote what is now called "v_tri_aspect_frobenius"

Definition at line 160 of file V_TriMetric.cpp.

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

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

{
    static const double normal_coeff = sqrt( 3. ) / 6.;

    // three vectors for each side
    VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                     coordinates[1][2] - coordinates[0][2] );

    VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
                     coordinates[2][2] - coordinates[1][2] );

    VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
                     coordinates[0][2] - coordinates[2][2] );

    double a1 = a.length();
    double b1 = b.length();
    double c1 = c.length();

    double hm = a1 > b1 ? a1 : b1;
    hm        = hm > c1 ? hm : c1;

    VerdictVector ab   = a * b;
    double denominator = ab.length();

    if( denominator < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;
    else
    {
        double aspect_ratio;
        aspect_ratio = normal_coeff * hm * ( a1 + b1 + c1 ) / denominator;

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

Calculates tri metric.

Condition number of the Jacobian matrix. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

The condition of a tri

Condition number of the jacobian matrix at any corner

Definition at line 412 of file V_TriMetric.cpp.

References compute_normal, VerdictVector::length(), VERDICT_DBL_MAX, VERDICT_MIN, VerdictVector::x(), VerdictVector::y(), and VerdictVector::z().

Referenced by moab::VerdictWrapper::quality_measure(), v_quad_condition(), and v_tri_shape().

{
    static const double rt3 = sqrt( 3.0 );

    VerdictVector v1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                      coordinates[1][2] - coordinates[0][2] );

    VerdictVector v2( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
                      coordinates[2][2] - coordinates[0][2] );

    VerdictVector tri_normal = v1 * v2;
    double areax2            = tri_normal.length();

    if( areax2 == 0.0 ) return (double)VERDICT_DBL_MAX;

    double condition = (double)( ( ( v1 % v1 ) + ( v2 % v2 ) - ( v1 % v2 ) ) / ( areax2 * rt3 ) );

    // check for inverted if we have access to the normal
    if( compute_normal )
    {
        // center of tri
        double point[3], surf_normal[3];
        point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
        point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
        point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;

        // dot product
        compute_normal( point, surf_normal );
        if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
            0 )
            return (double)VERDICT_DBL_MAX;
    }
    return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
}
C_FUNC_DEF double v_tri_distortion ( int  num_nodes,
double  coordinates[][3] 
)

Calculates tri metric.

{min(|J|)/actual area}*parent area, parent area = 1/2 for triangular element. Reference --- SDRC/IDEAS Simulation: Finite Element Modeling--User's Guide

The distortion of a tri

TODO: make a short definition of the distortion and comment below

Definition at line 579 of file V_TriMetric.cpp.

References GaussIntegration::calculate_derivative_at_nodes_2d_tri(), GaussIntegration::calculate_shape_function_2d_tri(), dot_product(), GaussIntegration::get_shape_func(), GaussIntegration::initialize(), VerdictVector::length(), maxNumberNodes, maxTotalNumberGaussPoints, VerdictVector::normalize(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{

    double distortion;
    int total_number_of_gauss_points = 0;
    VerdictVector aa, bb, cc, normal_at_point, xin;
    double element_area = 0.;

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

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

    VerdictVector tri_normal = aa * bb;

    int number_of_gauss_points = 0;
    if( num_nodes == 3 )
    {
        distortion = 1.0;
        return (double)distortion;
    }

    else if( num_nodes == 6 )
    {
        total_number_of_gauss_points = 6;
        number_of_gauss_points       = 6;
    }

    distortion = VERDICT_DBL_MAX;
    double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
    double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
    double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
    double weight[maxTotalNumberGaussPoints];

    // create an object of GaussIntegration
    int number_dims = 2;
    int is_tri      = 1;
    GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
    GaussIntegration::calculate_shape_function_2d_tri();
    GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );

    // calculate element area
    int ife, ja;
    for( ife = 0; ife < total_number_of_gauss_points; ife++ )
    {
        aa.set( 0.0, 0.0, 0.0 );
        bb.set( 0.0, 0.0, 0.0 );

        for( ja = 0; ja < num_nodes; ja++ )
        {
            xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
            aa += dndy1[ife][ja] * xin;
            bb += dndy2[ife][ja] * xin;
        }
        normal_at_point = aa * bb;
        double jacobian = normal_at_point.length();
        element_area += weight[ife] * jacobian;
    }

    element_area *= 0.8660254;
    double dndy1_at_node[maxNumberNodes][maxNumberNodes];
    double dndy2_at_node[maxNumberNodes][maxNumberNodes];

    GaussIntegration::calculate_derivative_at_nodes_2d_tri( dndy1_at_node, dndy2_at_node );

    VerdictVector normal_at_nodes[7];

    // evaluate normal at nodes and distortion values at nodes
    int jai = 0;
    for( ja = 0; ja < num_nodes; ja++ )
    {
        aa.set( 0.0, 0.0, 0.0 );
        bb.set( 0.0, 0.0, 0.0 );
        for( jai = 0; jai < num_nodes; jai++ )
        {
            xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
            aa += dndy1_at_node[ja][jai] * xin;
            bb += dndy2_at_node[ja][jai] * xin;
        }
        normal_at_nodes[ja] = aa * bb;
        normal_at_nodes[ja].normalize();
    }

    // determine if element is flat
    bool flat_element = true;
    double dot_product;

    for( ja = 0; ja < num_nodes; ja++ )
    {
        dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
        if( fabs( dot_product ) < 0.99 )
        {
            flat_element = false;
            break;
        }
    }

    // take into consideration of the thickness of the element
    double thickness, thickness_gauss;
    double distrt;
    // get_tri_thickness(tri, element_area, thickness );
    thickness = 0.001 * sqrt( element_area );

    // set thickness gauss point location
    double zl = 0.5773502691896;
    if( flat_element ) zl = 0.0;

    int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
    double thickness_z;

    // loop on integration points
    int igz;
    for( ife = 0; ife < total_number_of_gauss_points; ife++ )
    {
        // loop on the thickness direction gauss points
        for( igz = 0; igz < no_gauss_pts_z; igz++ )
        {
            zl          = -zl;
            thickness_z = zl * thickness / 2.0;

            aa.set( 0.0, 0.0, 0.0 );
            bb.set( 0.0, 0.0, 0.0 );
            cc.set( 0.0, 0.0, 0.0 );

            for( ja = 0; ja < num_nodes; ja++ )
            {
                xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
                xin += thickness_z * normal_at_nodes[ja];
                aa += dndy1[ife][ja] * xin;
                bb += dndy2[ife][ja] * xin;
                thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
                cc += thickness_gauss * normal_at_nodes[ja];
            }

            normal_at_point = aa * bb;
            distrt          = cc % normal_at_point;
            if( distrt < distortion ) distortion = distrt;
        }
    }

    // loop through nodal points
    for( ja = 0; ja < num_nodes; ja++ )
    {
        for( igz = 0; igz < no_gauss_pts_z; igz++ )
        {
            zl          = -zl;
            thickness_z = zl * thickness / 2.0;

            aa.set( 0.0, 0.0, 0.0 );
            bb.set( 0.0, 0.0, 0.0 );
            cc.set( 0.0, 0.0, 0.0 );

            for( jai = 0; jai < num_nodes; jai++ )
            {
                xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
                xin += thickness_z * normal_at_nodes[ja];
                aa += dndy1_at_node[ja][jai] * xin;
                bb += dndy2_at_node[ja][jai] * xin;
                if( jai == ja )
                    thickness_gauss = thickness / 2.0;
                else
                    thickness_gauss = 0.;
                cc += thickness_gauss * normal_at_nodes[jai];
            }
        }

        normal_at_point      = aa * bb;
        double sign_jacobian = ( tri_normal % normal_at_point ) > 0 ? 1. : -1.;
        distrt               = sign_jacobian * ( cc % normal_at_point );

        if( distrt < distortion ) distortion = distrt;
    }
    if( element_area * thickness != 0 )
        distortion *= 1. / ( element_area * thickness );
    else
        distortion *= 1.;

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

Calculates tri metric.

edge ratio Reference --- P. P. Pebay & T. J. Baker, Analysis of Triangle Quality Measures, AMS Math. Comp., 2003, 72(244):1817-1839

the edge ratio of a triangle

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

Definition at line 76 of file V_TriMetric.cpp.

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

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

{

    // three vectors for each side
    VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                     coordinates[1][2] - coordinates[0][2] );

    VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
                     coordinates[2][2] - coordinates[1][2] );

    VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
                     coordinates[0][2] - coordinates[2][2] );

    double a2 = a.length_squared();
    double b2 = b.length_squared();
    double c2 = c.length_squared();

    double m2, M2;
    if( a2 < b2 )
    {
        if( b2 < c2 )
        {
            m2 = a2;
            M2 = c2;
        }
        else  // b2 <= a2
        {
            if( a2 < c2 )
            {
                m2 = a2;
                M2 = b2;
            }
            else  // c2 <= a2
            {
                m2 = c2;
                M2 = b2;
            }
        }
    }
    else  // b2 <= a2
    {
        if( a2 < c2 )
        {
            m2 = b2;
            M2 = c2;
        }
        else  // c2 <= a2
        {
            if( b2 < c2 )
            {
                m2 = b2;
                M2 = a2;
            }
            else  // c2 <= b2
            {
                m2 = c2;
                M2 = a2;
            }
        }
    }

    if( m2 < VERDICT_DBL_MIN )
        return (double)VERDICT_DBL_MAX;
    else
    {
        double edge_ratio;
        edge_ratio = sqrt( M2 / m2 );

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

Calculates tri metric.

Maximum included angle in triangle

The maximum angle of a tri

The maximum angle of a tri is the maximum angle between two adjacents sides out of all three corners of the triangle.

Definition at line 358 of file V_TriMetric.cpp.

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

Referenced by moab::VerdictWrapper::quality_measure(), v_quad_maximum_angle(), and v_quad_quality().

{

    // vectors for all the sides
    VerdictVector sides[4];
    sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                  coordinates[1][2] - coordinates[0][2] );
    sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
                  coordinates[2][2] - coordinates[1][2] );
    sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
                  coordinates[2][2] - coordinates[0][2] );

    // in case we need to find the interior angle
    // between sides 0 and 1
    sides[3] = -sides[1];

    // calculate the lengths squared of the sides
    double sides_lengths[3];
    sides_lengths[0] = sides[0].length_squared();
    sides_lengths[1] = sides[1].length_squared();
    sides_lengths[2] = sides[2].length_squared();

    if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 ) { return 0.0; }

    // using the law of sines, we know that the maximum
    // angle is opposite of the longest side

    // find the longest side
    int short_side = 0;
    if( sides_lengths[1] > sides_lengths[0] ) short_side = 1;
    if( sides_lengths[2] > sides_lengths[short_side] ) short_side = 2;

    // from the longest side, calculate the angle of the
    // opposite angle
    double max_angle;
    if( short_side == 0 ) { max_angle = sides[2].interior_angle( sides[1] ); }
    else if( short_side == 1 )
    {
        max_angle = sides[0].interior_angle( sides[2] );
    }
    else
    {
        max_angle = sides[0].interior_angle( sides[3] );
    }

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

Calculates tri metric.

Minimum included angle in triangle

The minimum angle of a tri

The minimum angle of a tri is the minimum angle between two adjacents sides out of all three corners of the triangle.

Definition at line 303 of file V_TriMetric.cpp.

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

Referenced by moab::VerdictWrapper::quality_measure(), v_quad_minimum_angle(), and v_quad_quality().

{

    // vectors for all the sides
    VerdictVector sides[4];
    sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                  coordinates[1][2] - coordinates[0][2] );
    sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
                  coordinates[2][2] - coordinates[1][2] );
    sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
                  coordinates[2][2] - coordinates[0][2] );

    // in case we need to find the interior angle
    // between sides 0 and 1
    sides[3] = -sides[1];

    // calculate the lengths squared of the sides
    double sides_lengths[3];
    sides_lengths[0] = sides[0].length_squared();
    sides_lengths[1] = sides[1].length_squared();
    sides_lengths[2] = sides[2].length_squared();

    if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 ) return 0.0;

    // using the law of sines, we know that the minimum
    // angle is opposite of the shortest side

    // find the shortest side
    int short_side = 0;
    if( sides_lengths[1] < sides_lengths[0] ) short_side = 1;
    if( sides_lengths[2] < sides_lengths[short_side] ) short_side = 2;

    // from the shortest side, calculate the angle of the
    // opposite angle
    double min_angle;
    if( short_side == 0 ) { min_angle = sides[2].interior_angle( sides[1] ); }
    else if( short_side == 1 )
    {
        min_angle = sides[0].interior_angle( sides[2] );
    }
    else
    {
        min_angle = sides[0].interior_angle( sides[3] );
    }

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

Calculates quality metrics for triangle elements.

tri_quality for calculating multiple tri metrics at once

using this method is generally faster than using the individual method multiple times.

Definition at line 768 of file V_TriMetric.cpp.

References TriMetricVals::area, TriMetricVals::aspect_frobenius, TriMetricVals::aspect_ratio, compute_normal, TriMetricVals::condition, determinant(), TriMetricVals::distortion, TriMetricVals::edge_ratio, interior_angle(), length(), VerdictVector::length(), TriMetricVals::maximum_angle, TriMetricVals::minimum_angle, TriMetricVals::radius_ratio, TriMetricVals::relative_size_squared, TriMetricVals::scaled_jacobian, VerdictVector::set(), TriMetricVals::shape, TriMetricVals::shape_and_size, size, V_TRI_AREA, V_TRI_ASPECT_FROBENIUS, V_TRI_CONDITION, V_TRI_DISTORTION, v_tri_distortion(), V_TRI_EDGE_RATIO, v_tri_edge_ratio(), v_tri_get_weight(), V_TRI_MAXIMUM_ANGLE, V_TRI_MINIMUM_ANGLE, V_TRI_RADIUS_RATIO, v_tri_radius_ratio(), V_TRI_RELATIVE_SIZE_SQUARED, V_TRI_SCALED_JACOBIAN, V_TRI_SHAPE, V_TRI_SHAPE_AND_SIZE, VERDICT_DBL_MAX, VERDICT_MAX, VERDICT_MIN, VerdictVector::x(), VerdictVector::y(), and VerdictVector::z().

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

{

    memset( metric_vals, 0, sizeof( TriMetricVals ) );

    // for starts, lets set up some basic and common information

    /*  node numbers and side numbers used below

               2
               ++
              /  \
           2 /    \ 1
            /      \
           /        \
         0 ---------+ 1
               0
    */

    // vectors for each side
    VerdictVector sides[3];
    sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                  coordinates[1][2] - coordinates[0][2] );
    sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
                  coordinates[2][2] - coordinates[1][2] );
    sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
                  coordinates[2][2] - coordinates[0][2] );
    VerdictVector tri_normal = sides[0] * sides[2];
    // if we have access to normal information, check to see if the
    // element is inverted.  If we don't have the normal information
    // that we need for this, assume the element is not inverted.
    // This flag will be used for condition number, jacobian, shape,
    // and size and shape.
    bool is_inverted = false;
    if( compute_normal )
    {
        // center of tri
        double point[3], surf_normal[3];
        point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
        point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
        point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
        // dot product
        compute_normal( point, surf_normal );
        if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
            0 )
            is_inverted = true;
    }

    // lengths squared of each side
    double sides_lengths_squared[3];
    sides_lengths_squared[0] = sides[0].length_squared();
    sides_lengths_squared[1] = sides[1].length_squared();
    sides_lengths_squared[2] = sides[2].length_squared();

    // if we are doing angle calcuations
    if( metrics_request_flag & ( V_TRI_MINIMUM_ANGLE | V_TRI_MAXIMUM_ANGLE ) )
    {
        // which is short and long side
        int short_side = 0, long_side = 0;

        if( sides_lengths_squared[1] < sides_lengths_squared[0] ) short_side = 1;
        if( sides_lengths_squared[2] < sides_lengths_squared[short_side] ) short_side = 2;

        if( sides_lengths_squared[1] > sides_lengths_squared[0] ) long_side = 1;
        if( sides_lengths_squared[2] > sides_lengths_squared[long_side] ) long_side = 2;

        // calculate the minimum angle of the tri
        if( metrics_request_flag & V_TRI_MINIMUM_ANGLE )
        {
            if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
            { metric_vals->minimum_angle = 0.0; }
            else if( short_side == 0 )
                metric_vals->minimum_angle = (double)sides[2].interior_angle( sides[1] );
            else if( short_side == 1 )
                metric_vals->minimum_angle = (double)sides[0].interior_angle( sides[2] );
            else
                metric_vals->minimum_angle = (double)sides[0].interior_angle( -sides[1] );
        }

        // calculate the maximum angle of the tri
        if( metrics_request_flag & V_TRI_MAXIMUM_ANGLE )
        {
            if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
            { metric_vals->maximum_angle = 0.0; }
            else if( long_side == 0 )
                metric_vals->maximum_angle = (double)sides[2].interior_angle( sides[1] );
            else if( long_side == 1 )
                metric_vals->maximum_angle = (double)sides[0].interior_angle( sides[2] );
            else
                metric_vals->maximum_angle = (double)sides[0].interior_angle( -sides[1] );
        }
    }

    // calculate the area of the tri
    // the following metrics depend on area
    if( metrics_request_flag &
        ( V_TRI_AREA | V_TRI_SCALED_JACOBIAN | V_TRI_SHAPE | V_TRI_RELATIVE_SIZE_SQUARED | V_TRI_SHAPE_AND_SIZE ) )
    { metric_vals->area = (double)( ( sides[0] * sides[2] ).length() * 0.5 ); }

    // calculate the aspect ratio
    if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )
    {
        // sum the lengths squared
        double srms = sides_lengths_squared[0] + sides_lengths_squared[1] + sides_lengths_squared[2];

        // calculate once and reuse
        static const double twoTimesRootOf3 = 2 * sqrt( 3.0 );

        double div = ( twoTimesRootOf3 * ( ( sides[0] * sides[2] ).length() ) );

        if( div == 0.0 )
            metric_vals->aspect_frobenius = (double)VERDICT_DBL_MAX;
        else
            metric_vals->aspect_frobenius = (double)( srms / div );
    }

    // calculate the scaled jacobian
    if( metrics_request_flag & V_TRI_SCALED_JACOBIAN )
    {
        // calculate once and reuse
        static const double twoOverRootOf3 = 2 / sqrt( 3.0 );
        // use the area from above

        double tmp = tri_normal.length() * twoOverRootOf3;

        // now scale it by the lengths of the sides
        double min_scaled_jac = VERDICT_DBL_MAX;
        double temp_scaled_jac;
        for( int i = 0; i < 3; i++ )
        {
            if( sides_lengths_squared[i % 3] == 0.0 || sides_lengths_squared[( i + 2 ) % 3] == 0.0 )
                temp_scaled_jac = 0.0;
            else
                temp_scaled_jac =
                    tmp / sqrt( sides_lengths_squared[i % 3] ) / sqrt( sides_lengths_squared[( i + 2 ) % 3] );
            if( temp_scaled_jac < min_scaled_jac ) min_scaled_jac = temp_scaled_jac;
        }
        // multiply by -1 if the normals are in opposite directions
        if( is_inverted ) { min_scaled_jac *= -1; }
        metric_vals->scaled_jacobian = (double)min_scaled_jac;
    }

    // calculate the condition number
    if( metrics_request_flag & V_TRI_CONDITION )
    {
        // calculate once and reuse
        static const double rootOf3 = sqrt( 3.0 );
        // if it is inverted, the condition number is considered to be infinity.
        if( is_inverted ) { metric_vals->condition = VERDICT_DBL_MAX; }
        else
        {
            double area2x = ( sides[0] * sides[2] ).length();
            if( area2x == 0.0 )
                metric_vals->condition = (double)( VERDICT_DBL_MAX );
            else
                metric_vals->condition = (double)( ( sides[0] % sides[0] + sides[2] % sides[2] - sides[0] % sides[2] ) /
                                                   ( area2x * rootOf3 ) );
        }
    }

    // calculate the shape
    if( metrics_request_flag & V_TRI_SHAPE || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
    {
        // if element is inverted, shape is zero.  We don't need to
        // calculate anything.
        if( is_inverted ) { metric_vals->shape = 0.0; }
        else
        {  // otherwise, we calculate the shape
            // calculate once and reuse
            static const double rootOf3 = sqrt( 3.0 );
            // reuse area from before
            double area2x = metric_vals->area * 2;
            // dot products
            double dots[3] = { sides[0] % sides[0], sides[2] % sides[2], sides[0] % sides[2] };

            // add the dots
            double sum_dots = dots[0] + dots[1] - dots[2];

            // then the finale
            if( sum_dots == 0.0 )
                metric_vals->shape = 0.0;
            else
                metric_vals->shape = (double)( rootOf3 * area2x / sum_dots );
        }
    }

    // calculate relative size squared
    if( metrics_request_flag & V_TRI_RELATIVE_SIZE_SQUARED || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
    {
        // get weights
        double w11, w21, w12, w22;
        v_tri_get_weight( w11, w21, w12, w22 );
        // get the determinant
        double detw = determinant( w11, w21, w12, w22 );
        // use the area from above and divide with the determinant
        if( metric_vals->area == 0.0 || detw == 0.0 )
            metric_vals->relative_size_squared = 0.0;
        else
        {
            double size = metric_vals->area * 2.0 / detw;
            // square the size
            size *= size;
            // value ranges between 0 to 1
            metric_vals->relative_size_squared = (double)VERDICT_MIN( size, 1.0 / size );
        }
    }

    // calculate shape and size
    if( metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
    { metric_vals->shape_and_size = metric_vals->relative_size_squared * metric_vals->shape; }

    // calculate distortion
    if( metrics_request_flag & V_TRI_DISTORTION ) metric_vals->distortion = v_tri_distortion( num_nodes, coordinates );

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

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

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

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

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

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

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

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

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

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

    if( metrics_request_flag & V_TRI_EDGE_RATIO ) { metric_vals->edge_ratio = v_tri_edge_ratio( 3, coordinates ); }
    if( metrics_request_flag & V_TRI_RADIUS_RATIO )
    { metric_vals->radius_ratio = v_tri_radius_ratio( 3, coordinates ); }
    if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )  // there is no V_TRI_ASPECT_RATIO !
    { metric_vals->aspect_ratio = v_tri_radius_ratio( 3, coordinates ); }
}
C_FUNC_DEF double v_tri_radius_ratio ( int  ,
double  coordinates[][3] 
)

Calculates tri metric.

radius ratio Reference --- P. P. Pebay & T. J. Baker, Analysis of Triangle Quality Measures, AMS Math. Comp., 2003, 72(244):1817-1839

the radius ratio of a triangle

NB (P. Pebay 01/13/07): CR / (3.0*IR) where CR is the circumradius and IR is the inradius

this quality metric is also known to VERDICT, for tetrahedral elements only, a the "aspect beta"

Definition at line 206 of file V_TriMetric.cpp.

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

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

{

    // three vectors for each side
    VerdictVector a( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
                     coordinates[1][2] - coordinates[0][2] );

    VerdictVector b( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
                     coordinates[2][2] - coordinates[1][2] );

    VerdictVector c( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
                     coordinates[0][2] - coordinates[2][2] );

    double a2 = a.length_squared();
    double b2 = b.length_squared();
    double c2 = c.length_squared();

    VerdictVector ab   = a * b;
    double denominator = ab.length_squared();

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

    double radius_ratio;
    radius_ratio = .25 * a2 * b2 * c2 * ( a2 + b2 + c2 ) / denominator;

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

Calculates tri metric.

Min( J, 1/J ), where J is determinant of weighted Jacobian matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

The relative size of a tri

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

Definition at line 525 of file V_TriMetric.cpp.

References determinant(), VerdictVector::length(), VerdictVector::set(), size, v_tri_get_weight(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    double w11, w21, w12, w22;

    VerdictVector xxi, xet, tri_normal;

    v_tri_get_weight( w11, w21, w12, w22 );

    double detw = determinant( w11, w21, w12, w22 );

    if( detw == 0.0 ) return 0.0;

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

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

    tri_normal = xxi * xet;

    double deta = tri_normal.length();
    if( deta == 0.0 || detw == 0.0 ) return 0.0;

    double size = pow( deta / detw, 2 );

    double rel_size = VERDICT_MIN( size, 1.0 / size );

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

Calculates tri metric.

Minimum Jacobian divided by the lengths of 2 edge vectors. Reference --- P. Knupp, Achieving Finite Element Mesh Quality via Optimization of the Jacobian Matrix Norm and Associated Quantities, Intl. J. Numer. Meth. Engng. 2000, 48:1165-1185.

The scaled jacobian of a tri

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

Definition at line 452 of file V_TriMetric.cpp.

References compute_normal, moab::cross(), moab::GeomUtil::first(), length(), VerdictVector::length(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, VERDICT_MIN, VerdictVector::x(), VerdictVector::y(), and VerdictVector::z().

Referenced by moab::VerdictWrapper::quality_measure(), v_quad_quality(), and v_quad_scaled_jacobian().

{
    static const double detw = 2. / sqrt( 3.0 );
    VerdictVector first, second;
    double jacobian;

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

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

    edge[2].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
                 coordinates[2][2] - coordinates[1][2] );
    first  = edge[1] - edge[0];
    second = edge[2] - edge[0];

    VerdictVector cross = first * second;
    jacobian            = cross.length();

    double max_edge_length_product;
    max_edge_length_product =
        VERDICT_MAX( edge[0].length() * edge[1].length(),
                     VERDICT_MAX( edge[1].length() * edge[2].length(), edge[0].length() * edge[2].length() ) );

    if( max_edge_length_product < VERDICT_DBL_MIN ) return (double)0.0;

    jacobian *= detw;
    jacobian /= max_edge_length_product;

    if( compute_normal )
    {
        // center of tri
        double point[3], surf_normal[3];
        point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
        point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
        point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;

        // dot product
        compute_normal( point, surf_normal );
        if( ( cross.x() * surf_normal[0] + cross.y() * surf_normal[1] + cross.z() * surf_normal[2] ) < 0 )
            jacobian *= -1;
    }

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

Calculates tri metric.

2/Condition number of weighted Jacobian matrix. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

The shape of a tri

2 / condition number of weighted jacobian matrix

Definition at line 506 of file V_TriMetric.cpp.

References v_tri_condition(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.

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

{
    double condition = v_tri_condition( num_nodes, coordinates );

    double shape;
    if( condition <= VERDICT_DBL_MIN )
        shape = VERDICT_DBL_MAX;
    else
        shape = ( 1 / condition );

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

Calculates tri metric.

Product of Shape and Relative Size. Reference --- P. Knupp, Algebraic Mesh Quality Metrics for Unstructured Initial Meshes, submitted for publication.

The shape and size of a tri

Product of the Shape and Relative Size

Definition at line 561 of file V_TriMetric.cpp.

References size, v_tri_relative_size_squared(), v_tri_shape(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.

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

{
    double size, shape;

    size  = v_tri_relative_size_squared( num_nodes, coordinates );
    shape = v_tri_shape( num_nodes, coordinates );

    double shape_and_size = size * shape;

    if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
    return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
}
C_FUNC_DEF void v_wedge_quality ( int  num_nodes,
double  coordinates[][3],
unsigned int  metrics_request_flag,
struct WedgeMetricVals metric_vals 
)

Calculates quality metrics for wedge elements.

Definition at line 102 of file V_WedgeMetric.cpp.

References V_WEDGE_VOLUME, v_wedge_volume(), and WedgeMetricVals::volume.

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

    if( metrics_request_flag & V_WEDGE_VOLUME ) metric_vals->volume = v_wedge_volume( num_nodes, coordinates );
}
C_FUNC_DEF double v_wedge_volume ( int  num_nodes,
double  coordinates[][3] 
)

Calculates wedge volume.

calculate the volume of a wedge

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

Definition at line 54 of file V_WedgeMetric.cpp.

References VerdictVector::set().

Referenced by moab::VerdictWrapper::all_quality_measures(), moab::VerdictWrapper::quality_measure(), and v_wedge_quality().

{

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

    if( num_nodes == 6 )
    {

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

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

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

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

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

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

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

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

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

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

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

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

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

    return (double)volume;
}
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines