All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
plane.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 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 "plane.H"
27 #include "tensor.H"
28 
29 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
30 
31 void Foam::plane::calcPntAndVec(const scalarList& C)
32 {
33  normal_ = vector(C[0], C[1], C[2]);
34 
35  const scalar magNormal = mag(normal_);
36 
37  if (magNormal == 0)
38  {
40  << "Plane normal has zero length"
41  << abort(FatalError);
42  }
43 
44  normal_ /= magNormal;
45 
46  if (magNormal < mag(C[3])*vSmall)
47  {
49  << "Plane is too far from the origin"
50  << abort(FatalError);
51  }
52 
53  point_ = - C[3]/magNormal*normal_;
54 }
55 
56 
57 void Foam::plane::calcPntAndVec
58 (
59  const point& point1,
60  const point& point2,
61  const point& point3
62 )
63 {
64  point_ = (point1 + point2 + point3)/3;
65  vector line12 = point1 - point2;
66  vector line23 = point2 - point3;
67 
68  if
69  (
70  mag(line12) < vSmall
71  || mag(line23) < vSmall
72  || mag(point3-point1) < vSmall
73  )
74  {
76  << "Bad points:" << point1 << ' ' << point2 << ' ' << point3
77  << abort(FatalError);
78  }
79 
80  normal_ = line12 ^ line23;
81  scalar magUnitVector(mag(normal_));
82 
83  if (magUnitVector < vSmall)
84  {
86  << "Plane normal defined with zero length" << nl
87  << "Bad points:" << point1 << ' ' << point2 << ' ' << point3
88  << abort(FatalError);
89  }
90 
91  normal_ /= magUnitVector;
92 }
93 
94 
95 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
96 
97 Foam::plane::plane(const vector& normalVector)
98 :
99  normal_(normalVector),
100  point_(Zero)
101 {
102  scalar magUnitVector(mag(normal_));
103 
104  if (magUnitVector > vSmall)
105  {
106  normal_ /= magUnitVector;
107  }
108  else
109  {
111  << "plane normal has zero length. basePoint:" << point_
112  << abort(FatalError);
113  }
114 }
115 
116 
117 Foam::plane::plane(const point& basePoint, const vector& normalVector)
118 :
119  normal_(normalVector),
120  point_(basePoint)
121 {
122  scalar magUnitVector(mag(normal_));
123 
124  if (magUnitVector > vSmall)
125  {
126  normal_ /= magUnitVector;
127  }
128  else
129  {
131  << "plane normal has zero length. basePoint:" << point_
132  << abort(FatalError);
133  }
134 }
135 
136 
138 {
139  calcPntAndVec(C);
140 }
141 
142 
144 (
145  const point& a,
146  const point& b,
147  const point& c
148 )
149 {
150  calcPntAndVec(a, b, c);
151 }
152 
153 
155 :
156  normal_(Zero),
157  point_(Zero)
158 {
159  const word planeType(dict.lookup("planeType"));
160 
161  if (planeType == "planeEquation")
162  {
163  const dictionary& subDict = dict.subDict("planeEquationDict");
164  scalarList C(4);
165 
166  C[0] = subDict.lookup<scalar>("a");
167  C[1] = subDict.lookup<scalar>("b");
168  C[2] = subDict.lookup<scalar>("c");
169  C[3] = subDict.lookup<scalar>("d");
170 
171  calcPntAndVec(C);
172 
173  }
174  else if (planeType == "embeddedPoints")
175  {
176  const dictionary& subDict = dict.subDict("embeddedPointsDict");
177 
178  point point1(subDict.lookup("point1"));
179  point point2(subDict.lookup("point2"));
180  point point3(subDict.lookup("point3"));
181 
182  calcPntAndVec(point1, point2, point3);
183  }
184  else if (planeType == "pointAndNormal")
185  {
186  const dictionary& subDict = dict.subDict("pointAndNormalDict");
187 
188  point_ =
189  subDict.found("basePoint")
190  ? subDict.lookup("basePoint")
191  : subDict.lookup("point");
192 
193  normal_ =
194  subDict.found("normalVector")
195  ? subDict.lookup("normalVector")
196  : subDict.lookup("normal");
197 
198  normal_ /= mag(normal_);
199  }
200  else
201  {
203  << "Invalid plane type: " << planeType << nl
204  << "Valid options include: planeEquation, embeddedPoints and "
205  << "pointAndNormal"
206  << abort(FatalIOError);
207  }
208 }
209 
210 
212 :
213  normal_(is),
214  point_(is)
215 {
216  scalar magUnitVector(mag(normal_));
217 
218  if (magUnitVector > vSmall)
219  {
220  normal_ /= magUnitVector;
221  }
222  else
223  {
225  << "plane normal has zero length. basePoint:" << point_
226  << abort(FatalError);
227  }
228 }
229 
230 
231 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
232 
234 {
235  return normal_;
236 }
237 
238 
240 {
241  return point_;
242 }
243 
244 
246 {
248 
249  scalar magX = mag(normal_.x());
250  scalar magY = mag(normal_.y());
251  scalar magZ = mag(normal_.z());
252 
253  if (magX > magY)
254  {
255  if (magX > magZ)
256  {
257  C[0] = 1;
258  C[1] = normal_.y()/normal_.x();
259  C[2] = normal_.z()/normal_.x();
260  }
261  else
262  {
263  C[0] = normal_.x()/normal_.z();
264  C[1] = normal_.y()/normal_.z();
265  C[2] = 1;
266  }
267  }
268  else
269  {
270  if (magY > magZ)
271  {
272  C[0] = normal_.x()/normal_.y();
273  C[1] = 1;
274  C[2] = normal_.z()/normal_.y();
275  }
276  else
277  {
278  C[0] = normal_.x()/normal_.z();
279  C[1] = normal_.y()/normal_.z();
280  C[2] = 1;
281  }
282  }
283 
284  C[3] = - C[0] * point_.x()
285  - C[1] * point_.y()
286  - C[2] * point_.z();
287 
288  return C;
289 }
290 
291 
293 {
294  // Perturb base point
295  const point& refPt = refPoint();
296 
297  // ax + by + cz + d = 0
299 
300  const scalar perturbX = refPt.x() + 1e-3;
301  const scalar perturbY = refPt.y() + 1e-3;
302  const scalar perturbZ = refPt.z() + 1e-3;
303 
304  if (mag(planeCoeffs[2]) < small)
305  {
306  if (mag(planeCoeffs[1]) < small)
307  {
308  const scalar x =
309  -1.0
310  *(
311  planeCoeffs[3]
312  + planeCoeffs[1]*perturbY
313  + planeCoeffs[2]*perturbZ
314  )/planeCoeffs[0];
315 
316  return point
317  (
318  x,
319  perturbY,
320  perturbZ
321  );
322  }
323 
324  const scalar y =
325  -1.0
326  *(
327  planeCoeffs[3]
328  + planeCoeffs[0]*perturbX
329  + planeCoeffs[2]*perturbZ
330  )/planeCoeffs[1];
331 
332  return point
333  (
334  perturbX,
335  y,
336  perturbZ
337  );
338  }
339  else
340  {
341  const scalar z =
342  -1.0
343  *(
344  planeCoeffs[3]
345  + planeCoeffs[0]*perturbX
346  + planeCoeffs[1]*perturbY
347  )/planeCoeffs[2];
348 
349  return point
350  (
351  perturbX,
352  perturbY,
353  z
354  );
355  }
356 }
357 
358 
360 {
361  return p - normal_*((p - point_) & normal_);
362 }
363 
364 
365 Foam::scalar Foam::plane::distance(const point& p) const
366 {
367  return mag((p - point_) & normal_);
368 }
369 
370 
371 Foam::scalar Foam::plane::normalIntersect
372 (
373  const point& pnt0,
374  const vector& dir
375 ) const
376 {
377  const scalar num = (point_ - pnt0) & normal_;
378  const scalar den = dir & normal_;
379 
380  return mag(den) > mag(num)*vSmall ? num/den : vGreat;
381 }
382 
383 
385 {
386  // Mathworld plane-plane intersection. Assume there is a point on the
387  // intersection line with z=0 and solve the two plane equations
388  // for that (now 2x2 equation in x and y)
389  // Better: use either z=0 or x=0 or y=0.
390 
391  const vector& n1 = normal();
392  const vector& n2 = plane2.normal();
393 
394  const point& p1 = refPoint();
395  const point& p2 = plane2.refPoint();
396 
397  scalar n1p1 = n1&p1;
398  scalar n2p2 = n2&p2;
399 
400  vector dir = n1 ^ n2;
401 
402  // Determine zeroed out direction (can be x,y or z) by looking at which
403  // has the largest component in dir.
404  scalar magX = mag(dir.x());
405  scalar magY = mag(dir.y());
406  scalar magZ = mag(dir.z());
407 
408  direction iZero, i1, i2;
409 
410  if (magX > magY)
411  {
412  if (magX > magZ)
413  {
414  iZero = 0;
415  i1 = 1;
416  i2 = 2;
417  }
418  else
419  {
420  iZero = 2;
421  i1 = 0;
422  i2 = 1;
423  }
424  }
425  else
426  {
427  if (magY > magZ)
428  {
429  iZero = 1;
430  i1 = 2;
431  i2 = 0;
432  }
433  else
434  {
435  iZero = 2;
436  i1 = 0;
437  i2 = 1;
438  }
439  }
440 
441  vector pt;
442 
443  pt[iZero] = 0;
444  pt[i1] = (n2[i2]*n1p1 - n1[i2]*n2p2) / (n1[i1]*n2[i2] - n2[i1]*n1[i2]);
445  pt[i2] = (n2[i1]*n1p1 - n1[i1]*n2p2) / (n1[i2]*n2[i1] - n1[i1]*n2[i2]);
446 
447  return ray(pt, dir);
448 }
449 
450 
452 (
453  const plane& plane2,
454  const plane& plane3
455 ) const
456 {
458  FixedList<scalar, 4> coeffs2(plane2.planeCoeffs());
459  FixedList<scalar, 4> coeffs3(plane3.planeCoeffs());
460 
461  tensor a
462  (
463  coeffs1[0],coeffs1[1],coeffs1[2],
464  coeffs2[0],coeffs2[1],coeffs2[2],
465  coeffs3[0],coeffs3[1],coeffs3[2]
466  );
467 
468  vector b(coeffs1[3],coeffs2[3],coeffs3[3]);
469 
470  return (inv(a) & (-b));
471 }
472 
473 
475 {
476  const scalar angle((p - point_) & normal_);
477 
478  return (angle < 0 ? FLIP : NORMAL);
479 }
480 
481 
483 {
484  const vector mirroredPtDir = p - nearestPoint(p);
485 
486  if ((normal() & mirroredPtDir) > 0)
487  {
488  return p - 2.0*distance(p)*normal();
489  }
490  else
491  {
492  return p + 2.0*distance(p)*normal();
493  }
494 }
495 
496 
498 {
499  writeEntry(os, "planeType", "pointAndNormal");
500  os << indent << "pointAndNormalDict" << nl
502  writeEntry(os, "point", point_);
503  writeEntry(os, "normal", normal_);
504  os << decrIndent << indent << token::END_BLOCK << endl;
505 }
506 
507 
508 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
509 
510 bool Foam::operator==(const plane& a, const plane& b)
511 {
512  if (a.point_ == b.point_ && a.normal_ == b.normal_)
513  {
514  return true;
515  }
516  else
517  {
518  return false;
519  }
520 }
521 
522 bool Foam::operator!=(const plane& a, const plane& b)
523 {
524  return !(a == b);
525 }
526 
527 
528 // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
529 
531 {
532  os << a.normal_ << token::SPACE << a.point_;
533 
534  return os;
535 }
536 
537 
538 // ************************************************************************* //
dictionary dict
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:667
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:221
A direction and a reference point.
Definition: plane.H:72
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:54
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
point nearestPoint(const point &p) const
Return nearest point in the plane for the given point.
Definition: plane.C:359
uint8_t direction
Definition: direction.H:45
An Istream is an abstract base class for all input systems (streams, files, token lists etc)...
Definition: Istream.H:57
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
ray planeIntersect(const plane &) const
Return the cutting line between this plane and another.
Definition: plane.C:384
point planePlaneIntersect(const plane &, const plane &) const
Return the cutting point between this plane and two other planes.
Definition: plane.C:452
const Cmpt & z() const
Definition: VectorI.H:87
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
plane(const vector &normalVector)
Construct from normal vector through the origin.
Definition: plane.C:97
const Cmpt & y() const
Definition: VectorI.H:81
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:934
side sideOfPlane(const point &p) const
Return the side of the plane that the point is on.
Definition: plane.C:474
scalar y
scalar normalIntersect(const point &pnt0, const vector &dir) const
Return cut coefficient for plane and line defined by.
Definition: plane.C:372
const dimensionedScalar & b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
point aPoint() const
Return a point on the plane.
Definition: plane.C:292
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
A class for handling words, derived from string.
Definition: word.H:59
const point & refPoint() const
Return or return plane base point.
Definition: plane.C:239
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
scalar distance(const point &p) const
Return distance from the given point to the plane.
Definition: plane.C:365
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Cmpt & x() const
Definition: VectorI.H:75
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
static const char nl
Definition: Ostream.H:260
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:235
void writeEntry(Ostream &os, const HashTable< T, Key, Hash > &ht)
Definition: HashTableIO.C:96
point mirror(const point &p) const
Mirror the supplied point in the plane. Return the mirrored point.
Definition: plane.C:482
side
Side of the plane.
Definition: plane.H:65
const vector & normal() const
Return plane normal.
Definition: plane.C:233
vector point
Point is a vector.
Definition: point.H:41
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:331
Ostream & operator<<(Ostream &, const ensightPart &)
void writeDict(Ostream &) const
Write to dictionary.
Definition: plane.C:497
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
FixedList< scalar, 4 > planeCoeffs() const
Return coefficients for the.
Definition: plane.C:245
volScalarField & p
bool operator!=(const particle &, const particle &)
Definition: particle.C:1204
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:228
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:812
IOerror FatalIOError