GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/geometry/distance.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 128 129 99.2%
Functions: 103 106 97.2%
Branches: 108 206 52.4%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Geometry
10 * \brief Helper functions for distance queries
11 */
12 #ifndef DUMUX_GEOMETRY_DISTANCE_HH
13 #define DUMUX_GEOMETRY_DISTANCE_HH
14
15 #include <dune/common/fvector.hh>
16 #include <dune/geometry/quadraturerules.hh>
17
18 #include <dumux/common/math.hh>
19 #include <dumux/geometry/boundingboxtree.hh>
20
21 namespace Dumux {
22
23 /*!
24 * \ingroup Geometry
25 * \brief Compute the average distance from a point to a geometry by integration
26 */
27 template<class Geometry>
28 static inline typename Geometry::ctype
29 33058 averageDistancePointGeometry(const typename Geometry::GlobalCoordinate& p,
30 const Geometry& geometry,
31 std::size_t integrationOrder = 2)
32 {
33 33058 typename Geometry::ctype avgDist = 0.0;
34 33058 const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
35
3/3
✓ Branch 0 taken 246984 times.
✓ Branch 1 taken 48353 times.
✓ Branch 2 taken 2185 times.
 297522 for (const auto& qp : quad)
36 1287360 avgDist += (geometry.global(qp.position())-p).two_norm()*qp.weight()*geometry.integrationElement(qp.position());
37 33058 return avgDist/geometry.volume();
38 }
39
40 /*!
41 * \ingroup Geometry
42 * \brief Compute the squared distance from a point to a line through the points a and b
43 */
44 template<class Point>
45 static inline typename Point::value_type
46 2880 squaredDistancePointLine(const Point& p, const Point& a, const Point& b)
47 {
48 2880 const auto ab = b - a;
49
2/2
✓ Branch 0 taken 3600 times.
✓ Branch 1 taken 1440 times.
 12960 const auto t = (p - a)*ab/ab.two_norm2();
50 2880 auto proj = a;
51 2880 proj.axpy(t, ab);
52 2880 return (proj - p).two_norm2();
53 }
54
55 /*!
56 * \ingroup Geometry
57 * \brief Compute the squared distance from a point to a line given by a geometry with two corners
58 * \note We currently lack the a representation of a line geometry. This convenience function
59 * assumes a segment geometry (with two corners) is passed which represents a line geometry.
60 */
61 template<class Geometry>
62 static inline typename Geometry::ctype
63 squaredDistancePointLine(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
64 {
65 static_assert(Geometry::mydimension == 1, "Geometry has to be a line");
66 const auto& a = geometry.corner(0);
67 const auto& b = geometry.corner(1);
68 return squaredDistancePointLine(p, a, b);
69 }
70
71 /*!
72 * \ingroup Geometry
73 * \brief Compute the distance from a point to a line through the points a and b
74 */
75 template<class Point>
76 static inline typename Point::value_type
77 1440 distancePointLine(const Point& p, const Point& a, const Point& b)
78
24/48
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 60 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 60 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 60 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 60 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 60 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 60 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 60 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 60 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 60 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 60 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 60 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 60 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 60 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 60 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 60 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 60 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 60 times.
✗ Branch 56 not taken.
✓ Branch 58 taken 60 times.
✗ Branch 59 not taken.
✓ Branch 61 taken 60 times.
✗ Branch 62 not taken.
✓ Branch 64 taken 60 times.
✗ Branch 65 not taken.
✓ Branch 67 taken 60 times.
✗ Branch 68 not taken.
✓ Branch 70 taken 60 times.
✗ Branch 71 not taken.
 1440 { using std::sqrt; return sqrt(squaredDistancePointLine(p, a, b)); }
79
80 /*!
81 * \ingroup Geometry
82 * \brief Compute the distance from a point to a line given by a geometry with two corners
83 * \note We currently lack the a representation of a line geometry. This convenience function
84 * assumes a segment geometry (with two corners) is passed which represents a line geometry.
85 */
86 template<class Geometry>
87 static inline typename Geometry::ctype
88 distancePointLine(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
89 { using std::sqrt; return sqrt(squaredDistancePointLine(p, geometry)); }
90
91 /*!
92 * \ingroup Geometry
93 * \brief Compute the squared distance from a point to the segment connecting the points a and b
94 */
95 template<class Point>
96 static inline typename Point::value_type
97 4809417 squaredDistancePointSegment(const Point& p, const Point& a, const Point& b)
98 {
99 4809417 const auto ab = b - a;
100 4809417 const auto ap = p - a;
101 4809417 const auto t = ap*ab;
102
103
2/2
✓ Branch 0 taken 1687387 times.
✓ Branch 1 taken 722405 times.
 4809417 if (t <= 0.0)
104 4809417 return ap.two_norm2();
105
106 3364607 const auto lengthSq = ab.two_norm2();
107
2/2
✓ Branch 0 taken 855335 times.
✓ Branch 1 taken 832052 times.
 3364607 if (t >= lengthSq)
108 3421340 return (b - p).two_norm2();
109
110 1653937 auto proj = a;
111 3307874 proj.axpy(t/lengthSq, ab);
112 3307874 return (proj - p).two_norm2();
113 }
114
115 /*!
116 * \ingroup Geometry
117 * \brief Compute the squared distance from a point to a given segment geometry
118 */
119 template<class Geometry>
120 static inline typename Geometry::ctype
121 11726 squaredDistancePointSegment(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
122 {
123 static_assert(Geometry::mydimension == 1, "Geometry has to be a segment");
124 11726 const auto& a = geometry.corner(0);
125 11726 const auto& b = geometry.corner(1);
126 11726 return squaredDistancePointSegment(p, a, b);
127 }
128
129 /*!
130 * \ingroup Geometry
131 * \brief Compute the distance from a point to the segment connecting the points a and b
132 */
133 template<class Point>
134 static inline typename Point::value_type
135 3040 distancePointSegment(const Point& p, const Point& a, const Point& b)
136 3040 { using std::sqrt; return sqrt(squaredDistancePointSegment(p, a, b)); }
137
138 /*!
139 * \ingroup Geometry
140 * \brief Compute the distance from a point to a given segment geometry
141 */
142 template<class Geometry>
143 static inline typename Geometry::ctype
144 28 distancePointSegment(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
145 28 { using std::sqrt; return sqrt(squaredDistancePointSegment(p, geometry)); }
146
147 /*!
148 * \ingroup Geometry
149 * \brief Compute the shortest squared distance from a point to the triangle connecting the points a, b and c
150 * See https://www.iquilezles.org/www/articles/triangledistance/triangledistance.htm.
151 */
152 template<class Point>
153 static inline typename Point::value_type
154 845441 squaredDistancePointTriangle(const Point& p, const Point& a, const Point& b, const Point& c)
155 {
156 static_assert(Point::dimension == 3, "Only works in 3D");
157 845441 const auto ab = b - a;
158 845441 const auto bc = c - b;
159 845441 const auto ca = a - c;
160 845441 const auto normal = crossProduct(ab, ca);
161
162 845441 const auto ap = p - a;
163 845441 const auto bp = p - b;
164 845441 const auto cp = p - c;
165
166 1690882 const auto sum = sign(crossProduct(ab, normal)*ap)
167 1690882 + sign(crossProduct(bc, normal)*bp)
168
2/2
✓ Branch 1 taken 799023 times.
✓ Branch 2 taken 46418 times.
 1690882 + sign(crossProduct(ca, normal)*cp);
169
170 // if there is no orthogonal projection
171 // (point is outside the infinite prism implied by the triangle and its surface normal)
172 // compute distance to the edges (codim-1 facets)
173
2/2
✓ Branch 0 taken 799023 times.
✓ Branch 1 taken 46418 times.
 845441 if (sum < 2.0)
174 {
175 using std::min;
176 799023 return min({squaredDistancePointSegment(p, a, b),
177 799023 squaredDistancePointSegment(p, a, c),
178 799023 squaredDistancePointSegment(p, b, c)});
179 }
180 // compute distance via orthogonal projection
181 else
182 {
183 46418 const auto tmp = normal*ap;
184 92836 return tmp*tmp / (normal*normal);
185 }
186 }
187
188 /*!
189 * \ingroup Geometry
190 * \brief Compute the shortest squared distance from a point to a given triangle geometry
191 */
192 template<class Geometry>
193 static inline typename Geometry::ctype
194
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 213155 times.
 428154 squaredDistancePointTriangle(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
195 {
196 static_assert(Geometry::coorddimension == 3, "Only works in 3D");
197 static_assert(Geometry::mydimension == 2, "Geometry has to be a triangle");
198
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 214077 times.
 428154 assert(geometry.corners() == 3);
199 428154 const auto& a = geometry.corner(0);
200 428154 const auto& b = geometry.corner(1);
201 428154 const auto& c = geometry.corner(2);
202 428154 return squaredDistancePointTriangle(p, a, b, c);
203 }
204
205 /*!
206 * \ingroup Geometry
207 * \brief Compute the shortest distance from a point to the triangle connecting the points a, b and c
208 */
209 template<class Point>
210 static inline typename Point::value_type
211 480 distancePointTriangle(const Point& p, const Point& a, const Point& b, const Point& c)
212
8/16
✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 60 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 60 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 60 times.
✗ Branch 15 not taken.
✓ Branch 18 taken 60 times.
✗ Branch 19 not taken.
✓ Branch 22 taken 60 times.
✗ Branch 23 not taken.
✓ Branch 26 taken 60 times.
✗ Branch 27 not taken.
✓ Branch 30 taken 60 times.
✗ Branch 31 not taken.
 480 { using std::sqrt; return sqrt(squaredDistancePointTriangle(p, a, b, c)); }
213
214 /*!
215 * \ingroup Geometry
216 * \brief Compute the shortest distance from a point to a given triangle geometry
217 */
218 template<class Geometry>
219 static inline typename Geometry::ctype
220 240 distancePointTriangle(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
221
4/8
✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 60 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 60 times.
✗ Branch 11 not taken.
✓ Branch 14 taken 60 times.
✗ Branch 15 not taken.
 240 { using std::sqrt; return sqrt(squaredDistancePointTriangle(p, geometry)); }
222
223 /*!
224 * \ingroup Geometry
225 * \brief Compute the shortest squared distance from a point to a given polygon geometry
226 * \note We only support triangles and quadrilaterals so far.
227 */
228 template<class Geometry>
229 static inline typename Geometry::ctype
230
2/2
✓ Branch 0 taken 213155 times.
✓ Branch 1 taken 314487 times.
 1058537 squaredDistancePointPolygon(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
231 {
232
2/2
✓ Branch 0 taken 213837 times.
✓ Branch 1 taken 315421 times.
 1058516 if (geometry.corners() == 3)
233 427674 return squaredDistancePointTriangle(p, geometry);
234
1/2
✓ Branch 0 taken 315421 times.
✗ Branch 1 not taken.
 630842 else if (geometry.corners() == 4)
235 {
236 630863 const auto& a = geometry.corner(0);
237 630863 const auto& b = geometry.corner(1);
238 630863 const auto& c = geometry.corner(2);
239 630863 const auto& d = geometry.corner(3);
240
241 using std::min;
242
2/2
✓ Branch 0 taken 115601 times.
✓ Branch 1 taken 199841 times.
 630863 return min(squaredDistancePointTriangle(p, a, b, d),
243 1261705 squaredDistancePointTriangle(p, a, d, c));
244 }
245 else
246 ✗ DUNE_THROW(Dune::NotImplemented, "Polygon with " << geometry.corners() << " corners not supported");
247 }
248
249 /*!
250 * \ingroup Geometry
251 * \brief Compute the shortest distance from a point to a given polygon geometry
252 * \note We only support triangles and quadrilaterals so far.
253 */
254 template<class Geometry>
255 static inline typename Geometry::ctype
256 21 distancePointPolygon(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
257 21 { using std::sqrt; return sqrt(squaredDistancePointPolygon(p, geometry)); }
258
259 /*!
260 * \ingroup Geometry
261 * \brief Compute the average distance from a segment to a geometry by integration
262 */
263 template<class Geometry>
264 static inline typename Geometry::ctype
265 380 averageDistanceSegmentGeometry(const typename Geometry::GlobalCoordinate& a,
266 const typename Geometry::GlobalCoordinate& b,
267 const Geometry& geometry,
268 std::size_t integrationOrder = 2)
269 {
270 380 typename Geometry::ctype avgDist = 0.0;
271 380 const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
272
2/2
✓ Branch 0 taken 3040 times.
✓ Branch 1 taken 380 times.
 3420 for (const auto& qp : quad)
273 9120 avgDist += distancePointSegment(geometry.global(qp.position()), a, b)*qp.weight()*geometry.integrationElement(qp.position());
274 380 return avgDist/geometry.volume();
275 }
276
277 /*!
278 * \ingroup Geometry
279 * \brief Compute the shortest distance between two points
280 */
281 template<class ctype, int dimWorld>
282 static inline ctype distance(const Dune::FieldVector<ctype, dimWorld>& a,
283 const Dune::FieldVector<ctype, dimWorld>& b)
284 { return (a-b).two_norm(); }
285
286 /*!
287 * \ingroup Geometry
288 * \brief Compute the shortest squared distance between two points
289 */
290 template<class ctype, int dimWorld>
291 960 static inline ctype squaredDistance(const Dune::FieldVector<ctype, dimWorld>& a,
292 const Dune::FieldVector<ctype, dimWorld>& b)
293 960 { return (a-b).two_norm2(); }
294
295
296 namespace Detail {
297
298 // helper struct to compute distance between two geometries, specialized below
299 template<class Geo1, class Geo2,
300 int dimWorld = Geo1::coorddimension,
301 int dim1 = Geo1::mydimension, int dim2 = Geo2::mydimension>
302 struct GeometrySquaredDistance
303 {
304 static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
305 static auto distance(const Geo1& geo1, const Geo2& geo2)
306 {
307 DUNE_THROW(Dune::NotImplemented, "Geometry distance computation not implemented for dimworld = "
308 << dimWorld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
309 }
310 };
311
312 // distance point-point
313 template<class Geo1, class Geo2, int dimWorld>
314 struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 0, 0>
315 {
316 static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
317 840 static auto distance(const Geo1& geo1, const Geo2& geo2)
318 840 { return Dumux::squaredDistance(geo1.corner(0), geo2.corner(0)); }
319 };
320
321 // distance segment-point
322 template<class Geo1, class Geo2, int dimWorld>
323 struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 1, 0>
324 {
325 static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
326 1800 static auto distance(const Geo1& geo1, const Geo2& geo2)
327 1260 { return squaredDistancePointSegment(geo2.corner(0), geo1); }
328 };
329
330 // distance point-segment
331 template<class Geo1, class Geo2, int dimWorld>
332 struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 0, 1>
333 {
334 static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
335 static auto distance(const Geo1& geo1, const Geo2& geo2)
336 { return squaredDistancePointSegment(geo1.corner(0), geo2); }
337 };
338
339 // distance point-polygon
340 template<class Geo1, class Geo2, int dimWorld>
341 struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 0, 2>
342 {
343 static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
344 720 static inline auto distance(const Geo1& geo1, const Geo2& geo2)
345 720 { return squaredDistancePointPolygon(geo1.corner(0), geo2); }
346 };
347
348 // distance polygon-point
349 template<class Geo1, class Geo2, int dimWorld>
350 struct GeometrySquaredDistance<Geo1, Geo2, dimWorld, 2, 0>
351 {
352 static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
353 720 static inline auto distance(const Geo1& geo1, const Geo2& geo2)
354 720 { return squaredDistancePointPolygon(geo2.corner(0), geo1); }
355 };
356
357 } // end namespace Detail
358
359 /*!
360 * \ingroup Geometry
361 * \brief Compute the shortest distance between two geometries
362 */
363 template<class Geo1, class Geo2>
364 2880 static inline auto squaredDistance(const Geo1& geo1, const Geo2& geo2)
365 2880 { return Detail::GeometrySquaredDistance<Geo1, Geo2>::distance(geo1, geo2); }
366
367 /*!
368 * \ingroup Geometry
369 * \brief Compute the shortest distance between two geometries
370 */
371 template<class Geo1, class Geo2>
372
8/16
✓ Branch 5 taken 60 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 60 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 60 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 60 times.
✗ Branch 15 not taken.
✓ Branch 22 taken 60 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 60 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 60 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 60 times.
✗ Branch 32 not taken.
 4440 static inline auto distance(const Geo1& geo1, const Geo2& geo2)
373
20/40
✓ Branch 1 taken 60 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 60 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 60 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 60 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 60 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 60 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 60 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 60 times.
✗ Branch 23 not taken.
✓ Branch 25 taken 60 times.
✗ Branch 26 not taken.
✓ Branch 28 taken 60 times.
✗ Branch 29 not taken.
✓ Branch 31 taken 60 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 60 times.
✗ Branch 35 not taken.
✓ Branch 37 taken 60 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 60 times.
✗ Branch 41 not taken.
✓ Branch 43 taken 60 times.
✗ Branch 44 not taken.
✓ Branch 46 taken 60 times.
✗ Branch 47 not taken.
✓ Branch 49 taken 60 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 60 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 60 times.
✗ Branch 56 not taken.
✓ Branch 58 taken 60 times.
✗ Branch 59 not taken.
 4560 { using std::sqrt; return sqrt(squaredDistance(geo1, geo2)); }
374
375
376 //! Distance implementation details
377 namespace Detail {
378
379 /*!
380 * \ingroup Geometry
381 * \brief Compute the closest entity in an AABB tree (index and shortest squared distance) recursively
382 * \note specialization for geometries with dimension larger than points
383 */
384 template<class EntitySet, class ctype, int dimworld,
385 typename std::enable_if_t<(EntitySet::Entity::Geometry::mydimension > 0), int> = 0>
386 3641202 void closestEntity(const Dune::FieldVector<ctype, dimworld>& point,
387 const BoundingBoxTree<EntitySet>& tree,
388 std::size_t node,
389 ctype& minSquaredDistance,
390 std::size_t& eIdx)
391 {
392 // Get the bounding box for the current node
393
2/2
✓ Branch 0 taken 2340332 times.
✓ Branch 1 taken 1302094 times.
 7219864 const auto& bBox = tree.getBoundingBoxNode(node);
394
395 // If bounding box is outside radius, then don't search further
396 7219864 const auto squaredDistance = squaredDistancePointBoundingBox(
397 point, tree.getBoundingBoxCoordinates(node)
398 );
399
400 // do not continue if the AABB is further away than the current minimum
401
2/2
✓ Branch 0 taken 2340332 times.
✓ Branch 1 taken 1302094 times.
 7219864 if (squaredDistance > minSquaredDistance) return;
402
403 // If the bounding box is a leaf, update distance and index with primitive test
404
2/2
✓ Branch 0 taken 536513 times.
✓ Branch 1 taken 1803819 times.
 4644652 else if (tree.isLeaf(bBox, node))
405 {
406 1065990 const std::size_t entityIdx = bBox.child1;
407 2668493 const auto squaredDistance = [&]{
408 536513 const auto geometry = tree.entitySet().entity(entityIdx).geometry();
409 if constexpr (EntitySet::Entity::Geometry::mydimension == 2)
410
2/4
✓ Branch 1 taken 336 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 336 times.
✗ Branch 4 not taken.
 528298 return squaredDistancePointPolygon(point, geometry);
411 else if constexpr (EntitySet::Entity::Geometry::mydimension == 1)
412
2/4
✓ Branch 1 taken 145 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158 times.
✗ Branch 5 not taken.
 7915 return squaredDistancePointSegment(point, geometry);
413 else
414 DUNE_THROW(Dune::NotImplemented, "squaredDistance to entity with dim>2");
415 1066326 }();
416
417
2/2
✓ Branch 0 taken 50032 times.
✓ Branch 1 taken 486481 times.
 1065990 if (squaredDistance < minSquaredDistance)
418 {
419 93128 eIdx = entityIdx;
420 93128 minSquaredDistance = squaredDistance;
421 }
422 }
423
424 // Check the children nodes recursively
425 else
426 {
427 3578662 closestEntity(point, tree, bBox.child0, minSquaredDistance, eIdx);
428 3578662 closestEntity(point, tree, bBox.child1, minSquaredDistance, eIdx);
429 }
430 }
431
432 /*!
433 * \ingroup Geometry
434 * \brief Compute the closest entity in an AABB tree (index and shortest squared distance) recursively
435 * \note specialization for point geometries (point cloud AABB tree)
436 */
437 template<class EntitySet, class ctype, int dimworld,
438 typename std::enable_if_t<(EntitySet::Entity::Geometry::mydimension == 0), int> = 0>
439 10520702 void closestEntity(const Dune::FieldVector<ctype, dimworld>& point,
440 const BoundingBoxTree<EntitySet>& tree,
441 std::size_t node,
442 ctype& minSquaredDistance,
443 std::size_t& eIdx)
444 {
445 // Get the bounding box for the current node
446
2/2
✓ Branch 0 taken 3108638 times.
✓ Branch 1 taken 7463176 times.
 20978864 const auto& bBox = tree.getBoundingBoxNode(node);
447
448 // If the bounding box is a leaf, update distance and index with primitive test
449
2/2
✓ Branch 0 taken 3108638 times.
✓ Branch 1 taken 7463176 times.
 20978864 if (tree.isLeaf(bBox, node))
450 {
451 6141746 const std::size_t entityIdx = bBox.child1;
452 6141746 const auto& p = tree.entitySet().entity(entityIdx).geometry().corner(0);
453 6141746 const auto squaredDistance = (p-point).two_norm2();
454
455
2/2
✓ Branch 0 taken 1411151 times.
✓ Branch 1 taken 1697487 times.
 6141746 if (squaredDistance < minSquaredDistance)
456 {
457 2762632 eIdx = entityIdx;
458 2762632 minSquaredDistance = squaredDistance;
459 }
460 }
461
462 // Check the children nodes recursively
463 else
464 {
465 // If bounding box is outside radius, then don't search further
466 14837118 const auto squaredDistance = squaredDistancePointBoundingBox(
467 point, tree.getBoundingBoxCoordinates(node)
468 );
469
470 // do not continue if the AABB is further away than the current minimum
471
2/2
✓ Branch 0 taken 5268513 times.
✓ Branch 1 taken 2194663 times.
 14837118 if (squaredDistance > minSquaredDistance) return;
472
473 10458162 closestEntity(point, tree, bBox.child0, minSquaredDistance, eIdx);
474 10458162 closestEntity(point, tree, bBox.child1, minSquaredDistance, eIdx);
475 }
476 }
477
478 } // end namespace Detail
479
480 /*!
481 * \ingroup Geometry
482 * \brief Compute the closest entity in an AABB tree to a point (index and shortest squared distance)
483 * \param point the point
484 * \param tree the AABB tree
485 * \param minSquaredDistance conservative estimate of the minimum distance
486 *
487 * \note it is important that if an estimate is provided for minSquaredDistance to choose
488 * this estimate to be larger than the expected result. If the minSquaredDistance is smaller
489 * or equal to the result the returned entity index will be zero and the distance is equal
490 * to the estimate. However, this can also be the correct result. When in doubt, use the
491 * default parameter value.
492 */
493 template<class EntitySet, class ctype, int dimworld>
494 139152 std::pair<ctype, std::size_t> closestEntity(const Dune::FieldVector<ctype, dimworld>& point,
495 const BoundingBoxTree<EntitySet>& tree,
496 ctype minSquaredDistance = std::numeric_limits<ctype>::max())
497 {
498 139152 std::size_t eIdx = 0;
499 139152 Detail::closestEntity(point, tree, tree.numBoundingBoxes() - 1, minSquaredDistance, eIdx);
500 using std::sqrt;
501 139152 return { minSquaredDistance, eIdx };
502 }
503
504 /*!
505 * \ingroup Geometry
506 * \brief Compute the shortest squared distance to entities in an AABB tree
507 */
508 template<class EntitySet, class ctype, int dimworld>
509 34788 ctype squaredDistance(const Dune::FieldVector<ctype, dimworld>& point,
510 const BoundingBoxTree<EntitySet>& tree,
511 ctype minSquaredDistance = std::numeric_limits<ctype>::max())
512 {
513 34788 return closestEntity(point, tree, minSquaredDistance).first;
514 }
515
516 /*!
517 * \ingroup Geometry
518 * \brief Compute the shortest distance to entities in an AABB tree
519 */
520 template<class EntitySet, class ctype, int dimworld>
521 ctype distance(const Dune::FieldVector<ctype, dimworld>& point,
522 const BoundingBoxTree<EntitySet>& tree,
523 ctype minSquaredDistance = std::numeric_limits<ctype>::max())
524 {
525 using std::sqrt;
526 return sqrt(squaredDistance(point, tree, minSquaredDistance));
527 }
528
529 } // end namespace Dumux
530
531 #endif
532