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);
        
    class EvalSet
    {
  public:
      EvalFcn evalFcn;
        
      ReverseEvalFcn reverseEvalFcn;
        
      JacobianFcn jacobianFcn;
        
      IntegrateFcn integrateFcn;

      InitFcn initFcn;

      InsideFcn insideFcn;

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

      EvalSet(EvalFcn eval, ReverseEvalFcn rev, JacobianFcn jacob, IntegrateFcn integ, InitFcn initf, InsideFcn insidef)
              : evalFcn(eval), reverseEvalFcn(rev), jacobianFcn(jacob), integrateFcn(integ), initFcn(initf), insideFcn(insidef)
          {}

      static ErrorCode get_eval_set(Interface *mb, EntityHandle eh, EvalSet &eval_set);
      
      static ErrorCode get_eval_set(EntityType tp, unsigned int num_vertices, EvalSet &eval_set);
      
      EvalSet &operator=(const EvalSet &eval) {
        evalFcn = eval.evalFcn;
        reverseEvalFcn = eval.reverseEvalFcn;
        jacobianFcn = eval.jacobianFcn;
        integrateFcn = eval.integrateFcn;
        initFcn = eval.initFcn;
        insideFcn = eval.insideFcn;
        return *this;
      }

      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 int inside_function(const double *params, const int ndims, const double tol);
    };

    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);
    }
      
    class ElemEvaluator {
  public:
      ElemEvaluator(Interface *impl, EntityHandle ent = 0, Tag tag = 0, int tagged_ent_dim = -1);

      ErrorCode eval(const double *params, double *result, int num_tuples = -1) const;
        
      ErrorCode reverse_eval(const double *posn, double iter_tol, double inside_tol, double *params, 
                             int *is_inside = NULL) const;
        
      ErrorCode jacobian(const double *params, double *result) const;
        
      ErrorCode integrate(double *result) const;
      
      int inside(const double *params, const double tol) const;

      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);
      
      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);
      
      ErrorCode set_eval_set(EntityType tp, const EvalSet &eval_set);
      
      ErrorCode set_eval_set(const EntityHandle eh);
      
      inline EvalSet get_eval_set(EntityType tp) {return evalSets[tp];}

      inline ErrorCode set_ent_handle(EntityHandle ent);

      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:

      Interface *mbImpl;
      
      EntityHandle entHandle;

      EntityType entType;

      int entDim;
            
      int numVerts;

      const EntityHandle *vertHandles;
      
      CartVect vertPos[CN::MAX_NODES_PER_ELEMENT];
      
      Tag tagHandle;

      bool tagCoords;
      
      int numTuples;

      int taggedEntDim;

      std::vector<unsigned char> tagSpace;
      
      EvalSet evalSets[MBMAXTYPE];

      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 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.reserve(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);
    }
        
    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);
    }
        
    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*/