diff --git a/dumux/geometry/boundingboxtree.hh b/dumux/geometry/boundingboxtree.hh
index 85888341c424477bc2f1ff3c7d8101223999f34f..b9afa8ce463dbaf254e688859a77cf0cd625fd42 100644
--- a/dumux/geometry/boundingboxtree.hh
+++ b/dumux/geometry/boundingboxtree.hh
@@ -411,6 +411,25 @@ inline bool intersectsBoundingBoxBoundingBox(const ctypea* a, const ctypeb* b)
return b[0] - eps0 <= a[1] && a[0] <= b[1] + eps0;
}
+/*!
+ * \brief Compute squared distance between point and bounding box
+ * \param point The point
+ * \param b Pointer to bounding box coordinates
+ */
+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
diff --git a/dumux/geometry/distance.hh b/dumux/geometry/distance.hh
index 84b2b614973d7e01c896a2829b4f46ac5d519991..faa1b7840411765a62f1e487d6d4922da90436e2 100644
--- a/dumux/geometry/distance.hh
+++ b/dumux/geometry/distance.hh
@@ -24,7 +24,9 @@
#include <dune/common/fvector.hh>
#include <dune/geometry/quadraturerules.hh>
+
#include <dumux/common/math.hh>
+#include <dumux/geometry/boundingboxtree.hh>
namespace Dumux {
@@ -380,6 +382,160 @@ template<class Geo1, class Geo2>
static inline auto distance(const Geo1& geo1, const Geo2& geo2)
{ using std::sqrt; return sqrt(squaredDistance(geo1, geo2)); }
+
+//! Distance implementation details
+namespace Detail {
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the closest entity in an AABB tree (index and shortest squared distance) recursively
+ * \note specialization for geometries with dimension larger than points
+ */
+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);
+ }
+}
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the closest entity in an AABB tree (index and shortest squared distance) recursively
+ * \note specialization for point geometries (point cloud AABB tree)
+ */
+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
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the closest entity in an AABB tree to a point (index and shortest squared distance)
+ * \param point the point
+ * \param tree the AABB tree
+ * \param minSquaredDistance conservative estimate of the minimum distance
+ *
+ * \note it is important that if an estimate is provided for minSquaredDistance to choose
+ * this estimate to be larger than the expected result. If the minSquaredDistance is smaller
+ * or equal to the result the returned entity index will be zero and the distance is equal
+ * to the estimate. However, this can also be the correct result. When in doubt, use the
+ * default parameter value.
+ */
+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 };
+}
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the shortest squared distance to entities in an AABB tree
+ */
+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;
+}
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the shortest distance to entities in an AABB tree
+ */
+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