[geometry] Implement polyhedron-polygon and segment-polygon intersections in 3D

dumux/common/geometry/boundingboxtreeintersection.hh

@@ -45,12 +45,13 @@ class BoundingBoxTreeIntersection
     using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
 
 public:
+    template<class Corners>
     explicit BoundingBoxTreeIntersection(std::size_t a,
                                          std::size_t b,
-                                         std::vector<GlobalPosition>&& c)
+                                         Corners&& c)
     : a_(a)
     , b_(b)
-    , corners_(std::move(c))
+    , corners_(c.begin(), c.end())
     {
         static_assert(int(EntitySet0::dimensionworld) == int(EntitySet1::dimensionworld),
           "Can only store intersections of entity sets with the same world dimension");

dumux/common/geometry/geometryintersection.hh

@@ -23,14 +23,16 @@
 #ifndef DUMUX_GEOMETRY_INTERSECTION_HH
 #define DUMUX_GEOMETRY_INTERSECTION_HH
 
+#include <dune/common/version.hh>
 #include <dune/common/exceptions.hh>
 #include <dune/common/promotiontraits.hh>
 #include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/math.hh>
+#include <dumux/common/geometry/intersectspointgeometry.hh>
+#include <dumux/common/geometry/grahamconvexhull.hh>
 
-namespace Dumux
-{
+namespace Dumux {
 
 /*!
  * \ingroup Common
@@ -46,11 +48,8 @@ template
 class GeometryIntersection
 {
 public:
-    using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
-    using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
-    using IntersectionType = std::vector<GlobalPosition>;
-
     //! Determine if the two geometries intersect and compute the intersection corners
+    template<class IntersectionType>
     static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
     {
         static_assert(dimworld == Geometry2::coorddimension, "Can only intersect geometries of same coordinate dimension");
@@ -59,19 +58,18 @@ public:
     }
 };
 
-//! Geometry collision detection with 3d and 1d geometry in 3d space
+//! polyhedron--segment intersection in 3d space
 template <class Geometry1, class Geometry2>
 class GeometryIntersection<Geometry1, Geometry2, 3, 3, 1>
 {
     enum { dimworld = 3 };
     enum { dim1 = 3 };
     enum { dim2 = 1 };
-    enum { dimis = dim2 }; // the intersection dimension
 
 public:
     using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
-    using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
-    using IntersectionType = std::vector<GlobalPosition>;
+    using Point = Dune::FieldVector<ctype, dimworld>;
+    using IntersectionType = std::array<std::vector<Point>, 1>;
 
 private:
     static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
@@ -80,7 +78,8 @@ private:
 public:
     /*!
      *  \brief Colliding segment and convex polyhedron
-     *  \note Algorithm from "Real-Time Collision Detection" by Christer Ericson
+     *  \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson,
+     *        published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc.
      *        Basis is the theorem that for any two non-intersecting convex polyhedrons
      *        a separating plane exists.
      * \param geo1/geo2 The geometries to intersect
@@ -116,7 +115,8 @@ public:
                     facets = {{1, 2, 0}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}};
                     break;
                 default:
-                    DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type " << geo1.type() << ", "<< geo1.corners() << " corners.");
+                    DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
+                                   << geo1.type() << ", "<< geo1.corners() << " corners.");
             }
 
             return facets;
@@ -168,7 +168,258 @@ public:
         }
         // If we made it until here an intersection exists. We also export
         // the intersections geometry now s(t) = a + t(b-a) in [tfirst, tlast]
-        intersection = {geo2.global(tfirst), geo2.global(tlast)};
+        intersection = {std::vector<Point>({geo2.global(tfirst), geo2.global(tlast)})};
+        return true;
+    }
+};
+
+//! polyhedron--polygon intersection in 3d space
+template <class Geometry1, class Geometry2>
+class GeometryIntersection<Geometry1, Geometry2, 3, 3, 2>
+{
+    enum { dimworld = 3 };
+    enum { dim1 = 3 };
+    enum { dim2 = 2 };
+
+public:
+    using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
+    using Point = Dune::FieldVector<ctype, dimworld>;
+    using IntersectionType = std::vector<std::array<Point, 3>>;
+
+private:
+    static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
+    using ReferenceElementsGeo1 = typename Dune::ReferenceElements<ctype, dim1>;
+    using ReferenceElementsGeo2 = typename Dune::ReferenceElements<ctype, dim2>;
+
+public:
+    /*!
+     * \brief Colliding segment and convex polyhedron
+     * \note First we find the vertex candidates for the intersection region as follows:
+     *       Add triangle vertices that are inside the tetrahedron
+     *       Add tetrahedron vertices that are inside the triangle
+     *       Add all intersection points of tetrahedron edges (codim 2) with the triangle (codim 0) (6*1 tests)
+     *       Add all intersection points of triangle edges (codim 1) with tetrahedron faces (codim 1) (3*4 tests)
+     *       Remove duplicate points from the list
+     *       Compute the convex hull polygon
+     *       Return a triangulation of that polygon as intersection
+     * \param geo1/geo2 The geometries to intersect
+     * \param intersection A triangulation of the intersection polygon
+     */
+    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
+    {
+        static_assert(int(dimworld) == int(Geometry2::coorddimension),
+                      "Can only collide geometries of same coordinate dimension");
+
+        // the candidate intersection points
+        std::vector<Point> points; points.reserve(10);
+
+        // add 3d geometry corners that are inside the 2d geometry
+        for (int i = 0; i < geo1.corners(); ++i)
+            if (intersectsPointGeometry(geo1.corner(i), geo2))
+                points.emplace_back(geo1.corner(i));
+
+        // add 2d geometry corners that are inside the 3d geometry
+        for (int i = 0; i < geo2.corners(); ++i)
+            if (intersectsPointGeometry(geo2.corner(i), geo1))
+                points.emplace_back(geo2.corner(i));
+
+        // get some geometry types
+        using PolyhedronFaceGeometry = Dune::MultiLinearGeometry<ctype, 2, dimworld>;
+        using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>;
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON,2,6)
+        const auto referenceElement1 = ReferenceElementsGeo1::general(geo1.type());
+        const auto referenceElement2 = ReferenceElementsGeo2::general(geo2.type());
+#else
+        const auto& referenceElement1 = ReferenceElementsGeo1::general(geo1.type());
+        const auto& referenceElement2 = ReferenceElementsGeo2::general(geo2.type());
+#endif
+
+        // add intersection points of all polyhedron edges (codim dim-1) with the polygon
+        for (int i = 0; i < referenceElement1.size(dim1-1); ++i)
+        {
+            const auto localEdgeGeom = referenceElement1.template geometry<dim1-1>(i);
+            const auto p = geo1.global(localEdgeGeom.corner(0));
+            const auto q = geo1.global(localEdgeGeom.corner(1));
+            const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q});
+
+            using PolySegTest = GeometryIntersection<Geometry2, SegGeometry>;
+            typename PolySegTest::IntersectionType intersection;
+            if (PolySegTest::template intersection<2>(geo2, segGeo, intersection))
+                points.emplace_back(intersection[0]);
+        }
+
+        // add intersection points of all polygon faces (codim 1) with the polyhedron faces
+        for (int i = 0; i < referenceElement1.size(1); ++i)
+        {
+            const auto faceGeo = [&]()
+            {
+                const auto localFaceGeo = referenceElement1.template geometry<1>(i);
+                if (localFaceGeo.corners() == 4)
+                {
+                    const auto a = geo1.global(localFaceGeo.corner(0));
+                    const auto b = geo1.global(localFaceGeo.corner(1));
+                    const auto c = geo1.global(localFaceGeo.corner(2));
+                    const auto d = geo1.global(localFaceGeo.corner(3));
+                    return PolyhedronFaceGeometry(Dune::GeometryTypes::cube(2), std::vector<Point>{a, b, c, d});
+                }
+                else
+                {
+                    const auto a = geo1.global(localFaceGeo.corner(0));
+                    const auto b = geo1.global(localFaceGeo.corner(1));
+                    const auto c = geo1.global(localFaceGeo.corner(2));
+                    return PolyhedronFaceGeometry(Dune::GeometryTypes::simplex(2), std::vector<Point>{a, b, c});
+                }
+            }();
+
+            for (int j = 0; j < referenceElement2.size(1); ++j)
+            {
+                const auto localEdgeGeom = referenceElement2.template geometry<1>(j);
+                const auto p = geo2.global(localEdgeGeom.corner(0));
+                const auto q = geo2.global(localEdgeGeom.corner(1));
+                const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q});
+
+                using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry>;
+                typename PolySegTest::IntersectionType intersection;
+                if (PolySegTest::template intersection<2>(faceGeo, segGeo, intersection))
+                    points.emplace_back(intersection[0]);
+            }
+        }
+
+        // return if no intersection points were found
+        if (points.empty()) return false;
+
+        // remove duplicates
+        const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_;
+        std::sort(points.begin(), points.end(), [&eps](const auto& a, const auto& b) -> bool
+        {
+            return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : (abs(a[1]-b[1]) > eps ? a[1] < b[1] : (a[2] < b[2])));
+        });
+
+        auto removeIt = std::unique(points.begin(), points.end(), [&eps](const auto& a, const auto&b)
+        {
+            return (b-a).two_norm() < eps;
+        });
+
+        points.erase(removeIt, points.end());
+
+        // return false if we don't have more than three unique points
+        if (points.size() < 3) return false;
+
+        // compute convex hull
+        const auto convexHull = grahamConvexHull2d3d(points);
+
+        // the intersections are the triangulation of the convex hull polygon
+        intersection = triangulateConvexHull(convexHull);
+        return true;
+    }
+};
+
+//! polygon--segment intersection in 3d space
+template <class Geometry1, class Geometry2>
+class GeometryIntersection<Geometry1, Geometry2, 3, 2, 1>
+{
+    enum { dimworld = 3 };
+    enum { dim1 = 2 };
+    enum { dim2 = 1 };
+
+public:
+    using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
+    using Point = Dune::FieldVector<ctype, dimworld>;
+    using IntersectionType = std::vector<Point>;
+
+private:
+    static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
+    using ReferenceElements = typename Dune::ReferenceElements<ctype, dim1>;
+
+public:
+    /*!
+     *  \brief Colliding segment and convex polyhedron
+     *  \note Algorithm based on the one from "Real-Time Collision Detection" by Christer Ericson,
+     *        published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. (Chapter 5.3.6)
+     * \param geo1/geo2 The geometries to intersect
+     * \param intersection If the geometries collide intersection holds the corner points of
+     *        the intersection object in global coordinates.
+     */
+    template<int dimIntersection>
+    static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& is)
+    {
+        if (dimIntersection != 2)
+            DUNE_THROW(Dune::NotImplemented, "Only simplex intersections are currently implemented!");
+
+        static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
+
+        const auto p = geo2.corner(0);
+        const auto q = geo2.corner(1);
+
+        const auto a = geo1.corner(0);
+        const auto b = geo1.corner(1);
+        const auto c = geo1.corner(2);
+
+        if (geo1.corners() == 3)
+            return intersection<dimIntersection>(a, b, c, p, q, is);
+
+        else if (geo1.corners() == 4)
+        {
+            const auto d = geo1.corner(3);
+            if (intersection<dimIntersection>(a, b, d, p, q, is)) return true;
+            else if (intersection<dimIntersection>(a, d, c, p, q, is)) return true;
+            else return false;
+        }
+
+        else
+            DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
+                          << geo1.type() << ", "<< geo1.corners() << " corners.");
+    }
+
+    // triangle--segment intersection with points as input
+    template<int dimIntersection>
+    static bool intersection(const Point& a, const Point& b, const Point& c,
+                             const Point& p, const Point& q,
+                             IntersectionType& is)
+    {
+        if (dimIntersection != 2)
+            DUNE_THROW(Dune::NotImplemented, "Only simplex intersections are currently implemented!");
+
+        const auto ab = b - a;
+        const auto ac = c - a;
+        const auto qp = p - q;
+
+        // compute the triangle normal that defines the triangle plane
+        const auto normal = crossProduct(ab, ac);
+
+        // compute the denominator
+        // if denom is 0 the segment is parallel and we can return
+        const auto denom = normal*qp;
+        const auto eps = eps_*ab.two_norm2()*qp.two_norm();
+        using std::abs;
+        if (abs(denom) < eps)
+            return false;
+
+        // compute intersection t value of pq with plane of triangle.
+        // a segment intersects if and only if 0 <= t <= 1.
+        const auto ap = p - a;
+        const auto t = (ap*normal)/denom;
+        if (t < 0.0 - eps_) return false;
+        if (t > 1.0 + eps_) return false;
+
+        // compute the barycentric coordinates and check if the intersection point
+        // is inside the bounds of the triangle
+        const auto e = crossProduct(qp, ap);
+        const auto v = (ac*e)/denom;
+        if (v < -eps_ || v > 1.0 + eps_) return false;
+        const auto w = -(ab*e)/denom;
+        if (w < -eps_ || v + w > 1.0 + eps_) return false;
+
+        // Now we are sure there is an intersection points
+        // Perform delayed division compute the last barycentric coordinate component
+        const auto u = 1.0 - v - w;
+
+        Point ip(0.0);
+        ip.axpy(u, a);
+        ip.axpy(v, b);
+        ip.axpy(w, c);
+        is = {ip};
         return true;
     }
 };

dumux/common/geometry/intersectingentities.hh

@@ -145,7 +145,10 @@ void intersectingEntities(const BoundingBoxTree<EntitySet0>& treeA,
                                                             std::decay_t<decltype(geometryB)> >;
         typename IntersectionAlgorithm::IntersectionType intersection;
         if (IntersectionAlgorithm::intersection(geometryA, geometryB, intersection))
-            intersections.emplace_back(eIdxA, eIdxB, std::move(intersection));
+        {
+            for (int i = 0; i < intersection.size(); ++i)
+                intersections.emplace_back(eIdxA, eIdxB, std::move(intersection[i]));
+        }
     }
 
     // if we reached the leaf in treeA, just continue in treeB