GenericBoundingBoxTree.h

// Copyright (C) 2013 Anders Logg
//
// This file is part of DOLFIN.
//
// DOLFIN is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// DOLFIN is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
//
// First added:  2013-04-23
// Last changed: 2014-02-06

#ifndef __GENERIC_BOUNDING_BOX_TREE_H
#define __GENERIC_BOUNDING_BOX_TREE_H

#include <memory>
#include <sstream>
#include <set>
#include <vector>
#include <dolfin/geometry/Point.h>

namespace dolfin
{

  // Forward declarations
  class Mesh;
  class MeshEntity;


  class GenericBoundingBoxTree
  {
  public:

    GenericBoundingBoxTree();

    virtual ~GenericBoundingBoxTree() {}

    static std::shared_ptr<GenericBoundingBoxTree> create(unsigned int dim);

    void build(const Mesh& mesh, std::size_t tdim);

    void build(const std::vector<Point>& points);

    std::vector<unsigned int>
    compute_collisions(const Point& point) const;

    std::pair<std::vector<unsigned int>, std::vector<unsigned int> >
    compute_collisions(const GenericBoundingBoxTree& tree) const;

    std::vector<unsigned int>
    compute_entity_collisions(const Point& point,
                              const Mesh& mesh) const;

    std::vector<unsigned int>
    compute_process_collisions(const Point& point) const;

    std::pair<std::vector<unsigned int>, std::vector<unsigned int> >
    compute_entity_collisions(const GenericBoundingBoxTree& tree,
                              const Mesh& mesh_A, const Mesh& mesh_B) const;

    unsigned int compute_first_collision(const Point& point) const;

    unsigned int compute_first_entity_collision(const Point& point,
                                              const Mesh& mesh) const;

    std::pair<unsigned int, double> compute_closest_entity(const Point& point,
                                                           const Mesh& mesh) const;

    std::pair<unsigned int, double> compute_closest_point(const Point& point) const;

    std::string str(bool verbose=false);

  protected:

    struct BBox
    {
      unsigned int child_0;
      unsigned int child_1;
    };

    std::size_t _tdim;

    std::vector<BBox> _bboxes;

    std::vector<double> _bbox_coordinates;

    mutable std::shared_ptr<GenericBoundingBoxTree> _point_search_tree;

    std::shared_ptr<GenericBoundingBoxTree> _global_tree;

    void clear();

    //--- Recursive build functions ---

    unsigned int _build(const std::vector<double>& leaf_bboxes,
                        const std::vector<unsigned int>::iterator& begin,
                        const std::vector<unsigned int>::iterator& end,
                        std::size_t gdim);

    unsigned int _build(const std::vector<Point>& points,
                        const std::vector<unsigned int>::iterator& begin,
                        const std::vector<unsigned int>::iterator& end,
                        std::size_t gdim);

    //--- Recursive search functions ---

    // Note that these functions are made static for consistency as
    // some of them need to deal with more than tree.

    // Compute collisions with point (recursive)
    static void
    _compute_collisions(const GenericBoundingBoxTree& tree,
                        const Point& point,
                        unsigned int node,
                        std::vector<unsigned int>& entities,
                        const Mesh* mesh);

    // Compute collisions with tree (recursive)
    static void
    _compute_collisions(const GenericBoundingBoxTree& A,
                        const GenericBoundingBoxTree& B,
                        unsigned int node_A,
                        unsigned int node_B,
                        std::vector<unsigned int>& entities_A,
                        std::vector<unsigned int>& entities_B,
                        const Mesh* mesh_A,
                        const Mesh* mesh_B);

    // Compute first collision (recursive)
    static unsigned int
    _compute_first_collision(const GenericBoundingBoxTree& tree,
                             const Point& point,
                             unsigned int node);

    // Compute first entity collision (recursive)
    static unsigned int
    _compute_first_entity_collision(const GenericBoundingBoxTree& tree,
                                    const Point& point,
                                    unsigned int node,
                                    const Mesh& mesh);

    // Compute closest entity (recursive)
    static void _compute_closest_entity(const GenericBoundingBoxTree& tree,
                                        const Point& point,
                                        unsigned int node,
                                        const Mesh& mesh,
                                        unsigned int& closest_entity,
                                        double& R2);

    // Compute closest point (recursive)
    static void
    _compute_closest_point(const GenericBoundingBoxTree& tree,
                           const Point& point,
                           unsigned int node,
                           unsigned int& closest_point,
                           double& R2);

    //--- Utility functions ---

    void build_point_search_tree(const Mesh& mesh) const;

    void compute_bbox_of_entity(double* b,
                                const MeshEntity& entity,
                                std::size_t gdim) const;

    void sort_points(std::size_t axis,
                     const std::vector<Point>& points,
                     const std::vector<unsigned int>::iterator& begin,
                     const std::vector<unsigned int>::iterator& middle,
                     const std::vector<unsigned int>::iterator& end);

    inline unsigned int add_bbox(const BBox& bbox,
                                 const double* b,
                                 std::size_t gdim)
    {
      // Add bounding box
      _bboxes.push_back(bbox);

      // Add bounding box coordinates
      for (std::size_t i = 0; i < 2*gdim; ++i)
        _bbox_coordinates.push_back(b[i]);

      return _bboxes.size() - 1;
    }

    inline const BBox& get_bbox(unsigned int node) const
    {
      return _bboxes[node];
    }

    inline unsigned int num_bboxes() const
    {
      return _bboxes.size();
    }

    inline unsigned int add_point(const BBox& bbox,
                                  const Point& point,
                                  std::size_t gdim)
    {
      // Add bounding box
      _bboxes.push_back(bbox);

      // Add point coordinates (twice)
      const double* x = point.coordinates();
      for (std::size_t i = 0; i < gdim; ++i)
        _bbox_coordinates.push_back(x[i]);
      for (std::size_t i = 0; i < gdim; ++i)
        _bbox_coordinates.push_back(x[i]);

      return _bboxes.size() - 1;
    }

    inline bool is_leaf(const BBox& bbox, unsigned int node) const
    {
      // Leaf nodes are marked by setting child_0 equal to the node itself
      return bbox.child_0 == node;
    }

    struct less_x_point
    {
      const std::vector<Point>& points;

      less_x_point(const std::vector<Point>& points): points(points) {}

      inline bool operator()(unsigned int i, unsigned int j)
      {
        const double* pi = points[i].coordinates();
        const double* pj = points[j].coordinates();
        return pi[0] < pj[0];
      }
    };

    struct less_y_point
    {
      const std::vector<Point>& points;

      less_y_point(const std::vector<Point>& points): points(points) {}

      inline bool operator()(unsigned int i, unsigned int j)
      {
        const double* pi = points[i].coordinates();
        const double* pj = points[j].coordinates();
        return pi[1] < pj[1];
      }
    };

    struct less_z_point
    {
      const std::vector<Point>& points;

      less_z_point(const std::vector<Point>& points): points(points) {}

      inline bool operator()(unsigned int i, unsigned int j)
      {
        const double* pi = points[i].coordinates();
        const double* pj = points[j].coordinates();
        return pi[2] < pj[2];
      }
    };

    //--- Dimension-dependent functions to be implemented by subclass ---

    virtual std::size_t gdim() const = 0;

    virtual const double* get_bbox_coordinates(unsigned int node) const = 0;

    virtual bool
    point_in_bbox(const double* x, unsigned int node) const = 0;

    virtual bool
    bbox_in_bbox(const double* a, unsigned int node) const = 0;

    virtual double
    compute_squared_distance_bbox(const double* x, unsigned int node) const = 0;

    virtual double
    compute_squared_distance_point(const double* x, unsigned int node) const = 0;

    virtual void
    compute_bbox_of_bboxes(double* bbox,
                           std::size_t& axis,
                           const std::vector<double>& leaf_bboxes,
                           const std::vector<unsigned int>::iterator& begin,
                           const std::vector<unsigned int>::iterator& end) = 0;

    virtual void
    compute_bbox_of_points(double* bbox,
                           std::size_t& axis,
                           const std::vector<Point>& points,
                           const std::vector<unsigned int>::iterator& begin,
                           const std::vector<unsigned int>::iterator& end) = 0;

    virtual void
    sort_bboxes(std::size_t axis,
                const std::vector<double>& leaf_bboxes,
                const std::vector<unsigned int>::iterator& begin,
                const std::vector<unsigned int>::iterator& middle,
                const std::vector<unsigned int>::iterator& end) = 0;

    void tree_print(std::stringstream& s, unsigned int i);


  };

}

#endif