Cppcheck report - MOAB: src/moab/ElemEvaluator.hpp

#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 );
<--- Class 'ElemEvaluator' has a constructor with 1 argument that is not explicit. [+]
Class 'ElemEvaluator' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
<--- Class 'ElemEvaluator' has a constructor with 1 argument that is not explicit. [+]
Class 'ElemEvaluator' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
<--- Class 'ElemEvaluator' has a constructor with 1 argument that is not explicit. [+]
Class 'ElemEvaluator' has a constructor with 1 argument that is not explicit. Such constructors should in general be explicit for type safety reasons. Using the explicit keyword in the constructor means some mistakes when using the class can be avoided.
    ~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,<--- Casting between unsigned char * and const double * which have an incompatible binary data representation.<--- Casting between unsigned char * and const double * which have an incompatible binary data representation.<--- Casting between unsigned char * and const double * which have an incompatible binary data representation.
        ( -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;<--- Variable 'rval' is assigned a value that is never used.<--- Variable 'rval' is assigned a value that is never used.<--- Variable 'rval' is assigned a value that is never used.
    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] ),<--- Casting between unsigned char * and const double * which have an incompatible binary data representation.<--- Casting between unsigned char * and const double * which have an incompatible binary data representation.<--- Casting between unsigned char * and const double * which have an incompatible binary data representation.
                                                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*/