GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/geometry/intersectspointsimplex.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 65 66 98.5%
Functions: 6 6 100.0%
Branches: 108 116 93.1%
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 Detect if a point intersects a simplex (including boundary)
11 */
12 #ifndef DUMUX_GEOMETRY_INTERSECTS_POINT_SIMPLEX_HH
13 #define DUMUX_GEOMETRY_INTERSECTS_POINT_SIMPLEX_HH
14
15 #include <cmath>
16 #include <dune/common/fvector.hh>
17 #include <dumux/common/math.hh>
18
19 namespace Dumux {
20
21 /*!
22 * \ingroup Geometry
23 * \brief Find out whether a point is inside the tetrahedron (p0, p1, p2, p3) (dimworld is 3)
24 */
25 template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3), int> = 0>
26 1913269 bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
27 const Dune::FieldVector<ctype, dimworld>& p0,
28 const Dune::FieldVector<ctype, dimworld>& p1,
29 const Dune::FieldVector<ctype, dimworld>& p2,
30 const Dune::FieldVector<ctype, dimworld>& p3)
31 {
32 // Algorithm from http://www.blackpawn.com/texts/pointinpoly/
33 // See also "Real-Time Collision Detection" by Christer Ericson.
34 using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
35 static constexpr ctype eps_ = 1.0e-7;
36
37 // put the tetrahedron points in an array
38 1913269 const GlobalPosition *p[4] = {&p0, &p1, &p2, &p3};
39
40 // iterate over all faces
41
2/2
✓ Branch 0 taken 4697621 times.
✓ Branch 1 taken 329312 times.
 5026933 for (int i = 0; i < 4; ++i)
42 {
43 // compute all the vectors from vertex (local index 0) to the other points
44 4697621 const GlobalPosition v1 = *p[(i + 1)%4] - *p[i];
45 4697621 const GlobalPosition v2 = *p[(i + 2)%4] - *p[i];
46 4697621 const GlobalPosition v3 = *p[(i + 3)%4] - *p[i];
47 4697621 const GlobalPosition v = point - *p[i];
48 // compute the normal to the facet (cross product)
49 4697621 GlobalPosition n1 = crossProduct(v1, v2);
50 9395242 n1 /= n1.two_norm();
51 // find out on which side of the plane v and v3 are
52 const auto t1 = n1.dot(v);
53 const auto t2 = n1.dot(v3);
54 // If the point is not exactly on the plane the
55 // points have to be on the same side
56 4697621 const auto eps = eps_ * v1.two_norm();
57
10/10
✓ Branch 0 taken 2588287 times.
✓ Branch 1 taken 2109334 times.
✓ Branch 2 taken 2311842 times.
✓ Branch 3 taken 276445 times.
✓ Branch 4 taken 1583957 times.
✓ Branch 5 taken 2837219 times.
✓ Branch 6 taken 1583957 times.
✓ Branch 7 taken 2837219 times.
✓ Branch 8 taken 1583957 times.
✓ Branch 9 taken 2837219 times.
 4697621 if ((t1 > eps || t1 < -eps) && std::signbit(t1) != std::signbit(t2))
58 1583957 return false;
59 }
60 return true;
61 }
62
63 /*!
64 * \ingroup Geometry
65 * \brief Find out whether a point is inside the triangle (p0, p1, p2) (dimworld is 3)
66 */
67 template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3), int> = 0>
68 411165 bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
69 const Dune::FieldVector<ctype, dimworld>& p0,
70 const Dune::FieldVector<ctype, dimworld>& p1,
71 const Dune::FieldVector<ctype, dimworld>& p2)
72 {
73 // adapted from the algorithm from from "Real-Time Collision Detection" by Christer Ericson,
74 // published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. (Chapter 5.4.2)
75 411165 constexpr ctype eps_ = 1.0e-7;
76
77 // compute the normal of the triangle
78 411165 const auto v1 = p0 - p2;
79 822330 auto n = crossProduct(v1, p1 - p0);
80 411165 const ctype nnorm = n.two_norm();
81 411165 const ctype eps4 = eps_*nnorm*nnorm; // compute an epsilon for later
82 411165 n /= nnorm; // normalize
83
84 // first check if we are in the plane of the triangle
85 // if not we can return early
86 using std::abs;
87 411165 auto x = p0 - point;
88 411165 x /= x.two_norm(); // normalize
89
90
4/4
✓ Branch 0 taken 154692 times.
✓ Branch 1 taken 256473 times.
✓ Branch 2 taken 154692 times.
✓ Branch 3 taken 256473 times.
 822330 if (abs(x*n) > eps_)
91 return false;
92
93 // translate the triangle so that 'point' is the origin
94 154692 const auto a = p0 - point;
95 154692 const auto b = p1 - point;
96 154692 const auto c = p2 - point;
97
98 // compute the normal vectors for triangles P->A->B and P->B->C
99 154692 const auto u = crossProduct(b, c);
100 154692 const auto v = crossProduct(c, a);
101
102 // they have to point in the same direction or be orthogonal
103
2/2
✓ Branch 1 taken 102121 times.
✓ Branch 2 taken 52571 times.
 309384 if (u*v < 0.0 - eps4)
104 return false;
105
106 // compute the normal vector for triangle P->C->A
107 102121 const auto w = crossProduct(a, b);
108
109 // it also has to point in the same direction or be orthogonal
110
2/2
✓ Branch 1 taken 80354 times.
✓ Branch 2 taken 21767 times.
 204242 if (u*w < 0.0 - eps4)
111 return false;
112
113 // check if point is on the line of one of the edges
114
2/2
✓ Branch 0 taken 6583 times.
✓ Branch 1 taken 73771 times.
 80354 if (u.two_norm2() < eps4)
115 6583 return b*c < 0.0 + eps_*nnorm;
116
2/2
✓ Branch 0 taken 4357 times.
✓ Branch 1 taken 69414 times.
 73771 if (v.two_norm2() < eps4)
117 4357 return a*c < 0.0 + eps_*nnorm;
118
119 // now the point must be in the triangle (or on the faces)
120 return true;
121 }
122
123 /*!
124 * \ingroup Geometry
125 * \brief Find out whether a point is inside the triangle (p0, p1, p2) (dimworld is 2)
126 */
127 template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 2), int> = 0>
128 2722601 bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
129 const Dune::FieldVector<ctype, dimworld>& p0,
130 const Dune::FieldVector<ctype, dimworld>& p1,
131 const Dune::FieldVector<ctype, dimworld>& p2)
132 {
133 static constexpr ctype eps_ = 1.0e-7;
134
135 // Use barycentric coordinates
136
10/10
✓ Branch 0 taken 1747069 times.
✓ Branch 1 taken 975532 times.
✓ Branch 2 taken 1747069 times.
✓ Branch 3 taken 975532 times.
✓ Branch 4 taken 1747069 times.
✓ Branch 5 taken 975532 times.
✓ Branch 6 taken 1747069 times.
✓ Branch 7 taken 975532 times.
✓ Branch 8 taken 1747069 times.
✓ Branch 9 taken 975532 times.
 13613005 const ctype A = 0.5*(-p1[1]*p2[0] + p0[1]*(p2[0] - p1[0])
137
10/10
✓ Branch 0 taken 1747069 times.
✓ Branch 1 taken 975532 times.
✓ Branch 2 taken 1747069 times.
✓ Branch 3 taken 975532 times.
✓ Branch 4 taken 1747069 times.
✓ Branch 5 taken 975532 times.
✓ Branch 6 taken 1747069 times.
✓ Branch 7 taken 975532 times.
✓ Branch 8 taken 1747069 times.
✓ Branch 9 taken 975532 times.
 13613005 +p1[0]*p2[1] + p0[0]*(p1[1] - p2[1]));
138 2722601 const ctype sign = std::copysign(1.0, A);
139
10/10
✓ Branch 0 taken 1747069 times.
✓ Branch 1 taken 975532 times.
✓ Branch 2 taken 1747069 times.
✓ Branch 3 taken 975532 times.
✓ Branch 4 taken 1747069 times.
✓ Branch 5 taken 975532 times.
✓ Branch 6 taken 1747069 times.
✓ Branch 7 taken 975532 times.
✓ Branch 8 taken 1747069 times.
✓ Branch 9 taken 975532 times.
 13613005 const ctype s = sign*(p0[1]*p2[0] + point[0]*(p2[1]-p0[1])
140
10/10
✓ Branch 0 taken 1747069 times.
✓ Branch 1 taken 975532 times.
✓ Branch 2 taken 1747069 times.
✓ Branch 3 taken 975532 times.
✓ Branch 4 taken 1747069 times.
✓ Branch 5 taken 975532 times.
✓ Branch 6 taken 1747069 times.
✓ Branch 7 taken 975532 times.
✓ Branch 8 taken 1747069 times.
✓ Branch 9 taken 975532 times.
 13613005 -p0[0]*p2[1] + point[1]*(p0[0]-p2[0]));
141
10/10
✓ Branch 0 taken 1747069 times.
✓ Branch 1 taken 975532 times.
✓ Branch 2 taken 1747069 times.
✓ Branch 3 taken 975532 times.
✓ Branch 4 taken 1747069 times.
✓ Branch 5 taken 975532 times.
✓ Branch 6 taken 1747069 times.
✓ Branch 7 taken 975532 times.
✓ Branch 8 taken 1747069 times.
✓ Branch 9 taken 975532 times.
 13613005 const ctype t = sign*(p0[0]*p1[1] + point[0]*(p0[1]-p1[1])
142
10/10
✓ Branch 0 taken 1747069 times.
✓ Branch 1 taken 975532 times.
✓ Branch 2 taken 1747069 times.
✓ Branch 3 taken 975532 times.
✓ Branch 4 taken 1747069 times.
✓ Branch 5 taken 975532 times.
✓ Branch 6 taken 1747069 times.
✓ Branch 7 taken 975532 times.
✓ Branch 8 taken 1747069 times.
✓ Branch 9 taken 975532 times.
 13613005 -p0[1]*p1[0] + point[1]*(p1[0]-p0[0]));
143 2722601 const ctype eps = sign*A*eps_;
144
145 2722601 return (s > -eps
146
2/2
✓ Branch 0 taken 1056605 times.
✓ Branch 1 taken 690464 times.
 1747069 && t > -eps
147
4/4
✓ Branch 0 taken 1747069 times.
✓ Branch 1 taken 975532 times.
✓ Branch 2 taken 406658 times.
✓ Branch 3 taken 649947 times.
 3779206 && (s + t) < 2*A*sign + eps);
148 }
149
150 /*!
151 * \ingroup Geometry
152 * \brief Find out whether a point is inside the interval (p0, p1) (dimworld is 2 or 3)
153 * \note We assume the given interval has non-zero length and use it to scale the epsilon
154 */
155 template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3 || dimworld == 2), int> = 0>
156 1676577 bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
157 const Dune::FieldVector<ctype, dimworld>& p0,
158 const Dune::FieldVector<ctype, dimworld>& p1)
159 {
160 using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
161 static constexpr ctype eps_ = 1.0e-7;
162
163 // compute the vectors between p0 and the other points
164 1676577 const GlobalPosition v1 = p1 - p0;
165 1676577 const GlobalPosition v2 = point - p0;
166
167 1676577 const ctype v1norm = v1.two_norm();
168 1676577 const ctype v2norm = v2.two_norm();
169
170 // early exit if point and p0 are the same
171
2/2
✓ Branch 0 taken 1650047 times.
✓ Branch 1 taken 26434 times.
 1676577 if (v2norm < v1norm*eps_)
172 return true;
173
174 // early exit if the point is outside the segment (no epsilon in the
175 // first statement because we already did the above equality check)
176
4/4
✓ Branch 0 taken 1279453 times.
✓ Branch 1 taken 370594 times.
✓ Branch 2 taken 438830 times.
✓ Branch 3 taken 840623 times.
 1650135 if (v1.dot(v2) < 0.0 || v2norm > v1norm*(1.0 + eps_))
177 return false;
178
179 // If the area spanned by the 2 vectors is zero, the points are colinear.
180 // If that is the case, the given point is on the segment.
181 438886 const auto n = crossProduct(v1, v2);
182 438886 const auto eps2 = v1norm*v1norm*eps_;
183 if constexpr (dimworld == 3)
184 688052 return n.two_norm2() < eps2*eps2;
185 else
186 {
187 using std::abs;
188 189720 return abs(n) < eps2;
189 }
190 }
191
192 /*!
193 * \ingroup Geometry
194 * \brief Find out whether a point is inside the interval (p0, p1) (dimworld is 1)
195 */
196 template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 1), int> = 0>
197 20 bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
198 const Dune::FieldVector<ctype, dimworld>& p0,
199 const Dune::FieldVector<ctype, dimworld>& p1)
200 {
201 static constexpr ctype eps_ = 1.0e-7;
202
203 // sort the interval so interval[1] is the end and interval[0] the start
204
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 20 times.
 40 const ctype *interval[2] = {&p0[0], &p1[0]};
205
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 20 times.
 20 if (*interval[0] > *interval[1])
206 ✗ std::swap(interval[0], interval[1]);
207
208 20 const ctype v1 = point[0] - *interval[0];
209 20 const ctype v2 = *interval[1] - *interval[0]; // always positive
210
211 // the coordinates are the same
212 using std::abs;
213
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 8 times.
 40 if (abs(v1) < v2*eps_)
214 return true;
215
216 // the point doesn't coincide with p0
217 // so if p0 and p1 are equal it's not inside
218
1/2
✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
 12 if (v2 < 1.0e-30)
219 return false;
220
221 // the point is inside if the length is
222 // smaller than the interval length and the
223 // sign of v1 & v2 are the same
224 using std::signbit;
225
4/8
✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 12 times.
 24 return (!signbit(v1) && abs(v1) < v2*(1.0 + eps_));
226 }
227
228 } // end namespace Dumux
229
230 #endif
231