[geometryisection] introduce policy classes for intersections

dumux/common/geometry/intersectingentities.hh

@@ -24,12 +24,15 @@
 
 #include <cmath>
 #include <type_traits>
+
 #include <dune/common/fvector.hh>
+
 #include <dumux/common/math.hh>
 #include <dumux/common/geometry/boundingboxtree.hh>
 #include <dumux/common/geometry/intersectspointgeometry.hh>
 #include <dumux/common/geometry/geometryintersection.hh>
 #include <dumux/common/geometry/boundingboxtreeintersection.hh>
+#include <dumux/common/geometry/triangulation.hh>
 
 namespace Dumux {
 
@@ -149,13 +152,25 @@ void intersectingEntities(const BoundingBoxTree<EntitySet0>& treeA,
         const auto geometryA = treeA.entitySet().entity(eIdxA).geometry();
         const auto geometryB = treeB.entitySet().entity(eIdxB).geometry();
 
-        using IntersectionAlgorithm = GeometryIntersection< std::decay_t<decltype(geometryA)>,
-                                                            std::decay_t<decltype(geometryB)> >;
-        typename IntersectionAlgorithm::IntersectionType intersection;
+        using GeometryA = std::decay_t<decltype(geometryA)>;
+        using GeometryB = std::decay_t<decltype(geometryB)>;
+        using Policy = IntersectionPolicy::DefaultPolicy<GeometryA, GeometryB>;
+        using IntersectionAlgorithm = GeometryIntersection<GeometryA, GeometryB, Policy>;
+        using Intersection = typename IntersectionAlgorithm::Intersection;
+        Intersection intersection;
+
         if (IntersectionAlgorithm::intersection(geometryA, geometryB, intersection))
         {
-            for (int i = 0; i < intersection.size(); ++i)
-                intersections.emplace_back(eIdxA, eIdxB, std::move(intersection[i]));
+            static constexpr int dimIntersection = Policy::dimIntersection;
+
+            if (dimIntersection >= 2)
+            {
+                const auto triangulation = triangulate<2, dimworld>(intersection);
+                for (unsigned int i = 0; i < triangulation.size(); ++i)
+                    intersections.emplace_back(eIdxA, eIdxB, std::move(triangulation[i]));
+            }
+            else
+                intersections.emplace_back(eIdxA, eIdxB, intersection);
         }
     }
 

dumux/multidomain/facet/box/fluxvariablescache.hh

@@ -88,27 +88,26 @@ public:
             const auto& geometry = element.geometry();
             auto distVec = scvf.unitOuterNormal();
             distVec *= -1.0;
-            distVec *= diameter(geometry)*5.0; // make sure to be long enough
+            distVec *= diameter(geometry)*5.0; // make sure segment will be long enough
 
-            auto p1 = scvf.ipGlobal();
+            const auto p1 = scvf.ipGlobal();
             auto p2 = scvf.ipGlobal();
             p2 += distVec;
 
             static constexpr int dimWorld = GridView::dimensionworld;
             using Segment = Dune::MultiLinearGeometry<Scalar, 1, dimWorld>;
 
-            std::vector<GlobalPosition> corners({p1, p2});
-            Segment segment(Dune::GeometryTypes::line, corners);
+            using Policy = IntersectionPolicy::SegmentPolicy<typename GridView::ctype, dimWorld>;
+            using IntersectionAlgorithm = GeometryIntersection<typename Element::Geometry, Segment, Policy>;
+            typename IntersectionAlgorithm::Intersection intersection;
 
-            using Intersection = GeometryIntersection<typename Element::Geometry, Segment>;
-            typename Intersection::IntersectionType intersection;
-
-            if (!Intersection::intersection(geometry, segment, intersection))
+            Segment segment(Dune::GeometryTypes::line, std::vector<GlobalPosition>({p1, p2}));
+            if (!IntersectionAlgorithm::intersection(geometry, segment, intersection))
                 DUNE_THROW(Dune::InvalidStateException, "Could not compute interior integration point");
 
             // use center of intersection as integration point
-            ipGlobalInside_ = intersection[0][0];
-            ipGlobalInside_ += intersection[0][1];
+            ipGlobalInside_ = intersection[0];
+            ipGlobalInside_ += intersection[1];
             ipGlobalInside_ /= 2.0;
 
             fvGeometry.feLocalBasis().evaluateFunction(geometry.local(ipGlobalInside_), shapeValuesInside_);

test/common/geometry/test_2d3d_intersection.cc

@@ -11,6 +11,7 @@
 #include <dune/geometry/multilineargeometry.hh>
 
 #include <dumux/common/geometry/geometryintersection.hh>
+#include <dumux/common/geometry/triangulation.hh>
 #include <test/common/geometry/writetriangulation.hh>
 
 template<int dimworld = 3>
@@ -21,20 +22,22 @@ void testSegTriangle(const Dune::FieldVector<double, dimworld>& a,
                      const Dune::FieldVector<double, dimworld>& p,
                      bool expectIntersection = true)
 {
+    using GlobalPosition = Dune::FieldVector<double, dimworld>;
+    using CornerStorage = std::vector<GlobalPosition>;
     using TriGeometry = Dune::MultiLinearGeometry<double, 2, dimworld>;
     using SegGeometry = Dune::MultiLinearGeometry<double, 1, dimworld>;
-    using Test = Dumux::GeometryIntersection<TriGeometry, SegGeometry>;
+    using Policy = Dumux::IntersectionPolicy::PointPolicy<double, dimworld>;
+    using Test = Dumux::GeometryIntersection<TriGeometry, SegGeometry, Policy>;
     typename Test::IntersectionType intersection;
 
-    if (Test::template intersection<2>(a, b, c, p, q, intersection))
+    const auto tria = TriGeometry(Dune::GeometryTypes::triangle, CornerStorage({a, b, c}));
+    const auto seg = SegGeometry(Dune::GeometryTypes::line, CornerStorage({q, p}));
+    if (Test::intersection(tria, seg, intersection))
     {
-        if (intersection.size() != 1)
-            DUNE_THROW(Dune::InvalidStateException, "Found more than one intersection poin!");
-
         if (expectIntersection)
-            std::cout << "Found intersection point: " << intersection[0] << std::endl;
+            std::cout << "Found intersection point: " << intersection << std::endl;
         else
-            DUNE_THROW(Dune::InvalidStateException, "Found false positive: " << intersection[0]);
+            DUNE_THROW(Dune::InvalidStateException, "Found false positive: " << intersection);
     }
     else
     {
@@ -74,7 +77,7 @@ int main(int argc, char* argv[]) try
     using Geometry3D = Dune::MultiLinearGeometry<double, 3, dimworld>;
     using Geometry2D = Dune::MultiLinearGeometry<double, 2, dimworld>;
     using Test = Dumux::GeometryIntersection<Geometry3D, Geometry2D>;
-    typename Test::IntersectionType intersections;
+    typename Test::Intersection intersectionPolygon;
 
     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},
@@ -93,20 +96,22 @@ int main(int argc, char* argv[]) try
     Geometry2D quad(Dune::GeometryTypes::cube(dimworld-1), quadCorners);
     Geometry2D tri(Dune::GeometryTypes::simplex(dimworld-1), triCorners);
 
-    if (Test::intersection(cube, quad, intersections))
+    if (Test::intersection(cube, quad, intersectionPolygon))
     {
-        Dumux::writeVTKPolyDataTriangle(intersections, "quad_intersections");
-        if (intersections.size() != 4)
-            DUNE_THROW(Dune::InvalidStateException, "Should be 4 intersections!");
+        const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
+        Dumux::writeVTKPolyDataTriangle(triangulation, "quad_intersections");
+        if (triangulation.size() != 4)
+            DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 4 intersections!");
     }
     else
         DUNE_THROW(Dune::InvalidStateException, "No intersections found!");
 
-    if (Test::intersection(cube, tri, intersections))
+    if (Test::intersection(cube, tri, intersectionPolygon))
     {
-        Dumux::writeVTKPolyDataTriangle(intersections, "tri_intersections");
-        if (intersections.size() != 6)
-            DUNE_THROW(Dune::InvalidStateException, "Should be 4 intersections!");
+        const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
+        Dumux::writeVTKPolyDataTriangle(triangulation, "tri_intersections");
+        if (triangulation.size() != 6)
+            DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 6 intersections!");
     }
     else
         DUNE_THROW(Dune::InvalidStateException, "No intersections found!");