MOAB: Mesh Oriented datABase  (version 5.2.1)
quadratic_hex_map.hpp
Go to the documentation of this file.
00001 #ifndef MOAB_QUADRATIC_HEX_HPP
00002 #define MOAB_QUADRATIC_HEX_HPP
00003 
00004 #include "moab/Matrix3.hpp"
00005 #include "moab/CartVect.hpp"
00006 #include <sstream>
00007 #include <iomanip>
00008 #include <iostream>
00009 
00010 namespace moab
00011 {
00012 
00013 namespace element_utility
00014 {
00015 
00016     namespace
00017     {
00018 
00019         double SH( const int i, const double xi )
00020         {
00021             switch( i )
00022             {
00023                 case -1:
00024                     return ( xi * xi - xi ) / 2;
00025                 case 0:
00026                     return 1 - xi * xi;
00027                 case 1:
00028                     return ( xi * xi + xi ) / 2;
00029                 default:
00030                     return 0.;
00031             }
00032         }
00033         double DSH( const int i, const double xi )
00034         {
00035             switch( i )
00036             {
00037                 case -1:
00038                     return xi - 0.5;
00039                 case 0:
00040                     return -2 * xi;
00041                 case 1:
00042                     return xi + 0.5;
00043                 default:
00044                     return 0.;
00045             }
00046         }
00047 
00048     }  // namespace
00049 
00050     template < typename _Matrix >
00051     class Quadratic_hex_map
00052     {
00053       public:
00054         typedef _Matrix Matrix;
00055 
00056       private:
00057         typedef Quadratic_hex_map< Matrix > Self;
00058 
00059       public:
00060         // Constructor
00061         Quadratic_hex_map() {}
00062         // Copy constructor
00063         Quadratic_hex_map( const Self& f ) {}
00064 
00065       public:
00066         // Natural coordinates
00067         template < typename Moab, typename Entity_handle, typename Points, typename Point >
00068         std::pair< bool, Point > operator()( const Moab& /* moab */, const Entity_handle& /* h */, const Points& v,
00069                                              const Point& p, const double tol = 1.e-6 ) const
00070         {
00071             Point result( 3, 0.0 );
00072             bool point_found = solve_inverse( p, result, v, tol ) && is_contained( result, tol );
00073             return std::make_pair( point_found, result );
00074         }
00075 
00076       private:
00077         // This is a hack to avoid a .cpp file and C++11
00078         // reference_points(i,j) will be a 1 or -1;
00079         // This should unroll..
00080         inline double reference_points( const std::size_t& i, const std::size_t& j ) const
00081         {
00082             const double rpts[27][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 },  // reference_points nodes: 0-7
00083                                          { -1, 1, -1 },                                // mid-edge nodes: 8-19
00084                                          { -1, -1, 1 },                 // center-face nodes 20-25  center node  26
00085                                          { 1, -1, 1 },                  //
00086                                          { 1, 1, 1 },    { -1, 1, 1 },  //                    4   ----- 19   -----  7
00087                                          { 0, -1, -1 },                 //                .   |                 .   |
00088                                          { 1, 0, -1 },                  //            16         25         18      |
00089                                          { 0, 1, -1 },                  //         .          |          .          |
00090                                          { -1, 0, -1 },                 //      5   ----- 17   -----  6             |
00091                                          { -1, -1, 0 },                 //      |            12       | 23         15
00092                                          { 1, -1, 0 },                  //      |                     |             |
00093                                          { 1, 1, 0 },                   //      |     20      |  26   |     22      |
00094                                          { -1, 1, 0 },                  //      |                     |             |
00095                                          { 0, -1, 1 },                  //     13         21  |      14             |
00096                                          { 1, 0, 1 },                   //      |             0   ----- 11   -----  3
00097                                          { 0, 1, 1 },                   //      |         .           |         .
00098                                          { -1, 0, 1 },                  //      |      8         24   |     10
00099                                          { 0, -1, 0 },                  //      |  .                  |  .
00100                                          { 1, 0, 0 },                   //      1   -----  9   -----  2
00101                                          { 0, 1, 0 },                   //
00102                                          { -1, 0, 0 },   { 0, 0, -1 },  { 0, 0, 1 },  { 0, 0, 0 } };
00103             return rpts[i][j];
00104         }
00105 
00106         template < typename Point >
00107         bool is_contained( const Point& p, const double tol ) const
00108         {
00109             // just look at the box+tol here
00110             return ( p[0] >= -1. - tol ) && ( p[0] <= 1. + tol ) && ( p[1] >= -1. - tol ) && ( p[1] <= 1. + tol ) &&
00111                    ( p[2] >= -1. - tol ) && ( p[2] <= 1. + tol );
00112         }
00113 
00114         template < typename Point, typename Points >
00115         bool solve_inverse( const Point& x, Point& xi, const Points& points, const double tol = 1.e-6 ) const
00116         {
00117             const double error_tol_sqr = tol * tol;
00118             Point delta( 3, 0.0 );
00119             xi = delta;
00120             evaluate( xi, points, delta );
00121             vec_subtract( delta, x );
00122             std::size_t num_iterations = 0;
00123 #ifdef QUADRATIC_HEX_DEBUG
00124             std::stringstream ss;
00125             ss << "Point: ";
00126             ss << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
00127             ss << "Hex: ";
00128             for( int i = 0; i < 8; ++i )
00129             {
00130                 ss << points[i][0] << ", " << points[i][1] << ", " << points[i][2] << std::endl;
00131             }
00132             ss << std::endl;
00133 #endif
00134             while( normsq( delta ) > error_tol_sqr )
00135             {
00136 #ifdef QUADRATIC_HEX_DEBUG
00137                 ss << "Iter #: " << num_iterations << " Err: " << sqrt( normsq( delta ) ) << " Iterate: ";
00138                 ss << xi[0] << ", " << xi[1] << ", " << xi[2] << std::endl;
00139 #endif
00140                 if( ++num_iterations >= 5 ) { return false; }
00141                 Matrix J;
00142                 jacobian( xi, points, J );
00143                 double det = moab::Matrix::determinant3( J );
00144                 if( fabs( det ) < 1.e-10 )
00145                 {
00146 #ifdef QUADRATIC_HEX_DEBUG
00147                     std::cerr << ss.str();
00148 #endif
00149 #ifndef QUADRATIC_HEX_DEBUG
00150                     std::cerr << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
00151 #endif
00152                     std::cerr << "inverse solve failure: det: " << det << std::endl;
00153                     exit( -1 );
00154                 }
00155                 vec_subtract( xi, moab::Matrix::inverse( J, 1.0 / det ) * delta );
00156                 evaluate( xi, points, delta );
00157                 vec_subtract( delta, x );
00158             }
00159             return true;
00160         }
00161 
00162         template < typename Point, typename Points >
00163         Point& evaluate( const Point& p, const Points& points, Point& f ) const
00164         {
00165             typedef typename Points::value_type Vector;
00166             Vector result;
00167             for( int i = 0; i < 3; ++i )
00168             {
00169                 result[i] = 0;
00170             }
00171             for( unsigned i = 0; i < 27; ++i )
00172             {
00173                 const double sh = SH( reference_points( i, 0 ), p[0] ) * SH( reference_points( i, 1 ), p[1] ) *
00174                                   SH( reference_points( i, 2 ), p[2] );
00175                 result += sh * points[i];
00176             }
00177             for( int i = 0; i < 3; ++i )
00178             {
00179                 f[i] = result[i];
00180             }
00181             return f;
00182         }
00183         template < typename Point, typename Field >
00184         double evaluate_scalar_field( const Point& p, const Field& field ) const
00185         {
00186             double x = 0.0;
00187             for( int i = 0; i < 27; i++ )
00188             {
00189                 const double sh = SH( reference_points( i, 0 ), p[0] ) * SH( reference_points( i, 1 ), p[1] ) *
00190                                   SH( reference_points( i, 2 ), p[2] );
00191                 x += sh * field[i];
00192             }
00193             return x;
00194         }
00195         template < typename Field, typename Points >
00196         double integrate_scalar_field( const Points& p, const Field& field_values ) const
00197         {
00198             // TODO: gaussian integration , probably 2x2x2
00199             return 0.;
00200         }
00201 
00202         template < typename Point, typename Points >
00203         Matrix& jacobian( const Point& p, const Points& /* points */, Matrix& J ) const
00204         {
00205             J = Matrix( 0.0 );
00206             for( int i = 0; i < 27; i++ )
00207             {
00208                 const double sh[3]  = { SH( reference_points( i, 0 ), p[0] ), SH( reference_points( i, 1 ), p[1] ),
00209                                        SH( reference_points( i, 2 ), p[2] ) };
00210                 const double dsh[3] = { DSH( reference_points( i, 0 ), p[0] ), DSH( reference_points( i, 1 ), p[1] ),
00211                                         DSH( reference_points( i, 2 ), p[2] ) };
00212                 for( int j = 0; j < 3; j++ )
00213                 {
00214                     // dxj/dr first column
00215                     J( j, 0 ) += dsh[0] * sh[1] * sh[2] * reference_points( i, j );
00216                     J( j, 1 ) += sh[0] * dsh[1] * sh[2] * reference_points( i, j );  // dxj/ds
00217                     J( j, 2 ) += sh[0] * sh[1] * dsh[2] * reference_points( i, j );  // dxj/dt
00218                 }
00219             }
00220             return J;
00221         }
00222 
00223       private:
00224     };  // Class Quadratic_hex_map
00225 
00226 }  // namespace element_utility
00227 
00228 }  // namespace moab
00229 #endif  // MOAB_QUADRATIC_HEX_nPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines