67 centres[0] = 0.5*(point1_ + point2_);
70 if (radius1_ > radius2_)
87 auto&
pts = tpts.ref();
96 void Foam::searchableCone::findNearestAndNormal
99 const scalar nearestDistSqr,
104 vector v(sample - point1_);
107 const scalar parallel = (v & unitDir_);
110 v -= parallel*unitDir_;
112 const scalar magV =
mag(v);
116 point disk1Point(point1_ +
min(
max(magV, innerRadius1_), radius1_)*v);
117 vector disk1Normal(-unitDir_);
120 point disk2Point(point2_ +
min(
max(magV, innerRadius2_), radius2_)*v);
121 vector disk2Normal(unitDir_);
125 point nearCone(point::uniform(GREAT));
126 vector normalCone(1, 0, 0);
129 point iCnearCone(point::uniform(GREAT));
130 vector iCnormalCone(1, 0, 0);
132 point projPt1 = point1_+ radius1_*v;
133 point projPt2 = point2_+ radius2_*v;
135 vector p1 = (projPt1 - point1_);
136 if (
mag(p1) > ROOTVSMALL)
144 vector a = (sample - projPt1);
146 if (
mag(a) <= ROOTVSMALL)
149 a = (sample - projPt2);
152 nearCone = (a &
b)*
b + projPt2;
160 nearCone = (a &
b)*
b + projPt1;
167 if (innerRadius1_ > 0 || innerRadius2_ > 0)
170 point iCprojPt1 = point1_+ innerRadius1_*v;
171 point iCprojPt2 = point2_+ innerRadius2_*v;
180 vector iCa(sample - iCprojPt1);
182 if (
mag(iCa) <= ROOTVSMALL)
184 iCa = (sample - iCprojPt2);
187 iCnearCone = (iCa & iCb)*iCb+iCprojPt2;
189 vector b1 = (iCp1 & iCb)*iCb;
195 iCnearCone = (iCa & iCb)*iCb+iCprojPt1;
198 vector b1 = (iCp1 & iCb)*iCb;
208 FixedList<scalar, 4> dist;
209 dist[0] = sample.distSqr(nearCone);
210 dist[1] = sample.distSqr(disk1Point);
211 dist[2] = sample.distSqr(disk2Point);
212 dist[3] = sample.distSqr(iCnearCone);
214 const label minI =
findMin(dist);
224 vector v1(nearCone-point1_);
225 scalar para = (v1 & unitDir_);
228 const scalar magV1 =
mag(v1);
231 if (para < 0.0 && magV1 >= radius1_)
235 nearCone = disk1Point;
237 else if (para < 0.0 && magV1 < radius1_)
240 nearCone = disk1Point;
241 normalCone = disk1Normal;
243 else if (para > magDir_ && magV1 >= radius2_)
247 nearCone = disk2Point;
249 else if (para > magDir_ && magV1 < radius2_)
252 nearCone = disk2Point;
253 normalCone = disk2Normal;
256 info.setPoint(nearCone);
257 nearNormal = normalCone;
261 info.setPoint(disk1Point);
262 nearNormal = disk1Normal;
266 info.setPoint(disk2Point);
267 nearNormal = disk2Normal;
272 vector v1(iCnearCone-point1_);
273 scalar para = (v1 & unitDir_);
277 const scalar magV1 =
mag(v1);
280 if (para < 0.0 && magV1 >= innerRadius1_)
282 iCnearCone = disk1Point;
284 else if (para < 0.0 && magV1 < innerRadius1_)
286 iCnearCone = disk1Point;
287 iCnormalCone = disk1Normal;
289 else if (para > magDir_ && magV1 >= innerRadius2_)
291 iCnearCone = disk2Point;
293 else if (para > magDir_ && magV1 < innerRadius2_)
295 iCnearCone = disk2Point;
296 iCnormalCone = disk2Normal;
299 info.setPoint(iCnearCone);
300 nearNormal = iCnormalCone;
304 if (info.point().distSqr(sample) < nearestDistSqr)
312 Foam::scalar Foam::searchableCone::radius2
314 const searchableCone& cone,
318 const vector x = (pt-cone.point1_) ^ cone.unitDir_;
325 void Foam::searchableCone::findLineAll
327 const searchableCone& cone,
328 const scalar innerRadius1,
329 const scalar innerRadius2,
339 vector point1Start(start-cone.point1_);
340 vector point2Start(start-cone.point2_);
344 scalar s1 = point1Start & (cone.unitDir_);
345 scalar s2 = point1End & (cone.unitDir_);
347 if ((s1 < 0.0 && s2 < 0.0) || (s1 > cone.magDir_ && s2 > cone.magDir_))
354 scalar magV =
mag(V);
355 if (magV < ROOTVSMALL)
371 scalar tNear = VGREAT;
372 scalar tFar = VGREAT;
374 scalar radius_sec = cone.radius1_;
378 scalar
s = (V & unitDir_);
382 tPoint2 = -(point2Start&(cone.unitDir_))/
s;
384 if (tPoint2 < tPoint1)
386 std::swap(tPoint1, tPoint2);
388 if (tPoint1 > magV || tPoint2 < 0)
393 if (tPoint1 >= 0 && tPoint1 <= magV)
395 scalar r2 = radius2(cone, start+tPoint1*V);
396 vector p1 = (start+tPoint1*V-point1_);
397 vector p2 = (start+tPoint1*V-point2_);
398 radius_sec = cone.radius1_;
399 scalar inC_radius_sec = innerRadius1_;
401 if (
mag(p2&(cone.unitDir_)) <
mag(p1&(cone.unitDir_)))
403 radius_sec = cone.radius2_;
404 inC_radius_sec = innerRadius2_;
407 if (r2 <=
sqr(radius_sec) && r2 >=
sqr(inC_radius_sec))
412 if (tPoint2 >= 0 && tPoint2 <= magV)
415 vector p1 = (start+tPoint2*V-cone.point1_);
416 vector p2 = (start+tPoint2*V-cone.point2_);
417 radius_sec = cone.radius1_;
418 scalar inC_radius_sec = innerRadius1_;
419 if (
mag(p2&(cone.unitDir_)) <
mag(p1&(cone.unitDir_)))
421 radius_sec = cone.radius2_;
422 inC_radius_sec = innerRadius2_;
424 scalar r2 = radius2(cone, start+tPoint2*V);
425 if (r2 <=
sqr(radius_sec) && r2 >=
sqr(inC_radius_sec))
452 scalar deltaRadius = cone.radius2_-cone.radius1_;
453 if (
mag(deltaRadius) <= ROOTVSMALL)
455 vector point1Start(start-cone.point1_);
456 const vector x = point1Start ^ cone.unitDir_;
457 const vector y = V ^ cone.unitDir_;
458 const scalar d =
sqr(0.5*(cone.radius1_ + cone.radius2_));
466 vector va = cone.unitDir_;
473 cone.unitDir_*cone.radius1_*
mag(cone.point2_-cone.point1_)
477 scalar l2 =
sqr(deltaRadius)+
sqr(cone.magDir_);
478 scalar sqrCosAlpha =
sqr(cone.magDir_)/l2;
479 scalar sqrSinAlpha =
sqr(deltaRadius)/l2;
483 vector p1 = (delP-(delP&va)*va);
485 a = sqrCosAlpha*((v1-
p*va)&(v1-
p*va))-sqrSinAlpha*
p*
p;
487 2.0*sqrCosAlpha*(a1&p1)
488 -2.0*sqrSinAlpha*(v1&va)*(delP&va);
491 *((delP-(delP&va)*va)&(delP-(delP&va)*va))
492 -sqrSinAlpha*
sqr(delP&va);
495 const scalar disc =
b*
b-4.0*a*
c;
505 else if (disc < ROOTVSMALL)
508 if (
mag(a) > ROOTVSMALL)
512 if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
541 if (
mag(a) > ROOTVSMALL)
543 scalar sqrtDisc =
sqrt(disc);
545 t1 = (-
b - sqrtDisc)/(2.0*a);
546 t2 = (-
b + sqrtDisc)/(2.0*a);
552 if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
565 if (t2>=0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
590 if (tNear>=0 && tNear <= magV)
592 near.hitPoint(start+tNear*V);
596 far.hitPoint(start+tFar*V);
600 else if (tFar>=0 && tFar <= magV)
602 near.hitPoint(start+tFar*V);
608 void Foam::searchableCone::insertHit
612 List<pointIndexHit>& info,
616 scalar smallDistSqr = SMALL*
magSqr(
end-start);
618 scalar hitMagSqr = hit.hitPoint().distSqr(start);
622 scalar d2 = info[i].hitPoint().distSqr(start);
624 if (d2 > hitMagSqr+smallDistSqr)
627 label sz = info.size();
629 for (label j = sz; j > i; --j)
636 else if (d2 > hitMagSqr-smallDistSqr)
643 label sz = info.size();
684 if (radius2_ >= radius1_)
705 Foam::searchableCone::searchableCone
709 const scalar radius1,
710 const scalar innerRadius1,
712 const scalar radius2,
713 const scalar innerRadius2
719 innerRadius1_(innerRadius1),
722 innerRadius2_(innerRadius2),
723 magDir_(
mag(point2_-point1_)),
724 unitDir_((point2_-point1_)/magDir_)
726 bounds() = calcBounds();
730 Foam::searchableCone::searchableCone
738 radius1_(
dict.
get<scalar>(
"radius1")),
739 innerRadius1_(
dict.getOrDefault<scalar>(
"innerRadius1", 0)),
741 radius2_(
dict.
get<scalar>(
"radius2")),
742 innerRadius2_(
dict.getOrDefault<scalar>(
"innerRadius2", 0)),
743 magDir_(
mag(point2_-point1_)),
744 unitDir_((point2_-point1_)/magDir_)
746 bounds() = calcBounds();
754 if (regions_.empty())
757 regions_[0] =
"region0";
767 List<pointIndexHit>& info
775 findNearestAndNormal(
samples[i], nearestDistSqr[i], info[i], normal);
784 List<pointIndexHit>& info
787 info.setSize(start.size());
803 if (!info[i].hit() &&
b.hit())
810 if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
814 searchableCone innerCone
830 newEnd = info[i].point();
866 List<pointIndexHit>& info
869 info.setSize(start.size());
885 if (!info[i].hit() &&
b.hit())
890 if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
920 if (!info[i].hit() &&
b.hit())
930 void Foam::searchableCone::findLineAll
934 List<List<pointIndexHit>>& info
937 info.setSize(start.size());
981 if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
1013 insertHit(start[i],
end[i], info[i], near);
1017 insertHit(start[i],
end[i], info[i], far);
1026 const List<pointIndexHit>& info,
1030 region.setSize(info.size());
1049 findNearestAndNormal
1064 List<volumeType>& volType
1078 const scalar parallel = (v & unitDir_);
1081 if (parallel < 0 || parallel > magDir_)
1086 const scalar radius_sec =
1087 radius1_ + parallel * (radius2_-radius1_)/magDir_;
1089 const scalar radius_sec_inner =
1090 innerRadius1_ + parallel * (innerRadius2_-innerRadius1_)/magDir_;
1093 v -= parallel*unitDir_;
1095 if (
mag(v) >= radius_sec_inner &&
mag(v) <= radius_sec)
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
void size(const label n)
Older name for setAddressableSize.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
A list of keyword definitions, which are a keyword followed by a number of values (eg...
virtual tmp< pointField > points() const
Get the points that define the surface.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
virtual void rename(const word &newName)
Rename the object.
scalarField samples(nIntervals, Zero)
dimensionedScalar sqrt(const dimensionedScalar &ds)
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
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.
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find nearest intersection on line from start to end.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
virtual const wordList & regions() const
Names of regions.
virtual void getRegion(const List< pointIndexHit > &, labelList ®ion) const
From a set of points and indices get the region.
Macros for easy insertion into run-time selection tables.
addNamedToRunTimeSelectionTable(topoSetCellSource, badQualityToCell, word, badQuality)
virtual void getVolumeType(const pointField &points, List< volumeType > &volType) const
Determine type (inside/outside/mixed) for point.
#define forAll(list, i)
Loop across all elements in list.
label findMin(const ListType &input, label start=0)
Linear search for the index of the min element, similar to std::min_element but for lists and returns...
vectorField pointField
pointField is a vectorField.
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
void setSize(const label n)
Alias for resize()
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual void findNearest(const pointField &sample, const scalarField &nearestDistSqr, List< pointIndexHit > &) const
Find nearest point on cylinder.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
A location inside the volume.
virtual void boundingSpheres(pointField ¢res, scalarField &radiusSqr) const
Get bounding spheres (centre and radius squared), one per element.
A location outside the volume.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
virtual void getNormal(const List< pointIndexHit > &, vectorField &normal) const
From a set of points and indices get the normal.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
defineTypeNameAndDebug(combustionModel, 0)
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
vector point
Point is a vector.
const dimensionedScalar c
Speed of light in a vacuum.
virtual void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find any intersection on line from start to end.
virtual tmp< pointField > coordinates() const
Get representative set of element coordinates.
List< label > labelList
A List of labels.
A class for managing temporary objects.
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))
Defines the attributes of an object for which implicit objectRegistry management is supported...
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
static constexpr const zero Zero
Global zero (0)