From a8ef9e9742ba0ba4edeb925ec1e3c403aca73c77 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 14 Sep 2018 09:34:14 +0200
Subject: [PATCH] ENH: make cuttingPlane cell walker an algorithm

- takes two general actions:
  1. orient edge in canonical direction (positive gradient) and detect
     any edge intersection.
  2. edge intersection alpha (0-1)

- refactor into a cuttingSurfaceBase intermediate class with the
  actions as templated parameters rather than function pointers. This
  allows the use of lambda functions with captures from the caller.
---
 src/sampling/Make/files                       |   2 +-
 src/sampling/surface/cutting/cuttingPlane.C   |  87 ++++----
 src/sampling/surface/cutting/cuttingPlane.H   |  82 +-------
 .../surface/cutting/cuttingSurfaceBase.C      | 109 ++++++++++
 .../surface/cutting/cuttingSurfaceBase.H      | 187 ++++++++++++++++++
 ...neWalk.C => cuttingSurfaceBaseTemplates.C} |  43 +---
 6 files changed, 348 insertions(+), 162 deletions(-)
 create mode 100644 src/sampling/surface/cutting/cuttingSurfaceBase.C
 create mode 100644 src/sampling/surface/cutting/cuttingSurfaceBase.H
 rename src/sampling/surface/cutting/{cuttingPlaneWalk.C => cuttingSurfaceBaseTemplates.C} (91%)

diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 7d3665f9a4..5e6a2e2ac6 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -21,7 +21,7 @@ sampledSet/shortestPath/shortestPathSet.C
 
 surface/cutting/cuttingPlane.C
 surface/cutting/cuttingPlaneCuts.C
-surface/cutting/cuttingPlaneWalk.C
+surface/cutting/cuttingSurfaceBase.C
 surface/distanceSurface/distanceSurface.C
 surface/isoSurface/isoSurface.C
 surface/isoSurface/isoSurfaceCell.C
diff --git a/src/sampling/surface/cutting/cuttingPlane.C b/src/sampling/surface/cutting/cuttingPlane.C
index 46f41212d9..30d6e29fe7 100644
--- a/src/sampling/surface/cutting/cuttingPlane.C
+++ b/src/sampling/surface/cutting/cuttingPlane.C
@@ -25,11 +25,6 @@ License
 
 #include "cuttingPlane.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-int Foam::cuttingPlane::debug(Foam::debug::debugSwitch("cuttingPlane", 0));
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::cuttingPlane::cuttingPlane(const plane& pln)
@@ -89,6 +84,9 @@ void Foam::cuttingPlane::performCut
     bitSet&& cellIdLabels
 )
 {
+    const plane& pln = *this;
+    const pointField& pts = mesh.points();
+
     MeshStorage::clear();
     meshCells_.clear();
 
@@ -107,57 +105,40 @@ void Foam::cuttingPlane::performCut
     // - some ambiguity when plane is exactly between cells
     const label nFaceCuts = calcCellCuts(mesh, sides, cellCuts);
 
-    // Find closed loop from cell cuts
-    walkCellCuts(mesh, cellCuts, sides, triangulate, nFaceCuts);
-}
 
+    // Walk cell cuts to create faces
 
-void Foam::cuttingPlane::performCut
-(
-    const primitiveMesh& mesh,
-    const bool triangulate,
-    const bitSet& cellIdLabels
-)
-{
-    bitSet subsetCells(cellIdLabels);
-
-    performCut(mesh, triangulate, std::move(subsetCells));
-}
-
-
-void Foam::cuttingPlane::performCut
-(
-    const primitiveMesh& mesh,
-    const bool triangulate,
-    const labelUList& cellIdLabels
-)
-{
-    bitSet subsetCells;
-
-    if (notNull(cellIdLabels))
-    {
-        // Pre-populate with restriction
-        subsetCells.resize(mesh.nCells());
-        subsetCells.set(cellIdLabels);
-    }
-
-    performCut(mesh, triangulate, std::move(subsetCells));
-}
-
-
-void Foam::cuttingPlane::remapFaces(const labelUList& faceMap)
-{
-    if (notNull(faceMap) && !faceMap.empty())
-    {
-        MeshStorage::remapFaces(faceMap);
-
-        List<label> remappedCells(faceMap.size());
-        forAll(faceMap, facei)
+    // Action #1:
+    // - Orient edge so it points in the positive normal direction.
+    // - Edge/plane intersection when the sign changes
+    const auto edgeOrientIntersect =
+        [=](edge& e) -> bool
         {
-            remappedCells[facei] = meshCells_[faceMap[facei]];
-        }
-        meshCells_.transfer(remappedCells);
-    }
+            if (sides[e.last()] < sides[e.first()])
+            {
+                e.flip();
+            }
+
+            return sides[e.first()] != sides[e.last()];
+        };
+
+    // Action #2:
+    // - The edge intersection alpha
+    const auto edgeAlphaIntersect =
+        [=](const edge& e) -> scalar
+        {
+            return pln.lineIntersect(e.line(pts));
+        };
+
+    walkCellCuts
+    (
+        mesh,
+        cellCuts,
+        edgeOrientIntersect,
+        edgeAlphaIntersect,
+        triangulate,
+        nFaceCuts
+    );
 }
 
 
diff --git a/src/sampling/surface/cutting/cuttingPlane.H b/src/sampling/surface/cutting/cuttingPlane.H
index 305441f1a2..11667e67f4 100644
--- a/src/sampling/surface/cutting/cuttingPlane.H
+++ b/src/sampling/surface/cutting/cuttingPlane.H
@@ -39,18 +39,13 @@ SourceFiles
 #define cuttingPlane_H
 
 #include "plane.H"
-#include "faceList.H"
-#include "bitSet.H"
-#include "MeshedSurface.H"
+#include "cuttingSurfaceBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Forward declarations
-class primitiveMesh;
-
 /*---------------------------------------------------------------------------*\
                         Class cuttingPlane Declaration
 \*---------------------------------------------------------------------------*/
@@ -58,18 +53,8 @@ class primitiveMesh;
 class cuttingPlane
 :
     public plane,
-    public MeshedSurface<face>
+    public cuttingSurfaceBase
 {
-    //- Private typedef for convenience
-    typedef MeshedSurface<face> MeshStorage;
-
-
-    // Private data
-
-        //- List of cells cut by the plane
-        labelList meshCells_;
-
-
     // Private Member Functions
 
         //- Determine cut cells, possibly restricted to a list of cells
@@ -88,21 +73,6 @@ class cuttingPlane
             bitSet& cellCuts
         );
 
-        //- Walk the cell cuts to create faces
-        //
-        //  \param cellCuts [in] The cells to walk.
-        //  \param sides [in] For each mesh point, the encoded side of the
-        //      plane (0=BACK, 1=ONPLANE, 2=FRONT), which is used to determine
-        //      the edge cuts
-        void walkCellCuts
-        (
-            const primitiveMesh& mesh,
-            const bitSet& cellCuts,
-            const PackedList<2>& sides,
-            const bool triangulate,
-            const label nFaceCuts = 0
-        );
-
 
 protected:
 
@@ -114,41 +84,21 @@ protected:
 
     // Protected Member Functions
 
-        //- Cut mesh with existing plane, restricted to a list of cells
+        //- Cut mesh, restricted to a list of cells
+        using cuttingSurfaceBase::performCut;
+
+        //- Cut mesh, restricted to a list of cells
         //  Reclaim memory for cellIdLabels
-        void performCut
+        virtual void performCut
         (
             const primitiveMesh& mesh,
             const bool triangulate,
             bitSet&& cellIdLabels
         );
 
-        //- Cut mesh with existing plane, restricted to a list of cells
-        void performCut
-        (
-            const primitiveMesh& mesh,
-            const bool triangulate,
-            const bitSet& cellIdLabels = bitSet()
-        );
-
-        //- Cut mesh with existing plane, restricted to a list of cells
-        void performCut
-        (
-            const primitiveMesh& mesh,
-            const bool triangulate,
-            const labelUList& cellIdLabels
-        );
-
-        //- Remap action on triangulation or cleanup
-        virtual void remapFaces(const labelUList& faceMap);
-
 
 public:
 
-    //- Debug information
-    static int debug;
-
-
     // Constructors
 
         //- Construct from plane and mesh reference,
@@ -190,24 +140,6 @@ public:
             return *this;
         }
 
-        //- The mesh cells cut by the plane
-        const labelList& meshCells() const
-        {
-            return meshCells_;
-        }
-
-        //- The mesh cells cut by the plane
-        labelList& meshCells()
-        {
-            return meshCells_;
-        }
-
-        //- Have any cells been cut?
-        bool cut() const
-        {
-            return meshCells_.size();
-        }
-
 
     // Member Operators
 
diff --git a/src/sampling/surface/cutting/cuttingSurfaceBase.C b/src/sampling/surface/cutting/cuttingSurfaceBase.C
new file mode 100644
index 0000000000..699f367da6
--- /dev/null
+++ b/src/sampling/surface/cutting/cuttingSurfaceBase.C
@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "cuttingSurfaceBase.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+int Foam::cuttingSurfaceBase::debug
+(
+    Foam::debug::debugSwitch("cuttingSurfaceBase", 0)
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cuttingSurfaceBase::cuttingSurfaceBase()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::cuttingSurfaceBase::performCut
+(
+    const primitiveMesh& mesh,
+    const bool triangulate,
+    const bitSet& cellIdLabels
+)
+{
+    bitSet subsetCells(cellIdLabels);
+
+    performCut(mesh, triangulate, std::move(subsetCells));
+}
+
+
+void Foam::cuttingSurfaceBase::performCut
+(
+    const primitiveMesh& mesh,
+    const bool triangulate,
+    const labelUList& cellIdLabels
+)
+{
+    bitSet subsetCells;
+
+    if (notNull(cellIdLabels))
+    {
+        // Pre-populate with restriction
+        subsetCells.resize(mesh.nCells());
+        subsetCells.set(cellIdLabels);
+    }
+
+    performCut(mesh, triangulate, std::move(subsetCells));
+}
+
+
+void Foam::cuttingSurfaceBase::remapFaces(const labelUList& faceMap)
+{
+    if (notNull(faceMap) && !faceMap.empty())
+    {
+        MeshStorage::remapFaces(faceMap);
+
+        List<label> remappedCells(faceMap.size());
+        forAll(faceMap, facei)
+        {
+            remappedCells[facei] = meshCells_[faceMap[facei]];
+        }
+        meshCells_.transfer(remappedCells);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::cuttingSurfaceBase::operator=(const cuttingSurfaceBase& rhs)
+{
+    if (this == &rhs)
+    {
+        FatalErrorInFunction
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+
+    static_cast<MeshStorage&>(*this) = rhs;
+    meshCells_ = rhs.meshCells();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/cutting/cuttingSurfaceBase.H b/src/sampling/surface/cutting/cuttingSurfaceBase.H
new file mode 100644
index 0000000000..f645070f02
--- /dev/null
+++ b/src/sampling/surface/cutting/cuttingSurfaceBase.H
@@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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::cuttingSurfaceBase
+
+Description
+    Base for creating a MeshedSurface by performing some type of cell
+    cutting/intersection.
+
+    No attempt at resolving degenerate cases.
+    Since the cut faces can be quite ugly, they will often be triangulated.
+
+SourceFiles
+    cuttingSurfaceBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cuttingSurfaceBase_H
+#define cuttingSurfaceBase_H
+
+#include "bitSet.H"
+#include "faceList.H"
+#include "MeshedSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+class primitiveMesh;
+
+/*---------------------------------------------------------------------------*\
+                        Class cuttingSurfaceBase Declaration
+\*---------------------------------------------------------------------------*/
+
+class cuttingSurfaceBase
+:
+    public MeshedSurface<face>
+{
+protected:
+
+    //- Typedef for convenience
+    typedef MeshedSurface<face> MeshStorage;
+
+
+    // Protected Data
+
+        //- List of the cells cut
+        labelList meshCells_;
+
+
+    // Protected Member Functions
+
+        //- Walk cell cuts to create faces
+        //
+        //  \tparam EdgeOrientIntersect
+        //      Parameter (edge&), returns bool.
+        //      Orient edge for a consistent positive gradient.
+        //      Checks for edge intersection (true|false).
+        //
+        //  \tparam EdgeAlphaIntersect
+        //      Parameter (const edge&), returns scalar.
+        //      Determine alpha [0-1] for an intersecting edge.
+        //      No guarantees when used with non-intersecting edges.
+        //
+        //  \param cellCuts [in] The cells to walk.
+        template<class EdgeOrientIntersect, class EdgeAlphaIntersect>
+        void walkCellCuts
+        (
+            const primitiveMesh& mesh,
+            const bitSet& cellCuts,
+            const EdgeOrientIntersect& edgeOrientIntersect,
+            const EdgeAlphaIntersect&  edgeAlphaIntersect,
+            const bool triangulate,
+            label nFaceCuts = 0
+        );
+
+
+        //- Cut mesh, restricted to a list of cells
+        virtual void performCut
+        (
+            const primitiveMesh& mesh,
+            const bool triangulate,
+            const labelUList& cellIdLabels
+        );
+
+        //- Cut mesh, restricted to a list of cells
+        virtual void performCut
+        (
+            const primitiveMesh& mesh,
+            const bool triangulate,
+            const bitSet& cellSelectionMask = bitSet()
+        );
+
+        //- Cut mesh, restricted to a list of cells
+        //  Reclaim memory for cellSelectionMask
+        virtual void performCut
+        (
+            const primitiveMesh& mesh,
+            const bool triangulate,
+            bitSet&& cellSelectionMask
+        ) = 0;
+
+        //- Remap action on triangulation or cleanup
+        virtual void remapFaces(const labelUList& faceMap);
+
+
+public:
+
+    //- Debug information
+    static int debug;
+
+
+    // Constructors
+
+        //- Construct null
+        cuttingSurfaceBase();
+
+
+    //- Destructors
+    virtual ~cuttingSurfaceBase() = default;
+
+
+    // Member Functions
+
+        //- The mesh cells cut
+        const labelList& meshCells() const
+        {
+            return meshCells_;
+        }
+
+        //- The mesh cells cut
+        labelList& meshCells()
+        {
+            return meshCells_;
+        }
+
+        //- Have any cells been cut?
+        bool cut() const
+        {
+            return meshCells_.size();
+        }
+
+
+    // Member Operators
+
+        //- Copy assignment
+        void operator=(const cuttingSurfaceBase& rhs);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "cuttingSurfaceBaseTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/cutting/cuttingPlaneWalk.C b/src/sampling/surface/cutting/cuttingSurfaceBaseTemplates.C
similarity index 91%
rename from src/sampling/surface/cutting/cuttingPlaneWalk.C
rename to src/sampling/surface/cutting/cuttingSurfaceBaseTemplates.C
index 6265094cba..9249d6f7dc 100644
--- a/src/sampling/surface/cutting/cuttingPlaneWalk.C
+++ b/src/sampling/surface/cutting/cuttingSurfaceBaseTemplates.C
@@ -23,40 +23,19 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "cuttingPlane.H"
+#include "volFields.H"
 #include "edgeHashes.H"
 #include "HashOps.H"
 
-// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-namespace Foam
-{
-    // Check edge/plane intersection based on crossings ... trivial check.
-    // Orients the edge (first,last) points in the positive normal direction
-    inline bool intersectEdgeOrient(const PackedList<2>& sides, edge& e)
-    {
-        if (sides[e.first()] == sides[e.last()])
-        {
-            return false;
-        }
-        else if (sides[e.last()] < sides[e.first()])
-        {
-            e.flip();
-        }
-
-        return true;
-    }
-
-} // End namespace Foam
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::cuttingPlane::walkCellCuts
+template<class EdgeOrientIntersect, class EdgeAlphaIntersect>
+void Foam::cuttingSurfaceBase::walkCellCuts
 (
     const primitiveMesh& mesh,
     const bitSet& cellCuts,
-    const PackedList<2>& sides,
+    const EdgeOrientIntersect& edgeOrientIntersect,
+    const EdgeAlphaIntersect&  edgeAlphaIntersect,
     const bool triangulate,
     label nFaceCuts
 )
@@ -66,7 +45,6 @@ void Foam::cuttingPlane::walkCellCuts
     const cellList& cells = mesh.cells();
     const pointField& points = mesh.points();
 
-
     // Dynamic lists to handle triangulation and/or missed cuts etc
     const label nCellCuts = cellCuts.count();
 
@@ -156,8 +134,8 @@ void Foam::cuttingPlane::walkCellCuts
             {
                 edge e(f.faceEdge(fp));
 
-                // Action #1: orient edge and detect intersection
-                if (!intersectEdgeOrient(sides, e))
+                // Action #1: Orient edge (+ve gradient) and detect intersect
+                if (!edgeOrientIntersect(e))
                 {
                     continue;
                 }
@@ -191,8 +169,7 @@ void Foam::cuttingPlane::walkCellCuts
                 const point& p1 = points[e.last()];
 
                 // Action #2: edge cut alpha
-                const scalar alpha =
-                    this->lineIntersect(linePointRef(p0, p1));
+                const scalar alpha = edgeAlphaIntersect(e);
 
                 if (alpha < SMALL)
                 {
@@ -253,7 +230,7 @@ void Foam::cuttingPlane::walkCellCuts
 
 
         // Handling cuts between two cells
-        // After the previous intersectEdgeOrient call, the edge oriented
+        // After the previous edgeIntersectAndOrient call, the edge is oriented
         // according to the gradient.
         // If we only ever cut at the same edge end we know that we have
         // a cut coinciding with a cell face.