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 delta
scalar typDim() const
Typical dimension length,height,width.
Definition: treeBoundBoxI.H:57
direction subOctant(const point &pt) const
Returns octant number given point and the calculated midpoint.
Definition: treeBoundBoxI.H:75
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
A 1D vector of objects of type <T> with a fixed size <Size>.
Definition: FixedList.H:54
point corner(const direction) const
Corner point given octant.
Definition: treeBoundBoxI.H:63
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
dimensionedScalar sqrt(const dimensionedScalar &ds)
const Cmpt & z() const
Definition: VectorI.H:87
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
void searchOrder(const point &pt, FixedList< direction, 8 > &octantOrder) const
Calculates optimal order to look for nearest to point.
const Cmpt & y() const
Definition: VectorI.H:81
scalar avgDim() const
Average length/height/width dimension.
Definition: boundBoxI.H:114
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))
treeBoundBox()
Construct null setting points to zero.
Definition: treeBoundBoxI.H:31
point midpoint() const
The midpoint of the bounding box.
Definition: boundBoxI.H:78
const Cmpt & x() const
Definition: VectorI.H:75
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:60
vector point
Point is a vector.
Definition: point.H:41
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
static Vector< scalar > uniform(const scalar &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:168
dimensioned< scalar > mag(const dimensioned< Type > &)
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:54
treeBoundBox extend(const scalar s) const
Return asymetrically extended bounding box, with guaranteed.