 |
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
|
1 #include "TestUtil.hpp"
2 #include "ElemUtil.hpp"
3 #include <iostream>
4
5 using namespace moab;
6
7 void test_hex_nat_coords();
8
9 int main()
10 {
11 int rval = 0;
12 rval += RUN_TEST( test_hex_nat_coords );
13 return rval;
14 }
15
16 const CartVect cube_corners[8] = { CartVect( 0, 0, 0 ), CartVect( 1, 0, 0 ), CartVect( 1, 1, 0 ), CartVect( 0, 1, 0 ),
17 CartVect( 0, 0, 1 ), CartVect( 1, 0, 1 ), CartVect( 1, 1, 1 ), CartVect( 0, 1, 1 ) };
18
19 const CartVect hex_corners[8] = { CartVect( 1.0e0, 0.0e0, 0.0e0 ), CartVect( 1.0e0, 1.0e0, 0.3e0 ),
20 CartVect( 0.0e0, 2.0e0, 0.6e0 ), CartVect( 0.2e0, 1.1e0, 0.4e0 ),
21 CartVect( 1.5e0, 0.3e0, 1.0e0 ), CartVect( 1.5e0, 1.3e0, 1.0e0 ),
22 CartVect( 0.5e0, 2.3e0, 1.0e0 ), CartVect( 0.7e0, 1.4e0, 1.0e0 ) };
23
24 /** shape function for trilinear hex */
25 CartVect hex_map( const CartVect& xi, const CartVect* corners )
26 {
27 CartVect x( 0.0 );
28 x += ( 1 - xi[0] ) * ( 1 - xi[1] ) * ( 1 - xi[2] ) * corners[0];
29 x += ( 1 + xi[0] ) * ( 1 - xi[1] ) * ( 1 - xi[2] ) * corners[1];
30 x += ( 1 + xi[0] ) * ( 1 + xi[1] ) * ( 1 - xi[2] ) * corners[2];
31 x += ( 1 - xi[0] ) * ( 1 + xi[1] ) * ( 1 - xi[2] ) * corners[3];
32 x += ( 1 - xi[0] ) * ( 1 - xi[1] ) * ( 1 + xi[2] ) * corners[4];
33 x += ( 1 + xi[0] ) * ( 1 - xi[1] ) * ( 1 + xi[2] ) * corners[5];
34 x += ( 1 + xi[0] ) * ( 1 + xi[1] ) * ( 1 + xi[2] ) * corners[6];
35 x += ( 1 - xi[0] ) * ( 1 + xi[1] ) * ( 1 + xi[2] ) * corners[7];
36 return x *= 0.125;
37 }
38
39 static void hex_bounding_box( const CartVect* corners, CartVect& min, CartVect& max )
40 {
41 min = max = corners[0];
42 for( unsigned i = 1; i < 8; ++i )
43 for( unsigned d = 0; d < 3; ++d )
44 if( corners[i][d] < min[d] )
45 min[d] = corners[i][d];
46 else if( corners[i][d] > max[d] )
47 max[d] = corners[i][d];
48 }
49
50 static bool in_range( const CartVect& xi )
51 {
52 return fabs( xi[0] ) <= 1 && fabs( xi[1] ) <= 1 && fabs( xi[2] ) <= 1;
53 }
54
55 void test_hex_nat_coords()
56 {
57 CartVect xi, result_xi;
58 bool valid;
59 // rename EPS to EPS1 because of conflict with definition of EPS in types.h
60 // types.h is now included in the header.
61 const double EPS1 = 1e-6;
62
63 // first test with cube because it's easier to debug failures
64 std::vector< CartVect > cube_corners2;
65 std::copy( cube_corners, cube_corners + 8, std::back_inserter( cube_corners2 ) );
66 Element::LinearHex hex( cube_corners2 );
67 for( xi[0] = -1; xi[0] <= 1; xi[0] += 0.2 )
68 {
69 for( xi[1] = -1; xi[1] <= 1; xi[1] += 0.2 )
70 {
71 for( xi[2] = -1; xi[2] <= 1; xi[2] += 0.2 )
72 {
73 const CartVect pt = hex_map( xi, cube_corners );
74 result_xi = hex.ievaluate( pt, EPS1 / 10 );
75 double dum = EPS1 / 10;
76 valid = hex.inside_nat_space( result_xi, dum );
77 CHECK( valid );
78 CHECK_REAL_EQUAL( xi[0], result_xi[0], EPS1 );
79 CHECK_REAL_EQUAL( xi[1], result_xi[1], EPS1 );
80 CHECK_REAL_EQUAL( xi[2], result_xi[2], EPS1 );
81 }
82 }
83 }
84
85 // now test with distorted hex
86 std::vector< CartVect > hex_corners2;
87 std::copy( hex_corners, hex_corners + 8, std::back_inserter( hex_corners2 ) );
88 Element::LinearHex hex2( hex_corners2 );
89 for( xi[0] = -1; xi[0] <= 1; xi[0] += 0.2 )
90 {
91 for( xi[1] = -1; xi[1] <= 1; xi[1] += 0.2 )
92 {
93 for( xi[2] = -1; xi[2] <= 1; xi[2] += 0.2 )
94 {
95 const CartVect pt = hex_map( xi, hex_corners );
96 result_xi = hex2.ievaluate( pt, EPS1 / 10 );
97 double dum = EPS1 / 10;
98 valid = hex2.inside_nat_space( result_xi, dum );
99 CHECK( valid );
100 CHECK_REAL_EQUAL( xi[0], result_xi[0], EPS1 );
101 CHECK_REAL_EQUAL( xi[1], result_xi[1], EPS1 );
102 CHECK_REAL_EQUAL( xi[2], result_xi[2], EPS1 );
103 }
104 }
105 }
106
107 // test points outside of element
108 CartVect x, min, max;
109 hex_bounding_box( cube_corners, min, max );
110 for( x[0] = -1; x[0] <= 2; x[0] += 0.4 )
111 {
112 for( x[1] = -1; x[1] <= 2; x[1] += 0.4 )
113 {
114 for( x[2] = -1; x[2] <= 2; x[2] += 0.4 )
115 {
116 bool in_box = x[0] >= min[0] && x[0] <= max[0] && x[1] >= min[1] && x[1] <= max[1] && x[2] >= min[2] &&
117 x[2] <= max[2];
118 if( in_box ) continue;
119 result_xi = hex.ievaluate( x, EPS1 / 10 );
120 double dum = EPS1 / 10;
121 valid = hex.inside_nat_space( result_xi, dum );
122
123 // std::cout << (valid ? 'y' : 'n');
124 CHECK( !valid || !in_range( result_xi ) );
125 }
126 }
127 }
128 // std::cout << std::endl;
129
130 hex_bounding_box( hex_corners, min, max );
131 for( x[0] = -1; x[0] <= 3; x[0] += 0.5 )
132 {
133 for( x[1] = -2; x[1] <= 4; x[1] += 0.5 )
134 {
135 for( x[2] = -1; x[2] <= 2; x[2] += 0.4 )
136 {
137 bool in_box = x[0] >= min[0] && x[0] <= max[0] && x[1] >= min[1] && x[1] <= max[1] && x[2] >= min[2] &&
138 x[2] <= max[2];
139 if( in_box ) continue;
140 try
141 {
142 result_xi = hex2.ievaluate( x, EPS1 / 10 );
143 }
144 catch( Element::Map::EvaluationError& /* err */ )
145 {
146 valid = false;
147 }
148 // std::cout << (valid ? 'y' : 'n');
149 CHECK( !valid || !in_range( result_xi ) );
150 }
151 }
152 }
153 // std::cout << std::endl;
154 }
1.9.1.