Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

applications/utilities/mesh/generation/blockMesh/blockPoints.C

@@ -41,23 +41,23 @@ namespace Foam
 void block::blockPoints()
 {
     // set local variables for mesh specification
-    label ni = blockDef_.n().x();
-    label nj = blockDef_.n().y();
-    label nk = blockDef_.n().z();
+    const label ni = blockDef_.n().x();
+    const label nj = blockDef_.n().y();
+    const label nk = blockDef_.n().z();
 
-    point start = blockDef_.points()[blockDef_.blockShape()[0]];
-    point xEnd = blockDef_.points()[blockDef_.blockShape()[1]];
-    point xyEnd = blockDef_.points()[blockDef_.blockShape()[2]];
-    point yEnd = blockDef_.points()[blockDef_.blockShape()[3]];
+    const point p000 = blockDef_.points()[blockDef_.blockShape()[0]];
+    const point p100 = blockDef_.points()[blockDef_.blockShape()[1]];
+    const point p110 = blockDef_.points()[blockDef_.blockShape()[2]];
+    const point p010 = blockDef_.points()[blockDef_.blockShape()[3]];
 
-    point zEnd = blockDef_.points()[blockDef_.blockShape()[4]];
-    point xzEnd = blockDef_.points()[blockDef_.blockShape()[5]];
-    point xyzEnd = blockDef_.points()[blockDef_.blockShape()[6]];
-    point yzEnd = blockDef_.points()[blockDef_.blockShape()[7]];
+    const point p001 = blockDef_.points()[blockDef_.blockShape()[4]];
+    const point p101 = blockDef_.points()[blockDef_.blockShape()[5]];
+    const point p111 = blockDef_.points()[blockDef_.blockShape()[6]];
+    const point p011 = blockDef_.points()[blockDef_.blockShape()[7]];
 
-    // set reference to the list of edge point and weighting factors
-    const List<List<point> >& edgePoints = blockDef_.blockEdgePoints();
-    const scalarListList& edgeWeights = blockDef_.blockEdgeWeights();
+    // list of edge point and weighting factors
+    const List<List<point> >& p = blockDef_.blockEdgePoints();
+    const scalarListList& w = blockDef_.blockEdgeWeights();
 
     // generate vertices
 
@@ -69,193 +69,138 @@ void block::blockPoints()
             {
                 label vertexNo = vtxLabel(i, j, k);
 
-                vector edgex1 = start*(1.0 - edgeWeights[0][i])
-                      + xEnd*edgeWeights[0][i];
+                // points on edges
+                vector edgex1 = p000 + (p100 - p000)*w[0][i];
+                vector edgex2 = p010 + (p110 - p010)*w[1][i];
+                vector edgex3 = p011 + (p111 - p011)*w[2][i];
+                vector edgex4 = p001 + (p101 - p001)*w[3][i];
 
-                vector edgex2 = yEnd*(1.0 - edgeWeights[1][i])
-                      + xyEnd*edgeWeights[1][i];
+                vector edgey1 = p000 + (p010 - p000)*w[4][j];
+                vector edgey2 = p100 + (p110 - p100)*w[5][j];
+                vector edgey3 = p101 + (p111 - p101)*w[6][j];
+                vector edgey4 = p001 + (p011 - p001)*w[7][j];
 
-                vector edgex3 = yzEnd*(1.0 - edgeWeights[2][i])
-                      + xyzEnd*edgeWeights[2][i];
-
-                vector edgex4 = zEnd*(1.0 - edgeWeights[3][i])
-                      + xzEnd*edgeWeights[3][i];
-
-
-
-                vector edgey1 = start*(1.0 - edgeWeights[4][j])
-                      + yEnd*edgeWeights[4][j];
-
-                vector edgey2 = xEnd*(1.0 - edgeWeights[5][j])
-                      + xyEnd*edgeWeights[5][j];
-
-                vector edgey3 = xzEnd*(1.0 - edgeWeights[6][j])
-                      + xyzEnd*edgeWeights[6][j];
-
-                vector edgey4 = zEnd*(1.0 - edgeWeights[7][j])
-                      + yzEnd*edgeWeights[7][j];
-
-
-
-                vector edgez1 = start*(1.0 - edgeWeights[8][k])
-                      + zEnd*edgeWeights[8][k];
-
-                vector edgez2 = xEnd*(1.0 - edgeWeights[9][k])
-                      + xzEnd*edgeWeights[9][k];
-
-                vector edgez3 = xyEnd*(1.0 - edgeWeights[10][k])
-                      + xyzEnd*edgeWeights[10][k];
-
-                vector edgez4 = yEnd*(1.0 - edgeWeights[11][k])
-                      + yzEnd*edgeWeights[11][k];
+                vector edgez1 = p000 + (p001 - p000)*w[8][k];
+                vector edgez2 = p100 + (p101 - p100)*w[9][k];
+                vector edgez3 = p110 + (p111 - p110)*w[10][k];
+                vector edgez4 = p010 + (p011 - p010)*w[11][k];
 
                 // calculate the importance factors for all edges
                 // x - direction
                 scalar impx1 =
                 (
-                    (1.0 - edgeWeights[0][i])
-                        *(1.0 - edgeWeights[4][j])
-                        *(1.0 - edgeWeights[8][k])
-                  + edgeWeights[0][i]
-                        *(1.0 - edgeWeights[5][j])
-                        *(1.0 - edgeWeights[9][k])
+                    (1.0 - w[0][i])*(1.0 - w[4][j])*(1.0 - w[8][k])
+                  + w[0][i]*(1.0 - w[5][j])*(1.0 - w[9][k])
                 );
 
                 scalar impx2 =
                 (
-                    (1.0 - edgeWeights[1][i])
-                        *edgeWeights[4][j]
-                        *(1.0 - edgeWeights[11][k])
-                  + edgeWeights[1][i]
-                        *edgeWeights[5][j]
-                        *(1.0 - edgeWeights[10][k])
+                    (1.0 - w[1][i])*w[4][j]*(1.0 - w[11][k])
+                  + w[1][i]*w[5][j]*(1.0 - w[10][k])
                 );
 
-               scalar impx3 =
-               (
-                   (1.0 - edgeWeights[2][i])
-                       *edgeWeights[7][j]
-                       *edgeWeights[11][k]
-                 + edgeWeights[2][i]
-                       *edgeWeights[6][j]
-                       *edgeWeights[10][k]
-               );
+                scalar impx3 =
+                (
+                     (1.0 - w[2][i])*w[7][j]*w[11][k]
+                   + w[2][i]*w[6][j]*w[10][k]
+                );
 
 
-               scalar impx4 =
-               (
-                   (1.0 - edgeWeights[3][i])
-                       *(1.0 - edgeWeights[7][j])
-                       *edgeWeights[8][k]
-                 + edgeWeights[3][i]
-                       *(1.0 - edgeWeights[6][j])
-                       *edgeWeights[9][k]
-               );
+                scalar impx4 =
+                (
+                    (1.0 - w[3][i])*(1.0 - w[7][j])*w[8][k]
+                  + w[3][i]*(1.0 - w[6][j])*w[9][k]
+                );
 
+                scalar magImpx = impx1 + impx2 + impx3 + impx4;
+                impx1 /= magImpx;
+                impx2 /= magImpx;
+                impx3 /= magImpx;
+                impx4 /= magImpx;
 
 
                 // y - direction
                 scalar impy1 =
                 (
-                    (1.0 - edgeWeights[4][j])
-                        *(1.0 - edgeWeights[0][i])
-                        *(1.0 - edgeWeights[8][k])
-                  + edgeWeights[4][j]
-                        *(1.0 - edgeWeights[1][i])
-                        *(1.0 - edgeWeights[11][k])
+                    (1.0 - w[4][j])*(1.0 - w[0][i])*(1.0 - w[8][k])
+                  + w[4][j]*(1.0 - w[1][i])*(1.0 - w[11][k])
                 );
 
                 scalar impy2 =
                 (
-                    (1.0 - edgeWeights[5][j])
-                        *edgeWeights[0][i]
-                        *(1.0 - edgeWeights[9][k])
-                  + edgeWeights[5][j]
-                        *edgeWeights[1][i]
-                        *(1.0 - edgeWeights[10][k])
+                    (1.0 - w[5][j])*w[0][i]*(1.0 - w[9][k])
+                  + w[5][j]*w[1][i]*(1.0 - w[10][k])
                 );
 
                 scalar impy3 =
                 (
-                    (1.0 - edgeWeights[6][j])
-                        *edgeWeights[3][i]
-                        *edgeWeights[9][k]
-                  + edgeWeights[6][j]
-                        *edgeWeights[2][i]
-                        *edgeWeights[10][k]
+                    (1.0 - w[6][j])*w[3][i]*w[9][k]
+                  + w[6][j]*w[2][i]*w[10][k]
                 );
 
                 scalar impy4 =
                 (
-                    (1.0 - edgeWeights[7][j])
-                        *(1.0 - edgeWeights[3][i])
-                        *edgeWeights[8][k]
-                  + edgeWeights[7][j]
-                        *(1.0 - edgeWeights[2][i])
-                        *edgeWeights[11][k]
+                    (1.0 - w[7][j])*(1.0 - w[3][i])*w[8][k]
+                  + w[7][j]*(1.0 - w[2][i])*w[11][k]
                 );
 
+                scalar magImpy = impy1 + impy2 + impy3 + impy4;
+                impy1 /= magImpy;
+                impy2 /= magImpy;
+                impy3 /= magImpy;
+                impy4 /= magImpy;
 
 
                 // z - direction
                 scalar impz1 =
                 (
-                    (1.0 - edgeWeights[8][k])
-                        *(1.0 - edgeWeights[0][i])
-                        *(1.0 - edgeWeights[4][j])
-                  + edgeWeights[8][k]
-                        *(1.0 - edgeWeights[3][i])
-                        *(1.0 - edgeWeights[7][j])
+                    (1.0 - w[8][k])*(1.0 - w[0][i])*(1.0 - w[4][j])
+                  + w[8][k]*(1.0 - w[3][i])*(1.0 - w[7][j])
                 );
 
                 scalar impz2 =
                 (
-                    (1.0 - edgeWeights[9][k])
-                        *edgeWeights[0][i]
-                        *(1.0 - edgeWeights[5][j])
-                  + edgeWeights[9][k]
-                        *edgeWeights[3][i]
-                        *(1.0 - edgeWeights[6][j])
+                    (1.0 - w[9][k])*w[0][i]*(1.0 - w[5][j])
+                  + w[9][k]*w[3][i]*(1.0 - w[6][j])
                 );
 
                 scalar impz3 =
                 (
-                    (1.0 - edgeWeights[10][k])
-                        *edgeWeights[1][i]
-                        *edgeWeights[5][j]
-                  + edgeWeights[10][k]
-                        *edgeWeights[2][i]
-                        *edgeWeights[6][j]
+                    (1.0 - w[10][k])*w[1][i]*w[5][j]
+                  + w[10][k]*w[2][i]*w[6][j]
                 );
 
                 scalar impz4 =
                 (
-                    (1.0 - edgeWeights[11][k])
-                        *(1.0 - edgeWeights[1][i])
-                        *edgeWeights[4][j]
-                  + edgeWeights[11][k]
-                        *(1.0 - edgeWeights[2][i])
-                        *edgeWeights[7][j]
+                    (1.0 - w[11][k])*(1.0 - w[1][i])*w[4][j]
+                  + w[11][k]*(1.0 - w[2][i])*w[7][j]
                 );
 
+                scalar magImpz = impz1 + impz2 + impz3 + impz4;
+                impz1 /= magImpz;
+                impz2 /= magImpz;
+                impz3 /= magImpz;
+                impz4 /= magImpz;
+
                 // calculate the correction vectors
-                vector corx1 = impx1*(edgePoints[0][i] - edgex1);
-                vector corx2 = impx2*(edgePoints[1][i] - edgex2);
-                vector corx3 = impx3*(edgePoints[2][i] - edgex3);
-                vector corx4 = impx4*(edgePoints[3][i] - edgex4);
+                vector corx1 = impx1*(p[0][i] - edgex1);
+                vector corx2 = impx2*(p[1][i] - edgex2);
+                vector corx3 = impx3*(p[2][i] - edgex3);
+                vector corx4 = impx4*(p[3][i] - edgex4);
 
-                vector cory1 = impy1*(edgePoints[4][j] - edgey1);
-                vector cory2 = impy2*(edgePoints[5][j] - edgey2);
-                vector cory3 = impy3*(edgePoints[6][j] - edgey3);
-                vector cory4 = impy4*(edgePoints[7][j] - edgey4);
+                vector cory1 = impy1*(p[4][j] - edgey1);
+                vector cory2 = impy2*(p[5][j] - edgey2);
+                vector cory3 = impy3*(p[6][j] - edgey3);
+                vector cory4 = impy4*(p[7][j] - edgey4);
 
-                vector corz1 = impz1*(edgePoints[8][k] - edgez1);
-                vector corz2 = impz2*(edgePoints[9][k] - edgez2);
-                vector corz3 = impz3*(edgePoints[10][k] - edgez3);
-                vector corz4 = impz4*(edgePoints[11][k] - edgez4);
+                vector corz1 = impz1*(p[8][k] - edgez1);
+                vector corz2 = impz2*(p[9][k] - edgez2);
+                vector corz3 = impz3*(p[10][k] - edgez3);
+                vector corz4 = impz4*(p[11][k] - edgez4);
 
 
                 // multiply by the importance factor
+
                 // x - direction
                 edgex1 *= impx1;
                 edgex2 *= impx2;
@@ -285,7 +230,6 @@ void block::blockPoints()
                 vertices_[vertexNo] += corx1 + corx2 + corx3 + corx4;
                 vertices_[vertexNo] += cory1 + cory2 + cory3 + cory4;
                 vertices_[vertexNo] += corz1 + corz2 + corz3 + corz4;
-
             }
         }
     }

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -40,7 +40,7 @@ surfaceFormat vtk;
 // 1] vertex values determined from neighbouring cell-centre values
 // 2] face values determined using the current face interpolation scheme
 //    for the field (linear, gamma, etc.)
-interpolationScheme cellPointFace;
+interpolationScheme cellPoint;
 
 // Fields to sample.
 fields
@@ -155,21 +155,25 @@ surfaces
         // Optional: whether to leave as faces (=default) or triangulate
     }
 
-/* not yet (re)implemented --
+    interpolatedIso
+    {
+        // Iso surface for interpolated values only
+        type            isoSurface;
+        isoField        rho;
+        isoValue        0.5;
+        interpolate     true;
+        //regularise      false;    //optional: do not simplify
+    }
     constantIso
     {
-        name            iso;
-        field           rho;
-        value           0.5;
+        // 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
     }
-    someIso
-    {
-        type            iso;
-        field           rho;
-        value           0.5;
-        interpolate     true;
-    }
- */
 );
 
 

applications/utilities/surface/surfaceCheck/surfaceCheck.C

@@ -36,7 +36,7 @@ License
 using namespace Foam;
 
 // Does face use valid vertices?
-bool validTri(const triSurface& surf, const label faceI)
+bool validTri(const bool verbose, const triSurface& surf, const label faceI)
 {
     // Simple check on indices ok.
 
@@ -49,20 +49,21 @@ bool validTri(const triSurface& surf, const label faceI)
      || (f[2] < 0) || (f[2] >= surf.points().size())
     )
     {
-        //WarningIn("validTri(const triSurface&, const label)")
-        //    << "triangle " << faceI << " vertices " << f
-        //    << " uses point indices outside point range 0.."
-        //    << surf.points().size()-1 << endl;
+        WarningIn("validTri(const triSurface&, const label)")
+            << "triangle " << faceI << " vertices " << f
+            << " uses point indices outside point range 0.."
+            << surf.points().size()-1 << endl;
 
         return false;
     }
 
     if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
     {
-        //WarningIn("validTri(const triSurface&, const label)")
-        //    << "triangle " << faceI
-        //    << " uses non-unique vertices " << f
-        //    << endl;
+        WarningIn("validTri(const triSurface&, const label)")
+            << "triangle " << faceI
+            << " uses non-unique vertices " << f
+            << " coords:" << f.points(surf.points())
+            << endl;
         return false;
     }
 
@@ -91,11 +92,12 @@ bool validTri(const triSurface& surf, const label faceI)
          && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
         )
         {
-            //WarningIn("validTri(const triSurface&, const label)")
-            //    << "triangle " << faceI << " vertices " << f
-            //    << " has the same vertices as triangle " << nbrFaceI
-            //    << " vertices " << nbrF
-            //    << endl;
+            WarningIn("validTri(const triSurface&, const label)")
+                << "triangle " << faceI << " vertices " << f
+                << " has the same vertices as triangle " << nbrFaceI
+                << " vertices " << nbrF
+                << " coords:" << f.points(surf.points())
+                << endl;
 
             return false;
         }
@@ -170,9 +172,11 @@ int main(int argc, char *argv[])
     argList::validArgs.clear();
     argList::validArgs.append("surface file");
     argList::validOptions.insert("noSelfIntersection", "");
+    argList::validOptions.insert("verbose", "");
     argList args(argc, argv);
 
     bool checkSelfIntersection = !args.options().found("noSelfIntersection");
+    bool verbose = args.options().found("verbose");
 
     fileName surfFileName(args.additionalArgs()[0]);
     Pout<< "Reading surface from " << surfFileName << " ..." << nl << endl;
@@ -232,7 +236,7 @@ int main(int argc, char *argv[])
 
         forAll(surf, faceI)
         {
-            if (!validTri(surf, faceI))
+            if (!validTri(verbose, surf, faceI))
             {
                 illegalFaces.append(faceI);
             }

applications/utilities/surface/surfaceSubset/surfaceSubsetDict

@@ -42,4 +42,7 @@ surface
 // Extend selection with edge neighbours
 addFaceNeighbours no;
 
+// Invert selection
+invertSelection false;
+
 // ************************************************************************* //

src/OpenFOAM/containers/Lists/SortableList/SortableList.C

@@ -32,7 +32,7 @@ License
 
 // Construct from List
 template <class Type>
-Foam::SortableList<Type>::SortableList(const List<Type>& values)
+Foam::SortableList<Type>::SortableList(const UList<Type>& values)
 :
     List<Type>(values),
     indices_(values.size())

src/OpenFOAM/containers/Lists/SortableList/SortableList.H

@@ -88,7 +88,7 @@ public:
 
         //- Construct from List, sorting the elements.
         //  Starts with indices set to index in argument
-        explicit SortableList(const List<Type>&);
+        explicit SortableList(const UList<Type>&);
 
         //- Construct from tranferred List, sorting the elements.
         //  Starts with indices set to index in argument

src/OpenFOAM/fields/Fields/Field/Field.C

@@ -575,20 +575,6 @@ void Field<Type>::replace
 }
 
 
-template<class Type>
-void Field<Type>::transfer(Field<Type>& f)
-{
-    List<Type>::transfer(f);
-}
-
-
-template<class Type>
-void Field<Type>::transfer(List<Type>& lst)
-{
-    List<Type>::transfer(lst);
-}
-
-
 template<class Type>
 tmp<Field<Type> > Field<Type>::T() const
 {

src/OpenFOAM/fields/Fields/Field/Field.H

@@ -297,12 +297,6 @@ public:
         //- Replace a component field of the field
         void replace(const direction, const cmptType&);
 
-        //- Transfer the contents of the argument Field into this Field
-        void transfer(Field<Type>&);
-
-        //- Transfer the contents of the argument List into this Field
-        void transfer(List<Type>&);
-
         //- Return the field transpose (only defined for second rank tensors)
         tmp<Field<Type> > T() const;
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointFaces.C

@@ -38,6 +38,11 @@ const labelListList& primitiveMesh::pointFaces() const
 {
     if (!pfPtr_)
     {
+        if (debug)
+        {
+            Pout<< "primitiveMesh::pointFaces() : "
+                << "calculating pointFaces" << endl;
+        }
         // Invert faces()
         pfPtr_ = new labelListList(nPoints());
         invertManyToMany(nPoints(), faces(), *pfPtr_);

src/Pstream/mpi/IPread.C

@@ -193,15 +193,13 @@ void Foam::IPstream::waitRequests()
 {
     if (IPstream_outstandingRequests_.size() > 0)
     {
-        List<MPI_Status> status(IPstream_outstandingRequests_.size());
-
         if
         (
             MPI_Waitall
             (
                 IPstream_outstandingRequests_.size(),
                 IPstream_outstandingRequests_.begin(),
-                status.begin()
+                MPI_STATUSES_IGNORE
             )
         )
         {
@@ -231,9 +229,7 @@ bool Foam::IPstream::finishedRequest(const label i)
     }
 
     int flag;
-    MPI_Status status;
-
-    MPI_Test(&IPstream_outstandingRequests_[i], &flag, &status);
+    MPI_Test(&IPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE);
 
     return flag != 0;
 }

src/Pstream/mpi/OPwrite.C

@@ -131,15 +131,13 @@ void Foam::OPstream::waitRequests()
 {
     if (OPstream_outstandingRequests_.size() > 0)
     {
-        List<MPI_Status> status(OPstream_outstandingRequests_.size());
-
         if
         (
             MPI_Waitall
             (
                 OPstream_outstandingRequests_.size(),
                 OPstream_outstandingRequests_.begin(),
-                status.begin()
+                MPI_STATUSES_IGNORE
             )
         )
         {
@@ -169,9 +167,7 @@ bool Foam::OPstream::finishedRequest(const label i)
     }
 
     int flag;
-    MPI_Status status;
-
-    MPI_Test(&OPstream_outstandingRequests_[i], &flag, &status);
+    MPI_Test(&OPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE);
 
     return flag != 0;
 }

src/Pstream/mpi/Pstream.C

@@ -157,8 +157,6 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
 
     if (Pstream::nProcs() <= Pstream::nProcsSimpleSum)
     {
-        MPI_Status status;
-
         if (Pstream::master())
         {
             for
@@ -180,7 +178,7 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
                         Pstream::procID(slave),
                         Pstream::msgType(),
                         MPI_COMM_WORLD,
-                        &status
+                        MPI_STATUS_IGNORE
                     )
                 )
                 {
@@ -260,7 +258,7 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
                     Pstream::procID(Pstream::masterNo()),
                     Pstream::msgType(),
                     MPI_COMM_WORLD,
-                    &status
+                    MPI_STATUS_IGNORE
                 )
             )
             {
@@ -279,8 +277,6 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
         Value = sum;
 
         /*
-        MPI_Status status;
-
         int myProcNo = Pstream::myProcNo();
         int nProcs = Pstream::nProcs();
 
@@ -314,7 +310,7 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
                         Pstream::procID(childProcId),
                         Pstream::msgType(),
                         MPI_COMM_WORLD,
-                        &status
+                        MPI_STATUS_IGNORE
                     )
                 )
                 {
@@ -370,7 +366,7 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
                     Pstream::procID(parentId),
                     Pstream::msgType(),
                     MPI_COMM_WORLD,
-                    &status
+                    MPI_STATUS_IGNORE
                 )
             )
             {

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -2875,6 +2875,8 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
     }
 
     // 2. Extend to 2:1. I don't understand yet why this is not done
+    // 2. Extend to 2:1. For non-cube cells the scalar distance does not work
+    // so make sure it at least provides 2:1.
     PackedList<1> refineCell(mesh_.nCells(), 0);
     forAll(allCellInfo, cellI)
     {
@@ -2887,6 +2889,25 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
     }
     faceConsistentRefinement(true, refineCell);
 
+    while (true)
+    {
+        label nChanged = faceConsistentRefinement(true, refineCell);
+
+        reduce(nChanged, sumOp<label>());
+
+        if (debug)
+        {
+            Pout<< "hexRef8::consistentSlowRefinement2 : Changed " << nChanged
+                << " refinement levels due to 2:1 conflicts."
+                << endl;
+        }
+
+        if (nChanged == 0)
+        {
+            break;
+        }
+    }
+
     // 3. Convert back to labelList.
     label nRefined = 0;
 

src/finiteVolume/Make/files

@@ -297,15 +297,4 @@ cfdTools/general/SRF/SRFModel/SRFModel/newSRFModel.C
 cfdTools/general/SRF/SRFModel/rpm/rpm.C
 cfdTools/general/SRF/derivedFvPatchFields/SRFVelocityFvPatchVectorField/SRFVelocityFvPatchVectorField.C
 
-fvMeshCutSurface = fvMesh/fvMeshCutSurface
-
-meshCut = $(fvMeshCutSurface)/meshCut
-$(meshCut)/meshCutSurface.C
-$(meshCut)/cellAddressing.C
-
-edgeCuts = $(fvMeshCutSurface)/edgeCuts
-$(edgeCuts)/meshEdgeCuts.C
-$(edgeCuts)/faceDecompIsoSurfaceCuts.C
-$(edgeCuts)/cellDecompIsoSurfaceCuts.C
-
 LIB = $(FOAM_LIBBIN)/libfiniteVolume

src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C

@@ -33,15 +33,68 @@ void Foam::setRefCell
     const volScalarField& field,
     const dictionary& dict,
     label& refCelli,
-    scalar& refValue
+    scalar& refValue,
+    bool forceReference
 )
 {
-    if (field.needReference())
+    if (field.needReference() || forceReference)
     {
         word refCellName = field.name() + "RefCell";
+        word refPointName = field.name() + "RefPoint";
+
         word refValueName = field.name() + "RefValue";
 
-        refCelli = readLabel(dict.lookup(refCellName));
+        if (dict.found(refCellName))
+        {
+            if (Pstream::master())
+            {
+                refCelli = readLabel(dict.lookup(refCellName));
+            }
+            else
+            {
+                refCelli = -1;
+            }
+        }
+        else if (dict.found(refPointName))
+        {
+            point refPointi(dict.lookup(refPointName));
+            refCelli = field.mesh().findCell(refPointi);
+            label hasRef = (refCelli >= 0 ? 1 : 0);
+            label sumHasRef = returnReduce<label>(hasRef, sumOp<label>());
+            if (sumHasRef != 1)
+            {
+                FatalErrorIn
+                (
+                    "void Foam::setRefCell"
+                     "("
+                     "    const volScalarField&,"
+                     "    const dictionary&,"
+                     "    label& scalar&,"
+                     "    bool"
+                     ")"
+                )
+                  << "Unable to set reference cell for field " << field.name()
+                  << nl << "    Reference point " << refPointName
+                  << " found on multiple domains" << nl << abort(FatalError);
+            }
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "void Foam::setRefCell"
+                 "("
+                 "    const volScalarField&,"
+                 "    const dictionary&,"
+                 "    label& scalar&,"
+                 "    bool"
+                 ")"
+            )
+              << "Unable to set reference cell for field" << field.name() << nl
+              << "    Please supply either " << refCellName
+              << " or " << refPointName << nl << abort(FatalError);
+        }
+
         refValue = readScalar(dict.lookup(refValueName));
     }
 }

src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H

@@ -52,7 +52,8 @@ void setRefCell
     const volScalarField& field,
     const dictionary& dict,
     label& refCelli,
-    scalar& refValue
+    scalar& refValue,
+    bool forceReference = false
 );
 
 }

src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C

@@ -234,6 +234,16 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs()
                  (
                     mpp.samplePatch()
                  );
+            if (patchID < 0)
+            {
+                FatalErrorIn
+                (
+                    "void directMappedFixedValueFvPatchField<Type>::"
+                    "updateCoeffs()"
+                )<< "Unable to find sample patch " << mpp.samplePatch()
+                 << " for patch " << this->patch().name() << nl
+                 << abort(FatalError);
+            }
             typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
             const word& fieldName = this->dimensionedInternalField().name();
             const fieldType& sendField =

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C

@@ -483,7 +483,7 @@ void Foam::fvMatrix<Type>::setReference
 {
     if (psi_.needReference() || forceReference)
     {
-        if (Pstream::master())
+        if (cell >= 0)
         {
             source()[cell] += diag()[cell]*value;
             diag()[cell] += diag()[cell];

src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompCuts.H

@@ -1,161 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::cellDecompCuts
-
-Description
-    Container for cuts of edges of (implicit) tet decomposition.
-    Used to collect data for meshCut.
-
-    As much as possible, cuts are defined using mesh information:
-      - cut (exactly) through mesh vertex
-      - cut (exactly) through cell centre
-
-      - cut through mesh edge. Both edge label and position on edge given.
-
-      - cut through tet pyramidEdge (edge between vertex and cell centre).
-        Edge and position on edge given.
-
-      - cut through diagonalEdge (edge between vertices of a face)
-        Edge and position on edge given.
-
-SourceFiles
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef cellDecompCuts_H
-#define cellDecompCuts_H
-
-#include "meshEdgeCuts.H"
-#include "pyramidEdge.H"
-#include "diagonalEdge.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-
-/*---------------------------------------------------------------------------*\
-                           Class cellDecompCuts Declaration
-\*---------------------------------------------------------------------------*/
-
-class cellDecompCuts
-:
-    public meshEdgeCuts
-{
-
-protected:
-
-        labelList meshCellCentres_;
-
-        List<pyramidEdge> pyrEdges_;
-        scalarField pyrEdgeWeights_;
-
-        List<diagonalEdge> diagEdges_;
-        scalarField diagEdgeWeights_;
-
-    // Private Member Functions
-
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        cellDecompCuts
-        (
-            const primitiveMesh& mesh,
-            const labelList& cells,
-
-            const labelList& meshVerts,
-            const labelList& meshCellCentres,
-
-            const labelList& meshEdges,
-            const scalarField& meshEdgeWeights,
-
-            const List<pyramidEdge>& pyrEdges,
-            const scalarField& pyrEdgeWeights,
-
-            const List<diagonalEdge>& diagEdges,
-            const scalarField& diagEdgeWeights
-        )
-        :
-            meshEdgeCuts(mesh, cells, meshVerts, meshEdges, meshEdgeWeights),
-            meshCellCentres_(meshCellCentres),
-            pyrEdges_(pyrEdges),
-            pyrEdgeWeights_(pyrEdgeWeights),
-            diagEdges_(diagEdges),
-            diagEdgeWeights_(diagEdgeWeights)
-        {}
-
-
-    // Member Functions
-
-        const labelList& meshCellCentres() const
-        {
-            return meshCellCentres_;
-        }
-
-        const List<pyramidEdge>& pyrEdges() const
-        {
-            return pyrEdges_;
-        }
-
-        const scalarField& pyrEdgeWeights() const
-        {
-            return pyrEdgeWeights_;
-        }
-
-        const List<diagonalEdge>& diagEdges() const
-        {
-            return diagEdges_;
-        }
-
-        const scalarField& diagEdgeWeights() const
-        {
-            return diagEdgeWeights_;
-        }
-
-        label size() const
-        {
-            return
-                meshEdgeCuts::size() + meshCellCentres_.size()
-              + pyrEdges_.size() + diagEdges_.size();
-        }
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompIsoSurfaceCuts.C

@@ -1,267 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
-
-Description
-
-\*---------------------------------------------------------------------------*/
-
-#include "cellDecompIsoSurfaceCuts.H"
-#include "volPointInterpolation.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::cellDecompIsoSurfaceCuts::constructEdgeCuts
-(
-    const volScalarField& volField,
-    const pointScalarField& ptField,
-    const scalar isoVal,
-    const scalar tol
-)
-{
-    const primitiveMesh& mesh = volField.mesh();
-
-    // Intermediate storage for labels of cut edges and cut points
-    label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
-
-    labelHashSet cutCells(nBFaces);
-    labelHashSet cutVerts(nBFaces);
-    labelHashSet cutCellCentres(nBFaces);
-
-    DynamicList<label> meshCutEdges(nBFaces);
-    DynamicList<scalar> meshEdgeWeights(nBFaces);
-
-    DynamicList<pyramidEdge> cutPyrEdges(nBFaces);
-    DynamicList<scalar> pyrEdgeWeights(nBFaces);
-
-    DynamicList<diagonalEdge> cutDiagEdges(nBFaces);
-    DynamicList<scalar> diagEdgeWeights(nBFaces);
-
-
-    // Walk over faces
-    forAll(mesh.faces(), faceI)
-    {
-        const face& f = mesh.faces()[faceI];
-
-        label ownerI = mesh.faceOwner()[faceI];
-        label neighbourI = -1;
-        if (mesh.isInternalFace(faceI))
-        {
-            neighbourI = mesh.faceNeighbour()[faceI];
-        }
-
-        scalar weight;
-
-        //
-        // Check face diagonal. Simple triangulation from f[0]
-        //
-
-        for(label fp = 2; fp < f.size() - 1; fp++)
-        {
-            if (crosses(isoVal, ptField[f[0]], ptField[f[fp]], weight))
-            {
-                mark(ownerI, cutCells);
-                if (neighbourI != -1)
-                {
-                    mark(neighbourI, cutCells);
-                }
-
-                if (weight < tol)
-                {
-                    mark(f[0], cutVerts);
-                }
-                else if (weight > 1-tol)
-                {
-                    mark(f[fp], cutVerts);
-                }
-                else
-                {
-                    cutDiagEdges.append(diagonalEdge(faceI, 0, fp));
-                    diagEdgeWeights.append(weight);
-                }
-            }
-        }
-    }
-    
-
-    // Walk over all mesh points
-    forAll(mesh.points(), pointI)
-    {
-        const labelList& myCells = mesh.pointCells()[pointI];
-
-        forAll(myCells, myCellI)
-        {
-            label cellI = myCells[myCellI];
-
-            //
-            // Check pyramid edge crossing
-            //
-
-            scalar weight;
-
-            if (crosses(isoVal, ptField[pointI], volField[cellI], weight))
-            {
-                mark(cellI, cutCells);
-
-                if (weight < tol)
-                {
-                    mark(pointI, cutVerts);
-                }
-                else if (weight > 1-tol)
-                {
-                    mark(cellI, cutCellCentres);
-                }
-                else
-                {
-                    cutPyrEdges.append(pyramidEdge(pointI, cellI));
-                    pyrEdgeWeights.append(weight);
-                }
-            }
-        }
-    }            
-
-
-    // Walk over all mesh edges
-    forAll(mesh.edges(), edgeI)
-    {
-        const edge& e = mesh.edges()[edgeI];
-
-        scalar weight;
-
-        if (crosses(isoVal, ptField[e.start()], ptField[e.end()], weight))
-        {
-            mark(mesh.edgeCells()[edgeI], cutCells);
-
-            if (weight < tol)
-            {
-                mark(e.start(), cutVerts);
-            }
-            else if (weight > 1-tol)
-            {
-                mark(e.end(), cutVerts);
-            }
-            else
-            {
-                meshCutEdges.append(edgeI);
-                meshEdgeWeights.append(weight);
-            }
-        }
-    }
-
-    meshCutEdges.shrink();
-    meshEdgeWeights.shrink();
-
-    cutPyrEdges.shrink();
-    pyrEdgeWeights.shrink();
-
-    cutDiagEdges.shrink();
-    diagEdgeWeights.shrink();
-
-
-    // Tranfer lists to cellDecompCuts
-
-    cells_ = cutCells.toc();
-
-    meshVerts_ = cutVerts.toc();
-    meshCellCentres_ = cutCellCentres.toc();
-
-    meshEdges_.transfer(meshCutEdges);
-    meshEdgeWeights_.transfer(meshEdgeWeights);
-
-    pyrEdges_.transfer(cutPyrEdges);
-    pyrEdgeWeights_.transfer(pyrEdgeWeights);
-
-    diagEdges_.transfer(cutDiagEdges);
-    diagEdgeWeights_.transfer(diagEdgeWeights);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::cellDecompIsoSurfaceCuts::cellDecompIsoSurfaceCuts
-(
-    const volScalarField& volField,
-    const pointScalarField& ptField,
-    const scalar isoVal,
-    const scalar tol
-)
-:
-    cellDecompCuts
-    (
-        volField.mesh(),
-        labelList(),
-
-        labelList(),            // mesh vertices
-        labelList(),            // mesh cell centres
-
-        labelList(),            // mesh edges
-        scalarField(),
-
-        List<pyramidEdge>(),    // pyramid edges
-        scalarField(),
-
-        List<diagonalEdge>(),   // face diagonal edges
-        scalarField()
-    )
-{
-    constructEdgeCuts(volField, ptField, isoVal, tol);
-}
-
-
-// Construct from interpolator and field (does interpolation)
-Foam::cellDecompIsoSurfaceCuts::cellDecompIsoSurfaceCuts
-(
-    const volScalarField& volField,
-    const volPointInterpolation& pInterp,
-    const scalar isoVal,
-    const scalar tol
-)
-:
-    cellDecompCuts
-    (
-        volField.mesh(),
-        labelList(),
-
-        labelList(),            // mesh vertices
-        labelList(),            // mesh cell centres
-
-        labelList(),            // mesh edges
-        scalarField(),
-
-        List<pyramidEdge>(),    // pyramid edges
-        scalarField(),
-
-        List<diagonalEdge>(),   // face diagonal edges
-        scalarField()
-    )
-{
-    // Get field on vertices
-    tmp<pointScalarField> tptField = pInterp.interpolate(volField);
-    const pointScalarField& ptField = tptField();
-
-
-    constructEdgeCuts(volField, ptField, isoVal, tol);
-}
-
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompIsoSurfaceCuts.H

@@ -1,101 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::cellDecompIsoSurfaceCuts
-
-Description
-    cellCuts resulting from isoSurface reconstruction on cellDecomp tet mesh
-
-SourceFiles
-    cellDecompIsoSurfaceCuts.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef cellDecompIsoSurfaceCuts_H
-#define cellDecompIsoSurfaceCuts_H
-
-#include "cellDecompCuts.H"
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class volPointInterpolation;
-
-/*---------------------------------------------------------------------------*\
-                           Class cellDecompIsoSurfaceCuts Declaration
-\*---------------------------------------------------------------------------*/
-
-class cellDecompIsoSurfaceCuts
-:
-    public cellDecompCuts
-{
-    // Private Member Functions
-
-        //- Do all the hard work: determine edge crossings
-        void constructEdgeCuts
-        (
-            const volScalarField&,
-            const pointScalarField&,
-            const scalar,
-            const scalar
-        );
-
-public:
-
-    // Constructors
-
-        //- Construct from field on cells and field on vertices
-        cellDecompIsoSurfaceCuts
-        (
-            const volScalarField&,
-            const pointScalarField&,
-            const scalar isoValue,
-            const scalar tol
-        );
-
-        //- Construct from interpolator and field (does interpolation)
-        cellDecompIsoSurfaceCuts
-        (
-            const volScalarField&,
-            const volPointInterpolation&,
-            const scalar isoValue,
-            const scalar tol
-        );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompCuts.H

@@ -1,200 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::faceDecompCuts
-
-Description
-    Container for cuts of edges of (implicit) tet decomposition. Used to
-    collect data for meshCut. As much as possible cuts are defined using
-    mesh information:
-
-    - cut (exactly) through mesh vertex
-    -   ,,              ,,       cell centre
-    -   ,,              ,,       face centre
-
-    - cut through mesh edge. Both edge label and position on edge given.
-
-    - cut through tet pyramidEdge (edge between vertex and cell centre).
-      Edge and position on edge given.
-    - cut through tet faceEdge (edge between vertex and face centre).
-      Edge and position on edge given.
-    - cut through tet centreEdge (edge between face centre and cell centre).
-      Edge and position on edge given.
-
-SourceFiles
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef faceDecompCuts_H
-#define faceDecompCuts_H
-
-#include "meshEdgeCuts.H"
-#include "pyramidEdge.H"
-#include "faceEdge.H"
-#include "centreEdge.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-
-/*---------------------------------------------------------------------------*\
-                           Class faceDecompCuts Declaration
-\*---------------------------------------------------------------------------*/
-
-class faceDecompCuts
-:
-    public meshEdgeCuts
-{
-
-protected:
-
-        //- List of faces whose centre is exactly cut
-        labelList meshFaceCentres_;
-
-        //- List of cells whose centre is exactly cut
-        labelList meshCellCentres_;
-
-        //- List of pyramid edge descriptions
-        List<pyramidEdge> pyrEdges_;
-
-        //- For each pyramid edge description the position of the cut
-        scalarField pyrEdgeWeights_;
-
-        //- List of face edge descriptions
-        List<faceEdge> faceEdges_;
-
-        //- For each face edge description the position of the cut
-        scalarField faceEdgeWeights_;
-
-        //- List of centre edge descriptions
-        List<centreEdge> centreEdges_;
-
-        //- For each centre edge description the position of the cut
-        scalarField centreEdgeWeights_;
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        faceDecompCuts
-        (
-            const primitiveMesh& mesh,
-            const labelList& cells,
-
-            const labelList& meshVerts,
-            const labelList& meshFaceCentres,
-            const labelList& meshCellCentres,
-
-            const labelList& meshEdges,
-            const scalarField& meshEdgeWeights,
-
-            const List<pyramidEdge>& pyrEdges,
-            const scalarField& pyrEdgeWeights,
-
-            const List<faceEdge>& faceEdges,
-            const scalarField& faceEdgeWeights,
-
-            const List<centreEdge>& centreEdges,
-            const scalarField& centreEdgeWeights
-        )
-        :
-            meshEdgeCuts(mesh, cells, meshVerts, meshEdges, meshEdgeWeights),
-            meshFaceCentres_(meshFaceCentres),
-            meshCellCentres_(meshCellCentres),
-            pyrEdges_(pyrEdges),
-            pyrEdgeWeights_(pyrEdgeWeights),
-            faceEdges_(faceEdges),
-            faceEdgeWeights_(faceEdgeWeights),
-            centreEdges_(centreEdges),
-            centreEdgeWeights_(centreEdgeWeights)
-        {}
-
-
-    // Member Functions
-
-        const labelList& meshFaceCentres() const
-        {
-            return meshFaceCentres_;
-        }
-
-        const labelList& meshCellCentres() const
-        {
-            return meshCellCentres_;
-        }
-
-        const List<pyramidEdge>& pyrEdges() const
-        {
-            return pyrEdges_;
-        }
-
-        const scalarField& pyrEdgeWeights() const
-        {
-            return pyrEdgeWeights_;
-        }
-
-        const List<centreEdge>& centreEdges() const
-        {
-            return centreEdges_;
-        }
-
-        const scalarField& centreEdgeWeights() const
-        {
-            return centreEdgeWeights_;
-        }
-
-        const List<faceEdge>& faceEdges() const
-        {
-            return faceEdges_;
-        }
-
-        const scalarField& faceEdgeWeights() const
-        {
-            return faceEdgeWeights_;
-        }
-
-        label size() const
-        {
-            return
-                meshEdgeCuts::size()
-              + meshFaceCentres_.size() + meshCellCentres_.size()
-              + pyrEdges_.size() + centreEdges_.size() + faceEdges_.size();       
-        }
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompIsoSurfaceCuts.C

@@ -1,342 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
-
-Description
-
-\*---------------------------------------------------------------------------*/
-
-#include "faceDecompIsoSurfaceCuts.H"
-#include "volPointInterpolation.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::faceDecompIsoSurfaceCuts::constructEdgeCuts
-(
-    const volScalarField& volField,
-    const pointScalarField& ptField,
-    const scalar isoVal,
-    const scalar tol
-)
-{
-    const primitiveMesh& mesh = volField.mesh();
-
-    // Intermediate storage for labels of cut edges and cut points
-    label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
-
-    labelHashSet cutCells(nBFaces);
-    labelHashSet cutVerts(nBFaces);
-    labelHashSet cutFaceCentres(nBFaces);
-    labelHashSet cutCellCentres(nBFaces);
-
-    DynamicList<label> meshCutEdges(nBFaces);
-    DynamicList<scalar> meshEdgeWeights(nBFaces);
-
-    DynamicList<pyramidEdge> cutPyrEdges(nBFaces);
-    DynamicList<scalar> pyrEdgeWeights(nBFaces);
-
-    DynamicList<faceEdge> cutFaceEdges(nBFaces);
-    DynamicList<scalar> faceEdgeWeights(nBFaces);
-
-    DynamicList<centreEdge> cutCentreEdges(nBFaces);
-    DynamicList<scalar> centreEdgeWeights(nBFaces);
-
-
-    // Walk over faces
-    forAll(mesh.faces(), faceI)
-    {
-        const face& f = mesh.faces()[faceI];
-
-        // Face centre value
-        scalar fcVal = 0;
-
-        forAll(f, fp)
-        {
-            fcVal += ptField[f[fp]];
-        }
-        fcVal /= f.size();
-
-
-        label ownerI = mesh.faceOwner()[faceI];
-        label neighbourI = -1;
-        if (mesh.isInternalFace(faceI))
-        {
-            neighbourI = mesh.faceNeighbour()[faceI];
-        }
-
-        scalar weight;
-
-        //
-        // Check faceCentre-cellCentre crossing for owner and neigbour
-        //
-
-        if (crosses(isoVal, fcVal, volField[ownerI], weight))
-        {
-            mark(ownerI, cutCells);
-
-            if (weight < tol)
-            {
-                mark(faceI, cutFaceCentres);
-            }
-            else if (weight > 1-tol)
-            {
-                mark(ownerI, cutCellCentres);
-            }
-            else
-            {
-                cutCentreEdges.append(centreEdge(faceI, true));
-                centreEdgeWeights.append(weight);
-            }
-        }
-
-        if (neighbourI != -1)
-        {
-            if (crosses(isoVal, fcVal, volField[neighbourI], weight))
-            {
-                mark(neighbourI, cutCells);
-
-                if (weight < tol)
-                {
-                    mark(faceI, cutFaceCentres);
-                }
-                else if (weight > 1-tol)
-                {
-                    mark(neighbourI, cutCellCentres);
-                }
-                else
-                {
-                    cutCentreEdges.append(centreEdge(faceI, false));
-                    centreEdgeWeights.append(weight);
-                }
-            }
-        }
-
-
-        forAll(f, fp)
-        {
-            //
-            // Check face internal (faceEdge) edge crossing
-            //
-
-            if (crosses(isoVal, ptField[f[fp]], fcVal, weight))
-            {
-                mark(ownerI, cutCells);
-                if (neighbourI != -1)
-                {
-                    mark(neighbourI, cutCells);
-                }
-
-                if (weight < tol)
-                {
-                    mark(f[fp], cutVerts);
-                }
-                else if (weight > 1-tol)
-                {
-                    mark(faceI, cutFaceCentres);
-                }
-                else
-                {
-                    cutFaceEdges.append(faceEdge(faceI, fp));
-                    faceEdgeWeights.append(weight);
-                }
-            }
-        }
-    }
-    
-
-    // Walk over all mesh points
-    forAll(mesh.points(), pointI)
-    {
-        const labelList& myCells = mesh.pointCells()[pointI];
-
-        forAll(myCells, myCellI)
-        {
-            label cellI = myCells[myCellI];
-
-            //
-            // Check pyramid edge crossing
-            //
-
-            scalar weight;
-
-            if (crosses(isoVal, ptField[pointI], volField[cellI], weight))
-            {
-                mark(cellI, cutCells);
-
-                if (weight < tol)
-                {
-                    mark(pointI, cutVerts);
-                }
-                else if (weight > 1-tol)
-                {
-                    mark(cellI, cutCellCentres);
-                }
-                else
-                {
-                    cutPyrEdges.append(pyramidEdge(pointI, cellI));
-                    pyrEdgeWeights.append(weight);
-                }
-            }
-        }
-    }            
-
-
-    // Walk over all mesh edges
-    forAll(mesh.edges(), edgeI)
-    {
-        const edge& e = mesh.edges()[edgeI];
-
-        scalar weight;
-
-        if (crosses(isoVal, ptField[e.start()], ptField[e.end()], weight))
-        {
-            mark(mesh.edgeCells()[edgeI], cutCells);
-
-            if (weight < tol)
-            {
-                mark(e.start(), cutVerts);
-            }
-            else if (weight > 1-tol)
-            {
-                mark(e.end(), cutVerts);
-            }
-            else
-            {
-                meshCutEdges.append(edgeI);
-                meshEdgeWeights.append(weight);
-            }
-        }
-    }
-
-    meshCutEdges.shrink();
-    meshEdgeWeights.shrink();
-
-    cutPyrEdges.shrink();
-    pyrEdgeWeights.shrink();
-
-    cutFaceEdges.shrink();
-    faceEdgeWeights.shrink();
-
-    cutCentreEdges.shrink();
-    centreEdgeWeights.shrink();
-
-
-
-    // Tranfer lists to faceDecompCuts
-
-    cells_ = cutCells.toc();
-    meshVerts_ = cutVerts.toc();
-    meshFaceCentres_ = cutFaceCentres.toc();
-    meshCellCentres_ = cutCellCentres.toc();
-
-    meshEdges_.transfer(meshCutEdges);
-    meshEdgeWeights_.transfer(meshEdgeWeights);
-
-    pyrEdges_.transfer(cutPyrEdges);
-    pyrEdgeWeights_.transfer(pyrEdgeWeights);
-
-    faceEdges_.transfer(cutFaceEdges);
-    faceEdgeWeights_.transfer(faceEdgeWeights);
-
-    centreEdges_.transfer(cutCentreEdges);
-    centreEdgeWeights_.transfer(centreEdgeWeights);
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::faceDecompIsoSurfaceCuts::faceDecompIsoSurfaceCuts
-(
-    const volScalarField& volField,
-    const pointScalarField& ptField,
-    const scalar isoVal,
-    const scalar tol
-)
-:
-    faceDecompCuts
-    (
-        volField.mesh(),
-        labelList(),
-
-        labelList(),
-        labelList(),
-        labelList(),
-
-        labelList(),            // mesh edges
-        scalarField(),
-
-        List<pyramidEdge>(),    // pyramid edges
-        scalarField(),
-
-        List<faceEdge>(),       // face edges
-        scalarField(),
-
-        List<centreEdge>(),     // cell centre edges
-        scalarField()
-    )
-{
-    constructEdgeCuts(volField, ptField, isoVal, tol);
-}
-
-
-// Construct from interpolator and field (does interpolation)
-Foam::faceDecompIsoSurfaceCuts::faceDecompIsoSurfaceCuts
-(
-    const volScalarField& volField,
-    const volPointInterpolation& pInterp,
-    const scalar isoVal,
-    const scalar tol
-)
-:
-    faceDecompCuts
-    (
-        volField.mesh(),
-        labelList(),
-
-        labelList(),
-        labelList(),
-        labelList(),
-
-        labelList(),            // mesh edges
-        scalarField(),
-
-        List<pyramidEdge>(),    // pyramid edges
-        scalarField(),
-
-        List<faceEdge>(),       // face edges
-        scalarField(),
-
-        List<centreEdge>(),     // cell centre edges
-        scalarField()
-    )
-{
-    // Get field on vertices
-    tmp<pointScalarField> tptField = pInterp.interpolate(volField);
-    const pointScalarField& ptField = tptField();
-
-
-    constructEdgeCuts(volField, ptField, isoVal, tol);
-}
-
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompIsoSurfaceCuts.H

@@ -1,102 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::faceDecompIsoSurfaceCuts
-
-Description
-    List of triangles based on isoSurface reconstruction on faceDecomp mesh.
-
-SourceFiles
-    faceDecompIsoSurfaceCuts.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef faceDecompIsoSurfaceCuts_H
-#define faceDecompIsoSurfaceCuts_H
-
-#include "faceDecompCuts.H"
-#include "labelHashSet.H"
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class volPointInterpolation;
-
-/*---------------------------------------------------------------------------*\
-                           Class faceDecompIsoSurfaceCuts Declaration
-\*---------------------------------------------------------------------------*/
-
-class faceDecompIsoSurfaceCuts
-:
-    public faceDecompCuts
-{
-    // Private Member Functions
-
-        //- Do all the hard work: determine edge crossings
-        void constructEdgeCuts
-        (
-            const volScalarField&,
-            const pointScalarField&,
-            const scalar isoVal,
-            const scalar tol
-        );
-
-public:
-
-    // Constructors
-
-        //- Construct from field on cells and field on vertices
-        faceDecompIsoSurfaceCuts
-        (
-            const volScalarField&,
-            const pointScalarField&,
-            const scalar isoVal,
-            const scalar tol
-        );
-
-        //- Construct from interpolator and field (does interpolation)
-        faceDecompIsoSurfaceCuts
-        (
-            const volScalarField&,
-            const volPointInterpolation&,
-            const scalar isoVal,
-            const scalar tol
-        );
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/meshEdgeCuts.C

@@ -1,109 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
-
-Description
-
-\*---------------------------------------------------------------------------*/
-
-#include "meshEdgeCuts.H"
-
-// * * * * * * * * * * * * * Private Static Functions  * * * * * * * * * * * //
-
-void Foam::meshEdgeCuts::mark(const label elem, labelHashSet& markedElems)
-{
-    labelHashSet::const_iterator iter = markedElems.find(elem);
-
-    if (iter == markedElems.end())
-    {
-        markedElems.insert(elem);
-    }
-}
-
-
-void Foam::meshEdgeCuts::mark
-(
-    const labelList& elems,
-    labelHashSet& markedElems
-)
-{
-    forAll(elems, elemI)
-    {
-        mark(elems[elemI], markedElems);
-    }
-}
-
-
-bool Foam::meshEdgeCuts::crosses
-(
-    const scalar isoVal,
-    const scalar val0,
-    const scalar val1,
-    scalar& weight
-)
-{
-    if
-    ( 
-        ((val0 <= isoVal) && (val1 >= isoVal))
-     || ((val1 <= isoVal) && (val0 >= isoVal))
-    )
-    {
-        if (val0 == val1)
-        {
-            return false;
-        }
-        else
-        {
-            weight = (isoVal - val0)/((val1 - val0) + VSMALL);
-
-            return true;
-        }
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-Foam::meshEdgeCuts::meshEdgeCuts
-(
-    const primitiveMesh& mesh,
-    const labelList& cells,
-    const labelList& meshVerts,
-    const labelList& meshEdges,
-    const scalarField& meshEdgeWeights
-)
-:
-    mesh_(mesh),
-    cells_(cells),
-    meshVerts_(meshVerts),
-    meshEdges_(meshEdges),
-    meshEdgeWeights_(meshEdgeWeights)
-{}
-
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/meshEdgeCuts.H

@@ -1,161 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::meshEdgeCuts
-
-Description
-    Container for cuts of edges of mesh. Cuts either exactly through existing 
-    vertices or through edges.
-
-SourceFiles
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef meshEdgeCuts_H
-#define meshEdgeCuts_H
-
-#include "labelList.H"
-#include "labelHashSet.H"
-#include "scalarField.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class primitiveMesh;
-
-/*---------------------------------------------------------------------------*\
-                           Class meshEdgeCuts Declaration
-\*---------------------------------------------------------------------------*/
-
-class meshEdgeCuts
-{
-    // Private data
-
-        const primitiveMesh& mesh_;
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        meshEdgeCuts(const meshEdgeCuts&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const meshEdgeCuts&);
-
-
-protected:
-
-        //- List of cells containing the cuts
-        labelList cells_;
-
-        //- Points exactly cut by cuts
-        labelList meshVerts_;
-
-        //- List of edge labels cut
-        labelList meshEdges_;
-
-        //- Positions on edges
-        scalarField meshEdgeWeights_;
-
-
-    // Protected Static Functions
-
-        //- Mark element in a hashSet
-        static void mark(const label elem, labelHashSet& markedElems);
-
-        //- Mark list of elements in a hashSet
-        static void mark(const labelList& elems, labelHashSet& markedElems);
-
-        //- Return true and set weight if linear interpolation between
-        //  val0 and val1 crosses isoVal. weight=1 if isoVal==val1
-        static bool crosses
-        (
-            const scalar isoVal,
-            const scalar val0,
-            const scalar val1,
-            scalar& weight
-        );
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        meshEdgeCuts
-        (
-            const primitiveMesh& mesh,
-            const labelList& cells,
-            const labelList& meshVerts,
-            const labelList& meshEdges,
-            const scalarField& meshEdgeWeights
-        );
-
-
-    // Member Functions
-
-        const primitiveMesh& mesh() const
-        {
-            return mesh_;
-        }
-
-        const labelList& cells() const
-        {
-            return cells_;
-        }
-
-        const labelList& meshVerts() const
-        {
-            return meshVerts_;
-        }
-
-        const labelList& meshEdges() const
-        {
-            return meshEdges_;
-        }
-
-        const scalarField& meshEdgeWeights() const
-        {
-            return meshEdgeWeights_;
-        }
-
-        label size() const
-        {
-            return meshVerts_.size() + meshEdges_.size();
-        }
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/cellAddressing.C

@@ -1,132 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
-
-Description
-
-\*---------------------------------------------------------------------------*/
-
-#include "cellAddressing.H"
-#include "labelHashSet.H"
-#include "cellModel.H"
-
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-Foam::labelList Foam::cellAddressing::makeIdentity(const label size)
-{
-    labelList result(size);
-    forAll(result, resultI)
-    {
-        result[resultI] = resultI;
-    }
-    return result;
-}
-
-
-Foam::label Foam::cellAddressing::findEdge(const edge& e)
-{
-    forAll(edges_, edgeI)
-    {
-        if (edges_[edgeI] == e)
-        {
-            return edgeI;
-        }
-    }
-    FatalErrorIn
-    (
-        "cellAddressing::findEdge(const edge&)"
-    )   << "Problem: cannot find edge " << e << " in edges " << edges_
-        << exit(FatalError);
-
-    return -1;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-Foam::cellAddressing::cellAddressing(const Foam::cellModel& model)
-:
-    modelFaces_(model.modelFaces()),
-    edges_(model.edges(makeIdentity(model.nPoints()))),
-    edgeFaces_(model.nEdges()),
-    faceEdges_(model.nFaces()),
-    pointEdges_(model.nPoints())
-{
-    List<labelHashSet > dynEdgeFaces(edgeFaces_.size());
-
-    forAll(faceEdges_, faceI)
-    {
-        DynamicList<label> fEdges;
-
-        const face& f = modelFaces_[faceI];
-
-        forAll(f, fp)
-        {
-            label edgeI = findEdge(edge(f[fp], f[(fp+1)%f.size()]));
-
-            // Append to faceEdges
-            fEdges.append(edgeI);
-
-            // Append to edgefaces
-            labelHashSet& eFaces = dynEdgeFaces[edgeI];
-
-            if (!eFaces.found(faceI))
-            {
-                eFaces.insert(faceI);
-            }
-        }
-        fEdges.shrink();
-
-        // Copy into faceEdges_
-        faceEdges_[faceI].transfer(fEdges.shrink());
-    }
-
-    // Convert dynEdgeFaces into edgeFaces_
-
-    forAll(dynEdgeFaces, edgeI)
-    {
-        edgeFaces_[edgeI] = dynEdgeFaces[edgeI].toc();
-    }
-
-
-    // Convert edges_ into pointEdges
-    List<DynamicList<label> > dynPointEdges(model.nPoints());
-
-    forAll(edges_, edgeI)
-    {
-        const edge& e = edges_[edgeI];
-
-        dynPointEdges[e.start()].append(edgeI);
-        dynPointEdges[e.end()].append(edgeI);
-    }
-
-    forAll(dynPointEdges, pointI)
-    {
-        pointEdges_[pointI].transfer(dynPointEdges[pointI].shrink());
-    }
-}
-
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/cellAddressing.H

@@ -1,145 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::cellAddressing
-
-Description
-    Additional addressing for a cellModel. (cellModel themselves only
-    contain cell-points). Used in cellFeatures class.
-
-SourceFiles
-    cellAddressing.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef cellAddressing_H
-#define cellAddressing_H
-
-#include "cellModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class cellModel;
-
-/*---------------------------------------------------------------------------*\
-                           Class cellAddressing Declaration
-\*---------------------------------------------------------------------------*/
-
-class cellAddressing
-{
-    // Private data
-
-        faceList modelFaces_;
-
-        edgeList edges_;
-
-        labelListList edgeFaces_;
-
-        labelListList faceEdges_;
-
-        labelListList pointEdges_;
-
-    // Private static functions        
-
-        static labelList makeIdentity(const label size);
-
-    // Private static functions        
-
-        label findEdge(const edge& e); 
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        cellAddressing(const cellModel& model);
-
-    // Destructor
-
-    // Member Functions
-
-        // Access
-
-            inline const faceList& modelFaces() const
-            {
-                return modelFaces_;
-            }
-
-            inline const edgeList& edges() const
-            {
-                return edges_;
-            }
-
-            inline const labelListList& edgeFaces() const
-            {
-                return edgeFaces_;
-            }
-
-            inline const labelListList& faceEdges() const
-            {
-                return faceEdges_;
-            }
-
-            inline const labelListList& pointEdges() const
-            {
-                return pointEdges_;
-            }
-
-
-        // Check
-
-        // Edit
-
-        // Write
-
-
-    // Member Operators
-
-    // Friend Functions
-
-    // Friend Operators
-
-    // IOstream Operators
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-//#include "cellAddressingI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurface.H

@@ -1,209 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::meshCutSurface
-
-Description
-    Description of surface resulting from cuts through mesh. Region on
-    labelledTris is the cell the triangle is resulting from.
-
-SourceFiles
-    meshCutSurface.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef meshCutSurface_H
-#define meshCutSurface_H
-
-#include "labelList.H"
-#include "HashTable.H"
-#include "pointField.H"
-#include "triSurface.H"
-#include "boolList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-class primitiveMesh;
-class cellAddressing;
-class faceDecompCuts;
-class cellDecompCuts;
-
-/*---------------------------------------------------------------------------*\
-                           Class meshCutSurface Declaration
-\*---------------------------------------------------------------------------*/
-
-class meshCutSurface
-:
-    public triSurface
-{
-
-    // Static data members
-
-        //- Symbolic names of tet vertices
-        static label FP0;
-        static label FC;
-        static label FP1;
-        static label CC;
-
-        //- Symbolic names of tet edges
-        static label FP0_FC;
-        static label FP0_FP1;
-        static label FP0_CC;
-        static label FC_CC;
-        static label FC_FP1;
-        static label FP1_CC;
-
-
-    // Private Static Functions
-
-        //- Find key in hashtable. Returns value or -1 if not found
-        template<class T, class Hash>
-        static label find
-        (
-            const HashTable<label, T, Hash>& table,
-            const T& key
-        );
-
-        //- Find edge in subset edgeLabels of edges.
-        static label findEdge
-        (
-            const edgeList& edges,
-            const labelList& edgeLabels,
-            const edge& e
-        );
-
-        //- Find the other (i.e. not cutEdgeI) cut-edge on faceI of cell model
-        static label findCutEdge
-        (
-            const cellAddressing& model,
-            const labelList& tetEdgeCuts,
-            const label faceI,
-            const label cutEdgeI
-        );
-
-
-    // Private Member Functions
-
-        //- Remove duplicate faces (resulting from 3 meshVerts cut - is seen
-        //  from both sides)
-        static triSurface removeDuplicates(const triSurface& surf);
-
-
-public:
-
-        //- Generic tet cutting routine.
-        //  Generate triangles given the cut vertices and edges of the tet.
-        //  (Triangles can be incorrectly oriented)
-        void cutTet
-        (
-            const cellAddressing& model,
-            const labelList& tetVertCuts,
-            const labelList& tetEdgeCuts,
-            const label cellI,
-            List<labelledTri>& tris,
-            label& triI
-        );
-
-        //- Tet cutting when cuts are only through edges.
-        void cutTetThroughEdges
-        (
-            const cellAddressing& model,
-            const labelList& tetEdgeCuts,
-            const label nEdgeCuts,
-            const label cellI,
-            List<labelledTri>& tris,
-            label& triI
-        );
-
-        //- Tet cutting when cuts are through corner vertices (and possibly
-        //  edges)
-        void cutTetThroughVerts
-        (
-            const cellAddressing& model,
-            const labelList& tetVertCuts,
-            const label nVertCuts,
-            const labelList& tetEdgeCuts,
-            const label nEdgeCuts,
-            const label cellI,
-            List<labelledTri>& tris,
-            label& triI
-        );
-
-
-    // Static functions
-
-        //- Interpolate onto points given edge weights, cell centre values,
-        //  face centre values and vertex values
-        template <class T>
-        static tmp<Field<T> > interpolate
-        (
-            const faceDecompCuts& cuts,
-            const Field<T>& vField,
-            const Field<T>& fField,
-            const Field<T>& pField
-        );
-
-        //- Interpolate onto points given edge weights, cell centre values,
-        //  and vertex values
-        template <class T>
-        static tmp<Field<T> > interpolate
-        (
-            const cellDecompCuts& cuts,
-            const Field<T>& vField,
-            const Field<T>& pField
-        );
-
-
-    // Constructors
-
-        //- Construct from cuts on faceDecomp tet decomposition
-        //  (introduces cellcentres and facecentres)
-        meshCutSurface(const faceDecompCuts& cuts);
-
-        //- Construct from cuts on cellDecomp tet decomposition
-        //  (introduces cellcentres but not facecentres)
-        meshCutSurface(const cellDecompCuts& cuts);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#ifdef NoRepository
-#   include "meshCutSurfaceInterpolate.C"
-#endif
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurfaceInterpolate.C

@@ -1,174 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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 "meshCutSurface.H"
-#include "pyramidEdge.H"
-#include "faceEdge.H"
-#include "centreEdge.H"
-#include "diagonalEdge.H"
-#include "faceDecompCuts.H"
-#include "cellDecompCuts.H"
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template <class T>
-Foam::tmp<Foam::Field<T> > Foam::meshCutSurface::interpolate
-(
-    const faceDecompCuts& cuts,
-    const Field<T>& vField,
-    const Field<T>& fField,
-    const Field<T>& pField
-)
-{
-    const primitiveMesh& mesh = cuts.mesh();
-
-    tmp<Field<T> > tvals(new Field<T>(cuts.size()));
-    Field<T>& vals = tvals();
-
-    label pointI = 0;
-
-    forAll(cuts.meshVerts(), cutVertI)
-    {
-        vals[pointI++] = pField[cuts.meshVerts()[cutVertI]];
-    }
-
-    forAll(cuts.meshFaceCentres(), cutFaceCentreI)
-    {
-        label faceI = cuts.meshFaceCentres()[cutFaceCentreI];
-
-        vals[pointI++] = fField[faceI];
-    }
-
-    forAll(cuts.meshCellCentres(), cutCellCentreI)
-    {
-        label cellI = cuts.meshCellCentres()[cutCellCentreI];
-
-        vals[pointI++] = vField[cellI];
-    }
-
-    forAll(cuts.meshEdges(), meshCutSurfaceEdgeI)
-    {
-        label edgeI = cuts.meshEdges()[meshCutSurfaceEdgeI];
-
-        const edge& e = mesh.edges()[edgeI];
-
-        scalar weight = cuts.meshEdgeWeights()[meshCutSurfaceEdgeI];
-
-        vals[pointI++] = weight*pField[e.start()] + (1-weight)*pField[e.end()];
-    }
-
-    forAll(cuts.pyrEdges(), edgeI)
-    {
-        const pyramidEdge& e = cuts.pyrEdges()[edgeI];
-
-        scalar weight = cuts.pyrEdgeWeights()[edgeI];
-
-        vals[pointI++] = e.interpolate(vField, pField, weight);
-    }
-
-    forAll(cuts.centreEdges(), edgeI)
-    {
-        const centreEdge& e = cuts.centreEdges()[edgeI];
-
-        scalar weight = cuts.centreEdgeWeights()[edgeI];
-
-        vals[pointI++] = e.interpolate(mesh, vField, fField, weight);
-    }
-
-    forAll(cuts.faceEdges(), edgeI)
-    {
-        const faceEdge& e = cuts.faceEdges()[edgeI];
-
-        scalar weight = cuts.faceEdgeWeights()[edgeI];
-
-        vals[pointI++] = e.interpolate(mesh, fField, pField, weight);
-    }
-
-    return tvals;
-}
-
-
-template <class T>
-Foam::tmp<Foam::Field<T> > Foam::meshCutSurface::interpolate
-(
-    const cellDecompCuts& cuts,
-    const Field<T>& vField,
-    const Field<T>& pField
-)
-{
-    const primitiveMesh& mesh = cuts.mesh();
-
-    tmp<Field<T> > tvals(new Field<T>(cuts.size()));
-    Field<T>& vals = tvals();
-
-    label pointI = 0;
-
-    forAll(cuts.meshVerts(), cutVertI)
-    {
-        vals[pointI++] = pField[cuts.meshVerts()[cutVertI]];
-    }
-
-    forAll(cuts.meshCellCentres(), cutCellCentreI)
-    {
-        label cellI = cuts.meshCellCentres()[cutCellCentreI];
-
-        vals[pointI++] = vField[cellI];
-    }
-
-    forAll(cuts.meshEdges(), meshCutSurfaceEdgeI)
-    {
-        label edgeI = cuts.meshEdges()[meshCutSurfaceEdgeI];
-
-        const edge& e = mesh.edges()[edgeI];
-
-        scalar weight = cuts.meshEdgeWeights()[meshCutSurfaceEdgeI];
-
-        vals[pointI++] = weight*pField[e.start()] + (1-weight)*pField[e.end()];
-    }
-
-    forAll(cuts.pyrEdges(), edgeI)
-    {
-        const pyramidEdge& e = cuts.pyrEdges()[edgeI];
-
-        scalar weight = cuts.pyrEdgeWeights()[edgeI];
-
-        vals[pointI++] = e.interpolate(vField, pField, weight);
-    }
-
-    forAll(cuts.diagEdges(), edgeI)
-    {
-        const diagonalEdge& e = cuts.diagEdges()[edgeI];
-
-        scalar weight = cuts.diagEdgeWeights()[edgeI];
-
-        vals[pointI++] = e.interpolate(mesh, pField, weight);
-    }
-
-    return tvals;
-}
-
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/centreEdge.H

@@ -1,187 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::centreEdge
-
-Description
-    Implicit description of faceCentre to cellCentre edge
-    (from tet decomposition).
-
-    Stores
-      - face label
-      - whether in owner or neighbour of face
-
-SourceFiles
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef centreEdge_H
-#define centreEdge_H
-
-#include "label.H"
-#include "primitiveMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-
-/*---------------------------------------------------------------------------*\
-                           Class centreEdge Declaration
-\*---------------------------------------------------------------------------*/
-
-class centreEdge
-{
-    // Private data
-
-        //- face label
-        label faceLabel_;
-
-        //- is in owner or neighbour of face
-        bool inOwner_;
-
-public:
-
-    // Public classes
-
-        //- Hash function
-        class centreEdgeHash
-        {
-
-        public:
-
-            centreEdgeHash()
-            {}
-
-            //- simple hashing function of labels
-            label operator()(const centreEdge& ce, const label tableSize) const
-            {
-                return (ce.faceLabel() + ce.inOwner()) % tableSize;
-            }
-        };
-
-
-    // Constructors
-
-        //- Construct null
-        inline centreEdge()
-        :
-            faceLabel_(-1),
-            inOwner_(-1)
-        {}
-
-        //- Construct from components
-        inline centreEdge(const label faceLabel, const bool inOwner)
-        :
-            faceLabel_(faceLabel),
-            inOwner_(inOwner)
-        {}
-
-    // Member Functions
-
-        label faceLabel() const
-        {   
-            return faceLabel_;
-        }
-
-        bool inOwner() const
-        {   
-            return inOwner_;
-        }
-
-        template <class T>
-        T interpolate
-        (
-            const primitiveMesh& mesh,
-            const Field<T>& cellField,
-            const Field<T>& faceField,
-            const scalar weight
-        ) const
-        {
-            label cellI;
-            if (inOwner_)
-            {
-                cellI = mesh.faceOwner()[faceLabel_];
-            }
-            else
-            {
-                cellI = mesh.faceNeighbour()[faceLabel_];
-            }
-
-            return
-                (1-weight)*faceField[faceLabel_]
-              + weight*cellField[cellI];
-        }
-
-        point coord(const primitiveMesh& mesh, const scalar weight) const
-        {
-            return interpolate
-            (
-                mesh,
-                mesh.cellCentres(),
-                mesh.faceCentres(),
-                weight
-            );
-        }
-
-
-    // Member Operators
-
-        bool operator==(const centreEdge& ce) const
-        {
-            return
-                (faceLabel() == ce.faceLabel())
-             && (inOwner() == ce.inOwner());
-        }
-
-
-    // IOstream Operators
-
-        inline friend Ostream& operator<<(Ostream& os, const centreEdge& ce)
-        {
-            os  << token::BEGIN_LIST
-                << ce.faceLabel_ << token::SPACE
-                << ce.inOwner_
-                << token::END_LIST;
-
-            // Check state of Ostream
-            os.check("Ostream& operator<<(Ostream&, const centreEdge&)");
-
-            return os;
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/diagonalEdge.H

@@ -1,199 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::diagonalEdge
-
-Description
-    Implicit description of diagonal edge (from tet decomposition).
-    Diangonal edge is edge between two vertices of same face.
-
-    Stores
-        - face label
-        - indices in face
-
-SourceFiles
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef diagonalEdge_H
-#define diagonalEdge_H
-
-#include "label.H"
-#include "primitiveMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class diagonalEdge Declaration
-\*---------------------------------------------------------------------------*/
-
-class diagonalEdge
-{
-    // Private data
-
-        //- face label
-        label faceLabel_;
-
-        //- index in face
-        label faceIndex0_;
-
-        //- index in face
-        label faceIndex1_;
-
-public:
-
-    // Public classes
-
-        //- Hash function
-        class diagonalEdgeHash
-        {
-
-        public:
-
-            diagonalEdgeHash()
-            {}
-
-            //- simple hashing function of labels. Commutative on face indices
-            //  since edge is same.
-            label operator()(const diagonalEdge& de, const label tableSize)
-                const
-            {
-                // Note: mag since multiply might overflow and produce
-                // negative numbers
-                return
-                    mag
-                    (
-                        de.faceLabel()
-                      + (de.faceLabel()*de.faceLabel())
-                      + de.faceIndex0() * de.faceIndex1()
-                      + (de.faceIndex0() + de.faceIndex1())
-                    )
-                  % tableSize;
-            }
-        };
-
-
-    // Constructors
-
-        //- Construct null
-        inline diagonalEdge()
-        :
-            faceLabel_(-1),
-            faceIndex0_(-1),
-            faceIndex1_(-1)
-        {}
-
-        //- Construct from components
-        inline diagonalEdge
-        (
-            const label faceLabel,
-            const label faceIndex0,
-            const label faceIndex1
-        )
-        :
-            faceLabel_(faceLabel),
-            faceIndex0_(faceIndex0),
-            faceIndex1_(faceIndex1)
-        {}
-
-    // Member Functions
-
-        label faceLabel() const
-        {   
-            return faceLabel_;
-        }
-
-        label faceIndex0() const
-        {   
-            return faceIndex0_;
-        }
-
-        label faceIndex1() const
-        {   
-            return faceIndex1_;
-        }
-
-        template <class T>
-        T interpolate
-        (
-            const primitiveMesh& mesh,
-            const Field<T>& vertField,
-            const scalar weight
-        ) const
-        {
-            const face& f = mesh.faces()[faceLabel_];
-
-            return
-                (1-weight)*vertField[f[faceIndex0_]]
-              + weight*vertField[f[faceIndex1_]];
-        }
-
-        point coord(const primitiveMesh& mesh, const scalar weight) const
-        {
-            return interpolate(mesh, mesh.points(), weight);
-        }
-
-
-    // Member Operators
-
-        bool operator==(const diagonalEdge& de) const
-        {
-            return
-                (faceLabel() == de.faceLabel())
-             && (faceIndex0() == de.faceIndex0())
-             && (faceIndex1() == de.faceIndex1());
-        }
-
-
-    // IOstream Operators
-
-        inline friend Ostream& operator<<(Ostream& os, const diagonalEdge& de)
-        {
-            os  << token::BEGIN_LIST
-                << de.faceLabel_ << token::SPACE
-                << de.faceIndex0_ << token::SPACE
-                << de.faceIndex1_
-                << token::END_LIST;
-
-            // Check state of Ostream
-            os.check("Ostream& operator<<(Ostream&, const diagonalEdge&)");
-
-            return os;
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/faceEdge.H

@@ -1,186 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::faceEdge
-
-Description
-    Implicit description of face edge (from tet decomposition).
-    Face edge is edge between vertex of face and face centre.
-
-    Stores
-        - face label
-        - index in face (gives the mesh vertex)
-
-
-SourceFiles
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef faceEdge_H
-#define faceEdge_H
-
-#include "label.H"
-#include "primitiveMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-
-/*---------------------------------------------------------------------------*\
-                           Class faceEdge Declaration
-\*---------------------------------------------------------------------------*/
-
-class faceEdge
-{
-    // Private data
-
-        //- face label
-        label faceLabel_;
-
-        //- index in face
-        label faceIndex_;
-
-public:
-
-    // Public classes
-
-        //- Hash function
-        class faceEdgeHash
-        {
-
-        public:
-
-            faceEdgeHash()
-            {}
-
-            //- simple hashing function of labels
-            label operator()(const faceEdge& fe, const label tableSize)
-                const
-            {
-                // Note: mag since multiply might overflow and produce
-                // negative numbers
-                return
-                    mag
-                    (
-                        fe.faceLabel()
-                      + (fe.faceLabel() * fe.faceLabel())
-                      + fe.faceIndex()
-                      + (fe.faceIndex() * fe.faceIndex())
-                    )
-                  % tableSize;
-            }
-        };
-
-
-    // Constructors
-
-        //- Construct null
-        inline faceEdge()
-        :
-            faceLabel_(-1),
-            faceIndex_(-1)
-        {}
-
-        //- Construct from components
-        inline faceEdge(const label faceLabel, const label faceIndex)
-        :
-            faceLabel_(faceLabel),
-            faceIndex_(faceIndex)
-        {}
-
-    // Member Functions
-
-        label faceLabel() const
-        {
-            return faceLabel_;
-        }
-
-        label faceIndex() const
-        {
-            return faceIndex_;
-        }
-
-        template <class T>
-        T interpolate
-        (
-            const primitiveMesh& mesh,
-            const Field<T>& faceField,
-            const Field<T>& vertField,
-            const scalar weight
-        ) const
-        {
-            const face& f = mesh.faces()[faceLabel_];
-
-            return
-                (1-weight)*vertField[f[faceIndex_]]
-              + weight*faceField[faceLabel_];
-        }
-
-
-        point coord(const primitiveMesh& mesh, const scalar weight) const
-        {
-            return interpolate(mesh, mesh.faceCentres(), mesh.points(), weight);
-        }
-
-
-    // Member Operators
-
-        bool operator==(const faceEdge& fe) const
-        {
-            return
-                (faceLabel() == fe.faceLabel())
-             && (faceIndex() == fe.faceIndex());
-        }
-
-
-    // IOstream Operators
-
-        inline friend Ostream& operator<<(Ostream& os, const faceEdge& fe)
-        {
-            os  << token::BEGIN_LIST
-                << fe.faceLabel_ << token::SPACE
-                << fe.faceIndex_
-                << token::END_LIST;
-
-            // Check state of Ostream
-            os.check("Ostream& operator<<(Ostream&, const faceEdge&)");
-
-            return os;
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/pyramidEdge.H

@@ -1,185 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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::pyramidEdge
-
-Description
-    Implicit description of pyramid edge (from tet decomposition).
-    Pyramid edge is edge between mesh vertex and cell centre.
-
-    Stores
-        - vertex label
-        - cell label
-
-SourceFiles
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef pyramidEdge_H
-#define pyramidEdge_H
-
-#include "label.H"
-#include "primitiveMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-
-/*---------------------------------------------------------------------------*\
-                           Class pyramidEdge Declaration
-\*---------------------------------------------------------------------------*/
-
-class pyramidEdge
-{
-    // Private data
-
-        //- vertex label
-        label vertexLabel_;
-
-        //- cell label
-        label cellLabel_;
-
-public:
-
-    // Public classes
-
-        //- Hash function
-        class pyramidEdgeHash
-        {
-
-        public:
-
-            pyramidEdgeHash()
-            {}
-
-            //- simple hashing function of labels
-            label operator()(const pyramidEdge& pe, const label tableSize) const
-            {
-                // Note: mag since multiply might overflow and produce
-                // negative numbers
-                return
-                    mag
-                    (
-                        pe.vertexLabel()
-                      + (pe.vertexLabel() * pe.vertexLabel())
-                      + pe.cellLabel()
-                      + (pe.cellLabel() * pe.cellLabel())
-                    )
-                  % tableSize;
-            }
-        };
-
-
-    // Constructors
-
-        //- Construct null
-        inline pyramidEdge()
-        :
-            vertexLabel_(-1),
-            cellLabel_(-1)
-        {}
-
-        //- Construct from components
-        inline pyramidEdge
-        (
-            const label vertexLabel,
-            const label cellLabel
-        )
-        :
-            vertexLabel_(vertexLabel),
-            cellLabel_(cellLabel)
-        {}
-
-    // Member Functions
-
-        label vertexLabel() const
-        {   
-            return vertexLabel_;
-        }
-
-        label cellLabel() const
-        {   
-            return cellLabel_;
-        }
-
-        template <class T>
-        T interpolate
-        (
-            const Field<T>& cellField,
-            const Field<T>& vertField,
-            const scalar weight
-        ) const
-        {
-            return
-                (1-weight)*vertField[vertexLabel_]
-              + weight*cellField[cellLabel_];
-        }
-
-
-        point coord(const primitiveMesh& mesh, const scalar weight) const
-        {
-            return interpolate(mesh.cellCentres(), mesh.points(), weight);
-        }
-
-
-    // Member Operators
-
-        bool operator==(const pyramidEdge& pe) const
-        {
-            return
-                (vertexLabel() == pe.vertexLabel())
-             && (cellLabel() == pe.cellLabel());
-        }
-
-
-    // IOstream Operators
-
-        inline friend Ostream& operator<<(Ostream& os, const pyramidEdge& pe)
-        {
-            os  << token::BEGIN_LIST
-                << pe.vertexLabel_ << token::SPACE
-                << pe.cellLabel_
-                << token::END_LIST;
-
-            // Check state of Ostream
-            os.check("Ostream& operator<<(Ostream&, const pyramidEdge&)");
-
-            return os;
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/sampling/Make/files

@@ -24,6 +24,9 @@ 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
 sampledSurface/sampledSurfacesFunctionObject/sampledSurfacesFunctionObject.C

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -0,0 +1,342 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "distanceSurface.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(distanceSurface, 0);
+    addToRunTimeSelectionTable(sampledSurface, distanceSurface, word);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::distanceSurface::createGeometry() const
+{
+    // Clear any stored topo
+    facesPtr_.clear();
+
+    // Distance to cell centres
+    scalarField cellDistance(mesh().nCells());
+    {
+        List<pointIndexHit> nearest;
+        surfPtr_().findNearest
+        (
+            mesh().cellCentres(),
+            scalarField(mesh().nCells(), GREAT),
+            nearest
+        );
+
+        if (signed_)
+        {
+            vectorField normal;
+            surfPtr_().getNormal(nearest, normal);
+
+            forAll(cellDistance, cellI)
+            {
+                vector d(mesh().cellCentres()[cellI]-nearest[cellI].hitPoint());
+
+                if ((d&normal[cellI]) > 0)
+                {
+                    cellDistance[cellI] = Foam::mag(d);
+                }
+                else
+                {
+                    cellDistance[cellI] = -Foam::mag(d);
+                }
+            }
+        }
+        else
+        {
+            forAll(cellDistance, cellI)
+            {
+                cellDistance[cellI] = Foam::mag
+                (
+                    nearest[cellI].hitPoint()
+                  - mesh().cellCentres()[cellI]
+                );
+            }
+        }
+    }
+
+    // Distance to points
+    scalarField pointDistance(mesh().nPoints());
+    {
+        List<pointIndexHit> nearest;
+        surfPtr_().findNearest
+        (
+            mesh().points(),
+            scalarField(mesh().nPoints(), GREAT),
+            nearest
+        );
+        forAll(pointDistance, pointI)
+        {
+            pointDistance[pointI] = Foam::mag
+            (
+                nearest[pointI].hitPoint()
+              - mesh().points()[pointI]
+            );
+        }
+    }
+
+    //- Direct from cell field and point field.
+    const isoSurfaceCell iso
+    (
+        mesh(),
+        cellDistance,
+        pointDistance,
+        distance_,
+        regularise_
+    );
+
+    ////- From point field and interpolated cell.
+    //scalarField cellAvg(mesh().nCells(), scalar(0.0));
+    //labelField nPointCells(mesh().nCells(), 0);
+    //{
+    //    for (label pointI = 0; pointI < mesh().nPoints(); pointI++)
+    //    {
+    //        const labelList& pCells = mesh().pointCells(pointI);
+    //
+    //        forAll(pCells, i)
+    //        {
+    //            label cellI = pCells[i];
+    //
+    //            cellAvg[cellI] += pointDistance[pointI];
+    //            nPointCells[cellI]++;
+    //        }
+    //    }
+    //}
+    //forAll(cellAvg, cellI)
+    //{
+    //    cellAvg[cellI] /= nPointCells[cellI];
+    //}
+    //
+    //const isoSurface iso
+    //(
+    //    mesh(),
+    //    cellAvg,
+    //    pointDistance,
+    //    distance_,
+    //    regularise_
+    //);
+
+
+    const_cast<distanceSurface&>(*this).triSurface::operator=(iso);
+    meshCells_ = iso.meshCells();
+
+    if (debug)
+    {
+        print(Pout);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::distanceSurface::distanceSurface
+(
+    const word& name,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    sampledSurface(name, mesh, dict),
+    surfPtr_
+    (
+        searchableSurface::New
+        (
+            dict.lookup("surfaceType"),
+            IOobject
+            (
+                dict.lookupOrDefault("surfaceName", name),  // name
+                mesh.time().constant(),                     // directory
+                "triSurface",                               // instance
+                mesh.time(),                                // registry
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            ),
+            dict
+        )
+    ),
+    distance_(readScalar(dict.lookup("distance"))),
+    signed_(readBool(dict.lookup("signed"))),
+    regularise_(dict.lookupOrDefault("regularise", 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;
+//        }
+//    }
+    correct(true);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::distanceSurface::~distanceSurface()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::distanceSurface::correct(const bool meshChanged)
+{
+    // Only change of mesh changes plane - zone restriction gets lost
+    if (meshChanged)
+    {
+        createGeometry();
+    }
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::distanceSurface::sample
+(
+    const volScalarField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::distanceSurface::sample
+(
+    const volVectorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::distanceSurface::sample
+(
+    const volSphericalTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::distanceSurface::sample
+(
+    const volSymmTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::distanceSurface::sample
+(
+    const volTensorField& vField
+) const
+{
+    return sampleField(vField);
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::distanceSurface::interpolate
+(
+    const interpolation<scalar>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::vectorField>
+Foam::distanceSurface::interpolate
+(
+    const interpolation<vector>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+Foam::tmp<Foam::sphericalTensorField>
+Foam::distanceSurface::interpolate
+(
+    const interpolation<sphericalTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::symmTensorField>
+Foam::distanceSurface::interpolate
+(
+    const interpolation<symmTensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+Foam::tmp<Foam::tensorField>
+Foam::distanceSurface::interpolate
+(
+    const interpolation<tensor>& interpolator
+) const
+{
+    return interpolateField(interpolator);
+}
+
+
+void Foam::distanceSurface::print(Ostream& os) const
+{
+    os  << "distanceSurface: " << name() << " :"
+        << "  surface:" << surfPtr_().name()
+        << "  distance:" << distance_
+        << "  faces:" << faces().size()
+        << "  points:" << points().size();
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/distanceSurface/distanceSurface.H

@@ -0,0 +1,237 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::distanceSurface
+
+Description
+    A sampledSurface defined by a distance to a surface.
+
+SourceFiles
+    distanceSurface.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef distanceSurface_H
+#define distanceSurface_H
+
+#include "sampledSurface.H"
+#include "triSurface.H"
+#include "searchableSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class distanceSurface Declaration
+\*---------------------------------------------------------------------------*/
+
+class distanceSurface
+:
+    public sampledSurface,
+    public triSurface
+{
+    // Private data
+
+        //- Surface
+        const autoPtr<searchableSurface> surfPtr_;
+
+        //- distance value
+        const scalar distance_;
+
+        //- signed distance
+        const bool signed_;
+
+        //- Whether to coarsen
+        const Switch regularise_;
+
+        //- 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("distanceSurface");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        distanceSurface
+        (
+            const word& name,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~distanceSurface();
+
+
+    // 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 "distanceSurfaceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C

@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "distanceSurface.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::distanceSurface::sampleField
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vField
+) const
+{
+    return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
+}
+
+
+template <class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::distanceSurface::interpolateField
+(
+    const interpolation<Type>& interpolator
+) const
+{
+    // 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/isoSurface.H

@@ -27,7 +27,20 @@ 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",
+    G.M. Treece, R.W. Prager and A.H. Gee.
+
+    Note:
+    - in parallel the regularisation (coarsening) always takes place
+    and slightly different surfaces will be created compared to non-parallel.
+    The surface will still be continuous though!
+    - 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 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
     isoSurface.C
@@ -38,13 +51,17 @@ SourceFiles
 #define isoSurface_H
 
 #include "triSurface.H"
+#include "labelPair.H"
+#include "pointIndexHit.H"
+#include "PackedList.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class polyMesh;
+class fvMesh;
 
 /*---------------------------------------------------------------------------*\
                        Class isoSurface Declaration
@@ -56,18 +73,39 @@ class isoSurface
 {
     // 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_;
-
-        //- Reference to cell values
-        const scalarField& cellValues_;
-
-        //- Reference to point values
-        const scalarField& pointValues_;
+        const fvMesh& mesh_;
 
         //- Isosurface value
         const scalar iso_;
 
+        //- When to merge points
+        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_;
@@ -78,31 +116,210 @@ class isoSurface
 
     // Private Member Functions
 
-        template<class T>
-        static T vertexInterp
+        //- Get location of iso value as fraction inbetween s0,s1
+        scalar isoFraction
         (
-            const scalar iso,
-            const T& p0,
-            const T& p1,
             const scalar s0,
             const scalar s1
+        ) const;
+
+        //- Set faceCutType,cellCutType.
+        void calcCutTypes
+        (
+            const volScalarField& cVals,
+            const scalarField& pVals
         );
 
-        template<class T>
-        static void vertexInterp
+        static labelPair findCommonPoints
         (
-            const scalar iso,
+            const labelledTri&,
+            const labelledTri&
+        );
+
+        static point calcCentre(const triSurface&);
+
+        static pointIndexHit collapseSurface
+        (
+            pointField& localPoints,
+            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 labelList& boundaryRegion,
+            const volScalarField& cVals,
+            const scalarField& pVals,
+            DynamicList<point>& snappedPoints,
+            labelList& snappedCc
+        ) const;
+
+        //- Determine per point whether all near cuts can be snapped to single
+        //  point.
+        void calcSnappedPoint
+        (
+            const PackedList<1>& isBoundaryPoint,
+            const labelList& boundaryRegion,
+            const volScalarField& cVals,
+            const scalarField& pVals,
+            DynamicList<point>& snappedPoints,
+            labelList& snappedPoint
+        ) const;
+
+        //- Generate single point by interpolation or snapping
+        template<class Type>
+        Type generatePoint
+        (
+            const DynamicList<Type>& snappedPoints,
+
             const scalar s0,
+            const Type& p0,
+            const label p0Index,
+
             const scalar s1,
+            const Type& p1,
+            const label p1Index
+        ) const;
+
+        template<class Type>
+        void 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,
 
-            const T& p0,
-            const T& p1,
-            const T& p2,
-            const T& p3,
+            DynamicList<Type>& points
+        ) const;
 
-            DynamicList<T>& points
+        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
+        (
+            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
         );
 
 public:
@@ -113,13 +330,16 @@ 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 scalarField& cellValues,
-            const scalarField& pointValues,
-            const scalar iso
+            const volScalarField& cellIsoVals,
+            const scalarField& pointIsoVals,
+            const scalar iso,
+            const bool regularise,
+            const scalar mergeTol = 1E-6    // fraction of bounding box
         );
 
 
@@ -132,27 +352,22 @@ public:
         }
 
         //- For every unmerged triangle point the point in the triSurface
-        const labelList triPointMergeMap() const
+        const labelList& triPointMergeMap() const
         {
             return triPointMergeMap_;
         }
 
-
-        //- sample field on faces
+        //- Interpolates cCoords,pCoords. Takes the original fields
+        //  used to create the iso surface.
         template <class Type>
-        tmp<Field<Type> > sample
+        tmp<Field<Type> > interpolate
         (
-            const Field<Type>& sampleCellValues
+            const volScalarField& cellIsoVals,
+            const scalarField& pointIsoVals,
+            const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+            const Field<Type>& pCoords
         ) const;
 
-        //- interpolate field to points
-        template <class Type>
-        tmp<Field<Type> >
-        interpolate
-        (
-            const Field<Type>& sampleCellValues,
-            const Field<Type>& samplePointValues
-        ) const;
 };
 
 

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

@@ -26,69 +26,90 @@ License
 
 #include "isoSurface.H"
 #include "polyMesh.H"
-#include "tetMatcher.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class Type>
-Type Foam::isoSurface::vertexInterp
+Type Foam::isoSurface::generatePoint
 (
-    const scalar iso,
-    const Type& p0,
-    const Type& p1,
+    const DynamicList<Type>& snappedPoints,
+
     const scalar s0,
-    const scalar s1
-)
+    const Type& p0,
+    const label p0Index,
+
+    const scalar s1,
+    const Type& p1,
+    const label p1Index
+) const
 {
     scalar d = s1-s0;
 
     if (mag(d) > VSMALL)
     {
-        return (iso-s0)/d*p1 + (s1-iso)/d*p0;
+        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
     {
-        return 0.5*(p0+p1);
+        scalar s = 0.4999;
+
+        return s*p1 + (1.0-s)*p0;
     }
 }
 
 
-// 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
+void Foam::isoSurface::generateTriPoints
 (
-    const scalar iso,
-    const scalar s0,
-    const scalar s1,
-    const scalar s2,
-    const scalar s3,
+    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)
+    if (s0 < iso_)
     {
         triIndex |= 1;
     }
-    if (s1 < iso)
+    if (s1 < iso_)
     {
         triIndex |= 2;
     }
-    if (s2 < iso)
+    if (s2 < iso_)
     {
         triIndex |= 4;
     }
-    if (s3 < iso)
+    if (s3 < iso_)
     {
         triIndex |= 8;
     }
@@ -102,29 +123,29 @@ void Foam::isoSurface::vertexInterp
 
         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));
+            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(vertexInterp(iso,p1,p0,s1,s0));
-            points.append(vertexInterp(iso,p1,p3,s1,s3));
-            points.append(vertexInterp(iso,p1,p2,s1,s2));
+            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:
         {
-            const Type tp1 = vertexInterp(iso,p0,p2,s0,s2);
-            const Type tp2 = vertexInterp(iso,p1,p3,s1,s3);
+            Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+            Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
 
-            points.append(vertexInterp(iso,p0,p3,s0,s3));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
             points.append(tp1);
             points.append(tp2);
             points.append(tp2);
-            points.append(vertexInterp(iso,p1,p2,s1,s2));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
             points.append(tp1);
         }
         break;
@@ -132,23 +153,23 @@ void Foam::isoSurface::vertexInterp
         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));
+            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:
         {
-            const Type tp0 = vertexInterp(iso,p0,p1,s0,s1);
-            const Type tp1 = vertexInterp(iso,p2,p3,s2,s3);
+            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(vertexInterp(iso,p0,p3,s0,s3));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
             points.append(tp0);
-            points.append(vertexInterp(iso,p1,p2,s1,s2));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
             points.append(tp1);
         }
         break;
@@ -156,133 +177,305 @@ void Foam::isoSurface::vertexInterp
         case 0x09:
         case 0x06:
         {
-            const Type tp0 = vertexInterp(iso,p0,p1,s0,s1);
-            const Type tp1 = vertexInterp(iso,p2,p3,s2,s3);
+            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(vertexInterp(iso,p1,p3,s1,s3));
+            points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
             points.append(tp1);
             points.append(tp0);
-            points.append(vertexInterp(iso,p0,p2,s0,s2));
+            points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
             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));
+            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::tmp<Foam::Field<Type> >
-Foam::isoSurface::sample(const Field<Type>& vField) const
+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
 {
-    return tmp<Field<Type> >(new Field<Type>(vField, meshCells()));
+    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 Field<Type>& sampleCellValues,
-    const Field<Type>& samplePointValues
+    const volScalarField& cVals,
+    const scalarField& pVals,
+    const GeometricField<Type, fvPatchField, volMesh>& cCoords,
+    const Field<Type>& pCoords
 ) const
 {
-    tetMatcher tet;
+    DynamicList<Type> triPoints(nCutCells_);
+    DynamicList<label> triMeshCells(nCutCells_);
 
-    DynamicList<Type> triValues;
+    // Dummy snap data
+    DynamicList<Type> snappedPoints;
+    labelList snappedCc(mesh_.nCells(), -1);
+    labelList snappedPoint(mesh_.nPoints(), -1);
 
-    // Note: in same order as construction of triSurface
-    label oldCellI = -1;
-    forAll(meshCells_, triI)
-    {
-        label cellI = meshCells_[triI];
+    generateTriPoints
+    (
+        cVals,
+        pVals,
 
-        if (cellI != oldCellI)
-        {
-            oldCellI = cellI;
+        cCoords,
+        pCoords,
 
-            const cell& cFaces = mesh_.cells()[cellI];
+        snappedPoints,
+        snappedCc,
+        snappedPoint,
 
-            if (tet.isA(mesh_, cellI))
-            {
-                // For tets don't do cell-centre decomposition, just use the
-                // tet points and values
+        triPoints,
+        triMeshCells
+    );
 
-                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)
+    forAll(triPoints, i)
     {
-        values[triPointMergeMap_[i]] = triValues[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,116 +143,90 @@ 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 (average_)
         {
-            if (debug)
+            //- From point field and interpolated cell.
+            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);
             {
-                Info<< "sampledIsoSurface::createGeometry() : lookup "
-                    << isoField_ << endl;
-            }
+                for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
+                {
+                    const labelList& pCells = fvm.pointCells(pointI);
 
-            cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_);
+                    forAll(pCells, i)
+                    {
+                        label cellI = pCells[i];
+
+                        cellAvg[cellI] += (*pointFieldPtr_)[pointI];
+                        nPointCells[cellI]++;
+                    }
+                }
+            }
+            forAll(cellAvg, cellI)
+            {
+                cellAvg[cellI] /= nPointCells[cellI];
+            }
+            // Give value to calculatedFvPatchFields
+            cellAvg.correctBoundaryConditions();
+
+            surfPtr_.reset
+            (
+                new isoSurface
+                (
+                    cellAvg,
+                    *pointFieldPtr_,
+                    isoVal_,
+                    regularise_
+                )
+            );
         }
         else
         {
-            // Bit of a hack. Read field and store.
-
-            if (debug)
-            {
-                Info<< "sampledIsoSurface::createGeometry() : reading "
-                    << isoField_ << " from time " <<fvm.time().timeName()
-                    << endl;
-            }
-
-            readFieldPtr_.reset
+            //- Direct from cell field and point field.
+            surfPtr_.reset
             (
-                new volScalarField
+                new isoSurface
                 (
-                    IOobject
-                    (
-                        isoField_,
-                        fvm.time().timeName(),
-                        fvm,
-                        IOobject::MUST_READ,
-                        IOobject::NO_WRITE,
-                        false
-                    ),
-                    fvm
+                    *volFieldPtr_,
+                    *pointFieldPtr_,
+                    isoVal_,
+                    regularise_
                 )
             );
-
-            cellFldPtr = readFieldPtr_.operator->();
-        }
-        const volScalarField& cellFld = *cellFldPtr;
-
-        tmp<pointScalarField> pointFld
-        (
-            volPointInterpolation::New(fvm).interpolate(cellFld)
-        );
-
-        //- Direct from cell field and point field. Gives bad continuity.
-        //const isoSurface iso
-        //(
-        //    fvm,
-        //    cellFld.internalField(),
-        //    pointFld().internalField(),
-        //    isoVal_
-        //);
-
-        //- 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 isoSurface iso
-        (
-            fvm,
-            cellAvg,
-            pointFld().internalField(),
-            isoVal_
-        );
-
-
-        const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
-        meshCells_ = iso.meshCells();
-        triPointMergeMap_ = iso.triPointMergeMap();
 
         if (debug)
         {
             Pout<< "sampledIsoSurface::createGeometry() : constructed iso:"
                 << nl
+                << "    regularise     : " << regularise_ << nl
+                << "    average        : " << average_ << nl
                 << "    isoField       : " << isoField_ << nl
                 << "    isoValue       : " << isoVal_ << nl
-                << "    unmerged points: " << triPointMergeMap_.size() << nl
                 << "    points         : " << points().size() << nl
-                << "    tris           : " << triSurface::size() << nl
-                << "    cut cells      : " << meshCells_.size() << endl;
+                << "    tris           : " << surface().size() << nl
+                << "    cut cells      : " << surface().meshCells().size()
+                << endl;
         }
     }
 }
@@ -176,12 +244,27 @@ Foam::sampledIsoSurface::sampledIsoSurface
     sampledSurface(name, mesh, dict),
     isoField_(dict.lookup("isoField")),
     isoVal_(readScalar(dict.lookup("isoValue"))),
+    regularise_(dict.lookupOrDefault("regularise", true)),
+    average_(dict.lookupOrDefault("average", false)),
     zoneName_(word::null),
+    surfPtr_(NULL),
     facesPtr_(NULL),
     storedTimeIndex_(-1),
-    meshCells_(0),
-    triPointMergeMap_(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_))
 //    {
@@ -209,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();
     }
 }
@@ -317,9 +401,9 @@ void Foam::sampledIsoSurface::print(Ostream& os) const
 {
     os  << "sampledIsoSurface: " << name() << " :"
         << "  field:" << isoField_
-        << "  value:" << isoVal_
-        << "  faces:" << faces().size()
-        << "  points:" << points().size();
+        << "  value:" << isoVal_;
+        //<< "  faces:" << faces().size()       // note: possibly no geom yet
+        //<< "  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,20 +52,27 @@ namespace Foam
 
 class sampledIsoSurface
 :
-    public sampledSurface,
-    public triSurface
+    public sampledSurface
 {
     // Private data
 
         //- Field to get isoSurface of
-        word isoField_;
+        const word isoField_;
 
         //- iso value
-        scalar isoVal_;
+        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_;
 
+        mutable autoPtr<isoSurface> surfPtr_;
+
         //- triangles converted to faceList
         mutable autoPtr<faceList> facesPtr_;
 
@@ -75,15 +82,21 @@ class sampledIsoSurface
             //- Time at last call
             mutable label storedTimeIndex_;
 
-            //- For every triangle the original cell in mesh
-            mutable labelList meshCells_;
+            //- Cached volfield
+            mutable autoPtr<volScalarField> storedVolFieldPtr_;
+            mutable const volScalarField* volFieldPtr_;
+
+            //- Cached pointfield
+            mutable autoPtr<pointScalarField> storedPointFieldPtr_;
+            mutable const pointScalarField* pointFieldPtr_;
 
-            //- For every unmerged triangle point the point in the triSurface
-            mutable labelList triPointMergeMap_;
 
 
     // Private Member Functions
 
+        //- Get fields needed to recreate iso surface.
+        void getIsoFields() const;
+
         //- Create iso surface (if time has changed)
         void createGeometry() const;
 
@@ -127,7 +140,7 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return triSurface::points();
+            return surface().points();
         }
 
         //- Faces of surface
@@ -135,7 +148,7 @@ public:
         {
             if (!facesPtr_.valid())
             {
-                const triSurface& s = *this;
+                const triSurface& s = surface();
 
                 facesPtr_.reset(new faceList(s.size()));
 
@@ -147,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, meshCells_));
+    return tmp<Field<Type> >(new Field<Type>(vField, surface().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
-    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;
+    // Sample.
+    return surface().interpolate
+    (
+        *volFieldPtr_,
+        *pointFieldPtr_,
+        volFld,
+        pointFld()
+    );
 }