Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ElemUtilTest.cpp
Go to the documentation of this file.
1 #include "TestUtil.hpp"
2 #include "ElemUtil.hpp"
3 #include <iostream>
4 
5 using namespace moab;
6 
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 
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 }