ENH: support different iso-algorithms (#1374)

- add an 'isoAlgorithm' keyword to distance surface and cutting plane
  to advance further testing of the isoSurfaceTopo algorithm.

  Does not yet handle the full spectrum of bound boxes, cellZones etc.

src/sampling/Make/files

@@ -29,6 +29,7 @@ surface/cutting/cuttingSurfaceBase.C
 surface/cutting/cuttingSurfaceBaseSelection.C
 surface/distanceSurface/distanceSurface.C
 surface/isoSurface/isoSurface.C
+surface/isoSurface/isoSurfaceBase.C
 surface/isoSurface/isoSurfaceCell.C
 surface/isoSurface/isoSurfaceTopo.C
 surface/thresholdCellFaces/thresholdCellFaces.C

src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H

@@ -37,25 +37,27 @@ Usage
     (
         surface1
         {
-            type        dDistanceSurface;
+            type        distanceSurface;
         }
     );
     \endverbatim
 
     Where the sub-entries comprise:
     \table
-        Property | Description                             | Required | Default
-        type     | distanceSurface                         | yes      |
-        distance | distance from surface                   | yes      |
-        signed   | apply sign for +ve/-ve distance         | yes      |
-        cell     | use isoSurfaceCell algorithm            | no       | true
-        regularise | point snapping                        | yes      |
-        average  | cell values from averaged point values  | no       | false
-        bounds   | limit with bounding box                 | no       |
-        surfaceType | type of surface                      | yes      |
-        surfaceName | name of surface in triSurface/       | no  | dict name
+        Property | Description                              | Required | Default
+        type     | distanceSurface                          | yes |
+        distance | Distance from surface                    | yes |
+        signed   | Use sign when distance is positive       | partly |
+        isoAlgorithm | (cell/topo/point)                    | no  | cell
+        regularise | Point snapping for iso-surface         | no  | true
+        average  | Cell values from averaged point values   | no  | false
+        bounds   | Limit with bounding box                  | no  |
+        surfaceType | Type of surface                       | yes |
+        surfaceName | Name of surface in \c triSurface/     | no  | dict name
     \endtable
 
+Note
+    For compatibility, the keyword 'cell' (as a bool) is accepted
 
 SourceFiles
     sampledDistanceSurface.C

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C

@@ -102,7 +102,9 @@ void Foam::sampledCuttingPlane::createGeometry()
     }
 
     // Clear any stored topologies
-    isoSurfPtr_.ptr();
+    isoSurfPtr_.clear();
+    isoSurfCellPtr_.clear();
+    isoSurfTopoPtr_.clear();
     pointDistance_.clear();
     cellDistancePtr_.clear();
 
@@ -291,28 +293,54 @@ void Foam::sampledCuttingPlane::createGeometry()
         pDist.write();
     }
 
-    //- Direct from cell field and point field.
-    isoSurfPtr_.reset
-    (
-        new isoSurface
+
+    // Direct from cell field and point field.
+    if (isoAlgo_ == isoSurfaceBase::ALGO_CELL)
+    {
+        isoSurfCellPtr_.reset
         (
-            cellDistance,
-            pointDistance_,
-            0.0,
-            regularise_,
-            bounds_,
-            mergeTol_
-        )
-        //new isoSurfaceCell
-        //(
-        //    mesh,
-        //    cellDistance,
-        //    pointDistance_,
-        //    0.0,
-        //    regularise_,
-        //    mergeTol_
-        //)
-    );
+            new isoSurfaceCell
+            (
+                fvm,
+                cellDistance,
+                pointDistance_,
+                0,
+                regularise_,
+                bounds_,
+                mergeTol_
+            )
+        );
+    }
+    else if (isoAlgo_ == isoSurfaceBase::ALGO_TOPO)
+    {
+        isoSurfTopoPtr_.reset
+        (
+            new isoSurfaceTopo
+            (
+                fvm,
+                cellDistance,
+                pointDistance_,
+                0,
+                (regularise_ ? isoSurfaceTopo::DIAGCELL : isoSurfaceTopo::NONE),
+                bounds_
+            )
+        );
+    }
+    else
+    {
+        isoSurfPtr_.reset
+        (
+            new isoSurface
+            (
+                cellDistance,
+                pointDistance_,
+                0,
+                regularise_,
+                bounds_,
+                mergeTol_
+            )
+        );
+    }
 
     if (debug)
     {
@@ -337,12 +365,23 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
     mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
     regularise_(dict.lookupOrDefault("regularise", true)),
     average_(dict.lookupOrDefault("average", false)),
+    isoAlgo_
+    (
+        isoSurfaceBase::algorithmNames.getOrDefault
+        (
+            "isoAlgorithm",
+            dict,
+            isoSurfaceBase::ALGO_POINT
+        )
+    ),
     zoneNames_(),
     exposedPatchName_(),
     needsUpdate_(true),
     subMeshPtr_(nullptr),
     cellDistancePtr_(nullptr),
-    isoSurfPtr_(nullptr)
+    isoSurfPtr_(nullptr),
+    isoSurfCellPtr_(nullptr),
+    isoSurfTopoPtr_(nullptr)
 {
     if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
     {

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H

@@ -27,7 +27,7 @@ Class
     Foam::sampledCuttingPlane
 
 Description
-    A sampledSurface defined by a plane using the iso-surface algorithm
+    A sampledSurface defined by a plane using an iso-surface algorithm
     to \a cut the mesh.
 
     This is often embedded as part of a sampled surfaces function object.
@@ -51,15 +51,16 @@ Usage
 
     Where the sub-entries comprise:
     \table
-        Property | Description                             | Required | Default
-        type     | cuttingPlane                            | yes      |
-        planeType | plane description (pointAndNormal etc) | yes      |
-        mergeTol | tolerance for merging points            | no       | 1e-6
-        regularise | point snapping                        | no       | true
-        bounds   | limit with bounding box                 | no       |
-        zone     | limit to cell zone (name or regex)      | no       |
-        zones    | limit to cell zones (names, regexs)     | no       |
-        exposedPatchName | name for zone subset            | partly   |
+        Property | Description                              | Required | Default
+        type     | cuttingPlane                             | yes |
+        planeType | plane description (pointAndNormal etc)  | yes |
+        mergeTol | tolerance for merging points             | no  | 1e-6
+        isoAlgorithm | (cell/topo/point)                    | no  | point
+        regularise | point snapping                         | no  | true
+        bounds   | limit with bounding box                  | no  |
+        zone     | limit to cell zone (name or regex)       | no  |
+        zones    | limit to cell zones (names, regexs)      | no  |
+        exposedPatchName | name for zone subset             | partly |
     \endtable
 
 Note
@@ -77,10 +78,11 @@ SourceFiles
 #define sampledCuttingPlane_H
 
 #include "sampledSurface.H"
-#include "isoSurface.H"
-//#include "isoSurfaceCell.H"
 #include "plane.H"
 #include "fvMeshSubset.H"
+#include "isoSurface.H"
+#include "isoSurfaceCell.H"
+#include "isoSurfaceTopo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -112,6 +114,9 @@ class sampledCuttingPlane
         //- Whether to recalculate cell values as average of point values
         const bool average_;
 
+        //- The iso-surface algorithm type
+        const isoSurfaceBase::algorithmType isoAlgo_;
+
         //- The zone or zones in which cutting is to occur
         wordRes zoneNames_;
 
@@ -131,10 +136,15 @@ class sampledCuttingPlane
         //- Distance to points
         scalarField pointDistance_;
 
-        //- Constructed iso surface
-        //autoPtr<isoSurfaceCell> isoSurfCellPtr_;
+        //- Constructed iso surface (ALGO_POINT)
         autoPtr<isoSurface> isoSurfPtr_;
 
+        //- Constructed iso surface (ALGO_CELL)
+        autoPtr<isoSurfaceCell> isoSurfCellPtr_;
+
+        //- Constructed iso surface (ALGO_TOPO)
+        autoPtr<isoSurfaceTopo> isoSurfTopoPtr_;
+
 
     // Private Member Functions
 
@@ -162,6 +172,13 @@ class sampledCuttingPlane
             const interpolation<Type>& interpolator
         ) const;
 
+        //- Interpolates cCoords,pCoords.
+        template<class Type>
+        tmp<Field<Type>> isoSurfaceInterpolate
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
 
 public:
 
@@ -186,12 +203,6 @@ public:
 
     // Member Functions
 
-        //const isoSurfaceCell& surface() const
-        const isoSurface& surface() const
-        {
-            return *isoSurfPtr_;
-        }
-
         //- Does the surface need an update?
         virtual bool needsUpdate() const;
 
@@ -241,21 +252,59 @@ public:
         }
 
 
+        //- The underlying surface
+        const meshedSurface& surface() const
+        {
+            if (isoSurfCellPtr_)
+            {
+                return *isoSurfCellPtr_;
+            }
+            else if (isoSurfTopoPtr_)
+            {
+                return *isoSurfTopoPtr_;
+            }
+            return *isoSurfPtr_;
+        }
+
         //- The underlying surface
         meshedSurface& surface()
         {
+            if (isoSurfCellPtr_)
+            {
+                return *isoSurfCellPtr_;
+            }
+            else if (isoSurfTopoPtr_)
+            {
+                return *isoSurfTopoPtr_;
+            }
             return *isoSurfPtr_;
         }
 
         //- For each face, the original cell in mesh
         const labelList& meshCells() const
         {
+            if (isoSurfCellPtr_)
+            {
+                return isoSurfCellPtr_->meshCells();
+            }
+            else if (isoSurfTopoPtr_)
+            {
+                return isoSurfTopoPtr_->meshCells();
+            }
             return isoSurfPtr_->meshCells();
         }
 
         //- For each face, the original cell in mesh
         labelList& meshCells()
         {
+            if (isoSurfCellPtr_)
+            {
+                return isoSurfCellPtr_->meshCells();
+            }
+            else if (isoSurfTopoPtr_)
+            {
+                return isoSurfTopoPtr_->meshCells();
+            }
             return isoSurfPtr_->meshCells();
         }
 

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2018-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -67,7 +67,7 @@ Foam::sampledCuttingPlane::sampleOnPoints
         auto tpointFld =
             volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
 
-        return surface().interpolate
+        return this->isoSurfaceInterpolate
         (
             (average_ ? pointAverage(tpointFld())() : volSubFld),
             tpointFld()
@@ -78,7 +78,7 @@ Foam::sampledCuttingPlane::sampleOnPoints
     auto tpointFld =
         volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
 
-    return surface().interpolate
+    return this->isoSurfaceInterpolate
     (
         (average_ ? pointAverage(tpointFld())() : volFld),
         tpointFld()
@@ -86,4 +86,27 @@ Foam::sampledCuttingPlane::sampleOnPoints
 }
 
 
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::sampledCuttingPlane::isoSurfaceInterpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+    const Field<Type>& pCoords
+) const
+{
+    if (isoSurfCellPtr_)
+    {
+        return isoSurfCellPtr_->interpolate(cCoords, pCoords);
+    }
+    else if (isoSurfTopoPtr_)
+    {
+        return isoSurfTopoPtr_->interpolate(cCoords, pCoords);
+    }
+    else
+    {
+        return isoSurfPtr_->interpolate(cCoords, pCoords);
+    }
+}
+
+
 // ************************************************************************* //

src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.H

@@ -87,12 +87,12 @@ class sampledThresholdCellFaces
     public sampledSurface,
     public MeshedSurface<face>
 {
-    //- Private typedefs for convenience
+        //- Private typedefs for convenience
         typedef MeshedSurface<face> MeshStorage;
 
     // Private data
 
-        //- Field to get isoSurface of
+        //- Threshold field
         const word fieldName_;
 
         //- Threshold value

src/sampling/surface/distanceSurface/distanceSurface.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -40,6 +40,48 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    static isoSurfaceBase::algorithmType getIsoAlgorithm(const dictionary& dict)
+    {
+        // Previously 'cell' (bool), but now 'isoAlgorithm' (enum)
+
+        // Default (bool) for 1906 and earlier
+        bool useCell = true;
+
+        // Default (enum) after 1906
+        isoSurfaceBase::algorithmType algo = isoSurfaceBase::ALGO_CELL;
+
+        if
+        (
+            !isoSurfaceBase::algorithmNames.readIfPresent
+            (
+                "isoAlgorithm", dict, algo
+            )
+            // When above fails, use 'compat' to also get upgrade messages
+         && dict.readIfPresentCompat
+            (
+                "isoAlgorithm", {{"cell", 1906}}, useCell
+            )
+        )
+        {
+            return
+            (
+                useCell
+              ? isoSurfaceBase::ALGO_CELL
+              : isoSurfaceBase::ALGO_POINT
+            );
+        }
+
+        return algo;
+    }
+
+
+} // End namespace Foam
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::distanceSurface::distanceSurface
@@ -57,7 +99,7 @@ Foam::distanceSurface::distanceSurface
             dict.get<word>("surfaceType"),
             IOobject
             (
-                dict.lookupOrDefault("surfaceName", defaultSurfaceName),
+                dict.getOrDefault("surfaceName", defaultSurfaceName),
                 mesh.time().constant(), // directory
                 "triSurface",           // instance
                 mesh.time(),            // registry
@@ -72,11 +114,12 @@ Foam::distanceSurface::distanceSurface
     (
         distance_ < 0 || equal(distance_, Zero) || dict.get<bool>("signed")
     ),
-    cell_(dict.lookupOrDefault("cell", true)),
-    regularise_(dict.lookupOrDefault("regularise", true)),
-    bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)),
+    regularise_(dict.getOrDefault("regularise", true)),
+    isoAlgo_(getIsoAlgorithm(dict)),
+    bounds_(dict.getOrDefault("bounds", boundBox::invertedBox)),
+    isoSurfPtr_(nullptr),
     isoSurfCellPtr_(nullptr),
-    isoSurfPtr_(nullptr)
+    isoSurfTopoPtr_(nullptr)
 {}
 
 
@@ -88,7 +131,7 @@ Foam::distanceSurface::distanceSurface
     const word& surfaceName,
     const scalar distance,
     const bool signedDistance,
-    const bool cell,
+    const isoSurfaceBase::algorithmType algo,
     const bool regularise,
     const boundBox& bounds
 )
@@ -116,11 +159,12 @@ Foam::distanceSurface::distanceSurface
     (
         signedDistance || distance_ < 0 || equal(distance_, Zero)
     ),
-    cell_(cell),
     regularise_(regularise),
+    isoAlgo_(algo),
     bounds_(bounds),
+    isoSurfPtr_(nullptr),
     isoSurfCellPtr_(nullptr),
-    isoSurfPtr_(nullptr)
+    isoSurfTopoPtr_(nullptr)
 {}
 
 
@@ -134,8 +178,9 @@ void Foam::distanceSurface::createGeometry()
     }
 
     // Clear any stored topologies
-    isoSurfCellPtr_.clear();
     isoSurfPtr_.clear();
+    isoSurfCellPtr_.clear();
+    isoSurfTopoPtr_.clear();
 
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh_);
 
@@ -165,7 +210,11 @@ void Foam::distanceSurface::createGeometry()
     // to limit the extent of open edges.
 
     const bool isZeroDist  = equal(distance_, Zero);
-    const bool filterCells = (cell_ && isZeroDist);
+    const bool filterCells =
+    (
+        isZeroDist
+     && isoAlgo_ != isoSurfaceBase::ALGO_POINT
+    );
 
     bitSet ignoreCells;
     if (filterCells)
@@ -209,6 +258,9 @@ void Foam::distanceSurface::createGeometry()
                     cellBb.clear();
                     cellBb.add(fvm.points(), fvm.cellPoints(i));
 
+                    // Expand slightly to catch corners
+                    cellBb.inflate(0.1);
+
                     if (!cellBb.contains(nearest[i].hitPoint()))
                     {
                         ignoreCells.set(i);
@@ -349,7 +401,7 @@ void Foam::distanceSurface::createGeometry()
 
 
     // Direct from cell field and point field.
-    if (cell_)
+    if (isoAlgo_ == isoSurfaceBase::ALGO_CELL)
     {
         isoSurfCellPtr_.reset
         (
@@ -366,6 +418,22 @@ void Foam::distanceSurface::createGeometry()
             )
         );
     }
+    else if (isoAlgo_ == isoSurfaceBase::ALGO_TOPO)
+    {
+        isoSurfTopoPtr_.reset
+        (
+            new isoSurfaceTopo
+            (
+                fvm,
+                cellDistance,
+                pointDistance_,
+                distance_,
+                (regularise_ ? isoSurfaceTopo::DIAGCELL : isoSurfaceTopo::NONE),
+                bounds_,
+                ignoreCells
+            )
+        );
+    }
     else
     {
         isoSurfPtr_.reset

src/sampling/surface/distanceSurface/distanceSurface.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -28,20 +28,20 @@ Class
 
 Description
     A surface defined by a distance from an input searchable surface.
-    Uses an isoSurfaceCell or an isoSurface algorithm for constructing the
+    Uses an iso-surface algorithm (cell, topo, point) for constructing the
     distance surface.
 
 Usage
     Dictionary controls:
     \table
-        Property | Description                             | Required | Default
-        distance | distance from surface                   | yes      |
-        signed   | Use sign when distance is positive      | partly   |
-        cell     | use isoSurfaceCell algorithm            | no       | true
-        regularise | Point snapping for iso-surface        | yes      |
-        bounds   | Limit with bounding box                 | no       |
-        surfaceType | Type of surface                      | yes      |
-        surfaceName | Name of surface in \c triSurface/    | no  | dict name
+        Property | Description                              | Required | Default
+        distance | distance from surface                    | yes |
+        signed   | Use sign when distance is positive       | partly |
+        isoAlgorithm | (cell/topo/point)                    | no  | cell
+        regularise | Point snapping for iso-surface         | no  | true
+        bounds   | Limit with bounding box                  | no  |
+        surfaceType | Type of surface                       | yes |
+        surfaceName | Name of surface in \c triSurface/     | no  | dict name
     \endtable
 
 Note
@@ -53,6 +53,8 @@ Note
       surface. The resulting surface elements will not, however, contain
       partial cell coverage.
 
+    For compatibility, the keyword 'cell' (as a bool) is accepted
+
 SourceFiles
     distanceSurface.C
 
@@ -63,8 +65,9 @@ SourceFiles