searchableCone.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) 2015-2021 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "searchableCone.H"
30
31// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32
33namespace Foam
34{
37 (
40 dict
41 );
43 (
46 dict,
47 cone
48 );
49}
50
51
52// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53
55{
56 return tmp<pointField>::New(1, 0.5*(point1_ + point2_));
57}
58
59
61(
62 pointField& centres,
63 scalarField& radiusSqr
64) const
65{
66 centres.setSize(1);
67 centres[0] = 0.5*(point1_ + point2_);
68
69 radiusSqr.setSize(1);
70 if (radius1_ > radius2_)
71 {
72 radiusSqr[0] = Foam::magSqr(point1_-centres[0]) + Foam::sqr(radius1_);
73 }
74 else
75 {
76 radiusSqr[0] = Foam::magSqr(point2_-centres[0]) + Foam::sqr(radius2_);
77 }
78
79 // Add a bit to make sure all points are tested inside
80 radiusSqr += Foam::sqr(SMALL);
81}
82
83
85{
86 auto tpts = tmp<pointField>::New(2);
87 auto& pts = tpts.ref();
88
89 pts[0] = point1_;
90 pts[1] = point2_;
91
92 return tpts;
93}
94
95
96void Foam::searchableCone::findNearestAndNormal
97(
98 const point& sample,
99 const scalar nearestDistSqr,
100 pointIndexHit& info,
101 vector& nearNormal
102) const
103{
104 vector v(sample - point1_);
105
106 // Decompose sample-point1 into normal and parallel component
107 const scalar parallel = (v & unitDir_);
108
109 // Remove the parallel component and normalise
110 v -= parallel*unitDir_;
111
112 const scalar magV = mag(v);
113 v.normalise();
114
115 // Nearest and normal on disk at point1
116 point disk1Point(point1_ + min(max(magV, innerRadius1_), radius1_)*v);
117 vector disk1Normal(-unitDir_);
118
119 // Nearest and normal on disk at point2
120 point disk2Point(point2_ + min(max(magV, innerRadius2_), radius2_)*v);
121 vector disk2Normal(unitDir_);
122
123 // Nearest and normal on cone. Initialise to far-away point so if not
124 // set it picks one of the disk points
125 point nearCone(point::uniform(GREAT));
126 vector normalCone(1, 0, 0);
127
128 // Nearest and normal on inner cone. Initialise to far-away point
129 point iCnearCone(point::uniform(GREAT));
130 vector iCnormalCone(1, 0, 0);
131
132 point projPt1 = point1_+ radius1_*v;
133 point projPt2 = point2_+ radius2_*v;
134
135 vector p1 = (projPt1 - point1_);
136 if (mag(p1) > ROOTVSMALL)
137 {
138 p1 /= mag(p1);
139
140 // Find vector along the two end of cone
141 const vector b = normalised(projPt2 - projPt1);
142
143 // Find the vector along sample pt and pt at one end of cone
144 vector a = (sample - projPt1);
145
146 if (mag(a) <= ROOTVSMALL)
147 {
148 // Exception: sample on disk1. Redo with projPt2.
149 a = (sample - projPt2);
150
151 // Find normal unitvector
152 nearCone = (a & b)*b + projPt2;
153
154 vector b1 = (p1 & b)*b;
155 normalCone = normalised(p1 - b1);
156 }
157 else
158 {
159 // Find nearest point on cone surface
160 nearCone = (a & b)*b + projPt1;
161
162 // Find projection along surface of cone
163 vector b1 = (p1 & b)*b;
164 normalCone = normalised(p1 - b1);
165 }
166
167 if (innerRadius1_ > 0 || innerRadius2_ > 0)
168 {
169 // Same for inner radius
170 point iCprojPt1 = point1_+ innerRadius1_*v;
171 point iCprojPt2 = point2_+ innerRadius2_*v;
172
173 const vector iCp1 = normalised(iCprojPt1 - point1_);
174
175 // Find vector along the two end of cone
176 const vector iCb = normalised(iCprojPt2 - iCprojPt1);
177
178
179 // Find the vector along sample pt and pt at one end of conde
180 vector iCa(sample - iCprojPt1);
181
182 if (mag(iCa) <= ROOTVSMALL)
183 {
184 iCa = (sample - iCprojPt2);
185
186 // Find normal unitvector
187 iCnearCone = (iCa & iCb)*iCb+iCprojPt2;
188
189 vector b1 = (iCp1 & iCb)*iCb;
190 iCnormalCone = normalised(iCp1 - b1);
191 }
192 else
193 {
194 // Find nearest point on cone surface
195 iCnearCone = (iCa & iCb)*iCb+iCprojPt1;
196
197 // Find projection along surface of cone
198 vector b1 = (iCp1 & iCb)*iCb;
199 iCnormalCone = normalised(iCp1 - b1);
200 }
201 }
202 }
203
204
205 // Select nearest out of the 4 points (outer cone, disk1, disk2, inner
206 // cone)
207
208 FixedList<scalar, 4> dist;
209 dist[0] = magSqr(nearCone-sample);
210 dist[1] = magSqr(disk1Point-sample);
211 dist[2] = magSqr(disk2Point-sample);
212 dist[3] = magSqr(iCnearCone-sample);
213
214 const label minI = findMin(dist);
215
216
217 // Snap the point to the corresponding surface
218
219 if (minI == 0) // Outer cone
220 {
221 // Closest to (infinite) outer cone. See if needs clipping to end disks
222
223 {
224 vector v1(nearCone-point1_);
225 scalar para = (v1 & unitDir_);
226 // Remove the parallel component and normalise
227 v1 -= para*unitDir_;
228 const scalar magV1 = mag(v1);
229 v1 = v1/magV1;
230
231 if (para < 0.0 && magV1 >= radius1_)
232 {
233 // Near point 1. Set point to intersection of disk and cone.
234 // Keep normal from cone.
235 nearCone = disk1Point;
236 }
237 else if (para < 0.0 && magV1 < radius1_)
238 {
239 // On disk1
240 nearCone = disk1Point;
241 normalCone = disk1Normal;
242 }
243 else if (para > magDir_ && magV1 >= radius2_)
244 {
245 // Near point 2. Set point to intersection of disk and cone.
246 // Keep normal from cone.
247 nearCone = disk2Point;
248 }
249 else if (para > magDir_ && magV1 < radius2_)
250 {
251 // On disk2
252 nearCone = disk2Point;
253 normalCone = disk2Normal;
254 }
255 }
256 info.setPoint(nearCone);
257 nearNormal = normalCone;
258 }
259 else if (minI == 1) // Near to disk1
260 {
261 info.setPoint(disk1Point);
262 nearNormal = disk1Normal;
263 }
264 else if (minI == 2) // Near to disk2
265 {
266 info.setPoint(disk2Point);
267 nearNormal = disk2Normal;
268 }
269 else // Near to inner cone
270 {
271 {
272 vector v1(iCnearCone-point1_);
273 scalar para = (v1 & unitDir_);
274 // Remove the parallel component and normalise
275 v1 -= para*unitDir_;
276
277 const scalar magV1 = mag(v1);
278 v1 = v1/magV1;
279
280 if (para < 0.0 && magV1 >= innerRadius1_)
281 {
282 iCnearCone = disk1Point;
283 }
284 else if (para < 0.0 && magV1 < innerRadius1_)
285 {
286 iCnearCone = disk1Point;
287 iCnormalCone = disk1Normal;
288 }
289 else if (para > magDir_ && magV1 >= innerRadius2_)
290 {
291 iCnearCone = disk2Point;
292 }
293 else if (para > magDir_ && magV1 < innerRadius2_)
294 {
295 iCnearCone = disk2Point;
296 iCnormalCone = disk2Normal;
297 }
298 }
299 info.setPoint(iCnearCone);
300 nearNormal = iCnormalCone;
301 }
302
303
304 if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
305 {
306 info.setHit();
307 info.setIndex(0);
308 }
309}
310
311
312Foam::scalar Foam::searchableCone::radius2
313(
314 const searchableCone& cone,
315 const point& pt
316)
317{
318 const vector x = (pt-cone.point1_) ^ cone.unitDir_;
319 return x&x;
320}
321
322
323// From mrl.nyu.edu/~dzorin/rend05/lecture2.pdf,
324// Ray Tracing II, Infinite cone ray intersection.
325void Foam::searchableCone::findLineAll
326(
327 const searchableCone& cone,
328 const scalar innerRadius1,
329 const scalar innerRadius2,
330 const point& start,
331 const point& end,
332 pointIndexHit& near,
333 pointIndexHit& far
334) const
335{
336 near.setMiss();
337 far.setMiss();
338
339 vector point1Start(start-cone.point1_);
340 vector point2Start(start-cone.point2_);
341 vector point1End(end-cone.point1_);
342
343 // Quick rejection of complete vector outside endcaps
344 scalar s1 = point1Start & (cone.unitDir_);
345 scalar s2 = point1End & (cone.unitDir_);
346
347 if ((s1 < 0.0 && s2 < 0.0) || (s1 > cone.magDir_ && s2 > cone.magDir_))
348 {
349 return;
350 }
351
352 // Line as P = start+t*V where V is unit vector and t=[0..mag(end-start)]
353 vector V(end-start);
354 scalar magV = mag(V);
355 if (magV < ROOTVSMALL)
356 {
357 return;
358 }
359 V /= magV;
360
361
362 // We now get the nearest intersections to start. This can either be
363 // the intersection with the end plane or with the cylinder side.
364
365 // Get the two points (expressed in t) on the end planes. This is to
366 // clip any cylinder intersection against.
367 scalar tPoint1;
368 scalar tPoint2;
369
370 // Maintain the two intersections with the endcaps
371 scalar tNear = VGREAT;
372 scalar tFar = VGREAT;
373
374 scalar radius_sec = cone.radius1_;
375
376 {
377 // Find dot product: mag(s)>VSMALL suggest that it is greater
378 scalar s = (V & unitDir_);
379 if (mag(s) > VSMALL)
380 {
381 tPoint1 = -s1/s;
382 tPoint2 = -(point2Start&(cone.unitDir_))/s;
383
384 if (tPoint2 < tPoint1)
385 {
386 std::swap(tPoint1, tPoint2);
387 }
388 if (tPoint1 > magV || tPoint2 < 0)
389 {
390 return;
391 }
392 // See if the points on the endcaps are actually inside the cylinder
393 if (tPoint1 >= 0 && tPoint1 <= magV)
394 {
395 scalar r2 = radius2(cone, start+tPoint1*V);
396 vector p1 = (start+tPoint1*V-point1_);
397 vector p2 = (start+tPoint1*V-point2_);
398 radius_sec = cone.radius1_;
399 scalar inC_radius_sec = innerRadius1_;
400
401 if (mag(p2&(cone.unitDir_)) < mag(p1&(cone.unitDir_)))
402 {
403 radius_sec = cone.radius2_;
404 inC_radius_sec = innerRadius2_;
405 }
406
407 if (r2 <= sqr(radius_sec) && r2 >= sqr(inC_radius_sec))
408 {
409 tNear = tPoint1;
410 }
411 }
412 if (tPoint2 >= 0 && tPoint2 <= magV)
413 {
414 // Decompose sample-point1 into normal and parallel component
415 vector p1 = (start+tPoint2*V-cone.point1_);
416 vector p2 = (start+tPoint2*V-cone.point2_);
417 radius_sec = cone.radius1_;
418 scalar inC_radius_sec = innerRadius1_;
419 if (mag(p2&(cone.unitDir_)) < mag(p1&(cone.unitDir_)))
420 {
421 radius_sec = cone.radius2_;
422 inC_radius_sec = innerRadius2_;
423 }
424 scalar r2 = radius2(cone, start+tPoint2*V);
425 if (r2 <= sqr(radius_sec) && r2 >= sqr(inC_radius_sec))
426 {
427 // Check if already have a near hit from point1
428 if (tNear <= 1)
429 {
430 tFar = tPoint2;
431 }
432 else
433 {
434 tNear = tPoint2;
435 }
436 }
437 }
438 }
439 else
440 {
441 // Vector perpendicular to cylinder. Check for outside already done
442 // above so just set tpoint to allow all.
443 tPoint1 = -VGREAT;
444 tPoint2 = VGREAT;
445 }
446 }
447
448
449 // Second order equation of the form a*t^2 + b*t + c
450 scalar a, b, c;
451
452 scalar deltaRadius = cone.radius2_-cone.radius1_;
453 if (mag(deltaRadius) <= ROOTVSMALL)
454 {
455 vector point1Start(start-cone.point1_);
456 const vector x = point1Start ^ cone.unitDir_;
457 const vector y = V ^ cone.unitDir_;
458 const scalar d = sqr(0.5*(cone.radius1_ + cone.radius2_));
459
460 a = (y&y);
461 b = 2*(x&y);
462 c = (x&x)-d;
463 }
464 else
465 {
466 vector va = cone.unitDir_;
467 vector v1 = normalised(end-start);
468 scalar p = (va&v1);
469 vector a1 = (v1-p*va);
470
471 // Determine the end point of the cone
472 point pa =
473 cone.unitDir_*cone.radius1_*mag(cone.point2_-cone.point1_)
474 /(-deltaRadius)
475 +cone.point1_;
476
477 scalar l2 = sqr(deltaRadius)+sqr(cone.magDir_);
478 scalar sqrCosAlpha = sqr(cone.magDir_)/l2;
479 scalar sqrSinAlpha = sqr(deltaRadius)/l2;
480
481
482 vector delP(start-pa);
483 vector p1 = (delP-(delP&va)*va);
484
485 a = sqrCosAlpha*((v1-p*va)&(v1-p*va))-sqrSinAlpha*p*p;
486 b =
487 2.0*sqrCosAlpha*(a1&p1)
488 -2.0*sqrSinAlpha*(v1&va)*(delP&va);
489 c =
490 sqrCosAlpha
491 *((delP-(delP&va)*va)&(delP-(delP&va)*va))
492 -sqrSinAlpha*sqr(delP&va);
493 }
494
495 const scalar disc = b*b-4.0*a*c;
496
497 scalar t1 = -VGREAT;
498 scalar t2 = VGREAT;
499
500 if (disc < 0)
501 {
502 // Fully outside
503 return;
504 }
505 else if (disc < ROOTVSMALL)
506 {
507 // Single solution
508 if (mag(a) > ROOTVSMALL)
509 {
510 t1 = -b/(2.0*a);
511
512 if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
513 {
514 // Valid. Insert sorted.
515 if (t1 < tNear)
516 {
517 tFar = tNear;
518 tNear = t1;
519 }
520 else if (t1 < tFar)
521 {
522 tFar = t1;
523 }
524 }
525 else
526 {
527 return;
528 }
529 }
530 else
531 {
532 // Aligned with axis. Check if outside radius
533 if (c > 0.0)
534 {
535 return;
536 }
537 }
538 }
539 else
540 {
541 if (mag(a) > ROOTVSMALL)
542 {
543 scalar sqrtDisc = sqrt(disc);
544
545 t1 = (-b - sqrtDisc)/(2.0*a);
546 t2 = (-b + sqrtDisc)/(2.0*a);
547 if (t2 < t1)
548 {
549 std::swap(t1, t2);
550 }
551
552 if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
553 {
554 // Valid. Insert sorted.
555 if (t1 < tNear)
556 {
557 tFar = tNear;
558 tNear = t1;
559 }
560 else if (t1 < tFar)
561 {
562 tFar = t1;
563 }
564 }
565 if (t2>=0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
566 {
567 // Valid. Insert sorted.
568 if (t2 < tNear)
569 {
570 tFar = tNear;
571 tNear = t2;
572 }
573 else if (t2 < tFar)
574 {
575 tFar = t2;
576 }
577 }
578 }
579 else
580 {
581 // Aligned with axis. Check if outside radius
582 if (c > 0.0)
583 {
584 return;
585 }
586 }
587 }
588
589 // Check tNear, tFar
590 if (tNear>=0 && tNear <= magV)
591 {
592 near.setPoint(start+tNear*V);
593 near.setHit();
594 near.setIndex(0);
595 if (tFar <= magV)
596 {
597 far.setPoint(start+tFar*V);
598 far.setHit();
599 far.setIndex(0);
600 }
601 }
602 else if (tFar>=0 && tFar <= magV)
603 {
604 near.setPoint(start+tFar*V);
605 near.setHit();
606 near.setIndex(0);
607 }
608}
609
610
611void Foam::searchableCone::insertHit
612(
613 const point& start,
614 const point& end,
615 List<pointIndexHit>& info,
616 const pointIndexHit& hit
617) const
618{
619 scalar smallDistSqr = SMALL*magSqr(end-start);
620
621 scalar hitMagSqr = magSqr(hit.hitPoint()-start);
622
623 forAll(info, i)
624 {
625 scalar d2 = magSqr(info[i].hitPoint()-start);
626
627 if (d2 > hitMagSqr+smallDistSqr)
628 {
629 // Insert at i.
630 label sz = info.size();
631 info.setSize(sz+1);
632 for (label j = sz; j > i; --j)
633 {
634 info[j] = info[j-1];
635 }
636 info[i] = hit;
637 return;
638 }
639 else if (d2 > hitMagSqr-smallDistSqr)
640 {
641 // hit is same point as info[i].
642 return;
643 }
644 }
645 // Append
646 label sz = info.size();
647 info.setSize(sz+1);
648 info[sz] = hit;
649}
650
651
652Foam::boundBox Foam::searchableCone::calcBounds() const
653{
654 // Adapted from
655 // http://www.gamedev.net/community/forums
656 // /topic.asp?topic_id=338522&forum_id=20&gforum_id=0
657
658 // Let cylinder have end points A,B and radius r,
659
660 // Bounds in direction X (same for Y and Z) can be found as:
661 // Let A.X<B.X (otherwise swap points)
662 // Good approximate lowest bound is A.X-r and highest is B.X+r (precise for
663 // capsule). At worst, in one direction it can be larger than needed by 2*r.
664
665 // Accurate bounds for cylinder is
666 // A.X-kx*r, B.X+kx*r
667 // where
668 // kx=sqrt(((A.Y-B.Y)^2+(A.Z-B.Z)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
669
670 // similar thing for Y and Z
671 // (i.e.
672 // ky=sqrt(((A.X-B.X)^2+(A.Z-B.Z)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
673 // kz=sqrt(((A.X-B.X)^2+(A.Y-B.Y)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
674 // )
675
676 // How derived: geometric reasoning. Bounds of cylinder is same as for 2
677 // circles centered on A and B. This sqrt thingy gives sine of angle between
678 // axis and direction, used to find projection of radius.
679
680 vector kr
681 (
682 sqrt(sqr(unitDir_.y()) + sqr(unitDir_.z())),
683 sqrt(sqr(unitDir_.x()) + sqr(unitDir_.z())),
684 sqrt(sqr(unitDir_.x()) + sqr(unitDir_.y()))
685 );
686
687 if (radius2_ >= radius1_)
688 {
689 kr *= radius2_;
690 }
691 else
692 {
693 kr *= radius1_;
694 }
695
696 point min = point1_ - kr;
697 point max = point1_ + kr;
698
699 min = ::Foam::min(min, point2_ - kr);
700 max = ::Foam::max(max, point2_ + kr);
701
702 return boundBox(min, max);
703}
704
705
706// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
707
709(
710 const IOobject& io,
711 const point& point1,
712 const scalar radius1,
713 const scalar innerRadius1,
714 const point& point2,
715 const scalar radius2,
716 const scalar innerRadius2
717)
718:
720 point1_(point1),
721 radius1_(radius1),
722 innerRadius1_(innerRadius1),
723 point2_(point2),
724 radius2_(radius2),
725 innerRadius2_(innerRadius2),
726 magDir_(mag(point2_-point1_)),
727 unitDir_((point2_-point1_)/magDir_)
728{
729 bounds() = calcBounds();
730}
731
732
734(
735 const IOobject& io,
736 const dictionary& dict
737)
738:
740 point1_(dict.get<point>("point1")),
741 radius1_(dict.get<scalar>("radius1")),
742 innerRadius1_(dict.getOrDefault<scalar>("innerRadius1", 0)),
743 point2_(dict.get<point>("point2")),
744 radius2_(dict.get<scalar>("radius2")),
745 innerRadius2_(dict.getOrDefault<scalar>("innerRadius2", 0)),
746 magDir_(mag(point2_-point1_)),
747 unitDir_((point2_-point1_)/magDir_)
748{
749 bounds() = calcBounds();
750}
751
752
753// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
754
756{
757 if (regions_.empty())
758 {
759 regions_.setSize(1);
760 regions_[0] = "region0";
761 }
762 return regions_;
763}
764
765
767(
768 const pointField& samples,
769 const scalarField& nearestDistSqr,
771) const
772{
773 info.setSize(samples.size());
774
775 forAll(samples, i)
776 {
777 vector normal;
778 findNearestAndNormal(samples[i], nearestDistSqr[i], info[i], normal);
779 }
780}
781
782
784(
785 const pointField& start,
786 const pointField& end,
788) const
789{
790 info.setSize(start.size());
791
792 forAll(start, i)
793 {
794 // Pick nearest intersection. If none intersected take second one.
796 findLineAll
797 (
798 *this,
799 innerRadius1_,
800 innerRadius2_,
801 start[i],
802 end[i],
803 info[i],
804 b
805 );
806 if (!info[i].hit() && b.hit())
807 {
808 info[i] = b;
809 }
810 }
811
812
813 if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
814 {
815 IOobject io(*this);
816 io.rename(name()+"Inner");
817 searchableCone innerCone
818 (
819 io,
820 point1_,
821 innerRadius1_,
822 0.0,
823 point2_,
824 innerRadius2_,
825 0.0
826 );
827
828 forAll(info, i)
829 {
830 point newEnd;
831 if (info[i].hit())
832 {
833 newEnd = info[i].hitPoint();
834 }
835 else
836 {
837 newEnd = end[i];
838 }
839 pointIndexHit near;
840 pointIndexHit far;
841 findLineAll
842 (
843 innerCone,
844 innerRadius1_,
845 innerRadius2_,
846 start[i],
847 newEnd,
848 near,
849 far
850 );
851
852 if (near.hit())
853 {
854 info[i] = near;
855 }
856 else if (far.hit())
857 {
858 info[i] = far;
859 }
860 }
861 }
862}
863
864
866(
867 const pointField& start,
868 const pointField& end,
870) const
871{
872 info.setSize(start.size());
873
874 forAll(start, i)
875 {
876 // Discard far intersection
878 findLineAll
879 (
880 *this,
881 innerRadius1_,
882 innerRadius2_,
883 start[i],
884 end[i],
885 info[i],
886 b
887 );
888 if (!info[i].hit() && b.hit())
889 {
890 info[i] = b;
891 }
892 }
893 if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
894 {
895 IOobject io(*this);
896 io.rename(name()+"Inner");
897 searchableCone cone
898 (
899 io,
900 point1_,
901 innerRadius1_,
902 0.0,
903 point2_,
904 innerRadius2_,
905 0.0
906 );
907
908 forAll(info, i)
909 {
910 if (!info[i].hit())
911 {
913 findLineAll
914 (
915 cone,
916 innerRadius1_,
917 innerRadius2_,
918 start[i],
919 end[i],
920 info[i],
921 b
922 );
923 if (!info[i].hit() && b.hit())
924 {
925 info[i] = b;
926 }
927 }
928 }
929 }
930}
931
932
933void Foam::searchableCone::findLineAll
934(
935 const pointField& start,
936 const pointField& end,
938) const
939{
940 info.setSize(start.size());
941
942 forAll(start, i)
943 {
944 pointIndexHit near, far;
945 findLineAll
946 (
947 *this,
948 innerRadius1_,
949 innerRadius2_,
950 start[i],
951 end[i],
952 near,
953 far
954 );
955
956 if (near.hit())
957 {
958 if (far.hit())
959 {
960 info[i].setSize(2);
961 info[i][0] = near;
962 info[i][1] = far;
963 }
964 else
965 {
966 info[i].setSize(1);
967 info[i][0] = near;
968 }
969 }
970 else
971 {
972 if (far.hit())
973 {
974 info[i].setSize(1);
975 info[i][0] = far;
976 }
977 else
978 {
979 info[i].clear();
980 }
981 }
982 }
983
984 if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
985 {
986 IOobject io(*this);
987 io.rename(name()+"Inner");
988 searchableCone cone
989 (
990 io,
991 point1_,
992 innerRadius1_,
993 0.0,
994 point2_,
995 innerRadius2_,
996 0.0
997 );
998
999 forAll(info, i)
1000 {
1001 pointIndexHit near;
1002 pointIndexHit far;
1003 findLineAll
1004 (
1005 cone,
1006 innerRadius1_,
1007 innerRadius2_,
1008 start[i],
1009 end[i],
1010 near,
1011 far
1012 );
1013
1014 if (near.hit())
1015 {
1016 insertHit(start[i], end[i], info[i], near);
1017 }
1018 if (far.hit())
1019 {
1020 insertHit(start[i], end[i], info[i], far);
1021 }
1022 }
1023 }
1024}
1025
1026
1028(
1029 const List<pointIndexHit>& info,
1030 labelList& region
1031) const
1032{
1033 region.setSize(info.size());
1034 region = 0;
1035}
1036
1037
1039(
1040 const List<pointIndexHit>& info,
1041 vectorField& normal
1042) const
1043{
1044 normal.setSize(info.size());
1045 normal = Zero;
1046
1047 forAll(info, i)
1048 {
1049 if (info[i].hit())
1050 {
1051 pointIndexHit nearInfo;
1052 findNearestAndNormal
1053 (
1054 info[i].hitPoint(),
1055 Foam::sqr(GREAT),
1056 nearInfo,
1057 normal[i]
1058 );
1059 }
1060 }
1061}
1062
1063
1065(
1066 const pointField& points,
1067 List<volumeType>& volType
1068) const
1069{
1070 volType.setSize(points.size());
1071
1072 forAll(points, pointi)
1073 {
1074 const point& pt = points[pointi];
1075
1076 volType[pointi] = volumeType::OUTSIDE;
1077
1078 vector v(pt - point1_);
1079
1080 // Decompose sample-point1 into normal and parallel component
1081 const scalar parallel = (v & unitDir_);
1082
1083 // Quick rejection. Left of point1 endcap, or right of point2 endcap
1084 if (parallel < 0 || parallel > magDir_)
1085 {
1086 continue;
1087 }
1088
1089 const scalar radius_sec =
1090 radius1_ + parallel * (radius2_-radius1_)/magDir_;
1091
1092 const scalar radius_sec_inner =
1093 innerRadius1_ + parallel * (innerRadius2_-innerRadius1_)/magDir_;
1094
1095 // Remove the parallel component
1096 v -= parallel*unitDir_;
1097
1098 if (mag(v) >= radius_sec_inner && mag(v) <= radius_sec)
1099 {
1100 volType[pointi] = volumeType::INSIDE;
1101 }
1102 }
1103}
1104
1105
1106// ************************************************************************* //
scalar y
Macros for easy insertion into run-time selection tables.
#define addNamedToRunTimeSelectionTable(baseType, thisType, argNames, lookupName)
Add to construction table with 'lookupName' as the key.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Minimal example by using system/controlDict.functions:
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
virtual void rename(const word &newName)
Rename the object.
Definition: IOobject.H:497
void setSize(const label n)
Alias for resize()
Definition: List.H:218
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
void setHit() noexcept
Set the hit status on.
const point_type & rawPoint() const noexcept
The point, no checks. Same as point()
void setIndex(const label index) noexcept
Set the index.
void setPoint(const point_type &p)
Set the point.
bool hit() const noexcept
Is there a hit?
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Particle-size distribution model wherein random samples are drawn from the doubly-truncated uniform p...
Definition: uniform.H:164
Searching on (optionally hollow) cone.
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find nearest intersection on line from start to end.
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find any intersection on line from start to end.
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
virtual void getVolumeType(const pointField &points, List< volumeType > &volType) const
Determine type (inside/outside/mixed) for point.
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 void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const
Find nearest point on cylinder.
virtual tmp< pointField > coordinates() const
Get representative set of element coordinates.
virtual tmp< pointField > points() const
Get the points that define the surface.
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
virtual const boundBox & bounds() const
Return const reference to boundBox.
A class for managing temporary objects.
Definition: tmp.H:65
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
@ 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
volScalarField & p
const pointField & points
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))
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
const dimensionedScalar c
Speed of light in a vacuum.
Namespace for OpenFOAM.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
label findMin(const ListType &input, label start=0)
vector point
Point is a vector.
Definition: point.H:43
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dictionary dict
volScalarField & b
Definition: createFields.H:27
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
scalarField samples(nIntervals, Zero)