#include <spectral_hex_map.hpp>

Public Types
typedef _Matrix	Matrix
Public Member Functions
	Spectral_hex_map ()
	Spectral_hex_map (int order)
	Spectral_hex_map (const Self &f)
template<typename Moab , typename Entity_handle , typename Points , typename Point >
std::pair< bool, Point >	operator() (const Moab &, const Entity_handle &, const Points &v, const Point &p, const double tol=1.e-6)
Private Types
typedef Spectral_hex_map< Matrix >	Self
Private Member Functions
void	initialize_spectral_hex (int order)
void	free_data ()
void	set_gl_points (double x, double y, double *z)
template<typename Point >
bool	is_contained (const Point &p, const double tol) const
template<typename Point , typename Points >
bool	solve_inverse (const Point &x, Point &xi, const Points &points, const double tol=1.e-6)
template<typename Point , typename Points >
Point &	evaluate (const Point &p, const Points &, Point &f)
template<typename Point , typename Field >
double	evaluate_scalar_field (const Point &p, const Field &field) const
template<typename Points , typename Field >
double	integrate_scalar_field (const Points &p, const Field &field) const
template<typename Point , typename Points >
Matrix &	jacobian (const Point &, const Points &, Matrix &J)
Private Attributes
bool	_init
int	_n
real *	_z [3]
lagrange_data	_ld [3]
opt_data_3	_data
real *	_odwork
real *	_xyz [3]

Detailed Description

template<typename _Matrix>
class moab::element_utility::Spectral_hex_map< _Matrix >

Definition at line 18 of file spectral_hex_map.hpp.

Member Typedef Documentation

template<typename _Matrix>

typedef _Matrix moab::element_utility::Spectral_hex_map< _Matrix >::Matrix

Definition at line 20 of file spectral_hex_map.hpp.

template<typename _Matrix>

typedef Spectral_hex_map< Matrix> moab::element_utility::Spectral_hex_map< _Matrix >::Self [private]

Definition at line 22 of file spectral_hex_map.hpp.

Constructor & Destructor Documentation

template<typename _Matrix>

moab::element_utility::Spectral_hex_map< _Matrix >::Spectral_hex_map ( ) [inline]

Definition at line 25 of file spectral_hex_map.hpp.

{};

template<typename _Matrix>

moab::element_utility::Spectral_hex_map< _Matrix >::Spectral_hex_map ( int order ) [inline]

Definition at line 26 of file spectral_hex_map.hpp.

{ initialize_spectral_hex( order); }

template<typename _Matrix>

moab::element_utility::Spectral_hex_map< _Matrix >::Spectral_hex_map ( const Self & f ) [inline]

Definition at line 28 of file spectral_hex_map.hpp.

{}

Member Function Documentation

template<typename _Matrix>

template<typename Point , typename Points >

Point& moab::element_utility::Spectral_hex_map< _Matrix >::evaluate	(	const Point &	p,
		const Points &	,
		Point &	f
	)		`[inline, private]`

Definition at line 138 of file spectral_hex_map.hpp.

                                                                   {
      for (int d = 0; d < 3; ++d) { lagrange_0(&_ld[ d], p[ 0]); }
      for (int d = 0; d < 3; ++d) {
        f[ d] = tensor_i3( _ld[ 0].J, _ld[ 0].n,
        _ld[1].J, _ld[1].n,
        _ld[2].J, _ld[2].n,
        _xyz[ d],
        _odwork);
      }
      return f;
    }

template<typename _Matrix>

template<typename Point , typename Field >

double moab::element_utility::Spectral_hex_map< _Matrix >::evaluate_scalar_field	(	const Point &	p,
		const Field &	field
	)		const `[inline, private]`

Definition at line 151 of file spectral_hex_map.hpp.

                                                                              {
     int d;
     for(d=0; d<3; d++){ lagrange_0(&_ld[d], p[d]); }
     return tensor_i3( _ld[0].J,_ld[0].n,
                   _ld[1].J,_ld[1].n,
                       _ld[2].J,_ld[2].n,
                   field, _odwork);
   }

template<typename _Matrix>

void moab::element_utility::Spectral_hex_map< _Matrix >::free_data ( ) [inline, private]

Definition at line 44 of file spectral_hex_map.hpp.

                    {
       for(int d=0; d<3; d++){
         free(_z[d]);
         lagrange_free(&_ld[d]);
       }
       opt_free_3(&_data);
       free(_odwork);
     }

template<typename _Matrix>

void moab::element_utility::Spectral_hex_map< _Matrix >::initialize_spectral_hex ( int order ) [inline, private]

Definition at line 30 of file spectral_hex_map.hpp.

                                            {
    if (_init && _n==order){ return; }
    if( _init && _n != order){ free_data();}
    _init = true;
    _n = order;
    for( int d = 0; d < 3; d++){
        lobatto_nodes(_z[ d], _n);
        lagrange_setup(&_ld[ d], _z[ d], _n);
    }
    opt_alloc_3(&_data, _ld);
    std::size_t nf = _n*_n, ne = _n, nw = 2*_n*_n + 3*_n;
    _odwork = tmalloc(real, 6*nf + 9*ne + nw);
    }

template<typename _Matrix>

template<typename Points , typename Field >

double moab::element_utility::Spectral_hex_map< _Matrix >::integrate_scalar_field	(	const Points &	p,
		const Field &	field
	)		const `[inline, private]`

Definition at line 160 of file spectral_hex_map.hpp.

                                              {  
   // 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];
   double xi[3];
   //triple loop; the most inner loop is in r direction, then s, then t
   double integral = 0.;
   //double volume = 0;
   int index=0; // used fr the inner loop
   for (int k=0; k<_n; k++ ) {
     xi[2]=_ld[2].z[k];
     //double wk= _ld[2].w[k];
     for (int j=0; j<_n; j++) {
       xi[1]=_ld[1].z[j];
       //double wj= _ld[1].w[j];
       for (int i=0; i<_n; i++) {
         xi[0]=_ld[0].z[i];
         //double wi= _ld[0].w[i];
         opt_vol_set_intp_3((opt_data_3 *)&_data,xi); // cast away const-ness
         double wk= _ld[2].J[k];
         double wj= _ld[1].J[j];
         double wi= _ld[0].J[i];
         Matrix3 J(0.);
         for (int n = 0; n < 8; n++)
           J(n/3, n%3) = _data.jac[n];
         double bm = wk*wj*wi* J.determinant();
         integral+= bm*field[index++];
         //volume +=bm;
       }
     }
   }
   //std::cout << "volume: " << volume << "\n";
   return integral;
 }

template<typename _Matrix>

template<typename Point >

bool moab::element_utility::Spectral_hex_map< _Matrix >::is_contained	(	const Point &	p,
		const double	tol
	)		const `[inline, private]`

Definition at line 78 of file spectral_hex_map.hpp.

                                                               {
     //just look at the box+tol here
     return ( p[0]>=-1.-tol) && (p[0]<=1.+tol) &&
            ( p[1]>=-1.-tol) && (p[1]<=1.+tol) &&
            ( p[2]>=-1.-tol) && (p[2]<=1.+tol);
    }

template<typename _Matrix>

template<typename Point , typename Points >

Matrix& moab::element_utility::Spectral_hex_map< _Matrix >::jacobian	(	const Point &	,
		const Points &	,
		Matrix &	J
	)		`[inline, private]`

Definition at line 199 of file spectral_hex_map.hpp.

                                                                   {
    real x[ 3];
    for (int i = 0; i < 3; ++i) { _data.elx[ i] = _xyz[ i]; }
    opt_vol_set_intp_3(& _data, x);
    for (int i = 0; i < 9; ++i) { J(i%3, i/3) = _data.jac[ i]; }
    return J;
  }

template<typename _Matrix>

template<typename Moab , typename Entity_handle , typename Points , typename Point >

std::pair< bool, Point> moab::element_utility::Spectral_hex_map< _Matrix >::operator()	(	const Moab &	,
		const Entity_handle &	,
		const Points &	v,
		const Point &	p,
		const double	tol = `1.e-6`
	)		`[inline]`

Definition at line 57 of file spectral_hex_map.hpp.

                                              {
        Point result(3, 0.0);
        /*
        moab.tag_get_by_ptr(_xm1Tag, &eh, 1,(const void **) &_xyz[ 0] );
        moab.tag_get_by_ptr(_ym1Tag, &eh, 1,(const void **) &_xyz[ 1] );
        moab.tag_get_by_ptr(_zm1Tag, &eh, 1,(const void **) &_xyz[ 2] );
        */
        bool point_found = solve_inverse( p, result, v, tol) &&
                  is_contained( result, tol);
        return std::make_pair( point_found, result);
    }

template<typename _Matrix>

void moab::element_utility::Spectral_hex_map< _Matrix >::set_gl_points	(	double *	x,
		double *	y,
		double *	z
	)		`[inline, private]`

Definition at line 74 of file spectral_hex_map.hpp.

                                                           {
        _xyz[ 0] = x; _xyz[ 1] = y; _xyz[ 2] = z;
    }

template<typename _Matrix>

template<typename Point , typename Points >

bool moab::element_utility::Spectral_hex_map< _Matrix >::solve_inverse	(	const Point &	x,
		Point &	xi,
		const Points &	points,
		const double	tol = `1.e-6`
	)		`[inline, private]`

Definition at line 86 of file spectral_hex_map.hpp.

                                    {
      const double error_tol_sqr = tol*tol;
      Point delta(3,0.0);
      xi = delta;
      evaluate( xi, points, delta);
      vec_subtract( delta, x);
      std::size_t num_iterations=0;
      #ifdef SPECTRAL_HEX_DEBUG
    std::stringstream ss;
    ss << "Point: "; 
       ss << x[ 0 ] << ", " << x[ 1] 
          << ", " << x [ 2] << std::endl;
    ss << "Hex: ";
    for(int i = 0; i < 8; ++i){
            ss << points[ i][ 0] << ", " << points[ i][ 1] << ", "
           << points[ i][ 2] << std::endl;
    }
    ss << std::endl;
      #endif
      while ( normsq( delta) > error_tol_sqr) {
    #ifdef SPECTRAL_HEX_DEBUG
    ss << "Iter #: "  << num_iterations 
       << " Err: " << sqrt( normsq( delta)) << " Iterate: ";
    ss << xi[ 0 ] << ", " << xi[ 1] 
        << ", " << xi[ 2] << std::endl;
    #endif
    if( ++num_iterations >= 5){ return false; }
        Matrix J;
    jacobian( xi, points, J);
        double det = moab::Matrix::determinant3( J);
        if (fabs(det) < 1.e-10){
        #ifdef SPECTRAL_HEX_DEBUG
            std::cerr << ss.str();
        #endif
        #ifndef SPECTRAL_HEX_DEBUG
        std::cerr << x[ 0 ] << ", " << x[ 1] 
              << ", " << x [ 2] << std::endl;
        #endif
        std::cerr << "inverse solve failure: det: " << det << std::endl;
        exit( -1);
    }
        vec_subtract( xi, moab::Matrix::inverse(J, 1.0/det) * delta);
        evaluate( xi, points, delta);
    vec_subtract( delta, x);
      }
       return true;
    }