FreeFOAM The Cross-Platform CFD Toolkit
treeDataCell.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) 1991-2010 OpenCFD Ltd.
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 "treeDataCell.H"
27 #include "indexedOctree.H"
28 #include <OpenFOAM/primitiveMesh.H>
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
33 
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
37 Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label cellI) const
38 {
39  const cellList& cells = mesh_.cells();
40  const faceList& faces = mesh_.faces();
41  const pointField& points = mesh_.points();
42 
43  treeBoundBox cellBb
44  (
45  vector(GREAT, GREAT, GREAT),
46  vector(-GREAT, -GREAT, -GREAT)
47  );
48 
49  const cell& cFaces = cells[cellI];
50 
51  forAll(cFaces, cFaceI)
52  {
53  const face& f = faces[cFaces[cFaceI]];
54 
55  forAll(f, fp)
56  {
57  const point& p = points[f[fp]];
58 
59  cellBb.min() = min(cellBb.min(), p);
60  cellBb.max() = max(cellBb.max(), p);
61  }
62  }
63  return cellBb;
64 }
65 
66 
67 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
68 
69 // Construct from components
71 (
72  const bool cacheBb,
73  const primitiveMesh& mesh,
74  const labelList& cellLabels
75 )
76 :
77  mesh_(mesh),
78  cellLabels_(cellLabels),
79  cacheBb_(cacheBb)
80 {
81  if (cacheBb_)
82  {
83  bbs_.setSize(cellLabels_.size());
84 
85  forAll(cellLabels_, i)
86  {
87  bbs_[i] = calcCellBb(cellLabels_[i]);
88  }
89  }
90 }
91 
92 
94 (
95  const bool cacheBb,
96  const primitiveMesh& mesh
97 )
98 :
99  mesh_(mesh),
100  cellLabels_(identity(mesh_.nCells())),
101  cacheBb_(cacheBb)
102 {
103  if (cacheBb_)
104  {
105  bbs_.setSize(cellLabels_.size());
106 
107  forAll(cellLabels_, i)
108  {
109  bbs_[i] = calcCellBb(cellLabels_[i]);
110  }
111  }
112 }
113 
114 
115 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
116 
118 {
119  pointField cc(cellLabels_.size());
120 
121  forAll(cellLabels_, i)
122  {
123  cc[i] = mesh_.cellCentres()[cellLabels_[i]];
124  }
125 
126  return cc;
127 }
128 
129 
130 // Check if any point on shape is inside cubeBb.
132 (
133  const label index,
134  const treeBoundBox& cubeBb
135 ) const
136 {
137  if (cacheBb_)
138  {
139  return cubeBb.overlaps(bbs_[index]);
140  }
141  else
142  {
143  return cubeBb.overlaps(calcCellBb(cellLabels_[index]));
144  }
145 }
146 
147 
148 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
149 // nearestPoint.
151 (
152  const labelList& indices,
153  const point& sample,
154 
155  scalar& nearestDistSqr,
156  label& minIndex,
157  point& nearestPoint
158 ) const
159 {
160  forAll(indices, i)
161  {
162  label index = indices[i];
163  label cellI = cellLabels_[index];
164  scalar distSqr = magSqr(sample - mesh_.cellCentres()[cellI]);
165 
166  if (distSqr < nearestDistSqr)
167  {
168  nearestDistSqr = distSqr;
169  minIndex = index;
170  nearestPoint = mesh_.cellCentres()[cellI];
171  }
172  }
173 }
174 
175 
177 (
178  const label index,
179  const point& start,
180  const point& end,
181  point& intersectionPoint
182 ) const
183 {
184  // Do quick rejection test
185  if (cacheBb_)
186  {
187  const treeBoundBox& cellBb = bbs_[index];
188 
189  if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
190  {
191  // start and end in same block outside of cellBb.
192  return false;
193  }
194  }
195  else
196  {
197  const treeBoundBox cellBb = calcCellBb(cellLabels_[index]);
198 
199  if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
200  {
201  // start and end in same block outside of cellBb.
202  return false;
203  }
204  }
205 
206 
207  // Do intersection with all faces of cell
208  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
209 
210  // Disable picking up intersections behind us.
211  scalar oldTol = intersection::setPlanarTol(0.0);
212 
213  const cell& cFaces = mesh_.cells()[cellLabels_[index]];
214 
215  const vector dir(end - start);
216  scalar minDistSqr = magSqr(dir);
217  bool hasMin = false;
218 
219  forAll(cFaces, i)
220  {
221  const face& f = mesh_.faces()[cFaces[i]];
222 
223  pointHit inter = f.ray
224  (
225  start,
226  dir,
227  mesh_.points(),
229  );
230 
231  if (inter.hit() && sqr(inter.distance()) <= minDistSqr)
232  {
233  // Note: no extra test on whether intersection is in front of us
234  // since using half_ray AND zero tolerance. (note that tolerance
235  // is used to look behind us)
236  minDistSqr = sqr(inter.distance());
237  intersectionPoint = inter.hitPoint();
238  hasMin = true;
239  }
240  }
241 
242  // Restore picking tolerance
244 
245  return hasMin;
246 }
247 
248 
249 // ************************ vim: set sw=4 sts=4 et: ************************ //