diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index bd0abc24a9be4c307f0476f391e3ef3183699781..23b11884e2b53c9f95d3d942587bc0e350e256ef 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -65,6 +65,7 @@ searchableSurface = searchableSurface
$(searchableSurface)/searchableBox.C
$(searchableSurface)/searchableRotatedBox.C
$(searchableSurface)/searchableCylinder.C
+$(searchableSurface)/searchableCone.C
$(searchableSurface)/searchableDisk.C
$(searchableSurface)/searchablePlane.C
$(searchableSurface)/searchablePlate.C
diff --git a/src/meshTools/searchableSurface/searchableCone.C b/src/meshTools/searchableSurface/searchableCone.C
new file mode 100644
index 0000000000000000000000000000000000000000..c1047f8ec76f677e3ad088b4137c0634c4bf838b
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableCone.C
@@ -0,0 +1,1126 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "searchableCone.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchableCone, 0);
+addToRunTimeSelectionTable(searchableSurface, searchableCone, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::tmp<Foam::pointField> Foam::searchableCone::coordinates() const
+{
+ tmp<pointField> tCtrs(new pointField(1, 0.5*(point1_ + point2_)));
+
+ return tCtrs;
+}
+
+
+void Foam::searchableCone::boundingSpheres
+(
+ pointField& centres,
+ scalarField& radiusSqr
+) const
+{
+ centres.setSize(1);
+ centres[0] = 0.5*(point1_ + point2_);
+
+ radiusSqr.setSize(1);
+ if (radius1_ > radius2_)
+ {
+ radiusSqr[0] = Foam::magSqr(point1_-centres[0]) + Foam::sqr(radius1_);
+ }
+ else
+ {
+ radiusSqr[0] = Foam::magSqr(point2_-centres[0]) + Foam::sqr(radius2_);
+ }
+
+ // Add a bit to make sure all points are tested inside
+ radiusSqr += Foam::sqr(SMALL);
+}
+
+
+Foam::tmp<Foam::pointField> Foam::searchableCone::points() const
+{
+ tmp<pointField> tPts(new pointField(2));
+ pointField& pts = tPts();
+
+ pts[0] = point1_;
+ pts[1] = point2_;
+
+ return tPts;
+}
+
+
+void Foam::searchableCone::findNearest
+(
+ const point& sample,
+ const scalar nearestDistSqr,
+ pointIndexHit& info,
+ vector& nearNormal
+) const
+{
+ vector v(sample - point1_);
+
+ // Decompose sample-point1 into normal and parallel component
+ scalar parallel = (v & unitDir_);
+
+ // Remove the parallel component and normalise
+ v -= parallel*unitDir_;
+
+ scalar magV = mag(v);
+ if (magV < ROOTVSMALL)
+ {
+ v = vector::zero;
+ }
+ else
+ {
+ v /= magV;
+ }
+
+ // Nearest and normal on disk at point1
+ point disk1Point(point1_ + min(max(magV, innerRadius1_), radius1_)*v);
+ vector disk1Normal(-unitDir_);
+
+ // Nearest and normal on disk at point2
+ point disk2Point(point2_ + min(max(magV, innerRadius2_), radius2_)*v);
+ vector disk2Normal(unitDir_);
+
+ // Nearest and normal on cone. Initialise to far-away point so if not
+ // set it picks one of the disk points
+ point nearCone(GREAT, GREAT, GREAT);
+ vector normalCone(1, 0, 0);
+
+ // Nearest and normal on inner cone. Initialise to far-away point
+ point iCnearCone(GREAT, GREAT, GREAT);
+ vector iCnormalCone(1, 0, 0);
+
+ point projPt1 = point1_+ radius1_*v;
+ point projPt2 = point2_+ radius2_*v;
+
+ vector p1 = (projPt1 - point1_);
+ if (mag(p1) > ROOTVSMALL)
+ {
+ p1 /= mag(p1);
+
+ // Find vector along the two end of cone
+ vector b(projPt2 - projPt1);
+ scalar magS = mag(b);
+ b /= magS;
+
+ // Find the vector along sample pt and pt at one end of conde
+ vector a(sample - projPt1);
+
+ if (mag(a) <= ROOTVSMALL)
+ {
+ // Exception: sample on disk1. Redo with projPt2.
+ vector a(sample - projPt2);
+ // Find normal unitvector
+ nearCone = (a & b)*b+projPt2;
+ vector b1 = (p1 & b)*b;
+ normalCone = p1 - b1;
+ normalCone /= mag(normalCone);
+ }
+ else
+ {
+ // Find neartest point on cone surface
+ nearCone = (a & b)*b+projPt1;
+ // Find projection along surface of cone
+ vector b1 = (p1 & b)*b;
+ normalCone = p1 - b1;
+ normalCone /= mag(normalCone);
+ }
+
+ if (innerRadius1_ > 0 && innerRadius2_ > 0)
+ {
+ // Same for inner radius
+ point iCprojPt1 = point1_+ innerRadius1_*v;
+ point iCprojPt2 = point2_+ innerRadius2_*v;
+
+ vector iCp1 = (iCprojPt1 - point1_);
+ iCp1 /= mag(iCp1);
+
+ // Find vector along the two end of cone
+ vector iCb(iCprojPt2 - iCprojPt1);
+ magS = mag(iCb);
+ iCb /= magS;
+
+
+ // Find the vector along sample pt and pt at one end of conde
+ vector iCa(sample - iCprojPt1);
+
+ if (mag(iCa) <= ROOTVSMALL)
+ {
+ vector iCa(sample - iCprojPt2);
+
+ // Find normal unitvector
+ iCnearCone = (iCa & iCb)*iCb+iCprojPt2;
+ vector b1 = (iCp1 & iCb)*iCb;
+ iCnormalCone = iCp1 - b1;
+ iCnormalCone /= mag(iCnormalCone);
+ }
+ else
+ {
+ // Find nearest point on cone surface
+ iCnearCone = (iCa & iCb)*iCb+iCprojPt1;
+ // Find projection along surface of cone
+ vector b1 = (iCp1 & iCb)*iCb;
+ iCnormalCone = iCp1 - b1;
+ iCnormalCone /= mag(iCnormalCone);
+ }
+ }
+ }
+
+
+ // Select nearest out of the 4 points (outer cone, disk1, disk2, inner
+ // cone)
+
+ FixedList<scalar, 4> dist;
+ dist[0] = magSqr(nearCone-sample);
+ dist[1] = magSqr(disk1Point-sample);
+ dist[2] = magSqr(disk2Point-sample);
+ dist[3] = magSqr(iCnearCone-sample);
+
+ label minI = findMin(dist);
+
+
+ // Snap the point to the corresponding surface
+
+ if (minI == 0) // Outer cone
+ {
+ // Closest to (infinite) outer cone. See if needs clipping to end disks
+
+ {
+ vector v1(nearCone-point1_);
+ scalar para = (v1 & unitDir_);
+ // Remove the parallel component and normalise
+ v1 -= para*unitDir_;
+ scalar magV1 = mag(v1);
+ v1 = v1/magV1;
+ if (para < 0.0 && magV1 >= radius1_)
+ {
+ // Near point 1. Set point to intersection of disk and cone.
+ // Keep normal from cone.
+ nearCone = disk1Point;
+ }
+ else if (para < 0.0 && magV1 < radius1_)
+ {
+ // On disk1
+ nearCone = disk1Point;
+ normalCone = disk1Normal;
+ }
+ else if (para > magDir_ && magV1 >= radius2_)
+ {
+ // Near point 2. Set point to intersection of disk and cone.
+ // Keep normal from cone.
+ nearCone = disk2Point;
+ }
+ else if (para > magDir_ && magV1 < radius2_)
+ {
+ // On disk2
+ nearCone = disk2Point;
+ normalCone = disk2Normal;
+ }
+ }
+ info.setPoint(nearCone);
+ nearNormal = normalCone;
+ }
+ else if (minI == 1) // Near to disk1
+ {
+ info.setPoint(disk1Point);
+ nearNormal = disk1Normal;
+ }
+ else if (minI == 2) // Near to disk2
+ {
+ info.setPoint(disk2Point);
+ nearNormal = disk2Normal;
+ }
+ else // Near to inner cone
+ {
+ {
+ vector v1(iCnearCone-point1_);
+ scalar para = (v1 & unitDir_);
+ // Remove the parallel component and normalise
+ v1 -= para*unitDir_;
+ scalar magV1 = mag(v1);
+ v1 = v1/magV1;
+
+ if (para < 0.0 && magV1 >= innerRadius1_)
+ {
+ iCnearCone = disk1Point;
+ }
+ else if (para < 0.0 && magV1 < innerRadius1_)
+ {
+ iCnearCone = disk1Point;
+ iCnormalCone = disk1Normal;
+ }
+ else if (para > magDir_ && magV1 >= innerRadius2_)
+ {
+ iCnearCone = disk2Point;
+ }
+ else if (para > magDir_ && magV1 < innerRadius2_)
+ {
+ iCnearCone = disk2Point;
+ iCnormalCone = disk2Normal;
+ }
+ }
+ info.setPoint(iCnearCone);
+ nearNormal = iCnormalCone;
+ }
+
+
+ if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
+ {
+ info.setHit();
+ info.setIndex(0);
+ }
+}
+
+
+Foam::scalar Foam::searchableCone::radius2
+(
+ const searchableCone& cone,
+ const point& pt
+)
+{
+ const vector x = (pt-cone.point1_) ^ cone.unitDir_;
+ return x&x;
+}
+
+
+// From http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 -
+// intersection of cylinder with ray
+void Foam::searchableCone::findLineAll
+(
+ const searchableCone& cone,
+ const scalar innerRadius1,
+ const scalar innerRadius2,
+ const point& start,
+ const point& end,
+ pointIndexHit& near,
+ pointIndexHit& far
+) const
+{
+ near.setMiss();
+ far.setMiss();
+
+ vector point1Start(start-cone.point1_);
+ vector point2Start(start-cone.point2_);
+ vector point1End(end-cone.point1_);
+
+ // Quick rejection of complete vector outside endcaps
+ scalar s1 = point1Start&(cone.unitDir_);
+ scalar s2 = point1End&(cone.unitDir_);
+
+ if ((s1 < 0.0 && s2 < 0.0) || (s1 > cone.magDir_ && s2 > cone.magDir_))
+ {
+ return;
+ }
+
+ // Line as P = start+t*V where V is unit vector and t=[0..mag(end-start)]
+ vector V(end-start);
+ scalar magV = mag(V);
+ if (magV < ROOTVSMALL)
+ {
+ return;
+ }
+ V /= magV;
+
+
+ // We now get the nearest intersections to start. This can either be
+ // the intersection with the end plane or with the cylinder side.
+
+ // Get the two points (expressed in t) on the end planes. This is to
+ // clip any cylinder intersection against.
+ scalar tPoint1;
+ scalar tPoint2;
+
+ // Maintain the two intersections with the endcaps
+ scalar tNear = VGREAT;
+ scalar tFar = VGREAT;
+
+ scalar radius_sec = cone.radius1_;
+
+ {
+ // Find dot product: mag(s)>VSMALL suggest that it is greater
+ scalar s = (V&unitDir_);
+ if (mag(s) > VSMALL)
+ {
+ tPoint1 = -s1/s;
+ tPoint2 = -(point2Start&(cone.unitDir_))/s;
+
+ if (tPoint2 < tPoint1)
+ {
+ Swap(tPoint1, tPoint2);
+ }
+ if (tPoint1 > magV || tPoint2 < 0)
+ {
+ return;
+ }
+ // See if the points on the endcaps are actually inside the cylinder
+ if (tPoint1 >= 0 && tPoint1 <= magV)
+ {
+ scalar r2 = radius2(cone, start+tPoint1*V);
+ vector p1 = (start+tPoint1*V-point1_);
+ vector p2 = (start+tPoint1*V-point2_);
+ radius_sec = cone.radius1_;
+ scalar inC_radius_sec = innerRadius1_;
+
+ if (mag(p2&(cone.unitDir_)) < mag(p1&(cone.unitDir_)))
+ {
+ radius_sec = cone.radius2_;
+ inC_radius_sec = innerRadius2_;
+ }
+
+ if (r2 <= sqr(radius_sec) && r2 >= sqr(inC_radius_sec))
+ {
+ tNear = tPoint1;
+ }
+ }
+ if (tPoint2 >= 0 && tPoint2 <= magV)
+ {
+ // Decompose sample-point1 into normal and parallel component
+ vector p1 = (start+tPoint2*V-cone.point1_);
+ vector p2 = (start+tPoint2*V-cone.point2_);
+ radius_sec = cone.radius1_;
+ scalar inC_radius_sec = innerRadius1_;
+ if (mag(p2&(cone.unitDir_)) < mag(p1&(cone.unitDir_)))
+ {
+ radius_sec = cone.radius2_;
+ inC_radius_sec = innerRadius2_;
+ }
+ scalar r2 = radius2(cone, start+tPoint2*V);
+ if (r2 <= sqr(radius_sec) && r2 >= sqr(inC_radius_sec))
+ {
+ // Check if already have a near hit from point1
+ if (tNear <= 1)
+ {
+ tFar = tPoint2;
+ }
+ else
+ {
+ tNear = tPoint2;
+ }
+ }
+ }
+ }
+ else
+ {
+ // Vector perpendicular to cylinder. Check for outside already done
+ // above so just set tpoint to allow all.
+ tPoint1 = -VGREAT;
+ tPoint2 = VGREAT;
+ }
+ }
+
+
+ vector va = cone.unitDir_;
+ point pa =
+ cone.unitDir_*cone.radius1_*mag(cone.point2_-cone.point1_)
+ /(cone.radius1_-cone.radius2_)
+ +cone.point1_;
+
+ scalar deltaRadius = cone.radius2_-cone.radius1_;
+ scalar l2 = sqr(deltaRadius)+sqr(cone.magDir_);
+ scalar sqrCosAlpha = sqr(cone.magDir_)/l2;
+ scalar sqrSinAlpha = sqr(deltaRadius)/l2;
+
+
+ vector delP(start-pa);
+ vector v1(end-start);
+ v1 = v1/mag(v1);
+
+ scalar p = (va&v1);
+ vector a1 = (v1-p*va);
+ vector p1 = (delP-(delP&va)*va);
+ const scalar a = sqrCosAlpha*((v1-p*va)&(v1-p*va))-sqrSinAlpha*p*p;
+ const scalar b =
+ 2.0*sqrCosAlpha*(a1&p1)
+ -2.0*sqrSinAlpha*(v1&va)*(delP&va);
+ const scalar c =
+ sqrCosAlpha
+ *((delP-(delP&va)*va)&(delP-(delP&va)*va))
+ -sqrSinAlpha*sqr(delP&va);
+ const scalar disc = b*b-4.0*a*c;
+
+ scalar t1 = -VGREAT;
+ scalar t2 = VGREAT;
+
+ if (disc < 0)
+ {
+ // Fully outside
+ return;
+ }
+ else if (disc < ROOTVSMALL)
+ {
+ // Single solution
+ if (mag(a) > ROOTVSMALL)
+ {
+ t1 = -b/(2.0*a);
+
+ if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
+ {
+ // Valid. Insert sorted.
+ if (t1 < tNear)
+ {
+ tFar = tNear;
+ tNear = t1;
+ }
+ else if (t1 < tFar)
+ {
+ tFar = t1;
+ }
+ }
+ else
+ {
+ return;
+ }
+ }
+ else
+ {
+ // Aligned with axis. Check if outside radius
+ if (c > 0.0)
+ {
+ return;
+ }
+ }
+ }
+ else
+ {
+ if (mag(a) > ROOTVSMALL)
+ {
+ scalar sqrtDisc = sqrt(disc);
+
+ t1 = (-b - sqrtDisc)/(2.0*a);
+ t2 = (-b + sqrtDisc)/(2.0*a);
+ if (t2 < t1)
+ {
+ Swap(t1, t2);
+ }
+
+ if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
+ {
+ // Valid. Insert sorted.
+ if (t1 < tNear)
+ {
+ tFar = tNear;
+ tNear = t1;
+ }
+ else if (t1 < tFar)
+ {
+ tFar = t1;
+ }
+ }
+ if (t2>=0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
+ {
+ // Valid. Insert sorted.
+ if (t2 < tNear)
+ {
+ tFar = tNear;
+ tNear = t2;
+ }
+ else if (t2 < tFar)
+ {
+ tFar = t2;
+ }
+ }
+ }
+ else
+ {
+ // Aligned with axis. Check if outside radius
+ if (c > 0.0)
+ {
+ return;
+ }
+ }
+ }
+
+ // Check tNear, tFar
+ if (tNear>=0 && tNear <= magV)
+ {
+ near.setPoint(start+tNear*V);
+ near.setHit();
+ near.setIndex(0);
+ if (tFar <= magV)
+ {
+ far.setPoint(start+tFar*V);
+ far.setHit();
+ far.setIndex(0);
+ }
+ }
+ else if (tFar>=0 && tFar <= magV)
+ {
+ near.setPoint(start+tFar*V);
+ near.setHit();
+ near.setIndex(0);
+ }
+}
+
+
+void Foam::searchableCone::insertHit
+(
+ const point& start,
+ const point& end,
+ List<pointIndexHit>& info,
+ const pointIndexHit& hit
+) const
+{
+ scalar smallDistSqr = SMALL*magSqr(end-start);
+
+ scalar hitMagSqr = magSqr(hit.hitPoint()-start);
+
+ forAll(info, i)
+ {
+ scalar d2 = magSqr(info[i].hitPoint()-start);
+
+ if (d2 > hitMagSqr+smallDistSqr)
+ {
+ // Insert at i.
+ label sz = info.size();
+ info.setSize(sz+1);
+ for (label j = sz; j > i; --j)
+ {
+ info[j] = info[j-1];
+ }
+ info[i] = hit;
+ return;
+ }
+ else if (d2 > hitMagSqr-smallDistSqr)
+ {
+ // hit is same point as info[i].
+ return;
+ }
+ }
+ // Append
+ label sz = info.size();
+ info.setSize(sz+1);
+ info[sz] = hit;
+}
+
+
+Foam::boundBox Foam::searchableCone::calcBounds() const
+{
+ // Adapted from
+ // http://www.gamedev.net/community/forums
+ // /topic.asp?topic_id=338522&forum_id=20&gforum_id=0
+
+ // Let cylinder have end points A,B and radius r,
+
+ // Bounds in direction X (same for Y and Z) can be found as:
+ // Let A.X<B.X (otherwise swap points)
+ // Good approximate lowest bound is A.X-r and highest is B.X+r (precise for
+ // capsule). At worst, in one direction it can be larger than needed by 2*r.
+
+ // Accurate bounds for cylinder is
+ // A.X-kx*r, B.X+kx*r
+ // where
+ // kx=sqrt(((A.Y-B.Y)^2+(A.Z-B.Z)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
+
+ // similar thing for Y and Z
+ // (i.e.
+ // ky=sqrt(((A.X-B.X)^2+(A.Z-B.Z)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
+ // kz=sqrt(((A.X-B.X)^2+(A.Y-B.Y)^2)/((A.X-B.X)^2+(A.Y-B.Y)^2+(A.Z-B.Z)^2))
+ // )
+
+ // How derived: geometric reasoning. Bounds of cylinder is same as for 2
+ // circles centered on A and B. This sqrt thingy gives sine of angle between
+ // axis and direction, used to find projection of radius.
+
+ vector kr
+ (
+ sqrt(sqr(unitDir_.y()) + sqr(unitDir_.z())),
+ sqrt(sqr(unitDir_.x()) + sqr(unitDir_.z())),
+ sqrt(sqr(unitDir_.x()) + sqr(unitDir_.y()))
+ );
+
+ if (radius2_ >= radius1_)
+ {
+ kr *= radius2_;
+ }
+ else
+ {
+ kr *= radius1_;
+ }
+
+ point min = point1_ - kr;
+ point max = point1_ + kr;
+
+ min = ::Foam::min(min, point2_ - kr);
+ max = ::Foam::max(max, point2_ + kr);
+
+ return boundBox(min, max);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::searchableCone::searchableCone
+(
+ const IOobject& io,
+ const point& point1,
+ const scalar radius1,
+ const scalar innerRadius1,
+ const point& point2,
+ const scalar radius2,
+ const scalar innerRadius2
+)
+:
+ searchableSurface(io),
+ point1_(point1),
+ radius1_(radius1),
+ innerRadius1_(innerRadius1),
+ point2_(point2),
+ radius2_(radius2),
+ innerRadius2_(innerRadius2),
+ magDir_(mag(point2_-point1_)),
+ unitDir_((point2_-point1_)/magDir_)
+{
+ bounds() = calcBounds();
+}
+
+
+Foam::searchableCone::searchableCone
+(
+ const IOobject& io,
+ const dictionary& dict
+)
+:
+ searchableSurface(io),
+ point1_(dict.lookup("point1")),
+ radius1_(readScalar(dict.lookup("radius1"))),
+ innerRadius1_(dict.lookupOrDefault("innerRadius1", 0.0)),
+ point2_(dict.lookup("point2")),
+ radius2_(readScalar(dict.lookup("radius2"))),
+ innerRadius2_(dict.lookupOrDefault("innerRadius2", 0.0)),
+ magDir_(mag(point2_-point1_)),
+ unitDir_((point2_-point1_)/magDir_)
+{
+ bounds() = calcBounds();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::searchableCone::~searchableCone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchableCone::regions() const
+{
+ if (regions_.empty())
+ {
+ regions_.setSize(1);
+ regions_[0] = "region0";
+ }
+ return regions_;
+}
+
+
+void Foam::searchableCone::findNearest
+(
+ const pointField& samples,
+ const scalarField& nearestDistSqr,
+ List<pointIndexHit>& info
+) const
+{
+ info.setSize(samples.size());
+ forAll(samples, i)
+ {
+ vector normal;
+ findNearest(samples[i], nearestDistSqr[i], info[i], normal);
+ }
+}
+
+
+void Foam::searchableCone::findLine
+(
+ const pointField& start,
+ const pointField& end,
+ List<pointIndexHit>& info
+) const
+{
+ info.setSize(start.size());
+
+ forAll(start, i)
+ {
+ // Pick nearest intersection. If none intersected take second one.
+ pointIndexHit b;
+ findLineAll
+ (
+ *this,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ info[i],
+ b
+ );
+ if (!info[i].hit() && b.hit())
+ {
+ info[i] = b;
+ }
+ }
+
+
+ if (innerRadius1_ > 0.0 && innerRadius2_ > 0.0)
+ {
+ IOobject io(*this);
+ io.rename(name()+"Inner");
+ searchableCone innerCone
+ (
+ io,
+ point1_,
+ innerRadius1_,
+ 0.0,
+ point2_,
+ innerRadius2_,
+ 0.0
+ );
+
+ forAll(info, i)
+ {
+ point newEnd;
+ if (info[i].hit())
+ {
+ newEnd = info[i].hitPoint();
+ }
+ else
+ {
+ newEnd = end[i];
+ }
+ pointIndexHit near;
+ pointIndexHit far;
+ findLineAll
+ (
+ innerCone,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ newEnd,
+ near,
+ far
+ );
+
+ if (near.hit())
+ {
+ info[i] = near;
+ }
+ else if (far.hit())
+ {
+ info[i] = far;
+ }
+ }
+ }
+}
+
+
+void Foam::searchableCone::findLineAny
+(
+ const pointField& start,
+ const pointField& end,
+ List<pointIndexHit>& info
+) const
+{
+ info.setSize(start.size());
+
+ forAll(start, i)
+ {
+ // Discard far intersection
+ pointIndexHit b;
+ findLineAll
+ (
+ *this,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ info[i],
+ b
+ );
+ if (!info[i].hit() && b.hit())
+ {
+ info[i] = b;
+ }
+ }
+ if (innerRadius1_ > 0.0 && innerRadius2_ > 0.0)
+ {
+ IOobject io(*this);
+ io.rename(name()+"Inner");
+ searchableCone cone
+ (
+ io,
+ point1_,
+ innerRadius1_,
+ 0.0,
+ point2_,
+ innerRadius2_,
+ 0.0
+ );
+
+ forAll(info, i)
+ {
+ if (!info[i].hit())
+ {
+ pointIndexHit b;
+ findLineAll
+ (
+ cone,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ info[i],
+ b
+ );
+ if (!info[i].hit() && b.hit())
+ {
+ info[i] = b;
+ }
+ }
+ }
+ }
+}
+
+
+void Foam::searchableCone::findLineAll
+(
+ const pointField& start,
+ const pointField& end,
+ List<List<pointIndexHit> >& info
+) const
+{
+ info.setSize(start.size());
+
+ forAll(start, i)
+ {
+ pointIndexHit near, far;
+ findLineAll
+ (
+ *this,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ near,
+ far
+ );
+
+ if (near.hit())
+ {
+ if (far.hit())
+ {
+ info[i].setSize(2);
+ info[i][0] = near;
+ info[i][1] = far;
+ }
+ else
+ {
+ info[i].setSize(1);
+ info[i][0] = near;
+ }
+ }
+ else
+ {
+ if (far.hit())
+ {
+ info[i].setSize(1);
+ info[i][0] = far;
+ }
+ else
+ {
+ info[i].clear();
+ }
+ }
+ }
+
+ if (innerRadius1_ > 0.0 && innerRadius2_ > 0.0)
+ {
+ IOobject io(*this);
+ io.rename(name()+"Inner");
+ searchableCone cone
+ (
+ io,
+ point1_,
+ innerRadius1_,
+ 0.0,
+ point2_,
+ innerRadius2_,
+ 0.0
+ );
+
+ forAll(info, i)
+ {
+ pointIndexHit near;
+ pointIndexHit far;
+ findLineAll
+ (
+ cone,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ near,
+ far
+ );
+
+ if (near.hit())
+ {
+ insertHit(start[i], end[i], info[i], near);
+ }
+ if (far.hit())
+ {
+ insertHit(start[i], end[i], info[i], far);
+ }
+ }
+ }
+}
+
+
+void Foam::searchableCone::getRegion
+(
+ const List<pointIndexHit>& info,
+ labelList& region
+) const
+{
+ region.setSize(info.size());
+ region = 0;
+}
+
+
+void Foam::searchableCone::getNormal
+(
+ const List<pointIndexHit>& info,
+ vectorField& normal
+) const
+{
+ normal.setSize(info.size());
+ normal = vector::zero;
+
+ forAll(info, i)
+ {
+ if (info[i].hit())
+ {
+ pointIndexHit nearInfo;
+ findNearest
+ (
+ info[i].hitPoint(),
+ Foam::sqr(GREAT),
+ nearInfo,
+ normal[i]
+ );
+ }
+ }
+}
+
+
+void Foam::searchableCone::getVolumeType
+(
+ const pointField& points,
+ List<volumeType>& volType
+) const
+{
+ volType.setSize(points.size());
+ volType = volumeType::INSIDE;
+
+ scalar magDir = magDir_;
+
+ bool swap = false;
+ bool swapInner = false;
+
+ if (radius1_ > radius2_)
+ {
+ swap = true;
+ }
+
+ if (innerRadius1_ > innerRadius2_)
+ {
+ swapInner = true;
+ }
+
+ forAll(points, pointI)
+ {
+ const point& pt = points[pointI];
+
+ vector v(pt - point1_);
+
+ // Decompose sample-point1 into normal and parallel component
+ scalar parallel = v & unitDir_;
+ scalar comp = parallel;
+ scalar compInner = parallel;
+ if (swap)
+ {
+ comp = -(magDir-parallel);
+ }
+
+ if (swapInner)
+ {
+ compInner = -(magDir-parallel);
+ }
+
+ scalar radius_sec = radius1_+comp*(radius2_-radius1_)/magDir_;
+
+ scalar radius_sec_inner =
+ innerRadius1_
+ +compInner*(innerRadius2_-innerRadius1_)/magDir_;
+
+ if (parallel < 0)
+ {
+ // Left of point1 endcap
+ volType[pointI] = volumeType::OUTSIDE;
+ }
+ else if (parallel > magDir_)
+ {
+ // Right of point2 endcap
+ volType[pointI] = volumeType::OUTSIDE;
+ }
+ else
+ {
+ // Remove the parallel component
+ v -= parallel*unitDir_;
+ if (mag(v) >= radius_sec_inner && mag(v) <= radius_sec)
+ {
+ volType[pointI] = volumeType::INSIDE;
+ }
+ else
+ {
+ volType[pointI] = volumeType::OUTSIDE;
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurface/searchableCone.H b/src/meshTools/searchableSurface/searchableCone.H
new file mode 100644
index 0000000000000000000000000000000000000000..6e477bc4057eabdb083fe4d762653df3f83ef8dd
--- /dev/null
+++ b/src/meshTools/searchableSurface/searchableCone.H
@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::searchableCone
+
+Description
+ Searching on (optionally hollow) cone.
+
+ \heading Function object usage
+ \table
+ Property | Description | Required | Default value
+ point1 | coordinate of endpoint | yes |
+ radius1 | radius at point1 | yes | yes
+ innerRadius1 | inner radius at point1 | no |
+ point2 | coordinate of endpoint | yes |
+ radius2 | radius at point2 | yes | yes
+ innerRadius2 | inner radius at point2 | no |
+ \endtable
+
+Note
+ Initial implementation, might suffer from robustness (e.g. radius1==radius2)
+
+SourceFiles
+ searchableCone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableCone_H
+#define searchableCone_H
+
+#include "searchableSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class searchableCone Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableCone
+:
+ public searchableSurface
+{
+ // Private Member Data
+
+ //- 'Left' point
+ const point point1_;
+
+ //- Outer radius at point1
+ const scalar radius1_;
+
+ //- Inner radius at point1
+ const scalar innerRadius1_;
+
+
+ //- 'Right' point
+ const point point2_;
+
+ //- Outer radius at point2
+ const scalar radius2_;
+
+ //- Inner radius at point2
+ const scalar innerRadius2_;
+
+
+ //- Length of vector point2-point1
+ const scalar magDir_;
+
+ //- Normalised vector point2-point1
+ const vector unitDir_;
+
+
+ //- Names of regions
+ mutable wordList regions_;
+
+
+ // Private Member Functions
+
+ //- Find nearest point on cylinder.
+ void findNearest
+ (
+ const point& sample,
+ const scalar nearestDistSqr,
+ pointIndexHit & nearInfo,
+ vector& normal
+ ) const;
+
+ //- Determine radial coordinate (squared)
+ static scalar radius2(const searchableCone& cone, const point& pt);
+
+ //- Find both intersections with cone. innerRadii supplied externally.
+ void findLineAll
+ (
+ const searchableCone& cone,
+ const scalar innerRadius1,
+ const scalar innerRadius2,
+ const point& start,
+ const point& end,
+ pointIndexHit& near,
+ pointIndexHit& far
+ ) const;
+
+ //- Insert a hit if it differs (by a tolerance) from the existing ones
+ void insertHit
+ (
+ const point& start,
+ const point& end,
+ List<pointIndexHit>& info,
+ const pointIndexHit& hit
+ ) const;
+
+ //- Return the boundBox of the cylinder
+ boundBox calcBounds() const;
+
+ //- Disallow default bitwise copy construct
+ searchableCone(const searchableCone&);
+
+ //- Disallow default bitwise assignment
+ void operator=(const searchableCone&);
+
+
+public:
+
+ //- Runtime type information
+ TypeName("searchableCone");
+
+
+ // Constructors
+
+ //- Construct from components
+ searchableCone
+ (
+ const IOobject& io,
+ const point& point1,
+ const scalar radius1,
+ const scalar innerRadius1,
+ const point& point2,
+ const scalar radius2,
+ const scalar innerRadius2
+ );
+
+ //- Construct from dictionary (used by searchableSurface)
+ searchableCone
+ (
+ const IOobject& io,
+ const dictionary& dict
+ );
+
+
+ //- Destructor
+
+ virtual ~searchableCone();
+
+
+ // Member Functions
+
+ virtual const wordList& regions() const;
+
+ //- Whether supports volume type below
+ virtual bool hasVolumeType() const
+ {
+ return true;
+ }
+
+ //- Range of local indices that can be returned.
+ virtual label size() const
+ {
+ return 1;
+ }
+
+ //- Get representative set of element coordinates
+ // Usually the element centres (should be of length size()).
+ virtual tmp<pointField> coordinates() const;
+
+ //- Get bounding spheres (centre and radius squared), one per element.
+ // Any point on element is guaranteed to be inside.
+ virtual void boundingSpheres
+ (
+ pointField& centres,
+ scalarField& radiusSqr
+ ) const;
+
+ //- Get the points that define the surface.
+ virtual tmp<pointField> points() const;
+
+ //- Does any part of the surface overlap the supplied bound box?
+ virtual bool overlaps(const boundBox& bb) const
+ {
+ notImplemented
+ (
+ "searchableCone::overlaps(const boundBox&) const"
+ );
+
+ return false;
+ }
+
+
+ // Multiple point queries.
+
+ //- Find nearest point on cylinder
+ virtual void findNearest
+ (
+ const pointField& sample,
+ const scalarField& nearestDistSqr,
+ List<pointIndexHit>&
+ ) const;
+
+ //- Find nearest intersection on line from start to end
+ virtual void findLine
+ (
+ const pointField& start,
+ const pointField& end,
+ List<pointIndexHit>&
+ ) const;
+
+ //- Find any intersection on line from start to end
+ virtual void findLineAny
+ (
+ const pointField& start,
+ const pointField& end,
+ List<pointIndexHit>&
+ ) const;
+
+ //- Find all intersections in order from start to end
+ virtual void findLineAll
+ (
+ const pointField& start,
+ const pointField& end,
+ List<List<pointIndexHit> >&
+ ) const;
+
+ //- From a set of points and indices get the region
+ virtual void getRegion
+ (
+ const List<pointIndexHit>&,
+ labelList& region
+ ) const;
+
+ //- From a set of points and indices get the normal
+ virtual void getNormal
+ (
+ const List<pointIndexHit>&,
+ vectorField& normal
+ ) const;
+
+ //- Determine type (inside/outside/mixed) for point. unknown if
+ // cannot be determined (e.g. non-manifold surface)
+ virtual void getVolumeType
+ (
+ const pointField&,
+ List<volumeType>&
+ ) const;
+
+
+ // regIOobject implementation
+
+ bool writeData(Ostream&) const
+ {
+ notImplemented("searchableCone::writeData(Ostream&) const");
+ return false;
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //