SpectralHex.cpp

Go to the documentation of this file.

#include "moab/LocalDiscretization/SpectralHex.hpp"
#include "moab/Forward.hpp"
 
namespace moab
{
 
// SpectralHex
 
// filescope for static member data that is cached
int SpectralHex::_n;
real* SpectralHex::_z[3];
lagrange_data SpectralHex::_ld[3];
opt_data_3 SpectralHex::_data;
real* SpectralHex::_odwork;
 
bool SpectralHex::_init = false;
 
void SpectralHex::initFcn( const double* verts, const int nverts, double*& work )
{
    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;
    // triplicates! n is the same in all directions !!!
    for( int d = 0; d < 3; d++ )
    {
        _z[d] = tmalloc( double, _n );
        lobatto_nodes( _z[d], _n );
        lagrange_setup( &_ld[d], _z[d], _n );
    }
    opt_alloc_3( &_data, _ld );
 
    unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
    _odwork = tmalloc( double, 6 * nf + 9 * ne + nw );
}
void SpectralHex::freedata()
{
    for( int d = 0; d < 3; d++ )
    {
        free( _z[d] );
        lagrange_free( &_ld[d] );
    }
    opt_free_3( &_data );
    free( _odwork );
}
 
void SpectralHex::set_gl_points( double* x, double* y, double* z )
{
    _xyz[0] = x;
    _xyz[1] = y;
    _xyz[2] = z;
}
CartVect SpectralHex::evaluate( const CartVect& params ) const
{
    // piece that we shouldn't want to cache
    int d = 0;
    for( d = 0; d < 3; d++ )
    {
        lagrange_0( &_ld[d], params[d] );
    }
    CartVect result;
    for( d = 0; d < 3; d++ )
    {
        result[d] = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n,
                               _xyz[d],  // this is the "field"
                               _odwork );
    }
    return result;
}
// replicate the functionality of hex_findpt
bool SpectralHex::evaluate_reverse( CartVect const& xyz,
                                    CartVect& params,
                                    double iter_tol,
                                    const double inside_tol,
                                    const CartVect& init ) const
{
    params = init;
 
    // find nearest point
    double x_star[3];
    xyz.get( x_star );
 
    double r[3] = { 0, 0, 0 };  // initial guess for parametric coords
    unsigned c  = opt_no_constraints_3;
    double dist = opt_findpt_3( &_data, (const double**)_xyz, x_star, r, &c );
    // if it did not converge, get out with throw...
    if( dist > 0.9e+30 ) return false;
    // 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 is_inside( params, inside_tol );
}
Matrix3 SpectralHex::jacobian( const CartVect& params ) const
{
    real x_i[3];
    params.get( x_i );
    // set the positions of GL nodes, before evaluations
    _data.elx[0] = _xyz[0];
    _data.elx[1] = _xyz[1];
    _data.elx[2] = _xyz[2];
    opt_vol_set_intp_3( &_data, x_i );
    Matrix3 J( 0. );
    // it is organized differently
    J( 0, 0 ) = _data.jac[0];  // dx/dr
    J( 0, 1 ) = _data.jac[1];  // dx/ds
    J( 0, 2 ) = _data.jac[2];  // dx/dt
    J( 1, 0 ) = _data.jac[3];  // dy/dr
    J( 1, 1 ) = _data.jac[4];  // dy/ds
    J( 1, 2 ) = _data.jac[5];  // dy/dt
    J( 2, 0 ) = _data.jac[6];  // dz/dr
    J( 2, 1 ) = _data.jac[7];  // dz/ds
    J( 2, 2 ) = _data.jac[8];  // dz/dt
    return J;
}
void SpectralHex::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 < 3; d++ )
    {
        lagrange_0( &_ld[d], params[d] );
    }
 
    *eval = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, field, _odwork );
}
void SpectralHex::integrate_vector( const double* field_values, int num_tuples, double* integral ) const
{
    // set the position of GL points
    // set the positions of GL nodes, before evaluations
    _data.elx[0] = _xyz[0];
    _data.elx[1] = _xyz[1];
    _data.elx[2] = _xyz[2];
    real params[3];
    // triple loop; the most inner loop is in r direction, then s, then t
    for( int l = 0; l < num_tuples; l++ )
        integral[l] = 0.0;
    // double volume = 0;
    int index = 0;  // used fr the inner loop
    for( int k = 0; k < _n; k++ )
    {
        params[2] = _ld[2].z[k];
        // double wk= _ld[2].w[k];
        for( int j = 0; j < _n; j++ )
        {
            params[1] = _ld[1].z[j];
            // double wj= _ld[1].w[j];
            for( int i = 0; i < _n; i++ )
            {
                params[0] = _ld[0].z[i];
                // double wi= _ld[0].w[i];
                opt_vol_set_intp_3( &_data, params );
                double wk = _ld[2].J[k];
                double wj = _ld[1].J[j];
                double wi = _ld[0].J[i];
                Matrix3 J( 0. );
                // it is organized differently
                J( 0, 0 ) = _data.jac[0];  // dx/dr
                J( 0, 1 ) = _data.jac[1];  // dx/ds
                J( 0, 2 ) = _data.jac[2];  // dx/dt
                J( 1, 0 ) = _data.jac[3];  // dy/dr
                J( 1, 1 ) = _data.jac[4];  // dy/ds
                J( 1, 2 ) = _data.jac[5];  // dy/dt
                J( 2, 0 ) = _data.jac[6];  // dz/dr
                J( 2, 1 ) = _data.jac[7];  // dz/ds
                J( 2, 2 ) = _data.jac[8];  // dz/dt
                double bm = wk * wj * wi * J.determinant();
                for( int l = 0; l < num_tuples; l++ )
                    integral[l] += bm * field_values[num_tuples * index + l];
                // volume +=bm;
            }
        }
    }
    // std::cout << "volume: " << volume << "\n";
}
 
int SpectralHex::insideFcn( const double* params, const int ndim, const double tol )
{
    return EvalSet::inside( params, ndim, tol );
}
 
// SpectralHex
 
}  // namespace moab