SpatialLocator.hpp

Go to the documentation of this file.

/**\file SpatialLocator.hpp
 * \class moab::SpatialLocator
 * \brief Tool to facilitate spatial location of a point in a mesh
 *
 * SpatialLocator facilitates searching for points in or performing ray traces on collections of
 * mesh entities in 2D or 3D.  This searching is facilitated by a tree-based decomposition of the
 * mesh.  Various types of trees are implemented in MOAB and can be used by this tool, but by
 * default it uses AdaptiveKDTree (see child classes of Tree for which others are available).
 * Parallel and serial searching are both supported.
 *
 * SpatialLocator can either cache the search results for points or pass back this information in
 * arguments. Cached information is kept in locTable, indexed in the same order as search points
 * passed in.  This information consists of the entity handle containing the point and the
 * parametric coordinates inside that element. Information about the points searched, e.g. the
 * entities from which those points are derived, can be stored in the calling application if
 * desired.
 *
 * In parallel, there is a separation between the proc deciding which points to search for (the
 * "target" proc), and the proc locating the point in its local mesh (the "source" proc).  On the
 * source proc, location information is cached in locTable, as in the serial case.  By default, this
 * location information (handle and parametric coords) is not returned to the target proc, since it
 * would be of no use there.  Instead, the rank of the source proc locating the point, and the index
 * of that location info in the source proc's locTable, is returned; this information is stored on
 * the target proc in this class's parLocTable variable.  Again, information about the points
 * searched should be stored in the calling application, if desired.
 *
 * This class uses the ElemEvaluator class for specification and evaluation of basis functions (used
 * for computing parametric coords within an entity).  See documentation and examples for that class
 * for usage information.
 *
 */
 
#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 <cmath>
 
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 MOAB_HAVE_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
 
    /** \brief Return the MOAB interface associated with this locator
     */
    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 );
    }
 
    /** \brief Get spatial locator times object */
    SpatialLocatorTimes& sl_times()
    {
        return myTimes;
    }
 
    /** \brief Get spatial locator times object */
    const SpatialLocatorTimes& sl_times() const
    {
        return myTimes;
    }
 
  private:
#ifdef MOAB_HAVE_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;
 
#if 0
        /* given an ijk location in the intermediate partition, return the proc rank for that location
         */
      inline int proc_from_ijk(const int *ijk) const;
#endif
 
    /* 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
 
    /** Create a tree
     * Tree type depends on what's in myElems: if empty or all vertices, creates a kdtree,
     * otherwise creates a BVHTree.
     */
    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 MOAB_HAVE_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 );
    ;
}
 
#if 0
    inline int SpatialLocator::proc_from_ijk(const int *ijk) const
    {
      return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
    }
#endif
 
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