SpectralHex.cpp

Go to the documentation of this file.

#include "moab/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
{
  double 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];
  double 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