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
513inline Foam::searchableSphere::componentOrder
514Foam::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]; }
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)
687void 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);
747 vector dir(scalePoint(end) - relStart);
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);
801 order_(getOrdering(radii_))
817 dict.getCompat<
vector>(
"origin", {{
"centre", -1806}}),
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_;
937void Foam::searchableSphere::findNearest
948 info[i] = findNearest(
samples[i], nearestDistSqr[i]);
969 findLineAll(start[i], end[i], info[i],
b);
971 if (!info[i].hit() &&
b.hit())
995 findLineAll(start[i], end[i], info[i],
b);
997 if (!info[i].hit() &&
b.hit())
1005void Foam::searchableSphere::findLineAll
1012 info.resize(start.
size());
1018 findLineAll(start[i], end[i], near, far);
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();
1108 if (order_.shape == shapeType::SPHERE)
1112 const scalar rad2 =
sqr(radius());
Macros for easy insertion into run-time selection tables.
#define addNamedToRunTimeSelectionTable(baseType, thisType, argNames, lookupName)
Add to construction table with 'lookupName' as the key.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
const uniformDimensionedVectorField & g
bool surfacePoint() const
Part of a surface point pair.
Minimal example by using system/controlDict.functions:
Defines the attributes of an object for which implicit objectRegistry management is supported,...
void resize(const label len)
Adjust allocated size of list.
int overlaps
Flag to control which overlap calculations are performed.
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?
void size(const label n)
Older name for setAddressableSize.
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
A bounding box defined in terms of min/max extrema points.
bool valid() const
Bounding box is non-inverted.
bool overlaps(const boundBox &bb) const
Overlaps/touches boundingBox?
const point & min() const
Minimum describing the bounding box.
const point & max() const
Maximum describing the bounding box.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
ITstream & lookup(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
static constexpr direction nComponents
Number of components in bool is 1.
Searching on general spheroid.
const vector & radii() const noexcept
The radii of the spheroid.
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find first intersection on segment from start to end.
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Return any intersection on segment from start to end.
virtual void boundingSpheres(pointField ¢res, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
virtual void getVolumeType(const pointField &points, List< volumeType > &volType) const
Determine type (inside/outside/mixed) for point.
vector surfaceNormal(const scalar theta, const scalar phi) const
Surface normal on the sphere at given location.
const point & centre() const noexcept
The centre (origin) of the sphere.
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
virtual void getRegion(const List< pointIndexHit > &, labelList ®ion) const
From a set of points and indices get the region.
virtual const wordList & regions() const
Names of regions.
@ PROLATE
Prolate (major > mezzo = minor)
@ GENERAL
General spheroid.
@ SPHERE
Sphere (all components equal)
@ OBLATE
Oblate (major = mezzo > minor)
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
virtual const boundBox & bounds() const
Return const reference to boundBox.
Specialization of rigidBody to construct a sphere given the mass and radius.
A token holds an item read from Istream.
bool isNumber() const noexcept
Token is LABEL, FLOAT or DOUBLE.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
@ OUTSIDE
A location outside the volume.
@ INSIDE
A location inside the volume.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
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))
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
#define InfoInFunction
Report an information message using Foam::Info.
static unsigned getOctant(const point &p)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
dimensionedScalar y0(const dimensionedScalar &ds)
static vector getRadius(const word &name, const dictionary &dict)
dimensionedScalar sin(const dimensionedScalar &ds)
static scalar vectorMag(const scalar x, const scalar y)
vector point
Point is a vector.
static constexpr scalar tolCloseness
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool equal(const T &s1, const T &s2)
Compare two values for equality.
dimensionedScalar y1(const dimensionedScalar &ds)
static void applyOctant(point &p, unsigned octant)
static scalar findRootEllipsoidDistance(const scalar r0, const scalar r1, const scalar z0, const scalar z1, const scalar z2, scalar g)
dimensionedScalar sqrt(const dimensionedScalar &ds)
static constexpr int maxIters
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
static scalar findRootEllipseDistance(const scalar r0, const scalar z0, const scalar z1, scalar g)
static scalar distanceToEllipse(const scalar e0, const scalar e1, const scalar y0, const scalar y1, scalar &x0, scalar &x1)
static constexpr const zero Zero
Global zero (0)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
static scalar vectorMagSqr(const scalar x, const scalar y)
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)
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.
scalarField samples(nIntervals, Zero)