treeBoundBoxI.H
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-2018 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 "Random.H"
28 
29 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
30 
32 :
33  boundBox()
34 {}
35 
36 
38 :
39  boundBox(min, max)
40 {}
41 
42 
44 :
45  boundBox(bb)
46 {}
47 
48 
50 :
51  boundBox(is)
52 {}
53 
54 
55 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
56 
57 inline Foam::scalar Foam::treeBoundBox::typDim() const
58 {
59  return avgDim();
60 }
61 
62 
64 {
65  return point
66  (
67  (octant & octantBit::rightHalf) ? max().x() : min().x(),
68  (octant & octantBit::topHalf) ? max().y() : min().y(),
69  (octant & octantBit::frontHalf) ? max().z() : min().z()
70  );
71 }
72 
73 
74 // Returns octant in which point resides. Reverse of subBbox.
76 {
77  return subOctant(midpoint(), pt);
78 }
79 
80 
81 // Returns octant in which point resides. Reverse of subBbox.
82 // Precalculated midpoint
84 (
85  const point& mid,
86  const point& pt
87 )
88 {
89  direction octant = 0;
90 
91  if (pt.x() > mid.x())
92  {
93  octant |= octantBit::rightHalf;
94  }
95 
96  if (pt.y() > mid.y())
97  {
98  octant |= octantBit::topHalf;
99  }
100 
101  if (pt.z() > mid.z())
102  {
103  octant |= octantBit::frontHalf;
104  }
105 
106  return octant;
107 }
108 
109 
110 // Returns octant in which point resides. Reverse of subBbox.
111 // Flags point exactly on edge.
113 (
114  const point& pt,
115  bool& onEdge
116 ) const
117 {
118  return subOctant(midpoint(), pt, onEdge);
119 }
120 
121 
122 // Returns octant in which point resides. Reverse of subBbox.
123 // Precalculated midpoint
125 (
126  const point& mid,
127  const point& pt,
128  bool& onEdge
129 )
130 {
131  direction octant = 0;
132  onEdge = false;
133 
134  if (pt.x() > mid.x())
135  {
136  octant |= octantBit::rightHalf;
137  }
138  else if (pt.x() == mid.x())
139  {
140  onEdge = true;
141  }
142 
143  if (pt.y() > mid.y())
144  {
145  octant |= octantBit::topHalf;
146  }
147  else if (pt.y() == mid.y())
148  {
149  onEdge = true;
150  }
151 
152  if (pt.z() > mid.z())
153  {
154  octant |= octantBit::frontHalf;
155  }
156  else if (pt.z() == mid.z())
157  {
158  onEdge = true;
159  }
160 
161  return octant;
162 }
163 
164 
165 // Returns octant in which intersection resides.
166 // Precalculated midpoint. If the point is on the dividing line between
167 // the octants the direction vector determines which octant to use
168 // (i.e. in which octant the point would be if it were moved along dir)
170 (
171  const point& mid,
172  const vector& dir,
173  const point& pt,
174  bool& onEdge
175 )
176 {
177  direction octant = 0;
178  onEdge = false;
179 
180  if (pt.x() > mid.x())
181  {
182  octant |= octantBit::rightHalf;
183  }
184  else if (pt.x() == mid.x())
185  {
186  onEdge = true;
187  if (dir.x() > 0)
188  {
189  octant |= octantBit::rightHalf;
190  }
191  }
192 
193  if (pt.y() > mid.y())
194  {
195  octant |= octantBit::topHalf;
196  }
197  else if (pt.y() == mid.y())
198  {
199  onEdge = true;
200  if (dir.y() > 0)
201  {
202  octant |= octantBit::topHalf;
203  }
204  }
205 
206  if (pt.z() > mid.z())
207  {
208  octant |= octantBit::frontHalf;
209  }
210  else if (pt.z() == mid.z())
211  {
212  onEdge = true;
213  if (dir.z() > 0)
214  {
215  octant |= octantBit::frontHalf;
216  }
217  }
218 
219  return octant;
220 }
221 
222 
223 // Returns reference to octantOrder which defines the
224 // order to do the search.
226 (
227  const point& pt,
228  FixedList<direction,8>& octantOrder
229 ) const
230 {
231  vector dist = midpoint() - pt;
232 
233  direction octant = 0;
234 
235  if (dist.x() < 0)
236  {
237  octant |= octantBit::rightHalf;
238  dist.x() *= -1;
239  }
240 
241  if (dist.y() < 0)
242  {
243  octant |= octantBit::topHalf;
244  dist.y() *= -1;
245  }
246 
247  if (dist.z() < 0)
248  {
249  octant |= octantBit::frontHalf;
250  dist.z() *= -1;
251  }
252 
253  direction min = 0;
254  direction mid = 0;
255  direction max = 0;
256 
257  if (dist.x() < dist.y())
258  {
259  if (dist.y() < dist.z())
260  {
261  min = octantBit::rightHalf;
262  mid = octantBit::topHalf;
263  max = octantBit::frontHalf;
264  }
265  else if (dist.z() < dist.x())
266  {
267  min = octantBit::frontHalf;
268  mid = octantBit::rightHalf;
269  max = octantBit::topHalf;
270  }
271  else
272  {
273  min = octantBit::rightHalf;
274  mid = octantBit::frontHalf;
275  max = octantBit::topHalf;
276  }
277  }
278  else
279  {
280  if (dist.z() < dist.y())
281  {
282  min = octantBit::frontHalf;
283  mid = octantBit::topHalf;
284  max = octantBit::rightHalf;
285  }
286  else if (dist.x() < dist.z())
287  {
288  min = octantBit::topHalf;
289  mid = octantBit::rightHalf;
290  max = octantBit::frontHalf;
291  }
292  else
293  {
294  min = octantBit::topHalf;
295  mid = octantBit::frontHalf;
296  max = octantBit::rightHalf;
297  }
298  }
299 
300  // Primary subOctant
301  octantOrder[0] = octant;
302  // subOctants joined to the primary by faces.
303  octantOrder[1] = octant ^ min;
304  octantOrder[2] = octant ^ mid;
305  octantOrder[3] = octant ^ max;
306  // subOctants joined to the primary by edges.
307  octantOrder[4] = octantOrder[1] ^ mid;
308  octantOrder[5] = octantOrder[1] ^ max;
309  octantOrder[6] = octantOrder[2] ^ max;
310  // subOctants joined to the primary by corner.
311  octantOrder[7] = octantOrder[4] ^ max;
312 }
313 
314 
316 {
317  // Numbers that don't approximate rational fractions with which to make the
318  // box asymmetric. These are between one and two.
319  static const vector a((sqrt(5.0) + 1)/2, sqrt(2.0), (sqrt(13.0) - 1)/2);
320  static const vector b(a.y(), a.z(), a.x());
321 
322  treeBoundBox bb(*this);
323 
324  const scalar delta = s*Foam::mag(bb.span());
325  bb.min() -= Foam::max(delta*a, vector::uniform(rootVSmall));
326  bb.max() += Foam::max(delta*b, vector::uniform(rootVSmall));
327 
328  return bb;
329 }
330 
331 
332 // ************************************************************************* //
scalar y
scalar delta
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:78
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
Definition: Istream.H:60
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:168
const Cmpt & z() const
Definition: VectorI.H:87
const Cmpt & y() const
Definition: VectorI.H:81
const Cmpt & x() const
Definition: VectorI.H:75
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:59
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:54
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:60
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:84
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:90
void searchOrder(const point &pt, FixedList< direction, 8 > &octantOrder) const
Calculates optimal order to look for nearest to point.
treeBoundBox extend(const scalar s) const
Return asymetrically extended bounding box, with guaranteed.
point corner(const direction) const
Corner point given octant.
Definition: treeBoundBoxI.H:63
direction subOctant(const point &pt) const
Returns octant number given point and the calculated midpoint.
Definition: treeBoundBoxI.H:75
scalar typDim() const
Typical dimension length,height,width.
Definition: treeBoundBoxI.H:57
treeBoundBox()
Construct null setting points to zero.
Definition: treeBoundBoxI.H:31
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.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
volScalarField & b
Definition: createFields.H:27
vector point
Point is a vector.
Definition: point.H:41
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
uint8_t direction
Definition: direction.H:45