FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2011 OpenCFD Ltd.
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 
34 defineTypeNameAndDebug(searchableCylinder, 0);
35 addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
36 
37 }
38 
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
43 {
44  pointField ctrs(1, 0.5*(point1_ + point2_));
45 
46  return ctrs;
47 }
48 
49 
50 Foam::pointIndexHit Foam::searchableCylinder::findNearest
51 (
52  const point& sample,
53  const scalar nearestDistSqr
54 ) const
55 {
56  pointIndexHit info(false, sample, -1);
57 
58  vector v(sample - point1_);
59 
60  // Decompose sample-point1 into normal and parallel component
61  scalar parallel = (v & unitDir_);
62 
63  // Remove the parallel component and normalise
64  v -= parallel*unitDir_;
65  scalar magV = mag(v);
66 
67  if (magV < ROOTVSMALL)
68  {
69  v = vector::zero;
70  }
71  else
72  {
73  v /= magV;
74  }
75 
76  if (parallel <= 0)
77  {
78  // nearest is at point1 end cap. Clip to radius.
79  info.setPoint(point1_ + min(magV, radius_)*v);
80  }
81  else if (parallel >= magDir_)
82  {
83  // nearest is at point2 end cap. Clip to radius.
84  info.setPoint(point2_ + min(magV, radius_)*v);
85  }
86  else
87  {
88  // inbetween endcaps. Might either be nearer endcaps or cylinder wall
89 
90  // distance to endpoint: parallel or parallel-magDir
91  // distance to cylinder wall: magV-radius_
92 
93  // Nearest cylinder point
94  point cylPt;
95  if (magV < ROOTVSMALL)
96  {
97  // Point exactly on centre line. Take any point on wall.
98  vector e1 = point(1,0,0) ^ unitDir_;
99  scalar magE1 = mag(e1);
100  if (magE1 < SMALL)
101  {
102  e1 = point(0,1,0) ^ unitDir_;
103  magE1 = mag(e1);
104  }
105  e1 /= magE1;
106  cylPt = sample + radius_*e1;
107  }
108  else
109  {
110  cylPt = sample + (radius_-magV)*v;
111  }
112 
113  if (parallel < 0.5*magDir_)
114  {
115  // Project onto p1 endcap
116  point end1Pt = point1_ + min(magV, radius_)*v;
117 
118  if (magSqr(sample-cylPt) < magSqr(sample-end1Pt))
119  {
120  info.setPoint(cylPt);
121  }
122  else
123  {
124  info.setPoint(end1Pt);
125  }
126  }
127  else
128  {
129  // Project onto p2 endcap
130  point end2Pt = point2_ + min(magV, radius_)*v;
131 
132  if (magSqr(sample-cylPt) < magSqr(sample-end2Pt))
133  {
134  info.setPoint(cylPt);
135  }
136  else
137  {
138  info.setPoint(end2Pt);
139  }
140  }
141  }
142 
143  if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
144  {
145  info.setHit();
146  info.setIndex(0);
147  }
148 
149  return info;
150 }
151 
152 
153 Foam::scalar Foam::searchableCylinder::radius2(const point& pt) const
154 {
155  const vector x = (pt-point1_) ^ unitDir_;
156  return x&x;
157 }
158 
159 
160 // From http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 -
161 // intersection of cylinder with ray
162 void Foam::searchableCylinder::findLineAll
163 (
164  const point& start,
165  const point& end,
166  pointIndexHit& near,
167  pointIndexHit& far
168 ) const
169 {
170  near.setMiss();
171  far.setMiss();
172 
173  vector point1Start(start-point1_);
174  vector point2Start(start-point2_);
175  vector point1End(end-point1_);
176 
177  // Quick rejection of complete vector outside endcaps
178  scalar s1 = point1Start&unitDir_;
179  scalar s2 = point1End&unitDir_;
180 
181  if ((s1 < 0 && s2 < 0) || (s1 > magDir_ && s2 > magDir_))
182  {
183  return;
184  }
185 
186  // Line as P = start+t*V where V is unit vector and t=[0..mag(end-start)]
187  vector V(end-start);
188  scalar magV = mag(V);
189  if (magV < ROOTVSMALL)
190  {
191  return;
192  }
193  V /= magV;
194 
195 
196  // We now get the nearest intersections to start. This can either be
197  // the intersection with the end plane or with the cylinder side.
198 
199  // Get the two points (expressed in t) on the end planes. This is to
200  // clip any cylinder intersection against.
201  scalar tPoint1;
202  scalar tPoint2;
203 
204  // Maintain the two intersections with the endcaps
205  scalar tNear = VGREAT;
206  scalar tFar = VGREAT;
207 
208  {
209  scalar s = (V&unitDir_);
210  if (mag(s) > VSMALL)
211  {
212  tPoint1 = -s1/s;
213  tPoint2 = -(point2Start&unitDir_)/s;
214  if (tPoint2 < tPoint1)
215  {
216  Swap(tPoint1, tPoint2);
217  }
218  if (tPoint1 > magV || tPoint2 < 0)
219  {
220  return;
221  }
222 
223  // See if the points on the endcaps are actually inside the cylinder
224  if (tPoint1 >= 0 && tPoint1 <= magV)
225  {
226  if (radius2(start+tPoint1*V) <= sqr(radius_))
227  {
228  tNear = tPoint1;
229  }
230  }
231  if (tPoint2 >= 0 && tPoint2 <= magV)
232  {
233  if (radius2(start+tPoint2*V) <= sqr(radius_))
234  {
235  // Check if already have a near hit from point1
236  if (tNear <= 1)
237  {
238  tFar = tPoint2;
239  }
240  else
241  {
242  tNear = tPoint2;
243  }
244  }
245  }
246  }
247  else
248  {
249  // Vector perpendicular to cylinder. Check for outside already done
250  // above so just set tpoint to allow all.
251  tPoint1 = -VGREAT;
252  tPoint2 = VGREAT;
253  }
254  }
255 
256 
257  const vector x = point1Start ^ unitDir_;
258  const vector y = V ^ unitDir_;
259  const scalar d = sqr(radius_);
260 
261  // Second order equation of the form a*t^2 + b*t + c
262  const scalar a = (y&y);
263  const scalar b = 2*(x&y);
264  const scalar c = (x&x)-d;
265 
266  const scalar disc = b*b-4*a*c;
267 
268  scalar t1 = -VGREAT;
269  scalar t2 = VGREAT;
270 
271  if (disc < 0)
272  {
273  // Fully outside
274  return;
275  }
276  else if (disc < ROOTVSMALL)
277  {
278  // Single solution
279  if (mag(a) > ROOTVSMALL)
280  {
281  t1 = -b/(2*a);
282 
283  //Pout<< "single solution t:" << t1
284  // << " for start:" << start << " end:" << end
285  // << " c:" << c << endl;
286 
287  if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
288  {
289  // valid. Insert sorted.
290  if (t1 < tNear)
291  {
292  tFar = tNear;
293  tNear = t1;
294  }
295  else if (t1 < tFar)
296  {
297  tFar = t1;
298  }
299  }
300  else
301  {
302  return;
303  }
304  }
305  else
306  {
307  // Aligned with axis. Check if outside radius
308  //Pout<< "small discriminant:" << disc
309  // << " for start:" << start << " end:" << end
310  // << " magV:" << magV
311  // << " c:" << c << endl;
312  if (c > 0)
313  {
314  return;
315  }
316  }
317  }
318  else
319  {
320  if (mag(a) > ROOTVSMALL)
321  {
322  scalar sqrtDisc = sqrt(disc);
323 
324  t1 = (-b - sqrtDisc)/(2*a);
325  t2 = (-b + sqrtDisc)/(2*a);
326  if (t2 < t1)
327  {
328  Swap(t1, t2);
329  }
330 
331  if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
332  {
333  // valid. Insert sorted.
334  if (t1 < tNear)
335  {
336  tFar = tNear;
337  tNear = t1;
338  }
339  else if (t1 < tFar)
340  {
341  tFar = t1;
342  }
343  }
344  if (t2 >= 0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
345  {
346  // valid. Insert sorted.
347  if (t2 < tNear)
348  {
349  tFar = tNear;
350  tNear = t2;
351  }
352  else if (t2 < tFar)
353  {
354  tFar = t2;
355  }
356  }
357  //Pout<< "two solutions t1:" << t1 << " t2:" << t2
358  // << " for start:" << start << " end:" << end
359  // << " magV:" << magV
360  // << " c:" << c << endl;
361  }
362  else
363  {
364  // Aligned with axis. Check if outside radius
365  //Pout<< "large discriminant:" << disc
366  // << " small a:" << a
367  // << " for start:" << start << " end:" << end
368  // << " magV:" << magV
369  // << " c:" << c << endl;
370  if (c > 0)
371  {
372  return;
373  }
374  }
375  }
376 
377  // Check tNear, tFar
378  if (tNear >= 0 && tNear <= magV)
379  {
380  near.setPoint(start+tNear*V);
381  near.setHit();
382  near.setIndex(0);
383 
384  if (tFar <= magV)
385  {
386  far.setPoint(start+tFar*V);
387  far.setHit();
388  far.setIndex(0);
389  }
390  }
391  else if (tFar >= 0 && tFar <= magV)
392  {
393  near.setPoint(start+tFar*V);
394  near.setHit();
395  near.setIndex(0);
396  }
397 }
398 
399 
400 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
401 
402 Foam::searchableCylinder::searchableCylinder
403 (
404  const IOobject& io,
405  const point& point1,
406  const point& point2,
407  const scalar radius
408 )
409 :
410  searchableSurface(io),
411  point1_(point1),
412  point2_(point2),
413  magDir_(mag(point2_-point1_)),
414  unitDir_((point2_-point1_)/magDir_),
415  radius_(radius)
416 {}
417 
418 
419 Foam::searchableCylinder::searchableCylinder
420 (
421  const IOobject& io,
422  const dictionary& dict
423 )
424 :
425  searchableSurface(io),
426  point1_(dict.lookup("point1")),
427  point2_(dict.lookup("point2")),
428  magDir_(mag(point2_-point1_)),
429  unitDir_((point2_-point1_)/magDir_),
430  radius_(readScalar(dict.lookup("radius")))
431 {}
432 
433 
434 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
435 
437 {}
438 
439 
440 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
441 
443 {
444  if (regions_.empty())
445  {
446  regions_.setSize(1);
447  regions_[0] = "region0";
448  }
449  return regions_;
450 }
451 
452 
453 void Foam::searchableCylinder::findNearest
454 (
455  const pointField& samples,
456  const scalarField& nearestDistSqr,
457  List<pointIndexHit>& info
458 ) const
459 {
460  info.setSize(samples.size());
461 
462  forAll(samples, i)
463  {
464  info[i] = findNearest(samples[i], nearestDistSqr[i]);
465  }
466 }
467 
468 
470 (
471  const pointField& start,
472  const pointField& end,
473  List<pointIndexHit>& info
474 ) const
475 {
476  info.setSize(start.size());
477 
478  forAll(start, i)
479  {
480  // Pick nearest intersection. If none intersected take second one.
482  findLineAll(start[i], end[i], info[i], b);
483  if (!info[i].hit() && b.hit())
484  {
485  info[i] = b;
486  }
487  }
488 }
489 
490 
492 (
493  const pointField& start,
494  const pointField& end,
495  List<pointIndexHit>& info
496 ) const
497 {
498  info.setSize(start.size());
499 
500  forAll(start, i)
501  {
502  // Discard far intersection
504  findLineAll(start[i], end[i], info[i], b);
505  if (!info[i].hit() && b.hit())
506  {
507  info[i] = b;
508  }
509  }
510 }
511 
512 
513 void Foam::searchableCylinder::findLineAll
514 (
515  const pointField& start,
516  const pointField& end,
517  List<List<pointIndexHit> >& info
518 ) const
519 {
520  info.setSize(start.size());
521 
522  forAll(start, i)
523  {
524  pointIndexHit near, far;
525  findLineAll(start[i], end[i], near, far);
526 
527  if (near.hit())
528  {
529  if (far.hit())
530  {
531  info[i].setSize(2);
532  info[i][0] = near;
533  info[i][1] = far;
534  }
535  else
536  {
537  info[i].setSize(1);
538  info[i][0] = near;
539  }
540  }
541  else
542  {
543  if (far.hit())
544  {
545  info[i].setSize(1);
546  info[i][0] = far;
547  }
548  else
549  {
550  info[i].clear();
551  }
552  }
553  }
554 }
555 
556 
558 (
559  const List<pointIndexHit>& info,
560  labelList& region
561 ) const
562 {
563  region.setSize(info.size());
564  region = 0;
565 }
566 
567 
569 (
570  const List<pointIndexHit>& info,
572 ) const
573 {
574  normal.setSize(info.size());
575  normal = vector::zero;
576 
577  forAll(info, i)
578  {
579  if (info[i].hit())
580  {
581  vector v(info[i].hitPoint() - point1_);
582 
583  // Decompose sample-point1 into normal and parallel component
584  scalar parallel = v & unitDir_;
585 
586  if (parallel < 0)
587  {
588  normal[i] = -unitDir_;
589  }
590  else if (parallel > magDir_)
591  {
592  normal[i] = -unitDir_;
593  }
594  else
595  {
596  // Remove the parallel component
597  v -= parallel*unitDir_;
598  normal[i] = v/mag(v);
599  }
600  }
601  }
602 }
603 
604 
606 (
607  const pointField& points,
608  List<volumeType>& volType
609 ) const
610 {
611  volType.setSize(points.size());
612  volType = INSIDE;
613 
614  forAll(points, pointI)
615  {
616  const point& pt = points[pointI];
617 
618  vector v(pt - point1_);
619 
620  // Decompose sample-point1 into normal and parallel component
621  scalar parallel = v & unitDir_;
622 
623  if (parallel < 0)
624  {
625  // left of point1 endcap
626  volType[pointI] = OUTSIDE;
627  }
628  else if (parallel > magDir_)
629  {
630  // right of point2 endcap
631  volType[pointI] = OUTSIDE;
632  }
633  else
634  {
635  // Remove the parallel component
636  v -= parallel*unitDir_;
637 
638  if (mag(v) > radius_)
639  {
640  volType[pointI] = OUTSIDE;
641  }
642  else
643  {
644  volType[pointI] = INSIDE;
645  }
646  }
647  }
648 }
649 
650 
651 // ************************ vim: set sw=4 sts=4 et: ************************ //