 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
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 ( V_HEX_SKEW + V_HEX_TAPER ) |
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 ( V_TET_SHAPE + V_TET_RELATIVE_SIZE_SQUARED + V_TET_SHAPE_AND_SIZE ) |
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 ( V_QUAD_MAX_EDGE_RATIO + V_QUAD_SKEW + V_QUAD_TAPER ) |
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 ( V_TRI_SHAPE + V_TRI_SHAPE_AND_SIZE + V_TRI_RELATIVE_SIZE_SQUARED ) |
| #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. | |
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 C_FUNC_DEF |
| #define V_EDGE_LENGTH |
Definition at line 465 of file verdict.h.
Referenced by edge_quality().
| #define V_HEX_ALGEBRAIC |
| #define V_HEX_ALL |
Definition at line 340 of file verdict.h.
Referenced by moab::VerdictWrapper::all_quality_measures().
| #define V_HEX_CONDITION |
Definition at line 329 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_DIAGNOSTIC |
| #define V_HEX_DIAGONAL |
Definition at line 325 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_DIMENSION |
Definition at line 326 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_DISTORTION |
Definition at line 337 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_EDGE_RATIO |
Definition at line 338 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_JACOBIAN |
Definition at line 330 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_MAX_ASPECT_FROBENIUS |
| #define V_HEX_MAX_EDGE_RATIO |
Definition at line 320 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_MED_ASPECT_FROBENIUS |
Definition at line 339 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_ODDY |
Definition at line 327 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_RELATIVE_SIZE_SQUARED |
Definition at line 334 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_ROBINSON ( V_HEX_SKEW + V_HEX_TAPER ) |
| #define V_HEX_SCALED_JACOBIAN |
Definition at line 331 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_SHAPE |
Definition at line 333 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_SHAPE_AND_SIZE |
Definition at line 335 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_SHEAR |
Definition at line 332 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_SHEAR_AND_SIZE |
Definition at line 336 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_SKEW |
Definition at line 321 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_STRETCH |
Definition at line 324 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_TAPER |
Definition at line 322 of file verdict.h.
Referenced by v_hex_quality().
| #define V_HEX_TRADITIONAL |
| #define V_HEX_VOLUME |
Definition at line 323 of file verdict.h.
Referenced by v_hex_quality().
| #define V_KNIFE_VOLUME |
Definition at line 399 of file verdict.h.
Referenced by v_knife_quality().
| #define V_PYRAMID_VOLUME |
Definition at line 387 of file verdict.h.
Referenced by v_pyramid_quality().
| #define V_QUAD_ALGEBRAIC |
| #define V_QUAD_ALL |
Definition at line 428 of file verdict.h.
Referenced by moab::VerdictWrapper::all_quality_measures().
| #define V_QUAD_AREA |
Definition at line 409 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_ASPECT_RATIO |
Definition at line 424 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_CONDITION |
Definition at line 414 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_DIAGNOSTIC |
| #define V_QUAD_DISTORTION |
Definition at line 422 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_EDGE_RATIO |
Definition at line 423 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_JACOBIAN |
Definition at line 415 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_MAX_ASPECT_FROBENIUS |
Definition at line 427 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_MAX_EDGE_RATIO |
Definition at line 405 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_MAXIMUM_ANGLE |
Definition at line 412 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_MED_ASPECT_FROBENIUS |
Definition at line 426 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_MINIMUM_ANGLE |
Definition at line 411 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_ODDY |
Definition at line 413 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_RADIUS_RATIO |
Definition at line 425 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_RELATIVE_SIZE_SQUARED |
Definition at line 419 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_ROBINSON ( V_QUAD_MAX_EDGE_RATIO + V_QUAD_SKEW + V_QUAD_TAPER ) |
| #define V_QUAD_SCALED_JACOBIAN |
Definition at line 416 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_SHAPE |
Definition at line 418 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_SHAPE_AND_SIZE |
Definition at line 420 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_SHEAR |
Definition at line 417 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_SHEAR_AND_SIZE |
Definition at line 421 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_SKEW |
Definition at line 406 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_STRETCH |
Definition at line 410 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_TAPER |
Definition at line 407 of file verdict.h.
Referenced by v_quad_quality().
| #define V_QUAD_TRADITIONAL |
| #define V_QUAD_WARPAGE |
Definition at line 408 of file verdict.h.
Referenced by v_quad_quality().
| #define V_TET_ALGEBRAIC ( V_TET_SHAPE + V_TET_RELATIVE_SIZE_SQUARED + V_TET_SHAPE_AND_SIZE ) |
| #define V_TET_ALL |
Definition at line 374 of file verdict.h.
Referenced by moab::VerdictWrapper::all_quality_measures().
| #define V_TET_ASPECT_BETA |
Definition at line 359 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_ASPECT_FROBENIUS |
Definition at line 371 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_ASPECT_GAMMA |
Definition at line 360 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_ASPECT_RATIO |
Definition at line 370 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_COLLAPSE_RATIO |
Definition at line 373 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_CONDITION |
Definition at line 362 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_DIAGNOSTIC |
| #define V_TET_DISTORTION |
Definition at line 368 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_EDGE_RATIO |
| #define V_TET_JACOBIAN |
Definition at line 363 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_MINIMUM_ANGLE |
Definition at line 372 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_RADIUS_RATIO |
Definition at line 358 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_RELATIVE_SIZE_SQUARED |
Definition at line 366 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_SCALED_JACOBIAN |
Definition at line 364 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_SHAPE |
Definition at line 365 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_SHAPE_AND_SIZE |
Definition at line 367 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TET_TRADITIONAL |
| #define V_TET_VOLUME |
Definition at line 361 of file verdict.h.
Referenced by v_tet_quality().
| #define V_TRI_ALGEBRAIC ( V_TRI_SHAPE + V_TRI_SHAPE_AND_SIZE + V_TRI_RELATIVE_SIZE_SQUARED ) |
| #define V_TRI_ALL |
Definition at line 456 of file verdict.h.
Referenced by moab::VerdictWrapper::all_quality_measures().
| #define V_TRI_AREA |
Definition at line 445 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_ASPECT_FROBENIUS |
Definition at line 444 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_CONDITION |
Definition at line 448 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_DIAGNOSTIC |
| #define V_TRI_DISTORTION |
Definition at line 453 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_EDGE_RATIO |
Definition at line 455 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_MAXIMUM_ANGLE |
Definition at line 447 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_MINIMUM_ANGLE |
Definition at line 446 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_RADIUS_RATIO |
Definition at line 454 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_RELATIVE_SIZE_SQUARED |
Definition at line 451 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_SCALED_JACOBIAN |
Definition at line 449 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_SHAPE |
Definition at line 450 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_SHAPE_AND_SIZE |
Definition at line 452 of file verdict.h.
Referenced by v_tri_quality().
| #define V_TRI_TRADITIONAL |
| #define V_WEDGE_VOLUME |
Definition at line 393 of file verdict.h.
Referenced by v_wedge_quality().
| #define VERDICT_DBL_MAX 1.0E+30 |
Definition at line 37 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_DBL_MIN 1.0E-30 |
Definition at line 36 of file verdict.h.
Referenced by condition_comp(), oddy_comp(), safe_ratio(), skew_x(), v_hex_edge_ratio(), v_hex_quality(), v_hex_relative_size_squared(), v_hex_scaled_jacobian(), v_hex_shape(), v_hex_shear(), v_hex_skew(), v_quad_aspect_ratio(), v_quad_condition(), v_quad_edge_ratio(), 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_shear(), 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_edge_ratio(), v_tet_minimum_angle(), v_tet_quality(), v_tet_radius_ratio(), v_tet_relative_size_squared(), v_tet_scaled_jacobian(), v_tet_shape(), v_tri_aspect_ratio(), v_tri_edge_ratio(), v_tri_radius_ratio(), v_tri_scaled_jacobian(), and v_tri_shape().
| #define VERDICT_PI 3.1415926535897932384626 |
Definition at line 38 of file verdict.h.
Referenced by interior_angle(), VerdictVector::interior_angle(), VerdictVector::scale_angle(), v_quad_maximum_angle(), v_quad_minimum_angle(), v_quad_quality(), and VerdictVector::vector_angle().
| #define VERDICT_VERSION 120 |
| typedef int( * ComputeNormal)(double point[3], double normal[3]) |
| typedef double( * VerdictFunction)(int, double[][3]) |
| 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.
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 1371 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 771 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 791 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 2228 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 529 of file V_HexMetric.cpp.
References 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 1382 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 1243 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 643 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 1164 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 1014 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 2651 of file V_HexMetric.cpp.
References calc_hex_efg(), HexMetricVals::condition, condition_comp(), diag_length(), HexMetricVals::diagonal, HexMetricVals::dimension, HexMetricVals::distortion, HexMetricVals::edge_ratio, 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 2090 of file V_HexMetric.cpp.
References 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 1507 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 1931 of file V_HexMetric.cpp.
References 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 2198 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 1733 of file V_HexMetric.cpp.
References 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 2214 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 674 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 753 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 704 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 732 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 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 828 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 1052 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 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 870 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 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(), 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 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 734 of file V_QuadMetric.cpp.
References moab::angle(), 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 788 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 1258 of file V_QuadMetric.cpp.
References QuadMetricVals::area, QuadMetricVals::aspect_ratio, QuadMetricVals::condition, determinant(), QuadMetricVals::distortion, QuadMetricVals::edge_ratio, 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 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 988 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 888 of file V_QuadMetric.cpp.
References 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 944 of file V_QuadMetric.cpp.
References 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 1020 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 929 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 1037 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 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 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().
{
verdict_hex_size = 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().
{
verdict_quad_size = size;
}
| C_FUNC_DEF void v_set_tet_size | ( | double | size | ) |
Sets average size (volume) of tet, needed for v_tet_relative_size(...)
set the average volume of a tet
Definition at line 37 of file V_TetMetric.cpp.
References size, and verdict_tet_size.
Referenced by moab::VerdictWrapper::set_size().
{
verdict_tet_size = size;
}
| C_FUNC_DEF void v_set_tri_normal_func | ( | ComputeNormal | func | ) |
Sets fuction pointer to calculate tri normal wrt surface.
Definition at line 63 of file V_TriMetric.cpp.
References compute_normal.
{
compute_normal = func;
}
| C_FUNC_DEF void v_set_tri_size | ( | double | size | ) |
Sets average size (area) of tri, needed for v_tri_relative_size(...)
sets the average area of a tri
Definition at line 58 of file V_TriMetric.cpp.
References size, and verdict_tri_size.
Referenced by moab::VerdictWrapper::set_size().
{
verdict_tri_size = size;
}
| C_FUNC_DEF double v_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(), 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 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 VerdictVector::length(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_quality().
{
// 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 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 765 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 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 VerdictVector::length(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tet_quality().
{
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 862 of file V_TetMetric.cpp.
References TetMetricVals::aspect_beta, TetMetricVals::aspect_frobenius, TetMetricVals::aspect_gamma, TetMetricVals::aspect_ratio, TetMetricVals::collapse_ratio, TetMetricVals::condition, TetMetricVals::distortion, 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, 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 727 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 752 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 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 421 of file V_TriMetric.cpp.
References compute_normal, VerdictVector::length(), VERDICT_DBL_MAX, VERDICT_MIN, VerdictVector::x(), VerdictVector::y(), and VerdictVector::z().
Referenced by moab::VerdictWrapper::quality_measure(), v_quad_condition(), and v_tri_shape().
{
static const double rt3 = sqrt( 3.0 );
VerdictVector v1( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
coordinates[1][2] - coordinates[0][2] );
VerdictVector v2( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
coordinates[2][2] - coordinates[0][2] );
VerdictVector tri_normal = v1 * v2;
double areax2 = tri_normal.length();
if( areax2 == 0.0 ) return (double)VERDICT_DBL_MAX;
double condition = (double)( ( ( v1 % v1 ) + ( v2 % v2 ) - ( v1 % v2 ) ) / ( areax2 * rt3 ) );
// check for inverted if we have access to the normal
if( compute_normal )
{
// center of tri
double point[3], surf_normal[3];
point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
// dot product
compute_normal( point, surf_normal );
if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
0 )
return (double)VERDICT_DBL_MAX;
}
return (double)VERDICT_MIN( condition, VERDICT_DBL_MAX );
}
| C_FUNC_DEF double v_tri_distortion | ( | int | num_nodes, |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
{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 588 of file V_TriMetric.cpp.
References GaussIntegration::calculate_derivative_at_nodes_2d_tri(), GaussIntegration::calculate_shape_function_2d_tri(), dot_product(), GaussIntegration::get_shape_func(), GaussIntegration::initialize(), VerdictVector::length(), maxNumberNodes, maxTotalNumberGaussPoints, VerdictVector::normalize(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tri_quality().
{
double distortion;
int total_number_of_gauss_points = 0;
VerdictVector aa, bb, cc, normal_at_point, xin;
double element_area = 0.;
aa.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
coordinates[1][2] - coordinates[0][2] );
bb.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
coordinates[2][2] - coordinates[0][2] );
VerdictVector tri_normal = aa * bb;
int number_of_gauss_points = 0;
if( num_nodes == 3 )
{
distortion = 1.0;
return (double)distortion;
}
else if( num_nodes == 6 )
{
total_number_of_gauss_points = 6;
number_of_gauss_points = 6;
}
distortion = VERDICT_DBL_MAX;
double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
double weight[maxTotalNumberGaussPoints];
// create an object of GaussIntegration
int number_dims = 2;
int is_tri = 1;
GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
GaussIntegration::calculate_shape_function_2d_tri();
GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], weight );
// calculate element area
int ife, ja;
for( ife = 0; ife < total_number_of_gauss_points; ife++ )
{
aa.set( 0.0, 0.0, 0.0 );
bb.set( 0.0, 0.0, 0.0 );
for( ja = 0; ja < num_nodes; ja++ )
{
xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
aa += dndy1[ife][ja] * xin;
bb += dndy2[ife][ja] * xin;
}
normal_at_point = aa * bb;
double jacobian = normal_at_point.length();
element_area += weight[ife] * jacobian;
}
element_area *= 0.8660254;
double dndy1_at_node[maxNumberNodes][maxNumberNodes];
double dndy2_at_node[maxNumberNodes][maxNumberNodes];
GaussIntegration::calculate_derivative_at_nodes_2d_tri( dndy1_at_node, dndy2_at_node );
VerdictVector normal_at_nodes[7];
// evaluate normal at nodes and distortion values at nodes
int jai = 0;
for( ja = 0; ja < num_nodes; ja++ )
{
aa.set( 0.0, 0.0, 0.0 );
bb.set( 0.0, 0.0, 0.0 );
for( jai = 0; jai < num_nodes; jai++ )
{
xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
aa += dndy1_at_node[ja][jai] * xin;
bb += dndy2_at_node[ja][jai] * xin;
}
normal_at_nodes[ja] = aa * bb;
normal_at_nodes[ja].normalize();
}
// determine if element is flat
bool flat_element = true;
double dot_product;
for( ja = 0; ja < num_nodes; ja++ )
{
dot_product = normal_at_nodes[0] % normal_at_nodes[ja];
if( fabs( dot_product ) < 0.99 )
{
flat_element = false;
break;
}
}
// take into consideration of the thickness of the element
double thickness, thickness_gauss;
double distrt;
// get_tri_thickness(tri, element_area, thickness );
thickness = 0.001 * sqrt( element_area );
// set thickness gauss point location
double zl = 0.5773502691896;
if( flat_element ) zl = 0.0;
int no_gauss_pts_z = ( flat_element ) ? 1 : 2;
double thickness_z;
// loop on integration points
int igz;
for( ife = 0; ife < total_number_of_gauss_points; ife++ )
{
// loop on the thickness direction gauss points
for( igz = 0; igz < no_gauss_pts_z; igz++ )
{
zl = -zl;
thickness_z = zl * thickness / 2.0;
aa.set( 0.0, 0.0, 0.0 );
bb.set( 0.0, 0.0, 0.0 );
cc.set( 0.0, 0.0, 0.0 );
for( ja = 0; ja < num_nodes; ja++ )
{
xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
xin += thickness_z * normal_at_nodes[ja];
aa += dndy1[ife][ja] * xin;
bb += dndy2[ife][ja] * xin;
thickness_gauss = shape_function[ife][ja] * thickness / 2.0;
cc += thickness_gauss * normal_at_nodes[ja];
}
normal_at_point = aa * bb;
distrt = cc % normal_at_point;
if( distrt < distortion ) distortion = distrt;
}
}
// loop through nodal points
for( ja = 0; ja < num_nodes; ja++ )
{
for( igz = 0; igz < no_gauss_pts_z; igz++ )
{
zl = -zl;
thickness_z = zl * thickness / 2.0;
aa.set( 0.0, 0.0, 0.0 );
bb.set( 0.0, 0.0, 0.0 );
cc.set( 0.0, 0.0, 0.0 );
for( jai = 0; jai < num_nodes; jai++ )
{
xin.set( coordinates[jai][0], coordinates[jai][1], coordinates[jai][2] );
xin += thickness_z * normal_at_nodes[ja];
aa += dndy1_at_node[ja][jai] * xin;
bb += dndy2_at_node[ja][jai] * xin;
if( jai == ja )
thickness_gauss = thickness / 2.0;
else
thickness_gauss = 0.;
cc += thickness_gauss * normal_at_nodes[jai];
}
}
normal_at_point = aa * bb;
double sign_jacobian = ( tri_normal % normal_at_point ) > 0 ? 1. : -1.;
distrt = sign_jacobian * ( cc % normal_at_point );
if( distrt < distortion ) distortion = distrt;
}
if( element_area * thickness != 0 )
distortion *= 1. / ( element_area * thickness );
else
distortion *= 1.;
if( distortion > 0 ) return (double)VERDICT_MIN( distortion, VERDICT_DBL_MAX );
return (double)VERDICT_MAX( distortion, -VERDICT_DBL_MAX );
}
| C_FUNC_DEF double v_tri_edge_ratio | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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 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 361 of file V_TriMetric.cpp.
References VerdictVector::interior_angle(), VerdictVector::length_squared(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), v_quad_maximum_angle(), and v_quad_quality().
{
// vectors for all the sides
VerdictVector sides[4];
sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
coordinates[1][2] - coordinates[0][2] );
sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
coordinates[2][2] - coordinates[1][2] );
sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
coordinates[2][2] - coordinates[0][2] );
// in case we need to find the interior angle
// between sides 0 and 1
sides[3] = -sides[1];
// calculate the lengths squared of the sides
double sides_lengths[3];
sides_lengths[0] = sides[0].length_squared();
sides_lengths[1] = sides[1].length_squared();
sides_lengths[2] = sides[2].length_squared();
if( sides_lengths[0] == 0.0 || sides_lengths[1] == 0.0 || sides_lengths[2] == 0.0 )
{
return 0.0;
}
// using the law of sines, we know that the maximum
// angle is opposite of the longest side
// find the longest side
int short_side = 0;
if( sides_lengths[1] > sides_lengths[0] ) short_side = 1;
if( sides_lengths[2] > sides_lengths[short_side] ) short_side = 2;
// from the longest side, calculate the angle of the
// opposite angle
double max_angle;
if( short_side == 0 )
{
max_angle = sides[2].interior_angle( sides[1] );
}
else if( short_side == 1 )
{
max_angle = sides[0].interior_angle( sides[2] );
}
else
{
max_angle = sides[0].interior_angle( sides[3] );
}
if( max_angle > 0 ) return (double)VERDICT_MIN( max_angle, VERDICT_DBL_MAX );
return (double)VERDICT_MAX( max_angle, -VERDICT_DBL_MAX );
}
| C_FUNC_DEF double v_tri_minimum_angle | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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 777 of file V_TriMetric.cpp.
References TriMetricVals::area, TriMetricVals::aspect_frobenius, TriMetricVals::aspect_ratio, compute_normal, TriMetricVals::condition, determinant(), TriMetricVals::distortion, TriMetricVals::edge_ratio, interior_angle(), length(), VerdictVector::length(), TriMetricVals::maximum_angle, TriMetricVals::minimum_angle, TriMetricVals::radius_ratio, TriMetricVals::relative_size_squared, TriMetricVals::scaled_jacobian, VerdictVector::set(), TriMetricVals::shape, TriMetricVals::shape_and_size, size, V_TRI_AREA, V_TRI_ASPECT_FROBENIUS, V_TRI_CONDITION, V_TRI_DISTORTION, v_tri_distortion(), V_TRI_EDGE_RATIO, v_tri_edge_ratio(), v_tri_get_weight(), V_TRI_MAXIMUM_ANGLE, V_TRI_MINIMUM_ANGLE, V_TRI_RADIUS_RATIO, v_tri_radius_ratio(), V_TRI_RELATIVE_SIZE_SQUARED, V_TRI_SCALED_JACOBIAN, V_TRI_SHAPE, V_TRI_SHAPE_AND_SIZE, VERDICT_DBL_MAX, VERDICT_MAX, VERDICT_MIN, VerdictVector::x(), VerdictVector::y(), and VerdictVector::z().
Referenced by moab::VerdictWrapper::all_quality_measures().
{
memset( metric_vals, 0, sizeof( TriMetricVals ) );
// for starts, lets set up some basic and common information
/* node numbers and side numbers used below
2
++
/ \
2 / \ 1
/ \
/ \
0 ---------+ 1
0
*/
// vectors for each side
VerdictVector sides[3];
sides[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
coordinates[1][2] - coordinates[0][2] );
sides[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
coordinates[2][2] - coordinates[1][2] );
sides[2].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
coordinates[2][2] - coordinates[0][2] );
VerdictVector tri_normal = sides[0] * sides[2];
// if we have access to normal information, check to see if the
// element is inverted. If we don't have the normal information
// that we need for this, assume the element is not inverted.
// This flag will be used for condition number, jacobian, shape,
// and size and shape.
bool is_inverted = false;
if( compute_normal )
{
// center of tri
double point[3], surf_normal[3];
point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
// dot product
compute_normal( point, surf_normal );
if( ( tri_normal.x() * surf_normal[0] + tri_normal.y() * surf_normal[1] + tri_normal.z() * surf_normal[2] ) <
0 )
is_inverted = true;
}
// lengths squared of each side
double sides_lengths_squared[3];
sides_lengths_squared[0] = sides[0].length_squared();
sides_lengths_squared[1] = sides[1].length_squared();
sides_lengths_squared[2] = sides[2].length_squared();
// if we are doing angle calcuations
if( metrics_request_flag & ( V_TRI_MINIMUM_ANGLE | V_TRI_MAXIMUM_ANGLE ) )
{
// which is short and long side
int short_side = 0, long_side = 0;
if( sides_lengths_squared[1] < sides_lengths_squared[0] ) short_side = 1;
if( sides_lengths_squared[2] < sides_lengths_squared[short_side] ) short_side = 2;
if( sides_lengths_squared[1] > sides_lengths_squared[0] ) long_side = 1;
if( sides_lengths_squared[2] > sides_lengths_squared[long_side] ) long_side = 2;
// calculate the minimum angle of the tri
if( metrics_request_flag & V_TRI_MINIMUM_ANGLE )
{
if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
{
metric_vals->minimum_angle = 0.0;
}
else if( short_side == 0 )
metric_vals->minimum_angle = (double)sides[2].interior_angle( sides[1] );
else if( short_side == 1 )
metric_vals->minimum_angle = (double)sides[0].interior_angle( sides[2] );
else
metric_vals->minimum_angle = (double)sides[0].interior_angle( -sides[1] );
}
// calculate the maximum angle of the tri
if( metrics_request_flag & V_TRI_MAXIMUM_ANGLE )
{
if( sides_lengths_squared[0] == 0.0 || sides_lengths_squared[1] == 0.0 || sides_lengths_squared[2] == 0.0 )
{
metric_vals->maximum_angle = 0.0;
}
else if( long_side == 0 )
metric_vals->maximum_angle = (double)sides[2].interior_angle( sides[1] );
else if( long_side == 1 )
metric_vals->maximum_angle = (double)sides[0].interior_angle( sides[2] );
else
metric_vals->maximum_angle = (double)sides[0].interior_angle( -sides[1] );
}
}
// calculate the area of the tri
// the following metrics depend on area
if( metrics_request_flag &
( V_TRI_AREA | V_TRI_SCALED_JACOBIAN | V_TRI_SHAPE | V_TRI_RELATIVE_SIZE_SQUARED | V_TRI_SHAPE_AND_SIZE ) )
{
metric_vals->area = (double)( ( sides[0] * sides[2] ).length() * 0.5 );
}
// calculate the aspect ratio
if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS )
{
// sum the lengths squared
double srms = sides_lengths_squared[0] + sides_lengths_squared[1] + sides_lengths_squared[2];
// calculate once and reuse
static const double twoTimesRootOf3 = 2 * sqrt( 3.0 );
double div = ( twoTimesRootOf3 * ( ( sides[0] * sides[2] ).length() ) );
if( div == 0.0 )
metric_vals->aspect_frobenius = (double)VERDICT_DBL_MAX;
else
metric_vals->aspect_frobenius = (double)( srms / div );
}
// calculate the scaled jacobian
if( metrics_request_flag & V_TRI_SCALED_JACOBIAN )
{
// calculate once and reuse
static const double twoOverRootOf3 = 2 / sqrt( 3.0 );
// use the area from above
double tmp = tri_normal.length() * twoOverRootOf3;
// now scale it by the lengths of the sides
double min_scaled_jac = VERDICT_DBL_MAX;
double temp_scaled_jac;
for( int i = 0; i < 3; i++ )
{
if( sides_lengths_squared[i % 3] == 0.0 || sides_lengths_squared[( i + 2 ) % 3] == 0.0 )
temp_scaled_jac = 0.0;
else
temp_scaled_jac =
tmp / sqrt( sides_lengths_squared[i % 3] ) / sqrt( sides_lengths_squared[( i + 2 ) % 3] );
if( temp_scaled_jac < min_scaled_jac ) min_scaled_jac = temp_scaled_jac;
}
// multiply by -1 if the normals are in opposite directions
if( is_inverted )
{
min_scaled_jac *= -1;
}
metric_vals->scaled_jacobian = (double)min_scaled_jac;
}
// calculate the condition number
if( metrics_request_flag & V_TRI_CONDITION )
{
// calculate once and reuse
static const double rootOf3 = sqrt( 3.0 );
// if it is inverted, the condition number is considered to be infinity.
if( is_inverted )
{
metric_vals->condition = VERDICT_DBL_MAX;
}
else
{
double area2x = ( sides[0] * sides[2] ).length();
if( area2x == 0.0 )
metric_vals->condition = (double)( VERDICT_DBL_MAX );
else
metric_vals->condition = (double)( ( sides[0] % sides[0] + sides[2] % sides[2] - sides[0] % sides[2] ) /
( area2x * rootOf3 ) );
}
}
// calculate the shape
if( metrics_request_flag & V_TRI_SHAPE || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
{
// if element is inverted, shape is zero. We don't need to
// calculate anything.
if( is_inverted )
{
metric_vals->shape = 0.0;
}
else
{ // otherwise, we calculate the shape
// calculate once and reuse
static const double rootOf3 = sqrt( 3.0 );
// reuse area from before
double area2x = metric_vals->area * 2;
// dot products
double dots[3] = { sides[0] % sides[0], sides[2] % sides[2], sides[0] % sides[2] };
// add the dots
double sum_dots = dots[0] + dots[1] - dots[2];
// then the finale
if( sum_dots == 0.0 )
metric_vals->shape = 0.0;
else
metric_vals->shape = (double)( rootOf3 * area2x / sum_dots );
}
}
// calculate relative size squared
if( metrics_request_flag & V_TRI_RELATIVE_SIZE_SQUARED || metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
{
// get weights
double w11, w21, w12, w22;
v_tri_get_weight( w11, w21, w12, w22 );
// get the determinant
double detw = determinant( w11, w21, w12, w22 );
// use the area from above and divide with the determinant
if( metric_vals->area == 0.0 || detw == 0.0 )
metric_vals->relative_size_squared = 0.0;
else
{
double size = metric_vals->area * 2.0 / detw;
// square the size
size *= size;
// value ranges between 0 to 1
metric_vals->relative_size_squared = (double)VERDICT_MIN( size, 1.0 / size );
}
}
// calculate shape and size
if( metrics_request_flag & V_TRI_SHAPE_AND_SIZE )
{
metric_vals->shape_and_size = metric_vals->relative_size_squared * metric_vals->shape;
}
// calculate distortion
if( metrics_request_flag & V_TRI_DISTORTION ) metric_vals->distortion = v_tri_distortion( num_nodes, coordinates );
// take care of any over-flow problems
if( metric_vals->aspect_frobenius > 0 )
metric_vals->aspect_frobenius = (double)VERDICT_MIN( metric_vals->aspect_frobenius, VERDICT_DBL_MAX );
else
metric_vals->aspect_frobenius = (double)VERDICT_MAX( metric_vals->aspect_frobenius, -VERDICT_DBL_MAX );
if( metric_vals->area > 0 )
metric_vals->area = (double)VERDICT_MIN( metric_vals->area, VERDICT_DBL_MAX );
else
metric_vals->area = (double)VERDICT_MAX( metric_vals->area, -VERDICT_DBL_MAX );
if( metric_vals->minimum_angle > 0 )
metric_vals->minimum_angle = (double)VERDICT_MIN( metric_vals->minimum_angle, VERDICT_DBL_MAX );
else
metric_vals->minimum_angle = (double)VERDICT_MAX( metric_vals->minimum_angle, -VERDICT_DBL_MAX );
if( metric_vals->maximum_angle > 0 )
metric_vals->maximum_angle = (double)VERDICT_MIN( metric_vals->maximum_angle, VERDICT_DBL_MAX );
else
metric_vals->maximum_angle = (double)VERDICT_MAX( metric_vals->maximum_angle, -VERDICT_DBL_MAX );
if( metric_vals->condition > 0 )
metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
else
metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
if( metric_vals->shape > 0 )
metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
else
metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
if( metric_vals->scaled_jacobian > 0 )
metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
else
metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
if( metric_vals->relative_size_squared > 0 )
metric_vals->relative_size_squared = (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
else
metric_vals->relative_size_squared =
(double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
if( metric_vals->shape_and_size > 0 )
metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
else
metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
if( metric_vals->distortion > 0 )
metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
else
metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
if( metrics_request_flag & V_TRI_EDGE_RATIO )
{
metric_vals->edge_ratio = v_tri_edge_ratio( 3, coordinates );
}
if( metrics_request_flag & V_TRI_RADIUS_RATIO )
{
metric_vals->radius_ratio = v_tri_radius_ratio( 3, coordinates );
}
if( metrics_request_flag & V_TRI_ASPECT_FROBENIUS ) // there is no V_TRI_ASPECT_RATIO !
{
metric_vals->aspect_ratio = v_tri_radius_ratio( 3, coordinates );
}
}
| C_FUNC_DEF double v_tri_radius_ratio | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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 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 534 of file V_TriMetric.cpp.
References determinant(), VerdictVector::length(), VerdictVector::set(), size, v_tri_get_weight(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tri_shape_and_size().
{
double w11, w21, w12, w22;
VerdictVector xxi, xet, tri_normal;
v_tri_get_weight( w11, w21, w12, w22 );
double detw = determinant( w11, w21, w12, w22 );
if( detw == 0.0 ) return 0.0;
xxi.set( coordinates[0][0] - coordinates[1][0], coordinates[0][1] - coordinates[1][1],
coordinates[0][2] - coordinates[1][2] );
xet.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
coordinates[0][2] - coordinates[2][2] );
tri_normal = xxi * xet;
double deta = tri_normal.length();
if( deta == 0.0 || detw == 0.0 ) return 0.0;
double size = pow( deta / detw, 2 );
double rel_size = VERDICT_MIN( size, 1.0 / size );
if( rel_size > 0 ) return (double)VERDICT_MIN( rel_size, VERDICT_DBL_MAX );
return (double)VERDICT_MAX( rel_size, -VERDICT_DBL_MAX );
}
| C_FUNC_DEF double v_tri_scaled_jacobian | ( | int | , |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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 461 of file V_TriMetric.cpp.
References compute_normal, moab::cross(), moab::GeomUtil::first(), length(), VerdictVector::length(), VerdictVector::set(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, VERDICT_MIN, VerdictVector::x(), VerdictVector::y(), and VerdictVector::z().
Referenced by moab::VerdictWrapper::quality_measure(), v_quad_quality(), and v_quad_scaled_jacobian().
{
static const double detw = 2. / sqrt( 3.0 );
VerdictVector first, second;
double jacobian;
VerdictVector edge[3];
edge[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
coordinates[1][2] - coordinates[0][2] );
edge[1].set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
coordinates[2][2] - coordinates[0][2] );
edge[2].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
coordinates[2][2] - coordinates[1][2] );
first = edge[1] - edge[0];
second = edge[2] - edge[0];
VerdictVector cross = first * second;
jacobian = cross.length();
double max_edge_length_product;
max_edge_length_product =
VERDICT_MAX( edge[0].length() * edge[1].length(),
VERDICT_MAX( edge[1].length() * edge[2].length(), edge[0].length() * edge[2].length() ) );
if( max_edge_length_product < VERDICT_DBL_MIN ) return (double)0.0;
jacobian *= detw;
jacobian /= max_edge_length_product;
if( compute_normal )
{
// center of tri
double point[3], surf_normal[3];
point[0] = ( coordinates[0][0] + coordinates[1][0] + coordinates[2][0] ) / 3;
point[1] = ( coordinates[0][1] + coordinates[1][1] + coordinates[2][1] ) / 3;
point[2] = ( coordinates[0][2] + coordinates[1][2] + coordinates[2][2] ) / 3;
// dot product
compute_normal( point, surf_normal );
if( ( cross.x() * surf_normal[0] + cross.y() * surf_normal[1] + cross.z() * surf_normal[2] ) < 0 )
jacobian *= -1;
}
if( jacobian > 0 ) return (double)VERDICT_MIN( jacobian, VERDICT_DBL_MAX );
return (double)VERDICT_MAX( jacobian, -VERDICT_DBL_MAX );
}
| C_FUNC_DEF double v_tri_shape | ( | int | num_nodes, |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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 515 of file V_TriMetric.cpp.
References v_tri_condition(), VERDICT_DBL_MAX, VERDICT_DBL_MIN, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure(), and v_tri_shape_and_size().
{
double condition = v_tri_condition( num_nodes, coordinates );
double shape;
if( condition <= VERDICT_DBL_MIN )
shape = VERDICT_DBL_MAX;
else
shape = ( 1 / condition );
if( shape > 0 ) return (double)VERDICT_MIN( shape, VERDICT_DBL_MAX );
return (double)VERDICT_MAX( shape, -VERDICT_DBL_MAX );
}
| C_FUNC_DEF double v_tri_shape_and_size | ( | int | num_nodes, |
| double | coordinates[][3] | ||
| ) |
Calculates tri metric.
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 570 of file V_TriMetric.cpp.
References size, v_tri_relative_size_squared(), v_tri_shape(), VERDICT_DBL_MAX, VERDICT_MAX, and VERDICT_MIN.
Referenced by moab::VerdictWrapper::quality_measure().
{
double size, shape;
size = v_tri_relative_size_squared( num_nodes, coordinates );
shape = v_tri_shape( num_nodes, coordinates );
double shape_and_size = size * shape;
if( shape_and_size > 0 ) return (double)VERDICT_MIN( shape_and_size, VERDICT_DBL_MAX );
return (double)VERDICT_MAX( shape_and_size, -VERDICT_DBL_MAX );
}
| 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;
}
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.