makegeometry.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-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_GEOMETRY_MAKE_GEOMETRY_HH
#define DUMUX_GEOMETRY_MAKE_GEOMETRY_HH
 
#include <vector>
#include <array>
#include <limits>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/exceptions.hh>
#include <dune/geometry/multilineargeometry.hh>
#include <dumux/common/math.hh>
#include <dumux/geometry/grahamconvexhull.hh>
 
namespace Dumux {
 
template<class CoordScalar>
bool pointsAreCoplanar(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points, const CoordScalar scale)
{
    if (points.size() != 4)
        DUNE_THROW(Dune::InvalidStateException, "Check only works for 4 points!");
 
    // (see "Real-Time Collision Detection" by Christer Ericson)
    Dune::FieldMatrix<CoordScalar, 4, 4> M;
    for (int i = 0; i < 3; ++i )
        M[i] = {points[0][i], points[1][i], points[2][i], points[3][i]};
    M[3] = {1.0*scale, 1.0*scale, 1.0*scale, 1.0*scale};
 
    using std::abs;
    return abs(M.determinant()) < 1.5e-7*scale*scale*scale*scale;
}
 
template<class CoordScalar>
bool pointsAreCoplanar(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
{
    Dune::FieldVector<CoordScalar, 3> bBoxMin(std::numeric_limits<CoordScalar>::max());
    Dune::FieldVector<CoordScalar, 3> bBoxMax(std::numeric_limits<CoordScalar>::lowest());
    for (const auto& p : points)
    {
        for (int i=0; i<3; i++)
        {
            using std::min;
            using std::max;
            bBoxMin[i] = min(bBoxMin[i], p[i]);
            bBoxMax[i] = max(bBoxMax[i], p[i]);
        }
    }
 
    const auto size = (bBoxMax - bBoxMin).two_norm();
 
    return pointsAreCoplanar(points, size);
}
 
template<class CoordScalar>
std::vector<Dune::FieldVector<CoordScalar, 3>> getReorderedPoints(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
{
    std::array<int, 4> tmp;
    return getReorderedPoints(points, tmp);
}
 
template<class CoordScalar>
std::vector<Dune::FieldVector<CoordScalar, 3>> getReorderedPoints(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points,
                                                             std::array<int, 4>& orientations)
{
    if(points.size() == 4)
    {
        auto& p0 = points[0];
        auto& p1 = points[1];
        auto& p2 = points[2];
        auto& p3 = points[3];
 
        // check if the points define a proper quadrilateral
        const auto normal = crossProduct((p1 - p0), (p2 - p0));
 
        orientations = { getOrientation(p0, p3, p2, normal),
                         getOrientation(p0, p3, p1, normal),
                         getOrientation(p2, p1, p0, normal),
                         getOrientation(p2, p1, p3, normal) };
 
 
        // check if the points follow the dune ordering (see http://www.dcs.gla.ac.uk/~pat/52233/slides/Geometry1x1.pdf)
        const bool diagonalsIntersect = (orientations[0] != orientations[1]) && (orientations[2] != orientations[3]);
 
        // the points conform with the dune ordering
        if(diagonalsIntersect)
            return points;
 
        // the points do not conform with the dune ordering, re-order
        using GlobalPosition = Dune::FieldVector<CoordScalar, 3>;
        if(!diagonalsIntersect && orientations[0] == 1)
            return std::vector<GlobalPosition>{p1, p0, p2, p3};
        else if(!diagonalsIntersect && orientations[0] == -1)
            return std::vector<GlobalPosition>{p3, p1, p0, p2};
        else
            DUNE_THROW(Dune::InvalidStateException, "Could not reorder points");
    }
    else
        DUNE_THROW(Dune::NotImplemented, "Reorder for " << points.size() << " points.");
}
 
template<class CoordScalar, bool enableSanityCheck = true>
auto makeDuneQuadrilaterial(const std::vector<Dune::FieldVector<CoordScalar, 3>>& points)
{
    if (points.size() != 4)
        DUNE_THROW(Dune::InvalidStateException, "A quadrilateral needs 4 corner points!");
 
    using GlobalPosition = Dune::FieldVector<CoordScalar, 3>;
    static constexpr auto coordDim = GlobalPosition::dimension;
    static constexpr auto dim = coordDim-1;
    using GeometryType = Dune::MultiLinearGeometry<CoordScalar, dim, coordDim>;
 
    // if no sanity check if desired, use the given points directly to construct the geometry
    if (!enableSanityCheck)
        return GeometryType(Dune::GeometryTypes::quadrilateral, points);
 
    // compute size
    Dune::FieldVector<CoordScalar, 3> bBoxMin(std::numeric_limits<CoordScalar>::max());
    Dune::FieldVector<CoordScalar, 3> bBoxMax(std::numeric_limits<CoordScalar>::lowest());
    for (const auto& p : points)
    {
        for (int i = 0; i < 3; i++)
        {
            using std::min;
            using std::max;
            bBoxMin[i] = min(bBoxMin[i], p[i]);
            bBoxMax[i] = max(bBoxMax[i], p[i]);
        }
    }
 
    const auto size = (bBoxMax - bBoxMin).two_norm();
 
    // otherwise, perform a number of checks and corrections
    if (!pointsAreCoplanar(points, size))
        DUNE_THROW(Dune::InvalidStateException, "Points do not lie within a plane");
 
    auto corners = grahamConvexHull<2>(points);
    if (corners.size() != 4)
        DUNE_THROW(Dune::InvalidStateException, "Points do not span a strictly convex polygon!");
 
    // make sure points conform with dune ordering
    std::array<int, 4> orientations;
    corners = getReorderedPoints(corners, orientations);
 
    if (std::any_of(orientations.begin(), orientations.end(), [](auto i){ return i == 0; }))
        DUNE_THROW(Dune::InvalidStateException, "More than two points lie on the same line.");
 
    const auto quadrilateral = GeometryType(Dune::GeometryTypes::quadrilateral, corners);
 
    const auto eps = 1e-7;
    if (quadrilateral.volume() < eps*size*size)
        DUNE_THROW(Dune::InvalidStateException, "Something went wrong, geometry has area of zero");
 
    return quadrilateral;
}
 
 
 
} // end namespace Dumux
 
#endif