Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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,
00082                              EntityHandle* ents,
00083                              double* params,
00084                              int* is_inside            = NULL,
00085                              const double rel_iter_tol = 1.0e-10,
00086                              const double abs_iter_tol = 1.0e-10,
00087                              const double inside_tol   = 1.0e-6 );
00088 
00089     /* locate a set of points */
00090     ErrorCode locate_points( const double* pos,
00091                              int num_points,
00092                              EntityHandle* ents,
00093                              double* params,
00094                              int* is_inside            = NULL,
00095                              const double rel_iter_tol = 1.0e-10,
00096                              const double abs_iter_tol = 1.0e-10,
00097                              const double inside_tol   = 1.0e-6 );
00098 
00099     /* locate a set of vertices or entity centroids, storing results on TupleList in this class
00100      * Locate a set of vertices or entity centroids, storing the detailed results in member
00101      * variable (TupleList) locTable (see comments on locTable for structure of that tuple).
00102      */
00103     ErrorCode locate_points( Range& ents,
00104                              const double rel_iter_tol = 1.0e-10,
00105                              const double abs_iter_tol = 1.0e-10,
00106                              const double inside_tol   = 1.0e-6 );
00107 
00108     /* locate a set of points, storing results on TupleList in this class
00109      * Locate a set of points, storing the detailed results in member variable (TupleList) locTable
00110      * (see comments on locTable for structure of that tuple).
00111      */
00112     ErrorCode locate_points( const double* pos,
00113                              int num_points,
00114                              const double rel_iter_tol = 1.0e-10,
00115                              const double abs_iter_tol = 1.0e-10,
00116                              const double inside_tol   = 1.0e-6 );
00117 
00118     /* Count the number of located points in locTable
00119      * Return the number of entries in locTable that have non-zero entity handles, which
00120      * represents the number of points in targetEnts that were inside one element in sourceEnts
00121      *
00122      */
00123     int local_num_located();
00124 
00125     /* Count the number of located points in parLocTable
00126      * Return the number of entries in parLocTable that have a non-negative index in on a remote
00127      * proc in parLocTable, which gives the number of points located in at least one element in a
00128      * remote proc's sourceEnts.
00129      */
00130     int remote_num_located();
00131 
00132 #ifdef MOAB_HAVE_MPI
00133     /* locate a set of vertices or entity centroids, storing results on TupleList in this class
00134      * Locate a set of vertices or entity centroids, storing the detailed results in member
00135      * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
00136      * structure of those tuples).
00137      */
00138     ErrorCode par_locate_points( ParallelComm* pc,
00139                                  Range& vertices,
00140                                  const double rel_iter_tol = 1.0e-10,
00141                                  const double abs_iter_tol = 1.0e-10,
00142                                  const double inside_tol   = 1.0e-6 );
00143 
00144     /* locate a set of points, storing results on TupleList in this class
00145      * Locate a set of points, storing the detailed results in member
00146      * variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
00147      * structure of those tuples).
00148      */
00149     ErrorCode par_locate_points( ParallelComm* pc,
00150                                  const double* pos,
00151                                  int num_points,
00152                                  const double rel_iter_tol = 1.0e-10,
00153                                  const double abs_iter_tol = 1.0e-10,
00154                                  const double inside_tol   = 1.0e-6 );
00155 #endif
00156 
00157     /** \brief Return the MOAB interface associated with this locator
00158      */
00159     Interface* moab()
00160     {
00161         return mbImpl;
00162     }
00163 
00164     /* locate a point */
00165     ErrorCode locate_point( const double* pos,
00166                             EntityHandle& ent,
00167                             double* params,
00168                             int* is_inside            = NULL,
00169                             const double rel_iter_tol = 1.0e-10,
00170                             const double abs_iter_tol = 1.0e-10,
00171                             const double inside_tol   = 1.0e-6 );
00172 
00173     /* return the tree */
00174     Tree* get_tree()
00175     {
00176         return myTree;
00177     }
00178 
00179     /* get the locTable
00180      */
00181     TupleList& loc_table()
00182     {
00183         return locTable;
00184     }
00185 
00186     /* get the locTable
00187      */
00188     const TupleList& loc_table() const
00189     {
00190         return locTable;
00191     }
00192 
00193     /* get the parLocTable
00194      */
00195     TupleList& par_loc_table()
00196     {
00197         return parLocTable;
00198     }
00199 
00200     /* get the parLocTable
00201      */
00202     const TupleList& par_loc_table() const
00203     {
00204         return parLocTable;
00205     }
00206 
00207     /* get elemEval */
00208     ElemEvaluator* elem_eval()
00209     {
00210         return elemEval;
00211     }
00212 
00213     /* get elemEval */
00214     const ElemEvaluator* elem_eval() const
00215     {
00216         return elemEval;
00217     }
00218 
00219     /* set elemEval */
00220     void elem_eval( ElemEvaluator* eval )
00221     {
00222         elemEval = eval;
00223         if( myTree ) myTree->set_eval( eval );
00224     }
00225 
00226     /** \brief Get spatial locator times object */
00227     SpatialLocatorTimes& sl_times()
00228     {
00229         return myTimes;
00230     }
00231 
00232     /** \brief Get spatial locator times object */
00233     const SpatialLocatorTimes& sl_times() const
00234     {
00235         return myTimes;
00236     }
00237 
00238   private:
00239 #ifdef MOAB_HAVE_MPI
00240     /* MPI_ReduceAll source mesh bounding boxes to get global source mesh bounding box
00241      */
00242     ErrorCode initialize_intermediate_partition( ParallelComm* pc );
00243 
00244     /* for a given point in space, compute its ijk location in the intermediate decomposition;
00245      * tolerance is used only to test proximity to global box extent, not for local box test
00246      */
00247     inline ErrorCode get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const;
00248 
00249 #if 0
00250         /* given an ijk location in the intermediate partition, return the proc rank for that location
00251          */
00252       inline int proc_from_ijk(const int *ijk) const;
00253 #endif
00254 
00255     /* given a point in space, return the proc responsible for that point from the intermediate
00256      * decomp; no tolerances applied here, so first proc in lexicographic ijk ordering is returned
00257      */
00258     inline int proc_from_point( const double* pos, const double abs_iter_tol ) const;
00259 
00260     /* register my source mesh with intermediate decomposition procs
00261      */
00262     ErrorCode register_src_with_intermediate_procs( ParallelComm* pc, double abs_iter_tol, TupleList& TLreg_o );
00263 
00264 #endif
00265 
00266     /** Create a tree
00267      * Tree type depends on what's in myElems: if empty or all vertices, creates a kdtree,
00268      * otherwise creates a BVHTree.
00269      */
00270     void create_tree();
00271 
00272     /* MOAB instance */
00273     Interface* mbImpl;
00274 
00275     /* elements being located */
00276     Range myElems;
00277 
00278     /* dimension of entities in locator */
00279     int myDim;
00280 
00281     /* tree used for location */
00282     Tree* myTree;
00283 
00284     /* element evaluator */
00285     ElemEvaluator* elemEval;
00286 
00287     /* whether I created the tree or not (determines whether to delete it or not on destruction) */
00288     bool iCreatedTree;
00289 
00290     /* \brief local locations table
00291      * This table stores detailed local location data results from locate_points, that is, location
00292      * data for points located on the local mesh.  Data is stored in a TupleList, where each tuple
00293      * consists of (p_i, hs_ul, r[3]_d), where p_i = (int) proc from which request for this point
00294      * was made (0 if serial) hs_ul = (unsigned long) source entity containing the point r[3]_d =
00295      * (double) parametric coordinates of the point in hs
00296      */
00297     TupleList locTable;
00298 
00299     /* \brief parallel locations table
00300      * This table stores information about points located on a local or remote processor.  For
00301      * points located on this processor's local mesh, detailed location data is stored in locTable.
00302      * For points located on remote processors, more communication is required to retrieve specific
00303      * location data (usually that information isn't useful on this processor).
00304      *
00305      * The tuple structure of this TupleList is (p_i, ri_i), where:
00306      *   p_i = (int) processor rank containing this point
00307      *   ri_i = (int) index into locTable on remote proc containing this point's location
00308      * information The indexing of parLocTable corresponds to that of the points/entities passed in.
00309      */
00310     TupleList parLocTable;
00311 
00312     /* \brief Local bounding box, duplicated from tree
00313      */
00314     BoundBox localBox;
00315 
00316     /* \brief Global bounding box, used in parallel spatial location
00317      */
00318     BoundBox globalBox;
00319 
00320     /* \brief Regional delta xyz, used in parallel spatial location
00321      */
00322     CartVect regDeltaXYZ;
00323 
00324     /* \brief Number of regions in each of 3 directions
00325      */
00326     int regNums[3];
00327 
00328     /* \brief Map from source processor to bounding box of that proc's source mesh
00329      *
00330      */
00331     std::map< int, BoundBox > srcProcBoxes;
00332 
00333     /* \brief Timing object for spatial location
00334      */
00335     SpatialLocatorTimes myTimes;
00336 
00337     /* \brief Timer object to manage overloaded search functions
00338      */
00339     CpuTimer myTimer;
00340 
00341     /* \brief Flag to manage initialization of timer for overloaded search functions
00342      */
00343     bool timerInitialized;
00344 };
00345 
00346 inline SpatialLocator::~SpatialLocator()
00347 {
00348     if( iCreatedTree && myTree ) delete myTree;
00349 }
00350 
00351 inline ErrorCode SpatialLocator::locate_point( const double* pos,
00352                                                EntityHandle& ent,
00353                                                double* params,
00354                                                int* is_inside,
00355                                                const double rel_iter_tol,
00356                                                const double abs_iter_tol,
00357                                                const double inside_tol )
00358 {
00359     return locate_points( pos, 1, &ent, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol );
00360 }
00361 
00362 #ifdef MOAB_HAVE_MPI
00363 inline ErrorCode SpatialLocator::get_point_ijk( const CartVect& point, const double abs_iter_tol, int* ijk ) const
00364 {
00365     for( int i = 0; i < 3; i++ )
00366     {
00367         if( point[i] < globalBox.bMin[i] - abs_iter_tol || point[i] > globalBox.bMax[i] + abs_iter_tol )
00368             ijk[i] = -1;
00369         else
00370         {
00371             ijk[i] = point[i] - globalBox.bMin[i] / regDeltaXYZ[i];
00372             if( ijk[i] >= regNums[i] && point[i] <= globalBox.bMax[i] + abs_iter_tol ) ijk[i] = regNums[i] - 1;
00373         }
00374     }
00375 
00376     return ( ijk[0] >= 0 && ijk[1] >= 0 && ijk[2] >= 0 ? MB_SUCCESS : MB_FAILURE );
00377     ;
00378 }
00379 
00380 #if 0
00381     inline int SpatialLocator::proc_from_ijk(const int *ijk) const
00382     {
00383       return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
00384     }
00385 #endif
00386 
00387 inline int SpatialLocator::proc_from_point( const double* pos, const double abs_iter_tol ) const
00388 {
00389     int ijk[3];
00390     ErrorCode rval = get_point_ijk( CartVect( pos ), abs_iter_tol, ijk );
00391     if( MB_SUCCESS != rval ) return -1;
00392 
00393     return ijk[2] * regNums[0] * regNums[1] + ijk[1] * regNums[0] + ijk[0];
00394 }
00395 #endif
00396 
00397 }  // namespace moab
00398 
00399 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines