triangulation.hh

Go to the documentation of this file.

/*****************************************************************************
 *   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/>.   *
 *****************************************************************************/
#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 {
 
struct MidPointPolicy {};
 
struct DelaunayPolicy {};
 
template<int dim, int dimWorld>
using DefaultPolicy = std::conditional_t< dim >= 2, MidPointPolicy, DelaunayPolicy >;
 
} // end namespace TriangulationPolicy
 
template<int dim, int dimWorld, class ctype>
using Triangulation = std::vector< std::array<Dune::FieldVector<ctype, dimWorld>, dim+1> >;
 
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, 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, dimWorld>;
    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;
}
 
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, dimWorld>;
 
    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]} });
 
    assert(points.size() > 1);
    DUNE_THROW(Dune::NotImplemented, "1d triangulation for point cloud size > 2");
}
 
} // end namespace Dumux
 
# endif