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