rmt isosurface correction

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -155,21 +155,25 @@ surfaces
         // Optional: whether to leave as faces (=default) or triangulate
     }
 
-    constantIso
-    {
-        name            isoSurface;
-        isoField        rho;
-        isoValue        0.5;
-        interpolate     false
-    }
     interpolatedIso
     {
-        name            isoSurface;
+        // Iso surface for interpolated values only
+        type            isoSurface;
         isoField        rho;
         isoValue        0.5;
         interpolate     true;
         //regularise      false;    //optional: do not simplify
     }
+    constantIso
+    {
+        // Iso surface for constant values. Guarantees triangles to not
+        // cross cells.
+        type            isoSurfaceCell;
+        isoField        rho;
+        isoValue        0.5;
+        interpolate     false;
+        //regularise      false;    //optional: do not simplify
+    }
 );
 
 

src/sampling/Make/files

@@ -24,6 +24,8 @@ sampledSurface/patch/sampledPatch.C
 sampledSurface/plane/sampledPlane.C
 sampledSurface/isoSurface/isoSurface.C
 sampledSurface/isoSurface/sampledIsoSurface.C
+sampledSurface/isoSurface/isoSurfaceCell.C
+sampledSurface/isoSurface/sampledIsoSurfaceCell.C
 sampledSurface/distanceSurface/distanceSurface.C
 sampledSurface/sampledSurface/sampledSurface.C
 sampledSurface/sampledSurfaces/sampledSurfaces.C

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -30,6 +30,7 @@ License
 #include "volPointInterpolation.H"
 #include "addToRunTimeSelectionTable.H"
 #include "fvMesh.H"
+#include "isoSurfaceCell.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -110,7 +111,7 @@ void Foam::distanceSurface::createGeometry() const
     }
 
     //- Direct from cell field and point field.
-    const isoSurface iso
+    const isoSurfaceCell iso
     (
         mesh(),
         cellDistance,

src/sampling/sampledSurface/isoSurface/isoSurface.H

@@ -27,9 +27,7 @@ Class
 
 Description
     A surface formed by the iso value.
-    After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke
+    After "Regularised Marching Tetrahedra: improved iso-surface extraction",
-    and
-    "Regularised Marching Tetrahedra: improved iso-surface extraction",
     G.M. Treece, R.W. Prager and A.H. Gee.
 
     Note:
@@ -39,7 +37,9 @@ Description
     - does tets without using cell centres/cell values. Not tested.
     - regularisation can create duplicate triangles/non-manifold surfaces.
     Current handling of those is bit ad-hoc for now and not perfect.
-    Uses geometric merge.
+    - uses geometric merge with fraction of bounding box as distance.
+    - triangles can be between two cell centres so constant sampling
+      does not make sense.
     - does not do 2D correctly, creates non-flat iso surface.
 
 SourceFiles
@@ -54,13 +54,14 @@ SourceFiles
 #include "labelPair.H"
 #include "pointIndexHit.H"
 #include "PackedList.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class polyMesh;
+class fvMesh;
 
 /*---------------------------------------------------------------------------*\
                        Class isoSurface Declaration
@@ -88,7 +89,7 @@ class isoSurface
 
 
         //- Reference to mesh
-        const polyMesh& mesh_;
+        const fvMesh& mesh_;
 
         //- Isosurface value
         const scalar iso_;
@@ -97,6 +98,15 @@ class isoSurface
         const scalar mergeDistance_;
 
 
+        //- Whether face might be cut
+        List<cellCutType> faceCutType_;
+
+        //- Whether cell might be cut
+        List<cellCutType> cellCutType_;
+
+        //- Estimated number of cut cells
+        label nCutCells_;
+
         //- For every triangle the original cell in mesh
         labelList meshCells_;
 
@@ -113,16 +123,12 @@ class isoSurface
             const scalar s1
         ) const;
 
-        //List<triFace> triangulate(const face& f) const;
+        //- Set faceCutType,cellCutType.
-
+        void calcCutTypes
-        //- Determine whether cell is cut
-        cellCutType calcCutType
         (
-            const PackedList<1>&,
+            const volScalarField& cVals,
-            const scalarField& cellValues,
+            const scalarField& pVals
-            const scalarField& pointValues,
+        );
-            const label
-        ) const;
 
         static labelPair findCommonPoints
         (
@@ -138,79 +144,115 @@ class isoSurface
             DynamicList<labelledTri, 64>& localTris
         );
 
+        void getNeighbour
+        (
+            const labelList& boundaryRegion,
+            const volScalarField& cVals,
+            const label cellI,
+            const label faceI,
+            scalar& nbrValue,
+            point& nbrPoint
+        ) const;
+
         //- Determine per cc whether all near cuts can be snapped to single
         //  point.
         void calcSnappedCc
         (
-            const PackedList<1>& isTet,
+            const labelList& boundaryRegion,
-            const scalarField& cVals,
+            const volScalarField& cVals,
             const scalarField& pVals,
-            const List<cellCutType>& cellCutType,
+            DynamicList<point>& snappedPoints,
-            DynamicList<point>& triPoints,
             labelList& snappedCc
         ) const;
 
-        //- Generate triangles for face connected to pointI
-        void genPointTris
-        (
-            const scalarField& cellValues,
-            const scalarField& pointValues,
-            const label pointI,
-            const face& f,
-            const label cellI,
-            DynamicList<point, 64>& localTriPoints
-        ) const;
-
-        //- Generate triangles for tet connected to pointI
-        void genPointTris
-        (
-            const scalarField& pointValues,
-            const label pointI,
-            const label cellI,
-            DynamicList<point, 64>& localTriPoints
-        ) const;
-
         //- Determine per point whether all near cuts can be snapped to single
         //  point.
         void calcSnappedPoint
         (
             const PackedList<1>& isBoundaryPoint,
-            const PackedList<1>& isTet,
+            const labelList& boundaryRegion,
-            const scalarField& cVals,
+            const volScalarField& cVals,
             const scalarField& pVals,
-            const List<cellCutType>& cellCutType,
+            DynamicList<point>& snappedPoints,
-            DynamicList<point>& triPoints,
             labelList& snappedPoint
         ) const;
 
         //- Generate single point by interpolation or snapping
-        point generatePoint
+        template<class Type>
+        Type generatePoint
         (
-            const DynamicList<point>& snappedPoints,
+            const DynamicList<Type>& snappedPoints,
+
             const scalar s0,
-            const point& p0,
+            const Type& p0,
             const label p0Index,
+
             const scalar s1,
-            const point& p1,
+            const Type& p1,
             const label p1Index
         ) const;
 
+        template<class Type>
         void generateTriPoints
         (
-            const DynamicList<point>& snapped,
+            const DynamicList<Type>& snapped,
+
             const scalar s0,
-            const point& p0,
+            const Type& p0,
             const label p0Index,
+
             const scalar s1,
-            const point& p1,
+            const Type& p1,
             const label p1Index,
+
             const scalar s2,
-            const point& p2,
+            const Type& p2,
             const label p2Index,
+
             const scalar s3,
-            const point& p3,
+            const Type& p3,
             const label p3Index,
-            DynamicList<point>& points
+
+            DynamicList<Type>& points
+        ) const;
+
+        template<class Type>
+        label generateTriPoints
+        (
+            const volScalarField& cVals,
+            const scalarField& pVals,
+
+            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+            const Field<Type>& pCoords,
+
+            const DynamicList<Type>& snappedPoints,
+            const labelList& snappedCc,
+            const labelList& snappedPoint,
+            const label faceI,
+
+            const scalar neiVal,
+            const Type& neiPt,
+            const label neiSnap,
+
+            DynamicList<Type>& triPoints,
+            DynamicList<label>& triMeshCells
+        ) const;
+
+        template<class Type>
+        void generateTriPoints
+        (
+            const volScalarField& cVals,
+            const scalarField& pVals,
+
+            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+            const Field<Type>& pCoords,
+
+            const DynamicList<Type>& snappedPoints,
+            const labelList& snappedCc,
+            const labelList& snappedPoint,
+
+            DynamicList<Type>& triPoints,
+            DynamicList<label>& triMeshCells
         ) const;
 
         triSurface stitchTriPoints
@@ -276,12 +318,10 @@ class isoSurface
         (
             const triSurface& s,
             const labelList& newToOldFaces,
+            labelList& oldToNewPoints,
             labelList& newToOldPoints
         );
 
-        //- Combine all triangles inside a cell into a minimal triangulation
-        void combineCellTriangles();
-
 public:
 
     //- Runtime type information
@@ -290,12 +330,13 @@ public:
 
     // Constructors
 
-        //- Construct from dictionary
+        //- Construct from cell values and point values. Uses boundaryField
+        //  for boundary values. Requires on coupled patchfields to be set
+        //  to the opposite cell value.
         isoSurface
         (
-            const polyMesh& mesh,
+            const volScalarField& cellIsoVals,
-            const scalarField& cellValues,
+            const scalarField& pointIsoVals,
-            const scalarField& pointValues,
             const scalar iso,
             const bool regularise,
             const scalar mergeTol = 1E-6    // fraction of bounding box
@@ -311,10 +352,22 @@ public:
         }
 
         //- For every unmerged triangle point the point in the triSurface
-        const labelList triPointMergeMap() const
+        const labelList& triPointMergeMap() const
         {
             return triPointMergeMap_;
         }
+
+        //- Interpolates cCoords,pCoords. Takes the original fields
+        //  used to create the iso surface.
+        template <class Type>
+        tmp<Field<Type> > interpolate
+        (
+            const volScalarField& cellIsoVals,
+            const scalarField& pointIsoVals,
+            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+            const Field<Type>& pCoords
+        ) const;
+
 };
 
 
@@ -324,6 +377,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+#   include "isoSurfaceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H

@@ -0,0 +1,325 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::isoSurfaceCell
+
+Description
+    A surface formed by the iso value.
+    After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke
+    and
+    "Regularised Marching Tetrahedra: improved iso-surface extraction",
+    G.M. Treece, R.W. Prager and A.H. Gee.
+
+    See isoSurface. This is a variant which does tetrahedrisation from
+    triangulation of face to cell centre instead of edge of face to two
+    neighbouring cell centres. This gives much lower quality triangles
+    but they are local to a cell.
+
+SourceFiles
+    isoSurfaceCell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoSurfaceCell_H
+#define isoSurfaceCell_H
+
+#include "triSurface.H"
+#include "labelPair.H"
+#include "pointIndexHit.H"
+#include "PackedList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                       Class isoSurfaceCell Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoSurfaceCell
+:
+    public triSurface
+{
+    // Private data
+
+        enum segmentCutType
+        {
+            NEARFIRST,      // intersection close to e.first()
+            NEARSECOND,     //  ,,                   e.second()
+            NOTNEAR         // no intersection
+        };
+
+        enum cellCutType
+        {
+            NOTCUT,     // not cut
+            SPHERE,     // all edges to cell centre cut
+            CUT         // normal cut
+        };
+
+
+        //- Reference to mesh
+        const polyMesh& mesh_;
+
+        //- isoSurfaceCell value
+        const scalar iso_;
+
+        //- When to merge points
+        const scalar mergeDistance_;
+
+
+        //- For every triangle the original cell in mesh
+        labelList meshCells_;
+
+        //- For every unmerged triangle point the point in the triSurface
+        labelList triPointMergeMap_;
+
+
+    // Private Member Functions
+
+        //- Get location of iso value as fraction inbetween s0,s1
+        scalar isoFraction
+        (
+            const scalar s0,
+            const scalar s1
+        ) const;
+
+        //List<triFace> triangulate(const face& f) const;
+
+        //- Determine whether cell is cut
+        cellCutType calcCutType
+        (
+            const PackedList<1>&,
+            const scalarField& cellValues,
+            const scalarField& pointValues,
+            const label
+        ) const;
+
+        static labelPair findCommonPoints
+        (
+            const labelledTri&,
+            const labelledTri&
+        );
+
+        static point calcCentre(const triSurface&);
+
+        static pointIndexHit collapseSurface
+        (
+            pointField& localPoints,
+            DynamicList<labelledTri, 64>& localTris
+        );
+
+        //- Determine per cc whether all near cuts can be snapped to single
+        //  point.
+        void calcSnappedCc
+        (
+            const PackedList<1>& isTet,
+            const scalarField& cVals,
+            const scalarField& pVals,
+            const List<cellCutType>& cellCutType,
+            DynamicList<point>& triPoints,
+            labelList& snappedCc
+        ) const;
+
+        //- Generate triangles for face connected to pointI
+        void genPointTris
+        (
+            const scalarField& cellValues,
+            const scalarField& pointValues,
+            const label pointI,
+            const face& f,
+            const label cellI,
+            DynamicList<point, 64>& localTriPoints
+        ) const;
+
+        //- Generate triangles for tet connected to pointI
+        void genPointTris
+        (
+            const scalarField& pointValues,
+            const label pointI,
+            const label cellI,
+            DynamicList<point, 64>& localTriPoints
+        ) const;
+
+        //- Determine per point whether all near cuts can be snapped to single
+        //  point.
+        void calcSnappedPoint
+        (
+            const PackedList<1>& isBoundaryPoint,
+            const PackedList<1>& isTet,
+            const scalarField& cVals,
+            const scalarField& pVals,
+            const List<cellCutType>& cellCutType,
+            DynamicList<point>& triPoints,
+            labelList& snappedPoint
+        ) const;
+
+        //- Generate single point by interpolation or snapping
+        point generatePoint
+        (
+            const DynamicList<point>& snappedPoints,
+            const scalar s0,
+            const point& p0,
+            const label p0Index,
+            const scalar s1,
+            const point& p1,
+            const label p1Index
+        ) const;
+
+        void generateTriPoints
+        (
+            const DynamicList<point>& snapped,
+            const scalar s0,
+            const point& p0,
+            const label p0Index,
+            const scalar s1,
+            const point& p1,
+            const label p1Index,
+            const scalar s2,
+            const point& p2,
+            const label p2Index,
+            const scalar s3,
+            const point& p3,
+            const label p3Index,
+            DynamicList<point>& points
+        ) const;
+
+        triSurface stitchTriPoints
+        (
+            const bool checkDuplicates,
+            const List<point>& triPoints,
+            labelList& triPointReverseMap,  // unmerged to merged point
+            labelList& triMap               // merged to unmerged triangle
+        ) const;
+
+        //- Check single triangle for (topological) validity
+        static bool validTri(const triSurface&, const label);
+
+        //- Determine edge-face addressing
+        void calcAddressing
+        (
+            const triSurface& surf,
+            List<FixedList<label, 3> >& faceEdges,
+            labelList& edgeFace0,
+            labelList& edgeFace1,
+            Map<labelList>& edgeFacesRest
+        ) const;
+
+        //- Determine orientation
+        static void walkOrientation
+        (
+            const triSurface& surf,
+            const List<FixedList<label, 3> >& faceEdges,
+            const labelList& edgeFace0,
+            const labelList& edgeFace1,
+            const label seedTriI,
+            labelList& flipState
+        );
+
+        //- Orient surface
+        static void orientSurface
+        (
+            triSurface&,
+            const List<FixedList<label, 3> >& faceEdges,
+            const labelList& edgeFace0,
+            const labelList& edgeFace1,
+            const Map<labelList>& edgeFacesRest
+        );
+
+        //- Is triangle (given by 3 edges) not fully connected?
+        static bool danglingTriangle
+        (
+            const FixedList<label, 3>& fEdges,
+            const labelList& edgeFace1
+        );
+
+        //- Mark all non-fully connected triangles
+        static label markDanglingTriangles
+        (
+            const List<FixedList<label, 3> >& faceEdges,
+            const labelList& edgeFace0,
+            const labelList& edgeFace1,
+            const Map<labelList>& edgeFacesRest,
+            boolList& keepTriangles
+        );
+
+        static triSurface subsetMesh
+        (
+            const triSurface& s,
+            const labelList& newToOldFaces,
+            labelList& oldToNewPoints,
+            labelList& newToOldPoints
+        );
+
+        //- Combine all triangles inside a cell into a minimal triangulation
+        void combineCellTriangles();
+
+public:
+
+    //- Runtime type information
+    TypeName("isoSurfaceCell");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        isoSurfaceCell
+        (
+            const polyMesh& mesh,
+            const scalarField& cellValues,
+            const scalarField& pointValues,
+            const scalar iso,
+            const bool regularise,
+            const scalar mergeTol = 1E-6    // fraction of bounding box
+        );
+
+
+    // Member Functions
+
+        //- For every face original cell in mesh
+        const labelList& meshCells() const
+        {
+            return meshCells_;
+        }
+
+        //- For every unmerged triangle point the point in the triSurface
+        const labelList triPointMergeMap() const
+        {
+            return triPointMergeMap_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C

@@ -0,0 +1,485 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "isoSurface.H"
+#include "polyMesh.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::isoSurface::generatePoint
+(
+    const DynamicList<Type>& snappedPoints,
+
+    const scalar s0,
+    const Type& p0,
+    const label p0Index,
+
+    const scalar s1,
+    const Type& p1,
+    const label p1Index
+) const
+{
+    scalar d = s1-s0;
+
+    if (mag(d) > VSMALL)
+    {
+        scalar s = (iso_-s0)/d;
+
+        if (s >= 0.5 && s <= 1 && p1Index != -1)
+        {
+            return snappedPoints[p1Index];
+        }
+        else if (s >= 0.0 && s <= 0.5 && p0Index != -1)
+        {
+            return snappedPoints[p0Index];
+        }
+        else
+        {
+            return s*p1 + (1.0-s)*p0;
+        }
+    }
+    else
+    {
+        scalar s = 0.4999;
+
+        return s*p1 + (1.0-s)*p0;
+    }
+}
+
+
+template<class Type>
+void Foam::isoSurface::generateTriPoints
+(
+    const DynamicList<Type>& snapped,
+
+    const scalar s0,
+    const Type& p0,
+    const label p0Index,
+
+    const scalar s1,
+    const Type& p1,
+    const label p1Index,
+
+    const scalar s2,
+    const Type& p2,
+    const label p2Index,
+
+    const scalar s3,
+    const Type& p3,
+    const label p3Index,
+
+    DynamicList<Type>& points
+) const
+{
+    int triIndex = 0;
+    if (s0 < iso_)
+    {
+        triIndex |= 1;
+    }
+    if (s1 < iso_)
+    {
+        triIndex |= 2;
+    }
+    if (s2 < iso_)
+    {
+        triIndex |= 4;
+    }
+    if (s3 < iso_)
+    {
+        triIndex |= 8;
+    }
+
+    /* Form the vertices of the triangles for each case */
+    switch (triIndex)
+    {
+        case 0x00:
+        case 0x0F:
+        break;
+
+        case 0x0E:
+        case 0x01:
+            points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+        break;
+
+        case 0x0D:
+        case 0x02:
+            points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+        break;
+
+        case 0x0C:
+        case 0x03:
+        {
+            Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append(tp1);
+            points.append(tp2);
+            points.append(tp2);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x0B:
+        case 0x04:
+        {
+            points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
+            points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+        }
+        break;
+
+        case 0x0A:
+        case 0x05:
+        {
+            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            points.append(tp0);
+            points.append(tp1);
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+            points.append(tp0);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x09:
+        case 0x06:
+        {
+            Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+            Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+            points.append(tp0);
+            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
+            points.append(tp1);
+            points.append(tp0);
+            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x07:
+        case 0x08:
+            points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
+            points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
+            points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+        break;
+    }
+}
+
+
+template<class Type>
+Foam::label Foam::isoSurface::generateTriPoints
+(
+    const volScalarField& cVals,
+    const scalarField& pVals,
+
+    const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+    const Field<Type>& pCoords,
+
+    const DynamicList<Type>& snappedPoints,
+    const labelList& snappedCc,
+    const labelList& snappedPoint,
+    const label faceI,
+
+    const scalar neiVal,
+    const Type& neiPt,
+    const label neiSnap,
+
+    DynamicList<Type>& triPoints,
+    DynamicList<label>& triMeshCells
+) const
+{
+    label own = mesh_.faceOwner()[faceI];
+
+    label oldNPoints = triPoints.size();
+
+    const face& f = mesh_.faces()[faceI];
+
+    forAll(f, fp)
+    {
+        label pointI = f[fp];
+        label nextPointI = f[f.fcIndex(fp)];
+
+        generateTriPoints
+        (
+            snappedPoints,
+
+            pVals[pointI],
+            pCoords[pointI],
+            snappedPoint[pointI],
+
+            pVals[nextPointI],
+            pCoords[nextPointI],
+            snappedPoint[nextPointI],
+
+            cVals[own],
+            cCoords[own],
+            snappedCc[own],
+
+            neiVal,
+            neiPt,
+            neiSnap,
+
+            triPoints
+        );
+    }
+
+    // Every three triPoints is a triangle
+    label nTris = (triPoints.size()-oldNPoints)/3;
+    for (label i = 0; i < nTris; i++)
+    {
+        triMeshCells.append(own);
+    }
+
+    return nTris;
+}
+
+
+template<class Type>
+void Foam::isoSurface::generateTriPoints
+(
+    const volScalarField& cVals,
+    const scalarField& pVals,
+
+    const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+    const Field<Type>& pCoords,
+
+    const DynamicList<Type>& snappedPoints,
+    const labelList& snappedCc,
+    const labelList& snappedPoint,
+
+    DynamicList<Type>& triPoints,
+    DynamicList<label>& triMeshCells
+) const
+{
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+
+    // Determine neighbouring snap status
+    labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            forAll(pp, i)
+            {
+                neiSnappedCc[faceI-mesh_.nInternalFaces()] =
+                    snappedCc[own[faceI]];
+                faceI++;
+            }
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh_, neiSnappedCc, false);
+
+
+
+    // Generate triangle points
+
+    triPoints.clear();
+    triMeshCells.clear();
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        if (faceCutType_[faceI] != NOTCUT)
+        {
+            generateTriPoints
+            (
+                cVals,
+                pVals,
+
+                cCoords,
+                pCoords,
+
+                snappedPoints,
+                snappedCc,
+                snappedPoint,
+                faceI,
+
+                cVals[nei[faceI]],
+                cCoords[nei[faceI]],
+                snappedCc[nei[faceI]],
+
+                triPoints,
+                triMeshCells
+            );
+        }
+    }
+
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            if (refCast<const processorPolyPatch>(pp).owner())
+            {
+                label faceI = pp.start();
+
+                forAll(pp, i)
+                {
+                    if (faceCutType_[faceI] != NOTCUT)
+                    {
+                        generateTriPoints
+                        (
+                            cVals,
+                            pVals,
+
+                            cCoords,
+                            pCoords,
+
+                            snappedPoints,
+                            snappedCc,
+                            snappedPoint,
+                            faceI,
+
+                            cVals.boundaryField()[patchI][i],
+                            cCoords.boundaryField()[patchI][i],
+                            neiSnappedCc[faceI-mesh_.nInternalFaces()],
+
+                            triPoints,
+                            triMeshCells
+                        );
+                    }
+                    faceI++;
+                }
+            }
+        }
+        else
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                if (faceCutType_[faceI] != NOTCUT)
+                {
+                    generateTriPoints
+                    (
+                        cVals,
+                        pVals,
+
+                        cCoords,
+                        pCoords,
+
+                        snappedPoints,
+                        snappedCc,
+                        snappedPoint,
+                        faceI,
+
+                        cVals.boundaryField()[patchI][i],
+                        cCoords.boundaryField()[patchI][i],
+                        -1, // fc not snapped
+
+                        triPoints,
+                        triMeshCells
+                    );
+                }
+                faceI++;
+            }
+        }
+    }
+
+    triPoints.shrink();
+    triMeshCells.shrink();
+}
+
+
+//template <class Type>
+//Foam::tmp<Foam::Field<Type> >
+//Foam::isoSurface::sample(const Field<Type>& vField) const
+//{
+//    return tmp<Field<Type> >(new Field<Type>(vField, meshCells()));
+//}
+//
+//
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::isoSurface::interpolate
+(
+    const volScalarField& cVals,
+    const scalarField& pVals,
+    const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+    const Field<Type>& pCoords
+) const
+{
+    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<label> triMeshCells(nCutCells_);
+
+    // Dummy snap data
+    DynamicList<Type> snappedPoints;
+    labelList snappedCc(mesh_.nCells(), -1);
+    labelList snappedPoint(mesh_.nPoints(), -1);
+
+    generateTriPoints
+    (
+        cVals,
+        pVals,
+
+        cCoords,
+        pCoords,
+
+        snappedPoints,
+        snappedCc,
+        snappedPoint,
+
+        triPoints,
+        triMeshCells
+    );
+
+
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    Field<Type>& values = tvalues();
+
+    forAll(triPoints, i)
+    {
+        label mergedPointI = triPointMergeMap_[i];
+
+        if (mergedPointI >= 0)
+        {
+            values[mergedPointI] = triPoints[i];
+        }
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C

@@ -41,6 +41,100 @@ namespace Foam
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::sampledIsoSurface::getIsoFields() const
+{
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+    // Get volField
+    // ~~~~~~~~~~~~
+
+    if (fvm.foundObject<volScalarField>(isoField_))
+    {
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : lookup "
+                << isoField_ << endl;
+        }
+        storedVolFieldPtr_.clear();
+        volFieldPtr_ = &fvm.lookupObject<volScalarField>(isoField_);
+    }
+    else
+    {
+        // Bit of a hack. Read field and store.
+
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : reading "
+                << isoField_ << " from time " <<fvm.time().timeName()
+                << endl;
+        }
+
+        storedVolFieldPtr_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    isoField_,
+                    fvm.time().timeName(),
+                    fvm,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                fvm
+            )
+        );
+        volFieldPtr_ = storedVolFieldPtr_.operator->();
+    }
+
+
+
+    // Get pointField
+    // ~~~~~~~~~~~~~~
+
+    word pointFldName = "volPointInterpolate(" + isoField_ + ')';
+
+    if (fvm.foundObject<pointScalarField>(pointFldName))
+    {
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : lookup "
+                << pointFldName << endl;
+        }
+        storedPointFieldPtr_.clear();
+        pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
+    }
+    else
+    {
+        // Not in registry. Interpolate.
+
+        if (debug)
+        {
+            Info<< "sampledIsoSurface::getIsoField() : interpolating "
+                << pointFldName << endl;
+        }
+
+        storedPointFieldPtr_.reset
+        (
+            volPointInterpolation::New(fvm).interpolate(*volFieldPtr_).ptr()
+        );
+        pointFieldPtr_ = storedPointFieldPtr_.operator->();
+    }
+
+    if (debug)
+    {
+        Info<< "sampledIsoSurface::getIsoField() : obtained volField "
+            << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
+            << " max:" << max(*volFieldPtr_).value() << endl;
+        Info<< "sampledIsoSurface::getIsoField() : obtained pointField "
+            << pointFieldPtr_->name()
+            << " min:" << gMin(pointFieldPtr_->internalField())
+            << " max:" << gMax(pointFieldPtr_->internalField()) << endl;
+    }
+}
+
+
 void Foam::sampledIsoSurface::createGeometry() const
 {
     const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
@@ -49,66 +143,29 @@ void Foam::sampledIsoSurface::createGeometry() const
     {
         storedTimeIndex_ = fvm.time().timeIndex();
 
+        getIsoFields();
+
         // Clear any stored topo
+        surfPtr_.clear();
         facesPtr_.clear();
 
-        // Optionally read volScalarField
-        autoPtr<volScalarField> readFieldPtr_;
-
-        // 1. see if field in database
-        // 2. see if field can be read
-        const volScalarField* cellFldPtr = NULL;
-        if (fvm.foundObject<volScalarField>(isoField_))
-        {
-            if (debug)
-            {
-                Info<< "sampledIsoSurface::createGeometry() : lookup "
-                    << isoField_ << endl;
-            }
-
-            cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_);
-        }
-        else
-        {
-            // Bit of a hack. Read field and store.
-
-            if (debug)
-            {
-                Info<< "sampledIsoSurface::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));
+            volScalarField cellAvg
+            (
+                IOobject
+                (
+                    "cellAvg",
+                    fvm.time().timeName(),
+                    fvm.time(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                fvm,
+                dimensionedScalar("zero", dimless, scalar(0.0))
+            );
             labelField nPointCells(fvm.nCells(), 0);
             {
                 for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
@@ -119,7 +176,7 @@ void Foam::sampledIsoSurface::createGeometry() const
                     {
                         label cellI = pCells[i];
 
-                        cellAvg[cellI] += pointFld().internalField()[pointI];
+                        cellAvg[cellI] += (*pointFieldPtr_)[pointI];
                         nPointCells[cellI]++;
                     }
                 }
@@ -128,33 +185,33 @@ void Foam::sampledIsoSurface::createGeometry() const
             {
                 cellAvg[cellI] /= nPointCells[cellI];
             }
+            // Give value to calculatedFvPatchFields
+            cellAvg.correctBoundaryConditions();
 
-            const isoSurface iso
+            surfPtr_.reset
             (
-                fvm,
+                new isoSurface
-                cellAvg,
+                (
-                pointFld().internalField(),
+                    cellAvg,
-                isoVal_,
+                    *pointFieldPtr_,
-                regularise_
+                    isoVal_,
+                    regularise_
+                )
             );
-
-            const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
-            meshCells_ = iso.meshCells();
         }
         else
         {
-            //- Direct from cell field and point field. Gives bad continuity.
+            //- Direct from cell field and point field.
-            const isoSurface iso
+            surfPtr_.reset
             (
-                fvm,
+                new isoSurface
-                cellFld.internalField(),
+                (
-                pointFld().internalField(),
+                    *volFieldPtr_,
-                isoVal_,
+                    *pointFieldPtr_,
-                regularise_
+                    isoVal_,
+                    regularise_
+                )
             );
-
-            const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
-            meshCells_ = iso.meshCells();
         }
 
 
@@ -167,8 +224,9 @@ void Foam::sampledIsoSurface::createGeometry() const
                 << "    isoField       : " << isoField_ << nl
                 << "    isoValue       : " << isoVal_ << nl
                 << "    points         : " << points().size() << nl
-                << "    tris           : " << triSurface::size() << nl
+                << "    tris           : " << surface().size() << nl
-                << "    cut cells      : " << meshCells_.size() << endl;
+                << "    cut cells      : " << surface().meshCells().size()
+                << endl;
         }
     }
 }
@@ -187,12 +245,26 @@ Foam::sampledIsoSurface::sampledIsoSurface
     isoField_(dict.lookup("isoField")),
     isoVal_(readScalar(dict.lookup("isoValue"))),
     regularise_(dict.lookupOrDefault("regularise", true)),
-    average_(dict.lookupOrDefault("average", true)),
+    average_(dict.lookupOrDefault("average", false)),
     zoneName_(word::null),
+    surfPtr_(NULL),
     facesPtr_(NULL),
     storedTimeIndex_(-1),
-    meshCells_(0)
+    storedVolFieldPtr_(NULL),
+    volFieldPtr_(NULL),
+    storedPointFieldPtr_(NULL),
+    pointFieldPtr_(NULL)
 {
+    if (!sampledSurface::interpolate())
+    {
+        FatalErrorIn
+        (
+            "sampledIsoSurface::sampledIsoSurface"
+            "(const word&, const polyMesh&, const dictionary&)"
+        )   << "Non-interpolated iso surface not supported since triangles"
+            << " span across cells." << exit(FatalError);
+    }
+
 //    label zoneId = -1;
 //    if (dict.readIfPresent("zone", zoneName_))
 //    {
@@ -220,6 +292,7 @@ void Foam::sampledIsoSurface::correct(const bool meshChanged)
     // Only change of mesh changes plane - zone restriction gets lost
     if (meshChanged)
     {
+        surfPtr_.clear();
         facesPtr_.clear();
     }
 }
@@ -328,9 +401,9 @@ void Foam::sampledIsoSurface::print(Ostream& os) const
 {
     os  << "sampledIsoSurface: " << name() << " :"
         << "  field:" << isoField_
-        << "  value:" << isoVal_
+        << "  value:" << isoVal_;
-        << "  faces:" << faces().size()
+        //<< "  faces:" << faces().size()       // note: possibly no geom yet
-        << "  points:" << points().size();
+        //<< "  points:" << points().size();
 }
 
 

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H

@@ -39,7 +39,7 @@ SourceFiles
 #define sampledIsoSurface_H
 
 #include "sampledSurface.H"
-#include "triSurface.H"
+#include "isoSurface.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,8 +52,7 @@ namespace Foam
 
 class sampledIsoSurface
 :
-    public sampledSurface,
+    public sampledSurface
-    public triSurface
 {
     // Private data
 
@@ -72,6 +71,8 @@ class sampledIsoSurface
         //- zone name (if restricted to zones)
         word zoneName_;
 
+        mutable autoPtr<isoSurface> surfPtr_;
+
         //- triangles converted to faceList
         mutable autoPtr<faceList> facesPtr_;
 
@@ -81,12 +82,21 @@ class sampledIsoSurface
             //- Time at last call
             mutable label storedTimeIndex_;
 
-            //- For every triangle the original cell in mesh
+            //- Cached volfield
-            mutable labelList meshCells_;
+            mutable autoPtr<volScalarField> storedVolFieldPtr_;
+            mutable const volScalarField* volFieldPtr_;
+
+            //- Cached pointfield
+            mutable autoPtr<pointScalarField> storedPointFieldPtr_;
+            mutable const pointScalarField* pointFieldPtr_;
+
 
 
     // Private Member Functions
 
+        //- Get fields needed to recreate iso surface.
+        void getIsoFields() const;
+
         //- Create iso surface (if time has changed)
         void createGeometry() const;
 
@@ -130,7 +140,7 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return triSurface::points();
+            return surface().points();
         }
 
         //- Faces of surface
@@ -138,7 +148,7 @@ public:
         {
             if (!facesPtr_.valid())
             {
-                const triSurface& s = *this;
+                const triSurface& s = surface();
 
                 facesPtr_.reset(new faceList(s.size()));
 
@@ -150,11 +160,19 @@ public:
             return facesPtr_;
         }
 
-
         //- Correct for mesh movement and/or field changes
         virtual void correct(const bool meshChanged);
 
 
+        const isoSurface& surface() const
+        {
+            return surfPtr_();
+        }
+
+        //- Lookup or read isoField. Sets volFieldPtr_ and pointFieldPtr_.
+        void getIsoField();
+
+
         //- sample field on surface
         virtual tmp<scalarField> sample
         (

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C

@@ -0,0 +1,344 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledIsoSurfaceCell.H"
+#include "dictionary.H"
+#include "volFields.H"
+#include "volPointInterpolation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "isoSurfaceCell.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sampledIsoSurfaceCell, 0);
+    addNamedToRunTimeSelectionTable(sampledSurface, sampledIsoSurfaceCell, word, isoSurfaceCell);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sampledIsoSurfaceCell::createGeometry() const
+{
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+    if (fvm.time().timeIndex() != storedTimeIndex_)
+    {
+        storedTimeIndex_ = fvm.time().timeIndex();
+
+        // Clear any stored topo
+        facesPtr_.clear();
+
+        // Optionally read volScalarField
+        autoPtr<volScalarField> readFieldPtr_;
+
+        // 1. see if field in database
+        // 2. see if field can be read
+        const volScalarField* cellFldPtr = NULL;
+        if (fvm.foundObject<volScalarField>(isoField_))
+        {
+            if (debug)
+            {
+                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();
+        }
+
+
+        if (debug)
+        {
+            Pout<< "sampledIsoSurfaceCell::createGeometry() : 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;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    isoField_(dict.lookup("isoField")),
+    isoVal_(readScalar(dict.lookup("isoValue"))),
+    regularise_(dict.lookupOrDefault("regularise", true)),
+    average_(dict.lookupOrDefault("average", true)),
+    zoneName_(word::null),
+    facesPtr_(NULL),
+    storedTimeIndex_(-1),
+    meshCells_(0)
+{
+//    label zoneId = -1;
+//    if (dict.readIfPresent("zone", zoneName_))
+//    {
+//        zoneId = mesh.cellZones().findZoneID(zoneName_);
+//        if (debug && zoneId < 0)
+//        {
+//            Info<< "cellZone \"" << zoneName_
+//                << "\" not found - using entire mesh"
+//                << endl;
+//        }
+//    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfaceCell::~sampledIsoSurfaceCell()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::sampledIsoSurfaceCell::correct(const bool meshChanged)
+{
+    // Only change of mesh changes plane - zone restriction gets lost
+    if (meshChanged)
+    {
+        facesPtr_.clear();
+    }
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volScalarField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volVectorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volSphericalTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volSymmTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurfaceCell::sample
+(
+    const volTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurfaceCell::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+void Foam::sampledIsoSurfaceCell::print(Ostream& os) const
+{
+    os  << "sampledIsoSurfaceCell: " << name() << " :"
+        << "  field:" << isoField_
+        << "  value:" << isoVal_;
+        //<< "  faces:" << faces().size()   // possibly no geom yet
+        //<< "  points:" << points().size();
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H

@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::sampledIsoSurfaceCell
+
+Description
+    A sampledSurface defined by a surface of iso value. Always triangulated.
+    To be used in sampleSurfaces / functionObjects. Recalculates iso surface
+    only if time changes.
+
+SourceFiles
+    sampledIsoSurfaceCell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledIsoSurfaceCell_H
+#define sampledIsoSurfaceCell_H
+
+#include "sampledSurface.H"
+#include "triSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class sampledIsoSurfaceCell Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledIsoSurfaceCell
+:
+    public sampledSurface,
+    public triSurface
+{
+    // Private data
+
+        //- Field to get isoSurface of
+        const word isoField_;
+
+        //- iso value
+        const scalar isoVal_;
+
+        //- Whether to coarse
+        const Switch regularise_;
+
+        //- Whether to recalculate cell values as average of point values
+        const Switch average_;
+
+        //- zone name (if restricted to zones)
+        word zoneName_;
+
+        //- triangles converted to faceList
+        mutable autoPtr<faceList> facesPtr_;
+
+
+        // Recreated for every isoSurface
+
+            //- Time at last call
+            mutable label storedTimeIndex_;
+
+            //- For every triangle the original cell in mesh
+            mutable labelList meshCells_;
+
+
+    // Private Member Functions
+
+        //- Create iso surface (if time has changed)
+        void createGeometry() const;
+
+        //- sample field on faces
+        template <class Type>
+        tmp<Field<Type> > sampleField
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vField
+        ) const;
+
+
+        template <class Type>
+        tmp<Field<Type> >
+        interpolateField(const interpolation<Type>&) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sampledIsoSurfaceCell");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        sampledIsoSurfaceCell
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~sampledIsoSurfaceCell();
+
+
+    // Member Functions
+
+        //- Points of surface
+        virtual const pointField& points() const
+        {
+            return triSurface::points();
+        }
+
+        //- Faces of surface
+        virtual const faceList& faces() const
+        {
+            if (!facesPtr_.valid())
+            {
+                const triSurface& s = *this;
+
+                facesPtr_.reset(new faceList(s.size()));
+
+                forAll(s, i)
+                {
+                    facesPtr_()[i] = s[i].triFaceFace();
+                }
+            }
+            return facesPtr_;
+        }
+
+
+        //- Correct for mesh movement and/or field changes
+        virtual void correct(const bool meshChanged);
+
+
+        //- sample field on surface
+        virtual tmp<scalarField> sample
+        (
+            const volScalarField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<vectorField> sample
+        (
+            const volVectorField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<sphericalTensorField> sample
+        (
+            const volSphericalTensorField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<symmTensorField> sample
+        (
+            const volSymmTensorField&
+        ) const;
+
+        //- sample field on surface
+        virtual tmp<tensorField> sample
+        (
+            const volTensorField&
+        ) const;
+
+
+        //- interpolate field on surface
+        virtual tmp<scalarField> interpolate
+        (
+            const interpolation<scalar>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<vectorField> interpolate
+        (
+            const interpolation<vector>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<sphericalTensorField> interpolate
+        (
+            const interpolation<sphericalTensor>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<symmTensorField> interpolate
+        (
+            const interpolation<symmTensor>&
+        ) const;
+
+        //- interpolate field on surface
+        virtual tmp<tensorField> interpolate
+        (
+            const interpolation<tensor>&
+        ) const;
+
+        //- Write
+        virtual void print(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "sampledIsoSurfaceCellTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "sampledIsoSurfaceCell.H"
+#include "isoSurface.H"
+#include "volFieldsFwd.H"
+#include "pointFields.H"
+#include "volPointInterpolation.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::sampledIsoSurfaceCell::sampleField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+) const
+{
+    // Recreate geometry if time has changed
+    createGeometry();
+
+    return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
+}
+
+
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::sampledIsoSurfaceCell::interpolateField
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    // Recreate geometry if time has changed
+    createGeometry();
+
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    Field<Type>& values = tvalues();
+
+    boolList pointDone(points().size(), false);
+
+    forAll(faces(), cutFaceI)
+    {
+        const face& f = faces()[cutFaceI];
+
+        forAll(f, faceVertI)
+        {
+            label pointI = f[faceVertI];
+
+            if (!pointDone[pointI])
+            {
+                values[pointI] = interpolator.interpolate
+                (
+                    points()[pointI],
+                    meshCells_[cutFaceI]
+                );
+                pointDone[pointI] = true;
+            }
+        }
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sampledIsoSurface.H"
-#include "isoSurface.H"
 #include "volFieldsFwd.H"
 #include "pointFields.H"
 #include "volPointInterpolation.H"
@@ -41,8 +40,7 @@ Foam::sampledIsoSurface::sampleField
 {
     // Recreate geometry if time has changed
     createGeometry();
-
+    return tmp<Field<Type> >(new Field<Type>(vField, surface().meshCells()));
-    return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
 }
 
 
@@ -53,36 +51,28 @@ Foam::sampledIsoSurface::interpolateField
     const interpolation<Type>& interpolator
 ) const
 {
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+    // Get fields to sample. Assume volPointInterpolation!
+    const GeometricField<Type, fvPatchField, volMesh>& volFld =
+        interpolator.psi();
+
+    tmp<GeometricField<Type, pointPatchField, pointMesh> > pointFld
+    (
+        volPointInterpolation::New(fvm).interpolate(volFld)
+    );
+
     // Recreate geometry if time has changed
     createGeometry();
 
-    // One value per point
+    // Sample.
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    return surface().interpolate
-    Field<Type>& values = tvalues();
+    (
-
+        *volFieldPtr_,
-    boolList pointDone(points().size(), false);
+        *pointFieldPtr_,
-
+        volFld,
-    forAll(faces(), cutFaceI)
+        pointFld()
-    {
+    );
-        const face& f = faces()[cutFaceI];
-
-        forAll(f, faceVertI)
-        {
-            label pointI = f[faceVertI];
-
-            if (!pointDone[pointI])
-            {
-                values[pointI] = interpolator.interpolate
-                (
-                    points()[pointI],
-                    meshCells_[cutFaceI]
-                );
-                pointDone[pointI] = true;
-            }
-        }
-    }
-
-    return tvalues;
 }