| 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 | 411693 | 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 | 411693 | constexpr ctype eps_ = 1.0e-7; | |
| 76 | |||
| 77 | // compute the normal of the triangle | ||
| 78 | 411693 | const auto v1 = p0 - p2; | |
| 79 | 823386 | auto n = crossProduct(v1, p1 - p0); | |
| 80 | 411693 | const ctype nnorm = n.two_norm(); | |
| 81 | 411693 | const ctype eps4 = eps_*nnorm*nnorm; // compute an epsilon for later | |
| 82 | 411693 | 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 | 411693 | auto x = p0 - point; | |
| 88 | 411693 | x /= x.two_norm(); // normalize | |
| 89 | |||
| 90 |
4/4✓ Branch 0 taken 155220 times.
✓ Branch 1 taken 256473 times.
✓ Branch 2 taken 155220 times.
✓ Branch 3 taken 256473 times.
|
823386 | if (abs(x*n) > eps_) |
| 91 | return false; | ||
| 92 | |||
| 93 | // translate the triangle so that 'point' is the origin | ||
| 94 | 155220 | const auto a = p0 - point; | |
| 95 | 155220 | const auto b = p1 - point; | |
| 96 | 155220 | const auto c = p2 - point; | |
| 97 | |||
| 98 | // compute the normal vectors for triangles P->A->B and P->B->C | ||
| 99 | 155220 | const auto u = crossProduct(b, c); | |
| 100 | 155220 | 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 53099 times.
|
310440 | 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 |