intersectingentities.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_INTERSECTING_ENTITIES_HH
#define DUMUX_GEOMETRY_INTERSECTING_ENTITIES_HH
 
#include <cmath>
#include <type_traits>
#include <vector>
#include <algorithm>
#include <limits>
 
#include <dune/common/fvector.hh>
 
#include <dumux/common/math.hh>
#include <dumux/geometry/boundingboxtree.hh>
#include <dumux/geometry/geometrytraits.hh>
#include <dumux/geometry/intersectspointgeometry.hh>
#include <dumux/geometry/geometryintersection.hh>
#include <dumux/geometry/triangulation.hh>
 
namespace Dumux {

template<int dimworld, class CoordTypeA, class CoordTypeB = CoordTypeA>
class IntersectionInfo
{
public:
    using ctype = typename Dune::PromotionTraits<CoordTypeA, CoordTypeB>::PromotedType;
    static constexpr int dimensionworld = dimworld;
    using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
 
    template<class Corners>
    explicit IntersectionInfo(std::size_t a, std::size_t b, Corners&& c)
    : a_(a)
    , b_(b)
    , corners_(c.begin(), c.end())
    {}

    std::size_t first() const
    { return a_; }

    std::size_t second() const
    { return b_; }

    const std::vector<GlobalPosition>& corners() const
    { return corners_; }

    bool cornersMatch(const std::vector<GlobalPosition>& otherCorners) const
    {
        if (otherCorners.size() != corners_.size())
            return false;
 
        using std::max;
        ctype eps2 = std::numeric_limits<ctype>::min();
        for (int i = 1; i < corners_.size(); ++i)
            eps2 = max(eps2, (corners_[i] - corners_[0]).two_norm2());
 
        // We use a base epsilon of 1.5e-7 for comparisons of lengths.
        // Since here we compare squared lengths, we multiply by its square.
        eps2 *= 1.5e-7*1.5e-7;
 
        for (int i = 0; i < corners_.size(); ++i)
            // early return if none of the other corners are equal to this corner
            if (std::none_of(otherCorners.begin(),
                             otherCorners.end(),
                             [&] (const auto& other) { return (corners_[i] - other).two_norm2() < eps2; }))
                return false;
 
        return true;
    }
 
private:
    std::size_t a_, b_; 
    std::vector<GlobalPosition> corners_; 
};

template<class EntitySet, class ctype, int dimworld>
inline std::vector<std::size_t>
intersectingEntities(const Dune::FieldVector<ctype, dimworld>& point,
                     const BoundingBoxTree<EntitySet>& tree,
                     bool isCartesianGrid = false)
{
    // Call the recursive find function to find candidates
    std::vector<std::size_t> entities;
    intersectingEntities(point, tree, tree.numBoundingBoxes() - 1, entities, isCartesianGrid);
    return entities;
}

template<class EntitySet, class ctype, int dimworld>
void intersectingEntities(const Dune::FieldVector<ctype, dimworld>& point,
                          const BoundingBoxTree<EntitySet>& tree,
                          std::size_t node,
                          std::vector<std::size_t>& entities,
                          bool isCartesianGrid = false)
{
    // Get the bounding box for the current node
    const auto& bBox = tree.getBoundingBoxNode(node);
 
    // if the point is not in the bounding box we can stop
    if (!intersectsPointBoundingBox(point, tree.getBoundingBoxCoordinates(node))) return;
 
    // now we know it's inside
    // if the box is a leaf do the primitive test.
    else if (tree.isLeaf(bBox, node))
    {
        const std::size_t entityIdx = bBox.child1;
        using Geometry = typename EntitySet::Entity::Geometry;
        // for structured cube grids or axis-aligned geometries the bounding box
        // test is already exact, so we can skip the primitive test
        if (isCartesianGrid || GeometryTraits<Geometry>::axisAligned)
            entities.push_back(entityIdx);
        else
        {
            const auto geometry = tree.entitySet().entity(entityIdx).geometry();
            // if the primitive is positive it intersects the actual geometry, add the entity to the list
            if (intersectsPointGeometry(point, geometry))
                entities.push_back(entityIdx);
        }
    }
 
    // No leaf. Check both children nodes.
    else
    {
        intersectingEntities(point, tree, bBox.child0, entities, isCartesianGrid);
        intersectingEntities(point, tree, bBox.child1, entities, isCartesianGrid);
    }
}

template<class EntitySet, class ctype>
void intersectingBoxes(const ctype* box,
                       const BoundingBoxTree<EntitySet>& tree,
                       std::size_t node,
                       std::vector<std::size_t>& entities)
{
    // if the query box doesn't overlap the node's box we can stop
    static constexpr int dimworld = EntitySet::dimensionworld;
    if (!intersectsBoundingBoxBoundingBox<dimworld>(box, tree.getBoundingBoxCoordinates(node)))
        return;
 
    // if the box is a leaf the bounding box overlap is the result
    const auto& bBox = tree.getBoundingBoxNode(node);
    if (tree.isLeaf(bBox, node))
        entities.push_back(bBox.child1);
 
    // no leaf, check both children nodes
    else
    {
        intersectingBoxes(box, tree, bBox.child0, entities);
        intersectingBoxes(box, tree, bBox.child1, entities);
    }
}

template<class EntitySet, class ctype>
inline std::vector<std::size_t>
intersectingBoxes(const ctype* box, const BoundingBoxTree<EntitySet>& tree)
{
    std::vector<std::size_t> entities;
    intersectingBoxes(box, tree, tree.numBoundingBoxes() - 1, entities);
    return entities;
}

template<class EntitySet0, class EntitySet1>
void intersectingBoxes(const BoundingBoxTree<EntitySet0>& treeA,
                       const BoundingBoxTree<EntitySet1>& treeB,
                       std::size_t nodeA, std::size_t nodeB,
                       std::vector<std::pair<std::size_t, std::size_t>>& pairs)
{
    // if the two bounding boxes of the current nodes don't overlap we can stop searching
    static constexpr int dimworld = EntitySet0::dimensionworld;
    if (!intersectsBoundingBoxBoundingBox<dimworld>(treeA.getBoundingBoxCoordinates(nodeA),
                                                    treeB.getBoundingBoxCoordinates(nodeB)))
        return;
 
    const auto& bBoxA = treeA.getBoundingBoxNode(nodeA);
    const auto& bBoxB = treeB.getBoundingBoxNode(nodeB);
    const bool isLeafA = treeA.isLeaf(bBoxA, nodeA);
    const bool isLeafB = treeB.isLeaf(bBoxB, nodeB);
 
    // both leaves: their boxes overlap, record the entity index pair
    if (isLeafA && isLeafB)
        pairs.emplace_back(bBoxA.child1, bBoxB.child1);
 
    // descend into the tree that hasn't reached a leaf yet, larger tree (bigger node number) first
    else if (isLeafA)
    {
        intersectingBoxes(treeA, treeB, nodeA, bBoxB.child0, pairs);
        intersectingBoxes(treeA, treeB, nodeA, bBoxB.child1, pairs);
    }
    else if (isLeafB)
    {
        intersectingBoxes(treeA, treeB, bBoxA.child0, nodeB, pairs);
        intersectingBoxes(treeA, treeB, bBoxA.child1, nodeB, pairs);
    }
    else if (nodeA > nodeB)
    {
        intersectingBoxes(treeA, treeB, bBoxA.child0, nodeB, pairs);
        intersectingBoxes(treeA, treeB, bBoxA.child1, nodeB, pairs);
    }
    else
    {
        intersectingBoxes(treeA, treeB, nodeA, bBoxB.child0, pairs);
        intersectingBoxes(treeA, treeB, nodeA, bBoxB.child1, pairs);
    }
}

template<class EntitySet0, class EntitySet1>
inline std::vector<std::pair<std::size_t, std::size_t>>
intersectingBoxes(const BoundingBoxTree<EntitySet0>& treeA, const BoundingBoxTree<EntitySet1>& treeB)
{
    static_assert(int(EntitySet0::dimensionworld) == int(EntitySet1::dimensionworld),
        "Can only intersect bounding box trees of same world dimension");
 
    std::vector<std::pair<std::size_t, std::size_t>> pairs;
    intersectingBoxes(treeA, treeB, treeA.numBoundingBoxes() - 1, treeB.numBoundingBoxes() - 1, pairs);
    return pairs;
}

template<class Geometry, class EntitySet>
inline std::vector<IntersectionInfo<Geometry::coorddimension, typename Geometry::ctype, typename EntitySet::ctype>>
intersectingEntities(const Geometry& geometry,
                     const BoundingBoxTree<EntitySet>& tree)
{
    using IP = typename IntersectionPolicy::DefaultPolicy<Geometry, typename EntitySet::Entity::Geometry>;
    return intersectingEntities(geometry, tree, IP{});
}

template<class Geometry, class EntitySet, class IntersectionPolicy>
inline std::vector<IntersectionInfo<Geometry::coorddimension, typename Geometry::ctype, typename EntitySet::ctype>>
intersectingEntities(const Geometry& geometry,
                     const BoundingBoxTree<EntitySet>& tree,
                     IntersectionPolicy intersectionPolicy)
{
    // check if the world dimensions match
    static_assert(int(Geometry::coorddimension) == int(EntitySet::dimensionworld),
        "Can only intersect geometry and bounding box tree of same world dimension");
 
    // Create data structure for return type
    std::vector<IntersectionInfo<Geometry::coorddimension, typename Geometry::ctype, typename EntitySet::ctype>> intersections;
    using ctype = typename IntersectionInfo<Geometry::coorddimension, typename Geometry::ctype, typename EntitySet::ctype>::ctype;
    static constexpr int dimworld = Geometry::coorddimension;
 
    // compute the bounding box of the given geometry
    std::array<ctype, 2*Geometry::coorddimension> bBox;
    ctype* xMin = bBox.data(); ctype* xMax = xMin + Geometry::coorddimension;
 
    // Get coordinates of first vertex
    auto corner = geometry.corner(0);
    for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
        xMin[dimIdx] = xMax[dimIdx] = corner[dimIdx];
 
    // Compute the min and max over the remaining vertices
    for (std::size_t cornerIdx = 1; cornerIdx < geometry.corners(); ++cornerIdx)
    {
        corner = geometry.corner(cornerIdx);
        for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
        {
            using std::max;
            using std::min;
            xMin[dimIdx] = min(xMin[dimIdx], corner[dimIdx]);
            xMax[dimIdx] = max(xMax[dimIdx], corner[dimIdx]);
        }
    }
 
    // Call the recursive find function to find candidates
    intersectingEntities(geometry, tree,
                         bBox, tree.numBoundingBoxes() - 1,
                         intersections, intersectionPolicy);
 
    return intersections;
}

template<class Geometry, class EntitySet>
void intersectingEntities(const Geometry& geometry,
                          const BoundingBoxTree<EntitySet>& tree,
                          const std::array<typename Geometry::ctype, 2*Geometry::coorddimension>& bBox,
                          std::size_t nodeIdx,
                          std::vector<IntersectionInfo<Geometry::coorddimension,
                                                       typename Geometry::ctype,
                                                       typename EntitySet::ctype>>& intersections)
{
    using IP = typename IntersectionPolicy::DefaultPolicy<Geometry, typename EntitySet::Entity::Geometry>;
    intersectingEntities(geometry, tree, bBox, nodeIdx, intersections, IP{});
}

template<class Geometry, class EntitySet, class IntersectionPolicy>
void intersectingEntities(const Geometry& geometry,
                          const BoundingBoxTree<EntitySet>& tree,
                          const std::array<typename Geometry::ctype, 2*Geometry::coorddimension>& bBox,
                          std::size_t nodeIdx,
                          std::vector<IntersectionInfo<Geometry::coorddimension,
                                                       typename Geometry::ctype,
                                                       typename EntitySet::ctype>>& intersections,
                          IntersectionPolicy intersectionPolicy)
{
    // if the two bounding boxes don't intersect we can stop searching
    static constexpr int dimworld = Geometry::coorddimension;
    if (!intersectsBoundingBoxBoundingBox<dimworld>(bBox.data(), tree.getBoundingBoxCoordinates(nodeIdx)))
        return;
 
    // get node info for current bounding box node
    const auto& bBoxNode = tree.getBoundingBoxNode(nodeIdx);
 
    // if the box is a leaf do the primitive test.
    if (tree.isLeaf(bBoxNode, nodeIdx))
    {
        // eIdxA is always 0 since we intersect with exactly one geometry
        const auto eIdxA = 0;
        const auto eIdxB = bBoxNode.child1;
 
        const auto geometryTree = tree.entitySet().entity(eIdxB).geometry();
        using GeometryTree = std::decay_t<decltype(geometryTree)>;
        using IntersectionAlgorithm = GeometryIntersection<Geometry, GeometryTree, IntersectionPolicy>;
        using Intersection = typename IntersectionAlgorithm::Intersection;
        Intersection intersection;
 
        if (IntersectionAlgorithm::intersection(geometry, geometryTree, intersection))
        {
            static constexpr int dimIntersection = IntersectionPolicy::dimIntersection;
            if constexpr (dimIntersection >= 2)
            {
                const auto triangulation = triangulate<dimIntersection, dimworld>(intersection);
                for (unsigned int i = 0; i < triangulation.size(); ++i)
                    intersections.emplace_back(eIdxA, eIdxB, std::move(triangulation[i]));
            }
            else
                intersections.emplace_back(eIdxA, eIdxB, intersection);
        }
    }
 
    // No leaf. Check both children nodes.
    else
    {
        intersectingEntities(geometry, tree, bBox, bBoxNode.child0, intersections);
        intersectingEntities(geometry, tree, bBox, bBoxNode.child1, intersections);
    }
}

template<class EntitySet0, class EntitySet1>
inline std::vector<IntersectionInfo<EntitySet0::dimensionworld, typename EntitySet0::ctype, typename EntitySet1::ctype>>
intersectingEntities(const BoundingBoxTree<EntitySet0>& treeA,
                     const BoundingBoxTree<EntitySet1>& treeB)
{
    using IP = typename IntersectionPolicy::DefaultPolicy<typename EntitySet0::Entity::Geometry, typename EntitySet1::Entity::Geometry>;
    return intersectingEntities(treeA, treeB, IP{});
}

template<class EntitySet0, class EntitySet1, class IntersectionPolicy>
inline std::vector<IntersectionInfo<EntitySet0::dimensionworld, typename EntitySet0::ctype, typename EntitySet1::ctype>>
intersectingEntities(const BoundingBoxTree<EntitySet0>& treeA,
                     const BoundingBoxTree<EntitySet1>& treeB,
                     IntersectionPolicy intersectionPolicy)
{
    // check if the world dimensions match
    static_assert(int(EntitySet0::dimensionworld) == int(EntitySet1::dimensionworld),
        "Can only intersect bounding box trees of same world dimension");
 
    // Create data structure for return type
    std::vector<IntersectionInfo<EntitySet0::dimensionworld, typename EntitySet0::ctype, typename EntitySet1::ctype>> intersections;
 
    // Call the recursive find function to find candidates
    intersectingEntities(treeA, treeB,
                         treeA.numBoundingBoxes() - 1,
                         treeB.numBoundingBoxes() - 1,
                         intersections, intersectionPolicy);
 
    return intersections;
}

template<class EntitySet0, class EntitySet1>
void intersectingEntities(const BoundingBoxTree<EntitySet0>& treeA,
                          const BoundingBoxTree<EntitySet1>& treeB,
                          std::size_t nodeA, std::size_t nodeB,
                          std::vector<IntersectionInfo<EntitySet0::dimensionworld,
                                                       typename EntitySet0::ctype,
                                                       typename EntitySet1::ctype>>& intersections)
{
    using IP = typename IntersectionPolicy::DefaultPolicy<typename EntitySet0::Entity::Geometry, typename EntitySet1::Entity::Geometry>;
    intersectingEntities(treeA, treeB, nodeA, nodeB, intersections, IP{});
}

template<class EntitySet0, class EntitySet1, class IntersectionPolicy>
void intersectingEntities(const BoundingBoxTree<EntitySet0>& treeA,
                          const BoundingBoxTree<EntitySet1>& treeB,
                          std::size_t nodeA, std::size_t nodeB,
                          std::vector<IntersectionInfo<EntitySet0::dimensionworld,
                                                       typename EntitySet0::ctype,
                                                       typename EntitySet1::ctype>>& intersections,
                          IntersectionPolicy intersectionPolicy)
{
    // Get the bounding box for the current node
    const auto& bBoxA = treeA.getBoundingBoxNode(nodeA);
    const auto& bBoxB = treeB.getBoundingBoxNode(nodeB);
 
    // if the two bounding boxes of the current nodes don't intersect we can stop searching
    static constexpr int dimworld = EntitySet0::dimensionworld;
    if (!intersectsBoundingBoxBoundingBox<dimworld>(treeA.getBoundingBoxCoordinates(nodeA),
                                                    treeB.getBoundingBoxCoordinates(nodeB)))
        return;
 
    // Check if we have a leaf in treeA or treeB
    const bool isLeafA = treeA.isLeaf(bBoxA, nodeA);
    const bool isLeafB = treeB.isLeaf(bBoxB, nodeB);
 
    // If both boxes are leaves do the primitive test
    if (isLeafA && isLeafB)
    {
        const auto eIdxA = bBoxA.child1;
        const auto eIdxB = bBoxB.child1;
 
        const auto geometryA = treeA.entitySet().entity(eIdxA).geometry();
        const auto geometryB = treeB.entitySet().entity(eIdxB).geometry();
 
        using GeometryA = std::decay_t<decltype(geometryA)>;
        using GeometryB = std::decay_t<decltype(geometryB)>;
        using IntersectionAlgorithm = GeometryIntersection<GeometryA, GeometryB, IntersectionPolicy>;
        using Intersection = typename IntersectionAlgorithm::Intersection;
 
        if (Intersection intersection; IntersectionAlgorithm::intersection(geometryA, geometryB, intersection))
        {
            static constexpr int dimIntersection = IntersectionPolicy::dimIntersection;
 
            // intersection is returned as a point cloud for dim >= 2
            // so we have to triangulate first
            if constexpr (dimIntersection >= 2)
            {
                const auto triangulation = triangulate<dimIntersection, dimworld>(intersection);
                for (unsigned int i = 0; i < triangulation.size(); ++i)
                    intersections.emplace_back(eIdxA, eIdxB, std::move(triangulation[i]));
            }
            else
                intersections.emplace_back(eIdxA, eIdxB, intersection);
        }
    }
 
    // if we reached the leaf in treeA, just continue in treeB
    else if (isLeafA)
    {
        intersectingEntities(treeA, treeB, nodeA, bBoxB.child0, intersections);
        intersectingEntities(treeA, treeB, nodeA, bBoxB.child1, intersections);
    }
 
    // if we reached the leaf in treeB, just continue in treeA
    else if (isLeafB)
    {
        intersectingEntities(treeA, treeB, bBoxA.child0, nodeB, intersections);
        intersectingEntities(treeA, treeB, bBoxA.child1, nodeB, intersections);
    }
 
    // we know now that both trees didn't reach the leaf yet so
    // we continue with the larger tree first (bigger node number)
    else if (nodeA > nodeB)
    {
        intersectingEntities(treeA, treeB, bBoxA.child0, nodeB, intersections);
        intersectingEntities(treeA, treeB, bBoxA.child1, nodeB, intersections);
    }
    else
    {
        intersectingEntities(treeA, treeB, nodeA, bBoxB.child0, intersections);
        intersectingEntities(treeA, treeB, nodeA, bBoxB.child1, intersections);
    }
}

template<class ctype, int dimworld>
inline std::size_t intersectingEntityCartesianGrid(const Dune::FieldVector<ctype, dimworld>& point,
                                                   const Dune::FieldVector<ctype, dimworld>& min,
                                                   const Dune::FieldVector<ctype, dimworld>& max,
                                                   const std::array<int, std::size_t(dimworld)>& cells)
{
    std::size_t index = 0;
    for (int i = 0; i < dimworld; ++i)
    {
        using std::clamp; using std::floor;
        ctype dimOffset = clamp<ctype>(floor((point[i]-min[i])*cells[i]/(max[i]-min[i])), 0.0, cells[i]-1);
        for (int j = 0; j < i; ++j)
            dimOffset *= cells[j];
        index += static_cast<std::size_t>(dimOffset);
    }
    return index;
}
 
} // end namespace Dumux
 
#endif