00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "IOstreams.H"
00030 #include "pointHit.H"
00031
00032
00033
00034 namespace Foam
00035 {
00036
00037
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
00051 is.readBegin("line");
00052
00053 is >> a_ >> b_;
00054
00055
00056 is.readEnd("line");
00057
00058
00059 is.check("line::line(Istream& is)");
00060 }
00061
00062
00063
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
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
00155 if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
00156 {
00157
00158 thisPt = start() + a*s;
00159 edgePt = edge.start() + b*t;
00160 }
00161 else
00162 {
00163
00164
00165
00166
00167 pointHit this0(nearestDist(edge.start()));
00168 pointHit this1(nearestDist(edge.end()));
00169 scalar thisDist = min(this0.distance(), this1.distance());
00170
00171
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
00207
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
00230 edgePt = maxEdgePt;
00231 thisPt = minThisPt;
00232 }
00233 else if (maxEdge < maxThis)
00234 {
00235
00236 edgePt = maxEdgePt;
00237 thisPt = nearestDist(edgePt).rawPoint();
00238 }
00239 else
00240 {
00241
00242 if (minEdge < minThis)
00243 {
00244
00245
00246 thisPt = minThisPt;
00247 edgePt = edge.nearestDist(thisPt).rawPoint();
00248 }
00249 else if (minEdge < maxThis)
00250 {
00251
00252 edgePt = minEdgePt;
00253 thisPt = nearestDist(edgePt).rawPoint();
00254 }
00255 else
00256 {
00257
00258 edgePt = minEdgePt;
00259 thisPt = maxThisPt;
00260 }
00261 }
00262 }
00263
00264 return Foam::mag(thisPt - edgePt);
00265 }
00266
00267
00268
00269
00270 template<class point, class pointRef>
00271 inline Istream& operator>>(Istream& is, line<point, pointRef>& l)
00272 {
00273
00274 is.readBegin("line");
00275
00276 is >> l.a_ >> l.b_;
00277
00278
00279 is.readEnd("line");
00280
00281
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 }
00299
00300