MOAB: Mesh Oriented datABase  (version 5.2.1)
SpectralHex.cpp
Go to the documentation of this file.
00001 #include "moab/LocalDiscretization/SpectralHex.hpp"
00002 #include "moab/Forward.hpp"
00003 
00004 namespace moab
00005 {
00006 
00007 // SpectralHex
00008 
00009 // filescope for static member data that is cached
00010 int SpectralHex::_n;
00011 real* SpectralHex::_z[3];
00012 lagrange_data SpectralHex::_ld[3];
00013 opt_data_3 SpectralHex::_data;
00014 real* SpectralHex::_odwork;
00015 
00016 bool SpectralHex::_init = false;
00017 
00018 void SpectralHex::initFcn( const double* verts, const int nverts, double*& work )
00019 {
00020     if( _init && _n == order ) return;
00021     if( _init && _n != order )
00022     {
00023         // TODO: free data cached
00024         freedata();
00025     }
00026     // compute stuff that depends only on order
00027     _init = true;
00028     _n    = order;
00029     // triplicates! n is the same in all directions !!!
00030     for( int d = 0; d < 3; d++ )
00031     {
00032         _z[d] = tmalloc( double, _n );
00033         lobatto_nodes( _z[d], _n );
00034         lagrange_setup( &_ld[d], _z[d], _n );
00035     }
00036     opt_alloc_3( &_data, _ld );
00037 
00038     unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
00039     _odwork = tmalloc( double, 6 * nf + 9 * ne + nw );
00040 }
00041 void SpectralHex::freedata()
00042 {
00043     for( int d = 0; d < 3; d++ )
00044     {
00045         free( _z[d] );
00046         lagrange_free( &_ld[d] );
00047     }
00048     opt_free_3( &_data );
00049     free( _odwork );
00050 }
00051 
00052 void SpectralHex::set_gl_points( double* x, double* y, double* z )
00053 {
00054     _xyz[0] = x;
00055     _xyz[1] = y;
00056     _xyz[2] = z;
00057 }
00058 CartVect SpectralHex::evaluate( const CartVect& params ) const
00059 {
00060     // piece that we shouldn't want to cache
00061     int d = 0;
00062     for( d = 0; d < 3; d++ )
00063     {
00064         lagrange_0( &_ld[d], params[d] );
00065     }
00066     CartVect result;
00067     for( d = 0; d < 3; d++ )
00068     {
00069         result[d] = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n,
00070                                _xyz[d],  // this is the "field"
00071                                _odwork );
00072     }
00073     return result;
00074 }
00075 // replicate the functionality of hex_findpt
00076 bool SpectralHex::evaluate_reverse( CartVect const& xyz, CartVect& params, double iter_tol, const double inside_tol,
00077                                     const CartVect& init ) const
00078 {
00079     params = init;
00080 
00081     // find nearest point
00082     double x_star[3];
00083     xyz.get( x_star );
00084 
00085     double r[3] = { 0, 0, 0 };  // initial guess for parametric coords
00086     unsigned c  = opt_no_constraints_3;
00087     double dist = opt_findpt_3( &_data, (const double**)_xyz, x_star, r, &c );
00088     // if it did not converge, get out with throw...
00089     if( dist > 0.9e+30 ) return false;
00090     // c tells us if we landed inside the element or exactly on a face, edge, or node
00091     // also, dist shows the distance to the computed point.
00092     // copy parametric coords back
00093     params = r;
00094 
00095     return is_inside( params, inside_tol );
00096 }
00097 Matrix3 SpectralHex::jacobian( const CartVect& params ) const
00098 {
00099     real x_i[3];
00100     params.get( x_i );
00101     // set the positions of GL nodes, before evaluations
00102     _data.elx[0] = _xyz[0];
00103     _data.elx[1] = _xyz[1];
00104     _data.elx[2] = _xyz[2];
00105     opt_vol_set_intp_3( &_data, x_i );
00106     Matrix3 J( 0. );
00107     // it is organized differently
00108     J( 0, 0 ) = _data.jac[0];  // dx/dr
00109     J( 0, 1 ) = _data.jac[1];  // dx/ds
00110     J( 0, 2 ) = _data.jac[2];  // dx/dt
00111     J( 1, 0 ) = _data.jac[3];  // dy/dr
00112     J( 1, 1 ) = _data.jac[4];  // dy/ds
00113     J( 1, 2 ) = _data.jac[5];  // dy/dt
00114     J( 2, 0 ) = _data.jac[6];  // dz/dr
00115     J( 2, 1 ) = _data.jac[7];  // dz/ds
00116     J( 2, 2 ) = _data.jac[8];  // dz/dt
00117     return J;
00118 }
00119 void SpectralHex::evaluate_vector( const CartVect& params, const double* field, int num_tuples, double* eval ) const
00120 {
00121     // piece that we shouldn't want to cache
00122     int d;
00123     for( d = 0; d < 3; d++ )
00124     {
00125         lagrange_0( &_ld[d], params[d] );
00126     }
00127 
00128     *eval = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, field, _odwork );
00129 }
00130 void SpectralHex::integrate_vector( const double* field_values, int num_tuples, double* integral ) const
00131 {
00132     // set the position of GL points
00133     // set the positions of GL nodes, before evaluations
00134     _data.elx[0] = _xyz[0];
00135     _data.elx[1] = _xyz[1];
00136     _data.elx[2] = _xyz[2];
00137     real params[3];
00138     // triple loop; the most inner loop is in r direction, then s, then t
00139     for( int l = 0; l < num_tuples; l++ )
00140         integral[l] = 0.0;
00141     // double volume = 0;
00142     int index = 0;  // used fr the inner loop
00143     for( int k = 0; k < _n; k++ )
00144     {
00145         params[2] = _ld[2].z[k];
00146         // double wk= _ld[2].w[k];
00147         for( int j = 0; j < _n; j++ )
00148         {
00149             params[1] = _ld[1].z[j];
00150             // double wj= _ld[1].w[j];
00151             for( int i = 0; i < _n; i++ )
00152             {
00153                 params[0] = _ld[0].z[i];
00154                 // double wi= _ld[0].w[i];
00155                 opt_vol_set_intp_3( &_data, params );
00156                 double wk = _ld[2].J[k];
00157                 double wj = _ld[1].J[j];
00158                 double wi = _ld[0].J[i];
00159                 Matrix3 J( 0. );
00160                 // it is organized differently
00161                 J( 0, 0 ) = _data.jac[0];  // dx/dr
00162                 J( 0, 1 ) = _data.jac[1];  // dx/ds
00163                 J( 0, 2 ) = _data.jac[2];  // dx/dt
00164                 J( 1, 0 ) = _data.jac[3];  // dy/dr
00165                 J( 1, 1 ) = _data.jac[4];  // dy/ds
00166                 J( 1, 2 ) = _data.jac[5];  // dy/dt
00167                 J( 2, 0 ) = _data.jac[6];  // dz/dr
00168                 J( 2, 1 ) = _data.jac[7];  // dz/ds
00169                 J( 2, 2 ) = _data.jac[8];  // dz/dt
00170                 double bm = wk * wj * wi * J.determinant();
00171                 for( int l = 0; l < num_tuples; l++ )
00172                     integral[l] += bm * field_values[num_tuples * index + l];
00173                 // volume +=bm;
00174             }
00175         }
00176     }
00177     // std::cout << "volume: " << volume << "\n";
00178 }
00179 
00180 int SpectralHex::insideFcn( const double* params, const int ndim, const double tol )
00181 {
00182     return EvalSet::inside( params, ndim, tol );
00183 }
00184 
00185 // SpectralHex
00186 
00187 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines