boundingboxtree.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_BOUNDINGBOXTREE_HH
#define DUMUX_GEOMETRY_BOUNDINGBOXTREE_HH
 
#include <vector>
#include <array>
#include <algorithm>
#include <memory>
#include <numeric>
#include <limits>
#include <type_traits>
#include <iostream>
 
#include <dune/common/promotiontraits.hh>
#include <dune/common/fvector.hh>
 
#ifndef DOXYGEN
namespace Dumux::Detail {
 
template<typename ctype>
inline constexpr ctype minimumBaseEpsilon = 10.0*std::numeric_limits<ctype>::epsilon();
 
}  // namespace Dumux::Detail
#endif  // DOXYGEN
 
namespace Dumux {

template <int dimworld, class ctype, class Geometry>
  requires requires(const Geometry& g) { g.corner(0); g.corners(); }
inline void computeGeometryBoundingBox(ctype* b, const Geometry& geometry)
{
    // get the bounding box coordinates
    ctype* xMin = b;
    ctype* xMax = b + dimworld;
 
    // 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]);
        }
    }
}

template <int dimworld, class ctype, class Entity>
  requires requires(const Entity& e) { e.geometry(); }
inline void computeEntityBoundingBox(ctype* b, const Entity& entity)
{ computeGeometryBoundingBox<dimworld>(b, entity.geometry()); }

template <class GeometricEntitySet>
class BoundingBoxTree
{
    static constexpr int dimworld = GeometricEntitySet::dimensionworld;
    using ctype = typename GeometricEntitySet::ctype;

    struct BoundingBoxNode
    {
        std::size_t child0;
        std::size_t child1;
    };
 
public:
    using EntitySet = GeometricEntitySet;

    BoundingBoxTree() = default;

    explicit BoundingBoxTree(std::shared_ptr<const GeometricEntitySet> set)
    { build(set); }

    void build(std::shared_ptr<const GeometricEntitySet> set)
    {
        // set the pointer to the entity set
        entitySet_ = set;
 
        // clear all internal data
        boundingBoxNodes_.clear();
        boundingBoxCoordinates_.clear();
 
        // Create bounding boxes for all elements
        const auto numLeaves = set->size();
 
        // reserve enough space for the nodes and the coordinates
        const auto numNodes = 2*numLeaves - 1;
        boundingBoxNodes_.reserve(numNodes);
        boundingBoxCoordinates_.reserve(numNodes*2*dimworld);
 
        // create a vector for leaf boxes (min and max for all dims)
        std::vector<ctype> leafBoxes(2*dimworld*numLeaves);
 
        for (const auto& geometricEntity : *set)
            computeEntityBoundingBox<dimworld>(leafBoxes.data() + 2*dimworld*set->index(geometricEntity), geometricEntity);
 
        // create the leaf partition, the set of available indices (to be sorted)
        std::vector<std::size_t> leafPartition(numLeaves);
        std::iota(leafPartition.begin(), leafPartition.end(), 0);
 
        // Recursively build the bounding box tree
        build_(leafBoxes, leafPartition.begin(), leafPartition.end());
    }

    const EntitySet& entitySet() const
    { return *entitySet_; }


    const BoundingBoxNode& getBoundingBoxNode(std::size_t nodeIdx) const
    { return boundingBoxNodes_[nodeIdx]; }

    const ctype* getBoundingBoxCoordinates(std::size_t nodeIdx) const
    { return boundingBoxCoordinates_.data() + 2*dimworld*nodeIdx; }

    std::size_t numBoundingBoxes() const
    { return boundingBoxNodes_.size(); }

    bool isLeaf(const BoundingBoxNode& node, std::size_t nodeIdx) const
    { return node.child0 == nodeIdx; }
 
private:

    std::vector<BoundingBoxNode> boundingBoxNodes_;

    std::vector<ctype> boundingBoxCoordinates_;

    std::shared_ptr<const EntitySet> entitySet_;

    std::size_t build_(const std::vector<ctype>& leafBoxes,
                       const std::vector<std::size_t>::iterator& begin,
                       const std::vector<std::size_t>::iterator& end)
    {
        assert(begin < end);
 
        // If we reached the end of the recursion, i.e. only a leaf box is left
        if (end - begin == 1)
        {
            // Get the bounding box coordinates for the leaf
            const std::size_t leafNodeIdx = *begin;
            const auto beginCoords = leafBoxes.begin() + 2*dimworld*leafNodeIdx;
            const auto endCoords = beginCoords + 2*dimworld;
 
            // Store the data in the bounding box
            // leaf nodes are indicated by setting child0 to
            // the node itself and child1 to the index of the entity in the bounding box.
            return addBoundingBox_(BoundingBoxNode{numBoundingBoxes(), leafNodeIdx}, beginCoords, endCoords);
        }
 
        // Compute the bounding box of all bounding boxes in the range [begin, end]
        const auto bCoords = computeBBoxOfBBoxes_(leafBoxes, begin, end);
 
        // sort bounding boxes along the longest axis
        const auto axis = computeLongestAxis_(bCoords);
 
        // nth_element sorts the range to make sure that middle points to the coordinate median in axis direction
        // this is the most expensive part of the algorithm
        auto middle = begin + (end - begin)/2;
        std::nth_element(begin, middle, end, [&leafBoxes, axis](std::size_t i, std::size_t j)
                         {
                             const ctype* bi = leafBoxes.data() + 2*dimworld*i;
                             const ctype* bj = leafBoxes.data() + 2*dimworld*j;
                             return bi[axis] + bi[axis + dimworld] < bj[axis] + bj[axis + dimworld];
                         });
 
        // split the bounding boxes into two at the middle iterator and call build recursively, each
        // call resulting in a new node of this bounding box, i.e. the root will be added at the end of the process.
        return addBoundingBox_(BoundingBoxNode{build_(leafBoxes, begin, middle), build_(leafBoxes, middle, end)},
                               bCoords.begin(), bCoords.end());
    }

    template <class Iterator>
    std::size_t addBoundingBox_(BoundingBoxNode&& node,
                                const Iterator& coordBegin,
                                const Iterator& coordEnd)
    {
        // Add the bounding box
        boundingBoxNodes_.emplace_back(node);
 
        // Add the bounding box's coordinates
        boundingBoxCoordinates_.insert(boundingBoxCoordinates_.end(), coordBegin, coordEnd);
 
        // return the index of the new node
        return boundingBoxNodes_.size() - 1;
    }

    std::array<ctype, 2*dimworld>
    computeBBoxOfBBoxes_(const std::vector<ctype>& leafBoxes,
                         const std::vector<std::size_t>::iterator& begin,
                         const std::vector<std::size_t>::iterator& end)
    {
        std::array<ctype, 2*dimworld> bBoxCoords;
 
        // copy the iterator and get coordinates of first box
        auto it = begin;
        const auto* bFirst = leafBoxes.data() + 2*dimworld*(*it);
 
        for (int coordIdx = 0; coordIdx < 2*dimworld; ++coordIdx)
            bBoxCoords[coordIdx] = bFirst[coordIdx];
 
        // Compute min and max with the remaining boxes
        for (; it != end; ++it)
        {
            const auto* b = leafBoxes.data() + 2*dimworld*(*it);
            for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
                if (b[coordIdx] < bBoxCoords[coordIdx])
                    bBoxCoords[coordIdx] = b[coordIdx];
            for (int coordIdx = dimworld; coordIdx < 2*dimworld; ++coordIdx)
                if (b[coordIdx] > bBoxCoords[coordIdx])
                    bBoxCoords[coordIdx] = b[coordIdx];
        }
 
        return bBoxCoords;
    }

    std::size_t computeLongestAxis_(const std::array<ctype, 2*dimworld>& bCoords)
    {
        std::array<ctype, dimworld> axisLength;
        for (int coordIdx = 0; coordIdx < dimworld; ++coordIdx)
            axisLength[coordIdx] = bCoords[dimworld + coordIdx] - bCoords[coordIdx];
 
        return std::distance(axisLength.begin(), std::max_element(axisLength.begin(), axisLength.end()));
    }
};

template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 3, int> = 0>
inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
{
    static constexpr ctype eps_ = 1.0e-7;
 
    using std::max;
    const auto dx = b[3] - b[0];
    const auto dy = b[4] - b[1];
    const auto dz = b[5] - b[2];
    const ctype eps = max({dx, dy, dz})*eps_;
    return (b[0] - eps <= point[0] && point[0] <= b[3] + eps &&
            b[1] - eps <= point[1] && point[1] <= b[4] + eps &&
            b[2] - eps <= point[2] && point[2] <= b[5] + eps);
}

template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 2, int> = 0>
inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
{
    static constexpr ctype eps_ = 1.0e-7;
 
    using std::max;
    const auto dx = b[2] - b[0];
    const auto dy = b[3] - b[1];
    const ctype eps = max(dx, dy)*eps_;
    return (b[0] - eps <= point[0] && point[0] <= b[2] + eps &&
            b[1] - eps <= point[1] && point[1] <= b[3] + eps);
}

template<class ctype, int dimworld, typename std::enable_if_t<dimworld == 1, int> = 0>
inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
{
    static constexpr ctype eps_ = 1.0e-7;
    const ctype eps0 = eps_*(b[1] - b[0]);
    return b[0] - eps0 <= point[0] && point[0] <= b[1] + eps0;
}

template<class ctype, int dimworld>
inline bool intersectsPointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point,
                                       const Dune::FieldVector<ctype, dimworld>& min,
                                       const Dune::FieldVector<ctype, dimworld>& max)
{
    std::array<ctype, 2*dimworld> bBox;
    std::copy(min.begin(), min.end(), bBox.begin());
    std::copy(max.begin(), max.end(), bBox.begin()+dimworld);
    return intersectsPointBoundingBox(point, bBox.data());
}

template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 3, int> = 0>
inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
{
    using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
    static constexpr ctype eps_ = 1.0e-7;
    const ctype scale0 = std::max(b[3]-b[0], a[3]-a[0]);
    const ctype scale1 = std::max(b[4]-b[1], a[4]-a[1]);
    const ctype scale2 = std::max(b[5]-b[2], a[5]-a[2]);
    const ctype maxScale = std::max(scale0, std::max(scale1, scale2));
    const ctype minEps = maxScale*Detail::minimumBaseEpsilon<ctype>;
    const ctype eps0 = std::max(eps_*scale0, minEps);
    const ctype eps1 = std::max(eps_*scale1, minEps);
    const ctype eps2 = std::max(eps_*scale2, minEps);
    return (b[0] - eps0 <= a[3] && a[0] <= b[3] + eps0 &&
            b[1] - eps1 <= a[4] && a[1] <= b[4] + eps1 &&
            b[2] - eps2 <= a[5] && a[2] <= b[5] + eps2);
 
}

template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 2, int> = 0>
inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
{
    using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
    static constexpr ctype eps_ = 1.0e-7;
    const ctype scale0 = std::max(b[2]-b[0], a[2]-a[0]);
    const ctype scale1 = std::max(b[3]-b[1], a[3]-a[1]);
    const ctype maxScale = std::max(scale0, scale1);
    const ctype minEps = maxScale*Detail::minimumBaseEpsilon<ctype>;
    const ctype eps0 = std::max(eps_*scale0, minEps);
    const ctype eps1 = std::max(eps_*scale1, minEps);
    return (b[0] - eps0 <= a[2] && a[0] <= b[2] + eps0 &&
            b[1] - eps1 <= a[3] && a[1] <= b[3] + eps1);
}

template<int dimworld, class ctypea, class ctypeb, typename std::enable_if_t<dimworld == 1, int> = 0>
inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
{
    using ctype = typename Dune::PromotionTraits<ctypea, ctypeb>::PromotedType;
    static constexpr ctype eps_ = 1.0e-7;
    const ctype scale0 = std::max(b[1]-b[0], a[1]-a[0]);
    const ctype eps0 = std::max(eps_*scale0, Detail::minimumBaseEpsilon<ctype>*scale0);
    return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0;
}

template<int dimworld, class ctype>
inline ctype squaredDistancePointBoundingBox(const Dune::FieldVector<ctype, dimworld>& point, const ctype* b)
{
    ctype squaredDistance = 0.0;
    for (int d = 0; d < dimworld; ++d)
    {
        if (point[d] < b[d])
            squaredDistance += (point[d] - b[d])*(point[d] - b[d]);
        if (point[d] > b[d+dimworld])
            squaredDistance += (point[d] - b[d+dimworld])*(point[d] - b[d+dimworld]);
    }
    return squaredDistance;
}
 
} // end namespace Dumux
 
#endif