1 #ifndef PORTAGE_SEARCH_SEARCH_DIRECT_PRODUCT_H_ 2 #define PORTAGE_SEARCH_SEARCH_DIRECT_PRODUCT_H_ 12 #include "wonton/support/Point.h" 46 typename ID_t,
int D,
typename SourceMeshType,
typename TargetMeshType>
63 const TargetMeshType& target_mesh);
80 std::vector<ID_t>
operator() (
const int tgt_cell)
const;
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;
95 std::vector<ID_t> list_cells(
96 const std::array<int,D>& ilo,
const std::array<int,D>& ihi)
const;
101 const SourceMeshType& sourceMesh_;
102 const TargetMeshType& targetMesh_;
121 template<
int D,
typename SourceMeshType,
typename TargetMeshType>
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) {
144 typename ID_t,
int D,
typename SourceMeshType,
typename TargetMeshType>
153 Wonton::Point<D> tlo, thi;
154 targetMesh_.cell_get_bounds(tgt_cell, &tlo, &thi);
155 for (
int d = 0; d < D; ++d) {
158 auto eps = EPSILON * (thi[d] - tlo[d]);
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>();
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);
181 auto is_point_above = [
this,d] (
const double x,
const int pid) {
182 return (x < this->sourceMesh_.get_axis_point(d, pid));
184 auto itrlo = std::upper_bound(saxis_begin, saxis_end, tlo[d],
186 ilo[d] = itrlo - saxis_begin - 1;
189 ilo[d] = std::max(ilo[d], 0);
192 auto is_point_below = [
this,d] (
const int pid,
const double x) {
193 return (this->sourceMesh_.get_axis_point(d, pid) < x);
195 auto itrhi = std::lower_bound(itrlo, saxis_end, thi[d],
197 ihi[d] = itrhi - saxis_begin;
198 ihi[d] = std::min(ihi[d], sourceMesh_.num_axis_cells(d));
199 assert(ihi[d] > ilo[d]);
203 return list_cells(ilo, ihi);
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 {
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];
227 list[idx_start] = sourceMesh_.indices_to_cellid(indices);
234 typename ID_t,
int D,
typename SourceMeshType,
typename TargetMeshType>
238 const std::array<int,D>& ilo,
const std::array<int,D>& ihi)
const {
241 std::array<int,D> stepsize;
243 for(
int d = 0; d < D-1; ++d) {
244 stepsize[d+1] = stepsize[d] * (ihi[d] - ilo[d]);
248 int list_size = stepsize[D-1] * (ihi[D-1] - ilo[D-1]);
249 std::vector<ID_t> list(list_size);
252 std::array<int,D> indices;
253 fill_list_by_dim(list, D-1, 0, stepsize, ilo, ihi, indices);
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