SpectralQuad.cpp

Go to the documentation of this file.

#include "moab/LocalDiscretization/SpectralQuad.hpp"
#include "moab/Forward.hpp"
 
namespace moab
{
 
// filescope for static member data that is cached
int SpectralQuad::_n;
real* SpectralQuad::_z[2];
lagrange_data SpectralQuad::_ld[2];
opt_data_2 SpectralQuad::_data;
real* SpectralQuad::_odwork;
real* SpectralQuad::_glpoints;
bool SpectralQuad::_init = false;
 
SpectralQuad::SpectralQuad() : Map( 0 ) {}
// the preferred constructor takes pointers to GL blocked positions
SpectralQuad::SpectralQuad( int order, double* x, double* y, double* z ) : Map( 0 )
{
 Init( order );
 _xyz[0] = x;
 _xyz[1] = y;
 _xyz[2] = z;
}
SpectralQuad::SpectralQuad( int order ) : Map( 4 )
{
 Init( order );
 _xyz[0] = _xyz[1] = _xyz[2] = NULL;
}
SpectralQuad::~SpectralQuad()
{
 if( _init ) freedata();
 _init = false;
}
void SpectralQuad::Init( int order )
{
 if( _init && _n == order ) return;
 if( _init && _n != order )
 {
 // TODO: free data cached
 freedata();
 }
 // compute stuff that depends only on order
 _init = true;
 _n = order;
 // duplicates! n is the same in all directions !!!
 for( int d = 0; d < 2; d++ )
 {
 _z[d] = tmalloc( double, _n );
 lobatto_nodes( _z[d], _n );
 lagrange_setup( &_ld[d], _z[d], _n );
 }
 opt_alloc_2( &_data, _ld );
 
 unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
 _odwork = tmalloc( double, 6 * nf + 9 * ne + nw );
 _glpoints = tmalloc( double, 3 * nf );
}
 
void SpectralQuad::freedata()
{
 for( int d = 0; d < 2; d++ )
 {
 free( _z[d] );
 lagrange_free( &_ld[d] );
 }
 opt_free_2( &_data );
 free( _odwork );
 free( _glpoints );
}
 
void SpectralQuad::set_gl_points( double* x, double* y, double* z )
{
 _xyz[0] = x;
 _xyz[1] = y;
 _xyz[2] = z;
}
CartVect SpectralQuad::evalFcn( const double* params,
 const double* field,
 const int ndim,
 const int num_tuples,
 double* work,
 double* result )
{
 // piece that we shouldn't want to cache
 int d = 0;
 for( d = 0; d < 2; d++ )
 {
 lagrange_0( &_ld[d], params[d] );
 }
 CartVect result;
 for( d = 0; d < 3; d++ )
 {
 result[d] = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _xyz[d], _odwork );
 }
 return result;
}
// replicate the functionality of hex_findpt
bool SpectralQuad::reverseEvalFcn( const double* posn,
 const double* verts,
 const int nverts,
 const int ndim,
 const double iter_tol,
 const double inside_tol,
 double* work,
 double* params,
 int* is_inside )
{
 params = init;
 
 // find nearest point
 double x_star[3];
 xyz.get( x_star );
 
 double r[2] = { 0, 0 }; // initial guess for parametric coords
 unsigned c = opt_no_constraints_3;
 double dist = opt_findpt_2( &_data, (const double**)_xyz, x_star, r, &c );
 // if it did not converge, get out with throw...
 if( dist > 0.9e+30 )
 {
 std::vector< CartVect > dummy;
 throw Map::EvaluationError( CartVect( x_star ), dummy );
 }
 // c tells us if we landed inside the element or exactly on a face, edge, or node
 // also, dist shows the distance to the computed point.
 // copy parametric coords back
 params = r;
 
 return insideFcn( params, 2, inside_tol );
}
 
Matrix3 SpectralQuad::jacobian( const double* params,
 const double* verts,
 const int nverts,
 const int ndim,
 double* work,
 double* result )
{
 // not implemented
 Matrix3 J( 0. );
 return J;
}
 
void SpectralQuad::evaluate_vector( const CartVect& params, const double* field, int num_tuples, double* eval ) const
{
 // piece that we shouldn't want to cache
 int d;
 for( d = 0; d < 2; d++ )
 {
 lagrange_0( &_ld[d], params[d] );
 }
 
 *eval = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, field, _odwork );
}
void SpectralQuad::integrate_vector( const double* field,
 const double* verts,
 const int nverts,
 const int ndim,
 const int num_tuples,
 double* work,
 double* result )
{
 // not implemented
}
 
int SpectralQuad::insideFcn( const double* params, const int ndim, const double tol )
{
 return EvalSet::inside( params, ndim, tol );
}
 
// something we don't do for spectral hex; we do it here because
// we do not store the position of gl points in a tag yet
void SpectralQuad::compute_gl_positions()
{
 // will need to use shape functions on a simple linear quad to compute gl points
 // so we know the position of gl points in parametric space, we will just compute those
 // from the 3d vertex position (corner nodes of the quad), using simple mapping
 assert( this->vertex.size() == 4 );
 static double corner_params[4][2] = { { -1., -1. }, { 1., -1. }, { 1., 1. }, { -1., 1. } };
 // we will use the cached lobatto nodes in parametric space _z[d] (the same in both directions)
 int indexGL = 0;
 int n2 = _n * _n;
 for( int i = 0; i < _n; i++ )
 {
 double csi = _z[0][i];
 for( int j = 0; j < _n; j++ )
 {
 double eta = _z[1][j]; // we could really use the same _z[0] array of lobatto nodes
 CartVect pos( 0.0 );
 for( int k = 0; k < 4; k++ )
 {
 const double N_k = ( 1 + csi * corner_params[k][0] ) * ( 1 + eta * corner_params[k][1] );
 pos += N_k * vertex[k];
 }
 pos *= 0.25; // these are x, y, z of gl points; reorder them
 _glpoints[indexGL] = pos[0]; // x
 _glpoints[indexGL + n2] = pos[1]; // y
 _glpoints[indexGL + 2 * n2] = pos[2]; // z
 indexGL++;
 }
 }
 // now, we can set the _xyz pointers to internal memory allocated to these!
 _xyz[0] = &( _glpoints[0] );
 _xyz[1] = &( _glpoints[n2] );
 _xyz[2] = &( _glpoints[2 * n2] );
}
void SpectralQuad::get_gl_points( double*& x, double*& y, double*& z, int& size )
{
 x = (double*)_xyz[0];
 y = (double*)_xyz[1];
 z = (double*)_xyz[2];
 size = _n * _n;
}
 
} // namespace moab