From 10c6a07fba102041b8d230dbdeb686f3c40f98b6 Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Thu, 29 Jul 2021 11:42:09 +0200
Subject: [PATCH] [geometry][distance] Add convenience function for squared
 distance

---
 dumux/geometry/distance.hh | 32 ++++++++++++++++++++------------
 1 file changed, 20 insertions(+), 12 deletions(-)

diff --git a/dumux/geometry/distance.hh b/dumux/geometry/distance.hh
index 20530e4599..a9f110bf8a 100644
--- a/dumux/geometry/distance.hh
+++ b/dumux/geometry/distance.hh
@@ -303,7 +303,7 @@ namespace Detail {
 template<class Geo1, class Geo2,
          int dimWorld = Geo1::coorddimension,
          int dim1 = Geo1::mydimension, int dim2 = Geo2::mydimension>
-struct GeometryDistance
+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)
@@ -315,58 +315,66 @@ struct GeometryDistance
 
 // distance point-point
 template<class Geo1, class Geo2, int dimWorld>
-struct GeometryDistance<Geo1, Geo2, dimWorld, 0, 0>
+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::distance(geo1.corner(0), geo2.corner(0)); }
+    { return Dumux::squaredDistance(geo1.corner(0), geo2.corner(0)); }
 };
 
 // distance segment-point
 template<class Geo1, class Geo2, int dimWorld>
-struct GeometryDistance<Geo1, Geo2, dimWorld, 1, 0>
+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 distancePointSegment(geo2.corner(0), geo1); }
+    { return squaredDistancePointSegment(geo2.corner(0), geo1); }
 };
 
 // distance point-segment
 template<class Geo1, class Geo2, int dimWorld>
-struct GeometryDistance<Geo1, Geo2, dimWorld, 0, 1>
+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 distancePointSegment(geo1.corner(0), geo2); }
+    { return squaredDistancePointSegment(geo1.corner(0), geo2); }
 };
 
 // distance point-polygon
 template<class Geo1, class Geo2, int dimWorld>
-struct GeometryDistance<Geo1, Geo2, dimWorld, 0, 2>
+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 distancePointPolygon(geo1.corner(0), geo2); }
+    { return squaredDistancePointPolygon(geo1.corner(0), geo2); }
 };
 
 // distance polygon-point
 template<class Geo1, class Geo2, int dimWorld>
-struct GeometryDistance<Geo1, Geo2, dimWorld, 2, 0>
+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 distancePointPolygon(geo2.corner(0), geo1); }
+    { return squaredDistancePointPolygon(geo2.corner(0), geo1); }
 };
 
 } // end namespace Detail
 
+/*!
+ * \ingroup Geometry
+ * \brief Compute the shortest distance between two geometries
+ */
+template<class Geo1, class Geo2>
+static inline auto squaredDistance(const Geo1& geo1, const Geo2& geo2)
+{ return Detail::GeometrySquaredDistance<Geo1, Geo2>::distance(geo1, geo2); }
+
 /*!
  * \ingroup Geometry
  * \brief Compute the shortest distance between two geometries
  */
 template<class Geo1, class Geo2>
 static inline auto distance(const Geo1& geo1, const Geo2& geo2)
-{ return Detail::GeometryDistance<Geo1, Geo2>::distance(geo1, geo2); }
+{ using std::sqrt; return sqrt(squaredDistance(geo1, geo2)); }
 
 } // end namespace Dumux
 
-- 
GitLab