1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "searchableCylinder.H"
28 #include "addToRunTimeSelectionTable.H"
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 defineTypeNameAndDebug(searchableCylinder, 0);
36 addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 Foam::pointField Foam::searchableCylinder::coordinates() const
45 pointField ctrs(1, 0.5*(point1_ + point2_));
51 Foam::pointIndexHit Foam::searchableCylinder::findNearest
54 const scalar nearestDistSqr
57 pointIndexHit info(false, sample, -1);
59 vector v(sample - point1_);
61 // Decompose sample-point1 into normal and parallel component
62 scalar parallel = (v & unitDir_);
64 // Remove the parallel component
65 v -= parallel*unitDir_;
71 // nearest is at point1 end cap. Clip to radius.
72 if (magV < ROOTVSMALL)
74 info.setPoint(point1_);
78 info.setPoint(point1_ + min(magV, radius_)*v/magV);
81 else if (parallel >= magDir_)
83 // nearest is at point2 end cap
84 if (magV < ROOTVSMALL)
86 info.setPoint(point2_);
90 info.setPoint(point2_ + min(magV, radius_)*v/magV);
95 // inbetween endcaps. Might either be nearer endcaps or cylinder wall
97 if (magV < ROOTVSMALL)
99 if (parallel < 0.5*magDir_)
101 info.setPoint(point1_);
105 info.setPoint(point2_);
112 // distance to endpoint: parallel or parallel-magDir
113 // distance to cylinder wall: magV-radius_
115 // Nearest cylinder point
116 point cylPt = sample + (radius_-magV)*n;
118 if (parallel < 0.5*magDir_)
120 // Project onto p1 endcap
121 point end1Pt = point1_ + min(magV, radius_)*n;
123 if (magSqr(sample-cylPt) < magSqr(sample-end1Pt))
125 info.setPoint(cylPt);
129 info.setPoint(end1Pt);
134 // Project onto p2 endcap
135 point end2Pt = point2_ + min(magV, radius_)*n;
137 if (magSqr(sample-cylPt) < magSqr(sample-end2Pt))
139 info.setPoint(cylPt);
143 info.setPoint(end2Pt);
149 if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
159 Foam::scalar Foam::searchableCylinder::radius2(const point& pt) const
161 const vector x = (pt-point1_) ^ unitDir_;
166 // From http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 -
167 // intersection of cylinder with ray
168 void Foam::searchableCylinder::findLineAll
179 vector point1Start(start-point1_);
180 vector point2Start(start-point2_);
181 vector point1End(end-point1_);
183 // Quick rejection of complete vector outside endcaps
184 scalar s1 = point1Start&unitDir_;
185 scalar s2 = point1End&unitDir_;
187 if ((s1 < 0 && s2 < 0) || (s1 > magDir_ && s2 > magDir_))
192 // Line as P = start+t*V where V is unit vector and t=[0..mag(end-start)]
194 scalar magV = mag(V);
195 if (magV < ROOTVSMALL)
202 // We now get the nearest intersections to start. This can either be
203 // the intersection with the end plane or with the cylinder side.
205 // Get the two points (expressed in t) on the end planes. This is to
206 // clip any cylinder intersection against.
210 // Maintain the two intersections with the endcaps
211 scalar tNear = VGREAT;
212 scalar tFar = VGREAT;
215 scalar s = (V&unitDir_);
219 tPoint2 = -(point2Start&unitDir_)/s;
220 if (tPoint2 < tPoint1)
222 Swap(tPoint1, tPoint2);
224 if (tPoint1 > magV || tPoint2 < 0)
229 // See if the points on the endcaps are actually inside the cylinder
230 if (tPoint1 >= 0 && tPoint1 <= magV)
232 if (radius2(start+tPoint1*V) <= sqr(radius_))
237 if (tPoint2 >= 0 && tPoint2 <= magV)
239 if (radius2(start+tPoint2*V) <= sqr(radius_))
241 // Check if already have a near hit from point1
255 // Vector perpendicular to cylinder. Check for outside already done
256 // above so just set tpoint to allow all.
263 const vector x = point1Start ^ unitDir_;
264 const vector y = V ^ unitDir_;
265 const scalar d = sqr(radius_);
267 // Second order equation of the form a*t^2 + b*t + c
268 const scalar a = (y&y);
269 const scalar b = 2*(x&y);
270 const scalar c = (x&x)-d;
272 const scalar disc = b*b-4*a*c;
282 else if (disc < ROOTVSMALL)
285 if (mag(a) > ROOTVSMALL)
289 //Pout<< "single solution t:" << t1
290 // << " for start:" << start << " end:" << end
291 // << " c:" << c << endl;
293 if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
295 // valid. Insert sorted.
313 // Aligned with axis. Check if outside radius
314 //Pout<< "small discriminant:" << disc
315 // << " for start:" << start << " end:" << end
316 // << " magV:" << magV
317 // << " c:" << c << endl;
326 if (mag(a) > ROOTVSMALL)
328 scalar sqrtDisc = sqrt(disc);
330 t1 = (-b - sqrtDisc)/(2*a);
331 t2 = (-b + sqrtDisc)/(2*a);
337 if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
339 // valid. Insert sorted.
350 if (t2 >= 0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
352 // valid. Insert sorted.
363 //Pout<< "two solutions t1:" << t1 << " t2:" << t2
364 // << " for start:" << start << " end:" << end
365 // << " magV:" << magV
366 // << " c:" << c << endl;
370 // Aligned with axis. Check if outside radius
371 //Pout<< "large discriminant:" << disc
372 // << " small a:" << a
373 // << " for start:" << start << " end:" << end
374 // << " magV:" << magV
375 // << " c:" << c << endl;
384 if (tNear >= 0 && tNear <= magV)
386 near.setPoint(start+tNear*V);
392 far.setPoint(start+tFar*V);
397 else if (tFar >= 0 && tFar <= magV)
399 near.setPoint(start+tFar*V);
406 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
408 Foam::searchableCylinder::searchableCylinder
416 searchableSurface(io),
419 magDir_(mag(point2_-point1_)),
420 unitDir_((point2_-point1_)/magDir_),
425 Foam::searchableCylinder::searchableCylinder
428 const dictionary& dict
431 searchableSurface(io),
432 point1_(dict.lookup("point1")),
433 point2_(dict.lookup("point2")),
434 magDir_(mag(point2_-point1_)),
435 unitDir_((point2_-point1_)/magDir_),
436 radius_(readScalar(dict.lookup("radius")))
440 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
442 Foam::searchableCylinder::~searchableCylinder()
446 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
448 const Foam::wordList& Foam::searchableCylinder::regions() const
450 if (regions_.empty())
453 regions_[0] = "region0";
459 void Foam::searchableCylinder::findNearest
461 const pointField& samples,
462 const scalarField& nearestDistSqr,
463 List<pointIndexHit>& info
466 info.setSize(samples.size());
470 info[i] = findNearest(samples[i], nearestDistSqr[i]);
475 void Foam::searchableCylinder::findLine
477 const pointField& start,
478 const pointField& end,
479 List<pointIndexHit>& info
482 info.setSize(start.size());
486 // Pick nearest intersection. If none intersected take second one.
488 findLineAll(start[i], end[i], info[i], b);
489 if (!info[i].hit() && b.hit())
497 void Foam::searchableCylinder::findLineAny
499 const pointField& start,
500 const pointField& end,
501 List<pointIndexHit>& info
504 info.setSize(start.size());
508 // Discard far intersection
510 findLineAll(start[i], end[i], info[i], b);
511 if (!info[i].hit() && b.hit())
519 void Foam::searchableCylinder::findLineAll
521 const pointField& start,
522 const pointField& end,
523 List<List<pointIndexHit> >& info
526 info.setSize(start.size());
530 pointIndexHit near, far;
531 findLineAll(start[i], end[i], near, far);
563 void Foam::searchableCylinder::getRegion
565 const List<pointIndexHit>& info,
569 region.setSize(info.size());
574 void Foam::searchableCylinder::getNormal
576 const List<pointIndexHit>& info,
580 normal.setSize(info.size());
581 normal = vector::zero;
587 vector v(info[i].hitPoint() - point1_);
589 // Decompose sample-point1 into normal and parallel component
590 scalar parallel = v & unitDir_;
594 normal[i] = -unitDir_;
596 else if (parallel > magDir_)
598 normal[i] = -unitDir_;
602 // Remove the parallel component
603 v -= parallel*unitDir_;
604 normal[i] = v/mag(v);
611 void Foam::searchableCylinder::getVolumeType
613 const pointField& points,
614 List<volumeType>& volType
617 volType.setSize(points.size());
620 forAll(points, pointI)
622 const point& pt = points[pointI];
624 vector v(pt - point1_);
626 // Decompose sample-point1 into normal and parallel component
627 scalar parallel = v & unitDir_;
631 // left of point1 endcap
632 volType[pointI] = OUTSIDE;
634 else if (parallel > magDir_)
636 // right of point2 endcap
637 volType[pointI] = OUTSIDE;
641 // Remove the parallel component
642 v -= parallel*unitDir_;
644 if (mag(v) > radius_)
646 volType[pointI] = OUTSIDE;
650 volType[pointI] = INSIDE;
657 // ************************************************************************* //