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