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 -------------------------------------------------------------------------------
10 License
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 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(searchableCone, 0);
37  (
38  searchableSurface,
39  searchableCone,
40  dict
41  );
43  (
44  searchableSurface,
45  searchableCone,
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 
96 void 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 
312 Foam::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.
325 void 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 
611 void 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 
652 Foam::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 
708 Foam::searchableCone::searchableCone
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 :
719  searchableSurface(io),
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 
733 Foam::searchableCone::searchableCone
734 (
735  const IOobject& io,
736  const dictionary& dict
737 )
738 :
739  searchableSurface(io),
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,
770  List<pointIndexHit>& info
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,
787  List<pointIndexHit>& info
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,
869  List<pointIndexHit>& info
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 
933 void 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 
1029  const List<pointIndexHit>& info,
1030  labelList& region
1031 ) const
1032 {
1033  region.setSize(info.size());
1034  region = 0;
1035 }
1036 
1037 
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 
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 // ************************************************************************* //
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::searchableCone::findLineAny
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find any intersection on line from start to end.
Definition: searchableCone.C:866
Foam::VectorSpace< Vector< Cmpt >, Cmpt, 3 >::uniform
static Vector< Cmpt > uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:164
Foam::PointIndexHit::setIndex
void setIndex(const label index) noexcept
Set the index.
Definition: PointIndexHit.H:211
s
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))
Definition: gmvOutputSpray.H:25
Foam::searchableCone::coordinates
virtual tmp< pointField > coordinates() const
Get representative set of element coordinates.
Definition: searchableCone.C:54
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::searchableCone::findNearest
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const
Find nearest point on cylinder.
Definition: searchableCone.C:767
Foam::searchableCone::regions
virtual const wordList & regions() const
Names of regions.
Definition: searchableCone.C:755
Foam::searchableCone
Searching on (optionally hollow) cone.
Definition: searchableCone.H:108
Foam::dictionary::get
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:107
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::PointIndexHit
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:52
Foam::PointIndexHit::hit
bool hit() const noexcept
Is there a hit?
Definition: PointIndexHit.H:130
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::Field< vector >
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::searchableSurface
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
Definition: searchableSurface.H:69
Foam::addNamedToRunTimeSelectionTable
addNamedToRunTimeSelectionTable(topoSetCellSource, badQualityToCell, word, badQuality)
Foam::searchableCone::points
virtual tmp< pointField > points() const
Get the points that define the surface.
Definition: searchableCone.C:84
Foam::searchableCone::getVolumeType
virtual void getVolumeType(const pointField &points, List< volumeType > &volType) const
Determine type (inside/outside/mixed) for point.
Definition: searchableCone.C:1065
samples
scalarField samples(nIntervals, Zero)
searchableCone.H
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::PointIndexHit::rawPoint
const point_type & rawPoint() const noexcept
The point, no checks. Same as point()
Definition: PointIndexHit.H:178
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::IOobject::rename
virtual void rename(const word &newName)
Rename.
Definition: IOobject.H:506
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::searchableCone::findLine
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find nearest intersection on line from start to end.
Definition: searchableCone.C:784
Foam::normalised
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:487
Foam::PointIndexHit::setPoint
void setPoint(const point_type &p)
Set the point.
Definition: PointIndexHit.H:205
Foam::searchableCone::getRegion
virtual void getRegion(const List< pointIndexHit > &, labelList &region) const
From a set of points and indices get the region.
Definition: searchableCone.C:1028
Foam::searchableCone::boundingSpheres
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
Definition: searchableCone.C:61
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::pointIndexHit
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
Foam::Vector< scalar >
Foam::List< word >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::PointIndexHit::setHit
void setHit() noexcept
Set the hit status on.
Definition: PointIndexHit.H:193
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::searchableCone::getNormal
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
Definition: searchableCone.C:1039
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
Foam::tmp::New
static tmp< T > New(Args &&... args)
Construct tmp of T with forwarding arguments.
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::volumeType::INSIDE
A location inside the volume.
Definition: volumeType.H:68
Foam::boundBox
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:63
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
Foam::findMin
label findMin(const ListType &input, label start=0)
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::volumeType::OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
y
scalar y
Definition: LISASMDCalcMethod1.H:14
sample
Minimal example by using system/controlDict.functions: