diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh
index e3199636d76d4bd51c547223d0331a588aa5d298..69b5c894d3cd76dbf8e265d0708f5ce3e3c28595 100644
--- a/dumux/common/geometry/geometryintersection.hh
+++ b/dumux/common/geometry/geometryintersection.hh
@@ -32,6 +32,7 @@
#include <dumux/common/math.hh>
#include <dumux/common/geometry/intersectspointgeometry.hh>
#include <dumux/common/geometry/grahamconvexhull.hh>
+#include <dumux/common/geometry/boundingboxtree.hh>
namespace Dumux {
namespace IntersectionPolicy {
@@ -214,6 +215,126 @@ public:
}
};
+/*!
+ * \ingroup Geometry
+ * \brief A class for segment--segment intersection in 2d space
+ */
+template <class Geometry1, class Geometry2, class Policy>
+class GeometryIntersection<Geometry1, Geometry2, Policy, 2, 1, 1>
+{
+ enum { dimworld = 2 };
+ enum { dim1 = 1 };
+ enum { dim2 = 1 };
+
+ // base epsilon for floating point comparisons
+ static constexpr typename Policy::ctype eps_ = 1.5e-7;
+
+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;
+
+ /*!
+ * \brief Colliding two segments
+ * \param geo1/geo2 The geometries to intersect
+ * \param intersection The intersection point
+ * \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");
+
+ const auto v1 = geo1.corner(1) - geo1.corner(0);
+ const auto v2 = geo2.corner(1) - geo2.corner(0);
+ const auto ac = geo2.corner(0) - geo1.corner(0);
+
+ auto n2 = Point({-1.0*v2[1], v2[0]});
+ n2 /= n2.two_norm();
+
+ // compute distance of first corner on geo2 to line1
+ const auto dist12 = n2*ac;
+
+ // first check parallel segments
+ using std::abs;
+ using std::sqrt;
+
+ const auto v1Norm2 = v1.two_norm2();
+ const auto eps = eps_*sqrt(v1Norm2);
+ const auto eps2 = eps_*v1Norm2;
+
+ const auto sp = n2*v1;
+ if (abs(sp) < eps)
+ {
+ if (abs(dist12) > eps)
+ return false;
+
+ // point intersection can only be one of the corners
+ if ( (v1*v2) < 0.0 )
+ {
+ if ( ac.two_norm2() < eps2 )
+ { intersection = geo2.corner(0); return true; }
+
+ if ( (geo2.corner(1) - geo1.corner(1)).two_norm2() < eps2 )
+ { intersection = geo2.corner(1); return true; }
+ }
+ else
+ {
+ if ( (geo2.corner(1) - geo1.corner(0)).two_norm2() < eps2 )
+ { intersection = geo2.corner(1); return true; }
+
+ if ( (geo2.corner(0) - geo1.corner(1)).two_norm2() < eps2 )
+ { intersection = geo2.corner(0); return true; }
+ }
+
+ // no intersection
+ return false;
+ }
+
+ // intersection point on v1 in local coords
+ const auto t1 = dist12 / sp;
+
+ // check if the local coords are valid
+ if (t1 < -1.0*eps_ || t1 > 1.0 + eps_)
+ return false;
+
+ // compute global coordinates
+ auto isPoint = geo1.global(t1);
+
+ // check if point is in bounding box of v2
+ const auto c = geo2.corner(0);
+ const auto d = geo2.corner(1);
+
+ using std::min; using std::max;
+ std::array<ctype, 4> bBox({ min(c[0], d[0]), min(c[1], d[1]), max(c[0], d[0]), max(c[1], d[1]) });
+
+ if ( intersectsPointBoundingBox(isPoint, bBox.data()) )
+ {
+ intersection = std::move(isPoint);
+ return true;
+ }
+
+ return false;
+ }
+
+ /*!
+ * \brief Colliding two segments
+ * \param geo1/geo2 The geometries to intersect
+ * \param intersection Container to store 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, "segment-segment intersection detection for segment-like intersections");
+ }
+};
+
/*!
* \ingroup Geometry
* \brief A class for polygon--segment intersection in 2d space
@@ -363,6 +484,155 @@ 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));
+
+ const auto numPoints1 = points.size();
+ if (numPoints1 != geo1.corners())
+ {
+ // 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;
+
+ if (points.size() - numPoints1 != geo2.corners())
+ {
+ 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<2>(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
@@ -551,7 +821,7 @@ public:
* 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
+ * \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>
@@ -658,7 +928,7 @@ public:
if (points.size() < 3) return false;
// intersection polygon is convex hull of above points
- intersection = grahamConvexHull2d3d(points);
+ intersection = grahamConvexHull<2>(points);
assert(!intersection.empty());
return true;
@@ -675,7 +945,7 @@ public:
* 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
+ * \param intersection Container to store the intersection result
* \todo implement overloads for segment or point-like intersections
*/
template<class P = Policy, std::enable_if_t<P::dimIntersection != 2, int> = 0>
diff --git a/dumux/common/geometry/grahamconvexhull.hh b/dumux/common/geometry/grahamconvexhull.hh
index 463b61f90714689ac3069ea1f9132aa165d69e47..9a26d4e44fbfca38abb34703a69bf81926acc923 100644
--- a/dumux/common/geometry/grahamconvexhull.hh
+++ b/dumux/common/geometry/grahamconvexhull.hh
@@ -26,6 +26,7 @@
#include <vector>
#include <array>
#include <algorithm>
+#include <iterator>
#include <dune/common/deprecated.hh>
#include <dune/common/exceptions.hh>
@@ -60,13 +61,13 @@ int getOrientation(const Dune::FieldVector<ctype, 3>& a,
* \ingroup Geometry
* \brief Compute the points making up the convex hull around the given set of unordered points
* \note We assume that all points are coplanar and there are no indentical points in the list
+ * \note This is the overload if we are not allowed to write into the given points vector
*/
-template<class ctype>
-std::vector<Dune::FieldVector<ctype, 3>>
-grahamConvexHull2d3d(const std::vector<Dune::FieldVector<ctype, 3>>& points)
+template<int dim, class Points>
+Points grahamConvexHull(const Points& points)
{
auto copyPoints = points;
- return grahamConvexHull2d3d(copyPoints);
+ return grahamConvexHull<dim>(copyPoints);
}
/*!
@@ -76,9 +77,10 @@ grahamConvexHull2d3d(const std::vector<Dune::FieldVector<ctype, 3>>& points)
* \note This algorithm changes the order of the given points a bit
* as they are unordered anyway this shouldn't matter too much
*/
-template<class ctype>
+template<int dim, class ctype,
+ std::enable_if_t<(dim==2), int> = 0>
std::vector<Dune::FieldVector<ctype, 3>>
-grahamConvexHull2d3d(std::vector<Dune::FieldVector<ctype, 3>>& points)
+grahamConvexHull(std::vector<Dune::FieldVector<ctype, 3>>& points)
{
using Point = Dune::FieldVector<ctype, 3>;
std::vector<Point> convexHull;
@@ -175,6 +177,49 @@ grahamConvexHull2d3d(std::vector<Dune::FieldVector<ctype, 3>>& points)
return convexHull;
}
+/*!
+ * \ingroup Geometry
+ * \brief Compute the points making up the convex hull around the given set of unordered points
+ * \note This is the specialization for 2d space. Here, we make use of the generic implementation
+ * for the case of coplanar points in 3d space (a more efficient implementation could be provided).
+ */
+template<int dim, class ctype,
+ std::enable_if_t<(dim==2), int> = 0>
+std::vector<Dune::FieldVector<ctype, 2>>
+grahamConvexHull(const std::vector<Dune::FieldVector<ctype, 2>>& points)
+{
+ std::vector<Dune::FieldVector<ctype, 3>> points3D;
+ points3D.reserve(points.size());
+ std::transform(points.begin(), points.end(), std::back_inserter(points3D),
+ [](const auto& p) { return Dune::FieldVector<ctype, 3>({p[0], p[1], 0.0}); });
+
+ const auto result3D = grahamConvexHull<2>(points3D);
+
+ std::vector<Dune::FieldVector<ctype, 2>> result2D;
+ result2D.reserve(result3D.size());
+ std::transform(result3D.begin(), result3D.end(), std::back_inserter(result2D),
+ [](const auto& p) { return Dune::FieldVector<ctype, 2>({p[0], p[1]}); });
+
+ return result2D;
+}
+
+// deprecated interfaces
+#ifndef DOXYGEN
+template<class ctype>
+std::vector<Dune::FieldVector<ctype, 3>>
+DUNE_DEPRECATED_MSG("Use grahamConvexHull<dim> with dim as template argument. Will be removed after 3.1")
+grahamConvexHull2d3d(const std::vector<Dune::FieldVector<ctype, 3>>& points)
+{
+ auto copyPoints = points;
+ return grahamConvexHull<2>(copyPoints);
+}
+
+template<class ctype>
+std::vector<Dune::FieldVector<ctype, 3>>
+DUNE_DEPRECATED_MSG("Use grahamConvexHull<dim> with dim as template argument. Will be removed after 3.1")
+grahamConvexHull2d3d(std::vector<Dune::FieldVector<ctype, 3>>& points)
+{ return grahamConvexHull<2>(points); }
+
/*!
* \ingroup Geometry
* \brief Triangulate area given points of the convex hull
@@ -183,9 +228,10 @@ grahamConvexHull2d3d(std::vector<Dune::FieldVector<ctype, 3>>& points)
*/
template<class ctype>
std::vector<std::array<Dune::FieldVector<ctype, 3>, 3> >
-DUNE_DEPRECATED_MSG("Please use triangulate")
+DUNE_DEPRECATED_MSG("Please use triangulate. Will be removed after 3.1")
triangulateConvexHull(const std::vector<Dune::FieldVector<ctype, 3>>& convexHull)
{ return triangulate<2, 3>(convexHull); }
+#endif
} // end namespace Dumux
diff --git a/dumux/common/geometry/triangulation.hh b/dumux/common/geometry/triangulation.hh
index 8bc76f992dcdbbc7231c3a97ebea344cdf00dacd..2b04870cd6e7357f7d9ebc18a1fd2323caac982a 100644
--- a/dumux/common/geometry/triangulation.hh
+++ b/dumux/common/geometry/triangulation.hh
@@ -58,19 +58,19 @@ using Triangulation = std::vector< std::array<Dune::FieldVector<ctype, dimWorld>
* \tparam dimWorld The dimension of the coordinate space
* \tparam Policy Specifies the algorithm to be used for triangulation
*
- * \note this specialization is for 2d triangulations in 3d space using mid point policy
+ * \note this specialization is for 2d triangulations using mid point policy
* \note Assumes all points of the convex hull are coplanar
* \note This inserts a mid point and connects all corners with that point to triangles
*/
template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>,
class RandomAccessContainer,
std::enable_if_t< std::is_same<Policy, TriangulationPolicy::MidPointPolicy>::value
- && dim == 2 && dimWorld == 3, int> = 0 >
+ && dim == 2, int> = 0 >
inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type>
triangulate(const RandomAccessContainer& convexHull)
{
using ctype = typename RandomAccessContainer::value_type::value_type;
- using Point = Dune::FieldVector<ctype, 3>;
+ using Point = Dune::FieldVector<ctype, dimWorld>;
using Triangle = std::array<Point, 3>;
static_assert(std::is_same<typename RandomAccessContainer::value_type, Point>::value,
diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
index d6fb00cb610ac724531ef6311484431b8b0c0f29..28f60e2606f2dd9c02bb7d07fc883b172b730e77 100644
--- a/test/common/geometry/CMakeLists.txt
+++ b/test/common/geometry/CMakeLists.txt
@@ -1,6 +1,8 @@
dune_add_test(SOURCES test_0d1d_intersection.cc LABELS unit)
+dune_add_test(SOURCES test_1d1d_intersection.cc LABELS unit)
dune_add_test(SOURCES test_1d3d_intersection.cc LABELS unit)
dune_add_test(SOURCES test_1d2d_intersection.cc LABELS unit)
+dune_add_test(SOURCES test_2d2d_intersection.cc LABELS unit)
dune_add_test(SOURCES test_2d3d_intersection.cc LABELS unit)
dune_add_test(SOURCES test_graham_convex_hull.cc LABELS unit)
dune_add_test(SOURCES test_makegeometry.cc LABELS unit)
diff --git a/test/common/geometry/test_1d1d_intersection.cc b/test/common/geometry/test_1d1d_intersection.cc
new file mode 100644
index 0000000000000000000000000000000000000000..787f4176d0b99297c6441f6522c9dd8185de0df3
--- /dev/null
+++ b/test/common/geometry/test_1d1d_intersection.cc
@@ -0,0 +1,131 @@
+#include <config.h>
+
+#include <iostream>
+#include <algorithm>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/fvector.hh>
+
+#include <dumux/common/geometry/geometryintersection.hh>
+
+#ifndef DOXYGEN
+template<int dimWorld>
+Dune::FieldVector<double, dimWorld>
+makePosition( double x, double y )
+{
+ if (dimWorld == 2)
+ return {x, y};
+ else if (dimWorld == 3)
+ return {x, y, 1.0};
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Invalid dimworld");
+}
+
+template<int dimWorld>
+Dune::MultiLinearGeometry<double, 1, dimWorld>
+makeLine( Dune::FieldVector<double, dimWorld>&& source,
+ Dune::FieldVector<double, dimWorld>&& target)
+{
+ using CornerStorage = std::vector<Dune::FieldVector<double, dimWorld>>;
+ return { Dune::GeometryTypes::line,
+ CornerStorage({std::move(source), std::move(target)}) };
+}
+
+template<int dimworld>
+bool testPointIntersection(const Dune::MultiLinearGeometry<double, 1, dimworld>& seg1,
+ const Dune::MultiLinearGeometry<double, 1, dimworld>& seg2,
+ bool foundExpected)
+{
+ using Segment = Dune::MultiLinearGeometry<double, 1, dimworld>;
+ using Policy = Dumux::IntersectionPolicy::PointPolicy<double, dimworld>;
+ using Algorithm = Dumux::GeometryIntersection<Segment, Segment, Policy>;
+
+ typename Algorithm::Intersection is;
+ bool found = Algorithm::intersection(seg1, seg2, is);
+
+ if (!found && foundExpected)
+ std::cerr << "Failed detecting intersection with ";
+ else if (found && foundExpected)
+ std::cout << "Found intersection at " << is << " with ";
+ else if (found && !foundExpected)
+ std::cerr << "Found false positive: intersection at " << is << " with ";
+ else if (!found && !foundExpected)
+ std::cout << "No intersection with ";
+
+ std::cout << seg1.corner(0) << " - " << seg1.corner(1) << " / "
+ << seg2.corner(0) << " - " << seg2.corner(1) << std::endl;
+
+ return (found == foundExpected);
+}
+
+template<int dimWorld>
+void testPointIntersections(std::vector<bool>& returns)
+{
+ // test intersection points of segments
+ for (auto scaling : {1.0, 1e3, 1e12, 1e-12})
+ {
+ std::cout << "Test point intersections with scaling " << scaling << std::endl;
+
+ const auto seg1 = makeLine( makePosition<dimWorld>(0.0, 0.0), makePosition<dimWorld>(0.0, scaling*1.0) );
+ const auto seg2 = makeLine( makePosition<dimWorld>(0.0, 0.0), makePosition<dimWorld>(scaling*1.0, 0.0) );
+ const auto seg3 = makeLine( makePosition<dimWorld>(scaling*0.5, 0.0), makePosition<dimWorld>(scaling*1.0, 0.0) );
+ const auto seg4 = makeLine( makePosition<dimWorld>(-1.0*scaling*0.5, scaling*0.5), makePosition<dimWorld>(scaling*1.0, 0.0) );
+ const auto seg5 = makeLine( makePosition<dimWorld>(0.0, scaling*1.0), makePosition<dimWorld>(scaling*1.0, 0.0) );
+ const auto seg6 = makeLine( makePosition<dimWorld>(-1.0*scaling*0.5, scaling*0.5), makePosition<dimWorld>(0.0, scaling*0.5) );
+ const auto seg7 = makeLine( makePosition<dimWorld>(-1.0*scaling*0.5, scaling*0.5), makePosition<dimWorld>(-1.0*scaling*0.1, scaling*0.5) );
+
+ // segments with same orientation touching in source or target or not at all
+ const auto seg8 = makeLine( makePosition<dimWorld>(0.0, 0.0), makePosition<dimWorld>(0.0, -1.0*scaling*1.0) );
+ const auto seg9 = makeLine( makePosition<dimWorld>(0.0, scaling*1.0), makePosition<dimWorld>(0.0, scaling*2.0) );
+ const auto seg10 = makeLine( makePosition<dimWorld>(0.0, -1.0*0.01*scaling), makePosition<dimWorld>(0.0, -1.0*scaling*1.0) );
+ const auto seg11 = makeLine( makePosition<dimWorld>(0.0, scaling*1.01), makePosition<dimWorld>(0.0, scaling*2.0) );
+
+ returns.push_back(testPointIntersection(seg1, seg2, true));
+ returns.push_back(testPointIntersection(seg1, seg3, false));
+ returns.push_back(testPointIntersection(seg1, seg4, true));
+ returns.push_back(testPointIntersection(seg1, seg5, true));
+ returns.push_back(testPointIntersection(seg1, seg6, true));
+ returns.push_back(testPointIntersection(seg1, seg7, false));
+ returns.push_back(testPointIntersection(seg1, seg8, true));
+ returns.push_back(testPointIntersection(seg1, seg9, true));
+ returns.push_back(testPointIntersection(seg1, seg10, false));
+ returns.push_back(testPointIntersection(seg1, seg11, false));
+
+ std::cout << std::endl;
+ }
+}
+
+#endif
+
+int main (int argc, char *argv[]) try
+{
+ // maybe initialize mpi
+ Dune::MPIHelper::instance(argc, argv);
+
+ // collect returns to determine exit code
+ std::vector<bool> returns;
+
+ // test for dimWorld = 2
+ testPointIntersections<2>(returns);
+
+ // TODO: implement and test for dimWorld = 3
+ // testPointIntersections<3>(returns);
+
+ // TODO: implement and test segment intersections
+
+ // determine the exit code
+ if (std::any_of(returns.begin(), returns.end(), [](bool i){ return !i; }))
+ return 1;
+
+ std::cout << "All tests passed!" << std::endl;
+
+ return 0;
+}
+// //////////////////////////////////
+// Error handler
+// /////////////////////////////////
+catch (const Dune::Exception& e) {
+ std::cout << e << std::endl;
+ return 1;
+}
diff --git a/test/common/geometry/test_2d2d_intersection.cc b/test/common/geometry/test_2d2d_intersection.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1d54d848230bd4869bdaeb6e69229a2bb06532fd
--- /dev/null
+++ b/test/common/geometry/test_2d2d_intersection.cc
@@ -0,0 +1,155 @@
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/geometry/multilineargeometry.hh>
+
+#include <dumux/common/geometry/geometryintersection.hh>
+
+#ifndef DOXYGEN
+template<int dimWorld>
+Dune::FieldVector<double, dimWorld>
+makePosition( double x, double y )
+{
+ if (dimWorld == 2)
+ return {x, y};
+ else if (dimWorld == 3)
+ return {x, y, 1.0};
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Invalid dimworld");
+}
+
+template<int dimWorld>
+Dune::MultiLinearGeometry<double, 2, dimWorld>
+makeTriangle( Dune::FieldVector<double, dimWorld>&& a,
+ Dune::FieldVector<double, dimWorld>&& b,
+ Dune::FieldVector<double, dimWorld>&& c )
+{
+ using CornerStorage = std::vector<Dune::FieldVector<double, dimWorld>>;
+ return { Dune::GeometryTypes::triangle,
+ CornerStorage({std::move(a), std::move(b), std::move(c)}) };
+}
+
+template<int dimWorld>
+Dune::MultiLinearGeometry<double, 2, dimWorld>
+makeQuadrilateral( Dune::FieldVector<double, dimWorld>&& a,
+ Dune::FieldVector<double, dimWorld>&& b,
+ Dune::FieldVector<double, dimWorld>&& c,
+ Dune::FieldVector<double, dimWorld>&& d )
+{
+ using CornerStorage = std::vector<Dune::FieldVector<double, dimWorld>>;
+ return { Dune::GeometryTypes::quadrilateral,
+ CornerStorage({std::move(a), std::move(b), std::move(c), std::move(d)}) };
+}
+
+template<int dimworld, class Polygon1, class Polygon2>
+bool testPolygonIntersection(const Polygon1& pol1,
+ const Polygon2& pol2,
+ bool expectIntersection)
+{
+ using Policy = Dumux::IntersectionPolicy::PolygonPolicy<double, dimworld>;
+ using Algorithm = Dumux::GeometryIntersection<Polygon1, Polygon2, Policy>;
+
+ typename Algorithm::IntersectionType intersection;
+ const bool found = Algorithm::intersection(pol1, pol2, intersection);
+ if (found && !expectIntersection)
+ std::cout << "Found false positive!" << std::endl;
+ else if (!found && expectIntersection)
+ std::cout << "Could not find intersection!" << std::endl;
+ else if (found && expectIntersection)
+ std::cout << "Found intersection" << std::endl;
+ else
+ std::cout << "No intersection" << std::endl;
+
+ return (found == expectIntersection);
+}
+
+template<int dimWorld>
+void testPolygonIntersections(std::vector<bool>& returns)
+{
+ for (auto scaling : {1.0, 1e3, 1e12, 1e-12})
+ {
+ std::cout << "Test with scaling " << scaling << std::endl;
+ const auto tria1 = makeTriangle( makePosition<dimWorld>(0.0, 0.0),
+ makePosition<dimWorld>(0.0, 1.0*scaling),
+ makePosition<dimWorld>(1.0*scaling, 1.0*scaling) );
+ const auto tria2 = makeTriangle( makePosition<dimWorld>(1.0*scaling, 0.0),
+ makePosition<dimWorld>(0.0, 0.0),
+ makePosition<dimWorld>(0.0, 1.0*scaling) );
+ const auto tria3 = makeTriangle( makePosition<dimWorld>(1.0*scaling, 0.0),
+ makePosition<dimWorld>(0.0, 0.5*scaling),
+ makePosition<dimWorld>(1.0*scaling, 1.0*scaling) );
+ const auto tria4 = makeTriangle( makePosition<dimWorld>(1.0*scaling, 0.0),
+ makePosition<dimWorld>(1.0*scaling, 1.0*scaling),
+ makePosition<dimWorld>(2.0*scaling, 1.0*scaling) );
+ const auto tria5 = makeTriangle( makePosition<dimWorld>(1.0*scaling, 0.0),
+ makePosition<dimWorld>(0.5*scaling, 0.5*scaling),
+ makePosition<dimWorld>(2.0*scaling, 1.0*scaling) );
+ const auto tria6 = makeTriangle( makePosition<dimWorld>(1.0*scaling, 0.0),
+ makePosition<dimWorld>(0.501*scaling, 0.501*scaling),
+ makePosition<dimWorld>(2.0*scaling, 1.0*scaling) );
+ const auto tria7 = makeTriangle( makePosition<dimWorld>(1.0*scaling, 0.0),
+ makePosition<dimWorld>(0.499*scaling, 0.501*scaling),
+ makePosition<dimWorld>(2.0*scaling, 1.0*scaling) );
+
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, tria2, true));
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, tria3, true));
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, tria4, false)); // touches in point
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, tria5, false)); // touches in edge
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, tria6, false)); // little bit outside triangle
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, tria7, true)); // little bit inside triangle
+
+ // test some quadrilaterals
+ const auto quad1 = makeQuadrilateral( makePosition<dimWorld>(0.0, 0.0),
+ makePosition<dimWorld>(1.0*scaling, 0.0),
+ makePosition<dimWorld>(0.0, 1.0*scaling),
+ makePosition<dimWorld>(1.0*scaling, 1.0*scaling) );
+ const auto quad2 = makeQuadrilateral( makePosition<dimWorld>(1.0*scaling, 0.0),
+ makePosition<dimWorld>(2.0*scaling, 0.0),
+ makePosition<dimWorld>(1.0*scaling, 2.0*scaling),
+ makePosition<dimWorld>(2.0*scaling, 2.0*scaling) );
+ const auto quad3 = makeQuadrilateral( makePosition<dimWorld>(-1.0*scaling, 0.0),
+ makePosition<dimWorld>(0.0, 0.0),
+ makePosition<dimWorld>(-1.0*scaling, 1.0*scaling),
+ makePosition<dimWorld>(0.0, 1.0*scaling) );
+ const auto quad4 = makeQuadrilateral( makePosition<dimWorld>(0.5*scaling, 0.0),
+ makePosition<dimWorld>(0.5*scaling, 0.501*scaling),
+ makePosition<dimWorld>(1.0*scaling, 0.0),
+ makePosition<dimWorld>(1.0*scaling, 1.0*scaling) );
+
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, quad1, true));
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, quad2, false)); // touches in point
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, quad3, false)); // touches in edge
+ returns.push_back(testPolygonIntersection<dimWorld>(tria1, quad4, true)); // small overlap region
+
+ std::cout << std::endl;
+ }
+}
+
+#endif
+
+int main(int argc, char* argv[]) try
+{
+ std::vector<bool> returns;
+ testPolygonIntersections<2>(returns);
+
+ // TODO: implement and test intersections in 3d
+ // TODO: implement and test point and segment intersections
+
+ // determine the exit code
+ if (std::any_of(returns.begin(), returns.end(), [](bool i){ return !i; }))
+ return 1;
+
+ std::cout << "All tests passed!" << std::endl;
+
+ return 0;
+}
+// //////////////////////////////////
+// Error handler
+// /////////////////////////////////
+catch (const Dune::Exception& e) {
+ std::cout << e << std::endl;
+ return 1;
+}
diff --git a/test/common/geometry/test_graham_convex_hull.cc b/test/common/geometry/test_graham_convex_hull.cc
index 5a663e4e41c51845c9fe60725af7d08f2d1186fe..a59029cc12325393b8006727d82baa7f400cedbb 100644
--- a/test/common/geometry/test_graham_convex_hull.cc
+++ b/test/common/geometry/test_graham_convex_hull.cc
@@ -97,6 +97,8 @@ void writeVTKPolyData(const std::vector<Dune::FieldVector<double, 3>>& p,
int main(int argc, char* argv[]) try
{
+ using namespace Dumux;
+
using Point = Dune::FieldVector<double, 3>;
// create 100 random points
std::vector<Point> points(100);
@@ -109,30 +111,30 @@ int main(int argc, char* argv[]) try
writeCSV(points, "points");
Dune::Timer timer;
- auto convexHullPoints = Dumux::grahamConvexHull2d3d(points);
+ auto convexHullPoints = grahamConvexHull<2>(points);
std::cout << "Computed convex hull of " << points.size() << " points in "
<< timer.elapsed() << " seconds." << std::endl;
writeVTKPolyData(convexHullPoints, "convexhull");
Dune::Timer timer2;
- auto triangles = Dumux::triangulate<2, 3>(convexHullPoints);
+ auto triangles = triangulate<2, 3>(convexHullPoints);
std::cout << "Computed triangulation of convex hull with " << convexHullPoints.size()
<< " points in " << timer2.elapsed() << " seconds." << std::endl;
- Dumux::writeVTKPolyDataTriangle(triangles, "triangulation");
+ writeVTKPolyDataTriangle(triangles, "triangulation");
// some specific tests for corner cases with colinear points
points = {{0.0,0.0,0.0}, {0.0,0.0,1.0}, {0.0,0.0,2.0}, {0.0,0.0,3.0}};
- auto hull = Dumux::grahamConvexHull2d3d(points);
+ auto hull = grahamConvexHull<2>(points);
if (!(hull.empty())) DUNE_THROW(Dune::InvalidStateException, "False positive for finding a convex hull!");
points = {{0.0,0.0,0.0}, {0.0,0.0,1.0}, {0.0,0.0,2.0}, {0.0,0.0,3.0}, {0.0,0.0,4.0}, {2.0,3.0,3.0}};
- hull = Dumux::grahamConvexHull2d3d(points);
+ hull = grahamConvexHull<2>(points);
if (hull.empty()) DUNE_THROW(Dune::InvalidStateException, "Didn't find convex hull!");
points = {{2.0,3.0,3.0}, {0.0,0.0,4.0}, {0.0,0.0,0.0}, {0.0,0.0,1.0}, {0.0,0.0,2.0}, {0.0,0.0,3.0}};
- hull = Dumux::grahamConvexHull2d3d(points);
+ hull = grahamConvexHull<2>(points);
if (hull.empty()) DUNE_THROW(Dune::InvalidStateException, "Didn't find convex hull!");
return 0;
diff --git a/test/common/geometry/writetriangulation.hh b/test/common/geometry/writetriangulation.hh
index c106e285b76b3e2bb17429ef78b4141be86146dc..5b8bfa5d02bc72f0363c28327830c5c16299857e 100644
--- a/test/common/geometry/writetriangulation.hh
+++ b/test/common/geometry/writetriangulation.hh
@@ -42,8 +42,17 @@ void writeVTKPolyDataTriangle(const TriangleVector& triangles,
<< " <DataArray type=\"Float32\" Name=\"Coordinates\" NumberOfComponents=\"3\" format=\"ascii\">\n";
for (const auto& t : triangles)
+ {
for (const auto& p : t)
+ {
fout << p << " ";
+ if (p.size() == 1)
+ fout << "0.0 0.0 ";
+ else if (p.size() == 2)
+ fout << "0.0 ";
+ }
+
+ }
fout << '\n';
fout << " </DataArray>\n"