From cdb0a8a8d517b26f91a551f3c3af335fa16dd8fb Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 21 Oct 2008 14:06:10 +0100
Subject: [PATCH] iso surface

---
 .../sampledSurface/isoSurface/isoSurface.C    | 291 +++++++++++++++++
 .../sampledSurface/isoSurface/isoSurface.H    | 173 +++++++++++
 .../isoSurface/isoSurfaceTemplates.C          | 292 ++++++++++++++++++
 .../isoSurface/sampledIsoSurface.C            | 204 ++++++++++++
 .../isoSurface/sampledIsoSurface.H            | 232 ++++++++++++++
 .../isoSurface/sampledIsoSurfaceTemplates.C   | 144 +++++++++
 6 files changed, 1336 insertions(+)
 create mode 100644 src/sampling/sampledSurface/isoSurface/isoSurface.C
 create mode 100644 src/sampling/sampledSurface/isoSurface/isoSurface.H
 create mode 100644 src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
 create mode 100644 src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
 create mode 100644 src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
 create mode 100644 src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C

diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
new file mode 100644
index 0000000000..45be2c1334
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -0,0 +1,291 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "dictionary.H"
+#include "polyMesh.H"
+#include "volFields.H"
+#include "mergePoints.H"
+#include "tetMatcher.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(isoSurface, 0);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::isoSurface::isoSurface
+(
+    const polyMesh& mesh,
+    const scalarField& cellValues,
+    const scalarField& pointValues,
+    const scalar iso
+)
+:
+    mesh_(mesh),
+    cellValues_(cellValues),
+    pointValues_(pointValues),
+    iso_(iso)
+{
+    const pointField& cellCentres = mesh.cellCentres();
+
+    tetMatcher tet;
+
+    DynamicList<point> triPoints;
+    DynamicList<label> triMeshCells;
+
+    forAll(mesh.cells(), cellI)
+    {
+        label oldNPoints = triPoints.size();
+
+        const cell& cFaces = mesh.cells()[cellI];
+
+        if (tet.isA(mesh, cellI))
+        {
+            // For tets don't do cell-centre decomposition, just use the
+            // tet points and values
+
+            const face& f0 = mesh.faces()[cFaces[0]];
+
+            // Get the other point
+            const face& f1 = mesh.faces()[cFaces[1]];
+            label oppositeI = -1;
+            forAll(f1, fp)
+            {
+                oppositeI = f1[fp];
+
+                if (findIndex(f0, oppositeI) == -1)
+                {
+                    break;
+                }
+            }
+
+            vertexInterp
+            (
+                iso,
+                pointValues[f0[0]],
+                pointValues[f0[1]],
+                pointValues[f0[2]],
+                pointValues[oppositeI],
+
+                mesh.points()[f0[0]],
+                mesh.points()[f0[1]],
+                mesh.points()[f0[2]],
+                mesh.points()[oppositeI],
+
+                triPoints
+            );
+        }
+        else
+        {
+            const cell& cFaces = mesh.cells()[cellI];
+
+            forAll(cFaces, cFaceI)
+            {
+                label faceI = cFaces[cFaceI];
+                const face& f = mesh.faces()[faceI];
+
+                // Do a tetrahedrisation. Each face to cc becomes pyr.
+                // Each pyr gets split into two tets by diagionalisation
+                // of face. So
+                // - f[0], f[1], f[2], cc
+                // - f[0], f[2], f[3], cc
+
+                for(label fp = 1; fp < f.size() - 1; fp++)
+                {
+                    vertexInterp
+                    (
+                        iso,
+                        pointValues[f[0]],
+                        pointValues[f[fp]],
+                        pointValues[f[f.fcIndex(fp)]],
+                        cellValues[cellI],
+
+                        mesh.points()[f[0]],
+                        mesh.points()[f[fp]],
+                        mesh.points()[f[f.fcIndex(fp)]],
+                        cellCentres[cellI],
+
+                        triPoints
+                    );
+                }
+            }
+        }
+
+
+        // Every three triPoints is a cell
+        label nCells = (triPoints.size()-oldNPoints)/3;
+        for (label i = 0; i < nCells; i++)
+        {
+            triMeshCells.append(cellI);
+        }
+    }
+
+    triPoints.shrink();
+    triMeshCells.shrink();
+    meshCells_.transfer(triMeshCells);
+
+    pointField newPoints;
+    mergePoints(triPoints, SMALL, false, triPointMergeMap_, newPoints);
+
+    DynamicList<labelledTri> tris(meshCells_.size());
+    forAll(meshCells_, triI)
+    {
+        tris.append
+        (
+            labelledTri
+            (
+                triPointMergeMap_[3*triI],
+                triPointMergeMap_[3*triI+1],
+                triPointMergeMap_[3*triI+2],
+                0
+            )
+        );
+    }
+
+    //- 1.5.x:
+    //tris.shrink();
+    //triSurface::operator=
+    //(
+    //    triSurface(tris, geometricSurfacePatchList(0), newPoints)
+    //);
+
+    triSurface::operator=
+    (
+        triSurface(tris, geometricSurfacePatchList(0), newPoints, true)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//Foam::tmp<Foam::scalarField>
+//Foam::isoSurface::sample
+//(
+//    const volScalarField& vField
+//) const
+//{
+//    return sampleField(vField);
+//}
+//
+//
+//Foam::tmp<Foam::vectorField>
+//Foam::isoSurface::sample
+//(
+//    const volVectorField& vField
+//) const
+//{
+//    return sampleField(vField);
+//}
+//
+//Foam::tmp<Foam::sphericalTensorField>
+//Foam::isoSurface::sample
+//(
+//    const volSphericalTensorField& vField
+//) const
+//{
+//    return sampleField(vField);
+//}
+//
+//
+//Foam::tmp<Foam::symmTensorField>
+//Foam::isoSurface::sample
+//(
+//    const volSymmTensorField& vField
+//) const
+//{
+//    return sampleField(vField);
+//}
+//
+//
+//Foam::tmp<Foam::tensorField>
+//Foam::isoSurface::sample
+//(
+//    const volTensorField& vField
+//) const
+//{
+//    return sampleField(vField);
+//}
+//
+//
+//Foam::tmp<Foam::scalarField>
+//Foam::isoSurface::interpolate
+//(
+//    const interpolation<scalar>& interpolator
+//) const
+//{
+//    return interpolateField(interpolator);
+//}
+//
+//
+//Foam::tmp<Foam::vectorField>
+//Foam::isoSurface::interpolate
+//(
+//    const interpolation<vector>& interpolator
+//) const
+//{
+//    return interpolateField(interpolator);
+//}
+//
+//Foam::tmp<Foam::sphericalTensorField>
+//Foam::isoSurface::interpolate
+//(
+//    const interpolation<sphericalTensor>& interpolator
+//) const
+//{
+//    return interpolateField(interpolator);
+//}
+//
+//
+//Foam::tmp<Foam::symmTensorField>
+//Foam::isoSurface::interpolate
+//(
+//    const interpolation<symmTensor>& interpolator
+//) const
+//{
+//    return interpolateField(interpolator);
+//}
+//
+//
+//Foam::tmp<Foam::tensorField>
+//Foam::isoSurface::interpolate
+//(
+//    const interpolation<tensor>& interpolator
+//) const
+//{
+//    return interpolateField(interpolator);
+//}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
new file mode 100644
index 0000000000..ddd6854909
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::isoSurface
+
+Description
+    A surface formed by the iso value.
+    After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke
+
+SourceFiles
+    isoSurface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoSurface_H
+#define isoSurface_H
+
+#include "triSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                       Class isoSurface Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoSurface
+:
+    public triSurface
+{
+    // Private data
+
+        //- Reference to mesh
+        const polyMesh& mesh_;
+
+        //- Reference to cell values
+        const scalarField& cellValues_;
+
+        //- Reference to point values
+        const scalarField& pointValues_;
+
+        //- Isosurface value
+        const scalar iso_;
+
+
+        //- 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
+
+        template<class T>
+        static T vertexInterp
+        (
+            const scalar iso,
+            const T& p0,
+            const T& p1,
+            const scalar s0,
+            const scalar s1
+        );
+
+        template<class T>
+        static void vertexInterp
+        (
+            const scalar iso,
+            const scalar s0,
+            const scalar s1,
+            const scalar s2,
+            const scalar s3,
+
+            const T& p0,
+            const T& p1,
+            const T& p2,
+            const T& p3,
+
+            DynamicList<T>& points
+        );
+
+public:
+
+    //- Runtime type information
+    TypeName("isoSurface");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        isoSurface
+        (
+            const polyMesh& mesh,
+            const scalarField& cellValues,
+            const scalarField& pointValues,
+            const scalar iso
+        );
+
+
+    // 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_;
+        }
+
+
+        //- sample field on faces
+        template <class Type>
+        tmp<Field<Type> > sample
+        (
+            const Field<Type>& sampleCellValues
+        ) const;
+
+        //- interpolate field to points
+        template <class Type>
+        tmp<Field<Type> >
+        interpolate
+        (
+            const Field<Type>& sampleCellValues,
+            const Field<Type>& samplePointValues
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "isoSurfaceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
new file mode 100644
index 0000000000..e1500ca6e5
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -0,0 +1,292 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "tetMatcher.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::isoSurface::vertexInterp
+(
+    const scalar iso,
+    const Type& p0,
+    const Type& p1,
+    const scalar s0,
+    const scalar s1
+)
+{
+    scalar d = s1-s0;
+
+    if (mag(d) > VSMALL)
+    {
+        return (iso-s0)/d*p1 + (s1-iso)/d*p0;
+    }
+    else
+    {
+        return 0.5*(p0+p1);
+    }
+}
+
+
+// After "Polygonising A Scalar Field Using Tetrahedrons"
+// by Paul Bourke
+// Get value consistent with uncompacted triangle points.
+// Given tet corner sample values s0..s3 interpolate the corresponding
+// values p0..p3 to construct the surface corresponding to sample value iso.
+template<class Type>
+void Foam::isoSurface::vertexInterp
+(
+    const scalar iso,
+    const scalar s0,
+    const scalar s1,
+    const scalar s2,
+    const scalar s3,
+
+    const Type& p0,
+    const Type& p1,
+    const Type& p2,
+    const Type& p3,
+
+    DynamicList<Type>& points
+)
+{
+    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(vertexInterp(iso,p0,p1,s0,s1));
+            points.append(vertexInterp(iso,p0,p2,s0,s2));
+            points.append(vertexInterp(iso,p0,p3,s0,s3));
+        break;
+
+        case 0x0D:
+        case 0x02:
+            points.append(vertexInterp(iso,p1,p0,s1,s0));
+            points.append(vertexInterp(iso,p1,p3,s1,s3));
+            points.append(vertexInterp(iso,p1,p2,s1,s2));
+        break;
+
+        case 0x0C:
+        case 0x03:
+        {
+            const Type tp1 = vertexInterp(iso,p0,p2,s0,s2);
+            const Type tp2 = vertexInterp(iso,p1,p3,s1,s3);
+
+            points.append(vertexInterp(iso,p0,p3,s0,s3));
+            points.append(tp1);
+            points.append(tp2);
+            points.append(tp2);
+            points.append(vertexInterp(iso,p1,p2,s1,s2));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x0B:
+        case 0x04:
+        {
+            points.append(vertexInterp(iso,p2,p0,s2,s0));
+            points.append(vertexInterp(iso,p2,p1,s2,s1));
+            points.append(vertexInterp(iso,p2,p3,s2,s3));
+        }
+        break;
+
+        case 0x0A:
+        case 0x05:
+        {
+            const Type tp0 = vertexInterp(iso,p0,p1,s0,s1);
+            const Type tp1 = vertexInterp(iso,p2,p3,s2,s3);
+
+            points.append(tp0);
+            points.append(tp1);
+            points.append(vertexInterp(iso,p0,p3,s0,s3));
+            points.append(tp0);
+            points.append(vertexInterp(iso,p1,p2,s1,s2));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x09:
+        case 0x06:
+        {
+            const Type tp0 = vertexInterp(iso,p0,p1,s0,s1);
+            const Type tp1 = vertexInterp(iso,p2,p3,s2,s3);
+
+            points.append(tp0);
+            points.append(vertexInterp(iso,p1,p3,s1,s3));
+            points.append(tp1);
+            points.append(tp0);
+            points.append(vertexInterp(iso,p0,p2,s0,s2));
+            points.append(tp1);
+        }
+        break;
+
+        case 0x07:
+        case 0x08:
+            points.append(vertexInterp(iso,p3,p0,s3,s0));
+            points.append(vertexInterp(iso,p3,p2,s3,s2));
+            points.append(vertexInterp(iso,p3,p1,s3,s1));
+        break;
+    }
+}
+
+
+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 Field<Type>& sampleCellValues,
+    const Field<Type>& samplePointValues
+) const
+{
+    tetMatcher tet;
+
+    DynamicList<Type> triValues;
+
+    // Note: in same order as construction of triSurface
+    label oldCellI = -1;
+    forAll(meshCells_, triI)
+    {
+        label cellI = meshCells_[triI];
+
+        if (cellI != oldCellI)
+        {
+            oldCellI = cellI;
+
+            const cell& cFaces = mesh_.cells()[cellI];
+
+            if (tet.isA(mesh_, cellI))
+            {
+                // For tets don't do cell-centre decomposition, just use the
+                // tet points and values
+
+                const face& f0 = mesh_.faces()[cFaces[0]];
+
+                // Get the other point
+                const face& f1 = mesh_.faces()[cFaces[1]];
+                label oppositeI = -1;
+                forAll(f1, fp)
+                {
+                    oppositeI = f1[fp];
+
+                    if (findIndex(f0, oppositeI) == -1)
+                    {
+                        break;
+                    }
+                }
+
+                vertexInterp
+                (
+                    iso_,
+                    pointValues_[f0[0]],
+                    pointValues_[f0[1]],
+                    pointValues_[f0[2]],
+                    pointValues_[oppositeI],
+
+                    samplePointValues[f0[0]],
+                    samplePointValues[f0[1]],
+                    samplePointValues[f0[2]],
+                    samplePointValues[oppositeI],
+
+                    triValues
+                );
+            }
+            else
+            {
+                forAll(cFaces, cFaceI)
+                {
+                    label faceI = cFaces[cFaceI];
+                    const face& f = mesh_.faces()[faceI];
+
+                    for(label fp = 1; fp < f.size() - 1; fp++)
+                    {
+                        vertexInterp
+                        (
+                            iso_,
+                            pointValues_[f[0]],
+                            pointValues_[f[fp]],
+                            pointValues_[f[f.fcIndex(fp)]],
+                            cellValues_[cellI],
+
+                            samplePointValues[f[0]],
+                            samplePointValues[f[fp]],
+                            samplePointValues[f[f.fcIndex(fp)]],
+                            sampleCellValues[cellI],
+
+                            triValues
+                        );
+                    }
+                }
+            }
+        }
+    }
+
+    // One value per point
+    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    Field<Type>& values = tvalues();
+
+    forAll(triValues, i)
+    {
+        values[triPointMergeMap_[i]] = triValues[i];
+    }
+
+    return tvalues;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
new file mode 100644
index 0000000000..53344e6282
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "sampledIsoSurface.H"
+#include "dictionary.H"
+#include "volFields.H"
+#include "volPointInterpolation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sampledIsoSurface, 0);
+    addNamedToRunTimeSelectionTable(sampledSurface, sampledIsoSurface, word, isoSurface);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurface::sampledIsoSurface
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    isoField_(dict.lookup("isoField")),
+    isoVal_(readScalar(dict.lookup("isoValue"))),
+    zoneName_(word::null),
+    facesPtr_(NULL),
+    storedTimeIndex_(-1),
+    meshCells_(0),
+    triPointMergeMap_(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::sampledIsoSurface::~sampledIsoSurface()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::sampledIsoSurface::correct(const bool meshChanged)
+{
+    // Only change of mesh changes plane - zone restriction gets lost
+    if (meshChanged)
+    {
+        facesPtr_.clear();
+    }
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurface::sample
+(
+    const volScalarField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurface::sample
+(
+    const volVectorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurface::sample
+(
+    const volSphericalTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurface::sample
+(
+    const volSymmTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurface::sample
+(
+    const volTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::sampledIsoSurface::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::sampledIsoSurface::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::sampledIsoSurface::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::sampledIsoSurface::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::sampledIsoSurface::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+void Foam::sampledIsoSurface::print(Ostream& os) const
+{
+    os  << "sampledIsoSurface: " << name() << " :"
+        << "  field:" << isoField_
+        << "  value:" << isoVal_
+        << "  faces:" << faces().size()
+        << "  points:" << points().size();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
new file mode 100644
index 0000000000..57445082f6
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -0,0 +1,232 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::sampledIsoSurface
+
+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
+    sampledIsoSurface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledIsoSurface_H
+#define sampledIsoSurface_H
+
+#include "sampledSurface.H"
+#include "triSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class sampledIsoSurface Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledIsoSurface
+:
+    public sampledSurface,
+    public triSurface
+{
+    // Private data
+
+        //- Field to get isoSurface of
+        word isoField_;
+
+        //- iso value
+        scalar isoVal_;
+
+        //- 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_;
+
+            //- For every unmerged triangle point the point in the triSurface
+            mutable labelList triPointMergeMap_;
+
+
+    // Private Member Functions
+
+        //- 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("sampledIsoSurface");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        sampledIsoSurface
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~sampledIsoSurface();
+
+
+    // 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 "sampledIsoSurfaceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
new file mode 100644
index 0000000000..755fc95939
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "sampledIsoSurface.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::sampledIsoSurface::sampleField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+) const
+{
+    const fvMesh& fvm = vField.mesh();
+
+    if (fvm.time().timeIndex() != storedTimeIndex_)
+    {
+        storedTimeIndex_ = fvm.time().timeIndex();
+
+        //- Clear any stored topo
+        facesPtr_.clear();
+
+        const volScalarField& cellFld =
+            fvm.lookupObject<volScalarField>(isoField_);
+
+        tmp<pointScalarField> pointFld
+        (
+            volPointInterpolation::New(fvm).interpolate(cellFld)
+        );
+
+        const isoSurface iso
+        (
+            fvm,
+            cellFld.internalField(),
+            pointFld().internalField(),
+            isoVal_
+        );
+
+        const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
+        meshCells_ = iso.meshCells();
+        triPointMergeMap_ = iso.triPointMergeMap();
+    }
+
+    return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
+}
+
+
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::sampledIsoSurface::interpolateField
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+    if (fvm.time().timeIndex() != storedTimeIndex_)
+    {
+        //- Clear any stored topo
+        facesPtr_.clear();
+
+        storedTimeIndex_ = fvm.time().timeIndex();
+
+        const volScalarField& cellFld =
+            fvm.lookupObject<volScalarField>(isoField_);
+
+        tmp<pointScalarField> pointFld
+        (
+            volPointInterpolation::New(fvm).interpolate(cellFld)
+        );
+
+        const isoSurface iso
+        (
+            fvm,
+            cellFld.internalField(),
+            pointFld().internalField(),
+            isoVal_
+        );
+
+        const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
+        meshCells_ = iso.meshCells();
+        triPointMergeMap_ = iso.triPointMergeMap();
+    }
+
+
+    // 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;
+}
+
+
+// ************************************************************************* //