Go to the documentation of this file.
88 if (
p.x() < 0) { octant |= 1; }
89 if (
p.y() < 0) { octant |= 2; }
90 if (
p.z() < 0) { octant |= 4; }
99 if (octant & 1) {
p.x() = -
p.x(); }
100 if (octant & 2) {
p.y() = -
p.y(); }
101 if (octant & 4) {
p.z() = -
p.z(); }
171 const scalar n0 = r0*z0;
174 scalar s1 = (
g < 0 ? 0 :
vectorMag(n0, z1) - 1);
204 <<
"Located root in " << nIters <<
" iterations" <<
endl;
222 const scalar n0 = r0*z0;
223 const scalar n1 = r1*z1;
226 scalar s1 = (
g < 0 ? 0 :
vectorMag(n0, n1, z2) - 1);
256 <<
"root at " <<
s <<
" found in " << nIters
257 <<
" iterations" <<
endl;
268 const scalar e0,
const scalar e1,
270 const scalar
y0,
const scalar
y1,
272 scalar& x0, scalar& x1
279 const scalar numer0 = e0*
y0;
280 const scalar denom0 =
sqr(e0) -
sqr(e1);
284 const scalar xde0 = numer0/denom0;
309 const scalar z0 =
y0 / e0;
310 const scalar z1 =
y1 / e1;
311 scalar eval =
sqr(z0) +
sqr(z1);
336 const scalar r0 =
sqr(e0 / e1);
341 x0 = r0 *
y0 / (sbar + r0);
342 x1 =
y1 / (sbar + 1);
345 eval =
sqr(x0/e0) +
sqr(x1/e1);
365 <<
"Programming/logic error" <<
nl
376 const scalar e0,
const scalar e1,
const scalar e2,
378 const scalar
y0,
const scalar
y1,
const scalar y2,
380 scalar& x0, scalar& x1, scalar& x2
387 const scalar numer0 = e0*
y0;
388 const scalar numer1 = e1*
y1;
389 const scalar denom0 =
sqr(e0) -
sqr(e2);
390 const scalar denom1 =
sqr(e1) -
sqr(e2);
392 if (numer0 < denom0 && numer1 < denom1)
394 const scalar xde0 = numer0/denom0;
395 const scalar xde1 = numer1/denom1;
397 const scalar disc = (1 -
sqr(xde0) -
sqr(xde1));
438 const scalar z0 =
y0/e0;
439 const scalar z1 =
y1/e1;
440 const scalar z2 = y2/e2;
470 const scalar r0 =
sqr(e0/e2);
471 const scalar r1 =
sqr(e1/e2);
476 x0 = r0*
y0/(sbar+r0);
477 x1 = r1*
y1/(sbar+r1);
502 <<
"Programming/logic error" <<
nl
513 inline Foam::searchableSphere::componentOrder
514 Foam::searchableSphere::getOrdering(
const vector& radii)
519 if (
radii[cmpt] <= 0)
522 <<
"Radii must be positive, non-zero: " <<
radii <<
endl
528 std::array<uint8_t, 3> idx{0, 1, 2};
536 [&](uint8_t a, uint8_t
b){
return radii[a] >
radii[
b]; }
539 componentOrder order{idx[0], idx[1], idx[2], shapeType::GENERAL};
543 order.shape = shapeType::SPHERE;
547 order.shape = shapeType::OBLATE;
551 order.shape = shapeType::PROLATE;
563 const scalar nearestDistSqr
570 if (order_.shape == shapeType::SPHERE)
574 const scalar magN =
mag(
n);
578 if (nearestDistSqr >=
sqr(magN - radii_[0]))
583 + (magN < ROOTVSMALL ?
vector(radii_[0],0,0) : (radii_[0]*
n/magN))
601 const unsigned octant =
getOctant(relPt);
610 point& near = info.rawPoint();
613 if (order_.shape == shapeType::OBLATE)
617 const scalar axisDist =
hypot(relPt[order_.major], relPt[order_.mezzo]);
624 radii_[order_.major], radii_[order_.minor],
625 axisDist, relPt[order_.minor],
626 nearAxis, near[order_.minor]
630 nearAxis /= (axisDist + VSMALL);
632 near[order_.major] = relPt[order_.major] * nearAxis;
633 near[order_.mezzo] = relPt[order_.mezzo] * nearAxis;
636 else if (order_.shape == shapeType::PROLATE)
640 const scalar axisDist =
hypot(relPt[order_.mezzo], relPt[order_.minor]);
647 radii_[order_.major], radii_[order_.minor],
648 relPt[order_.major], axisDist,
649 near[order_.major], nearAxis
653 nearAxis /= (axisDist + VSMALL);
656 near[order_.mezzo] = relPt[order_.mezzo] * nearAxis;
657 near[order_.minor] = relPt[order_.minor] * nearAxis;
663 radii_[order_.major], radii_[order_.mezzo], radii_[order_.minor],
664 relPt[order_.major], relPt[order_.mezzo], relPt[order_.minor],
665 near[order_.major], near[order_.mezzo], near[order_.minor]
677 if (distSqr <= nearestDistSqr)
687 void Foam::searchableSphere::findLineAll
698 if (order_.shape == shapeType::SPHERE)
701 const scalar magSqrDir =
magSqr(dir);
703 if (magSqrDir > ROOTVSMALL)
707 const vector relStart(start - origin_);
709 const scalar v = -(relStart & dir);
711 const scalar disc =
sqr(radius()) - (
magSqr(relStart) -
sqr(v));
717 const scalar nearParam = v - d;
718 const scalar farParam = v + d;
720 if (nearParam >= 0 &&
sqr(nearParam) <= magSqrDir)
722 near.hitPoint(start + nearParam*dir, 0);
725 if (farParam >= 0 &&
sqr(farParam) <= magSqrDir)
727 far.hitPoint(start + farParam*dir, 0);
745 const point relStart = scalePoint(start);
748 const scalar magSqrDir =
magSqr(dir);
750 if (magSqrDir > ROOTVSMALL)
754 const scalar v = -(relStart & dir);
756 const scalar disc = scalar(1) - (
magSqr(relStart) -
sqr(v));
762 const scalar nearParam = v - d;
763 const scalar farParam = v + d;
765 if (nearParam >= 0 &&
sqr(nearParam) <= magSqrDir)
767 near.hitPoint(unscalePoint(relStart + nearParam*dir), 0);
769 if (farParam >= 0 &&
sqr(farParam) <= magSqrDir)
771 far.hitPoint(unscalePoint(relStart + farParam*dir), 0);
780 Foam::searchableSphere::searchableSphere
791 Foam::searchableSphere::searchableSphere
801 order_(getOrdering(radii_))
803 bounds().min() = (centre() - radii_);
804 bounds().max() = (centre() + radii_);
808 Foam::searchableSphere::searchableSphere
833 origin_.x() + radii_.x() *
cos(theta)*
sin(
phi),
834 origin_.y() + radii_.y() *
sin(theta)*
sin(
phi),
835 origin_.z() + radii_.z() *
cos(
phi)
859 if (order_.shape == shapeType::SPHERE)
881 const scalar d0 = bb.
min()[dir] - origin_[dir];
882 const scalar d1 = bb.
max()[dir] - origin_[dir];
884 if ((d0 > 0) == (d1 > 0))
911 if (regions_.empty())
914 regions_.first() =
"region0";
929 centres[0] = origin_;
937 void Foam::searchableSphere::findNearest
948 info[i] = findNearest(
samples[i], nearestDistSqr[i]);
960 info.
resize(start.size());
969 findLineAll(start[i],
end[i], info[i],
b);
971 if (!info[i].hit() &&
b.hit())
986 info.
resize(start.size());
995 findLineAll(start[i],
end[i], info[i],
b);
997 if (!info[i].hit() &&
b.hit())
1005 void Foam::searchableSphere::findLineAll
1012 info.
resize(start.size());
1018 findLineAll(start[i],
end[i], near, far);
1056 region.
resize(info.size());
1067 normal.resize(info.size());
1073 if (order_.shape == shapeType::SPHERE)
1076 normal[i] =
normalised(info[i].hitPoint() - origin_);
1083 normal[i] = scalePoint(info[i].hitPoint());
1085 normal[i].x() /= radii_.x();
1086 normal[i].y() /= radii_.y();
1087 normal[i].z() /= radii_.z();
1088 normal[i].normalise();
1108 if (order_.shape == shapeType::SPHERE)
1112 const scalar rad2 =
sqr(radius());
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
Defines the attributes of an object for which implicit objectRegistry management is supported,...
#define InfoInFunction
Report an information message using Foam::Info.
A class for handling words, derived from Foam::string.
static Vector< scalar > uniform(const scalar &s)
Return a VectorSpace with all elements = s.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
void resize(const label len)
Adjust allocated size of list.
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find first intersection on segment from start to end.
dimensionedScalar y1(const dimensionedScalar &ds)
static constexpr const zero Zero
Global zero (0)
virtual void getRegion(const List< pointIndexHit > &, labelList ®ion) const
From a set of points and indices get the region.
static constexpr int maxIters
dimensionedScalar sin(const dimensionedScalar &ds)
static void applyOctant(point &p, unsigned octant)
virtual const wordList & regions() const
Names of regions.
Ostream & endl(Ostream &os)
Add newline and flush stream.
A token holds an item read from Istream.
static scalar vectorMagSqr(const scalar x, const scalar y)
const point & max() const
Maximum describing the bounding box.
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
T getCompat(const word &keyword, std::initializer_list< std::pair< const char *, int >> compat, enum keyType::option=keyType::REGEX) const
static scalar distanceToEllipse(const scalar e0, const scalar e1, const scalar y0, const scalar y1, scalar &x0, scalar &x1)
static vector getRadius(const word &name, const dictionary &dict)
point surfacePoint(const scalar theta, const scalar phi) const
A point on the sphere at given location.
const point & min() const
Minimum describing the bounding box.
#define forAll(list, i)
Loop across all elements in list.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
bool overlaps(const boundBox &bb) const
Overlaps/touches boundingBox?
Specialization of rigidBody to construct a sphere given the mass and radius.
const vector & radii() const noexcept
The radii of the spheroid.
static unsigned getOctant(const point &p)
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
bool hit() const noexcept
Is there a hit?
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
static scalar vectorMag(const scalar x, const scalar y)
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
addNamedToRunTimeSelectionTable(topoSetCellSource, badQualityToCell, word, badQuality)
vector surfaceNormal(const scalar theta, const scalar phi) const
Surface normal on the sphere at given location.
dimensionedScalar y0(const dimensionedScalar &ds)
scalarField samples(nIntervals, Zero)
virtual void boundingSpheres(pointField ¢res, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Return any intersection on segment from start to end.
ITstream & lookup(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Macros for easy insertion into run-time selection tables.
virtual void getVolumeType(const pointField &points, List< volumeType > &volType) const
Determine type (inside/outside/mixed) for point.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
const uniformDimensionedVectorField & g
Vector< scalar > vector
A scalar version of the templated Vector.
errorManipArg< error, int > exit(error &err, const int errNo=1)
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
static scalar distanceToEllipsoid(const scalar e0, const scalar e1, const scalar e2, const scalar y0, const scalar y1, const scalar y2, scalar &x0, scalar &x1, scalar &x2)
static scalar findRootEllipseDistance(const scalar r0, const scalar z0, const scalar z1, scalar g)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
const dimensionedScalar e
Elementary charge.
void clear()
Clear the list, i.e. set size to zero.
static scalar findRootEllipsoidDistance(const scalar r0, const scalar r1, const scalar z0, const scalar z1, const scalar z2, scalar g)
A location inside the volume.
A bounding box defined in terms of min/max extrema points.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
vector point
Point is a vector.
bool valid() const
Bounding box is non-inverted.
virtual bool overlaps(const boundBox &bb) const
Does any part of the surface overlap the supplied bound box?
bool equal(const T &s1, const T &s2)
Compare two values for equality.
defineTypeNameAndDebug(combustionModel, 0)
static constexpr direction nComponents
Number of components in this vector space.
A location outside the volume.
static constexpr scalar tolCloseness
Minimal example by using system/controlDict.functions:
dimensionedScalar cos(const dimensionedScalar &ds)
Searching on general spheroid.