FreeFOAM The Cross-Platform CFD Toolkit
plane.C
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 "plane.H"
27 #include <OpenFOAM/tensor.H>
28 
29 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
30 
31 // Calculate base point and unit normal vector from plane equation
32 void Foam::plane::calcPntAndVec(const scalarList& C)
33 {
34  if (mag(C[0]) > VSMALL)
35  {
36  basePoint_ = vector((-C[3]/C[0]), 0, 0);
37  }
38  else
39  {
40  if (mag(C[1]) > VSMALL)
41  {
42  basePoint_ = vector(0, (-C[3]/C[1]), 0);
43  }
44  else
45  {
46  if (mag(C[2]) > VSMALL)
47  {
48  basePoint_ = vector(0, 0, (-C[3]/C[2]));
49  }
50  else
51  {
52  FatalErrorIn("void plane::calcPntAndVec(const scalarList&)")
53  << "At least one plane coefficient must have a value"
54  << abort(FatalError);
55  }
56  }
57  }
58 
59  unitVector_ = vector(C[0], C[1], C[2]);
60  scalar magUnitVector(mag(unitVector_));
61 
62  if (magUnitVector < VSMALL)
63  {
64  FatalErrorIn("void plane::calcPntAndVec(const scalarList&)")
65  << "Plane normal defined with zero length"
66  << abort(FatalError);
67  }
68 
69  unitVector_ /= magUnitVector;
70 }
71 
72 
73 void Foam::plane::calcPntAndVec
74 (
75  const point& point1,
76  const point& point2,
77  const point& point3
78 )
79 {
80  basePoint_ = (point1 + point2 + point3)/3;
81  vector line12 = point1 - point2;
82  vector line23 = point2 - point3;
83 
84  if
85  (
86  mag(line12) < VSMALL
87  || mag(line23) < VSMALL
88  || mag(point3-point1) < VSMALL
89  )
90  {
92  (
93  "void plane::calcPntAndVec\n"
94  "(\n"
95  " const point&,\n"
96  " const point&,\n"
97  " const point&\n"
98  ")\n"
99  ) << "Bad points." << abort(FatalError);
100  }
101 
102  unitVector_ = line12 ^ line23;
103  scalar magUnitVector(mag(unitVector_));
104 
105  if (magUnitVector < VSMALL)
106  {
108  (
109  "void plane::calcPntAndVec\n"
110  "(\n"
111  " const point&,\n"
112  " const point&,\n"
113  " const point&\n"
114  ")\n"
115  ) << "Plane normal defined with zero length"
116  << abort(FatalError);
117  }
118 
119  unitVector_ /= magUnitVector;
120 }
121 
122 
123 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
124 
125 // Construct from normal vector through the origin
126 Foam::plane::plane(const vector& normalVector)
127 :
128  unitVector_(normalVector),
129  basePoint_(vector::zero)
130 {}
131 
132 
133 // Construct from point and normal vector
134 Foam::plane::plane(const point& basePoint, const vector& normalVector)
135 :
136  unitVector_(normalVector),
137  basePoint_(basePoint)
138 {
139  scalar magUnitVector(mag(unitVector_));
140 
141  if (magUnitVector > VSMALL)
142  {
143  unitVector_ /= magUnitVector;
144  }
145  else
146  {
147  FatalErrorIn("plane::plane(const point&, const vector&)")
148  << "plane normal has zero length"
149  << abort(FatalError);
150  }
151 }
152 
153 
154 // Construct from plane equation
156 {
157  calcPntAndVec(C);
158 }
159 
160 
161 // Construct from three points
163 (
164  const point& a,
165  const point& b,
166  const point& c
167 )
168 {
169  calcPntAndVec(a, b, c);
170 }
171 
172 
173 // Construct from dictionary
175 :
176  unitVector_(vector::zero),
177  basePoint_(point::zero)
178 {
179  word planeType(dict.lookup("planeType"));
180 
181  if (planeType == "planeEquation")
182  {
183  const dictionary& subDict = dict.subDict("planeEquationDict");
184  scalarList C(4);
185 
186  C[0] = readScalar(subDict.lookup("a"));
187  C[1] = readScalar(subDict.lookup("b"));
188  C[2] = readScalar(subDict.lookup("c"));
189  C[3] = readScalar(subDict.lookup("d"));
190 
191  calcPntAndVec(C);
192 
193  }
194  else if (planeType == "embeddedPoints")
195  {
196  const dictionary& subDict = dict.subDict("embeddedPoints");
197 
198  point point1(subDict.lookup("point1"));
199  point point2(subDict.lookup("point2"));
200  point point3(subDict.lookup("point3"));
201 
202  calcPntAndVec(point1, point2, point3);
203  }
204  else if (planeType == "pointAndNormal")
205  {
206  const dictionary& subDict = dict.subDict("pointAndNormalDict");
207 
208  basePoint_ = subDict.lookup("basePoint");
209  unitVector_ = subDict.lookup("normalVector");
210  unitVector_ /= mag(unitVector_);
211  }
212  else
213  {
215  (
216  "plane::plane(const dictionary&)",
217  dict
218  )
219  << "Invalid plane type: " << planeType
220  << abort(FatalIOError);
221  }
222 }
223 
224 
225 // Construct from Istream. Assumes point and normal vector.
227 :
228  unitVector_(is),
229  basePoint_(is)
230 {
231  scalar magUnitVector(mag(unitVector_));
232 
233  if (magUnitVector > VSMALL)
234  {
235  unitVector_ /= magUnitVector;
236  }
237  else
238  {
239  FatalErrorIn("plane::plane(Istream& is)")
240  << "plane normal has zero length"
241  << abort(FatalError);
242  }
243 }
244 
245 
246 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
247 
248 // Return plane normal vector
250 {
251  return unitVector_;
252 }
253 
254 
255 // Return plane base point
257 {
258  return basePoint_;
259 }
260 
261 
262 // Return coefficients for plane equation: ax + by + cz + d = 0
264 {
266 
267  scalar magX = mag(unitVector_.x());
268  scalar magY = mag(unitVector_.y());
269  scalar magZ = mag(unitVector_.z());
270 
271  if (magX > magY)
272  {
273  if (magX > magZ)
274  {
275  C[0] = 1;
276  C[1] = unitVector_.y()/unitVector_.x();
277  C[2] = unitVector_.z()/unitVector_.x();
278  }
279  else
280  {
281  C[0] = unitVector_.x()/unitVector_.z();
282  C[1] = unitVector_.y()/unitVector_.z();
283  C[2] = 1;
284  }
285  }
286  else
287  {
288  if (magY > magZ)
289  {
290  C[0] = unitVector_.x()/unitVector_.y();
291  C[1] = 1;
292  C[2] = unitVector_.z()/unitVector_.y();
293  }
294  else
295  {
296  C[0] = unitVector_.x()/unitVector_.z();
297  C[1] = unitVector_.y()/unitVector_.z();
298  C[2] = 1;
299  }
300  }
301 
302  C[3] = - C[0] * basePoint_.x()
303  - C[1] * basePoint_.y()
304  - C[2] * basePoint_.z();
305 
306  return C;
307 }
308 
309 
310 // Return nearest point in the plane for the given point
312 {
313  return p - unitVector_*((p - basePoint_) & unitVector_);
314 }
315 
316 
317 // Return distance from the given point to the plane
318 Foam::scalar Foam::plane::distance(const point& p) const
319 {
320  return mag((p - basePoint_) & unitVector_);
321 }
322 
323 
324 // Cutting point for plane and line defined by origin and direction
325 Foam::scalar Foam::plane::normalIntersect
326 (
327  const point& pnt0,
328  const vector& dir
329 ) const
330 {
331  scalar denom = stabilise((dir & unitVector_), VSMALL);
332 
333  return ((basePoint_ - pnt0) & unitVector_)/denom;
334 }
335 
336 
337 // Cutting line of two planes
339 {
340  // Mathworld plane-plane intersection. Assume there is a point on the
341  // intersection line with z=0 and solve the two plane equations
342  // for that (now 2x2 equation in x and y)
343  // Better: use either z=0 or x=0 or y=0.
344 
345  const vector& n1 = normal();
346  const vector& n2 = plane2.normal();
347 
348  const point& p1 = refPoint();
349  const point& p2 = plane2.refPoint();
350 
351  scalar n1p1 = n1&p1;
352  scalar n2p2 = n2&p2;
353 
354  vector dir = n1 ^ n2;
355 
356  // Determine zeroed out direction (can be x,y or z) by looking at which
357  // has the largest component in dir.
358  scalar magX = mag(dir.x());
359  scalar magY = mag(dir.y());
360  scalar magZ = mag(dir.z());
361 
362  direction iZero, i1, i2;
363 
364  if (magX > magY)
365  {
366  if (magX > magZ)
367  {
368  iZero = 0;
369  i1 = 1;
370  i2 = 2;
371  }
372  else
373  {
374  iZero = 2;
375  i1 = 0;
376  i2 = 1;
377  }
378  }
379  else
380  {
381  if (magY > magZ)
382  {
383  iZero = 1;
384  i1 = 2;
385  i2 = 0;
386  }
387  else
388  {
389  iZero = 2;
390  i1 = 0;
391  i2 = 1;
392  }
393  }
394 
395  vector pt;
396 
397  pt[iZero] = 0;
398  pt[i1] = (n2[i2]*n1p1 - n1[i2]*n2p2) / (n1[i1]*n2[i2] - n2[i1]*n1[i2]);
399  pt[i2] = (n2[i1]*n1p1 - n1[i1]*n2p2) / (n1[i2]*n2[i1] - n1[i1]*n2[i2]);
400 
401  return ray(pt, dir);
402 }
403 
404 
405 // Cutting point of three planes
407 (
408  const plane& plane2,
409  const plane& plane3
410 ) const
411 {
412  FixedList<scalar, 4> coeffs1(planeCoeffs());
413  FixedList<scalar, 4> coeffs2(plane2.planeCoeffs());
414  FixedList<scalar, 4> coeffs3(plane3.planeCoeffs());
415 
416  tensor a
417  (
418  coeffs1[0],coeffs1[1],coeffs1[2],
419  coeffs2[0],coeffs2[1],coeffs2[2],
420  coeffs3[0],coeffs3[1],coeffs3[2]
421  );
422 
423  vector b(coeffs1[3],coeffs2[3],coeffs3[3]);
424 
425  return (inv(a) & (-b));
426 }
427 
428 
430 {
431  os.writeKeyword("planeType") << "pointAndNormal"
432  << token::END_STATEMENT << nl;
433  os << indent << "pointAndNormalDict" << nl
435  os.writeKeyword("basePoint") << basePoint_ << token::END_STATEMENT << nl;
436  os.writeKeyword("normalVector") << unitVector_ << token::END_STATEMENT
437  << nl;
438  os << decrIndent << indent << token::END_BLOCK << endl;
439 }
440 
441 
442 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
443 
444 bool Foam::operator==(const plane& a, const plane& b)
445 {
446  if (a.basePoint_ == b.basePoint_ && a.unitVector_ == b.unitVector_)
447  {
448  return true;
449  }
450  else
451  {
452  return false;
453  }
454 }
455 
456 bool Foam::operator!=(const plane& a, const plane& b)
457 {
458  return !(a == b);
459 }
460 
461 
462 // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
463 
465 {
466  os << a.unitVector_ << token::SPACE << a.basePoint_;
467 
468  return os;
469 }
470 
471 
472 // ************************ vim: set sw=4 sts=4 et: ************************ //