#include <ElemEvaluator.hpp>

Public Member Functions
	EvalSet ()
	Bare constructor.
	EvalSet (EvalFcn eval, ReverseEvalFcn rev, JacobianFcn jacob, IntegrateFcn integ, InitFcn initf, InsideFcn insidef)
	Constructor.
EvalSet &	operator= (const EvalSet &eval)
	Operator=.
Static Public Member Functions
static ErrorCode	get_eval_set (Interface *mb, EntityHandle eh, EvalSet &eval_set)
	Given an entity handle, get an appropriate eval set, based on type & #vertices.
static ErrorCode	get_eval_set (EntityType tp, unsigned int num_vertices, EvalSet &eval_set)
	Given type & #vertices, get an appropriate eval set.
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)
	Common function to do reverse evaluation based on evaluation and jacobian functions.
static int	inside_function (const double *params, const int ndims, const double tol)
	Common function that returns true if params is in [-1,1]^ndims.
Public Attributes
EvalFcn	evalFcn
	Forward-evaluation of field at parametric coordinates.
ReverseEvalFcn	reverseEvalFcn
	Reverse-evaluation of parametric coordinates at physical space position.
JacobianFcn	jacobianFcn
	Evaluate the jacobian at a specified parametric position.
IntegrateFcn	integrateFcn
	Forward-evaluation of field at parametric coordinates.
InitFcn	initFcn
	Initialization function for an element.
InsideFcn	insideFcn
	Function that returns whether or not the parameters are inside the natural space of the element.

Detailed Description

Definition at line 31 of file ElemEvaluator.hpp.

Constructor & Destructor Documentation

moab::EvalSet::EvalSet ( ) [inline]

Bare constructor.

Definition at line 53 of file ElemEvaluator.hpp.

: evalFcn(NULL), reverseEvalFcn(NULL), jacobianFcn(NULL), integrateFcn(NULL), initFcn(NULL), insideFcn(NULL) {}

moab::EvalSet::EvalSet	(	EvalFcn	eval,
		ReverseEvalFcn	rev,
		JacobianFcn	jacob,
		IntegrateFcn	integ,
		InitFcn	initf,
		InsideFcn	insidef
	)		`[inline]`

Constructor.

Definition at line 56 of file ElemEvaluator.hpp.

              : evalFcn(eval), reverseEvalFcn(rev), jacobianFcn(jacob), integrateFcn(integ), initFcn(initf), insideFcn(insidef)
          {}

Member Function Documentation

ErrorCode moab::EvalSet::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]`

Common function to do reverse evaluation based on evaluation and jacobian functions.

Definition at line 18 of file ElemEvaluator.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);

        // initialize to center of element
      *cvparams = CartVect(-.4);
  
      CartVect new_pos;
        // evaluate that first guess to get a new position
      ErrorCode rval = (*eval)(cvparams->array(), verts, ndim, 
                               3, // hardwire to num_tuples to 3 since the field is coords
                               work, new_pos.array());
      if (MB_SUCCESS != rval) return rval;
      
        // residual is diff between old and new pos; need to minimize that
      CartVect res = new_pos - *cvposn;
      Matrix3 J;
      int dum, *tmp_inside = (inside ? inside : &dum);

      int iters=0;
        // while |res| larger than tol
      while (res % res > error_tol_sqr) {
        if(++iters>10) {
            // 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_FAILURE;
        }

          // get jacobian at current params
        rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
        double det = J.determinant();
        if (det < std::numeric_limits<double>::epsilon()) {
          *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 -= J.inverse(1.0/det) * res;

          // get the new forward-evaluated position, and its difference from the target pt
        rval = (*eval)(params, verts, ndim, 
                       3, // hardwire to num_tuples to 3 since the field is coords
                       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::EvalSet::get_eval_set	(	Interface *	mb,
		EntityHandle	eh,
		EvalSet &	eval_set
	)		`[inline, static]`

Given an entity handle, get an appropriate eval set, based on type & #vertices.

Definition at line 87 of file ElemEvaluator.hpp.

    {
      int nv;
      EntityType tp = mb->type_from_handle(eh);
      const EntityHandle *connect;
      std::vector<EntityHandle> dum_vec;
      ErrorCode rval = mb->get_connectivity(eh, connect, nv, false, &dum_vec);
      if (MB_SUCCESS != rval) return rval;
      
      return get_eval_set(tp, nv, eval_set);
    }

ErrorCode moab::EvalSet::get_eval_set	(	EntityType	tp,
		unsigned int	num_vertices,
		EvalSet &	eval_set
	)		`[static]`

Given type & #vertices, get an appropriate eval set.

Definition at line 90 of file ElemEvaluator.cpp.

    {
      switch (tp) {
        case MBEDGE:
            break;
        case MBTRI:
            if (LinearTri::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
            break;
        case MBQUAD:
            if (LinearQuad::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
//            if (SpectralQuad::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
            break;
        case MBTET:
            if (LinearTet::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
            break;
        case MBHEX:
            if (LinearHex::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
            if (QuadraticHex::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
//            if (SpectralHex::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
            break;
        default:
            break;
      }

      return MB_NOT_IMPLEMENTED;
    }

int moab::EvalSet::inside_function	(	const double *	params,
		const int	ndims,
		const double	tol
	)		`[static]`

Common function that returns true if params is in [-1,1]^ndims.

Definition at line 80 of file ElemEvaluator.cpp.

    {
      if (params[0] >= -1-tol && params[0] <= 1+tol &&
          (ndims < 2 || (params[1] >= -1-tol && params[1] <= 1+tol)) &&
          (ndims < 3 || (params[2] >= -1-tol && params[2] <= 1+tol))) 
        return true;
      else return false;
    }

EvalSet& moab::EvalSet::operator= ( const EvalSet & eval ) [inline]

Operator=.

Definition at line 67 of file ElemEvaluator.hpp.

                                              {
        evalFcn = eval.evalFcn;
        reverseEvalFcn = eval.reverseEvalFcn;
        jacobianFcn = eval.jacobianFcn;
        integrateFcn = eval.integrateFcn;
        initFcn = eval.initFcn;
        insideFcn = eval.insideFcn;
        return *this;
      }