55 centres[0] = 0.5*(point1_ + point2_);
80 const scalar nearestDistSqr
85 vector v(sample - point1_);
88 scalar parallel = (v & unitDir_);
91 v -= parallel*unitDir_;
94 if (magV < ROOTVSMALL)
108 else if (parallel >= magDir_)
122 if (magV < ROOTVSMALL)
126 scalar magE1 =
mag(e1);
129 e1 =
point(0,1,0) ^ unitDir_;
133 cylPt = sample + radius_*e1;
137 cylPt = sample + (radius_-magV)*v;
140 if (parallel < 0.5*magDir_)
143 point end1Pt = point1_ +
min(magV, radius_)*v;
157 point end2Pt = point2_ +
min(magV, radius_)*v;
180 Foam::scalar Foam::searchableCylinder::radius2(
const point& pt)
const 182 const vector x = (pt-point1_) ^ unitDir_;
189 void Foam::searchableCylinder::findLineAll
200 vector point1Start(start-point1_);
201 vector point2Start(start-point2_);
202 vector point1End(end-point1_);
205 scalar s1 = point1Start&unitDir_;
206 scalar s2 = point1End&unitDir_;
208 if ((s1 < 0 && s2 < 0) || (s1 > magDir_ && s2 > magDir_))
215 scalar magV =
mag(V);
216 if (magV < ROOTVSMALL)
232 scalar tNear = VGREAT;
233 scalar tFar = VGREAT;
236 scalar
s = (V&unitDir_);
240 tPoint2 = -(point2Start&unitDir_)/s;
241 if (tPoint2 < tPoint1)
243 Swap(tPoint1, tPoint2);
245 if (tPoint1 > magV || tPoint2 < 0)
251 if (tPoint1 >= 0 && tPoint1 <= magV)
253 if (radius2(start+tPoint1*V) <=
sqr(radius_))
258 if (tPoint2 >= 0 && tPoint2 <= magV)
260 if (radius2(start+tPoint2*V) <=
sqr(radius_))
284 const vector x = point1Start ^ unitDir_;
286 const scalar d =
sqr(radius_);
289 const scalar a = (y&
y);
290 const scalar
b = 2*(x&
y);
291 const scalar
c = (x&
x)-d;
293 const scalar disc = b*b-4*a*
c;
303 else if (disc < ROOTVSMALL)
306 if (
mag(a) > ROOTVSMALL)
314 if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
347 if (
mag(a) > ROOTVSMALL)
349 scalar sqrtDisc =
sqrt(disc);
351 t1 = (-b - sqrtDisc)/(2*a);
352 t2 = (-b + sqrtDisc)/(2*a);
358 if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
371 if (t2 >= 0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
405 if (tNear >= 0 && tNear <= magV)
418 else if (tFar >= 0 && tFar <= magV)
477 Foam::searchableCylinder::searchableCylinder
488 magDir_(
mag(point2_-point1_)),
489 unitDir_((point2_-point1_)/magDir_),
492 bounds() = calcBounds();
496 Foam::searchableCylinder::searchableCylinder
503 point1_(dict.
lookup(
"point1")),
504 point2_(dict.
lookup(
"point2")),
505 magDir_(
mag(point2_-point1_)),
506 unitDir_((point2_-point1_)/magDir_),
509 bounds() = calcBounds();
523 if (regions_.empty())
526 regions_[0] =
"region0";
532 void Foam::searchableCylinder::findNearest
543 info[i] = findNearest(samples[i], nearestDistSqr[i]);
561 findLineAll(start[i], end[i], info[i], b);
562 if (!info[i].hit() && b.
hit())
583 findLineAll(start[i], end[i], info[i], b);
584 if (!info[i].hit() && b.
hit())
592 void Foam::searchableCylinder::findLineAll
604 findLineAll(start[i], end[i], near, far);
660 vector v(info[i].hitPoint() - point1_);
663 scalar parallel = (v & unitDir_);
666 v -= parallel*unitDir_;
667 scalar magV =
mag(v);
671 if ((magV-radius_) <
mag(parallel))
674 normal[i] = -unitDir_;
681 else if (parallel <= 0.5*magDir_)
684 if (magV >= radius_ || (radius_-magV) < parallel)
691 normal[i] = -unitDir_;
694 else if (parallel <= magDir_)
697 if (magV >= radius_ || (radius_-magV) < (magDir_-parallel))
704 normal[i] = unitDir_;
709 if ((magV-radius_) < (parallel-magDir_))
712 normal[i] = unitDir_;
735 const point& pt = points[pointi];
740 scalar parallel = v & unitDir_;
747 else if (parallel > magDir_)
755 v -= parallel*unitDir_;
757 if (
mag(v) > radius_)
#define forAll(list, i)
Loop across all elements in list.
virtual void boundingSpheres(pointField ¢res, 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...
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.
void setIndex(const label index)
dimensionedScalar sqrt(const dimensionedScalar &ds)
A bounding box defined in terms of the points at its extremities.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
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.
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.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
void clear()
Clear the list, i.e. set size to zero.
virtual ~searchableCylinder()
Destructor.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
const Point & rawPoint() const
Return point with no checking.
virtual void getRegion(const List< pointIndexHit > &, labelList ®ion) 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.
vector point
Point is a vector.
const dimensionedScalar c
Speed of light in a vacuum.
dimensioned< scalar > mag(const dimensioned< Type > &)
A class for managing temporary objects.
T & ref() const
Return non-const reference or generate a fatal error.
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...
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.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
virtual void getVolumeType(const pointField &, List< volumeType > &) const
Determine type (inside/outside/mixed) for point. unknown if.