treeBoundBox.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 "treeBoundBox.H"
27 #include "ListOps.H"
28 #include "OFstream.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 const Foam::scalar Foam::treeBoundBox::great(great);
33 
35 (
36  vector(-great, -great, -great),
37  vector(great, great, great)
38 );
39 
40 
42 (
43  vector(great, great, great),
44  vector(-great, -great, -great)
45 );
46 
47 
49 //- Skip documentation : local scope only
50 const Foam::label facesArray[6][4] =
51 {
52  {0, 4, 6, 2}, // left
53  {1, 3, 7, 5}, // right
54  {0, 1, 5, 4}, // bottom
55  {2, 6, 7, 3}, // top
56  {0, 2, 3, 1}, // back
57  {4, 5, 7, 6} // front
58 };
60 
61 
63 (
64  initListList<face, label, 6, 4>(facesArray)
65 );
66 
67 
69 //- Skip documentation : local scope only
70 const Foam::label edgesArray[12][2] =
71 {
72  {0, 1}, // 0
73  {1, 3},
74  {2, 3}, // 2
75  {0, 2},
76  {4, 5}, // 4
77  {5, 7},
78  {6, 7}, // 6
79  {4, 6},
80  {0, 4}, // 8
81  {1, 5},
82  {3, 7}, // 10
83  {2, 6}
84 };
86 
87 
89 (
90  // initListList<edge, label, 12, 2>(edgesArray)
91  calcEdges(edgesArray)
92 );
93 
94 
96 (
97  calcFaceNormals()
98 );
99 
100 
101 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
102 
103 Foam::edgeList Foam::treeBoundBox::calcEdges(const label edgesArray[12][2])
104 {
105  edgeList edges(12);
106  forAll(edges, edgeI)
107  {
108  edges[edgeI][0] = edgesArray[edgeI][0];
109  edges[edgeI][1] = edgesArray[edgeI][1];
110  }
111  return edges;
112 }
113 
114 
115 Foam::FixedList<Foam::vector, 6> Foam::treeBoundBox::calcFaceNormals()
116 {
117  FixedList<vector, 6> normals;
118  normals[faceId::left] = vector(-1, 0, 0);
119  normals[faceId::right] = vector( 1, 0, 0);
120  normals[faceId::bottom] = vector( 0, -1, 0);
121  normals[faceId::top] = vector( 0, 1, 0);
122  normals[faceId::back] = vector( 0, 0, -1);
123  normals[faceId::front] = vector( 0, 0, 1);
124  return normals;
125 }
126 
127 
128 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
129 
131 :
132  boundBox(points, false)
133 {
134  if (points.empty())
135  {
137  << "cannot find bounding box for zero-sized pointField, "
138  << "returning zero" << endl;
139 
140  return;
141  }
142 }
143 
144 
146 (
147  const UList<point>& points,
148  const labelUList& indices
149 )
150 :
151  boundBox(points, indices, false)
152 {
153  if (points.empty() || indices.empty())
154  {
156  << "cannot find bounding box for zero-sized pointField, "
157  << "returning zero" << endl;
158 
159  return;
160  }
161 }
162 
163 
164 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
165 
167 {
169  pointField& points = tPts.ref();
170 
171  forAll(points, octant)
172  {
173  points[octant] = corner(octant);
174  }
175 
176  return tPts;
177 }
178 
179 
181 {
182  return subBbox(midpoint(), octant);
183 }
184 
185 
187 (
188  const point& mid,
189  const direction octant
190 ) const
191 {
192  if (octant > 7)
193  {
195  << "octant should be [0..7]"
196  << abort(FatalError);
197  }
198 
199  // start with a copy of this bounding box and adjust limits accordingly
200  treeBoundBox subBb(*this);
201  point& bbMin = subBb.min();
202  point& bbMax = subBb.max();
203 
204  if (octant & octantBit::rightHalf)
205  {
206  bbMin.x() = mid.x(); // mid -> max
207  }
208  else
209  {
210  bbMax.x() = mid.x(); // min -> mid
211  }
212 
213  if (octant & octantBit::topHalf)
214  {
215  bbMin.y() = mid.y(); // mid -> max
216  }
217  else
218  {
219  bbMax.y() = mid.y(); // min -> mid
220  }
221 
222  if (octant & octantBit::frontHalf)
223  {
224  bbMin.z() = mid.z(); // mid -> max
225  }
226  else
227  {
228  bbMax.z() = mid.z(); // min -> mid
229  }
230 
231  return subBb;
232 }
233 
234 
236 (
237  const point& overallStart,
238  const vector& overallVec,
239  const point& start,
240  const point& end,
241  point& pt,
242  direction& ptOnFaces
243 ) const
244 {
245  // Sutherlands algorithm:
246  // loop
247  // - start = intersection of line with one of the planes bounding
248  // the bounding box
249  // - stop if start inside bb (return true)
250  // - stop if start and end in same 'half' (e.g. both above bb)
251  // (return false)
252  //
253  // Uses posBits to efficiently determine 'half' in which start and end
254  // point are.
255  //
256  // Note:
257  // - sets coordinate to exact position: e.g. pt.x() = min().x()
258  // since plane intersect routine might have truncation error.
259  // This makes sure that posBits tests 'inside'
260 
261  const direction endBits = posBits(end);
262  pt = start;
263 
264  // Allow maximum of 3 clips.
265  for (label i = 0; i < 4; ++i)
266  {
267  direction ptBits = posBits(pt);
268 
269  if (ptBits == 0)
270  {
271  // pt inside bb
272  ptOnFaces = faceBits(pt);
273  return true;
274  }
275 
276  if ((ptBits & endBits) != 0)
277  {
278  // pt and end in same block outside of bb
279  ptOnFaces = faceBits(pt);
280  return false;
281  }
282 
283  if (ptBits & faceBit::left)
284  {
285  // Intersect with plane V=min, n=-1,0,0
286  if (Foam::mag(overallVec.x()) > vSmall)
287  {
288  scalar s = (min().x() - overallStart.x())/overallVec.x();
289  pt.x() = min().x();
290  pt.y() = overallStart.y() + overallVec.y()*s;
291  pt.z() = overallStart.z() + overallVec.z()*s;
292  }
293  else
294  {
295  // Vector not in x-direction. But still intersecting bb planes.
296  // So must be close - just snap to bb.
297  pt.x() = min().x();
298  }
299  }
300  else if (ptBits & faceBit::right)
301  {
302  // Intersect with plane V=max, n=1,0,0
303  if (Foam::mag(overallVec.x()) > vSmall)
304  {
305  scalar s = (max().x() - overallStart.x())/overallVec.x();
306  pt.x() = max().x();
307  pt.y() = overallStart.y() + overallVec.y()*s;
308  pt.z() = overallStart.z() + overallVec.z()*s;
309  }
310  else
311  {
312  pt.x() = max().x();
313  }
314  }
315  else if (ptBits & faceBit::bottom)
316  {
317  // Intersect with plane V=min, n=0,-1,0
318  if (Foam::mag(overallVec.y()) > vSmall)
319  {
320  scalar s = (min().y() - overallStart.y())/overallVec.y();
321  pt.x() = overallStart.x() + overallVec.x()*s;
322  pt.y() = min().y();
323  pt.z() = overallStart.z() + overallVec.z()*s;
324  }
325  else
326  {
327  pt.x() = min().y();
328  }
329  }
330  else if (ptBits & faceBit::top)
331  {
332  // Intersect with plane V=max, n=0,1,0
333  if (Foam::mag(overallVec.y()) > vSmall)
334  {
335  scalar s = (max().y() - overallStart.y())/overallVec.y();
336  pt.x() = overallStart.x() + overallVec.x()*s;
337  pt.y() = max().y();
338  pt.z() = overallStart.z() + overallVec.z()*s;
339  }
340  else
341  {
342  pt.y() = max().y();
343  }
344  }
345  else if (ptBits & faceBit::back)
346  {
347  // Intersect with plane V=min, n=0,0,-1
348  if (Foam::mag(overallVec.z()) > vSmall)
349  {
350  scalar s = (min().z() - overallStart.z())/overallVec.z();
351  pt.x() = overallStart.x() + overallVec.x()*s;
352  pt.y() = overallStart.y() + overallVec.y()*s;
353  pt.z() = min().z();
354  }
355  else
356  {
357  pt.z() = min().z();
358  }
359  }
360  else if (ptBits & faceBit::front)
361  {
362  // Intersect with plane V=max, n=0,0,1
363  if (Foam::mag(overallVec.z()) > vSmall)
364  {
365  scalar s = (max().z() - overallStart.z())/overallVec.z();
366  pt.x() = overallStart.x() + overallVec.x()*s;
367  pt.y() = overallStart.y() + overallVec.y()*s;
368  pt.z() = max().z();
369  }
370  else
371  {
372  pt.z() = max().z();
373  }
374  }
375  }
376 
377  // Can end up here if the end point is on the edge of the boundBox
378  return true;
379 }
380 
381 
383 (
384  const point& start,
385  const point& end,
386  point& pt
387 ) const
388 {
389  direction ptBits;
390  return intersects(start, end-start, start, end, pt, ptBits);
391 }
392 
393 
394 bool Foam::treeBoundBox::contains(const vector& dir, const point& pt) const
395 {
396  // Compare all components against min and max of bb
397 
398  for (direction cmpt=0; cmpt<3; cmpt++)
399  {
400  if (pt[cmpt] < min()[cmpt])
401  {
402  return false;
403  }
404  else if (pt[cmpt] == min()[cmpt])
405  {
406  // On edge. Outside if direction points outwards.
407  if (dir[cmpt] < 0)
408  {
409  return false;
410  }
411  }
412 
413  if (pt[cmpt] > max()[cmpt])
414  {
415  return false;
416  }
417  else if (pt[cmpt] == max()[cmpt])
418  {
419  // On edge. Outside if direction points outwards.
420  if (dir[cmpt] > 0)
421  {
422  return false;
423  }
424  }
425  }
426 
427  // All components inside bb
428  return true;
429 }
430 
431 
433 {
434  direction faceBits = 0;
435  if (pt.x() == min().x())
436  {
437  faceBits |= faceBit::left;
438  }
439  else if (pt.x() == max().x())
440  {
441  faceBits |= faceBit::right;
442  }
443 
444  if (pt.y() == min().y())
445  {
446  faceBits |= faceBit::bottom;
447  }
448  else if (pt.y() == max().y())
449  {
450  faceBits |= faceBit::top;
451  }
452 
453  if (pt.z() == min().z())
454  {
455  faceBits |= faceBit::back;
456  }
457  else if (pt.z() == max().z())
458  {
459  faceBits |= faceBit::front;
460  }
461  return faceBits;
462 }
463 
464 
466 {
467  direction posBits = 0;
468 
469  if (pt.x() < min().x())
470  {
471  posBits |= faceBit::left;
472  }
473  else if (pt.x() > max().x())
474  {
475  posBits |= faceBit::right;
476  }
477 
478  if (pt.y() < min().y())
479  {
480  posBits |= faceBit::bottom;
481  }
482  else if (pt.y() > max().y())
483  {
484  posBits |= faceBit::top;
485  }
486 
487  if (pt.z() < min().z())
488  {
489  posBits |= faceBit::back;
490  }
491  else if (pt.z() > max().z())
492  {
493  posBits |= faceBit::front;
494  }
495  return posBits;
496 }
497 
498 
500 (
501  const point& pt,
502  point& nearest,
503  point& furthest
504 ) const
505 {
506  scalar nearX, nearY, nearZ;
507  scalar farX, farY, farZ;
508 
509  if (Foam::mag(min().x() - pt.x()) < Foam::mag(max().x() - pt.x()))
510  {
511  nearX = min().x();
512  farX = max().x();
513  }
514  else
515  {
516  nearX = max().x();
517  farX = min().x();
518  }
519 
520  if (Foam::mag(min().y() - pt.y()) < Foam::mag(max().y() - pt.y()))
521  {
522  nearY = min().y();
523  farY = max().y();
524  }
525  else
526  {
527  nearY = max().y();
528  farY = min().y();
529  }
530 
531  if (Foam::mag(min().z() - pt.z()) < Foam::mag(max().z() - pt.z()))
532  {
533  nearZ = min().z();
534  farZ = max().z();
535  }
536  else
537  {
538  nearZ = max().z();
539  farZ = min().z();
540  }
541 
542  nearest = point(nearX, nearY, nearZ);
543  furthest = point(farX, farY, farZ);
544 }
545 
546 
547 Foam::scalar Foam::treeBoundBox::maxDist(const point& pt) const
548 {
549  point near, far;
550  calcExtremities(pt, near, far);
551 
552  return Foam::mag(far - pt);
553 }
554 
555 
557 (
558  const point& pt,
559  const treeBoundBox& other
560 ) const
561 {
562  //
563  // Distance point <-> nearest and furthest away vertex of this
564  //
565 
566  point nearThis, farThis;
567 
568  // get nearest and furthest away vertex
569  calcExtremities(pt, nearThis, farThis);
570 
571  const scalar minDistThis =
572  sqr(nearThis.x() - pt.x())
573  + sqr(nearThis.y() - pt.y())
574  + sqr(nearThis.z() - pt.z());
575  const scalar maxDistThis =
576  sqr(farThis.x() - pt.x())
577  + sqr(farThis.y() - pt.y())
578  + sqr(farThis.z() - pt.z());
579 
580  //
581  // Distance point <-> other
582  //
583 
584  point nearOther, farOther;
585 
586  // get nearest and furthest away vertex
587  other.calcExtremities(pt, nearOther, farOther);
588 
589  const scalar minDistOther =
590  sqr(nearOther.x() - pt.x())
591  + sqr(nearOther.y() - pt.y())
592  + sqr(nearOther.z() - pt.z());
593  const scalar maxDistOther =
594  sqr(farOther.x() - pt.x())
595  + sqr(farOther.y() - pt.y())
596  + sqr(farOther.z() - pt.z());
597 
598  //
599  // Categorise
600  //
601  if (maxDistThis < minDistOther)
602  {
603  // All vertices of this are nearer to point than any vertex of other
604  return -1;
605  }
606  else if (minDistThis > maxDistOther)
607  {
608  // All vertices of this are further from point than any vertex of other
609  return 1;
610  }
611  else
612  {
613  // Mixed bag
614  return 0;
615  }
616 }
617 
618 
619 void Foam::treeBoundBox::writeOBJ(const fileName& fName) const
620 {
621  OFstream str(fName);
622 
623  Info<< "Dumping bounding box " << *this << " as lines to obj file "
624  << str.name() << endl;
625 
626  pointField bbPoints(points());
627 
628  forAll(bbPoints, i)
629  {
630  const point& pt = bbPoints[i];
631 
632  str<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
633  }
634 
636  {
637  const edge& e = treeBoundBox::edges[i];
638 
639  str<< "l " << e[0]+1 << ' ' << e[1]+1 << nl;
640  }
641 }
642 
643 
644 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
645 
647 {
648  return operator==
649  (
650  static_cast<const boundBox&>(a),
651  static_cast<const boundBox&>(b)
652  );
653 }
654 
655 
657 {
658  return !(a == b);
659 }
660 
661 
662 // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
663 
665 {
666  return os << static_cast<const boundBox&>(bb);
667 }
668 
669 
671 {
672  return is >> static_cast<boundBox&>(bb);
673 }
674 
675 
676 // ************************************************************************* //
static const FixedList< vector, 6 > faceNormals
Per face the unit normal.
Definition: treeBoundBox.H:181
boundBox()
Construct null, setting points to zero.
Definition: boundBoxI.H:32
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A class for handling file names.
Definition: fileName.H:79
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:54
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
point corner(const direction) const
Corner point given octant.
Definition: treeBoundBoxI.H:63
static const treeBoundBox greatBox
As per boundBox::greatBox, but with great instead of vGreat.
Definition: treeBoundBox.H:105
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedSymmTensor sqr(const dimensionedVector &dv)
treeBoundBox subBbox(const direction) const
Sub box given by octant number. Midpoint calculated.
Definition: treeBoundBox.C:180
Output to file stream.
Definition: OFstream.H:82
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
direction faceBits(const point &) const
Code position of point on bounding box faces.
Definition: treeBoundBox.C:432
static faceList faces()
Return faces with correct point order.
Definition: boundBox.C:167
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
const Cmpt & z() const
Definition: VectorI.H:87
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const fileName & name() const
Return the name of the stream.
Definition: OFstream.H:120
direction posBits(const point &) const
Position of point relative to bounding box.
Definition: treeBoundBox.C:465
const Cmpt & y() const
Definition: VectorI.H:81
Various functions to operate on Lists.
scalar y
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))
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
bool intersects(const point &overallStart, const vector &overallVec, const point &start, const point &end, point &pt, direction &ptBits) const
Intersects segment; set point to intersection position and face,.
Definition: treeBoundBox.C:236
List< edge > edgeList
Definition: edgeList.H:38
treeBoundBox()
Construct null setting points to zero.
Definition: treeBoundBoxI.H:31
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
tmp< fvMatrix< Type > > operator==(const fvMatrix< Type > &, const fvMatrix< Type > &)
Istream & operator>>(Istream &, directionInfo &)
static const edgeList edges
Edge to point addressing.
Definition: treeBoundBox.H:178
errorManip< error > abort(error &err)
Definition: errorManip.H:131
point midpoint() const
The midpoint of the bounding box.
Definition: boundBoxI.H:78
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
const Cmpt & x() const
Definition: VectorI.H:75
void calcExtremities(const point &pt, point &nearest, point &furthest) const
Calculate nearest and furthest (to point) vertex coords of.
Definition: treeBoundBox.C:500
bool contains(const vector &dir, const point &) const
Contains point (inside or on edge) and moving in direction.
Definition: treeBoundBox.C:394
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
scalar maxDist(const point &) const
Returns distance point to furthest away corner.
Definition: treeBoundBox.C:547
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
label distanceCmp(const point &, const treeBoundBox &other) const
Compare distance to point with other bounding box.
Definition: treeBoundBox.C:557
Ostream & operator<<(Ostream &, const ensightPart &)
static const treeBoundBox invertedBox
As per boundBox::invertedBox, but with great instead of vGreat.
Definition: treeBoundBox.H:108
tmp< pointField > points() const
Vertex coordinates. In octant coding.
Definition: treeBoundBox.C:166
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
A class for managing temporary objects.
Definition: PtrList.H:53
bool operator!=(const particle &, const particle &)
Definition: particle.C:1205
static const scalar great
The great value used for greatBox and invertedBox.
Definition: treeBoundBox.H:102
void writeOBJ(const fileName &fName) const
Definition: treeBoundBox.C:619