moab
ElemUtilTest.cpp
Go to the documentation of this file.
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   
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines