Merge commit 'OpenCFD/master' into olesenm

applications/utilities/mesh/manipulation/cellSet/cellSetDict

@@ -82,7 +82,7 @@ topoSetSources
     // Cells in cell zone
     zoneToCell
     {
-        name cellZone;      // name of cellZone
+        name ".*Zone";      // name of cellZone, wildcards allowed.
     }
 
     // values of field within certain range

applications/utilities/mesh/manipulation/faceSet/faceSetDict

@@ -56,13 +56,13 @@ topoSetSources
     // All faces of patch
     patchToFace
     {
-        name movingWall;
+        name ".*Wall";      // Name of patch, regular expressions allowed
     }
 
     // All faces of faceZone
     zoneToFace
     {
-        name faceZone1;     // Name of faceZone
+        name ".*Zone1";     // Name of faceZone, regular expressions allowed
     }
 
     // Faces with face centre within box

applications/utilities/mesh/manipulation/pointSet/pointSetDict

@@ -57,6 +57,12 @@ topoSetSources
        box   (0 0 0) (1 1 1);
     }
 
+    // All points in pointzone
+    zoneToPoint
+    {
+        name ".*Zone";      // name of pointZone, wildcards allowed.
+    }
+
     // Select based on surface
     surfaceToPoint
     {

src/meshTools/Make/files

@@ -55,6 +55,7 @@ indexedOctree/treeDataTriSurface.C
 searchableSurface = searchableSurface
 $(searchableSurface)/distributedTriSurfaceMesh.C
 $(searchableSurface)/searchableBox.C
+$(searchableSurface)/searchablePlane.C
 $(searchableSurface)/searchablePlate.C
 $(searchableSurface)/searchableSphere.C
 $(searchableSurface)/searchableSurface.C

src/meshTools/searchableSurface/searchablePlane.C

@@ -0,0 +1,231 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "searchablePlane.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SortableList.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchablePlane, 0);
+addToRunTimeSelectionTable(searchableSurface, searchablePlane, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::pointIndexHit Foam::searchablePlane::findLine
+(
+    const point& start,
+    const point& end
+) const
+{
+    pointIndexHit info(true, vector::zero, 0);
+
+    linePointRef l(start, end);
+
+    scalar t = lineIntersect(l);
+
+    if (t < 0 || t > 1)
+    {
+        info.setMiss();
+        info.setIndex(-1);
+    }
+    else
+    {
+        info.setPoint(start+t*l.vec());
+    }
+
+    return info;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchablePlane::searchablePlane
+(
+    const IOobject& io,
+    const point& basePoint,
+    const vector& normal
+)
+:
+    searchableSurface(io),
+    plane(basePoint, normal)
+{}
+
+
+Foam::searchablePlane::searchablePlane
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+:
+    searchableSurface(io),
+    plane(dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchablePlane::~searchablePlane()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchablePlane::regions() const
+{
+    if (regions_.empty())
+    {
+        regions_.setSize(1);
+        regions_[0] = "region0";
+    }
+    return regions_;
+}
+
+
+void Foam::searchablePlane::findNearest
+(
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(samples.size());
+
+    forAll(samples, i)
+    {
+        info[i].setPoint(nearestPoint(samples[i]));
+
+        if (magSqr(samples[i]-info[i].rawPoint()) > nearestDistSqr[i])
+        {
+            info[i].setIndex(-1);
+            info[i].setMiss();
+        }
+        else
+        {
+            info[i].setIndex(0);
+            info[i].setHit();
+        }
+    }
+}
+
+
+void Foam::searchablePlane::findLine
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    info.setSize(start.size());
+
+    forAll(start, i)
+    {
+        info[i] = findLine(start[i], end[i]);
+    }
+}
+
+
+void Foam::searchablePlane::findLineAny
+(
+    const pointField& start,
+    const pointField& end,
+    List<pointIndexHit>& info
+) const
+{
+    findLine(start, end, info);
+}
+
+
+void Foam::searchablePlane::findLineAll
+(
+    const pointField& start,
+    const pointField& end,
+    List<List<pointIndexHit> >& info
+) const
+{
+    List<pointIndexHit> nearestInfo;
+    findLine(start, end, nearestInfo);
+
+    info.setSize(start.size());
+    forAll(info, pointI)
+    {
+        if (nearestInfo[pointI].hit())
+        {
+            info[pointI].setSize(1);
+            info[pointI][0] = nearestInfo[pointI];
+        }
+        else
+        {
+            info[pointI].clear();
+        }
+    }
+}
+
+
+void Foam::searchablePlane::getRegion
+(
+    const List<pointIndexHit>& info,
+    labelList& region
+) const
+{
+    region.setSize(info.size());
+    region = 0;
+}
+
+
+void Foam::searchablePlane::getNormal
+(
+    const List<pointIndexHit>& info,
+    vectorField& n
+) const
+{
+    n.setSize(info.size());
+    n = normal();
+}
+
+
+void Foam::searchablePlane::getVolumeType
+(
+    const pointField& points,
+    List<volumeType>& volType
+) const
+{
+    FatalErrorIn
+    (
+        "searchableCollection::getVolumeType(const pointField&"
+        ", List<volumeType>&) const"
+    )   << "Volume type not supported for plane."
+        << exit(FatalError);
+}
+
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchablePlane.H

@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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::searchablePlane
+
+Description
+    Searching on plane. See plane.H
+
+SourceFiles
+    searchablePlane.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchablePlane_H
+#define searchablePlane_H
+
+#include "searchableSurface.H"
+#include "plane.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+/*---------------------------------------------------------------------------*\
+                           Class searchablePlane Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchablePlane
+:
+    public searchableSurface,
+    public plane
+{
+private:
+
+    // Private Member Data
+
+        mutable wordList regions_;
+
+
+    // Private Member Functions
+
+        pointIndexHit findLine
+        (
+            const point& start,
+            const point& end
+        ) const;
+
+
+        //- Disallow default bitwise copy construct
+        searchablePlane(const searchablePlane&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchablePlane&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchablePlane");
+
+
+    // Constructors
+
+        //- Construct from components
+        searchablePlane
+        (
+            const IOobject& io,
+            const point& basePoint,
+            const vector& normal
+        );
+
+        //- Construct from dictionary (used by searchableSurface)
+        searchablePlane
+        (
+            const IOobject& io,
+            const dictionary& dict
+        );
+
+    // Destructor
+
+        virtual ~searchablePlane();
+
+
+    // Member Functions
+
+        virtual const wordList& regions() const;
+
+        //- Whether supports volume type below
+        virtual bool hasVolumeType() const
+        {
+            return false;
+        }
+
+        //- Range of local indices that can be returned.
+        virtual label size() const
+        {
+            return 1;
+        }
+
+
+        // Multiple point queries.
+
+            virtual void findNearest
+            (
+                const pointField& sample,
+                const scalarField& nearestDistSqr,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLine
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            virtual void findLineAny
+            (
+                const pointField& start,
+                const pointField& end,
+                List<pointIndexHit>&
+            ) const;
+
+            //- Get all intersections in order from start to end.
+            virtual void findLineAll
+            (
+                const pointField& start,
+                const pointField& end,
+                List<List<pointIndexHit> >&
+            ) const;
+
+            //- From a set of points and indices get the region
+            virtual void getRegion
+            (
+                const List<pointIndexHit>&,
+                labelList& region
+            ) const;
+
+            //- From a set of points and indices get the normal
+            virtual void getNormal
+            (
+                const List<pointIndexHit>&,
+                vectorField& normal
+            ) const;
+
+            //- Determine type (inside/outside/mixed) for point. unknown if
+            //  cannot be determined (e.g. non-manifold surface)
+            virtual void getVolumeType
+            (
+                const pointField&,
+                List<volumeType>&
+            ) const;
+
+
+        // regIOobject implementation
+
+            bool writeData(Ostream&) const
+            {
+                notImplemented("searchablePlane::writeData(Ostream&) const");
+                return false;
+            }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/searchableSurface/searchablePlate.C

@@ -218,7 +218,16 @@ Foam::searchablePlate::searchablePlate
     origin_(origin),
     span_(span),
     normalDir_(calcNormal(span_))
-{}
+{
+    if (debug)
+    {
+        Info<< "searchablePlate::searchablePlate :"
+            << " origin:" << origin_
+            << " origin+span:" << origin_+span_
+            << " normal:" << vector::componentNames[normalDir_]
+            << endl;
+    }
+}
 
 
 Foam::searchablePlate::searchablePlate
@@ -231,7 +240,16 @@ Foam::searchablePlate::searchablePlate
     origin_(dict.lookup("origin")),
     span_(dict.lookup("span")),
     normalDir_(calcNormal(span_))
-{}
+{
+    if (debug)
+    {
+        Info<< "searchablePlate::searchablePlate :"
+            << " origin:" << origin_
+            << " origin+span:" << origin_+span_
+            << " normal:" << vector::componentNames[normalDir_]
+            << endl;
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

src/meshTools/searchableSurface/searchablePlate.H

@@ -27,7 +27,14 @@ Class
 
 Description
     Searching on finite plate. Plate has to be aligned with coordinate
-    axes!
+    axes.
+    Plate defined as origin and span. One of the components of span has
+    to be 0 which defines the normal direction. E.g.
+
+    span    = (Sx Sy 0)     // plate in x-y plane
+    origin  = (Ox Oy Oz)
+
+    now plane is from (Ox Oy Oz) to (Ox+Sx Oy+Sy Oz)
 
 SourceFiles
     searchablePlate.C

src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C

@@ -22,8 +22,6 @@ 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 "nbrToCell.H"
@@ -58,22 +56,46 @@ Foam::topoSetSource::addToUsageTable Foam::nbrToCell::usage_
 void Foam::nbrToCell::combine(topoSet& set, const bool add) const
 {
     const cellList& cells = mesh().cells();
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    boolList isCoupled(mesh_.nFaces()-mesh_.nInternalFaces(), false);
+
+    forAll(patches, patchI)
+    {   
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            forAll(pp, i)
+            {
+                isCoupled[faceI-mesh_.nInternalFaces()] = true;
+                faceI++;
+            }
+        }
+    }
 
     forAll(cells, cellI)
     {
-        const cell& cll = cells[cellI];
+        const cell& cFaces = cells[cellI];
 
-        label nInternalFaces = 0;
+        label nNbrCells = 0;
 
-        forAll(cll, i)
+        forAll(cFaces, i)
         {
-            if (mesh().isInternalFace(cll[i]))
+            label faceI = cFaces[i];
+
+            if (mesh_.isInternalFace(faceI))
             {
-                nInternalFaces++;
+                nNbrCells++;
+            }
+            else if (isCoupled[faceI-mesh_.nInternalFaces()])
+            {
+                nNbrCells++;
             }
         }
 
-        if (nInternalFaces <= minNbrs_)
+        if (nNbrCells <= minNbrs_)
         {
             addOrDelete(set, cellI, add);
         }

src/meshTools/sets/cellSources/nbrToCell/nbrToCell.H

@@ -27,7 +27,7 @@ Class
 
 Description
     A topoSetSource to select cells based on number of neighbouring cells
-    (i.e. number of internal faces)
+    (i.e. number of internal or coupled faces)
 
 SourceFiles
     nbrToCell.C

src/meshTools/sets/faceSources/cellToFace/cellToFace.C

@@ -22,14 +22,13 @@ 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 "cellToFace.H"
 #include "polyMesh.H"
 #include "cellSet.H"
 #include "Time.H"
+#include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -102,48 +101,58 @@ void Foam::cellToFace::combine(topoSet& set, const bool add) const
     {
         // Add all faces whose both neighbours are in set.
 
-        // Count number of cells using face.
+        label nInt = mesh_.nInternalFaces();
-        Map<label> numCells(loadedSet.size());
+        const labelList& own = mesh_.faceOwner();
+        const labelList& nei = mesh_.faceNeighbour();
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-        for
+
-        (
+        // Check all internal faces
-            cellSet::const_iterator iter = loadedSet.begin();
+        for (label faceI = 0; faceI < nInt; faceI++)
-            iter != loadedSet.end();
-            ++iter
-        )
         {
-            label cellI = iter.key();
+            if (loadedSet.found(own[faceI]) && loadedSet.found(nei[faceI]))
-
-            const labelList& cFaces = mesh_.cells()[cellI];
-
-            forAll(cFaces, cFaceI)
             {
-                label faceI = cFaces[cFaceI];
+                addOrDelete(set, faceI, add);
-
-                Map<label>::iterator fndFace = numCells.find(faceI);
-
-                if (fndFace == numCells.end())
-                {
-                    numCells.insert(faceI, 1);
-                }
-                else
-                {
-                    fndFace()++;
-                }
             }
         }
 
-        // Include faces that are referenced twice
+
-        for
+        // Get coupled cell status
-        (
+        boolList neiInSet(mesh_.nFaces()-nInt, false);
-            Map<label>::const_iterator iter = numCells.begin();
+
-            iter != numCells.end();
+        forAll(patches, patchI)
-            ++iter
-        )
         {
-            if (iter() == 2)
+            const polyPatch& pp = patches[patchI];
+
+            if (pp.coupled())
             {
-                addOrDelete(set, iter.key(), add);
+                label faceI = pp.start();
+                forAll(pp, i)
+                {
+                    neiInSet[faceI-nInt] = loadedSet.found(own[faceI]);
+                    faceI++;
+                }
+            }
+        }
+        syncTools::swapBoundaryFaceList(mesh_, neiInSet, false);
+
+
+        // Check all boundary faces
+        forAll(patches, patchI)
+        {
+            const polyPatch& pp = patches[patchI];
+
+            if (pp.coupled())
+            {
+                label faceI = pp.start();
+                forAll(pp, i)
+                {
+                    if (loadedSet.found(own[faceI]) && neiInSet[faceI-nInt])
+                    {
+                        addOrDelete(set, faceI, add);
+                    }
+                    faceI++;
+                }
             }
         }
     }

src/meshTools/sets/faceSources/cellToFace/cellToFace.H

@@ -30,6 +30,7 @@ Description
 
     Either all faces of cell or some other criterion.
     See implementation.
+    Note: when picking up coupled faces uses cells on neighbouring processors.
 
 SourceFiles
     cellToFace.C

src/sampling/sampledSurface/distanceSurface/distanceSurface.C

@@ -45,6 +45,11 @@ namespace Foam
 
 void Foam::distanceSurface::createGeometry()
 {
+    if (debug)
+    {
+        Pout<< "distanceSurface::createGeometry :updating geometry." << endl;
+    }
+
     // Clear any stored topologies
     facesPtr_.clear();
 
@@ -67,7 +72,7 @@ void Foam::distanceSurface::createGeometry()
                 false
             ),
             fvm,
-            dimensionedScalar("zero", dimless/dimTime, 0)
+            dimensionedScalar("zero", dimLength, 0)
         )
     );
     volScalarField& cellDistance = cellDistancePtr_();
@@ -157,6 +162,7 @@ void Foam::distanceSurface::createGeometry()
         }
     }
 
+
     // On processor patches the mesh.C() will already be the cell centre
     // on the opposite side so no need to swap cellDistance.
 
@@ -164,23 +170,70 @@ void Foam::distanceSurface::createGeometry()
     // Distance to points
     pointDistance_.setSize(fvm.nPoints());
     {
+        const pointField& pts = fvm.points();
+
         List<pointIndexHit> nearest;
         surfPtr_().findNearest
         (
-            fvm.points(),
+            pts,
-            scalarField(fvm.nPoints(), GREAT),
+            scalarField(pts.size(), GREAT),
             nearest
         );
-        forAll(pointDistance_, pointI)
+
+        if (signed_)
         {
-            pointDistance_[pointI] = Foam::mag
+            vectorField normal;
-            (
+            surfPtr_().getNormal(nearest, normal);
-                nearest[pointI].hitPoint()
+
-              - fvm.points()[pointI]
+            forAll(nearest, i)
-            );
+            {
+                vector d(pts[i]-nearest[i].hitPoint());
+
+                if ((d&normal[i]) > 0)
+                {
+                    pointDistance_[i] = Foam::mag(d);
+                }
+                else
+                {
+                    pointDistance_[i] = -Foam::mag(d);
+                }
+            }
+        }
+        else
+        {
+            forAll(nearest, i)
+            {
+                pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint());
+            }
         }
     }
 
+
+    if (debug)
+    {
+        Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
+        cellDistance.write();
+        pointScalarField pDist
+        (
+            IOobject
+            (
+                "pointDistance",
+                fvm.time().timeName(),
+                fvm.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            pointMesh::New(fvm),
+            dimensionedScalar("zero", dimLength, 0)
+        );
+        pDist.internalField() = pointDistance_;
+
+        Pout<< "Writing point distance:" << pDist.objectPath() << endl;
+        pDist.write();
+    }
+
+
     //- Direct from cell field and point field.
     isoSurfPtr_.reset
     (
@@ -196,6 +249,7 @@ void Foam::distanceSurface::createGeometry()
     if (debug)
     {
         print(Pout);
+        Pout<< endl;
     }
 }
 
@@ -264,6 +318,13 @@ bool Foam::distanceSurface::needsUpdate() const
 
 bool Foam::distanceSurface::expire()
 {
+    if (debug)
+    {
+        Pout<< "distanceSurface::expire :"
+            << " have-facesPtr_:" << facesPtr_.valid()
+            << " needsUpdate_:" << needsUpdate_ << endl;
+    }
+
     // Clear any stored topologies
     facesPtr_.clear();
 
@@ -280,6 +341,13 @@ bool Foam::distanceSurface::expire()
 
 bool Foam::distanceSurface::update()
 {
+    if (debug)
+    {
+        Pout<< "distanceSurface::update :"
+            << " have-facesPtr_:" << facesPtr_.valid()
+            << " needsUpdate_:" << needsUpdate_ << endl;
+    }
+
     if (!needsUpdate_)
     {
         return false;

src/sampling/sampledSurface/isoSurface/isoSurface.C

@@ -58,6 +58,31 @@ Foam::scalar Foam::isoSurface::isoFraction
 }
 
 
+bool Foam::isoSurface::isEdgeOfFaceCut
+(
+    const scalarField& pVals,
+    const face& f,
+    const bool ownLower,
+    const bool neiLower
+) const
+{
+    forAll(f, fp)
+    {
+        bool fpLower = (pVals[f[fp]] < iso_);
+        if
+        (
+            (fpLower != ownLower)
+         || (fpLower != neiLower)
+         || (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
+        )
+        {
+            return true;
+        }
+    }
+    return false;
+}
+
+
 // Determine for every face/cell whether it (possibly) generates triangles.
 void Foam::isoSurface::calcCutTypes
 (
@@ -84,22 +109,13 @@ void Foam::isoSurface::calcCutTypes
         }
         else
         {
-            // Mesh edge.
+            // See if any mesh edge is cut by looping over all the edges of the
+            // face.
             const face f = mesh_.faces()[faceI];
 
-            forAll(f, fp)
+            if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
             {
-                bool fpLower = (pVals[f[fp]] < iso_);
+                faceCutType_[faceI] = CUT;
-                if
-                (
-                    (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
-                 || (fpLower != ownLower)
-                 || (fpLower != neiLower)
-                )
-                {
-                    faceCutType_[faceI] = CUT;
-                    break;
-                }
             }
         }
     }
@@ -117,22 +133,13 @@ void Foam::isoSurface::calcCutTypes
             {
                 bool ownLower = (cVals[own[faceI]] < iso_);
 
-                // Mesh edge.
                 const face f = mesh_.faces()[faceI];
 
-                forAll(f, fp)
+                if (isEdgeOfFaceCut(pVals, f, ownLower, ownLower))
                 {
-                    bool fpLower = (pVals[f[fp]] < iso_);
+                    faceCutType_[faceI] = CUT;
-                    if
-                    (
-                        (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
-                     || (fpLower != ownLower)
-                    )
-                    {
-                        faceCutType_[faceI] = CUT;
-                        break;
-                    }
                 }
+
                 faceI++;
             }
         }
@@ -152,19 +159,9 @@ void Foam::isoSurface::calcCutTypes
                     // Mesh edge.
                     const face f = mesh_.faces()[faceI];
 
-                    forAll(f, fp)
+                    if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower))
                     {
-                        bool fpLower = (pVals[f[fp]] < iso_);
+                        faceCutType_[faceI] = CUT;
-                        if
-                        (
-                            (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_))
-                         || (fpLower != ownLower)
-                         || (fpLower != neiLower)
-                        )
-                        {
-                            faceCutType_[faceI] = CUT;
-                            break;
-                        }
                     }
                 }
                 faceI++;
@@ -1355,6 +1352,30 @@ Foam::isoSurface::isoSurface
     iso_(iso),
     mergeDistance_(mergeTol*mesh_.bounds().mag())
 {
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const labelList& own = mesh_.faceOwner();
+
+    // Check
+    forAll(patches, patchI)
+    {
+        if (isA<emptyPolyPatch>(patches[patchI]))
+        {
+            FatalErrorIn
+            (
+                "isoSurface::isoSurface\n"
+                "(\n"
+                "    const volScalarField& cVals,\n"
+                "    const scalarField& pVals,\n"
+                "    const scalar iso,\n"
+                "    const bool regularise,\n"
+                "    const scalar mergeTol\n"
+                ")\n"
+            )   << "Iso surfaces not supported on case with empty patches."
+                << exit(FatalError);
+        }
+    }
+
+
     // Determine if any cut through face/cell
     calcCutTypes(cVals, pVals);
 
@@ -1364,8 +1385,6 @@ Foam::isoSurface::isoSurface
     PackedList<1> isBoundaryPoint(mesh_.nPoints());
 
     labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    const labelList& own = mesh_.faceOwner();
 
     forAll(patches, patchI)
     {

src/sampling/sampledSurface/isoSurface/isoSurface.H

@@ -31,6 +31,9 @@ Description
     G.M. Treece, R.W. Prager and A.H. Gee.
 
     Note:
+    - not possible on patches of type 'empty'. There are no values on
+    'empty' patch fields so even the api would have to change 
+     (no volScalarField as argument). Too messy.
     - 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!
@@ -123,6 +126,15 @@ class isoSurface
             const scalar s1
         ) const;
 
+        //- Check if any edge of a face is cut
+        bool isEdgeOfFaceCut
+        (
+            const scalarField& pVals,
+            const face& f,
+            const bool ownLower,
+            const bool neiLower
+        ) const;
+
         //- Set faceCutType,cellCutType.
         void calcCutTypes
         (
@@ -217,7 +229,7 @@ class isoSurface
         ) const;
 
         template<class Type>
-        label generateTriPoints
+        label generateFaceTriPoints
         (
             const volScalarField& cVals,
             const scalarField& pVals,

src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C

@@ -200,7 +200,7 @@ void Foam::isoSurface::generateTriPoints
 
 
 template<class Type>
-Foam::label Foam::isoSurface::generateTriPoints
+Foam::label Foam::isoSurface::generateFaceTriPoints
 (
     const volScalarField& cVals,
     const scalarField& pVals,
@@ -288,6 +288,19 @@ void Foam::isoSurface::generateTriPoints
     const labelList& own = mesh_.faceOwner();
     const labelList& nei = mesh_.faceNeighbour();
 
+    if
+    (
+        (cVals.size() != mesh_.nCells())
+     || (pVals.size() != mesh_.nPoints())
+     || (cCoords.size() != mesh_.nCells())
+     || (pCoords.size() != mesh_.nPoints())
+     || (snappedCc.size() != mesh_.nCells())
+     || (snappedPoint.size() != mesh_.nPoints())
+    )
+    {
+        FatalErrorIn("isoSurface::generateTriPoints(..)")
+            << "Incorrect size." << abort(FatalError);
+    }
 
     // Determine neighbouring snap status
     labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
@@ -319,7 +332,7 @@ void Foam::isoSurface::generateTriPoints
     {
         if (faceCutType_[faceI] != NOTCUT)
         {
-            generateTriPoints
+            generateFaceTriPoints
             (
                 cVals,
                 pVals,
@@ -357,7 +370,7 @@ void Foam::isoSurface::generateTriPoints
                 {
                     if (faceCutType_[faceI] != NOTCUT)
                     {
-                        generateTriPoints
+                        generateFaceTriPoints
                         (
                             cVals,
                             pVals,
@@ -384,14 +397,14 @@ void Foam::isoSurface::generateTriPoints
         }
         else if (isA<emptyPolyPatch>(pp))
         {
-            // Assume zero-gradient.
+            // Assume zero-gradient. But what about coordinates?
             label faceI = pp.start();
 
             forAll(pp, i)
             {
                 if (faceCutType_[faceI] != NOTCUT)
                 {
-                    generateTriPoints
+                    generateFaceTriPoints
                     (
                         cVals,
                         pVals,
@@ -423,7 +436,7 @@ void Foam::isoSurface::generateTriPoints
             {
                 if (faceCutType_[faceI] != NOTCUT)
                 {
-                    generateTriPoints
+                    generateFaceTriPoints
                     (
                         cVals,
                         pVals,