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-2021 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.optionalSubDict("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.optionalSubDict("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.optionalSubDict("pointAndNormalDict");
187 
188  point_ =
190  (
191  {"point", "basePoint"}
192  );
193  normal_ =
194  normalised
195  (
197  (
198  {"normal", "normalVector"}
199  )
200  );
201  }
202  else
203  {
205  << "Invalid plane type: " << planeType << nl
206  << "Valid options include: planeEquation, embeddedPoints and "
207  << "pointAndNormal"
208  << abort(FatalIOError);
209  }
210 }
211 
212 
214 :
215  normal_(is),
216  point_(is)
217 {
218  scalar magUnitVector(mag(normal_));
219 
220  if (magUnitVector > vSmall)
221  {
222  normal_ /= magUnitVector;
223  }
224  else
225  {
227  << "plane normal has zero length. basePoint:" << point_
228  << abort(FatalError);
229  }
230 }
231 
232 
233 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
234 
236 {
237  return normal_;
238 }
239 
240 
242 {
243  return point_;
244 }
245 
246 
248 {
250 
251  scalar magX = mag(normal_.x());
252  scalar magY = mag(normal_.y());
253  scalar magZ = mag(normal_.z());
254 
255  if (magX > magY)
256  {
257  if (magX > magZ)
258  {
259  C[0] = 1;
260  C[1] = normal_.y()/normal_.x();
261  C[2] = normal_.z()/normal_.x();
262  }
263  else
264  {
265  C[0] = normal_.x()/normal_.z();
266  C[1] = normal_.y()/normal_.z();
267  C[2] = 1;
268  }
269  }
270  else
271  {
272  if (magY > magZ)
273  {
274  C[0] = normal_.x()/normal_.y();
275  C[1] = 1;
276  C[2] = normal_.z()/normal_.y();
277  }
278  else
279  {
280  C[0] = normal_.x()/normal_.z();
281  C[1] = normal_.y()/normal_.z();
282  C[2] = 1;
283  }
284  }
285 
286  C[3] = - C[0] * point_.x()
287  - C[1] * point_.y()
288  - C[2] * point_.z();
289 
290  return C;
291 }
292 
293 
295 {
296  // Perturb base point
297  const point& refPt = refPoint();
298 
299  // ax + by + cz + d = 0
301 
302  const scalar perturbX = refPt.x() + 1e-3;
303  const scalar perturbY = refPt.y() + 1e-3;
304  const scalar perturbZ = refPt.z() + 1e-3;
305 
306  if (mag(planeCoeffs[2]) < small)
307  {
308  if (mag(planeCoeffs[1]) < small)
309  {
310  const scalar x =
311  -1.0
312  *(
313  planeCoeffs[3]
314  + planeCoeffs[1]*perturbY
315  + planeCoeffs[2]*perturbZ
316  )/planeCoeffs[0];
317 
318  return point
319  (
320  x,
321  perturbY,
322  perturbZ
323  );
324  }
325 
326  const scalar y =
327  -1.0
328  *(
329  planeCoeffs[3]
330  + planeCoeffs[0]*perturbX
331  + planeCoeffs[2]*perturbZ
332  )/planeCoeffs[1];
333 
334  return point
335  (
336  perturbX,
337  y,
338  perturbZ
339  );
340  }
341  else
342  {
343  const scalar z =
344  -1.0
345  *(
346  planeCoeffs[3]
347  + planeCoeffs[0]*perturbX
348  + planeCoeffs[1]*perturbY
349  )/planeCoeffs[2];
350 
351  return point
352  (
353  perturbX,
354  perturbY,
355  z
356  );
357  }
358 }
359 
360 
362 {
363  return p - normal_*((p - point_) & normal_);
364 }
365 
366 
367 Foam::scalar Foam::plane::distance(const point& p) const
368 {
369  return mag((p - point_) & normal_);
370 }
371 
372 
373 Foam::scalar Foam::plane::normalIntersect
374 (
375  const point& pnt0,
376  const vector& dir
377 ) const
378 {
379  const scalar num = (point_ - pnt0) & normal_;
380  const scalar den = dir & normal_;
381 
382  return mag(den) > mag(num)*vSmall ? num/den : vGreat;
383 }
384 
385 
387 {
388  // Mathworld plane-plane intersection. Assume there is a point on the
389  // intersection line with z=0 and solve the two plane equations
390  // for that (now 2x2 equation in x and y)
391  // Better: use either z=0 or x=0 or y=0.
392 
393  const vector& n1 = normal();
394  const vector& n2 = plane2.normal();
395 
396  const point& p1 = refPoint();
397  const point& p2 = plane2.refPoint();
398 
399  scalar n1p1 = n1&p1;
400  scalar n2p2 = n2&p2;
401 
402  vector dir = n1 ^ n2;
403 
404  // Determine zeroed out direction (can be x,y or z) by looking at which
405  // has the largest component in dir.
406  scalar magX = mag(dir.x());
407  scalar magY = mag(dir.y());
408  scalar magZ = mag(dir.z());
409 
410  direction iZero, i1, i2;
411 
412  if (magX > magY)
413  {
414  if (magX > magZ)
415  {
416  iZero = 0;
417  i1 = 1;
418  i2 = 2;
419  }
420  else
421  {
422  iZero = 2;
423  i1 = 0;
424  i2 = 1;
425  }
426  }
427  else
428  {
429  if (magY > magZ)
430  {
431  iZero = 1;
432  i1 = 2;
433  i2 = 0;
434  }
435  else
436  {
437  iZero = 2;
438  i1 = 0;
439  i2 = 1;
440  }
441  }
442 
443  vector pt;
444 
445  pt[iZero] = 0;
446  pt[i1] = (n2[i2]*n1p1 - n1[i2]*n2p2) / (n1[i1]*n2[i2] - n2[i1]*n1[i2]);
447  pt[i2] = (n2[i1]*n1p1 - n1[i1]*n2p2) / (n1[i2]*n2[i1] - n1[i1]*n2[i2]);
448 
449  return ray(pt, dir);
450 }
451 
452 
454 (
455  const plane& plane2,
456  const plane& plane3
457 ) const
458 {
460  FixedList<scalar, 4> coeffs2(plane2.planeCoeffs());
461  FixedList<scalar, 4> coeffs3(plane3.planeCoeffs());
462 
463  tensor a
464  (
465  coeffs1[0],coeffs1[1],coeffs1[2],
466  coeffs2[0],coeffs2[1],coeffs2[2],
467  coeffs3[0],coeffs3[1],coeffs3[2]
468  );
469 
470  vector b(coeffs1[3],coeffs2[3],coeffs3[3]);
471 
472  return (inv(a) & (-b));
473 }
474 
475 
477 {
478  const scalar angle((p - point_) & normal_);
479 
480  return (angle < 0 ? FLIP : NORMAL);
481 }
482 
483 
485 {
486  const vector mirroredPtDir = p - nearestPoint(p);
487 
488  if ((normal() & mirroredPtDir) > 0)
489  {
490  return p - 2.0*distance(p)*normal();
491  }
492  else
493  {
494  return p + 2.0*distance(p)*normal();
495  }
496 }
497 
498 
500 {
501  writeEntry(os, "planeType", "pointAndNormal");
502  os << indent << "pointAndNormalDict" << nl
504  writeEntry(os, "point", point_);
505  writeEntry(os, "normal", normal_);
506  os << decrIndent << indent << token::END_BLOCK << endl;
507 }
508 
509 
510 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
511 
512 bool Foam::operator==(const plane& a, const plane& b)
513 {
514  if (a.point_ == b.point_ && a.normal_ == b.normal_)
515  {
516  return true;
517  }
518  else
519  {
520  return false;
521  }
522 }
523 
524 bool Foam::operator!=(const plane& a, const plane& b)
525 {
526  return !(a == b);
527 }
528 
529 
530 // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
531 
533 {
534  os << a.normal_ << token::SPACE << a.point_;
535 
536  return os;
537 }
538 
539 
540 // ************************************************************************* //
dictionary dict
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:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
point nearestPoint(const point &p) const
Return nearest point in the plane for the given point.
Definition: plane.C:361
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
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
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:386
point planePlaneIntersect(const plane &, const plane &) const
Return the cutting point between this plane and two other planes.
Definition: plane.C:454
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
side sideOfPlane(const point &p) const
Return the side of the plane that the point is on.
Definition: plane.C:476
scalar y
scalar normalIntersect(const point &pnt0, const vector &dir) const
Return cut coefficient for plane and line defined by.
Definition: plane.C:374
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
const dictionary & optionalSubDict(const word &) const
Find and return a sub-dictionary if found.
Definition: dictionary.C:1040
point aPoint() const
Return a point on the plane.
Definition: plane.C:294
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:241
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:367
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:484
side
Side of the plane.
Definition: plane.H:65
const vector & normal() const
Return plane normal.
Definition: plane.C:235
vector point
Point is a vector.
Definition: point.H:41
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:335
Ostream & operator<<(Ostream &, const ensightPart &)
void writeDict(Ostream &) const
Write to dictionary.
Definition: plane.C:499
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:247
volScalarField & p
bool operator!=(const particle &, const particle &)
Definition: particle.C:1205
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:228
ITstream & lookupBackwardsCompatible(const wordList &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream, trying a list of keywords.
Definition: dictionary.C:855
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844
IOerror FatalIOError