MOAB: Mesh Oriented datABase  (version 5.4.0)
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,
00077                                     CartVect& params,
00078                                     double iter_tol,
00079                                     const double inside_tol,
00080                                     const CartVect& init ) const
00081 {
00082     params = init;
00083 
00084     // find nearest point
00085     double x_star[3];
00086     xyz.get( x_star );
00087 
00088     double r[3] = { 0, 0, 0 };  // initial guess for parametric coords
00089     unsigned c  = opt_no_constraints_3;
00090     double dist = opt_findpt_3( &_data, (const double**)_xyz, x_star, r, &c );
00091     // if it did not converge, get out with throw...
00092     if( dist > 0.9e+30 ) return false;
00093     // c tells us if we landed inside the element or exactly on a face, edge, or node
00094     // also, dist shows the distance to the computed point.
00095     // copy parametric coords back
00096     params = r;
00097 
00098     return is_inside( params, inside_tol );
00099 }
00100 Matrix3 SpectralHex::jacobian( const CartVect& params ) const
00101 {
00102     real x_i[3];
00103     params.get( x_i );
00104     // set the positions of GL nodes, before evaluations
00105     _data.elx[0] = _xyz[0];
00106     _data.elx[1] = _xyz[1];
00107     _data.elx[2] = _xyz[2];
00108     opt_vol_set_intp_3( &_data, x_i );
00109     Matrix3 J( 0. );
00110     // it is organized differently
00111     J( 0, 0 ) = _data.jac[0];  // dx/dr
00112     J( 0, 1 ) = _data.jac[1];  // dx/ds
00113     J( 0, 2 ) = _data.jac[2];  // dx/dt
00114     J( 1, 0 ) = _data.jac[3];  // dy/dr
00115     J( 1, 1 ) = _data.jac[4];  // dy/ds
00116     J( 1, 2 ) = _data.jac[5];  // dy/dt
00117     J( 2, 0 ) = _data.jac[6];  // dz/dr
00118     J( 2, 1 ) = _data.jac[7];  // dz/ds
00119     J( 2, 2 ) = _data.jac[8];  // dz/dt
00120     return J;
00121 }
00122 void SpectralHex::evaluate_vector( const CartVect& params, const double* field, int num_tuples, double* eval ) const
00123 {
00124     // piece that we shouldn't want to cache
00125     int d;
00126     for( d = 0; d < 3; d++ )
00127     {
00128         lagrange_0( &_ld[d], params[d] );
00129     }
00130 
00131     *eval = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, field, _odwork );
00132 }
00133 void SpectralHex::integrate_vector( const double* field_values, int num_tuples, double* integral ) const
00134 {
00135     // set the position of GL points
00136     // set the positions of GL nodes, before evaluations
00137     _data.elx[0] = _xyz[0];
00138     _data.elx[1] = _xyz[1];
00139     _data.elx[2] = _xyz[2];
00140     real params[3];
00141     // triple loop; the most inner loop is in r direction, then s, then t
00142     for( int l = 0; l < num_tuples; l++ )
00143         integral[l] = 0.0;
00144     // double volume = 0;
00145     int index = 0;  // used fr the inner loop
00146     for( int k = 0; k < _n; k++ )
00147     {
00148         params[2] = _ld[2].z[k];
00149         // double wk= _ld[2].w[k];
00150         for( int j = 0; j < _n; j++ )
00151         {
00152             params[1] = _ld[1].z[j];
00153             // double wj= _ld[1].w[j];
00154             for( int i = 0; i < _n; i++ )
00155             {
00156                 params[0] = _ld[0].z[i];
00157                 // double wi= _ld[0].w[i];
00158                 opt_vol_set_intp_3( &_data, params );
00159                 double wk = _ld[2].J[k];
00160                 double wj = _ld[1].J[j];
00161                 double wi = _ld[0].J[i];
00162                 Matrix3 J( 0. );
00163                 // it is organized differently
00164                 J( 0, 0 ) = _data.jac[0];  // dx/dr
00165                 J( 0, 1 ) = _data.jac[1];  // dx/ds
00166                 J( 0, 2 ) = _data.jac[2];  // dx/dt
00167                 J( 1, 0 ) = _data.jac[3];  // dy/dr
00168                 J( 1, 1 ) = _data.jac[4];  // dy/ds
00169                 J( 1, 2 ) = _data.jac[5];  // dy/dt
00170                 J( 2, 0 ) = _data.jac[6];  // dz/dr
00171                 J( 2, 1 ) = _data.jac[7];  // dz/ds
00172                 J( 2, 2 ) = _data.jac[8];  // dz/dt
00173                 double bm = wk * wj * wi * J.determinant();
00174                 for( int l = 0; l < num_tuples; l++ )
00175                     integral[l] += bm * field_values[num_tuples * index + l];
00176                 // volume +=bm;
00177             }
00178         }
00179     }
00180     // std::cout << "volume: " << volume << "\n";
00181 }
00182 
00183 int SpectralHex::insideFcn( const double* params, const int ndim, const double tol )
00184 {
00185     return EvalSet::inside( params, ndim, tol );
00186 }
00187 
00188 // SpectralHex
00189 
00190 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines