search_kdtree.h

Go to the documentation of this file.

/*
This file is part of the Ristra portage project.
Please see the license file at the root of this repository, or at:
    https://github.com/laristra/portage/blob/master/LICENSE
*/

#ifndef PORTAGE_SEARCH_SEARCH_KDTREE_H_
#define PORTAGE_SEARCH_SEARCH_KDTREE_H_

#include <vector>
#include <memory>

// portage includes
#include "portage/support/portage.h"
#include "portage/search/BoundBox.h"
#include "portage/search/kdtree.h"
#include "wonton/support/Point.h"

namespace Portage {

template <int D, Entity_kind on_what,
          typename SourceMeshType, typename TargetMeshType>
class SearchKDTree {
 public:

  SearchKDTree() = delete;

  SearchKDTree(const SourceMeshType & source_mesh,
               const TargetMeshType & target_mesh)
      : sourceMesh_(source_mesh), targetMesh_(target_mesh)  {}

  std::vector<int> operator() (const int entityId) const {
    std::vector<int> candidates;
    std::cerr << "Search not implemented for generic entity kind" << std::endl;
    return candidates;
  }

 private:
  const SourceMeshType & sourceMesh_;
  const TargetMeshType & targetMesh_;
  std::shared_ptr<Portage::KDTree<D>> tree_;
};  // class SearchKDTree





template <int D, typename SourceMeshType, typename TargetMeshType>
class SearchKDTree<D, Entity_kind::CELL, SourceMeshType, TargetMeshType> {
 public:

  SearchKDTree() = delete;

  SearchKDTree(const SourceMeshType & source_mesh,
               const TargetMeshType & target_mesh)
      : sourceMesh_(source_mesh), targetMesh_(target_mesh)  {

    const int numCells = sourceMesh_.num_owned_cells();
    std::vector<Portage::IsotheticBBox<D>> bboxes;
    bboxes.reserve(numCells);

    // find bounding boxes for all cells
    for (int c = 0; c < numCells; ++c) {
      std::vector<Wonton::Point<D>> cell_coord;
      sourceMesh_.cell_get_coordinates(c, &cell_coord);
      Portage::IsotheticBBox<D> bb;
      for (const auto& cc : cell_coord) {
        bb.add(cc);
      }
      bboxes.emplace_back(bb);
    }

    // create the k-d tree
    tree_ = std::make_shared<Portage::KDTree<D>>(*Portage::KDTreeCreate(bboxes));

  }  // SearchKDTree::SearchKDTree

  std::vector<int> operator() (const int cellId) const {
    // find bounding box for target cell
    std::vector<Wonton::Point<D>> cell_coord;
    targetMesh_.cell_get_coordinates(cellId, &cell_coord);
    Portage::IsotheticBBox<D> bb;
    for (const auto& cc : cell_coord)
      bb.add(cc);

    // now see which sourceMesh cells have bounding boxes overlapping
    // with target cell, using the kdtree - since Portage::Intersect does
    // not take a shared_ptr, we have have to dereference and take
    // address of tree_
    std::vector<int> lcandidates;
    Portage::Intersect(bb, &(*tree_), lcandidates);

    std::vector<int> candidates(lcandidates.begin(), lcandidates.end());
    return candidates;
  }  // SearchKDTree::operator()

 private:
  const SourceMeshType & sourceMesh_;
  const TargetMeshType & targetMesh_;
  std::shared_ptr<Portage::KDTree<D>> tree_;
};  // class SearchKDTree (CELL specialization)





template <int D, typename SourceMeshType, typename TargetMeshType>
class SearchKDTree<D, Entity_kind::NODE, SourceMeshType, TargetMeshType> {
 public:

  SearchKDTree() = delete;

  SearchKDTree(const SourceMeshType & source_mesh,
               const TargetMeshType & target_mesh)
      : sourceMesh_(source_mesh), targetMesh_(target_mesh)  {

    const int numNodes = sourceMesh_.num_owned_nodes();
    std::vector<Portage::IsotheticBBox<D>> bboxes;
    bboxes.reserve(numNodes);

    // find bounding boxes for all cells
    for (int n = 0; n < numNodes; ++n) {
      std::vector<Wonton::Point<D>> dual_cell_coord;
      sourceMesh_.dual_cell_get_coordinates(n, &dual_cell_coord);
      Portage::IsotheticBBox<D> bb;
      for (const auto& cc : dual_cell_coord)
        bb.add(cc);
      bboxes.emplace_back(bb);
    }

    // create the k-d tree
    tree_ = std::make_shared<Portage::KDTree<D>>(*Portage::KDTreeCreate(bboxes));

  }  // SearchKDTree::SearchKDTree

  //  ~SearchKDTree() { if (tree_) delete tree_; }

  std::vector<int> operator() (const int nodeId) const {
    // find bounding box for dual cell of target node
    std::vector<Wonton::Point<D>> dual_cell_coord;
    targetMesh_.dual_cell_get_coordinates(nodeId, &dual_cell_coord);
    Portage::IsotheticBBox<D> bb;
    for (const auto& cc : dual_cell_coord)
      bb.add(cc);

    // now see which sourceMesh dual cells have bounding boxes
    // overlapping with dual cell of targetMesh, using the kdtree -
    // since Portage::Intersect does not take a shared_ptr, we have have to
    // dereference and take address of tree_
    std::vector<int> lcandidates;
    Portage::Intersect(bb, &(*tree_), lcandidates);

    std::vector<int> candidates(lcandidates.begin(), lcandidates.end());
    return candidates;
  }  // SearchKDTree::operator()

 private:
  const SourceMeshType & sourceMesh_;
  const TargetMeshType & targetMesh_;
  std::shared_ptr<Portage::KDTree<D>> tree_;
};  // class SearchKDTree (NODE specialization)

}  // namespace Portage

#endif  // PORTAGE_SEARCH_SEARCH_KDTREE_H_