diff --git a/dumux/common/geometry/boundingboxtree.hh b/dumux/common/geometry/boundingboxtree.hh
index 1b02e16d6b861beb5f222de2aae3d54c61812a67..a1cd4cf8b5b222ecd5f84891d2bccedfa43707af 100644
--- a/dumux/common/geometry/boundingboxtree.hh
+++ b/dumux/common/geometry/boundingboxtree.hh
@@ -351,9 +351,9 @@ inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
{
using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
static constexpr ctype eps_ = 1.0e-7;
- const ctype eps0 = eps_*(b[3] - b[0]);
- const ctype eps1 = eps_*(b[4] - b[1]);
- const ctype eps2 = eps_*(b[5] - b[2]);
+ const ctype eps0 = eps_*std::max(b[3]-b[0], a[3]-a[0]);
+ const ctype eps1 = eps_*std::max(b[4]-b[1], a[4]-a[1]);
+ const ctype eps2 = eps_*std::max(b[5]-b[2], a[5]-a[2]);
return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 &&
b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 &&
b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2);
@@ -370,8 +370,8 @@ inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
{
using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
static constexpr ctype eps_ = 1.0e-7;
- const ctype eps0 = eps_*(b[2] - b[0]);
- const ctype eps1 = eps_*(b[3] - b[1]);
+ const ctype eps0 = eps_*std::max(b[2]-b[0], a[2]-a[0]);
+ const ctype eps1 = eps_*std::max(b[3]-b[1], a[3]-a[1]);
return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 &&
b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1);
}
@@ -386,7 +386,7 @@ inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
{
using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
static constexpr ctype eps_ = 1.0e-1;
- const ctype eps0 = eps_*(b[1] - b[0]);
+ const ctype eps0 = eps_*std::max(b[1]-b[0], a[1]-a[0]);
return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0;
}
diff --git a/test/common/boundingboxtree/test_bboxtree.cc b/test/common/boundingboxtree/test_bboxtree.cc
index 9da087f54995fb259544db481d034238d4fabfc8..cff08c9c82ee59db25a054c7c35ef78427fe683f 100644
--- a/test/common/boundingboxtree/test_bboxtree.cc
+++ b/test/common/boundingboxtree/test_bboxtree.cc
@@ -68,13 +68,26 @@ public:
template <class OtherEntitySet, class OtherGridView>
int intersectTree(const Dumux::BoundingBoxTree<OtherEntitySet>& otherTree,
const OtherGridView& otherGridView,
- std::size_t expectedIntersections)
+ std::size_t expectedUniqueIntersections,
+ bool checkTotalIntersections = false,
+ std::size_t expectedIntersections = 0)
{
Dune::Timer timer;
const auto intersections = intersectingEntities(*tree_, otherTree);
- std::cout << "Computed tree intersections in " << timer.elapsed() << std::endl;
+ std::cout << "Computed " << intersections.size() << " tree intersections in " << timer.elapsed() << std::endl;
timer.reset();
+ if (checkTotalIntersections)
+ {
+ if (intersections.size() != expectedIntersections)
+ {
+ std::cerr << "BoundingBoxTree intersection failed: Expected "
+ << expectedIntersections << " (total) and got "
+ << intersections.size() << "!" <<std::endl;
+ return 1;
+ }
+ }
+
std::vector<std::vector<std::vector<GlobalPosition>>> map;
map.resize(otherGridView.size(0));
std::vector<std::vector<GlobalPosition>> uniqueIntersections;
@@ -96,10 +109,10 @@ public:
std::cout << "Found " << uniqueIntersections.size() << " unique intersections "
<< "in " << timer.elapsed() << std::endl;
- if (uniqueIntersections.size() != expectedIntersections)
+ if (uniqueIntersections.size() != expectedUniqueIntersections)
{
std::cerr << "BoundingBoxTree intersection failed: Expected "
- << expectedIntersections << " and got "
+ << expectedUniqueIntersections << " (unique) and got "
<< uniqueIntersections.size() << "!" <<std::endl;
return 1;
}
@@ -168,37 +181,93 @@ int main (int argc, char *argv[]) try
using NetworkGrid = Dune::FoamGrid<1, dimworld>;
using NetworkGridView = NetworkGrid::LeafGridView;
+ using EntitySet = Dumux::GridViewGeometricEntitySet<NetworkGridView, 0>;
+ Dumux::BoundingBoxTree<EntitySet> networkTree;
- std::cout << std::endl
- << "Intersect with other bounding box tree:" << std::endl
- << "***************************************"
- << std::endl;
+ {
+ std::cout << std::endl
+ << "Intersect with other bounding box tree:" << std::endl
+ << "***************************************"
+ << std::endl;
- auto networkGrid = std::shared_ptr<NetworkGrid>(Dune::GmshReader<NetworkGrid>::read("network1d.msh", false, false));
+ // create a network grid from gmsh
+ auto networkGrid = std::shared_ptr<NetworkGrid>(Dune::GmshReader<NetworkGrid>::read("network1d.msh", false, false));
- // scaling
- for (const auto& vertex : vertices(networkGrid->leafGridView()))
- {
- auto newPos = vertex.geometry().corner(0);
- newPos *= scaling;
- networkGrid->setPosition(vertex, newPos);
+ // scaling
+ for (const auto& vertex : vertices(networkGrid->leafGridView()))
+ {
+ auto newPos = vertex.geometry().corner(0);
+ newPos *= scaling;
+ networkGrid->setPosition(vertex, newPos);
+ }
+
+ // Dune::VTKWriter<NetworkGridView> lowDimVtkWriter(networkGrid->leafGridView());
+ // lowDimVtkWriter.write("network", Dune::VTK::ascii);
+
+ std::cout << "Constructed 1d network grid with " << networkGrid->leafGridView().size(0) << " elements." << std::endl;
+
+ // build the bulk grid bounding box tree
+ returns.push_back(test.build(grid->leafGridView()));
+
+ // build the network grid bounding box tree
+ networkTree.build(std::make_shared<EntitySet>(networkGrid->leafGridView()));
+
+ // intersect the two bounding box trees
+ returns.push_back(test.intersectTree(networkTree, networkGrid->leafGridView(), 20));
}
+ {
+ std::cout << std::endl
+ << "Intersect with other bounding box tree (2):" << std::endl
+ << "*******************************************"
+ << std::endl;
- // Dune::VTKWriter<NetworkGridView> lowDimVtkWriter(networkGrid->leafGridView());
- // lowDimVtkWriter.write("network", Dune::VTK::ascii);
+ // construct a line network grid
+ const GlobalPosition lowerLeftNW({0.5, 0.5, 0.0});
+ const GlobalPosition upperRightNW({0.5, 0.5, 1.0});
- std::cout << "Constructed 1d network grid with " << networkGrid->leafGridView().size(0) << " elements." << std::endl;
+ // make the grid (structured interval grid in dimworld space)
+ Dune::GridFactory<NetworkGrid> factory;
- // build the bulk grid bounding box tree
- returns.push_back(test.build(grid->leafGridView()));
+ #if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+ constexpr auto geomType = Dune::GeometryTypes::line;
+ #else
+ auto geomType = Dune::GeometryType(1);
+ #endif
- // build the network grid bounding box tree
- using EntitySet = Dumux::GridViewGeometricEntitySet<NetworkGridView, 0>;
- Dumux::BoundingBoxTree<EntitySet> networkTree;
- networkTree.build(std::make_shared<EntitySet>(networkGrid->leafGridView()));
+ // create a step vector
+ auto step = upperRightNW;
+ step -= lowerLeftNW, step /= numCellsX;
+
+ // create the vertices
+ auto globalPos = lowerLeftNW;
+ for (unsigned int vIdx = 0; vIdx <= numCellsX; vIdx++, globalPos += step)
+ factory.insertVertex(globalPos);
+
+ // create the cells
+ for(unsigned int eIdx = 0; eIdx < numCellsX; eIdx++)
+ factory.insertElement(geomType, {eIdx, eIdx+1});
- // intersect the two bounding box trees
- returns.push_back(test.intersectTree(networkTree, networkGrid->leafGridView(), 20));
+ auto networkGrid = std::shared_ptr<NetworkGrid>(factory.createGrid());
+
+ // scaling
+ for (const auto& vertex : vertices(networkGrid->leafGridView()))
+ {
+ auto newPos = vertex.geometry().corner(0);
+ newPos *= scaling;
+ networkGrid->setPosition(vertex, newPos);
+ }
+
+ std::cout << "Constructed 1d network grid with " << networkGrid->leafGridView().size(0) << " elements." << std::endl;
+
+ // build the bulk grid bounding box tree
+ returns.push_back(test.build(grid->leafGridView()));
+
+ // build the network grid bounding box tree
+ networkTree.build(std::make_shared<EntitySet>(networkGrid->leafGridView()));
+
+ // intersect the two bounding box trees
+ returns.push_back(test.intersectTree(networkTree, networkGrid->leafGridView(), 10, true, 40));
+ }
}
#endif
}