Merge commit 'OpenCFD/master' into merging

Conflicts:

	src/OpenFOAM/containers/misc/xferI.H
	src/OpenFOAM/fields/tmp/refCount.H
	src/OpenFOAM/memory/refCount/refCount.H
	src/OpenFOAM/memory/xfer/xferI.H
	src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.C
	src/finiteVolume/fvMesh/extendedStencil/faceStencil/faceEdgeCellStencil.H

applications/utilities/mesh/manipulation/createPatch/createPatch.C

@@ -625,10 +625,8 @@ int main(int argc, char *argv[])
                     << " as patch " << destPatchI
                     << " from " << patchDict << endl;
 
-                patchDict.remove("nFaces");
-                patchDict.add("nFaces", 0);
-                patchDict.remove("startFace");
-                patchDict.add("startFace", startFaceI);
+                patchDict.set("nFaces", 0);
+                patchDict.set("startFace", startFaceI);
 
                 // Add an empty patch.
                 allPatches.append

applications/utilities/parallelProcessing/reconstructParMesh/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

@@ -46,6 +46,8 @@ Description
 #include "mapAddedPolyMesh.H"
 #include "polyMeshAdder.H"
 #include "faceCoupleInfo.H"
+#include "fvMeshAdder.H"
+#include "polyTopoChange.H"
 
 using namespace Foam;
 
@@ -203,6 +205,78 @@ autoPtr<faceCoupleInfo> determineCoupledFaces
 }
 
 
+autoPtr<mapPolyMesh> mergeSharedPoints
+(
+    const scalar mergeDist,
+    polyMesh& mesh,
+    labelListList& pointProcAddressing
+)
+{
+    // Find out which sets of points get merged and create a map from
+    // mesh point to unique point.
+    Map<label> pointToMaster
+    (
+        fvMeshAdder::findSharedPoints
+        (
+            mesh,
+            mergeDist
+        )
+    );
+
+    Info<< "mergeSharedPoints : detected " << pointToMaster.size()
+        << " points that are to be merged." << endl;
+
+    if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
+    {
+        return autoPtr<mapPolyMesh>(NULL);
+    }
+
+    polyTopoChange meshMod(mesh);
+
+    fvMeshAdder::mergePoints(mesh, pointToMaster, meshMod);
+
+    // Change the mesh (no inflation). Note: parallel comms allowed.
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
+
+    // Update fields. No inflation, parallel sync.
+    mesh.updateMesh(map);
+
+    // pointProcAddressing give indices into the master mesh so adapt them
+    // for changed point numbering.
+
+    // Adapt constructMaps for merged points.
+    forAll(pointProcAddressing, procI)
+    {
+        labelList& constructMap = pointProcAddressing[procI];
+
+        forAll(constructMap, i)
+        {
+            label oldPointI = constructMap[i];
+
+            // New label of point after changeMesh.
+            label newPointI = map().reversePointMap()[oldPointI];
+
+            if (newPointI < -1)
+            {
+                constructMap[i] = -newPointI-2;
+            }
+            else if (newPointI >= 0)
+            {
+                constructMap[i] = newPointI;
+            }
+            else
+            {
+                FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
+                    << "Problem. oldPointI:" << oldPointI
+                    << " newPointI:" << newPointI << abort(FatalError);
+            }
+        }
+    }
+
+    return map;    
+}
+
+
 int main(int argc, char *argv[])
 {
     argList::noParallel();
@@ -214,7 +288,7 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 
-    Pout<< "This is an experimental tool which tries to merge"
+    Info<< "This is an experimental tool which tries to merge"
         << " individual processor" << nl
         << "meshes back into one master mesh. Use it if the original"
         << " master mesh has" << nl
@@ -246,7 +320,7 @@ int main(int argc, char *argv[])
     }
     scalar writeTol = Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
 
-    Pout<< "Merge tolerance : " << mergeTol << nl
+    Info<< "Merge tolerance : " << mergeTol << nl
         << "Write tolerance : " << writeTol << endl;
 
     if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
@@ -267,11 +341,11 @@ int main(int argc, char *argv[])
 
     if (fullMatch)
     {
-        Pout<< "Doing geometric matching on all boundary faces." << nl << endl;
+        Info<< "Doing geometric matching on all boundary faces." << nl << endl;
     }
     else
     {
-        Pout<< "Doing geometric matching on correct procBoundaries only."
+        Info<< "Doing geometric matching on correct procBoundaries only."
             << nl << "This assumes a correct decomposition." << endl;
     }
 
@@ -292,7 +366,7 @@ int main(int argc, char *argv[])
         nProcs++;
     }
 
-    Pout<< "Found " << nProcs << " processor directories" << nl << endl;
+    Info<< "Found " << nProcs << " processor directories" << nl << endl;
 
 
     // Read all databases.
@@ -300,7 +374,7 @@ int main(int argc, char *argv[])
 
     forAll (databases, procI)
     {
-        Pout<< "Reading database "
+        Info<< "Reading database "
             << args.caseName()/fileName(word("processor") + name(procI))
             << endl;
 
@@ -337,7 +411,7 @@ int main(int argc, char *argv[])
     }
 
     // Set master time
-    Pout<< "Setting master time to " << databases[0].timeName() << nl << endl;
+    Info<< "Setting master time to " << databases[0].timeName() << nl << endl;
     runTime.setTime(databases[0]);
 
 
@@ -373,7 +447,7 @@ int main(int argc, char *argv[])
                 << endl << exit(FatalError);
         }
 
-        Pout<< "Reading points from "
+        Info<< "Reading points from "
             << databases[procI].caseName()
             << " for time = " << databases[procI].timeName()
             << nl << endl;
@@ -403,7 +477,7 @@ int main(int argc, char *argv[])
     }
     const scalar mergeDist = mergeTol*mag(bb.max() - bb.min());
 
-    Pout<< "Overall mesh bounding box  : " << bb << nl
+    Info<< "Overall mesh bounding box  : " << bb << nl
         << "Relative tolerance         : " << mergeTol << nl
         << "Absolute matching distance : " << mergeDist << nl
         << endl;
@@ -422,7 +496,7 @@ int main(int argc, char *argv[])
 
     {
         // Construct empty mesh.
-        Pout<< "Constructing empty mesh to add to." << nl << endl;
+        Info<< "Constructing empty mesh to add to." << nl << endl;
         polyMesh masterMesh
         (
             IOobject
@@ -439,7 +513,7 @@ int main(int argc, char *argv[])
 
         for (label procI = 0; procI < nProcs; procI++)
         {
-            Pout<< "Reading mesh to add from "
+            Info<< "Reading mesh to add from "
                 << databases[procI].caseName()
                 << " for time = " << databases[procI].timeName()
                 << nl << endl;
@@ -475,7 +549,7 @@ int main(int argc, char *argv[])
 
 
             // Add elements to mesh
-            Pout<< "Adding to master mesh" << nl << endl;
+            Info<< "Adding to master mesh" << nl << endl;
 
             autoPtr<mapAddedPolyMesh> map = polyMeshAdder::add
             (
@@ -503,16 +577,19 @@ int main(int argc, char *argv[])
             renumber(map().addedPointMap(), pointProcAddressing[procI]);
             renumber(map().addedPatchMap(), boundaryProcAddressing[procI]);
 
-            Pout<< endl;
+            Info<< endl;
         }
 
+        // See if any points on the mastermesh have become connected
+        // because of connections through processor meshes.
+        mergeSharedPoints(mergeDist, masterMesh, pointProcAddressing);
 
         // Save some properties on the reconstructed mesh
         masterInternalFaces = masterMesh.nInternalFaces();
         masterOwner = masterMesh.faceOwner();
 
 
-        Pout<< "\nWriting merged mesh to "
+        Info<< "\nWriting merged mesh to "
             << runTime.path()/runTime.timeName()
             << nl << endl;
 
@@ -527,12 +604,12 @@ int main(int argc, char *argv[])
 
     // Write the addressing
 
-    Pout<< "Reconstructing the addressing from the processor meshes"
+    Info<< "Reconstructing the addressing from the processor meshes"
         << " to the newly reconstructed mesh" << nl << endl;
 
     forAll(databases, procI)
     {
-        Pout<< "Reading processor " << procI << " mesh from "
+        Info<< "Reading processor " << procI << " mesh from "
             << databases[procI].caseName() << endl;
 
         polyMesh procMesh
@@ -548,7 +625,7 @@ int main(int argc, char *argv[])
 
         // From processor point to reconstructed mesh point
 
-        Pout<< "Writing pointProcAddressing to "
+        Info<< "Writing pointProcAddressing to "
             << databases[procI].caseName()
               /procMesh.facesInstance()
               /polyMesh::meshSubDir
@@ -572,7 +649,7 @@ int main(int argc, char *argv[])
 
         // From processor face to reconstructed mesh face
 
-        Pout<< "Writing faceProcAddressing to "
+        Info<< "Writing faceProcAddressing to "
             << databases[procI].caseName()
               /procMesh.facesInstance()
               /polyMesh::meshSubDir
@@ -635,7 +712,7 @@ int main(int argc, char *argv[])
 
         // From processor cell to reconstructed mesh cell
 
-        Pout<< "Writing cellProcAddressing to "
+        Info<< "Writing cellProcAddressing to "
             << databases[procI].caseName()
               /procMesh.facesInstance()
               /polyMesh::meshSubDir
@@ -660,7 +737,7 @@ int main(int argc, char *argv[])
 
         // From processor patch to reconstructed mesh patch
 
-        Pout<< "Writing boundaryProcAddressing to "
+        Info<< "Writing boundaryProcAddressing to "
             << databases[procI].caseName()
               /procMesh.facesInstance()
               /polyMesh::meshSubDir
@@ -681,10 +758,10 @@ int main(int argc, char *argv[])
             boundaryProcAddressing[procI]
         ).write();
 
-        Pout<< endl;
+        Info<< endl;
     }
 
-    Pout<< "End.\n" << endl;
+    Info<< "End.\n" << endl;
 
     return 0;
 }

applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C

@@ -54,9 +54,11 @@ Foam::Field<T> Foam::channelIndex::collapse
 ) const
 {
     // Average and order
+    const Field<T> summedField(regionSum(cellField));
+
     Field<T> regionField
     (
-        regionSum(cellField)
+        summedField
       / regionCount_,
         sortMap_
     );

etc/settings.csh

@@ -112,7 +112,7 @@ case OpenFOAM:
     endif
 
     set compilerBin=$WM_COMPILER_DIR/bin
-    set compilerLib=$WM_COMPILER_DIR/lib$WM_COMPILER_LIB_ARCH:$WM_COMPILER_DIR/lib
+    set compilerLib=$WM_COMPILER_DIR/lib${WM_COMPILER_LIB_ARCH}:$WM_COMPILER_DIR/lib
     breaksw
 endsw
 

src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C

@@ -77,7 +77,7 @@ Foam::Istream& Foam::operator>>(Foam::Istream& is, FixedList<T, Size>& L)
         }
         else
         {
-            // Putback the openning bracket
+            // Putback the opening bracket
             is.putBack(firstToken);
         }
 
@@ -190,8 +190,9 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const FixedList<T, Size>& L)
 
         if (uniform)
         {
-            // Write size of list and start contents delimiter
-            os << token::BEGIN_BLOCK;
+            // Write size of list (so it is valid dictionary entry) and
+            // start contents delimiter
+            os << L.size() << token::BEGIN_BLOCK;
 
             // Write list contents
             os << L[0];
@@ -201,7 +202,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const FixedList<T, Size>& L)
         }
         else if (Size < 11 && contiguous<T>())
         {
-            // Write size of list and start contents delimiter
+            // Write start of contents delimiter
             os << token::BEGIN_LIST;
 
             // Write list contents
@@ -216,7 +217,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const FixedList<T, Size>& L)
         }
         else
         {
-            // Write size of list and start contents delimiter
+            // Write start of contents delimiter
             os << nl << token::BEGIN_LIST;
 
             // Write list contents

src/OpenFOAM/meshes/MeshObject/MeshObject.C

@@ -65,11 +65,11 @@ const Type& Foam::MeshObject<Mesh, Type>::New
 
 
 template<class Mesh, class Type>
-template<class Data>
+template<class Data1>
 const Type& Foam::MeshObject<Mesh, Type>::New
 (
     const Mesh& mesh,
-    const Data& d
+    const Data1& d
 )
 {
     if (!mesh.db().objectRegistry::foundObject<Type>(Type::typeName))
@@ -83,6 +83,69 @@ const Type& Foam::MeshObject<Mesh, Type>::New
 }
 
 
+template<class Mesh, class Type>
+template<class Data1, class Data2>
+const Type& Foam::MeshObject<Mesh, Type>::New
+(
+    const Mesh& mesh,
+    const Data1& d1,
+    const Data2& d2
+)
+{
+    if (!mesh.db().objectRegistry::foundObject<Type>(Type::typeName))
+    {
+        return store(new Type(mesh, d1, d2));
+    }
+    else
+    {
+        return mesh.db().objectRegistry::lookupObject<Type>(Type::typeName);
+    }
+}
+
+
+template<class Mesh, class Type>
+template<class Data1, class Data2, class Data3>
+const Type& Foam::MeshObject<Mesh, Type>::New
+(
+    const Mesh& mesh,
+    const Data1& d1,
+    const Data2& d2,
+    const Data3& d3
+)
+{
+    if (!mesh.db().objectRegistry::foundObject<Type>(Type::typeName))
+    {
+        return store(new Type(mesh, d1, d2, d3));
+    }
+    else
+    {
+        return mesh.db().objectRegistry::lookupObject<Type>(Type::typeName);
+    }
+}
+
+
+template<class Mesh, class Type>
+template<class Data1, class Data2, class Data3, class Data4>
+const Type& Foam::MeshObject<Mesh, Type>::New
+(
+    const Mesh& mesh,
+    const Data1& d1,
+    const Data2& d2,
+    const Data3& d3,
+    const Data4& d4
+)
+{
+    if (!mesh.db().objectRegistry::foundObject<Type>(Type::typeName))
+    {
+        return store(new Type(mesh, d3, d4));
+    }
+    else
+    {
+        return mesh.db().objectRegistry::lookupObject<Type>(Type::typeName);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
 
 template<class Mesh, class Type>

src/OpenFOAM/meshes/MeshObject/MeshObject.H

@@ -69,8 +69,39 @@ public:
 
         static const Type& New(const Mesh& mesh);
 
-        template<class Data>
-        static const Type& New(const Mesh& mesh, const Data& d);
+        template<class Data1>
+        static const Type& New
+        (
+            const Mesh& mesh,
+            const Data1& d
+        );
+
+        template<class Data1, class Data2>
+        static const Type& New
+        (
+            const Mesh& mesh,
+            const Data1&,
+            const Data2&
+        );
+
+        template<class Data1, class Data2, class Data3>
+        static const Type& New
+        (
+            const Mesh& mesh,
+            const Data1&,
+            const Data2&,
+            const Data3&
+        );
+
+        template<class Data1, class Data2, class Data3, class Data4>
+        static const Type& New
+        (
+            const Mesh& mesh,
+            const Data1&,
+            const Data2&,
+            const Data3&,
+            const Data4&
+        );
 
 
     // Destructor

src/OpenFOAM/meshes/meshTools/matchPoints.C

@@ -45,9 +45,19 @@ bool Foam::matchPoints
 
     bool fullMatch = true;
 
-    SortableList<scalar> pts0MagSqr(magSqr(pts0 - origin));
+    point compareOrigin = origin;
 
-    SortableList<scalar> pts1MagSqr(magSqr(pts1 - origin));
+    if (origin == point(VGREAT, VGREAT, VGREAT))
+    {
+        if (pts1.size() > 0)
+        {
+            compareOrigin = sum(pts1)/pts1.size();
+        }
+    }
+
+    SortableList<scalar> pts0MagSqr(magSqr(pts0 - compareOrigin));
+
+    SortableList<scalar> pts1MagSqr(magSqr(pts1 - compareOrigin));
 
     forAll(pts0MagSqr, i)
     {

src/OpenFOAM/meshes/meshTools/mergePoints.C

@@ -40,6 +40,16 @@ bool Foam::mergePoints
     const point& origin
 )
 {
+    point compareOrigin = origin;
+
+    if (origin == point(VGREAT, VGREAT, VGREAT))
+    {
+        if (points.size() > 0)
+        {
+            compareOrigin = sum(points)/points.size();
+        }
+    }
+
     // Create a old to new point mapping array
     pointMap.setSize(points.size());
     pointMap = -1;
@@ -56,7 +66,7 @@ bool Foam::mergePoints
     const scalar mergeTolSqr = sqr(mergeTol);
 
     // Sort points by magSqr
-    SortableList<scalar> sortedMagSqr(magSqr(points - origin));
+    SortableList<scalar> sortedMagSqr(magSqr(points - compareOrigin));
 
     bool hasMerged = false;
 

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C

@@ -271,4 +271,144 @@ Foam::mapDistribute::mapDistribute
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::mapDistribute::compact(const boolList& elemIsUsed)
+{
+    // 1. send back to sender. Have him delete the corresponding element
+    //    from the submap and do the same to the constructMap locally
+    //    (and in same order).
+
+    // Send elemIsUsed field to neighbour. Use nonblocking code from
+    // mapDistribute but in reverse order.
+    {
+        List<boolList> sendFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap_[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                boolList& subField = sendFields[domain];
+                subField.setSize(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = elemIsUsed[map[i]];
+                }
+
+                OPstream::write
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<const char*>(subField.begin()),
+                    subField.size()*sizeof(bool)
+                );
+            }
+        }
+
+        // Set up receives from neighbours
+
+        List<boolList> recvFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap_[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                recvFields[domain].setSize(map.size());
+                IPstream::read
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<char*>(recvFields[domain].begin()),
+                    recvFields[domain].size()*sizeof(bool)
+                );
+            }
+        }
+
+
+        // Set up 'send' to myself - write directly into recvFields
+
+        {
+            const labelList& map = constructMap_[Pstream::myProcNo()];
+
+            recvFields[Pstream::myProcNo()].setSize(map.size());
+            forAll(map, i)
+            {
+                recvFields[Pstream::myProcNo()][i] = elemIsUsed[map[i]];
+            }
+        }
+
+
+        // Wait for all to finish
+
+        OPstream::waitRequests();
+        IPstream::waitRequests();
+
+
+        // Compact out all submap entries that are referring to unused elements
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap_[domain];
+
+            labelList newMap(map.size());
+            label newI = 0;
+
+            forAll(map, i)
+            {
+                if (recvFields[domain][i])
+                {
+                    // So element is used on destination side
+                    newMap[newI++] = map[i];
+                }
+            }
+            if (newI < map.size())
+            {
+                newMap.setSize(newI);
+                subMap_[domain].transfer(newMap);
+            }
+        }
+    }
+
+
+    // 2. remove from construct map - since end-result (element in elemIsUsed)
+    //    not used.
+
+    label maxConstructIndex = -1;
+
+    for (label domain = 0; domain < Pstream::nProcs(); domain++)
+    {
+        const labelList& map = constructMap_[domain];
+
+        labelList newMap(map.size());
+        label newI = 0;
+
+        forAll(map, i)
+        {
+            label destinationI = map[i];
+
+            // Is element is used on destination side
+            if (elemIsUsed[destinationI])
+            {
+                maxConstructIndex = max(maxConstructIndex, destinationI);
+
+                newMap[newI++] = destinationI;
+            }
+        }
+        if (newI < map.size())
+        {
+            newMap.setSize(newI);
+            constructMap_[domain].transfer(newMap);
+        }
+    }
+
+    constructSize_ = maxConstructIndex+1;
+
+    // Clear the schedule (note:not necessary if nothing changed)
+    schedulePtr_.clear();
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H

@@ -51,6 +51,7 @@ SourceFiles
 #include "labelList.H"
 #include "labelPair.H"
 #include "Pstream.H"
+#include "boolList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -155,6 +156,12 @@ public:
 
         // Other
 
+            //- Compact maps. Gets per field a bool whether it is used (locally)
+            //  and works out itself what this side and sender side can remove
+            //  from maps.
+            void compact(const boolList& elemIsUsed);
+
+
             //- Distribute data. Note:schedule only used for Pstream::scheduled
             //  for now, all others just use send-to-all, receive-from-all.
             template<class T>
@@ -168,6 +175,21 @@ public:
                 List<T>&
             );
 
+            //- Distribute data. If multiple processors writing to same
+            //  position adds contributions using cop.
+            template<class T, class CombineOp>
+            static void distribute
+            (
+                const Pstream::commsTypes commsType,
+                const List<labelPair>& schedule,
+                const label constructSize,
+                const labelListList& subMap,
+                const labelListList& constructMap,
+                List<T>&,
+                const CombineOp& cop,
+                const T& nullValue
+            );
+
             //- Distribute data using default commsType.
             template<class T>
             void distribute(List<T>& fld) const

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C

@@ -268,25 +268,34 @@ void Foam::mapDistribute::distribute
         }
 
 
-        // Combine bits. Note that can reuse field storage
+        // Set up 'send' to myself
 
-        // Subset myself
-        const labelList& mySubMap = subMap[Pstream::myProcNo()];
-
-        List<T> subField(mySubMap.size());
-        forAll(mySubMap, i)
         {
-            subField[i] = field[mySubMap[i]];
+            const labelList& map = subMap[Pstream::myProcNo()];
+
+            List<T>& subField = sendFields[Pstream::myProcNo()];
+            subField.setSize(map.size());
+            forAll(map, i)
+            {
+                subField[i] = field[map[i]];
+            }
         }
 
+
+        // Combine bits. Note that can reuse field storage
+
         field.setSize(constructSize);
 
-        // Receive sub field from myself (subField)
-        const labelList& map = constructMap[Pstream::myProcNo()];
 
-        forAll(map, i)
+        // Receive sub field from myself (sendFields[Pstream::myProcNo()])
         {
-            field[map[i]] = subField[i];
+            const labelList& map = constructMap[Pstream::myProcNo()];
+            const List<T>& subField = sendFields[Pstream::myProcNo()];
+
+            forAll(map, i)
+            {
+                field[map[i]] = subField[i];
+            }
         }
 
 
@@ -339,4 +348,326 @@ void Foam::mapDistribute::distribute
 }
 
 
+// Distribute list.
+template<class T, class CombineOp>
+void Foam::mapDistribute::distribute
+(
+    const Pstream::commsTypes commsType,
+    const List<labelPair>& schedule,
+    const label constructSize,
+    const labelListList& subMap,
+    const labelListList& constructMap,
+    List<T>& field,
+    const CombineOp& cop,
+    const T& nullValue
+)
+{
+    if (commsType == Pstream::blocking)
+    {
+        // Since buffered sending can reuse the field to collect the
+        // received data.
+
+        // Send sub field to neighbour
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                List<T> subField(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = field[map[i]];
+                }
+                OPstream toNbr(Pstream::blocking, domain);
+                toNbr << subField;
+            }
+        }
+
+        // Subset myself
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = field[mySubMap[i]];
+        }
+
+        // Receive sub field from myself (subField)
+        const labelList& map = constructMap[Pstream::myProcNo()];
+
+        field.setSize(constructSize);
+        field = nullValue;
+
+        forAll(map, i)
+        {
+            cop(field[map[i]], subField[i]);
+        }
+
+        // Receive sub field from neighbour
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                IPstream fromNbr(Pstream::blocking, domain);
+                List<T> subField(fromNbr);
+
+                if (subField.size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << domain
+                        << " " << map.size() << " but received "
+                        << subField.size() << " elements."
+                        << abort(FatalError);
+                }
+
+                forAll(map, i)
+                {
+                    cop(field[map[i]], subField[i]);
+                }
+            }
+        }
+    }
+    else if (commsType == Pstream::scheduled)
+    {
+        // Need to make sure I don't overwrite field with received data
+        // since the data might need to be sent to another processor. So
+        // allocate a new field for the results.
+        List<T> newField(constructSize, nullValue);
+
+        // Subset myself
+        const labelList& mySubMap = subMap[Pstream::myProcNo()];
+
+        List<T> subField(mySubMap.size());
+        forAll(mySubMap, i)
+        {
+            subField[i] = field[mySubMap[i]];
+        }
+
+        // Receive sub field from myself (subField)
+        const labelList& map = constructMap[Pstream::myProcNo()];
+
+        forAll(map, i)
+        {
+            cop(newField[map[i]], subField[i]);
+        }
+
+        // Schedule will already have pruned 0-sized comms
+        forAll(schedule, i)
+        {
+            const labelPair& twoProcs = schedule[i];
+            label sendProc = twoProcs[0];
+            label recvProc = twoProcs[1];
+
+            if (Pstream::myProcNo() == sendProc)
+            {
+                // I am sender. Send to recvProc.
+                const labelList& map = subMap[recvProc];
+
+                List<T> subField(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = field[map[i]];
+                }
+
+                OPstream toNbr(Pstream::scheduled, recvProc);
+                toNbr << subField;
+            }
+            else
+            {
+                // I am receiver. Receive from sendProc.
+                IPstream fromNbr(Pstream::scheduled, sendProc);
+                List<T> subField(fromNbr);
+
+                const labelList& map = constructMap[sendProc];
+
+                if (subField.size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << sendProc
+                        << " " << map.size() << " but received "
+                        << subField.size() << " elements."
+                        << abort(FatalError);
+                }
+
+                forAll(map, i)
+                {
+                    cop(newField[map[i]], subField[i]);
+                }
+            }
+        }
+        field.transfer(newField);
+    }
+    else if (commsType == Pstream::nonBlocking)
+    {
+        if (!contiguous<T>())
+        {
+            FatalErrorIn
+            (
+                "template<class T>\n"
+                "void mapDistribute::distribute\n"
+                "(\n"
+                "    const Pstream::commsTypes commsType,\n"
+                "    const List<labelPair>& schedule,\n"
+                "    const label constructSize,\n"
+                "    const labelListList& subMap,\n"
+                "    const labelListList& constructMap,\n"
+                "    List<T>& field\n"
+                ")\n"
+            )   << "Non-blocking only supported for contiguous data."
+                << exit(FatalError);
+        }
+
+        // Set up sends to neighbours
+
+        List<List<T > > sendFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = subMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                List<T>& subField = sendFields[domain];
+                subField.setSize(map.size());
+                forAll(map, i)
+                {
+                    subField[i] = field[map[i]];
+                }
+
+                OPstream::write
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<const char*>(subField.begin()),
+                    subField.size()*sizeof(T)
+                );
+            }
+        }
+
+        // Set up receives from neighbours
+
+        List<List<T > > recvFields(Pstream::nProcs());
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                recvFields[domain].setSize(map.size());
+                IPstream::read
+                (
+                    Pstream::nonBlocking,
+                    domain,
+                    reinterpret_cast<char*>(recvFields[domain].begin()),
+                    recvFields[domain].size()*sizeof(T)
+                );
+            }
+        }
+
+        // Set up 'send' to myself
+
+        {
+            const labelList& map = subMap[Pstream::myProcNo()];
+
+            List<T>& subField = sendFields[Pstream::myProcNo()];
+            subField.setSize(map.size());
+            forAll(map, i)
+            {
+                subField[i] = field[map[i]];
+            }
+        }
+
+
+        // Combine bits. Note that can reuse field storage
+
+        field.setSize(constructSize);
+        field = nullValue;
+
+        // Receive sub field from myself (subField)
+        {
+            const labelList& map = constructMap[Pstream::myProcNo()];
+            const List<T>& subField = sendFields[Pstream::myProcNo()];
+
+            forAll(map, i)
+            {
+                cop(field[map[i]], subField[i]);
+            }
+        }
+
+
+        // Wait for all to finish
+
+        OPstream::waitRequests();
+        IPstream::waitRequests();
+        
+        // Collect neighbour fields
+
+        for (label domain = 0; domain < Pstream::nProcs(); domain++)
+        {
+            const labelList& map = constructMap[domain];
+
+            if (domain != Pstream::myProcNo() && map.size() > 0)
+            {
+                if (recvFields[domain].size() != map.size())
+                {
+                    FatalErrorIn
+                    (
+                        "template<class T>\n"
+                        "void mapDistribute::distribute\n"
+                        "(\n"
+                        "    const Pstream::commsTypes commsType,\n"
+                        "    const List<labelPair>& schedule,\n"
+                        "    const label constructSize,\n"
+                        "    const labelListList& subMap,\n"
+                        "    const labelListList& constructMap,\n"
+                        "    List<T>& field\n"
+                        ")\n"
+                    )   << "Expected from processor " << domain
+                        << " " << map.size() << " but received "
+                        << recvFields[domain].size() << " elements."
+                        << abort(FatalError);
+                }
+
+                forAll(map, i)
+                {
+                    cop(field[map[i]], recvFields[domain][i]);
+                }
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn("mapDistribute::distribute(..)")
+            << "Unknown communication schedule " << commsType
+            << abort(FatalError);
+    }
+}
+
+
 // ************************************************************************* //

src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C

@@ -279,27 +279,13 @@ Foam::layerParameters::layerParameters
 {
     const dictionary& layersDict = dict.subDict("layers");
 
-    forAllConstIter(dictionary, layersDict, iter)
+    forAll(boundaryMesh, patchI)
     {
-        const word& key = iter().keyword();
+        const word& patchName = boundaryMesh[patchI].name();
 
-        if (layersDict.isDict(key))
+        if (layersDict.found(patchName))
         {
-            label patchI = boundaryMesh.findPatchID(key);
-
-            if (patchI == -1)
-            {
-                FatalErrorIn
-                (
-                    "layerParameters::layerParameters"
-                    "(const dictionary&, const polyBoundaryMesh&)"
-                )   << "Specified illegal patch " << key
-                    << " in layer dictionary." << endl
-                    << "Valid patch names are " << boundaryMesh.names()
-                    << exit(FatalError);
-            }
-
-            const dictionary& layerDict = layersDict.subDict(key);
+            const dictionary& layerDict = layersDict.subDict(patchName);
 
             numLayers_[patchI] =
                 readLabel(layerDict.lookup("nSurfaceLayers"));

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -1914,27 +1914,28 @@ void Foam::meshRefinement::baffleAndSplitMesh
 
         if (debug)
         {
-            // Dump all these faces to a faceSet.
-            faceSet problemGeom(mesh_, "problemFacesGeom", 100);
-
-            const labelList facePatchGeom
-            (
-                markFacesOnProblemCellsGeometric
-                (
-                    motionDict,
-                    globalToPatch
-                )
-            );
-            forAll(facePatchGeom, faceI)
-            {
-                if (facePatchGeom[faceI] != -1)
-                {
-                    problemGeom.insert(faceI);
-                }
-            }
-            Pout<< "Dumping " << problemGeom.size()
-                << " problem faces to " << problemGeom.objectPath() << endl;
-            problemGeom.write();
+            //- Note: commented out since not properly parallel yet.
+            //// Dump all these faces to a faceSet.
+            //faceSet problemGeom(mesh_, "problemFacesGeom", 100);
+            //
+            //const labelList facePatchGeom
+            //(
+            //    markFacesOnProblemCellsGeometric
+            //    (
+            //        motionDict,
+            //        globalToPatch
+            //    )
+            //);
+            //forAll(facePatchGeom, faceI)
+            //{
+            //    if (facePatchGeom[faceI] != -1)
+            //    {
+            //        problemGeom.insert(faceI);
+            //    }
+            //}
+            //Pout<< "Dumping " << problemGeom.size()
+            //    << " problem faces to " << problemGeom.objectPath() << endl;
+            //problemGeom.write();
 
 
             faceSet problemTopo(mesh_, "problemFacesTopo", 100);

src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C

@@ -539,27 +539,6 @@ void Foam::refinementSurfaces::findHigherIntersection
         return;
     }
 
-
-    // Precalculate per surface whether it has a minlevelfield
-    UPtrList<triSurfaceLabelField> minLevelFields(surfaces_.size());
-    forAll(surfaces_, surfI)
-    {
-        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
-
-        if (isA<triSurfaceMesh>(geom))
-        {
-            const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
-            minLevelFields.set
-            (
-                surfI,
-               &const_cast<triSurfaceLabelField&>
-                (
-                    triMesh.lookupObject<triSurfaceLabelField>("minLevel")
-                )
-            );
-        }
-    }
-
     // Work arrays
     labelList hitMap(identity(start.size()));
     pointField p0(start);
@@ -568,7 +547,21 @@ void Foam::refinementSurfaces::findHigherIntersection
 
     forAll(surfaces_, surfI)
     {
-        allGeometry_[surfaces_[surfI]].findLineAny(p0, p1, hitInfo);
+        const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+        geom.findLineAny(p0, p1, hitInfo);
+
+        labelList minLevelField;
+        if (isA<triSurfaceMesh>(geom))
+        {
+            const triSurfaceMesh& triMesh = refCast<const triSurfaceMesh>(geom);
+            triMesh.getField
+            (
+                "minLevel",
+                hitInfo,
+                minLevelField
+            );
+        }
 
 
         // Copy all hits into arguments, continue with misses
@@ -581,14 +574,9 @@ void Foam::refinementSurfaces::findHigherIntersection
             if (hitInfo[hitI].hit())
             {
                 // Check if minLevelField for this surface.
-                if
-                (
-                    minLevelFields.set(surfI)
-                 && minLevelFields[surfI].size() > 0
-                )
+                if (minLevelField.size() > 0)
                 {
-                    minLocalLevel =
-                        minLevelFields[surfI][hitInfo[hitI].index()];
+                    minLocalLevel = minLevelField[hitI];
                 }
                 else
                 {
@@ -618,7 +606,7 @@ void Foam::refinementSurfaces::findHigherIntersection
         }
 
         // All done? Note that this decision should be synchronised
-        if (newI == 0)
+        if (returnReduce(newI, sumOp<label>()) == 0)
         {
             break;
         }

src/dynamicMesh/fvMeshAdder/fvMeshAdder.H

@@ -40,6 +40,7 @@ SourceFiles
 
 #include "polyMeshAdder.H"
 #include "fvPatchFieldsFwd.H"
+#include "fvsPatchFieldsFwd.H"
 #include "fvPatchFieldMapper.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C

@@ -586,8 +586,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::repatch
                 "fvMeshDistribute::repatch(const labelList&, labelListList&)"
             )   << "reverseFaceMap contains -1 at index:"
                 << index << endl
-                << "This means that the repatch operation was not just a shuffle?"
-                << abort(FatalError);
+                << "This means that the repatch operation was not just"
+                << " a shuffle?" << abort(FatalError);
         }
     }
 
@@ -622,9 +622,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
         )
     );
 
-    bool merged = pointToMaster.size() > 0;
-
-    if (!returnReduce(merged, orOp<bool>()))
+    if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
     {
         return autoPtr<mapPolyMesh>(NULL);
     }
@@ -639,14 +637,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
     // Update fields. No inflation, parallel sync.
     mesh_.updateMesh(map);
 
-    // Move mesh (since morphing does not do this)
-    if (map().hasMotionPoints())
-    {
-        mesh_.movePoints(map().preMotionPoints());
-    }
-
     // Adapt constructMaps for merged points.
-    // 1.4.1: use reversePointMap < -1 feature.
     forAll(constructPointMap, procI)
     {
         labelList& constructMap = constructPointMap[procI];
@@ -655,16 +646,38 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
         {
             label oldPointI = constructMap[i];
 
-            // See if merged into other point
-            Map<label>::const_iterator iter = pointToMaster.find(oldPointI);
+            label newPointI = map().reversePointMap()[oldPointI];
 
-            if (iter != pointToMaster.end())
+            if (newPointI < -1)
             {
-                oldPointI = iter();
+                constructMap[i] = -newPointI-2;
+            }
+            else if (newPointI >= 0)
+            {
+                constructMap[i] = newPointI;
+            }
+            else
+            {
+                FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
+                    << "Problem. oldPointI:" << oldPointI
+                    << " newPointI:" << newPointI << abort(FatalError);
             }
-
-            constructMap[i] = map().reversePointMap()[oldPointI];
         }
+        //- old: use pointToMaster map.
+        //forAll(constructMap, i)
+        //{
+        //    label oldPointI = constructMap[i];
+        //
+        //    // See if merged into other point
+        //    Map<label>::const_iterator iter = pointToMaster.find(oldPointI);
+        //
+        //    if (iter != pointToMaster.end())
+        //    {
+        //        oldPointI = iter();
+        //    }
+        //
+        //    constructMap[i] = map().reversePointMap()[oldPointI];
+        //}
     }
     return map;
 }
@@ -1267,8 +1280,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
             domainMesh.boundaryMesh()
         ).ptr();
     }
-    domainMesh.addFvPatches(patches);
-
+    // Add patches; no parallel comms
+    domainMesh.addFvPatches(patches, false);
 
     // Construct zones
     List<pointZone*> pZonePtrs(pointZoneNames.size());

src/dynamicMesh/polyMeshAdder/polyMeshAdder.C

@@ -1737,85 +1737,218 @@ Foam::autoPtr<Foam::mapAddedPolyMesh> Foam::polyMeshAdder::add
 Foam::Map<Foam::label> Foam::polyMeshAdder::findSharedPoints
 (
     const polyMesh& mesh,
-    const scalar mergeTol
+    const scalar mergeDist
 )
 {
     const labelList& sharedPointLabels = mesh.globalData().sharedPointLabels();
+    const labelList& sharedPointAddr = mesh.globalData().sharedPointAddr();
 
-    labelList sharedToMerged;
-    pointField mergedPoints;
-    bool hasMerged = Foam::mergePoints
-    (
-        pointField
-        (
-            IndirectList<point>
+    // Because of adding the missing pieces e.g. when redistributing a mesh
+    // it can be that there are multiple points on the same processor that
+    // refer to the same shared point.
+
+    // Invert point-to-shared addressing
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    Map<labelList> sharedToMesh(sharedPointLabels.size());
+
+    label nMultiple = 0;
+
+    forAll(sharedPointLabels, i)
+    {
+        label pointI = sharedPointLabels[i];
+
+        label sharedI = sharedPointAddr[i];
+
+        Map<labelList>::iterator iter = sharedToMesh.find(sharedI);
+
+        if (iter != sharedToMesh.end())
+        {
+            // sharedI already used by other point. Add this one.
+
+            nMultiple++;
+
+            labelList& connectedPointLabels = iter();
+
+            label sz = connectedPointLabels.size();
+
+            // Check just to make sure.
+            if (findIndex(connectedPointLabels, pointI) != -1)
+            {
+                FatalErrorIn("polyMeshAdder::findSharedPoints(..)")
+                    << "Duplicate point in sharedPoint addressing." << endl
+                    << "When trying to add point " << pointI << " on shared "
+                    << sharedI  << " with connected points "
+                    << connectedPointLabels
+                    << abort(FatalError);
+            }
+
+            connectedPointLabels.setSize(sz+1);
+            connectedPointLabels[sz] = pointI;
+        }
+        else
+        {
+            sharedToMesh.insert(sharedI, labelList(1, pointI));
+        }
+    }
+
+    
+    // Assign single master for every shared with multiple geometric points
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    Map<label> pointToMaster(nMultiple);
+
+    forAllConstIter(Map<labelList>, sharedToMesh, iter)
+    {
+        const labelList& connectedPointLabels = iter();
+
+        //Pout<< "For shared:" << iter.key()
+        //    << " found points:" << connectedPointLabels
+        //    << " at coords:"
+        //    <<  pointField(mesh.points(), connectedPointLabels) << endl;
+
+        if (connectedPointLabels.size() > 1)
+        {
+            const pointField connectedPoints
             (
                 mesh.points(),
-                sharedPointLabels
-            )()
-        ),
-        mergeTol,
-        false,
-        sharedToMerged,
-        mergedPoints
-    );
+                connectedPointLabels
+            );
 
-    // Find out which sets of points get merged and create a map from
-    // mesh point to unique point.
-
-    Map<label> pointToMaster(10*sharedToMerged.size());
-
-    if (hasMerged)
-    {
-        labelListList mergeSets
-        (
-            invertOneToMany
+            labelList toMergedPoints;
+            pointField mergedPoints;
+            bool hasMerged = Foam::mergePoints
             (
-                sharedToMerged.size(),
-                sharedToMerged
-            )
-        );
+                connectedPoints,
+                mergeDist,
+                false,
+                toMergedPoints,
+                mergedPoints
+            );
 
-        label nMergeSets = 0;
-
-        forAll(mergeSets, setI)
-        {
-            const labelList& mergeSet = mergeSets[setI];
-
-            if (mergeSet.size() > 1)
+            if (hasMerged)
             {
-                // Take as master the shared point with the lowest mesh
-                // point label. (rather arbitrarily - could use max or any other
-                // one of the points)
+                // Invert toMergedPoints
+                const labelListList mergeSets
+                (
+                    invertOneToMany
+                    (
+                        mergedPoints.size(),
+                        toMergedPoints
+                    )
+                );
 
-                nMergeSets++;
-
-                label masterI = labelMax;
-
-                forAll(mergeSet, i)
+                // Find master for valid merges
+                forAll(mergeSets, setI)
                 {
-                    label sharedI = mergeSet[i];
+                    const labelList& mergeSet = mergeSets[setI];
 
-                    masterI = min(masterI, sharedPointLabels[sharedI]);
-                }
+                    if (mergeSet.size() > 1)
+                    {
+                        // Pick lowest numbered point
+                        label masterPointI = labelMax;
 
-                forAll(mergeSet, i)
-                {
-                    label sharedI = mergeSet[i];
+                        forAll(mergeSet, i)
+                        {
+                            label pointI = connectedPointLabels[mergeSet[i]];
 
-                    pointToMaster.insert(sharedPointLabels[sharedI], masterI);
+                            masterPointI = min(masterPointI, pointI);
+                        }
+
+                        forAll(mergeSet, i)
+                        {
+                            label pointI = connectedPointLabels[mergeSet[i]];
+
+                            //Pout<< "Merging point " << pointI
+                            //    << " at " << mesh.points()[pointI]
+                            //    << " into master point "
+                            //    << masterPointI
+                            //    << " at " << mesh.points()[masterPointI]
+                            //    << endl;
+
+                            pointToMaster.insert(pointI, masterPointI);
+                        }
+                    }
                 }
             }
         }
-
-        //if (debug)
-        //{
-        //    Pout<< "polyMeshAdder : merging:"
-        //        << pointToMaster.size() << " into " << nMergeSets << " sets."
-        //        << endl;
-        //}
     }
 
+    //- Old: geometric merging. Causes problems for two close shared points.
+    //labelList sharedToMerged;
+    //pointField mergedPoints;
+    //bool hasMerged = Foam::mergePoints
+    //(
+    //    pointField
+    //    (
+    //        IndirectList<point>
+    //        (
+    //            mesh.points(),
+    //            sharedPointLabels
+    //        )()
+    //    ),
+    //    mergeDist,
+    //    false,
+    //    sharedToMerged,
+    //    mergedPoints
+    //);
+    //
+    //// Find out which sets of points get merged and create a map from
+    //// mesh point to unique point.
+    //
+    //Map<label> pointToMaster(10*sharedToMerged.size());
+    //
+    //if (hasMerged)
+    //{
+    //    labelListList mergeSets
+    //    (
+    //        invertOneToMany
+    //        (
+    //            sharedToMerged.size(),
+    //            sharedToMerged
+    //        )
+    //    );
+    //
+    //    label nMergeSets = 0;
+    //
+    //    forAll(mergeSets, setI)
+    //    {
+    //        const labelList& mergeSet = mergeSets[setI];
+    //
+    //        if (mergeSet.size() > 1)
+    //        {
+    //            // Take as master the shared point with the lowest mesh
+    //            // point label. (rather arbitrarily - could use max or
+    //            // any other one of the points)
+    //
+    //            nMergeSets++;
+    //
+    //            label masterI = labelMax;
+    //
+    //            forAll(mergeSet, i)
+    //            {
+    //                label sharedI = mergeSet[i];
+    //
+    //                masterI = min(masterI, sharedPointLabels[sharedI]);
+    //            }
+    //
+    //            forAll(mergeSet, i)
+    //            {
+    //                label sharedI = mergeSet[i];
+    //
+    //                pointToMaster.insert(sharedPointLabels[sharedI], masterI);
+    //            }
+    //        }
+    //    }
+    //
+    //    //if (debug)
+    //    //{
+    //    //    Pout<< "polyMeshAdder : merging:"
+    //    //        << pointToMaster.size() << " into " << nMergeSets
+    //    //        << " sets." << endl;
+    //    //}
+    //}
+
     return pointToMaster;
 }
 
@@ -1836,8 +1969,7 @@ void Foam::polyMeshAdder::mergePoints
         {
             if (iter() != pointI)
             {
-                //1.4.1: meshMod.removePoint(pointI, iter());
-                meshMod.setAction(polyRemovePoint(pointI));
+                meshMod.removePoint(pointI, iter());
             }
         }
     }

src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C

@@ -4765,69 +4765,70 @@ void Foam::hexRef8::checkRefinementLevels
     }
 
 
-    // Hanging points
-    {
-        // Any patches with points having only two edges.
-
-        boolList isHangingPoint(mesh_.nPoints(), false);
-
-        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            const labelList& meshPoints = pp.meshPoints();
-
-            forAll(meshPoints, i)
-            {
-                label pointI = meshPoints[i];
-
-                const labelList& pEdges = mesh_.pointEdges()[pointI];
-
-                if (pEdges.size() == 2)
-                {
-                    isHangingPoint[pointI] = true;
-                }
-            }
-        }
-
-        syncTools::syncPointList
-        (
-            mesh_,
-            isHangingPoint,
-            andEqOp<bool>(),        // only if all decide it is hanging point
-            true,                   // null
-            false                   // no separation
-        );
-
-        //OFstream str(mesh_.time().path()/"hangingPoints.obj");
-
-        label nHanging = 0;
-
-        forAll(isHangingPoint, pointI)
-        {
-            if (isHangingPoint[pointI])
-            {
-                nHanging++;
-
-                Pout<< "Hanging boundary point " << pointI
-                    << " at " << mesh_.points()[pointI]
-                    << endl;
-                //meshTools::writeOBJ(str, mesh_.points()[pointI]);
-            }
-        }
-
-        if (returnReduce(nHanging, sumOp<label>()) > 0)
-        {
-            FatalErrorIn
-            (
-                "hexRef8::checkRefinementLevels(const label)"
-            )   << "Detected a point used by two edges only (hanging point)"
-                << nl << "This is not allowed"
-                << abort(FatalError);
-        }
-    }
+    //- Gives problems after first splitting off inside mesher.
+    //// Hanging points
+    //{
+    //    // Any patches with points having only two edges.
+    //
+    //    boolList isHangingPoint(mesh_.nPoints(), false);
+    //
+    //    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    //
+    //    forAll(patches, patchI)
+    //    {
+    //        const polyPatch& pp = patches[patchI];
+    //
+    //        const labelList& meshPoints = pp.meshPoints();
+    //
+    //        forAll(meshPoints, i)
+    //        {
+    //            label pointI = meshPoints[i];
+    //
+    //            const labelList& pEdges = mesh_.pointEdges()[pointI];
+    //
+    //            if (pEdges.size() == 2)
+    //            {
+    //                isHangingPoint[pointI] = true;
+    //            }
+    //        }
+    //    }
+    //
+    //    syncTools::syncPointList
+    //    (
+    //        mesh_,
+    //        isHangingPoint,
+    //        andEqOp<bool>(),        // only if all decide it is hanging point
+    //        true,                   // null
+    //        false                   // no separation
+    //    );
+    //
+    //    //OFstream str(mesh_.time().path()/"hangingPoints.obj");
+    //
+    //    label nHanging = 0;
+    //
+    //    forAll(isHangingPoint, pointI)
+    //    {
+    //        if (isHangingPoint[pointI])
+    //        {
+    //            nHanging++;
+    //
+    //            Pout<< "Hanging boundary point " << pointI
+    //                << " at " << mesh_.points()[pointI]
+    //                << endl;
+    //            //meshTools::writeOBJ(str, mesh_.points()[pointI]);
+    //        }
+    //    }
+    //
+    //    if (returnReduce(nHanging, sumOp<label>()) > 0)
+    //    {
+    //        FatalErrorIn
+    //        (
+    //            "hexRef8::checkRefinementLevels(const label)"
+    //        )   << "Detected a point used by two edges only (hanging point)"
+    //            << nl << "This is not allowed"
+    //            << abort(FatalError);
+    //    }
+    //}
 }
 
 

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C

@@ -2617,6 +2617,13 @@ void Foam::polyTopoChange::removePoint
             << abort(FatalError);
     }
 
+    if (pointI == mergePointI)
+    {
+        FatalErrorIn("polyTopoChange::removePoint(const label, const label)")
+            << "Cannot remove/merge point " << pointI << " onto itself."
+            << abort(FatalError);
+    }
+
     points_[pointI] = greatPoint;
     pointMap_[pointI] = -1;
     if (mergePointI >= 0)

src/finiteVolume/Make/files

@@ -40,6 +40,24 @@ $(fvMeshMapper)/fvSurfaceMapper.C
 
 extendedStencil = fvMesh/extendedStencil
 $(extendedStencil)/extendedStencil.C
+$(extendedStencil)/extendedUpwindStencil.C
+$(extendedStencil)/extendedCentredStencil.C
+
+$(extendedStencil)/faceStencil/faceStencil.C
+$(extendedStencil)/faceStencil/faceEdgeCellStencil.C
+$(extendedStencil)/faceStencil/cellFaceCellStencil.C
+$(extendedStencil)/faceStencil/cellPointCellStencil.C
+$(extendedStencil)/faceStencil/cellEdgeCellStencil.C
+
+$(extendedStencil)/extendedStencilMeshObjects/centredCECStencilObject.C
+$(extendedStencil)/extendedStencilMeshObjects/centredCFCStencilObject.C
+$(extendedStencil)/extendedStencilMeshObjects/centredCPCStencilObject.C
+$(extendedStencil)/extendedStencilMeshObjects/centredFECStencilObject.C
+
+$(extendedStencil)/extendedStencilMeshObjects/upwindCECStencilObject.C
+$(extendedStencil)/extendedStencilMeshObjects/upwindCFCStencilObject.C
+$(extendedStencil)/extendedStencilMeshObjects/upwindCPCStencilObject.C
+$(extendedStencil)/extendedStencilMeshObjects/upwindFECStencilObject.C
 
 fvPatchFields = fields/fvPatchFields
 $(fvPatchFields)/fvPatchField/fvPatchFields.C
@@ -168,6 +186,8 @@ $(schemes)/harmonic/harmonic.C
 $(schemes)/localBlended/localBlended.C
 $(schemes)/localMax/localMax.C
 $(schemes)/localMin/localMin.C
+//$(schemes)/linearFit/linearFit.C
+//$(schemes)/linearFit/linearFitData.C
 $(schemes)/quadraticFit/quadraticFit.C
 $(schemes)/quadraticFit/quadraticFitData.C
 
@@ -251,8 +271,8 @@ $(snGradSchemes)/snGradScheme/snGradSchemes.C
 $(snGradSchemes)/correctedSnGrad/correctedSnGrads.C
 $(snGradSchemes)/limitedSnGrad/limitedSnGrads.C
 $(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C
-$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGradData.C
-$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C
+//$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGradData.C
+//$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C
 
 convectionSchemes = finiteVolume/convectionSchemes
 $(convectionSchemes)/convectionScheme/convectionSchemes.C

src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "mapDistribute.H"
+#include "extendedCentredStencil.H"
+#include "faceStencil.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::extendedCentredStencil::extendedCentredStencil(const faceStencil& stencil)
+:
+    extendedStencil(stencil.mesh())
+{
+    stencil_ = stencil;
+
+    // Calculate distribute map (also renumbers stencil)
+    mapPtr_ = calcDistributeMap(stencil.globalNumbering(), stencil_);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// Per face which elements of the stencil to keep.
+void Foam::extendedCentredStencil::compact()
+{
+    boolList isInStencil(map().constructSize(), false);
+
+    forAll(stencil_, faceI)
+    {
+        const labelList& stencilCells = stencil_[faceI];
+
+        forAll(stencilCells, i)
+        {
+            isInStencil[stencilCells[i]] = true;
+        }
+    }
+
+    mapPtr_().compact(isInStencil);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedCentredStencil.H

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::extendedCentredStencil
+
+Description
+
+SourceFiles
+    extendedCentredStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedCentredStencil_H
+#define extendedCentredStencil_H
+
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class faceStencil;
+
+/*---------------------------------------------------------------------------*\
+                           Class extendedCentredStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedCentredStencil
+:
+    public extendedStencil
+{
+    // Private data
+
+        //- Swap map for getting neigbouring data
+        autoPtr<mapDistribute> mapPtr_;
+
+        //- Per face the stencil.
+        labelListList stencil_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        extendedCentredStencil(const extendedCentredStencil&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const extendedCentredStencil&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from uncompacted face stencil
+        explicit extendedCentredStencil(const faceStencil&);
+
+
+    // Member Functions
+
+        //- Return reference to the parallel distribution map
+        const mapDistribute& map() const
+        {
+            return mapPtr_();
+        }
+
+        //- Return reference to the stencil
+        const labelListList& stencil() const
+        {
+            return stencil_;
+        }
+
+        //- After removing elements from the stencil adapt the schedule (map).
+        void compact();
+
+        //- Use map to get the data into stencil order
+        template<class T>
+        void collectData
+        (
+            const GeometricField<T, fvPatchField, volMesh>& fld,
+            List<List<T> >& stencilFld
+        ) const
+        {
+            extendedStencil::collectData(map(), stencil(), fld, stencilFld);
+        }
+
+        //- Sum vol field contributions to create face values
+        template<class Type>
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > weightedSum
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& fld,
+            const List<List<scalar> >& stencilWeights
+        ) const
+        {
+            return extendedStencil::weightedSum
+            (
+                map(),
+                stencil(),
+                fld,
+                stencilWeights
+            );
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencil.C

@@ -32,154 +32,83 @@ License
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Calculates per face a list of global cell/face indices.
-void Foam::extendedStencil::calcFaceStencils
+void Foam::extendedStencil::calcFaceStencil
 (
-    const polyMesh& mesh,
-    const globalIndex& globalNumbering
+    const labelListList& globalCellCells,
+    labelListList& faceStencil
 )
 {
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    const label nBnd = mesh.nFaces()-mesh.nInternalFaces();
-    const labelList& own = mesh.faceOwner();
-    const labelList& nei = mesh.faceNeighbour();
-
-
-    // Determine neighbouring global cell or boundary face
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelList neiGlobal(nBnd);
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-        label faceI = pp.start();
-
-        if (pp.coupled())
-        {
-            // For coupled faces get the cell on the other side
-            forAll(pp, i)
-            {
-                label bFaceI = faceI-mesh.nInternalFaces();
-                neiGlobal[bFaceI] = globalNumbering.toGlobal(own[faceI]);
-                faceI++;
-            }
-        }
-        else if (isA<emptyPolyPatch>(pp))
-        {
-            forAll(pp, i)
-            {
-                label bFaceI = faceI-mesh.nInternalFaces();
-                neiGlobal[bFaceI] = -1;
-                faceI++;
-            }
-        }
-        else
-        {
-            // For noncoupled faces get the boundary face.
-            forAll(pp, i)
-            {
-                label bFaceI = faceI-mesh.nInternalFaces();
-                neiGlobal[bFaceI] =
-                    globalNumbering.toGlobal(mesh.nCells()+bFaceI);
-                faceI++;
-            }
-        }
-    }
-    syncTools::swapBoundaryFaceList(mesh, neiGlobal, false);
-
-
-    // Determine cellCells in global numbering
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    labelListList globalCellCells(mesh.nCells());
-    forAll(globalCellCells, cellI)
-    {
-        const cell& cFaces = mesh.cells()[cellI];
-
-        labelList& cCells = globalCellCells[cellI];
-
-        cCells.setSize(cFaces.size());
-
-        // Collect neighbouring cells/faces
-        label nNbr = 0;
-        forAll(cFaces, i)
-        {
-            label faceI = cFaces[i];
-
-            if (mesh.isInternalFace(faceI))
-            {
-                label nbrCellI = own[faceI];
-                if (nbrCellI == cellI)
-                {
-                    nbrCellI = nei[faceI];
-                }
-                cCells[nNbr++] = globalNumbering.toGlobal(nbrCellI);
-            }
-            else
-            {
-                label nbrCellI = neiGlobal[faceI-mesh.nInternalFaces()];
-                if (nbrCellI != -1)
-                {
-                    cCells[nNbr++] = nbrCellI;
-                }
-            }
-        }
-        cCells.setSize(nNbr);
-    }
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
 
 
     // Determine neighbouring global cell Cells
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     labelListList neiGlobalCellCells(nBnd);
-    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
+    forAll(patches, patchI)
     {
-        neiGlobalCellCells[faceI-mesh.nInternalFaces()] =
-            globalCellCells[own[faceI]];
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                neiGlobalCellCells[faceI-mesh_.nInternalFaces()] =
+                    globalCellCells[own[faceI]];
+                faceI++;
+            }
+        }
     }
-    syncTools::swapBoundaryFaceList(mesh, neiGlobalCellCells, false);
+    syncTools::swapBoundaryFaceList(mesh_, neiGlobalCellCells, false);
 
 
 
     // Construct stencil in global numbering
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    stencil_.setSize(mesh.nFaces());
+    faceStencil.setSize(mesh_.nFaces());
 
-    labelHashSet faceStencil;
+    labelHashSet faceStencilSet;
 
-    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
     {
-        faceStencil.clear();
-        label globalOwn = globalNumbering.toGlobal(own[faceI]);
-        faceStencil.insert(globalOwn);
+        faceStencilSet.clear();
+
         const labelList& ownCCells = globalCellCells[own[faceI]];
+        label globalOwn = ownCCells[0];
+        // Insert cellCells
         forAll(ownCCells, i)
         {
-            faceStencil.insert(ownCCells[i]);
+            faceStencilSet.insert(ownCCells[i]);
         }
 
-        label globalNei = globalNumbering.toGlobal(nei[faceI]);
-        faceStencil.insert(globalNei);
         const labelList& neiCCells = globalCellCells[nei[faceI]];
+        label globalNei = neiCCells[0];
+        // Insert cellCells
         forAll(neiCCells, i)
         {
-            faceStencil.insert(neiCCells[i]);
+            faceStencilSet.insert(neiCCells[i]);
         }
 
         // Guarantee owner first, neighbour second.
-        stencil_[faceI].setSize(faceStencil.size());
+        faceStencil[faceI].setSize(faceStencilSet.size());
         label n = 0;
-        stencil_[faceI][n++] = globalOwn;
-        stencil_[faceI][n++] = globalNei;
-        forAllConstIter(labelHashSet, faceStencil, iter)
+        faceStencil[faceI][n++] = globalOwn;
+        faceStencil[faceI][n++] = globalNei;
+        forAllConstIter(labelHashSet, faceStencilSet, iter)
         {
             if (iter.key() != globalOwn && iter.key() != globalNei)
             {
-                stencil_[faceI][n++] = iter.key();
+                faceStencil[faceI][n++] = iter.key();
             }
         }
-        //Pout<< "internalface:" << faceI << " toc:" << faceStencil.toc()
-        //    << " stencil:" << stencil_[faceI] << endl;
+        //Pout<< "internalface:" << faceI << " toc:" << faceStencilSet.toc()
+        //    << " faceStencil:" << faceStencil[faceI] << endl;
     }
     forAll(patches, patchI)
     {
@@ -190,41 +119,40 @@ void Foam::extendedStencil::calcFaceStencils
         {
             forAll(pp, i)
             {
-                faceStencil.clear();
-                label globalOwn = globalNumbering.toGlobal(own[faceI]);
-                faceStencil.insert(globalOwn);
+                faceStencilSet.clear();
+
                 const labelList& ownCCells = globalCellCells[own[faceI]];
+                label globalOwn = ownCCells[0];
                 forAll(ownCCells, i)
                 {
-                    faceStencil.insert(ownCCells[i]);
+                    faceStencilSet.insert(ownCCells[i]);
                 }
-                // Get the coupled cell
-                label globalNei = neiGlobal[faceI-mesh.nInternalFaces()];
-                faceStencil.insert(globalNei);
+
                 // And the neighbours of the coupled cell
                 const labelList& neiCCells =
-                    neiGlobalCellCells[faceI-mesh.nInternalFaces()];
+                    neiGlobalCellCells[faceI-mesh_.nInternalFaces()];
+                label globalNei = neiCCells[0];
                 forAll(neiCCells, i)
                 {
-                    faceStencil.insert(neiCCells[i]);
+                    faceStencilSet.insert(neiCCells[i]);
                 }
 
                 // Guarantee owner first, neighbour second.
-                stencil_[faceI].setSize(faceStencil.size());
+                faceStencil[faceI].setSize(faceStencilSet.size());
                 label n = 0;
-                stencil_[faceI][n++] = globalOwn;
-                stencil_[faceI][n++] = globalNei;
-                forAllConstIter(labelHashSet, faceStencil, iter)
+                faceStencil[faceI][n++] = globalOwn;
+                faceStencil[faceI][n++] = globalNei;
+                forAllConstIter(labelHashSet, faceStencilSet, iter)
                 {
                     if (iter.key() != globalOwn && iter.key() != globalNei)
                     {
-                        stencil_[faceI][n++] = iter.key();
+                        faceStencil[faceI][n++] = iter.key();
                     }
                 }
 
                 //Pout<< "coupledface:" << faceI
-                //    << " toc:" << faceStencil.toc()
-                //    << " stencil:" << stencil_[faceI] << endl;
+                //    << " toc:" << faceStencilSet.toc()
+                //    << " faceStencil:" << faceStencil[faceI] << endl;
 
                 faceI++;
             }
@@ -233,31 +161,30 @@ void Foam::extendedStencil::calcFaceStencils
         {
             forAll(pp, i)
             {
-                faceStencil.clear();
-                label globalOwn = globalNumbering.toGlobal(own[faceI]);
-                faceStencil.insert(globalOwn);
+                faceStencilSet.clear();
+
                 const labelList& ownCCells = globalCellCells[own[faceI]];
+                label globalOwn = ownCCells[0];
                 forAll(ownCCells, i)
                 {
-                    faceStencil.insert(ownCCells[i]);
+                    faceStencilSet.insert(ownCCells[i]);
                 }
 
-
                 // Guarantee owner first, neighbour second.
-                stencil_[faceI].setSize(faceStencil.size());
+                faceStencil[faceI].setSize(faceStencilSet.size());
                 label n = 0;
-                stencil_[faceI][n++] = globalOwn;
-                forAllConstIter(labelHashSet, faceStencil, iter)
+                faceStencil[faceI][n++] = globalOwn;
+                forAllConstIter(labelHashSet, faceStencilSet, iter)
                 {
                     if (iter.key() != globalOwn)
                     {
-                        stencil_[faceI][n++] = iter.key();
+                        faceStencil[faceI][n++] = iter.key();
                     }
                 }
 
                 //Pout<< "boundaryface:" << faceI
-                //    << " toc:" << faceStencil.toc()
-                //    << " stencil:" << stencil_[faceI] << endl;
+                //    << " toc:" << faceStencilSet.toc()
+                //    << " faceStencil:" << faceStencil[faceI] << endl;
 
                 faceI++;
             }
@@ -266,34 +193,13 @@ void Foam::extendedStencil::calcFaceStencils
 }
 
 
-// Calculates extended stencil. This is per face
-// - owner
-// - cellCells of owner
-// - neighbour
-// - cellCells of neighbour
-// It comes in two parts:
-// - a map which collects/distributes all necessary data in a compact array
-// - the stencil (a labelList per face) which is a set of indices into this
-//   compact array.
-// The compact array is laid out as follows:
-// - first data for current processor (Pstream::myProcNo())
-//      - all cells
-//      - all boundary faces
-// - then per processor
-//      - all used cells and boundary faces
-void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
+Foam::autoPtr<Foam::mapDistribute> Foam::extendedStencil::calcDistributeMap
+(
+    const globalIndex& globalNumbering,
+    labelListList& faceStencil
+)
 {
-    const label nBnd = mesh.nFaces()-mesh.nInternalFaces();
-
-    // Global numbering for cells and boundary faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    globalIndex globalNumbering(mesh.nCells()+nBnd);
-
-
-    // Calculate stencil in global cell indices
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    calcFaceStencils(mesh, globalNumbering);
+    const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
 
 
     // Convert stencil to schedule
@@ -309,8 +215,8 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
     //    these are always all needed.
     List<Map<label> > globalToProc(Pstream::nProcs());
     {
-        const labelList& procPatchMap = mesh.globalData().procPatchMap();
-        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+        const labelList& procPatchMap = mesh_.globalData().procPatchMap();
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
         // Presize with (as estimate) size of patch to neighbour.
         forAll(procPatchMap, procI)
@@ -325,9 +231,9 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
         }
 
         // Collect all (non-local) globalcells/faces needed.
-        forAll(stencil_, faceI)
+        forAll(faceStencil, faceI)
         {
-            const labelList& stencilCells = stencil_[faceI];
+            const labelList& stencilCells = faceStencil[faceI];
 
             forAll(stencilCells, i)
             {
@@ -359,17 +265,17 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
         }
 
 
-        // forAll(globalToProc, procI)
-        // {
-        //     Pout<< "From processor:" << procI << " want cells/faces:" << endl;
-        //     forAllConstIter(Map<label>, globalToProc[procI], iter)
-        //     {
-        //         Pout<< "    global:" << iter.key()
-        //             << " local:" << globalNumbering.toLocal(procI, iter.key())
-        //             << endl;
-        //     }
-        //     Pout<< endl;
-        // }
+        //forAll(globalToProc, procI)
+        //{
+        //    Pout<< "From processor:" << procI << " want cells/faces:" << endl;
+        //    forAllConstIter(Map<label>, globalToProc[procI], iter)
+        //    {
+        //        Pout<< "    global:" << iter.key()
+        //            << " local:" << globalNumbering.toLocal(procI, iter.key())
+        //            << endl;
+        //    }
+        //    Pout<< endl;
+        //}
     }
 
 
@@ -379,19 +285,15 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
 
     labelList compactStart(Pstream::nProcs());
     compactStart[Pstream::myProcNo()] = 0;
-    label nCompact = mesh.nCells()+nBnd;
+    label nCompact = mesh_.nCells()+nBnd;
     forAll(compactStart, procI)
     {
         if (procI != Pstream::myProcNo())
         {
             compactStart[procI] = nCompact;
             nCompact += globalToProc[procI].size();
-
-            // Pout<< "Data wanted from " << procI << " starts at "
-            //     << compactStart[procI] << endl;
         }
     }
-    // Pout<< "Overall cells needed:" << nCompact << endl;
 
 
     // 3. Find out what to receive/send in compact addressing.
@@ -411,11 +313,6 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
                 i++;
             }
 
-            // Pout<< "From proc:" << procI
-            //     << " I need (globalcells):" << wantedGlobals
-            //     << " which are my compact:" << recvCompact[procI]
-            //     << endl;
-
             // Send the global cell numbers I need from procI
             OPstream str(Pstream::blocking, procI);
             str << wantedGlobals;
@@ -424,7 +321,7 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
         {
             recvCompact[procI] =
                 compactStart[procI]
-              + identity(mesh.nCells()+nBnd);
+              + identity(mesh_.nCells()+nBnd);
         }
     }
     labelListList sendCompact(Pstream::nProcs());
@@ -455,9 +352,9 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
     }
 
     // Convert stencil to compact numbering
-    forAll(stencil_, faceI)
+    forAll(faceStencil, faceI)
     {
-        labelList& stencilCells = stencil_[faceI];
+        labelList& stencilCells = faceStencil[faceI];
 
         forAll(stencilCells, i)
         {
@@ -476,10 +373,9 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
 
         }
     }
-    // Pout<< "***stencil_:" << stencil_ << endl;
 
     // Constuct map for distribution of compact data.
-    mapPtr_.reset
+    return autoPtr<mapDistribute>
     (
         new mapDistribute
         (
@@ -494,32 +390,10 @@ void Foam::extendedStencil::calcExtendedFaceStencil(const polyMesh& mesh)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::extendedStencil::extendedStencil
-(
-    const mapDistribute& map,
-    const labelListList& stencil
-)
+Foam::extendedStencil::extendedStencil(const polyMesh& mesh)
 :
-    mapPtr_
-    (
-        autoPtr<mapDistribute>
-        (
-            new mapDistribute
-            (
-                map.constructSize(),
-                map.subMap(),
-                map.constructMap()
-            )
-        )
-    ),
-    stencil_(stencil)
+    mesh_(mesh)
 {}
 
 
-Foam::extendedStencil::extendedStencil(const polyMesh& mesh)
-{
-    calcExtendedFaceStencil(mesh);
-}
-
-
 // ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencil.H

@@ -42,6 +42,7 @@ Description
 
 SourceFiles
     extendedStencil.C
+    extendedStencilTemplates.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -65,23 +66,34 @@ class globalIndex;
 
 class extendedStencil
 {
-    // Private data
+protected:
 
-        //- Swap map for getting neigbouring data
-        autoPtr<mapDistribute> mapPtr_;
+    // Protected data
 
-        //- Per face the stencil.
-        labelListList stencil_;
+        const polyMesh& mesh_;
 
 
+    // Protected Member Functions
+
+        //- Collect cell neighbours into extended stencil
+        void calcFaceStencil
+        (
+            const labelListList& globalCellCells,
+            labelListList& faceStencil
+        );
+
+        //- Calculate distribute map
+        autoPtr<mapDistribute> calcDistributeMap
+        (
+            const globalIndex& globalNumbering,
+            labelListList& faceStencil
+        );
+
+
+private:
+
     // Private Member Functions
 
-        void calcFaceStencils(const polyMesh&, const globalIndex&);
-
-        //- Calculate the stencil (but not weights)
-        void calcExtendedFaceStencil(const polyMesh&);
-
-
         //- Disallow default bitwise copy construct
         extendedStencil(const extendedStencil&);
 
@@ -93,43 +105,32 @@ public:
 
     // Constructors
 
-        //- Construct from components
-        extendedStencil(const mapDistribute& map, const labelListList&);
-
-        //- Construct from all cells and boundary faces
-        extendedStencil(const polyMesh&);
-
+        //- Construct from mesh
+        explicit extendedStencil(const polyMesh&);
 
 
     // Member Functions
 
-        //- Return reference to the parallel distribution map
-        const mapDistribute& map() const
-        {
-            return mapPtr_();
-        }
-
-        //- Return reference to the stencil
-        const labelListList& stencil() const
-        {
-            return stencil_;
-        }
-
         //- Use map to get the data into stencil order
         template<class T>
-        void collectData
+        static void collectData
         (
+            const mapDistribute& map,
+            const labelListList& stencil,
             const GeometricField<T, fvPatchField, volMesh>& fld,
             List<List<T> >& stencilFld
-        ) const;
+        );
 
-        //- Calculate weighted sum of vol field
+        //- Sum vol field contributions to create face values
         template<class Type>
-        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > weightedSum
+        static tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        weightedSum
         (
+            const mapDistribute& map,
+            const labelListList& stencil,
             const GeometricField<Type, fvPatchField, volMesh>& fld,
             const List<List<scalar> >& stencilWeights
-        ) const;
+        );
 
 };
 

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCECStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "centredCECStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(centredCECStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCECStencilObject.H

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::centredCECStencilObject
+
+Description
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef centredCECStencilObject_H
+#define centredCECStencilObject_H
+
+#include "extendedCentredStencil.H"
+#include "cellEdgeCellStencil.H"
+#include "MeshObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class centredCECStencilObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class centredCECStencilObject
+:
+    public MeshObject<fvMesh, centredCECStencilObject>,
+    public extendedCentredStencil
+{
+
+public:
+
+    TypeName("centredCECStencil");
+
+    // Constructors
+
+        //- Construct from uncompacted face stencil
+        explicit centredCECStencilObject
+        (
+            const fvMesh& mesh
+        )
+        :
+            MeshObject<fvMesh, centredCECStencilObject>(mesh),
+            extendedCentredStencil(cellEdgeCellStencil(mesh))
+        {}
+
+
+    // Destructor
+
+        virtual ~centredCECStencilObject()
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "centredCFCStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(centredCFCStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCFCStencilObject.H

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::centredCFCStencilObject
+
+Description
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef centredCFCStencilObject_H
+#define centredCFCStencilObject_H
+
+#include "extendedCentredStencil.H"
+#include "cellFaceCellStencil.H"
+#include "MeshObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class centredCFCStencilObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class centredCFCStencilObject
+:
+    public MeshObject<fvMesh, centredCFCStencilObject>,
+    public extendedCentredStencil
+{
+
+public:
+
+    TypeName("centredCFCStencil");
+
+    // Constructors
+
+        //- Construct from uncompacted face stencil
+        explicit centredCFCStencilObject
+        (
+            const fvMesh& mesh
+        )
+        :
+            MeshObject<fvMesh, centredCFCStencilObject>(mesh),
+            extendedCentredStencil(cellFaceCellStencil(mesh))
+        {}
+
+
+    // Destructor
+
+        virtual ~centredCFCStencilObject()
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCPCStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "centredCPCStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(centredCPCStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredCPCStencilObject.H

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::centredCPCStencilObject
+
+Description
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef centredCPCStencilObject_H
+#define centredCPCStencilObject_H
+
+#include "extendedCentredStencil.H"
+#include "cellPointCellStencil.H"
+#include "MeshObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class centredCPCStencilObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class centredCPCStencilObject
+:
+    public MeshObject<fvMesh, centredCPCStencilObject>,
+    public extendedCentredStencil
+{
+
+public:
+
+    TypeName("centredCFCStencil");
+
+    // Constructors
+
+        //- Construct from uncompacted face stencil
+        explicit centredCPCStencilObject
+        (
+            const fvMesh& mesh
+        )
+        :
+            MeshObject<fvMesh, centredCPCStencilObject>(mesh),
+            extendedCentredStencil(cellPointCellStencil(mesh))
+        {}
+
+
+    // Destructor
+
+        virtual ~centredCPCStencilObject()
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredFECStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "centredFECStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(centredFECStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/centredFECStencilObject.H

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::centredFECStencilObject
+
+Description
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef centredFECStencilObject_H
+#define centredFECStencilObject_H
+
+#include "extendedCentredStencil.H"
+#include "faceEdgeCellStencil.H"
+#include "MeshObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class centredFECStencilObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class centredFECStencilObject
+:
+    public MeshObject<fvMesh, centredFECStencilObject>,
+    public extendedCentredStencil
+{
+
+public:
+
+    TypeName("centredCFCStencil");
+
+    // Constructors
+
+        //- Construct from uncompacted face stencil
+        explicit centredFECStencilObject
+        (
+            const fvMesh& mesh
+        )
+        :
+            MeshObject<fvMesh, centredFECStencilObject>(mesh),
+            extendedCentredStencil(faceEdgeCellStencil(mesh))
+        {}
+
+
+    // Destructor
+
+        virtual ~centredFECStencilObject()
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCECStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "upwindCECStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(upwindCECStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCECStencilObject.H

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::upwindCECStencilObject
+
+Description
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef upwindCECStencilObject_H
+#define upwindCECStencilObject_H
+
+#include "extendedUpwindStencil.H"
+#include "cellEdgeCellStencil.H"
+#include "MeshObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class upwindCECStencilObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class upwindCECStencilObject
+:
+    public MeshObject<fvMesh, upwindCECStencilObject>,
+    public extendedUpwindStencil
+{
+
+public:
+
+    TypeName("upwindCFCStencil");
+
+    // Constructors
+
+        //- Construct from uncompacted face stencil
+        explicit upwindCECStencilObject
+        (
+            const fvMesh& mesh,
+            const scalar minOpposedness
+        )
+        :
+            MeshObject<fvMesh, upwindCECStencilObject>(mesh),
+            extendedUpwindStencil(cellEdgeCellStencil(mesh), minOpposedness)
+        {}
+
+
+    // Destructor
+
+        virtual ~upwindCECStencilObject()
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCFCStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "upwindCFCStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(upwindCFCStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCFCStencilObject.H

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::upwindCFCStencilObject
+
+Description
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef upwindCFCStencilObject_H
+#define upwindCFCStencilObject_H
+
+#include "extendedUpwindStencil.H"
+#include "cellFaceCellStencil.H"
+#include "MeshObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class upwindCFCStencilObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class upwindCFCStencilObject
+:
+    public MeshObject<fvMesh, upwindCFCStencilObject>,
+    public extendedUpwindStencil
+{
+
+public:
+
+    TypeName("upwindCFCStencil");
+
+    // Constructors
+
+        //- Construct from uncompacted face stencil
+        explicit upwindCFCStencilObject
+        (
+            const fvMesh& mesh,
+            const scalar minOpposedness
+        )
+        :
+            MeshObject<fvMesh, upwindCFCStencilObject>(mesh),
+            extendedUpwindStencil(cellFaceCellStencil(mesh), minOpposedness)
+        {}
+
+
+    // Destructor
+
+        virtual ~upwindCFCStencilObject()
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCPCStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "upwindCPCStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(upwindCPCStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindCPCStencilObject.H

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::upwindCPCStencilObject
+
+Description
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef upwindCPCStencilObject_H
+#define upwindCPCStencilObject_H
+
+#include "extendedUpwindStencil.H"
+#include "cellPointCellStencil.H"
+#include "MeshObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class upwindCPCStencilObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class upwindCPCStencilObject
+:
+    public MeshObject<fvMesh, upwindCPCStencilObject>,
+    public extendedUpwindStencil
+{
+
+public:
+
+    TypeName("upwindCFCStencil");
+
+    // Constructors
+
+        //- Construct from uncompacted face stencil
+        explicit upwindCPCStencilObject
+        (
+            const fvMesh& mesh,
+            const scalar minOpposedness
+        )
+        :
+            MeshObject<fvMesh, upwindCPCStencilObject>(mesh),
+            extendedUpwindStencil(cellPointCellStencil(mesh), minOpposedness)
+        {}
+
+
+    // Destructor
+
+        virtual ~upwindCPCStencilObject()
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindFECStencilObject.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "upwindFECStencilObject.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(upwindFECStencilObject, 0);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilMeshObjects/upwindFECStencilObject.H

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::upwindFECStencilObject
+
+Description
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef upwindFECStencilObject_H
+#define upwindFECStencilObject_H
+
+#include "extendedUpwindStencil.H"
+#include "faceEdgeCellStencil.H"
+#include "MeshObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class upwindFECStencilObject Declaration
+\*---------------------------------------------------------------------------*/
+
+class upwindFECStencilObject
+:
+    public MeshObject<fvMesh, upwindFECStencilObject>,
+    public extendedUpwindStencil
+{
+
+public:
+
+    TypeName("upwindCFCStencil");
+
+    // Constructors
+
+        //- Construct from uncompacted face stencil
+        explicit upwindFECStencilObject
+        (
+            const fvMesh& mesh,
+            const scalar minOpposedness
+        )
+        :
+            MeshObject<fvMesh, upwindFECStencilObject>(mesh),
+            extendedUpwindStencil(faceEdgeCellStencil(mesh), minOpposedness)
+        {}
+
+
+    // Destructor
+
+        virtual ~upwindFECStencilObject()
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedStencilTemplates.C

@@ -31,12 +31,14 @@ License
 template<class Type>
 void Foam::extendedStencil::collectData
 (
+    const mapDistribute& map,
+    const labelListList& stencil,
     const GeometricField<Type, fvPatchField, volMesh>& fld,
     List<List<Type> >& stencilFld
-) const
+)
 {
     // 1. Construct cell data in compact addressing
-    List<Type> compactFld(map().constructSize(), pTraits<Type>::zero);
+    List<Type> compactFld(map.constructSize(), pTraits<Type>::zero);
 
     // Insert my internal values
     forAll(fld, cellI)
@@ -56,14 +58,14 @@ void Foam::extendedStencil::collectData
     }
 
     // Do all swapping
-    map().distribute(compactFld);
+    map.distribute(compactFld);
 
     // 2. Pull to stencil
-    stencilFld.setSize(stencil_.size());
+    stencilFld.setSize(stencil.size());
 
-    forAll(stencil_, faceI)
+    forAll(stencil, faceI)
     {
-        const labelList& compactCells = stencil_[faceI];
+        const labelList& compactCells = stencil[faceI];
 
         stencilFld[faceI].setSize(compactCells.size());
 
@@ -79,15 +81,17 @@ template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
 Foam::extendedStencil::weightedSum
 (
+    const mapDistribute& map,
+    const labelListList& stencil,
     const GeometricField<Type, fvPatchField, volMesh>& fld,
     const List<List<scalar> >& stencilWeights
-) const
+)
 {
     const fvMesh& mesh = fld.mesh();
 
     // Collect internal and boundary values
     List<List<Type> > stencilFld;
-    collectData(fld, stencilFld);
+    collectData(map, stencil, fld, stencilFld);
 
     tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
     (
@@ -110,6 +114,7 @@ Foam::extendedStencil::weightedSum
     );
     GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsfCorr();
 
+    // Internal faces
     for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
     {
         const List<Type>& stField = stencilFld[faceI];
@@ -121,33 +126,30 @@ Foam::extendedStencil::weightedSum
         }
     }
 
-    // Coupled boundaries
-    /*
+    // Boundaries. Either constrained or calculated so assign value
+    // directly (instead of nicely using operator==)
     typename GeometricField<Type, fvsPatchField, surfaceMesh>::
         GeometricBoundaryField& bSfCorr = sf.boundaryField();
+
     forAll(bSfCorr, patchi)
     {
         fvsPatchField<Type>& pSfCorr = bSfCorr[patchi];
 
-        if (pSfCorr.coupled())
+        label faceI = pSfCorr.patch().patch().start();
+
+        forAll(pSfCorr, i)
         {
-            label faceI = pSfCorr.patch().patch().start();
+            const List<Type>& stField = stencilFld[faceI];
+            const List<scalar>& stWeight = stencilWeights[faceI];
 
-            forAll(pSfCorr, i)
+            forAll(stField, j)
             {
-                const List<Type>& stField = stencilFld[faceI];
-                const List<scalar>& stWeight = stencilWeights[faceI];
-
-                forAll(stField, j)
-                {
-                    pSfCorr[i] += stField[j]*stWeight[j];
-                }
-
-                faceI++;
+                pSfCorr[i] += stField[j]*stWeight[j];
             }
+
+            faceI++;
         }
     }
-    */
 
     return tsfCorr;
 }

src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencil.C

@@ -0,0 +1,429 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "extendedUpwindStencil.H"
+#include "faceStencil.H"
+#include "syncTools.H"
+#include "SortableList.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::extendedUpwindStencil::selectOppositeFaces
+(
+    const boolList& nonEmptyFace,
+    const scalar minOpposedness,
+    const label faceI,
+    const label cellI,
+    DynamicList<label>& oppositeFaces
+) const
+{
+    const vectorField& areas = mesh_.faceAreas();
+    const labelList& own = mesh_.faceOwner();
+    const cell& cFaces = mesh_.cells()[cellI];
+
+    SortableList<scalar> opposedness(cFaces.size(), -GREAT);
+
+    // Pick up all the faces that oppose this one.
+    forAll(cFaces, i)
+    {
+        label otherFaceI = cFaces[i];
+
+        if (otherFaceI != faceI && nonEmptyFace[otherFaceI])
+        {
+            if ((own[otherFaceI] == cellI) == (own[faceI] == cellI))
+            {
+                opposedness[i] = -(areas[otherFaceI] & areas[faceI]);
+            }
+            else
+            {
+                opposedness[i] = (areas[otherFaceI] & areas[faceI]);
+            }
+        }
+    }
+
+    label sz = opposedness.size();
+
+    oppositeFaces.clear();
+
+    scalar myAreaSqr = magSqr(areas[faceI]);
+
+    if (myAreaSqr > VSMALL)
+    {
+        forAll(opposedness, i)
+        {
+            opposedness[i] /= myAreaSqr;
+        }
+        // Sort in incrementing order
+        opposedness.sort();
+
+        // Pick largest no matter what
+        oppositeFaces.append(cFaces[opposedness.indices()[sz-1]]);
+
+        for (label i = sz-2; i >= 0; --i)
+        {
+            if (opposedness[i] < minOpposedness)
+            {
+                break;
+            }
+            oppositeFaces.append(cFaces[opposedness.indices()[i]]);
+        }
+    }
+    else
+    {
+        // Sort in incrementing order
+        opposedness.sort();
+
+        // Tiny face. Do what?
+        // Pick largest no matter what
+        oppositeFaces.append(cFaces[opposedness.indices()[sz-1]]);
+    }
+}
+
+
+void Foam::extendedUpwindStencil::transportStencil
+(
+    const boolList& nonEmptyFace,
+    const labelListList& faceStencil,
+    const scalar minOpposedness,
+    const label faceI,
+    const label cellI,
+    const bool stencilHasNeighbour,
+
+    DynamicList<label>& oppositeFaces,
+    labelHashSet& faceStencilSet,
+    labelList& transportedStencil
+) const
+{
+    label globalOwn = faceStencil[faceI][0];
+    label globalNei = -1;
+    if (stencilHasNeighbour && faceStencil[faceI].size() >= 2)
+    {
+        globalNei = faceStencil[faceI][1];
+    }
+
+
+    selectOppositeFaces
+    (
+        nonEmptyFace,
+        minOpposedness,
+        faceI,
+        cellI,
+        oppositeFaces
+    );
+
+    // Collect all stencils of oppositefaces
+    faceStencilSet.clear();
+    forAll(oppositeFaces, i)
+    {
+        const labelList& fStencil = faceStencil[oppositeFaces[i]];
+
+        forAll(fStencil, j)
+        {
+            label globalI = fStencil[j];
+
+            if (globalI != globalOwn && globalI != globalNei)
+            {
+                faceStencilSet.insert(globalI);
+            }
+        }
+    }
+
+    // Add my owner and neighbour first.
+    if (stencilHasNeighbour)
+    {
+        transportedStencil.setSize(faceStencilSet.size()+2);
+        label n = 0;
+        transportedStencil[n++] = globalOwn;
+        transportedStencil[n++] = globalNei;
+
+        forAllConstIter(labelHashSet, faceStencilSet, iter)
+        {
+            if (iter.key() != globalOwn && iter.key() != globalNei)
+            {
+                transportedStencil[n++] = iter.key();
+            }
+        }
+        if (n != transportedStencil.size())
+        {
+            FatalErrorIn("extendedUpwindStencil::transportStencil(..)")
+                << "problem:" << faceStencilSet
+                << abort(FatalError);
+        }
+    }
+    else
+    {
+        transportedStencil.setSize(faceStencilSet.size()+1);
+        label n = 0;
+        transportedStencil[n++] = globalOwn;
+
+        forAllConstIter(labelHashSet, faceStencilSet, iter)
+        {
+            if (iter.key() != globalOwn)
+            {
+                transportedStencil[n++] = iter.key();
+            }
+        }
+        if (n != transportedStencil.size())
+        {
+            FatalErrorIn("extendedUpwindStencil::transportStencil(..)")
+                << "problem:" << faceStencilSet
+                << abort(FatalError);
+        }
+    }
+}
+
+
+void Foam::extendedUpwindStencil::transportStencils
+(
+    const labelListList& faceStencil,
+    const scalar minOpposedness,
+    labelListList& ownStencil,
+    labelListList& neiStencil
+)
+{
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+    // Work arrays
+    DynamicList<label> oppositeFaces;
+    labelHashSet faceStencilSet;
+
+
+    // For quick detection of empty faces
+    boolList nonEmptyFace(mesh_.nFaces(), true);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (isA<emptyPolyPatch>(pp))
+        {
+            label faceI = pp.start();
+            forAll(pp, i)
+            {
+                nonEmptyFace[faceI++] = false;
+            }
+        }
+    }
+
+
+    // Do the owner side
+    // ~~~~~~~~~~~~~~~~~
+    // stencil is synchronised at entry so no need to swap.
+
+    ownStencil.setSize(mesh_.nFaces());
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        transportStencil
+        (
+            nonEmptyFace,
+            faceStencil,
+            minOpposedness,
+            faceI,
+            own[faceI],
+            true,                   //stencilHasNeighbour
+            oppositeFaces,
+            faceStencilSet,
+            ownStencil[faceI]
+        );
+    }
+    // Boundary faces
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        label faceI = pp.start();
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                transportStencil
+                (
+                    nonEmptyFace,
+                    faceStencil,
+                    minOpposedness,
+                    faceI,
+                    own[faceI],
+                    true,                   //stencilHasNeighbour
+
+                    oppositeFaces,
+                    faceStencilSet,
+                    ownStencil[faceI]
+                );
+                faceI++;
+            }
+        }
+        else if (!isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                // faceStencil does not contain neighbour
+                transportStencil
+                (
+                    nonEmptyFace,
+                    faceStencil,
+                    minOpposedness,
+                    faceI,
+                    own[faceI],
+                    false,                  //stencilHasNeighbour
+
+                    oppositeFaces,
+                    faceStencilSet,
+                    ownStencil[faceI]
+                );
+                faceI++;
+            }
+        }
+    }
+
+
+    // Swap coupled boundary stencil
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    
+    labelListList neiBndStencil(nBnd);
+    for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+    {
+        neiBndStencil[faceI-mesh_.nInternalFaces()] = ownStencil[faceI];
+    }
+    syncTools::swapBoundaryFaceList(mesh_, neiBndStencil, false);
+
+
+    // Do the neighbour side
+    // ~~~~~~~~~~~~~~~~~~~~~
+    // - internal faces : get opposite faces on neighbour side
+    // - boundary faces : empty
+    // - coupled faces  : in neiBndStencil
+
+    neiStencil.setSize(mesh_.nFaces());
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        transportStencil
+        (
+            nonEmptyFace,
+            faceStencil,
+            minOpposedness,
+            faceI,
+            nei[faceI],
+            true,                   //stencilHasNeighbour
+
+            oppositeFaces,
+            faceStencilSet,
+            neiStencil[faceI]
+        );
+    }
+
+    // Boundary faces
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        label faceI = pp.start();
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                neiStencil[faceI].transfer
+                (
+                    neiBndStencil[faceI-mesh_.nInternalFaces()]
+                );
+                faceI++;
+            }
+        }
+        else
+        {
+            // Boundary has empty neighbour stencil
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::extendedUpwindStencil::extendedUpwindStencil
+(
+    const faceStencil& stencil,
+    const scalar minOpposedness
+)
+:
+    extendedStencil(stencil.mesh())
+{
+    //forAll(stencil, faceI)
+    //{
+    //    const labelList& fCells = stencil[faceI];
+    //
+    //    Pout<< "Face:" << faceI << " at:" << mesh_.faceCentres()[faceI]
+    //        << endl;
+    //
+    //    forAll(fCells, i)
+    //    {
+    //        label globalI = fCells[i];
+    //
+    //        if (globalI < mesh_.nCells())
+    //        {
+    //            Pout<< "    cell:" << globalI
+    //                << " at:" << mesh_.cellCentres()[globalI] << endl;
+    //        }
+    //        else
+    //        {
+    //            label faceI = globalI-mesh_.nCells() + mesh_.nInternalFaces();
+    //
+    //            Pout<< "    boundary:" << faceI
+    //                << " at:" << mesh_.faceCentres()[faceI] << endl;
+    //        }
+    //    }
+    //}
+    //Pout<< endl << endl;
+
+
+    // Transport centred stencil to upwind/downwind face
+    transportStencils
+    (
+        stencil,
+        minOpposedness,
+        ownStencil_,
+        neiStencil_
+    );
+
+    ownMapPtr_ = calcDistributeMap
+    (
+        stencil.globalNumbering(),
+        ownStencil_
+    );
+
+    neiMapPtr_ = calcDistributeMap
+    (
+        stencil.globalNumbering(),
+        neiStencil_
+    );
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencil.H

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::extendedUpwindStencil
+
+Description
+
+SourceFiles
+    extendedUpwindStencil.C
+    extendedUpwindStencilTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef extendedUpwindStencil_H
+#define extendedUpwindStencil_H
+
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class faceStencil;
+
+/*---------------------------------------------------------------------------*\
+                           Class extendedUpwindStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class extendedUpwindStencil
+:
+    public extendedStencil
+{
+    // Private data
+
+        //- Swap map for getting neigbouring data
+        autoPtr<mapDistribute> ownMapPtr_;
+        autoPtr<mapDistribute> neiMapPtr_;
+
+        //- Per face the stencil.
+        labelListList ownStencil_;
+        labelListList neiStencil_;
+
+
+
+    // Private Member Functions
+
+        //- Find most 'opposite' faces of cell
+        void selectOppositeFaces
+        (
+            const boolList& nonEmptyFace,
+            const scalar minOpposedness,
+            const label faceI,
+            const label cellI,
+            DynamicList<label>& oppositeFaces
+        ) const;
+
+        //- Transport (centred) face stencil to 'opposite' face.
+        void transportStencil
+        (
+            const boolList& nonEmptyFace,
+            const labelListList& faceStencil,
+            const scalar minOpposedness,
+            const label faceI,
+            const label cellI,
+            const bool stencilHasNeighbour,
+
+            DynamicList<label>& oppositeFaces,
+            labelHashSet& faceStencilSet,
+            labelList& transportedStencil
+        ) const;
+
+        //- Transport (centred) face stencil to 'opposite' faces.
+        void transportStencils
+        (
+            const labelListList& faceStencil,
+            const scalar minOpposedness,
+            labelListList& ownStencil,
+            labelListList& neiStencil
+        );
+
+
+        //- Disallow default bitwise copy construct
+        extendedUpwindStencil(const extendedUpwindStencil&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const extendedUpwindStencil&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh and uncompacted face stencil
+        extendedUpwindStencil
+        (
+            const faceStencil&,
+            const scalar minOpposedness
+        );
+
+
+    // Member Functions
+
+        //- Return reference to the parallel distribution map
+        const mapDistribute& ownMap() const
+        {
+            return ownMapPtr_();
+        }
+
+        //- Return reference to the parallel distribution map
+        const mapDistribute& neiMap() const
+        {
+            return neiMapPtr_();
+        }
+
+        //- Return reference to the stencil
+        const labelListList& ownStencil() const
+        {
+            return ownStencil_;
+        }
+
+        //- Return reference to the stencil
+        const labelListList& neiStencil() const
+        {
+            return neiStencil_;
+        }
+
+        //- Sum vol field contributions to create face values
+        template<class Type>
+        tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > weightedSum
+        (
+            const surfaceScalarField& phi,
+            const GeometricField<Type, fvPatchField, volMesh>& fld,
+            const List<List<scalar> >& ownWeights,
+            const List<List<scalar> >& neiWeights
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "extendedUpwindStencilTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/extendedUpwindStencilTemplates.C

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh> >
+Foam::extendedUpwindStencil::weightedSum
+(
+    const surfaceScalarField& phi,
+    const GeometricField<Type, fvPatchField, volMesh>& fld,
+    const List<List<scalar> >& ownWeights,
+    const List<List<scalar> >& neiWeights
+) const
+{
+    const fvMesh& mesh = fld.mesh();
+
+    // Collect internal and boundary values
+    List<List<Type> > ownFld;
+    collectData(ownMap(), ownStencil(), fld, ownFld);
+    List<List<Type> > neiFld;
+    collectData(neiMap(), neiStencil(), fld, neiFld);
+
+    tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsfCorr
+    (
+        new GeometricField<Type, fvsPatchField, surfaceMesh>
+        (
+            IOobject
+            (
+                fld.name(),
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensioned<Type>
+            (
+                fld.name(),
+                fld.dimensions(),
+                pTraits<Type>::zero
+            )
+        )
+    );
+    GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsfCorr();
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        if (phi[faceI] > 0)
+        {
+            // Flux out of owner. Use upwind (= owner side) stencil.
+            const List<Type>& stField = ownFld[faceI];
+            const List<scalar>& stWeight = ownWeights[faceI];
+
+            forAll(stField, i)
+            {
+                sf[faceI] += stField[i]*stWeight[i];
+            }
+        }
+        else
+        {
+            const List<Type>& stField = neiFld[faceI];
+            const List<scalar>& stWeight = neiWeights[faceI];
+
+            forAll(stField, i)
+            {
+                sf[faceI] += stField[i]*stWeight[i];
+            }
+        }
+    }
+
+    // Boundaries. Either constrained or calculated so assign value
+    // directly (instead of nicely using operator==)
+    typename GeometricField<Type, fvsPatchField, surfaceMesh>::
+        GeometricBoundaryField& bSfCorr = sf.boundaryField();
+
+    forAll(bSfCorr, patchi)
+    {
+        fvsPatchField<Type>& pSfCorr = bSfCorr[patchi];
+
+        label faceI = pSfCorr.patch().patch().start();
+
+        forAll(pSfCorr, i)
+        {
+            if (phi[faceI] > 0)
+            {
+                // Flux out of owner. Use upwind (= owner side) stencil.
+                const List<Type>& stField = ownFld[faceI];
+                const List<scalar>& stWeight = ownWeights[faceI];
+
+                forAll(stField, j)
+                {
+                    pSfCorr[i] += stField[j]*stWeight[j];
+                }
+            }
+            else
+            {
+                const List<Type>& stField = neiFld[faceI];
+                const List<scalar>& stWeight = neiWeights[faceI];
+
+                forAll(stField, j)
+                {
+                    pSfCorr[i] += stField[j]*stWeight[j];
+                }
+            }
+            faceI++;
+        }
+    }
+
+    return tsfCorr;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellEdgeCellStencil.C

@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "cellEdgeCellStencil.H"
+#include "syncTools.H"
+//#include "meshTools.H"
+//#include "OFstream.H"
+//#include "Time.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Calculates per edge the neighbour data (= edgeCells)
+void Foam::cellEdgeCellStencil::calcEdgeBoundaryData
+(
+    const boolList& isValidBFace,
+    const labelList& boundaryEdges,
+    EdgeMap<labelList>& neiGlobal
+) const
+{
+    neiGlobal.resize(2*boundaryEdges.size());
+
+    labelHashSet edgeGlobals;
+
+    forAll(boundaryEdges, i)
+    {
+        label edgeI = boundaryEdges[i];
+
+        neiGlobal.insert
+        (
+            mesh().edges()[edgeI],
+            calcFaceCells
+            (
+                isValidBFace,
+                mesh().edgeFaces(edgeI),
+                edgeGlobals
+            )
+        );
+    }
+
+    syncTools::syncEdgeMap
+    (
+        mesh(),
+        neiGlobal,
+        unionEqOp(),
+        false           // apply separation
+    );
+}
+
+
+// Calculates per cell the neighbour data (= cell or boundary in global
+// numbering). First element is always cell itself!
+void Foam::cellEdgeCellStencil::calcCellStencil
+(
+    labelListList& globalCellCells
+) const
+{
+    // Calculate edges on coupled patches
+    labelList boundaryEdges
+    (
+        allCoupledFacesPatch()().meshEdges
+        (
+            mesh().edges(),
+            mesh().pointEdges()
+        )
+    );
+
+    //{
+    //    OFstream str(mesh().time().path()/"boundaryEdges.obj");
+    //    Pout<< "DUmping boundary edges to " << str.name() << endl;
+    //
+    //    label vertI = 0;
+    //    forAll(boundaryEdges, i)
+    //    {
+    //        label edgeI = boundaryEdges[i];
+    //        const edge& e = mesh().edges()[edgeI];
+    //        const point& p0 = mesh().points()[e[0]];
+    //        const point& p1 = mesh().points()[e[1]];
+    //
+    //        Pout<< "boundary edge " << edgeI << " between " << p0 << p1
+    //            << endl;
+    //
+    //        meshTools::writeOBJ(str, p0);
+    //        vertI++;
+    //        meshTools::writeOBJ(str, p1);
+    //        vertI++;
+    //        str << "l " << vertI-1 << ' ' << vertI << nl;
+    //    }
+    //}
+
+
+    // Mark boundary faces to be included in stencil (i.e. not coupled or empty)
+    boolList isValidBFace;
+    validBoundaryFaces(isValidBFace);
+
+
+    // Swap edgeCells for coupled edges. Note: use EdgeMap for now since we've
+    // got syncTools::syncEdgeMap for those. Should be replaced with Map and
+    // syncTools functionality to handle those.
+    EdgeMap<labelList> neiGlobal;
+    calcEdgeBoundaryData
+    (
+        isValidBFace,
+        boundaryEdges,
+        neiGlobal
+    );
+
+    globalCellCells.setSize(mesh().nCells());
+
+    // Do coupled edges first
+
+    forAll(boundaryEdges, i)
+    {
+        label edgeI = boundaryEdges[i];
+
+        const labelList& eGlobals = neiGlobal[mesh().edges()[edgeI]];
+
+        // Distribute to all edgeCells
+        const labelList& eCells = mesh().edgeCells(edgeI);
+
+        forAll(eCells, j)
+        {
+            label cellI = eCells[j];
+
+            // Insert pGlobals into globalCellCells
+            merge
+            (
+                globalNumbering().toGlobal(cellI),
+                eGlobals,
+                globalCellCells[cellI]
+            );
+        }
+    }
+    neiGlobal.clear();
+
+    // Do remaining edges cells
+    labelHashSet edgeGlobals;
+
+    for (label edgeI = 0; edgeI < mesh().nEdges(); edgeI++)
+    {
+        labelList eGlobals
+        (
+            calcFaceCells
+            (
+                isValidBFace,
+                mesh().edgeFaces(edgeI),
+                edgeGlobals
+            )
+        );
+
+        const labelList& eCells = mesh().edgeCells(edgeI);
+
+        forAll(eCells, j)
+        {
+            label cellI = eCells[j];
+
+            merge
+            (
+                globalNumbering().toGlobal(cellI),
+                eGlobals,
+                globalCellCells[cellI]
+            );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cellEdgeCellStencil::cellEdgeCellStencil(const polyMesh& mesh)
+:
+    faceStencil(mesh)
+{
+    // Calculate per cell the (edge) connected cells (in global numbering)
+    labelListList globalCellCells;
+    calcCellStencil(globalCellCells);
+
+    // Add stencils of neighbouring cells to create faceStencil
+    labelListList faceStencil;
+    calcFaceStencil(globalCellCells, faceStencil);
+
+    // Transfer to *this
+    transfer(faceStencil);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellEdgeCellStencil.H

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::cellEdgeCellStencil
+
+Description
+
+SourceFiles
+    cellEdgeCellStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellEdgeCellStencil_H
+#define cellEdgeCellStencil_H
+
+#include "faceStencil.H"
+#include "boolList.H"
+#include "labelHashSet.H"
+#include "Map.H"
+#include "EdgeMap.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class cellEdgeCellStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellEdgeCellStencil
+:
+    public faceStencil
+{
+    // Private Member Functions
+
+        //- Calculates per edge the neighbour data (= edgeCells)
+        void calcEdgeBoundaryData
+        (
+            const boolList& isValidBFace,
+            const labelList& boundaryEdges,
+            EdgeMap<labelList>& neiGlobal
+        ) const;
+
+        void calcCellStencil(labelListList& globalCellCells) const;
+
+
+        //- Disallow default bitwise copy construct
+        cellEdgeCellStencil(const cellEdgeCellStencil&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const cellEdgeCellStencil&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from all cells and boundary faces
+        explicit cellEdgeCellStencil(const polyMesh&);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.C

@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "cellFaceCellStencil.H"
+#include "syncTools.H"
+#include "SortableList.H"
+#include "emptyPolyPatch.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Calculates per face the neighbour data (= cell or boundary face)
+void Foam::cellFaceCellStencil::calcFaceBoundaryData
+(
+    labelList& neiGlobal
+) const
+{
+    const polyBoundaryMesh& patches = mesh().boundaryMesh();
+    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+    const labelList& own = mesh().faceOwner();
+
+    neiGlobal.setSize(nBnd);
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        label faceI = pp.start();
+
+        if (pp.coupled())
+        {
+            // For coupled faces get the cell on the other side
+            forAll(pp, i)
+            {
+                label bFaceI = faceI-mesh().nInternalFaces(); 
+                neiGlobal[bFaceI] = globalNumbering().toGlobal(own[faceI]);
+                faceI++;
+            }
+        }
+        else if (isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                label bFaceI = faceI-mesh().nInternalFaces(); 
+                neiGlobal[bFaceI] = -1;
+                faceI++;
+            }
+        }
+        else
+        {
+            // For noncoupled faces get the boundary face.
+            forAll(pp, i)
+            {
+                label bFaceI = faceI-mesh().nInternalFaces(); 
+                neiGlobal[bFaceI] =
+                    globalNumbering().toGlobal(mesh().nCells()+bFaceI);
+                faceI++;
+            }
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh(), neiGlobal, false);
+}
+
+
+// Calculates per cell the neighbour data (= cell or boundary in global
+// numbering). First element is always cell itself!
+void Foam::cellFaceCellStencil::calcCellStencil(labelListList& globalCellCells)
+ const
+{
+    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+    const labelList& own = mesh().faceOwner();
+    const labelList& nei = mesh().faceNeighbour();
+
+
+    // Calculate coupled neighbour (in global numbering)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList neiGlobal(nBnd);
+    calcFaceBoundaryData(neiGlobal);
+
+
+    // Determine cellCells in global numbering
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    globalCellCells.setSize(mesh().nCells());
+    forAll(globalCellCells, cellI)
+    {
+        const cell& cFaces = mesh().cells()[cellI];
+
+        labelList& cCells = globalCellCells[cellI];
+
+        cCells.setSize(cFaces.size()+1);
+
+        label nNbr = 0;
+
+        // Myself
+        cCells[nNbr++] = globalNumbering().toGlobal(cellI);
+
+        // Collect neighbouring cells/faces
+        forAll(cFaces, i)
+        {
+            label faceI = cFaces[i];
+
+            if (mesh().isInternalFace(faceI))
+            {
+                label nbrCellI = own[faceI];
+                if (nbrCellI == cellI)
+                {
+                    nbrCellI = nei[faceI];
+                }
+                cCells[nNbr++] = globalNumbering().toGlobal(nbrCellI);
+            }
+            else
+            {
+                label nbrCellI = neiGlobal[faceI-mesh().nInternalFaces()];
+                if (nbrCellI != -1)
+                {
+                    cCells[nNbr++] = nbrCellI;
+                }
+            }
+        }
+        cCells.setSize(nNbr);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cellFaceCellStencil::cellFaceCellStencil(const polyMesh& mesh)
+:
+    faceStencil(mesh)
+{
+    // Calculate per cell the (face) connected cells (in global numbering)
+    labelListList globalCellCells;
+    calcCellStencil(globalCellCells);
+
+    // Add stencils of neighbouring cells to create faceStencil
+    labelListList faceStencil;
+    calcFaceStencil(globalCellCells, faceStencil);
+
+    // Transfer to *this
+    transfer(faceStencil);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellFaceCellStencil.H

@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::cellFaceCellStencil
+
+Description
+
+SourceFiles
+    cellFaceCellStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellFaceCellStencil_H
+#define cellFaceCellStencil_H
+
+#include "faceStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class cellFaceCellStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellFaceCellStencil
+:
+    public faceStencil
+{
+    // Private Member Functions
+
+        void calcFaceBoundaryData(labelList& neiGlobal) const;
+
+        void calcCellStencil(labelListList& globalCellCells) const;
+
+        //- Disallow default bitwise copy construct
+        cellFaceCellStencil(const cellFaceCellStencil&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const cellFaceCellStencil&);
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh
+        explicit cellFaceCellStencil(const polyMesh& mesh);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellPointCellStencil.C

@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "cellPointCellStencil.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Calculates per point the neighbour data (= pointCells)
+void Foam::cellPointCellStencil::calcPointBoundaryData
+(
+    const boolList& isValidBFace,
+    const labelList& boundaryPoints,
+    Map<labelList>& neiGlobal
+) const
+{
+    neiGlobal.resize(2*boundaryPoints.size());
+
+    labelHashSet pointGlobals;
+
+    forAll(boundaryPoints, i)
+    {
+        label pointI = boundaryPoints[i];
+
+        neiGlobal.insert
+        (
+            pointI,
+            calcFaceCells
+            (
+                isValidBFace,
+                mesh().pointFaces()[pointI],
+                pointGlobals
+            )
+        );
+    }
+
+    syncTools::syncPointMap
+    (
+        mesh(),
+        neiGlobal,
+        unionEqOp(),
+        false           // apply separation
+    );
+}
+
+
+// Calculates per cell the neighbour data (= cell or boundary in global
+// numbering). First element is always cell itself!
+void Foam::cellPointCellStencil::calcCellStencil
+(
+    labelListList& globalCellCells
+) const
+{
+    // Calculate points on coupled patches
+    labelList boundaryPoints(allCoupledFacesPatch()().meshPoints());
+
+
+    // Mark boundary faces to be included in stencil (i.e. not coupled or empty)
+    boolList isValidBFace;
+    validBoundaryFaces(isValidBFace);
+
+
+    // Swap pointCells for coupled points
+    Map<labelList> neiGlobal;
+    calcPointBoundaryData
+    (
+        isValidBFace,
+        boundaryPoints,
+        neiGlobal
+    );
+
+    globalCellCells.setSize(mesh().nCells());
+
+    // Do coupled points first
+
+    forAll(boundaryPoints, i)
+    {
+        label pointI = boundaryPoints[i];
+
+        const labelList& pGlobals = neiGlobal[pointI];
+
+        // Distribute to all pointCells
+        const labelList& pCells = mesh().pointCells(pointI);
+
+        forAll(pCells, j)
+        {
+            label cellI = pCells[j];
+
+            // Insert pGlobals into globalCellCells
+            merge
+            (
+                globalNumbering().toGlobal(cellI),
+                pGlobals,
+                globalCellCells[cellI]
+            );
+        }
+    }
+    neiGlobal.clear();
+
+    // Do remaining points cells
+    labelHashSet pointGlobals;
+
+    for (label pointI = 0; pointI < mesh().nPoints(); pointI++)
+    {
+        labelList pGlobals
+        (
+            calcFaceCells
+            (
+                isValidBFace,
+                mesh().pointFaces()[pointI],
+                pointGlobals
+            )
+        );
+
+        const labelList& pCells = mesh().pointCells(pointI);
+
+        forAll(pCells, j)
+        {
+            label cellI = pCells[j];
+
+            merge
+            (
+                globalNumbering().toGlobal(cellI),
+                pGlobals,
+                globalCellCells[cellI]
+            );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cellPointCellStencil::cellPointCellStencil(const polyMesh& mesh)
+:
+    faceStencil(mesh)
+{
+    // Calculate per cell the (point) connected cells (in global numbering)
+    labelListList globalCellCells;
+    calcCellStencil(globalCellCells);
+
+    // Add stencils of neighbouring cells to create faceStencil
+    labelListList faceStencil;
+    calcFaceStencil(globalCellCells, faceStencil);
+
+    // Transfer to *this
+    transfer(faceStencil);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/cellPointCellStencil.H

@@ -0,0 +1,94 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::cellPointCellStencil
+
+Description
+
+SourceFiles
+    cellPointCellStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cellPointCellStencil_H
+#define cellPointCellStencil_H
+
+#include "faceStencil.H"
+#include "boolList.H"
+#include "labelHashSet.H"
+#include "Map.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class cellPointCellStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class cellPointCellStencil
+:
+    public faceStencil
+{
+    // Private Member Functions
+
+        //- Calculates per point the neighbour data (= pointCells)
+        void calcPointBoundaryData
+        (
+            const boolList& isValidBFace,
+            const labelList& boundaryPoints,
+            Map<labelList>& neiGlobal
+        ) const;
+
+        void calcCellStencil(labelListList& globalCellCells) const;
+
+
+        //- Disallow default bitwise copy construct
+        cellPointCellStencil(const cellPointCellStencil&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const cellPointCellStencil&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from all cells and boundary faces
+        explicit cellPointCellStencil(const polyMesh&);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/faceEdgeCellStencil.C

@@ -0,0 +1,445 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "faceEdgeCellStencil.H"
+#include "syncTools.H"
+#include "emptyPolyPatch.H"
+//#include "meshTools.H"
+//#include "OFstream.H"
+//#include "Time.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Calculates per edge the neighbour data (= edgeCells)
+void Foam::faceEdgeCellStencil::calcEdgeBoundaryData
+(
+    const boolList& isValidBFace,
+    const labelList& boundaryEdges,
+    EdgeMap<labelList>& neiGlobal
+) const
+{
+    neiGlobal.resize(2*boundaryEdges.size());
+
+    labelHashSet edgeGlobals;
+
+    forAll(boundaryEdges, i)
+    {
+        label edgeI = boundaryEdges[i];
+
+        neiGlobal.insert
+        (
+            mesh().edges()[edgeI],
+            calcFaceCells
+            (
+                isValidBFace,
+                mesh().edgeFaces(edgeI),
+                edgeGlobals
+            )
+        );
+    }
+
+    syncTools::syncEdgeMap
+    (
+        mesh(),
+        neiGlobal,
+        unionEqOp(),
+        false           // apply separation
+    );
+}
+
+
+// Calculates per face the edge connected data (= cell or boundary in global
+// numbering).
+void Foam::faceEdgeCellStencil::calcFaceStencil
+(
+    labelListList& faceStencil
+) const
+{
+    const polyBoundaryMesh& patches = mesh().boundaryMesh();
+    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+    const labelList& own = mesh().faceOwner();
+    const labelList& nei = mesh().faceNeighbour();
+
+
+
+    // Determine neighbouring global cell
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList neiGlobalCell(nBnd);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                neiGlobalCell[faceI-mesh().nInternalFaces()] =
+                    globalNumbering().toGlobal(own[faceI]);
+                faceI++;
+            }
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh(), neiGlobalCell, false);
+
+
+
+    // Determine on coupled edges the edge cells on the other side
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Calculate edges on coupled patches
+    labelList boundaryEdges
+    (
+        allCoupledFacesPatch()().meshEdges
+        (
+            mesh().edges(),
+            mesh().pointEdges()
+        )
+    );
+
+    // Mark boundary faces to be included in stencil (i.e. not coupled or empty)
+    boolList isValidBFace;
+    validBoundaryFaces(isValidBFace);
+
+    // Swap edgeCells for coupled edges. Note: use EdgeMap for now since we've
+    // got syncTools::syncEdgeMap for those. Should be replaced with Map and
+    // syncTools functionality to handle those.
+    EdgeMap<labelList> neiGlobal;
+    calcEdgeBoundaryData
+    (
+        isValidBFace,
+        boundaryEdges,
+        neiGlobal
+    );
+
+
+
+    // Construct stencil in global numbering
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    faceStencil.setSize(mesh().nFaces());
+
+    // Do coupled edges first
+
+    forAll(boundaryEdges, i)
+    {
+        label edgeI = boundaryEdges[i];
+
+        const labelList& eGlobals = neiGlobal[mesh().edges()[edgeI]];
+
+        // Distribute to all edgeFaces
+        const labelList& eFaces = mesh().edgeFaces(edgeI);
+
+        forAll(eFaces, j)
+        {
+            label faceI = eFaces[j];
+
+            // Insert eGlobals into faceStencil.
+            merge(-1, -1, eGlobals, faceStencil[faceI]);
+        }
+    }
+    neiGlobal.clear();
+
+
+    // Do remaining edges by looping over all faces
+
+    // Work arrays
+    DynamicList<label> fEdgesSet;
+    DynamicList<label> eFacesSet;
+    labelHashSet faceStencilSet;
+
+    for (label faceI = 0; faceI < mesh().nInternalFaces(); faceI++)
+    {
+        label globalOwn = globalNumbering().toGlobal(own[faceI]);
+        label globalNei = globalNumbering().toGlobal(nei[faceI]);
+
+        // Convert any existing faceStencil (from coupled edges) into
+        // set and operate on this.
+
+        faceStencilSet.clear();
+
+        // Insert all but global owner and neighbour
+        forAll(faceStencil[faceI], i)
+        {
+            label globalI = faceStencil[faceI][i];
+            if (globalI != globalOwn && globalI != globalNei)
+            {
+                faceStencilSet.insert(globalI);
+            }
+        }
+        faceStencil[faceI].clear();
+
+        // Collect all edge connected (internal) cells
+        const labelList& fEdges = mesh().faceEdges(faceI, fEdgesSet);
+
+        forAll(fEdges, i)
+        {
+            label edgeI = fEdges[i];
+
+            insertFaceCells
+            (
+                globalOwn,
+                globalNei,
+                isValidBFace,
+                mesh().edgeFaces(edgeI, eFacesSet),
+                faceStencilSet
+            );
+        }
+
+        // Extract, guarantee owner first, neighbour second.
+        faceStencil[faceI].setSize(faceStencilSet.size()+2);
+        label n = 0;
+        faceStencil[faceI][n++] = globalOwn;
+        faceStencil[faceI][n++] = globalNei;
+        forAllConstIter(labelHashSet, faceStencilSet, iter)
+        {
+            if (iter.key() == globalOwn || iter.key() == globalNei)
+            {
+                FatalErrorIn("faceEdgeCellStencil::calcFaceStencil(..)")
+                    << "problem:" << faceStencilSet
+                    << abort(FatalError);
+            }
+            faceStencil[faceI][n++] = iter.key();
+        }
+    }
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        label faceI = pp.start();
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                label globalOwn = globalNumbering().toGlobal(own[faceI]);
+                label globalNei = neiGlobalCell[faceI-mesh().nInternalFaces()];
+
+                // Convert any existing faceStencil (from coupled edges) into
+                // set and operate on this.
+
+                faceStencilSet.clear();
+
+                // Insert all but global owner and neighbour
+                forAll(faceStencil[faceI], i)
+                {
+                    label globalI = faceStencil[faceI][i];
+                    if (globalI != globalOwn && globalI != globalNei)
+                    {
+                        faceStencilSet.insert(globalI);
+                    }
+                }
+                faceStencil[faceI].clear();
+
+                // Collect all edge connected (internal) cells
+                const labelList& fEdges = mesh().faceEdges(faceI, fEdgesSet);
+
+                forAll(fEdges, i)
+                {
+                    label edgeI = fEdges[i];
+
+                    insertFaceCells
+                    (
+                        globalOwn,
+                        globalNei,
+                        isValidBFace,
+                        mesh().edgeFaces(edgeI, eFacesSet),
+                        faceStencilSet
+                    );
+                }
+
+                // Extract, guarantee owner first, neighbour second.
+                faceStencil[faceI].setSize(faceStencilSet.size()+2);
+                label n = 0;
+                faceStencil[faceI][n++] = globalOwn;
+                faceStencil[faceI][n++] = globalNei;
+                forAllConstIter(labelHashSet, faceStencilSet, iter)
+                {
+                    if (iter.key() == globalOwn || iter.key() == globalNei)
+                    {
+                        FatalErrorIn("faceEdgeCellStencil::calcFaceStencil(..)")
+                            << "problem:" << faceStencilSet
+                            << abort(FatalError);
+                    }
+                    faceStencil[faceI][n++] = iter.key();
+                }
+
+                if (n != faceStencil[faceI].size())
+                {
+                    FatalErrorIn("problem") << "n:" << n
+                        << " size:" << faceStencil[faceI].size()
+                        << abort(FatalError);
+                }
+
+                faceI++;
+            }
+        }
+        else if (!isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                label globalOwn = globalNumbering().toGlobal(own[faceI]);
+
+                // Convert any existing faceStencil (from coupled edges) into
+                // set and operate on this.
+
+                faceStencilSet.clear();
+
+                // Insert all but global owner and neighbour
+                forAll(faceStencil[faceI], i)
+                {
+                    label globalI = faceStencil[faceI][i];
+                    if (globalI != globalOwn)
+                    {
+                        faceStencilSet.insert(globalI);
+                    }
+                }
+                faceStencil[faceI].clear();
+
+                // Collect all edge connected (internal) cells
+                const labelList& fEdges = mesh().faceEdges(faceI, fEdgesSet);
+
+                forAll(fEdges, i)
+                {
+                    label edgeI = fEdges[i];
+
+                    insertFaceCells
+                    (
+                        globalOwn,
+                        -1,
+                        isValidBFace,
+                        mesh().edgeFaces(edgeI, eFacesSet),
+                        faceStencilSet
+                    );
+                }
+
+                // Extract, guarantee owner first, neighbour second.
+                faceStencil[faceI].setSize(faceStencilSet.size()+1);
+                label n = 0;
+                faceStencil[faceI][n++] = globalOwn;
+                forAllConstIter(labelHashSet, faceStencilSet, iter)
+                {
+                    if (iter.key() == globalOwn)
+                    {
+                        FatalErrorIn("faceEdgeCellStencil::calcFaceStencil(..)")
+                            << "problem:" << faceStencilSet
+                            << abort(FatalError);
+                    }
+                    faceStencil[faceI][n++] = iter.key();
+                }
+
+                faceI++;
+            }
+        }
+    }
+
+
+    for (label faceI = 0; faceI < mesh().nInternalFaces(); faceI++)
+    {
+        label globalOwn = globalNumbering().toGlobal(own[faceI]);
+        if (faceStencil[faceI][0] != globalOwn)
+        {
+            FatalErrorIn("faceEdgeCellStencil::calcFaceStencil(..)")
+                << "problem:" << faceStencil[faceI]
+                << " globalOwn:" << globalOwn
+                << abort(FatalError);
+        }
+        label globalNei = globalNumbering().toGlobal(nei[faceI]);
+        if (faceStencil[faceI][1] != globalNei)
+        {
+            FatalErrorIn("faceEdgeCellStencil::calcFaceStencil(..)")
+                << "problem:" << faceStencil[faceI]
+                << " globalNei:" << globalNei
+                << abort(FatalError);
+        }
+    }
+
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+
+                label globalOwn = globalNumbering().toGlobal(own[faceI]);
+                if (faceStencil[faceI][0] != globalOwn)
+                {
+                    FatalErrorIn("faceEdgeCellStencil::calcFaceStencil(..)")
+                        << "problem:" << faceStencil[faceI]
+                        << " globalOwn:" << globalOwn
+                        << abort(FatalError);
+                }
+                label globalNei = neiGlobalCell[faceI-mesh().nInternalFaces()];
+                if (faceStencil[faceI][1] != globalNei)
+                {
+                    FatalErrorIn("faceEdgeCellStencil::calcFaceStencil(..)")
+                        << "problem:" << faceStencil[faceI]
+                        << " globalNei:" << globalNei
+                        << abort(FatalError);
+                }
+            }
+        }
+        else if (!isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+
+                label globalOwn = globalNumbering().toGlobal(own[faceI]);
+                if (faceStencil[faceI][0] != globalOwn)
+                {
+                    FatalErrorIn("faceEdgeCellStencil::calcFaceStencil(..)")
+                        << "problem:" << faceStencil[faceI]
+                        << " globalOwn:" << globalOwn
+                        << abort(FatalError);
+                }
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faceEdgeCellStencil::faceEdgeCellStencil(const polyMesh& mesh)
+:
+    faceStencil(mesh)
+{
+    // Calculate per face the (edge) connected cells (in global numbering)
+    labelListList faceStencil;
+    calcFaceStencil(faceStencil);
+
+    // Transfer to *this
+    transfer(faceStencil);
+}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/faceEdgeCellStencil.H

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::faceEdgeCellStencil
+
+Description
+    All cells connected via edge to face.
+
+SourceFiles
+    faceEdgeCellStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceEdgeCellStencil_H
+#define faceEdgeCellStencil_H
+
+#include "faceStencil.H"
+#include "boolList.H"
+#include "labelHashSet.H"
+#include "Map.H"
+#include "EdgeMap.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class faceEdgeCellStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceEdgeCellStencil
+:
+    public faceStencil
+{
+    // Private Member Functions
+
+        //- Calculates per edge the neighbour data (= edgeCells)
+        void calcEdgeBoundaryData
+        (
+            const boolList& isValidBFace,
+            const labelList& boundaryEdges,
+            EdgeMap<labelList>& neiGlobal
+        ) const;
+
+        void calcFaceStencil(labelListList& faceStencil) const;
+
+
+        //- Disallow default bitwise copy construct
+        faceEdgeCellStencil(const faceEdgeCellStencil&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const faceEdgeCellStencil&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from all cells and boundary faces
+        explicit faceEdgeCellStencil(const polyMesh&);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/faceStencil.C

@@ -0,0 +1,512 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "faceStencil.H"
+#include "syncTools.H"
+#include "SortableList.H"
+#include "emptyPolyPatch.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Merge two list and guarantee global0,global1 are first.
+void Foam::faceStencil::merge
+(
+    const label global0,
+    const label global1,
+    const labelList& listA,
+    labelList& listB
+)
+{
+    sort(listB);
+
+    // See if global0, global1 already present in listB
+    label nGlobalInsert = 0;
+
+    if (global0 != -1)
+    {
+        label index0 = findSortedIndex(listB, global0);
+        if (index0 == -1)
+        {
+            nGlobalInsert++;
+        }
+    }
+
+    if (global1 != -1)
+    {
+        label index1 = findSortedIndex(listB, global1);
+        if (index1 == -1)
+        {
+            nGlobalInsert++;
+        }
+    }
+
+
+    // For all in listA see if they are present
+    label nInsert = 0;
+
+    forAll(listA, i)
+    {
+        label elem = listA[i];
+
+        if (elem != global0 && elem != global1)
+        {
+            if (findSortedIndex(listB, elem) == -1)
+            {
+                nInsert++;
+            }
+        }
+    }
+
+    // Extend B with nInsert and whether global0,global1 need to be inserted.
+    labelList result(listB.size() + nGlobalInsert + nInsert);
+
+    label resultI = 0;
+
+    // Insert global0,1 first
+    if (global0 != -1)
+    {
+        result[resultI++] = global0;
+    }
+    if (global1 != -1)
+    {
+        result[resultI++] = global1;
+    }
+
+
+    // Insert listB
+    forAll(listB, i)
+    {
+        label elem = listB[i];
+
+        if (elem != global0 && elem != global1)
+        {
+            result[resultI++] = elem;
+        }
+    }
+
+
+    // Insert listA
+    forAll(listA, i)
+    {
+        label elem = listA[i];
+
+        if (elem != global0 && elem != global1)
+        {
+            if (findSortedIndex(listB, elem) == -1)
+            {
+                result[resultI++] = elem;
+            }
+        }
+    }
+
+    if (resultI != result.size())
+    {
+        FatalErrorIn("faceStencil::merge(..)")
+            << "problem" << abort(FatalError);
+    }
+
+    listB.transfer(result);
+}
+
+
+// Merge two list and guarantee globalI is first.
+void Foam::faceStencil::merge
+(
+    const label globalI,
+    const labelList& pGlobals,
+    labelList& cCells
+)
+{
+    labelHashSet set;
+    forAll(cCells, i)
+    {
+        if (cCells[i] != globalI)
+        {
+            set.insert(cCells[i]);
+        }
+    }
+
+    forAll(pGlobals, i)
+    {
+        if (pGlobals[i] != globalI)
+        {
+            set.insert(pGlobals[i]);
+        }
+    }
+
+    cCells.setSize(set.size()+1);
+    label n = 0;
+    cCells[n++] = globalI;
+
+    forAllConstIter(labelHashSet, set, iter)
+    {
+        cCells[n++] = iter.key();
+    }
+}
+
+
+void Foam::faceStencil::validBoundaryFaces(boolList& isValidBFace) const
+{
+    const polyBoundaryMesh& patches = mesh().boundaryMesh();
+
+    isValidBFace.setSize(mesh().nFaces()-mesh().nInternalFaces(), true);
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled() || isA<emptyPolyPatch>(pp))
+        {
+            label bFaceI = pp.start()-mesh().nInternalFaces();
+            forAll(pp, i)
+            {
+                isValidBFace[bFaceI++] = false;
+            }
+        }
+    }
+}
+
+
+Foam::autoPtr<Foam::indirectPrimitivePatch>
+Foam::faceStencil::allCoupledFacesPatch() const
+{
+    const polyBoundaryMesh& patches = mesh().boundaryMesh();
+
+    label nCoupled = 0;
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            nCoupled += pp.size();
+        }
+    }
+    labelList coupledFaces(nCoupled);
+    nCoupled = 0;
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                coupledFaces[nCoupled++] = faceI++;
+            }
+        }
+    }
+
+    return autoPtr<indirectPrimitivePatch>
+    (
+        new indirectPrimitivePatch
+        (
+            IndirectList<face>
+            (
+                mesh().faces(),
+                coupledFaces
+            ),
+            mesh().points()
+        )
+    );
+}
+
+
+void Foam::faceStencil::unionEqOp::operator()
+(
+    labelList& x,
+    const labelList& y
+) const
+{
+    if (y.size() > 0)
+    {
+        if (x.size() == 0)
+        {
+            x = y;
+        }
+        else
+        {
+            labelHashSet set(x);
+            forAll(y, i)
+            {
+                set.insert(y[i]);
+            }
+            x = set.toc();
+        }
+    }
+}
+
+
+void Foam::faceStencil::insertFaceCells
+(
+    const label exclude0,
+    const label exclude1,
+    const boolList& isValidBFace,
+    const labelList& faceLabels,
+    labelHashSet& globals
+) const
+{
+    const labelList& own = mesh().faceOwner();
+    const labelList& nei = mesh().faceNeighbour();
+
+    forAll(faceLabels, i)
+    {
+        label faceI = faceLabels[i];
+
+        label globalOwn = globalNumbering().toGlobal(own[faceI]);
+        if (globalOwn != exclude0 && globalOwn != exclude1)
+        {
+            globals.insert(globalOwn);
+        }
+
+        if (mesh().isInternalFace(faceI))
+        {
+            label globalNei = globalNumbering().toGlobal(nei[faceI]);
+            if (globalNei != exclude0 && globalNei != exclude1)
+            {
+                globals.insert(globalNei);
+            }
+        }
+        else
+        {
+            label bFaceI = faceI-mesh().nInternalFaces();
+
+            if (isValidBFace[bFaceI])
+            {
+                label globalI = globalNumbering().toGlobal
+                (
+                    mesh().nCells()
+                  + bFaceI
+                );
+
+                if (globalI != exclude0 && globalI != exclude1)
+                {
+                    globals.insert(globalI);
+                }
+            }
+        }
+    }
+}
+
+
+Foam::labelList Foam::faceStencil::calcFaceCells
+(
+    const boolList& isValidBFace,
+    const labelList& faceLabels,
+    labelHashSet& globals
+) const
+{
+    globals.clear();
+
+    insertFaceCells
+    (
+        -1,
+        -1,
+        isValidBFace,
+        faceLabels,
+        globals
+    );
+
+    return globals.toc();
+}
+
+
+// Calculates per face a list of global cell/face indices.
+void Foam::faceStencil::calcFaceStencil
+(
+    const labelListList& globalCellCells,
+    labelListList& faceStencil
+) const
+{
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+    const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+
+    // Determine neighbouring global cell Cells
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelListList neiGlobalCellCells(nBnd);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                neiGlobalCellCells[faceI-mesh_.nInternalFaces()] =
+                    globalCellCells[own[faceI]];
+                faceI++;
+            }
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh_, neiGlobalCellCells, false);
+
+
+
+    // Construct stencil in global numbering
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    faceStencil.setSize(mesh_.nFaces());
+
+    labelHashSet faceStencilSet;
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        faceStencilSet.clear();
+
+        const labelList& ownCCells = globalCellCells[own[faceI]];
+        label globalOwn = ownCCells[0];
+        // Insert cellCells
+        forAll(ownCCells, i)
+        {
+            faceStencilSet.insert(ownCCells[i]);
+        }
+
+        const labelList& neiCCells = globalCellCells[nei[faceI]];
+        label globalNei = neiCCells[0];
+        // Insert cellCells
+        forAll(neiCCells, i)
+        {
+            faceStencilSet.insert(neiCCells[i]);
+        }
+
+        // Guarantee owner first, neighbour second.
+        faceStencil[faceI].setSize(faceStencilSet.size());
+        label n = 0;
+        faceStencil[faceI][n++] = globalOwn;
+        faceStencil[faceI][n++] = globalNei;
+        forAllConstIter(labelHashSet, faceStencilSet, iter)
+        {
+            if (iter.key() != globalOwn && iter.key() != globalNei)
+            {
+                faceStencil[faceI][n++] = iter.key();
+            }
+        }
+        //Pout<< "internalface:" << faceI << " toc:" << faceStencilSet.toc()
+        //    << " faceStencil:" << faceStencil[faceI] << endl;
+    }
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+        label faceI = pp.start();
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                faceStencilSet.clear();
+
+                const labelList& ownCCells = globalCellCells[own[faceI]];
+                label globalOwn = ownCCells[0];
+                forAll(ownCCells, i)
+                {
+                    faceStencilSet.insert(ownCCells[i]);
+                }
+
+                // And the neighbours of the coupled cell
+                const labelList& neiCCells =
+                    neiGlobalCellCells[faceI-mesh_.nInternalFaces()];
+                label globalNei = neiCCells[0];
+                forAll(neiCCells, i)
+                {
+                    faceStencilSet.insert(neiCCells[i]);
+                }
+
+                // Guarantee owner first, neighbour second.
+                faceStencil[faceI].setSize(faceStencilSet.size());
+                label n = 0;
+                faceStencil[faceI][n++] = globalOwn;
+                faceStencil[faceI][n++] = globalNei;
+                forAllConstIter(labelHashSet, faceStencilSet, iter)
+                {
+                    if (iter.key() != globalOwn && iter.key() != globalNei)
+                    {
+                        faceStencil[faceI][n++] = iter.key();
+                    }
+                }
+
+                //Pout<< "coupledface:" << faceI
+                //    << " toc:" << faceStencilSet.toc()
+                //    << " faceStencil:" << faceStencil[faceI] << endl;
+
+                faceI++;
+            }
+        }
+        else if (!isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                faceStencilSet.clear();
+
+                const labelList& ownCCells = globalCellCells[own[faceI]];
+                label globalOwn = ownCCells[0];
+                forAll(ownCCells, i)
+                {
+                    faceStencilSet.insert(ownCCells[i]);
+                }
+
+                // Guarantee owner first
+                faceStencil[faceI].setSize(faceStencilSet.size());
+                label n = 0;
+                faceStencil[faceI][n++] = globalOwn;
+                forAllConstIter(labelHashSet, faceStencilSet, iter)
+                {
+                    if (iter.key() != globalOwn)
+                    {
+                        faceStencil[faceI][n++] = iter.key();
+                    }
+                }
+
+                //Pout<< "boundaryface:" << faceI
+                //    << " toc:" << faceStencilSet.toc()
+                //    << " faceStencil:" << faceStencil[faceI] << endl;
+
+                faceI++;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::faceStencil::faceStencil(const polyMesh& mesh)
+:
+    mesh_(mesh),
+    globalNumbering_(mesh_.nCells()+mesh_.nFaces()-mesh_.nInternalFaces())
+{}
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/extendedStencil/faceStencil/faceStencil.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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::faceStencil
+
+Description
+    baseclass for extended face-cell addressing.
+
+SourceFiles
+    faceStencil.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceStencil_H
+#define faceStencil_H
+
+#include "globalIndex.H"
+#include "boolList.H"
+#include "labelHashSet.H"
+#include "indirectPrimitivePatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class faceStencil Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceStencil
+:
+    public labelListList
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        //- Global numbering for cells and boundary faces
+        const globalIndex globalNumbering_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise assignment
+        void operator=(const faceStencil&);
+
+        
+
+protected:
+
+        //- Merge two lists.
+        static void merge
+        (
+            const label,
+            const label,
+            const labelList&,
+            labelList&
+        );
+
+        //- Merge two lists.
+        static void merge(const label, const labelList&, labelList&);
+
+        //- Valid boundary faces (not empty and not coupled)
+        void validBoundaryFaces(boolList& isValidBFace) const;
+
+        //- Return patch of all coupled faces.
+        autoPtr<indirectPrimitivePatch> allCoupledFacesPatch() const;
+
+        //- Combine operator for labelLists
+        class unionEqOp
+        {
+            public:
+            void operator()( labelList& x, const labelList& y ) const;
+        };
+
+        //- Collect cell neighbours of faces in global numbering
+        void insertFaceCells
+        (
+            const label exclude0,
+            const label exclude1,
+            const boolList& nonEmptyFace,
+            const labelList& faceLabels,
+            labelHashSet& globals
+        ) const;
+
+        //- Collect cell neighbours of faces in global numbering
+        labelList calcFaceCells
+        (
+            const boolList& nonEmptyFace,
+            const labelList& faceLabels,
+            labelHashSet& globals
+        ) const;
+
+
+        //- Collect cell neighbours into extended stencil
+        void calcFaceStencil
+        (
+            const labelListList& globalCellCells,
+            labelListList& faceStencil
+        ) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from mesh
+        explicit faceStencil(const polyMesh&);
+
+
+    // Member Functions
+
+        const polyMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        //- Global numbering for cells and boundary faces
+        const globalIndex& globalNumbering() const
+        {
+            return globalNumbering_;
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C

@@ -1330,7 +1330,7 @@ void Foam::fvMeshSubset::setLargeCellSubset
 
 
     // Add the fvPatches
-    fvMeshSubsetPtr_().addFvPatches(newBoundary);
+    fvMeshSubsetPtr_().addFvPatches(newBoundary, syncPar);
 
     // Subset and add any zones
     subsetZones();

src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.C

@@ -0,0 +1,36 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "linearFit.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeSurfaceInterpolationScheme(linearFit);
+}
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFit.H

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+    linearFit
+
+Description
+    Linear fit interpolation scheme which applies an explicit correction to
+    linear.
+
+SourceFiles
+    linearFit.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linearFit_H
+#define linearFit_H
+
+#include "linear.H"
+#include "linearFitData.H"
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class linearFit Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class linearFit
+:
+    public linear<Type>
+{
+    // Private Data
+        const scalar centralWeight_;
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        linearFit(const linearFit&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const linearFit&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("linearFit");
+
+
+    // Constructors
+
+        //- Construct from mesh and Istream
+        linearFit(const fvMesh& mesh, Istream& is)
+        :
+            linear<Type>(mesh),
+            centralWeight_(readScalar(is))
+        {}
+
+
+        //- Construct from mesh, faceFlux and Istream
+        linearFit
+        (
+            const fvMesh& mesh,
+            const surfaceScalarField& faceFlux,
+            Istream& is
+        )
+        :
+            linear<Type>(mesh),
+            centralWeight_(readScalar(is))
+        {}
+
+
+    // Member Functions
+
+        //- Return true if this scheme uses an explicit correction
+        virtual bool corrected() const
+        {
+            return true;
+        }
+
+        //- Return the explicit correction to the face-interpolate
+        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        correction
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& vf
+        ) const
+        {
+            const fvMesh& mesh = this->mesh();
+
+            const linearFitData& cfd = linearFitData::New
+            (
+                mesh,
+                centralWeight_
+            );
+
+            const extendedStencil& stencil = cfd.stencil();
+            const List<scalarList>& f = cfd.fit();
+
+            return stencil.weightedSum(vf, f);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.C

@@ -0,0 +1,371 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 "linearFitData.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+#include "SVD.H"
+#include "syncTools.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(linearFitData, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+static int count = 0;
+
+Foam::linearFitData::linearFitData
+(
+    const fvMesh& mesh,
+    const scalar cWeight
+)
+:
+    MeshObject<fvMesh, linearFitData>(mesh),
+    centralWeight_(cWeight),
+#   ifdef SPHERICAL_GEOMETRY
+    dim_(2),
+#   else
+    dim_(mesh.nGeometricD()),
+#   endif
+    minSize_
+    (
+        dim_ == 1 ? 2 :
+        dim_ == 2 ? 3 :
+        dim_ == 3 ? 4 : 0
+    ),
+    stencil_(mesh),
+    fit_(mesh.nInternalFaces())
+{
+    if (debug)
+    {
+        Info << "Contructing linearFitData" << endl;
+    }
+
+    // check input
+    if (centralWeight_ < 1 - SMALL)
+    {
+        FatalErrorIn("linearFitData::linearFitData")
+            << "centralWeight requested = " << centralWeight_
+            << " should not be less than one"
+            << exit(FatalError);
+    }
+
+    if (minSize_ == 0)
+    {
+        FatalErrorIn("linearFitSnGradData")
+            << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
+    }
+
+    // store the polynomial size for each cell to write out
+    surfaceScalarField interpPolySize
+    (
+        IOobject
+        (
+            "linearFitInterpPolySize",
+            "constant",
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("linearFitInterpPolySize", dimless, scalar(0))
+    );
+
+    // Get the cell/face centres in stencil order.
+    // Centred face stencils no good for triangles of tets. Need bigger stencils
+    List<List<point> > stencilPoints(stencil_.stencil().size());
+    stencil_.collectData
+    (
+        mesh.C(),
+        stencilPoints
+    );
+
+    // find the fit coefficients for every face in the mesh
+
+    for(label faci = 0; faci < mesh.nInternalFaces(); faci++)
+    {
+        interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
+    }
+
+    Pout<< "count = " << count << endl;
+
+    if (debug)
+    {
+        Info<< "linearFitData::linearFitData() :"
+            << "Finished constructing polynomialFit data"
+            << endl;
+
+        interpPolySize.write();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::linearFitData::findFaceDirs
+(
+    vector& idir,        // value changed in return
+    vector& jdir,        // value changed in return
+    vector& kdir,        // value changed in return
+    const fvMesh& mesh,
+    const label faci
+)
+{
+    idir = mesh.Sf()[faci];
+    //idir = mesh.C()[mesh.neighbour()[faci]] - mesh.C()[mesh.owner()[faci]];
+    idir /= mag(idir);
+
+#   ifndef SPHERICAL_GEOMETRY
+    if (mesh.nGeometricD() <= 2) // find the normal direcion
+    {
+        if (mesh.directions()[0] == -1)
+        {
+            kdir = vector(1, 0, 0);
+        }
+        else if (mesh.directions()[1] == -1)
+        {
+            kdir = vector(0, 1, 0);
+        }
+        else
+        {
+            kdir = vector(0, 0, 1);
+        }
+    }
+    else // 3D so find a direction in the place of the face
+    {
+        const face& f = mesh.faces()[faci];
+        kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
+    }
+#   else
+    // Spherical geometry so kdir is the radial direction
+    kdir = mesh.Cf()[faci];
+#   endif
+
+    if (mesh.nGeometricD() == 3)
+    {
+        // Remove the idir component from kdir and normalise
+        kdir -= (idir & kdir)*idir;
+
+        scalar magk = mag(kdir);
+
+        if (magk < SMALL)
+        {
+            FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
+                << exit(FatalError);
+        }
+        else
+        {
+            kdir /= magk;
+        }
+    }
+
+    jdir = kdir ^ idir;
+}
+
+
+Foam::label Foam::linearFitData::calcFit
+(
+    const List<point>& C,
+    const label faci
+)
+{
+    vector idir(1,0,0);
+    vector jdir(0,1,0);
+    vector kdir(0,0,1);
+    findFaceDirs(idir, jdir, kdir, mesh(), faci);
+
+    scalarList wts(C.size(), scalar(1));
+    wts[0] = centralWeight_;
+    wts[1] = centralWeight_;
+
+    point p0 = mesh().faceCentres()[faci];
+    scalar scale = 0;
+
+    // calculate the matrix of the polynomial components
+    scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
+
+    for(label ip = 0; ip < C.size(); ip++)
+    {
+        const point& p = C[ip];
+
+        scalar px = (p - p0)&idir;
+        scalar py = (p - p0)&jdir;
+#       ifndef SPHERICAL_GEOMETRY
+        scalar pz = (p - p0)&kdir;
+#       else
+        scalar pz = mag(p) - mag(p0);
+#       endif
+
+        if (ip == 0)
+        {
+            scale = max(max(mag(px), mag(py)), mag(pz));
+        }
+
+        px /= scale;
+        py /= scale;
+        pz /= scale;
+
+        label is = 0;
+
+        B[ip][is++] = wts[0]*wts[ip];
+        B[ip][is++] = wts[0]*wts[ip]*px;
+
+        if (dim_ >= 2)
+        {
+            B[ip][is++] = wts[ip]*py;
+        }
+        if (dim_ == 3)
+        {
+            B[ip][is++] = wts[ip]*pz;
+        }
+    }
+
+    // Set the fit
+    label stencilSize = C.size();
+    fit_[faci].setSize(stencilSize);
+    scalarList singVals(minSize_);
+    label nSVDzeros = 0;
+
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& w =
+        mesh().surfaceInterpolation::weights();
+
+    bool goodFit = false;
+    for(int iIt = 0; iIt < 10 && !goodFit; iIt++)
+    {
+        SVD svd(B, SMALL);
+
+        scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0];
+        scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1];
+
+        //goodFit = (fit0 > 0 && fit1 > 0);
+
+        goodFit =
+            (mag(fit0 - w[faci])/w[faci] < 0.5)
+         && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < 0.5);
+
+        //scalar w0Err = fit0/w[faci];
+        //scalar w1Err = fit1/(1 - w[faci]);
+
+        //goodFit =
+        //    (w0Err > 0.5 && w0Err < 1.5)
+        // && (w1Err > 0.5 && w1Err < 1.5);
+
+        if (goodFit)
+        {
+            fit_[faci][0] = fit0;
+            fit_[faci][1] = fit1;
+
+            for(label i=2; i<stencilSize; i++)
+            {
+                fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
+            }
+
+            singVals = svd.S();
+            nSVDzeros = svd.nZeros();
+        }
+        else // (not good fit so increase weight in the centre and for linear)
+        {
+            wts[0] *= 10;
+            wts[1] *= 10;
+
+            for(label i = 0; i < B.n(); i++)
+            {
+                B[i][0] *= 10;
+                B[i][1] *= 10;
+            }
+
+            for(label j = 0; j < B.m(); j++)
+            {
+                B[0][j] *= 10;
+                B[1][j] *= 10;
+            }
+        }
+    }
+
+    // static const scalar L = 0.1;
+    // static const scalar R = 0.2;
+
+    // static const scalar beta = 1.0/(R - L);
+    // static const scalar alpha = R*beta;
+
+    if (goodFit)
+    {
+         if ((mag(fit_[faci][0] - w[faci])/w[faci] < 0.15)
+         && (mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]) < 0.15))
+         {
+             count++;
+             //Pout<< "fit " << mag(fit_[faci][0] - w[faci])/w[faci] << " " << mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]) << endl;
+         }
+
+        // scalar limiter =
+        // max
+        // (
+        //     min
+        //     (
+        //         min(alpha - beta*mag(fit_[faci][0] - w[faci])/w[faci], 1),
+        //         min(alpha - beta*mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]), 1)
+        //     ), 0
+        // );
+
+        // Remove the uncorrected linear coefficients
+        fit_[faci][0] -= w[faci];
+        fit_[faci][1] -= 1 - w[faci];
+
+        // if (limiter < 0.99)
+        // {
+        //     for(label i = 0; i < stencilSize; i++)
+        //     {
+        //         fit_[faci][i] *= limiter;
+        //     }
+        // }
+    }
+    else
+    {
+        Pout<< "Could not fit face " << faci
+            << " " << fit_[faci][0] << " " << w[faci]
+            << " " << fit_[faci][1] << " " << 1 - w[faci]<< endl;
+        fit_[faci] = 0;
+    }
+
+    return minSize_ - nSVDzeros;
+}
+
+
+bool Foam::linearFitData::movePoints()
+{
+    notImplemented("linearFitData::movePoints()");
+
+    return true;
+}
+
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearFit/linearFitData.H

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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
+    linearFitData
+
+Description
+    Data for the linear fit correction interpolation scheme
+
+SourceFiles
+    linearFitData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linearFitData_H
+#define linearFitData_H
+
+#include "MeshObject.H"
+#include "fvMesh.H"
+#include "extendedStencil.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class globalIndex;
+
+/*---------------------------------------------------------------------------*\
+                    Class linearFitData Declaration
+\*---------------------------------------------------------------------------*/
+
+class linearFitData
+:
+    public MeshObject<fvMesh, linearFitData>
+{
+    // Private data
+
+        //- weights for central stencil
+        const scalar centralWeight_;
+
+        //- dimensionality of the geometry
+        const label dim_;
+
+        //- minimum stencil size
+        const label minSize_;
+
+        //- Extended stencil addressing
+        extendedStencil stencil_;
+
+        //- For each cell in the mesh store the values which multiply the
+        //  values of the stencil to obtain the gradient for each direction
+        List<scalarList> fit_;
+
+
+    // Private member functions
+
+        //- Find the normal direction and i, j and k directions for face faci
+        static void findFaceDirs
+        (
+            vector& idir,        // value changed in return
+            vector& jdir,        // value changed in return
+            vector& kdir,        // value changed in return
+            const fvMesh& mesh,
+            const label faci
+        );
+
+        label calcFit(const List<point>&, const label faci);
+
+
+public:
+
+    TypeName("linearFitData");
+
+
+    // Constructors
+
+        explicit linearFitData
+        (
+            const fvMesh& mesh,
+            scalar cWeightDim
+        );
+
+
+    // Destructor
+
+        virtual ~linearFitData()
+        {}
+
+
+    // Member functions
+
+
+        //- Return reference to the stencil
+        const extendedStencil& stencil() const
+        {
+            return stencil_;
+        }
+
+        //- Return reference to fit coefficients
+        const List<scalarList>& fit() const
+        {
+            return fit_;
+        }
+
+        //- Delete the data when the mesh moves not implemented
+        virtual bool movePoints();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFit.H

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    quadraticFit
+    Foam::quadraticFit
 
 Description
     Quadratic fit interpolation scheme which applies an explicit correction to
@@ -37,9 +37,9 @@ SourceFiles
 #ifndef quadraticFit_H
 #define quadraticFit_H
 
-#include "linear.H"
 #include "quadraticFitData.H"
-#include "extendedStencil.H"
+#include "linear.H"
+#include "centredCFCStencilObject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,8 +56,16 @@ class quadraticFit
     public linear<Type>
 {
     // Private Data
+
+        //- Factor the fit is allowed to deviate from linear.
+        //  This limits the amount of high-order correction and increases
+        //  stability on bad meshes
+        const scalar linearLimitFactor_;
+
+        //- Weights for central stencil
         const scalar centralWeight_;
 
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -79,6 +87,7 @@ public:
         quadraticFit(const fvMesh& mesh, Istream& is)
         :
             linear<Type>(mesh),
+            linearLimitFactor_(readScalar(is)),
             centralWeight_(readScalar(is))
         {}
 
@@ -92,6 +101,7 @@ public:
         )
         :
             linear<Type>(mesh),
+            linearLimitFactor_(readScalar(is)),
             centralWeight_(readScalar(is))
         {}
 
@@ -113,13 +123,20 @@ public:
         {
             const fvMesh& mesh = this->mesh();
 
+            const extendedCentredStencil& stencil =
+            centredCFCStencilObject::New
+            (
+                mesh
+            );
+
             const quadraticFitData& cfd = quadraticFitData::New
             (
                 mesh,
+                stencil,
+                linearLimitFactor_,
                 centralWeight_
             );
 
-            const extendedStencil& stencil = cfd.stencil();
             const List<scalarList>& f = cfd.fit();
 
             return stencil.weightedSum(vf, f);

src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C

@@ -29,7 +29,7 @@ License
 #include "volFields.H"
 #include "SVD.H"
 #include "syncTools.H"
-
+#include "extendedCentredStencil.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -41,15 +41,16 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
 
-static int count = 0;
-
 Foam::quadraticFitData::quadraticFitData
 (
     const fvMesh& mesh,
+    const extendedCentredStencil& stencil,
+    const scalar linearLimitFactor,
     const scalar cWeight
 )
 :
     MeshObject<fvMesh, quadraticFitData>(mesh),
+    linearLimitFactor_(linearLimitFactor),
     centralWeight_(cWeight),
 #   ifdef SPHERICAL_GEOMETRY
     dim_(2),
@@ -60,9 +61,8 @@ Foam::quadraticFitData::quadraticFitData
     (
         dim_ == 1 ? 3 :
         dim_ == 2 ? 6 :
-        dim_ == 3 ? 9 : 0
+        dim_ == 3 ? 7 : 0
     ),
-    stencil_(mesh),
     fit_(mesh.nInternalFaces())
 {
     if (debug)
@@ -102,8 +102,8 @@ Foam::quadraticFitData::quadraticFitData
 
     // Get the cell/face centres in stencil order.
     // Centred face stencils no good for triangles of tets. Need bigger stencils
-    List<List<point> > stencilPoints(stencil_.stencil().size());
-    stencil_.collectData
+    List<List<point> > stencilPoints(mesh.nFaces());
+    stencil.collectData
     (
         mesh.C(),
         stencilPoints
@@ -116,8 +116,6 @@ Foam::quadraticFitData::quadraticFitData
         interpPolySize[faci] = calcFit(stencilPoints[faci], faci);
     }
 
-    Pout<< "count = " << count << endl;
-
     if (debug)
     {
         Info<< "quadraticFitData::quadraticFitData() :"
@@ -243,14 +241,14 @@ Foam::label Foam::quadraticFitData::calcFit
         {
             B[ip][is++] = wts[ip]*py;
             B[ip][is++] = wts[ip]*px*py;
-            B[ip][is++] = wts[ip]*sqr(py);
+            //B[ip][is++] = wts[ip]*sqr(py);
         }
         if (dim_ == 3)
         {
             B[ip][is++] = wts[ip]*pz;
             B[ip][is++] = wts[ip]*px*pz;
             //B[ip][is++] = wts[ip]*py*pz;
-            B[ip][is++] = wts[ip]*sqr(pz);
+            //B[ip][is++] = wts[ip]*sqr(pz);
         }
     }
 
@@ -274,15 +272,15 @@ Foam::label Foam::quadraticFitData::calcFit
         //goodFit = (fit0 > 0 && fit1 > 0);
 
         goodFit =
-            (mag(fit0 - w[faci])/w[faci] < 0.15)
-         && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < 0.15);
+            (mag(fit0 - w[faci])/w[faci] < linearLimitFactor_)
+         && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < linearLimitFactor_);
 
         //scalar w0Err = fit0/w[faci];
         //scalar w1Err = fit1/(1 - w[faci]);
 
         //goodFit =
-        //    (w0Err > 0.5 && w0Err < 1.5)
-        // && (w1Err > 0.5 && w1Err < 1.5);
+        //    (w0Err > linearLimitFactor_ && w0Err < (1 + linearLimitFactor_))
+        // && (w1Err > linearLimitFactor_ && w1Err < (1 + linearLimitFactor_));
 
         if (goodFit)
         {
@@ -324,13 +322,6 @@ Foam::label Foam::quadraticFitData::calcFit
 
     if (goodFit)
     {
-         if ((mag(fit_[faci][0] - w[faci])/w[faci] < 0.15)
-         && (mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]) < 0.15))
-         {
-             count++;
-             //Pout<< "fit " << mag(fit_[faci][0] - w[faci])/w[faci] << " " << mag(fit_[faci][1] - (1 - w[faci]))/(1 - w[faci]) << endl;
-         }
-
         // scalar limiter =
         // max
         // (
@@ -358,6 +349,7 @@ Foam::label Foam::quadraticFitData::calcFit
         Pout<< "Could not fit face " << faci
             << " " << fit_[faci][0] << " " << w[faci]
             << " " << fit_[faci][1] << " " << 1 - w[faci]<< endl;
+
         fit_[faci] = 0;
     }
 

src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.H

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    quadraticFitData
+    Foam::quadraticFitData
 
 Description
     Data for the quadratic fit correction interpolation scheme
@@ -38,14 +38,13 @@ SourceFiles
 
 #include "MeshObject.H"
 #include "fvMesh.H"
-#include "extendedStencil.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class globalIndex;
+class extendedCentredStencil;
 
 /*---------------------------------------------------------------------------*\
                     Class quadraticFitData Declaration
@@ -57,18 +56,20 @@ class quadraticFitData
 {
     // Private data
 
-        //- weights for central stencil
+        //- Factor the fit is allowed to deviate from linear.
+        //  This limits the amount of high-order correction and increases
+        //  stability on bad meshes
+        const scalar linearLimitFactor_;
+
+        //- Weights for central stencil
         const scalar centralWeight_;
 
-        //- dimensionality of the geometry
+        //- Dimensionality of the geometry
         const label dim_;
 
-        //- minimum stencil size
+        //- Minimum stencil size
         const label minSize_;
 
-        //- Extended stencil addressing
-        extendedStencil stencil_;
-
         //- For each cell in the mesh store the values which multiply the
         //  values of the stencil to obtain the gradient for each direction
         List<scalarList> fit_;
@@ -99,25 +100,19 @@ public:
         explicit quadraticFitData
         (
             const fvMesh& mesh,
-            scalar cWeightDim
+            const extendedCentredStencil& stencil,
+            const scalar linearLimitFactor,
+            const scalar centralWeight
         );
 
 
-    // Destructor
-
-        virtual ~quadraticFitData()
-        {}
+    //- Destructor
+    virtual ~quadraticFitData()
+    {}
 
 
     // Member functions
 
-
-        //- Return reference to the stencil
-        const extendedStencil& stencil() const
-        {
-            return stencil_;
-        }
-
         //- Return reference to fit coefficients
         const List<scalarList>& fit() const
         {

src/meshTools/searchableSurface/searchableSurfaceWithGaps.C

@@ -232,9 +232,9 @@ void Foam::searchableSurfaceWithGaps::findLine
         // do we accept the hit.
 
         List<pointIndexHit> plusInfo;
-        surface().findLine(start+offset0, end+offset0, plusInfo);
+        surface().findLine(compactStart+offset0, compactEnd+offset0, plusInfo);
         List<pointIndexHit> minInfo;
-        surface().findLine(start-offset0, end-offset0, minInfo);
+        surface().findLine(compactStart-offset0, compactEnd-offset0, minInfo);
 
         // Extract any hits
         forAll(plusInfo, i)
@@ -257,14 +257,28 @@ void Foam::searchableSurfaceWithGaps::findLine
                 compactStart[i] = compactStart[mapI];
                 compactEnd[i] = compactEnd[mapI];
                 compactMap[i] = compactMap[mapI];
+                offset0[i] = offset0[mapI];
+                offset1[i] = offset1[mapI];
             }
             compactStart.setSize(plusMissMap.size());
             compactEnd.setSize(plusMissMap.size());
             compactMap.setSize(plusMissMap.size());
+            offset0.setSize(plusMissMap.size());
+            offset1.setSize(plusMissMap.size());
 
 
-            surface().findLine(start+offset1, end+offset1, plusInfo);
-            surface().findLine(start-offset1, end-offset1, minInfo);
+            surface().findLine
+            (
+                compactStart+offset1,
+                compactEnd+offset1,
+                plusInfo
+            );
+            surface().findLine
+            (
+                compactStart-offset1,
+                compactEnd-offset1,
+                minInfo
+            );
 
             // Extract any hits
             forAll(plusInfo, i)

src/meshTools/searchableSurface/searchableSurfaceWithGaps.H

@@ -148,6 +148,8 @@ public:
 
         // Multiple point queries.
 
+            //- Find nearest on original surface. Note:does not use perturbation
+            //  and hence might be inconsistent with intersections.
             virtual void findNearest
             (
                 const pointField& sample,

src/sampling/probes/probes.C

@@ -240,7 +240,7 @@ bool Foam::probes::checkFieldTypes()
 
                     forAll(probeLocations_, probeI)
                     {
-                        *sPtr<< setw(w) << probeLocations_[probeI][cmpt];
+                        *sPtr<< ' ' << setw(w) << probeLocations_[probeI][cmpt];
                     }
                     *sPtr << endl;
                 }

src/sampling/probes/probesTemplates.C

@@ -109,7 +109,7 @@ void Foam::probes::sampleAndWrite
 
         forAll(values, probeI)
         {
-            probeStream << setw(w) << values[probeI];
+            probeStream << ' ' << setw(w) << values[probeI];
         }
         probeStream << endl;
     }

src/triSurface/triSurface/interfaces/STL/readSTLASCII.L

@@ -459,6 +459,7 @@ bool triSurface::readSTLASCII(const fileName& STLfileName)
         patches_[iter()].name() = iter.key();
     }
 
+    // Fill in the missing information in the patches
     setDefaultPatches();
 
     return true;

src/triSurface/triSurface/triSurface.C

@@ -631,15 +631,16 @@ surfacePatchList triSurface::calcPatches(labelList& faceMap) const
 
     faceMap = sortedRegion.indices();
 
-    // Compact regions
-
-    // Get last region
-    //label maxRegion = 0;                  // for compacted regions
+    // Extend regions
     label maxRegion = patches_.size()-1;    // for non-compacted regions
 
     if (faceMap.size() > 0)
     {
-        maxRegion = operator[](faceMap[faceMap.size() - 1]).region();
+        maxRegion = max
+        (
+            maxRegion,
+            operator[](faceMap[faceMap.size() - 1]).region()
+        );
     }
 
     // Get new region list

wmake/wmakeScheduler

@@ -163,18 +163,6 @@ colourPipe()
 }
 
 
-#
-# prefix message with [HOSTNAME]
-#
-prefixPipe()
-{
-   while read line
-   do
-      echo "[$@] $line"
-   done
-}
-
-
 colourIndex=0
 
 while :
@@ -210,11 +198,11 @@ do
                retval=$?
             else
                if [ "$host" = "$HOST" ]; then
-                  eval $* 2>&1 | prefixPipe "$host"
+                  eval $* 2>&1
                elif [ -n "$JOB_ID" ]; then
-                  qrsh -inherit -v PWD $host "$rcmd" | prefixPipe "$host"
+                  qrsh -inherit -v PWD $host "$rcmd"
                else
-                  ssh $host "$sourceFoam 2>/dev/null; cd $PWD && $rcmd" 2>&1 | prefixPipe "$host"
+                  ssh $host "$sourceFoam 2>/dev/null; cd $PWD && $rcmd" 2>&1
                fi
                retval=$?
             fi