point_locater.hpp

Go to the documentation of this file.

    1 /**
    2  * point_locater.hpp
    3  * Ryan H. Lewis
    4  * Copyright 2012
    5  */
    6 #include <vector>
    7 #include <iostream>
    8 #include <ctime>
    9  
   10 #include <sstream>
   11 #include <fstream>
   12  
   13 #ifndef POINT_LOCATER_HPP
   14 #define POINT_LOCATER_HPP
   15  
   16 namespace moab
   17 {
   18  
   19 template < typename _Tree, typename _Boxes >
   20 class Point_search
   21 {
   22  
   23     // public types
   24   public:
   25     typedef _Tree Tree;
   26     typedef _Boxes Boxes;
   27     //  typedef typename Tree::Elements::value_type Element;
   28     typedef int Error;
   29  
   30     // private types
   31   private:
   32     typedef Point_search< _Tree, _Boxes > Self;
   33     // public methods
   34   public:
   35     // Constructor
   36     Point_search( Tree& _tree, Boxes& _boxes ) : tree_( _tree ), boxes( _boxes ) {}
   37  
   38     // Copy constructor
   39     Point_search( Self& s ) : tree_( s.tree_ ), boxes( s.boxes ) {}
   40  
   41     // private functionality
   42   private:
   43     // TODO: deprecate this
   44     template < typename Point_map, typename List >
   45     void resolve_boxes( const Point_map& query_points, List& list )
   46     {
   47         /*
   48      typedef typename std::vector< bool> Bitmask;
   49      typedef typename Points::const_iterator Point;
   50      typedef typename Tree::const_element_iterator Element;
   51      typedef typename std::vector< std::size_t> Processor_list;
   52      typedef typename List::value_type Tuple;
   53      const Element & end = tree_.element_end();
   54      Bitmask located_points( query_points.size(), 0);
   55      std::size_t index=0;
   56      for( Point i = query_points.begin();
   57             i != query_points.end(); ++i,++index){
   58          const Element  & element = tree_.find( *i);
   59          if(element != end){
   60              located_points[ index] = 1;
   61          }
   62      }
   63      for(int i = 0; i < located_points.size(); ++i){
   64          if(!located_points[ i]){
   65              Processor_list processors;
   66              const Point & point = query_point.begin()+i;
   67              resolve_box_for_point( point, processors);
   68              for( std::size_t p = processors.begin();
   69                       p != processors.end(); ++p){
   70                  list.push_back( Tuple( *point, *p) );
   71              }
   72          }
   73      }
   74      */
   75     }
   76  
   77     // public functionality
   78   public:
   79     template < typename Point_map, typename Entities, typename Communicator >
   80     Error locate_points( Point_map& query_points, Entities& entities, Communicator& comm, double tol )
   81     {
   82         /*TODO: implement a parallel location algorithm here
   83         original algorithm:
   84         all-scatter/gather bounding boxes
   85         for each point perform box checks
   86         scatter/gather points
   87         perform location and solution interpolation
   88         //send results back
   89         ------------
   90         new algorithm
   91         Assume power num_proc = n^2;
   92         all-reduce bounding boxes
   93         divide box up into n x n grid
   94         each source processor can communicate its elements to grid proc
   95         each target processor then communicates its points to grid proc
   96         grid proc sends points to correct source proc.
   97         */
   98         return 0;
   99     }
  100  
  101     template < typename Points, typename Entities >
  102     Error locate_points( const Points& query_points, Entities& entities, double tol ) const
  103     {
  104         typedef typename Points::const_iterator Point_iterator;
  105         typedef typename Entities::value_type Result;
  106         Entities result;
  107         result.reserve( query_points.size() );
  108         for( Point_iterator i = query_points.begin(); i != query_points.end(); ++i )
  109         {
  110             Result h;
  111             tree_.find( *i, tol, h );
  112             result.push_back( h );
  113         }
  114         entities = result;
  115         return 0;
  116     }
  117  
  118     template < typename Points, typename Entities >
  119     Error bruteforce_locate_points( const Points& query_points, Entities& entities, double tol ) const
  120     {
  121         // TODO: this could be faster with caching, but of course this is
  122         // really just for testing
  123         typedef typename Points::const_iterator Point_iterator;
  124         typedef typename Entities::value_type::first_type Entity_handle;
  125         Entities result;
  126         result.reserve( query_points.size() );
  127         std::size_t count = 0;
  128         std::stringstream ss;
  129         typename Entities::iterator j = entities.begin();
  130         for( Point_iterator i = query_points.begin(); i != query_points.end(); ++i, ++j )
  131         {
  132             if( j->first == 0 )
  133             {
  134                 const Entity_handle h = tree_.bruteforce_find( *i, tol );
  135                 if( h == 0 )
  136                 {
  137                     ++count;
  138                     for( int k = 0; k < 3; ++k )
  139                     {
  140                         ss << ( *i )[k];
  141                         if( k < 2 )
  142                         {
  143                             ss << ", ";
  144                         }
  145                         else
  146                         {
  147                             ss << std::endl;
  148                         }
  149                     }
  150                 }
  151             }
  152         }
  153         std::ofstream out( "unlocated_pts" );
  154         out << ss.str();
  155         out.close();
  156         std::cout << count << " vertices are not contained in _any_ elements!" << std::endl;
  157         return 0;
  158     }
  159  
  160     // public accessor methods
  161   public:
  162     Tree& tree() const
  163     {
  164         return tree_;
  165     }
  166  
  167     // private data members
  168   private:
  169     const Tree& tree_;
  170     const Boxes& boxes;
  171 };  // class Point_search
  172  
  173 }  // namespace moab
  174  
  175 #endif  // POINT_LOCATER_HPP