GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/geometry/distance.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 119 121 98.3%
Functions: 95 108 88.0%
Branches: 110 198 55.6%
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-FileCopyrightInfo: 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 63931 const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
35
4/4
✓ Branch 0 taken 48353 times.
✓ Branch 1 taken 249169 times.
✓ Branch 2 taken 48353 times.
✓ Branch 3 taken 249169 times.
 628102 for (const auto& qp : quad)
36 1534344 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 1440 squaredDistancePointLine(const Point& p, const Point& a, const Point& b)
47 {
48 1440 const auto ab = b - a;
49 1440 const auto t = (p - a)*ab/ab.two_norm2();
50 1440 auto proj = a;
51 1440 proj.axpy(t, ab);
52 1440 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 ✗ distancePointLine(const Point& p, const Point& a, const Point& b)
78
24/48
✓ Branch 3 taken 60 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 60 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 60 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 60 times.
✗ Branch 16 not taken.
✓ Branch 19 taken 60 times.
✗ Branch 20 not taken.
✓ Branch 23 taken 60 times.
✗ Branch 24 not taken.
✓ Branch 27 taken 60 times.
✗ Branch 28 not taken.
✓ Branch 31 taken 60 times.
✗ Branch 32 not taken.
✓ Branch 35 taken 60 times.
✗ Branch 36 not taken.
✓ Branch 39 taken 60 times.
✗ Branch 40 not taken.
✓ Branch 43 taken 60 times.
✗ Branch 44 not taken.
✓ Branch 47 taken 60 times.
✗ Branch 48 not taken.
✓ Branch 51 taken 60 times.
✗ Branch 52 not taken.
✓ Branch 55 taken 60 times.
✗ Branch 56 not taken.
✓ Branch 59 taken 60 times.
✗ Branch 60 not taken.
✓ Branch 63 taken 60 times.
✗ Branch 64 not taken.
✓ Branch 67 taken 60 times.
✗ Branch 68 not taken.
✓ Branch 71 taken 60 times.
✗ Branch 72 not taken.
✓ Branch 75 taken 60 times.
✗ Branch 76 not taken.
✓ Branch 79 taken 60 times.
✗ Branch 80 not taken.
✓ Branch 83 taken 60 times.
✗ Branch 84 not taken.
✓ Branch 87 taken 60 times.
✗ Branch 88 not taken.
✓ Branch 91 taken 60 times.
✗ Branch 92 not taken.
✓ Branch 95 taken 60 times.
✗ Branch 96 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 4798761 squaredDistancePointSegment(const Point& p, const Point& a, const Point& b)
98 {
99 4798761 const auto ab = b - a;
100 4798761 const auto ap = p - a;
101 4798761 const auto t = ap*ab;
102
103
2/2
✓ Branch 0 taken 1679250 times.
✓ Branch 1 taken 721673 times.
 4798761 if (t <= 0.0)
104 return ap.two_norm2();
105
106 3355415 const auto lengthSq = ab.two_norm2();
107
2/2
✓ Branch 0 taken 854634 times.
✓ Branch 1 taken 824616 times.
 3355415 if (t >= lengthSq)
108 3418536 return (b - p).two_norm2();
109
110 1646147 auto proj = a;
111 1646147 proj.axpy(t/lengthSq, ab);
112 3292294 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 4114 squaredDistancePointSegment(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
122 {
123 static_assert(Geometry::mydimension == 1, "Geometry has to be a segment");
124 17938 const auto& a = geometry.corner(0);
125 11026 const auto& b = geometry.corner(1);
126 11026 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 distancePointSegment(const Point& p, const Point& a, const Point& b)
136 6080 { 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 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 845393 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 845393 const auto ab = b - a;
158 845393 const auto bc = c - b;
159 845393 const auto ca = a - c;
160 845393 const auto normal = crossProduct(ab, ca);
161
162 845393 const auto ap = p - a;
163 845393 const auto bp = p - b;
164 845393 const auto cp = p - c;
165
166 1690786 const auto sum = sign(crossProduct(ab, normal)*ap)
167 1690786 + sign(crossProduct(bc, normal)*bp)
168
2/2
✓ Branch 1 taken 46435 times.
✓ Branch 2 taken 798958 times.
 1690786 + 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 46435 times.
✓ Branch 1 taken 798958 times.
 845393 if (sum < 2.0)
174 {
175 using std::min;
176 2396874 return min({squaredDistancePointSegment(p, a, b),
177 798958 squaredDistancePointSegment(p, a, c),
178 798958 squaredDistancePointSegment(p, b, c)});
179 }
180 // compute distance via orthogonal projection
181 else
182 {
183 46435 const auto tmp = normal*ap;
184 92870 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 428134 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
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 214067 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 214067 times.
 856268 assert(geometry.corners() == 3);
199 428134 const auto& a = geometry.corner(0);
200 428134 const auto& b = geometry.corner(1);
201 428134 const auto& c = geometry.corner(2);
202 428134 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 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 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 1058485 squaredDistancePointPolygon(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
231 {
232
4/4
✓ Branch 0 taken 213827 times.
✓ Branch 1 taken 315408 times.
✓ Branch 2 taken 213827 times.
✓ Branch 3 taken 315408 times.
 2116955 if (geometry.corners() == 3)
233 427654 return squaredDistancePointTriangle(p, geometry);
234
2/4
✓ Branch 0 taken 315408 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 315408 times.
✗ Branch 3 not taken.
 1261647 else if (geometry.corners() == 4)
235 {
236 630831 const auto& a = geometry.corner(0);
237 630831 const auto& b = geometry.corner(1);
238 630831 const auto& c = geometry.corner(2);
239 630831 const auto& d = geometry.corner(3);
240
241 using std::min;
242
2/2
✓ Branch 0 taken 115577 times.
✓ Branch 1 taken 199846 times.
 630831 return min(squaredDistancePointTriangle(p, a, b, d),
243 1261647 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 distancePointPolygon(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
257 15 { 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 760 const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
272
4/4
✓ Branch 0 taken 380 times.
✓ Branch 1 taken 3040 times.
✓ Branch 2 taken 380 times.
✓ Branch 3 taken 3040 times.
 7220 for (const auto& qp : quad)
273 12160 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 480 static inline ctype squaredDistance(const Dune::FieldVector<ctype, dimWorld>& a,
292 const Dune::FieldVector<ctype, dimWorld>& b)
293 480 { 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 720 static auto distance(const Geo1& geo1, const Geo2& geo2)
318 1440 { 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 720 static auto distance(const Geo1& geo1, const Geo2& geo2)
327 1800 { 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 480 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 480 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 static inline auto squaredDistance(const Geo1& geo1, const Geo2& geo2)
365 4560 { 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 3360 static inline auto distance(const Geo1& geo1, const Geo2& geo2)
373
28/56
✓ 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 17 taken 60 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 60 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 60 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 60 times.
✗ Branch 27 not taken.
✓ Branch 29 taken 60 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 60 times.
✗ Branch 33 not taken.
✓ Branch 35 taken 60 times.
✗ Branch 36 not taken.
✓ Branch 38 taken 60 times.
✗ Branch 39 not taken.
✓ Branch 42 taken 60 times.
✗ Branch 43 not taken.
✓ Branch 46 taken 60 times.
✗ Branch 47 not taken.
✓ Branch 50 taken 60 times.
✗ Branch 51 not taken.
✓ Branch 54 taken 60 times.
✗ Branch 55 not taken.
✓ Branch 58 taken 60 times.
✗ Branch 59 not taken.
✓ Branch 62 taken 60 times.
✗ Branch 63 not taken.
✓ Branch 65 taken 60 times.
✗ Branch 66 not taken.
✓ Branch 68 taken 60 times.
✗ Branch 69 not taken.
✓ Branch 71 taken 60 times.
✗ Branch 72 not taken.
✓ Branch 74 taken 60 times.
✗ Branch 75 not taken.
✓ Branch 77 taken 60 times.
✗ Branch 78 not taken.
✓ Branch 80 taken 60 times.
✗ Branch 81 not taken.
✓ Branch 83 taken 60 times.
✗ Branch 84 not taken.
✓ Branch 86 taken 60 times.
✗ Branch 87 not taken.
✓ Branch 90 taken 60 times.
✗ Branch 91 not taken.
✓ Branch 94 taken 60 times.
✗ Branch 95 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 7215428 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 2337735 times.
✓ Branch 1 taken 1300243 times.
 7215428 const auto& bBox = tree.getBoundingBoxNode(node);
394
395 // If bounding box is outside radius, then don't search further
396
2/2
✓ Branch 0 taken 2337735 times.
✓ Branch 1 taken 1300243 times.
 7215428 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 2337735 times.
✓ Branch 1 taken 1300243 times.
 7215428 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 535790 times.
✓ Branch 1 taken 1801945 times.
 4642038 else if (tree.isLeaf(bBox, node))
405 {
406 1065244 const std::size_t entityIdx = bBox.child1;
407 535790 const auto squaredDistance = [&]{
408 535790 const auto geometry = tree.entitySet().entity(entityIdx).geometry();
409 if constexpr (EntitySet::Entity::Geometry::mydimension == 2)
410 528597 return squaredDistancePointPolygon(point, geometry);
411 else if constexpr (EntitySet::Entity::Geometry::mydimension == 1)
412 7515 return squaredDistancePointSegment(point, geometry);
413 else
414 DUNE_THROW(Dune::NotImplemented, "squaredDistance to entity with dim>2");
415 1601034 }();
416
417
2/2
✓ Branch 0 taken 49432 times.
✓ Branch 1 taken 486358 times.
 1065244 if (squaredDistance < minSquaredDistance)
418 {
419 92528 eIdx = entityIdx;
420 92528 minSquaredDistance = squaredDistance;
421 }
422 }
423
424 // Check the children nodes recursively
425 else
426 {
427 3576794 closestEntity(point, tree, bBox.child0, minSquaredDistance, eIdx);
428 3576794 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 10517322 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 3105778 times.
✓ Branch 1 taken 7459956 times.
 20972804 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 3105778 times.
✓ Branch 1 taken 7459956 times.
 20972804 if (tree.isLeaf(bBox, node))
450 {
451 6138906 const std::size_t entityIdx = bBox.child1;
452 12277812 const auto& p = tree.entitySet().entity(entityIdx).geometry().corner(0);
453 6220856 const auto squaredDistance = (p-point).two_norm2();
454
455
2/2
✓ Branch 0 taken 1408555 times.
✓ Branch 1 taken 1697223 times.
 6138906 if (squaredDistance < minSquaredDistance)
456 {
457 2760040 eIdx = entityIdx;
458 2760040 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
2/2
✓ Branch 0 taken 5265823 times.
✓ Branch 1 taken 2194133 times.
 14833898 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 5265823 times.
✓ Branch 1 taken 2194133 times.
 14833898 if (squaredDistance > minSquaredDistance) return;
472
473 10455482 closestEntity(point, tree, bBox.child0, minSquaredDistance, eIdx);
474 10455482 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 136352 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 136352 std::size_t eIdx = 0;
499 272704 Detail::closestEntity(point, tree, tree.numBoundingBoxes() - 1, minSquaredDistance, eIdx);
500 using std::sqrt;
501 272704 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 ctype squaredDistance(const Dune::FieldVector<ctype, dimworld>& point,
510 const BoundingBoxTree<EntitySet>& tree,
511 ctype minSquaredDistance = std::numeric_limits<ctype>::max())
512 {
513 34088 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