searchablePlate.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-2014 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 "searchablePlate.H"
28 #include "SortableList.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 defineTypeNameAndDebug(searchablePlate, 0);
36 addToRunTimeSelectionTable(searchableSurface, searchablePlate, dict);
37 
38 }
39 
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 Foam::direction Foam::searchablePlate::calcNormal(const point& span)
44 {
45  direction normalDir = 3;
46 
47  for (direction dir = 0; dir < vector::nComponents; dir++)
48  {
49  if (span[dir] < 0)
50  {
51  FatalErrorIn("searchablePlate::calcNormal()")
52  << "Span should have two positive and one zero entry. Now:"
53  << span << exit(FatalError);
54  }
55  else if (span[dir] < VSMALL)
56  {
57  if (normalDir == 3)
58  {
59  normalDir = dir;
60  }
61  else
62  {
63  // Multiple zero entries. Flag and exit.
64  normalDir = 3;
65  break;
66  }
67  }
68  }
69 
70  if (normalDir == 3)
71  {
72  FatalErrorIn("searchablePlate::calcNormal()")
73  << "Span should have two positive and one zero entry. Now:"
74  << span << exit(FatalError);
75  }
76 
77  return normalDir;
78 }
79 
80 
81 // Returns miss or hit with face (always 0)
82 Foam::pointIndexHit Foam::searchablePlate::findNearest
83 (
84  const point& sample,
85  const scalar nearestDistSqr
86 ) const
87 {
88  // For every component direction can be
89  // left of min, right of max or inbetween.
90  // - outside points: project first one x plane (either min().x()
91  // or max().x()), then onto y plane and finally z. You should be left
92  // with intersection point
93  // - inside point: find nearest side (compare to mid point). Project onto
94  // that.
95 
96  // Project point on plane.
97  pointIndexHit info(true, sample, 0);
98  info.rawPoint()[normalDir_] = origin_[normalDir_];
99 
100  // Clip to edges if outside
101  for (direction dir = 0; dir < vector::nComponents; dir++)
102  {
103  if (dir != normalDir_)
104  {
105  if (info.rawPoint()[dir] < origin_[dir])
106  {
107  info.rawPoint()[dir] = origin_[dir];
108  }
109  else if (info.rawPoint()[dir] > origin_[dir]+span_[dir])
110  {
111  info.rawPoint()[dir] = origin_[dir]+span_[dir];
112  }
113  }
114  }
115 
116  // Check if outside. Optimisation: could do some checks on distance already
117  // on components above
118  if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
119  {
120  info.setMiss();
121  info.setIndex(-1);
122  }
123 
124  return info;
125 }
126 
127 
128 Foam::pointIndexHit Foam::searchablePlate::findLine
129 (
130  const point& start,
131  const point& end
132 ) const
133 {
134  pointIndexHit info
135  (
136  true,
137  vector::zero,
138  0
139  );
140 
141  const vector dir(end-start);
142 
143  if (mag(dir[normalDir_]) < VSMALL)
144  {
145  info.setMiss();
146  info.setIndex(-1);
147  }
148  else
149  {
150  scalar t = (origin_[normalDir_]-start[normalDir_]) / dir[normalDir_];
151 
152  if (t < 0 || t > 1)
153  {
154  info.setMiss();
155  info.setIndex(-1);
156  }
157  else
158  {
159  info.rawPoint() = start+t*dir;
160  info.rawPoint()[normalDir_] = origin_[normalDir_];
161 
162  // Clip to edges
163  for (direction dir = 0; dir < vector::nComponents; dir++)
164  {
165  if (dir != normalDir_)
166  {
167  if (info.rawPoint()[dir] < origin_[dir])
168  {
169  info.setMiss();
170  info.setIndex(-1);
171  break;
172  }
173  else if (info.rawPoint()[dir] > origin_[dir]+span_[dir])
174  {
175  info.setMiss();
176  info.setIndex(-1);
177  break;
178  }
179  }
180  }
181  }
182  }
183 
184  // Debug
185  if (info.hit())
186  {
187  treeBoundBox bb(origin_, origin_+span_);
188  bb.min()[normalDir_] -= 1e-6;
189  bb.max()[normalDir_] += 1e-6;
190 
191  if (!bb.contains(info.hitPoint()))
192  {
193  FatalErrorIn("searchablePlate::findLine(..)")
194  << "bb:" << bb << endl
195  << "origin_:" << origin_ << endl
196  << "span_:" << span_ << endl
197  << "normalDir_:" << normalDir_ << endl
198  << "hitPoint:" << info.hitPoint()
199  << abort(FatalError);
200  }
201  }
202 
203  return info;
204 }
205 
206 
207 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
208 
209 Foam::searchablePlate::searchablePlate
210 (
211  const IOobject& io,
212  const point& origin,
213  const vector& span
214 )
215 :
216  searchableSurface(io),
217  origin_(origin),
218  span_(span),
219  normalDir_(calcNormal(span_))
220 {
221  if (debug)
222  {
223  Info<< "searchablePlate::searchablePlate :"
224  << " origin:" << origin_
225  << " origin+span:" << origin_+span_
226  << " normal:" << vector::componentNames[normalDir_]
227  << endl;
228  }
229 
230  bounds() = boundBox(origin_, origin_ + span_);
231 }
232 
233 
234 Foam::searchablePlate::searchablePlate
235 (
236  const IOobject& io,
237  const dictionary& dict
238 )
239 :
240  searchableSurface(io),
241  origin_(dict.lookup("origin")),
242  span_(dict.lookup("span")),
243  normalDir_(calcNormal(span_))
244 {
245  if (debug)
246  {
247  Info<< "searchablePlate::searchablePlate :"
248  << " origin:" << origin_
249  << " origin+span:" << origin_+span_
250  << " normal:" << vector::componentNames[normalDir_]
251  << endl;
252  }
253 
254  bounds() = boundBox(origin_, origin_ + span_);
255 }
256 
257 
258 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
259 
261 {}
262 
263 
264 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
265 
267 {
268  if (regions_.empty())
269  {
270  regions_.setSize(1);
271  regions_[0] = "region0";
272  }
273  return regions_;
274 }
275 
276 
278 {
279  return tmp<pointField>(new pointField(1, origin_ + 0.5*span_));
280 }
281 
282 
284 (
285  pointField& centres,
286  scalarField& radiusSqr
287 ) const
288 {
289  centres.setSize(1);
290  centres[0] = origin_ + 0.5*span_;
291 
292  radiusSqr.setSize(1);
293  radiusSqr[0] = Foam::magSqr(0.5*span_);
294 
295  // Add a bit to make sure all points are tested inside
296  radiusSqr += Foam::sqr(SMALL);
297 }
298 
299 
301 {
302  tmp<pointField> tPts(new pointField(4));
303  pointField& pts = tPts();
304 
305  pts[0] = origin_;
306  pts[2] = origin_ + span_;
307 
308  if (span_.x() < span_.y() && span_.x() < span_.z())
309  {
310  pts[1] = origin_ + point(0, span_.y(), 0);
311  pts[3] = origin_ + point(0, 0, span_.z());
312  }
313  else if (span_.y() < span_.z())
314  {
315  pts[1] = origin_ + point(span_.x(), 0, 0);
316  pts[3] = origin_ + point(0, 0, span_.z());
317  }
318  else
319  {
320  pts[1] = origin_ + point(span_.x(), 0, 0);
321  pts[3] = origin_ + point(0, span_.y(), 0);
322  }
323 
324  return tPts;
325 }
326 
327 
329 {
330  return
331  (
332  (origin_.x() + span_.x()) >= bb.min().x()
333  && origin_.x() <= bb.max().x()
334  && (origin_.y() + span_.y()) >= bb.min().y()
335  && origin_.y() <= bb.max().y()
336  && (origin_.z() + span_.z()) >= bb.min().z()
337  && origin_.z() <= bb.max().z()
338  );
339 }
340 
341 
342 void Foam::searchablePlate::findNearest
343 (
344  const pointField& samples,
345  const scalarField& nearestDistSqr,
346  List<pointIndexHit>& info
347 ) const
348 {
349  info.setSize(samples.size());
350 
351  forAll(samples, i)
352  {
353  info[i] = findNearest(samples[i], nearestDistSqr[i]);
354  }
355 }
356 
357 
358 void Foam::searchablePlate::findLine
359 (
360  const pointField& start,
361  const pointField& end,
362  List<pointIndexHit>& info
363 ) const
364 {
365  info.setSize(start.size());
366 
367  forAll(start, i)
368  {
369  info[i] = findLine(start[i], end[i]);
370  }
371 }
372 
373 
375 (
376  const pointField& start,
377  const pointField& end,
378  List<pointIndexHit>& info
379 ) const
380 {
381  findLine(start, end, info);
382 }
383 
384 
386 (
387  const pointField& start,
388  const pointField& end,
389  List<List<pointIndexHit> >& info
390 ) const
391 {
392  List<pointIndexHit> nearestInfo;
393  findLine(start, end, nearestInfo);
394 
395  info.setSize(start.size());
396  forAll(info, pointI)
397  {
398  if (nearestInfo[pointI].hit())
399  {
400  info[pointI].setSize(1);
401  info[pointI][0] = nearestInfo[pointI];
402  }
403  else
404  {
405  info[pointI].clear();
406  }
407  }
408 }
409 
410 
412 (
413  const List<pointIndexHit>& info,
414  labelList& region
415 ) const
416 {
417  region.setSize(info.size());
418  region = 0;
419 }
420 
421 
423 (
424  const List<pointIndexHit>& info,
425  vectorField& normal
426 ) const
427 {
428  normal.setSize(info.size());
429  normal = vector::zero;
430  forAll(normal, i)
431  {
432  normal[i][normalDir_] = 1.0;
433  }
434 }
435 
436 
438 (
439  const pointField& points,
440  List<volumeType>& volType
441 ) const
442 {
444  (
445  "searchableCollection::getVolumeType(const pointField&"
446  ", List<volumeType>&) const"
447  ) << "Volume type not supported for plate."
448  << exit(FatalError);
449 }
450 
451 
452 // ************************************************************************* //
unsigned char direction
Definition: direction.H:43
virtual bool overlaps(const boundBox &bb) const
Does any part of the surface overlap the supplied bound box?
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
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
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
dimensioned< scalar > mag(const dimensioned< Type > &)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static const char * componentNames[]
Definition: Vector.H:79
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
messageStream Info
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
virtual const wordList & regions() const
Names of regions.
Number of components in this vector space.
Definition: VectorSpace.H:88
void clear()
Clear the list, i.e. set size to zero.
Definition: List.C:379
const double e
Elementary charge.
Definition: doubleFloat.H:78
void setSize(const label)
Reset size of List.
Definition: List.C:318
const Cmpt & y() const
Definition: VectorI.H:71
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
virtual tmp< pointField > points() const
Get the points that define the surface.
#define forAll(list, i)
Definition: UList.H:421
virtual ~searchablePlate()
Destructor.
const Cmpt & x() const
Definition: VectorI.H:65
Macros for easy insertion into run-time selection tables.
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
const Cmpt & z() const
Definition: VectorI.H:77
errorManip< error > abort(error &err)
Definition: errorManip.H:131
virtual void findLineAll(const pointField &start, const pointField &end, List< List< pointIndexHit > > &) const
Get all intersections in order from start to end.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:55
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
error FatalError
static const Vector zero
Definition: Vector.H:80
virtual tmp< pointField > coordinates() const
Get representative set of element coordinates.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
virtual void getRegion(const List< pointIndexHit > &, labelList &region) const
From a set of points and indices get the region.
virtual void getVolumeType(const pointField &, List< volumeType > &) const
Determine type (inside/outside/mixed) for point. unknown if.
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Return any intersection on segment from start to end.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
A class for managing temporary objects.
Definition: PtrList.H:118
defineTypeNameAndDebug(combustionModel, 0)