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) 2011 OpenFOAM Foundation
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 "IOstreams.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template<class Point, class PointRef>
31 inline Foam::line<Point, PointRef>::line(const Point& start, const Point& end)
32 :
33  a_(start),
34  b_(end)
35 {}
36 
37 
38 template<class Point, class PointRef>
40 (
41  const UList<Point>& points,
42  const FixedList<label, 2>& indices
43 )
44 :
45  a_(points[indices[0]]),
46  b_(points[indices[1]])
47 {}
48 
49 
50 template<class Point, class PointRef>
52 {
53  is >> *this;
54 }
55 
56 
57 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
58 
59 template<class Point, class PointRef>
60 inline PointRef Foam::line<Point, PointRef>::start() const
61 {
62  return a_;
63 }
64 
65 template<class Point, class PointRef>
66 inline PointRef Foam::line<Point, PointRef>::end() const
67 {
68  return b_;
69 }
70 
71 
72 template<class Point, class PointRef>
74 {
75  return 0.5*(a_ + b_);
76 }
77 
78 
79 template<class Point, class PointRef>
80 inline Foam::scalar Foam::line<Point, PointRef>::mag() const
81 {
83 }
84 
85 
86 template<class Point, class PointRef>
88 {
89  return b_ - a_;
90 }
91 
92 
93 template<class Point, class PointRef>
95 (
96  const Point& p
97 ) const
98 {
99  Point v = vec();
100 
101  Point w(p - a_);
102 
103  scalar c1 = v & w;
104 
105  if (c1 <= 0)
106  {
107  return PointHit<Point>(false, a_, Foam::mag(p - a_), true);
108  }
109 
110  scalar c2 = v & v;
111 
112  if (c2 <= c1)
113  {
114  return PointHit<Point>(false, b_, Foam::mag(p - b_), true);
115  }
116 
117  scalar b = c1/c2;
118 
119  Point pb(a_ + b*v);
120 
121  return PointHit<Point>(true, pb, Foam::mag(p - pb), false);
122 }
123 
124 
125 template<class Point, class PointRef>
127 (
128  const line<Point, const Point&>& edge,
129  Point& thisPt,
130  Point& edgePt
131 ) const
132 {
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());
137 
138  Point crossab = a ^ b;
139  scalar magCrossSqr = magSqr(crossab);
140 
141  if (magCrossSqr > VSMALL)
142  {
143  scalar s = ((c ^ b) & crossab)/magCrossSqr;
144  scalar t = ((c ^ a) & crossab)/magCrossSqr;
145 
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!
157 
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());
162 
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());
167 
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).
200 
201  scalar edge0 = edge.start() & a;
202  scalar edge1 = edge.end() & a;
203  bool edgeOrder = edge0 < edge1;
204 
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());
209 
210  scalar this0 = start() & a;
211  scalar this1 = end() & a;
212  bool thisOrder = this0 < this1;
213 
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());
218 
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  }
255 
256  return Foam::mag(thisPt - edgePt);
257 }
258 
259 
260 // * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * //
261 
262 template<class Point, class PointRef>
263 inline Foam::Istream& Foam::operator>>
264 (
265  Istream& is,
267 )
268 {
269  is.readBegin("line");
270  is >> l.a_ >> l.b_;
271  is.readEnd("line");
272 
273  is.check("Istream& operator>>(Istream&, line&)");
274  return is;
275 }
276 
277 
278 template<class Point, class PointRef>
279 inline Foam::Ostream& Foam::operator<<
280 (
281  Ostream& os,
282  const line<Point, PointRef>& l
283 )
284 {
285  os << token::BEGIN_LIST
286  << l.a_ << token::SPACE
287  << l.b_
288  << token::END_LIST;
289  return os;
290 }
291 
292 
293 // ************************************************************************* //
PointHit< Point > nearestDist(const Point &p) const
Return nearest distance to line from a given point.
Definition: lineI.H:95
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
PointRef start() const
Return first vertex.
Definition: lineI.H:60
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject( name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE ))
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
This class describes the interaction of a face and a point. It carries the info of a successful hit a...
Definition: PointHit.H:51
line(const Point &start, const Point &end)
Construct from two points.
Definition: lineI.H:31
dimensioned< scalar > mag(const dimensioned< Type > &)
PointRef end() const
Return second vertex.
Definition: lineI.H:66
dimensioned< scalar > magSqr(const dimensioned< Type > &)
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
A line primitive.
Definition: line.H:56
virtual bool check(const char *operation) const
Check IOstream status for given operation.
Definition: IOstream.C:92
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
Istream & readBegin(const char *funcName)
Definition: Istream.C:88
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
Point centre() const
Return centre (centroid)
Definition: lineI.H:73
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
scalar mag() const
Return scalar magnitude.
Definition: lineI.H:80
const dimensionedScalar c
Speed of light in a vacuum.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
Point vec() const
Return start-end vector.
Definition: lineI.H:87
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
Istream & readEnd(const char *funcName)
Definition: Istream.C:105