blockMesh: Added projected vertices and edges

Patch contributed by Mattijs Janssens

    - Added projected vertices
    - Added projected edges
    - Change of blockEdges API (operate on list lambdas)
    - Change of blockFaces API (pass in blockDescriptor and blockFacei)
    - Added sphere7ProjectedEdges tutorial to demonstrate vertex and edge projection

src/mesh/blockMesh/Make/files

@@ -1,5 +1,6 @@
 blockVertices/blockVertex/blockVertex.C
 blockVertices/pointVertex/pointVertex.C
+blockVertices/projectVertex/projectVertex.C
 
 blockEdges/blockEdge/blockEdge.C
 blockEdges/lineDivide/lineDivide.C
@@ -11,6 +12,7 @@ blockEdges/BSplineEdge/BSpline.C
 blockEdges/BSplineEdge/BSplineEdge.C
 blockEdges/splineEdge/CatmullRomSpline.C
 blockEdges/splineEdge/splineEdge.C
+blockEdges/projectEdge/projectEdge.C
 
 blockFaces/blockFace/blockFace.C
 blockFaces/projectFace/projectFace.C

src/mesh/blockMesh/blockDescriptor/blockDescriptor.C

@@ -336,7 +336,12 @@ void Foam::blockDescriptor::correctFacePoints
     {
         if (curvedFaces_[blockFacei] != -1)
         {
-            faces_[curvedFaces_[blockFacei]].project(facePoints[blockFacei]);
+            faces_[curvedFaces_[blockFacei]].project
+            (
+                *this,
+                blockFacei,
+                facePoints[blockFacei]
+            );
         }
     }
 }

src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C

@@ -125,6 +125,20 @@ Foam::pointField Foam::blockEdge::appendEndPoints
 }
 
 
+Foam::tmp<Foam::pointField>
+Foam::blockEdge::position(const scalarList& lambdas) const
+{
+    tmp<pointField> tpoints(new pointField(lambdas.size()));
+    pointField& points = tpoints.ref();
+
+    forAll(lambdas, i)
+    {
+        points[i] = position(lambdas[i]);
+    }
+    return tpoints;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const blockEdge& p)

src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H

@@ -188,6 +188,10 @@ public:
         //  0 <= lambda <= 1
         virtual point position(const scalar) const = 0;
 
+        //- Return the point positions corresponding to the curve parameters
+        //  0 <= lambda <= 1
+        virtual tmp<pointField> position(const scalarList&) const;
+
         //- Return the length of the curve
         virtual scalar length() const = 0;
 

src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C

@@ -131,10 +131,7 @@ Foam::lineDivide::lineDivide
     }
 
     // Calculate the points
-    for (label i = 0; i <= nDiv; i++)
+    points_ = cedge.position(divisions_);
-    {
-        points_[i] = cedge.position(divisions_[i]);
-    }
 }
 
 

src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     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 "searchableSurfacesQueries.H"
+#include "projectEdge.H"
+#include "unitConversion.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointConstraint.H"
+#include "plane.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(projectEdge, 0);
+    addToRunTimeSelectionTable(blockEdge, projectEdge, Istream);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::projectEdge::findNearest
+(
+    const point& pt,
+    point& near,
+    pointConstraint& constraint
+) const
+{
+    if (surfaces_.size())
+    {
+        const scalar distSqr = magSqr(points_[end_]-points_[start_]);
+
+        pointField boundaryNear(1);
+        List<pointConstraint> boundaryConstraint(1);
+
+        searchableSurfacesQueries::findNearest
+        (
+            geometry_,
+            surfaces_,
+            pointField(1, pt),
+            scalarField(1, distSqr),
+            boundaryNear,
+            boundaryConstraint
+        );
+        near = boundaryNear[0];
+        constraint = boundaryConstraint[0];
+    }
+    else
+    {
+        near = pt;
+        constraint = pointConstraint();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::projectEdge::projectEdge
+(
+    const searchableSurfaces& geometry,
+    const pointField& points,
+    Istream& is
+)
+:
+    blockEdge(points, is),
+    geometry_(geometry)
+{
+    wordList names(is);
+    surfaces_.setSize(names.size());
+    forAll(names, i)
+    {
+        surfaces_[i] = geometry_.findSurfaceID(names[i]);
+
+        if (surfaces_[i] == -1)
+        {
+            FatalIOErrorInFunction(is)
+                << "Cannot find surface " << names[i] << " in geometry"
+                << exit(FatalIOError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::point Foam::projectEdge::position(const scalar lambda) const
+{
+    // Initial guess
+    const point start(points_[start_]+lambda*(points_[end_]-points_[start_]));
+
+    point near(start);
+
+    if (lambda >= SMALL && lambda < 1.0-SMALL)
+    {
+        pointConstraint constraint;
+        findNearest(start, near, constraint);
+    }
+
+    return near;
+}
+
+
+Foam::tmp<Foam::pointField>
+Foam::projectEdge::position(const scalarList& lambdas) const
+{
+    tmp<pointField> tpoints(new pointField(lambdas.size()));
+    pointField& points = tpoints.ref();
+
+    const point& startPt = points_[start_];
+    const point& endPt = points_[end_];
+    const vector d = endPt-startPt;
+
+    forAll(lambdas, i)
+    {
+        points[i] = startPt+lambdas[i]*d;
+    }
+
+    forAll(lambdas, i)
+    {
+        if (lambdas[i] >= SMALL && lambdas[i] < 1.0-SMALL)
+        {
+            pointConstraint constraint;
+            findNearest(points[i], points[i], constraint);
+        }
+    }
+
+    return tpoints;
+}
+
+
+// ************************************************************************* //

src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     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::projectEdge
+
+Description
+    Defines the edge from the projection onto a surface (single surface)
+    or intersection of two surfaces.
+
+SourceFiles
+    projectEdge.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef projectEdge_H
+#define projectEdge_H
+
+#include "blockEdge.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class pointConstraint;
+
+/*---------------------------------------------------------------------------*\
+                         Class projectEdge Declaration
+\*---------------------------------------------------------------------------*/
+
+class projectEdge
+:
+    public blockEdge
+{
+    // Private data
+
+        const searchableSurfaces& geometry_;
+
+        //- The indices of surfaces onto which the points are projected
+        labelList surfaces_;
+
+
+    // Private Member Functions
+
+        //- Single point find nearest
+        void findNearest(const point&, point& near, pointConstraint&) const;
+
+        //- Disallow default bitwise copy construct
+        projectEdge(const projectEdge&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const projectEdge&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("project");
+
+
+    // Constructors
+
+        //- Construct from Istream setting pointsList
+        projectEdge
+        (
+            const searchableSurfaces& geometry,
+            const pointField& points,
+            Istream&
+        );
+
+
+    //- Destructor
+    virtual ~projectEdge()
+    {}
+
+
+    // Member Functions
+
+        //- Return the point positions corresponding to the curve parameters
+        //  0 <= lambda <= 1
+        virtual point position(const scalar) const;
+
+        //- Return the point positions corresponding to the curve parameters
+        //  0 <= lambda <= 1
+        virtual tmp<pointField> position(const scalarList&) const;
+
+        //- Return the length of the curve
+        virtual scalar length() const
+        {
+            NotImplemented;
+            return 1;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/mesh/blockMesh/blockFaces/blockFace/blockFace.H

@@ -44,11 +44,11 @@ namespace Foam
 
 // Forward declaration of friend functions and operators
 
+class blockDescriptor;
 class blockFace;
 
 Ostream& operator<<(Ostream&, const blockFace&);
 
-
 /*---------------------------------------------------------------------------*\
                          Class blockFace Declaration
 \*---------------------------------------------------------------------------*/
@@ -138,7 +138,12 @@ public:
         //- Compare with the given block and block face
         inline bool compare(const face& vertices) const;
 
-        virtual void project(pointField& points) const = 0;
+        virtual void project
+        (
+            const blockDescriptor&,
+            const label blockFacei,
+            pointField& points
+        ) const = 0;
 
 
     // Ostream operator

src/mesh/blockMesh/blockFaces/projectFace/projectFace.C

@@ -26,6 +26,7 @@ License
 #include "projectFace.H"
 #include "unitConversion.H"
 #include "addToRunTimeSelectionTable.H"
+#include "blockDescriptor.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -80,7 +81,12 @@ Foam::blockFaces::projectFace::projectFace
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::blockFaces::projectFace::project(pointField& points) const
+void Foam::blockFaces::projectFace::project
+(
+    const blockDescriptor& desc,
+    const label blockFacei,
+    pointField& points
+) const
 {
     List<pointIndexHit> hits;
     scalarField nearestDistSqr

src/mesh/blockMesh/blockFaces/projectFace/projectFace.H

@@ -98,7 +98,12 @@ public:
     // Member Functions
 
         //- Project the given points onto the surface
-        virtual void project(pointField& points) const;
+        virtual void project
+        (
+            const blockDescriptor&,
+            const label blockFacei,
+            pointField& points
+        ) const;
 };
 
 

src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H

@@ -51,7 +51,9 @@ class pointVertex
 :
     public blockVertex
 {
-    // Private member data
+protected:
+
+    // Protected member data
 
         //- The vertex location
         point vertex_;

src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C

@@ -23,77 +23,76 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "projectFace.H"
+#include "projectVertex.H"
 #include "unitConversion.H"
 #include "addToRunTimeSelectionTable.H"
+#include "searchableSurfacesQueries.H"
+#include "pointConstraint.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(projectFace, 0);
+namespace blockVertices
-    addToRunTimeSelectionTable(blockVertex, projectFace, Istream);
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-const Foam::searchableSurface& Foam::projectFace::lookupSurface
-(
-    const searchableSurfaces& geometry,
-    Istream& is
-) const
 {
-    word name(is);
+    defineTypeNameAndDebug(projectVertex, 0);
-
+    addToRunTimeSelectionTable(blockVertex, projectVertex, Istream);
-    forAll(geometry, i)
-    {
-        if (geometry[i].name() == name)
-        {
-            return geometry[i];
 }
 }
 
-    FatalIOErrorInFunction(is)
-        << "Cannot find surface " << name << " in geometry"
-        << exit(FatalIOError);
-
-    return geometry[0];
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::projectFace::projectFace
+Foam::blockVertices::projectVertex::projectVertex
 (
     const searchableSurfaces& geometry,
     Istream& is
 )
 :
-    blockVertex(is),
+    pointVertex(geometry, is),
-    surface_(lookupSurface(geometry, is))
+    geometry_(geometry)
-{}
+{
+    wordList names(is);
+    surfaces_.setSize(names.size());
+    forAll(names, i)
+    {
+        surfaces_[i] = geometry_.findSurfaceID(names[i]);
+
+        if (surfaces_[i] == -1)
+        {
+            FatalIOErrorInFunction(is)
+                << "Cannot find surface " << names[i] << " in geometry"
+                << exit(FatalIOError);
+        }
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::projectFace::project(pointField& points) const
+Foam::blockVertices::projectVertex::operator point() const
 {
-    List<pointIndexHit> hits;
+    pointField start(1, pointVertex::operator point());
-    scalarField nearestDistSqr
-    (
-        points.size(),
-        magSqr(points[0] - points[points.size()-1])
-    );
-    surface_.findNearest(points, nearestDistSqr, hits);
 
-    forAll(hits, i)
+    pointField boundaryNear(start);
-    {
+    List<pointConstraint> boundaryConstraint;
-        if (hits[i].hit())
+
-        {
+
-            points[i] = hits[i].hitPoint();
+    // Note: how far do we need to search? Probably not further than
-        }
+    //       span of surfaces themselves.
-    }
+    boundBox bb(searchableSurfacesQueries::bounds(geometry_, surfaces_));
+
+    searchableSurfacesQueries::findNearest
+    (
+        geometry_,
+        surfaces_,
+        start,
+        scalarField(start.size(), magSqr(bb.span())),
+        boundaryNear,
+        boundaryConstraint
+    );
+
+    return boundaryNear[0];
 }
 
 

src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H

@@ -22,54 +22,52 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::projectFace
+    Foam::blockVertices::projectVertex
 
 Description
-    Projects the given set of face points onto the selected surface of the
+    Projects the vertex onto the selected surfacees of the
     geometry provided as a searchableSurfaces object.
 
 SourceFiles
-    projectFace.C
+    projectVertex.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef projectFace_H
+#ifndef blockVertices_projectVertex_H
-#define projectFace_H
+#define blockVertices_projectVertex_H
 
-#include "blockVertex.H"
+#include "pointVertex.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+namespace blockVertices
+{
 
 /*---------------------------------------------------------------------------*\
-                           Class projectFace Declaration
+                        Class projectVertex Declaration
 \*---------------------------------------------------------------------------*/
 
-class projectFace
+class projectVertex
 :
-    public blockVertex
+    public pointVertex
 {
     // Private data
 
-        //- The surface onto which the points are projected
+        const searchableSurfaces& geometry_;
-        const searchableSurface& surface_;
+
+        //- The indices of surfaces onto which the points are projected
+        labelList surfaces_;
 
 
     // Private Member Functions
 
-        const searchableSurface& lookupSurface
-        (
-            const searchableSurfaces& geometry,
-            Istream& is
-        ) const;
-
         //- Disallow default bitwise copy construct
-        projectFace(const projectFace&);
+        projectVertex(const projectVertex&);
 
         //- Disallow default bitwise assignment
-        void operator=(const projectFace&);
+        void operator=(const projectVertex&);
 
 
 public:
@@ -81,7 +79,7 @@ public:
     // Constructors
 
         //- Construct from Istream setting pointsList
-        projectFace
+        projectVertex
         (
             const searchableSurfaces& geometry,
             Istream&
@@ -89,19 +87,20 @@ public:
 
 
     //- Destructor
-    virtual ~projectFace()
+    virtual ~projectVertex()
     {}
 
 
     // Member Functions
 
         //- Project the given points onto the surface
-        virtual void project(pointField& points) const;
+        virtual operator point() const;
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace blockVertices
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/meshTools/searchableSurface/searchableSurfaces.C

@@ -419,24 +419,6 @@ Foam::boundBox Foam::searchableSurfaces::bounds() const
 }
 
 
-Foam::pointIndexHit Foam::searchableSurfaces::facesIntersection
-(
-    const scalar initDistSqr,
-    const scalar convergenceDistSqr,
-    const point& start
-) const
-{
-    return searchableSurfacesQueries::facesIntersection
-    (
-        *this,
-        allSurfaces_,
-        initDistSqr,
-        convergenceDistSqr,
-        start
-    );
-}
-
-
 bool Foam::searchableSurfaces::checkClosed(const bool report) const
 {
     if (report)

src/meshTools/searchableSurface/searchableSurfaces.H

@@ -213,18 +213,6 @@ public:
             //- Calculate bounding box
             boundBox bounds() const;
 
-
-        // Single point queries
-
-            //- Calculate point which is on a set of surfaces.
-            pointIndexHit facesIntersection
-            (
-                const scalar initialDistSqr,
-                const scalar convergenceDistSqr,
-                const point& start
-            ) const;
-
-
         // Checking
 
             //- Are all surfaces closed and manifold

src/meshTools/searchableSurface/searchableSurfacesQueries.C

@@ -28,6 +28,8 @@ License
 #include "OFstream.H"
 #include "meshTools.H"
 #include "DynamicField.H"
+#include "pointConstraint.H"
+#include "plane.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -39,288 +41,6 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::pointIndexHit Foam::searchableSurfacesQueries::tempFindNearest
-(
-    const searchableSurface& surf,
-    const point& pt,
-    const scalar initDistSqr
-)
-{
-    pointField onePoint(1, pt);
-    scalarField oneDist(1, initDistSqr);
-    List<pointIndexHit> oneHit(1);
-    surf.findNearest(onePoint, oneDist, oneHit);
-    return oneHit[0];
-}
-
-
-// Calculate sum of distance to surfaces.
-Foam::scalar Foam::searchableSurfacesQueries::sumDistSqr
-(
-    const PtrList<searchableSurface>& allSurfaces,
-    const labelList& surfacesToTest,
-    const scalar initDistSqr,
-    const point& pt
-)
-{
-    scalar sum = 0;
-
-    forAll(surfacesToTest, testI)
-    {
-        label surfI = surfacesToTest[testI];
-
-        pointIndexHit hit
-        (
-            tempFindNearest(allSurfaces[surfI], pt, initDistSqr)
-        );
-
-        // Note: make it fall over if not hit.
-        sum += magSqr(hit.hitPoint()-pt);
-    }
-    return sum;
-}
-
-
-// Reflects the point furthest away around the triangle centre by a factor fac.
-// (triangle centre is the average of all points but the ihi. pSum is running
-//  sum of all points)
-Foam::scalar Foam::searchableSurfacesQueries::tryMorphTet
-(
-    const PtrList<searchableSurface>& allSurfaces,
-    const labelList& surfacesToTest,
-    const scalar initDistSqr,
-    List<vector>& p,
-    List<scalar>& y,
-    vector& pSum,
-    const label ihi,
-    const scalar fac
-)
-{
-    scalar fac1 = (1.0-fac)/vector::nComponents;
-    scalar fac2 = fac1-fac;
-
-    vector ptry = pSum*fac1-p[ihi]*fac2;
-
-    scalar ytry = sumDistSqr(allSurfaces, surfacesToTest, initDistSqr, ptry);
-
-    if (ytry < y[ihi])
-    {
-        y[ihi] = ytry;
-        pSum += ptry - p[ihi];
-        p[ihi] = ptry;
-    }
-    return ytry;
-}
-
-
-bool Foam::searchableSurfacesQueries::morphTet
-(
-    const PtrList<searchableSurface>& allSurfaces,
-    const labelList& surfacesToTest,
-    const scalar initDistSqr,
-    const scalar convergenceDistSqr,
-    const label maxIter,
-    List<vector>& p,
-    List<scalar>& y
-)
-{
-    vector pSum = sum(p);
-
-    autoPtr<OFstream> str;
-    label vertI = 0;
-    if (debug)
-    {
-        wordList names(surfacesToTest.size());
-        forAll(surfacesToTest, i)
-        {
-            names[i] = allSurfaces[surfacesToTest[i]].name();
-        }
-        Pout<< "searchableSurfacesQueries::morphTet : intersection of "
-            << names << " starting from points:" << p << endl;
-        str.reset(new OFstream("track.obj"));
-        meshTools::writeOBJ(str(), p[0]);
-        vertI++;
-    }
-
-    for (label iter = 0; iter < maxIter; iter++)
-    {
-        // Get the indices of lowest, highest and second-highest values.
-        label ilo, ihi, inhi;
-        {
-            labelList sortedIndices;
-            sortedOrder(y, sortedIndices);
-            ilo  = sortedIndices[0];
-            ihi  = sortedIndices[sortedIndices.size()-1];
-            inhi = sortedIndices[sortedIndices.size()-2];
-        }
-
-        if (debug)
-        {
-            Pout<< "Iteration:" << iter
-                << " lowest:" << y[ilo] << " highest:" << y[ihi]
-                << " points:" << p << endl;
-
-            meshTools::writeOBJ(str(), p[ilo]);
-            vertI++;
-            str()<< "l " << vertI-1 << ' ' << vertI << nl;
-        }
-
-        if (y[ihi] < convergenceDistSqr)
-        {
-            // Get point on 0th surface.
-            Swap(p[0], p[ilo]);
-            Swap(y[0], y[ilo]);
-            return true;
-        }
-
-        // Reflection: point furthest away gets reflected.
-        scalar ytry = tryMorphTet
-        (
-            allSurfaces,
-            surfacesToTest,
-            10*y[ihi],             // search box.
-            p,
-            y,
-            pSum,
-            ihi,
-            -1.0
-        );
-
-        if (ytry <= y[ilo])
-        {
-            // If in right direction (y lower) expand by two.
-            ytry = tryMorphTet
-            (
-                allSurfaces,
-                surfacesToTest,
-                10*y[ihi],
-                p,
-                y,
-                pSum,
-                ihi,
-                2.0
-            );
-        }
-        else if (ytry >= y[inhi])
-        {
-            // If inside tet try contraction.
-
-            scalar ysave = y[ihi];
-
-            ytry = tryMorphTet
-            (
-                allSurfaces,
-                surfacesToTest,
-                10*y[ihi],
-                p,
-                y,
-                pSum,
-                ihi,
-                0.5
-            );
-
-            if (ytry >= ysave)
-            {
-                // Contract around lowest point.
-                forAll(p, i)
-                {
-                    if (i != ilo)
-                    {
-                        p[i] = 0.5*(p[i] + p[ilo]);
-                        y[i] = sumDistSqr
-                        (
-                            allSurfaces,
-                            surfacesToTest,
-                            y[ihi],
-                            p[i]
-                        );
-                    }
-                }
-                pSum = sum(p);
-            }
-        }
-    }
-
-    if (debug)
-    {
-        meshTools::writeOBJ(str(), p[0]);
-        vertI++;
-        str()<< "l " << vertI-1 << ' ' << vertI << nl;
-    }
-
-    // Failure to converge. Return best guess so far.
-    label ilo = findMin(y);
-    Swap(p[0], p[ilo]);
-    Swap(y[0], y[ilo]);
-    return false;
-}
-
-
-//// Get all intersections (in order) for single surface.
-//void Foam::searchableSurfacesQueries::findAllIntersections
-//(
-//    const searchableSurface& s,
-//    const pointField& start,
-//    const pointField& end,
-//    const vectorField& smallVec,
-//    List<List<pointIndexHit>>& surfaceHitInfo
-//)
-//{
-//    surfaceHitInfo.setSize(start.size());
-//
-//    // Current start point of vector
-//    pointField p0(start);
-//
-//    List<pointIndexHit> intersectInfo(start.size());
-//
-//    // For test whether finished doing vector.
-//    const vectorField dirVec(end-start);
-//    const scalarField magSqrDirVec(magSqr(dirVec));
-//
-//    while (true)
-//    {
-//        // Find first intersection. Synced.
-//        s.findLine(p0, end, intersectInfo);
-//
-//        label nHits = 0;
-//
-//        forAll(intersectInfo, i)
-//        {
-//            if (intersectInfo[i].hit())
-//            {
-//                nHits++;
-//
-//                label sz = surfaceHitInfo[i].size();
-//                surfaceHitInfo[i].setSize(sz+1);
-//                surfaceHitInfo[i][sz] = intersectInfo[i];
-//
-//                p0[i] = intersectInfo[i].hitPoint() + smallVec[i];
-//
-//                // If beyond endpoint set to endpoint so as not to pick up
-//                // any intersections. Could instead just filter out hits.
-//                if (((p0[i]-start[i])&dirVec[i]) > magSqrDirVec[i])
-//                {
-//                    p0[i] = end[i];
-//                }
-//            }
-//            else
-//            {
-//                // Set to endpoint to stop intersection test. See above.
-//                p0[i] = end[i];
-//            }
-//        }
-//
-//        // returnReduce(nHits) ?
-//        if (nHits == 0)
-//        {
-//            break;
-//        }
-//    }
-//}
-
-
-// Given current set of hits (allSurfaces, allInfo) merge in those coming from
-// surface surfI.
 void Foam::searchableSurfacesQueries::mergeHits
 (
     const point& start,
@@ -333,6 +53,9 @@ void Foam::searchableSurfacesQueries::mergeHits
     scalarList& allDistSqr
 )
 {
+    // Given current set of hits (allSurfaces, allInfo) merge in those coming
+    // from surface surfI.
+
     // Precalculate distances
     scalarList surfDistSqr(surfHits.size());
     forAll(surfHits, i)
@@ -532,7 +255,6 @@ void Foam::searchableSurfacesQueries::findAllIntersections
 }
 
 
-//Find intersections of edge nearest to both endpoints.
 void Foam::searchableSurfacesQueries::findNearestIntersection
 (
    const PtrList<searchableSurface>& allSurfaces,
@@ -620,7 +342,6 @@ void Foam::searchableSurfacesQueries::findNearestIntersection
 }
 
 
-// Find nearest. Return -1 or nearest point
 void Foam::searchableSurfacesQueries::findNearest
 (
     const PtrList<searchableSurface>& allSurfaces,
@@ -631,6 +352,8 @@ void Foam::searchableSurfacesQueries::findNearest
     List<pointIndexHit>& nearestInfo
 )
 {
+    // Find nearest. Return -1 or nearest point
+
     // Initialise
     nearestSurfaces.setSize(samples.size());
     nearestSurfaces = -1;
@@ -667,7 +390,6 @@ void Foam::searchableSurfacesQueries::findNearest
 }
 
 
-// Find nearest. Return -1 or nearest point
 void Foam::searchableSurfacesQueries::findNearest
 (
     const PtrList<searchableSurface>& allSurfaces,
@@ -679,6 +401,8 @@ void Foam::searchableSurfacesQueries::findNearest
     List<pointIndexHit>& nearestInfo
 )
 {
+    // Find nearest. Return -1 or nearest point
+
     if (regionIndices.empty())
     {
         findNearest
@@ -729,6 +453,105 @@ void Foam::searchableSurfacesQueries::findNearest
 }
 
 
+void Foam::searchableSurfacesQueries::findNearest
+(
+    const PtrList<searchableSurface>& allSurfaces,
+    const labelList& surfacesToTest,
+    const pointField& start,
+    const scalarField& distSqr,
+    pointField& near,
+    List<pointConstraint>& constraint,
+    const label nIter
+)
+{
+    // Multi-surface findNearest
+
+    vectorField normal;
+    List<pointIndexHit> info;
+
+    allSurfaces[surfacesToTest[0]].findNearest(start, distSqr, info);
+    allSurfaces[surfacesToTest[0]].getNormal(info, normal);
+
+    // Extract useful info from initial start point
+    near = start;
+    forAll(info, i)
+    {
+        if (info[i].hit())
+        {
+            near[i] = info[i].hitPoint();
+        }
+    }
+    constraint.setSize(near.size());
+
+
+    if (surfacesToTest.size() == 1)
+    {
+        constraint = pointConstraint();
+        forAll(info, i)
+        {
+            if (info[i].hit())
+            {
+                constraint[i].applyConstraint(normal[i]);
+            }
+        }
+    }
+    else if (surfacesToTest.size() >= 2)
+    {
+        // Work space
+        pointField near1;
+        vectorField normal1;
+
+        label surfi = 1;
+        for (label iter = 0; iter < nIter; iter++)
+        {
+            constraint = pointConstraint();
+            forAll(constraint, i)
+            {
+                if (info[i].hit())
+                {
+                    constraint[i].applyConstraint(normal[i]);
+                }
+            }
+
+            // Find intersection with next surface
+            const searchableSurface& s = allSurfaces[surfacesToTest[surfi]];
+            s.findNearest(near, distSqr, info);
+            s.getNormal(info, normal1);
+            near1.setSize(info.size());
+            forAll(info, i)
+            {
+                if (info[i].hit())
+                {
+                    near1[i] = info[i].hitPoint();
+                }
+            }
+
+            // Move to intersection
+            forAll(near, pointi)
+            {
+                if (info[pointi].hit())
+                {
+                    plane pl0(near[pointi], normal[pointi]);
+                    plane pl1(near1[pointi], normal1[pointi]);
+                    plane::ray r(pl0.planeIntersect(pl1));
+                    vector n = r.dir() / mag(r.dir());
+
+                    vector d(r.refPoint()-near[pointi]);
+                    d -= (d&n)*n;
+
+                    near[pointi] += d;
+                    normal[pointi] = normal1[pointi];
+                    constraint[pointi].applyConstraint(normal1[pointi]);
+                }
+            }
+
+            // Step to next surface
+            surfi = surfacesToTest.fcIndex(surfi);
+        }
+    }
+}
+
+
 void Foam::searchableSurfacesQueries::signedDistance
 (
     const PtrList<searchableSurface>& allSurfaces,
@@ -849,126 +672,4 @@ Foam::boundBox Foam::searchableSurfacesQueries::bounds
 }
 
 
-Foam::pointIndexHit Foam::searchableSurfacesQueries::facesIntersection
-(
-    const PtrList<searchableSurface>& allSurfaces,
-    const labelList& surfacesToTest,
-    const scalar initDistSqr,
-    const scalar convergenceDistSqr,
-    const point& start
-)
-{
-    // Get four starting points. Take these as the projection of the
-    // starting point onto the surfaces and the mid point
-    List<point> nearest(surfacesToTest.size()+1);
-
-    point sumNearest = Zero;
-
-    forAll(surfacesToTest, i)
-    {
-        pointIndexHit hit
-        (
-            tempFindNearest(allSurfaces[surfacesToTest[i]], start, initDistSqr)
-        );
-
-        if (hit.hit())
-        {
-            nearest[i] = hit.hitPoint();
-            sumNearest += nearest[i];
-        }
-        else
-        {
-            FatalErrorInFunction
-                << "Did not find point within distance "
-                << initDistSqr << " of starting point " << start
-                << " on surface "
-                << allSurfaces[surfacesToTest[i]].IOobject::name()
-                << abort(FatalError);
-        }
-    }
-
-    nearest.last() = sumNearest / surfacesToTest.size();
-
-
-    // Get the sum of distances (initial evaluation)
-    List<scalar> nearestDist(nearest.size());
-
-    forAll(nearestDist, i)
-    {
-        nearestDist[i] = sumDistSqr
-        (
-            allSurfaces,
-            surfacesToTest,
-            initDistSqr,
-            nearest[i]
-        );
-    }
-
-
-    //- Downhill Simplex method
-
-    bool converged = morphTet
-    (
-        allSurfaces,
-        surfacesToTest,
-        initDistSqr,
-        convergenceDistSqr,
-        2000,
-        nearest,
-        nearestDist
-    );
-
-
-    pointIndexHit intersection;
-
-    if (converged)
-    {
-        // Project nearest onto 0th surface.
-        intersection = tempFindNearest
-        (
-            allSurfaces[surfacesToTest[0]],
-            nearest[0],
-            nearestDist[0]
-        );
-    }
-
-    //if (!intersection.hit())
-    //{
-    //    // Restart
-    //    scalar smallDist = Foam::sqr(convergenceDistSqr);
-    //    nearest[0] = intersection.hitPoint();
-    //    nearest[1] = nearest[0];
-    //    nearest[1].x() += smallDist;
-    //    nearest[2] = nearest[0];
-    //    nearest[2].y() += smallDist;
-    //    nearest[3] = nearest[0];
-    //    nearest[3].z() += smallDist;
-    //
-    //    forAll(nearestDist, i)
-    //    {
-    //        nearestDist[i] = sumDistSqr
-    //        (
-    //            surfacesToTest,
-    //            initDistSqr,
-    //            nearest[i]
-    //        );
-    //    }
-    //
-    //    intersection = morphTet
-    //    (
-    //        allSurfaces,
-    //        surfacesToTest,
-    //        initDistSqr,
-    //        convergenceDistSqr,
-    //        1000,
-    //        nearest,
-    //        nearestDist
-    //    );
-    //}
-
-    return intersection;
-}
-
-
-
 // ************************************************************************* //

src/meshTools/searchableSurface/searchableSurfacesQueries.H

@@ -44,6 +44,7 @@ namespace Foam
 
 // Forward declaration of classes
 class plane;
+class pointConstraint;
 
 /*---------------------------------------------------------------------------*\
                            Class searchableSurfacesQueries Declaration
@@ -51,68 +52,8 @@ class plane;
 
 class searchableSurfacesQueries
 {
-    // Private data
-
     // Private Member Functions
 
-        //- Temporary wrapper around findNearest. Used in facesIntersection only
-        static pointIndexHit tempFindNearest
-        (
-            const searchableSurface&,
-            const point& pt,
-            const scalar initDistSqr
-        );
-
-        //- Calculate sum of distances to nearest point on surfaces. Is used
-        //  in minimisation to find intersection. Returns sum of (square of)
-        //  distances to the surfaces.
-        static scalar sumDistSqr
-        (
-            const PtrList<searchableSurface>&,
-            const labelList& surfacesToTest,
-            const scalar initialDistSqr,    // search box
-            const point& pt
-        );
-
-        //- Takes the tet (points p) and reflects the point with the
-        //  highest value around the centre (pSum). Checks if it gets closer
-        //  and updates p, y if so.
-        static scalar tryMorphTet
-        (
-            const PtrList<searchableSurface>&,
-            const labelList& surfacesToTest,
-            const scalar initialDistSqr,
-            List<vector>& p,
-            List<scalar>& y,
-            vector& pSum,
-            const label ihi,
-            const scalar fac
-        );
-
-        //- Downhill simplex method: find the point with min cumulative
-        //  distance to all surfaces. Does so by morphing a tet (points p).
-        //  Returns the point on the 0th surface or hit if not reached within
-        //  maxIters iterations.
-        static bool morphTet
-        (
-            const PtrList<searchableSurface>&,
-            const labelList& surfacesToTest,
-            const scalar initialDistSqr,
-            const scalar convergenceDistSqr,
-            const label maxIter,
-            List<vector>& p,
-            List<scalar>& y
-        );
-
-        //static void findAllIntersections
-        //(
-        //    const searchableSurface& s,
-        //    const pointField& start,
-        //    const pointField& end,
-        //    const vectorField& smallVec,
-        //    List<List<pointIndexHit>>&
-        //);
-
         static void mergeHits
         (
             const point& start,
@@ -158,7 +99,7 @@ public:
                 List<List<pointIndexHit>>& surfaceHits
             );
 
-            //Find intersections of edge nearest to both endpoints.
+            //- Find intersections of edge nearest to both endpoints.
             static void findNearestIntersection
             (
                 const PtrList<searchableSurface>& allSurfaces,
@@ -195,6 +136,21 @@ public:
                 List<pointIndexHit>& nearestInfo
             );
 
+            //- Find nearest points that are on all supplied surfaces
+            //  (nearest point if single surface; nearest intersection by
+            //   steepst descent if on multiple surfaces). Returns current
+            //   best guess). Wip.
+            static void findNearest
+            (
+                const PtrList<searchableSurface>& allSurfaces,
+                const labelList& surfacesToTest,
+                const pointField& start,
+                const scalarField& distSqr,
+                pointField& near,
+                List<pointConstraint>& constraint,
+                const label nIter = 20
+            );
+
             //- Find signed distance to nearest surface. Outside is positive.
             //  illegalHandling: how to handle non-inside or outside
             //      OUTSIDE : treat as outside
@@ -217,18 +173,6 @@ public:
                 const PtrList<searchableSurface>& allSurfaces,
                 const labelList& surfacesToTest
             );
-
-        // Single point queries
-
-            //- Calculate point which is on a set of surfaces. WIP.
-            static pointIndexHit facesIntersection
-            (
-                const PtrList<searchableSurface>& allSurfaces,
-                const labelList& surfacesToTest,
-                const scalar initDistSqr,
-                const scalar convergenceDistSqr,
-                const point& start
-            );
 };
 
 

tutorials/mesh/blockMesh/sphere7ProjectedEdges/Allrun

@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+
+#------------------------------------------------------------------------------

tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/blockMeshDict

@@ -0,0 +1,143 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+verbose no;
+
+geometry
+{
+    sphere
+    {
+        type searchableSphere;
+        centre (0 0 0);
+        radius 1;
+    }
+
+    innerSphere
+    {
+        type searchableSphere;
+        centre (0 0 0);
+        radius 0.5;
+    }
+}
+
+scale 1;
+
+n    10;
+
+// Rough approximation of corner points
+v    0.5773502;
+mv  -0.5773502;
+vh   0.2886751;
+mvh -0.2886751;
+
+vertices
+(
+    // Inner block points
+    project ($mvh $mvh $mvh) (innerSphere)
+    project ( $vh $mvh $mvh) (innerSphere)
+    project ( $vh  $vh $mvh) (innerSphere)
+    project ($mvh  $vh $mvh) (innerSphere)
+    project ($mvh $mvh  $vh) (innerSphere)
+    project ( $vh $mvh  $vh) (innerSphere)
+    project ( $vh  $vh  $vh) (innerSphere)
+    project ($mvh  $vh  $vh) (innerSphere)
+
+    // Outer block points
+    project ($mv $mv $mv) (sphere)
+    project ( $v $mv $mv) (sphere)
+    project ( $v  $v $mv) (sphere)
+    project ($mv  $v $mv) (sphere)
+    project ($mv $mv  $v) (sphere)
+    project ( $v $mv  $v) (sphere)
+    project ( $v  $v  $v) (sphere)
+    project ($mv  $v  $v) (sphere)
+);
+
+blocks
+(
+    hex ( 0  1  2  3  4  5  6  7) ($n $n $n) simpleGrading (1 1 1)
+    hex ( 9  8 12 13  1  0  4  5) ($n $n $n) simpleGrading (1 1 1)
+    hex (10  9 13 14  2  1  5  6) ($n $n $n) simpleGrading (1 1 1)
+    hex (11 10 14 15  3  2  6  7) ($n $n $n) simpleGrading (1 1 1)
+    hex ( 8 11 15 12  0  3  7  4) ($n $n $n) simpleGrading (1 1 1)
+    hex ( 8  9 10 11  0  1  2  3) ($n $n $n) simpleGrading (1 1 1)
+    hex (13 12 15 14  5  4  7  6) ($n $n $n) simpleGrading (1 1 1)
+);
+
+edges
+(
+    // Outer edges
+    project  8  9 (sphere)
+    project 10 11 (sphere)
+    project 14 15 (sphere)
+    project 12 13 (sphere)
+
+    project  8 11 (sphere)
+    project  9 10 (sphere)
+    project 13 14 (sphere)
+    project 12 15 (sphere)
+
+    project  8 12 (sphere)
+    project  9 13 (sphere)
+    project 10 14 (sphere)
+    project 11 15 (sphere)
+
+
+    // Inner edges
+    project  0  1 (innerSphere)
+    project  2  3 (innerSphere)
+    project  6  7 (innerSphere)
+    project  4  5 (innerSphere)
+
+    project  0  3 (innerSphere)
+    project  1  2 (innerSphere)
+    project  5  6 (innerSphere)
+    project  4  7 (innerSphere)
+
+    project  0  4 (innerSphere)
+    project  1  5 (innerSphere)
+    project  2  6 (innerSphere)
+    project  3  7 (innerSphere)
+);
+
+faces
+(
+    project ( 8 12 15 11) sphere
+    project (10 14 13  9) sphere
+    project ( 9 13 12  8) sphere
+    project (11 15 14 10) sphere
+    project ( 8 11 10  9) sphere
+    project (12 13 14 15) sphere
+);
+
+boundary
+(
+    walls
+    {
+        type wall;
+        faces
+        (
+            ( 8 12 15 11)
+            (10 14 13  9)
+            ( 9 13 12  8)
+            (11 15 14 10)
+            ( 8 11 10  9)
+            (12 13 14 15)
+        );
+    }
+);
+
+// ************************************************************************* //

tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/controlDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     blockMesh;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0;
+
+deltaT          0;
+
+writeControl    timeStep;
+
+writeInterval   1;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //

tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/fvSchemes

@@ -0,0 +1,37 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{}
+
+gradSchemes
+{}
+
+divSchemes
+{}
+
+laplacianSchemes
+{}
+
+interpolationSchemes
+{}
+
+snGradSchemes
+{}
+
+
+// ************************************************************************* //

tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/fvSolution

@@ -0,0 +1,18 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// ************************************************************************* //