Go to the documentation of this file.
87 if (
p.x() < 0) { octant |= 1; }
88 if (
p.y() < 0) { octant |= 2; }
89 if (
p.z() < 0) { octant |= 4; }
98 if (octant & 1) {
p.x() = -
p.x(); }
99 if (octant & 2) {
p.y() = -
p.y(); }
100 if (octant & 4) {
p.z() = -
p.z(); }
170 const scalar n0 = r0*z0;
173 scalar s1 = (
g < 0 ? 0 :
vectorMag(n0, z1) - 1);
203 <<
"Located root in " << nIters <<
" iterations" <<
endl;
221 const scalar n0 = r0*z0;
222 const scalar n1 = r1*z1;
225 scalar s1 = (
g < 0 ? 0 :
vectorMag(n0, n1, z2) - 1);
255 <<
"root at " <<
s <<
" found in " << nIters
256 <<
" iterations" <<
endl;
267 const scalar e0,
const scalar e1,
269 const scalar
y0,
const scalar
y1,
271 scalar& x0, scalar& x1
278 const scalar numer0 = e0*
y0;
279 const scalar denom0 =
sqr(e0) -
sqr(e1);
283 const scalar xde0 = numer0/denom0;
308 const scalar z0 =
y0 / e0;
309 const scalar z1 =
y1 / e1;
310 scalar eval =
sqr(z0) +
sqr(z1);
335 const scalar r0 =
sqr(e0 / e1);
340 x0 = r0 *
y0 / (sbar + r0);
341 x1 =
y1 / (sbar + 1);
344 eval =
sqr(x0/e0) +
sqr(x1/e1);
364 <<
"Programming/logic error" <<
nl
375 const scalar e0,
const scalar e1,
const scalar e2,
377 const scalar
y0,
const scalar
y1,
const scalar y2,
379 scalar& x0, scalar& x1, scalar& x2
386 const scalar numer0 = e0*
y0;
387 const scalar numer1 = e1*
y1;
388 const scalar denom0 =
sqr(e0) -
sqr(e2);
389 const scalar denom1 =
sqr(e1) -
sqr(e2);
391 if (numer0 < denom0 && numer1 < denom1)
393 const scalar xde0 = numer0/denom0;
394 const scalar xde1 = numer1/denom1;
396 const scalar disc = (1 -
sqr(xde0) -
sqr(xde1));
437 const scalar z0 =
y0/e0;
438 const scalar z1 =
y1/e1;
439 const scalar z2 = y2/e2;
469 const scalar r0 =
sqr(e0/e2);
470 const scalar r1 =
sqr(e1/e2);
475 x0 = r0*
y0/(sbar+r0);
476 x1 = r1*
y1/(sbar+r1);
501 <<
"Programming/logic error" <<
nl
512 inline Foam::searchableSphere::componentOrder
513 Foam::searchableSphere::getOrdering(
const vector& radii)
518 if (
radii[cmpt] <= 0)
521 <<
"Radii must be positive, non-zero: " <<
radii <<
endl
527 std::array<uint8_t, 3> idx{0, 1, 2};
535 [&](uint8_t a, uint8_t
b){
return radii[a] >
radii[
b]; }
538 componentOrder order{idx[0], idx[1], idx[2], shapeType::GENERAL};
542 order.shape = shapeType::SPHERE;
546 order.shape = shapeType::OBLATE;
550 order.shape = shapeType::PROLATE;
562 const scalar nearestDistSqr
569 if (order_.shape == shapeType::SPHERE)
572 const vector n(sample - origin_);
573 const scalar magN =
mag(
n);
577 if (nearestDistSqr >=
sqr(magN - radii_[0]))
582 + (magN < ROOTVSMALL ?
vector(radii_[0],0,0) : (radii_[0]*
n/magN))
597 vector relPt(sample - origin_);
600 const unsigned octant =
getOctant(relPt);
609 point& near = info.rawPoint();
612 if (order_.shape == shapeType::OBLATE)
616 const scalar axisDist =
hypot(relPt[order_.major], relPt[order_.mezzo]);
623 radii_[order_.major], radii_[order_.minor],
624 axisDist, relPt[order_.minor],
625 nearAxis, near[order_.minor]
629 nearAxis /= (axisDist + VSMALL);
631 near[order_.major] = relPt[order_.major] * nearAxis;
632 near[order_.mezzo] = relPt[order_.mezzo] * nearAxis;
635 else if (order_.shape == shapeType::PROLATE)
639 const scalar axisDist =
hypot(relPt[order_.mezzo], relPt[order_.minor]);
646 radii_[order_.major], radii_[order_.minor],
647 relPt[order_.major], axisDist,
648 near[order_.major], nearAxis
652 nearAxis /= (axisDist + VSMALL);
655 near[order_.mezzo] = relPt[order_.mezzo] * nearAxis;
656 near[order_.minor] = relPt[order_.minor] * nearAxis;
662 radii_[order_.major], radii_[order_.mezzo], radii_[order_.minor],
663 relPt[order_.major], relPt[order_.mezzo], relPt[order_.minor],
664 near[order_.major], near[order_.mezzo], near[order_.minor]
676 if (distSqr <= nearestDistSqr)
686 void Foam::searchableSphere::findLineAll
697 if (order_.shape == shapeType::SPHERE)
700 const scalar magSqrDir =
magSqr(dir);
702 if (magSqrDir > ROOTVSMALL)
706 const vector relStart(start - origin_);
708 const scalar v = -(relStart & dir);
710 const scalar disc =
sqr(radius()) - (
magSqr(relStart) -
sqr(v));
716 const scalar nearParam = v - d;
717 const scalar farParam = v + d;
719 if (nearParam >= 0 &&
sqr(nearParam) <= magSqrDir)
721 near.hitPoint(start + nearParam*dir, 0);
724 if (farParam >= 0 &&
sqr(farParam) <= magSqrDir)
726 far.hitPoint(start + farParam*dir, 0);
744 const point relStart = scalePoint(start);
747 const scalar magSqrDir =
magSqr(dir);
749 if (magSqrDir > ROOTVSMALL)
753 const scalar v = -(relStart & dir);
755 const scalar disc = scalar(1) - (
magSqr(relStart) -
sqr(v));
761 const scalar nearParam = v - d;
762 const scalar farParam = v + d;
764 if (nearParam >= 0 &&
sqr(nearParam) <= magSqrDir)
766 near.hitPoint(unscalePoint(relStart + nearParam*dir), 0);
768 if (farParam >= 0 &&
sqr(farParam) <= magSqrDir)
770 far.hitPoint(unscalePoint(relStart + farParam*dir), 0);
779 Foam::searchableSphere::searchableSphere
790 Foam::searchableSphere::searchableSphere
800 order_(getOrdering(radii_))
802 bounds().min() = (centre() - radii_);
803 bounds().max() = (centre() + radii_);
807 Foam::searchableSphere::searchableSphere
832 origin_.x() + radii_.x() *
cos(theta)*
sin(
phi),
833 origin_.y() + radii_.y() *
sin(theta)*
sin(
phi),
834 origin_.z() + radii_.z() *
cos(
phi)
858 if (order_.shape == shapeType::SPHERE)
880 const scalar d0 = bb.
min()[dir] - origin_[dir];
881 const scalar d1 = bb.
max()[dir] - origin_[dir];
883 if ((d0 > 0) == (d1 > 0))
910 if (regions_.empty())
913 regions_.first() =
"region0";
928 centres[0] = origin_;
936 void Foam::searchableSphere::findNearest
947 info[i] = findNearest(
samples[i], nearestDistSqr[i]);
959 info.
resize(start.size());
968 findLineAll(start[i],
end[i], info[i],
b);
970 if (!info[i].hit() &&
b.hit())
985 info.
resize(start.size());
994 findLineAll(start[i],
end[i], info[i],
b);
996 if (!info[i].hit() &&
b.hit())
1004 void Foam::searchableSphere::findLineAll
1011 info.
resize(start.size());
1017 findLineAll(start[i],
end[i], near, far);
1055 region.
resize(info.size());
1066 normal.resize(info.size());
1072 if (order_.shape == shapeType::SPHERE)
1075 normal[i] =
normalised(info[i].hitPoint() - origin_);
1082 normal[i] = scalePoint(info[i].hitPoint());
1084 normal[i].x() /= radii_.x();
1085 normal[i].y() /= radii_.y();
1086 normal[i].z() /= radii_.z();
1087 normal[i].normalise();
1107 if (order_.shape == shapeType::SPHERE)
1111 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< Cmpt > uniform(const Cmpt &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))
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.
word name(const complex &c)
Return string representation of complex.
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.
void resize(const label newSize)
Adjust allocated size of list.
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.
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
dimensionedScalar cos(const dimensionedScalar &ds)
Searching on general spheroid.