MOAB: Mesh Oriented datABase  (version 5.4.0)
spectral_hex_map.hpp
Go to the documentation of this file.
00001 #ifndef MOAB_SPECTRAL_HEX_HPP
00002 #define MOAB_SPECTRAL_HEX_HPP
00003 
00004 #include "moab/Matrix3.hpp"
00005 #include "moab/CartVect.hpp"
00006 #include "moab/FindPtFuncs.h"
00007 #include <sstream>
00008 #include <iomanip>
00009 #include <iostream>
00010 
00011 namespace moab
00012 {
00013 
00014 namespace element_utility
00015 {
00016 
00017     namespace
00018     {
00019     }  // namespace
00020 
00021     template < typename _Matrix >
00022     class Spectral_hex_map
00023     {
00024       public:
00025         typedef _Matrix Matrix;
00026 
00027       private:
00028         typedef Spectral_hex_map< Matrix > Self;
00029 
00030       public:
00031         // Constructor
00032         Spectral_hex_map(){};
00033         Spectral_hex_map( int order )
00034         {
00035             initialize_spectral_hex( order );
00036         }
00037         // Copy constructor
00038         Spectral_hex_map( const Self& f ) {}
00039 
00040       private:
00041         void initialize_spectral_hex( int order )
00042         {
00043             if( _init && _n == order )
00044             {
00045                 return;
00046             }
00047             if( _init && _n != order )
00048             {
00049                 free_data();
00050             }
00051             _init = true;
00052             _n    = order;
00053             for( int d = 0; d < 3; d++ )
00054             {
00055                 lobatto_nodes( _z[d], _n );
00056                 lagrange_setup( &_ld[d], _z[d], _n );
00057             }
00058             opt_alloc_3( &_data, _ld );
00059             std::size_t nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
00060             _odwork = tmalloc( real, 6 * nf + 9 * ne + nw );
00061         }
00062 
00063         void free_data()
00064         {
00065             for( int d = 0; d < 3; d++ )
00066             {
00067                 free( _z[d] );
00068                 lagrange_free( &_ld[d] );
00069             }
00070             opt_free_3( &_data );
00071             free( _odwork );
00072         }
00073 
00074       public:
00075         // Natural coordinates
00076         template < typename Moab, typename Entity_handle, typename Points, typename Point >
00077         std::pair< bool, Point > operator()( const Moab& /* moab */,
00078                                              const Entity_handle& /* h */,
00079                                              const Points& v,
00080                                              const Point& p,
00081                                              const double tol = 1.e-6 )
00082         {
00083             Point result( 3, 0.0 );
00084             /*
00085             moab.tag_get_by_ptr(_xm1Tag, &eh, 1,(const void **) &_xyz[ 0] );
00086             moab.tag_get_by_ptr(_ym1Tag, &eh, 1,(const void **) &_xyz[ 1] );
00087             moab.tag_get_by_ptr(_zm1Tag, &eh, 1,(const void **) &_xyz[ 2] );
00088             */
00089             bool point_found = solve_inverse( p, result, v, tol ) && is_contained( result, tol );
00090             return std::make_pair( point_found, result );
00091         }
00092 
00093       private:
00094         void set_gl_points( double* x, double* y, double* z )
00095         {
00096             _xyz[0] = x;
00097             _xyz[1] = y;
00098             _xyz[2] = z;
00099         }
00100         template < typename Point >
00101         bool is_contained( const Point& p, const double tol ) const
00102         {
00103             // just look at the box+tol here
00104             return ( p[0] >= -1. - tol ) && ( p[0] <= 1. + tol ) && ( p[1] >= -1. - tol ) && ( p[1] <= 1. + tol ) &&
00105                    ( p[2] >= -1. - tol ) && ( p[2] <= 1. + tol );
00106         }
00107 
00108         template < typename Point, typename Points >
00109         bool solve_inverse( const Point& x, Point& xi, const Points& points, const double tol = 1.e-6 )
00110         {
00111             const double error_tol_sqr = tol * tol;
00112             Point delta( 3, 0.0 );
00113             xi = delta;
00114             evaluate( xi, points, delta );
00115             vec_subtract( delta, x );
00116             std::size_t num_iterations = 0;
00117 #ifdef SPECTRAL_HEX_DEBUG
00118             std::stringstream ss;
00119             ss << "Point: ";
00120             ss << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
00121             ss << "Hex: ";
00122             for( int i = 0; i < 8; ++i )
00123             {
00124                 ss << points[i][0] << ", " << points[i][1] << ", " << points[i][2] << std::endl;
00125             }
00126             ss << std::endl;
00127 #endif
00128             while( normsq( delta ) > error_tol_sqr )
00129             {
00130 #ifdef SPECTRAL_HEX_DEBUG
00131                 ss << "Iter #: " << num_iterations << " Err: " << sqrt( normsq( delta ) ) << " Iterate: ";
00132                 ss << xi[0] << ", " << xi[1] << ", " << xi[2] << std::endl;
00133 #endif
00134                 if( ++num_iterations >= 5 )
00135                 {
00136                     return false;
00137                 }
00138                 Matrix J;
00139                 jacobian( xi, points, J );
00140                 double det = moab::Matrix::determinant3( J );
00141                 if( fabs( det ) < 1.e-10 )
00142                 {
00143 #ifdef SPECTRAL_HEX_DEBUG
00144                     std::cerr << ss.str();
00145 #endif
00146 #ifndef SPECTRAL_HEX_DEBUG
00147                     std::cerr << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
00148 #endif
00149                     std::cerr << "inverse solve failure: det: " << det << std::endl;
00150                     exit( -1 );
00151                 }
00152                 vec_subtract( xi, moab::Matrix::inverse( J, 1.0 / det ) * delta );
00153                 evaluate( xi, points, delta );
00154                 vec_subtract( delta, x );
00155             }
00156             return true;
00157         }
00158 
00159         template < typename Point, typename Points >
00160         Point& evaluate( const Point& p, const Points& /* points */, Point& f )
00161         {
00162             for( int d = 0; d < 3; ++d )
00163             {
00164                 lagrange_0( &_ld[d], p[0] );
00165             }
00166             for( int d = 0; d < 3; ++d )
00167             {
00168                 f[d] = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, _xyz[d], _odwork );
00169             }
00170             return f;
00171         }
00172 
00173         template < typename Point, typename Field >
00174         double evaluate_scalar_field( const Point& p, const Field& field ) const
00175         {
00176             int d;
00177             for( d = 0; d < 3; d++ )
00178             {
00179                 lagrange_0( &_ld[d], p[d] );
00180             }
00181             return tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, field, _odwork );
00182         }
00183         template < typename Points, typename Field >
00184         double integrate_scalar_field( const Points& p, const Field& field ) const
00185         {
00186             // set the position of GL points
00187             // set the positions of GL nodes, before evaluations
00188             _data.elx[0] = _xyz[0];
00189             _data.elx[1] = _xyz[1];
00190             _data.elx[2] = _xyz[2];
00191             double xi[3];
00192             // triple loop; the most inner loop is in r direction, then s, then t
00193             double integral = 0.;
00194             // double volume = 0;
00195             int index = 0;  // used fr the inner loop
00196             for( int k = 0; k < _n; k++ )
00197             {
00198                 xi[2] = _ld[2].z[k];
00199                 // double wk= _ld[2].w[k];
00200                 for( int j = 0; j < _n; j++ )
00201                 {
00202                     xi[1] = _ld[1].z[j];
00203                     // double wj= _ld[1].w[j];
00204                     for( int i = 0; i < _n; i++ )
00205                     {
00206                         xi[0] = _ld[0].z[i];
00207                         // double wi= _ld[0].w[i];
00208                         opt_vol_set_intp_3( (opt_data_3*)&_data, xi );  // cast away const-ness
00209                         double wk = _ld[2].J[k];
00210                         double wj = _ld[1].J[j];
00211                         double wi = _ld[0].J[i];
00212                         Matrix3 J( 0. );
00213                         for( int n = 0; n < 8; n++ )
00214                             J( n / 3, n % 3 ) = _data.jac[n];
00215                         double bm = wk * wj * wi * J.determinant();
00216                         integral += bm * field[index++];
00217                         // volume +=bm;
00218                     }
00219                 }
00220             }
00221             // std::cout << "volume: " << volume << "\n";
00222             return integral;
00223         }
00224 
00225         template < typename Point, typename Points >
00226         Matrix& jacobian( const Point& /* p */, const Points& /* points */, Matrix& J )
00227         {
00228             real x[3];
00229             for( int i = 0; i < 3; ++i )
00230             {
00231                 _data.elx[i] = _xyz[i];
00232             }
00233             opt_vol_set_intp_3( &_data, x );
00234             for( int i = 0; i < 9; ++i )
00235             {
00236                 J( i % 3, i / 3 ) = _data.jac[i];
00237             }
00238             return J;
00239         }
00240 
00241       private:
00242         bool _init;
00243         int _n;
00244         real* _z[3];
00245         lagrange_data _ld[3];
00246         opt_data_3 _data;
00247         real* _odwork;
00248         real* _xyz[3];
00249     };  // Class Spectral_hex_map
00250 
00251 }  // namespace element_utility
00252 
00253 }  // namespace moab
00254 #endif  // MOAB_SPECTRAL_HEX_nPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines