FreeFOAM The Cross-Platform CFD Toolkit
lineI.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include <OpenFOAM/IOstreams.H>
27 
28 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32 
33 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
34 
35 template<class Point, class PointRef>
36 inline line<Point, PointRef>::line(const Point& start, const Point& end)
37 :
38  a_(start),
39  b_(end)
40 {}
41 
42 
43 template<class Point, class PointRef>
45 {
46  // Read beginning of line point pair
47  is.readBegin("line");
48 
49  is >> a_ >> b_;
50 
51  // Read end of line point pair
52  is.readEnd("line");
53 
54  // Check state of Istream
55  is.check("line::line(Istream& is)");
56 }
57 
58 
59 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
60 
61 template<class Point, class PointRef>
62 inline PointRef line<Point, PointRef>::start() const
63 {
64  return a_;
65 }
66 
67 template<class Point, class PointRef>
68 inline PointRef line<Point, PointRef>::end() const
69 {
70  return b_;
71 }
72 
73 
74 template<class Point, class PointRef>
75 inline Point line<Point, PointRef>::centre() const
76 {
77  return 0.5*(a_ + b_);
78 }
79 
80 
81 template<class Point, class PointRef>
82 inline scalar line<Point, PointRef>::mag() const
83 {
84  return ::Foam::mag(vec());
85 }
86 
87 
88 template<class Point, class PointRef>
89 inline Point line<Point, PointRef>::vec() const
90 {
91  return b_ - a_;
92 }
93 
94 
95 template<class Point, class PointRef>
97 {
98  Point v = vec();
99 
100  Point w(p - a_);
101 
102  scalar c1 = v & w;
103 
104  if (c1 <= 0)
105  {
106  return PointHit<Point>(false, a_, Foam::mag(p - a_), true);
107  }
108 
109  scalar c2 = v & v;
110 
111  if (c2 <= c1)
112  {
113  return PointHit<Point>(false, b_, Foam::mag(p - b_), true);
114  }
115 
116  scalar b = c1/c2;
117 
118  Point pb(a_ + b*v);
119 
120  return PointHit<Point>(true, pb, Foam::mag(p - pb), false);
121 }
122 
123 
124 template<class Point, class PointRef>
126 (
128  Point& thisPt,
129  Point& edgePt
130 ) const
131 {
132  // From Mathworld Line-Line distance/(Gellert et al. 1989, p. 538).
133  Point a(end() - start());
134  Point b(edge.end() - edge.start());
135  Point c(edge.start() - start());
136 
137  Point crossab = a ^ b;
138  scalar magCrossSqr = magSqr(crossab);
139 
140  if (magCrossSqr > VSMALL)
141  {
142  scalar s = ((c ^ b) & crossab)/magCrossSqr;
143  scalar t = ((c ^ a) & crossab)/magCrossSqr;
144 
145  // Check for end points outside of range 0..1
146  if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
147  {
148  // Both inside range 0..1
149  thisPt = start() + a*s;
150  edgePt = edge.start() + b*t;
151  }
152  else
153  {
154  // Do brute force. Distance of everything to everything.
155  // Can quite possibly be improved!
156 
157  // From edge endpoints to *this
158  PointHit<Point> this0(nearestDist(edge.start()));
159  PointHit<Point> this1(nearestDist(edge.end()));
160  scalar thisDist = min(this0.distance(), this1.distance());
161 
162  // From *this to edge
163  PointHit<Point> edge0(edge.nearestDist(start()));
164  PointHit<Point> edge1(edge.nearestDist(end()));
165  scalar edgeDist = min(edge0.distance(), edge1.distance());
166 
167  if (thisDist < edgeDist)
168  {
169  if (this0.distance() < this1.distance())
170  {
171  thisPt = this0.rawPoint();
172  edgePt = edge.start();
173  }
174  else
175  {
176  thisPt = this1.rawPoint();
177  edgePt = edge.end();
178  }
179  }
180  else
181  {
182  if (edge0.distance() < edge1.distance())
183  {
184  thisPt = start();
185  edgePt = edge0.rawPoint();
186  }
187  else
188  {
189  thisPt = end();
190  edgePt = edge1.rawPoint();
191  }
192  }
193  }
194  }
195  else
196  {
197  // Parallel lines. Find overlap of both lines by projecting onto
198  // direction vector (now equal for both lines).
199 
200  scalar edge0 = edge.start() & a;
201  scalar edge1 = edge.end() & a;
202  bool edgeOrder = edge0 < edge1;
203 
204  scalar minEdge = (edgeOrder ? edge0 : edge1);
205  scalar maxEdge = (edgeOrder ? edge1 : edge0);
206  const Point& minEdgePt = (edgeOrder ? edge.start() : edge.end());
207  const Point& maxEdgePt = (edgeOrder ? edge.end() : edge.start());
208 
209  scalar this0 = start() & a;
210  scalar this1 = end() & a;
211  bool thisOrder = this0 < this1;
212 
213  scalar minThis = min(this0, this1);
214  scalar maxThis = max(this1, this0);
215  const Point& minThisPt = (thisOrder ? start() : end());
216  const Point& maxThisPt = (thisOrder ? end() : start());
217 
218  if (maxEdge < minThis)
219  {
220  // edge completely below *this
221  edgePt = maxEdgePt;
222  thisPt = minThisPt;
223  }
224  else if (maxEdge < maxThis)
225  {
226  // maxEdge inside interval of *this
227  edgePt = maxEdgePt;
228  thisPt = nearestDist(edgePt).rawPoint();
229  }
230  else
231  {
232  // maxEdge outside. Check if minEdge inside.
233  if (minEdge < minThis)
234  {
235  // Edge completely envelops this. Take any this point and
236  // determine nearest on edge.
237  thisPt = minThisPt;
238  edgePt = edge.nearestDist(thisPt).rawPoint();
239  }
240  else if (minEdge < maxThis)
241  {
242  // minEdge inside this interval.
243  edgePt = minEdgePt;
244  thisPt = nearestDist(edgePt).rawPoint();
245  }
246  else
247  {
248  // minEdge outside this interval
249  edgePt = minEdgePt;
250  thisPt = maxThisPt;
251  }
252  }
253  }
254 
255  return Foam::mag(thisPt - edgePt);
256 }
257 
258 
259 // * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * //
260 
261 template<class Point, class PointRef>
263 {
264  // Read beginning of line point pair
265  is.readBegin("line");
266 
267  is >> l.a_ >> l.b_;
268 
269  // Read end of line point pair
270  is.readEnd("line");
271 
272  // Check state of Ostream
273  is.check("Istream& operator>>(Istream&, line&)");
274 
275  return is;
276 }
277 
278 
279 template<class Point, class PointRef>
280 inline Ostream& operator<<(Ostream& os, const line<Point, PointRef>& l)
281 {
282  os << token::BEGIN_LIST << l.a_ << token::SPACE << l.b_ << token::END_LIST;
283  return os;
284 }
285 
286 
287 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
288 
289 } // End namespace Foam
290 
291 // ************************ vim: set sw=4 sts=4 et: ************************ //