From 53c7d3e9f34e57f371b9d9911132a4b487b3a0f1 Mon Sep 17 00:00:00 2001
From: Tim Jupe <tim.jupe@iws.uni-stuttgart.de>
Date: Wed, 1 Apr 2020 11:26:11 +0200
Subject: [PATCH] [test] Add new test for point intersections with a 3D
 geometry

---
 test/common/geometry/CMakeLists.txt           |   1 +
 .../common/geometry/test_0d3d_intersection.cc | 102 ++++++++++++++++++
 2 files changed, 103 insertions(+)
 create mode 100644 test/common/geometry/test_0d3d_intersection.cc

diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
index c2e2c5c6c7..2ea2a16f3f 100644
--- a/test/common/geometry/CMakeLists.txt
+++ b/test/common/geometry/CMakeLists.txt
@@ -1,4 +1,5 @@
 dumux_add_test(SOURCES test_0d1d_intersection.cc LABELS unit)
+dumux_add_test(SOURCES test_0d3d_intersection.cc LABELS unit)
 dumux_add_test(SOURCES test_1d1d_intersection.cc LABELS unit)
 dumux_add_test(SOURCES test_1d3d_intersection.cc LABELS unit)
 dumux_add_test(SOURCES test_1d2d_intersection.cc LABELS unit)
diff --git a/test/common/geometry/test_0d3d_intersection.cc b/test/common/geometry/test_0d3d_intersection.cc
new file mode 100644
index 0000000000..f75a161a42
--- /dev/null
+++ b/test/common/geometry/test_0d3d_intersection.cc
@@ -0,0 +1,102 @@
+#include <config.h>
+
+#include <iostream>
+#include <algorithm>
+
+#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/intersectspointgeometry.hh>
+
+#ifndef DOXYGEN
+template<int dimworld = 3, class Geometry>
+bool testIntersection(const Geometry& geo,
+                      const Dune::FieldVector<double, dimworld>& p,
+                      bool foundExpected = false)
+{
+    bool found = Dumux::intersectsPointGeometry(p, geo);
+    if (!found && foundExpected)
+        std::cerr << "Failed detecting intersection with " << p << std::endl;
+    else if (found && foundExpected)
+        std::cout << "Found intersection with " << p << std::endl;
+    else if (found && !foundExpected)
+        std::cerr << "Found false positive: intersection with " << p << std::endl;
+    else if (!found && !foundExpected)
+        std::cout << "No intersection with " << p << 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;
+    
+    // collect returns to determine exit code
+    std::vector<bool> returns;
+
+    for (auto scaling : {1.0, 1e3, 1e12, 1e-12})
+    {
+        std::cout << "Test with scaling " << scaling << std::endl;
+
+        using Points = std::vector<Dune::FieldVector<double, dimworld>>;
+        using Geo = Dune::MultiLinearGeometry<double, 3, dimworld>;
+    
+        // test tetrahedron-point intersections
+        std::cout << "test tetrahedron-point intersections" << std::endl;
+
+        auto cornersTetrahedron = Points({{0.0, 0.0, 0.0}, {1.0*scaling, 0.0, 0.0}, {0.0, 1.0*scaling, 0.0}, {0.0, 0.0, 1.0*scaling}});
+        auto tetrahedron = Geo(Dune::GeometryTypes::tetrahedron, cornersTetrahedron);
+
+        returns.push_back(testIntersection(tetrahedron, {0.0, 0.0, 0.0}, true));
+        returns.push_back(testIntersection(tetrahedron, {0.0, 0.0, 0.5*scaling}, true));
+        returns.push_back(testIntersection(tetrahedron, {0.25*scaling, 0.25*scaling, 0.5*scaling}, true));
+        returns.push_back(testIntersection(tetrahedron, {0.5*scaling, 0.5*scaling, 0.5*scaling}));
+        returns.push_back(testIntersection(tetrahedron, {1.01*scaling, 0.0, 0.0}));
+        returns.push_back(testIntersection(tetrahedron, {0.5*scaling, 0.0, 0.51*scaling}));
+
+        // test hexahedron-point intersections
+        std::cout << "test hexahedron-point intersections" << std::endl;
+
+        auto cornersHexahedron = Points({{0.0, 0.0, 0.0}, {1.0*scaling, 0.0, 0.0}, {0.0, 1.0*scaling, 0.0}, {1.0*scaling, 1.0*scaling, 0.0}, {0.0, 0.0, 1.0*scaling}, {1.0*scaling, 0.0, 1.0*scaling}, {0.0, 1.0*scaling, 1.0*scaling}, {1.0*scaling, 1.0*scaling, 1.0*scaling}});
+        auto hexahedron = Geo(Dune::GeometryTypes::hexahedron, cornersHexahedron);
+
+        returns.push_back(testIntersection(hexahedron, {0.0, 0.0, 0.0}, true));
+
+        // test pyramid-point intersections
+        std::cout << "test pyramid-point intersections" << std::endl;
+
+        auto cornersPyramid = Points({{0.0, 0.0, 0.0}, {1.0*scaling, 0.0, 0.0}, {0.0, 1.0*scaling, 0.0}, {1.0*scaling, 1.0*scaling, 0.0}, {0.5*scaling, 0.5*scaling, 1.0*scaling}});
+        auto pyramid = Geo(Dune::GeometryTypes::pyramid, cornersPyramid);
+
+        returns.push_back(testIntersection(pyramid, {0.0, 0.0, 0.0}, true));
+
+        // test prism-point intersections
+        std::cout << "test prism-point intersections" << std::endl;
+
+        auto cornersPrism = Points({{0.0, 0.0, 0.0}, {1.0*scaling, 0.0, 0.0}, {0.0, 1.0*scaling, 0.0}, {0.0, 0.0, 1.0*scaling}, {1.0*scaling, 0.0, 1.0*scaling}, {0.0, 1.0*scaling, 1.0*scaling}});
+        auto prism = Geo(Dune::GeometryTypes::prism, cornersPrism);
+
+        returns.push_back(testIntersection(prism, {0.0, 0.0, 0.0}, true));
+    }
+
+    // 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;
+}
-- 
GitLab