ElemEvaluator.hpp

Go to the documentation of this file.

#ifndef ELEM_EVALUATOR_HPP
#define ELEM_EVALUATOR_HPP
 
#include "moab/Interface.hpp"
#include "moab/CartVect.hpp"
#include "moab/Matrix3.hpp"
#include "moab/CN.hpp"
 
#include <vector>
 
namespace moab
{
 
typedef ErrorCode ( *EvalFcn )( const double* params,
 const double* field,
 const int ndim,
 const int num_tuples,
 double* work,
 double* result );
 
typedef ErrorCode ( *JacobianFcn )( const double* params,
 const double* verts,
 const int nverts,
 const int ndim,
 double* work,
 double* result );
 
typedef ErrorCode ( *IntegrateFcn )( const double* field,
 const double* verts,
 const int nverts,
 const int ndim,
 const int num_tuples,
 double* work,
 double* result );
 
typedef ErrorCode ( *InitFcn )( const double* verts, const int nverts, double*& work );
 
typedef int ( *InsideFcn )( const double* verts, const int ndims, const double tol );
 
typedef 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 );
 
typedef ErrorCode (
 *NormalFcn )( const int ientDim, const int facet, const int nverts, const double* verts, double normal[3] );
 
class EvalSet
{
 public:
 /** \brief Forward-evaluation of field at parametric coordinates */
 EvalFcn evalFcn;
 
 /** \brief Reverse-evaluation of parametric coordinates at physical space position */
 ReverseEvalFcn reverseEvalFcn;
 
 /** \brief Evaluate the normal at a local facet (edge/face for 2D/3D) */
 NormalFcn normalFcn;
 
 /** \brief Evaluate the jacobian at a specified parametric position */
 JacobianFcn jacobianFcn;
 
 /** \brief Forward-evaluation of field at parametric coordinates */
 IntegrateFcn integrateFcn;
 
 /** \brief Initialization function for an element */
 InitFcn initFcn;
 
 /** \brief Function that returns whether or not the parameters are inside the natural space of
 * the element */
 InsideFcn insideFcn;
 
 /** \brief Bare constructor */
 EvalSet()
 : evalFcn( NULL ), reverseEvalFcn( NULL ), normalFcn( NULL ), jacobianFcn( NULL ), integrateFcn( NULL ),
 initFcn( NULL ), insideFcn( NULL )
 {
 }
 
 /** \brief Constructor */
 EvalSet( EvalFcn eval,
 ReverseEvalFcn rev,
 NormalFcn normal,
 JacobianFcn jacob,
 IntegrateFcn integ,
 InitFcn initf,
 InsideFcn insidef )
 : evalFcn( eval ), reverseEvalFcn( rev ), normalFcn( normal ), jacobianFcn( jacob ), integrateFcn( integ ),
 initFcn( initf ), insideFcn( insidef )
 {
 }
 
 /** \brief Copy constructor */
 EvalSet( EvalSet const& eval )
 : evalFcn( eval.evalFcn ), reverseEvalFcn( eval.reverseEvalFcn ), normalFcn( eval.normalFcn ),
 jacobianFcn( eval.jacobianFcn ), integrateFcn( eval.integrateFcn ), initFcn( eval.initFcn ),
 insideFcn( eval.insideFcn )
 {
 }
 
 /** \brief Given an entity handle, get an appropriate eval set, based on type & #vertices */
 static ErrorCode get_eval_set( Interface* mb, EntityHandle eh, EvalSet& eval_set );
 
 /** \brief Given type & #vertices, get an appropriate eval set */
 static ErrorCode get_eval_set( EntityType tp, unsigned int num_vertices, EvalSet& eval_set );
 
 /** \brief Operator= */
 EvalSet& operator=( const EvalSet& eval )
 {
 evalFcn = eval.evalFcn;
 reverseEvalFcn = eval.reverseEvalFcn;
 normalFcn = eval.normalFcn;
 jacobianFcn = eval.jacobianFcn;
 integrateFcn = eval.integrateFcn;
 initFcn = eval.initFcn;
 insideFcn = eval.insideFcn;
 return *this;
 }
 
 /** \brief Common function to do reverse evaluation based on evaluation and jacobian functions
 */
 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 );
 /** \brief Common function that returns true if params is in [-1,1]^ndims */
 static int inside_function( const double* params, const int ndims, const double tol );
};
 
/** \brief Given an entity handle, get an appropriate eval set, based on type & #vertices */
inline ErrorCode EvalSet::get_eval_set( Interface* mb, EntityHandle eh, EvalSet& eval_set )
{
 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 );
}
 
/**\brief Class facilitating local discretization-related functions
 * \class ElemEvaluator
 * This class implements discretization-related functionality operating
 * on data in MOAB. A member of this class caches certain data about the element
 * it's currently operating on, but is not meant to be instantiated one-per-element,
 * but rather one-per-search (or other operation on a collection of elements).
 *
 * Actual discretization functionality is accessed through function pointers,
 * allowing applications to specialize the implementation of specific functions
 * while still using this class.
 *
 * This class depends on MOAB functionality for gathering entity-based data; the functions
 * it calls through function pointers depend only on POD (plain old data, or intrinsic data types).
 * This allows the use of other packages for serving these functions, without having to modify
 * them to get data through MOAB. This should also promote efficiency, since in many cases they
 * will be able to read data from its native storage locations.
 */
 
class ElemEvaluator
{
 public:
 /** \brief Constructor
 * \param impl MOAB instance
 * \param ent Entity handle to cache on the evaluator; if non-zero, calls set_ent_handle, which
 * does some other stuff. \param tag Tag to cache on the evaluator; if non-zero, calls
 * set_tag_handle, which does some other stuff too. \param tagged_ent_dim Dimension of entities
 * to be tagged to cache on the evaluator
 */
 ElemEvaluator( Interface* impl, EntityHandle ent = 0, Tag tag = 0, int tagged_ent_dim = -1 );
 ~ElemEvaluator();
 
 /** \brief Evaluate cached tag at a given parametric location within the cached entity
 * If evaluating coordinates, call set_tag(0, 0), which indicates coords instead of a tag.
 * \param params Parameters at which to evaluate field
 * \param result Result of evaluation
 * \param num_tuples Size of evaluated field, in number of values
 */
 ErrorCode eval( const double* params, double* result, int num_tuples = -1 ) const;
 
 /** \brief Reverse-evaluate the cached entity at a given physical position
 * \param posn Position at which to evaluate parameters
 * \param iter_tol Tolerance of reverse evaluation non-linear iteration, usually 10^-10 or so
 * \param inside_tol Tolerance of is_inside evaluation, usually 10^-6 or so
 * \param params Result of evaluation
 * \param is_inside If non-NULL, returns true of resulting parameters place the point inside the
 * element (in most cases, within [-1]*(dim)
 */
 ErrorCode reverse_eval( const double* posn,
 double iter_tol,
 double inside_tol,
 double* params,
 int* is_inside = NULL ) const;
 
 /**
 * \brief Evaluate the normal to a facet of an entity
 * \param ientDim Dimension of the facet. Should be (d-1) for d-dimensional entities
 * \param facet Local id of the facet w.r.t the entity
 * \param normal Returns the normal.
 */
 ErrorCode get_normal( const int ientDim, const int facet, double normal[3] ) const;
 
 /** \brief Evaluate the jacobian of the cached entity at a given parametric location
 * \param params Parameters at which to evaluate jacobian
 * \param result Result of evaluation, in the form of a 3x3 matrix, stored in column-major order
 */
 ErrorCode jacobian( const double* params, double* result ) const;
 
 /** \brief Integrate the cached tag over the cached entity
 * \param result Result of the integration
 */
 ErrorCode integrate( double* result ) const;
 
 /** \brief Return whether a physical position is inside the cached entity to within a tolerance
 * \param params Parameters at which to query the element
 * \param tol Tolerance, usually 10^-6 or so
 */
 int inside( const double* params, const double tol ) const;
 
 /** \brief Given a list of entities, return the entity the point is in, or none
 * This function reverse-evaluates the entities, returning the first entity containing the
 * point. If no entity contains the point, containing_ent is returned as 0 and params are
 * unchanged. This function returns something other than MB_SUCCESS only when queries go wrong
 * for some reason. num_evals, if non-NULL, is always incremented for each call to reverse_eval.
 * This function calls set_ent_handle for each entity before calling reverse_eval, so the
 * ElemEvaluator object is changed. \param entities Entities tested \param point Point tested,
 * must have 3 dimensions, even for edge and face entities \param iter_tol Tolerance for
 * non-linear reverse evaluation \param inside_tol Tolerance for is_inside test \param
 * containing_ent Entity containing the point, returned 0 if no entity \param params Parameters
 * of point in containing entity, unchanged if no containing entity \param num_evals If
 * non-NULL, incremented each time reverse_eval is called \return Returns non-success only if
 * evaulation failed for some reason (point not in element is NOT a reason for failure)
 */
 ErrorCode find_containing_entity( Range& entities,
 const double* point,
 const double iter_tol,
 const double inside_tol,
 EntityHandle& containing_ent,
 double* params,
 unsigned int* num_evals = NULL );
 
 /** \brief Given an entity set, return the contained entity the point is in, or none
 * This function reverse-evaluates the entities, returning the first entity containing the
 * point. If no entity contains the point, containing_ent is returned as 0 and params are
 * unchanged. This function returns something other than MB_SUCCESS only when queries go wrong
 * for some reason. num_evals, if non-NULL, is always incremented for each call to reverse_eval.
 * This function calls set_ent_handle for each entity before calling reverse_eval, so the
 * ElemEvaluator object is changed. \param ent_set Entity set containing the entities to be
 * tested \param point Point tested, must have 3 dimensions, even for edge and face entities
 * \param iter_tol Tolerance for non-linear reverse evaluation
 * \param inside_tol Tolerance for is_inside test
 * \param containing_ent Entity containing the point, returned 0 if no entity
 * \param params Parameters of point in containing entity, unchanged if no containing entity
 * \param num_evals If non-NULL, incremented each time reverse_eval is called
 * \return Returns non-success only if evaulation failed for some reason (point not in element
 * is NOT a reason for failure)
 */
 ErrorCode find_containing_entity( EntityHandle ent_set,
 const double* point,
 const double iter_tol,
 const double inside_tol,
 EntityHandle& containing_ent,
 double* params,
 unsigned int* num_evals = NULL );
 
 /** \brief Set the eval set for a given type entity
 * \param tp Entity type for which to set the eval set
 * \param eval_set Eval set object to use
 */
 ErrorCode set_eval_set( EntityType tp, const EvalSet& eval_set );
 
 /** \brief Set the eval set using a given entity handle
 * This function queries the entity type and number of vertices on the entity to decide which
 * type of shape function to use. NOTE: this function should only be called after setting a
 * valid MOAB instance on the evaluator \param eh Entity handle whose type and #vertices are
 * queried
 */
 ErrorCode set_eval_set( const EntityHandle eh );
 
 /** \brief Get the eval set for a given type entity */
 inline EvalSet get_eval_set( EntityType tp )
 {
 return evalSets[tp];
 }
 
 /** \brief Set entity handle & cache connectivty & vertex positions */
 inline ErrorCode set_ent_handle( EntityHandle ent );
 
 /** \brief Get entity handle for this ElemEval */
 inline EntityHandle get_ent_handle() const
 {
 return entHandle;
 }
 
 /* \brief Get vertex positions cached on this EE
 */
 inline double* get_vert_pos()
 {
 return vertPos[0].array();
 }
 
 /* \brief Get the vertex handles cached here */
 inline const EntityHandle* get_vert_handles() const
 {
 return vertHandles;
 }
 
 /* \brief Get the number of vertices for the cached entity */
 inline int get_num_verts() const
 {
 return numVerts;
 }
 
 /* \brief Get the tag handle cached on this entity */
 inline Tag get_tag_handle() const
 {
 return tagHandle;
 };
 
 /* \brief Set the tag handle to cache on this evaluator
 * To designate that vertex coordinates are the desired tag, pass in a tag handle of 0
 * and a tag dimension of 0.
 * \param tag Tag handle to cache, or 0 to cache vertex positions
 * \param tagged_ent_dim Dimension of entities tagged with this tag
 */
 inline ErrorCode set_tag_handle( Tag tag, int tagged_ent_dim = -1 );
 
 /* \brief Set the name of the tag to cache on this evaluator
 * To designate that vertex coordinates are the desired tag, pass in "COORDS" as the name
 * and a tag dimension of 0.
 * \param tag_name Tag handle to cache, or 0 to cache vertex positions
 * \param tagged_ent_dim Dimension of entities tagged with this tag
 */
 inline ErrorCode set_tag( const char* tag_name, int tagged_ent_dim = -1 );
 
 /* \brief Get the dimension of the entities on which tag is set */
 inline int get_tagged_ent_dim() const
 {
 return taggedEntDim;
 };
 
 /* \brief Set the dimension of entities having the tag */
 inline ErrorCode set_tagged_ent_dim( int dim );
 
 /* \brief Get work space, sometimes this is useful for evaluating data you don't want to set as
 * tags on entities Can't be const because most of the functions (evaluate, integrate, etc.)
 * take non-const work space *'s
 */
 inline double* get_work_space()
 {
 return workSpace;
 }
 
 /* \brief MOAB interface cached on this evaluator */
 inline Interface* get_moab()
 {
 return mbImpl;
 }
 
 private:
 /** \brief Interface */
 Interface* mbImpl;
 
 /** \brief Entity handle being evaluated */
 EntityHandle entHandle;
 
 /** \brief Entity type */
 EntityType entType;
 
 /** \brief Entity dimension */
 int entDim;
 
 /** \brief Number of vertices cached here */
 int numVerts;
 
 /** \brief Cached copy of vertex handle ptr */
 const EntityHandle* vertHandles;
 
 /** \brief Cached copy of vertex positions */
 CartVect vertPos[CN::MAX_NODES_PER_ELEMENT];
 
 /** \brief Tag being evaluated */
 Tag tagHandle;
 
 /** \brief Whether tag is coordinates or something else */
 bool tagCoords;
 
 /** \brief Number of values in this tag */
 int numTuples;
 
 /** \brief Dimension of entities from which to grab tag */
 int taggedEntDim;
 
 /** \brief Tag space */
 std::vector< unsigned char > tagSpace;
 
 /** \brief Evaluation methods for elements of various topologies */
 EvalSet evalSets[MBMAXTYPE];
 
 /** \brief Work space for element-specific data */
 double* workSpace;
 
}; // class ElemEvaluator
 
inline ElemEvaluator::ElemEvaluator( Interface* impl, EntityHandle ent, Tag tag, int tagged_ent_dim )
 : mbImpl( impl ), entHandle( 0 ), entType( MBMAXTYPE ), entDim( -1 ), numVerts( 0 ), vertHandles( NULL ),
 tagHandle( 0 ), tagCoords( false ), numTuples( 0 ), taggedEntDim( 0 ), workSpace( NULL )
{
 if( ent ) set_ent_handle( ent );
 if( tag ) set_tag_handle( tag, tagged_ent_dim );
}
 
inline ElemEvaluator::~ElemEvaluator()
{
 if( workSpace ) delete[] workSpace;
}
 
inline ErrorCode ElemEvaluator::set_ent_handle( EntityHandle ent )
{
 entHandle = ent;
 if( workSpace )
 {
 delete[] workSpace;
 workSpace = NULL;
 }
 
 entType = mbImpl->type_from_handle( ent );
 entDim = mbImpl->dimension_from_handle( ent );
 
 std::vector< EntityHandle > dum_vec;
 ErrorCode rval = mbImpl->get_connectivity( ent, vertHandles, numVerts, false, &dum_vec );
 if( MB_SUCCESS != rval ) return rval;
 
 if( !evalSets[entType].evalFcn ) EvalSet::get_eval_set( entType, numVerts, evalSets[entType] );
 
 rval = mbImpl->get_coords( vertHandles, numVerts, vertPos[0].array() );
 if( MB_SUCCESS != rval ) return rval;
 
 if( tagHandle )
 {
 rval = set_tag_handle( tagHandle );
 if( MB_SUCCESS != rval ) return rval;
 }
 if( evalSets[entType].initFcn ) return ( *evalSets[entType].initFcn )( vertPos[0].array(), numVerts, workSpace );
 return MB_SUCCESS;
}
 
inline ErrorCode ElemEvaluator::set_tag_handle( Tag tag, int tagged_ent_dim )
{
 ErrorCode rval = MB_SUCCESS;
 if( !tag && !tagged_ent_dim )
 {
 tagCoords = true;
 numTuples = 3;
 taggedEntDim = 0;
 tagHandle = 0;
 return rval;
 }
 else if( tagHandle != tag )
 {
 tagHandle = tag;
 rval = mbImpl->tag_get_length( tagHandle, numTuples );
 if( MB_SUCCESS != rval ) return rval;
 int sz;
 rval = mbImpl->tag_get_bytes( tag, sz );
 if( MB_SUCCESS != rval ) return rval;
 tagSpace.resize( CN::MAX_NODES_PER_ELEMENT * sz );
 tagCoords = false;
 }
 
 taggedEntDim = ( -1 == tagged_ent_dim ? 0 : tagged_ent_dim );
 
 if( entHandle )
 {
 if( 0 == taggedEntDim )
 {
 rval = mbImpl->tag_get_data( tagHandle, vertHandles, numVerts, &tagSpace[0] );
 if( MB_SUCCESS != rval ) return rval;
 }
 else if( taggedEntDim == entDim )
 {
 rval = mbImpl->tag_get_data( tagHandle, &entHandle, 1, &tagSpace[0] );
 if( MB_SUCCESS != rval ) return rval;
 }
 }
 
 return rval;
}
 
inline ErrorCode ElemEvaluator::set_tag( const char* tag_name, int tagged_ent_dim )
{
 ErrorCode rval = MB_SUCCESS;
 if( !tag_name || strlen( tag_name ) == 0 ) return MB_FAILURE;
 Tag tag;
 if( !strcmp( tag_name, "COORDS" ) )
 {
 tagCoords = true;
 taggedEntDim = 0;
 numTuples = 3;
 tagHandle = 0;
 // can return here, because vertex coords already cached when entity handle set
 return rval;
 }
 else
 {
 rval = mbImpl->tag_get_handle( tag_name, tag );
 if( MB_SUCCESS != rval ) return rval;
 
 if( tagHandle != tag )
 {
 tagHandle = tag;
 rval = mbImpl->tag_get_length( tagHandle, numTuples );
 if( MB_SUCCESS != rval ) return rval;
 int sz;
 rval = mbImpl->tag_get_bytes( tag, sz );
 if( MB_SUCCESS != rval ) return rval;
 tagSpace.reserve( CN::MAX_NODES_PER_ELEMENT * sz );
 tagCoords = false;
 }
 
 taggedEntDim = ( -1 == tagged_ent_dim ? entDim : tagged_ent_dim );
 }
 
 if( entHandle )
 {
 if( 0 == taggedEntDim )
 {
 rval = mbImpl->tag_get_data( tagHandle, vertHandles, numVerts, &tagSpace[0] );
 if( MB_SUCCESS != rval ) return rval;
 }
 else if( taggedEntDim == entDim )
 {
 rval = mbImpl->tag_get_data( tagHandle, &entHandle, 1, &tagSpace[0] );
 if( MB_SUCCESS != rval ) return rval;
 }
 }
 
 return rval;
}
 
inline ErrorCode ElemEvaluator::set_eval_set( EntityType tp, const EvalSet& eval_set )
{
 evalSets[tp] = eval_set;
 if( entHandle && evalSets[entType].initFcn )
 {
 ErrorCode rval = ( *evalSets[entType].initFcn )( vertPos[0].array(), numVerts, workSpace );
 if( MB_SUCCESS != rval ) return rval;
 }
 return MB_SUCCESS;
}
 
inline ErrorCode ElemEvaluator::eval( const double* params, double* result, int num_tuples ) const
{
 assert( entHandle && MBMAXTYPE != entType );
 return ( *evalSets[entType].evalFcn )(
 params, ( tagCoords ? (const double*)vertPos[0].array() : (const double*)&tagSpace[0] ), entDim,
 ( -1 == num_tuples ? numTuples : num_tuples ), workSpace, result );
}
 
inline ErrorCode ElemEvaluator::reverse_eval( const double* posn,
 const double iter_tol,
 const double inside_tol,
 double* params,
 int* ins ) const
{
 assert( entHandle && MBMAXTYPE != entType );
 return ( *evalSets[entType].reverseEvalFcn )( evalSets[entType].evalFcn, evalSets[entType].jacobianFcn,
 evalSets[entType].insideFcn, posn, vertPos[0].array(), numVerts,
 entDim, iter_tol, inside_tol, workSpace, params, ins );
}
 
/** \brief Evaluate the normal of the cached entity at a given facet */
inline ErrorCode ElemEvaluator::get_normal( const int ientDim, const int facet, double normal[] ) const
{
 assert( entHandle && MBMAXTYPE != entType );
 return ( *evalSets[entType].normalFcn )( ientDim, facet, numVerts, vertPos[0].array(), normal );
}
 
/** \brief Evaluate the jacobian of the cached entity at a given parametric location */
inline ErrorCode ElemEvaluator::jacobian( const double* params, double* result ) const
{
 assert( entHandle && MBMAXTYPE != entType );
 return ( *evalSets[entType].jacobianFcn )( params, vertPos[0].array(), numVerts, entDim, workSpace, result );
}
 
/** \brief Integrate the cached tag over the cached entity */
inline ErrorCode ElemEvaluator::integrate( double* result ) const
{
 assert( entHandle && MBMAXTYPE != entType && ( tagCoords || tagHandle ) );
 ErrorCode rval = MB_SUCCESS;
 if( !tagCoords )
 {
 if( 0 == taggedEntDim )
 rval = mbImpl->tag_get_data( tagHandle, vertHandles, numVerts, (void*)&tagSpace[0] );
 else
 rval = mbImpl->tag_get_data( tagHandle, &entHandle, 1, (void*)&tagSpace[0] );
 if( MB_SUCCESS != rval ) return rval;
 }
 return ( *evalSets[entType].integrateFcn )( ( tagCoords ? vertPos[0].array() : (const double*)&tagSpace[0] ),
 vertPos[0].array(), numVerts, entDim, numTuples, workSpace, result );
}
 
inline ErrorCode ElemEvaluator::find_containing_entity( EntityHandle ent_set,
 const double* point,
 const double iter_tol,
 const double inside_tol,
 EntityHandle& containing_ent,
 double* params,
 unsigned int* num_evals )
{
 assert( mbImpl->type_from_handle( ent_set ) == MBENTITYSET );
 Range entities;
 ErrorCode rval = mbImpl->get_entities_by_handle( ent_set, entities );
 if( MB_SUCCESS != rval )
 return rval;
 else
 return find_containing_entity( entities, point, iter_tol, inside_tol, containing_ent, params, num_evals );
}
 
inline int ElemEvaluator::inside( const double* params, const double tol ) const
{
 return ( *evalSets[entType].insideFcn )( params, entDim, tol );
}
 
inline ErrorCode ElemEvaluator::set_eval_set( const EntityHandle eh )
{
 EvalSet eset;
 ErrorCode rval = EvalSet::get_eval_set( mbImpl, eh, evalSets[mbImpl->type_from_handle( eh )] );
 return rval;
}
 
} // namespace moab
 
#endif /*MOAB_ELEM_EVALUATOR_HPP*/