OpenFOAM logo
Open Source CFD Toolkit

lineI.H

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2005 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software; you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by the
00013     Free Software Foundation; either version 2 of the License, or (at your
00014     option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM; if not, write to the Free Software Foundation,
00023     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
00024 
00025 Description
00026 
00027 \*---------------------------------------------------------------------------*/
00028 
00029 #include "IOstreams.H"
00030 #include "pointHit.H"
00031 
00032 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00033 
00034 namespace Foam
00035 {
00036 
00037 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00038 
00039 template<class Point, class PointRef>
00040 inline line<Point, PointRef>::line(const Point& start, const Point& end)
00041 :
00042     a_(start),
00043     b_(end)
00044 {}
00045 
00046 
00047 template<class Point, class PointRef>
00048 inline line<Point, PointRef>::line(Istream& is)
00049 {
00050     // Read beginning of line point pair
00051     is.readBegin("line");
00052 
00053     is >> a_ >> b_;
00054 
00055     // Read end of line point pair
00056     is.readEnd("line");
00057 
00058     // Check state of Istream
00059     is.check("line::line(Istream& is)");
00060 }
00061 
00062 
00063 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00064 
00065 template<class Point, class PointRef>
00066 inline PointRef line<Point, PointRef>::start() const
00067 {
00068     return a_;
00069 }
00070 
00071 template<class Point, class PointRef>
00072 inline PointRef line<Point, PointRef>::end() const
00073 {
00074     return b_;
00075 }
00076 
00077 
00078 template<class Point, class PointRef>
00079 inline Point line<Point, PointRef>::centre() const
00080 {
00081     return 0.5*(a_ + b_);
00082 }
00083 
00084 
00085 template<class Point, class PointRef>
00086 inline scalar line<Point, PointRef>::mag() const
00087 {
00088     return ::Foam::mag(vec());
00089 }
00090 
00091 
00092 template<class Point, class PointRef>
00093 inline vector line<Point, PointRef>::vec() const
00094 {
00095     return b_ - a_;
00096 }
00097 
00098 
00099 template<class Point, class PointRef>
00100 pointHit line<Point, PointRef>::nearestDist
00101 (
00102     const point& p
00103 ) const
00104 {
00105     vector v = vec();
00106 
00107     vector w(p - a_);
00108 
00109     scalar c1 = v & w;
00110 
00111     if (c1 <= 0)
00112     {
00113         return pointHit(false, a_, Foam::mag(p - a_), true);
00114     }
00115 
00116     scalar c2 = v & v;
00117 
00118     if (c2 <= c1)
00119     {
00120         return pointHit(false, b_, Foam::mag(p - b_), true);
00121     }
00122 
00123     scalar b = c1/c2;
00124 
00125     point pb(a_ + b*v);
00126 
00127     return pointHit(true, pb, Foam::mag(p - pb), false);
00128 }
00129 
00130 
00131 template<class Point, class PointRef>
00132 scalar line<Point, PointRef>::nearestDist
00133 (
00134     const line<point, const point&>& edge,
00135     point& thisPt,
00136     point& edgePt
00137 ) const
00138 {
00139     // From Mathworld Line-Line distance/(Gellert et al. 1989, p. 538).
00140     vector a(end() - start());
00141     vector b(edge.end() - edge.start());
00142     vector c(edge.start() - start());
00143 
00144     vector crossab = a ^ b;
00145     scalar magCrossSqr = magSqr(crossab);
00146 
00147     if (magCrossSqr > VSMALL)
00148     {
00149         vector crosscb = c ^ b;
00150 
00151         scalar s = ((c ^ b) & crossab)/magCrossSqr;
00152         scalar t = ((c ^ a) & crossab)/magCrossSqr;
00153 
00154         // Check for end points outside of range 0..1
00155         if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
00156         {
00157             // Both inside range 0..1
00158             thisPt = start() + a*s;
00159             edgePt = edge.start() + b*t;
00160         }
00161         else
00162         {
00163             // Do brute force. Distance of everything to everything.
00164             // Can quite possibly be improved!
00165 
00166             // From edge endpoints to *this
00167             pointHit this0(nearestDist(edge.start()));
00168             pointHit this1(nearestDist(edge.end()));
00169             scalar thisDist = min(this0.distance(), this1.distance());
00170 
00171             // From *this to edge
00172             pointHit edge0(edge.nearestDist(start()));
00173             pointHit edge1(edge.nearestDist(end()));
00174             scalar edgeDist = min(edge0.distance(), edge1.distance());
00175 
00176             if (thisDist < edgeDist)
00177             {
00178                 if (this0.distance() < this1.distance())
00179                 {
00180                     thisPt = this0.rawPoint();
00181                     edgePt = edge.start();
00182                 }
00183                 else
00184                 {
00185                     thisPt = this1.rawPoint();
00186                     edgePt = edge.end();
00187                 }
00188             }
00189             else
00190             {
00191                 if (edge0.distance() < edge1.distance())
00192                 {
00193                     thisPt = start();
00194                     edgePt = edge0.rawPoint();
00195                 }
00196                 else
00197                 {
00198                     thisPt = end();
00199                     edgePt = edge1.rawPoint();
00200                 }
00201             }
00202         }
00203     }
00204     else
00205     {
00206         // Parallel lines. Find overlap of both lines by projecting onto
00207         // direction vector (now equal for both lines).
00208 
00209         scalar edge0 = edge.start() & a;
00210         scalar edge1 = edge.end() & a;
00211         bool edgeOrder = edge0 < edge1;
00212 
00213         scalar minEdge = (edgeOrder ? edge0 : edge1);
00214         scalar maxEdge = (edgeOrder ? edge1 : edge0);
00215         const point& minEdgePt = (edgeOrder ? edge.start() : edge.end());
00216         const point& maxEdgePt = (edgeOrder ? edge.end() : edge.start());
00217 
00218         scalar this0 = start() & a;
00219         scalar this1 = end() & a;
00220         bool thisOrder = this0 < this1;
00221 
00222         scalar minThis = min(this0, this1);
00223         scalar maxThis = max(this1, this0);
00224         const point& minThisPt = (thisOrder ? start() : end());
00225         const point& maxThisPt = (thisOrder ? end() : start());
00226 
00227         if (maxEdge < minThis)
00228         {
00229             // edge completely below *this
00230             edgePt = maxEdgePt;
00231             thisPt = minThisPt;
00232         }
00233         else if (maxEdge < maxThis)
00234         {
00235             // maxEdge inside interval of *this
00236             edgePt = maxEdgePt;
00237             thisPt = nearestDist(edgePt).rawPoint();
00238         }
00239         else
00240         {
00241             // maxEdge outside. Check if minEdge inside.
00242             if (minEdge < minThis)
00243             {
00244                 // Edge completely envelops this. Take any this point and
00245                 // determine nearest on edge.
00246                 thisPt = minThisPt;
00247                 edgePt = edge.nearestDist(thisPt).rawPoint();
00248             }
00249             else if (minEdge < maxThis)
00250             {
00251                 // minEdge inside this interval.
00252                 edgePt = minEdgePt;
00253                 thisPt = nearestDist(edgePt).rawPoint();
00254             }
00255             else
00256             {
00257                 // minEdge outside this interval
00258                 edgePt = minEdgePt;
00259                 thisPt = maxThisPt;
00260             }
00261         }
00262     }
00263 
00264     return Foam::mag(thisPt - edgePt);
00265 }
00266 
00267 
00268 // * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
00269 
00270 template<class point, class pointRef>
00271 inline Istream& operator>>(Istream& is, line<point, pointRef>& l)
00272 {
00273     // Read beginning of line point pair
00274     is.readBegin("line");
00275 
00276     is >> l.a_ >> l.b_;
00277 
00278     // Read end of line point pair
00279     is.readEnd("line");
00280 
00281     // Check state of Ostream
00282     is.check("Istream& operator>>(Istream&, line&)");
00283 
00284     return is;
00285 }
00286 
00287 
00288 template<class Point, class PointRef>
00289 inline Ostream& operator<<(Ostream& os, const line<Point, PointRef>& l)
00290 {
00291     os << token::BEGIN_LIST << l.a_ << token::SPACE << l.b_ << token::END_LIST;
00292     return os;
00293 }
00294 
00295 
00296 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00297 
00298 } // End namespace Foam
00299 
00300 // ************************************************************************* //
For further information go to www.openfoam.org