From fc6a8b648ffd2e95ef95598bee4b2338afb5aa76 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 17 Aug 2017 10:00:14 +0100
Subject: [PATCH] ENH: shortestPath: new line sampling method ('set') that
 detects the      shortest path between two sets of points. See
 shortestPathSet.H.      Fixes #572.

---
 src/sampling/Make/files                       |   1 +
 .../sampledSet/sampledSet/sampledSet.C        |  72 ++++
 .../sampledSet/sampledSet/sampledSet.H        |   6 +-
 .../sampledSet/sampledSets/sampledSets.C      |  59 +--
 .../sampledSet/shortestPath/shortestPathSet.C | 344 ++++++++++++++++++
 .../sampledSet/shortestPath/shortestPathSet.H | 149 ++++++++
 6 files changed, 572 insertions(+), 59 deletions(-)
 create mode 100644 src/sampling/sampledSet/shortestPath/shortestPathSet.C
 create mode 100644 src/sampling/sampledSet/shortestPath/shortestPathSet.H

diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index de74ebb8de..3bb94cc39b 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -16,6 +16,7 @@ sampledSet/sampledSets/sampledSetsGrouping.C
 sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
 sampledSet/uniform/uniformSet.C
 sampledSet/array/arraySet.C
+sampledSet/shortestPath/shortestPathSet.C
 
 surface/cuttingPlane/cuttingPlane.C
 surface/isoSurface/isoSurface.C
diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.C b/src/sampling/sampledSet/sampledSet/sampledSet.C
index 78694160f3..a6da49aad3 100644
--- a/src/sampling/sampledSet/sampledSet/sampledSet.C
+++ b/src/sampling/sampledSet/sampledSet/sampledSet.C
@@ -29,6 +29,7 @@ License
 #include "meshSearch.H"
 #include "writer.H"
 #include "particle.H"
+#include "SortableList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -394,6 +395,77 @@ void Foam::sampledSet::setSamples
 }
 
 
+Foam::autoPtr<Foam::coordSet> Foam::sampledSet::gather
+(
+    labelList& indexSet
+) const
+{
+    // Combine sampleSet from processors. Sort by curveDist. Return
+    // ordering in indexSet.
+    // Note: only master results are valid
+
+    // Collect data from all processors
+    List<List<point>> gatheredPts(Pstream::nProcs());
+    gatheredPts[Pstream::myProcNo()] = *this;
+    Pstream::gatherList(gatheredPts);
+
+    List<labelList> gatheredSegments(Pstream::nProcs());
+    gatheredSegments[Pstream::myProcNo()] = segments();
+    Pstream::gatherList(gatheredSegments);
+
+    List<scalarList> gatheredDist(Pstream::nProcs());
+    gatheredDist[Pstream::myProcNo()] = curveDist();
+    Pstream::gatherList(gatheredDist);
+
+
+    // Combine processor lists into one big list.
+    List<point> allPts
+    (
+        ListListOps::combine<List<point>>
+        (
+            gatheredPts, accessOp<List<point>>()
+        )
+    );
+    labelList allSegments
+    (
+        ListListOps::combine<labelList>
+        (
+            gatheredSegments, accessOp<labelList>()
+        )
+    );
+    scalarList allCurveDist
+    (
+        ListListOps::combine<scalarList>
+        (
+            gatheredDist, accessOp<scalarList>()
+        )
+    );
+
+
+    if (Pstream::master() && allCurveDist.size() == 0)
+    {
+        WarningInFunction
+            << "Sample set " << name()
+            << " has zero points." << endl;
+    }
+
+    // Sort curveDist and use to fill masterSamplePts
+    SortableList<scalar> sortedDist(allCurveDist);
+    indexSet = sortedDist.indices();
+
+    return autoPtr<coordSet>
+    (
+        new coordSet
+        (
+            name(),
+            axis(),
+            List<point>(UIndirectList<point>(allPts, indexSet)),
+            sortedDist
+        )
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::sampledSet::sampledSet
diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.H b/src/sampling/sampledSet/sampledSet/sampledSet.H
index 1f146c9a86..172f9bfa65 100644
--- a/src/sampling/sampledSet/sampledSet/sampledSet.H
+++ b/src/sampling/sampledSet/sampledSet/sampledSet.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -266,6 +266,10 @@ public:
 
         //- Output for debugging
         Ostream& write(Ostream&) const;
+
+        //- Helper: gather onto master and sort. Return (on master) gathered set
+        //  and overall sort order
+        autoPtr<coordSet> gather(labelList& indexSet) const;
 };
 
 
diff --git a/src/sampling/sampledSet/sampledSets/sampledSets.C b/src/sampling/sampledSet/sampledSets/sampledSets.C
index 2115280372..b762de557a 100644
--- a/src/sampling/sampledSet/sampledSets/sampledSets.C
+++ b/src/sampling/sampledSet/sampledSets/sampledSets.C
@@ -70,67 +70,10 @@ void Foam::sampledSets::combineSampledSets
 
     forAll(sampledSets, setI)
     {
-        const sampledSet& samplePts = sampledSets[setI];
-
-        // Collect data from all processors
-        List<List<point>> gatheredPts(Pstream::nProcs());
-        gatheredPts[Pstream::myProcNo()] = samplePts;
-        Pstream::gatherList(gatheredPts);
-
-        List<labelList> gatheredSegments(Pstream::nProcs());
-        gatheredSegments[Pstream::myProcNo()] = samplePts.segments();
-        Pstream::gatherList(gatheredSegments);
-
-        List<scalarList> gatheredDist(Pstream::nProcs());
-        gatheredDist[Pstream::myProcNo()] = samplePts.curveDist();
-        Pstream::gatherList(gatheredDist);
-
-
-        // Combine processor lists into one big list.
-        List<point> allPts
-        (
-            ListListOps::combine<List<point>>
-            (
-                gatheredPts, accessOp<List<point>>()
-            )
-        );
-        labelList allSegments
-        (
-            ListListOps::combine<labelList>
-            (
-                gatheredSegments, accessOp<labelList>()
-            )
-        );
-        scalarList allCurveDist
-        (
-            ListListOps::combine<scalarList>
-            (
-                gatheredDist, accessOp<scalarList>()
-            )
-        );
-
-
-        if (Pstream::master() && allCurveDist.size() == 0)
-        {
-            WarningInFunction
-                << "Sample set " << samplePts.name()
-                << " has zero points." << endl;
-        }
-
-        // Sort curveDist and use to fill masterSamplePts
-        SortableList<scalar> sortedDist(allCurveDist);
-        indexSets[setI] = sortedDist.indices();
-
         masterSampledSets.set
         (
             setI,
-            new coordSet
-            (
-                samplePts.name(),
-                samplePts.axis(),
-                List<point>(UIndirectList<point>(allPts, indexSets[setI])),
-                sortedDist
-            )
+            sampledSets[setI].gather(indexSets[setI])
         );
     }
 }
diff --git a/src/sampling/sampledSet/shortestPath/shortestPathSet.C b/src/sampling/sampledSet/shortestPath/shortestPathSet.C
new file mode 100644
index 0000000000..d289f43b97
--- /dev/null
+++ b/src/sampling/sampledSet/shortestPath/shortestPathSet.C
@@ -0,0 +1,344 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "shortestPathSet.H"
+#include "meshSearch.H"
+#include "DynamicList.H"
+#include "topoDistanceData.H"
+#include "addToRunTimeSelectionTable.H"
+#include "FaceCellWave.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(shortestPathSet, 0);
+    addToRunTimeSelectionTable(sampledSet, shortestPathSet, word);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::label Foam::shortestPathSet::findMinFace
+(
+    const polyMesh& mesh,
+    const label cellI,
+    const List<topoDistanceData>& allFaceInfo,
+    const point& origin
+)
+{
+    const cell& cFaces2 = mesh.cells()[cellI];
+
+    // 1. Get topologically nearest face
+
+    label minDist = labelMax;
+    label minFaceI = -1;
+    forAll(cFaces2, i)
+    {
+        label faceI = cFaces2[i];
+        const topoDistanceData& info = allFaceInfo[faceI];
+        if (info.distance() < minDist)
+        {
+            minDist = info.distance();
+            minFaceI = faceI;
+        }
+    }
+
+    // 2. Check all faces with minDist for minimum distance to origin
+    scalar minDist2 = ROOTVGREAT;
+    forAll(cFaces2, i)
+    {
+        label faceI = cFaces2[i];
+        if (allFaceInfo[faceI].distance() == minDist)
+        {
+            scalar d2 = magSqr(mesh.faceCentres()[faceI]-origin);
+            if (d2 < minDist2)
+            {
+                minDist2 = d2;
+                minFaceI  = faceI;
+            }
+        }
+    }
+
+    return minFaceI;
+}
+
+
+void Foam::shortestPathSet::genSamples(const polyMesh& mesh)
+{
+    // Storage for sample points
+    DynamicList<point> samplingPts;
+    DynamicList<label> samplingCells;
+    DynamicList<label> samplingFaces;
+    DynamicList<label> samplingSegments;
+    DynamicList<scalar> samplingCurveDist;
+
+    forAll(insidePoints_, pointI)
+    {
+        label cell1I = mesh.findCell(insidePoints_[pointI]);
+
+        //
+        // Pass1: Set initial changed faces from cell1 (seed)
+        //
+
+        List<topoDistanceData> faceDist;
+        labelList cFaces1;
+
+        if (cell1I != -1)
+        {
+            cFaces1 = mesh.cells()[cell1I];
+            faceDist.setSize
+            (
+                cFaces1.size(),
+                topoDistanceData(123, 0)
+            );
+        }
+
+        List<topoDistanceData> allFaceInfo(mesh.nFaces());
+        List<topoDistanceData> allCellInfo(mesh.nCells());
+
+        // Walk through face-cell wave till all cells are reached
+        FaceCellWave
+        <
+            topoDistanceData
+        > wallDistCalc
+        (
+            mesh,
+            cFaces1,
+            faceDist,
+            allFaceInfo,
+            allCellInfo,
+            mesh.globalData().nTotalCells()+1   // max iterations
+        );
+
+        // Pass2: walk from outside points backwards. Note: could be done using
+        //        FaceCellWave as well but is overly complex since does
+        //        not allow logic comparing all faces of a cell.
+
+        const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+        // Get the target point
+        label cell2I = mesh.findCell(outsidePoints_[pointI]);
+
+        // The number of cells between cell1 and cell2 is the max number of
+        // iterations to search backward
+        label nPathPoints = 0;
+
+        if (cell2I != -1)
+        {
+            if (!allCellInfo[cell2I].valid(wallDistCalc.data()))
+            {
+                WarningInFunction
+                    << "Point " << outsidePoints_[pointI]
+                    << " not reachable by walk. Probably mesh has "
+                    << " island/regions. Skipped route detection." << endl;
+                return;
+            }
+
+            nPathPoints = allCellInfo[cell2I].distance();
+        }
+        reduce(nPathPoints, maxOp<label>());
+
+        // Start with given target cell and walk back
+        label frontCellI = cell2I;
+
+        while (nPathPoints--)
+        {
+            label frontFaceI = -1;
+
+            // Work within same processor
+            if (frontCellI != -1)
+            {
+                // Find face with lowest distance from seed
+                //   x  |  x  2  1  2  2  |  x  |  x
+                //  --- + --- + -1- + -2- + --- + ---
+                //   x  |  1  1  0  1  1  |  x  |  x
+                //  --- + --- + -1- + -2- + --- + ---
+                //   x  |  x  2  1  2  2  3  3  4  4
+                //  --- + --- + --- + -3- + -4- + -5-
+                //   x  |  x  3  2  3  3  4  4  5  5
+                // e.g. if we start from cell with value = 4, we have neighbour
+                // faces 4, 4, 5, 5. Choose 4 (least distance to seed)
+                // and continue...
+
+                frontFaceI = findMinFace
+                (
+                    mesh,
+                    frontCellI,
+                    allFaceInfo,
+                    outsidePoints_[pointI]
+                );
+
+                // Loop until we hit a boundary face
+                while (mesh.isInternalFace(frontFaceI))
+                {
+                    // Step to neighbouring cell
+                    label nbrCellI = mesh.faceOwner()[frontFaceI];
+                    if (nbrCellI == frontCellI)
+                    {
+                        nbrCellI = mesh.faceNeighbour()[frontFaceI];
+                    }
+
+                    if (nbrCellI == cell1I)
+                    {
+                        // Pout<< " Found connection seed cell!" << endl;
+                        frontCellI = -1;
+                        break;
+                    }
+
+                    frontCellI = nbrCellI;
+
+                    // Pick best face on cell
+                    frontFaceI = findMinFace
+                    (
+                        mesh,
+                        frontCellI,
+                        allFaceInfo,
+                        outsidePoints_[pointI]
+                    );
+
+                    // Set the sampling point
+                    samplingPts.append(mesh.cellCentres()[frontCellI]);
+                    samplingCells.append(frontCellI);
+                    samplingFaces.append(-1);
+                    samplingSegments.append(pointI);
+                    //Check if mag of distance is useful
+                    samplingCurveDist.append(nPathPoints);
+                }
+            }
+
+            // Situation 1: we found the destination cell (do nothing)
+            if (!returnReduce(frontCellI != -1, orOp<bool>()))
+            {
+                break;
+            }
+
+            // Situation 2: we're on a coupled patch and might need to
+            //              switch processor/cell
+            boolList isFront(mesh.nFaces()-mesh.nInternalFaces(), false);
+
+            if (frontFaceI != -1)
+            {
+                isFront[frontFaceI-mesh.nInternalFaces()] = true;
+            }
+            syncTools::swapBoundaryFaceList(mesh, isFront);
+
+            frontCellI = -1;
+            forAll(pbm, patchI)
+            {
+                const polyPatch& pp = pbm[patchI];
+                forAll(pp, i)
+                {
+                    label faceI = pp.start()+i;
+                    if (isFront[faceI-mesh.nInternalFaces()])
+                    {
+                        frontCellI = pp.faceCells()[i];
+                        break;
+                    }
+                }
+
+                if (frontCellI != -1)
+                {
+                    samplingPts.append(mesh.cellCentres()[frontCellI]);
+                    samplingCells.append(frontCellI);
+                    samplingFaces.append(-1);
+                    samplingSegments.append(pointI);
+                    samplingCurveDist.append(nPathPoints);
+                    break;
+                }
+            }
+        }
+    }
+    samplingPts.shrink();
+    samplingCells.shrink();
+    samplingFaces.shrink();
+    samplingSegments.shrink();
+    samplingCurveDist.shrink();
+
+    setSamples
+    (
+        samplingPts,
+        samplingCells,
+        samplingFaces,
+        samplingSegments,
+        samplingCurveDist
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::shortestPathSet::shortestPathSet
+(
+    const word& name,
+    const polyMesh& mesh,
+    const meshSearch& searchEngine,
+    const word& axis,
+    const pointField& insidePoints,
+    const pointField& outsidePoints
+)
+:
+    sampledSet(name, mesh, searchEngine, axis),
+    insidePoints_(insidePoints),
+    outsidePoints_(outsidePoints)
+{
+    genSamples(mesh);
+
+    if (debug)
+    {
+        write(Info);
+    }
+}
+
+
+Foam::shortestPathSet::shortestPathSet
+(
+    const word& name,
+    const polyMesh& mesh,
+    const meshSearch& searchEngine,
+    const dictionary& dict
+)
+:
+    sampledSet(name, mesh, searchEngine, dict),
+    insidePoints_(dict.lookup("insidePoints")),
+    outsidePoints_(dict.lookup("outsidePoints"))
+{
+    genSamples(mesh);
+
+    if (debug)
+    {
+        write(Info);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::shortestPathSet::~shortestPathSet()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSet/shortestPath/shortestPathSet.H b/src/sampling/sampledSet/shortestPath/shortestPathSet.H
new file mode 100644
index 0000000000..0cfa5b1598
--- /dev/null
+++ b/src/sampling/sampledSet/shortestPath/shortestPathSet.H
@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::shortestPathSet
+
+Description
+    Finds shortest path (in terms of cell centres) to walk on mesh from
+    any point in insidePoints to any point in outsidePoints.
+
+Usage
+    Example of function object specification:
+    \verbatim
+    leakFind
+    {
+        type            sets;
+
+        writeControl    timeStep;
+        interpolationScheme cell;
+        setFormat       vtk;
+
+        sets
+        (
+            leakFind
+            {
+                type    shortestPath;
+                insidePoints   ((0.08 -0.020  -0.005) (-0.05 -0.020  -0.005));
+                outsidePoints  ((-0.08 -0.020  -0.005)(0.05 -0.020  -0.005));
+                axis    xyz;
+            }
+        );
+
+        // Needs at least one field
+        fields          ( p );
+    }
+    \endverbatim
+
+SourceFiles
+    shortestPathSet.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef shortestPathSet_H
+#define shortestPathSet_H
+
+#include "sampledSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class topoDistanceData;
+
+/*---------------------------------------------------------------------------*\
+                       Class shortestPathSet Declaration
+\*---------------------------------------------------------------------------*/
+
+class shortestPathSet
+:
+    public sampledSet
+{
+    // Private data
+
+        //- Originating set of points
+        const pointField insidePoints_;
+
+        //- Destination set of points
+        const pointField outsidePoints_;
+
+
+    // Private Member Functions
+
+        //- Get face with least distance along route
+        static label findMinFace
+        (
+            const polyMesh& mesh,
+            const label cellI,
+            const List<topoDistanceData>& allFaceInfo,
+            const point& origin
+        );
+
+        //- Generate whole path
+        void genSamples(const polyMesh& mesh);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("shortestPath");
+
+
+    // Constructors
+
+        //- Construct from components
+        shortestPathSet
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const meshSearch& searchEngine,
+            const word& axis,
+            const pointField& insidePoints,
+            const pointField& outsidePoints
+        );
+
+        //- Construct from dictionary
+        shortestPathSet
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const meshSearch& searchEngine,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~shortestPathSet();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //