searchableCylinder.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 | Copyright (C) 2011-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "searchableCylinder.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineTypeNameAndDebug(searchableCylinder, 0);
34  addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
35 }
36 
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
41 {
42  tmp<pointField> tCtrs(new pointField(1, 0.5*(point1_ + point2_)));
43 
44  return tCtrs;
45 }
46 
47 
49 (
50  pointField& centres,
51  scalarField& radiusSqr
52 ) const
53 {
54  centres.setSize(1);
55  centres[0] = 0.5*(point1_ + point2_);
56 
57  radiusSqr.setSize(1);
58  radiusSqr[0] = Foam::magSqr(point1_-centres[0]) + Foam::sqr(radius_);
59 
60  // Add a bit to make sure all points are tested inside
61  radiusSqr += Foam::sqr(SMALL);
62 }
63 
64 
66 {
67  tmp<pointField> tPts(new pointField(2));
68  pointField& pts = tPts.ref();
69 
70  pts[0] = point1_;
71  pts[1] = point2_;
72 
73  return tPts;
74 }
75 
76 
77 Foam::pointIndexHit Foam::searchableCylinder::findNearest
78 (
79  const point& sample,
80  const scalar nearestDistSqr
81 ) const
82 {
83  pointIndexHit info(false, sample, -1);
84 
85  vector v(sample - point1_);
86 
87  // Decompose sample-point1 into normal and parallel component
88  scalar parallel = (v & unitDir_);
89 
90  // Remove the parallel component and normalise
91  v -= parallel*unitDir_;
92  scalar magV = mag(v);
93 
94  if (magV < ROOTVSMALL)
95  {
96  v = Zero;
97  }
98  else
99  {
100  v /= magV;
101  }
102 
103  if (parallel <= 0)
104  {
105  // nearest is at point1 end cap. Clip to radius.
106  info.setPoint(point1_ + min(magV, radius_)*v);
107  }
108  else if (parallel >= magDir_)
109  {
110  // nearest is at point2 end cap. Clip to radius.
111  info.setPoint(point2_ + min(magV, radius_)*v);
112  }
113  else
114  {
115  // inbetween endcaps. Might either be nearer endcaps or cylinder wall
116 
117  // distance to endpoint: parallel or parallel-magDir
118  // distance to cylinder wall: magV-radius_
119 
120  // Nearest cylinder point
121  point cylPt;
122  if (magV < ROOTVSMALL)
123  {
124  // Point exactly on centre line. Take any point on wall.
125  vector e1 = point(1,0,0) ^ unitDir_;
126  scalar magE1 = mag(e1);
127  if (magE1 < SMALL)
128  {
129  e1 = point(0,1,0) ^ unitDir_;
130  magE1 = mag(e1);
131  }
132  e1 /= magE1;
133  cylPt = sample + radius_*e1;
134  }
135  else
136  {
137  cylPt = sample + (radius_-magV)*v;
138  }
139 
140  if (parallel < 0.5*magDir_)
141  {
142  // Project onto p1 endcap
143  point end1Pt = point1_ + min(magV, radius_)*v;
144 
145  if (magSqr(sample-cylPt) < magSqr(sample-end1Pt))
146  {
147  info.setPoint(cylPt);
148  }
149  else
150  {
151  info.setPoint(end1Pt);
152  }
153  }
154  else
155  {
156  // Project onto p2 endcap
157  point end2Pt = point2_ + min(magV, radius_)*v;
158 
159  if (magSqr(sample-cylPt) < magSqr(sample-end2Pt))
160  {
161  info.setPoint(cylPt);
162  }
163  else
164  {
165  info.setPoint(end2Pt);
166  }
167  }
168  }
169 
170  if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
171  {
172  info.setHit();
173  info.setIndex(0);
174  }
175 
176  return info;
177 }
178 
179 
180 Foam::scalar Foam::searchableCylinder::radius2(const point& pt) const
181 {
182  const vector x = (pt-point1_) ^ unitDir_;
183  return x&x;
184 }
185 
186 
187 // From http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 -
188 // intersection of cylinder with ray
189 void Foam::searchableCylinder::findLineAll
190 (
191  const point& start,
192  const point& end,
193  pointIndexHit& near,
194  pointIndexHit& far
195 ) const
196 {
197  near.setMiss();
198  far.setMiss();
199 
200  vector point1Start(start-point1_);
201  vector point2Start(start-point2_);
202  vector point1End(end-point1_);
203 
204  // Quick rejection of complete vector outside endcaps
205  scalar s1 = point1Start&unitDir_;
206  scalar s2 = point1End&unitDir_;
207 
208  if ((s1 < 0 && s2 < 0) || (s1 > magDir_ && s2 > magDir_))
209  {
210  return;
211  }
212 
213  // Line as P = start+t*V where V is unit vector and t=[0..mag(end-start)]
214  vector V(end-start);
215  scalar magV = mag(V);
216  if (magV < ROOTVSMALL)
217  {
218  return;
219  }
220  V /= magV;
221 
222 
223  // We now get the nearest intersections to start. This can either be
224  // the intersection with the end plane or with the cylinder side.
225 
226  // Get the two points (expressed in t) on the end planes. This is to
227  // clip any cylinder intersection against.
228  scalar tPoint1;
229  scalar tPoint2;
230 
231  // Maintain the two intersections with the endcaps
232  scalar tNear = VGREAT;
233  scalar tFar = VGREAT;
234 
235  {
236  scalar s = (V&unitDir_);
237  if (mag(s) > VSMALL)
238  {
239  tPoint1 = -s1/s;
240  tPoint2 = -(point2Start&unitDir_)/s;
241  if (tPoint2 < tPoint1)
242  {
243  Swap(tPoint1, tPoint2);
244  }
245  if (tPoint1 > magV || tPoint2 < 0)
246  {
247  return;
248  }
249 
250  // See if the points on the endcaps are actually inside the cylinder
251  if (tPoint1 >= 0 && tPoint1 <= magV)
252  {
253  if (radius2(start+tPoint1*V) <= sqr(radius_))
254  {
255  tNear = tPoint1;
256  }
257  }
258  if (tPoint2 >= 0 && tPoint2 <= magV)
259  {
260  if (radius2(start+tPoint2*V) <= sqr(radius_))
261  {
262  // Check if already have a near hit from point1
263  if (tNear <= 1)
264  {
265  tFar = tPoint2;
266  }
267  else
268  {
269  tNear = tPoint2;
270  }
271  }
272  }
273  }
274  else
275  {
276  // Vector perpendicular to cylinder. Check for outside already done
277  // above so just set tpoint to allow all.
278  tPoint1 = -VGREAT;
279  tPoint2 = VGREAT;
280  }
281  }
282 
283 
284  const vector x = point1Start ^ unitDir_;
285  const vector y = V ^ unitDir_;
286  const scalar d = sqr(radius_);
287 
288  // Second order equation of the form a*t^2 + b*t + c
289  const scalar a = (y&y);
290  const scalar b = 2*(x&y);
291  const scalar c = (x&x)-d;
292 
293  const scalar disc = b*b-4*a*c;
294 
295  scalar t1 = -VGREAT;
296  scalar t2 = VGREAT;
297 
298  if (disc < 0)
299  {
300  // Fully outside
301  return;
302  }
303  else if (disc < ROOTVSMALL)
304  {
305  // Single solution
306  if (mag(a) > ROOTVSMALL)
307  {
308  t1 = -b/(2*a);
309 
310  //Pout<< "single solution t:" << t1
311  // << " for start:" << start << " end:" << end
312  // << " c:" << c << endl;
313 
314  if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
315  {
316  // valid. Insert sorted.
317  if (t1 < tNear)
318  {
319  tFar = tNear;
320  tNear = t1;
321  }
322  else if (t1 < tFar)
323  {
324  tFar = t1;
325  }
326  }
327  else
328  {
329  return;
330  }
331  }
332  else
333  {
334  // Aligned with axis. Check if outside radius
335  //Pout<< "small discriminant:" << disc
336  // << " for start:" << start << " end:" << end
337  // << " magV:" << magV
338  // << " c:" << c << endl;
339  if (c > 0)
340  {
341  return;
342  }
343  }
344  }
345  else
346  {
347  if (mag(a) > ROOTVSMALL)
348  {
349  scalar sqrtDisc = sqrt(disc);
350 
351  t1 = (-b - sqrtDisc)/(2*a);
352  t2 = (-b + sqrtDisc)/(2*a);
353  if (t2 < t1)
354  {
355  Swap(t1, t2);
356  }
357 
358  if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
359  {
360  // valid. Insert sorted.
361  if (t1 < tNear)
362  {
363  tFar = tNear;
364  tNear = t1;
365  }
366  else if (t1 < tFar)
367  {
368  tFar = t1;
369  }
370  }
371  if (t2 >= 0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
372  {
373  // valid. Insert sorted.
374  if (t2 < tNear)
375  {
376  tFar = tNear;
377  tNear = t2;
378  }
379  else if (t2 < tFar)
380  {
381  tFar = t2;
382  }
383  }
384  //Pout<< "two solutions t1:" << t1 << " t2:" << t2
385  // << " for start:" << start << " end:" << end
386  // << " magV:" << magV
387  // << " c:" << c << endl;
388  }
389  else
390  {
391  // Aligned with axis. Check if outside radius
392  //Pout<< "large discriminant:" << disc
393  // << " small a:" << a
394  // << " for start:" << start << " end:" << end
395  // << " magV:" << magV
396  // << " c:" << c << endl;
397  if (c > 0)
398  {
399  return;
400  }
401  }
402  }
403 
404  // Check tNear, tFar
405  if (tNear >= 0 && tNear <= magV)
406  {
407  near.setPoint(start+tNear*V);
408  near.setHit();
409  near.setIndex(0);
410 
411  if (tFar <= magV)
412  {
413  far.setPoint(start+tFar*V);
414  far.setHit();
415  far.setIndex(0);
416  }
417  }
418  else if (tFar >= 0 && tFar <= magV)
419  {
420  near.setPoint(start+tFar*V);
421  near.setHit();
422  near.setIndex(0);
423  }
424 }
425 
426 
427 Foam::boundBox Foam::searchableCylinder::calcBounds() const
428 {
429 
430  // Adapted from
431  // http://www.gamedev.net/community/forums
432  // /topic.asp?topic_id=338522&forum_id=20&gforum_id=0
433 
434  // Let cylinder have end points A,B and radius r,
435 
436  // Bounds in direction X (same for Y and Z) can be found as:
437  // Let A.X<B.X (otherwise swap points)
438  // Good approximate lowest bound is A.X-r and highest is B.X+r (precise for
439  // capsule). At worst, in one direction it can be larger than needed by 2*r.
440 
441  // Accurate bounds for cylinder is
442  // A.X-kx*r, B.X+kx*r
443  // where
444  // 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))
445 
446  // similar thing for Y and Z
447  // (i.e.
448  // 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))
449  // 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))
450  // )
451 
452  // How derived: geometric reasoning. Bounds of cylinder is same as for 2
453  // circles centered on A and B. This sqrt thingy gives sine of angle between
454  // axis and direction, used to find projection of radius.
455 
456  vector kr
457  (
458  sqrt(sqr(unitDir_.y()) + sqr(unitDir_.z())),
459  sqrt(sqr(unitDir_.x()) + sqr(unitDir_.z())),
460  sqrt(sqr(unitDir_.x()) + sqr(unitDir_.y()))
461  );
462 
463  kr *= radius_;
464 
465  point min = point1_ - kr;
466  point max = point1_ + kr;
467 
468  min = ::Foam::min(min, point2_ - kr);
469  max = ::Foam::max(max, point2_ + kr);
470 
471  return boundBox(min, max);
472 }
473 
474 
475 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
476 
477 Foam::searchableCylinder::searchableCylinder
478 (
479  const IOobject& io,
480  const point& point1,
481  const point& point2,
482  const scalar radius
483 )
484 :
485  searchableSurface(io),
486  point1_(point1),
487  point2_(point2),
488  magDir_(mag(point2_-point1_)),
489  unitDir_((point2_-point1_)/magDir_),
490  radius_(radius)
491 {
492  bounds() = calcBounds();
493 }
494 
495 
496 Foam::searchableCylinder::searchableCylinder
497 (
498  const IOobject& io,
499  const dictionary& dict
500 )
501 :
502  searchableSurface(io),
503  point1_(dict.lookup("point1")),
504  point2_(dict.lookup("point2")),
505  magDir_(mag(point2_-point1_)),
506  unitDir_((point2_-point1_)/magDir_),
507  radius_(readScalar(dict.lookup("radius")))
508 {
509  bounds() = calcBounds();
510 }
511 
512 
513 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
514 
516 {}
517 
518 
519 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
520 
522 {
523  if (regions_.empty())
524  {
525  regions_.setSize(1);
526  regions_[0] = "region0";
527  }
528  return regions_;
529 }
530 
531 
532 void Foam::searchableCylinder::findNearest
533 (
534  const pointField& samples,
535  const scalarField& nearestDistSqr,
536  List<pointIndexHit>& info
537 ) const
538 {
539  info.setSize(samples.size());
540 
541  forAll(samples, i)
542  {
543  info[i] = findNearest(samples[i], nearestDistSqr[i]);
544  }
545 }
546 
547 
549 (
550  const pointField& start,
551  const pointField& end,
552  List<pointIndexHit>& info
553 ) const
554 {
555  info.setSize(start.size());
556 
557  forAll(start, i)
558  {
559  // Pick nearest intersection. If none intersected take second one.
561  findLineAll(start[i], end[i], info[i], b);
562  if (!info[i].hit() && b.hit())
563  {
564  info[i] = b;
565  }
566  }
567 }
568 
569 
571 (
572  const pointField& start,
573  const pointField& end,
574  List<pointIndexHit>& info
575 ) const
576 {
577  info.setSize(start.size());
578 
579  forAll(start, i)
580  {
581  // Discard far intersection
583  findLineAll(start[i], end[i], info[i], b);
584  if (!info[i].hit() && b.hit())
585  {
586  info[i] = b;
587  }
588  }
589 }
590 
591 
592 void Foam::searchableCylinder::findLineAll
593 (
594  const pointField& start,
595  const pointField& end,
597 ) const
598 {
599  info.setSize(start.size());
600 
601  forAll(start, i)
602  {
603  pointIndexHit near, far;
604  findLineAll(start[i], end[i], near, far);
605 
606  if (near.hit())
607  {
608  if (far.hit())
609  {
610  info[i].setSize(2);
611  info[i][0] = near;
612  info[i][1] = far;
613  }
614  else
615  {
616  info[i].setSize(1);
617  info[i][0] = near;
618  }
619  }
620  else
621  {
622  if (far.hit())
623  {
624  info[i].setSize(1);
625  info[i][0] = far;
626  }
627  else
628  {
629  info[i].clear();
630  }
631  }
632  }
633 }
634 
635 
637 (
638  const List<pointIndexHit>& info,
639  labelList& region
640 ) const
641 {
642  region.setSize(info.size());
643  region = 0;
644 }
645 
646 
648 (
649  const List<pointIndexHit>& info,
650  vectorField& normal
651 ) const
652 {
653  normal.setSize(info.size());
654  normal = Zero;
655 
656  forAll(info, i)
657  {
658  if (info[i].hit())
659  {
660  vector v(info[i].hitPoint() - point1_);
661 
662  // Decompose sample-point1 into normal and parallel component
663  scalar parallel = (v & unitDir_);
664 
665  // Remove the parallel component and normalise
666  v -= parallel*unitDir_;
667  scalar magV = mag(v);
668 
669  if (parallel <= 0)
670  {
671  if ((magV-radius_) < mag(parallel))
672  {
673  // either above endcap (magV<radius) or outside but closer
674  normal[i] = -unitDir_;
675  }
676  else
677  {
678  normal[i] = v/magV;
679  }
680  }
681  else if (parallel <= 0.5*magDir_)
682  {
683  // See if endcap closer or sidewall
684  if (magV >= radius_ || (radius_-magV) < parallel)
685  {
686  normal[i] = v/magV;
687  }
688  else
689  {
690  // closer to endcap
691  normal[i] = -unitDir_;
692  }
693  }
694  else if (parallel <= magDir_)
695  {
696  // See if endcap closer or sidewall
697  if (magV >= radius_ || (radius_-magV) < (magDir_-parallel))
698  {
699  normal[i] = v/magV;
700  }
701  else
702  {
703  // closer to endcap
704  normal[i] = unitDir_;
705  }
706  }
707  else // beyond cylinder
708  {
709  if ((magV-radius_) < (parallel-magDir_))
710  {
711  // above endcap
712  normal[i] = unitDir_;
713  }
714  else
715  {
716  normal[i] = v/magV;
717  }
718  }
719  }
720  }
721 }
722 
723 
725 (
726  const pointField& points,
727  List<volumeType>& volType
728 ) const
729 {
730  volType.setSize(points.size());
731  volType = volumeType::INSIDE;
732 
733  forAll(points, pointi)
734  {
735  const point& pt = points[pointi];
736 
737  vector v(pt - point1_);
738 
739  // Decompose sample-point1 into normal and parallel component
740  scalar parallel = v & unitDir_;
741 
742  if (parallel < 0)
743  {
744  // left of point1 endcap
745  volType[pointi] = volumeType::OUTSIDE;
746  }
747  else if (parallel > magDir_)
748  {
749  // right of point2 endcap
750  volType[pointi] = volumeType::OUTSIDE;
751  }
752  else
753  {
754  // Remove the parallel component
755  v -= parallel*unitDir_;
756 
757  if (mag(v) > radius_)
758  {
759  volType[pointi] = volumeType::OUTSIDE;
760  }
761  else
762  {
763  volType[pointi] = volumeType::INSIDE;
764  }
765  }
766  }
767 }
768 
769 
770 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
virtual void boundingSpheres(pointField &centres, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
bool hit() const
Is there a hit.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
void setIndex(const label index)
dimensionedScalar sqrt(const dimensionedScalar &ds)
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
void setPoint(const Point &p)
Macros for easy insertion into run-time selection tables.
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Return any intersection on segment from start to end.
virtual tmp< pointField > points() const
Get the points that define the surface.
scalar y
virtual const wordList & regions() const
Names of regions.
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))
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
void Swap(T &a, T &b)
Definition: Swap.H:43
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
void clear()
Clear the list, i.e. set size to zero.
Definition: List.C:356
virtual ~searchableCylinder()
Destructor.
static const zero Zero
Definition: zero.H:91
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
const Point & rawPoint() const
Return point with no checking.
virtual void getRegion(const List< pointIndexHit > &, labelList &region) const
From a set of points and indices get the region.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void setSize(const label)
Reset size of List.
Definition: List.C:295
vector point
Point is a vector.
Definition: point.H:41
const dimensionedScalar c
Speed of light in a vacuum.
dimensioned< scalar > mag(const dimensioned< Type > &)
A class for managing temporary objects.
Definition: PtrList.H:54
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:174
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find first intersection on segment from start to end.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
virtual tmp< pointField > coordinates() const
Get representative set of element coordinates.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451
virtual void getVolumeType(const pointField &, List< volumeType > &) const
Determine type (inside/outside/mixed) for point. unknown if.