Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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,
00079                                 const double* field,
00080                                 const int ndim,
00081                                 const int num_tuples,
00082                                 double* work,
00083                                 double* result )
00084 {
00085     // piece that we shouldn't want to cache
00086     int d = 0;
00087     for( d = 0; d < 2; d++ )
00088     {
00089         lagrange_0( &_ld[d], params[d] );
00090     }
00091     CartVect result;
00092     for( d = 0; d < 3; d++ )
00093     {
00094         result[d] = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _xyz[d], _odwork );
00095     }
00096     return result;
00097 }
00098 // replicate the functionality of hex_findpt
00099 bool SpectralQuad::reverseEvalFcn( const double* posn,
00100                                    const double* verts,
00101                                    const int nverts,
00102                                    const int ndim,
00103                                    const double iter_tol,
00104                                    const double inside_tol,
00105                                    double* work,
00106                                    double* params,
00107                                    int* is_inside )
00108 {
00109     params = init;
00110 
00111     // find nearest point
00112     double x_star[3];
00113     xyz.get( x_star );
00114 
00115     double r[2] = { 0, 0 };  // initial guess for parametric coords
00116     unsigned c  = opt_no_constraints_3;
00117     double dist = opt_findpt_2( &_data, (const double**)_xyz, x_star, r, &c );
00118     // if it did not converge, get out with throw...
00119     if( dist > 0.9e+30 )
00120     {
00121         std::vector< CartVect > dummy;
00122         throw Map::EvaluationError( CartVect( x_star ), dummy );
00123     }
00124     // c tells us if we landed inside the element or exactly on a face, edge, or node
00125     // also, dist shows the distance to the computed point.
00126     // copy parametric coords back
00127     params = r;
00128 
00129     return insideFcn( params, 2, inside_tol );
00130 }
00131 
00132 Matrix3 SpectralQuad::jacobian( const double* params,
00133                                 const double* verts,
00134                                 const int nverts,
00135                                 const int ndim,
00136                                 double* work,
00137                                 double* result )
00138 {
00139     // not implemented
00140     Matrix3 J( 0. );
00141     return J;
00142 }
00143 
00144 void SpectralQuad::evaluate_vector( const CartVect& params, const double* field, int num_tuples, double* eval ) const
00145 {
00146     // piece that we shouldn't want to cache
00147     int d;
00148     for( d = 0; d < 2; d++ )
00149     {
00150         lagrange_0( &_ld[d], params[d] );
00151     }
00152 
00153     *eval = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, field, _odwork );
00154 }
00155 void SpectralQuad::integrate_vector( const double* field,
00156                                      const double* verts,
00157                                      const int nverts,
00158                                      const int ndim,
00159                                      const int num_tuples,
00160                                      double* work,
00161                                      double* result )
00162 {
00163     // not implemented
00164 }
00165 
00166 int SpectralQuad::insideFcn( const double* params, const int ndim, const double tol )
00167 {
00168     return EvalSet::inside( params, ndim, tol );
00169 }
00170 
00171 // something we don't do for spectral hex; we do it here because
00172 //       we do not store the position of gl points in a tag yet
00173 void SpectralQuad::compute_gl_positions()
00174 {
00175     // will need to use shape functions on a simple linear quad to compute gl points
00176     // so we know the position of gl points in parametric space, we will just compute those
00177     // from the 3d vertex position (corner nodes of the quad), using simple mapping
00178     assert( this->vertex.size() == 4 );
00179     static double corner_params[4][2] = { { -1., -1. }, { 1., -1. }, { 1., 1. }, { -1., 1. } };
00180     // we will use the cached lobatto nodes in parametric space _z[d] (the same in both directions)
00181     int indexGL = 0;
00182     int n2      = _n * _n;
00183     for( int i = 0; i < _n; i++ )
00184     {
00185         double csi = _z[0][i];
00186         for( int j = 0; j < _n; j++ )
00187         {
00188             double eta = _z[1][j];  // we could really use the same _z[0] array of lobatto nodes
00189             CartVect pos( 0.0 );
00190             for( int k = 0; k < 4; k++ )
00191             {
00192                 const double N_k = ( 1 + csi * corner_params[k][0] ) * ( 1 + eta * corner_params[k][1] );
00193                 pos += N_k * vertex[k];
00194             }
00195             pos *= 0.25;                           // these are x, y, z of gl points; reorder them
00196             _glpoints[indexGL]          = pos[0];  // x
00197             _glpoints[indexGL + n2]     = pos[1];  // y
00198             _glpoints[indexGL + 2 * n2] = pos[2];  // z
00199             indexGL++;
00200         }
00201     }
00202     // now, we can set the _xyz pointers to internal memory allocated to these!
00203     _xyz[0] = &( _glpoints[0] );
00204     _xyz[1] = &( _glpoints[n2] );
00205     _xyz[2] = &( _glpoints[2 * n2] );
00206 }
00207 void SpectralQuad::get_gl_points( double*& x, double*& y, double*& z, int& size )
00208 {
00209     x    = (double*)_xyz[0];
00210     y    = (double*)_xyz[1];
00211     z    = (double*)_xyz[2];
00212     size = _n * _n;
00213 }
00214 
00215 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines