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