From a5200767c9d9e4f27305ef5b8eb212e4cd78f47f Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 30 May 2019 12:32:30 +0200 Subject: [PATCH 01/10] [geometryintersection] add algo for 1d-1d intersections in 2d space --- dumux/common/geometry/geometryintersection.hh | 121 ++++++++++++++++++ 1 file changed, 121 insertions(+) diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh index e3199636d7..91588c948b 100644 --- a/dumux/common/geometry/geometryintersection.hh +++ b/dumux/common/geometry/geometryintersection.hh @@ -32,6 +32,7 @@ #include #include #include +#include namespace Dumux { namespace IntersectionPolicy { @@ -214,6 +215,126 @@ public: } }; +/*! + * \ingroup Geometry + * \brief A class for segment--segment intersection in 2d space + */ +template +class GeometryIntersection +{ + 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 = 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 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 = 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 -- GitLab From abfa13bd4ba44c0042db8057a3fb0d3b617ac471 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 30 May 2019 12:32:52 +0200 Subject: [PATCH 02/10] [geometryintersection] improve docu --- dumux/common/geometry/geometryintersection.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh index 91588c948b..d7c1336abf 100644 --- a/dumux/common/geometry/geometryintersection.hh +++ b/dumux/common/geometry/geometryintersection.hh @@ -672,7 +672,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 = 0> @@ -796,7 +796,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 = 0> -- GitLab From c4527cb1fd69069487a3aadf39fe47e861d2446c Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 30 May 2019 12:33:42 +0200 Subject: [PATCH 03/10] [geometryintersection] add test for 1d-1d intersections --- test/common/geometry/CMakeLists.txt | 1 + .../common/geometry/test_1d1d_intersection.cc | 131 ++++++++++++++++++ 2 files changed, 132 insertions(+) create mode 100644 test/common/geometry/test_1d1d_intersection.cc diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt index d6fb00cb61..6efe5cd972 100644 --- a/test/common/geometry/CMakeLists.txt +++ b/test/common/geometry/CMakeLists.txt @@ -1,4 +1,5 @@ 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_2d3d_intersection.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 0000000000..787f4176d0 --- /dev/null +++ b/test/common/geometry/test_1d1d_intersection.cc @@ -0,0 +1,131 @@ +#include + +#include +#include + +#include +#include +#include + +#include + +#ifndef DOXYGEN +template +Dune::FieldVector +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 +Dune::MultiLinearGeometry +makeLine( Dune::FieldVector&& source, + Dune::FieldVector&& target) +{ + using CornerStorage = std::vector>; + return { Dune::GeometryTypes::line, + CornerStorage({std::move(source), std::move(target)}) }; +} + +template +bool testPointIntersection(const Dune::MultiLinearGeometry& seg1, + const Dune::MultiLinearGeometry& seg2, + bool foundExpected) +{ + using Segment = Dune::MultiLinearGeometry; + using Policy = Dumux::IntersectionPolicy::PointPolicy; + using Algorithm = Dumux::GeometryIntersection; + + 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 +void testPointIntersections(std::vector& 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(0.0, 0.0), makePosition(0.0, scaling*1.0) ); + const auto seg2 = makeLine( makePosition(0.0, 0.0), makePosition(scaling*1.0, 0.0) ); + const auto seg3 = makeLine( makePosition(scaling*0.5, 0.0), makePosition(scaling*1.0, 0.0) ); + const auto seg4 = makeLine( makePosition(-1.0*scaling*0.5, scaling*0.5), makePosition(scaling*1.0, 0.0) ); + const auto seg5 = makeLine( makePosition(0.0, scaling*1.0), makePosition(scaling*1.0, 0.0) ); + const auto seg6 = makeLine( makePosition(-1.0*scaling*0.5, scaling*0.5), makePosition(0.0, scaling*0.5) ); + const auto seg7 = makeLine( makePosition(-1.0*scaling*0.5, scaling*0.5), makePosition(-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(0.0, 0.0), makePosition(0.0, -1.0*scaling*1.0) ); + const auto seg9 = makeLine( makePosition(0.0, scaling*1.0), makePosition(0.0, scaling*2.0) ); + const auto seg10 = makeLine( makePosition(0.0, -1.0*0.01*scaling), makePosition(0.0, -1.0*scaling*1.0) ); + const auto seg11 = makeLine( makePosition(0.0, scaling*1.01), makePosition(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 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; +} -- GitLab From 540f09650c410cb6991f9a35ec7157a9805b935a Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 30 May 2019 22:04:24 +0200 Subject: [PATCH 04/10] [grahamconvexhull] add impl for 2d hulls This actually makes use of the 2d in 3d algorithm at the moment. In the future a more efficient implementation for 2d hulls could be implemented. --- dumux/common/geometry/grahamconvexhull.hh | 25 +++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/dumux/common/geometry/grahamconvexhull.hh b/dumux/common/geometry/grahamconvexhull.hh index 463b61f907..9fff6845c7 100644 --- a/dumux/common/geometry/grahamconvexhull.hh +++ b/dumux/common/geometry/grahamconvexhull.hh @@ -175,6 +175,31 @@ grahamConvexHull2d3d(std::vector>& 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 +std::vector> +grahamConvexHull(const std::vector>& points) +{ + std::vector> tmp; + tmp.reserve(points.size()); + for (const auto& p : points) + tmp.emplace_back( Dune::FieldVector({p[0], p[1], 0.0}) ); + + auto result3d = grahamConvexHull2d3d(tmp); + + std::vector> result; + result.reserve(result3d.size()); + for (const auto& p : result3d) + result.emplace_back( Dune::FieldVector({p[0], p[1]}) ); + + return result; +} + /*! * \ingroup Geometry * \brief Triangulate area given points of the convex hull -- GitLab From e29c59990e1df5018bb49f26a319161c4d997f58 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 30 May 2019 22:09:12 +0200 Subject: [PATCH 05/10] [triangulation] avoid algorithm being specific to dimworld=3 --- dumux/common/geometry/triangulation.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dumux/common/geometry/triangulation.hh b/dumux/common/geometry/triangulation.hh index 8bc76f992d..2b04870cd6 100644 --- a/dumux/common/geometry/triangulation.hh +++ b/dumux/common/geometry/triangulation.hh @@ -58,19 +58,19 @@ using Triangulation = std::vector< std::array * \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, class RandomAccessContainer, std::enable_if_t< std::is_same::value - && dim == 2 && dimWorld == 3, int> = 0 > + && dim == 2, int> = 0 > inline Triangulation triangulate(const RandomAccessContainer& convexHull) { using ctype = typename RandomAccessContainer::value_type::value_type; - using Point = Dune::FieldVector; + using Point = Dune::FieldVector; using Triangle = std::array; static_assert(std::is_same::value, -- GitLab From 31db57ef2acdca056cc27db3fdf785e51448e9bf Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 30 May 2019 22:10:03 +0200 Subject: [PATCH 06/10] [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 GeometryIntersection +{ + 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; + using ReferenceElementsGeo2 = typename Dune::ReferenceElements; + +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 = 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 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; + using PointPolicy = IntersectionPolicy::PointPolicy; + for (int i = 0; i < referenceElement1.size(dim1-1); ++i) + { + const auto localEdgeGeom1 = referenceElement1.template geometry(i); + const auto edge1 = SegGeometry( Dune::GeometryTypes::line, + std::vector( {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(j); + const auto edge2 = SegGeometry( Dune::GeometryTypes::line, + std::vector( {geo2.global(localEdgeGeom2.corner(0)), + geo2.global(localEdgeGeom2.corner(1))} )); + + using EdgeTest = GeometryIntersection; + 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 = 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 = 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 From ac5baf27cf3d3d3a4fb0c3f1e4ddd9643e85e8b8 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Thu, 30 May 2019 22:10:58 +0200 Subject: [PATCH 07/10] [geointersection] add test for 2d-2d intersections --- test/common/geometry/CMakeLists.txt | 1 + .../common/geometry/test_2d2d_intersection.cc | 155 ++++++++++++++++++ 2 files changed, 156 insertions(+) create mode 100644 test/common/geometry/test_2d2d_intersection.cc diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt index 6efe5cd972..28f60e2606 100644 --- a/test/common/geometry/CMakeLists.txt +++ b/test/common/geometry/CMakeLists.txt @@ -2,6 +2,7 @@ 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_2d2d_intersection.cc b/test/common/geometry/test_2d2d_intersection.cc new file mode 100644 index 0000000000..1d54d84823 --- /dev/null +++ b/test/common/geometry/test_2d2d_intersection.cc @@ -0,0 +1,155 @@ +#include + +#include + +#include +#include +#include + +#include + +#ifndef DOXYGEN +template +Dune::FieldVector +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 +Dune::MultiLinearGeometry +makeTriangle( Dune::FieldVector&& a, + Dune::FieldVector&& b, + Dune::FieldVector&& c ) +{ + using CornerStorage = std::vector>; + return { Dune::GeometryTypes::triangle, + CornerStorage({std::move(a), std::move(b), std::move(c)}) }; +} + +template +Dune::MultiLinearGeometry +makeQuadrilateral( Dune::FieldVector&& a, + Dune::FieldVector&& b, + Dune::FieldVector&& c, + Dune::FieldVector&& d ) +{ + using CornerStorage = std::vector>; + return { Dune::GeometryTypes::quadrilateral, + CornerStorage({std::move(a), std::move(b), std::move(c), std::move(d)}) }; +} + +template +bool testPolygonIntersection(const Polygon1& pol1, + const Polygon2& pol2, + bool expectIntersection) +{ + using Policy = Dumux::IntersectionPolicy::PolygonPolicy; + using Algorithm = Dumux::GeometryIntersection; + + 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 +void testPolygonIntersections(std::vector& returns) +{ + for (auto scaling : {1.0, 1e3, 1e12, 1e-12}) + { + std::cout << "Test with scaling " << scaling << std::endl; + const auto tria1 = makeTriangle( makePosition(0.0, 0.0), + makePosition(0.0, 1.0*scaling), + makePosition(1.0*scaling, 1.0*scaling) ); + const auto tria2 = makeTriangle( makePosition(1.0*scaling, 0.0), + makePosition(0.0, 0.0), + makePosition(0.0, 1.0*scaling) ); + const auto tria3 = makeTriangle( makePosition(1.0*scaling, 0.0), + makePosition(0.0, 0.5*scaling), + makePosition(1.0*scaling, 1.0*scaling) ); + const auto tria4 = makeTriangle( makePosition(1.0*scaling, 0.0), + makePosition(1.0*scaling, 1.0*scaling), + makePosition(2.0*scaling, 1.0*scaling) ); + const auto tria5 = makeTriangle( makePosition(1.0*scaling, 0.0), + makePosition(0.5*scaling, 0.5*scaling), + makePosition(2.0*scaling, 1.0*scaling) ); + const auto tria6 = makeTriangle( makePosition(1.0*scaling, 0.0), + makePosition(0.501*scaling, 0.501*scaling), + makePosition(2.0*scaling, 1.0*scaling) ); + const auto tria7 = makeTriangle( makePosition(1.0*scaling, 0.0), + makePosition(0.499*scaling, 0.501*scaling), + makePosition(2.0*scaling, 1.0*scaling) ); + + returns.push_back(testPolygonIntersection(tria1, tria2, true)); + returns.push_back(testPolygonIntersection(tria1, tria3, true)); + returns.push_back(testPolygonIntersection(tria1, tria4, false)); // touches in point + returns.push_back(testPolygonIntersection(tria1, tria5, false)); // touches in edge + returns.push_back(testPolygonIntersection(tria1, tria6, false)); // little bit outside triangle + returns.push_back(testPolygonIntersection(tria1, tria7, true)); // little bit inside triangle + + // test some quadrilaterals + const auto quad1 = makeQuadrilateral( makePosition(0.0, 0.0), + makePosition(1.0*scaling, 0.0), + makePosition(0.0, 1.0*scaling), + makePosition(1.0*scaling, 1.0*scaling) ); + const auto quad2 = makeQuadrilateral( makePosition(1.0*scaling, 0.0), + makePosition(2.0*scaling, 0.0), + makePosition(1.0*scaling, 2.0*scaling), + makePosition(2.0*scaling, 2.0*scaling) ); + const auto quad3 = makeQuadrilateral( makePosition(-1.0*scaling, 0.0), + makePosition(0.0, 0.0), + makePosition(-1.0*scaling, 1.0*scaling), + makePosition(0.0, 1.0*scaling) ); + const auto quad4 = makeQuadrilateral( makePosition(0.5*scaling, 0.0), + makePosition(0.5*scaling, 0.501*scaling), + makePosition(1.0*scaling, 0.0), + makePosition(1.0*scaling, 1.0*scaling) ); + + returns.push_back(testPolygonIntersection(tria1, quad1, true)); + returns.push_back(testPolygonIntersection(tria1, quad2, false)); // touches in point + returns.push_back(testPolygonIntersection(tria1, quad3, false)); // touches in edge + returns.push_back(testPolygonIntersection(tria1, quad4, true)); // small overlap region + + std::cout << std::endl; + } +} + +#endif + +int main(int argc, char* argv[]) try +{ + std::vector 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; +} -- GitLab From 6a00eca4f9e9ac217965193a44d0baa2c5fb6b5f Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Fri, 31 May 2019 15:37:39 +0200 Subject: [PATCH 08/10] [geometry] Add template param dim to grahamConvexHull Rename grahamConvexHull2d3d to grahamConvexHull<2>, deprecate old interface. --- dumux/common/geometry/geometryintersection.hh | 4 +- dumux/common/geometry/grahamconvexhull.hh | 57 +++++++++++++------ .../geometry/test_graham_convex_hull.cc | 14 +++-- 3 files changed, 49 insertions(+), 26 deletions(-) diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh index 2502ec1a6d..5b00e12154 100644 --- a/dumux/common/geometry/geometryintersection.hh +++ b/dumux/common/geometry/geometryintersection.hh @@ -592,7 +592,7 @@ public: return false; // intersection polygon is convex hull of above points - intersection = grahamConvexHull(points); + intersection = grahamConvexHull<2>(points); assert(!intersection.empty()); return true; } @@ -921,7 +921,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; diff --git a/dumux/common/geometry/grahamconvexhull.hh b/dumux/common/geometry/grahamconvexhull.hh index 9fff6845c7..9a26d4e44f 100644 --- a/dumux/common/geometry/grahamconvexhull.hh +++ b/dumux/common/geometry/grahamconvexhull.hh @@ -26,6 +26,7 @@ #include #include #include +#include #include #include @@ -60,13 +61,13 @@ int getOrientation(const Dune::FieldVector& 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 -std::vector> -grahamConvexHull2d3d(const std::vector>& points) +template +Points grahamConvexHull(const Points& points) { auto copyPoints = points; - return grahamConvexHull2d3d(copyPoints); + return grahamConvexHull(copyPoints); } /*! @@ -76,9 +77,10 @@ grahamConvexHull2d3d(const std::vector>& 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 +template = 0> std::vector> -grahamConvexHull2d3d(std::vector>& points) +grahamConvexHull(std::vector>& points) { using Point = Dune::FieldVector; std::vector convexHull; @@ -181,25 +183,43 @@ grahamConvexHull2d3d(std::vector>& 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 +template = 0> std::vector> grahamConvexHull(const std::vector>& points) { - std::vector> tmp; - tmp.reserve(points.size()); - for (const auto& p : points) - tmp.emplace_back( Dune::FieldVector({p[0], p[1], 0.0}) ); + std::vector> points3D; + points3D.reserve(points.size()); + std::transform(points.begin(), points.end(), std::back_inserter(points3D), + [](const auto& p) { return Dune::FieldVector({p[0], p[1], 0.0}); }); - auto result3d = grahamConvexHull2d3d(tmp); + const auto result3D = grahamConvexHull<2>(points3D); - std::vector> result; - result.reserve(result3d.size()); - for (const auto& p : result3d) - result.emplace_back( Dune::FieldVector({p[0], p[1]}) ); + std::vector> result2D; + result2D.reserve(result3D.size()); + std::transform(result3D.begin(), result3D.end(), std::back_inserter(result2D), + [](const auto& p) { return Dune::FieldVector({p[0], p[1]}); }); - return result; + return result2D; } +// deprecated interfaces +#ifndef DOXYGEN +template +std::vector> +DUNE_DEPRECATED_MSG("Use grahamConvexHull with dim as template argument. Will be removed after 3.1") +grahamConvexHull2d3d(const std::vector>& points) +{ + auto copyPoints = points; + return grahamConvexHull<2>(copyPoints); +} + +template +std::vector> +DUNE_DEPRECATED_MSG("Use grahamConvexHull with dim as template argument. Will be removed after 3.1") +grahamConvexHull2d3d(std::vector>& points) +{ return grahamConvexHull<2>(points); } + /*! * \ingroup Geometry * \brief Triangulate area given points of the convex hull @@ -208,9 +228,10 @@ grahamConvexHull(const std::vector>& points) */ template std::vector, 3> > -DUNE_DEPRECATED_MSG("Please use triangulate") +DUNE_DEPRECATED_MSG("Please use triangulate. Will be removed after 3.1") triangulateConvexHull(const std::vector>& convexHull) { return triangulate<2, 3>(convexHull); } +#endif } // end namespace Dumux diff --git a/test/common/geometry/test_graham_convex_hull.cc b/test/common/geometry/test_graham_convex_hull.cc index 5a663e4e41..a59029cc12 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>& p, int main(int argc, char* argv[]) try { + using namespace Dumux; + using Point = Dune::FieldVector; // create 100 random points std::vector 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; -- GitLab From 38db0204831757f898a025c80df4ac44729cb66e Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Mon, 3 Jun 2019 17:20:24 +0200 Subject: [PATCH 09/10] [geometryisection] improve 2d-2d algo --- dumux/common/geometry/geometryintersection.hh | 65 ++++++++++--------- 1 file changed, 36 insertions(+), 29 deletions(-) diff --git a/dumux/common/geometry/geometryintersection.hh b/dumux/common/geometry/geometryintersection.hh index 5b00e12154..69b5c894d3 100644 --- a/dumux/common/geometry/geometryintersection.hh +++ b/dumux/common/geometry/geometryintersection.hh @@ -534,38 +534,45 @@ public: 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; - using PointPolicy = IntersectionPolicy::PointPolicy; - for (int i = 0; i < referenceElement1.size(dim1-1); ++i) + const auto numPoints1 = points.size(); + if (numPoints1 != geo1.corners()) { - const auto localEdgeGeom1 = referenceElement1.template geometry(i); - const auto edge1 = SegGeometry( Dune::GeometryTypes::line, - std::vector( {geo1.global(localEdgeGeom1.corner(0)), - geo1.global(localEdgeGeom1.corner(1))} )); + // 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)); - for (int j = 0; j < referenceElement2.size(dim2-1); ++j) + if (points.empty()) + return false; + + if (points.size() - numPoints1 != geo2.corners()) { - const auto localEdgeGeom2 = referenceElement2.template geometry(j); - const auto edge2 = SegGeometry( Dune::GeometryTypes::line, - std::vector( {geo2.global(localEdgeGeom2.corner(0)), - geo2.global(localEdgeGeom2.corner(1))} )); - - using EdgeTest = GeometryIntersection; - typename EdgeTest::Intersection edgeIntersection; - if (EdgeTest::intersection(edge1, edge2, edgeIntersection)) - points.emplace_back(edgeIntersection); + const auto referenceElement1 = ReferenceElementsGeo1::general(geo1.type()); + const auto referenceElement2 = ReferenceElementsGeo2::general(geo2.type()); + + // add intersections of edges + using SegGeometry = Dune::MultiLinearGeometry; + using PointPolicy = IntersectionPolicy::PointPolicy; + for (int i = 0; i < referenceElement1.size(dim1-1); ++i) + { + const auto localEdgeGeom1 = referenceElement1.template geometry(i); + const auto edge1 = SegGeometry( Dune::GeometryTypes::line, + std::vector( {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(j); + const auto edge2 = SegGeometry( Dune::GeometryTypes::line, + std::vector( {geo2.global(localEdgeGeom2.corner(0)), + geo2.global(localEdgeGeom2.corner(1))} )); + + using EdgeTest = GeometryIntersection; + typename EdgeTest::Intersection edgeIntersection; + if (EdgeTest::intersection(edge1, edge2, edgeIntersection)) + points.emplace_back(edgeIntersection); + } + } } } -- GitLab From a46222afce3b1a25cfaf3509985810c90edb5778 Mon Sep 17 00:00:00 2001 From: "Dennis.Glaeser" Date: Mon, 3 Jun 2019 17:22:10 +0200 Subject: [PATCH 10/10] [test][geom] fix triangulation vtk output for 2d grids --- test/common/geometry/writetriangulation.hh | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/test/common/geometry/writetriangulation.hh b/test/common/geometry/writetriangulation.hh index c106e285b7..5b8bfa5d02 100644 --- a/test/common/geometry/writetriangulation.hh +++ b/test/common/geometry/writetriangulation.hh @@ -42,8 +42,17 @@ void writeVTKPolyDataTriangle(const TriangleVector& triangles, << " \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 << " \n" -- GitLab