triangulation.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_GEOMETRY_TRIANGULATION_HH
#define DUMUX_GEOMETRY_TRIANGULATION_HH
 
#include <vector>
#include <array>
#include <algorithm>
#include <type_traits>
#include <tuple>
 
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
 
#include <dumux/common/math.hh>
#include <dumux/geometry/intersectspointsimplex.hh>
#include <dumux/geometry/grahamconvexhull.hh>
 
namespace Dumux {
namespace TriangulationPolicy {

struct MidPointPolicy {};

struct ConvexHullPolicy {};

struct DelaunayPolicy {};
 
#ifndef DOXYGEN
namespace Detail {
using DefaultDimPolicies = std::tuple<DelaunayPolicy, MidPointPolicy, ConvexHullPolicy>;
} // end namespace Detail
#endif

template<int dim, int dimWorld>
using DefaultPolicy = std::tuple_element_t<dim-1, Detail::DefaultDimPolicies>;
 
} // 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_v<Policy, TriangulationPolicy::DelaunayPolicy>
                            && 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_v<typename RandomAccessContainer::value_type, Point>,
                  "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");
}

template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>,
          class RandomAccessContainer,
          std::enable_if_t< std::is_same_v<Policy, TriangulationPolicy::MidPointPolicy>
                            && dim == 2, int> = 0 >
inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type>
triangulate(const RandomAccessContainer& convexHullPoints)
{
    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_v<typename RandomAccessContainer::value_type, Point>,
                  "Triangulation expects Dune::FieldVector as point type");
 
    if (convexHullPoints.size() < 3)
        DUNE_THROW(Dune::InvalidStateException, "Try to triangulate point cloud with less than 3 points!");
 
    if (convexHullPoints.size() == 3)
        return std::vector<Triangle>(1, {convexHullPoints[0], convexHullPoints[1], convexHullPoints[2]});
 
    Point midPoint(0.0);
    for (const auto& p : convexHullPoints)
        midPoint += p;
    midPoint /= convexHullPoints.size();
 
    std::vector<Triangle> triangulation;
    triangulation.reserve(convexHullPoints.size());
 
    for (std::size_t i = 0; i < convexHullPoints.size()-1; ++i)
        triangulation.emplace_back(Triangle{midPoint, convexHullPoints[i], convexHullPoints[i+1]});
 
    triangulation.emplace_back(Triangle{midPoint, convexHullPoints[convexHullPoints.size()-1], convexHullPoints[0]});
 
    return triangulation;
}

template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>,
          class RandomAccessContainer,
          std::enable_if_t< std::is_same_v<Policy, TriangulationPolicy::ConvexHullPolicy>
                            && dim == 2, int> = 0 >
inline Triangulation<dim, dimWorld, typename RandomAccessContainer::value_type::value_type>
triangulate(const RandomAccessContainer& points)
{
    const auto convexHullPoints = grahamConvexHull<2>(points);
    return triangulate<dim, dimWorld, TriangulationPolicy::MidPointPolicy>(convexHullPoints);
}

template< int dim, int dimWorld, class Policy = TriangulationPolicy::DefaultPolicy<dim, dimWorld>,
          class RandomAccessContainer,
          std::enable_if_t< std::is_same_v<Policy, TriangulationPolicy::ConvexHullPolicy>
                            && dim == 3, 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>;
    using Tetrahedron = std::array<Point, 4>;
 
    static_assert(std::is_same_v<typename RandomAccessContainer::value_type, Point>,
                  "Triangulation expects Dune::FieldVector as point type");
 
    const auto numPoints = points.size();
    if (numPoints < 4)
        DUNE_THROW(Dune::InvalidStateException, "Trying to create 3D triangulation of point cloud with less than 4 points!");
 
    if (numPoints == 4)
        return std::vector<Tetrahedron>(1, {points[0], points[1], points[2], points[3]});
 
    // compute the mid point of the point cloud (not the midpoint of the convex hull but this
    // should not matter too much for the applications we have in mind here)
    Point midPoint(0.0);
    Point lowerLeft(1e100);
    Point upperRight(-1e100);
    for (const auto& p : points)
    {
        midPoint += p;
        for (int i = 0; i < dimWorld; ++i)
        {
            using std::max; using std::min;
            lowerLeft[i] = min(p[i], lowerLeft[i]);
            upperRight[i] = max(p[i], upperRight[i]);
        }
    }
    midPoint /= numPoints;
 
    auto magnitude = 0.0;
    using std::max;
    for (int i = 0; i < dimWorld; ++i)
        magnitude = max(upperRight[i] - lowerLeft[i], magnitude);
    const auto eps = 1e-7*magnitude;
    const auto eps2 = eps*eps;
    const auto epsDist = 1e-8*magnitude;
 
    // reserve memory conservatively to avoid reallocation
    std::vector<Tetrahedron> triangulation;
    triangulation.reserve(numPoints);
 
    // make a buffer for storing coplanar points and indices
    std::vector<Point> coplanarPointBuffer;
    coplanarPointBuffer.reserve(std::min<std::size_t>(12, numPoints-1));
 
    // remember coplanar cluster planes
    // coplanar clusters are uniquely identified by their plane (outward unit normal and offset)
    // we only want to add each cluster once (when handling the first triangle in the cluster)
    std::vector<std::pair<Point, ctype>> coplanarClusters;
    coplanarClusters.reserve(numPoints/3);
 
    // brute force algorithm: Try all possible triangles and check
    // if they are triangles of the convex hull. This is achieved by checking if all
    // other points are in the same half-space (side)
    // the algorithm is O(n^4), so only do this for very small point clouds
    for (int i = 0; i < numPoints; ++i)
    {
        for (int j = i+1; j < numPoints; ++j)
        {
            for (int k = j+1; k < numPoints; ++k)
            {
                const auto pointI = points[i];
                const auto ab = points[j] - pointI;
                const auto ac = points[k] - pointI;
                const auto normal = crossProduct(ab, ac);
 
                // clear list of coplanar points w.r.t to triangle ijk
                coplanarPointBuffer.clear();
 
                // check if this triangle ijk is admissible which means
                // it is on the convex hull (all other points in the cloud are in the same half-space/side)
                const bool isAdmissible = [&]()
                {
                    // check if points are colinear and we can't form a triangle
                    // if so, skip this triangle
                    if (normal.two_norm2() < eps2*eps2)
                        return false;
 
                    int marker = 0; // 0 means undecided side (or coplanar)
                    for (int m = 0; m < numPoints; ++m)
                    {
                        if (m != i && m != j && m != k)
                        {
                            // check scalar product with surface normal to decide side
                            const auto ad = points[m] - pointI;
                            const auto sp = normal*ad;
 
                            // if the sign changes wrt the previous sign, the triangle is not part of the convex hull
                            // sp is a triple product, so divide by the normal to obtain a distance
                            using std::abs; using std::signbit;
                            const bool coplanar = abs(sp) < epsDist*normal.two_norm();
                            int newMarker = coplanar ? 0 : signbit(sp) ? -1 : 1;
 
                            // make decision for a side as soon as the next marker is != 0
                            // keep track of the previous marker
                            if (marker == 0 && newMarker != 0)
                                marker = newMarker;
 
                            // if marker flips, not all points are on one side
                            // zero marker (undecided side / coplanar point) shouldn't abort the process
                            if (newMarker != 0 && marker != newMarker)
                                return false;
 
                            // handle possible coplanar points
                            if (coplanar)
                                coplanarPointBuffer.push_back(points[m]);
                        }
                    }
 
                    // if there are coplanar points complete the cluster with i, j, and k
                    // and store the cluster information for future lookup
                    if (!coplanarPointBuffer.empty())
                    {
                        auto unitNormal = normal;
                        unitNormal /= unitNormal.two_norm();
                        if (unitNormal*(pointI - midPoint) < 0.0)
                            unitNormal *= -1.0;
                        const auto offset = unitNormal*pointI;
 
                        using std::abs;
                        if (std::find_if(
                                coplanarClusters.begin(), coplanarClusters.end(),
                                [&](const auto& c){ return unitNormal*c.first > 1.0 - 1e-6 && abs(offset - c.second) < epsDist; }
                            ) != coplanarClusters.end())
                        {
                            // this cluster has already been handled
                            coplanarPointBuffer.clear();
                            return false;
                        }
 
                        coplanarPointBuffer.insert(coplanarPointBuffer.end(), { points[i], points[j], points[k] });
                        coplanarClusters.emplace_back(std::make_pair(unitNormal, offset));
                    }
 
                    // we require that not all points are coplanar, so
                    // there will be always at least one non-coplanar point
                    // to check on which side the other points are.
                    // Hence, once we get here, the triangle or coplanar point cluster is part of the convex hull
                    return true;
                }();
 
                if (isAdmissible)
                {
                    // check if we have a cluster of coplanar points forming on of the
                    // faces of the convex hull, if yes, compute (2d) convex hull first and triangulate
                    if (!coplanarPointBuffer.empty())
                    {
                        const auto triangles = triangulate<2, 3, TriangulationPolicy::ConvexHullPolicy>(coplanarPointBuffer);
                        for (const auto& triangle : triangles)
                        {
                            const auto ab = triangle[1] - triangle[0];
                            const auto ac = triangle[2] - triangle[0];
                            const auto normal = crossProduct(ab, ac);
                            const auto am = midPoint - triangle[0];
                            const auto sp = normal*am;
                            using std::signbit;
                            const bool isBelow = signbit(sp);
                            if (isBelow)
                                triangulation.emplace_back(Tetrahedron{
                                    triangle[0], triangle[2], triangle[1], midPoint
                                });
                            else
                                triangulation.emplace_back(Tetrahedron{
                                    triangle[0], triangle[1], triangle[2], midPoint
                                });
                        }
                    }
                    else
                    {
                        const auto am = midPoint - pointI;
                        const auto sp = normal*am;
                        using std::signbit;
                        const bool isBelow = signbit(sp);
                        if (isBelow)
                            triangulation.emplace_back(Tetrahedron{
                                pointI, points[k], points[j], midPoint
                            });
                        else
                            triangulation.emplace_back(Tetrahedron{
                                pointI, points[j], points[k], midPoint
                            });
                    }
                }
            }
        }
    }
 
    // sanity check: if points are not coplanar, then using the mid point policy, we get at least 4 tetrahedrons
    if (triangulation.size() < 4)
        DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the triangulation!");
 
    return triangulation;
}
 
} // end namespace Dumux
 
#endif