searchableSurfacesQueries.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-2021 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 "ListOps.H"
28 #include "OFstream.H"
29 #include "meshTools.H"
30 #include "DynamicField.H"
31 #include "pointConstraint.H"
32 #include "plane.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(searchableSurfacesQueries, 0);
39 }
40 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 void Foam::searchableSurfacesQueries::mergeHits
45 (
46  const point& start,
47 
48  const label testI, // index of surface
49  const List<pointIndexHit>& surfHits, // hits on surface
50 
51  labelList& allSurfaces,
52  List<pointIndexHit>& allInfo,
53  scalarList& allDistSqr
54 )
55 {
56  // Given current set of hits (allSurfaces, allInfo) merge in those coming
57  // from surface surfI.
58 
59  // Precalculate distances
60  scalarList surfDistSqr(surfHits.size());
61  forAll(surfHits, i)
62  {
63  surfDistSqr[i] = magSqr(surfHits[i].hitPoint() - start);
64  }
65 
66  forAll(surfDistSqr, i)
67  {
68  label index = findLower(allDistSqr, surfDistSqr[i]);
69 
70  // Check if equal to lower.
71  if (index >= 0)
72  {
73  // Same. Do not count.
74  // Pout<< "point:" << surfHits[i].hitPoint()
75  // << " considered same as:" << allInfo[index].hitPoint()
76  // << " within tol:" << mergeDist
77  // << endl;
78  }
79  else
80  {
81  // Check if equal to higher
82  label next = index + 1;
83 
84  if (next < allDistSqr.size())
85  {
86  // Pout<< "point:" << surfHits[i].hitPoint()
87  // << " considered same as:" << allInfo[next].hitPoint()
88  // << " within tol:" << mergeDist
89  // << endl;
90  }
91  else
92  {
93  // Insert after index
94  label sz = allSurfaces.size();
95  allSurfaces.setSize(sz+1);
96  allInfo.setSize(allSurfaces.size());
97  allDistSqr.setSize(allSurfaces.size());
98  // Make space.
99  for (label j = sz-1; j > index; --j)
100  {
101  allSurfaces[j+1] = allSurfaces[j];
102  allInfo[j+1] = allInfo[j];
103  allDistSqr[j+1] = allDistSqr[j];
104  }
105  // Insert new value
106  allSurfaces[index+1] = testI;
107  allInfo[index+1] = surfHits[i];
108  allDistSqr[index+1] = surfDistSqr[i];
109  }
110  }
111  }
112 }
113 
114 
115 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
116 
117 // Find any intersection
119 (
120  const PtrList<searchableSurface>& allSurfaces,
121  const labelList& surfacesToTest,
122  const pointField& start,
123  const pointField& end,
124  labelList& hitSurfaces,
125  List<pointIndexHit>& hitInfo
126 )
127 {
128  hitSurfaces.setSize(start.size());
129  hitSurfaces = -1;
130  hitInfo.setSize(start.size());
131 
132  // Work arrays
133  labelList hitMap(identity(start.size()));
134  pointField p0(start);
135  pointField p1(end);
136  List<pointIndexHit> intersectInfo(start.size());
137 
138  forAll(surfacesToTest, testI)
139  {
140  // Do synchronised call to all surfaces.
141  allSurfaces[surfacesToTest[testI]].findLineAny(p0, p1, intersectInfo);
142 
143  // Copy all hits into arguments, continue with misses
144  label newI = 0;
145  forAll(intersectInfo, i)
146  {
147  if (intersectInfo[i].hit())
148  {
149  hitInfo[hitMap[i]] = intersectInfo[i];
150  hitSurfaces[hitMap[i]] = testI;
151  }
152  else
153  {
154  if (i != newI)
155  {
156  hitMap[newI] = hitMap[i];
157  p0[newI] = p0[i];
158  p1[newI] = p1[i];
159  }
160  newI++;
161  }
162  }
163 
164  // All done? Note that this decision should be synchronised
165  if (newI == 0)
166  {
167  break;
168  }
169 
170  // Trim and continue
171  hitMap.setSize(newI);
172  p0.setSize(newI);
173  p1.setSize(newI);
174  intersectInfo.setSize(newI);
175  }
176 }
177 
178 
180 (
181  const PtrList<searchableSurface>& allSurfaces,
182  const labelList& surfacesToTest,
183  const pointField& start,
184  const pointField& end,
185  labelListList& hitSurfaces,
186  List<List<pointIndexHit>>& hitInfo
187 )
188 {
189  // Note: maybe move the single-surface all intersections test into
190  // searchable surface? Some of the tolerance issues might be
191  // lessened.
192 
193  // 2. Currently calling searchableSurface::findLine with start==end
194  // is expected to find no intersection. Problem if it does.
195 
196  hitSurfaces.setSize(start.size());
197  hitInfo.setSize(start.size());
198 
199  if (surfacesToTest.empty())
200  {
201  return;
202  }
203 
204  // Test first surface
205  allSurfaces[surfacesToTest[0]].findLineAll(start, end, hitInfo);
206 
207  // Set hitSurfaces and distance
208  List<scalarList> hitDistSqr(hitInfo.size());
209  forAll(hitInfo, pointi)
210  {
211  const List<pointIndexHit>& pHits = hitInfo[pointi];
212 
213  labelList& pSurfaces = hitSurfaces[pointi];
214  pSurfaces.setSize(pHits.size());
215  pSurfaces = 0;
216 
217  scalarList& pDistSqr = hitDistSqr[pointi];
218  pDistSqr.setSize(pHits.size());
219  forAll(pHits, i)
220  {
221  pDistSqr[i] = magSqr(pHits[i].hitPoint() - start[pointi]);
222  }
223  }
224 
225 
226  if (surfacesToTest.size() > 1)
227  {
228  // Test the other surfaces and merge (according to distance from start).
229  for (label testI = 1; testI < surfacesToTest.size(); testI++)
230  {
231  List<List<pointIndexHit>> surfHits;
232  allSurfaces[surfacesToTest[testI]].findLineAll
233  (
234  start,
235  end,
236  surfHits
237  );
238 
239  forAll(surfHits, pointi)
240  {
241  mergeHits
242  (
243  start[pointi], // Current segment
244 
245  testI, // Surface and its hits
246  surfHits[pointi],
247 
248  hitSurfaces[pointi], // Merge into overall hit info
249  hitInfo[pointi],
250  hitDistSqr[pointi]
251  );
252  }
253  }
254  }
255 }
256 
257 
259 (
260  const PtrList<searchableSurface>& allSurfaces,
261  const labelList& surfacesToTest,
262  const pointField& start,
263  const pointField& end,
264  labelList& surface1,
265  List<pointIndexHit>& hit1,
266  labelList& surface2,
267  List<pointIndexHit>& hit2
268 )
269 {
270  // 1. intersection from start to end
271  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
272 
273  // Initialise arguments
274  surface1.setSize(start.size());
275  surface1 = -1;
276  hit1.setSize(start.size());
277 
278  // Current end of segment to test.
279  pointField nearest(end);
280  // Work array
281  List<pointIndexHit> nearestInfo(start.size());
282 
283  forAll(surfacesToTest, testI)
284  {
285  // See if any intersection between start and current nearest
286  allSurfaces[surfacesToTest[testI]].findLine
287  (
288  start,
289  nearest,
290  nearestInfo
291  );
292 
293  forAll(nearestInfo, pointi)
294  {
295  if (nearestInfo[pointi].hit())
296  {
297  hit1[pointi] = nearestInfo[pointi];
298  surface1[pointi] = testI;
299  nearest[pointi] = hit1[pointi].hitPoint();
300  }
301  }
302  }
303 
304 
305  // 2. intersection from end to last intersection
306  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
307 
308  // Find the nearest intersection from end to start. Note that we
309  // initialise to the first intersection (if any).
310  surface2 = surface1;
311  hit2 = hit1;
312 
313  // Set current end of segment to test.
314  forAll(nearest, pointi)
315  {
316  if (hit1[pointi].hit())
317  {
318  nearest[pointi] = hit1[pointi].hitPoint();
319  }
320  else
321  {
322  // Disable testing by setting to end.
323  nearest[pointi] = end[pointi];
324  }
325  }
326 
327  forAll(surfacesToTest, testI)
328  {
329  // See if any intersection between end and current nearest
330  allSurfaces[surfacesToTest[testI]].findLine(end, nearest, nearestInfo);
331 
332  forAll(nearestInfo, pointi)
333  {
334  if (nearestInfo[pointi].hit())
335  {
336  hit2[pointi] = nearestInfo[pointi];
337  surface2[pointi] = testI;
338  nearest[pointi] = hit2[pointi].hitPoint();
339  }
340  }
341  }
342 }
343 
344 
346 (
347  const PtrList<searchableSurface>& allSurfaces,
348  const labelList& surfacesToTest,
349  const pointField& samples,
350  const scalarField& nearestDistSqr,
351  labelList& nearestSurfaces,
352  List<pointIndexHit>& nearestInfo
353 )
354 {
355  // Find nearest. Return -1 or nearest point
356 
357  // Initialise
358  nearestSurfaces.setSize(samples.size());
359  nearestSurfaces = -1;
360  nearestInfo.setSize(samples.size());
361 
362  // Work arrays
363  scalarField minDistSqr(nearestDistSqr);
364  List<pointIndexHit> hitInfo(samples.size());
365 
366  forAll(surfacesToTest, testI)
367  {
368  allSurfaces[surfacesToTest[testI]].findNearest
369  (
370  samples,
371  minDistSqr,
372  hitInfo
373  );
374 
375  // Update minDistSqr and arguments
376  forAll(hitInfo, pointi)
377  {
378  if (hitInfo[pointi].hit())
379  {
380  minDistSqr[pointi] = magSqr
381  (
382  hitInfo[pointi].hitPoint()
383  - samples[pointi]
384  );
385  nearestInfo[pointi] = hitInfo[pointi];
386  nearestSurfaces[pointi] = testI;
387  }
388  }
389  }
390 }
391 
392 
394 (
395  const PtrList<searchableSurface>& allSurfaces,
396  const labelList& surfacesToTest,
397  const pointField& samples,
398  const scalarField& nearestDistSqr,
399  const labelList& regionIndices,
400  labelList& nearestSurfaces,
401  List<pointIndexHit>& nearestInfo
402 )
403 {
404  // Find nearest. Return -1 or nearest point
405 
406  if (regionIndices.empty())
407  {
408  findNearest
409  (
410  allSurfaces,
411  surfacesToTest,
412  samples,
413  nearestDistSqr,
414  nearestSurfaces,
415  nearestInfo
416  );
417  }
418 
419  // Initialise
420  nearestSurfaces.setSize(samples.size());
421  nearestSurfaces = -1;
422  nearestInfo.setSize(samples.size());
423 
424  // Work arrays
425  scalarField minDistSqr(nearestDistSqr);
426  List<pointIndexHit> hitInfo(samples.size());
427 
428  forAll(surfacesToTest, testI)
429  {
430  allSurfaces[surfacesToTest[testI]].findNearest
431  (
432  samples,
433  minDistSqr,
434  regionIndices,
435  hitInfo
436  );
437 
438  // Update minDistSqr and arguments
439  forAll(hitInfo, pointi)
440  {
441  if (hitInfo[pointi].hit())
442  {
443  minDistSqr[pointi] = magSqr
444  (
445  hitInfo[pointi].hitPoint()
446  - samples[pointi]
447  );
448  nearestInfo[pointi] = hitInfo[pointi];
449  nearestSurfaces[pointi] = testI;
450  }
451  }
452  }
453 }
454 
455 
457 (
458  const PtrList<searchableSurface>& allSurfaces,
459  const labelList& surfacesToTest,
460  const pointField& start,
461  const scalarField& distSqr,
462  pointField& near,
463  List<pointConstraint>& constraint,
464  const label nIter
465 )
466 {
467  // Multi-surface findNearest
468 
469  vectorField normal;
470  List<pointIndexHit> info;
471 
472  allSurfaces[surfacesToTest[0]].findNearest(start, distSqr, info);
473  allSurfaces[surfacesToTest[0]].getNormal(info, normal);
474 
475  // Extract useful info from initial start point
476  near = start;
477  forAll(info, i)
478  {
479  if (info[i].hit())
480  {
481  near[i] = info[i].hitPoint();
482  }
483  }
484  constraint.setSize(near.size());
485 
486 
487  if (surfacesToTest.size() == 1)
488  {
489  constraint = pointConstraint();
490  forAll(info, i)
491  {
492  if (info[i].hit())
493  {
494  constraint[i].applyConstraint(normal[i]);
495  }
496  }
497  }
498  else if (surfacesToTest.size() >= 2)
499  {
500  // Work space
501  pointField near1;
502  vectorField normal1;
503 
504  label surfi = 1;
505  for (label iter = 0; iter < nIter; iter++)
506  {
507  constraint = pointConstraint();
508  forAll(constraint, i)
509  {
510  if (info[i].hit())
511  {
512  constraint[i].applyConstraint(normal[i]);
513  }
514  }
515 
516  // Find intersection with next surface
517  const searchableSurface& s = allSurfaces[surfacesToTest[surfi]];
518  s.findNearest(near, distSqr, info);
519  s.getNormal(info, normal1);
520  near1.setSize(info.size());
521  forAll(info, i)
522  {
523  if (info[i].hit())
524  {
525  near1[i] = info[i].hitPoint();
526  }
527  }
528 
529  // Move to intersection
530  forAll(near, pointi)
531  {
532  if (info[pointi].hit())
533  {
534  plane pl0(near[pointi], normal[pointi]);
535  plane pl1(near1[pointi], normal1[pointi]);
536  plane::ray r(pl0.planeIntersect(pl1));
537  vector n = r.dir() / mag(r.dir());
538 
539  vector d(r.refPoint()-near[pointi]);
540  d -= (d&n)*n;
541 
542  near[pointi] += d;
543  normal[pointi] = normal1[pointi];
544  constraint[pointi].applyConstraint(normal1[pointi]);
545  }
546  }
547 
548  // Step to next surface
549  surfi = surfacesToTest.fcIndex(surfi);
550  }
551  }
552 }
553 
554 
556 (
557  const PtrList<searchableSurface>& allSurfaces,
558  const labelList& surfacesToTest,
559  const pointField& samples,
560  const scalarField& nearestDistSqr,
561  const volumeType illegalHandling,
562  labelList& nearestSurfaces,
563  scalarField& distance
564 )
565 {
566  // Initialise
567  distance.setSize(samples.size());
568  distance = -great;
569 
570  // Find nearest
571  List<pointIndexHit> nearestInfo;
572  findNearest
573  (
574  allSurfaces,
575  surfacesToTest,
576  samples,
577  nearestDistSqr,
578  nearestSurfaces,
579  nearestInfo
580  );
581 
582  // Determine sign of nearest. Sort by surface to do this.
583  DynamicField<point> surfPoints(samples.size());
584  DynamicList<label> surfIndices(samples.size());
585 
586  forAll(surfacesToTest, testI)
587  {
588  // Extract samples on this surface
589  surfPoints.clear();
590  surfIndices.clear();
591  forAll(nearestSurfaces, i)
592  {
593  if (nearestSurfaces[i] == testI)
594  {
595  surfPoints.append(samples[i]);
596  surfIndices.append(i);
597  }
598  }
599 
600  // Calculate sideness of these surface points
601  List<volumeType> volType;
602  allSurfaces[surfacesToTest[testI]].getVolumeType(surfPoints, volType);
603 
604  // Push back to original
605  forAll(volType, i)
606  {
607  label pointi = surfIndices[i];
608  scalar dist = mag(samples[pointi] - nearestInfo[pointi].hitPoint());
609 
610  volumeType vT = volType[i];
611 
612  if (vT == volumeType::outside)
613  {
614  distance[pointi] = dist;
615  }
616  else if (vT == volumeType::inside)
617  {
618  distance[i] = -dist;
619  }
620  else
621  {
622  switch (illegalHandling)
623  {
624  case volumeType::outside:
625  {
626  distance[pointi] = dist;
627  break;
628  }
629  case volumeType::inside:
630  {
631  distance[pointi] = -dist;
632  break;
633  }
634  default:
635  {
637  << "getVolumeType failure,"
638  << " neither INSIDE or OUTSIDE."
639  << " point:" << surfPoints[i]
640  << " surface:"
641  << allSurfaces[surfacesToTest[testI]].name()
642  << " volType:"
643  << volumeType::names[vT]
644  << exit(FatalError);
645  break;
646  }
647  }
648  }
649  }
650  }
651 }
652 
653 
655 (
656  const PtrList<searchableSurface>& allSurfaces,
657  const labelList& surfacesToTest
658 )
659 {
660  pointField bbPoints(2*surfacesToTest.size());
661 
662  forAll(surfacesToTest, testI)
663  {
664  const searchableSurface& surface(allSurfaces[surfacesToTest[testI]]);
665 
666  bbPoints[2*testI] = surface.bounds().min();
667 
668  bbPoints[2*testI + 1] = surface.bounds().max();
669  }
670 
671  return boundBox(bbPoints);
672 }
673 
674 
675 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
bool empty() const
Return true if the UList is empty (ie, size() is zero)
Definition: UListI.H:325
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A direction and a reference point.
Definition: plane.H:72
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
static void findNearest(const PtrList< searchableSurface > &, const labelList &surfacesToTest, const pointField &, const scalarField &nearestDistSqr, labelList &surfaces, List< pointIndexHit > &)
Find nearest. Return -1 (and a miss()) or surface and nearest.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
ray planeIntersect(const plane &) const
Return the cutting line between this plane and another.
Definition: plane.C:386
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
static void findAllIntersections(const PtrList< searchableSurface > &, const labelList &surfacesToTest, const pointField &start, const pointField &end, labelListList &surfaces, List< List< pointIndexHit >> &surfaceHits)
Find all intersections in order from start to end. Returns for.
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
static void signedDistance(const PtrList< searchableSurface > &allSurfaces, const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, const volumeType illegalHandling, labelList &nearestSurfaces, scalarField &distance)
Find signed distance to nearest surface. Outside is positive.
Various functions to operate on Lists.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const =0
static void findNearestIntersection(const PtrList< searchableSurface > &allSurfaces, const labelList &surfacesToTest, const pointField &start, const pointField &end, labelList &surface1, List< pointIndexHit > &hit1, labelList &surface2, List< pointIndexHit > &hit2)
Find intersections of edge nearest to both endpoints.
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
Dynamically sized Field.
Definition: DynamicField.H:49
List< label > labelList
A List of labels.
Definition: labelList.H:56
static void findAnyIntersection(const PtrList< searchableSurface > &, const labelList &surfacesToTest, const pointField &start, const pointField &end, labelList &surfaces, List< pointIndexHit > &)
Find any intersection. Return hit point information and.
Accumulates point constraints through successive applications of the applyConstraint function...
dimensioned< scalar > magSqr(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
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
vector point
Point is a vector.
Definition: point.H:41
label findLower(const ListType &, typename ListType::const_reference, const label stary, const BinaryOp &bop)
Find last element < given value in sorted list and return index,.
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
const boundBox & bounds() const
Return const reference to boundBox.
static const NamedEnum< volumeType, 4 > names
Definition: volumeType.H:80
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const =0
From a set of points and indices get the normal.
static boundBox bounds(const PtrList< searchableSurface > &allSurfaces, const labelList &surfacesToTest)
Find the boundBox of the selected surfaces.
Namespace for OpenFOAM.