search_direct_product.h
Go to the documentation of this file.
1 #ifndef PORTAGE_SEARCH_SEARCH_DIRECT_PRODUCT_H_
2 #define PORTAGE_SEARCH_SEARCH_DIRECT_PRODUCT_H_
3 
4 #include <cassert>
5 #include <cstdint>
6 #include <limits>
7 #include <algorithm>
8 #include <array>
9 #include <utility>
10 #include <vector>
11 
12 #include "wonton/support/Point.h"
14 
27 namespace Portage {
28 
45 template <
46  typename ID_t, int D, typename SourceMeshType, typename TargetMeshType>
48 
49  public:
50 
51  // ==========================================================================
52  // Constructors and destructors
53 
55  SearchDirectProductBase() = delete;
56 
62  SearchDirectProductBase(const SourceMeshType& source_mesh,
63  const TargetMeshType& target_mesh);
64 
67  delete;
68 
69  // ==========================================================================
70  // Callable operation
71 
80  std::vector<ID_t> operator() (const int tgt_cell) const;
81 
82  private:
83 
84  // ==========================================================================
85  // Private support methods
86 
88  void fill_list_by_dim(
89  std::vector<ID_t> &list,
90  const int d, const int idx_start, const std::array<int,D> &step_size,
91  const std::array<int,D> &ilo, const std::array<int,D> &ihi,
92  std::array<int,D> &indices) const;
93 
95  std::vector<ID_t> list_cells(
96  const std::array<int,D>& ilo, const std::array<int,D>& ihi) const;
97 
98  // ==========================================================================
99  // Class data
100 
101  const SourceMeshType& sourceMesh_;
102  const TargetMeshType& targetMesh_;
103 
104 }; // class SearchDirectProductBase
105 
106 
107 // ============================================================================
108 // Typedef
109 
121 template<int D, typename SourceMeshType, typename TargetMeshType>
122 using SearchDirectProduct =
124 
125 
126 // ============================================================================
127 // Constructors and destructors
128 
129 // ____________________________________________________________________________
130 // Constructor with meshes.
131 template <
132  typename ID_t, int D, typename SourceMeshType, typename TargetMeshType>
135  const SourceMeshType& source_mesh, const TargetMeshType& target_mesh)
136  : sourceMesh_(source_mesh), targetMesh_(target_mesh) {
137 }
138 
139 // ============================================================================
140 // Callable
141 
142 // Search for source cells that intersect a given target cell.
143 template <
144  typename ID_t, int D, typename SourceMeshType, typename TargetMeshType>
145 std::vector<ID_t>
147  operator() ( const int tgt_cell) const {
148 
149  // Tolerance for floating-point round-off
150  const auto EPSILON = 10. * std::numeric_limits<double>::epsilon();
151 
152  // Get the bounding box for the target mesh cell
153  Wonton::Point<D> tlo, thi;
154  targetMesh_.cell_get_bounds(tgt_cell, &tlo, &thi);
155  for (int d = 0; d < D; ++d) {
156  // allow for roundoff error: if the intersection is within epsilon, assume
157  // this is a precision issue rather than true physical overlap
158  auto eps = EPSILON * (thi[d] - tlo[d]);
159  tlo[d] += eps;
160  thi[d] -= eps;
161  }
162 
163  // quick check: does this cell intersect the source mesh at all?
164  // if not, exit early
165  Wonton::Point<D> sglo, sghi;
166  sourceMesh_.get_global_bounds(&sglo, &sghi);
167  for (int d = 0; d < D; ++d) {
168  assert(tlo[d] < thi[d]);
169  assert(sglo[d] < sghi[d]);
170  if (tlo[d] >= sghi[d] || thi[d] <= sglo[d])
171  return std::vector<ID_t>();
172  }
173 
174  // find which source cells overlap target cell, in each dimension
175  std::array<int,D> ilo, ihi;
176  for (int d = 0; d < D; ++d) {
177  auto saxis_begin = sourceMesh_.axis_point_begin(d);
178  auto saxis_end = sourceMesh_.axis_point_end(d);
179 
180  // find last axis point less than or equal to tlo
181  auto is_point_above = [this,d] (const double x, const int pid) {
182  return (x < this->sourceMesh_.get_axis_point(d, pid));
183  };
184  auto itrlo = std::upper_bound(saxis_begin, saxis_end, tlo[d],
185  is_point_above);
186  ilo[d] = itrlo - saxis_begin - 1;
187  // std::max(ilo[d],0) instead of assert(ilo[d] >= 0) to account for cells
188  // that only partially overlap the target cell
189  ilo[d] = std::max(ilo[d], 0);
190 
191  // find first axis point greater than or equal to thi
192  auto is_point_below = [this,d] (const int pid, const double x) {
193  return (this->sourceMesh_.get_axis_point(d, pid) < x);
194  };
195  auto itrhi = std::lower_bound(itrlo, saxis_end, thi[d],
196  is_point_below);
197  ihi[d] = itrhi - saxis_begin;
198  ihi[d] = std::min(ihi[d], sourceMesh_.num_axis_cells(d));
199  assert(ihi[d] > ilo[d]);
200  } // for d
201 
202  // Generate list of cells from lower and upper bounds, return list
203  return list_cells(ilo, ihi);
204 
205 } // operator()
206 
207 
208 // ============================================================================
209 // Private support methods
210 
211 // Loop over each dimension recursively and fill the list
212 template<
213  typename ID_t, int D, typename SourceMeshType, typename TargetMeshType>
216  std::vector<ID_t> &list,
217  const int d, const int idx_start, const std::array<int,D> &step_size,
218  const std::array<int,D> &ilo, const std::array<int,D> &ihi,
219  std::array<int,D> &indices) const {
220  if (d >= 0) {
221  int index = idx_start;
222  for (indices[d] = ilo[d]; indices[d] < ihi[d]; ++indices[d]) {
223  fill_list_by_dim(list, d-1, index, step_size, ilo, ihi, indices);
224  index += step_size[d];
225  }
226  } else {
227  list[idx_start] = sourceMesh_.indices_to_cellid(indices);
228  }
229 }
230 
231 
232 // List cells given index bounds in each dimension
233 template <
234  typename ID_t, int D, typename SourceMeshType, typename TargetMeshType>
235 std::vector<ID_t>
238  const std::array<int,D>& ilo, const std::array<int,D>& ihi) const {
239 
240  // Compute step sizes for recursion
241  std::array<int,D> stepsize;
242  stepsize[0] = 1;
243  for(int d = 0; d < D-1; ++d) {
244  stepsize[d+1] = stepsize[d] * (ihi[d] - ilo[d]);
245  }
246 
247  // Declare the list of cells to be returned
248  int list_size = stepsize[D-1] * (ihi[D-1] - ilo[D-1]);
249  std::vector<ID_t> list(list_size);
250 
251  // Recurse across dimensions
252  std::array<int,D> indices;
253  fill_list_by_dim(list, D-1, 0, stepsize, ilo, ihi, indices);
254 
255  return list;
256 
257 } // list_cells
258 
259 } // namespace Portage
260 
261 #endif // PORTAGE_SEARCH_SEARCH_DIRECT_PRODUCT_H_
double eps
Definition: test_tangram_intersect_2D.cc:20
SearchDirectProductBase()=delete
Default constructor (disabled)
double const epsilon
Numerical tolerance.
Definition: weight.h:34
std::vector< ID_t > operator()(const int tgt_cell) const
Search for source cells that intersect a given target cell.
Definition: search_direct_product.h:147
SearchDirectProductBase & operator=(const SearchDirectProductBase &)=delete
Assignment operator (disabled)
Definition of the SearchDirectProductBase class.
Definition: search_direct_product.h:47
Definition: coredriver.h:42