3.3.0
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
geometry/intersectspointsimplex.hh
Go to the documentation of this file.
1/*****************************************************************************
2 * See the file COPYING for full copying permissions. *
3 * *
4 * This program is free software: you can redistribute it and/or modify *
5 * it under the terms of the GNU General Public License as published by *
6 * the Free Software Foundation, either version 3 of the License, or *
7 * (at your option) any later version. *
8 * *
9 * This program is distributed in the hope that it will be useful, *
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
12 * GNU General Public License for more details. *
13 * *
14 * You should have received a copy of the GNU General Public License *
15 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
16 *****************************************************************************/
22#ifndef DUMUX_GEOMETRY_INTERSECTS_POINT_SIMPLEX_HH
23#define DUMUX_GEOMETRY_INTERSECTS_POINT_SIMPLEX_HH
24
25#include <cmath>
26#include <dune/common/fvector.hh>
27#include <dumux/common/math.hh>
28
29namespace Dumux {
30
35template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3), int> = 0>
36bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
37 const Dune::FieldVector<ctype, dimworld>& p0,
38 const Dune::FieldVector<ctype, dimworld>& p1,
39 const Dune::FieldVector<ctype, dimworld>& p2,
40 const Dune::FieldVector<ctype, dimworld>& p3)
41{
42 // Algorithm from http://www.blackpawn.com/texts/pointinpoly/
43 // See also "Real-Time Collision Detection" by Christer Ericson.
44 using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
45 static constexpr ctype eps_ = 1.0e-7;
46
47 // put the tetrahedron points in an array
48 const GlobalPosition *p[4] = {&p0, &p1, &p2, &p3};
49
50 // iterate over all faces
51 for (int i = 0; i < 4; ++i)
52 {
53 // compute all the vectors from vertex (local index 0) to the other points
54 const GlobalPosition v1 = *p[(i + 1)%4] - *p[i];
55 const GlobalPosition v2 = *p[(i + 2)%4] - *p[i];
56 const GlobalPosition v3 = *p[(i + 3)%4] - *p[i];
57 const GlobalPosition v = point - *p[i];
58 // compute the normal to the facet (cross product)
59 GlobalPosition n1 = crossProduct(v1, v2);
60 n1 /= n1.two_norm();
61 // find out on which side of the plane v and v3 are
62 const auto t1 = n1.dot(v);
63 const auto t2 = n1.dot(v3);
64 // If the point is not exactly on the plane the
65 // points have to be on the same side
66 const auto eps = eps_ * v1.two_norm();
67 if ((t1 > eps || t1 < -eps) && std::signbit(t1) != std::signbit(t2))
68 return false;
69 }
70 return true;
71}
72
77template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3), int> = 0>
78bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
79 const Dune::FieldVector<ctype, dimworld>& p0,
80 const Dune::FieldVector<ctype, dimworld>& p1,
81 const Dune::FieldVector<ctype, dimworld>& p2)
82{
83 // adapted from the algorithm from from "Real-Time Collision Detection" by Christer Ericson,
84 // published by Morgan Kaufmann Publishers, (c) 2005 Elsevier Inc. (Chapter 5.4.2)
85 constexpr ctype eps_ = 1.0e-7;
86
87 // compute the normal of the triangle
88 const auto v1 = p0 - p2;
89 auto n = crossProduct(v1, p1 - p0);
90 const ctype nnorm = n.two_norm();
91 const ctype eps4 = eps_*nnorm*nnorm; // compute an epsilon for later
92 n /= nnorm; // normalize
93
94 // first check if we are in the plane of the triangle
95 // if not we can return early
96 using std::abs;
97 auto x = p0 - point;
98 x /= x.two_norm(); // normalize
99
100 if (abs(x*n) > eps_)
101 return false;
102
103 // translate the triangle so that 'point' is the origin
104 const auto a = p0 - point;
105 const auto b = p1 - point;
106 const auto c = p2 - point;
107
108 // compute the normal vectors for triangles P->A->B and P->B->C
109 const auto u = crossProduct(b, c);
110 const auto v = crossProduct(c, a);
111
112 // they have to point in the same direction or be orthogonal
113 if (u*v < 0.0 - eps4)
114 return false;
115
116 // compute the normal vector for triangle P->C->A
117 const auto w = crossProduct(a, b);
118
119 // it also has to point in the same direction or be orthogonal
120 if (u*w < 0.0 - eps4)
121 return false;
122
123 // check if point is on the line of one of the edges
124 if (u.two_norm2() < eps4)
125 return b*c < 0.0 + eps_*nnorm;
126 if (v.two_norm2() < eps4)
127 return a*c < 0.0 + eps_*nnorm;
128
129 // now the point must be in the triangle (or on the faces)
130 return true;
131}
132
137template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 2), int> = 0>
138bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
139 const Dune::FieldVector<ctype, dimworld>& p0,
140 const Dune::FieldVector<ctype, dimworld>& p1,
141 const Dune::FieldVector<ctype, dimworld>& p2)
142{
143 static constexpr ctype eps_ = 1.0e-7;
144
145 // Use barycentric coordinates
146 const ctype A = 0.5*(-p1[1]*p2[0] + p0[1]*(p2[0] - p1[0])
147 +p1[0]*p2[1] + p0[0]*(p1[1] - p2[1]));
148 const ctype sign = std::copysign(1.0, A);
149 const ctype s = sign*(p0[1]*p2[0] + point[0]*(p2[1]-p0[1])
150 -p0[0]*p2[1] + point[1]*(p0[0]-p2[0]));
151 const ctype t = sign*(p0[0]*p1[1] + point[0]*(p0[1]-p1[1])
152 -p0[1]*p1[0] + point[1]*(p1[0]-p0[0]));
153 const ctype eps = sign*A*eps_;
154
155 return (s > -eps
156 && t > -eps
157 && (s + t) < 2*A*sign + eps);
158}
159
165template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 3 || dimworld == 2), int> = 0>
166bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
167 const Dune::FieldVector<ctype, dimworld>& p0,
168 const Dune::FieldVector<ctype, dimworld>& p1)
169{
170 using GlobalPosition = Dune::FieldVector<ctype, dimworld>;
171 static constexpr ctype eps_ = 1.0e-7;
172
173 // compute the vectors between p0 and the other points
174 const GlobalPosition v1 = p1 - p0;
175 const GlobalPosition v2 = point - p0;
176
177 const ctype v1norm = v1.two_norm();
178 const ctype v2norm = v2.two_norm();
179
180 // check if point and p0 are the same
181 if (v2norm < v1norm*eps_)
182 return true;
183
184 return (v1.dot(v2) > v1norm*v2norm*(1.0 - eps_) && v2norm < v1norm*(1.0 + eps_));
185}
186
191template<class ctype, int dimworld, typename std::enable_if_t<(dimworld == 1), int> = 0>
192bool intersectsPointSimplex(const Dune::FieldVector<ctype, dimworld>& point,
193 const Dune::FieldVector<ctype, dimworld>& p0,
194 const Dune::FieldVector<ctype, dimworld>& p1)
195{
196 static constexpr ctype eps_ = 1.0e-7;
197
198 // sort the interval so interval[1] is the end and interval[0] the start
199 const ctype *interval[2] = {&p0[0], &p1[0]};
200 if (*interval[0] > *interval[1])
201 std::swap(interval[0], interval[1]);
202
203 const ctype v1 = point[0] - *interval[0];
204 const ctype v2 = *interval[1] - *interval[0]; // always positive
205
206 // the coordinates are the same
207 using std::abs;
208 if (abs(v1) < v2*eps_)
209 return true;
210
211 // the point doesn't coincide with p0
212 // so if p0 and p1 are equal it's not inside
213 if (v2 < 1.0e-30)
214 return false;
215
216 // the point is inside if the length is
217 // smaller than the interval length and the
218 // sign of v1 & v2 are the same
219 using std::signbit;
220 return (!signbit(v1) && abs(v1) < v2*(1.0 + eps_));
221}
222
223} // end namespace Dumux
224
225#endif
Define some often used mathematical functions.
bool intersectsPointSimplex(const Dune::FieldVector< ctype, dimworld > &point, const Dune::FieldVector< ctype, dimworld > &p0, const Dune::FieldVector< ctype, dimworld > &p1, const Dune::FieldVector< ctype, dimworld > &p2, const Dune::FieldVector< ctype, dimworld > &p3)
Find out whether a point is inside a tetrahedron (p0, p1, p2, p3) (dimworld is 3)
Definition: geometry/intersectspointsimplex.hh:36
Dune::FieldVector< Scalar, 3 > crossProduct(const Dune::FieldVector< Scalar, 3 > &vec1, const Dune::FieldVector< Scalar, 3 > &vec2)
Cross product of two vectors in three-dimensional Euclidean space.
Definition: math.hh:640
constexpr int sign(const ValueType &value) noexcept
Sign or signum function.
Definition: math.hh:627
Definition: adapt.hh:29