distance.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_DISTANCE_HH
#define DUMUX_GEOMETRY_DISTANCE_HH
 
#include <dune/common/fvector.hh>
#include <dune/geometry/quadraturerules.hh>
 
#include <dumux/common/math.hh>
#include <dumux/geometry/boundingboxtree.hh>
 
namespace Dumux {

template<class Geometry>
static inline typename Geometry::ctype
averageDistancePointGeometry(const typename Geometry::GlobalCoordinate& p,
                             const Geometry& geometry,
                             std::size_t integrationOrder = 2)
{
    typename Geometry::ctype avgDist = 0.0;
    const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
    for (const auto& qp : quad)
        avgDist += (geometry.global(qp.position())-p).two_norm()*qp.weight()*geometry.integrationElement(qp.position());
    return avgDist/geometry.volume();
}

template<class Point>
static inline typename Point::value_type
squaredDistancePointLine(const Point& p, const Point& a, const Point& b)
{
    const auto ab = b - a;
    const auto t = (p - a)*ab/ab.two_norm2();
    auto proj = a;
    proj.axpy(t, ab);
    return (proj - p).two_norm2();
}

template<class Geometry>
static inline typename Geometry::ctype
squaredDistancePointLine(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
{
    static_assert(Geometry::mydimension == 1, "Geometry has to be a line");
    const auto& a = geometry.corner(0);
    const auto& b = geometry.corner(1);
    return squaredDistancePointLine(p, a, b);
}

template<class Point>
static inline typename Point::value_type
distancePointLine(const Point& p, const Point& a, const Point& b)
{ using std::sqrt; return sqrt(squaredDistancePointLine(p, a, b)); }

template<class Geometry>
static inline typename Geometry::ctype
distancePointLine(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
{ using std::sqrt; return sqrt(squaredDistancePointLine(p, geometry)); }

template<class Point>
static inline typename Point::value_type
squaredDistancePointSegment(const Point& p, const Point& a, const Point& b)
{
    const auto ab = b - a;
    const auto ap = p - a;
    const auto t = ap*ab;
 
    if (t <= 0.0)
        return ap.two_norm2();
 
    const auto lengthSq = ab.two_norm2();
    if (t >= lengthSq)
        return (b - p).two_norm2();
 
    auto proj = a;
    proj.axpy(t/lengthSq, ab);
    return (proj - p).two_norm2();
}

template<class Geometry>
static inline typename Geometry::ctype
squaredDistancePointSegment(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
{
    static_assert(Geometry::mydimension == 1, "Geometry has to be a segment");
    const auto& a = geometry.corner(0);
    const auto& b = geometry.corner(1);
    return squaredDistancePointSegment(p, a, b);
}

template<class Point>
static inline typename Point::value_type
distancePointSegment(const Point& p, const Point& a, const Point& b)
{ using std::sqrt; return sqrt(squaredDistancePointSegment(p, a, b)); }

template<class Geometry>
static inline typename Geometry::ctype
distancePointSegment(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
{ using std::sqrt; return sqrt(squaredDistancePointSegment(p, geometry)); }

template<class Point>
static inline typename Point::value_type
squaredDistancePointTriangle(const Point& p, const Point& a, const Point& b, const Point& c)
{
    static_assert(Point::dimension == 3, "Only works in 3D");
    const auto ab = b - a;
    const auto bc = c - b;
    const auto ca = a - c;
    const auto normal = crossProduct(ab, ca);
 
    const auto ap = p - a;
    const auto bp = p - b;
    const auto cp = p - c;
 
    const auto sum = sign(crossProduct(ab, normal)*ap)
                   + sign(crossProduct(bc, normal)*bp)
                   + sign(crossProduct(ca, normal)*cp);
 
    // if there is no orthogonal projection
    // (point is outside the infinite prism implied by the triangle and its surface normal)
    // compute distance to the edges (codim-1 facets)
    if (sum < 2.0)
    {
        using std::min;
        return min({squaredDistancePointSegment(p, a, b),
                    squaredDistancePointSegment(p, a, c),
                    squaredDistancePointSegment(p, b, c)});
    }
    // compute distance via orthogonal projection
    else
    {
        const auto tmp = normal*ap;
        return tmp*tmp / (normal*normal);
    }
}

template<class Geometry>
static inline typename Geometry::ctype
squaredDistancePointTriangle(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
{
    static_assert(Geometry::coorddimension == 3, "Only works in 3D");
    static_assert(Geometry::mydimension == 2, "Geometry has to be a triangle");
    assert(geometry.corners() == 3);
    const auto& a = geometry.corner(0);
    const auto& b = geometry.corner(1);
    const auto& c = geometry.corner(2);
    return squaredDistancePointTriangle(p, a, b, c);
}

template<class Point>
static inline typename Point::value_type
distancePointTriangle(const Point& p, const Point& a, const Point& b, const Point& c)
{ using std::sqrt; return sqrt(squaredDistancePointTriangle(p, a, b, c)); }

template<class Geometry>
static inline typename Geometry::ctype
distancePointTriangle(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
{ using std::sqrt; return sqrt(squaredDistancePointTriangle(p, geometry)); }

template<class Geometry>
static inline typename Geometry::ctype
squaredDistancePointPolygon(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
{
    if (geometry.corners() == 3)
        return squaredDistancePointTriangle(p, geometry);
    else if (geometry.corners() == 4)
    {
        const auto& a = geometry.corner(0);
        const auto& b = geometry.corner(1);
        const auto& c = geometry.corner(2);
        const auto& d = geometry.corner(3);
 
        using std::min;
        return min(squaredDistancePointTriangle(p, a, b, d),
                   squaredDistancePointTriangle(p, a, d, c));
    }
    else
        DUNE_THROW(Dune::NotImplemented, "Polygon with " << geometry.corners() << " corners not supported");
}

template<class Geometry>
static inline typename Geometry::ctype
distancePointPolygon(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
{ using std::sqrt; return sqrt(squaredDistancePointPolygon(p, geometry)); }

template<class Geometry>
static inline typename Geometry::ctype
averageDistanceSegmentGeometry(const typename Geometry::GlobalCoordinate& a,
                               const typename Geometry::GlobalCoordinate& b,
                               const Geometry& geometry,
                               std::size_t integrationOrder = 2)
{
    typename Geometry::ctype avgDist = 0.0;
    const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
    for (const auto& qp : quad)
        avgDist += distancePointSegment(geometry.global(qp.position()), a, b)*qp.weight()*geometry.integrationElement(qp.position());
    return avgDist/geometry.volume();
}

template<class Geometry>
static inline typename Geometry::ctype
averageDistanceLineGeometry(const typename Geometry::GlobalCoordinate& a,
                            const typename Geometry::GlobalCoordinate& b,
                            const Geometry& geometry,
                            std::size_t integrationOrder = 2)
{
    typename Geometry::ctype avgDist = 0.0;
    const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
    for (const auto& qp : quad)
        avgDist += distancePointLine(geometry.global(qp.position()), a, b)*qp.weight()*geometry.integrationElement(qp.position());
    return avgDist/geometry.volume();
}

template<class ctype, int dimWorld>
static inline ctype distance(const Dune::FieldVector<ctype, dimWorld>& a,
                             const Dune::FieldVector<ctype, dimWorld>& b)
{ return (a-b).two_norm(); }

template<class ctype, int dimWorld>
static inline ctype squaredDistance(const Dune::FieldVector<ctype, dimWorld>& a,
                                    const Dune::FieldVector<ctype, dimWorld>& b)
{ return (a-b).two_norm2(); }
 
 
#ifndef DOXYGEN
namespace Detail {
 
// helper struct to compute distance between two geometries, specialized below
template<class Geo1, class Geo2,
         int dimWorld = Geo1::coorddimension,
         int dim1 = Geo1::mydimension, int dim2 = Geo2::mydimension>
struct GeometrySquaredDistance
{
    static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
    static auto distance(const Geo1& geo1, const Geo2& geo2)
    {
        DUNE_THROW(Dune::NotImplemented, "Geometry distance computation not implemented for dimworld = "
                    << dimWorld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
    }
};
 
// distance point-point
template<class Geo1, class Geo2, int dimWorld>
struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 0, 0>
{
    static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
    static auto distance(const Geo1& geo1, const Geo2& geo2)
    { return Dumux::squaredDistance(geo1.corner(0), geo2.corner(0)); }
};
 
// distance segment-point
template<class Geo1, class Geo2, int dimWorld>
struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 1, 0>
{
    static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
    static auto distance(const Geo1& geo1, const Geo2& geo2)
    { return squaredDistancePointSegment(geo2.corner(0), geo1); }
};
 
// distance point-segment
template<class Geo1, class Geo2, int dimWorld>
struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 0, 1>
{
    static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
    static auto distance(const Geo1& geo1, const Geo2& geo2)
    { return squaredDistancePointSegment(geo1.corner(0), geo2); }
};
 
// distance point-polygon
template<class Geo1, class Geo2, int dimWorld>
struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 0, 2>
{
    static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
    static inline auto distance(const Geo1& geo1, const Geo2& geo2)
    { return squaredDistancePointPolygon(geo1.corner(0), geo2); }
};
 
// distance polygon-point
template<class Geo1, class Geo2, int dimWorld>
struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 2, 0>
{
    static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
    static inline auto distance(const Geo1& geo1, const Geo2& geo2)
    { return squaredDistancePointPolygon(geo2.corner(0), geo1); }
};
 
} // end namespace Detail
#endif // DOXYGEN

template<class Geo1, class Geo2>
static inline auto squaredDistance(const Geo1& geo1, const Geo2& geo2)
{ return Detail::GeometrySquaredDistance<Geo1, Geo2>::distance(geo1, geo2); }

template<class Geo1, class Geo2>
static inline auto distance(const Geo1& geo1, const Geo2& geo2)
{ using std::sqrt; return sqrt(squaredDistance(geo1, geo2)); }
 
 
#ifndef DOXYGEN
namespace Detail {

template<class EntitySet, class ctype, int dimworld,
         typename std::enable_if_t<(EntitySet::Entity::Geometry::mydimension > 0), int> = 0>
void closestEntity(const Dune::FieldVector<ctype, dimworld>& point,
                   const BoundingBoxTree<EntitySet>& tree,
                   std::size_t node,
                   ctype& minSquaredDistance,
                   std::size_t& eIdx)
{
    // Get the bounding box for the current node
    const auto& bBox = tree.getBoundingBoxNode(node);
 
    // If bounding box is outside radius, then don't search further
    const auto squaredDistance = squaredDistancePointBoundingBox(
        point, tree.getBoundingBoxCoordinates(node)
    );
 
    // do not continue if the AABB is further away than the current minimum
    if (squaredDistance > minSquaredDistance) return;
 
    // If the bounding box is a leaf, update distance and index with primitive test
    else if (tree.isLeaf(bBox, node))
    {
        const std::size_t entityIdx = bBox.child1;
        const auto squaredDistance = [&]{
            const auto geometry = tree.entitySet().entity(entityIdx).geometry();
            if constexpr (EntitySet::Entity::Geometry::mydimension == 2)
                return squaredDistancePointPolygon(point, geometry);
            else if constexpr (EntitySet::Entity::Geometry::mydimension == 1)
                return squaredDistancePointSegment(point, geometry);
            else
                DUNE_THROW(Dune::NotImplemented, "squaredDistance to entity with dim>2");
        }();
 
        if (squaredDistance < minSquaredDistance)
        {
            eIdx = entityIdx;
            minSquaredDistance = squaredDistance;
        }
    }
 
    // Check the children nodes recursively
    else
    {
        closestEntity(point, tree, bBox.child0, minSquaredDistance, eIdx);
        closestEntity(point, tree, bBox.child1, minSquaredDistance, eIdx);
    }
}

template<class EntitySet, class ctype, int dimworld,
         typename std::enable_if_t<(EntitySet::Entity::Geometry::mydimension == 0), int> = 0>
void closestEntity(const Dune::FieldVector<ctype, dimworld>& point,
                   const BoundingBoxTree<EntitySet>& tree,
                   std::size_t node,
                   ctype& minSquaredDistance,
                   std::size_t& eIdx)
{
    // Get the bounding box for the current node
    const auto& bBox = tree.getBoundingBoxNode(node);
 
    // If the bounding box is a leaf, update distance and index with primitive test
    if (tree.isLeaf(bBox, node))
    {
        const std::size_t entityIdx = bBox.child1;
        const auto& p = tree.entitySet().entity(entityIdx).geometry().corner(0);
        const auto squaredDistance = (p-point).two_norm2();
 
        if (squaredDistance < minSquaredDistance)
        {
            eIdx = entityIdx;
            minSquaredDistance = squaredDistance;
        }
    }
 
    // Check the children nodes recursively
    else
    {
        // If bounding box is outside radius, then don't search further
        const auto squaredDistance = squaredDistancePointBoundingBox(
            point, tree.getBoundingBoxCoordinates(node)
        );
 
        // do not continue if the AABB is further away than the current minimum
        if (squaredDistance > minSquaredDistance) return;
 
        closestEntity(point, tree, bBox.child0, minSquaredDistance, eIdx);
        closestEntity(point, tree, bBox.child1, minSquaredDistance, eIdx);
    }
}
 
} // end namespace Detail
#endif // DOXYGEN

template<class EntitySet, class ctype, int dimworld>
std::pair<ctype, std::size_t> closestEntity(const Dune::FieldVector<ctype, dimworld>& point,
                                            const BoundingBoxTree<EntitySet>& tree,
                                            ctype minSquaredDistance = std::numeric_limits<ctype>::max())
{
    std::size_t eIdx = 0;
    Detail::closestEntity(point, tree, tree.numBoundingBoxes() - 1, minSquaredDistance, eIdx);
    using std::sqrt;
    return { minSquaredDistance, eIdx };
}

template<class EntitySet, class ctype, int dimworld>
ctype squaredDistance(const Dune::FieldVector<ctype, dimworld>& point,
                      const BoundingBoxTree<EntitySet>& tree,
                      ctype minSquaredDistance = std::numeric_limits<ctype>::max())
{
    return closestEntity(point, tree, minSquaredDistance).first;
}

template<class EntitySet, class ctype, int dimworld>
ctype distance(const Dune::FieldVector<ctype, dimworld>& point,
               const BoundingBoxTree<EntitySet>& tree,
               ctype minSquaredDistance = std::numeric_limits<ctype>::max())
{
    using std::sqrt;
    return sqrt(squaredDistance(point, tree, minSquaredDistance));
}
 
} // end namespace Dumux
 
#endif