Commit 6a00eca4 authored by Timo Koch's avatar Timo Koch
Browse files

[geometry] Add template param dim to grahamConvexHull

Rename grahamConvexHull2d3d to grahamConvexHull<2>,
deprecate old interface.
parent ac5baf27
......@@ -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;
......
......@@ -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
......
......@@ -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;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment