From 18e9b2751eebf918bc2f689106793bce0294d0a7 Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Mon, 1 Apr 2019 09:56:46 +0200
Subject: [PATCH] [geometryisection] introduce policy classes for intersections
---
dumux/common/geometry/geometryintersection.hh | 688 ++++++++++++------
dumux/common/geometry/intersectingentities.hh | 25 +-
.../facet/box/fluxvariablescache.hh | 19 +-
.../common/geometry/test_2d3d_intersection.cc | 37 +-
4 files changed, 534 insertions(+), 235 deletions(-)
diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh
index 4235e2037b..ab77faece4 100644
--- a/dumux/common/geometry/geometryintersection.hh
+++ b/dumux/common/geometry/geometryintersection.hh
@@ -23,6 +23,7 @@
#ifndef DUMUX_GEOMETRY_INTERSECTION_HH
#define DUMUX_GEOMETRY_INTERSECTION_HH
+#include <dune/common/deprecated.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/geometry/referenceelements.hh>
@@ -33,6 +34,65 @@
#include <dumux/common/geometry/grahamconvexhull.hh>
namespace Dumux {
+namespace IntersectionPolicy {
+
+//! Policy structure for point-like intersections
+template<class ct, int dw>
+struct PointPolicy
+{
+ using ctype = ct;
+ using Point = Dune::FieldVector<ctype, dw>;
+
+ static constexpr int dimWorld = dw;
+ static constexpr int dimIntersection = 0;
+ using Intersection = Point;
+};
+
+//! Policy structure for segment-like intersections
+template<class ct, int dw>
+struct SegmentPolicy
+{
+ using ctype = ct;
+ using Point = Dune::FieldVector<ctype, dw>;
+ using Segment = std::array<Point, 2>;
+
+ static constexpr int dimWorld = dw;
+ static constexpr int dimIntersection = 1;
+ using Intersection = Segment;
+};
+
+//! Policy structure for polygon-like intersections
+template<class ct, int dw>
+struct PolygonPolicy
+{
+ using ctype = ct;
+ using Point = Dune::FieldVector<ctype, dw>;
+ using Polygon = std::vector<Point>;
+
+ static constexpr int dimWorld = dw;
+ static constexpr int dimIntersection = 2;
+ using Intersection = Polygon;
+};
+
+//! default policy chooser class
+template<class Geometry1, class Geometry2, int dim2 = Geometry2::mydimension>
+class DefaultPolicyChooser
+{
+ static constexpr int dimworld = Geometry1::coorddimension;
+ static_assert(dimworld == int(Geometry2::coorddimension), "Geometries must have the same coordinate dimension!");
+ static_assert(dim2 == 1 || dim2 == 2, "No chooser class specialization available for given dimensions");
+ using ctype = typename Dune::PromotionTraits<typename Geometry1::ctype, typename Geometry2::ctype>::PromotedType;
+public:
+ using type = std::conditional_t< dim2 == 2,
+ PolygonPolicy<ctype, Geometry1::coorddimension>,
+ SegmentPolicy<ctype, Geometry1::coorddimension> >;
+};
+
+//! Helper alias to define the default intersection policy
+template<class Geometry1, class Geometry2>
+using DefaultPolicy = typename DefaultPolicyChooser<Geometry1, Geometry2>::type;
+
+} // end namespace IntersectionPolicy
/*!
* \ingroup Geometry
@@ -42,16 +102,26 @@ namespace Dumux {
*/
template
<class Geometry1, class Geometry2,
- int dimworld = Geometry1::coorddimension,
- int dim1 = Geometry1::mydimension,
- int dim2 = Geometry2::mydimension>
+ class Policy = IntersectionPolicy::DefaultPolicy<Geometry1, Geometry2>,
+ int dimworld = Geometry1::coorddimension,
+ int dim1 = Geometry1::mydimension,
+ int dim2 = Geometry2::mydimension>
class GeometryIntersection
{
+ static constexpr int dimWorld = Policy::dimWorld;
+
public:
- //! 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)
+ 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;
+
+ //! Determine if the two geometries intersect and compute the intersection geometry
+ static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
{
+ static_assert(int(dim1) >= int(dim2), "Geometries must be ordered such that dim1 >= dim2");
static_assert(dimworld == Geometry2::coorddimension, "Can only intersect geometries of same coordinate dimension");
DUNE_THROW(Dune::NotImplemented, "Geometry intersection detection with intersection computation not implemented for dimworld = "
<< dimworld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
@@ -62,17 +132,20 @@ public:
* \ingroup Geometry
* \brief A class for polygon--segment intersection in 2d space
*/
-template <class Geometry1, class Geometry2>
-class GeometryIntersection<Geometry1, Geometry2, 2, 2, 1>
+template <class Geometry1, class Geometry2, class Policy>
+class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 2, 1>
{
enum { dimworld = 2 };
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::array<std::vector<Point>, 1>;
+ 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
@@ -84,110 +157,161 @@ public:
* \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.
+ * \note This overload is used when segment-like intersections are seeked.
*/
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, int> = 0>
+ static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection)
{
- const auto a = geo2.corner(0);
- const auto b = geo2.corner(1);
- const auto d = b - a;
+ static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
- const std::vector<std::array<int, 2>> edges = [&]()
+ ctype tfirst, tlast;
+ if (intersect_(geo1, geo2, tfirst, tlast))
{
- std::vector<std::array<int, 2>> edges;
- switch (geo1.corners())
- {
- case 4: // quadrilateral
- edges = {{0, 2}, {3, 1}, {1, 0}, {2, 3}};
- break;
- case 3: // triangle
- edges = {{1, 0}, {0, 2}, {2, 1}};
- break;
- default:
- DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
- << geo1.type() << " with "<< geo1.corners() << " corners.");
- }
+ // the intersection exists. Export the intersections geometry now:
+ // s(t) = a + t(b-a) in [tfirst, tlast]
+ intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)});
+ return true;
+ }
- return edges;
- } ();
+ return false;
+ }
- // The initial interval is the whole segment
- // afterward we start clipping the interval
- // by the lines decribed by the edges
- ctype tfirst = 0.0;
- ctype tlast = 1.0;
+ /*!
+ * \brief Colliding segment and convex polygon
+ * \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.
+ * \note This overload is used when point-like intersections are seeked.
+ */
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
+ static bool intersection(const Geometry1& geo1, const Geometry2& geo2, typename P::Intersection& intersection)
+ {
+ static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension");
- const auto center1 = geo1.center();
- for (const auto& e : edges)
+ ctype tfirst, tlast;
+ if (!intersect_(geo1, geo2, tfirst, tlast))
{
- // compute outer normal vector of the edge
- const auto c0 = geo1.corner(e[0]);
- const auto c1 = geo1.corner(e[1]);
- const auto edge = c1 - c0;
- const auto eps = eps_*c0.two_norm();
-
- Dune::FieldVector<ctype, dimworld> n({-1.0*edge[1], edge[0]});
- n /= n.two_norm();
-
- // orientation of the element is not known
- // make sure normal points outwards of element
- if ( n*(center1-c0) > 0.0 )
- n *= -1.0;
-
- const ctype denom = n*d;
- const ctype dist = n*(a-c0);
-
- // if denominator is zero the segment in parallel to the edge.
- // If the distance is positive there is no intersection
+ // if start & end point are same, the intersection is a point
using std::abs;
- if (abs(denom) < eps)
- {
- if (dist > eps)
- return false;
- }
- else // not parallel: compute line-line intersection
+ if ( tfirst > tlast - eps_ )
{
- const ctype t = -dist / denom;
- // if entering half space cut tfirst if t is larger
- using std::signbit;
- if (signbit(denom))
- {
- if (t > tfirst)
- tfirst = t;
- }
- // if exiting half space cut tlast if t is smaller
- else
- {
- if (t < tlast)
- tlast = t;
- }
- // there is no intersection of the interval is empty
- // use unscaled epsilon since t is in local coordinates
- if (tfirst > tlast - eps_)
- return false;
+ intersection = Intersection(geo2.global(tfirst));
+ return true;
}
}
- // 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 = {std::vector<Point>({geo2.global(tfirst), geo2.global(tlast)})};
- return true;
+
+ return false;
}
+
+private:
+ /*!
+ * \brief Obtain local coordinates of start/end point of the intersecting segment.
+ */
+ static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast)
+ {
+ const auto a = geo2.corner(0);
+ const auto b = geo2.corner(1);
+ const auto d = b - a;
+
+ const std::vector<std::array<int, 2>> edges = [&]()
+ {
+ std::vector<std::array<int, 2>> edges;
+ switch (geo1.corners())
+ {
+ case 4: // quadrilateral
+ edges = {{0, 2}, {3, 1}, {1, 0}, {2, 3}};
+ break;
+ case 3: // triangle
+ edges = {{1, 0}, {0, 2}, {2, 1}};
+ break;
+ default:
+ DUNE_THROW(Dune::NotImplemented, "Collision of segment and geometry of type "
+ << geo1.type() << " with "<< geo1.corners() << " corners.");
+ }
+
+ return edges;
+ } ();
+
+ // The initial interval is the whole segment
+ // afterward we start clipping the interval
+ // by the lines decribed by the edges
+ tfirst = 0.0;
+ tlast = 1.0;
+
+ const auto center1 = geo1.center();
+ for (const auto& e : edges)
+ {
+ // compute outer normal vector of the edge
+ const auto c0 = geo1.corner(e[0]);
+ const auto c1 = geo1.corner(e[1]);
+ const auto edge = c1 - c0;
+ const auto eps = eps_*c0.two_norm();
+
+ Dune::FieldVector<ctype, dimworld> n({-1.0*edge[1], edge[0]});
+ n /= n.two_norm();
+
+ // orientation of the element is not known
+ // make sure normal points outwards of element
+ if ( n*(center1-c0) > 0.0 )
+ n *= -1.0;
+
+ const ctype denom = n*d;
+ const ctype dist = n*(a-c0);
+
+ // if denominator is zero the segment in parallel to the edge.
+ // If the distance is positive there is no intersection
+ using std::abs;
+ if (abs(denom) < eps)
+ {
+ if (dist > eps)
+ return false;
+ }
+ else // not parallel: compute line-line intersection
+ {
+ const ctype t = -dist / denom;
+ // if entering half space cut tfirst if t is larger
+ using std::signbit;
+ if (signbit(denom))
+ {
+ if (t > tfirst)
+ tfirst = t;
+ }
+ // if exiting half space cut tlast if t is smaller
+ else
+ {
+ if (t < tlast)
+ tlast = t;
+ }
+ // there is no intersection of the interval is empty
+ // use unscaled epsilon since t is in local coordinates
+ if (tfirst > tlast - eps_)
+ return false;
+ }
+ }
+
+ // If we made it until here an intersection exists.
+ return true;
+ }
};
/*!
* \ingroup Geometry
* \brief A class for polyhedron--segment intersection in 3d space
*/
-template <class Geometry1, class Geometry2>
-class GeometryIntersection<Geometry1, Geometry2, 3, 3, 1>
+template <class Geometry1, class Geometry2, class Policy>
+class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 3, 1>
{
enum { dimworld = 3 };
enum { dim1 = 3 };
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::array<std::vector<Point>, 1>;
+ 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
@@ -203,108 +327,161 @@ public:
* \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.
+ * \note This overload is used when segment-like intersections are seeked.
*/
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
+ 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");
- const auto a = geo2.corner(0);
- const auto b = geo2.corner(1);
- const auto d = b - a;
-
- // The initial interval is the whole segment
- // afterward we start clipping the interval
- // by the planes decribed by the facet
- ctype tfirst = 0.0;
- ctype tlast = 1.0;
-
- const std::vector<std::vector<int>> facets = [&]()
+ ctype tfirst, tlast;
+ if (intersect_(geo1, geo2, tfirst, tlast))
{
- std::vector<std::vector<int>> facets;
- // sort facet corner so that normal n = (p1-p0)x(p2-p0) always points outwards
- switch (geo1.corners())
- {
- // todo compute intersection geometries
- case 8: // hexahedron
- facets = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
- {3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
- break;
- case 4: // tetrahedron
- facets = {{1, 0, 2}, {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.");
- }
-
- return facets;
- }();
-
- for (const auto& f : facets)
- {
- // compute normal vector by cross product
- const auto v0 = geo1.corner(f[1]) - geo1.corner(f[0]);
- const auto v1 = geo1.corner(f[2]) - geo1.corner(f[0]);
- const auto eps = eps_*v0.two_norm();
+ // Intersection exists. Export the intersections geometry now:
+ // s(t) = a + t(b-a) in [tfirst, tlast]
+ intersection = Intersection({geo2.global(tfirst), geo2.global(tlast)});
+ return true;
+ }
- auto n = crossProduct(v0, v1);
- n /= n.two_norm();
+ return false;
+ }
- const ctype denom = n*d;
- const ctype dist = n*(a-geo1.corner(f[0]));
+ /*!
+ * \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.
+ * Basis is the theorem that for any two non-intersecting convex polyhedrons
+ * a separating plane exists.
+ * \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.
+ * \note This overload is used when point-like intersections are seeked.
+ */
+ 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");
- // if denominator is zero the segment in parallel to
- // the plane. If the distance is positive there is no intersection
+ ctype tfirst, tlast;
+ if (intersect_(geo1, geo2, tfirst, tlast))
+ {
+ // if start & end point are same, the intersection is a point
using std::abs;
- if (abs(denom) < eps)
- {
- if (dist > eps)
- return false;
- }
- else // not parallel: compute line-plane intersection
+ if ( abs(tfirst - tlast) < eps_ )
{
- const ctype t = -dist / denom;
- // if entering half space cut tfirst if t is larger
- using std::signbit;
- if (signbit(denom))
- {
- if (t > tfirst)
- tfirst = t;
- }
- // if exiting half space cut tlast if t is smaller
- else
- {
- if (t < tlast)
- tlast = t;
- }
- // there is no intersection of the interval is empty
- // use unscaled epsilon since t is in local coordinates
- if (tfirst > tlast - eps_)
- return false;
+ intersection = Intersection(geo2.global(tfirst));
+ return true;
}
}
- // 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 = {std::vector<Point>({geo2.global(tfirst), geo2.global(tlast)})};
- return true;
+
+ return false;
}
+
+private:
+ /*!
+ * \brief Obtain local coordinates of start/end point of the intersecting segment.
+ */
+ static bool intersect_(const Geometry1& geo1, const Geometry2& geo2, ctype& tfirst, ctype& tlast)
+ {
+ const auto a = geo2.corner(0);
+ const auto b = geo2.corner(1);
+ const auto d = b - a;
+
+ // The initial interval is the whole segment
+ // afterward we start clipping the interval
+ // by the planes decribed by the facet
+ tfirst = 0.0;
+ tlast = 1.0;
+
+ const std::vector<std::vector<int>> facets = [&]()
+ {
+ std::vector<std::vector<int>> facets;
+ // sort facet corner so that normal n = (p1-p0)x(p2-p0) always points outwards
+ switch (geo1.corners())
+ {
+ // todo compute intersection geometries
+ case 8: // hexahedron
+ facets = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
+ {3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
+ break;
+ case 4: // tetrahedron
+ facets = {{1, 0, 2}, {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.");
+ }
+
+ return facets;
+ }();
+
+ for (const auto& f : facets)
+ {
+ // compute normal vector by cross product
+ const auto v0 = geo1.corner(f[1]) - geo1.corner(f[0]);
+ const auto v1 = geo1.corner(f[2]) - geo1.corner(f[0]);
+ const auto eps = eps_*v0.two_norm();
+
+ auto n = crossProduct(v0, v1);
+ n /= n.two_norm();
+
+ const ctype denom = n*d;
+ const ctype dist = n*(a-geo1.corner(f[0]));
+
+ // if denominator is zero the segment in parallel to
+ // the plane. If the distance is positive there is no intersection
+ using std::abs;
+ if (abs(denom) < eps)
+ {
+ if (dist > eps)
+ return false;
+ }
+ else // not parallel: compute line-plane intersection
+ {
+ const ctype t = -dist / denom;
+ // if entering half space cut tfirst if t is larger
+ using std::signbit;
+ if (signbit(denom))
+ {
+ if (t > tfirst)
+ tfirst = t;
+ }
+ // if exiting half space cut tlast if t is smaller
+ else
+ {
+ if (t < tlast)
+ tlast = t;
+ }
+ // there is no intersection of the interval is empty
+ // use unscaled epsilon since t is in local coordinates
+ if (tfirst > tlast - eps_)
+ return false;
+ }
+ }
+
+ // If we made it until here an intersection exists.
+ return true;
+ }
};
/*!
* \ingroup Geometry
* \brief A class for polyhedron--polygon intersection in 3d space
*/
-template <class Geometry1, class Geometry2>
-class GeometryIntersection<Geometry1, Geometry2, 3, 3, 2>
+template <class Geometry1, class Geometry2, class Policy>
+class GeometryIntersection<Geometry1, Geometry2, Policy, 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>>;
+ 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
@@ -324,11 +501,12 @@ public:
* Return a triangulation of that polygon as intersection
* \param geo1/geo2 The geometries to intersect
* \param intersection A triangulation of the intersection polygon
+ * \note This overload is used when polygon like intersections are seeked
*/
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& intersection)
+ 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");
+ 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);
@@ -350,6 +528,9 @@ public:
const auto referenceElement1 = ReferenceElementsGeo1::general(geo1.type());
const auto referenceElement2 = ReferenceElementsGeo2::general(geo2.type());
+ // intersection policy for point-like intersections (used later)
+ using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;
+
// add intersection points of all polyhedron edges (codim dim-1) with the polygon
for (int i = 0; i < referenceElement1.size(dim1-1); ++i)
{
@@ -358,10 +539,10 @@ public:
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 polySegIntersection;
- if (PolySegTest::template intersection<2>(geo2, segGeo, polySegIntersection))
- points.emplace_back(polySegIntersection[0]);
+ using PolySegTest = GeometryIntersection<Geometry2, SegGeometry, PointPolicy>;
+ typename PolySegTest::Intersection polySegIntersection;
+ if (PolySegTest::intersection(geo2, segGeo, polySegIntersection))
+ points.emplace_back(polySegIntersection);
}
// add intersection points of all polygon faces (codim 1) with the polyhedron faces
@@ -397,10 +578,10 @@ public:
const auto segGeo = SegGeometry(Dune::GeometryTypes::line, std::vector<Point>{p, q});
- using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry>;
- typename PolySegTest::IntersectionType polySegIntersection;
- if (PolySegTest::template intersection<2>(faceGeo, segGeo, polySegIntersection))
- points.emplace_back(polySegIntersection[0]);
+ using PolySegTest = GeometryIntersection<PolyhedronFaceGeometry, SegGeometry, PointPolicy>;
+ typename PolySegTest::Intersection polySegIntersection;
+ if (PolySegTest::intersection(faceGeo, segGeo, polySegIntersection))
+ points.emplace_back(polySegIntersection);
}
}
@@ -425,31 +606,53 @@ public:
// 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);
- assert(!convexHull.empty());
+ // intersection polygon is convex hull of above points
+ intersection = grahamConvexHull2d3d(points);
+ assert(!intersection.empty());
- // the intersections are the triangulation of the convex hull polygon
- intersection = triangulateConvexHull(convexHull);
return true;
}
+
+ /*!
+ * \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
+ * \todo implement overloads for segment or point-like intersections
+ */
+ 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");
+ DUNE_THROW(Dune::NotImplemented, "Polyhedron-polygon intersection detection only implemented for polygon-like intersections");
+ }
};
/*!
* \ingroup Geometry
* \brief A class for polygon--segment intersection in 3d space
*/
-template <class Geometry1, class Geometry2>
-class GeometryIntersection<Geometry1, Geometry2, 3, 2, 1>
+template <class Geometry1, class Geometry2, class Policy>
+class GeometryIntersection<Geometry1, Geometry2, Policy, 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>;
+ 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
@@ -463,13 +666,11 @@ public:
* \param geo1/geo2 The geometries to intersect
* \param is If the geometries collide, is holds the corner points of
* the intersection object in global coordinates.
+ * \note This overload is used when point-like intersections are seeked
*/
- template<int dimIntersection>
- static bool intersection(const Geometry1& geo1, const Geometry2& geo2, IntersectionType& is)
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
+ static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& 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);
@@ -480,14 +681,24 @@ public:
const auto c = geo1.corner(2);
if (geo1.corners() == 3)
- return intersection<dimIntersection>(a, b, c, p, q, is);
+ return intersect_<Policy>(a, b, c, p, q, is);
else if (geo1.corners() == 4)
{
+ Intersection is1, is2;
+ bool hasSegment1, hasSegment2;
+
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;
+ const bool intersects1 = intersect_<Policy>(a, b, d, p, q, is1, hasSegment1);
+ const bool intersects2 = intersect_<Policy>(a, d, c, p, q, is2, hasSegment2);
+
+ if (hasSegment1 || hasSegment2)
+ return false;
+
+ if (intersects1) { is = is1; return true; }
+ if (intersects2) { is = is2; return true; }
+
+ return false;
}
else
@@ -495,14 +706,45 @@ public:
<< geo1.type() << ", "<< geo1.corners() << " corners.");
}
- // triangle--segment intersection with points as input
- template<int dimIntersection>
+ /*!
+ * \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 is If the geometries collide, is holds the corner points of
+ * the intersection object in global coordinates.
+ * \note This overload is used when point-like intersections are seeked
+ */
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
+ DUNE_DEPRECATED_MSG("Please use intersection(triangle, segment, ...) instead")
static bool intersection(const Point& a, const Point& b, const Point& c,
const Point& p, const Point& q,
- IntersectionType& is)
+ Intersection& is)
{
- if (dimIntersection != 2)
- DUNE_THROW(Dune::NotImplemented, "Only simplex intersections are currently implemented!");
+ return intersect_<Policy>(a, b, c, p, q, is);
+ }
+
+private:
+ // triangle--segment intersection with points as input
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
+ static bool intersect_(const Point& a, const Point& b, const Point& c,
+ const Point& p, const Point& q,
+ Intersection& is)
+ {
+ bool hasSegment;
+ return intersect_(a, b, c, p, q, is, hasSegment);
+ }
+
+ /*!
+ * \brief triangle--segment point-like intersection with points as input. Also
+ * stores if a segment-like intersection was found in the provided boolean.
+ */
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0>
+ static bool intersect_(const Point& a, const Point& b, const Point& c,
+ const Point& p, const Point& q,
+ Intersection& is, bool& hasSegmentIs)
+ {
+ hasSegmentIs = false;
const auto ab = b - a;
const auto ac = c - a;
@@ -514,10 +756,48 @@ public:
// 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();
+ const auto abnorm2 = ab.two_norm2();
+ const auto eps = eps_*abnorm2*qp.two_norm();
using std::abs;
if (abs(denom) < eps)
+ {
+ const auto pa = a - p;
+ const auto denom2 = normal*pa;
+ if (abs(denom2) > eps_*pa.two_norm()*abnorm2)
+ return false;
+
+ // if we get here, we are in-plane. Check if a
+ // segment-like intersection with geometry 1 exists.
+ using SegmentPolicy = typename IntersectionPolicy::SegmentPolicy<ctype, dimworld>;
+ using Triangle = Dune::AffineGeometry<ctype, 2, dimworld>;
+ using Segment = Dune::AffineGeometry<ctype, 1, dimworld>;
+ using SegmentIntersectionAlgorithm = GeometryIntersection<Triangle, Segment, SegmentPolicy>;
+ using SegmentIntersectionType = typename SegmentIntersectionAlgorithm::Intersection;
+ SegmentIntersectionType segmentIs;
+
+ Triangle triangle(Dune::GeometryTypes::simplex(2), std::array<Point, 3>({a, b, c}));
+ Segment segment(Dune::GeometryTypes::simplex(1), std::array<Point, 2>({p, q}));
+ if (SegmentIntersectionAlgorithm::intersection(triangle, segment, segmentIs))
+ {
+ hasSegmentIs = true;
+ return false;
+ }
+
+ // We are in-plane. A point-like
+ // intersection can only be on the edges
+ for (const auto& ip : {p, q})
+ {
+ if ( intersectsPointSimplex(ip, a, b)
+ || intersectsPointSimplex(ip, b, c)
+ || intersectsPointSimplex(ip, c, a) )
+ {
+ is = ip;
+ return true;
+ }
+ }
+
return false;
+ }
// compute intersection t value of pq with plane of triangle.
// a segment intersects if and only if 0 <= t <= 1.
@@ -542,7 +822,7 @@ public:
ip.axpy(u, a);
ip.axpy(v, b);
ip.axpy(w, c);
- is = {ip};
+ is = ip;
return true;
}
};
diff --git a/dumux/common/geometry/intersectingentities.hh b/dumux/common/geometry/intersectingentities.hh
index f4990f46eb..eb7189e52e 100644
--- a/dumux/common/geometry/intersectingentities.hh
+++ b/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);
}
}
diff --git a/dumux/multidomain/facet/box/fluxvariablescache.hh b/dumux/multidomain/facet/box/fluxvariablescache.hh
index 23d2d78bcb..4074e68bc5 100644
--- a/dumux/multidomain/facet/box/fluxvariablescache.hh
+++ b/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_);
diff --git a/test/common/geometry/test_2d3d_intersection.cc b/test/common/geometry/test_2d3d_intersection.cc
index a8a681af49..3b2af1dad0 100644
--- a/test/common/geometry/test_2d3d_intersection.cc
+++ b/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!");
--
GitLab