distancefield.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_DISTANCE_FIELD_HH
#define DUMUX_GEOMETRY_DISTANCE_FIELD_HH
 
#include <memory>
#include <vector>
#include <utility>
 
#include <dumux/geometry/distance.hh>
#include <dumux/geometry/boundingboxtree.hh>
#include <dumux/geometry/geometricentityset.hh>
 
namespace Dumux {
 
template<class Geometry>
class AABBDistanceField
{
    using Point = typename Geometry::GlobalCoordinate;
    using Scalar = typename Geometry::ctype;
    using GeoSet = GeometriesEntitySet<Geometry>;
    using AABBTree = BoundingBoxTree<GeoSet>;
 
    class PointGeometry
    {
    public:
        static constexpr int mydimension = 0;
        static constexpr int coorddimension = Geometry::coorddimension;
        using ctype = typename Geometry::ctype;
 
        PointGeometry(Point&& point) : point_(std::move(point)) {}
        const Point& corner(std::size_t i) const { assert(i == 0); return point_; }
        std::size_t corners() const { return 1; }
    private:
        Point point_;
    };
 
    using PointGeoSet = GeometriesEntitySet<PointGeometry>;
    using AABBTreeMidPoints = BoundingBoxTree<PointGeoSet>;
 
public:
    AABBDistanceField(const std::vector<Geometry>& geometries)
    : tree_(std::make_unique<AABBTree>(std::make_shared<GeoSet>(geometries)))
    {
        std::vector<PointGeometry> points;
        points.reserve(geometries.size());
        for (const auto& geo : geometries)
        {
            auto center = geo.center();
            points.emplace_back(std::move(center));
        }
 
        pointTree_ = std::make_unique<AABBTreeMidPoints>(
            std::make_shared<PointGeoSet>(std::move(points))
        );
    }
 
    std::pair<Scalar, std::size_t> distanceAndIndex(const Point& p) const
    { return distanceAndIndex_(p); }
 
    Scalar distance(const Point& p) const
    { return distanceAndIndex_(p).first; }
 
private:
    std::pair<Scalar, std::size_t> distanceAndIndex_(const Point& p) const
    {
        // find a good first guess by checking the mid point tree
        const auto minSquaredDistanceEstimate = squaredDistance(p, *pointTree_);
 
        // find actual entity and distance to the entity's geometry
        // we choose the distance estimate a bit larger in case it actually is already the minimum distance
        const auto [squaredDistance, index] = closestEntity(p, *tree_, minSquaredDistanceEstimate*1.00001);
 
        using std::sqrt;
        return { sqrt(squaredDistance), index };
    }
 
    std::unique_ptr<AABBTree> tree_;
    std::unique_ptr<AABBTreeMidPoints> pointTree_;
};
 
template<class Geometry>
using DistanceField = AABBDistanceField<Geometry>;
 
} // end namespace Dumux
 
#endif