From 6a00eca4f9e9ac217965193a44d0baa2c5fb6b5f Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 31 May 2019 15:37:39 +0200
Subject: [PATCH] [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 <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;
@@ -181,25 +183,43 @@ grahamConvexHull2d3d(std::vector<Dune::FieldVector<ctype, 3>>& 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<class ctype>
+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>> tmp;
- tmp.reserve(points.size());
- for (const auto& p : points)
- tmp.emplace_back( Dune::FieldVector<ctype, 3>({p[0], p[1], 0.0}) );
+ 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}); });
- auto result3d = grahamConvexHull2d3d(tmp);
+ const auto result3D = grahamConvexHull<2>(points3D);
- std::vector<Dune::FieldVector<ctype, 2>> result;
- result.reserve(result3d.size());
- for (const auto& p : result3d)
- result.emplace_back( Dune::FieldVector<ctype, 2>({p[0], p[1]}) );
+ 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 result;
+ 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
@@ -208,9 +228,10 @@ grahamConvexHull(const std::vector<Dune::FieldVector<ctype, 2>>& 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/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<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;
--
GitLab