MOAB: Mesh Oriented datABase  (version 5.4.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 */,
00069                                              const Entity_handle& /* h */,
00070                                              const Points& v,
00071                                              const Point& p,
00072                                              const double tol = 1.e-6 ) const
00073         {
00074             Point result( 3, 0.0 );
00075             bool point_found = solve_inverse( p, result, v, tol ) && is_contained( result, tol );
00076             return std::make_pair( point_found, result );
00077         }
00078 
00079       private:
00080         // This is a hack to avoid a .cpp file and C++11
00081         // reference_points(i,j) will be a 1 or -1;
00082         // This should unroll..
00083         inline double reference_points( const std::size_t& i, const std::size_t& j ) const
00084         {
00085             const double rpts[27][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 },  // reference_points nodes: 0-7
00086                                          { -1, 1, -1 },                                // mid-edge nodes: 8-19
00087                                          { -1, -1, 1 },                 // center-face nodes 20-25  center node  26
00088                                          { 1, -1, 1 },                  //
00089                                          { 1, 1, 1 },    { -1, 1, 1 },  //                    4   ----- 19   -----  7
00090                                          { 0, -1, -1 },                 //                .   |                 .   |
00091                                          { 1, 0, -1 },                  //            16         25         18      |
00092                                          { 0, 1, -1 },                  //         .          |          .          |
00093                                          { -1, 0, -1 },                 //      5   ----- 17   -----  6             |
00094                                          { -1, -1, 0 },                 //      |            12       | 23         15
00095                                          { 1, -1, 0 },                  //      |                     |             |
00096                                          { 1, 1, 0 },                   //      |     20      |  26   |     22      |
00097                                          { -1, 1, 0 },                  //      |                     |             |
00098                                          { 0, -1, 1 },                  //     13         21  |      14             |
00099                                          { 1, 0, 1 },                   //      |             0   ----- 11   -----  3
00100                                          { 0, 1, 1 },                   //      |         .           |         .
00101                                          { -1, 0, 1 },                  //      |      8         24   |     10
00102                                          { 0, -1, 0 },                  //      |  .                  |  .
00103                                          { 1, 0, 0 },                   //      1   -----  9   -----  2
00104                                          { 0, 1, 0 },                   //
00105                                          { -1, 0, 0 },   { 0, 0, -1 },  { 0, 0, 1 },  { 0, 0, 0 } };
00106             return rpts[i][j];
00107         }
00108 
00109         template < typename Point >
00110         bool is_contained( const Point& p, const double tol ) const
00111         {
00112             // just look at the box+tol here
00113             return ( p[0] >= -1. - tol ) && ( p[0] <= 1. + tol ) && ( p[1] >= -1. - tol ) && ( p[1] <= 1. + tol ) &&
00114                    ( p[2] >= -1. - tol ) && ( p[2] <= 1. + tol );
00115         }
00116 
00117         template < typename Point, typename Points >
00118         bool solve_inverse( const Point& x, Point& xi, const Points& points, const double tol = 1.e-6 ) const
00119         {
00120             const double error_tol_sqr = tol * tol;
00121             Point delta( 3, 0.0 );
00122             xi = delta;
00123             evaluate( xi, points, delta );
00124             vec_subtract( delta, x );
00125             std::size_t num_iterations = 0;
00126 #ifdef QUADRATIC_HEX_DEBUG
00127             std::stringstream ss;
00128             ss << "Point: ";
00129             ss << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
00130             ss << "Hex: ";
00131             for( int i = 0; i < 8; ++i )
00132             {
00133                 ss << points[i][0] << ", " << points[i][1] << ", " << points[i][2] << std::endl;
00134             }
00135             ss << std::endl;
00136 #endif
00137             while( normsq( delta ) > error_tol_sqr )
00138             {
00139 #ifdef QUADRATIC_HEX_DEBUG
00140                 ss << "Iter #: " << num_iterations << " Err: " << sqrt( normsq( delta ) ) << " Iterate: ";
00141                 ss << xi[0] << ", " << xi[1] << ", " << xi[2] << std::endl;
00142 #endif
00143                 if( ++num_iterations >= 5 )
00144                 {
00145                     return false;
00146                 }
00147                 Matrix J;
00148                 jacobian( xi, points, J );
00149                 double det = moab::Matrix::determinant3( J );
00150                 if( fabs( det ) < 1.e-10 )
00151                 {
00152 #ifdef QUADRATIC_HEX_DEBUG
00153                     std::cerr << ss.str();
00154 #endif
00155 #ifndef QUADRATIC_HEX_DEBUG
00156                     std::cerr << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
00157 #endif
00158                     std::cerr << "inverse solve failure: det: " << det << std::endl;
00159                     exit( -1 );
00160                 }
00161                 vec_subtract( xi, moab::Matrix::inverse( J, 1.0 / det ) * delta );
00162                 evaluate( xi, points, delta );
00163                 vec_subtract( delta, x );
00164             }
00165             return true;
00166         }
00167 
00168         template < typename Point, typename Points >
00169         Point& evaluate( const Point& p, const Points& points, Point& f ) const
00170         {
00171             typedef typename Points::value_type Vector;
00172             Vector result;
00173             for( int i = 0; i < 3; ++i )
00174             {
00175                 result[i] = 0;
00176             }
00177             for( unsigned i = 0; i < 27; ++i )
00178             {
00179                 const double sh = SH( reference_points( i, 0 ), p[0] ) * SH( reference_points( i, 1 ), p[1] ) *
00180                                   SH( reference_points( i, 2 ), p[2] );
00181                 result += sh * points[i];
00182             }
00183             for( int i = 0; i < 3; ++i )
00184             {
00185                 f[i] = result[i];
00186             }
00187             return f;
00188         }
00189         template < typename Point, typename Field >
00190         double evaluate_scalar_field( const Point& p, const Field& field ) const
00191         {
00192             double x = 0.0;
00193             for( int i = 0; i < 27; i++ )
00194             {
00195                 const double sh = SH( reference_points( i, 0 ), p[0] ) * SH( reference_points( i, 1 ), p[1] ) *
00196                                   SH( reference_points( i, 2 ), p[2] );
00197                 x += sh * field[i];
00198             }
00199             return x;
00200         }
00201         template < typename Field, typename Points >
00202         double integrate_scalar_field( const Points& p, const Field& field_values ) const
00203         {
00204             // TODO: gaussian integration , probably 2x2x2
00205             return 0.;
00206         }
00207 
00208         template < typename Point, typename Points >
00209         Matrix& jacobian( const Point& p, const Points& /* points */, Matrix& J ) const
00210         {
00211             J = Matrix( 0.0 );
00212             for( int i = 0; i < 27; i++ )
00213             {
00214                 const double sh[3]  = { SH( reference_points( i, 0 ), p[0] ), SH( reference_points( i, 1 ), p[1] ),
00215                                        SH( reference_points( i, 2 ), p[2] ) };
00216                 const double dsh[3] = { DSH( reference_points( i, 0 ), p[0] ), DSH( reference_points( i, 1 ), p[1] ),
00217                                         DSH( reference_points( i, 2 ), p[2] ) };
00218                 for( int j = 0; j < 3; j++ )
00219                 {
00220                     // dxj/dr first column
00221                     J( j, 0 ) += dsh[0] * sh[1] * sh[2] * reference_points( i, j );
00222                     J( j, 1 ) += sh[0] * dsh[1] * sh[2] * reference_points( i, j );  // dxj/ds
00223                     J( j, 2 ) += sh[0] * sh[1] * dsh[2] * reference_points( i, j );  // dxj/dt
00224                 }
00225             }
00226             return J;
00227         }
00228 
00229       private:
00230     };  // Class Quadratic_hex_map
00231 
00232 }  // namespace element_utility
00233 
00234 }  // namespace moab
00235 #endif  // MOAB_QUADRATIC_HEX_nPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines