point_locater.hpp

Go to the documentation of this file.

/**
 * point_locater.hpp
 * Ryan H. Lewis
 * Copyright 2012
 */
#include <vector>
#include <iostream>
#include <ctime>
 
#include <sstream>
#include <fstream>
 
#ifndef POINT_LOCATER_HPP
#define POINT_LOCATER_HPP
 
namespace moab
{
 
template < typename _Tree, typename _Boxes >
class Point_search
{
 
    // public types
  public:
    typedef _Tree Tree;
    typedef _Boxes Boxes;
    //  typedef typename Tree::Elements::value_type Element;
    typedef int Error;
 
    // private types
  private:
    typedef Point_search< _Tree, _Boxes > Self;
    // public methods
  public:
    // Constructor
    Point_search( Tree& _tree, Boxes& _boxes ) : tree_( _tree ), boxes( _boxes ) {}
 
    // Copy constructor
    Point_search( Self& s ) : tree_( s.tree_ ), boxes( s.boxes ) {}
 
    // private functionality
  private:
    // TODO: deprecate this
    template < typename Point_map, typename List >
    void resolve_boxes( const Point_map& query_points, List& list )
    {
        /*
     typedef typename std::vector< bool> Bitmask;
     typedef typename Points::const_iterator Point;
     typedef typename Tree::const_element_iterator Element;
     typedef typename std::vector< std::size_t> Processor_list;
     typedef typename List::value_type Tuple;
     const Element & end = tree_.element_end();
     Bitmask located_points( query_points.size(), 0);
     std::size_t index=0;
     for( Point i = query_points.begin();
            i != query_points.end(); ++i,++index){
         const Element  & element = tree_.find( *i);
         if(element != end){
             located_points[ index] = 1;
         }
     }
     for(int i = 0; i < located_points.size(); ++i){
         if(!located_points[ i]){
             Processor_list processors;
             const Point & point = query_point.begin()+i;
             resolve_box_for_point( point, processors);
             for( std::size_t p = processors.begin();
                      p != processors.end(); ++p){
                 list.push_back( Tuple( *point, *p) );
             }
         }
     }
     */
    }
 
    // public functionality
  public:
    template < typename Point_map, typename Entities, typename Communicator >
    Error locate_points( Point_map& query_points, Entities& entities, Communicator& comm, double tol )
    {
        /*TODO: implement a parallel location algorithm here
        original algorithm:
        all-scatter/gather bounding boxes
        for each point perform box checks
        scatter/gather points
        perform location and solution interpolation
        //send results back
        ------------
        new algorithm
        Assume power num_proc = n^2;
        all-reduce bounding boxes
        divide box up into n x n grid
        each source processor can communicate its elements to grid proc
        each target processor then communicates its points to grid proc
        grid proc sends points to correct source proc.
        */
        return 0;
    }
 
    template < typename Points, typename Entities >
    Error locate_points( const Points& query_points, Entities& entities, double tol ) const
    {
        typedef typename Points::const_iterator Point_iterator;
        typedef typename Entities::value_type Result;
        Entities result;
        result.reserve( query_points.size() );
        for( Point_iterator i = query_points.begin(); i != query_points.end(); ++i )
        {
            Result h;
            tree_.find( *i, tol, h );
            result.push_back( h );
        }
        entities = result;
        return 0;
    }
 
    template < typename Points, typename Entities >
    Error bruteforce_locate_points( const Points& query_points, Entities& entities, double tol ) const
    {
        // TODO: this could be faster with caching, but of course this is
        // really just for testing
        typedef typename Points::const_iterator Point_iterator;
        typedef typename Entities::value_type::first_type Entity_handle;
        Entities result;
        result.reserve( query_points.size() );
        std::size_t count = 0;
        std::stringstream ss;
        typename Entities::iterator j = entities.begin();
        for( Point_iterator i = query_points.begin(); i != query_points.end(); ++i, ++j )
        {
            if( j->first == 0 )
            {
                const Entity_handle h = tree_.bruteforce_find( *i, tol );
                if( h == 0 )
                {
                    ++count;
                    for( int k = 0; k < 3; ++k )
                    {
                        ss << ( *i )[k];
                        if( k < 2 )
                        {
                            ss << ", ";
                        }
                        else
                        {
                            ss << std::endl;
                        }
                    }
                }
            }
        }
        std::ofstream out( "unlocated_pts" );
        out << ss.str();
        out.close();
        std::cout << count << " vertices are not contained in _any_ elements!" << std::endl;
        return 0;
    }
 
    // public accessor methods
  public:
    Tree& tree() const
    {
        return tree_;
    }
 
    // private data members
  private:
    const Tree& tree_;
    const Boxes& boxes;
};  // class Point_search
 
}  // namespace moab
 
#endif  // POINT_LOCATER_HPP