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