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/CMakeLists.txt b/test/common/CMakeLists.txt
index 8c8c1ee2e273e88e4d5050dac2158f6b83de742a..ea0bdb19633e37b326a1675139f36d13c0f6afd4 100644
--- a/test/common/CMakeLists.txt
+++ b/test/common/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory(boundingboxtree)
+add_subdirectory(geometry)
add_subdirectory(math)
add_subdirectory(parameters)
add_subdirectory(propertysystem)
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
}
diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..41465b6623d54694b5d64f7d47e63283cd2161ec
--- /dev/null
+++ b/test/common/geometry/CMakeLists.txt
@@ -0,0 +1,8 @@
+dune_add_test(SOURCES test_1d3d_intersection.cc)
+
+#install sources
+install(FILES
+test_1d3d_intersection.cc
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/common/geometry)
+
+set(CMAKE_BUILD_TYPE Release)
diff --git a/test/common/geometry/test_1d3d_intersection.cc b/test/common/geometry/test_1d3d_intersection.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8f7cf4d6ca249a33ce46dec3b612002d26056cfa
--- /dev/null
+++ b/test/common/geometry/test_1d3d_intersection.cc
@@ -0,0 +1,118 @@
+#include <config.h>
+
+#include <iostream>
+#include <algorithm>
+#include <initializer_list>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/fvector.hh>
+#include <dune/geometry/multilineargeometry.hh>
+
+#include <dumux/common/geometry/geometryintersection.hh>
+
+
+#ifndef DOXYGEN
+template<int dimworld = 3>
+Dune::MultiLinearGeometry<double, 1, dimworld>
+makeLine(std::initializer_list<Dune::FieldVector<double, dimworld>>&& c)
+{
+ return {Dune::GeometryTypes::line, c};
+}
+
+template<int dimworld = 3>
+bool testIntersection(const Dune::MultiLinearGeometry<double, dimworld, dimworld>& cube,
+ const Dune::MultiLinearGeometry<double, 1, dimworld>& line,
+ bool foundExpected = true)
+{
+ using Test = Dumux::GeometryIntersection<Dune::MultiLinearGeometry<double,dimworld,dimworld>,
+ Dune::MultiLinearGeometry<double,1,dimworld> >;
+ typename Test::IntersectionType intersection;
+ bool found = Test::intersection(cube, line, intersection);
+ if (!found && foundExpected)
+ std::cerr << "Failed detecting intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
+ else if (found && foundExpected)
+ std::cout << "Found intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
+ else if (found && !foundExpected)
+ std::cerr << "Found false positive: intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
+ else if (!found && !foundExpected)
+ std::cout << "No intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
+ return (found == foundExpected);
+}
+#endif
+
+int main (int argc, char *argv[]) try
+{
+ // maybe initialize mpi
+ Dune::MPIHelper::instance(argc, argv);
+
+ constexpr int dimworld = 3;
+ constexpr int dim = 3;
+
+ std::vector<Dune::FieldVector<double, dimworld>> cubeCorners({
+ {0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {1.0, 1.0, 0.0},
+ {0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}, {0.0, 1.0, 1.0}, {1.0, 1.0, 1.0}
+ });
+
+ Dune::MultiLinearGeometry<double, dim, dimworld>
+ cube(Dune::GeometryTypes::cube(dimworld), cubeCorners);
+
+
+
+ // collect returns to determine exit code
+ std::vector<bool> returns;
+
+ // the tests
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}})));
+
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 0.0}, {1.0, 0.0, 1.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 0.0}, {1.0, 1.0, 0.0}})));
+
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 0.0}, {0.0, 0.0, 0.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 0.0}, {0.0, 1.0, 1.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 0.0}, {1.0, 1.0, 0.0}})));
+
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 0.0}, {1.0, 1.0, 1.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 0.0}, {0.0, 1.0, 0.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 0.0}, {1.0, 0.0, 0.0}})));
+
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 1.0}, {0.0, 1.0, 1.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 1.0}, {0.0, 0.0, 0.0}})));
+
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 1.0}, {0.0, 0.0, 1.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 1.0}, {1.0, 0.0, 0.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 1.0}, {1.0, 1.0, 1.0}})));
+
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 1.0}, {0.0, 0.0, 1.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 1.0}, {0.0, 1.0, 0.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 1.0}, {1.0, 1.0, 1.0}})));
+
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 1.0}, {1.0, 1.0, 0.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 1.0}, {0.0, 1.0, 1.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 1.0}, {1.0, 0.0, 1.0}})));
+
+ returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}})));
+ returns.push_back(testIntersection(cube, makeLine({{0.5, 0.5, 0.5}, {0.5, 0.5, -2.0}})));
+
+ returns.push_back(testIntersection(cube, makeLine({{0.5, 0.5, 0.0}, {0.5, 0.5, -2.0}}), false));
+ returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 1.0}, {2.0, 2.0, 2.0}}), false));
+
+ // determine the exit code
+ if (std::any_of(returns.begin(), returns.end(), [](bool i){ return !i; }))
+ return 1;
+
+ std::cout << "All tests passed!" << std::endl;
+
+ return 0;
+}
+// //////////////////////////////////
+// Error handler
+// /////////////////////////////////
+catch (const Dune::Exception& e) {
+ std::cout << e << std::endl;
+ return 1;
+}