From 31db57ef2acdca056cc27db3fdf785e51448e9bf Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Thu, 30 May 2019 22:10:03 +0200
Subject: [PATCH] [geometryintersection] add algorithm for 2d-2d intersections
in 2d space
---
dumux/common/geometry/geometryintersection.hh | 142 ++++++++++++++++++
1 file changed, 142 insertions(+)
diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh
index d7c1336abf..2502ec1a6d 100644
--- a/dumux/common/geometry/geometryintersection.hh
+++ b/dumux/common/geometry/geometryintersection.hh
@@ -484,6 +484,148 @@ public:
}
};
+/*!
+ * \ingroup Geometry
+ * \brief A class for polygon--polygon intersection in 2d space
+ */
+template <class Geometry1, class Geometry2, class Policy>
+class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 2, 2>
+{
+ enum { dimworld = 2 };
+ enum { dim1 = 2 };
+ enum { dim2 = 2 };
+
+public:
+ using ctype = typename Policy::ctype;
+ using Point = typename Policy::Point;
+ using Intersection = typename Policy::Intersection;
+
+ //! Deprecated alias, will be removed after 3.1
+ using IntersectionType DUNE_DEPRECATED_MSG("Please use Intersection instead") = Intersection;
+
+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 two polygons
+ * \note First we find the vertex candidates for the intersection region as follows:
+ * Add polygon vertices that are inside the other polygon
+ * Add intersections of polygon edges
+ * 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 Container to store the corner points of the polygon (as convex hull)
+ * \note This overload is used when polygon like intersections are seeked
+ */
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, int> = 0>
+ static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& 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(6);
+
+ // add polygon1 corners that are inside polygon2
+ for (int i = 0; i < geo1.corners(); ++i)
+ if (intersectsPointGeometry(geo1.corner(i), geo2))
+ points.emplace_back(geo1.corner(i));
+
+ // add polygon2 corners that are inside polygon1
+ for (int i = 0; i < geo2.corners(); ++i)
+ if (intersectsPointGeometry(geo2.corner(i), geo1))
+ points.emplace_back(geo2.corner(i));
+
+ if (points.empty())
+ return false;
+
+ const auto referenceElement1 = ReferenceElementsGeo1::general(geo1.type());
+ const auto referenceElement2 = ReferenceElementsGeo2::general(geo2.type());
+
+ // add intersections of edges
+ using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>;
+ using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;
+ for (int i = 0; i < referenceElement1.size(dim1-1); ++i)
+ {
+ const auto localEdgeGeom1 = referenceElement1.template geometry<dim1-1>(i);
+ const auto edge1 = SegGeometry( Dune::GeometryTypes::line,
+ std::vector<Point>( {geo1.global(localEdgeGeom1.corner(0)),
+ geo1.global(localEdgeGeom1.corner(1))} ));
+
+ for (int j = 0; j < referenceElement2.size(dim2-1); ++j)
+ {
+ const auto localEdgeGeom2 = referenceElement2.template geometry<dim2-1>(j);
+ const auto edge2 = SegGeometry( Dune::GeometryTypes::line,
+ std::vector<Point>( {geo2.global(localEdgeGeom2.corner(0)),
+ geo2.global(localEdgeGeom2.corner(1))} ));
+
+ using EdgeTest = GeometryIntersection<SegGeometry, SegGeometry, PointPolicy>;
+ typename EdgeTest::Intersection edgeIntersection;
+ if (EdgeTest::intersection(edge1, edge2, edgeIntersection))
+ points.emplace_back(edgeIntersection);
+ }
+ }
+
+ 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
+ {
+ using std::abs;
+ return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : a[1] < b[1]);
+ });
+
+ 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 at least three unique points
+ if (points.size() < 3)
+ return false;
+
+ // intersection polygon is convex hull of above points
+ intersection = grahamConvexHull(points);
+ assert(!intersection.empty());
+ return true;
+ }
+
+ /*!
+ * \brief Colliding two polygons
+ * \param geo1/geo2 The geometries to intersect
+ * \param intersection Container to store the corners of intersection segment
+ * \note this overload is used when segment-like intersections are seeked
+ * \todo implement this query
+ */
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
+ static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
+ {
+ static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
+ DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for segment-like intersections");
+ }
+
+ /*!
+ * \brief Colliding two polygons
+ * \param geo1/geo2 The geometries to intersect
+ * \param intersection The intersection point
+ * \note this overload is used when point-like intersections are seeked
+ * \todo implement this query
+ */
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
+ static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
+ {
+ static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
+ DUNE_THROW(Dune::NotImplemented, "Polygon-polygon intersection detection for touching points");
+ }
+};
+
/*!
* \ingroup Geometry
* \brief A class for polyhedron--segment intersection in 3d space
--
GitLab