ENH: adjustments to searchable surfaces

- code reduction, documentation, code stubs for spheroid (#1901)

- make searchableSurfaceCollection available as 'collection'
  for consistency with other objects

applications/test/searchableSphere/Make/files

@@ -0,0 +1,3 @@
+Test-searchableSphere.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-searchableSphere

applications/test/searchableSphere/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfileFormats \
+    -lsurfMesh \
+    -lmeshTools

applications/test/searchableSphere/Test-searchableSphere.C

@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Application
+    Test-searchableSphere
+
+Description
+    Basic tests for searchable sphere
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "searchableSphere.H"
+#include "unitConversion.H"
+#include "Random.H"
+
+using namespace Foam;
+using namespace Foam::constant::mathematical;
+
+//- Evaluate using implicit form of the spheroid equation.
+scalar evaluate(const point& p, const searchableSphere& sph)
+{
+    return
+    (
+        sqr((p.x() - sph.centre().x()) / sph.radii().x())
+      + sqr((p.y() - sph.centre().y()) / sph.radii().y())
+      + sqr((p.z() - sph.centre().z()) / sph.radii().z())
+    );
+}
+
+
+void doTest1(const searchableSphere& sph)
+{
+    Info<< nl << "origin:" << sph.centre() << " radius:" << sph.radius();
+
+    if (sph.shape() == searchableSphere::SPHERE)
+    {
+        Info<< " [sphere]" << nl;
+    }
+    else
+    {
+        Info<< " radii:" << sph.radii() << nl;
+    }
+
+//     boundBox bb(point(0, 0, 0), point(30, 30, 30));
+//     Info<< "overlaps: " << Switch(sph.overlaps(bb)) << endl;
+
+    Pair<scalar> angles(-pi/2, pi/2);
+
+    point surfPt = sph.surfacePoint(angles.first(), angles.second());
+    vector norm = sph.surfaceNormal(angles.first(), angles.second());
+
+    Info<< "point at ("
+        << radToDeg(angles.first()) << ' '
+        << radToDeg(angles.second()) << ") deg" << nl;
+
+    Info<< "surface" << nl
+        << "    eval: " << evaluate(surfPt, sph) << nl
+        << "   point:" << surfPt << nl
+        << "  normal:" << norm << nl;
+
+    {
+        List<pointIndexHit> hits(1);
+        vectorField norms;
+
+        hits[0].hitPoint(surfPt, 0);
+
+        sph.getNormal(hits, norms);
+
+        Info<< "  normal:" << norms[0] << nl;
+    }
+
+
+    Random rndGen;
+
+    point testPt1 =
+    (
+        surfPt + 0.1*sph.radius()*norm
+      + 0.01*rndGen.sample01<vector>()
+    );
+
+
+    // Scale by max radius and shift by origin
+    const auto adjustPoint =
+        [&](point& p) -> void
+        {
+            p *= sph.radius();
+            p += sph.centre();
+        };
+
+
+    List<pointIndexHit> hits;
+    pointField query
+    ({
+        testPt1,
+        point(-2, -2, -2)
+    });
+
+    forAll(query, pointi)
+    {
+        if (pointi) adjustPoint(query[pointi]);
+    }
+
+    sph.findNearest
+    (
+        query,
+        scalarField(query.size(), GREAT),
+        hits
+    );
+
+    Info<< "query:" << nl;
+
+    forAll(hits, pointi)
+    {
+        // Should never miss
+        Info<< "  "
+            << query[pointi] << " => "
+            << hits[pointi].hitPoint() << nl;
+    }
+
+
+    pointField lineBeg
+    ({
+        point(-2, -2, -2),
+        point(0,0,0),
+        point(0,0,0),
+        point(0,0,0)
+    });
+    pointField lineEnd
+    ({
+        point(2, 2, 2),
+        point(2, 0, 0),
+        point(0, 2, 0),
+        point(0, 0, 4)
+    });
+
+    for (point& p : lineBeg)
+    {
+        adjustPoint(p);
+    }
+
+    for (point& p : lineEnd)
+    {
+        adjustPoint(p);
+    }
+
+    sph.findLineAny
+    (
+        lineBeg,
+        lineEnd,
+        hits
+    );
+
+    Info<< "lines:" << nl;
+
+    forAll(hits, pointi)
+    {
+        // Should never miss
+        Info<< "  "
+            << lineBeg[pointi] << " -> "
+            << lineEnd[pointi] << " = "
+            << hits[pointi].hitPoint() << nl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+    argList::noBanner();
+    #include "setRootCase.H"
+
+    // Use dummy Time for objectRegistry
+    autoPtr<Time> dummyTimePtr(Time::New());
+
+    const IOobject io
+    (
+        "sphere",
+        *dummyTimePtr,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE,
+        false  // do not register
+    );
+
+    Info<< "Testing searchable sphere" << endl;
+
+    doTest1
+    (
+        searchableSphere
+        (
+            io,
+            point(0.5, 0.5, 0.5),
+            scalar(2)
+        )
+    );
+
+    Info<< "\nDone\n" << endl;
+    return 0;
+}
+
+// ************************************************************************* //

src/meshTools/searchableSurfaces/searchableBox/searchableBox.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,7 +77,8 @@ void Foam::searchableBox::projectOntoCoordPlane
         FatalErrorInFunction
             << "Point on plane " << planePt
             << " is not on coordinate " << min()[dir]
-            << " nor " << max()[dir] << abort(FatalError);
+            << " nor " << max()[dir] << nl
+            << abort(FatalError);
     }
 }
 
@@ -135,31 +136,26 @@ Foam::pointIndexHit Foam::searchableBox::findNearest
     // using the three near distances. Project onto the nearest plane.
     if (!outside)
     {
-        vector dist(cmptMag(info.rawPoint() - near));
+        const vector dist(cmptMag(info.point() - near));
+
+        direction projNorm(vector::Z);
 
         if (dist.x() < dist.y())
         {
             if (dist.x() < dist.z())
             {
-                // Project onto x plane
+                projNorm = vector::X;
-                projectOntoCoordPlane(vector::X, near, info);
-            }
-            else
-            {
-                projectOntoCoordPlane(vector::Z, near, info);
             }
         }
         else
         {
             if (dist.y() < dist.z())
             {
-                projectOntoCoordPlane(vector::Y, near, info);
+                projNorm = vector::Y;
-            }
-            else
-            {
-                projectOntoCoordPlane(vector::Z, near, info);
             }
         }
+
+        projectOntoCoordPlane(projNorm, near, info);
     }
 
 
@@ -194,7 +190,7 @@ Foam::searchableBox::searchableBox
             << exit(FatalError);
     }
 
-    bounds() = static_cast<boundBox>(*this);
+    bounds() = static_cast<treeBoundBox>(*this);
 }
 
 
@@ -204,19 +200,12 @@ Foam::searchableBox::searchableBox
     const dictionary& dict
 )
 :
-    searchableSurface(io),
+    searchableBox
-    treeBoundBox(dict.get<point>("min"), dict.get<point>("max"))
+    (
-{
+        io,
-    if (!treeBoundBox::valid())
+        treeBoundBox(dict.get<point>("min"), dict.get<point>("max"))
-    {
+    )
-        FatalErrorInFunction
+{}
-            << "Illegal bounding box specification : "
-            << static_cast<const treeBoundBox>(*this) << nl
-            << exit(FatalError);
-    }
-
-    bounds() = static_cast<boundBox>(*this);
-}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -257,7 +246,7 @@ void Foam::searchableBox::boundingSpheres
 {
     centres.setSize(size());
     radiusSqr.setSize(size());
-    radiusSqr = 0.0;
+    radiusSqr = Zero;
 
     const pointField pts(treeBoundBox::points());
     const faceList& fcs = treeBoundBox::faces;
@@ -527,8 +516,7 @@ void Foam::searchableBox::findLineAll
     const scalarField magSqrDirVec(magSqr(dirVec));
     const vectorField smallVec
     (
-        ROOTSMALL*dirVec
+        ROOTSMALL*dirVec + vector::uniform(ROOTVSMALL)
-      + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
     );
 
     forAll(start, pointi)

src/meshTools/searchableSurfaces/searchableBox/searchableBox.H

@@ -33,11 +33,14 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                       | Required | Default
-        type        | box / searchableBox               | selector |
+        type        | box                               | selector |
         min         | minimum point for bounding box    | yes   |
         max         | maximum point for bounding box    | yes   |
     \endtable
 
+Note
+    Longer type name : \c searchableBox
+
 SourceFiles
     searchableBox.C
 
@@ -115,6 +118,7 @@ public:
             const dictionary& dict
         );
 
+
     //- Destructor
     virtual ~searchableBox() = default;
 

src/meshTools/searchableSurfaces/searchableCone/searchableCone.C

@@ -122,11 +122,11 @@ void Foam::searchableCone::findNearestAndNormal
 
     // 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);
+    point nearCone(point::uniform(GREAT));
     vector normalCone(1, 0, 0);
 
     // Nearest and normal on inner cone. Initialise to far-away point
-    point iCnearCone(GREAT, GREAT, GREAT);
+    point iCnearCone(point::uniform(GREAT));
     vector iCnormalCone(1, 0, 0);
 
     point projPt1 = point1_+ radius1_*v;

src/meshTools/searchableSurfaces/searchableCone/searchableCone.H

@@ -32,7 +32,7 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                       | Required | Default
-        type        | cone / searchableCone             | selector |
+        type        | cone                              | selector |
         point1      | coordinate of endpoint            | yes   |
         radius1     | radius at point1                  | yes   |
         innerRadius1| inner radius at point1            | no    | 0
@@ -44,6 +44,9 @@ Description
 Note
     Initial implementation, might suffer from robustness (e.g. radius1==radius2)
 
+Note
+    Longer type name : \c searchableCone
+
 SourceFiles
     searchableCone.C
 

src/meshTools/searchableSurfaces/searchableCylinder/searchableCylinder.H

@@ -33,12 +33,15 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                       | Required | Default
-        type        | cylinder / searchableCylinder     | selector |
+        type        | cylinder                          | selector |
         point1      | coordinate of endpoint            | yes   |
         point2      | coordinate of endpoint            | yes   |
         radius      | cylinder radius                   | yes   |
     \endtable
 
+Note
+    Longer type name : \c searchableCylinder
+
 SourceFiles
     searchableCylinder.C
 

src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.H

@@ -34,12 +34,15 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                       | Required | Default
-        type        | disk / searchableDisk             | selector |
+        type        | disk                              | selector |
         origin      | centre of disk                    | yes   |
         normal      | normal vector                     | yes   |
         radius      | disk radius                       | yes   |
     \endtable
 
+Note
+    Longer type name : \c searchableDisk
+
 SourceFiles
     searchableDisk.C
 

src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.H

@@ -32,7 +32,7 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                       | Required | Default
-        type        | extrudedCircle / searchableExtrudedCircle | selector |
+        type        | extrudedCircle                    | selector |
         file        | The name of the edge mesh         | yes   |
         radius      | Search radius around the edges    | yes   |
     \endtable
@@ -42,6 +42,9 @@ Note
     - Can not be used with snappyHexMesh since only implements nearest
       searching.
 
+Note
+    Longer type name : \c searchableExtrudedCircle
+
 SourceFiles
     searchableExtrudedCircle.C
 

src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.C

@@ -80,21 +80,19 @@ Foam::pointIndexHit Foam::searchablePlane::findLine
 
 Foam::boundBox Foam::searchablePlane::calcBounds() const
 {
-    point max(VGREAT, VGREAT, VGREAT);
+    boundBox bb(boundBox::greatBox);
 
     for (direction dir = 0; dir < vector::nComponents; dir++)
     {
         if (mag(normal()[dir]) - 1 < SMALL)
         {
-            max[dir] = 0;
+            bb.min()[dir] = 0;
-
+            bb.max()[dir] = 0;
             break;
         }
     }
 
-    point min = -max;
+    return bb;
-
-    return boundBox(min, max);
 }
 
 
@@ -146,10 +144,10 @@ void Foam::searchablePlane::boundingSpheres
     scalarField& radiusSqr
 ) const
 {
-    centres.setSize(1);
+    centres.resize(1);
-    centres[0] = origin();
+    radiusSqr.resize(1);
 
-    radiusSqr.setSize(1);
+    centres[0] = origin();
     radiusSqr[0] = Foam::sqr(GREAT);
 }
 

src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,9 +35,12 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                       | Required  | Default
-        type        | plane / searchablePlane           | selector |
+        type        | plane                             | selector |
     \endtable
 
+Note
+    Longer type name : \c searchablePlane
+
 SourceFiles
     searchablePlane.C
 
@@ -111,6 +115,7 @@ public:
             const dictionary& dict
         );
 
+
     //- Destructor
     virtual ~searchablePlane() = default;
 

src/meshTools/searchableSurfaces/searchablePlate/searchablePlate.C

@@ -61,30 +61,29 @@ Foam::direction Foam::searchablePlate::calcNormal(const point& span)
     {
         if (span[dir] < 0)
         {
-            FatalErrorInFunction
+            // Negative entry. Flag and exit.
-                << "Span should have two positive and one zero entry. Now:"
+            normalDir = 3;
-                << span << exit(FatalError);
+            break;
         }
         else if (span[dir] < VSMALL)
         {
-            if (normalDir == 3)
+            if (normalDir != 3)
-            {
-                normalDir = dir;
-            }
-            else
             {
                 // Multiple zero entries. Flag and exit.
                 normalDir = 3;
                 break;
             }
+
+            normalDir = dir;
         }
     }
 
     if (normalDir == 3)
     {
         FatalErrorInFunction
-            << "Span should have two positive and one zero entry. Now:"
+            << "Span should have two positive and one zero entry: "
-            << span << exit(FatalError);
+            << span << nl
+            << exit(FatalError);
     }
 
     return normalDir;
@@ -246,18 +245,8 @@ Foam::searchablePlate::searchablePlate
     const dictionary& dict
 )
 :
-    searchableSurface(io),
+    searchablePlate(io, dict.get<point>("origin"), dict.get<vector>("span"))
-    origin_(dict.get<point>("origin")),
+{}
-    span_(dict.get<vector>("span")),
-    normalDir_(calcNormal(span_))
-{
-    DebugInFunction
-        << " origin:" << origin_
-        << " origin+span:" << origin_+span_
-        << " normal:" << vector::componentNames[normalDir_] << nl;
-
-    bounds() = boundBox(origin_, origin_ + span_);
-}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -266,8 +255,8 @@ const Foam::wordList& Foam::searchablePlate::regions() const
 {
     if (regions_.empty())
     {
-        regions_.setSize(1);
+        regions_.resize(1);
-        regions_[0] = "region0";
+        regions_.first() = "region0";
     }
     return regions_;
 }
@@ -285,10 +274,10 @@ void Foam::searchablePlate::boundingSpheres
     scalarField& radiusSqr
 ) const
 {
-    centres.setSize(1);
+    centres.resize(1);
-    centres[0] = origin_ + 0.5*span_;
+    radiusSqr.resize(1);
 
-    radiusSqr.setSize(1);
+    centres[0] = origin_ + 0.5*span_;
     radiusSqr[0] = Foam::magSqr(0.5*span_);
 
     // Add a bit to make sure all points are tested inside
@@ -298,26 +287,25 @@ void Foam::searchablePlate::boundingSpheres
 
 Foam::tmp<Foam::pointField> Foam::searchablePlate::points() const
 {
-    auto tpts = tmp<pointField>::New(4);
+    auto tpts = tmp<pointField>::New(4, origin_);
     auto& pts = tpts.ref();
 
-    pts[0] = origin_;
+    pts[2] += span_;
-    pts[2] = origin_ + span_;
 
     if (span_.x() < span_.y() && span_.x() < span_.z())
     {
-        pts[1] = origin_ + point(0, span_.y(), 0);
+        pts[1].y() += span_.y();
-        pts[3] = origin_ + point(0, 0, span_.z());
+        pts[3].z() += span_.z();
     }
     else if (span_.y() < span_.z())
     {
-        pts[1] = origin_ + point(span_.x(), 0, 0);
+        pts[1].x() += span_.x();
-        pts[3] = origin_ + point(0, 0, span_.z());
+        pts[3].z() += span_.z();
     }
     else
     {
-        pts[1] = origin_ + point(span_.x(), 0, 0);
+        pts[1].x() += span_.x();
-        pts[3] = origin_ + point(0, span_.y(), 0);
+        pts[3].y() += span_.y();
     }
 
     return tpts;
@@ -326,15 +314,7 @@ Foam::tmp<Foam::pointField> Foam::searchablePlate::points() const
 
 bool Foam::searchablePlate::overlaps(const boundBox& bb) const
 {
-    return
+    return bb.overlaps(bounds());
-    (
-        (origin_.x() + span_.x()) >= bb.min().x()
-     && origin_.x() <= bb.max().x()
-     && (origin_.y() + span_.y()) >= bb.min().y()
-     && origin_.y() <= bb.max().y()
-     && (origin_.z() + span_.z()) >= bb.min().z()
-     && origin_.z() <= bb.max().z()
-    );
 }
 
 

src/meshTools/searchableSurfaces/searchablePlate/searchablePlate.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,11 +43,14 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                       | Required | Default
-        type        | plate / searchablePlate           | selector |
+        type        | plate                             | selector |
-        origin      | centre of the plate               | yes   |
+        origin      | The minimum corner of the plate   | yes |
-        span        | The plate dimensions              | yes   |
+        span        | The plate dimensions              | yes |
     \endtable
 
+Note
+    Longer type name : \c searchablePlate
+
 SourceFiles
     searchablePlate.C
 
@@ -72,12 +75,12 @@ class searchablePlate
 :
     public searchableSurface
 {
-private:
-
     // Private Member Data
 
+        //- The minimum corner of the plate
         const point origin_;
 
+        //- The span size of the plate
         const vector span_;
 
         //- Coordinate direction which is normal
@@ -144,6 +147,12 @@ public:
 
     // Member Functions
 
+        //- The centre (origin) of the plate
+        inline const point& centre() const noexcept
+        {
+             return origin_;
+        }
+
         //- Names of regions
         virtual const wordList& regions() const;
 

src/meshTools/searchableSurfaces/searchableRotatedBox/searchableRotatedBox.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014 OpenFOAM Foundation
-    Copyright (C) 2015-2018 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,7 +55,8 @@ namespace Foam
 Foam::searchableRotatedBox::searchableRotatedBox
 (
     const IOobject& io,
-    const dictionary& dict
+    const vector& span,
+    const coordSystem::cartesian& csys
 )
 :
     searchableSurface(io),
@@ -69,21 +70,36 @@ Foam::searchableRotatedBox::searchableRotatedBox
             io.db(),
             io.readOpt(),
             io.writeOpt(),
-            false      //io.registerObject(),
+            false      // never register
         ),
-        treeBoundBox(Zero, dict.get<vector>("span"))
+        treeBoundBox(Zero, span)
     ),
-    transform_
+    transform_(csys.origin(), csys.e3(), csys.e1())
-    (
-        dict.get<point>("origin"),
-        dict.get<vector>("e3"),
-        dict.get<vector>("e1")
-    )
 {
     points_ = transform_.globalPosition(box_.points());
 }
 
 
+Foam::searchableRotatedBox::searchableRotatedBox
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableRotatedBox
+    (
+        io,
+        dict.get<vector>("span"),
+        coordSystem::cartesian
+        (
+            dict.get<point>("origin"),
+            dict.get<vector>("e3"),
+            dict.get<vector>("e1")
+        )
+    )
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 const Foam::wordList& Foam::searchableRotatedBox::regions() const
@@ -125,14 +141,10 @@ bool Foam::searchableRotatedBox::overlaps(const boundBox& bb) const
         return false;
     }
 
-    // 2. Check if one or more face points inside
+    // 2. Check if any corner points inside
-    const faceList& fcs = treeBoundBox::faces;
+    if (bb.containsAny(points_))
-    forAll(fcs, faceI)
     {
-        if (bb.containsAny(points_))
+        return true;
-        {
-            return true;
-        }
     }
 
     // 3. Difficult case: all points are outside but connecting edges might
@@ -141,8 +153,7 @@ bool Foam::searchableRotatedBox::overlaps(const boundBox& bb) const
     const treeBoundBox treeBb(bb);
 
     // 3a. my edges through bb faces
-    const edgeList& edges = treeBoundBox::edges;
+    for (const edge& e : treeBoundBox::edges)
-    for (const edge& e : edges)
     {
         point inter;
         if (treeBb.intersects(points_[e[0]], points_[e[1]], inter))
@@ -155,11 +166,11 @@ bool Foam::searchableRotatedBox::overlaps(const boundBox& bb) const
 
     const pointField bbPoints(bb.points());
 
-    for (const face& f : fcs)
+    for (const face& f : treeBoundBox::faces)
     {
-        point fc = f.centre(points_);
+        const point fc = f.centre(points_);
 
-        for (const edge& e : edges)
+        for (const edge& e : treeBoundBox::edges)
         {
             pointHit inter = f.intersection
             (

src/meshTools/searchableSurfaces/searchableRotatedBox/searchableRotatedBox.H

@@ -46,13 +46,16 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                       | Required | Default
-        type        | rotatedBox / searchableRotatedBox | selector |
+        type        | rotatedBox                        | selector |
         span        | The box dimensions                | yes   |
         origin      | The box corner                    | yes   |
         e1          | Local x-axis of the box           | yes   |
         e3          | Local z-axis of the box           | yes   |
     \endtable
 
+Note
+    Longer type name : \c searchableRotatedBox
+
 SourceFiles
     searchableRotatedBox.C
 
@@ -78,8 +81,6 @@ class searchableRotatedBox
 :
     public searchableSurface
 {
-private:
-
     // Private Member Data
 
         //- Box in local coordinate system
@@ -109,6 +110,14 @@ public:
 
     // Constructors
 
+        //- Construct from components
+        searchableRotatedBox
+        (
+            const IOobject& io,
+            const vector& span,
+            const coordSystem::cartesian& csys = coordSystem::cartesian{}
+        );
+
         //- Construct from dictionary (used by searchableSurface)
         searchableRotatedBox
         (

src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,95 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+// General handling
+namespace Foam
+{
+
+// Dictionary entry with single scalar or vector quantity
+inline static vector getRadius(const word& name, const dictionary& dict)
+{
+    if (token(dict.lookup(name)).isNumber())
+    {
+        return vector::uniform(dict.get<scalar>(name));
+    }
+
+    return dict.get<vector>(name);
+}
+
+
+// Test point for negative components, return the sign-changes
+inline static unsigned getOctant(const point& p)
+{
+    unsigned octant = 0;
+
+    if (p.x() < 0) { octant |= 1; }
+    if (p.y() < 0) { octant |= 2; }
+    if (p.z() < 0) { octant |= 4; }
+
+    return octant;
+}
+
+
+// Apply sign-changes to point
+inline static void applyOctant(point& p, unsigned octant)
+{
+    if (octant & 1) { p.x() = -p.x(); }
+    if (octant & 2) { p.y() = -p.y(); }
+    if (octant & 4) { p.z() = -p.z(); }
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+inline Foam::searchableSphere::componentOrder
+Foam::searchableSphere::getOrdering(const vector& radii)
+{
+    #ifdef FULLDEBUG
+    for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
+    {
+        if (radii[cmpt] <= 0)
+        {
+            FatalErrorInFunction
+                << "Radii must be positive, non-zero: " << radii << endl
+                << exit(FatalError);
+        }
+    }
+    #endif
+
+    std::array<direction, 3> idx{0, 1, 2};
+
+    // Reverse sort by magnitude (largest first...)
+    // Radii are positive (checked above, or just always true)
+    std::stable_sort
+    (
+        idx.begin(),
+        idx.end(),
+        [&](direction a, direction b){ return radii[a] > radii[b]; }
+    );
+
+    componentOrder order{idx[0], idx[1], idx[2], shapeType::GENERAL};
+
+    if (equal(radii[order.major], radii[order.minor]))
+    {
+        order.shape = shapeType::SPHERE;
+    }
+    else if (equal(radii[order.major], radii[order.mezzo]))
+    {
+        order.shape = shapeType::OBLATE;
+    }
+    else if (equal(radii[order.mezzo], radii[order.minor]))
+    {
+        order.shape = shapeType::PROLATE;
+    }
+
+    return order;
+}
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 Foam::pointIndexHit Foam::searchableSphere::findNearest
@@ -60,23 +149,33 @@ Foam::pointIndexHit Foam::searchableSphere::findNearest
 {
     pointIndexHit info(false, sample, -1);
 
-    const vector n(sample - origin_);
+    // Handle special cases first
-    scalar magN = mag(n);
 
-    if (nearestDistSqr >= sqr(magN - radius_))
+    // if (order_.shape == shapeType::SPHERE)
+    if (true)
     {
-        if (magN < ROOTVSMALL)
+        // Point relative to origin, simultaneously the normal on the sphere
+        const vector n(sample - origin_);
+        const scalar magN = mag(n);
+
+        // It is a sphere, all radii are identical
+
+        if (nearestDistSqr >= sqr(magN - radii_[0]))
         {
-            info.rawPoint() = origin_ + vector(1,0,0)*radius_;
+            info.setPoint
+            (
+                origin_
+              + (magN < ROOTVSMALL ? vector(radii_[0],0,0) : (radii_[0]*n/magN))
+            );
+            info.setHit();
+            info.setIndex(0);
         }
-        else
+
-        {
+        return info;
-            info.rawPoint() = origin_ + n/magN*radius_;
-        }
-        info.setHit();
-        info.setIndex(0);
     }
 
+    //[code]
+
     return info;
 }
 
@@ -93,43 +192,47 @@ void Foam::searchableSphere::findLineAll
     near.setMiss();
     far.setMiss();
 
-    vector dir(end-start);
+    // if (order_.shape == shapeType::SPHERE)
-    scalar magSqrDir = magSqr(dir);
+    if (true)
-
-    if (magSqrDir > ROOTVSMALL)
     {
-        const vector toCentre(origin_ - start);
+        vector dir(end-start);
-        scalar magSqrToCentre = magSqr(toCentre);
+        const scalar magSqrDir = magSqr(dir);
 
-        dir /= Foam::sqrt(magSqrDir);
+        if (magSqrDir > ROOTVSMALL)
-
-        scalar v = (toCentre & dir);
-
-        scalar disc = sqr(radius_) - (magSqrToCentre - sqr(v));
-
-        if (disc >= 0)
         {
-            scalar d = Foam::sqrt(disc);
+            dir /= Foam::sqrt(magSqrDir);
 
-            scalar nearParam = v-d;
+            const vector relStart(start - origin_);
 
-            if (nearParam >= 0 && sqr(nearParam) <= magSqrDir)
+            const scalar v = -(relStart & dir);
+
+            const scalar disc = sqr(radius()) - (magSqr(relStart) - sqr(v));
+
+            if (disc >= 0)
             {
-                near.setHit();
+                const scalar d = Foam::sqrt(disc);
-                near.setPoint(start + nearParam*dir);
-                near.setIndex(0);
-            }
 
-            scalar farParam = v+d;
+                const scalar nearParam = v - d;
+                const scalar farParam = v + d;
 
-            if (farParam >= 0 && sqr(farParam) <= magSqrDir)
+                if (nearParam >= 0 && sqr(nearParam) <= magSqrDir)
-            {
+                {
-                far.setHit();
+                    near.setHit();
-                far.setPoint(start + farParam*dir);
+                    near.setPoint(start + nearParam*dir);
-                far.setIndex(0);
+                    near.setIndex(0);
+                }
+
+                if (farParam >= 0 && sqr(farParam) <= magSqrDir)
+                {
+                    far.setHit();
+                    far.setPoint(start + farParam*dir);
+                    far.setIndex(0);
+                }
             }
         }
     }
+
+    //[code]
 }
 
 
@@ -141,16 +244,26 @@ Foam::searchableSphere::searchableSphere
     const point& origin,
     const scalar radius
 )
+:
+    searchableSphere(io, origin, vector::uniform(radius))
+{}
+
+
+Foam::searchableSphere::searchableSphere
+(
+    const IOobject& io,
+    const point& origin,
+    const vector& radii
+)
 :
     searchableSurface(io),
     origin_(origin),
-    radius_(radius)
+    // radii_(radii),
+    radii_(vector::uniform(cmptMax(radii))) /* Transition */,
+    order_{getOrdering(radii_)}
 {
-    bounds() = boundBox
+    bounds().min() = (centre() - radii_);
-    (
+    bounds().max() = (centre() + radii_);
-        origin_ - radius_*vector::one,
-        origin_ + radius_*vector::one
-    );
 }
 
 
@@ -164,16 +277,61 @@ Foam::searchableSphere::searchableSphere
     (
         io,
         dict.getCompat<vector>("origin", {{"centre", -1806}}),
-        dict.get<scalar>("radius")
+        getRadius("radius", dict)
     )
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::point Foam::searchableSphere::surfacePoint
+(
+    const scalar theta,
+    const scalar phi
+) const
+{
+    return point
+    (
+        origin_.x() + radii_.x() * cos(theta)*sin(phi),
+        origin_.y() + radii_.y() * sin(theta)*sin(phi),
+        origin_.z() + radii_.z() * cos(phi)
+    );
+}
+
+
+Foam::vector Foam::searchableSphere::surfaceNormal
+(
+    const scalar theta,
+    const scalar phi
+) const
+{
+    // Normal is (x0/r0^2, x1/r1^2, x2/r2^2)
+
+    return vector
+    (
+        cos(theta)*sin(phi) / radii_.x(),
+        sin(theta)*sin(phi) / radii_.y(),
+        cos(phi) / radii_.z()
+    ).normalise();
+}
+
+
 bool Foam::searchableSphere::overlaps(const boundBox& bb) const
 {
-    return bb.overlaps(origin_, sqr(radius_));
+    // if (order_.shape == shapeType::SPHERE)
+    if (true)
+    {
+        return bb.overlaps(origin_, sqr(radius()));
+    }
+
+    if (!bb.valid())
+    {
+        return false;
+    }
+
+    //[code]
+
+    return true;
 }
 
 
@@ -188,7 +346,6 @@ const Foam::wordList& Foam::searchableSphere::regions() const
 }
 
 
-
 void Foam::searchableSphere::boundingSpheres
 (
     pointField& centres,
@@ -196,10 +353,10 @@ void Foam::searchableSphere::boundingSpheres
 ) const
 {
     centres.resize(1);
-    centres[0] = origin_;
-
     radiusSqr.resize(1);
-    radiusSqr[0] = Foam::sqr(radius_);
+
+    centres[0] = origin_;
+    radiusSqr[0] = Foam::sqr(radius());
 
     // Add a bit to make sure all points are tested inside
     radiusSqr += Foam::sqr(SMALL);
@@ -213,7 +370,7 @@ void Foam::searchableSphere::findNearest
     List<pointIndexHit>& info
 ) const
 {
-    info.setSize(samples.size());
+    info.resize(samples.size());
 
     forAll(samples, i)
     {
@@ -229,14 +386,17 @@ void Foam::searchableSphere::findLine
     List<pointIndexHit>& info
 ) const
 {
-    info.setSize(start.size());
+    info.resize(start.size());
 
     pointIndexHit b;
 
     forAll(start, i)
     {
-        // Pick nearest intersection. If none intersected take second one.
+        // Pick nearest intersection.
+        // If none intersected take second one.
+
         findLineAll(start[i], end[i], info[i], b);
+
         if (!info[i].hit() && b.hit())
         {
             info[i] = b;
@@ -252,14 +412,17 @@ void Foam::searchableSphere::findLineAny
     List<pointIndexHit>& info
 ) const
 {
-    info.setSize(start.size());
+    info.resize(start.size());
 
     pointIndexHit b;
 
     forAll(start, i)
     {
+        // Pick nearest intersection.
         // Discard far intersection
+
         findLineAll(start[i], end[i], info[i], b);
+
         if (!info[i].hit() && b.hit())
         {
             info[i] = b;
@@ -275,24 +438,25 @@ void Foam::searchableSphere::findLineAll
     List<List<pointIndexHit>>& info
 ) const
 {
-    info.setSize(start.size());
+    info.resize(start.size());
 
     forAll(start, i)
     {
         pointIndexHit near, far;
+
         findLineAll(start[i], end[i], near, far);
 
         if (near.hit())
         {
             if (far.hit())
             {
-                info[i].setSize(2);
+                info[i].resize(2);
                 info[i][0] = near;
                 info[i][1] = far;
             }
             else
             {
-                info[i].setSize(1);
+                info[i].resize(1);
                 info[i][0] = near;
             }
         }
@@ -300,7 +464,7 @@ void Foam::searchableSphere::findLineAll
         {
             if (far.hit())
             {
-                info[i].setSize(1);
+                info[i].resize(1);
                 info[i][0] = far;
             }
             else
@@ -318,7 +482,7 @@ void Foam::searchableSphere::getRegion
     labelList& region
 ) const
 {
-    region.setSize(info.size());
+    region.resize(info.size());
     region = 0;
 }
 
@@ -329,18 +493,18 @@ void Foam::searchableSphere::getNormal
     vectorField& normal
 ) const
 {
-    normal.setSize(info.size());
+    normal.resize(info.size());
-    normal = Zero;
 
     forAll(info, i)
     {
         if (info[i].hit())
         {
-            normal[i] = normalised(info[i].hitPoint() - origin_);
+            normal[i] = normalised(info[i].point() - origin_);
         }
         else
         {
             // Set to what?
+            normal[i] = Zero;
         }
     }
 }
@@ -352,20 +516,26 @@ void Foam::searchableSphere::getVolumeType
     List<volumeType>& volType
 ) const
 {
-    volType.setSize(points.size());
+    volType.resize(points.size());
 
-    const scalar rad2 = sqr(radius_);
+    // if (order_.shape == shapeType::SPHERE)
-
+    if (true)
-    forAll(points, pointi)
     {
-        const point& pt = points[pointi];
+        const scalar rad2 = sqr(radius());
 
-        volType[pointi] =
+        forAll(points, pointi)
-        (
+        {
-            (magSqr(pt - origin_) <= rad2)
+            const point& p = points[pointi];
-          ? volumeType::INSIDE : volumeType::OUTSIDE
+
-        );
+            volType[pointi] =
+            (
+                (magSqr(p - origin_) <= rad2)
+              ? volumeType::INSIDE : volumeType::OUTSIDE
+            );
+        }
     }
+
+    //[code]
 }
 
 

src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,12 +33,15 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                           | Required | Default
-        type        | sphere / searchableSphere             | selector |
+        type        | sphere                                | selector |
-        origin      | The origin (centre) of the sphere     | yes   |
+        origin      | The origin (centre) of the sphere     | yes |
-        radius      | The (outside) radius of sphere        | yes   |
+        radius      | The (outside) radius of sphere        | yes |
-        centre      | Alternative for 'origin'              | no    |
+        centre      | Alternative name for 'origin'         | no  |
     \endtable
 
+Note
+    Longer type name : \c searchableSphere
+
 SourceFiles
     searchableSphere.C
 
@@ -63,15 +66,44 @@ class searchableSphere
 :
     public searchableSurface
 {
+public:
+
+    // Public Types
+
+        //- The type of shape
+        enum shapeType : uint8_t
+        {
+            SPHERE = 0,     //!< Sphere (all components equal)
+            OBLATE = 1,     //!< Oblate (major = mezzo > minor)
+            PROLATE = 2,    //!< Prolate (major > mezzo = minor)
+            GENERAL = 3     //!< General spheroid
+        };
+
+
 private:
 
-    // Private Member Data
+    // Private Types
+
+        //- Component order (largest to smallest)
+        struct componentOrder
+        {
+            uint8_t major;      //!< Component with major radius
+            uint8_t mezzo;      //!< Component with intermediate radius
+            uint8_t minor;      //!< Component with minor radius
+            shapeType shape;    //!< The type of shape
+        };
+
+
+    // Private Data
 
         //- Centre point of the sphere
         const point origin_;
 
-        //- The outer radius of the sphere
+        //- The outer radii of the spheroid
-        const scalar radius_;
+        const vector radii_;
+
+        //- The canonical (sorted) order and shape
+        const struct componentOrder order_;
 
         //- Names of regions
         mutable wordList regions_;
@@ -79,6 +111,9 @@ private:
 
     // Private Member Functions
 
+        //- Determine sorted order and classify the shape
+        inline static componentOrder getOrdering(const vector& radii);
+
         //- Inherit findNearest from searchableSurface
         using searchableSurface::findNearest;
 
@@ -114,12 +149,20 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct a sphere from components
+        searchableSphere
+        (
+            const IOobject& io,
+            const point& origin,
+            const scalar radius
+        );
+
+        //- Construct a spheroid from components
         searchableSphere
         (
             const IOobject& io,
             const point& centre,
-            const scalar radius
+            const vector& radii
         );
 
         //- Construct from dictionary (used by searchableSurface)
@@ -136,6 +179,44 @@ public:
 
     // Member Functions
 
+    // Geometric
+
+        //- The centre (origin) of the sphere
+        const point& centre() const noexcept
+        {
+             return origin_;
+        }
+
+        //- The radius of the sphere, or major radius of the spheroid
+        scalar radius() const noexcept
+        {
+            return radii_[order_.major];
+        }
+
+        //- The radii of the spheroid
+        const vector& radii() const noexcept
+        {
+            return radii_;
+        }
+
+        //- The type of shape
+        enum shapeType shape() const noexcept
+        {
+            return shapeType::SPHERE;
+        }
+
+        //- A point on the sphere at given location
+        //  theta [-pi,pi], phi [0,pi]
+        point surfacePoint(const scalar theta, const scalar phi) const;
+
+        //- Surface normal on the sphere at given location
+        //  theta [-pi,pi], phi [0,pi]
+        vector surfaceNormal(const scalar theta, const scalar phi) const;
+
+
+
+    // Searching
+
         //- Names of regions
         virtual const wordList& regions() const;
 
@@ -228,7 +309,6 @@ public:
             ) const;
 
             //- Determine type (inside/outside/mixed) for point.
-            //  Unknown if cannot be determined (e.g. non-manifold surface)
             virtual void getVolumeType
             (
                 const pointField& points,
@@ -236,14 +316,14 @@ public:
             ) const;
 
 
-        // regIOobject implementation
+    // Output
-
-            bool writeData(Ostream&) const
-            {
-                NotImplemented;
-                return false;
-            }
 
+        // Implementation for regIOobject. NotImplemented
+        bool writeData(Ostream&) const
+        {
+            NotImplemented;
+            return false;
+        }
 };
 
 

src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,6 +43,13 @@ namespace Foam
         searchableSurfaceCollection,
         dict
     );
+    addNamedToRunTimeSelectionTable
+    (
+        searchableSurface,
+        searchableSurfaceCollection,
+        dict,
+        collection
+    );
 }
 
 

src/meshTools/searchableSurfaces/searchableSurfaceCollection/searchableSurfaceCollection.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015 OpenCFD Ltd.
+    Copyright (C) 2015-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,10 +34,13 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                       | Required | Default
-        type        | searchableSurfaceCollection       | selector |
+        type        | collection                        | selector |
         mergeSubRegions | boolean                       | yes   |
     \endtable
 
+Note
+    Longer type name : \c searchableSurfaceCollection
+
 SourceFiles
     searchableSurfaceCollection.C
 

tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict

@@ -44,7 +44,7 @@ geometry
 
     twoFridgeFreezers
     {
-        type searchableSurfaceCollection;
+        type collection;
 
         mergeSubRegions true;
 

tutorials/incompressible/pimpleFoam/laminar/cylinder2D/system/snappyHexMeshDict

@@ -22,7 +22,7 @@ geometry
 {
     cylinder
     {
-        type    searchableCylinder;
+        type    cylinder;
         point1  (0 0 -0.00375);
         point2  (0 0 0.00375);
         radius  0.06;

tutorials/incompressible/pimpleFoam/laminar/mixerVesselAMI2D/mixerVesselAMI2D-topologyChange/system/blockMeshDict.m4

@@ -1083,7 +1083,7 @@ boundary
 /* optional
         surface
         {
-            type            searchableCylinder; // none
+            type            cylinder; // none
             point1          (0 0 -1);
             point2          (0 0 1);
             radius          0.5;
@@ -1111,7 +1111,7 @@ boundary
 /* optional
         surface
         {
-            type            searchableCylinder; // none
+            type            cylinder; // none
             point1          (0 0 -1);
             point2          (0 0 1);
             radius          0.5;

tutorials/mesh/snappyHexMesh/iglooWithFridgesDirectionalRefinement/system/snappyHexMeshDict

@@ -30,21 +30,21 @@ geometry
 {
     igloo
     {
-        type searchableSphere;
+        type   sphere;
         origin (3 3 0);
         radius 4;
     }
 
     box1
     {
-        type searchableBox;
+        type box;
-        min (0 0 0);
+        min  (0 0 0);
-        max (1 1 1);
+        max  (1 1 1);
     }
 
     twoFridgeFreezers
     {
-        type searchableSurfaceCollection;
+        type collection;
 
         mergeSubRegions true;
 
@@ -87,13 +87,13 @@ geometry
     // For directional refinement
     dirRefineBox1
     {
-        type searchableBox;
+        type box;
         min (-10 -10 -10);
         max ( 10  10  0.5);
     }
     dirRefineBox2
     {
-        type searchableBox;
+        type box;
         min (-10 -10 -10);
         max ( 10  10  0.25);
     }

tutorials/mesh/snappyHexMesh/motorBike_leakDetection/system/snappyHexMeshDict

@@ -43,7 +43,7 @@ geometry
 
     refinementBox
     {
-        type searchableBox;
+        type box;
         min (-1.0 -0.7 0.0);
         max ( 8.0  0.7 2.5);
     }

tutorials/mesh/snappyHexMesh/opposite_walls/system/snappyHexMeshDict

@@ -40,7 +40,7 @@ geometry
     }
     all
     {
-        type searchableBox;
+        type box;
         min (-1000 -1000 -1000);
         max (1000 1000 1000);
     }