Merge branch 'feature/0d2d-intersection' into 'master'

Unit tests for point intersection with 2d geometries in  2d and 3d coorddim

Closes #856

See merge request !2131

test/common/geometry/CMakeLists.txt

 dumux_add_test(SOURCES test_0d1d_intersection.cc LABELS unit)
+dumux_add_test(SOURCES test_0d2d_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)

test/common/geometry/test_0d2d_intersection.cc

+#include <config.h>
+
+#include <iostream>
+#include <algorithm>
+#include <functional>
+#include <type_traits>
+
+#include <dune/common/classname.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/multilineargeometry.hh>
+
+#include <dumux/common/math.hh>
+
+#include "transformation.hh"
+#include "test_intersection.hh"
+
+namespace Dumux {
+
+template<int dimworld, class Transformation>
+void runIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
+{
+    using ctype = typename std::decay_t<decltype(transform({0.0}))>::value_type;
+    using Geo = Dune::MultiLinearGeometry<ctype, 2, dimworld>;
+    using Points2D = std::vector<Dune::FieldVector<ctype, 2>>;
+    using PointsDimWorld = std::vector<Dune::FieldVector<ctype, dimworld>>;
+
+    // test triangle-point intersections
+    if (verbose) std::cout << "\n  -- Test triangle-point intersections" << std::endl;
+
+    auto cornersTri2D = Points2D ({{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}});
+    auto cornersTri = PointsDimWorld(3);
+    std::transform(cornersTri2D.begin(), cornersTri2D.end(), cornersTri.begin(),
+                   [&](const auto& p) { return transform(p); });
+    auto triangle = Geo(Dune::GeometryTypes::triangle, cornersTri);
+
+    for (const auto& corner : cornersTri)
+        returns.push_back(testIntersection(triangle, corner, true, verbose));
+
+    returns.push_back(testIntersection(triangle, transform({0.25, 0.25}), true, verbose));
+    returns.push_back(testIntersection(triangle, transform({0.5, 0.5}), true, verbose));
+    returns.push_back(testIntersection(triangle, transform({0.0, 0.5}), true, verbose));
+    returns.push_back(testIntersection(triangle, transform({1.01, 0.0}), false, verbose));
+    returns.push_back(testIntersection(triangle, transform({0.5, 0.51}), false, verbose));
+    returns.push_back(testIntersection(triangle, transform({0.0, -0.01}), false, verbose));
+
+    // test quadrilateral-point intersections
+    if (verbose) std::cout << "\n  -- Test quadrilateral-point intersections" << std::endl;
+
+    auto cornersQuad2D = Points2D ({{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}, {1.0, 1.0}});
+    auto cornersQuad = PointsDimWorld(4);
+    std::transform(cornersQuad2D.begin(), cornersQuad2D.end(), cornersQuad.begin(),
+                   [&](const auto& p) { return transform(p); });
+    auto quadrilateral = Geo(Dune::GeometryTypes::quadrilateral, cornersQuad);
+
+    for (const auto& corner : cornersQuad)
+        returns.push_back(testIntersection(quadrilateral, corner, true, verbose));
+
+    returns.push_back(testIntersection(quadrilateral, transform({0.5, 0.5}), true, verbose));
+    returns.push_back(testIntersection(quadrilateral, transform({0.5, 0.0}), true, verbose));
+    returns.push_back(testIntersection(quadrilateral, transform({0.5, 1.0}), true, verbose));
+    returns.push_back(testIntersection(quadrilateral, transform({1.01, 1.0}), false, verbose));
+    returns.push_back(testIntersection(quadrilateral, transform({0.5, 1.01}), false, verbose));
+    returns.push_back(testIntersection(quadrilateral, transform({0.0, -0.01}), false, verbose));
+
+}
+
+template<class Transformation>
+void run3DIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
+{
+    // augment 2d vector by a third coordinate
+    using Point = std::decay_t<decltype(transform({0.0, 0.0, 0.0}))>;
+    const auto transform2D3D = [&](Dune::FieldVector<typename Point::value_type, 2> p2D){
+        return transform({p2D[0], p2D[1], 2.0});
+    };
+
+    runIntersectionTest<3>(returns, transform2D3D, verbose);
+}
+
+template<class Transformation>
+void run2DIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
+{ runIntersectionTest<2>(returns, transform, verbose); }
+
+} // end namespace Dumux
+
+int main (int argc, char *argv[]) try
+{
+    using namespace Dumux;
+
+    // collect returns to determine exit code
+    std::vector<bool> returns;
+    constexpr bool verbose = false;
+
+    {
+        using Vec = Dune::FieldVector<double, 2>;
+
+        std::cout << "Test for point intersection with 2d geometries in 2d corddim" << std::endl;
+        for (const double scaling : {1.0, 1e3, 1e12, 1e-12})
+            for (const auto& translation : {0.0, 1.0})
+                for (const double angle : {0.0, 0.2*M_PI, 0.5*M_PI, 0.567576567*M_PI, M_PI})
+                    run2DIntersectionTest(returns,
+                        make2DTransformation<double>(scaling, Vec(translation), angle, true), verbose);
+    }
+
+    {
+        using Vec = Dune::FieldVector<double, 3>;
+
+        std::cout << "Test for point intersection with 2d geometries in 3d corddim" << std::endl;
+        for (const double scaling : {1.0, 1e3, 1e12, 1e-12})
+            for (const auto& translation : {0.0, 1.0})
+                for (const double angle : {0.0, 0.2*M_PI, 0.5*M_PI, 0.567576567*M_PI, M_PI})
+                    for (const auto& rotAxis : {Vec(std::sqrt(3.0)/3.0), Vec({std::sqrt(2.0)/2.0, std::sqrt(2.0)/2.0, 0.0})})
+                        run3DIntersectionTest(returns,
+                            make3DTransformation<double>(scaling, Vec(translation), rotAxis, angle, true), verbose);
+    }
+
+    // determine the exit code
+    if (std::any_of(returns.begin(), returns.end(), std::logical_not<bool>{}))
+        return 1;
+
+    std::cout << "\n++++++++++++++++++++++\n"
+              << "All tests passed!"
+              << "\n++++++++++++++++++++++" << std::endl;
+
+    return 0;
+}
+// //////////////////////////////////
+//   Error handler
+// /////////////////////////////////
+catch (const Dune::Exception& e) {
+    std::cout << e << std::endl;
+    return 1;
+}

test/common/geometry/test_0d3d_intersection.cc

@@ -12,38 +12,11 @@
 #include <dune/geometry/multilineargeometry.hh>
 
 #include <dumux/common/math.hh>
-#include <dumux/common/geometry/intersectspointgeometry.hh>
 
-#ifndef DOXYGEN
-template<class Geometry>
-bool testIntersection(const Geometry& geo,
-                      const typename Geometry::GlobalCoordinate& p,
-                      bool foundExpected, bool verbose)
-{
-    bool found = Dumux::intersectsPointGeometry(p, geo);
-    if (!found && foundExpected)
-    {
-        std::cerr << "  Failed detecting intersection of " << geo.type();
-        for (int i = 0; i < geo.corners(); ++i)
-            std::cerr << " (" << geo.corner(i) << ")";
-        std::cerr << " with point: " << p << std::endl;
-    }
-    else if (found && !foundExpected)
-    {
-        std::cerr << "  Found false positive: intersection of " << geo.type();
-        for (int i = 0; i < geo.corners(); ++i)
-            std::cerr << " (" << geo.corner(i) << ")";
-        std::cerr << " with point: " << p << std::endl;
-    }
-    if (verbose)
-    {
-        if (found && foundExpected)
-            std::cout << "  Found intersection with " << p << std::endl;
-        else if (!found && !foundExpected)
-            std::cout << "  No intersection with " << p << std::endl;
-    }
-    return (found == foundExpected);
-}
+#include "transformation.hh"
+#include "test_intersection.hh"
+
+namespace Dumux {
 
 template<class Transformation>
 void runIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
@@ -57,7 +30,7 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
 
     auto cornersTet = Points({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}});
     std::transform(cornersTet.begin(), cornersTet.end(), cornersTet.begin(),
-                   [&transform](const auto& p) { return transform(p); });
+                   [&](const auto& p) { return transform(p); });
     const auto tetrahedron = Geo(Dune::GeometryTypes::tetrahedron, cornersTet);
 
     for (const auto& corner : cornersTet)
@@ -74,7 +47,7 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
     auto cornersHex = Points({{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}});
     std::transform(cornersHex.begin(), cornersHex.end(), cornersHex.begin(),
-                   [&transform](const auto& p) { return transform(p); });
+                   [&](const auto& p) { return transform(p); });
     auto hexahedron = Geo(Dune::GeometryTypes::hexahedron, cornersHex);
 
     for (const auto& corner : cornersHex)
@@ -91,7 +64,7 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
     auto cornersPyramid = Points({{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.5, 0.5, 1.0}});
     std::transform(cornersPyramid.begin(), cornersPyramid.end(), cornersPyramid.begin(),
-                   [&transform](const auto& p) { return transform(p); });
+                   [&](const auto& p) { return transform(p); });
     auto pyramid = Geo(Dune::GeometryTypes::pyramid, cornersPyramid);
 
     for (const auto& corner : cornersPyramid)
@@ -111,7 +84,7 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
     auto cornersPrism = Points({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0},
                                 {0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}, {0.0, 1.0, 1.0}});
     std::transform(cornersPrism.begin(), cornersPrism.end(), cornersPrism.begin(),
-                   [&transform](const auto& p) { return transform(p); });
+                   [&](const auto& p) { return transform(p); });
     auto prism = Geo(Dune::GeometryTypes::prism, cornersPrism);
 
     for (const auto& corner : cornersPrism)
@@ -124,34 +97,12 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
     returns.push_back(testIntersection(prism, transform({0.25, 0.25, 1.0001}), false, verbose));
 }
 
-template<class ctype>
-auto createTransformation(const ctype scale,
-                          const Dune::FieldVector<ctype, 3>& translate,
-                          const Dune::FieldVector<ctype, 3>& rotationAxis,
-                          const ctype rotationAngle)
-{
-    std::cout << "Intersection test with transformation:"
-              << " ctype: " << Dune::className<ctype>()
-              << ", scaling: " << scale
-              << ", translation: " << translate
-              << ", rotationAxis: " << rotationAxis
-              << ", rotationAngle: " << rotationAngle << std::endl;
-    const ctype sinAngle = std::sin(rotationAngle);
-    const ctype cosAngle = std::cos(rotationAngle);
-    return [=](Dune::FieldVector<ctype, 3> p){
-        p *= scale;
-        p.axpy(scale, translate);
-        auto tp = p;
-        tp *= cosAngle;
-        tp.axpy(sinAngle, Dumux::crossProduct(rotationAxis, p));
-        return tp.axpy((1.0-cosAngle)*(rotationAxis*p), rotationAxis);
-    };
-}
-
-#endif
+} // end namespace Dumux
 
 int main (int argc, char *argv[]) try
 {
+    using namespace Dumux;
+
     // collect returns to determine exit code
     std::vector<bool> returns;
     constexpr bool verbose = false;
@@ -162,7 +113,7 @@ int main (int argc, char *argv[]) try
             for (const double angle : {0.0, 0.2*M_PI, 0.5*M_PI, 0.567576567*M_PI, M_PI})
                 for (const auto& rotAxis : {Vec(std::sqrt(3.0)/3.0), Vec({std::sqrt(2.0)/2.0, std::sqrt(2.0)/2.0, 0.0})})
                     runIntersectionTest(returns,
-                        createTransformation<double>(scaling, Vec(translation), rotAxis, angle), verbose);
+                        make3DTransformation<double>(scaling, Vec(translation), rotAxis, angle, true), verbose);
 
     // determine the exit code
     if (std::any_of(returns.begin(), returns.end(), std::logical_not<bool>{}))

test/common/geometry/test_intersection.hh

+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief Intersection test helpers
+ */
+#ifndef DUMUX_TEST_INTERSECTION_HH
+#define DUMUX_TEST_INTERSECTION_HH
+
+#include <dumux/common/geometry/intersectspointgeometry.hh>
+
+namespace Dumux {
+
+/*!
+ * \file
+ * \ingroup Common
+ * \brief Test intersection overload for point and geometry
+ */
+template<class Geometry>
+bool testIntersection(const Geometry& geo,
+                      const typename Geometry::GlobalCoordinate& p,
+                      bool foundExpected, bool verbose)
+{
+    const bool found = intersectsPointGeometry(p, geo);
+    if (!found && foundExpected)
+    {
+        std::cerr << "  Failed detecting intersection of " << geo.type();
+        for (int i = 0; i < geo.corners(); ++i)
+            std::cerr << " (" << geo.corner(i) << ")";
+        std::cerr << " with point: " << p << std::endl;
+    }
+    else if (found && !foundExpected)
+    {
+        std::cerr << "  Found false positive: intersection of " << geo.type();
+        for (int i = 0; i < geo.corners(); ++i)
+            std::cerr << " (" << geo.corner(i) << ")";
+        std::cerr << " with point: " << p << std::endl;
+    }
+    if (verbose)
+    {
+        if (found && foundExpected)
+            std::cout << "  Found intersection with " << p << std::endl;
+        else if (!found && !foundExpected)
+            std::cout << "  No intersection with " << p << std::endl;
+    }
+    return (found == foundExpected);
+}
+
+} // end namespace Dumux
+
+#endif

test/common/geometry/transformation.hh

+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 3 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief Create a transformation function (translation + rotation)
+ */
+#ifndef DUMUX_TEST_GEOMETRY_TRANSFORMATION_HH
+#define DUMUX_TEST_GEOMETRY_TRANSFORMATION_HH
+
+#include <iostream>
+
+#include <dune/common/classname.hh>
+#include <dune/common/fvector.hh>
+
+#include <dumux/common/math.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Common
+ * \brief Create a transformation function (translation) in 1D
+ */
+template<class ctype>
+auto make1DTransformation(const ctype scale,
+                          const Dune::FieldVector<ctype, 1>& translate,
+                          bool verbose = true)
+{
+    if (verbose)
+        std::cout << "Created 1D transformation with"
+                  << " ctype: " << Dune::className<ctype>()
+                  << ", scaling: " << scale
+                  << ", translation: " << translate << std::endl;
+
+    return [=](Dune::FieldVector<ctype, 1> p){
+        p *= scale;
+        p.axpy(scale, translate);
+        auto tp = p;
+        return p;
+    };
+}
+
+/*!
+ * \ingroup Common
+ * \brief Create a transformation function (translation + rotation) in 2D
+ * \note https://en.wikipedia.org/wiki/Rotation_matrix
+ */
+template<class ctype>
+auto make2DTransformation(const ctype scale,
+                          const Dune::FieldVector<ctype, 2>& translate,
+                          const ctype rotationAngle,
+                          bool verbose = true)
+{
+    if (verbose)
+        std::cout << "Created 2D transformation with"
+                  << " ctype: " << Dune::className<ctype>()
+                  << ", scaling: " << scale
+                  << ", translation: " << translate
+                  << ", rotationAngle: " << rotationAngle << std::endl;
+
+    using std::sin; using std::cos;
+    const ctype sinAngle = sin(rotationAngle);
+    const ctype cosAngle = cos(rotationAngle);
+    return [=](Dune::FieldVector<ctype, 2> p){
+        p *= scale;
+        p.axpy(scale, translate);
+        auto tp = p;
+        tp[0] = p[0]*cosAngle-p[1]*sinAngle;
+        tp[1] = p[0]*sinAngle+p[1]*cosAngle;
+        return tp;
+    };
+}
+
+/*!
+ * \ingroup Common
+ * \brief Create a transformation function (translation + rotation) in 3D
+ * \note https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
+ */
+template<class ctype>
+auto make3DTransformation(const ctype scale,
+                          const Dune::FieldVector<ctype, 3>& translate,
+                          const Dune::FieldVector<ctype, 3>& rotationAxis,
+                          const ctype rotationAngle,
+                          bool verbose = true)
+{
+    if (verbose)
+        std::cout << "Created transformation with"
+                  << " ctype: " << Dune::className<ctype>()
+                  << ", scaling: " << scale
+                  << ", translation: " << translate
+                  << ", rotationAxis: " << rotationAxis
+                  << ", rotationAngle: " << rotationAngle << std::endl;
+
+    using std::sin; using std::cos;
+    const ctype sinAngle = sin(rotationAngle);
+    const ctype cosAngle = cos(rotationAngle);
+    return [=](Dune::FieldVector<ctype, 3> p){
+        p *= scale;
+        p.axpy(scale, translate);
+        auto tp = p;
+        tp *= cosAngle;
+        tp.axpy(sinAngle, Dumux::crossProduct({rotationAxis}, p));
+        return tp.axpy((1.0-cosAngle)*(rotationAxis*p), rotationAxis);
+    };
+}
+
+} // end namespace Dumux
+
+#endif