SpatialLocator.hpp

Go to the documentation of this file.

#ifndef MOAB_SPATIALLOCATOR_HPP
#define MOAB_SPATIALLOCATOR_HPP

#include "moab/Types.hpp"
#include "moab/Tree.hpp"
#include "moab/Range.hpp"
#include "moab/TupleList.hpp"
#include "moab/BoundBox.hpp"
#include "moab/SpatialLocatorTimes.hpp"
#include "moab/CpuTimer.hpp"

#include <string>
#include <vector>
#include <map>
#include <math.h>

namespace moab {

    class Interface;
    class ElemEvaluator;
    class ParallelComm;

    class SpatialLocator
    {
  public:
        /* constructor */
      SpatialLocator(Interface *impl, Range &elems, Tree *tree = NULL, ElemEvaluator *eval = NULL);

        /* destructor */
      virtual ~SpatialLocator();

        /* add elements to be searched */
      ErrorCode add_elems(Range &elems);
      
        /* get bounding box of this locator */
      BoundBox &local_box() {return localBox;}
      
        /* get bounding box of this locator */
      const BoundBox &local_box() const {return localBox;}
      
        /* locate a set of vertices, Range variant */
      ErrorCode locate_points(Range &vertices,
                              EntityHandle *ents, double *params, int *is_inside = NULL,
                              const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
                              const double inside_tol = 1.0e-6);
      
        /* locate a set of points */
      ErrorCode locate_points(const double *pos, int num_points,
                              EntityHandle *ents, double *params, int *is_inside = NULL,
                              const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
                              const double inside_tol = 1.0e-6);
      
        /* locate a set of vertices or entity centroids, storing results on TupleList in this class
         * Locate a set of vertices or entity centroids, storing the detailed results in member 
         * variable (TupleList) locTable (see comments on locTable for structure of that tuple).
         */
      ErrorCode locate_points(Range &ents,
                              const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
                              const double inside_tol = 1.0e-6);
      
        /* locate a set of points, storing results on TupleList in this class
         * Locate a set of points, storing the detailed results in member variable (TupleList) locTable
         * (see comments on locTable for structure of that tuple).
         */
      ErrorCode locate_points(const double *pos, int num_points,
                              const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
                              const double inside_tol = 1.0e-6);

        /* Count the number of located points in locTable
         * Return the number of entries in locTable that have non-zero entity handles, which
         * represents the number of points in targetEnts that were inside one element in sourceEnts
         *
         */
      int local_num_located();

        /* Count the number of located points in parLocTable
         * Return the number of entries in parLocTable that have a non-negative index in on a remote
         * proc in parLocTable, which gives the number of points located in at least one element in a
         * remote proc's sourceEnts.
         */
      int remote_num_located();

#ifdef USE_MPI      
        /* locate a set of vertices or entity centroids, storing results on TupleList in this class
         * Locate a set of vertices or entity centroids, storing the detailed results in member 
         * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for 
         * structure of those tuples).
         */
      ErrorCode par_locate_points(ParallelComm *pc,
                                  Range &vertices,
                                  const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
                                  const double inside_tol = 1.0e-6);
      
        /* locate a set of points, storing results on TupleList in this class
         * Locate a set of points, storing the detailed results in member 
         * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for 
         * structure of those tuples).
         */
      ErrorCode par_locate_points(ParallelComm *pc,
                                  const double *pos, int num_points,
                                  const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
                                  const double inside_tol = 1.0e-6);
#endif

      Interface* moab() { return mbImpl; }

        /* locate a point */
      ErrorCode locate_point(const double *pos, 
                             EntityHandle &ent, double *params, int *is_inside = NULL,
                              const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
                              const double inside_tol = 1.0e-6);

        /* return the tree */
      Tree *get_tree() {return myTree;}

        /* get the locTable
         */
      TupleList &loc_table() {return locTable;}
      
        /* get the locTable
         */
      const TupleList &loc_table() const {return locTable;}
      
        /* get the parLocTable
         */
      TupleList &par_loc_table() {return parLocTable;}
      
        /* get the parLocTable
         */
      const TupleList &par_loc_table() const {return parLocTable;}

        /* get elemEval */
      ElemEvaluator *elem_eval() {return elemEval;}
      
        /* get elemEval */
      const ElemEvaluator *elem_eval() const {return elemEval;}
      
        /* set elemEval */
      void elem_eval(ElemEvaluator *eval) {elemEval = eval; if (myTree) myTree->set_eval(eval);}

      SpatialLocatorTimes &sl_times() {return myTimes;}
      
      const SpatialLocatorTimes &sl_times() const {return myTimes;}
        
  private:

#ifdef USE_MPI
        /* MPI_ReduceAll source mesh bounding boxes to get global source mesh bounding box
         */
      ErrorCode initialize_intermediate_partition(ParallelComm *pc);
      
        /* for a given point in space, compute its ijk location in the intermediate decomposition; tolerance is
         * used only to test proximity to global box extent, not for local box test
         */
      inline ErrorCode get_point_ijk(const CartVect &point, const double abs_iter_tol, int *ijk) const;

        /* given an ijk location in the intermediate partition, return the proc rank for that location 
         */
      inline int proc_from_ijk(const int *ijk) const;

        /* given a point in space, return the proc responsible for that point from the intermediate decomp; no tolerances
         * applied here, so first proc in lexicographic ijk ordering is returned
         */
      inline int proc_from_point(const double *pos, const double abs_iter_tol) const;
      
        /* register my source mesh with intermediate decomposition procs
         */
      ErrorCode register_src_with_intermediate_procs(ParallelComm *pc, double abs_iter_tol, TupleList &TLreg_o);
      
#endif

      void create_tree();
      
        /* MOAB instance */
      Interface* mbImpl;

        /* elements being located */
      Range myElems;

        /* dimension of entities in locator */
      int myDim;
      
        /* tree used for location */
      Tree *myTree;
      
        /* element evaluator */
      ElemEvaluator *elemEval;

        /* whether I created the tree or not (determines whether to delete it or not on destruction) */
      bool iCreatedTree;

        /* \brief local locations table
         * This table stores detailed local location data results from locate_points, that is, location data
         * for points located on the local mesh.  Data is stored
         * in a TupleList, where each tuple consists of (p_i, hs_ul, r[3]_d), where
         *   p_i = (int) proc from which request for this point was made (0 if serial)
         *   hs_ul = (unsigned long) source entity containing the point
         *   r[3]_d = (double) parametric coordinates of the point in hs 
         */
      TupleList locTable;

        /* \brief parallel locations table
         * This table stores information about points located on a local or remote processor.  For 
         * points located on this processor's local mesh, detailed location data is stored in locTable.
         * For points located on remote processors, more communication is required to retrieve specific
         * location data (usually that information isn't useful on this processor).
         *
         * The tuple structure of this TupleList is (p_i, ri_i), where:
         *   p_i = (int) processor rank containing this point
         *   ri_i = (int) index into locTable on remote proc containing this point's location information
         * The indexing of parLocTable corresponds to that of the points/entities passed in.
         */
      TupleList parLocTable;

        /* \brief Local bounding box, duplicated from tree
         */
      BoundBox localBox;

        /* \brief Global bounding box, used in parallel spatial location
         */
      BoundBox globalBox;

        /* \brief Regional delta xyz, used in parallel spatial location
         */
      CartVect regDeltaXYZ;

        /* \brief Number of regions in each of 3 directions
         */
      int regNums[3];

        /* \brief Map from source processor to bounding box of that proc's source mesh
         *
         */
      std::map<int, BoundBox> srcProcBoxes;

        /* \brief Timing object for spatial location
         */
      SpatialLocatorTimes myTimes;

        /* \brief Timer object to manage overloaded search functions
         */
      CpuTimer myTimer;
      
        /* \brief Flag to manage initialization of timer for overloaded search functions
         */
      bool timerInitialized;
    };

    inline SpatialLocator::~SpatialLocator() 
    {
      if (iCreatedTree && myTree) delete myTree;
    }
    
    inline ErrorCode SpatialLocator::locate_point(const double *pos, 
                                                  EntityHandle &ent, double *params, int *is_inside, 
                                                  const double rel_iter_tol, const double abs_iter_tol,
                                                  const double inside_tol)
    {
      return locate_points(pos, 1, &ent, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol);
    }

#ifdef USE_MPI
    inline ErrorCode SpatialLocator::get_point_ijk(const CartVect &point, const double abs_iter_tol, int *ijk) const
    {
      for (int i = 0; i < 3; i++) {
        if (point[i] < globalBox.bMin[i]-abs_iter_tol || point[i] > globalBox.bMax[i]+abs_iter_tol)
          ijk[i] = -1;
        else {
          ijk[i] = point[i] - globalBox.bMin[i] / regDeltaXYZ[i];
          if (ijk[i] >= regNums[i] && point[i] <= globalBox.bMax[i]+abs_iter_tol)
            ijk[i] = regNums[i]-1;
        }
      }
      
      return (ijk[0] >= 0 && ijk[1] >= 0 && ijk[2] >= 0 ? MB_SUCCESS : MB_FAILURE);;
    }

    inline int SpatialLocator::proc_from_ijk(const int *ijk) const
    {
      return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
    }
    
    inline int SpatialLocator::proc_from_point(const double *pos, const double abs_iter_tol) const
    {
      int ijk[3];
      ErrorCode rval = get_point_ijk(CartVect(pos), abs_iter_tol, ijk);
      if (MB_SUCCESS != rval) return -1;
      
      return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
    }
#endif
    
} // namespace moab 

#endif