new management for sampledSurface(s)

//- Does the surface need an update?
        virtual bool needsUpdate() const = 0;

        //- Mark the surface as needing an update.
        //  May also free up unneeded data.
        //  Return false if surface was already marked as expired.
        virtual bool expire() = 0;

        //- Update the surface as required.
        //  Do nothing (and return false) if no update was required
        virtual bool update() = 0;

    The constructors for the derived classes should generally start in a
    'expired' condition (ie, needsUpdate() == true) and rely on a subsequent
    call to the update() method to complete the initialization. Delaying the
    final construction as late as possible allows the construction of
    surfaces that may depend on intermediate calculation results (eg,
    iso-surfaces) and also avoids the unnecessary reconstruction of surfaces
    between sampling intervals.

    It is the responsibility of the caller to ensure that the surface
    update() is called before the surface is used.  The update() method
    implementation should do nothing when the surface is already up-to-date.

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -109,26 +109,24 @@ sets
 //      plane : values on plane defined by point, normal.
 //      patch : values on patch.
 //
-// 1] planes are triangulated by default
+// 1] patches are not triangulated by default
-// 2] patches are not triangulated by default
+// 2] planes are always triangulated
+// 3] iso-surfaces are always triangulated
 surfaces
 (
     constantPlane
     {
-        type            plane;
+        type            plane;    // always triangulated
         basePoint       (0.0501 0.0501 0.005);
         normalVector    (0.1 0.1 1);
 
         //- Optional: restrict to a particular zone
         // zoneName        zone1;
-
-        // Optional: whether to leave as faces or triangulate (=default)
-        triangulate     false;
     }
 
     interpolatedPlane
     {
-        type            plane;
+        type            plane;    // always triangulated
         // make plane relative to the coordinateSystem (Cartesian)
         coordinateSystem
         {
@@ -136,7 +134,6 @@ surfaces
         }
         basePoint       (0 0 0);
         normalVector    (0.1 0.1 1);
-        triangulate     false;
         interpolate     true;
     }
 
@@ -145,6 +142,7 @@ surfaces
         type            patch;
         patchName       movingWall;
         // Optional: whether to leave as faces (=default) or triangulate
+        // triangulate     false;
     }
 
     movingWall_interpolated
@@ -153,26 +151,27 @@ surfaces
         patchName       movingWall;
         interpolate     true;
         // Optional: whether to leave as faces (=default) or triangulate
+        // triangulate     false;
     }
 
     interpolatedIso
     {
         // Iso surface for interpolated values only
-        type            isoSurface;
+        type            isoSurface;    // always triangulated
         isoField        rho;
         isoValue        0.5;
         interpolate     true;
-        //regularise      false;    //optional: do not simplify
+        // regularise      false;    // Optional: do not simplify
     }
     constantIso
     {
-        // Iso surface for constant values. Guarantees triangles to not
+        // Iso surface for constant values.
-        // cross cells.
+        // Triangles guaranteed not to cross cells.
-        type            isoSurfaceCell;
+        type            isoSurfaceCell;    // always triangulated
         isoField        rho;
         isoValue        0.5;
         interpolate     false;
-        //regularise      false;    //optional: do not simplify
+        // regularise      false;    // Optional: do not simplify
     }
 );
 

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -31,7 +31,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
 #include "isoSurface.H"
-//#include "isoSurfaceCell.H"
+// #include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -45,10 +45,9 @@ namespace Foam
 
 void Foam::distanceSurface::createGeometry()
 {
-    // Clear any stored topo
+    // Clear any stored topologies
     facesPtr_.clear();
 
-
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
 
     // Distance to cell centres
@@ -232,21 +231,20 @@ Foam::distanceSurface::distanceSurface
     signed_(readBool(dict.lookup("signed"))),
     regularise_(dict.lookupOrDefault("regularise", true)),
     zoneName_(word::null),
+    needsUpdate_(true),
     isoSurfPtr_(NULL),
     facesPtr_(NULL)
 {
-//    label zoneId = -1;
+//    dict.readIfPresent("zone", zoneName_);
-//    if (dict.readIfPresent("zone", zoneName_))
+//
+//    if (debug && zoneName_.size())
 //    {
-//        zoneId = mesh.cellZones().findZoneID(zoneName_);
+//        if (mesh.cellZones().findZoneID(zoneName_) < 0)
-//        if (debug && zoneId < 0)
 //        {
 //            Info<< "cellZone \"" << zoneName_
-//                << "\" not found - using entire mesh"
+//                << "\" not found - using entire mesh" << endl;
-//                << endl;
 //        }
 //    }
-    correct(true);
 }
 
 
@@ -258,13 +256,39 @@ Foam::distanceSurface::~distanceSurface()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::distanceSurface::correct(const bool meshChanged)
+bool Foam::distanceSurface::needsUpdate() const
 {
-    // Only change of mesh changes plane - zone restriction gets lost
+    return needsUpdate_;
-    if (meshChanged)
+}
+
+
+bool Foam::distanceSurface::expire()
+{
+    // Clear any stored topologies
+    facesPtr_.clear();
+
+    // already marked as expired
+    if (needsUpdate_)
     {
-        createGeometry();
+        return false;
     }
+
+    needsUpdate_ = true;
+    return true;
+}
+
+
+bool Foam::distanceSurface::update()
+{
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    createGeometry();
+
+    needsUpdate_ = false;
+    return true;
 }
 
 

src/sampling/sampledSurface/distanceSurface/distanceSurface.H

@@ -71,6 +71,9 @@ class distanceSurface
         //- zone name (if restricted to zones)
         word zoneName_;
 
+        //- Track if the surface needs an update
+        mutable bool needsUpdate_;
+
 
         //- Distance to cell centres
         autoPtr<volScalarField> cellDistancePtr_;
@@ -127,6 +130,18 @@ public:
 
     // Member Functions
 
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surface as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if surface was already marked as expired.
+        virtual bool expire();
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was needed
+        virtual bool update();
+
         //- Points of surface
         virtual const pointField& points() const
         {
@@ -151,9 +166,6 @@ public:
         }
 
 
-        //- Correct for mesh movement and/or field changes
-        virtual void correct(const bool meshChanged);
-
         const isoSurface& surface() const
         {
             return isoSurfPtr_();

src/sampling/sampledSurface/isoSurface/isoSurface.C

@@ -512,7 +512,7 @@ void Foam::isoSurface::calcSnappedCc
                     (
                         false,              // do not check for duplicate tris
                         localTriPoints,
-                        triPointReverseMap,  
+                        triPointReverseMap,
                         triMap
                     )
                 );
@@ -696,7 +696,7 @@ void Foam::isoSurface::calcSnappedPoint
                 (
                     false,                  // do not check for duplicate tris
                     localTriPoints,
-                    triPointReverseMap,  
+                    triPointReverseMap,
                     triMap
                 )
             );
@@ -1084,7 +1084,7 @@ void Foam::isoSurface::walkOrientation
             forAll(fEdges, fp)
             {
                 label edgeI = fEdges[fp];
-    
+
                 // my points:
                 label p0 = tri[fp];
                 label p1 = tri[tri.fcIndex(fp)];
@@ -1121,7 +1121,7 @@ void Foam::isoSurface::walkOrientation
 
         changedFaces.transfer(newChangedFaces);
     }
-}    
+}
 
 
 void Foam::isoSurface::orientSurface
@@ -1146,7 +1146,7 @@ void Foam::isoSurface::orientSurface
         for
         (
             ;
-            seedTriI < surf.size() && flipState[seedTriI] != -1; 
+            seedTriI < surf.size() && flipState[seedTriI] != -1;
             seedTriI++
         )
         {}
@@ -1355,7 +1355,7 @@ Foam::isoSurface::isoSurface
 :
     mesh_(cVals.mesh()),
     iso_(iso),
-    mergeDistance_(mergeTol*mag(mesh_.bounds().max()-mesh_.bounds().min()))
+    mergeDistance_(mergeTol*mesh_.bounds().mag())
 {
     // Determine if any cut through face/cell
     calcCutTypes(cVals, pVals);

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C

@@ -1409,7 +1409,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
 :
     mesh_(mesh),
     iso_(iso),
-    mergeDistance_(mergeTol*mag(mesh.bounds().max()-mesh.bounds().min()))
+    mergeDistance_(mergeTol*mesh.bounds().mag())
 {
     // Determine if cell is tet
     PackedList<1> isTet(mesh_.nCells());

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C

@@ -176,100 +176,105 @@ void Foam::sampledIsoSurface::getIsoFields() const
 }
 
 
-void Foam::sampledIsoSurface::createGeometry() const
+bool Foam::sampledIsoSurface::updateGeometry() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
 
-    if (fvm.time().timeIndex() != storedTimeIndex_)
+    // no update needed
+    if (fvm.time().timeIndex() == prevTimeIndex_)
     {
-        storedTimeIndex_ = fvm.time().timeIndex();
+        return false;
+    }
 
-        getIsoFields();
 
-        // Clear any stored topo
+    prevTimeIndex_ = fvm.time().timeIndex();
-        surfPtr_.clear();
+    getIsoFields();
-        facesPtr_.clear();
 
-        if (average_)
+    // Clear any stored topo
-        {
+    surfPtr_.clear();
-            //- From point field and interpolated cell.
+    facesPtr_.clear();
-            volScalarField cellAvg
+
+    if (average_)
+    {
+        //- From point field and interpolated cell.
+        volScalarField cellAvg
+        (
+            IOobject
             (
-                IOobject
+                "cellAvg",
-                (
+                fvm.time().timeName(),
-                    "cellAvg",
+                fvm.time(),
-                    fvm.time().timeName(),
+                IOobject::NO_READ,
-                    fvm.time(),
+                IOobject::NO_WRITE,
-                    IOobject::NO_READ,
+                false
-                    IOobject::NO_WRITE,
+            ),
-                    false
+            fvm,
-                ),
+            dimensionedScalar("zero", dimless, scalar(0.0))
-                fvm,
+        );
-                dimensionedScalar("zero", dimless, scalar(0.0))
+        labelField nPointCells(fvm.nCells(), 0);
-            );
+        {
-            labelField nPointCells(fvm.nCells(), 0);
+            for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
             {
-                for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
+                const labelList& pCells = fvm.pointCells(pointI);
+
+                forAll(pCells, i)
                 {
-                    const labelList& pCells = fvm.pointCells(pointI);
+                    label cellI = pCells[i];
 
-                    forAll(pCells, i)
+                    cellAvg[cellI] += (*pointFieldPtr_)[pointI];
-                    {
+                    nPointCells[cellI]++;
-                        label cellI = pCells[i];
-
-                        cellAvg[cellI] += (*pointFieldPtr_)[pointI];
-                        nPointCells[cellI]++;
-                    }
                 }
             }
-            forAll(cellAvg, cellI)
-            {
-                cellAvg[cellI] /= nPointCells[cellI];
-            }
-            // Give value to calculatedFvPatchFields
-            cellAvg.correctBoundaryConditions();
-
-            surfPtr_.reset
-            (
-                new isoSurface
-                (
-                    cellAvg,
-                    *pointFieldPtr_,
-                    isoVal_,
-                    regularise_
-                )
-            );
         }
-        else
+        forAll(cellAvg, cellI)
         {
-            //- Direct from cell field and point field.
+            cellAvg[cellI] /= nPointCells[cellI];
-            surfPtr_.reset
+        }
+        // Give value to calculatedFvPatchFields
+        cellAvg.correctBoundaryConditions();
+
+        surfPtr_.reset
+        (
+            new isoSurface
             (
-                new isoSurface
+                cellAvg,
-                (
+                *pointFieldPtr_,
-                    *volFieldPtr_,
+                isoVal_,
-                    *pointFieldPtr_,
+                regularise_
-                    isoVal_,
+            )
-                    regularise_
+        );
-                )
-            );
-        }
-
-
-        if (debug)
-        {
-            Pout<< "sampledIsoSurface::createGeometry() : constructed iso:"
-                << nl
-                << "    regularise     : " << regularise_ << nl
-                << "    average        : " << average_ << nl
-                << "    isoField       : " << isoField_ << nl
-                << "    isoValue       : " << isoVal_ << nl
-                << "    points         : " << points().size() << nl
-                << "    tris           : " << surface().size() << nl
-                << "    cut cells      : " << surface().meshCells().size()
-                << endl;
-        }
     }
+    else
+    {
+        //- Direct from cell field and point field.
+        surfPtr_.reset
+        (
+            new isoSurface
+            (
+                *volFieldPtr_,
+                *pointFieldPtr_,
+                isoVal_,
+                regularise_
+            )
+        );
+    }
+
+
+    if (debug)
+    {
+        Pout<< "sampledIsoSurface::updateGeometry() : constructed iso:"
+            << nl
+            << "    regularise     : " << regularise_ << nl
+            << "    average        : " << average_ << nl
+            << "    isoField       : " << isoField_ << nl
+            << "    isoValue       : " << isoVal_ << nl
+            << "    points         : " << points().size() << nl
+            << "    tris           : " << surface().size() << nl
+            << "    cut cells      : " << surface().meshCells().size()
+            << endl;
+    }
+
+    return true;
 }
 
 
@@ -290,7 +295,7 @@ Foam::sampledIsoSurface::sampledIsoSurface
     zoneName_(word::null),
     surfPtr_(NULL),
     facesPtr_(NULL),
-    storedTimeIndex_(-1),
+    prevTimeIndex_(-1),
     storedVolFieldPtr_(NULL),
     volFieldPtr_(NULL),
     storedPointFieldPtr_(NULL),
@@ -306,15 +311,14 @@ Foam::sampledIsoSurface::sampledIsoSurface
             << " span across cells." << exit(FatalError);
     }
 
-//    label zoneId = -1;
+//    dict.readIfPresent("zone", zoneName_);
-//    if (dict.readIfPresent("zone", zoneName_))
+//
+//    if (debug && zoneName_.size())
 //    {
-//        zoneId = mesh.cellZones().findZoneID(zoneName_);
+//        if (mesh.cellZones().findZoneID(zoneName_) < 0)
-//        if (debug && zoneId < 0)
 //        {
 //            Info<< "cellZone \"" << zoneName_
-//                << "\" not found - using entire mesh"
+//                << "\" not found - using entire mesh" << endl;
-//                << endl;
 //        }
 //    }
 }
@@ -328,14 +332,34 @@ Foam::sampledIsoSurface::~sampledIsoSurface()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::sampledIsoSurface::correct(const bool meshChanged)
+bool Foam::sampledIsoSurface::needsUpdate() const
 {
-    // Only change of mesh changes plane - zone restriction gets lost
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-    if (meshChanged)
+
+    return fvm.time().timeIndex() != prevTimeIndex_;
+}
+
+
+bool Foam::sampledIsoSurface::expire()
+{
+    surfPtr_.clear();
+    facesPtr_.clear();
+
+    // already marked as expired
+    if (prevTimeIndex_ == -1)
     {
-        surfPtr_.clear();
+        return false;
-        facesPtr_.clear();
     }
+
+    // force update
+    prevTimeIndex_ = -1;
+    return true;
+}
+
+
+bool Foam::sampledIsoSurface::update()
+{
+    return updateGeometry();
 }
 
 

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H

@@ -79,8 +79,8 @@ class sampledIsoSurface
 
         // Recreated for every isoSurface
 
-            //- Time at last call
+            //- Time at last call, also track it surface needs an update
-            mutable label storedTimeIndex_;
+            mutable label prevTimeIndex_;
 
             //- Cached volfield
             mutable autoPtr<volScalarField> storedVolFieldPtr_;
@@ -98,7 +98,8 @@ class sampledIsoSurface
         void getIsoFields() const;
 
         //- Create iso surface (if time has changed)
-        void createGeometry() const;
+        //  Do nothing (and return false) if no update was needed
+        bool updateGeometry() const;
 
         //- sample field on faces
         template <class Type>
@@ -137,6 +138,19 @@ public:
 
     // Member Functions
 
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surface as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if surface was already marked as expired.
+        virtual bool expire();
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was needed
+        virtual bool update();
+
+
         //- Points of surface
         virtual const pointField& points() const
         {
@@ -160,9 +174,6 @@ public:
             return facesPtr_;
         }
 
-        //- Correct for mesh movement and/or field changes
-        virtual void correct(const bool meshChanged);
-
 
         const isoSurface& surface() const
         {

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C

@@ -42,142 +42,147 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::sampledIsoSurfaceCell::createGeometry() const
+bool Foam::sampledIsoSurfaceCell::updateGeometry() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
 
-    if (fvm.time().timeIndex() != storedTimeIndex_)
+    // no update needed
+    if (fvm.time().timeIndex() == prevTimeIndex_)
     {
-        storedTimeIndex_ = fvm.time().timeIndex();
+        return false;
+    }
 
-        // Clear any stored topo
+    prevTimeIndex_ = fvm.time().timeIndex();
-        facesPtr_.clear();
 
-        // Optionally read volScalarField
+    // Clear any stored topo
-        autoPtr<volScalarField> readFieldPtr_;
+    facesPtr_.clear();
 
-        // 1. see if field in database
+    // Optionally read volScalarField
-        // 2. see if field can be read
+    autoPtr<volScalarField> readFieldPtr_;
-        const volScalarField* cellFldPtr = NULL;
+
-        if (fvm.foundObject<volScalarField>(isoField_))
+    // 1. see if field in database
+    // 2. see if field can be read
+    const volScalarField* cellFldPtr = NULL;
+    if (fvm.foundObject<volScalarField>(isoField_))
+    {
+        if (debug)
         {
-            if (debug)
+            Info<< "sampledIsoSurfaceCell::updateGeometry() : lookup "
-            {
+                << isoField_ << endl;
-                Info<< "sampledIsoSurfaceCell::createGeometry() : lookup "
-                    << isoField_ << endl;
-            }
-
-            cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_);
-        }
-        else
-        {
-            // Bit of a hack. Read field and store.
-
-            if (debug)
-            {
-                Info<< "sampledIsoSurfaceCell::createGeometry() : reading "
-                    << isoField_ << " from time " <<fvm.time().timeName()
-                    << endl;
-            }
-
-            readFieldPtr_.reset
-            (
-                new volScalarField
-                (
-                    IOobject
-                    (
-                        isoField_,
-                        fvm.time().timeName(),
-                        fvm,
-                        IOobject::MUST_READ,
-                        IOobject::NO_WRITE,
-                        false
-                    ),
-                    fvm
-                )
-            );
-
-            cellFldPtr = readFieldPtr_.operator->();
-        }
-        const volScalarField& cellFld = *cellFldPtr;
-
-        tmp<pointScalarField> pointFld
-        (
-            volPointInterpolation::New(fvm).interpolate(cellFld)
-        );
-
-        if (average_)
-        {
-            //- From point field and interpolated cell.
-            scalarField cellAvg(fvm.nCells(), scalar(0.0));
-            labelField nPointCells(fvm.nCells(), 0);
-            {
-                for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
-                {
-                    const labelList& pCells = fvm.pointCells(pointI);
-
-                    forAll(pCells, i)
-                    {
-                        label cellI = pCells[i];
-
-                        cellAvg[cellI] += pointFld().internalField()[pointI];
-                        nPointCells[cellI]++;
-                    }
-                }
-            }
-            forAll(cellAvg, cellI)
-            {
-                cellAvg[cellI] /= nPointCells[cellI];
-            }
-
-            const isoSurfaceCell iso
-            (
-                fvm,
-                cellAvg,
-                pointFld().internalField(),
-                isoVal_,
-                regularise_
-            );
-
-            const_cast<sampledIsoSurfaceCell&>
-            (
-                *this
-            ).triSurface::operator=(iso);
-            meshCells_ = iso.meshCells();
-        }
-        else
-        {
-            //- Direct from cell field and point field. Gives bad continuity.
-            const isoSurfaceCell iso
-            (
-                fvm,
-                cellFld.internalField(),
-                pointFld().internalField(),
-                isoVal_,
-                regularise_
-            );
-
-            const_cast<sampledIsoSurfaceCell&>
-            (
-                *this
-            ).triSurface::operator=(iso);
-            meshCells_ = iso.meshCells();
         }
 
+        cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_);
+    }
+    else
+    {
+        // Bit of a hack. Read field and store.
 
         if (debug)
         {
-            Pout<< "sampledIsoSurfaceCell::createGeometry() : constructed iso:"
+            Info<< "sampledIsoSurfaceCell::updateGeometry() : reading "
-                << nl
+                << isoField_ << " from time " <<fvm.time().timeName()
-                << "    regularise     : " << regularise_ << nl
+                << endl;
-                << "    average        : " << average_ << nl
-                << "    isoField       : " << isoField_ << nl
-                << "    isoValue       : " << isoVal_ << nl
-                << "    points         : " << points().size() << nl
-                << "    tris           : " << triSurface::size() << nl
-                << "    cut cells      : " << meshCells_.size() << endl;
         }
+
+        readFieldPtr_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    isoField_,
+                    fvm.time().timeName(),
+                    fvm,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                fvm
+            )
+        );
+
+        cellFldPtr = readFieldPtr_.operator->();
     }
+    const volScalarField& cellFld = *cellFldPtr;
+
+    tmp<pointScalarField> pointFld
+    (
+        volPointInterpolation::New(fvm).interpolate(cellFld)
+    );
+
+    if (average_)
+    {
+        //- From point field and interpolated cell.
+        scalarField cellAvg(fvm.nCells(), scalar(0.0));
+        labelField nPointCells(fvm.nCells(), 0);
+        {
+            for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
+            {
+                const labelList& pCells = fvm.pointCells(pointI);
+
+                forAll(pCells, i)
+                {
+                    label cellI = pCells[i];
+
+                    cellAvg[cellI] += pointFld().internalField()[pointI];
+                    nPointCells[cellI]++;
+                }
+            }
+        }
+        forAll(cellAvg, cellI)
+        {
+            cellAvg[cellI] /= nPointCells[cellI];
+        }
+
+        const isoSurfaceCell iso
+        (
+            fvm,
+            cellAvg,
+            pointFld().internalField(),
+            isoVal_,
+            regularise_
+        );
+
+        const_cast<sampledIsoSurfaceCell&>
+        (
+            *this
+        ).triSurface::operator=(iso);
+        meshCells_ = iso.meshCells();
+    }
+    else
+    {
+        //- Direct from cell field and point field. Gives bad continuity.
+        const isoSurfaceCell iso
+        (
+            fvm,
+            cellFld.internalField(),
+            pointFld().internalField(),
+            isoVal_,
+            regularise_
+        );
+
+        const_cast<sampledIsoSurfaceCell&>
+        (
+            *this
+        ).triSurface::operator=(iso);
+        meshCells_ = iso.meshCells();
+    }
+
+
+    if (debug)
+    {
+        Pout<< "sampledIsoSurfaceCell::updateGeometry() : constructed iso:"
+            << nl
+            << "    regularise     : " << regularise_ << nl
+            << "    average        : " << average_ << nl
+            << "    isoField       : " << isoField_ << nl
+            << "    isoValue       : " << isoVal_ << nl
+            << "    points         : " << points().size() << nl
+            << "    tris           : " << triSurface::size() << nl
+            << "    cut cells      : " << meshCells_.size() << endl;
+    }
+
+    return true;
 }
 
 
@@ -197,18 +202,17 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
     average_(dict.lookupOrDefault("average", true)),
     zoneName_(word::null),
     facesPtr_(NULL),
-    storedTimeIndex_(-1),
+    prevTimeIndex_(-1),
     meshCells_(0)
 {
-//    label zoneId = -1;
+//    dict.readIfPresent("zone", zoneName_);
-//    if (dict.readIfPresent("zone", zoneName_))
+//
+//    if (debug && zoneName_.size())
 //    {
-//        zoneId = mesh.cellZones().findZoneID(zoneName_);
+//        if (mesh.cellZones().findZoneID(zoneName_) < 0)
-//        if (debug && zoneId < 0)
 //        {
 //            Info<< "cellZone \"" << zoneName_
-//                << "\" not found - using entire mesh"
+//                << "\" not found - using entire mesh" << endl;
-//                << endl;
 //        }
 //    }
 }
@@ -222,13 +226,33 @@ Foam::sampledIsoSurfaceCell::~sampledIsoSurfaceCell()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::sampledIsoSurfaceCell::correct(const bool meshChanged)
+bool Foam::sampledIsoSurfaceCell::needsUpdate() const
 {
-    // Only change of mesh changes plane - zone restriction gets lost
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
-    if (meshChanged)
+
+    return fvm.time().timeIndex() != prevTimeIndex_;
+}
+
+
+bool Foam::sampledIsoSurfaceCell::expire()
+{
+    facesPtr_.clear();
+
+    // already marked as expired
+    if (prevTimeIndex_ == -1)
     {
-        facesPtr_.clear();
+        return false;
     }
+    
+    // force update
+    prevTimeIndex_ = -1;
+    return true;
+}
+
+
+bool Foam::sampledIsoSurfaceCell::update()
+{
+    return updateGeometry();
 }
 
 

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H

@@ -47,7 +47,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                       Class sampledIsoSurfaceCell Declaration
+                    Class sampledIsoSurfaceCell Declaration
 \*---------------------------------------------------------------------------*/
 
 class sampledIsoSurfaceCell
@@ -78,8 +78,8 @@ class sampledIsoSurfaceCell
 
         // Recreated for every isoSurface
 
-            //- Time at last call
+            //- Time at last call, also track it surface needs an update
-            mutable label storedTimeIndex_;
+            mutable label prevTimeIndex_;
 
             //- For every triangle the original cell in mesh
             mutable labelList meshCells_;
@@ -88,7 +88,8 @@ class sampledIsoSurfaceCell
     // Private Member Functions
 
         //- Create iso surface (if time has changed)
-        void createGeometry() const;
+        //  Do nothing (and return false) if no update was needed
+        bool updateGeometry() const;
 
         //- sample field on faces
         template <class Type>
@@ -127,6 +128,19 @@ public:
 
     // Member Functions
 
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surface as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if surface was already marked as expired.
+        virtual bool expire();
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was needed
+        virtual bool update();
+
+
         //- Points of surface
         virtual const pointField& points() const
         {
@@ -151,10 +165,6 @@ public:
         }
 
 
-        //- Correct for mesh movement and/or field changes
-        virtual void correct(const bool meshChanged);
-
-
         //- sample field on surface
         virtual tmp<scalarField> sample
         (

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C

@@ -40,7 +40,7 @@ Foam::sampledIsoSurfaceCell::sampleField
 ) const
 {
     // Recreate geometry if time has changed
-    createGeometry();
+    updateGeometry();
 
     return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
 }
@@ -54,7 +54,7 @@ Foam::sampledIsoSurfaceCell::interpolateField
 ) const
 {
     // Recreate geometry if time has changed
-    createGeometry();
+    updateGeometry();
 
     // One value per point
     tmp<Field<Type> > tvalues(new Field<Type>(points().size()));

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C

@@ -39,7 +39,8 @@ Foam::sampledIsoSurface::sampleField
 ) const
 {
     // Recreate geometry if time has changed
-    createGeometry();
+    updateGeometry();
+
     return tmp<Field<Type> >(new Field<Type>(vField, surface().meshCells()));
 }
 
@@ -67,8 +68,9 @@ Foam::sampledIsoSurface::interpolateField
 
     // Get pointers to sampling field (both original and interpolated one)
     getIsoFields();
+
     // Recreate geometry if time has changed
-    createGeometry();
+    updateGeometry();
 
     // Sample.
     return surface().interpolate

src/sampling/sampledSurface/patch/sampledPatch.C

@@ -40,17 +40,83 @@ namespace Foam
     addNamedToRunTimeSelectionTable(sampledSurface, sampledPatch, word, patch);
 }
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::sampledPatch::createGeometry()
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledPatch::sampledPatch
+(
+    const word& name,
+    const polyMesh& mesh,
+    const word& patchName,
+    const bool triangulate
+)
+:
+    sampledSurface(name, mesh),
+    patchName_(patchName),
+    triangulate_(triangulate),
+    needsUpdate_(true),
+    patchFaceLabels_(0)
+{}
+
+
+Foam::sampledPatch::sampledPatch
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    patchName_(dict.lookup("patchName")),
+    triangulate_(dict.lookupOrDefault("triangulate", false)),
+    needsUpdate_(true),
+    patchFaceLabels_(0)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledPatch::~sampledPatch()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::sampledPatch::needsUpdate() const
 {
+    return needsUpdate_;
+}
+
+
+bool Foam::sampledPatch::expire()
+{
+    // already marked as expired
+    if (needsUpdate_)
+    {
+        return false;
+    }
+
     sampledSurface::clearGeom();
     MeshStorage::clear();
     patchFaceLabels_.clear();
 
-    if (patchIndex_ != -1)
+    needsUpdate_ = true;
+    return true;
+}
+
+
+bool Foam::sampledPatch::update()
+{
+    if (!needsUpdate_)
     {
-        const polyPatch& p = mesh().boundaryMesh()[patchIndex_];
+        return false;
+    }
+
+    label patchI = mesh().boundaryMesh().findPatchID(patchName_);
+
+    if (patchI != -1)
+    {
+        const polyPatch& p = mesh().boundaryMesh()[patchI];
         this->storedPoints() = p.localPoints();
         this->storedFaces()  = p.localFaces();
 
@@ -76,60 +142,9 @@ void Foam::sampledPatch::createGeometry()
         print(Pout);
         Pout << endl;
     }
-}
 
-
+    needsUpdate_ = false;
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+    return true;
-
-Foam::sampledPatch::sampledPatch
-(
-    const word& name,
-    const polyMesh& mesh,
-    const word& patchName,
-    const bool triangulate
-)
-:
-    sampledSurface(name, mesh),
-    patchName_(patchName),
-    patchIndex_(mesh.boundaryMesh().findPatchID(patchName_)),
-    triangulate_(triangulate),
-    patchFaceLabels_(0)
-{
-    createGeometry();
-}
-
-
-Foam::sampledPatch::sampledPatch
-(
-    const word& name,
-    const polyMesh& mesh,
-    const dictionary& dict
-)
-:
-    sampledSurface(name, mesh, dict),
-    patchName_(dict.lookup("patchName")),
-    patchIndex_(mesh.boundaryMesh().findPatchID(patchName_)),
-    triangulate_(dict.lookupOrDefault("triangulate", false)),
-    patchFaceLabels_(0)
-{
-    createGeometry();
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledPatch::~sampledPatch()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::sampledPatch::correct(const bool meshChanged)
-{
-    if (meshChanged)
-    {
-        createGeometry();
-    }
 }
 
 
@@ -143,13 +158,6 @@ void Foam::sampledPatch::remapFaces
     if (&faceMap && faceMap.size())
     {
         MeshStorage::remapFaces(faceMap);
-//
-//        List<label> newCutCells(faceMap.size());
-//        forAll(faceMap, faceI)
-//        {
-//            newCutCells[faceI] = cutCells_[faceMap[faceI]];
-//        }
-//        cutCells_.transfer(newCutCells);
     }
 }
 

src/sampling/sampledSurface/patch/sampledPatch.H

@@ -50,30 +50,28 @@ namespace Foam
 
 class sampledPatch
 :
-    public PrimitiveMeshedSurface<face>,
+    public BasicMeshedSurface<face>,
     public sampledSurface
 {
     //- Private typedefs for convenience
-        typedef PrimitiveMeshedSurface<face> MeshStorage;
+        typedef BasicMeshedSurface<face> MeshStorage;
 
     // Private data
 
         //- Name of patch
         const word patchName_;
 
-        //- Index of patch in boundaryMesh
-        const label patchIndex_;
-
         //- Triangulated faces or keep faces as is
         bool triangulate_;
 
+        //- Track if the surface needs an update
+        mutable bool needsUpdate_;
+
         //- Local patch face labels
         labelList patchFaceLabels_;
 
     // Private Member Functions
 
-        //- Do all to construct geometry.
-        void createGeometry();
 
         //- sample field on faces
         template <class Type>
@@ -88,7 +86,7 @@ class sampledPatch
         interpolateField(const interpolation<Type>&) const;
 
 
-        //- remap action on triangulation or cleanup 
+        //- remap action on triangulation or cleanup
         virtual void remapFaces(const UList<label>& faceMap);
 
 public:
@@ -124,6 +122,19 @@ public:
 
     // Member Functions
 
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surface as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if surface was already marked as expired.
+        virtual bool expire();
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was needed
+        virtual bool update();
+
+
         const word patchName() const
         {
             return patchName_;
@@ -131,7 +142,7 @@ public:
 
         label patchIndex() const
         {
-            return patchIndex_;
+            return mesh().boundaryMesh().findPatchID(patchName_);
         }
 
         const labelList& patchFaceLabels() const
@@ -152,10 +163,6 @@ public:
         }
 
 
-        //- Correct for mesh movement and/or field changes
-        virtual void correct(const bool meshChanged);
-
-
         //- sample field on surface
         virtual tmp<scalarField> sample
         (

src/sampling/sampledSurface/plane/sampledPlane.C

@@ -39,33 +39,6 @@ namespace Foam
     addNamedToRunTimeSelectionTable(sampledSurface, sampledPlane, word, plane);
 }
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::sampledPlane::createGeometry
-(
-    const polyMesh& mesh,
-    const label zoneId
-)
-{
-    sampledSurface::clearGeom();
-
-    if (zoneId < 0)
-    {
-        reCut(mesh);
-    }
-    else
-    {
-        reCut(mesh, mesh.cellZones()[zoneId]);
-    }
-
-    if (debug)
-    {
-        print(Pout);
-        Pout << endl;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::sampledPlane::sampledPlane
@@ -78,21 +51,17 @@ Foam::sampledPlane::sampledPlane
 :
     sampledSurface(name, mesh),
     cuttingPlane(planeDesc),
-    zoneName_(zoneName)
+    zoneName_(zoneName),
+    needsUpdate_(true)
 {
-    label zoneId = -1;
+    if (debug && zoneName_.size())
-    if (zoneName_.size())
     {
-        zoneId = mesh.cellZones().findZoneID(zoneName_);
+        if (mesh.cellZones().findZoneID(zoneName_) < 0)
-        if (debug && zoneId < 0)
         {
             Info<< "cellZone \"" << zoneName_
-                << "\" not found - using entire mesh"
+                << "\" not found - using entire mesh" << endl;
-                << endl;
         }
     }
-
-    createGeometry(mesh, zoneId);
 }
 
 
@@ -105,9 +74,9 @@ Foam::sampledPlane::sampledPlane
 :
     sampledSurface(name, mesh, dict),
     cuttingPlane(plane(dict.lookup("basePoint"), dict.lookup("normalVector"))),
-    zoneName_(word::null)
+    zoneName_(word::null),
+    needsUpdate_(true)
 {
-
     // make plane relative to the coordinateSystem (Cartesian)
     // allow lookup from global coordinate systems
     if (dict.found("coordinateSystem"))
@@ -121,21 +90,17 @@ Foam::sampledPlane::sampledPlane
         static_cast<plane&>(*this) = plane(base, norm);
     }
 
+    dict.readIfPresent("zone", zoneName_);
 
-    label zoneId = -1;
+    if (debug && zoneName_.size())
-    if (dict.readIfPresent("zone", zoneName_))
     {
-        zoneId = mesh.cellZones().findZoneID(zoneName_);
+        if (mesh.cellZones().findZoneID(zoneName_) < 0)
-        if (debug && zoneId < 0)
         {
             Info<< "cellZone \"" << zoneName_
-                << "\" not found - using entire mesh"
+                << "\" not found - using entire mesh" << endl;
-                << endl;
         }
     }
 
-
-    createGeometry(mesh, zoneId);
 }
 
 
@@ -147,19 +112,59 @@ Foam::sampledPlane::~sampledPlane()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::sampledPlane::correct(const bool meshChanged)
+bool Foam::sampledPlane::needsUpdate() const
 {
-    // Only change of mesh changes plane - zone restriction gets lost
+    return needsUpdate_;
-    if (meshChanged)
+}
-    {
-        label zoneId = -1;
-        if (zoneName_.size())
-        {
-            zoneId = mesh().cellZones().findZoneID(zoneName_);
-        }
 
-        createGeometry(mesh(), zoneId);
+
+bool Foam::sampledPlane::expire()
+{
+    // already marked as expired
+    if (needsUpdate_)
+    {
+        return false;
     }
+
+    sampledSurface::clearGeom();
+
+    needsUpdate_ = true;
+    return true;
+}
+
+
+bool Foam::sampledPlane::update()
+{
+    if (!needsUpdate_)
+    {
+        return false;
+    }
+
+    sampledSurface::clearGeom();
+
+    label zoneId = -1;
+    if (zoneName_.size())
+    {
+        zoneId = mesh().cellZones().findZoneID(zoneName_);
+    }
+
+    if (zoneId < 0)
+    {
+        reCut(mesh());
+    }
+    else
+    {
+        reCut(mesh(), mesh().cellZones()[zoneId]);
+    }
+
+    if (debug)
+    {
+        print(Pout);
+        Pout << endl;
+    }
+
+    needsUpdate_ = false;
+    return true;
 }
 
 
@@ -182,6 +187,7 @@ Foam::sampledPlane::sample
     return sampleField(vField);
 }
 
+
 Foam::tmp<Foam::sphericalTensorField>
 Foam::sampledPlane::sample
 (

src/sampling/sampledSurface/plane/sampledPlane.H

@@ -58,14 +58,10 @@ class sampledPlane
         //- zone name (if restricted to zones)
         word zoneName_;
 
-    // Private Member Functions
+        //- Track if the surface needs an update
+        mutable bool needsUpdate_;
 
-        //- Do all to create geometry.
+    // Private Member Functions
-        void createGeometry
-        (
-            const polyMesh& mesh,
-            const label zoneId
-        );
 
         //- sample field on faces
         template <class Type>
@@ -113,6 +109,19 @@ public:
 
     // Member Functions
 
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surface as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if surface was already marked as expired.
+        virtual bool expire();
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was needed
+        virtual bool update();
+
+
         //- Points of surface
         virtual const pointField& points() const
         {
@@ -131,10 +140,6 @@ public:
             return cuttingPlane::cutCells();
         }
 
-        //- Correct for mesh movement and/or field changes
-        virtual void correct(const bool meshChanged);
-
-
         //- sample field on surface
         virtual tmp<scalarField> sample
         (

src/sampling/sampledSurface/sampledSurface/sampledSurface.H

@@ -28,6 +28,19 @@ Class
 Description
     An abstract class for surfaces with sampling.
 
+    The constructors for the derived classes should generally start in a
+    'expired' condition (ie, needsUpdate() == true) and rely on a
+    subsequent call to the update() method to complete the initialization.
+    Delaying the final construction as late as possible allows the
+    construction of surfaces that may depend on intermediate calculation
+    results (eg, iso-surfaces) and also avoids the unnecessary
+    reconstruction of surfaces between sampling intervals.
+
+    It is the responsibility of the caller to ensure that the surface
+    update() is called before the surface is used.  The update() method
+    implementation should do nothing when the surface is already
+    up-to-date.
+
 SourceFiles
     sampledSurface.C
     sampledSurfaceTemplates.C
@@ -237,15 +250,25 @@ public:
             return interpolate_;
         }
 
+        //- Does the surface need an update?
+        virtual bool needsUpdate() const = 0;
+
+        //- Mark the surface as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if surface was already marked as expired.
+        virtual bool expire() = 0;
+
+        //- Update the surface as required.
+        //  Do nothing (and return false) if no update was required
+        virtual bool update() = 0;
+
+
         //- Points of surface
         virtual const pointField& points() const = 0;
 
         //- Faces of surface
         virtual const faceList& faces() const = 0;
 
-        //- Correct for mesh movement and/or field changes
-        virtual void correct(const bool meshChanged) = 0;
-
         //- Return face area vectors
         virtual const vectorField& Sf() const;
 

src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C

@@ -60,11 +60,13 @@ namespace Foam
         }
     };
 
+
     defineTypeNameAndDebug(sampledSurfaces, 0);
 }
 
-bool Foam::sampledSurfaces::verbose_ = false;
 
+bool Foam::sampledSurfaces::verbose_(false);
+Foam::scalar Foam::sampledSurfaces::mergeTol_(1e-10);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -155,18 +157,213 @@ bool Foam::sampledSurfaces::checkFieldTypes()
 }
 
 
-void Foam::sampledSurfaces::mergeSurfaces()
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledSurfaces::sampledSurfaces
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    PtrList<sampledSurface>(),
+    name_(name),
+    mesh_(refCast<const fvMesh>(obr)),
+    loadFromFiles_(loadFromFiles),
+    outputPath_(fileName::null),
+    fieldNames_(),
+    interpolationScheme_(word::null),
+    writeFormat_(word::null),
+    mergeList_(),
+    scalarFields_(),
+    vectorFields_(),
+    sphericalTensorFields_(),
+    symmTensorFields_(),
+    tensorFields_()
 {
-    if (!Pstream::parRun())
+    if (Pstream::parRun())
     {
-        return;
+        outputPath_ = mesh_.time().path()/".."/name_;
+    }
+    else
+    {
+        outputPath_ = mesh_.time().path()/name_;
     }
 
-    // Merge close points (1E-10 of mesh bounding box)
+    read(dict);
-    const scalar mergeTol = 1e-10;
+}
 
-    const boundBox& bb = mesh_.globalData().bb();
+
-    scalar mergeDim = mergeTol * mag(bb.max() - bb.min());
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledSurfaces::~sampledSurfaces()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+void Foam::sampledSurfaces::verbose(const bool verbosity)
+{
+    verbose_ = verbosity;
+}
+
+
+void Foam::sampledSurfaces::execute()
+{
+    // Do nothing - only valid on write
+}
+
+
+void Foam::sampledSurfaces::write()
+{
+    if (size() && checkFieldTypes())
+    {
+        // finalize surfaces, merge points etc.
+        update();
+
+        sampleAndWrite(scalarFields_);
+        sampleAndWrite(vectorFields_);
+        sampleAndWrite(sphericalTensorFields_);
+        sampleAndWrite(symmTensorFields_);
+        sampleAndWrite(tensorFields_);
+    }
+}
+
+
+void Foam::sampledSurfaces::read(const dictionary& dict)
+{
+    fieldNames_ = wordList(dict.lookup("fields"));
+
+    interpolationScheme_ = dict.lookupOrDefault<word>
+    (
+        "interpolationScheme",
+        "cell"
+    );
+    writeFormat_ = dict.lookupOrDefault<word>
+    (
+        "surfaceFormat",
+        "null"
+    );
+
+    PtrList<sampledSurface> newList
+    (
+        dict.lookup("surfaces"),
+        sampledSurface::iNew(mesh_)
+    );
+
+    transfer(newList);
+
+    if (Pstream::parRun())
+    {
+        mergeList_.setSize(size());
+    }
+
+    // ensure all surfaces and merge information are expired
+    expire();
+
+    if (Pstream::master() && debug)
+    {
+        Pout<< "sample fields:" << fieldNames_ << nl
+            << "sample surfaces:" << nl << "(" << nl;
+
+        forAll(*this, surfI)
+        {
+            Pout<< "  " << operator[](surfI) << endl;
+        }
+        Pout<< ")" << endl;
+    }
+}
+
+
+void Foam::sampledSurfaces::updateMesh(const mapPolyMesh&)
+{
+    expire();
+}
+
+
+void Foam::sampledSurfaces::movePoints(const pointField&)
+{
+    expire();
+}
+
+
+void Foam::sampledSurfaces::readUpdate(const polyMesh::readUpdateState state)
+{
+    if (state != polyMesh::UNCHANGED)
+    {
+        expire();
+    }
+}
+
+
+bool Foam::sampledSurfaces::needsUpdate() const
+{
+    forAll(*this, surfI)
+    {
+        if (operator[](surfI).needsUpdate())
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+bool Foam::sampledSurfaces::expire()
+{
+    bool justExpired = false;
+
+    forAll(*this, surfI)
+    {
+        if (operator[](surfI).expire())
+        {
+            justExpired = true;
+        }
+
+        // clear merge information
+        if (Pstream::parRun())
+        {
+            mergeList_[surfI].clear();
+        }
+    }
+
+    // reset interpolation
+    pointMesh::Delete(mesh_);
+    volPointInterpolation::Delete(mesh_);
+
+    // true if any surfaces just expired
+    return justExpired;
+}
+
+
+bool Foam::sampledSurfaces::update()
+{
+    bool updated = false;
+
+    if (!needsUpdate())
+    {
+        return updated;
+    }
+
+    // serial: quick and easy, no merging required
+    if (!Pstream::parRun())
+    {
+        forAll(*this, surfI)
+        {
+            if (operator[](surfI).update())
+            {
+                updated = true;
+            }
+        }
+
+        return updated;
+    }
+
+    // dimension as fraction of mesh bounding box
+    scalar mergeDim = mergeTol_ * mesh_.globalData().bb().mag();
 
     if (Pstream::master() && debug)
     {
@@ -174,11 +371,20 @@ void Foam::sampledSurfaces::mergeSurfaces()
             << mergeDim << " meter" << endl;
     }
 
-    mergeList_.setSize(size());
     forAll(*this, surfI)
     {
         sampledSurface& s = operator[](surfI);
 
+        if (s.update())
+        {
+            updated = true;
+        }
+        else
+        {
+            continue;
+        }
+
+
         // Collect points from all processors
         List<pointField> gatheredPoints(Pstream::nProcs());
         gatheredPoints[Pstream::myProcNo()] = s.points();
@@ -253,147 +459,8 @@ void Foam::sampledSurfaces::mergeSurfaces()
             }
         }
     }
-}
 
-
+    return updated;
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::sampledSurfaces::sampledSurfaces
-(
-    const word& name,
-    const objectRegistry& obr,
-    const dictionary& dict,
-    const bool loadFromFiles
-)
-:
-    PtrList<sampledSurface>(),
-    name_(name),
-    mesh_(refCast<const fvMesh>(obr)),
-    loadFromFiles_(loadFromFiles),
-    outputPath_(fileName::null),
-    fieldNames_(),
-    interpolationScheme_(word::null),
-    writeFormat_(word::null),
-    mergeList_(),
-    scalarFields_(),
-    vectorFields_(),
-    sphericalTensorFields_(),
-    symmTensorFields_(),
-    tensorFields_()
-{
-    if (Pstream::parRun())
-    {
-        outputPath_ = mesh_.time().path()/".."/name_;
-    }
-    else
-    {
-        outputPath_ = mesh_.time().path()/name_;
-    }
-
-    read(dict);
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::sampledSurfaces::~sampledSurfaces()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::sampledSurfaces::verbose(const bool verbosity)
-{
-    verbose_ = verbosity;
-}
-
-
-void Foam::sampledSurfaces::execute()
-{
-    // Do nothing - only valid on write
-}
-
-
-void Foam::sampledSurfaces::write()
-{
-    if (size() && checkFieldTypes())
-    {
-        sampleAndWrite(scalarFields_);
-        sampleAndWrite(vectorFields_);
-        sampleAndWrite(sphericalTensorFields_);
-        sampleAndWrite(symmTensorFields_);
-        sampleAndWrite(tensorFields_);
-    }
-}
-
-
-void Foam::sampledSurfaces::read(const dictionary& dict)
-{
-    fieldNames_ = wordList(dict.lookup("fields"));
-
-    interpolationScheme_ = "cell";
-    dict.readIfPresent("interpolationScheme", interpolationScheme_);
-
-    writeFormat_ = "null";
-    dict.readIfPresent("surfaceFormat", writeFormat_);
-
-
-    PtrList<sampledSurface> newList
-    (
-        dict.lookup("surfaces"),
-        sampledSurface::iNew(mesh_)
-    );
-
-    transfer(newList);
-    mergeSurfaces();
-
-    if (Pstream::master() && debug)
-    {
-        Pout<< "sample fields:" << fieldNames_ << nl
-            << "sample surfaces:" << nl << "(" << nl;
-
-        forAll(*this, surfI)
-        {
-            Pout << "  " << operator[](surfI) << endl;
-        }
-        Pout << ")" << endl;
-    }
-}
-
-
-void Foam::sampledSurfaces::correct()
-{
-    forAll(*this, surfI)
-    {
-        operator[](surfI).correct(true);
-    }
-
-    // reset interpolation
-    pointMesh::Delete(mesh_);
-    volPointInterpolation::Delete(mesh_);
-
-    mergeSurfaces();
-}
-
-
-void Foam::sampledSurfaces::updateMesh(const mapPolyMesh&)
-{
-    correct();
-}
-
-
-void Foam::sampledSurfaces::movePoints(const pointField&)
-{
-    correct();
-}
-
-
-void Foam::sampledSurfaces::readUpdate(const polyMesh::readUpdateState state)
-{
-    if (state != polyMesh::UNCHANGED)
-    {
-        correct();
-    }
 }
 
 

src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.H

@@ -116,6 +116,14 @@ class sampledSurfaces
             pointField points;
             faceList   faces;
             labelList  pointsMap;
+
+            //- Clear all storage
+            void clear()
+            {
+                points.clear();
+                faces.clear();
+                pointsMap.clear();
+            }
         };
 
 
@@ -124,6 +132,8 @@ class sampledSurfaces
         //- output verbosity
         static bool verbose_;
 
+        //- Tolerance for merging points (fraction of mesh bounding box)
+        static scalar mergeTol_;
 
     // Private data
 
@@ -174,9 +184,6 @@ class sampledSurfaces
         //- Classify field types, return true if nFields > 0
         bool checkFieldTypes();
 
-        //- Merge points on surfaces
-        void mergeSurfaces();
-
         //- Find the fields in the list of the given type, return count
         template<class Type>
         label grep
@@ -228,7 +235,20 @@ public:
 
     // Member Functions
 
-        //- Return name of the set of probes
+        //- Does any of the surfaces need an update?
+        virtual bool needsUpdate() const;
+
+        //- Mark the surfaces as needing an update.
+        //  May also free up unneeded data.
+        //  Return false if all surfaces were already marked as expired.
+        virtual bool expire();
+
+        //- Update the surfaces as required and merge surface points (parallel).
+        //  Return false if no surfaces required an update.
+        virtual bool update();
+
+
+        //- Return name of the set of surfaces
         virtual const word& name() const
         {
             return name_;
@@ -243,20 +263,18 @@ public:
         //- Sample and write
         virtual void write();
 
-        //- Read the sampledSurfaces
+        //- Read the sampledSurfaces dictionary
         virtual void read(const dictionary&);
 
-        //- Correct for mesh changes
+        //- Update for changes of mesh - expires the surfaces
-        void correct();
-
-        //- Update for changes of mesh
         virtual void updateMesh(const mapPolyMesh&);
 
-        //- Update for mesh point-motion
+        //- Update for mesh point-motion - expires the surfaces
         virtual void movePoints(const pointField&);
 
-        //- Update for changes of mesh due to readUpdate
+        //- Update for changes of mesh due to readUpdate - expires the surfaces
         virtual void readUpdate(const polyMesh::readUpdateState state);
+
 };
 
 

src/sampling/sampledSurface/sampledSurfaces/sampledSurfacesTemplates.C

@@ -161,7 +161,7 @@ void Foam::sampledSurfaces::sampleAndWrite
 }
 
 
-template <class Type>
+template<class Type>
 void Foam::sampledSurfaces::sampleAndWrite
 (
     fieldGroup<Type>& fields