Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
Macros | Functions | Variables
V_TetMetric.cpp File Reference
#include "moab/verdict.h"
#include "verdict_defines.hpp"
#include "VerdictVector.hpp"
#include "V_GaussIntegration.hpp"
#include <memory.h>
+ Include dependency graph for V_TetMetric.cpp:

Go to the source code of this file.

Macros

#define VERDICT_EXPORTS
 

Functions

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

Variables

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

Macro Definition Documentation

◆ VERDICT_EXPORTS

#define VERDICT_EXPORTS

Definition at line 23 of file V_TetMetric.cpp.

Function Documentation

◆ get_weight()

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

get the weights based on the average size of a tet

Definition at line 46 of file V_TetMetric.cpp.

47 {
48  static const double rt3 = sqrt( 3.0 );
49  static const double root_of_2 = sqrt( 2.0 );
50 
51  w1.set( 1, 0, 0 );
52  w2.set( 0.5, 0.5 * rt3, 0 );
53  w3.set( 0.5, rt3 / 6.0, root_of_2 / rt3 );
54 
55  double scale = pow( 6. * verdict_tet_size / determinant( w1, w2, w3 ), 0.3333333333333 );
56 
57  w1 *= scale;
58  w2 *= scale;
59  w3 *= scale;
60 
61  return 1;
62 }

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

Referenced by v_tet_quality(), and v_tet_relative_size_squared().

◆ v_set_tet_size()

C_FUNC_DEF void v_set_tet_size ( double  size)

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

set the average volume of a tet

Definition at line 37 of file V_TetMetric.cpp.

38 {
39  verdict_tet_size = size;
40 }

References size, and verdict_tet_size.

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

◆ v_tet_aspect_beta()

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

Calculates the radius ratio metric of a positively oriented tet.

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

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

Definition at line 265 of file V_TetMetric.cpp.

266 {
267 
268  // Determine side vectors
269  VerdictVector side[6];
270 
271  side[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
272  coordinates[1][2] - coordinates[0][2] );
273 
274  side[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
275  coordinates[2][2] - coordinates[1][2] );
276 
277  side[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
278  coordinates[0][2] - coordinates[2][2] );
279 
280  side[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
281  coordinates[3][2] - coordinates[0][2] );
282 
283  side[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
284  coordinates[3][2] - coordinates[1][2] );
285 
286  side[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
287  coordinates[3][2] - coordinates[2][2] );
288 
289  VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
290  side[2].length_squared() * ( side[3] * side[0] ) +
291  side[0].length_squared() * ( side[3] * side[2] );
292 
293  double area_sum;
294  area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
295  ( side[3] * side[2] ).length() ) *
296  0.5;
297 
298  double volume = v_tet_volume( 4, coordinates );
299 
300  if( volume < VERDICT_DBL_MIN )
301  return (double)VERDICT_DBL_MAX;
302  else
303  {
304  double radius_ratio;
305  radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );
306 
307  return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
308  }
309 }

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

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

◆ v_tet_aspect_frobenius()

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

Calculates tet aspect frobenius metric.

The aspect frobenius of a tet

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

Definition at line 426 of file V_TetMetric.cpp.

427 {
428  static const double normal_exp = 1. / 3.;
429 
430  VerdictVector ab, ac, ad;
431 
432  ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
433  coordinates[1][2] - coordinates[0][2] );
434 
435  ac.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
436  coordinates[2][2] - coordinates[0][2] );
437 
438  ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
439  coordinates[3][2] - coordinates[0][2] );
440 
441  double denominator = ab % ( ac * ad );
442  denominator *= denominator;
443  denominator *= 2.;
444  denominator = 3. * pow( denominator, normal_exp );
445 
446  if( denominator < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
447 
448  double u[3];
449  ab.get_xyz( u );
450  double v[3];
451  ac.get_xyz( v );
452  double w[3];
453  ad.get_xyz( w );
454 
455  double numerator = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
456  numerator += v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
457  numerator += w[0] * w[0] + w[1] * w[1] + w[2] * w[2];
458  numerator *= 1.5;
459  numerator -= v[0] * u[0] + v[1] * u[1] + v[2] * u[2];
460  numerator -= w[0] * u[0] + w[1] * u[1] + w[2] * u[2];
461  numerator -= w[0] * v[0] + w[1] * v[1] + w[2] * v[2];
462 
463  double aspect_frobenius = numerator / denominator;
464 
465  if( aspect_frobenius > 0 ) return (double)VERDICT_MIN( aspect_frobenius, VERDICT_DBL_MAX );
466  return (double)VERDICT_MAX( aspect_frobenius, -VERDICT_DBL_MAX );
467 }

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

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

◆ v_tet_aspect_gamma()

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

Calculates tet aspect gamma metric.

the aspect gamma of a tet

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

Definition at line 381 of file V_TetMetric.cpp.

382 {
383 
384  // Determine side vectors
385  VerdictVector side0, side1, side2, side3, side4, side5;
386 
387  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
388  coordinates[1][2] - coordinates[0][2] );
389 
390  side1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
391  coordinates[2][2] - coordinates[1][2] );
392 
393  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
394  coordinates[0][2] - coordinates[2][2] );
395 
396  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
397  coordinates[3][2] - coordinates[0][2] );
398 
399  side4.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
400  coordinates[3][2] - coordinates[1][2] );
401 
402  side5.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
403  coordinates[3][2] - coordinates[2][2] );
404 
405  double volume = fabs( v_tet_volume( 4, coordinates ) );
406 
407  if( volume < VERDICT_DBL_MIN )
408  return (double)VERDICT_DBL_MAX;
409  else
410  {
411  double srms = sqrt( ( side0.length_squared() + side1.length_squared() + side2.length_squared() +
412  side3.length_squared() + side4.length_squared() + side5.length_squared() ) /
413  6.0 );
414 
415  double aspect_ratio_gamma = pow( srms, 3 ) / ( 8.48528137423857 * volume );
416  return (double)aspect_ratio_gamma;
417  }
418 }

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

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

◆ v_tet_aspect_ratio()

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

Calculates tet aspect ratio metric.

The aspect ratio of a tet

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

Definition at line 318 of file V_TetMetric.cpp.

319 {
320  static const double normal_coeff = sqrt( 6. ) / 12.;
321 
322  // Determine side vectors
323  VerdictVector ab, bc, ac, ad, bd, cd;
324 
325  ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
326  coordinates[1][2] - coordinates[0][2] );
327 
328  ac.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
329  coordinates[2][2] - coordinates[0][2] );
330 
331  ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
332  coordinates[3][2] - coordinates[0][2] );
333 
334  double detTet = ab % ( ac * ad );
335 
336  if( detTet < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
337 
338  bc.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
339  coordinates[2][2] - coordinates[1][2] );
340 
341  bd.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
342  coordinates[3][2] - coordinates[1][2] );
343 
344  cd.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
345  coordinates[3][2] - coordinates[2][2] );
346 
347  double ab2 = ab.length_squared();
348  double bc2 = bc.length_squared();
349  double ac2 = ac.length_squared();
350  double ad2 = ad.length_squared();
351  double bd2 = bd.length_squared();
352  double cd2 = cd.length_squared();
353 
354  double A = ab2 > bc2 ? ab2 : bc2;
355  double B = ac2 > ad2 ? ac2 : ad2;
356  double C = bd2 > cd2 ? bd2 : cd2;
357  double D = A > B ? A : B;
358  double hm = D > C ? sqrt( D ) : sqrt( C );
359 
360  bd = ab * bc;
361  A = bd.length();
362  bd = ab * ad;
363  B = bd.length();
364  bd = ac * ad;
365  C = bd.length();
366  bd = bc * cd;
367  D = bd.length();
368 
369  double aspect_ratio;
370  aspect_ratio = normal_coeff * hm * ( A + B + C + D ) / fabs( detTet );
371 
372  if( aspect_ratio > 0 ) return (double)VERDICT_MIN( aspect_ratio, VERDICT_DBL_MAX );
373  return (double)VERDICT_MAX( aspect_ratio, -VERDICT_DBL_MAX );
374 }

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

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

◆ v_tet_collapse_ratio()

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

Calculates tet collapse ratio metric.

The collapse ratio of a tet

Definition at line 526 of file V_TetMetric.cpp.

527 {
528  // Determine side vectors
529  VerdictVector e01, e02, e03, e12, e13, e23;
530 
531  e01.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
532  coordinates[1][2] - coordinates[0][2] );
533 
534  e02.set( coordinates[2][0] - coordinates[0][0], coordinates[2][1] - coordinates[0][1],
535  coordinates[2][2] - coordinates[0][2] );
536 
537  e03.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
538  coordinates[3][2] - coordinates[0][2] );
539 
540  e12.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
541  coordinates[2][2] - coordinates[1][2] );
542 
543  e13.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
544  coordinates[3][2] - coordinates[1][2] );
545 
546  e23.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
547  coordinates[3][2] - coordinates[2][2] );
548 
549  double l[6];
550  l[0] = e01.length();
551  l[1] = e02.length();
552  l[2] = e03.length();
553  l[3] = e12.length();
554  l[4] = e13.length();
555  l[5] = e23.length();
556 
557  // Find longest edge for each bounding triangle of tetrahedron
558  double l012 = l[4] > l[0] ? l[4] : l[0];
559  l012 = l[1] > l012 ? l[1] : l012;
560  double l031 = l[0] > l[2] ? l[0] : l[2];
561  l031 = l[3] > l031 ? l[3] : l031;
562  double l023 = l[2] > l[1] ? l[2] : l[1];
563  l023 = l[5] > l023 ? l[5] : l023;
564  double l132 = l[4] > l[3] ? l[4] : l[3];
565  l132 = l[5] > l132 ? l[5] : l132;
566 
567  // Compute collapse ratio for each vertex/triangle pair
568  VerdictVector N;
569  double h, magN;
570  double cr;
571  double crMin;
572 
573  N = e01 * e02;
574  magN = N.length();
575  h = ( e03 % N ) / magN; // height of vertex 3 above 0-1-2
576  crMin = h / l012; // ratio of height to longest edge of 0-1-2
577 
578  N = e03 * e01;
579  magN = N.length();
580  h = ( e02 % N ) / magN; // height of vertex 2 above 0-3-1
581  cr = h / l031; // ratio of height to longest edge of 0-3-1
582  if( cr < crMin ) crMin = cr;
583 
584  N = e02 * e03;
585  magN = N.length();
586  h = ( e01 % N ) / magN; // height of vertex 1 above 0-2-3
587  cr = h / l023; // ratio of height to longest edge of 0-2-3
588  if( cr < crMin ) crMin = cr;
589 
590  N = e12 * e13;
591  magN = N.length();
592  h = ( e01 % N ) / magN; // height of vertex 0 above 1-3-2
593  cr = h / l132; // ratio of height to longest edge of 1-3-2
594  if( cr < crMin ) crMin = cr;
595 
596  if( crMin < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
597  if( crMin > 0 ) return (double)VERDICT_MIN( crMin, VERDICT_DBL_MAX );
598  return (double)VERDICT_MAX( crMin, -VERDICT_DBL_MAX );
599 }

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

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

◆ v_tet_condition()

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

Calculates tet condition metric.

the condition of a tet

condition number of the jacobian matrix at any corner

Definition at line 629 of file V_TetMetric.cpp.

630 {
631 
632  double condition, term1, term2, det;
633  double rt3 = sqrt( 3.0 );
634  double rt6 = sqrt( 6.0 );
635 
636  VerdictVector side0, side2, side3;
637 
638  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
639  coordinates[1][2] - coordinates[0][2] );
640 
641  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
642  coordinates[0][2] - coordinates[2][2] );
643 
644  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
645  coordinates[3][2] - coordinates[0][2] );
646 
647  VerdictVector c_1, c_2, c_3;
648 
649  c_1 = side0;
650  c_2 = ( -2 * side2 - side0 ) / rt3;
651  c_3 = ( 3 * side3 + side2 - side0 ) / rt6;
652 
653  term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
654  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 );
655  det = c_1 % ( c_2 * c_3 );
656 
657  if( det <= VERDICT_DBL_MIN )
658  return VERDICT_DBL_MAX;
659  else
660  condition = sqrt( term1 * term2 ) / ( 3.0 * det );
661 
662  return (double)condition;
663 }

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

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

◆ v_tet_distortion()

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

Calculates tet distortion metric.

the distortion of a tet

Definition at line 765 of file V_TetMetric.cpp.

766 {
767 
768  double distortion = VERDICT_DBL_MAX;
769  int number_of_gauss_points = 0;
770  if( num_nodes == 4 )
771  // for linear tet, the distortion is always 1 because
772  // straight edge tets are the target shape for tet
773  return 1.0;
774 
775  else if( num_nodes == 10 )
776  // use four integration points for quadratic tet
777  number_of_gauss_points = 4;
778 
779  int number_dims = 3;
780  int total_number_of_gauss_points = number_of_gauss_points;
781  // use is_tri=1 to indicate this is for tet in 3D
782  int is_tri = 1;
783 
784  double shape_function[maxTotalNumberGaussPoints][maxNumberNodes];
785  double dndy1[maxTotalNumberGaussPoints][maxNumberNodes];
786  double dndy2[maxTotalNumberGaussPoints][maxNumberNodes];
787  double dndy3[maxTotalNumberGaussPoints][maxNumberNodes];
788  double weight[maxTotalNumberGaussPoints];
789 
790  // create an object of GaussIntegration for tet
791  GaussIntegration::initialize( number_of_gauss_points, num_nodes, number_dims, is_tri );
792  GaussIntegration::calculate_shape_function_3d_tet();
793  GaussIntegration::get_shape_func( shape_function[0], dndy1[0], dndy2[0], dndy3[0], weight );
794 
795  // vector xxi is the derivative vector of coordinates w.r.t local xi coordinate in the
796  // computation space
797  // vector xet is the derivative vector of coordinates w.r.t local et coordinate in the
798  // computation space
799  // vector xze is the derivative vector of coordinates w.r.t local ze coordinate in the
800  // computation space
801  VerdictVector xxi, xet, xze, xin;
802 
803  double jacobian, minimum_jacobian;
804  double element_volume = 0.0;
805  minimum_jacobian = VERDICT_DBL_MAX;
806 
807  // calculate element volume
808  int ife, ja;
809  for( ife = 0; ife < total_number_of_gauss_points; ife++ )
810  {
811  xxi.set( 0.0, 0.0, 0.0 );
812  xet.set( 0.0, 0.0, 0.0 );
813  xze.set( 0.0, 0.0, 0.0 );
814 
815  for( ja = 0; ja < num_nodes; ja++ )
816  {
817  xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
818  xxi += dndy1[ife][ja] * xin;
819  xet += dndy2[ife][ja] * xin;
820  xze += dndy3[ife][ja] * xin;
821  }
822 
823  // determinant
824  jacobian = xxi % ( xet * xze );
825  if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
826 
827  element_volume += weight[ife] * jacobian;
828  } // element_volume is 6 times the actual volume
829 
830  // loop through all nodes
831  double dndy1_at_node[maxNumberNodes][maxNumberNodes];
832  double dndy2_at_node[maxNumberNodes][maxNumberNodes];
833  double dndy3_at_node[maxNumberNodes][maxNumberNodes];
834 
835  GaussIntegration::calculate_derivative_at_nodes_3d_tet( dndy1_at_node, dndy2_at_node, dndy3_at_node );
836  int node_id;
837  for( node_id = 0; node_id < num_nodes; node_id++ )
838  {
839  xxi.set( 0.0, 0.0, 0.0 );
840  xet.set( 0.0, 0.0, 0.0 );
841  xze.set( 0.0, 0.0, 0.0 );
842 
843  for( ja = 0; ja < num_nodes; ja++ )
844  {
845  xin.set( coordinates[ja][0], coordinates[ja][1], coordinates[ja][2] );
846  xxi += dndy1_at_node[node_id][ja] * xin;
847  xet += dndy2_at_node[node_id][ja] * xin;
848  xze += dndy3_at_node[node_id][ja] * xin;
849  }
850 
851  jacobian = xxi % ( xet * xze );
852  if( minimum_jacobian > jacobian ) minimum_jacobian = jacobian;
853  }
854  distortion = minimum_jacobian / element_volume;
855 
856  return (double)distortion;
857 }

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

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

◆ v_tet_edge_ratio()

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

Calculates tet edge ratio metric.

the edge ratio of a tet

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

Definition at line 71 of file V_TetMetric.cpp.

72 {
73  VerdictVector a, b, c, d, e, f;
74 
75  a.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
76  coordinates[1][2] - coordinates[0][2] );
77 
78  b.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
79  coordinates[2][2] - coordinates[1][2] );
80 
81  c.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
82  coordinates[0][2] - coordinates[2][2] );
83 
84  d.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
85  coordinates[3][2] - coordinates[0][2] );
86 
87  e.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
88  coordinates[3][2] - coordinates[1][2] );
89 
90  f.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
91  coordinates[3][2] - coordinates[2][2] );
92 
93  double a2 = a.length_squared();
94  double b2 = b.length_squared();
95  double c2 = c.length_squared();
96  double d2 = d.length_squared();
97  double e2 = e.length_squared();
98  double f2 = f.length_squared();
99 
100  double m2, M2, mab, mcd, mef, Mab, Mcd, Mef;
101 
102  if( a2 < b2 )
103  {
104  mab = a2;
105  Mab = b2;
106  }
107  else // b2 <= a2
108  {
109  mab = b2;
110  Mab = a2;
111  }
112  if( c2 < d2 )
113  {
114  mcd = c2;
115  Mcd = d2;
116  }
117  else // d2 <= c2
118  {
119  mcd = d2;
120  Mcd = c2;
121  }
122  if( e2 < f2 )
123  {
124  mef = e2;
125  Mef = f2;
126  }
127  else // f2 <= e2
128  {
129  mef = f2;
130  Mef = e2;
131  }
132 
133  m2 = mab < mcd ? mab : mcd;
134  m2 = m2 < mef ? m2 : mef;
135 
136  if( m2 < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
137 
138  M2 = Mab > Mcd ? Mab : Mcd;
139  M2 = M2 > Mef ? M2 : Mef;
140 
141  double edge_ratio = sqrt( M2 / m2 );
142 
143  if( edge_ratio > 0 ) return (double)VERDICT_MIN( edge_ratio, VERDICT_DBL_MAX );
144  return (double)VERDICT_MAX( edge_ratio, -VERDICT_DBL_MAX );
145 }

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

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

◆ v_tet_jacobian()

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

Calculates tet jacobian.

the jacobian of a tet

TODO

Definition at line 670 of file V_TetMetric.cpp.

671 {
672  VerdictVector side0, side2, side3;
673 
674  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
675  coordinates[1][2] - coordinates[0][2] );
676 
677  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
678  coordinates[0][2] - coordinates[2][2] );
679 
680  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
681  coordinates[3][2] - coordinates[0][2] );
682 
683  return (double)( side3 % ( side2 * side0 ) );
684 }

References VerdictVector::set().

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

◆ v_tet_minimum_angle()

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

Calculates tet minimum dihedral angle.

The minimum angle of a tet

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

Definition at line 475 of file V_TetMetric.cpp.

476 {
477  static const double normal_coeff = 180. * .3183098861837906715377675267450287;
478 
479  // Determine side vectors
480  VerdictVector ab, bc, ad, cd;
481 
482  ab.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
483  coordinates[1][2] - coordinates[0][2] );
484 
485  ad.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
486  coordinates[3][2] - coordinates[0][2] );
487 
488  bc.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
489  coordinates[2][2] - coordinates[1][2] );
490 
491  cd.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
492  coordinates[3][2] - coordinates[2][2] );
493 
494  VerdictVector abc = ab * bc;
495  double nabc = abc.length();
496  VerdictVector abd = ab * ad;
497  double nabd = abd.length();
498  VerdictVector acd = ad * cd;
499  double nacd = acd.length();
500  VerdictVector bcd = bc * cd;
501  double nbcd = bcd.length();
502 
503  double alpha = acos( ( abc % abd ) / ( nabc * nabd ) );
504  double beta = acos( ( abc % acd ) / ( nabc * nacd ) );
505  double gamma = acos( ( abc % bcd ) / ( nabc * nbcd ) );
506  double delta = acos( ( abd % acd ) / ( nabd * nacd ) );
507  double epsilon = acos( ( abd % bcd ) / ( nabd * nbcd ) );
508  double zeta = acos( ( acd % bcd ) / ( nacd * nbcd ) );
509 
510  alpha = alpha < beta ? alpha : beta;
511  alpha = alpha < gamma ? alpha : gamma;
512  alpha = alpha < delta ? alpha : delta;
513  alpha = alpha < epsilon ? alpha : epsilon;
514  alpha = alpha < zeta ? alpha : zeta;
515  alpha *= normal_coeff;
516 
517  if( alpha < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
518 
519  if( alpha > 0 ) return (double)VERDICT_MIN( alpha, VERDICT_DBL_MAX );
520  return (double)VERDICT_MAX( alpha, -VERDICT_DBL_MAX );
521 }

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

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

◆ v_tet_quality()

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

Calculates quality metrics for tetrahedral elements.

the quality metrics of a tet

Definition at line 862 of file V_TetMetric.cpp.

866 {
867 
868  memset( metric_vals, 0, sizeof( TetMetricVals ) );
869 
870  /*
871 
872  node numbers and edge numbers below
873 
874 
875 
876  3
877  + edge 0 is node 0 to 1
878  +|+ edge 1 is node 1 to 2
879  3/ | \5 edge 2 is node 0 to 2
880  / 4| \ edge 3 is node 0 to 3
881  0 - -|- + 2 edge 4 is node 1 to 3
882  \ | + edge 5 is node 2 to 3
883  0\ | /1
884  +|/ edge 2 is behind edge 4
885  1
886 
887 
888  */
889 
890  // lets start with making the vectors
891  VerdictVector edges[6];
892  edges[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
893  coordinates[1][2] - coordinates[0][2] );
894 
895  edges[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
896  coordinates[2][2] - coordinates[1][2] );
897 
898  edges[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
899  coordinates[0][2] - coordinates[2][2] );
900 
901  edges[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
902  coordinates[3][2] - coordinates[0][2] );
903 
904  edges[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
905  coordinates[3][2] - coordinates[1][2] );
906 
907  edges[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
908  coordinates[3][2] - coordinates[2][2] );
909 
910  // common numbers
911  static const double root_of_2 = sqrt( 2.0 );
912 
913  // calculate the jacobian
914  static const int do_jacobian = V_TET_JACOBIAN | V_TET_VOLUME | V_TET_ASPECT_BETA | V_TET_ASPECT_GAMMA |
915  V_TET_SHAPE | V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE |
916  V_TET_SCALED_JACOBIAN | V_TET_CONDITION;
917  if( metrics_request_flag & do_jacobian )
918  {
919  metric_vals->jacobian = (double)( edges[3] % ( edges[2] * edges[0] ) );
920  }
921 
922  // calculate the volume
923  if( metrics_request_flag & V_TET_VOLUME )
924  {
925  metric_vals->volume = (double)( metric_vals->jacobian / 6.0 );
926  }
927 
928  // calculate aspect ratio
929  if( metrics_request_flag & V_TET_ASPECT_BETA )
930  {
931  double surface_area = ( ( edges[2] * edges[0] ).length() + ( edges[3] * edges[0] ).length() +
932  ( edges[4] * edges[1] ).length() + ( edges[3] * edges[2] ).length() ) *
933  0.5;
934 
935  VerdictVector numerator = edges[3].length_squared() * ( edges[2] * edges[0] ) +
936  edges[2].length_squared() * ( edges[3] * edges[0] ) +
937  edges[0].length_squared() * ( edges[3] * edges[2] );
938 
939  double volume = metric_vals->jacobian / 6.0;
940 
941  if( volume < VERDICT_DBL_MIN )
942  metric_vals->aspect_beta = (double)( VERDICT_DBL_MAX );
943  else
944  metric_vals->aspect_beta = (double)( numerator.length() * surface_area / ( 108 * volume * volume ) );
945  }
946 
947  // calculate the aspect gamma
948  if( metrics_request_flag & V_TET_ASPECT_GAMMA )
949  {
950  double volume = fabs( metric_vals->jacobian / 6.0 );
951  if( fabs( volume ) < VERDICT_DBL_MIN )
952  metric_vals->aspect_gamma = VERDICT_DBL_MAX;
953  else
954  {
955  double srms = sqrt( ( edges[0].length_squared() + edges[1].length_squared() + edges[2].length_squared() +
956  edges[3].length_squared() + edges[4].length_squared() + edges[5].length_squared() ) /
957  6.0 );
958 
959  // cube the srms
960  srms *= ( srms * srms );
961  metric_vals->aspect_gamma = (double)( srms / ( 8.48528137423857 * volume ) );
962  }
963  }
964 
965  // calculate the shape of the tet
966  if( metrics_request_flag & ( V_TET_SHAPE | V_TET_SHAPE_AND_SIZE ) )
967  {
968  // if the jacobian is non-positive, the shape is 0
969  if( metric_vals->jacobian < VERDICT_DBL_MIN )
970  {
971  metric_vals->shape = (double)0.0;
972  }
973  else
974  {
975  static const double two_thirds = 2.0 / 3.0;
976  double num = 3.0 * pow( root_of_2 * metric_vals->jacobian, two_thirds );
977  double den = 1.5 * ( edges[0] % edges[0] + edges[2] % edges[2] + edges[3] % edges[3] ) -
978  ( edges[0] % -edges[2] + -edges[2] % edges[3] + edges[3] % edges[0] );
979 
980  if( den < VERDICT_DBL_MIN )
981  metric_vals->shape = (double)0.0;
982  else
983  metric_vals->shape = (double)VERDICT_MAX( num / den, 0 );
984  }
985  }
986 
987  // calculate the relative size of the tet
988  if( metrics_request_flag & ( V_TET_RELATIVE_SIZE_SQUARED | V_TET_SHAPE_AND_SIZE ) )
989  {
990  VerdictVector w1, w2, w3;
991  get_weight( w1, w2, w3 );
992  double avg_vol = ( w1 % ( w2 * w3 ) ) / 6;
993 
994  if( avg_vol < VERDICT_DBL_MIN )
995  metric_vals->relative_size_squared = 0.0;
996  else
997  {
998  double tmp = metric_vals->jacobian / ( 6 * avg_vol );
999  if( tmp < VERDICT_DBL_MIN )
1000  metric_vals->relative_size_squared = 0.0;
1001  else
1002  {
1003  tmp *= tmp;
1004  metric_vals->relative_size_squared = (double)VERDICT_MIN( tmp, 1 / tmp );
1005  }
1006  }
1007  }
1008 
1009  // calculate the shape and size
1010  if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
1011  {
1012  metric_vals->shape_and_size = (double)( metric_vals->shape * metric_vals->relative_size_squared );
1013  }
1014 
1015  // calculate the scaled jacobian
1016  if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
1017  {
1018  // find out which node the normalized jacobian can be calculated at
1019  // and it will be the smaller than at other nodes
1020  double length_squared[4] = { edges[0].length_squared() * edges[2].length_squared() * edges[3].length_squared(),
1021  edges[0].length_squared() * edges[1].length_squared() * edges[4].length_squared(),
1022  edges[1].length_squared() * edges[2].length_squared() * edges[5].length_squared(),
1023  edges[3].length_squared() * edges[4].length_squared() *
1024  edges[5].length_squared() };
1025 
1026  int which_node = 0;
1027  if( length_squared[1] > length_squared[which_node] ) which_node = 1;
1028  if( length_squared[2] > length_squared[which_node] ) which_node = 2;
1029  if( length_squared[3] > length_squared[which_node] ) which_node = 3;
1030 
1031  // find the scaled jacobian at this node
1032  double length_product = sqrt( length_squared[which_node] );
1033  if( length_product < fabs( metric_vals->jacobian ) ) length_product = fabs( metric_vals->jacobian );
1034 
1035  if( length_product < VERDICT_DBL_MIN )
1036  metric_vals->scaled_jacobian = (double)VERDICT_DBL_MAX;
1037  else
1038  metric_vals->scaled_jacobian = (double)( root_of_2 * metric_vals->jacobian / length_product );
1039  }
1040 
1041  // calculate the condition number
1042  if( metrics_request_flag & V_TET_CONDITION )
1043  {
1044  static const double root_of_3 = sqrt( 3.0 );
1045  static const double root_of_6 = sqrt( 6.0 );
1046 
1047  VerdictVector c_1, c_2, c_3;
1048  c_1 = edges[0];
1049  c_2 = ( -2 * edges[2] - edges[0] ) / root_of_3;
1050  c_3 = ( 3 * edges[3] + edges[2] - edges[0] ) / root_of_6;
1051 
1052  double term1 = c_1 % c_1 + c_2 % c_2 + c_3 % c_3;
1053  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 );
1054 
1055  double det = c_1 % ( c_2 * c_3 );
1056 
1057  if( det <= VERDICT_DBL_MIN )
1058  metric_vals->condition = (double)VERDICT_DBL_MAX;
1059  else
1060  metric_vals->condition = (double)( sqrt( term1 * term2 ) / ( 3.0 * det ) );
1061  }
1062 
1063  // calculate the distortion
1064  if( metrics_request_flag & V_TET_DISTORTION )
1065  {
1066  metric_vals->distortion = v_tet_distortion( num_nodes, coordinates );
1067  }
1068 
1069  // check for overflow
1070  if( metrics_request_flag & V_TET_ASPECT_BETA )
1071  {
1072  if( metric_vals->aspect_beta > 0 )
1073  metric_vals->aspect_beta = (double)VERDICT_MIN( metric_vals->aspect_beta, VERDICT_DBL_MAX );
1074  metric_vals->aspect_beta = (double)VERDICT_MAX( metric_vals->aspect_beta, -VERDICT_DBL_MAX );
1075  }
1076 
1077  if( metrics_request_flag & V_TET_ASPECT_GAMMA )
1078  {
1079  if( metric_vals->aspect_gamma > 0 )
1080  metric_vals->aspect_gamma = (double)VERDICT_MIN( metric_vals->aspect_gamma, VERDICT_DBL_MAX );
1081  metric_vals->aspect_gamma = (double)VERDICT_MAX( metric_vals->aspect_gamma, -VERDICT_DBL_MAX );
1082  }
1083 
1084  if( metrics_request_flag & V_TET_VOLUME )
1085  {
1086  if( metric_vals->volume > 0 ) metric_vals->volume = (double)VERDICT_MIN( metric_vals->volume, VERDICT_DBL_MAX );
1087  metric_vals->volume = (double)VERDICT_MAX( metric_vals->volume, -VERDICT_DBL_MAX );
1088  }
1089 
1090  if( metrics_request_flag & V_TET_CONDITION )
1091  {
1092  if( metric_vals->condition > 0 )
1093  metric_vals->condition = (double)VERDICT_MIN( metric_vals->condition, VERDICT_DBL_MAX );
1094  metric_vals->condition = (double)VERDICT_MAX( metric_vals->condition, -VERDICT_DBL_MAX );
1095  }
1096 
1097  if( metrics_request_flag & V_TET_JACOBIAN )
1098  {
1099  if( metric_vals->jacobian > 0 )
1100  metric_vals->jacobian = (double)VERDICT_MIN( metric_vals->jacobian, VERDICT_DBL_MAX );
1101  metric_vals->jacobian = (double)VERDICT_MAX( metric_vals->jacobian, -VERDICT_DBL_MAX );
1102  }
1103 
1104  if( metrics_request_flag & V_TET_SCALED_JACOBIAN )
1105  {
1106  if( metric_vals->scaled_jacobian > 0 )
1107  metric_vals->scaled_jacobian = (double)VERDICT_MIN( metric_vals->scaled_jacobian, VERDICT_DBL_MAX );
1108  metric_vals->scaled_jacobian = (double)VERDICT_MAX( metric_vals->scaled_jacobian, -VERDICT_DBL_MAX );
1109  }
1110 
1111  if( metrics_request_flag & V_TET_SHAPE )
1112  {
1113  if( metric_vals->shape > 0 ) metric_vals->shape = (double)VERDICT_MIN( metric_vals->shape, VERDICT_DBL_MAX );
1114  metric_vals->shape = (double)VERDICT_MAX( metric_vals->shape, -VERDICT_DBL_MAX );
1115  }
1116 
1117  if( metrics_request_flag & V_TET_RELATIVE_SIZE_SQUARED )
1118  {
1119  if( metric_vals->relative_size_squared > 0 )
1120  metric_vals->relative_size_squared =
1121  (double)VERDICT_MIN( metric_vals->relative_size_squared, VERDICT_DBL_MAX );
1122  metric_vals->relative_size_squared =
1123  (double)VERDICT_MAX( metric_vals->relative_size_squared, -VERDICT_DBL_MAX );
1124  }
1125 
1126  if( metrics_request_flag & V_TET_SHAPE_AND_SIZE )
1127  {
1128  if( metric_vals->shape_and_size > 0 )
1129  metric_vals->shape_and_size = (double)VERDICT_MIN( metric_vals->shape_and_size, VERDICT_DBL_MAX );
1130  metric_vals->shape_and_size = (double)VERDICT_MAX( metric_vals->shape_and_size, -VERDICT_DBL_MAX );
1131  }
1132 
1133  if( metrics_request_flag & V_TET_DISTORTION )
1134  {
1135  if( metric_vals->distortion > 0 )
1136  metric_vals->distortion = (double)VERDICT_MIN( metric_vals->distortion, VERDICT_DBL_MAX );
1137  metric_vals->distortion = (double)VERDICT_MAX( metric_vals->distortion, -VERDICT_DBL_MAX );
1138  }
1139 
1140  if( metrics_request_flag & V_TET_ASPECT_RATIO ) metric_vals->aspect_ratio = v_tet_aspect_ratio( 4, coordinates );
1141 
1142  if( metrics_request_flag & V_TET_ASPECT_FROBENIUS )
1143  metric_vals->aspect_frobenius = v_tet_aspect_frobenius( 4, coordinates );
1144 
1145  if( metrics_request_flag & V_TET_MINIMUM_ANGLE ) metric_vals->minimum_angle = v_tet_minimum_angle( 4, coordinates );
1146 
1147  if( metrics_request_flag & V_TET_COLLAPSE_RATIO )
1148  metric_vals->collapse_ratio = v_tet_collapse_ratio( 4, coordinates );
1149 
1150  if( metrics_request_flag & V_TET_RADIUS_RATIO ) metric_vals->radius_ratio = v_tet_radius_ratio( 4, coordinates );
1151 }

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

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

◆ v_tet_radius_ratio()

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

Calculates tet radius ratio metric.

The radius ratio of a tet

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

Definition at line 209 of file V_TetMetric.cpp.

210 {
211 
212  // Determine side vectors
213  VerdictVector side[6];
214 
215  side[0].set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
216  coordinates[1][2] - coordinates[0][2] );
217 
218  side[1].set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
219  coordinates[2][2] - coordinates[1][2] );
220 
221  side[2].set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
222  coordinates[0][2] - coordinates[2][2] );
223 
224  side[3].set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
225  coordinates[3][2] - coordinates[0][2] );
226 
227  side[4].set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
228  coordinates[3][2] - coordinates[1][2] );
229 
230  side[5].set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
231  coordinates[3][2] - coordinates[2][2] );
232 
233  VerdictVector numerator = side[3].length_squared() * ( side[2] * side[0] ) +
234  side[2].length_squared() * ( side[3] * side[0] ) +
235  side[0].length_squared() * ( side[3] * side[2] );
236 
237  double area_sum;
238  area_sum = ( ( side[2] * side[0] ).length() + ( side[3] * side[0] ).length() + ( side[4] * side[1] ).length() +
239  ( side[3] * side[2] ).length() ) *
240  0.5;
241 
242  double volume = v_tet_volume( 4, coordinates );
243 
244  if( fabs( volume ) < VERDICT_DBL_MIN )
245  return (double)VERDICT_DBL_MAX;
246  else
247  {
248  double radius_ratio;
249  radius_ratio = numerator.length() * area_sum / ( 108 * volume * volume );
250 
251  return (double)VERDICT_MIN( radius_ratio, VERDICT_DBL_MAX );
252  }
253 }

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

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

◆ v_tet_relative_size_squared()

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

Calculates tet relative size metric.

the relative size of a tet

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

Definition at line 727 of file V_TetMetric.cpp.

728 {
729  double size;
730  VerdictVector w1, w2, w3;
731  get_weight( w1, w2, w3 );
732  double avg_volume = ( w1 % ( w2 * w3 ) ) / 6.0;
733 
734  double volume = v_tet_volume( 4, coordinates );
735 
736  if( avg_volume < VERDICT_DBL_MIN )
737  return 0.0;
738  else
739  {
740  size = volume / avg_volume;
741  if( size <= VERDICT_DBL_MIN ) return 0.0;
742  if( size > 1 ) size = (double)( 1 ) / size;
743  }
744  return (double)( size * size );
745 }

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

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

◆ v_tet_scaled_jacobian()

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

Calculates tet scaled jacobian.

the scaled jacobian of a tet

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

Definition at line 153 of file V_TetMetric.cpp.

154 {
155 
156  VerdictVector side0, side1, side2, side3, side4, side5;
157 
158  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
159  coordinates[1][2] - coordinates[0][2] );
160 
161  side1.set( coordinates[2][0] - coordinates[1][0], coordinates[2][1] - coordinates[1][1],
162  coordinates[2][2] - coordinates[1][2] );
163 
164  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
165  coordinates[0][2] - coordinates[2][2] );
166 
167  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
168  coordinates[3][2] - coordinates[0][2] );
169 
170  side4.set( coordinates[3][0] - coordinates[1][0], coordinates[3][1] - coordinates[1][1],
171  coordinates[3][2] - coordinates[1][2] );
172 
173  side5.set( coordinates[3][0] - coordinates[2][0], coordinates[3][1] - coordinates[2][1],
174  coordinates[3][2] - coordinates[2][2] );
175 
176  double jacobi;
177 
178  jacobi = side3 % ( side2 * side0 );
179 
180  // products of lengths squared of each edge attached to a node.
181  double length_squared[4] = { side0.length_squared() * side2.length_squared() * side3.length_squared(),
182  side0.length_squared() * side1.length_squared() * side4.length_squared(),
183  side1.length_squared() * side2.length_squared() * side5.length_squared(),
184  side3.length_squared() * side4.length_squared() * side5.length_squared() };
185  int which_node = 0;
186  if( length_squared[1] > length_squared[which_node] ) which_node = 1;
187  if( length_squared[2] > length_squared[which_node] ) which_node = 2;
188  if( length_squared[3] > length_squared[which_node] ) which_node = 3;
189 
190  double length_product = sqrt( length_squared[which_node] );
191  if( length_product < fabs( jacobi ) ) length_product = fabs( jacobi );
192 
193  if( length_product < VERDICT_DBL_MIN ) return (double)VERDICT_DBL_MAX;
194 
195  static const double root_of_2 = sqrt( 2.0 );
196 
197  return (double)( root_of_2 * jacobi / length_product );
198 }

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

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

◆ v_tet_shape()

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

Calculates tet shape metric.

the shape of a tet

3/ condition number of weighted jacobian matrix

Definition at line 691 of file V_TetMetric.cpp.

692 {
693 
694  static const double two_thirds = 2.0 / 3.0;
695  static const double root_of_2 = sqrt( 2.0 );
696 
697  VerdictVector edge0, edge2, edge3;
698 
699  edge0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
700  coordinates[1][2] - coordinates[0][2] );
701 
702  edge2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
703  coordinates[0][2] - coordinates[2][2] );
704 
705  edge3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
706  coordinates[3][2] - coordinates[0][2] );
707 
708  double jacobian = edge3 % ( edge2 * edge0 );
709  if( jacobian < VERDICT_DBL_MIN )
710  {
711  return (double)0.0;
712  }
713  double num = 3 * pow( root_of_2 * jacobian, two_thirds );
714  double den =
715  1.5 * ( edge0 % edge0 + edge2 % edge2 + edge3 % edge3 ) - ( edge0 % -edge2 + -edge2 % edge3 + edge3 % edge0 );
716 
717  if( den < VERDICT_DBL_MIN ) return (double)0.0;
718 
719  return (double)VERDICT_MAX( num / den, 0 );
720 }

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

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

◆ v_tet_shape_and_size()

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

Calculates tet shape-size metric.

the shape and size of a tet

Product of the shape and relative size

Definition at line 752 of file V_TetMetric.cpp.

753 {
754 
755  double shape, size;
756  shape = v_tet_shape( num_nodes, coordinates );
757  size = v_tet_relative_size_squared( num_nodes, coordinates );
758 
759  return (double)( shape * size );
760 }

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

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

◆ v_tet_volume()

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

Calculates tet volume.

the volume of a tet

1/6 * jacobian at a corner node

Definition at line 606 of file V_TetMetric.cpp.

607 {
608 
609  // Determine side vectors
610  VerdictVector side0, side2, side3;
611 
612  side0.set( coordinates[1][0] - coordinates[0][0], coordinates[1][1] - coordinates[0][1],
613  coordinates[1][2] - coordinates[0][2] );
614 
615  side2.set( coordinates[0][0] - coordinates[2][0], coordinates[0][1] - coordinates[2][1],
616  coordinates[0][2] - coordinates[2][2] );
617 
618  side3.set( coordinates[3][0] - coordinates[0][0], coordinates[3][1] - coordinates[0][1],
619  coordinates[3][2] - coordinates[0][2] );
620 
621  return (double)( ( side3 % ( side2 * side0 ) ) / 6.0 );
622 }

References VerdictVector::set().

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

Variable Documentation

◆ verdict_tet_size

double verdict_tet_size = 0

the average volume of a tet

Definition at line 32 of file V_TetMetric.cpp.

Referenced by get_weight(), and v_set_tet_size().