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