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