moab
Static Public Member Functions | Static Protected Attributes
moab::LinearTet Class Reference

#include <LinearTet.hpp>

List of all members.

Static Public Member Functions

static ErrorCode evalFcn (const double *params, const double *field, const int ndim, const int num_tuples, double *work, double *result)
 Forward-evaluation of field at parametric coordinates.
static ErrorCode reverseEvalFcn (EvalFcn eval, JacobianFcn jacob, InsideFcn ins, 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)
 Reverse-evaluation of parametric coordinates at physical space position.
static ErrorCode jacobianFcn (const double *params, const double *verts, const int nverts, const int ndim, double *work, double *result)
 Evaluate the jacobian at a specified parametric position.
static ErrorCode integrateFcn (const double *field, const double *verts, const int nverts, const int ndim, const int num_tuples, double *work, double *result)
 Forward-evaluation of field at parametric coordinates.
static ErrorCode initFcn (const double *verts, const int nverts, double *&work)
 Initialize this EvalSet.
static int insideFcn (const double *params, const int ndim, const double tol)
 Function that returns whether or not the parameters are inside the natural space of the element.
static ErrorCode evaluate_reverse (EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f, 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 *inside)
static EvalSet eval_set ()
static bool compatible (EntityType tp, int numv, EvalSet &eset)

Static Protected Attributes

static const double corner [4][3]

Detailed Description

Definition at line 10 of file LinearTet.hpp.

Member Function Documentation

static bool moab::LinearTet::compatible ( EntityType  tp,
int  numv,
EvalSet &  eset 
) [inline, static]

Definition at line 47 of file LinearTet.hpp.

      {
        if (tp == MBTET && numv >= 4) {
          eset = eval_set();
          return true;
        }
        else return false;
      }
static EvalSet moab::LinearTet::eval_set ( ) [inline, static]

Definition at line 42 of file LinearTet.hpp.

      {
        return EvalSet(evalFcn, reverseEvalFcn, jacobianFcn, integrateFcn, initFcn, insideFcn);
      }
ErrorCode moab::LinearTet::evalFcn ( const double *  params,
const double *  field,
const int  ndim,
const int  num_tuples,
double *  work,
double *  result 
) [static]

Forward-evaluation of field at parametric coordinates.

Definition at line 35 of file LinearTet.cpp.

                                                            {
      assert(params && field && num_tuples > 0);
      std::vector<double> f0(num_tuples);
      std::copy(field, field+num_tuples, f0.begin());
      std::copy(field, field+num_tuples, result);

      for (unsigned i = 1; i < 4; ++i) {
        double p = 0.5*(params[i-1] + 1); // transform from -1 <= p <= 1 to 0 <= p <= 1
        for (int j = 0; j < num_tuples; j++)
          result[j] += (field[i*num_tuples+j]-f0[j])*p;
      }

      return MB_SUCCESS;
    }
ErrorCode moab::LinearTet::evaluate_reverse ( EvalFcn  eval,
JacobianFcn  jacob,
InsideFcn  inside_f,
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 *  inside 
) [static]

Definition at line 92 of file LinearTet.cpp.

                                                                                     {
        // TODO: should differentiate between epsilons used for
        // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
        // right now, fix the tolerance used for NR
      const double error_tol_sqr = iter_tol*iter_tol;
      CartVect *cvparams = reinterpret_cast<CartVect*>(params);
      const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);

        // find best initial guess to improve convergence
      CartVect tmp_params[] = {CartVect(-1,-1,-1), CartVect(1,-1,-1), CartVect(-1,1,-1), CartVect(-1,-1,1)};
      double resl = HUGE;
      CartVect new_pos, tmp_pos;
      ErrorCode rval;
      for (unsigned int i = 0; i < 4; i++) {
        rval = (*eval)(tmp_params[i].array(), verts, ndim, ndim, work, tmp_pos.array());
        if (MB_SUCCESS != rval) return rval;
        double tmp_resl = (tmp_pos-*cvposn).length_squared();
        if (tmp_resl < resl) {
          *cvparams = tmp_params[i];
          new_pos = tmp_pos;
          resl = tmp_resl;
        }        
      }

        // residual is diff between old and new pos; need to minimize that
      CartVect res = new_pos - *cvposn;
      Matrix3 J;
      rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
      double det = J.determinant();
      assert(det > std::numeric_limits<double>::epsilon());
      Matrix3 Ji = J.inverse(1.0/det);

      int iters=0;
        // while |res| larger than tol
      int dum, *tmp_inside = (inside ? inside : &dum);
      while (res % res > error_tol_sqr) {
        if(++iters>25) {
            // if we haven't converged but we're outside, that's defined as success
          *tmp_inside = (*inside_f)(params, ndim, inside_tol);
          if (!(*tmp_inside)) return MB_SUCCESS;
          else return MB_INDEX_OUT_OF_RANGE;
        }
        
          // new params tries to eliminate residual
        *cvparams -= Ji * res;

          // get the new forward-evaluated position, and its difference from the target pt
        rval = (*eval)(params, verts, ndim, ndim, work, new_pos.array());
        if (MB_SUCCESS != rval) return rval;
        res = new_pos - *cvposn;
      }

      if (inside)
        *inside = (*inside_f)(params, ndim, inside_tol);

      return MB_SUCCESS;
    }// Map::evaluate_reverse()
ErrorCode moab::LinearTet::initFcn ( const double *  verts,
const int  nverts,
double *&  work 
) [static]

Initialize this EvalSet.

Definition at line 13 of file LinearTet.cpp.

                                                                                {
        // allocate work array as: 
        // work[0..8] = T
        // work[9..17] = Tinv
        // work[18] = detT
        // work[19] = detTinv
      assert(!work && verts);
      work = new double[20];
      Matrix3 *T = reinterpret_cast<Matrix3*>(work),
          *Tinv = reinterpret_cast<Matrix3*>(work+9);
      double *detT = work+18, *detTinv = work+19;
      
      *T = Matrix3(verts[1*3+0]-verts[0*3+0],verts[2*3+0]-verts[0*3+0],verts[3*3+0]-verts[0*3+0],
                   verts[1*3+1]-verts[0*3+1],verts[2*3+1]-verts[0*3+1],verts[3*3+1]-verts[0*3+1],
                   verts[1*3+2]-verts[0*3+2],verts[2*3+2]-verts[0*3+2],verts[3*3+2]-verts[0*3+2]);
      *Tinv = T->inverse();
      *detT = T->determinant();
      *detTinv = (0.0 == *detT ? HUGE : 1.0 / *detT);

      return MB_SUCCESS;
    }
int moab::LinearTet::insideFcn ( const double *  params,
const int  ndim,
const double  tol 
) [static]

Function that returns whether or not the parameters are inside the natural space of the element.

Definition at line 85 of file LinearTet.cpp.

    {
      return (params[0] >= -1.0-tol && params[1] >= -1.0-tol && params[2] >= -1.0-tol && 
              params[0] + params[1] + params[2] <= 1.0+tol);
      
    }
ErrorCode moab::LinearTet::integrateFcn ( const double *  field,
const double *  verts,
const int  nverts,
const int  ndim,
const int  num_tuples,
double *  work,
double *  result 
) [static]

Forward-evaluation of field at parametric coordinates.

Definition at line 51 of file LinearTet.cpp.

    {
      assert(field && num_tuples > 0);
      std::fill(result, result+num_tuples, 0.0);
      for(int i = 0; i < nverts; ++i) {
        for (int j = 0; j < num_tuples; j++)
          result[j] += field[i*num_tuples+j];
      }
      double tmp = work[18]/24.0;
      for (int i = 0; i < num_tuples; i++) result[i] *= tmp;

      return MB_SUCCESS;
    }
ErrorCode moab::LinearTet::jacobianFcn ( const double *  params,
const double *  verts,
const int  nverts,
const int  ndim,
double *  work,
double *  result 
) [static]

Evaluate the jacobian at a specified parametric position.

Definition at line 66 of file LinearTet.cpp.

    {
        // jacobian is cached in work array
      assert(work);
      std::copy(work, work+9, result);
      return MB_SUCCESS;
    }
ErrorCode moab::LinearTet::reverseEvalFcn ( EvalFcn  eval,
JacobianFcn  jacob,
InsideFcn  ins,
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 
) [static]

Reverse-evaluation of parametric coordinates at physical space position.

Definition at line 75 of file LinearTet.cpp.

    {
      assert(posn && verts);
      return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, 
                              work, params, is_inside);
    }

Member Data Documentation

const double moab::LinearTet::corner [static, protected]
Initial value:
 { {0,0,0},
                                             {1,0,0},
                                             {0,1,0},
                                             {0,0,1}}

Definition at line 58 of file LinearTet.hpp.

The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines