 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
#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... | |
| #define VERDICT_EXPORTS |
Definition at line 23 of file V_TetMetric.cpp.
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
| 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().
1.9.1.