searchableSurfacesQueries.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2015-2018 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
30#include "ListOps.H"
31#include "OFstream.H"
32#include "meshTools.H"
33#include "DynamicField.H"
34#include "pointConstraint.H"
35#include "plane.H"
36
37// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38
39namespace Foam
40{
42}
43
44
45// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46
47void Foam::searchableSurfacesQueries::mergeHits
48(
49 const point& start,
50
51 const label testI, // index of surface
52 const List<pointIndexHit>& surfHits, // hits on surface
53
54 labelList& allSurfaces,
55 List<pointIndexHit>& allInfo,
56 scalarList& allDistSqr
57)
58{
59 // Given current set of hits (allSurfaces, allInfo) merge in those coming
60 // from surface surfI.
61
62 // Precalculate distances
63 scalarList surfDistSqr(surfHits.size());
64 forAll(surfHits, i)
65 {
66 surfDistSqr[i] = magSqr(surfHits[i].hitPoint() - start);
67 }
68
69 forAll(surfDistSqr, i)
70 {
71 label index = findLower(allDistSqr, surfDistSqr[i]);
72
73 // Check if equal to lower.
74 if (index >= 0)
75 {
76 // Same. Do not count.
77 //Pout<< "point:" << surfHits[i].hitPoint()
78 // << " considered same as:" << allInfo[index].hitPoint()
79 // << " within tol:" << mergeDist
80 // << endl;
81 }
82 else
83 {
84 // Check if equal to higher
85 label next = index + 1;
86
87 if (next < allDistSqr.size())
88 {
89 //Pout<< "point:" << surfHits[i].hitPoint()
90 // << " considered same as:" << allInfo[next].hitPoint()
91 // << " within tol:" << mergeDist
92 // << endl;
93 }
94 else
95 {
96 // Insert after index
97 label sz = allSurfaces.size();
98 allSurfaces.setSize(sz+1);
99 allInfo.setSize(allSurfaces.size());
100 allDistSqr.setSize(allSurfaces.size());
101 // Make space.
102 for (label j = sz-1; j > index; --j)
103 {
104 allSurfaces[j+1] = allSurfaces[j];
105 allInfo[j+1] = allInfo[j];
106 allDistSqr[j+1] = allDistSqr[j];
107 }
108 // Insert new value
109 allSurfaces[index+1] = testI;
110 allInfo[index+1] = surfHits[i];
111 allDistSqr[index+1] = surfDistSqr[i];
112 }
113 }
114 }
115}
116
117
118// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
119
120// Find any intersection
122(
123 const PtrList<searchableSurface>& allSurfaces,
124 const labelList& surfacesToTest,
125 const pointField& start,
126 const pointField& end,
127 labelList& hitSurfaces,
128 List<pointIndexHit>& hitInfo
129)
130{
131 hitSurfaces.setSize(start.size());
132 hitSurfaces = -1;
133 hitInfo.setSize(start.size());
134
135 // Work arrays
136 labelList hitMap(identity(start.size()));
137 pointField p0(start);
138 pointField p1(end);
139 List<pointIndexHit> intersectInfo(start.size());
140
141 forAll(surfacesToTest, testI)
142 {
143 // Do synchronised call to all surfaces.
144 allSurfaces[surfacesToTest[testI]].findLineAny(p0, p1, intersectInfo);
145
146 // Copy all hits into arguments, continue with misses
147 label newI = 0;
148 forAll(intersectInfo, i)
149 {
150 if (intersectInfo[i].hit())
151 {
152 hitInfo[hitMap[i]] = intersectInfo[i];
153 hitSurfaces[hitMap[i]] = testI;
154 }
155 else
156 {
157 if (i != newI)
158 {
159 hitMap[newI] = hitMap[i];
160 p0[newI] = p0[i];
161 p1[newI] = p1[i];
162 }
163 newI++;
164 }
165 }
166
167 // All done? Note that this decision should be synchronised
168 if (newI == 0)
169 {
170 break;
171 }
172
173 // Trim and continue
174 hitMap.setSize(newI);
175 p0.setSize(newI);
176 p1.setSize(newI);
177 intersectInfo.setSize(newI);
178 }
179}
180
181
183(
184 const PtrList<searchableSurface>& allSurfaces,
185 const labelList& surfacesToTest,
186 const pointField& start,
187 const pointField& end,
188 labelListList& hitSurfaces,
189 List<List<pointIndexHit>>& hitInfo
190)
191{
192 // Note: maybe move the single-surface all intersections test into
193 // searchable surface? Some of the tolerance issues might be
194 // lessened.
195
196 // 2. Currently calling searchableSurface::findLine with start==end
197 // is expected to find no intersection. Problem if it does.
198
199 hitSurfaces.setSize(start.size());
200 hitInfo.setSize(start.size());
201
202 if (surfacesToTest.empty())
203 {
204 return;
205 }
206
207 // Test first surface
208 allSurfaces[surfacesToTest[0]].findLineAll(start, end, hitInfo);
209
210 // Set hitSurfaces and distance
211 List<scalarList> hitDistSqr(hitInfo.size());
212 forAll(hitInfo, pointi)
213 {
214 const List<pointIndexHit>& pHits = hitInfo[pointi];
215
216 labelList& pSurfaces = hitSurfaces[pointi];
217 pSurfaces.setSize(pHits.size());
218 pSurfaces = 0;
219
220 scalarList& pDistSqr = hitDistSqr[pointi];
221 pDistSqr.setSize(pHits.size());
222 forAll(pHits, i)
223 {
224 pDistSqr[i] = magSqr(pHits[i].hitPoint() - start[pointi]);
225 }
226 }
227
228
229 if (surfacesToTest.size() > 1)
230 {
231 // Test the other surfaces and merge (according to distance from start).
232 for (label testI = 1; testI < surfacesToTest.size(); testI++)
233 {
234 List<List<pointIndexHit>> surfHits;
235 allSurfaces[surfacesToTest[testI]].findLineAll
236 (
237 start,
238 end,
239 surfHits
240 );
241
242 forAll(surfHits, pointi)
243 {
244 mergeHits
245 (
246 start[pointi], // Current segment
247
248 testI, // Surface and its hits
249 surfHits[pointi],
250
251 hitSurfaces[pointi], // Merge into overall hit info
252 hitInfo[pointi],
253 hitDistSqr[pointi]
254 );
255 }
256 }
257 }
258}
259
260
262(
263 const PtrList<searchableSurface>& allSurfaces,
264 const labelList& surfacesToTest,
265 const pointField& start,
266 const pointField& end,
267 labelList& surface1,
269 labelList& surface2,
271)
272{
273 // 1. intersection from start to end
274 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
275
276 // Initialize arguments
277 surface1.setSize(start.size());
278 surface1 = -1;
279 hit1.setSize(start.size());
280
281 // Current end of segment to test.
282 pointField nearest(end);
283 // Work array
284 List<pointIndexHit> nearestInfo(start.size());
285
286 forAll(surfacesToTest, testI)
287 {
288 // See if any intersection between start and current nearest
289 allSurfaces[surfacesToTest[testI]].findLine
290 (
291 start,
292 nearest,
293 nearestInfo
294 );
295
296 forAll(nearestInfo, pointi)
297 {
298 if (nearestInfo[pointi].hit())
299 {
300 hit1[pointi] = nearestInfo[pointi];
301 surface1[pointi] = testI;
302 nearest[pointi] = hit1[pointi].hitPoint();
303 }
304 }
305 }
306
307
308 // 2. intersection from end to last intersection
309 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
310
311 // Find the nearest intersection from end to start. Note that we
312 // initialize to the first intersection (if any).
313 surface2 = surface1;
314 hit2 = hit1;
315
316 // Set current end of segment to test.
317 forAll(nearest, pointi)
318 {
319 if (hit1[pointi].hit())
320 {
321 nearest[pointi] = hit1[pointi].hitPoint();
322 }
323 else
324 {
325 // Disable testing by setting to end.
326 nearest[pointi] = end[pointi];
327 }
328 }
329
330 forAll(surfacesToTest, testI)
331 {
332 // See if any intersection between end and current nearest
333 allSurfaces[surfacesToTest[testI]].findLine(end, nearest, nearestInfo);
334
335 forAll(nearestInfo, pointi)
336 {
337 if (nearestInfo[pointi].hit())
338 {
339 hit2[pointi] = nearestInfo[pointi];
340 surface2[pointi] = testI;
341 nearest[pointi] = hit2[pointi].hitPoint();
342 }
343 }
344 }
345}
346
347
349(
350 const PtrList<searchableSurface>& allSurfaces,
351 const labelList& surfacesToTest,
352 const pointField& samples,
353 const scalarField& nearestDistSqr,
354 labelList& nearestSurfaces,
355 List<pointIndexHit>& nearestInfo
356)
357{
358 // Find nearest. Return -1 or nearest point
359
360 if (samples.size() != nearestDistSqr.size())
361 {
362 FatalErrorInFunction << "Inconsistent sizes. samples:" << samples.size()
363 << " search-radius:" << nearestDistSqr.size()
364 << exit(FatalError);
365 }
366
367 // Initialise
368 nearestSurfaces.setSize(samples.size());
369 nearestSurfaces = -1;
370 nearestInfo.setSize(samples.size());
371
372 // Work arrays
373 scalarField minDistSqr(nearestDistSqr);
375
376 forAll(surfacesToTest, testI)
377 {
378 allSurfaces[surfacesToTest[testI]].findNearest
379 (
380 samples,
381 minDistSqr,
382 hitInfo
383 );
384
385 // Update minDistSqr and arguments
386 forAll(hitInfo, pointi)
387 {
388 if (hitInfo[pointi].hit())
389 {
390 minDistSqr[pointi] = magSqr
391 (
392 hitInfo[pointi].hitPoint()
393 - samples[pointi]
394 );
395 nearestInfo[pointi] = hitInfo[pointi];
396 nearestSurfaces[pointi] = testI;
397 }
398 }
399 }
400}
401
402
404(
405 const PtrList<searchableSurface>& allSurfaces,
406 const labelList& surfacesToTest,
407 const labelListList& regionIndices,
408
409 const pointField& samples,
410 const scalarField& nearestDistSqr,
411
412 labelList& nearestSurfaces,
413 List<pointIndexHit>& nearestInfo
414)
415{
416 // Find nearest. Return -1 or nearest point
417
418 if (samples.size() != nearestDistSqr.size())
419 {
420 FatalErrorInFunction << "Inconsistent sizes. samples:" << samples.size()
421 << " search-radius:" << nearestDistSqr.size()
422 << exit(FatalError);
423 }
424
425
426 if (regionIndices.empty())
427 {
428 findNearest
429 (
430 allSurfaces,
431 surfacesToTest,
432 samples,
433 nearestDistSqr,
434 nearestSurfaces,
435 nearestInfo
436 );
437 }
438
439 // Initialise
440 nearestSurfaces.setSize(samples.size());
441 nearestSurfaces = -1;
442 nearestInfo.setSize(samples.size());
443
444 // Work arrays
445 scalarField minDistSqr(nearestDistSqr);
447
448 forAll(surfacesToTest, testI)
449 {
450 allSurfaces[surfacesToTest[testI]].findNearest
451 (
452 samples,
453 minDistSqr,
454 regionIndices[testI],
455 hitInfo
456 );
457
458 // Update minDistSqr and arguments
459 forAll(hitInfo, pointi)
460 {
461 if (hitInfo[pointi].hit())
462 {
463 minDistSqr[pointi] = magSqr
464 (
465 hitInfo[pointi].hitPoint()
466 - samples[pointi]
467 );
468 nearestInfo[pointi] = hitInfo[pointi];
469 nearestSurfaces[pointi] = testI;
470 }
471 }
472 }
473}
474
475
477(
478 const PtrList<searchableSurface>& allSurfaces,
479 const labelList& surfacesToTest,
480 const pointField& start,
481 const scalarField& distSqr,
482 pointField& near,
483 List<pointConstraint>& constraint,
484 const label nIter
485)
486{
487 // Multi-surface findNearest
488
489
490 if (start.size() != distSqr.size())
491 {
492 FatalErrorInFunction << "Inconsistent sizes. samples:" << start.size()
493 << " search-radius:" << distSqr.size()
494 << exit(FatalError);
495 }
496
497
498 vectorField normal;
500
501 allSurfaces[surfacesToTest[0]].findNearest(start, distSqr, info);
502 allSurfaces[surfacesToTest[0]].getNormal(info, normal);
503
504 // Extract useful info from initial start point
505 near = start;
506 forAll(info, i)
507 {
508 if (info[i].hit())
509 {
510 near[i] = info[i].hitPoint();
511 }
512 }
513
514 // Store normal as constraint
515 constraint.setSize(near.size());
516 constraint = pointConstraint();
517 forAll(constraint, i)
518 {
519 if (info[i].hit())
520 {
521 constraint[i].applyConstraint(normal[i]);
522 }
523 }
524
525 if (surfacesToTest.size() >= 2)
526 {
527 // Work space
528 //pointField near1;
529 vectorField normal1;
530
531 label surfi = 1;
532 for (label iter = 0; iter < nIter; iter++)
533 {
534 // Find nearest on next surface
535 const searchableSurface& s = allSurfaces[surfacesToTest[surfi]];
536
537 // Update: info, normal1
538 s.findNearest(near, distSqr, info);
539 s.getNormal(info, normal1);
540
541 // Move to intersection of
542 // - previous surface(s) : near+normal
543 // - current surface : info+normal1
544 forAll(near, i)
545 {
546 if (info[i].hit())
547 {
548 if (normal[i] != vector::zero)
549 {
550 // Have previous hit. Find intersection
551 if (mag(normal[i]&normal1[i]) < 1.0-1e-6)
552 {
553 plane pl0(near[i], normal[i], false);
554 plane pl1(info[i].hitPoint(), normal1[i], false);
555
556 plane::ray r(pl0.planeIntersect(pl1));
557 vector n = r.dir() / mag(r.dir());
558
559 // Calculate vector to move onto intersection line
560 vector d(r.refPoint()-near[i]);
562
563 // Trim the max distance
564 scalar magD = mag(d);
565 if (magD > SMALL)
566 {
567 scalar maxDist = Foam::sqrt(distSqr[i]);
568 if (magD > maxDist)
569 {
570 // Clip
571 d /= magD;
572 d *= maxDist;
573 }
574
575 near[i] += d;
576 normal[i] = normal1[i];
577 constraint[i].applyConstraint(normal1[i]);
578 }
579 }
580 }
581 else
582 {
583 // First hit
584 near[i] = info[i].hitPoint();
585 normal[i] = normal1[i];
586 constraint[i].applyConstraint(normal1[i]);
587 }
588 }
589 }
590
591 // Step to next surface
592 surfi = surfacesToTest.fcIndex(surfi);
593 }
594 }
595}
596
597
599(
600 const PtrList<searchableSurface>& allSurfaces,
601 const labelList& surfacesToTest,
602 const pointField& samples,
603 const scalarField& nearestDistSqr,
604 const volumeType illegalHandling,
605 labelList& nearestSurfaces,
607)
608{
609 // Initialise
610 distance.setSize(samples.size());
611 distance = -GREAT;
612
613 // Find nearest
614 List<pointIndexHit> nearestInfo;
615 findNearest
616 (
617 allSurfaces,
618 surfacesToTest,
619 samples,
620 nearestDistSqr,
621 nearestSurfaces,
622 nearestInfo
623 );
624
625 // Determine sign of nearest. Sort by surface to do this.
626 DynamicField<point> surfPoints(samples.size());
627 DynamicList<label> surfIndices(samples.size());
628
629 forAll(surfacesToTest, testI)
630 {
631 // Extract samples on this surface
632 surfPoints.clear();
633 surfIndices.clear();
634 forAll(nearestSurfaces, i)
635 {
636 if (nearestSurfaces[i] == testI)
637 {
638 surfPoints.append(samples[i]);
639 surfIndices.append(i);
640 }
641 }
642
643 // Calculate sideness of these surface points
644 List<volumeType> volType;
645 allSurfaces[surfacesToTest[testI]].getVolumeType(surfPoints, volType);
646
647 // Push back to original
648 forAll(volType, i)
649 {
650 label pointi = surfIndices[i];
651 scalar dist = mag(samples[pointi] - nearestInfo[pointi].hitPoint());
652
653 volumeType vT = volType[i];
654
655 if (vT == volumeType::OUTSIDE)
656 {
657 distance[pointi] = dist;
658 }
659 else if (vT == volumeType::INSIDE)
660 {
661 distance[i] = -dist;
662 }
663 else
664 {
665 switch (illegalHandling)
666 {
668 {
669 distance[pointi] = dist;
670 break;
671 }
673 {
674 distance[pointi] = -dist;
675 break;
676 }
677 default:
678 {
680 << "getVolumeType failure,"
681 << " neither INSIDE or OUTSIDE."
682 << " point:" << surfPoints[i]
683 << " surface:"
684 << allSurfaces[surfacesToTest[testI]].name()
685 << " volType:" << vT.str()
686 << exit(FatalError);
687 break;
688 }
689 }
690 }
691 }
692 }
693}
694
695
697(
698 const PtrList<searchableSurface>& allSurfaces,
699 const labelUList& surfacesToTest
700)
701{
703
704 for (const label surfi : surfacesToTest)
705 {
706 bb.add(allSurfaces[surfi].bounds());
707 }
708
709 return bb;
710}
711
712
713// ************************************************************************* //
Various functions to operate on Lists.
label n
Dynamically sized Field.
Definition: DynamicField.H:65
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
void append(const T &val)
Append an element at the end of the list.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:391
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
void setSize(const label n)
Alias for resize()
Definition: List.H:218
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
label fcIndex(const label i) const noexcept
Definition: UListI.H:60
Vector< Cmpt > & removeCollinear(const Vector< Cmpt > &unitVec)
Definition: VectorI.H:142
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
static const boundBox invertedBox
A large inverted boundBox: min/max == +/- ROOTVGREAT.
Definition: boundBox.H:86
void add(const boundBox &bb)
Extend to include the second box.
Definition: boundBoxI.H:191
A reference point and direction.
Definition: plane.H:110
const point & refPoint() const noexcept
Definition: plane.H:122
const vector & dir() const noexcept
Definition: plane.H:127
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Definition: plane.H:95
ray planeIntersect(const plane &plane2) const
Return the cutting line between this plane and another.
Definition: plane.C:311
Accumulates point constraints through successive applications of the applyConstraint function.
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
A collection of tools for searchableSurfaces.
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.
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.
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.
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.
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.
static boundBox bounds(const PtrList< searchableSurface > &allSurfaces, const labelUList &surfacesToTest)
Find the boundBox of the selected surfaces.
An enumeration wrapper for classification of a location as being inside/outside of a volume.
Definition: volumeType.H:61
const word & str() const
The string representation of the volume type enumeration.
Definition: volumeType.C:62
@ OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
@ INSIDE
A location inside the volume.
Definition: volumeType.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const volScalarField & p0
Definition: EEqn.H:36
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Namespace for OpenFOAM.
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:64
label findLower(const ListType &input, const T &val, const label start, const ComparePredicate &comp)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
A non-counting (dummy) refCount.
Definition: refCount.H:59
scalarField samples(nIntervals, Zero)