searchableBox.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-2024 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 "searchableBox.H"
28 #include "SortableList.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
36 }
37 
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 void Foam::searchableBox::projectOntoCoordPlane
42 (
43  const direction dir,
44  const point& planePt,
45  pointIndexHit& info
46 ) const
47 {
48  // Set point
49  info.rawPoint()[dir] = planePt[dir];
50  // Set face
51  if (planePt[dir] == min()[dir])
52  {
53  info.setIndex(dir*2);
54  }
55  else if (planePt[dir] == max()[dir])
56  {
57  info.setIndex(dir*2+1);
58  }
59  else
60  {
62  << "Point on plane " << planePt
63  << " is not on coordinate " << min()[dir]
64  << " nor " << max()[dir] << abort(FatalError);
65  }
66 }
67 
68 
69 // Returns miss or hit with face (0..5) and region(always 0)
70 Foam::pointIndexHit Foam::searchableBox::findNearest
71 (
72  const point& bbMid,
73  const point& sample,
74  const scalar nearestDistSqr
75 ) const
76 {
77  // Point can be inside or outside. For every component direction can be
78  // left of min, right of max or in between.
79  // - outside points: project first one x plane (either min().x()
80  // or max().x()), then onto y plane and finally z. You should be left
81  // with intersection point
82  // - inside point: find nearest side (compare to mid point). Project onto
83  // that.
84 
85  // The face is set to the last projected face.
86 
87 
88  // Outside point projected onto cube. Assume faces 0..5.
89  pointIndexHit info(true, sample, -1);
90  bool outside = false;
91 
92  // (for internal points) per direction what nearest cube side is
93  point near;
94 
95  for (direction dir = 0; dir < vector::nComponents; dir++)
96  {
97  if (info.rawPoint()[dir] < min()[dir])
98  {
99  projectOntoCoordPlane(dir, min(), info);
100  outside = true;
101  }
102  else if (info.rawPoint()[dir] > max()[dir])
103  {
104  projectOntoCoordPlane(dir, max(), info);
105  outside = true;
106  }
107  else if (info.rawPoint()[dir] > bbMid[dir])
108  {
109  near[dir] = max()[dir];
110  }
111  else
112  {
113  near[dir] = min()[dir];
114  }
115  }
116 
117 
118  // For outside points the info will be correct now. Handle inside points
119  // using the three near distances. Project onto the nearest plane.
120  if (!outside)
121  {
122  vector dist(cmptMag(info.rawPoint() - near));
123 
124  if (dist.x() < dist.y())
125  {
126  if (dist.x() < dist.z())
127  {
128  // Project onto x plane
129  projectOntoCoordPlane(vector::X, near, info);
130  }
131  else
132  {
133  projectOntoCoordPlane(vector::Z, near, info);
134  }
135  }
136  else
137  {
138  if (dist.y() < dist.z())
139  {
140  projectOntoCoordPlane(vector::Y, near, info);
141  }
142  else
143  {
144  projectOntoCoordPlane(vector::Z, near, info);
145  }
146  }
147  }
148 
149 
150  // Check if outside. Optimisation: could do some checks on distance already
151  // on components above
152  if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
153  {
154  info.setMiss();
155  info.setIndex(-1);
156  }
157 
158  return info;
159 }
160 
161 
162 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
163 
165 (
166  const IOobject& io,
167  const treeBoundBox& bb
168 )
169 :
170  searchableSurface(io),
171  treeBoundBox(bb)
172 {
173  if (!contains(midpoint()))
174  {
176  << "Illegal bounding box specification : "
177  << static_cast<const treeBoundBox>(*this) << exit(FatalError);
178  }
179 
180  bounds() = static_cast<boundBox>(*this);
181 }
182 
183 
185 (
186  const IOobject& io,
187  const dictionary& dict
188 )
189 :
190  searchableSurface(io),
192 {
193  if (!contains(midpoint()))
194  {
196  << "Illegal bounding box specification : "
197  << static_cast<const treeBoundBox>(*this) << exit(FatalError);
198  }
199 
200  bounds() = static_cast<boundBox>(*this);
201 }
202 
203 
204 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
205 
207 {}
208 
209 
210 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
211 
213 {
214  if (regions_.empty())
215  {
216  regions_.setSize(1);
217  regions_[0] = "region0";
218  }
219  return regions_;
220 }
221 
222 
224 {
226  pointField& ctrs = tCtrs.ref();
227 
228  const pointField pts(treeBoundBox::points());
229  const faceList& fcs = treeBoundBox::faces;
230 
231  forAll(fcs, i)
232  {
233  ctrs[i] = fcs[i].centre(pts);
234  }
235 
236  return tCtrs;
237 }
238 
239 
241 (
242  pointField& centres,
243  scalarField& radiusSqr
244 ) const
245 {
246  centres.setSize(size());
247  radiusSqr.setSize(size());
248  radiusSqr = 0.0;
249 
250  const pointField pts(treeBoundBox::points());
251  const faceList& fcs = treeBoundBox::faces;
252 
253  forAll(fcs, i)
254  {
255  const face& f = fcs[i];
256 
257  centres[i] = f.centre(pts);
258  forAll(f, fp)
259  {
260  const point& pt = pts[f[fp]];
261 
262  radiusSqr[i] = Foam::max
263  (
264  radiusSqr[i],
265  Foam::magSqr(pt-centres[i])
266  );
267  }
268  }
269 
270  // Add a bit to make sure all points are tested inside
271  radiusSqr += Foam::sqr(small);
272 }
273 
274 
276 {
277  return treeBoundBox::points();
278 }
279 
280 
281 Foam::pointIndexHit Foam::searchableBox::findNearest
282 (
283  const point& sample,
284  const scalar nearestDistSqr
285 ) const
286 {
287  return findNearest(midpoint(), sample, nearestDistSqr);
288 }
289 
290 
292 (
293  const point& sample,
294  const scalar nearestDistSqr
295 ) const
296 {
297  const point bbMid(midpoint());
298 
299  // Outside point projected onto cube. Assume faces 0..5.
300  pointIndexHit info(true, sample, -1);
301  bool outside = false;
302 
303  // (for internal points) per direction what nearest cube side is
304  point near;
305 
306  for (direction dir = 0; dir < vector::nComponents; dir++)
307  {
308  if (info.rawPoint()[dir] < min()[dir])
309  {
310  projectOntoCoordPlane(dir, min(), info);
311  outside = true;
312  }
313  else if (info.rawPoint()[dir] > max()[dir])
314  {
315  projectOntoCoordPlane(dir, max(), info);
316  outside = true;
317  }
318  else if (info.rawPoint()[dir] > bbMid[dir])
319  {
320  near[dir] = max()[dir];
321  }
322  else
323  {
324  near[dir] = min()[dir];
325  }
326  }
327 
328 
329  // For outside points the info will be correct now. Handle inside points
330  // using the three near distances. Project onto the nearest two planes.
331  if (!outside)
332  {
333  // Get the per-component distance to nearest wall
334  vector dist(cmptMag(info.rawPoint() - near));
335 
336  SortableList<scalar> sortedDist(3);
337  sortedDist[0] = dist[0];
338  sortedDist[1] = dist[1];
339  sortedDist[2] = dist[2];
340  sortedDist.sort();
341 
342  // Project onto nearest
343  projectOntoCoordPlane(sortedDist.indices()[0], near, info);
344  // Project onto second nearest
345  projectOntoCoordPlane(sortedDist.indices()[1], near, info);
346  }
347 
348 
349  // Check if outside. Optimisation: could do some checks on distance already
350  // on components above
351  if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
352  {
353  info.setMiss();
354  info.setIndex(-1);
355  }
356 
357  return info;
358 }
359 
360 
361 Foam::pointIndexHit Foam::searchableBox::findNearest
362 (
363  const linePointRef& ln,
364  treeBoundBox& tightest,
365  point& linePoint
366 ) const
367 {
369  return pointIndexHit();
370 }
371 
372 
374 (
375  const point& start,
376  const point& end
377 ) const
378 {
379  pointIndexHit info(false, start, -1);
380 
381  bool foundInter;
382 
383  if (posBits(start) == 0)
384  {
385  if (posBits(end) == 0)
386  {
387  // Both start and end inside.
388  foundInter = false;
389  }
390  else
391  {
392  // end is outside. Clip to bounding box.
393  foundInter = intersects(end, start, info.rawPoint());
394  }
395  }
396  else
397  {
398  // start is outside. Clip to bounding box.
399  foundInter = intersects(start, end, info.rawPoint());
400  }
401 
402 
403  // Classify point
404  if (foundInter)
405  {
406  info.setHit();
407 
408  for (direction dir = 0; dir < vector::nComponents; dir++)
409  {
410  if (info.rawPoint()[dir] == min()[dir])
411  {
412  info.setIndex(2*dir);
413  break;
414  }
415  else if (info.rawPoint()[dir] == max()[dir])
416  {
417  info.setIndex(2*dir+1);
418  break;
419  }
420  }
421 
422  if (info.index() == -1)
423  {
425  << "point " << info.rawPoint()
426  << " on segment " << start << end
427  << " should be on face of " << *this
428  << " but it isn't." << abort(FatalError);
429  }
430  }
431 
432  return info;
433 }
434 
435 
437 (
438  const point& start,
439  const point& end
440 ) const
441 {
442  return findLine(start, end);
443 }
444 
445 
446 void Foam::searchableBox::findNearest
447 (
448  const pointField& samples,
449  const scalarField& nearestDistSqr,
450  List<pointIndexHit>& info
451 ) const
452 {
453  info.setSize(samples.size());
454 
455  const point bbMid(midpoint());
456 
457  forAll(samples, i)
458  {
459  info[i] = findNearest(bbMid, samples[i], nearestDistSqr[i]);
460  }
461 }
462 
463 
465 (
466  const pointField& start,
467  const pointField& end,
468  List<pointIndexHit>& info
469 ) const
470 {
471  info.setSize(start.size());
472 
473  forAll(start, i)
474  {
475  info[i] = findLine(start[i], end[i]);
476  }
477 }
478 
479 
481 (
482  const pointField& start,
483  const pointField& end,
484  List<pointIndexHit>& info
485 ) const
486 {
487  info.setSize(start.size());
488 
489  forAll(start, i)
490  {
491  info[i] = findLineAny(start[i], end[i]);
492  }
493 }
494 
495 
497 (
498  const pointField& start,
499  const pointField& end,
501 ) const
502 {
503  info.setSize(start.size());
504 
505  // Work array
507 
508  // Tolerances:
509  // To find all intersections we add a small vector to the last intersection
510  // This is chosen such that
511  // - it is significant (small is smallest representative relative tolerance;
512  // we need something bigger since we're doing calculations)
513  // - if the start-end vector is zero we still progress
514  const vectorField dirVec(end-start);
515  const scalarField magSqrDirVec(magSqr(dirVec));
516  const vectorField smallVec
517  (
518  rootSmall*dirVec
519  + vector(rootVSmall,rootVSmall,rootVSmall)
520  );
521 
522  forAll(start, pointi)
523  {
524  // See if any intersection between pt and end
525  pointIndexHit inter = findLine(start[pointi], end[pointi]);
526 
527  if (inter.hit())
528  {
529  hits.clear();
530  hits.append(inter);
531 
532  point pt = inter.hitPoint() + smallVec[pointi];
533 
534  while (((pt-start[pointi])&dirVec[pointi]) <= magSqrDirVec[pointi])
535  {
536  // See if any intersection between pt and end
537  pointIndexHit inter = findLine(pt, end[pointi]);
538 
539  // Check for not hit or hit same face as before (can happen
540  // if vector along surface of face)
541  if
542  (
543  !inter.hit()
544  || (inter.index() == hits.last().index())
545  )
546  {
547  break;
548  }
549  hits.append(inter);
550 
551  pt = inter.hitPoint() + smallVec[pointi];
552  }
553 
554  info[pointi].transfer(hits);
555  }
556  else
557  {
558  info[pointi].clear();
559  }
560  }
561 }
562 
563 
565 (
566  const List<pointIndexHit>& info,
567  labelList& region
568 ) const
569 {
570  region.setSize(info.size());
571  region = 0;
572 }
573 
574 
576 (
577  const List<pointIndexHit>& info,
578  vectorField& normal
579 ) const
580 {
581  normal.setSize(info.size());
582  normal = Zero;
583 
584  forAll(info, i)
585  {
586  if (info[i].hit())
587  {
588  normal[i] = treeBoundBox::faceNormals[info[i].index()];
589  }
590  else
591  {
592  // Set to what?
593  }
594  }
595 }
596 
597 
599 (
600  const pointField& points,
601  List<volumeType>& volType
602 ) const
603 {
604  volType.setSize(points.size());
605  volType = volumeType::inside;
606 
607  forAll(points, pointi)
608  {
609  const point& pt = points[pointi];
610 
611  for (direction dir = 0; dir < vector::nComponents; dir++)
612  {
613  if (pt[dir] < min()[dir] || pt[dir] > max()[dir])
614  {
615  volType[pointi] = volumeType::outside;
616  break;
617  }
618  }
619  }
620 }
621 
622 
623 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:78
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
InfoProxy< IOobject > info() const
Return info proxy.
Definition: IOobject.H:486
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
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
const Point & hitPoint() const
Return hit point.
const Point & rawPoint() const
Return point with no checking.
label index() const
Return index.
bool hit() const
Is there a hit.
void setIndex(const label index)
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: SortableList.H:55
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:96
void sort()
(stable) sort the list (if changed after construction time)
Definition: SortableList.C:112
T & last()
Return the last element of the list.
Definition: UListI.H:128
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:105
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
point midpoint() const
The midpoint of the bounding box.
Definition: boundBoxI.H:84
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:76
A line primitive.
Definition: line.H:71
Surface geometry with a rectangular box shape, aligned with the coordinate axes, which can be used wi...
Definition: searchableBox.H:82
virtual ~searchableBox()
Destructor.
virtual void getVolumeType(const pointField &, List< volumeType > &) const
Determine type (inside/outside/mixed) for point. unknown if.
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
searchableBox(const IOobject &io, const treeBoundBox &bb)
Construct from components.
pointIndexHit findNearestOnEdge(const point &sample, const scalar nearestDistSqr) const
Calculate nearest point on edge.
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
pointIndexHit findLine(const point &start, const point &end) const
Find nearest intersection of line between start and end.
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.
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
virtual tmp< pointField > points() const
Get the points that define the surface.
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
const boundBox & bounds() const
Return const reference to boundBox.
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:181
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:90
bool contains(const point &) const
Contains point or other bounding box?
Definition: boundBoxI.H:176
static const FixedList< vector, 6 > faceNormals
Per face the unit normal.
Definition: treeBoundBox.H:181
static const faceList faces
Face to point addressing.
Definition: treeBoundBox.H:175
tmp< pointField > points() const
Vertex coordinates. In octant coding.
Definition: treeBoundBox.C:159
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:381
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
const pointField & points
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
addToRunTimeSelectionTable(polyPatch, mergedCyclicPolyPatch, word)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
errorManip< error > abort(error &err)
Definition: errorManip.H:131
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
vector point
Point is a vector.
Definition: point.H:41
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
dimensionSet cmptMag(const dimensionSet &)
Definition: dimensionSet.C:296
error FatalError
dimensioned< scalar > magSqr(const dimensioned< Type > &)
bool ln(const fileName &src, const fileName &dst)
Create a softlink. dst should not exist. Returns true if successful.
Definition: POSIX.C:908
uint8_t direction
Definition: direction.H:45
labelList f(nPoints)
dictionary dict
scalarField samples(nIntervals, 0)