Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
ElemEvaluator.hpp
Go to the documentation of this file.
00001 #ifndef ELEM_EVALUATOR_HPP
00002 #define ELEM_EVALUATOR_HPP
00003 
00004 #include "moab/Interface.hpp"
00005 #include "moab/CartVect.hpp"
00006 #include "moab/Matrix3.hpp"
00007 #include "moab/CN.hpp"
00008 
00009 #include <vector>
00010 
00011 namespace moab
00012 {
00013 
00014 typedef ErrorCode ( *EvalFcn )( const double* params,
00015                                 const double* field,
00016                                 const int ndim,
00017                                 const int num_tuples,
00018                                 double* work,
00019                                 double* result );
00020 
00021 typedef ErrorCode ( *JacobianFcn )( const double* params,
00022                                     const double* verts,
00023                                     const int nverts,
00024                                     const int ndim,
00025                                     double* work,
00026                                     double* result );
00027 
00028 typedef ErrorCode ( *IntegrateFcn )( const double* field,
00029                                      const double* verts,
00030                                      const int nverts,
00031                                      const int ndim,
00032                                      const int num_tuples,
00033                                      double* work,
00034                                      double* result );
00035 
00036 typedef ErrorCode ( *InitFcn )( const double* verts, const int nverts, double*& work );
00037 
00038 typedef int ( *InsideFcn )( const double* verts, const int ndims, const double tol );
00039 
00040 typedef ErrorCode ( *ReverseEvalFcn )( EvalFcn eval,
00041                                        JacobianFcn jacob,
00042                                        InsideFcn ins,
00043                                        const double* posn,
00044                                        const double* verts,
00045                                        const int nverts,
00046                                        const int ndim,
00047                                        const double iter_tol,
00048                                        const double inside_tol,
00049                                        double* work,
00050                                        double* params,
00051                                        int* is_inside );
00052 
00053 typedef ErrorCode (
00054     *NormalFcn )( const int ientDim, const int facet, const int nverts, const double* verts, double normal[3] );
00055 
00056 class EvalSet
00057 {
00058   public:
00059     /** \brief Forward-evaluation of field at parametric coordinates */
00060     EvalFcn evalFcn;
00061 
00062     /** \brief Reverse-evaluation of parametric coordinates at physical space position */
00063     ReverseEvalFcn reverseEvalFcn;
00064 
00065     /** \brief Evaluate the normal at a local facet (edge/face for 2D/3D) */
00066     NormalFcn normalFcn;
00067 
00068     /** \brief Evaluate the jacobian at a specified parametric position */
00069     JacobianFcn jacobianFcn;
00070 
00071     /** \brief Forward-evaluation of field at parametric coordinates */
00072     IntegrateFcn integrateFcn;
00073 
00074     /** \brief Initialization function for an element */
00075     InitFcn initFcn;
00076 
00077     /** \brief Function that returns whether or not the parameters are inside the natural space of
00078      * the element */
00079     InsideFcn insideFcn;
00080 
00081     /** \brief Bare constructor */
00082     EvalSet()
00083         : evalFcn( NULL ), reverseEvalFcn( NULL ), normalFcn( NULL ), jacobianFcn( NULL ), integrateFcn( NULL ),
00084           initFcn( NULL ), insideFcn( NULL )
00085     {
00086     }
00087 
00088     /** \brief Constructor */
00089     EvalSet( EvalFcn eval,
00090              ReverseEvalFcn rev,
00091              NormalFcn normal,
00092              JacobianFcn jacob,
00093              IntegrateFcn integ,
00094              InitFcn initf,
00095              InsideFcn insidef )
00096         : evalFcn( eval ), reverseEvalFcn( rev ), normalFcn( normal ), jacobianFcn( jacob ), integrateFcn( integ ),
00097           initFcn( initf ), insideFcn( insidef )
00098     {
00099     }
00100 
00101     /** \brief Copy constructor */
00102     EvalSet( EvalSet const& eval )
00103         : evalFcn( eval.evalFcn ), reverseEvalFcn( eval.reverseEvalFcn ), normalFcn( eval.normalFcn ),
00104           jacobianFcn( eval.jacobianFcn ), integrateFcn( eval.integrateFcn ), initFcn( eval.initFcn ),
00105           insideFcn( eval.insideFcn )
00106     {
00107     }
00108 
00109     /** \brief Given an entity handle, get an appropriate eval set, based on type & #vertices */
00110     static ErrorCode get_eval_set( Interface* mb, EntityHandle eh, EvalSet& eval_set );
00111 
00112     /** \brief Given type & #vertices, get an appropriate eval set */
00113     static ErrorCode get_eval_set( EntityType tp, unsigned int num_vertices, EvalSet& eval_set );
00114 
00115     /** \brief Operator= */
00116     EvalSet& operator=( const EvalSet& eval )
00117     {
00118         evalFcn        = eval.evalFcn;
00119         reverseEvalFcn = eval.reverseEvalFcn;
00120         normalFcn      = eval.normalFcn;
00121         jacobianFcn    = eval.jacobianFcn;
00122         integrateFcn   = eval.integrateFcn;
00123         initFcn        = eval.initFcn;
00124         insideFcn      = eval.insideFcn;
00125         return *this;
00126     }
00127 
00128     /** \brief Common function to do reverse evaluation based on evaluation and jacobian functions
00129      */
00130     static ErrorCode evaluate_reverse( EvalFcn eval,
00131                                        JacobianFcn jacob,
00132                                        InsideFcn inside_f,
00133                                        const double* posn,
00134                                        const double* verts,
00135                                        const int nverts,
00136                                        const int ndim,
00137                                        const double iter_tol,
00138                                        const double inside_tol,
00139                                        double* work,
00140                                        double* params,
00141                                        int* inside );
00142     /** \brief Common function that returns true if params is in [-1,1]^ndims */
00143     static int inside_function( const double* params, const int ndims, const double tol );
00144 };
00145 
00146 /** \brief Given an entity handle, get an appropriate eval set, based on type & #vertices */
00147 inline ErrorCode EvalSet::get_eval_set( Interface* mb, EntityHandle eh, EvalSet& eval_set )
00148 {
00149     int nv;
00150     EntityType tp = mb->type_from_handle( eh );
00151     const EntityHandle* connect;
00152     std::vector< EntityHandle > dum_vec;
00153     ErrorCode rval = mb->get_connectivity( eh, connect, nv, false, &dum_vec );
00154     if( MB_SUCCESS != rval ) return rval;
00155 
00156     return get_eval_set( tp, nv, eval_set );
00157 }
00158 
00159 /**\brief Class facilitating local discretization-related functions
00160  * \class ElemEvaluator
00161  * This class implements discretization-related functionality operating
00162  * on data in MOAB.  A member of this class caches certain data about the element
00163  * it's currently operating on, but is not meant to be instantiated one-per-element,
00164  * but rather one-per-search (or other operation on a collection of elements).
00165  *
00166  * Actual discretization functionality is accessed through function pointers,
00167  * allowing applications to specialize the implementation of specific functions
00168  * while still using this class.
00169  *
00170  * This class depends on MOAB functionality for gathering entity-based data; the functions
00171  * it calls through function pointers depend only on POD (plain old data, or intrinsic data types).
00172  * This allows the use of other packages for serving these functions, without having to modify
00173  * them to get data through MOAB.  This should also promote efficiency, since in many cases they
00174  * will be able to read data from its native storage locations.
00175  */
00176 
00177 class ElemEvaluator
00178 {
00179   public:
00180     /** \brief Constructor
00181      * \param impl MOAB instance
00182      * \param ent Entity handle to cache on the evaluator; if non-zero, calls set_ent_handle, which
00183      * does some other stuff. \param tag Tag to cache on the evaluator; if non-zero, calls
00184      * set_tag_handle, which does some other stuff too. \param tagged_ent_dim Dimension of entities
00185      * to be tagged to cache on the evaluator
00186      */
00187     ElemEvaluator( Interface* impl, EntityHandle ent = 0, Tag tag = 0, int tagged_ent_dim = -1 );
00188     ~ElemEvaluator();
00189 
00190     /** \brief Evaluate cached tag at a given parametric location within the cached entity
00191      * If evaluating coordinates, call set_tag(0, 0), which indicates coords instead of a tag.
00192      * \param params Parameters at which to evaluate field
00193      * \param result Result of evaluation
00194      * \param num_tuples Size of evaluated field, in number of values
00195      */
00196     ErrorCode eval( const double* params, double* result, int num_tuples = -1 ) const;
00197 
00198     /** \brief Reverse-evaluate the cached entity at a given physical position
00199      * \param posn Position at which to evaluate parameters
00200      * \param iter_tol Tolerance of reverse evaluation non-linear iteration, usually 10^-10 or so
00201      * \param inside_tol Tolerance of is_inside evaluation, usually 10^-6 or so
00202      * \param params Result of evaluation
00203      * \param is_inside If non-NULL, returns true of resulting parameters place the point inside the
00204      * element (in most cases, within [-1]*(dim)
00205      */
00206     ErrorCode reverse_eval( const double* posn,
00207                             double iter_tol,
00208                             double inside_tol,
00209                             double* params,
00210                             int* is_inside = NULL ) const;
00211 
00212     /**
00213      * \brief Evaluate the normal to a facet of an entity
00214      * \param ientDim Dimension of the facet. Should be (d-1) for d-dimensional entities
00215      * \param facet Local id of the facet w.r.t the entity
00216      * \param normal Returns the normal.
00217      */
00218     ErrorCode get_normal( const int ientDim, const int facet, double normal[3] ) const;
00219 
00220     /** \brief Evaluate the jacobian of the cached entity at a given parametric location
00221      * \param params Parameters at which to evaluate jacobian
00222      * \param result Result of evaluation, in the form of a 3x3 matrix, stored in column-major order
00223      */
00224     ErrorCode jacobian( const double* params, double* result ) const;
00225 
00226     /** \brief Integrate the cached tag over the cached entity
00227      * \param result Result of the integration
00228      */
00229     ErrorCode integrate( double* result ) const;
00230 
00231     /** \brief Return whether a physical position is inside the cached entity to within a tolerance
00232      * \param params Parameters at which to query the element
00233      * \param tol Tolerance, usually 10^-6 or so
00234      */
00235     int inside( const double* params, const double tol ) const;
00236 
00237     /** \brief Given a list of entities, return the entity the point is in, or none
00238      * This function reverse-evaluates the entities, returning the first entity containing the
00239      * point. If no entity contains the point, containing_ent is returned as 0 and params are
00240      * unchanged. This function returns something other than MB_SUCCESS only when queries go wrong
00241      * for some reason. num_evals, if non-NULL, is always incremented for each call to reverse_eval.
00242      * This function calls set_ent_handle for each entity before calling reverse_eval, so the
00243      * ElemEvaluator object is changed. \param entities Entities tested \param point Point tested,
00244      * must have 3 dimensions, even for edge and face entities \param iter_tol Tolerance for
00245      * non-linear reverse evaluation \param inside_tol Tolerance for is_inside test \param
00246      * containing_ent Entity containing the point, returned 0 if no entity \param params Parameters
00247      * of point in containing entity, unchanged if no containing entity \param num_evals If
00248      * non-NULL, incremented each time reverse_eval is called \return Returns non-success only if
00249      * evaulation failed for some reason (point not in element is NOT a reason for failure)
00250      */
00251     ErrorCode find_containing_entity( Range& entities,
00252                                       const double* point,
00253                                       const double iter_tol,
00254                                       const double inside_tol,
00255                                       EntityHandle& containing_ent,
00256                                       double* params,
00257                                       unsigned int* num_evals = NULL );
00258 
00259     /** \brief Given an entity set, return the contained entity the point is in, or none
00260      * This function reverse-evaluates the entities, returning the first entity containing the
00261      * point. If no entity contains the point, containing_ent is returned as 0 and params are
00262      * unchanged. This function returns something other than MB_SUCCESS only when queries go wrong
00263      * for some reason. num_evals, if non-NULL, is always incremented for each call to reverse_eval.
00264      * This function calls set_ent_handle for each entity before calling reverse_eval, so the
00265      * ElemEvaluator object is changed. \param ent_set Entity set containing the entities to be
00266      * tested \param point Point tested, must have 3 dimensions, even for edge and face entities
00267      * \param iter_tol Tolerance for non-linear reverse evaluation
00268      * \param inside_tol Tolerance for is_inside test
00269      * \param containing_ent Entity containing the point, returned 0 if no entity
00270      * \param params Parameters of point in containing entity, unchanged if no containing entity
00271      * \param num_evals If non-NULL, incremented each time reverse_eval is called
00272      * \return Returns non-success only if evaulation failed for some reason (point not in element
00273      * is NOT a reason for failure)
00274      */
00275     ErrorCode find_containing_entity( EntityHandle ent_set,
00276                                       const double* point,
00277                                       const double iter_tol,
00278                                       const double inside_tol,
00279                                       EntityHandle& containing_ent,
00280                                       double* params,
00281                                       unsigned int* num_evals = NULL );
00282 
00283     /** \brief Set the eval set for a given type entity
00284      * \param tp Entity type for which to set the eval set
00285      * \param eval_set Eval set object to use
00286      */
00287     ErrorCode set_eval_set( EntityType tp, const EvalSet& eval_set );
00288 
00289     /** \brief Set the eval set using a given entity handle
00290      * This function queries the entity type and number of vertices on the entity to decide which
00291      * type of shape function to use. NOTE: this function should only be called after setting a
00292      * valid MOAB instance on the evaluator \param eh Entity handle whose type and #vertices are
00293      * queried
00294      */
00295     ErrorCode set_eval_set( const EntityHandle eh );
00296 
00297     /** \brief Get the eval set for a given type entity */
00298     inline EvalSet get_eval_set( EntityType tp )
00299     {
00300         return evalSets[tp];
00301     }
00302 
00303     /** \brief Set entity handle & cache connectivty & vertex positions */
00304     inline ErrorCode set_ent_handle( EntityHandle ent );
00305 
00306     /** \brief Get entity handle for this ElemEval */
00307     inline EntityHandle get_ent_handle() const
00308     {
00309         return entHandle;
00310     }
00311 
00312     /* \brief Get vertex positions cached on this EE
00313      */
00314     inline double* get_vert_pos()
00315     {
00316         return vertPos[0].array();
00317     }
00318 
00319     /* \brief Get the vertex handles cached here */
00320     inline const EntityHandle* get_vert_handles() const
00321     {
00322         return vertHandles;
00323     }
00324 
00325     /* \brief Get the number of vertices for the cached entity */
00326     inline int get_num_verts() const
00327     {
00328         return numVerts;
00329     }
00330 
00331     /* \brief Get the tag handle cached on this entity */
00332     inline Tag get_tag_handle() const
00333     {
00334         return tagHandle;
00335     };
00336 
00337     /* \brief Set the tag handle to cache on this evaluator
00338      * To designate that vertex coordinates are the desired tag, pass in a tag handle of 0
00339      * and a tag dimension of 0.
00340      * \param tag Tag handle to cache, or 0 to cache vertex positions
00341      * \param tagged_ent_dim Dimension of entities tagged with this tag
00342      */
00343     inline ErrorCode set_tag_handle( Tag tag, int tagged_ent_dim = -1 );
00344 
00345     /* \brief Set the name of the tag to cache on this evaluator
00346      * To designate that vertex coordinates are the desired tag, pass in "COORDS" as the name
00347      * and a tag dimension of 0.
00348      * \param tag_name Tag handle to cache, or 0 to cache vertex positions
00349      * \param tagged_ent_dim Dimension of entities tagged with this tag
00350      */
00351     inline ErrorCode set_tag( const char* tag_name, int tagged_ent_dim = -1 );
00352 
00353     /* \brief Get the dimension of the entities on which tag is set */
00354     inline int get_tagged_ent_dim() const
00355     {
00356         return taggedEntDim;
00357     };
00358 
00359     /* \brief Set the dimension of entities having the tag */
00360     inline ErrorCode set_tagged_ent_dim( int dim );
00361 
00362     /* \brief Get work space, sometimes this is useful for evaluating data you don't want to set as
00363      * tags on entities Can't be const because most of the functions (evaluate, integrate, etc.)
00364      * take non-const work space *'s
00365      */
00366     inline double* get_work_space()
00367     {
00368         return workSpace;
00369     }
00370 
00371     /* \brief MOAB interface cached on this evaluator */
00372     inline Interface* get_moab()
00373     {
00374         return mbImpl;
00375     }
00376 
00377   private:
00378     /** \brief Interface */
00379     Interface* mbImpl;
00380 
00381     /** \brief Entity handle being evaluated */
00382     EntityHandle entHandle;
00383 
00384     /** \brief Entity type */
00385     EntityType entType;
00386 
00387     /** \brief Entity dimension */
00388     int entDim;
00389 
00390     /** \brief Number of vertices cached here */
00391     int numVerts;
00392 
00393     /** \brief Cached copy of vertex handle ptr */
00394     const EntityHandle* vertHandles;
00395 
00396     /** \brief Cached copy of vertex positions */
00397     CartVect vertPos[CN::MAX_NODES_PER_ELEMENT];
00398 
00399     /** \brief Tag being evaluated */
00400     Tag tagHandle;
00401 
00402     /** \brief Whether tag is coordinates or something else */
00403     bool tagCoords;
00404 
00405     /** \brief Number of values in this tag */
00406     int numTuples;
00407 
00408     /** \brief Dimension of entities from which to grab tag */
00409     int taggedEntDim;
00410 
00411     /** \brief Tag space */
00412     std::vector< unsigned char > tagSpace;
00413 
00414     /** \brief Evaluation methods for elements of various topologies */
00415     EvalSet evalSets[MBMAXTYPE];
00416 
00417     /** \brief Work space for element-specific data */
00418     double* workSpace;
00419 
00420 };  // class ElemEvaluator
00421 
00422 inline ElemEvaluator::ElemEvaluator( Interface* impl, EntityHandle ent, Tag tag, int tagged_ent_dim )
00423     : mbImpl( impl ), entHandle( 0 ), entType( MBMAXTYPE ), entDim( -1 ), numVerts( 0 ), vertHandles( NULL ),
00424       tagHandle( 0 ), tagCoords( false ), numTuples( 0 ), taggedEntDim( 0 ), workSpace( NULL )
00425 {
00426     if( ent ) set_ent_handle( ent );
00427     if( tag ) set_tag_handle( tag, tagged_ent_dim );
00428 }
00429 
00430 inline ElemEvaluator::~ElemEvaluator()
00431 {
00432     if( workSpace ) delete[] workSpace;
00433 }
00434 
00435 inline ErrorCode ElemEvaluator::set_ent_handle( EntityHandle ent )
00436 {
00437     entHandle = ent;
00438     if( workSpace )
00439     {
00440         delete[] workSpace;
00441         workSpace = NULL;
00442     }
00443 
00444     entType = mbImpl->type_from_handle( ent );
00445     entDim  = mbImpl->dimension_from_handle( ent );
00446 
00447     std::vector< EntityHandle > dum_vec;
00448     ErrorCode rval = mbImpl->get_connectivity( ent, vertHandles, numVerts, false, &dum_vec );
00449     if( MB_SUCCESS != rval ) return rval;
00450 
00451     if( !evalSets[entType].evalFcn ) EvalSet::get_eval_set( entType, numVerts, evalSets[entType] );
00452 
00453     rval = mbImpl->get_coords( vertHandles, numVerts, vertPos[0].array() );
00454     if( MB_SUCCESS != rval ) return rval;
00455 
00456     if( tagHandle )
00457     {
00458         rval = set_tag_handle( tagHandle );
00459         if( MB_SUCCESS != rval ) return rval;
00460     }
00461     if( evalSets[entType].initFcn ) return ( *evalSets[entType].initFcn )( vertPos[0].array(), numVerts, workSpace );
00462     return MB_SUCCESS;
00463 }
00464 
00465 inline ErrorCode ElemEvaluator::set_tag_handle( Tag tag, int tagged_ent_dim )
00466 {
00467     ErrorCode rval = MB_SUCCESS;
00468     if( !tag && !tagged_ent_dim )
00469     {
00470         tagCoords    = true;
00471         numTuples    = 3;
00472         taggedEntDim = 0;
00473         tagHandle    = 0;
00474         return rval;
00475     }
00476     else if( tagHandle != tag )
00477     {
00478         tagHandle = tag;
00479         rval      = mbImpl->tag_get_length( tagHandle, numTuples );
00480         if( MB_SUCCESS != rval ) return rval;
00481         int sz;
00482         rval = mbImpl->tag_get_bytes( tag, sz );
00483         if( MB_SUCCESS != rval ) return rval;
00484         tagSpace.resize( CN::MAX_NODES_PER_ELEMENT * sz );
00485         tagCoords = false;
00486     }
00487 
00488     taggedEntDim = ( -1 == tagged_ent_dim ? 0 : tagged_ent_dim );
00489 
00490     if( entHandle )
00491     {
00492         if( 0 == taggedEntDim )
00493         {
00494             rval = mbImpl->tag_get_data( tagHandle, vertHandles, numVerts, &tagSpace[0] );
00495             if( MB_SUCCESS != rval ) return rval;
00496         }
00497         else if( taggedEntDim == entDim )
00498         {
00499             rval = mbImpl->tag_get_data( tagHandle, &entHandle, 1, &tagSpace[0] );
00500             if( MB_SUCCESS != rval ) return rval;
00501         }
00502     }
00503 
00504     return rval;
00505 }
00506 
00507 inline ErrorCode ElemEvaluator::set_tag( const char* tag_name, int tagged_ent_dim )
00508 {
00509     ErrorCode rval = MB_SUCCESS;
00510     if( !tag_name || strlen( tag_name ) == 0 ) return MB_FAILURE;
00511     Tag tag;
00512     if( !strcmp( tag_name, "COORDS" ) )
00513     {
00514         tagCoords    = true;
00515         taggedEntDim = 0;
00516         numTuples    = 3;
00517         tagHandle    = 0;
00518         // can return here, because vertex coords already cached when entity handle set
00519         return rval;
00520     }
00521     else
00522     {
00523         rval = mbImpl->tag_get_handle( tag_name, tag );
00524         if( MB_SUCCESS != rval ) return rval;
00525 
00526         if( tagHandle != tag )
00527         {
00528             tagHandle = tag;
00529             rval      = mbImpl->tag_get_length( tagHandle, numTuples );
00530             if( MB_SUCCESS != rval ) return rval;
00531             int sz;
00532             rval = mbImpl->tag_get_bytes( tag, sz );
00533             if( MB_SUCCESS != rval ) return rval;
00534             tagSpace.reserve( CN::MAX_NODES_PER_ELEMENT * sz );
00535             tagCoords = false;
00536         }
00537 
00538         taggedEntDim = ( -1 == tagged_ent_dim ? entDim : tagged_ent_dim );
00539     }
00540 
00541     if( entHandle )
00542     {
00543         if( 0 == taggedEntDim )
00544         {
00545             rval = mbImpl->tag_get_data( tagHandle, vertHandles, numVerts, &tagSpace[0] );
00546             if( MB_SUCCESS != rval ) return rval;
00547         }
00548         else if( taggedEntDim == entDim )
00549         {
00550             rval = mbImpl->tag_get_data( tagHandle, &entHandle, 1, &tagSpace[0] );
00551             if( MB_SUCCESS != rval ) return rval;
00552         }
00553     }
00554 
00555     return rval;
00556 }
00557 
00558 inline ErrorCode ElemEvaluator::set_eval_set( EntityType tp, const EvalSet& eval_set )
00559 {
00560     evalSets[tp] = eval_set;
00561     if( entHandle && evalSets[entType].initFcn )
00562     {
00563         ErrorCode rval = ( *evalSets[entType].initFcn )( vertPos[0].array(), numVerts, workSpace );
00564         if( MB_SUCCESS != rval ) return rval;
00565     }
00566     return MB_SUCCESS;
00567 }
00568 
00569 inline ErrorCode ElemEvaluator::eval( const double* params, double* result, int num_tuples ) const
00570 {
00571     assert( entHandle && MBMAXTYPE != entType );
00572     return ( *evalSets[entType].evalFcn )(
00573         params, ( tagCoords ? (const double*)vertPos[0].array() : (const double*)&tagSpace[0] ), entDim,
00574         ( -1 == num_tuples ? numTuples : num_tuples ), workSpace, result );
00575 }
00576 
00577 inline ErrorCode ElemEvaluator::reverse_eval( const double* posn,
00578                                               const double iter_tol,
00579                                               const double inside_tol,
00580                                               double* params,
00581                                               int* ins ) const
00582 {
00583     assert( entHandle && MBMAXTYPE != entType );
00584     return ( *evalSets[entType].reverseEvalFcn )( evalSets[entType].evalFcn, evalSets[entType].jacobianFcn,
00585                                                   evalSets[entType].insideFcn, posn, vertPos[0].array(), numVerts,
00586                                                   entDim, iter_tol, inside_tol, workSpace, params, ins );
00587 }
00588 
00589 /** \brief Evaluate the normal of the cached entity at a given facet */
00590 inline ErrorCode ElemEvaluator::get_normal( const int ientDim, const int facet, double normal[] ) const
00591 {
00592     assert( entHandle && MBMAXTYPE != entType );
00593     return ( *evalSets[entType].normalFcn )( ientDim, facet, numVerts, vertPos[0].array(), normal );
00594 }
00595 
00596 /** \brief Evaluate the jacobian of the cached entity at a given parametric location */
00597 inline ErrorCode ElemEvaluator::jacobian( const double* params, double* result ) const
00598 {
00599     assert( entHandle && MBMAXTYPE != entType );
00600     return ( *evalSets[entType].jacobianFcn )( params, vertPos[0].array(), numVerts, entDim, workSpace, result );
00601 }
00602 
00603 /** \brief Integrate the cached tag over the cached entity */
00604 inline ErrorCode ElemEvaluator::integrate( double* result ) const
00605 {
00606     assert( entHandle && MBMAXTYPE != entType && ( tagCoords || tagHandle ) );
00607     ErrorCode rval = MB_SUCCESS;
00608     if( !tagCoords )
00609     {
00610         if( 0 == taggedEntDim )
00611             rval = mbImpl->tag_get_data( tagHandle, vertHandles, numVerts, (void*)&tagSpace[0] );
00612         else
00613             rval = mbImpl->tag_get_data( tagHandle, &entHandle, 1, (void*)&tagSpace[0] );
00614         if( MB_SUCCESS != rval ) return rval;
00615     }
00616     return ( *evalSets[entType].integrateFcn )( ( tagCoords ? vertPos[0].array() : (const double*)&tagSpace[0] ),
00617                                                 vertPos[0].array(), numVerts, entDim, numTuples, workSpace, result );
00618 }
00619 
00620 inline ErrorCode ElemEvaluator::find_containing_entity( EntityHandle ent_set,
00621                                                         const double* point,
00622                                                         const double iter_tol,
00623                                                         const double inside_tol,
00624                                                         EntityHandle& containing_ent,
00625                                                         double* params,
00626                                                         unsigned int* num_evals )
00627 {
00628     assert( mbImpl->type_from_handle( ent_set ) == MBENTITYSET );
00629     Range entities;
00630     ErrorCode rval = mbImpl->get_entities_by_handle( ent_set, entities );
00631     if( MB_SUCCESS != rval )
00632         return rval;
00633     else
00634         return find_containing_entity( entities, point, iter_tol, inside_tol, containing_ent, params, num_evals );
00635 }
00636 
00637 inline int ElemEvaluator::inside( const double* params, const double tol ) const
00638 {
00639     return ( *evalSets[entType].insideFcn )( params, entDim, tol );
00640 }
00641 
00642 inline ErrorCode ElemEvaluator::set_eval_set( const EntityHandle eh )
00643 {
00644     EvalSet eset;
00645     ErrorCode rval = EvalSet::get_eval_set( mbImpl, eh, evalSets[mbImpl->type_from_handle( eh )] );
00646     return rval;
00647 }
00648 
00649 }  // namespace moab
00650 
00651 #endif /*MOAB_ELEM_EVALUATOR_HPP*/
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines