53 const scalar nearestDistSqr
58 vector v(sample - point1_);
61 scalar parallel = (v & unitDir_);
64 v -= parallel*unitDir_;
67 if (magV < ROOTVSMALL)
79 info.setPoint(point1_ +
min(magV, radius_)*v);
81 else if (parallel >= magDir_)
84 info.setPoint(point2_ +
min(magV, radius_)*v);
95 if (magV < ROOTVSMALL)
99 scalar magE1 =
mag(e1);
102 e1 =
point(0,1,0) ^ unitDir_;
106 cylPt = sample + radius_*e1;
110 cylPt = sample + (radius_-magV)*v;
113 if (parallel < 0.5*magDir_)
116 point end1Pt = point1_ +
min(magV, radius_)*v;
120 info.setPoint(cylPt);
124 info.setPoint(end1Pt);
130 point end2Pt = point2_ +
min(magV, radius_)*v;
134 info.setPoint(cylPt);
138 info.setPoint(end2Pt);
143 if (
magSqr(sample - info.rawPoint()) < nearestDistSqr)
153 Foam::scalar Foam::searchableCylinder::radius2(
const point& pt)
const
155 const vector x = (pt-point1_) ^ unitDir_;
162 void Foam::searchableCylinder::findLineAll
173 vector point1Start(start-point1_);
174 vector point2Start(start-point2_);
175 vector point1End(end-point1_);
178 scalar s1 = point1Start&unitDir_;
179 scalar s2 = point1End&unitDir_;
181 if ((s1 < 0 && s2 < 0) || (s1 > magDir_ && s2 > magDir_))
188 scalar magV =
mag(V);
189 if (magV < ROOTVSMALL)
205 scalar tNear = VGREAT;
206 scalar tFar = VGREAT;
209 scalar s = (V&unitDir_);
213 tPoint2 = -(point2Start&unitDir_)/s;
214 if (tPoint2 < tPoint1)
216 Swap(tPoint1, tPoint2);
218 if (tPoint1 > magV || tPoint2 < 0)
224 if (tPoint1 >= 0 && tPoint1 <= magV)
226 if (radius2(start+tPoint1*V) <=
sqr(radius_))
231 if (tPoint2 >= 0 && tPoint2 <= magV)
233 if (radius2(start+tPoint2*V) <=
sqr(radius_))
257 const vector x = point1Start ^ unitDir_;
259 const scalar
d =
sqr(radius_);
262 const scalar a = (y&
y);
263 const scalar
b = 2*(x&
y);
264 const scalar c = (x&x)-d;
266 const scalar disc = b*b-4*a*c;
276 else if (disc < ROOTVSMALL)
279 if (
mag(a) > ROOTVSMALL)
287 if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
320 if (
mag(a) > ROOTVSMALL)
322 scalar sqrtDisc =
sqrt(disc);
324 t1 = (-b - sqrtDisc)/(2*a);
325 t2 = (-b + sqrtDisc)/(2*a);
331 if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
344 if (t2 >= 0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
378 if (tNear >= 0 && tNear <= magV)
380 near.setPoint(start+tNear*V);
386 far.setPoint(start+tFar*V);
391 else if (tFar >= 0 && tFar <= magV)
393 near.setPoint(start+tFar*V);
402 Foam::searchableCylinder::searchableCylinder
410 searchableSurface(io),
413 magDir_(
mag(point2_-point1_)),
414 unitDir_((point2_-point1_)/magDir_),
419 Foam::searchableCylinder::searchableCylinder
426 point1_(dict.
lookup(
"point1")),
427 point2_(dict.
lookup(
"point2")),
428 magDir_(
mag(point2_-point1_)),
429 unitDir_((point2_-point1_)/magDir_),
444 if (regions_.empty())
447 regions_[0] =
"region0";
453 void Foam::searchableCylinder::findNearest
464 info[i] = findNearest(samples[i], nearestDistSqr[i]);
482 findLineAll(start[i], end[i], info[i], b);
483 if (!info[i].hit() && b.
hit())
504 findLineAll(start[i], end[i], info[i], b);
505 if (!info[i].hit() && b.
hit())
513 void Foam::searchableCylinder::findLineAll
525 findLineAll(start[i], end[i], near, far);
581 vector v(info[i].hitPoint() - point1_);
584 scalar parallel = v & unitDir_;
588 normal[i] = -unitDir_;
590 else if (parallel > magDir_)
592 normal[i] = -unitDir_;
597 v -= parallel*unitDir_;
598 normal[i] = v/
mag(v);
616 const point& pt = points[pointI];
621 scalar parallel = v & unitDir_;
626 volType[pointI] = OUTSIDE;
628 else if (parallel > magDir_)
631 volType[pointI] = OUTSIDE;
636 v -= parallel*unitDir_;
638 if (
mag(v) > radius_)
640 volType[pointI] = OUTSIDE;
644 volType[pointI] = INSIDE;