moab
|
00001 #include "TestUtil.hpp" 00002 #include "ElemUtil.hpp" 00003 #include <iostream> 00004 00005 using namespace moab; 00006 00007 void test_hex_nat_coords(); 00008 00009 int main() 00010 { 00011 int rval = 0; 00012 rval += RUN_TEST(test_hex_nat_coords); 00013 return rval; 00014 } 00015 00016 const CartVect cube_corners[8] = { CartVect( 0, 0, 0 ), 00017 CartVect( 1, 0, 0 ), 00018 CartVect( 1, 1, 0 ), 00019 CartVect( 0, 1, 0 ), 00020 CartVect( 0, 0, 1 ), 00021 CartVect( 1, 0, 1 ), 00022 CartVect( 1, 1, 1 ), 00023 CartVect( 0, 1, 1 ) }; 00024 00025 00026 const CartVect hex_corners[8] = { CartVect( 1.0e0, 0.0e0, 0.0e0 ), 00027 CartVect( 1.0e0, 1.0e0, 0.3e0 ), 00028 CartVect( 0.0e0, 2.0e0, 0.6e0 ), 00029 CartVect( 0.2e0, 1.1e0, 0.4e0 ), 00030 CartVect( 1.5e0, 0.3e0, 1.0e0 ), 00031 CartVect( 1.5e0, 1.3e0, 1.0e0 ), 00032 CartVect( 0.5e0, 2.3e0, 1.0e0 ), 00033 CartVect( 0.7e0, 1.4e0, 1.0e0 ) }; 00034 00036 CartVect hex_map( const CartVect& xi, const CartVect* corners ) 00037 { 00038 CartVect x(0.0); 00039 x += (1 - xi[0]) * (1 - xi[1]) * (1 - xi[2]) * corners[0]; 00040 x += (1 + xi[0]) * (1 - xi[1]) * (1 - xi[2]) * corners[1]; 00041 x += (1 + xi[0]) * (1 + xi[1]) * (1 - xi[2]) * corners[2]; 00042 x += (1 - xi[0]) * (1 + xi[1]) * (1 - xi[2]) * corners[3]; 00043 x += (1 - xi[0]) * (1 - xi[1]) * (1 + xi[2]) * corners[4]; 00044 x += (1 + xi[0]) * (1 - xi[1]) * (1 + xi[2]) * corners[5]; 00045 x += (1 + xi[0]) * (1 + xi[1]) * (1 + xi[2]) * corners[6]; 00046 x += (1 - xi[0]) * (1 + xi[1]) * (1 + xi[2]) * corners[7]; 00047 return x *= 0.125; 00048 } 00049 00050 static void hex_bounding_box( const CartVect* corners, CartVect& min, CartVect& max ) 00051 { 00052 min = max = corners[0]; 00053 for (unsigned i = 1; i < 8; ++i) 00054 for (unsigned d = 0; d < 3; ++d) 00055 if (corners[i][d] < min[d]) 00056 min[d] = corners[i][d]; 00057 else if (corners[i][d] > max[d]) 00058 max[d] = corners[i][d]; 00059 } 00060 00061 static bool in_range( const CartVect& xi ) 00062 { return fabs(xi[0]) <= 1 00063 && fabs(xi[1]) <= 1 00064 && fabs(xi[2]) <= 1; 00065 } 00066 00067 void test_hex_nat_coords() 00068 { 00069 CartVect xi, result_xi; 00070 bool valid; 00071 // rename EPS to EPS1 because of conflict with definition of EPS in types.h 00072 // types.h is now included in the header. 00073 const double EPS1 = 1e-6; 00074 00075 // first test with cube because it's easier to debug failures 00076 std::vector<CartVect> cube_corners2; 00077 std::copy(cube_corners, cube_corners+8, std::back_inserter(cube_corners2)); 00078 Element::LinearHex hex(cube_corners2); 00079 for (xi[0] = -1; xi[0] <= 1; xi[0] += 0.2) { 00080 for (xi[1] = -1; xi[1] <= 1; xi[1] += 0.2) { 00081 for (xi[2] = -1; xi[2] <= 1; xi[2] += 0.2) { 00082 const CartVect pt = hex_map(xi, cube_corners); 00083 result_xi = hex.ievaluate(pt, EPS1/10); 00084 double dum = EPS1/10; 00085 valid = hex.inside_nat_space(result_xi, dum); 00086 CHECK(valid); 00087 CHECK_REAL_EQUAL( xi[0], result_xi[0], EPS1 ); 00088 CHECK_REAL_EQUAL( xi[1], result_xi[1], EPS1 ); 00089 CHECK_REAL_EQUAL( xi[2], result_xi[2], EPS1 ); 00090 } 00091 } 00092 } 00093 00094 // now test with distorted hex 00095 std::vector<CartVect> hex_corners2; 00096 std::copy(hex_corners, hex_corners+8, std::back_inserter(hex_corners2)); 00097 Element::LinearHex hex2(hex_corners2); 00098 for (xi[0] = -1; xi[0] <= 1; xi[0] += 0.2) { 00099 for (xi[1] = -1; xi[1] <= 1; xi[1] += 0.2) { 00100 for (xi[2] = -1; xi[2] <= 1; xi[2] += 0.2) { 00101 const CartVect pt = hex_map(xi, hex_corners); 00102 result_xi = hex2.ievaluate(pt, EPS1/10); 00103 double dum = EPS1/10; 00104 valid = hex2.inside_nat_space(result_xi, dum); 00105 CHECK(valid); 00106 CHECK_REAL_EQUAL( xi[0], result_xi[0], EPS1 ); 00107 CHECK_REAL_EQUAL( xi[1], result_xi[1], EPS1 ); 00108 CHECK_REAL_EQUAL( xi[2], result_xi[2], EPS1 ); 00109 } 00110 } 00111 } 00112 00113 // test points outside of element 00114 CartVect x, min, max; 00115 hex_bounding_box( cube_corners, min, max ); 00116 for (x[0] = -1; x[0] <= 2; x[0] += 0.4) { 00117 for (x[1] = -1; x[1] <= 2; x[1] += 0.4) { 00118 for (x[2] = -1; x[2] <= 2; x[2] += 0.4) { 00119 bool in_box = x[0] >= min[0] && x[0] <= max[0] 00120 && x[1] >= min[1] && x[1] <= max[1] 00121 && x[2] >= min[2] && x[2] <= max[2]; 00122 if (in_box) 00123 continue; 00124 result_xi = hex.ievaluate(x, EPS1/10); 00125 double dum = EPS1/10; 00126 valid = hex.inside_nat_space(result_xi, dum); 00127 00128 //std::cout << (valid ? 'y' : 'n'); 00129 CHECK(!valid || !in_range(result_xi)); 00130 } 00131 } 00132 } 00133 //std::cout << std::endl; 00134 00135 hex_bounding_box( hex_corners, min, max ); 00136 for (x[0] = -1; x[0] <= 3; x[0] += 0.5) { 00137 for (x[1] = -2; x[1] <= 4; x[1] += 0.5) { 00138 for (x[2] = -1; x[2] <= 2; x[2] += 0.4) { 00139 bool in_box = x[0] >= min[0] && x[0] <= max[0] 00140 && x[1] >= min[1] && x[1] <= max[1] 00141 && x[2] >= min[2] && x[2] <= max[2]; 00142 if (in_box) 00143 continue; 00144 try { 00145 result_xi = hex2.ievaluate(x, EPS1/10); 00146 } 00147 catch (Element::Map::EvaluationError err) { 00148 valid = false; 00149 } 00150 //std::cout << (valid ? 'y' : 'n'); 00151 CHECK(!valid || !in_range(result_xi)); 00152 } 00153 } 00154 } 00155 //std::cout << std::endl; 00156 } 00157