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