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) 2011-2013 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 "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)*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  bounds() = boundBox
137  (
138  centre_ - radius_*vector::one,
139  centre_ + radius_*vector::one
140  );
141 }
142 
143 
144 Foam::searchableSphere::searchableSphere
145 (
146  const IOobject& io,
147  const dictionary& dict
148 )
149 :
150  searchableSurface(io),
151  centre_(dict.lookup("centre")),
152  radius_(readScalar(dict.lookup("radius")))
153 {
154  bounds() = boundBox
155  (
156  centre_ - radius_*vector::one,
157  centre_ + radius_*vector::one
158  );
159 }
160 
161 
162 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
163 
165 {}
166 
167 
168 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
169 
171 {
172  return bb.overlaps(centre_, sqr(radius_));
173 }
174 
175 
177 {
178  if (regions_.empty())
179  {
180  regions_.setSize(1);
181  regions_[0] = "region0";
182  }
183  return regions_;
184 }
185 
186 
187 
189 (
190  pointField& centres,
191  scalarField& radiusSqr
192 ) const
193 {
194  centres.setSize(1);
195  centres[0] = centre_;
196 
197  radiusSqr.setSize(1);
198  radiusSqr[0] = Foam::sqr(radius_);
199 
200  // Add a bit to make sure all points are tested inside
201  radiusSqr += Foam::sqr(SMALL);
202 }
203 
204 
205 void Foam::searchableSphere::findNearest
206 (
207  const pointField& samples,
208  const scalarField& nearestDistSqr,
209  List<pointIndexHit>& info
210 ) const
211 {
212  info.setSize(samples.size());
213 
214  forAll(samples, i)
215  {
216  info[i] = findNearest(samples[i], nearestDistSqr[i]);
217  }
218 }
219 
220 
222 (
223  const pointField& start,
224  const pointField& end,
225  List<pointIndexHit>& info
226 ) const
227 {
228  info.setSize(start.size());
229 
231 
232  forAll(start, i)
233  {
234  // Pick nearest intersection. If none intersected take second one.
235  findLineAll(start[i], end[i], info[i], b);
236  if (!info[i].hit() && b.hit())
237  {
238  info[i] = b;
239  }
240  }
241 }
242 
243 
245 (
246  const pointField& start,
247  const pointField& end,
248  List<pointIndexHit>& info
249 ) const
250 {
251  info.setSize(start.size());
252 
254 
255  forAll(start, i)
256  {
257  // Discard far intersection
258  findLineAll(start[i], end[i], info[i], b);
259  if (!info[i].hit() && b.hit())
260  {
261  info[i] = b;
262  }
263  }
264 }
265 
266 
267 void Foam::searchableSphere::findLineAll
268 (
269  const pointField& start,
270  const pointField& end,
271  List<List<pointIndexHit> >& info
272 ) const
273 {
274  info.setSize(start.size());
275 
276  forAll(start, i)
277  {
278  pointIndexHit near, far;
279  findLineAll(start[i], end[i], near, far);
280 
281  if (near.hit())
282  {
283  if (far.hit())
284  {
285  info[i].setSize(2);
286  info[i][0] = near;
287  info[i][1] = far;
288  }
289  else
290  {
291  info[i].setSize(1);
292  info[i][0] = near;
293  }
294  }
295  else
296  {
297  if (far.hit())
298  {
299  info[i].setSize(1);
300  info[i][0] = far;
301  }
302  else
303  {
304  info[i].clear();
305  }
306  }
307  }
308 }
309 
310 
312 (
313  const List<pointIndexHit>& info,
314  labelList& region
315 ) const
316 {
317  region.setSize(info.size());
318  region = 0;
319 }
320 
321 
323 (
324  const List<pointIndexHit>& info,
325  vectorField& normal
326 ) const
327 {
328  normal.setSize(info.size());
329  normal = vector::zero;
330 
331  forAll(info, i)
332  {
333  if (info[i].hit())
334  {
335  normal[i] = info[i].hitPoint() - centre_;
336 
337  normal[i] /= mag(normal[i])+VSMALL;
338  }
339  else
340  {
341  // Set to what?
342  }
343  }
344 }
345 
346 
348 (
349  const pointField& points,
350  List<volumeType>& volType
351 ) const
352 {
353  volType.setSize(points.size());
354  volType = volumeType::INSIDE;
355 
356  forAll(points, pointI)
357  {
358  const point& pt = points[pointI];
359 
360  if (magSqr(pt - centre_) <= sqr(radius_))
361  {
362  volType[pointI] = volumeType::INSIDE;
363  }
364  else
365  {
366  volType[pointI] = volumeType::OUTSIDE;
367  }
368  }
369 }
370 
371 
372 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
vector point
Point is a vector.
Definition: point.H:41
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find first intersection on segment from start to end.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
virtual void getRegion(const List< pointIndexHit > &, labelList &region) const
From a set of points and indices get the region.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Return any intersection on segment from start to end.
bool hit() const
Is there a hit.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
virtual void getVolumeType(const pointField &, List< volumeType > &) const
Determine type (inside/outside/mixed) for point. unknown if.
Namespace for OpenFOAM.
label n
virtual ~searchableSphere()
Destructor.
void clear()
Clear the list, i.e. set size to zero.
Definition: List.C:379
void setSize(const label)
Reset size of List.
Definition: List.C:318
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
virtual bool overlaps(const boundBox &bb) const
Does any part of the surface overlap the supplied bound box?
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
virtual const wordList & regions() const
Names of regions.
Macros for easy insertion into run-time selection tables.
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:55
static const Vector one
Definition: Vector.H:81
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
static const Vector zero
Definition: Vector.H:80
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
dimensionedSymmTensor sqr(const dimensionedVector &dv)
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
Definition: boundBoxI.H:120
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)