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