FreeFOAM The Cross-Platform CFD Toolkit
searchableSphere.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 "searchableSphere.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 
34 defineTypeNameAndDebug(searchableSphere, 0);
35 addToRunTimeSelectionTable(searchableSurface, searchableSphere, dict);
36 
37 }
38 
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 Foam::pointIndexHit Foam::searchableSphere::findNearest
43 (
44  const point& sample,
45  const scalar nearestDistSqr
46 ) const
47 {
48  pointIndexHit info(false, sample, -1);
49 
50  const vector n(sample-centre_);
51  scalar magN = mag(n);
52 
53  if (nearestDistSqr > sqr(magN-radius_))
54  {
55  if (magN < ROOTVSMALL)
56  {
57  info.rawPoint() = centre_ + vector(1,0,0)/magN*radius_;
58  }
59  else
60  {
61  info.rawPoint() = centre_ + n/magN*radius_;
62  }
63  info.setHit();
64  info.setIndex(0);
65  }
66 
67  return info;
68 }
69 
70 
71 // From Graphics Gems - intersection of sphere with ray
72 void Foam::searchableSphere::findLineAll
73 (
74  const point& start,
75  const point& end,
76  pointIndexHit& near,
77  pointIndexHit& far
78 ) const
79 {
80  near.setMiss();
81  far.setMiss();
82 
83  vector dir(end-start);
84  scalar magSqrDir = magSqr(dir);
85 
86  if (magSqrDir > ROOTVSMALL)
87  {
88  const vector toCentre(centre_-start);
89  scalar magSqrToCentre = magSqr(toCentre);
90 
91  dir /= Foam::sqrt(magSqrDir);
92 
93  scalar v = (toCentre & dir);
94 
95  scalar disc = sqr(radius_) - (magSqrToCentre - sqr(v));
96 
97  if (disc >= 0)
98  {
99  scalar d = Foam::sqrt(disc);
100 
101  scalar nearParam = v-d;
102 
103  if (nearParam >= 0 && sqr(nearParam) <= magSqrDir)
104  {
105  near.setHit();
106  near.setPoint(start + nearParam*dir);
107  near.setIndex(0);
108  }
109 
110  scalar farParam = v+d;
111 
112  if (farParam >= 0 && sqr(farParam) <= magSqrDir)
113  {
114  far.setHit();
115  far.setPoint(start + farParam*dir);
116  far.setIndex(0);
117  }
118  }
119  }
120 }
121 
122 
123 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
124 
125 Foam::searchableSphere::searchableSphere
126 (
127  const IOobject& io,
128  const point& centre,
129  const scalar radius
130 )
131 :
132  searchableSurface(io),
133  centre_(centre),
134  radius_(radius)
135 {}
136 
137 
138 Foam::searchableSphere::searchableSphere
139 (
140  const IOobject& io,
141  const dictionary& dict
142 )
143 :
144  searchableSurface(io),
145  centre_(dict.lookup("centre")),
146  radius_(readScalar(dict.lookup("radius")))
147 {}
148 
149 
150 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
151 
153 {}
154 
155 
156 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
157 
159 {
160  if (regions_.empty())
161  {
162  regions_.setSize(1);
163  regions_[0] = "region0";
164  }
165  return regions_;
166 }
167 
168 
169 
170 void Foam::searchableSphere::findNearest
171 (
172  const pointField& samples,
173  const scalarField& nearestDistSqr,
174  List<pointIndexHit>& info
175 ) const
176 {
177  info.setSize(samples.size());
178 
179  forAll(samples, i)
180  {
181  info[i] = findNearest(samples[i], nearestDistSqr[i]);
182  }
183 }
184 
185 
187 (
188  const pointField& start,
189  const pointField& end,
190  List<pointIndexHit>& info
191 ) const
192 {
193  info.setSize(start.size());
194 
196 
197  forAll(start, i)
198  {
199  // Pick nearest intersection. If none intersected take second one.
200  findLineAll(start[i], end[i], info[i], b);
201  if (!info[i].hit() && b.hit())
202  {
203  info[i] = b;
204  }
205  }
206 }
207 
208 
210 (
211  const pointField& start,
212  const pointField& end,
213  List<pointIndexHit>& info
214 ) const
215 {
216  info.setSize(start.size());
217 
219 
220  forAll(start, i)
221  {
222  // Discard far intersection
223  findLineAll(start[i], end[i], info[i], b);
224  if (!info[i].hit() && b.hit())
225  {
226  info[i] = b;
227  }
228  }
229 }
230 
231 
232 void Foam::searchableSphere::findLineAll
233 (
234  const pointField& start,
235  const pointField& end,
236  List<List<pointIndexHit> >& info
237 ) const
238 {
239  info.setSize(start.size());
240 
241  forAll(start, i)
242  {
243  pointIndexHit near, far;
244  findLineAll(start[i], end[i], near, far);
245 
246  if (near.hit())
247  {
248  if (far.hit())
249  {
250  info[i].setSize(2);
251  info[i][0] = near;
252  info[i][1] = far;
253  }
254  else
255  {
256  info[i].setSize(1);
257  info[i][0] = near;
258  }
259  }
260  else
261  {
262  if (far.hit())
263  {
264  info[i].setSize(1);
265  info[i][0] = far;
266  }
267  else
268  {
269  info[i].clear();
270  }
271  }
272  }
273 }
274 
275 
277 (
278  const List<pointIndexHit>& info,
279  labelList& region
280 ) const
281 {
282  region.setSize(info.size());
283  region = 0;
284 }
285 
286 
288 (
289  const List<pointIndexHit>& info,
291 ) const
292 {
293  normal.setSize(info.size());
294  normal = vector::zero;
295 
296  forAll(info, i)
297  {
298  if (info[i].hit())
299  {
300  normal[i] = info[i].hitPoint() - centre_;
301 
302  normal[i] /= mag(normal[i])+VSMALL;
303  }
304  else
305  {
306  // Set to what?
307  }
308  }
309 }
310 
311 
313 (
314  const pointField& points,
315  List<volumeType>& volType
316 ) const
317 {
318  volType.setSize(points.size());
319  volType = INSIDE;
320 
321  forAll(points, pointI)
322  {
323  const point& pt = points[pointI];
324 
325  if (magSqr(pt - centre_) <= sqr(radius_))
326  {
327  volType[pointI] = INSIDE;
328  }
329  else
330  {
331  volType[pointI] = OUTSIDE;
332  }
333  }
334 }
335 
336 
337 // ************************ vim: set sw=4 sts=4 et: ************************ //