MOAB: Mesh Oriented datABase  (version 5.2.1)
SpectralQuad.cpp
Go to the documentation of this file.
00001 #include "moab/LocalDiscretization/SpectralQuad.hpp"
00002 #include "moab/Forward.hpp"
00003 
00004 namespace moab
00005 {
00006 
00007 // filescope for static member data that is cached
00008 int SpectralQuad::_n;
00009 real* SpectralQuad::_z[2];
00010 lagrange_data SpectralQuad::_ld[2];
00011 opt_data_2 SpectralQuad::_data;
00012 real* SpectralQuad::_odwork;
00013 real* SpectralQuad::_glpoints;
00014 bool SpectralQuad::_init = false;
00015 
00016 SpectralQuad::SpectralQuad() : Map( 0 ) {}
00017 // the preferred constructor takes pointers to GL blocked positions
00018 SpectralQuad::SpectralQuad( int order, double* x, double* y, double* z ) : Map( 0 )
00019 {
00020     Init( order );
00021     _xyz[0] = x;
00022     _xyz[1] = y;
00023     _xyz[2] = z;
00024 }
00025 SpectralQuad::SpectralQuad( int order ) : Map( 4 )
00026 {
00027     Init( order );
00028     _xyz[0] = _xyz[1] = _xyz[2] = NULL;
00029 }
00030 SpectralQuad::~SpectralQuad()
00031 {
00032     if( _init ) freedata();
00033     _init = false;
00034 }
00035 void SpectralQuad::Init( int order )
00036 {
00037     if( _init && _n == order ) return;
00038     if( _init && _n != order )
00039     {
00040         // TODO: free data cached
00041         freedata();
00042     }
00043     // compute stuff that depends only on order
00044     _init = true;
00045     _n    = order;
00046     // duplicates! n is the same in all directions !!!
00047     for( int d = 0; d < 2; d++ )
00048     {
00049         _z[d] = tmalloc( double, _n );
00050         lobatto_nodes( _z[d], _n );
00051         lagrange_setup( &_ld[d], _z[d], _n );
00052     }
00053     opt_alloc_2( &_data, _ld );
00054 
00055     unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
00056     _odwork   = tmalloc( double, 6 * nf + 9 * ne + nw );
00057     _glpoints = tmalloc( double, 3 * nf );
00058 }
00059 
00060 void SpectralQuad::freedata()
00061 {
00062     for( int d = 0; d < 2; d++ )
00063     {
00064         free( _z[d] );
00065         lagrange_free( &_ld[d] );
00066     }
00067     opt_free_2( &_data );
00068     free( _odwork );
00069     free( _glpoints );
00070 }
00071 
00072 void SpectralQuad::set_gl_points( double* x, double* y, double* z )
00073 {
00074     _xyz[0] = x;
00075     _xyz[1] = y;
00076     _xyz[2] = z;
00077 }
00078 CartVect SpectralQuad::evalFcn( const double* params, const double* field, const int ndim, const int num_tuples,
00079                                 double* work, double* result )
00080 {
00081     // piece that we shouldn't want to cache
00082     int d = 0;
00083     for( d = 0; d < 2; d++ )
00084     {
00085         lagrange_0( &_ld[d], params[d] );
00086     }
00087     CartVect result;
00088     for( d = 0; d < 3; d++ )
00089     {
00090         result[d] = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _xyz[d], _odwork );
00091     }
00092     return result;
00093 }
00094 // replicate the functionality of hex_findpt
00095 bool SpectralQuad::reverseEvalFcn( const double* posn, const double* verts, const int nverts, const int ndim,
00096                                    const double iter_tol, const double inside_tol, double* work, double* params,
00097                                    int* is_inside )
00098 {
00099     params = init;
00100 
00101     // find nearest point
00102     double x_star[3];
00103     xyz.get( x_star );
00104 
00105     double r[2] = { 0, 0 };  // initial guess for parametric coords
00106     unsigned c  = opt_no_constraints_3;
00107     double dist = opt_findpt_2( &_data, (const double**)_xyz, x_star, r, &c );
00108     // if it did not converge, get out with throw...
00109     if( dist > 0.9e+30 )
00110     {
00111         std::vector< CartVect > dummy;
00112         throw Map::EvaluationError( CartVect( x_star ), dummy );
00113     }
00114     // c tells us if we landed inside the element or exactly on a face, edge, or node
00115     // also, dist shows the distance to the computed point.
00116     // copy parametric coords back
00117     params = r;
00118 
00119     return insideFcn( params, 2, inside_tol );
00120 }
00121 
00122 Matrix3 SpectralQuad::jacobian( const double* params, const double* verts, const int nverts, const int ndim,
00123                                 double* work, double* result )
00124 {
00125     // not implemented
00126     Matrix3 J( 0. );
00127     return J;
00128 }
00129 
00130 void SpectralQuad::evaluate_vector( const CartVect& params, const double* field, int num_tuples, double* eval ) const
00131 {
00132     // piece that we shouldn't want to cache
00133     int d;
00134     for( d = 0; d < 2; d++ )
00135     {
00136         lagrange_0( &_ld[d], params[d] );
00137     }
00138 
00139     *eval = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, field, _odwork );
00140 }
00141 void SpectralQuad::integrate_vector( const double* field, const double* verts, const int nverts, const int ndim,
00142                                      const int num_tuples, double* work, double* result )
00143 {
00144     // not implemented
00145 }
00146 
00147 int SpectralQuad::insideFcn( const double* params, const int ndim, const double tol )
00148 {
00149     return EvalSet::inside( params, ndim, tol );
00150 }
00151 
00152 // something we don't do for spectral hex; we do it here because
00153 //       we do not store the position of gl points in a tag yet
00154 void SpectralQuad::compute_gl_positions()
00155 {
00156     // will need to use shape functions on a simple linear quad to compute gl points
00157     // so we know the position of gl points in parametric space, we will just compute those
00158     // from the 3d vertex position (corner nodes of the quad), using simple mapping
00159     assert( this->vertex.size() == 4 );
00160     static double corner_params[4][2] = { { -1., -1. }, { 1., -1. }, { 1., 1. }, { -1., 1. } };
00161     // we will use the cached lobatto nodes in parametric space _z[d] (the same in both directions)
00162     int indexGL = 0;
00163     int n2      = _n * _n;
00164     for( int i = 0; i < _n; i++ )
00165     {
00166         double csi = _z[0][i];
00167         for( int j = 0; j < _n; j++ )
00168         {
00169             double eta = _z[1][j];  // we could really use the same _z[0] array of lobatto nodes
00170             CartVect pos( 0.0 );
00171             for( int k = 0; k < 4; k++ )
00172             {
00173                 const double N_k = ( 1 + csi * corner_params[k][0] ) * ( 1 + eta * corner_params[k][1] );
00174                 pos += N_k * vertex[k];
00175             }
00176             pos *= 0.25;                           // these are x, y, z of gl points; reorder them
00177             _glpoints[indexGL]          = pos[0];  // x
00178             _glpoints[indexGL + n2]     = pos[1];  // y
00179             _glpoints[indexGL + 2 * n2] = pos[2];  // z
00180             indexGL++;
00181         }
00182     }
00183     // now, we can set the _xyz pointers to internal memory allocated to these!
00184     _xyz[0] = &( _glpoints[0] );
00185     _xyz[1] = &( _glpoints[n2] );
00186     _xyz[2] = &( _glpoints[2 * n2] );
00187 }
00188 void SpectralQuad::get_gl_points( double*& x, double*& y, double*& z, int& size )
00189 {
00190     x    = (double*)_xyz[0];
00191     y    = (double*)_xyz[1];
00192     z    = (double*)_xyz[2];
00193     size = _n * _n;
00194 }
00195 
00196 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines