ENH: improved handling of bounded sampled planes (issue #714)

- now warn about the following:
  * the bounding box does not overlap wih the global mesh
  * plane does not intersect the (valid) bounding box
  * plane does not intersect the global mesh

- add bounding to the "plane" variant of a sampled plane.

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C

@@ -114,7 +114,7 @@ void Foam::sampledCuttingPlane::createGeometry()
                 false
             ),
             mesh,
-            dimensionedScalar("zero", dimLength, 0)
+            dimensionedScalar(dimLength)
         )
     );
     volScalarField& cellDistance = cellDistancePtr_();
@@ -212,7 +212,7 @@ void Foam::sampledCuttingPlane::createGeometry()
                 false
             ),
             pointMesh::New(mesh),
-            dimensionedScalar("zero", dimLength, 0)
+            dimensionedScalar(dimLength)
         );
         pDist.primitiveFieldRef() = pointDistance_;
 
@@ -220,7 +220,6 @@ void Foam::sampledCuttingPlane::createGeometry()
         pDist.write();
     }
 
-
     //- Direct from cell field and point field.
     isoSurfPtr_.reset
     (
@@ -244,6 +243,43 @@ void Foam::sampledCuttingPlane::createGeometry()
         //)
     );
 
+    // Verify specified bounding box
+    if (!bounds_.empty())
+    {
+        // Bounding box does not overlap with (global) mesh!
+        if (!bounds_.overlaps(mesh.bounds()))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Bounds " << bounds_
+                << " do not overlap the mesh bounding box " << mesh.bounds()
+                << nl << endl;
+        }
+
+        // Plane does not intersect the bounding box
+        if (!bounds_.intersects(plane_))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Plane "<< plane_ << " does not intersect the bounds "
+                << bounds_
+                << nl << endl;
+        }
+    }
+
+    // Plane does not intersect the (global) mesh!
+    if (!mesh.bounds().intersects(plane_))
+    {
+        WarningInFunction
+            << nl
+            << name() << " : "
+            << "Plane "<< plane_ << " does not intersect the mesh bounds "
+            << mesh.bounds()
+            << nl << endl;
+    }
+
     if (debug)
     {
         print(Pout);
@@ -297,12 +333,6 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledCuttingPlane::~sampledCuttingPlane()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::sampledCuttingPlane::needsUpdate() const

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H

@@ -133,7 +133,7 @@ public:
 
 
     //- Destructor
-    virtual ~sampledCuttingPlane();
+    virtual ~sampledCuttingPlane() = default;
 
 
     // Member Functions

src/sampling/sampledSurface/sampledPlane/sampledPlane.C

@@ -52,12 +52,13 @@ Foam::sampledPlane::sampledPlane
     sampledSurface(name, mesh),
     cuttingPlane(planeDesc),
     zoneKey_(zoneKey),
+    bounds_(),
     triangulate_(triangulate),
     needsUpdate_(true)
 {
     if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1)
     {
-        Info<< "cellZone " << zoneKey_
+        Info<< "cellZone(s) " << zoneKey_
             << " not found - using entire mesh" << endl;
     }
 }
@@ -72,7 +73,8 @@ Foam::sampledPlane::sampledPlane
 :
     sampledSurface(name, mesh, dict),
     cuttingPlane(plane(dict)),
-    zoneKey_(keyType::null),
+    zoneKey_(dict.lookupOrDefault<keyType>("zone", keyType::null)),
+    bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)),
     triangulate_(dict.lookupOrDefault("triangulate", true)),
     needsUpdate_(true)
 {
@@ -82,29 +84,21 @@ Foam::sampledPlane::sampledPlane
     {
         coordinateSystem cs(mesh, dict.subDict("coordinateSystem"));
 
-        point  base = cs.globalPosition(planeDesc().refPoint());
+        const point  base = cs.globalPosition(planeDesc().refPoint());
-        vector norm = cs.globalVector(planeDesc().normal());
+        const vector norm = cs.globalVector(planeDesc().normal());
 
         // Assign the plane description
         static_cast<plane&>(*this) = plane(base, norm);
     }
 
-    dict.readIfPresent("zone", zoneKey_);
-
     if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1)
     {
-        Info<< "cellZone " << zoneKey_
+        Info<< "cellZone(s) " << zoneKey_
             << " not found - using entire mesh" << endl;
     }
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledPlane::~sampledPlane()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::sampledPlane::needsUpdate() const
@@ -137,9 +131,88 @@ bool Foam::sampledPlane::update()
 
     sampledSurface::clearGeom();
 
+    const plane& pln = static_cast<const plane&>(*this);
+
+    // Verify specified bounding box
+    if (!bounds_.empty())
+    {
+        // Bounding box does not overlap with (global) mesh!
+        if (!bounds_.overlaps(mesh().bounds()))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Bounds " << bounds_
+                << " do not overlap the mesh bounding box " << mesh().bounds()
+                << nl << endl;
+        }
+
+        // Plane does not intersect the bounding box
+        if (!bounds_.intersects(pln))
+        {
+            WarningInFunction
+                << nl
+                << name() << " : "
+                << "Plane "<< pln << " does not intersect the bounds "
+                << bounds_
+                << nl << endl;
+        }
+    }
+
+    // Plane does not intersect the (global) mesh!
+    if (!mesh().bounds().intersects(pln))
+    {
+        WarningInFunction
+            << nl
+            << name() << " : "
+            << "Plane "<< pln << " does not intersect the mesh bounds "
+            << mesh().bounds()
+            << nl << endl;
+    }
+
     labelList selectedCells = mesh().cellZones().findMatching(zoneKey_).used();
 
-    if (returnReduce(selectedCells.empty(), andOp<bool>()))
+    bool fullMesh = returnReduce(selectedCells.empty(), andOp<bool>());
+
+    if (!bounds_.empty())
+    {
+        const auto& cellCentres = static_cast<const fvMesh&>(mesh()).C();
+
+        if (fullMesh)
+        {
+            const label len = mesh().nCells();
+
+            selectedCells.setSize(len);
+
+            label count = 0;
+            for (label celli=0; celli < len; ++celli)
+            {
+                if (bounds_.contains(cellCentres[celli]))
+                {
+                    selectedCells[count++] = celli;
+                }
+            }
+
+            selectedCells.setSize(count);
+        }
+        else
+        {
+            label count = 0;
+            for (const label celli : selectedCells)
+            {
+                if (bounds_.contains(cellCentres[celli]))
+                {
+                    selectedCells[count++] = celli;
+                }
+            }
+
+            selectedCells.setSize(count);
+        }
+
+        fullMesh = false;
+    }
+
+    if (fullMesh)
     {
         reCut(mesh(), triangulate_);
     }

src/sampling/sampledSurface/sampledPlane/sampledPlane.H

@@ -59,7 +59,10 @@ class sampledPlane
     // Private data
 
         //- If restricted to zones, name of this zone or a regular expression
-        keyType zoneKey_;
+        const keyType zoneKey_;
+
+        //- Optional bounding box to trim triangles against
+        const boundBox bounds_;
 
         //- Triangulated faces or keep faces as is
         const bool triangulate_;
@@ -110,7 +113,7 @@ public:
 
 
     //- Destructor
-    virtual ~sampledPlane();
+    virtual ~sampledPlane() = default;
 
 
     // Member Functions

src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C

@@ -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) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,8 +45,10 @@ Foam::sampledPlane::interpolateField
     const interpolation<Type>& interpolator
 ) const
 {
-    // One value per point
+    // One value per point.
-    tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
+    // Initialize with Zero to handle missed/degenerate faces
+
+    tmp<Field<Type>> tvalues(new Field<Type>(points().size(), Zero));
     Field<Type>& values = tvalues.ref();
 
     boolList pointDone(points().size(), false);
@@ -54,19 +56,18 @@ Foam::sampledPlane::interpolateField
     forAll(faces(), cutFacei)
     {
         const face& f = faces()[cutFacei];
+        const label celli = meshCells()[cutFacei];
 
-        forAll(f, faceVertI)
+        for (const label pointi : f)
         {
-            label pointi = f[faceVertI];
-
             if (!pointDone[pointi])
             {
+                pointDone[pointi] = true;
                 values[pointi] = interpolator.interpolate
                 (
                     points()[pointi],
-                    meshCells()[cutFacei]
+                    celli
                 );
-                pointDone[pointi] = true;
             }
         }
     }

src/sampling/surface/cuttingPlane/cuttingPlane.C

@@ -33,8 +33,8 @@ License
 // Set values for what is close to zero and what is considered to
 // be positive (and not just rounding noise)
 //! \cond localScope
-const Foam::scalar zeroish  = Foam::SMALL;
+static const Foam::scalar zeroish  = Foam::SMALL;
-const Foam::scalar positive = Foam::SMALL * 1E3;
+static const Foam::scalar positive = Foam::SMALL * 1E3;
 //! \endcond
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -49,31 +49,25 @@ void Foam::cuttingPlane::calcCutCells
     const labelListList& cellEdges = mesh.cellEdges();
     const edgeList& edges = mesh.edges();
 
-    label listSize = cellEdges.size();
+    const label len =
-    if (notNull(cellIdLabels))
+        (notNull(cellIdLabels) ? cellIdLabels.size() : cellEdges.size());
-    {
-        listSize = cellIdLabels.size();
-    }
 
-    meshCells_.setSize(listSize);
+    meshCells_.setSize(len);
     label cutcelli(0);
 
     // Find the cut cells by detecting any cell that uses points with
     // opposing dotProducts.
-    for (label listI = 0; listI < listSize; ++listI)
+    for (label listi = 0; listi < len; ++listi)
     {
-        label celli = listI;
+        const label celli =
-
+            (notNull(cellIdLabels) ? cellIdLabels[listi] : listi);
-        if (notNull(cellIdLabels))
-        {
-            celli = cellIdLabels[listI];
-        }
 
         const labelList& cEdges = cellEdges[celli];
         label nCutEdges = 0;
-        forAll(cEdges, i)
+
+        for (const label edgei : cEdges)
         {
-            const edge& e = edges[cEdges[i]];
+            const edge& e = edges[edgei];
 
             if
             (
@@ -81,9 +75,7 @@ void Foam::cuttingPlane::calcCutCells
              || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
             )
             {
-                nCutEdges++;
+                if (++nCutEdges > 2)
-
-                if (nCutEdges > 2)
                 {
                     meshCells_[cutcelli++] = celli;
                     break;
@@ -129,7 +121,7 @@ void Foam::cuttingPlane::intersectEdges
             const point& p0 = points[e[0]];
             const point& p1 = points[e[1]];
 
-            scalar alpha = lineIntersect(linePointRef(p0, p1));
+            const scalar alpha = lineIntersect(linePointRef(p0, p1));
 
             if (alpha < zeroish)
             {
@@ -186,10 +178,8 @@ bool Foam::cuttingPlane::walkCell
         // If so should e.g. decompose the cells on both faces and redo
         // the calculation.
 
-        forAll(fEdges, i)
+        for (const label edge2I : fEdges)
         {
-            label edge2I = fEdges[i];
-
             if (edge2I != edgeI && edgePoint[edge2I] != -1)
             {
                 nextEdgeI = edge2I;
@@ -212,7 +202,7 @@ bool Foam::cuttingPlane::walkCell
 
         edgeI = nextEdgeI;
 
-        nIter++;
+        ++nIter;
 
         if (nIter > 1000)
         {
@@ -232,8 +222,7 @@ bool Foam::cuttingPlane::walkCell
     {
         return true;
     }
-    else
+
-    {
     WarningInFunction
         << "Did not find closed walk along surface of cell " << celli
         << " starting from edge " << startEdgeI << nl
@@ -242,7 +231,6 @@ bool Foam::cuttingPlane::walkCell
 
     return false;
 }
-}
 
 
 void Foam::cuttingPlane::walkCellCuts
@@ -261,19 +249,15 @@ void Foam::cuttingPlane::walkCellCuts
     // scratch space for calculating the face vertices
     DynamicList<label> faceVerts(10);
 
-    forAll(meshCells_, i)
+    for (const label celli : meshCells_)
     {
-        label celli = meshCells_[i];
-
         // Find the starting edge to walk from.
         const labelList& cEdges = mesh.cellEdges()[celli];
 
         label startEdgeI = -1;
 
-        forAll(cEdges, cEdgeI)
+        for (const label edgeI : cEdges)
         {
-            label edgeI = cEdges[cEdgeI];
-
             if (edgePoint[edgeI] != -1)
             {
                 startEdgeI = edgeI;
@@ -309,7 +293,7 @@ void Foam::cuttingPlane::walkCellCuts
                 f.flip();
             }
 
-            // the cut faces are usually quite ugly, so optionally triangulate
+            // The cut faces can be quite ugly, so optionally triangulate
             if (triangulate)
             {
                 label nTri = f.triangles(cutPoints, dynCutFaces);
@@ -326,9 +310,19 @@ void Foam::cuttingPlane::walkCellCuts
         }
     }
 
+    // No cuts? Then no need for any of this information
+    if (dynCutCells.empty())
+    {
+        this->storedPoints().clear();
+        this->storedFaces().clear();
+        meshCells_.clear();
+    }
+    else
+    {
         this->storedFaces().transfer(dynCutFaces);
         meshCells_.transfer(dynCutCells);
     }
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -353,7 +347,6 @@ Foam::cuttingPlane::cuttingPlane
 }
 
 
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::cuttingPlane::reCut
@@ -391,6 +384,7 @@ void Foam::cuttingPlane::remapFaces
     {
         MeshStorage::remapFaces(faceMap);
 
+        // Renumber
         List<label> newCutCells(faceMap.size());
         forAll(faceMap, facei)
         {

src/sampling/surface/cuttingPlane/cuttingPlane.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) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -116,7 +116,7 @@ protected:
     // Constructors
 
         //- Construct plane description without cutting
-        cuttingPlane(const plane&);
+        cuttingPlane(const plane& pln);
 
 
     // Protected Member Functions
@@ -141,8 +141,8 @@ public:
         //  possibly restricted to a list of cells
         cuttingPlane
         (
-            const plane&,
+            const plane& pln,
-            const primitiveMesh&,
+            const primitiveMesh& mesh,
             const bool triangulate,
             const labelUList& cellIdLabels = labelUList::null()
         );
@@ -178,7 +178,8 @@ public:
 
     // Member Operators
 
-        void operator=(const cuttingPlane&);
+        //- Copy assignment
+        void operator=(const cuttingPlane& rhs);
 };