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(); }
156 static constexpr
int maxIters = 100;
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)
517 for (
direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
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)
573 const vector n(sample - origin_);
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))
597 vector relPt(sample - origin_);
600 const unsigned octant =
getOctant(relPt);
609 point& near = info.point();
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_))
807 Foam::searchableSphere::searchableSphere
816 dict.getCompat<
vector>(
"origin", {{
"centre", -1806}}),
832 origin_.x() + radii_.x() *
cos(theta)*
sin(
phi),
833 origin_.y() + radii_.y() *
sin(theta)*
sin(
phi),
858 if (order_.shape == shapeType::SPHERE)
878 for (
direction dir = 0; dir < vector::nComponents; ++dir)
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]);
956 List<pointIndexHit>& info
959 info.resize(start.size());
968 findLineAll(start[i],
end[i], info[i],
b);
970 if (!info[i].hit() &&
b.hit())
982 List<pointIndexHit>& info
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
1008 List<List<pointIndexHit>>& info
1011 info.resize(start.size());
1017 findLineAll(start[i],
end[i], near, far);
1051 const List<pointIndexHit>& info,
1055 region.resize(info.size());
1072 if (order_.shape == shapeType::SPHERE)
1082 normal[i] = scalePoint(info[i].
point());
1084 normal[i].x() /= radii_.x();
1085 normal[i].y() /= radii_.y();
1086 normal[i].z() /= radii_.z();
1102 List<volumeType>& volType
1107 if (order_.shape == shapeType::SPHERE)
1111 const scalar rad2 =
sqr(radius());
static scalar vectorMagSqr(const scalar x, const scalar y)
ITstream & lookup(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return an entry data stream. FatalIOError if not found, or not a stream. ...
void size(const label n)
Older name for setAddressableSize.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void resize(const label len)
Adjust allocated size of list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Return any intersection on segment from start to end.
bool equal(const T &a, const T &b)
Compare two values for equality.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
static scalar findRootEllipseDistance(const scalar r0, const scalar z0, const scalar z1, scalar g)
vector surfaceNormal(const scalar theta, const scalar phi) const
Surface normal on the sphere at given location.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
A token holds an item read from Istream.
scalarField samples(nIntervals, Zero)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
static scalar findRootEllipsoidDistance(const scalar r0, const scalar r1, const scalar z0, const scalar z1, const scalar z2, scalar g)
dimensionedScalar y0(const dimensionedScalar &ds)
A bounding box defined in terms of min/max extrema points.
This class describes the interaction of an object (often a face) and a point. It carries the info of ...
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
const point & min() const noexcept
Minimum describing the bounding box.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T. FatalIOError if not found, or if the number of tokens is incorrect.
virtual const boundBox & bounds() const
Return const reference to boundBox.
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.
Macros for easy insertion into run-time selection tables.
addNamedToRunTimeSelectionTable(topoSetCellSource, badQualityToCell, word, badQuality)
#define forAll(list, i)
Loop across all elements in list.
const point & centre() const noexcept
The centre (origin) of the sphere.
const point & max() const noexcept
Maximum describing the bounding box.
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
virtual bool overlaps(const boundBox &bb) const
Does any part of the surface overlap the supplied bound box?
dimensionedScalar cos(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
static constexpr scalar tolCloseness
bool valid() const
Bounding box is non-inverted - same as good().
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual void getRegion(const List< pointIndexHit > &, labelList ®ion) const
From a set of points and indices get the region.
static vector getRadius(const word &name, const dictionary &dict)
virtual void getVolumeType(const pointField &points, List< volumeType > &volType) const
Determine type (inside/outside/mixed) for point.
point surfacePoint(const scalar theta, const scalar phi) const
A point on the sphere at given location.
A location inside the volume.
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
A location outside the volume.
dimensionedScalar y1(const dimensionedScalar &ds)
const uniformDimensionedVectorField & g
dimensionedScalar sin(const dimensionedScalar &ds)
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
static scalar vectorMag(const scalar x, const scalar y)
defineTypeNameAndDebug(combustionModel, 0)
Searching on general spheroid.
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
vector point
Point is a vector.
static scalar distanceToEllipse(const scalar e0, const scalar e1, const scalar y0, const scalar y1, scalar &x0, scalar &x1)
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)
virtual const wordList & regions() const
Names of regions.
virtual void boundingSpheres(pointField ¢res, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
static constexpr int maxIters
List< label > labelList
A List of labels.
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, IOobject::NO_REGISTER)
void normalise()
Inplace normalise this field. Generally a no-op except for vector fields.
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))
bool overlaps(const boundBox &bb) const
Overlaps/touches boundingBox?
dimensionedScalar hypot(const dimensionedScalar &x, const dimensionedScalar &y)
static unsigned getOctant(const point &p)
Defines the attributes of an object for which implicit objectRegistry management is supported...
static void applyOctant(point &p, unsigned octant)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
#define InfoInFunction
Report an information message using Foam::Info.
static constexpr const zero Zero
Global zero (0)