Initial commit for version 2.0.x patch release
[OpenFOAM-2.0.x.git] / src / OpenFOAM / meshes / primitiveShapes / line / lineI.H
blobae857900009374f7587b128a661950730854cf74
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
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.
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.
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/>.
24 \*---------------------------------------------------------------------------*/
26 #include "IOstreams.H"
28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
30 template<class Point, class PointRef>
31 inline Foam::line<Point, PointRef>::line(const Point& start, const Point& end)
33     a_(start),
34     b_(end)
38 template<class Point, class PointRef>
39 inline Foam::line<Point, PointRef>::line
41     const UList<Point>& points,
42     const FixedList<label, 2>& indices
45     a_(points[indices[0]]),
46     b_(points[indices[1]])
50 template<class Point, class PointRef>
51 inline Foam::line<Point, PointRef>::line(Istream& is)
53     is  >> *this;
57 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
59 template<class Point, class PointRef>
60 inline PointRef Foam::line<Point, PointRef>::start() const
62     return a_;
65 template<class Point, class PointRef>
66 inline PointRef Foam::line<Point, PointRef>::end() const
68     return b_;
72 template<class Point, class PointRef>
73 inline Point Foam::line<Point, PointRef>::centre() const
75     return 0.5*(a_ + b_);
79 template<class Point, class PointRef>
80 inline Foam::scalar Foam::line<Point, PointRef>::mag() const
82     return ::Foam::mag(vec());
86 template<class Point, class PointRef>
87 inline Point Foam::line<Point, PointRef>::vec() const
89     return b_ - a_;
93 template<class Point, class PointRef>
94 Foam::PointHit<Point> Foam::line<Point, PointRef>::nearestDist
96     const Point& p
97 ) const
99     Point v = vec();
101     Point w(p - a_);
103     scalar c1 = v & w;
105     if (c1 <= 0)
106     {
107         return PointHit<Point>(false, a_, Foam::mag(p - a_), true);
108     }
110     scalar c2 = v & v;
112     if (c2 <= c1)
113     {
114         return PointHit<Point>(false, b_, Foam::mag(p - b_), true);
115     }
117     scalar b = c1/c2;
119     Point pb(a_ + b*v);
121     return PointHit<Point>(true, pb, Foam::mag(p - pb), false);
125 template<class Point, class PointRef>
126 Foam::scalar Foam::line<Point, PointRef>::nearestDist
128     const line<Point, const Point&>& edge,
129     Point& thisPt,
130     Point& edgePt
131 ) const
133     // From Mathworld Line-Line distance/(Gellert et al. 1989, p. 538).
134     Point a(end() - start());
135     Point b(edge.end() - edge.start());
136     Point c(edge.start() - start());
138     Point crossab = a ^ b;
139     scalar magCrossSqr = magSqr(crossab);
141     if (magCrossSqr > VSMALL)
142     {
143         scalar s = ((c ^ b) & crossab)/magCrossSqr;
144         scalar t = ((c ^ a) & crossab)/magCrossSqr;
146         // Check for end points outside of range 0..1
147         if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
148         {
149             // Both inside range 0..1
150             thisPt = start() + a*s;
151             edgePt = edge.start() + b*t;
152         }
153         else
154         {
155             // Do brute force. Distance of everything to everything.
156             // Can quite possibly be improved!
158             // From edge endpoints to *this
159             PointHit<Point> this0(nearestDist(edge.start()));
160             PointHit<Point> this1(nearestDist(edge.end()));
161             scalar thisDist = min(this0.distance(), this1.distance());
163             // From *this to edge
164             PointHit<Point> edge0(edge.nearestDist(start()));
165             PointHit<Point> edge1(edge.nearestDist(end()));
166             scalar edgeDist = min(edge0.distance(), edge1.distance());
168             if (thisDist < edgeDist)
169             {
170                 if (this0.distance() < this1.distance())
171                 {
172                     thisPt = this0.rawPoint();
173                     edgePt = edge.start();
174                 }
175                 else
176                 {
177                     thisPt = this1.rawPoint();
178                     edgePt = edge.end();
179                 }
180             }
181             else
182             {
183                 if (edge0.distance() < edge1.distance())
184                 {
185                     thisPt = start();
186                     edgePt = edge0.rawPoint();
187                 }
188                 else
189                 {
190                     thisPt = end();
191                     edgePt = edge1.rawPoint();
192                 }
193             }
194         }
195     }
196     else
197     {
198         // Parallel lines. Find overlap of both lines by projecting onto
199         // direction vector (now equal for both lines).
201         scalar edge0 = edge.start() & a;
202         scalar edge1 = edge.end() & a;
203         bool edgeOrder = edge0 < edge1;
205         scalar minEdge = (edgeOrder ? edge0 : edge1);
206         scalar maxEdge = (edgeOrder ? edge1 : edge0);
207         const Point& minEdgePt = (edgeOrder ? edge.start() : edge.end());
208         const Point& maxEdgePt = (edgeOrder ? edge.end() : edge.start());
210         scalar this0 = start() & a;
211         scalar this1 = end() & a;
212         bool thisOrder = this0 < this1;
214         scalar minThis = min(this0, this1);
215         scalar maxThis = max(this1, this0);
216         const Point& minThisPt = (thisOrder ? start() : end());
217         const Point& maxThisPt = (thisOrder ? end() : start());
219         if (maxEdge < minThis)
220         {
221             // edge completely below *this
222             edgePt = maxEdgePt;
223             thisPt = minThisPt;
224         }
225         else if (maxEdge < maxThis)
226         {
227             // maxEdge inside interval of *this
228             edgePt = maxEdgePt;
229             thisPt = nearestDist(edgePt).rawPoint();
230         }
231         else
232         {
233             // maxEdge outside. Check if minEdge inside.
234             if (minEdge < minThis)
235             {
236                 // Edge completely envelops this. Take any this point and
237                 // determine nearest on edge.
238                 thisPt = minThisPt;
239                 edgePt = edge.nearestDist(thisPt).rawPoint();
240             }
241             else if (minEdge < maxThis)
242             {
243                 // minEdge inside this interval.
244                 edgePt = minEdgePt;
245                 thisPt = nearestDist(edgePt).rawPoint();
246             }
247             else
248             {
249                 // minEdge outside this interval
250                 edgePt = minEdgePt;
251                 thisPt = maxThisPt;
252             }
253         }
254     }
256     return Foam::mag(thisPt - edgePt);
260 // * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
262 template<class Point, class PointRef>
263 inline Foam::Istream& Foam::operator>>
265     Istream& is,
266     line<Point, PointRef>& l
269     is.readBegin("line");
270     is  >> l.a_ >> l.b_;
271     is.readEnd("line");
273     is.check("Istream& operator>>(Istream&, line&)");
274     return is;
278 template<class Point, class PointRef>
279 inline Foam::Ostream& Foam::operator<<
281     Ostream& os,
282     const line<Point, PointRef>& l
285     os  << token::BEGIN_LIST
286         << l.a_ << token::SPACE
287         << l.b_
288         << token::END_LIST;
289     return os;
293 // ************************************************************************* //