Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "LineSegment.h"
11 :
12 : #include "JsonIO.h"
13 :
14 : #include "libmesh/plane.h"
15 : #include "libmesh/vector_value.h"
16 :
17 : using libMesh::Point;
18 : using libMesh::Real;
19 : using libMesh::RealVectorValue;
20 :
21 1067 : LineSegment::LineSegment(const Point & p0, const Point & p1) : _p0(p0), _p1(p1) {}
22 :
23 : bool
24 793 : LineSegment::closest_point(const Point & p, bool clamp_to_segment, Point & closest_p) const
25 : {
26 793 : Point p0_p = p - _p0;
27 793 : Point p0_p1 = _p1 - _p0;
28 793 : Real p0_p1_2 = p0_p1.norm_sq();
29 793 : Real perp = p0_p(0) * p0_p1(0) + p0_p(1) * p0_p1(1) + p0_p(2) * p0_p1(2);
30 793 : Real t = perp / p0_p1_2;
31 793 : bool on_segment = true;
32 :
33 793 : if (t < 0.0 || t > 1.0)
34 268 : on_segment = false;
35 :
36 793 : if (clamp_to_segment)
37 : {
38 218 : if (t < 0.0)
39 32 : t = 0.0;
40 186 : else if (t > 1.0)
41 66 : t = 1.0;
42 : }
43 :
44 793 : closest_p = _p0 + p0_p1 * t;
45 793 : return on_segment;
46 : }
47 :
48 : Point
49 218 : LineSegment::closest_point(const Point & p) const
50 : {
51 218 : Point closest_p;
52 218 : closest_point(p, true, closest_p);
53 218 : return closest_p;
54 : }
55 :
56 : bool
57 66 : LineSegment::closest_normal_point(const Point & p, Point & closest_p) const
58 : {
59 66 : return closest_point(p, false, closest_p);
60 : }
61 :
62 : bool
63 509 : LineSegment::contains_point(const Point & p) const
64 : {
65 509 : Point closest_p;
66 1018 : return closest_point(p, false, closest_p) && closest_p.absolute_fuzzy_equals(p);
67 : }
68 :
69 : bool
70 624 : LineSegment::intersect(const libMesh::Plane & pl, Point & intersect_p) const
71 : {
72 : /**
73 : * There are three cases in 3D for intersection of a line and a plane
74 : * Case 1: The line is parallel to the plane - No intersection
75 : * Numerator = non-zero
76 : * Denominator = zero
77 : *
78 : * Case 2: The line is within the plane - Inf intersection
79 : * Numerator = zero
80 : * Denominator = zero
81 : *
82 : * Case 3: The line intersects the plane at a single point
83 : * Denominator = non-zero
84 : */
85 :
86 624 : Point pl0 = pl.get_planar_point();
87 624 : RealVectorValue N = pl.unit_normal(_p0);
88 624 : RealVectorValue I = (_p1 - _p0).unit();
89 :
90 624 : Real numerator = (pl0 - _p0) * N;
91 624 : Real denominator = I * N;
92 :
93 : // The Line is parallel to the plane
94 624 : if (std::abs(denominator) < 1.e-10)
95 : {
96 : // The Line is on the plane
97 462 : if (std::abs(numerator) < 1.e-10)
98 : {
99 : // The solution is not unique so we'll just pick an end point for now
100 10 : intersect_p = _p0;
101 10 : return true;
102 : }
103 452 : return false;
104 : }
105 :
106 162 : Real d = numerator / denominator;
107 :
108 : // Make sure we haven't moved off the line segment!
109 162 : if (d + libMesh::TOLERANCE < 0 || d - libMesh::TOLERANCE > (_p1 - _p0).norm())
110 52 : return false;
111 :
112 110 : intersect_p = d * I + _p0;
113 :
114 110 : return true;
115 : }
116 :
117 : bool
118 439 : LineSegment::intersect(const LineSegment & l, Point & intersect_p) const
119 : {
120 : /**
121 : * First check for concurance:
122 : *
123 : *
124 : * | x1 y1 z1 1 |
125 : * | x2 y2 z2 1 | = (x3 - x1) * [(x2-x1) x (x4-x3)] = 0
126 : * | x3 y3 z3 1 |
127 : * | x4 y4 z4 1 |
128 : *
129 : *
130 : * Solve:
131 : * x = _p0 + (_p1 - _p0)*s
132 : * x = l.p0 + (l._p1 - l.p0)*t
133 : *
134 : * where
135 : * a = _p1 - _p0
136 : * b = l._p1 - l._p0
137 : * c = l._p0 - _p0
138 : *
139 : * s = (c x b) * (a x b) / | a x b |^2
140 : */
141 439 : RealVectorValue a = _p1 - _p0;
142 439 : RealVectorValue b = l._p1 - l._p0;
143 439 : RealVectorValue c = l._p0 - _p0;
144 :
145 439 : RealVectorValue v = a.cross(b);
146 :
147 439 : const auto tol = 1.e-10;
148 :
149 : // Parallel lines check
150 439 : if (v.norm() < tol)
151 : {
152 : // Parallel but not collinear: no intersection
153 215 : if (c.cross(a).norm() >= tol)
154 185 : return false;
155 :
156 : // Collinear: segments overlap if any endpoint lies on the other segment
157 42 : const bool overlap = this->contains_point(l._p0) || this->contains_point(l._p1) ||
158 42 : l.contains_point(_p0) || l.contains_point(_p1);
159 :
160 30 : if (!overlap)
161 6 : return false;
162 :
163 : // Pick a deterministic intersection point on the overlap
164 24 : if (this->contains_point(l._p0))
165 18 : intersect_p = l._p0;
166 6 : else if (this->contains_point(l._p1))
167 4 : intersect_p = l._p1;
168 2 : else if (l.contains_point(_p0))
169 2 : intersect_p = _p0;
170 : else
171 0 : intersect_p = _p1;
172 :
173 24 : return true;
174 : }
175 :
176 : // Check that the lines are coplanar
177 224 : Real concur = c * (a.cross(b));
178 224 : if (std::abs(concur) > 1.e-10)
179 16 : return false;
180 :
181 208 : Real s = (c.cross(b) * v) / (v * v);
182 208 : Real t = (c.cross(a) * v) / (v * v);
183 :
184 : // if s and t are between 0 and 1 then the Line Segments intersect
185 : // TODO: We could handle other case of clamping to the end of Line
186 : // Segements if we want to here
187 :
188 208 : if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
189 : {
190 147 : intersect_p = _p0 + s * a;
191 147 : return true;
192 : }
193 61 : return false;
194 :
195 : /**
196 : * Parameteric Equation of lines
197 : * _p0 + t(v0) = l._p0 + u(v1)
198 : *
199 : * Case 1: Parallel Lines
200 : * v0 x v1 == 0
201 : *
202 : * Case 1a: Collinear Lines
203 : * v0 x v1 == 0
204 : * (l._p0 - _p0) x (_p1 - _p0) == 0
205 : *
206 : * Case 2: Intersecting Lines
207 : * 0 <= t <= 1
208 : * 0 <= u <= 1
209 : *
210 : *
211 : * Case 1: The lines do not intersect
212 : * vleft cross vright = non-zero
213 : *
214 : * Case 2: The lines are co-linear
215 : * vleft cross vright = zero
216 : * vleft (Denominator) = zero
217 : *
218 : * Case 3: The line intersect at a single point
219 : * vleft cross vright = zero
220 : * vleft (Denominator) = non-zero
221 : RealVectorValue v0 = _p1 - _p0;
222 : RealVectorValue v1 = l._p1 - l._p0;
223 : RealVectorValue v2 = l._p0 - _p0;
224 :
225 : RealVectorValue vbot = v0.cross(v1);
226 : RealVectorValue vtop = v2.cross(v1);
227 :
228 : RealVectorValue crossed = vleft.cross(vright);
229 :
230 : // Case 1: No intersection
231 : if (std::abs(vleft.cross(vright).size()) > 1.e-10)
232 : return false;
233 :
234 : // Case 2: Co-linear (just return one of the end points)
235 : if (std::abs(vleft.size()) < 1.e-10)
236 : {
237 : intersect_p = _p0;
238 : return true;
239 : }
240 :
241 : // Case 3:
242 :
243 : //TODO: We could detect whether the Line Segments actually overlap
244 : // instead of whether the Lines are co-linear
245 :
246 : Real a = vright.size()/vleft.size();
247 : intersect_p = _p0 + a*v0;
248 : return true;
249 : */
250 : }
251 :
252 : bool
253 0 : LineSegment::intersect(const LineSegment & line_segment) const
254 : {
255 0 : Point p;
256 0 : return intersect(line_segment, p);
257 : }
258 :
259 : void
260 0 : LineSegment::set(const Point & p0, const Point & p1)
261 : {
262 0 : setStart(p0);
263 0 : setEnd(p1);
264 0 : }
265 :
266 : void
267 0 : dataStore(std::ostream & stream, LineSegment & l, void * context)
268 : {
269 0 : dataStore(stream, l.start(), context);
270 0 : dataStore(stream, l.end(), context);
271 0 : }
272 :
273 : void
274 0 : dataLoad(std::istream & stream, LineSegment & l, void * context)
275 : {
276 0 : Point p0;
277 0 : dataLoad(stream, p0, context);
278 0 : Point p1;
279 0 : dataLoad(stream, p1, context);
280 0 : l.set(p0, p1);
281 0 : }
282 :
283 : void
284 0 : to_json(nlohmann::json & json, const LineSegment & l)
285 : {
286 0 : to_json(json["start"], l.start());
287 0 : to_json(json["end"], l.end());
288 0 : }
289 :
290 : Point
291 0 : LineSegment::normal() const
292 : {
293 0 : const auto tangent = _p1 - _p0;
294 :
295 : // Rotate 90 degrees counter-clockwise (2D)
296 0 : Point n(-tangent(1), tangent(0), 0.0);
297 0 : n /= n.norm();
298 :
299 0 : return n;
300 : }
301 :
302 : Ball
303 0 : LineSegment::computeBoundingBall() const
304 : {
305 0 : const Point center = 0.5 * (_p0 + _p1);
306 0 : const Real radius = 0.5 * (_p1 - _p0).norm();
307 0 : return Ball(center, radius);
308 : }
|
