diff --git a/dumux/common/geometry/CMakeLists.txt b/dumux/common/geometry/CMakeLists.txt index 5d6fc9ee409ce63a66e8b755d9508369a87b67ca..fd5972a41baebb36223d8a7cb7fbdd7756423c19 100644 --- a/dumux/common/geometry/CMakeLists.txt +++ b/dumux/common/geometry/CMakeLists.txt @@ -9,4 +9,5 @@ intersectingentities.hh intersectspointgeometry.hh intersectspointsimplex.hh makegeometry.hh +triangulation.hh DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/common/geometry) diff --git a/dumux/common/geometry/triangulation.hh b/dumux/common/geometry/triangulation.hh new file mode 100644 index 0000000000000000000000000000000000000000..bfb6cde114e5786ac8c710b80f1089919652e389 --- /dev/null +++ b/dumux/common/geometry/triangulation.hh @@ -0,0 +1,134 @@ +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup Geometry + * \brief Functionality to triangulate point clouds + */ +#ifndef DUMUX_TRIANGULATION_HH +#define DUMUX_TRIANGULATION_HH + +#include <vector> +#include <array> +#include <algorithm> +#include <type_traits> + +#include <dune/common/exceptions.hh> +#include <dune/common/fvector.hh> + +#include <dumux/common/math.hh> + +namespace Dumux { +namespace TriangulationPolicy { + +//! Policy that expects a point cloud that represents a convex +//! hull. Inserts the mid point and connects it to the points of the hull. +struct MidPointPolicy {}; + +//! Delaunay-type triangulations. +struct DelaunayPolicy {}; + +template<int dim, int dimWorld> +using DefaultPolicy = std::conditional_t< dim >= 2, MidPointPolicy, DelaunayPolicy >; + +} // end namespace TriangulationPolicy + +//! The data type to store triangulations +template<int dim, int dimWorld, class ctype> +using Triangulation = std::vector< std::array<Dune::FieldVector<ctype, dimWorld>, dim+1> >; + +/*! + * \ingroup Geometry + * \brief Triangulate area given points of a convex hull + * \tparam dim Specifies the dimension of the resulting triangulation (2 or 3) + * \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 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 > +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 Triangle = std::array<Point, 3>; + + static_assert(std::is_same<typename RandomAccessContainer::value_type, Point>::value, + "Triangulation expects Dune::FieldVector as point type"); + + if (convexHull.size() < 3) + DUNE_THROW(Dune::InvalidStateException, "Try to triangulate point cloud with less than 3 points!"); + + if (convexHull.size() == 3) + return std::vector<Triangle>(1, {convexHull[0], convexHull[1], convexHull[2]}); + + Point midPoint(0.0); + for (const auto p : convexHull) + midPoint += p; + midPoint /= convexHull.size(); + + std::vector<Triangle> triangulation; + triangulation.reserve(convexHull.size()); + + for (std::size_t i = 0; i < convexHull.size()-1; ++i) + triangulation.emplace_back(Triangle{midPoint, convexHull[i], convexHull[i+1]}); + + triangulation.emplace_back(Triangle{midPoint, convexHull[convexHull.size()-1], convexHull[0]}); + + return triangulation; +} + +/*! + * \ingroup Geometry + * \brief Triangulate area given points of a convex hull + * \tparam dim Specifies the dimension of the resulting triangulation (2 or 3) + * \tparam dimWorld The dimension of the coordinate space + * \tparam Policy Specifies the algorithm to be used for triangulation + * + * \note this specialization is for 1d discretization using segments + * \note sorts the points and connects them via segments + */ +template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>, + class RandomAccessContainer, + std::enable_if_t< std::is_same<Policy, TriangulationPolicy::DelaunayPolicy>::value + && dim == 1, int> = 0 > +inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type> +triangulate(const RandomAccessContainer& points) +{ + using ctype = typename RandomAccessContainer::value_type::value_type; + using Point = Dune::FieldVector<ctype, 3>; + + static_assert(std::is_same<typename RandomAccessContainer::value_type, Point>::value, + "Triangulation expects Dune::FieldVector as point type"); + + if (points.size() == 2) + return Triangulation<dim, dimWorld, ctype>({points[0], points[1]}); + + //! \todo sort points and create polyline + assert(points.size() > 1); + DUNE_THROW(Dune::NotImplemented, "1d triangulation for point cloud size > 2"); +} + +} // end namespace Dumux + +# endif