diff --git a/test/common/boundingboxtree/test_bboxtree.cc b/test/common/boundingboxtree/test_bboxtree.cc
index 007eca06256ed446fc5ff87b9b76cec08c265181..9c6066f7f27119c87195502ac28950a4035db2a6 100644
--- a/test/common/boundingboxtree/test_bboxtree.cc
+++ b/test/common/boundingboxtree/test_bboxtree.cc
@@ -83,8 +83,8 @@ public:
return 0;
}
- template <class OtherEntitySet, class OtherGridView>
- int intersectTree(const Dumux::BoundingBoxTree<OtherEntitySet>& otherTree,
+ template <class OtherGridView>
+ int intersectTree(const Dumux::BoundingBoxTree<GridViewGeometricEntitySet<OtherGridView>>& otherTree,
const OtherGridView& otherGridView,
std::size_t expectedUniqueIntersections,
bool checkTotalIntersections = false,
@@ -142,6 +142,37 @@ public:
return 0;
}
+ template <class GeometryType, class ExpectedIntersectingElements>
+ int intersectTree(const Dumux::BoundingBoxTree<GeometriesEntitySet<GeometryType>>& otherTree,
+ const ExpectedIntersectingElements& expectedIntersectingElements)
+ {
+ Dune::Timer timer;
+ const auto intersections = intersectingEntities(otherTree, *tree_);
+ std::cout << "Computed " << intersections.size() << " tree intersections in " << timer.elapsed() << std::endl;
+
+ std::unordered_map<std::size_t, std::vector<std::size_t>> sortedResults;
+
+ for (const auto& i : intersections)
+ {
+ sortedResults[i.first()].push_back(i.second());
+ std::sort(sortedResults[i.first()].begin(), sortedResults[i.first()].end());
+ sortedResults[i.first()].erase(std::unique(sortedResults[i.first()].begin(), sortedResults[i.first()].end()), sortedResults[i.first()].end());
+ }
+
+ for (int i = 0; i < expectedIntersectingElements.size(); ++i)
+ {
+ if (expectedIntersectingElements.at(i) != sortedResults[i])
+ {
+ std::cout << "Error for geometry type " << otherTree.entitySet().entity(i).geometry().type() << std::endl;
+ for (int j = 0; j < expectedIntersectingElements.at(i).size(); ++j)
+ if (expectedIntersectingElements.at(i)[j] != sortedResults[i][j])
+ std::cout << "expected index " << expectedIntersectingElements.at(i)[j] << ", got " << sortedResults[i][j] << std::endl;
+ return 1;
+ }
+ }
+ return 0;
+ }
+
private:
std::shared_ptr<BoundingBoxTree> tree_;
};
@@ -198,6 +229,97 @@ int main (int argc, char *argv[]) try
Dune::AffineGeometry<double, 2, WORLD_DIMENSION> geometry(Dune::GeometryTypes::simplex(2), corners);
returns.push_back(test.intersectGeometry(geometry, 2145)); // (33*33/2 - 33/2)*4 + 33
#endif
+
+ // test intersection of grid with 1D geometries (lines)
+ if (dimWorld > 1)
+ {
+ using GeometryType = Dune::MultiLinearGeometry<Scalar, 1, dimWorld>;
+ const std::vector<GlobalPosition> cornersDiagonal{lowerLeft, upperRight};
+ std::vector<GlobalPosition> cornersVertical{GlobalPosition(0.4*scaling), GlobalPosition(0.4*scaling)};
+ cornersVertical[1][dimWorld-1] = 0.6*scaling;
+ GeometryType diagonalLine(Dune::GeometryTypes::line, cornersDiagonal);
+ GeometryType verticalLine(Dune::GeometryTypes::line, cornersVertical);
+ std::vector<GeometryType> geometries{std::move(diagonalLine), std::move(verticalLine)};
+ using GeometriesEntitySet = Dumux::GeometriesEntitySet<GeometryType>;
+ GeometriesEntitySet entitySet(std::move(geometries));
+ Dumux::BoundingBoxTree<GeometriesEntitySet> geometriesTree(std::make_shared<GeometriesEntitySet>(entitySet));
+ std::unordered_map<std::size_t, std::vector<std::size_t>> expectedIntersectingElements;
+
+ if (dimWorld == 2)
+ {
+ expectedIntersectingElements[0] = {0, 34, 68, 102, 136, 170, 204, 238, 272, 306, 340, 374, 408,
+ 442, 476, 510, 544, 578, 612, 646, 680, 714, 748, 782, 816,
+ 850, 884, 918, 952, 986, 1020, 1054, 1088};
+ expectedIntersectingElements[1] = {442, 475, 508, 541, 574, 607, 640};
+ returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements));
+ }
+ else
+ {
+ expectedIntersectingElements[0] = {0, 1123, 2246, 3369, 4492, 5615, 6738, 7861, 8984, 10107, 11230,
+ 12353, 13476, 14599, 15722, 16845, 17968, 19091, 20214, 21337, 22460,
+ 23583, 24706, 25829, 26952, 28075, 29198, 30321, 31444, 32567,
+ 33690, 34813, 35936};
+ expectedIntersectingElements[1] = {14599 , 15688, 16777, 17866, 18955, 20044, 21133};
+ returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements));
+ }
+ }
+
+ // test intersection with 2D geometry (triangle)
+ if (dimWorld > 1)
+ {
+ using GeometryType = Dune::MultiLinearGeometry<Scalar, 2, dimWorld>;
+ const GlobalPosition corner0 = lowerLeft;
+ GlobalPosition corner1 = upperRight; corner1*= 0.1;
+ GlobalPosition corner2 = corner1; corner2[dimWorld-1] = lowerLeft[dimWorld-1];
+ std::vector<GlobalPosition> cornersTriangle{corner0, corner1, corner2};
+ const GeometryType triangle(Dune::GeometryTypes::simplex(2), cornersTriangle);
+ using GeometriesEntitySet = Dumux::GeometriesEntitySet<GeometryType>;
+ GeometriesEntitySet entitySet({triangle});
+ Dumux::BoundingBoxTree<GeometriesEntitySet> geometriesTree(std::make_shared<GeometriesEntitySet>(entitySet));
+ std::unordered_map<std::size_t, std::vector<std::size_t>> expectedIntersectingElements;
+
+ if (dimWorld == 2)
+ {
+ expectedIntersectingElements[0] = {0, 1, 2, 3, 34, 35, 36, 68, 69, 102};
+ returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements));
+ }
+ else
+ {
+ expectedIntersectingElements[0] = {0, 34, 68, 102, 1123, 1157, 1191, 2246, 2280, 3369};
+ returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements));
+ }
+ }
+
+ // test intersection with 2D geometry (quadrilateral)
+ if (dimWorld > 1)
+ {
+ using GeometryType = Dune::AxisAlignedCubeGeometry<Scalar, 2, dimWorld>;
+ GlobalPosition lowerLeftCube = upperRight; lowerLeftCube *= 0.4;
+ GlobalPosition upperRightCube = lowerLeftCube; upperRightCube[0] += 0.2*scaling; upperRightCube[1] += 0.2*scaling;
+ GeometryType cube(lowerLeftCube, upperRightCube);
+ using GeometriesEntitySet = Dumux::GeometriesEntitySet<GeometryType>;
+ GeometriesEntitySet entitySet(std::vector<GeometryType>{cube});
+ Dumux::BoundingBoxTree<GeometriesEntitySet> geometriesTree(std::make_shared<GeometriesEntitySet>(entitySet));
+ std::unordered_map<std::size_t, std::vector<std::size_t>> expectedIntersectingElements;
+
+ if (dimWorld == 2)
+ {
+ expectedIntersectingElements[0] = {442, 443, 444, 445, 446, 447, 448, 475, 476, 477, 478, 479, 480, 481,
+ 508, 509, 510, 511, 512, 513, 514, 541, 542, 543, 544, 545, 546, 547,
+ 574, 575, 576, 577, 578, 579, 580, 607, 608, 609, 610, 611, 612, 613,
+ 640, 641, 642, 643, 644, 645, 646};
+ returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements));
+ }
+ else
+ {
+ expectedIntersectingElements[0] = {14599, 14600, 14601, 14602, 14603, 14604, 14605, 14632, 14633, 14634,
+ 14635, 14636, 14637, 14638, 14665, 14666, 14667, 14668, 14669,14670,
+ 14671, 14698, 14699, 14700, 14701, 14702, 14703, 14704, 14731, 14732,
+ 14733, 14734, 14735, 14736, 14737, 14764, 14765, 14766, 14767,14768,
+ 14769, 14770, 14797, 14798, 14799, 14800, 14801, 14802, 14803};
+ returns.push_back(test.intersectTree(geometriesTree, expectedIntersectingElements));
+ }
+ }
}
#if HAVE_DUNE_FOAMGRID && WORLD_DIMENSION == 3