MOAB: Mesh Oriented datABase  (version 5.2.1)
SpatialLocator.hpp
Go to the documentation of this file.
00001 /**\file SpatialLocator.hpp
00002  * \class moab::SpatialLocator
00003  * \brief Tool to facilitate spatial location of a point in a mesh
00004  *
00005  * SpatialLocator facilitates searching for points in or performing ray traces on collections of
00006  * mesh entities in 2D or 3D.  This searching is facilitated by a tree-based decomposition of the
00007  * mesh.  Various types of trees are implemented in MOAB and can be used by this tool, but by
00008  * default it uses AdaptiveKDTree (see child classes of Tree for which others are available).
00009  * Parallel and serial searching are both supported.
00010  *
00011  * SpatialLocator can either cache the search results for points or pass back this information in
00012  * arguments. Cached information is kept in locTable, indexed in the same order as search points
00013  * passed in.  This information consists of the entity handle containing the point and the
00014  * parametric coordinates inside that element. Information about the points searched, e.g. the
00015  * entities from which those points are derived, can be stored in the calling application if
00016  * desired.
00017  *
00018  * In parallel, there is a separation between the proc deciding which points to search for (the
00019  * "target" proc), and the proc locating the point in its local mesh (the "source" proc).  On the
00020  * source proc, location information is cached in locTable, as in the serial case.  By default, this
00021  * location information (handle and parametric coords) is not returned to the target proc, since it
00022  * would be of no use there.  Instead, the rank of the source proc locating the point, and the index
00023  * of that location info in the source proc's locTable, is returned; this information is stored on
00024  * the target proc in this class's parLocTable variable.  Again, information about the points
00025  * searched should be stored in the calling application, if desired.
00026  *
00027  * This class uses the ElemEvaluator class for specification and evaluation of basis functions (used
00028  * for computing parametric coords within an entity).  See documentation and examples for that class
00029  * for usage information.
00030  *
00031  */
00032 
00033 #ifndef MOAB_SPATIALLOCATOR_HPP
00034 #define MOAB_SPATIALLOCATOR_HPP
00035 
00036 #include "moab/Types.hpp"
00037 #include "moab/Tree.hpp"
00038 #include "moab/Range.hpp"
00039 #include "moab/TupleList.hpp"
00040 #include "moab/BoundBox.hpp"
00041 #include "moab/SpatialLocatorTimes.hpp"
00042 #include "moab/CpuTimer.hpp"
00043 
00044 #include <string>
00045 #include <vector>
00046 #include <map>
00047 #include <cmath>
00048 
00049 namespace moab
00050 {
00051 
00052 class Interface;
00053 class ElemEvaluator;
00054 class ParallelComm;
00055 
00056 class SpatialLocator
00057 {
00058   public:
00059     /* constructor */
00060     SpatialLocator( Interface* impl, Range& elems, Tree* tree = NULL, ElemEvaluator* eval = NULL );
00061 
00062     /* destructor */
00063     virtual ~SpatialLocator();
00064 
00065     /* add elements to be searched */
00066     ErrorCode add_elems( Range& elems );
00067 
00068     /* get bounding box of this locator */
00069     BoundBox& local_box()
00070     {
00071         return localBox;
00072     }
00073 
00074     /* get bounding box of this locator */
00075     const BoundBox& local_box() const
00076     {
00077         return localBox;
00078     }
00079 
00080     /* locate a set of vertices, Range variant */
00081     ErrorCode locate_points( Range& vertices, EntityHandle* ents, double* params, int* is_inside = NULL,
00082                              const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
00083                              const double inside_tol = 1.0e-6 );
00084 
00085     /* locate a set of points */
00086     ErrorCode locate_points( const double* pos, int num_points, EntityHandle* ents, double* params,
00087                              int* is_inside = NULL, const double rel_iter_tol = 1.0e-10,
00088                              const double abs_iter_tol = 1.0e-10, const double inside_tol = 1.0e-6 );
00089 
00090     /* locate a set of vertices or entity centroids, storing results on TupleList in this class
00091      * Locate a set of vertices or entity centroids, storing the detailed results in member
00092      * variable (TupleList) locTable (see comments on locTable for structure of that tuple).
00093      */
00094     ErrorCode locate_points( Range& ents, const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
00095                              const double inside_tol = 1.0e-6 );
00096 
00097     /* locate a set of points, storing results on TupleList in this class
00098      * Locate a set of points, storing the detailed results in member variable (TupleList) locTable
00099      * (see comments on locTable for structure of that tuple).
00100      */
00101     ErrorCode locate_points( const double* pos, int num_points, const double rel_iter_tol = 1.0e-10,
00102                              const double abs_iter_tol = 1.0e-10, const double inside_tol = 1.0e-6 );
00103 
00104     /* Count the number of located points in locTable
00105      * Return the number of entries in locTable that have non-zero entity handles, which
00106      * represents the number of points in targetEnts that were inside one element in sourceEnts
00107      *
00108      */
00109     int local_num_located();
00110 
00111     /* Count the number of located points in parLocTable
00112      * Return the number of entries in parLocTable that have a non-negative index in on a remote
00113      * proc in parLocTable, which gives the number of points located in at least one element in a
00114      * remote proc's sourceEnts.
00115      */
00116     int remote_num_located();
00117 
00118 #ifdef MOAB_HAVE_MPI
00119     /* locate a set of vertices or entity centroids, storing results on TupleList in this class
00120      * Locate a set of vertices or entity centroids, storing the detailed results in member
00121      * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
00122      * structure of those tuples).
00123      */
00124     ErrorCode par_locate_points( ParallelComm* pc, Range& vertices, const double rel_iter_tol = 1.0e-10,
00125                                  const double abs_iter_tol = 1.0e-10, const double inside_tol = 1.0e-6 );
00126 
00127     /* locate a set of points, storing results on TupleList in this class
00128      * Locate a set of points, storing the detailed results in member
00129      * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
00130      * structure of those tuples).
00131      */
00132     ErrorCode par_locate_points( ParallelComm* pc, const double* pos, int num_points,
00133                                  const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
00134                                  const double inside_tol = 1.0e-6 );
00135 #endif
00136 
00137     /** \brief Return the MOAB interface associated with this locator
00138      */
00139     Interface* moab()
00140     {
00141         return mbImpl;
00142     }
00143 
00144     /* locate a point */
00145     ErrorCode locate_point( const double* pos, EntityHandle& ent, double* params, int* is_inside = NULL,
00146                             const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
00147                             const double inside_tol = 1.0e-6 );
00148 
00149     /* return the tree */
00150     Tree* get_tree()
00151     {
00152         return myTree;
00153     }
00154 
00155     /* get the locTable
00156      */
00157     TupleList& loc_table()
00158     {
00159         return locTable;
00160     }
00161 
00162     /* get the locTable
00163      */
00164     const TupleList& loc_table() const
00165     {
00166         return locTable;
00167     }
00168 
00169     /* get the parLocTable
00170      */
00171     TupleList& par_loc_table()
00172     {
00173         return parLocTable;
00174     }
00175 
00176     /* get the parLocTable
00177      */
00178     const TupleList& par_loc_table() const
00179     {
00180         return parLocTable;
00181     }
00182 
00183     /* get elemEval */
00184     ElemEvaluator* elem_eval()
00185     {
00186         return elemEval;
00187     }
00188 
00189     /* get elemEval */
00190     const ElemEvaluator* elem_eval() const
00191     {
00192         return elemEval;
00193     }
00194 
00195     /* set elemEval */
00196     void elem_eval( ElemEvaluator* eval )
00197     {
00198         elemEval = eval;
00199         if( myTree ) myTree->set_eval( eval );
00200     }
00201 
00202     /** \brief Get spatial locator times object */
00203     SpatialLocatorTimes& sl_times()
00204     {
00205         return myTimes;
00206     }
00207 
00208     /** \brief Get spatial locator times object */
00209     const SpatialLocatorTimes& sl_times() const
00210     {
00211         return myTimes;
00212     }
00213 
00214   private:
00215 #ifdef MOAB_HAVE_MPI
00216     /* MPI_ReduceAll source mesh bounding boxes to get global source mesh bounding box
00217      */
00218     ErrorCode initialize_intermediate_partition( ParallelComm* pc );
00219 
00220     /* for a given point in space, compute its ijk location in the intermediate decomposition;
00221      * tolerance is used only to test proximity to global box extent, not for local box test
00222      */
00223     inline ErrorCode get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const;
00224 
00225 #if 0
00226         /* given an ijk location in the intermediate partition, return the proc rank for that location
00227          */
00228       inline int proc_from_ijk(const int *ijk) const;
00229 #endif
00230 
00231     /* given a point in space, return the proc responsible for that point from the intermediate
00232      * decomp; no tolerances applied here, so first proc in lexicographic ijk ordering is returned
00233      */
00234     inline int proc_from_point( const double* pos, const double abs_iter_tol ) const;
00235 
00236     /* register my source mesh with intermediate decomposition procs
00237      */
00238     ErrorCode register_src_with_intermediate_procs( ParallelComm* pc, double abs_iter_tol, TupleList& TLreg_o );
00239 
00240 #endif
00241 
00242     /** Create a tree
00243      * Tree type depends on what's in myElems: if empty or all vertices, creates a kdtree,
00244      * otherwise creates a BVHTree.
00245      */
00246     void create_tree();
00247 
00248     /* MOAB instance */
00249     Interface* mbImpl;
00250 
00251     /* elements being located */
00252     Range myElems;
00253 
00254     /* dimension of entities in locator */
00255     int myDim;
00256 
00257     /* tree used for location */
00258     Tree* myTree;
00259 
00260     /* element evaluator */
00261     ElemEvaluator* elemEval;
00262 
00263     /* whether I created the tree or not (determines whether to delete it or not on destruction) */
00264     bool iCreatedTree;
00265 
00266     /* \brief local locations table
00267      * This table stores detailed local location data results from locate_points, that is, location
00268      * data for points located on the local mesh.  Data is stored in a TupleList, where each tuple
00269      * consists of (p_i, hs_ul, r[3]_d), where p_i = (int) proc from which request for this point
00270      * was made (0 if serial) hs_ul = (unsigned long) source entity containing the point r[3]_d =
00271      * (double) parametric coordinates of the point in hs
00272      */
00273     TupleList locTable;
00274 
00275     /* \brief parallel locations table
00276      * This table stores information about points located on a local or remote processor.  For
00277      * points located on this processor's local mesh, detailed location data is stored in locTable.
00278      * For points located on remote processors, more communication is required to retrieve specific
00279      * location data (usually that information isn't useful on this processor).
00280      *
00281      * The tuple structure of this TupleList is (p_i, ri_i), where:
00282      *   p_i = (int) processor rank containing this point
00283      *   ri_i = (int) index into locTable on remote proc containing this point's location
00284      * information The indexing of parLocTable corresponds to that of the points/entities passed in.
00285      */
00286     TupleList parLocTable;
00287 
00288     /* \brief Local bounding box, duplicated from tree
00289      */
00290     BoundBox localBox;
00291 
00292     /* \brief Global bounding box, used in parallel spatial location
00293      */
00294     BoundBox globalBox;
00295 
00296     /* \brief Regional delta xyz, used in parallel spatial location
00297      */
00298     CartVect regDeltaXYZ;
00299 
00300     /* \brief Number of regions in each of 3 directions
00301      */
00302     int regNums[3];
00303 
00304     /* \brief Map from source processor to bounding box of that proc's source mesh
00305      *
00306      */
00307     std::map< int, BoundBox > srcProcBoxes;
00308 
00309     /* \brief Timing object for spatial location
00310      */
00311     SpatialLocatorTimes myTimes;
00312 
00313     /* \brief Timer object to manage overloaded search functions
00314      */
00315     CpuTimer myTimer;
00316 
00317     /* \brief Flag to manage initialization of timer for overloaded search functions
00318      */
00319     bool timerInitialized;
00320 };
00321 
00322 inline SpatialLocator::~SpatialLocator()
00323 {
00324     if( iCreatedTree && myTree ) delete myTree;
00325 }
00326 
00327 inline ErrorCode SpatialLocator::locate_point( const double* pos, EntityHandle& ent, double* params, int* is_inside,
00328                                                const double rel_iter_tol, const double abs_iter_tol,
00329                                                const double inside_tol )
00330 {
00331     return locate_points( pos, 1, &ent, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol );
00332 }
00333 
00334 #ifdef MOAB_HAVE_MPI
00335 inline ErrorCode SpatialLocator::get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const
00336 {
00337     for( int i = 0; i < 3; i++ )
00338     {
00339         if( point[i] < globalBox.bMin[i] - abs_iter_tol || point[i] > globalBox.bMax[i] + abs_iter_tol )
00340             ijk[i] = -1;
00341         else
00342         {
00343             ijk[i] = point[i] - globalBox.bMin[i] / regDeltaXYZ[i];
00344             if( ijk[i] >= regNums[i] && point[i] <= globalBox.bMax[i] + abs_iter_tol ) ijk[i] = regNums[i] - 1;
00345         }
00346     }
00347 
00348     return ( ijk[0] >= 0 && ijk[1] >= 0 && ijk[2] >= 0 ? MB_SUCCESS : MB_FAILURE );
00349     ;
00350 }
00351 
00352 #if 0
00353     inline int SpatialLocator::proc_from_ijk(const int *ijk) const
00354     {
00355       return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
00356     }
00357 #endif
00358 
00359 inline int SpatialLocator::proc_from_point( const double* pos, const double abs_iter_tol ) const
00360 {
00361     int ijk[3];
00362     ErrorCode rval = get_point_ijk( CartVect( pos ), abs_iter_tol, ijk );
00363     if( MB_SUCCESS != rval ) return -1;
00364 
00365     return ijk[2] * regNums[0] * regNums[1] + ijk[1] * regNums[0] + ijk[0];
00366 }
00367 #endif
00368 
00369 }  // namespace moab
00370 
00371 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines