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