From 12910bf8ba7914db916389e955860d1ddaf47101 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 18 Oct 2010 18:03:57 +0100
Subject: [PATCH 01/17] LienLeschzinerLowRe: Corrected name of nut

---
 .../RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C               | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
index 486eccaa98..4a76228d23 100644
--- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
+++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
@@ -172,7 +172,7 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
     (
         IOobject
         (
-            "epsilon",
+            "nut",
             runTime_.timeName(),
             mesh_,
             IOobject::NO_READ,

From 82eedb5b3cf3b761098203d753e6e0abe2955c85 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 19 Oct 2010 09:55:33 +0100
Subject: [PATCH 02/17] BUG: Time : set up file monitoring only if
 runTimeModifiable

---
 src/OpenFOAM/db/Time/Time.C | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C
index 487576586b..ff1247dcc7 100644
--- a/src/OpenFOAM/db/Time/Time.C
+++ b/src/OpenFOAM/db/Time/Time.C
@@ -245,13 +245,13 @@ Foam::Time::Time
     readLibs_(controlDict_, "libs"),
     functionObjects_(*this)
 {
+    setControls();
+
     // Time objects not registered so do like objectRegistry::checkIn ourselves.
     if (runTimeModifiable_)
     {
         controlDict_.watchIndex() = addWatch(controlDict_.filePath());
     }
-
-    setControls();
 }
 
 
@@ -307,14 +307,18 @@ Foam::Time::Time
     readLibs_(controlDict_, "libs"),
     functionObjects_(*this)
 {
+    setControls();
 
     // Time objects not registered so do like objectRegistry::checkIn ourselves.
     if (runTimeModifiable_)
     {
-        controlDict_.watchIndex() = addWatch(controlDict_.filePath());
+        // File might not exist yet.
+        fileName f(controlDict_.filePath());
+        if (f != fileName::null)
+        {
+            controlDict_.watchIndex() = addWatch(f);
+        }
     }
-
-    setControls();
 }
 
 

From a5ee1669b4d6e2abd0758fbb2382574cb7948751 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 19 Oct 2010 09:56:29 +0100
Subject: [PATCH 03/17] ENH: fileMonitor : continue if inotify_init failed

---
 src/OSspecific/POSIX/fileMonitor.C | 62 ++++++++++++++++++++++--------
 1 file changed, 45 insertions(+), 17 deletions(-)

diff --git a/src/OSspecific/POSIX/fileMonitor.C b/src/OSspecific/POSIX/fileMonitor.C
index 982d5e28dc..a6befd5908 100644
--- a/src/OSspecific/POSIX/fileMonitor.C
+++ b/src/OSspecific/POSIX/fileMonitor.C
@@ -39,7 +39,7 @@ Class
 #else
 #   include <sys/inotify.h>
 #   include <sys/ioctl.h>
-
+#   include <errno.h>
 #   define EVENT_SIZE  ( sizeof (struct inotify_event) )
 #   define EVENT_LEN   (EVENT_SIZE + 16)
 #   define EVENT_BUF_LEN     ( 1024 * EVENT_LEN )
@@ -144,7 +144,7 @@ namespace Foam
 
 #else
         //- File descriptor for the inotify instance
-        int fd;
+        int inotifyFd_;
 
         //- Current watchIDs and corresponding directory id
         DynamicList<label> dirWatches_;
@@ -153,23 +153,39 @@ namespace Foam
         //- initialise inotify
         inline fileMonitorWatcher(const label sz = 20)
         :
-            fd(inotify_init()),
+            inotifyFd_(inotify_init()),
             dirWatches_(sz),
             dirFiles_(sz)
-        {}
+        {
+            if (inotifyFd_ < 0)
+            {
+                WarningIn("fileMonitorWatcher(const label)")
+                    << "Failed allocating an inotify descriptor : "
+                    << string(strerror(errno)) << endl
+                    << "    Please increase the number of allowable "
+                    << "inotify instances" << endl
+                    << "    (/proc/sys/fs/inotify/max_user_instances on Linux)"
+                    << endl
+                    << "    Continuing without additional file monitoring."
+                    << endl;
+            }
+        }
 
         //- remove all watches
         inline ~fileMonitorWatcher()
         {
-            forAll(dirWatches_, i)
+            if (inotifyFd_ >= 0)
             {
-                if (dirWatches_[i] >= 0)
+                forAll(dirWatches_, i)
                 {
-                    if (inotify_rm_watch(fd, int(dirWatches_[i])))
+                    if (dirWatches_[i] >= 0)
                     {
-                        WarningIn("fileMonitor::~fileMonitor()")
-                            << "Failed deleting directory watch "
-                            << dirWatches_[i] << endl;
+                        if (inotify_rm_watch(inotifyFd_, int(dirWatches_[i])))
+                        {
+                            WarningIn("fileMonitor::~fileMonitor()")
+                                << "Failed deleting directory watch "
+                                << dirWatches_[i] << endl;
+                        }
                     }
                 }
             }
@@ -177,10 +193,15 @@ namespace Foam
 
         inline bool addWatch(const label watchFd, const fileName& fName)
         {
+            if (inotifyFd_ < 0)
+            {
+                return false;
+            }
+
             // Add/retrieve watch on directory containing file
             label dirWatchID = inotify_add_watch
             (
-                fd,
+                inotifyFd_,
                 fName.path().c_str(),
                 IN_CLOSE_WRITE
             );
@@ -189,7 +210,8 @@ namespace Foam
             {
                 FatalErrorIn("addWatch(const label, const fileName&)")
                     << "Failed adding watch " << watchFd
-                    << " to directory " << fName
+                    << " to directory " << fName << " due to "
+                    << string(strerror(errno))
                     << exit(FatalError);
             }
 
@@ -209,6 +231,11 @@ namespace Foam
 
         inline bool removeWatch(const label watchFd)
         {
+            if (inotifyFd_ < 0)
+            {
+                return false;
+            }
+
             dirWatches_[watchFd] = -1;
             return true;
         }
@@ -263,11 +290,11 @@ void Foam::fileMonitor::checkFiles() const
         fd_set fdSet;
         // Add notify descriptor to select fd_set
         FD_ZERO(&fdSet);
-        FD_SET(watcher_->fd, &fdSet);
+        FD_SET(watcher_->inotifyFd_, &fdSet);
 
         int ready = select
         (
-            watcher_->fd+1,     // num filedescriptors in fdSet
+            watcher_->inotifyFd_+1,     // num filedescriptors in fdSet
             &fdSet,             // fdSet with only inotifyFd
             NULL,               // No writefds
             NULL,               // No errorfds
@@ -280,15 +307,15 @@ void Foam::fileMonitor::checkFiles() const
                 << "Problem in issuing select."
                 << abort(FatalError);
         }
-        else if (FD_ISSET(watcher_->fd, &fdSet))
+        else if (FD_ISSET(watcher_->inotifyFd_, &fdSet))
         {
             // Read events
-            ssize_t nBytes = read(watcher_->fd, buffer, EVENT_BUF_LEN);
+            ssize_t nBytes = read(watcher_->inotifyFd_, buffer, EVENT_BUF_LEN);
 
             if (nBytes < 0)
             {
                 FatalErrorIn("fileMonitor::updateStates(const fileName&)")
-                    << "read of " << watcher_->fd
+                    << "read of " << watcher_->inotifyFd_
                     << " failed with " << label(nBytes)
                     << abort(FatalError);
             }
@@ -364,6 +391,7 @@ Foam::label Foam::fileMonitor::addWatch(const fileName& fName)
     label watchFd;
 
     label sz = freeWatchFds_.size();
+
     if (sz)
     {
         watchFd = freeWatchFds_[sz-1];

From 51b8a7367c03d79b276a55e4b02f616b413c9db2 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 19 Oct 2010 17:40:36 +0100
Subject: [PATCH 04/17] ENH: MeshWave: allow 1 cells cases.

---
 .../utilities/mesh/advanced/selectCells/selectCells.C         | 2 +-
 src/dynamicMesh/meshCut/directions/directions.C               | 3 ++-
 src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C       | 4 ++--
 .../inverseFaceDistance/inverseFaceDistanceDiffusivity.C      | 2 +-
 src/meshTools/cellClassification/cellClassification.C         | 2 +-
 src/meshTools/cellDist/patchWave/patchDataWave.C              | 2 +-
 src/meshTools/cellDist/patchWave/patchWave.C                  | 2 +-
 .../decompositionMethods/structuredDecomp/structuredDecomp.C  | 2 +-
 src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C  | 2 +-
 9 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/applications/utilities/mesh/advanced/selectCells/selectCells.C b/applications/utilities/mesh/advanced/selectCells/selectCells.C
index a298c4207d..feab5ff4c7 100644
--- a/applications/utilities/mesh/advanced/selectCells/selectCells.C
+++ b/applications/utilities/mesh/advanced/selectCells/selectCells.C
@@ -298,7 +298,7 @@ label selectOutsideCells
         mesh,
         outsideFaces.shrink(),
         outsideFacesInfo.shrink(),
-        mesh.nCells()  // max iterations
+        mesh.globalData().nTotalCells()+1  // max iterations
     );
 
     // Now regionCalc should hold info on cells that are reachable from
diff --git a/src/dynamicMesh/meshCut/directions/directions.C b/src/dynamicMesh/meshCut/directions/directions.C
index 5f5fb220af..7024c2ffbc 100644
--- a/src/dynamicMesh/meshCut/directions/directions.C
+++ b/src/dynamicMesh/meshCut/directions/directions.C
@@ -32,6 +32,7 @@ License
 #include "meshTools.H"
 #include "hexMatcher.H"
 #include "Switch.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -206,7 +207,7 @@ Foam::vectorField Foam::directions::propagateDirection
         mesh,
         changedFaces,
         changedFacesInfo,
-        mesh.nCells()
+        mesh.globalData().nTotalCells()+1
     );
 
     const List<directionInfo>& cellInfo = directionCalc.allCellInfo();
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
index 9482861bde..c9919351fc 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C
@@ -2332,7 +2332,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
         seedFacesInfo.clear();
 
         // Iterate until no change. Now 2:1 face difference should be satisfied
-        levelCalc.iterate(mesh_.globalData().nTotalFaces());
+        levelCalc.iterate(mesh_.globalData().nTotalFaces()+1);
 
 
         // Now check point-connected cells (face-connected cells already ok):
@@ -2836,7 +2836,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2
         seedFacesInfo,
         allFaceInfo,
         allCellInfo,
-        mesh_.globalData().nTotalCells()
+        mesh_.globalData().nTotalCells()+1
     );
 
 
diff --git a/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C b/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C
index 72531d5cfd..02db9ad941 100644
--- a/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C
+++ b/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C
@@ -115,7 +115,7 @@ void Foam::inverseFaceDistanceDiffusivity::correct()
         mesh,
         changedFaces,
         faceDist,
-        mesh.globalData().nTotalCells() // max iterations
+        mesh.globalData().nTotalCells()+1 // max iterations
     );
 
     const List<wallPoint>& faceInfo = waveInfo.allFaceInfo();
diff --git a/src/meshTools/cellClassification/cellClassification.C b/src/meshTools/cellClassification/cellClassification.C
index a723665477..daf1de6ac0 100644
--- a/src/meshTools/cellClassification/cellClassification.C
+++ b/src/meshTools/cellClassification/cellClassification.C
@@ -322,7 +322,7 @@ void Foam::cellClassification::markCells
         changedFaces,                       // Labels of changed faces
         changedFacesInfo,                   // Information on changed faces
         cellInfoList,                       // Information on all cells
-        mesh_.globalData().nTotalCells()    // max iterations
+        mesh_.globalData().nTotalCells()+1  // max iterations
     );
 
     // Get information out of cellInfoList
diff --git a/src/meshTools/cellDist/patchWave/patchDataWave.C b/src/meshTools/cellDist/patchWave/patchDataWave.C
index bb6da1750d..988b8a0be5 100644
--- a/src/meshTools/cellDist/patchWave/patchDataWave.C
+++ b/src/meshTools/cellDist/patchWave/patchDataWave.C
@@ -225,7 +225,7 @@ void Foam::patchDataWave<TransferType>::correct()
         mesh(),
         changedFaces,
         faceDist,
-        mesh().globalData().nTotalCells() // max iterations
+        mesh().globalData().nTotalCells()+1 // max iterations
     );
 
 
diff --git a/src/meshTools/cellDist/patchWave/patchWave.C b/src/meshTools/cellDist/patchWave/patchWave.C
index 85118da58d..826430b003 100644
--- a/src/meshTools/cellDist/patchWave/patchWave.C
+++ b/src/meshTools/cellDist/patchWave/patchWave.C
@@ -186,7 +186,7 @@ void Foam::patchWave::correct()
         mesh(),
         changedFaces,
         faceDist,
-        mesh().globalData().nTotalCells()   // max iterations
+        mesh().globalData().nTotalCells()+1 // max iterations
     );
 
 
diff --git a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
index 32f86c82c8..04338dc46b 100644
--- a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C
@@ -148,7 +148,7 @@ Foam::labelList Foam::structuredDecomp::decompose
         patchData,
         faceData,
         cellData,
-        mesh.globalData().nTotalCells()
+        mesh.globalData().nTotalCells()+1
     );
 
     // And extract
diff --git a/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C b/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C
index af250a25f7..d97395d864 100644
--- a/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C
+++ b/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C
@@ -137,7 +137,7 @@ void smoothDelta::calcDelta()
         changedFacesInfo,
         faceDeltaData,
         cellDeltaData,
-        mesh_.globalData().nTotalCells()    // max iterations
+        mesh_.globalData().nTotalCells()+1   // max iterations
     );
 
     forAll(delta_, cellI)

From 17c42465b1ad9ffc0f3017fb72e5e1b431af8f1d Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 19 Oct 2010 17:45:51 +0100
Subject: [PATCH 05/17] ENH: selectCells : add one offset to number of cells

---
 applications/utilities/mesh/advanced/selectCells/selectCells.C | 1 +
 1 file changed, 1 insertion(+)

diff --git a/applications/utilities/mesh/advanced/selectCells/selectCells.C b/applications/utilities/mesh/advanced/selectCells/selectCells.C
index feab5ff4c7..827449c02c 100644
--- a/applications/utilities/mesh/advanced/selectCells/selectCells.C
+++ b/applications/utilities/mesh/advanced/selectCells/selectCells.C
@@ -55,6 +55,7 @@ Description
 #include "edgeStats.H"
 #include "treeDataTriSurface.H"
 #include "indexedOctree.H"
+#include "globalMeshData.H"
 
 using namespace Foam;
 

From 0d470d124e88f1cebb76e110241abba358a2c690 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 19 Oct 2010 17:46:21 +0100
Subject: [PATCH 06/17] ENH: foamConfigurePaths : check that sed has done
 something

---
 bin/tools/foamConfigurePaths | 72 ++++++++++++++++++++++++------------
 1 file changed, 48 insertions(+), 24 deletions(-)

diff --git a/bin/tools/foamConfigurePaths b/bin/tools/foamConfigurePaths
index 823c816d51..b820bece64 100755
--- a/bin/tools/foamConfigurePaths
+++ b/bin/tools/foamConfigurePaths
@@ -34,10 +34,10 @@ usage() {
     cat<<USAGE
 
 usage: ${0##*/}
-  --foamInstall dir       specify installation directory (e.g. /opt)
-  --projectName name      specify project name (e.g. openfoam170)
-  --archOption  arch      specify architecture option (only 32 or 64 applicable)
-  --paraviewInstall dir   specify ParaView_DIR (e.g. /opt/paraviewopenfoam380)
+  --foamInstall dir   		specify installation directory (e.g. /opt)
+  --projectName name 		specify project name (e.g. openfoam170)
+  --archOption  arch		specify architecture option (only 32 or 64 applicable)
+  --paraviewInstall dir 	specify ParaView_DIR (e.g. /opt/paraviewopenfoam380)
 
 * hardcode paths to installation
 
@@ -46,6 +46,17 @@ USAGE
 }
 
 
+# Function to do replacement on file. Checks if any replacement has been done.
+# inlineSed <file> <sedCommand> <description>
+_inlineSed()
+{
+    backup=`tempfile`
+    cp $1 $backup
+    sed -i -e "$2" $1
+    cmp $1 $backup || usage "Failed : $3"
+}
+
+
 [ -f etc/bashrc -a -f etc/settings.sh ] || usage "Please run from top-level directory of installation"
 
 unset foamInstall projectName archOption paraviewInstall
@@ -59,30 +70,40 @@ do
         ;;
     --foamInstall)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        foamInstall="$2"
-        echo "Replacing foamInstall setting by $foamInstall"
-        sed -i -e '/^[^#]/s@foamInstall=.*@foamInstall='"$foamInstall@" etc/bashrc
+	foamInstall="$2"
+        echo "** foamInstall:$foamInstall"
+
+        _inlineSed \
+            etc/bashrc \
+            '/^[^#]/s@foamInstall=.*@foamInstall='"$foamInstall@" \
+            "Replacing foamInstall setting by $foamInstall"
         shift 2
         ;;
    --projectName)
-        [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        projectName="$2"
-        echo "Replacing WM_PROJECT_DIR setting by $projectName"
-        sed -i -e '/^[^#]/s@WM_PROJECT_DIR=.*@WM_PROJECT_DIR=$WM_PROJECT_INST_DIR/'"$projectName@" etc/bashrc
-        shift 2
-        ;;
+       	[ "$#" -ge 2 ] || usage "'$1' option requires an argument"
+	projectName="$2"
+        _inlineSed \
+            etc/bashrc \
+            '/^[^#]/s@WM_PROJECT_DIR=.*@WM_PROJECT_DIR=$WM_PROJECT_INST_DIR/'"$projectName@" \
+            "Replacing WM_PROJECT_DIR setting by $projectName"
+       	shift 2
+       	;;
     --archOption)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        archOption="$2"
-        echo "Replacing WM_ARCH_OPTION setting by $archOption"
-        sed -i -e '/^[^#]/s@: ${WM_ARCH_OPTION:=64}@WM_ARCH_OPTION='"$archOption@" etc/bashrc
+	archOption="$2"
+        _inlineSed \
+            etc/bashrc \
+            '/^[^#]/s@: ${WM_ARCH_OPTION:=64}@WM_ARCH_OPTION='"$archOption@" \
+             "Replacing WM_ARCH_OPTION setting by $archOption"
         shift 2
         ;;
     --paraviewInstall)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-        paraviewInstall="$2"
-        echo "Replacing ParaView_DIR setting by $paraviewInstall"
-        sed -i -e '/^[^#]/s@ParaView_DIR=.*@ParaView_DIR='"$paraviewInstall@" etc/apps/paraview3/bashrc
+	paraviewInstall="$2"
+        _inlineSed \
+            etc/apps/paraview3/bashrc \
+            '/^[^#]/s@ParaView_DIR=.*@ParaView_DIR='"$paraviewInstall@" \
+             "Replacing ParaView_DIR setting by $paraviewInstall"
         shift 2
         ;;
     *)
@@ -97,11 +118,14 @@ done
 #sed -i -e 's@WM_PROJECT=.*@WM_PROJECT='"$projectName@" etc/bashrc
 
 # Replace the WM_MPLIB always
-echo "Replacing WM_MPLIB setting by SYSTEMOPENMPI"
-sed -i -e '/^[^#]/s@: ${WM_MPLIB:=.*}@WM_MPLIB=SYSTEMOPENMPI@' etc/bashrc
-
+_inlineSed \
+    etc/bashrc \
+    '/^[^#]/s@: ${WM_MPLIB:=.*}@WM_MPLIB=SYSTEMOPENMPI@' \
+    "Replacing WM_MPLIB setting by SYSTEMOPENMPI"
 # Replace the compilerInstall always
-echo "Replacing compilerInstall setting by system"
-sed -i -e '/^[^#]/s@: ${compilerInstall:=.*}@compilerInstall=system@' etc/settings.sh
+_inlineSed \
+    etc/settings.sh \
+    '/^[^#]/s@: ${compilerInstall:=.*}@compilerInstall=system@' \
+    "Replacing compilerInstall setting by system"
 
 #------------------------------------------------------------------------------

From 497daa30611baf45a4eb9d92960d72acbec9613c Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Tue, 19 Oct 2010 17:47:33 +0100
Subject: [PATCH 07/17] STYLE: foamConfigurePaths : expanded tab

---
 bin/tools/foamConfigurePaths | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)
 mode change 100755 => 100644 bin/tools/foamConfigurePaths

diff --git a/bin/tools/foamConfigurePaths b/bin/tools/foamConfigurePaths
old mode 100755
new mode 100644
index b820bece64..eb85ed9d97
--- a/bin/tools/foamConfigurePaths
+++ b/bin/tools/foamConfigurePaths
@@ -34,10 +34,10 @@ usage() {
     cat<<USAGE
 
 usage: ${0##*/}
-  --foamInstall dir   		specify installation directory (e.g. /opt)
-  --projectName name 		specify project name (e.g. openfoam170)
-  --archOption  arch		specify architecture option (only 32 or 64 applicable)
-  --paraviewInstall dir 	specify ParaView_DIR (e.g. /opt/paraviewopenfoam380)
+  --foamInstall dir         specify installation directory (e.g. /opt)
+  --projectName name        specify project name (e.g. openfoam170)
+  --archOption  arch        specify architecture option (only 32 or 64 applicable)
+  --paraviewInstall dir     specify ParaView_DIR (e.g. /opt/paraviewopenfoam380)
 
 * hardcode paths to installation
 
@@ -70,7 +70,7 @@ do
         ;;
     --foamInstall)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-	foamInstall="$2"
+    foamInstall="$2"
         echo "** foamInstall:$foamInstall"
 
         _inlineSed \
@@ -80,17 +80,17 @@ do
         shift 2
         ;;
    --projectName)
-       	[ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-	projectName="$2"
+        [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
+    projectName="$2"
         _inlineSed \
             etc/bashrc \
             '/^[^#]/s@WM_PROJECT_DIR=.*@WM_PROJECT_DIR=$WM_PROJECT_INST_DIR/'"$projectName@" \
             "Replacing WM_PROJECT_DIR setting by $projectName"
-       	shift 2
-       	;;
+        shift 2
+        ;;
     --archOption)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-	archOption="$2"
+    archOption="$2"
         _inlineSed \
             etc/bashrc \
             '/^[^#]/s@: ${WM_ARCH_OPTION:=64}@WM_ARCH_OPTION='"$archOption@" \
@@ -99,7 +99,7 @@ do
         ;;
     --paraviewInstall)
         [ "$#" -ge 2 ] || usage "'$1' option requires an argument"
-	paraviewInstall="$2"
+    paraviewInstall="$2"
         _inlineSed \
             etc/apps/paraview3/bashrc \
             '/^[^#]/s@ParaView_DIR=.*@ParaView_DIR='"$paraviewInstall@" \

From f240122ab080fd467bf6a4102684ac3fb991675d Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Wed, 20 Oct 2010 15:16:36 +0100
Subject: [PATCH 08/17] COMP: Using scotch 5.1.10a include files to match libs.

---
 src/parallel/decompose/ptscotchDecomp/Make/options | 2 +-
 src/parallel/decompose/scotchDecomp/Make/options   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/parallel/decompose/ptscotchDecomp/Make/options b/src/parallel/decompose/ptscotchDecomp/Make/options
index 488f5ef3bb..f3fa4b92b5 100644
--- a/src/parallel/decompose/ptscotchDecomp/Make/options
+++ b/src/parallel/decompose/ptscotchDecomp/Make/options
@@ -3,7 +3,7 @@ sinclude $(RULES)/mplib$(WM_MPLIB)
 
 EXE_INC = \
     $(PFLAGS) $(PINC) \
-    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/include \
+    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1.10a/include \
     -I/usr/include/scotch \
     -I../decompositionMethods/lnInclude
 
diff --git a/src/parallel/decompose/scotchDecomp/Make/options b/src/parallel/decompose/scotchDecomp/Make/options
index 2f01188ac6..10374b63c3 100644
--- a/src/parallel/decompose/scotchDecomp/Make/options
+++ b/src/parallel/decompose/scotchDecomp/Make/options
@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/include \
+    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1.10a/include \
     -I/usr/include/scotch \
     -I../decompositionMethods/lnInclude
 

From 0ad7b2cda88af24a648ae44cf1a203da7cef7071 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 21 Oct 2010 14:07:45 +0100
Subject: [PATCH 09/17] STYLE: foamToEnsight : indentation

---
 .../foamToEnsight/ensightField.C              | 172 +++++++++---------
 1 file changed, 86 insertions(+), 86 deletions(-)

diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
index 7ffa7a4ca5..5434209711 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C
@@ -575,99 +575,99 @@ void ensightPointField
     }
 
 
-     label ensightPatchI = eMesh.patchPartOffset();
+    label ensightPatchI = eMesh.patchPartOffset();
  
-     forAll(allPatchNames, patchi)
-     {
-         const word& patchName = allPatchNames[patchi];
+    forAll(allPatchNames, patchi)
+    {
+        const word& patchName = allPatchNames[patchi];
  
-         eMesh.barrier();
+        eMesh.barrier();
  
-         if (patchNames.empty() || patchNames.found(patchName))
-         {
-             const fvPatch& p = mesh.boundary()[patchi];
-             if
-             (
-                 returnReduce(p.size(), sumOp<label>())
-               > 0
-             )
-             {
-                 // Renumber the patch points/faces into unique points
-                 labelList pointToGlobal;
-                 labelList uniqueMeshPointLabels;
-                 autoPtr<globalIndex> globalPointsPtr =
-                 mesh.globalData().mergePoints
-                 (
-                     p.patch().meshPoints(),
-                     p.patch().meshPointMap(),
-                     pointToGlobal,
-                     uniqueMeshPointLabels
-                 );
+        if (patchNames.empty() || patchNames.found(patchName))
+        {
+            const fvPatch& p = mesh.boundary()[patchi];
+            if
+            (
+                returnReduce(p.size(), sumOp<label>())
+              > 0
+            )
+            {
+                // Renumber the patch points/faces into unique points
+                labelList pointToGlobal;
+                labelList uniqueMeshPointLabels;
+                autoPtr<globalIndex> globalPointsPtr =
+                mesh.globalData().mergePoints
+                (
+                    p.patch().meshPoints(),
+                    p.patch().meshPointMap(),
+                    pointToGlobal,
+                    uniqueMeshPointLabels
+                );
  
-                 if (Pstream::master())
-                 {
-                     ensightFile.writePartHeader(ensightPatchI);
-                 }
+                if (Pstream::master())
+                {
+                    ensightFile.writePartHeader(ensightPatchI);
+                }
  
-                 writeField
-                 (
-                     "coordinates",
-                     Field<Type>(pf.internalField(), uniqueMeshPointLabels),
-                     ensightFile
-                 );
+                writeField
+                (
+                    "coordinates",
+                    Field<Type>(pf.internalField(), uniqueMeshPointLabels),
+                    ensightFile
+                );
  
-                 ensightPatchI++;
-             }
-         }
-     }
+                ensightPatchI++;
+            }
+        }
+    }
  
-     // write faceZones, if requested
-     if (faceZoneNames.size())
-     {
-         forAllConstIter(wordHashSet, faceZoneNames, iter)
-         {
-             const word& faceZoneName = iter.key();
- 
-             eMesh.barrier();
- 
-             label zoneID = mesh.faceZones().findZoneID(faceZoneName);
- 
-             const faceZone& fz = mesh.faceZones()[zoneID];
- 
-             if (returnReduce(fz().nPoints(), sumOp<label>()) > 0)
-             {
-                 // Renumber the faceZone points/faces into unique points
-                 labelList pointToGlobal;
-                 labelList uniqueMeshPointLabels;
-                 autoPtr<globalIndex> globalPointsPtr =
-                 mesh.globalData().mergePoints
-                 (
-                     fz().meshPoints(),
-                     fz().meshPointMap(),
-                     pointToGlobal,
-                     uniqueMeshPointLabels
-                 );
- 
-                 if (Pstream::master())
-                 {
-                     ensightFile.writePartHeader(ensightPatchI);
-                 }
- 
-                 writeField
-                 (
-                     "coordinates",
-                     Field<Type>
-                     (
-                         pf.internalField(),
-                         uniqueMeshPointLabels
-                     ),
-                     ensightFile
-                 );
- 
-                 ensightPatchI++;
-             }
-         }
-     }
+    // write faceZones, if requested
+    if (faceZoneNames.size())
+    {
+        forAllConstIter(wordHashSet, faceZoneNames, iter)
+        {
+            const word& faceZoneName = iter.key();
+
+            eMesh.barrier();
+
+            label zoneID = mesh.faceZones().findZoneID(faceZoneName);
+
+            const faceZone& fz = mesh.faceZones()[zoneID];
+
+            if (returnReduce(fz().nPoints(), sumOp<label>()) > 0)
+            {
+                // Renumber the faceZone points/faces into unique points
+                labelList pointToGlobal;
+                labelList uniqueMeshPointLabels;
+                autoPtr<globalIndex> globalPointsPtr =
+                mesh.globalData().mergePoints
+                (
+                    fz().meshPoints(),
+                    fz().meshPointMap(),
+                    pointToGlobal,
+                    uniqueMeshPointLabels
+                );
+
+                if (Pstream::master())
+                {
+                    ensightFile.writePartHeader(ensightPatchI);
+                }
+
+                writeField
+                (
+                    "coordinates",
+                    Field<Type>
+                    (
+                        pf.internalField(),
+                        uniqueMeshPointLabels
+                    ),
+                    ensightFile
+                );
+
+                ensightPatchI++;
+            }
+        }
+    }
 
     if (Pstream::master())
     {

From 7e0f1e2a017ab10242eea726520106191f9c0d7a Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 21 Oct 2010 14:19:10 +0100
Subject: [PATCH 10/17] ENH: allocate fileMonitor only if wanted

---
 src/OSspecific/POSIX/fileMonitor.C | 24 +++++++++++++++---------
 src/OpenFOAM/db/Time/Time.C        | 13 ++++++++-----
 src/OpenFOAM/db/Time/Time.H        |  2 +-
 src/OpenFOAM/db/Time/TimeIO.C      |  5 +----
 4 files changed, 25 insertions(+), 19 deletions(-)

diff --git a/src/OSspecific/POSIX/fileMonitor.C b/src/OSspecific/POSIX/fileMonitor.C
index a6befd5908..4779142778 100644
--- a/src/OSspecific/POSIX/fileMonitor.C
+++ b/src/OSspecific/POSIX/fileMonitor.C
@@ -159,15 +159,21 @@ namespace Foam
         {
             if (inotifyFd_ < 0)
             {
-                WarningIn("fileMonitorWatcher(const label)")
-                    << "Failed allocating an inotify descriptor : "
-                    << string(strerror(errno)) << endl
-                    << "    Please increase the number of allowable "
-                    << "inotify instances" << endl
-                    << "    (/proc/sys/fs/inotify/max_user_instances on Linux)"
-                    << endl
-                    << "    Continuing without additional file monitoring."
-                    << endl;
+                static bool hasWarned = false;
+                if (!hasWarned)
+                {
+                    hasWarned = true;
+                    WarningIn("fileMonitorWatcher(const label)")
+                        << "Failed allocating an inotify descriptor : "
+                        << string(strerror(errno)) << endl
+                        << "    Please increase the number of allowable "
+                        << "inotify instances" << endl
+                        << "    (/proc/sys/fs/inotify/max_user_instances"
+                        << " on Linux)" << endl
+                        << "    or switch off runTimeModifiable." << endl
+                        << "    Continuing without additional file monitoring."
+                        << endl;
+                }
             }
         }
 
diff --git a/src/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C
index ff1247dcc7..5a5fbb8fa7 100644
--- a/src/OpenFOAM/db/Time/Time.C
+++ b/src/OpenFOAM/db/Time/Time.C
@@ -250,6 +250,7 @@ Foam::Time::Time
     // Time objects not registered so do like objectRegistry::checkIn ourselves.
     if (runTimeModifiable_)
     {
+        monitorPtr_.reset(new fileMonitor());
         controlDict_.watchIndex() = addWatch(controlDict_.filePath());
     }
 }
@@ -312,6 +313,8 @@ Foam::Time::Time
     // Time objects not registered so do like objectRegistry::checkIn ourselves.
     if (runTimeModifiable_)
     {
+        monitorPtr_.reset(new fileMonitor());
+
         // File might not exist yet.
         fileName f(controlDict_.filePath());
         if (f != fileName::null)
@@ -392,18 +395,18 @@ Foam::Time::~Time()
 
 Foam::label Foam::Time::addWatch(const fileName& fName) const
 {
-    return monitor_.addWatch(fName);
+    return monitorPtr_().addWatch(fName);
 }
 
 
 bool Foam::Time::removeWatch(const label watchIndex) const
 {
-    return monitor_.removeWatch(watchIndex);
+    return monitorPtr_().removeWatch(watchIndex);
 }
 
 const Foam::fileName& Foam::Time::getFile(const label watchIndex) const
 {
-    return monitor_.getFile(watchIndex);
+    return monitorPtr_().getFile(watchIndex);
 }
 
 
@@ -412,13 +415,13 @@ Foam::fileMonitor::fileState Foam::Time::getState
     const label watchFd
 ) const
 {
-    return monitor_.getState(watchFd);
+    return monitorPtr_().getState(watchFd);
 }
 
 
 void Foam::Time::setUnmodified(const label watchFd) const
 {
-    monitor_.setUnmodified(watchFd);
+    monitorPtr_().setUnmodified(watchFd);
 }
 
 
diff --git a/src/OpenFOAM/db/Time/Time.H b/src/OpenFOAM/db/Time/Time.H
index c644c757be..f9e38b58ee 100644
--- a/src/OpenFOAM/db/Time/Time.H
+++ b/src/OpenFOAM/db/Time/Time.H
@@ -73,7 +73,7 @@ class Time
     // Private data
 
         //- file-change monitor for all registered files
-        mutable fileMonitor monitor_;
+        mutable autoPtr<fileMonitor> monitorPtr_;
 
         //- The controlDict
         IOdictionary controlDict_;
diff --git a/src/OpenFOAM/db/Time/TimeIO.C b/src/OpenFOAM/db/Time/TimeIO.C
index e7c38a6aa6..86db2e1c2f 100644
--- a/src/OpenFOAM/db/Time/TimeIO.C
+++ b/src/OpenFOAM/db/Time/TimeIO.C
@@ -211,10 +211,7 @@ void Foam::Time::readModifiedObjects()
         // valid filePath).
         // Note: requires same ordering in objectRegistries on different
         // processors!
-        monitor_.updateStates(Pstream::parRun());
-
-//Pout<< "Time : runTimeModifiable_ and watchIndex:"
-//    << controlDict_.watchIndex() << endl;
+        monitorPtr_().updateStates(Pstream::parRun());
 
         // Time handling is special since controlDict_ is the one dictionary
         // that is not registered to any database.

From a0517f0ea9b8db5149e3acd152c45bb8e21d540f Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 21 Oct 2010 18:38:58 +0100
Subject: [PATCH 11/17] BUG: fan b.c. was using phi on slave side.

---
 .../constraint/cyclic/cyclicFvPatchField.C    | 17 ++++++++++++
 .../constraint/cyclic/cyclicFvPatchField.H    |  5 +++-
 .../derived/fan/fanFvPatchField.H             | 12 ++++++++-
 .../pimpleFoam/t-junction-with-fan/0/U        | 26 +++++++++----------
 .../pimpleFoam/t-junction-with-fan/0/epsilon  | 18 ++++++++-----
 .../pimpleFoam/t-junction-with-fan/0/k        | 12 ++++-----
 .../pimpleFoam/t-junction-with-fan/0/nuTilda  | 10 +++----
 .../pimpleFoam/t-junction-with-fan/0/nut      | 18 ++++++++-----
 .../pimpleFoam/t-junction-with-fan/0/p        | 26 +++++++++----------
 .../t-junction-with-fan/system/controlDict    |  3 +--
 10 files changed, 93 insertions(+), 54 deletions(-)

diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C
index 0010ab4a0a..c3949a6ac0 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C
@@ -171,6 +171,23 @@ tmp<Field<Type> > cyclicFvPatchField<Type>::patchNeighbourField() const
 }
 
 
+template<class Type>
+const cyclicFvPatchField<Type>& cyclicFvPatchField<Type>::neighbourPatchField()
+const
+{
+    const GeometricField<Type, fvPatchField, volMesh>& fld =
+    static_cast<const GeometricField<Type, fvPatchField, volMesh>&>
+    (
+        this->internalField()
+    );
+
+    return refCast<const cyclicFvPatchField<Type> >
+    (
+        fld.boundaryField()[this->cyclicPatch().neighbPatchID()]
+    );
+}
+
+
 template<class Type>
 void cyclicFvPatchField<Type>::updateInterfaceMatrix
 (
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H
index 654bd0831a..198dbab741 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H
@@ -150,9 +150,12 @@ public:
 
         // Evaluation functions
 
-            //- Return neighbour coupled given internal cell data
+            //- Return neighbour coupled internal cell data
             tmp<Field<Type> > patchNeighbourField() const;
 
+            //- Return reference to neighbour patchField
+            const cyclicFvPatchField<Type>& neighbourPatchField() const;
+
             //- Update result field based on interface functionality
             virtual void updateInterfaceMatrix
             (
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
index 73014e048a..31c51c999d 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H
@@ -140,7 +140,17 @@ public:
             //- Return the "jump" across the patch.
             virtual tmp<Field<Type> > jump() const
             {
-                return jump_;
+                if (this->cyclicPatch().owner())
+                {
+                    return jump_;
+                }
+                else
+                {
+                    return refCast<const fanFvPatchField<Type> >
+                    (
+                        this->neighbourPatchField()
+                    ).jump();
+                }
             }
 
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U
index 51ef892fb8..f1006ceb5e 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U
@@ -15,46 +15,46 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 1 -1 0 0 0 0 ];
+dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform ( 0 0 0 );
+internalField   uniform (0 0 0);
 
 boundaryField
 {
     inlet
     {
         type            pressureInletOutletVelocity;
-        value           uniform ( 0 0 0 );
+        value           uniform (0 0 0);
     }
     outlet1
     {
         type            inletOutlet;
-        inletValue      uniform ( 0 0 0 );
-        value           uniform ( 0 0 0 );
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
     }
     outlet2
     {
         type            inletOutlet;
-        inletValue      uniform ( 0 0 0 );
-        value           uniform ( 0 0 0 );
+        inletValue      uniform (0 0 0);
+        value           uniform (0 0 0);
     }
     baffles
     {
         type            fixedValue;
-        value           uniform ( 0 0 0 );
+        value           uniform (0 0 0);
     }
     fan_half0
     {
         type            cyclic;
     }
+    fan_half1
+    {
+        type            cyclic;
+    }
     defaultFaces
     {
         type            fixedValue;
-        value           uniform ( 0 0 0 );
-    }
-    fan_half1
-    {
-        type            cyclic;
+        value           uniform (0 0 0);
     }
 }
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon
index dc1eb9da3e..5c6c53664d 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 2 -3 0 0 0 0 ];
+dimensions      [0 2 -3 0 0 0 0];
 
 internalField   uniform 1;
 
@@ -42,21 +42,27 @@ boundaryField
     baffles
     {
         type            epsilonWallFunction;
-        value           uniform 1;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
     }
     fan_half0
     {
         type            cyclic;
     }
+    fan_half1
+    {
+        type            cyclic;
+    }
     defaultFaces
     {
         type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
         value           uniform 1;
     }
-    fan_half1
-    {
-        type            cyclic;
-    }
 }
 
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k
index 06aeff0f1f..18cb53d519 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 2 -2 0 0 0 0 ];
+dimensions      [0 2 -2 0 0 0 0];
 
 internalField   uniform 1;
 
@@ -42,21 +42,21 @@ boundaryField
     baffles
     {
         type            kqRWallFunction;
-        value           uniform 1;
+        value           uniform 0;
     }
     fan_half0
     {
         type            cyclic;
     }
+    fan_half1
+    {
+        type            cyclic;
+    }
     defaultFaces
     {
         type            kqRWallFunction;
         value           uniform 1;
     }
-    fan_half1
-    {
-        type            cyclic;
-    }
 }
 
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda
index 699a951041..ffc7ac906f 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 2 -1 0 0 0 0 ];
+dimensions      [0 2 -1 0 0 0 0];
 
 internalField   uniform 0;
 
@@ -41,14 +41,14 @@ boundaryField
     {
         type            cyclic;
     }
-    defaultFaces
-    {
-        type            zeroGradient;
-    }
     fan_half1
     {
         type            cyclic;
     }
+    defaultFaces
+    {
+        type            zeroGradient;
+    }
 }
 
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut
index 1a5afb9689..6ed8f42dfc 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 2 -1 0 0 0 0 ];
+dimensions      [0 2 -1 0 0 0 0];
 
 internalField   uniform 0;
 
@@ -39,21 +39,27 @@ boundaryField
     baffles
     {
         type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
         value           uniform 0;
     }
     fan_half0
     {
         type            cyclic;
     }
-    defaultFaces
-    {
-        type            nutkWallFunction;
-        value           uniform 0;
-    }
     fan_half1
     {
         type            cyclic;
     }
+    defaultFaces
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
 }
 
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p
index c77ad85725..66f5c1d712 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 2 -2 0 0 0 0 ];
+dimensions      [0 2 -2 0 0 0 0];
 
 internalField   uniform 100000;
 
@@ -24,14 +24,12 @@ boundaryField
     inlet
     {
         type            timeVaryingTotalPressure;
-        p0              100040;
-        outOfBounds     clamp;
-        fileName        "$FOAM_CASE/constant/p0vsTime";
-        U               U;
-        phi             phi;
         rho             none;
         psi             none;
         gamma           1;
+        p0              100040;
+        fileName        "$FOAM_CASE/constant/p0vsTime";
+        outOfBounds     clamp;
         value           uniform 100040;
     }
     outlet1
@@ -52,19 +50,19 @@ boundaryField
     {
         type            fan;
         patchType       cyclic;
-        f               2 ( 50 -0.1 );
-        value           $internalField;
-    }
-    defaultFaces
-    {
-        type            zeroGradient;
+        f               2(100 -0.1);
+        value           uniform 0;
     }
     fan_half1
     {
         type            fan;
         patchType       cyclic;
-        f               2 ( 50 -0.1 );
-        value           $internalField;
+        f               2(100 -0.1);
+        value           uniform 0;
+    }
+    defaultFaces
+    {
+        type            zeroGradient;
     }
 }
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict
index 9dc31d5561..63b2669008 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict
@@ -28,7 +28,6 @@ endTime         1;
 deltaT          0.001;
 
 writeControl    adjustableRunTime;
-
 writeInterval   0.1;
 
 purgeWrite      0;
@@ -47,7 +46,7 @@ runTimeModifiable true;
 
 adjustTimeStep  yes;
 
-maxCo           5;
+maxCo           3;
 
 libs
 (

From e44359e19011486fb9cefefe2f321f2be3167a6a Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 22 Oct 2010 11:50:28 +0100
Subject: [PATCH 12/17] BUG: cyclics in mesh conversion : have 'old' style
 cyclics

---
 .../fluent3DMeshToFoam/fluent3DMeshToFoam.L   |    4 +-
 .../mesh/conversion/kivaToFoam/kivaToFoam.C   |    2 +-
 .../mesh/conversion/kivaToFoam/readKivaGrid.H |    2 +-
 .../mesh/conversion/sammToFoam/readBoundary.C |    4 +-
 .../conversion/star3ToFoam/readBoundary.C     |    4 +-
 src/OpenFOAM/Make/files                       |    1 +
 .../constraint/oldCyclic/oldCyclicPolyPatch.C | 1298 +++++++++++++++++
 .../constraint/oldCyclic/oldCyclicPolyPatch.H |  318 ++++
 8 files changed, 1625 insertions(+), 8 deletions(-)
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
 create mode 100644 src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H

diff --git a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
index 5fd94ef38d..6f0f7c976d 100644
--- a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
+++ b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L
@@ -43,7 +43,7 @@ Description
 #include "polyMeshZipUpCells.H"
 #include "wallPolyPatch.H"
 #include "symmetryPolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "Swap.H"
 #include "IFstream.H"
 #include "readHexLabel.H"
@@ -900,7 +900,7 @@ int main(int argc, char *argv[])
     fluentToFoamType.insert("interface", polyPatch::typeName);
     fluentToFoamType.insert("internal", polyPatch::typeName);
     fluentToFoamType.insert("solid", polyPatch::typeName);
-    fluentToFoamType.insert("fan", cyclicPolyPatch::typeName);
+    fluentToFoamType.insert("fan", oldCyclicPolyPatch::typeName);
     fluentToFoamType.insert("radiator", polyPatch::typeName);
     fluentToFoamType.insert("porous-jump", polyPatch::typeName);
 
diff --git a/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C b/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C
index ab764547a3..066fab386f 100644
--- a/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C
+++ b/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C
@@ -41,7 +41,7 @@ Description
 #include "wallPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "wedgePolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "unitConversion.H"
 
 using namespace Foam;
diff --git a/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H b/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H
index 04b36020ea..1a559559e9 100644
--- a/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H
+++ b/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H
@@ -196,7 +196,7 @@ const word* kivaPatchTypes[nBCs] =
     &polyPatch::typeName,
     &polyPatch::typeName,
     &symmetryPolyPatch::typeName,
-    &cyclicPolyPatch::typeName
+    &oldCyclicPolyPatch::typeName
 };
 
 enum patchTypeNames
diff --git a/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C b/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C
index 53dffd39c6..d01a114111 100644
--- a/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C
+++ b/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C
@@ -29,7 +29,7 @@ Description
 #include "sammMesh.H"
 #include "Time.H"
 #include "wallPolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "preservePatchTypes.H"
 #include "IFstream.H"
@@ -208,7 +208,7 @@ void sammMesh::readBoundary()
             {
                 // incorrect. should be cyclicPatch but this
                 // requires info on connected faces.
-                patchTypes_[patchLabel] = cyclicPolyPatch::typeName;
+                patchTypes_[patchLabel] = oldCyclicPolyPatch::typeName;
             }
             else
             {
diff --git a/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C b/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C
index a96ee6e96b..75a1860006 100644
--- a/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C
+++ b/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C
@@ -29,7 +29,7 @@ Description
 #include "starMesh.H"
 #include "Time.H"
 #include "wallPolyPatch.H"
-#include "cyclicPolyPatch.H"
+#include "oldCyclicPolyPatch.H"
 #include "symmetryPolyPatch.H"
 #include "preservePatchTypes.H"
 #include "IFstream.H"
@@ -206,7 +206,7 @@ void starMesh::readBoundary()
             {
                 // incorrect. should be cyclicPatch but this
                 // requires info on connected faces.
-                patchTypes_[patchLabel] = cyclicPolyPatch::typeName;
+                patchTypes_[patchLabel] = oldCyclicPolyPatch::typeName;
             }
             else
             {
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 97513b339d..09fa5faa94 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -346,6 +346,7 @@ $(basicPolyPatches)/generic/genericPolyPatch.C
 constraintPolyPatches = $(polyPatches)/constraint
 $(constraintPolyPatches)/cyclic/cyclicPolyPatch.C
 $(constraintPolyPatches)/cyclicSlip/cyclicSlipPolyPatch.C
+$(constraintPolyPatches)/oldCyclic/oldCyclicPolyPatch.C
 $(constraintPolyPatches)/empty/emptyPolyPatch.C
 $(constraintPolyPatches)/nonuniformTransformCyclic/nonuniformTransformCyclicPolyPatch.C
 $(constraintPolyPatches)/processorCyclic/processorCyclicPolyPatch.C
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
new file mode 100644
index 0000000000..2dcbfef30e
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
@@ -0,0 +1,1298 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "oldCyclicPolyPatch.H"
+#include "addToRunTimeSelectionTable.H"
+#include "polyBoundaryMesh.H"
+#include "polyMesh.H"
+#include "demandDrivenData.H"
+#include "OFstream.H"
+#include "patchZones.H"
+#include "matchPoints.H"
+#include "Time.H"
+#include "transformList.H"
+#include "cyclicPolyPatch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(oldCyclicPolyPatch, 0);
+
+    addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, word);
+    addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, dictionary);
+
+
+template<>
+const char* NamedEnum<oldCyclicPolyPatch::transformType, 3>::names[] =
+{
+    "unknown",
+    "rotational",
+    "translational"
+};
+
+const NamedEnum<oldCyclicPolyPatch::transformType, 3>
+    oldCyclicPolyPatch::transformTypeNames;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::pointField Foam::oldCyclicPolyPatch::calcFaceCentres
+(
+    const UList<face>& faces,
+    const pointField& points
+)
+{
+    pointField ctrs(faces.size());
+
+    forAll(faces, faceI)
+    {
+        ctrs[faceI] = faces[faceI].centre(points);
+    }
+
+    return ctrs;
+}
+
+
+Foam::pointField Foam::oldCyclicPolyPatch::getAnchorPoints
+(
+    const UList<face>& faces,
+    const pointField& points
+)
+{
+    pointField anchors(faces.size());
+
+    forAll(faces, faceI)
+    {
+        anchors[faceI] = points[faces[faceI][0]];
+    }
+
+    return anchors;
+}
+
+
+Foam::label Foam::oldCyclicPolyPatch::findMaxArea
+(
+    const pointField& points,
+    const faceList& faces
+)
+{
+    label maxI = -1;
+    scalar maxAreaSqr = -GREAT;
+
+    forAll(faces, faceI)
+    {
+        scalar areaSqr = magSqr(faces[faceI].normal(points));
+
+        if (areaSqr > maxAreaSqr)
+        {
+            maxAreaSqr = areaSqr;
+            maxI = faceI;
+        }
+    }
+    return maxI;
+}
+
+
+// Get geometric zones of patch by looking at normals.
+// Method 1: any edge with sharpish angle is edge between two halves.
+//           (this will handle e.g. wedge geometries).
+//           Also two fully disconnected regions will be handled this way.
+// Method 2: sort faces into two halves based on face normal.
+bool Foam::oldCyclicPolyPatch::getGeometricHalves
+(
+    const primitivePatch& pp,
+    labelList& half0ToPatch,
+    labelList& half1ToPatch
+) const
+{
+    // Calculate normals
+    const vectorField& faceNormals = pp.faceNormals();
+
+    // Find edges with sharp angles.
+    boolList regionEdge(pp.nEdges(), false);
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+
+    label nRegionEdges = 0;
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        // Check manifold edges for sharp angle.
+        // (Non-manifold already handled by patchZones)
+        if (eFaces.size() == 2)
+        {
+            if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
+            {
+                regionEdge[edgeI] = true;
+
+                nRegionEdges++;
+            }
+        }
+    }
+
+
+    // For every face determine zone it is connected to (without crossing
+    // any regionEdge)
+    patchZones ppZones(pp, regionEdge);
+
+    if (debug)
+    {
+        Pout<< "oldCyclicPolyPatch::getGeometricHalves : "
+            << "Found " << nRegionEdges << " edges on patch " << name()
+            << " where the cos of the angle between two connected faces"
+            << " was less than " << featureCos_ << nl
+            << "Patch divided by these and by single sides edges into "
+            << ppZones.nZones() << " parts." << endl;
+
+
+        // Dumping zones to obj files.
+
+        labelList nZoneFaces(ppZones.nZones());
+
+        for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
+        {
+            OFstream stream
+            (
+                boundaryMesh().mesh().time().path()
+               /name()+"_zone_"+Foam::name(zoneI)+".obj"
+            );
+            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing zone "
+                << zoneI << " face centres to OBJ file " << stream.name()
+                << endl;
+
+            labelList zoneFaces(findIndices(ppZones, zoneI));
+
+            forAll(zoneFaces, i)
+            {
+                writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
+            }
+
+            nZoneFaces[zoneI] = zoneFaces.size();
+        }
+    }
+
+
+    if (ppZones.nZones() == 2)
+    {
+        half0ToPatch = findIndices(ppZones, 0);
+        half1ToPatch = findIndices(ppZones, 1);
+    }
+    else
+    {
+        if (debug)
+        {
+            Pout<< "oldCyclicPolyPatch::getGeometricHalves :"
+                << " falling back to face-normal comparison" << endl;
+        }
+        label n0Faces = 0;
+        half0ToPatch.setSize(pp.size());
+
+        label n1Faces = 0;
+        half1ToPatch.setSize(pp.size());
+
+        // Compare to face 0 normal.
+        forAll(faceNormals, faceI)
+        {
+            if ((faceNormals[faceI] & faceNormals[0]) > 0)
+            {
+                half0ToPatch[n0Faces++] = faceI;
+            }
+            else
+            {
+                half1ToPatch[n1Faces++] = faceI;
+            }
+        }
+        half0ToPatch.setSize(n0Faces);
+        half1ToPatch.setSize(n1Faces);
+
+        if (debug)
+        {
+            Pout<< "oldCyclicPolyPatch::getGeometricHalves :"
+                << " Number of faces per zone:("
+                << n0Faces << ' ' << n1Faces << ')' << endl;
+        }
+    }
+
+    if (half0ToPatch.size() != half1ToPatch.size())
+    {
+        fileName casePath(boundaryMesh().mesh().time().path());
+
+        // Dump halves
+        {
+            fileName nm0(casePath/name()+"_half0_faces.obj");
+            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0"
+                << " faces to OBJ file " << nm0 << endl;
+            writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
+
+            fileName nm1(casePath/name()+"_half1_faces.obj");
+            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1"
+                << " faces to OBJ file " << nm1 << endl;
+            writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
+        }
+
+        // Dump face centres
+        {
+            OFstream str0(casePath/name()+"_half0.obj");
+            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0"
+                << " face centres to OBJ file " << str0.name() << endl;
+
+            forAll(half0ToPatch, i)
+            {
+                writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
+            }
+
+            OFstream str1(casePath/name()+"_half1.obj");
+            Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1"
+                << " face centres to OBJ file " << str1.name() << endl;
+            forAll(half1ToPatch, i)
+            {
+                writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
+            }
+        }
+
+        SeriousErrorIn
+        (
+            "oldCyclicPolyPatch::getGeometricHalves"
+            "(const primitivePatch&, labelList&, labelList&) const"
+        )   << "Patch " << name() << " gets decomposed in two zones of"
+            << "inequal size: " << half0ToPatch.size()
+            << " and " << half1ToPatch.size() << endl
+            << "This means that the patch is either not two separate regions"
+            << " or one region where the angle between the different regions"
+            << " is not sufficiently sharp." << endl
+            << "Please adapt the featureCos setting." << endl
+            << "Continuing with incorrect face ordering from now on!" << endl;
+
+        return false;
+    }
+    else
+    {
+        return true;
+    }
+}
+
+
+// Given a split of faces into left and right half calculate the centres
+// and anchor points. Transform the left points so they align with the
+// right ones.
+void Foam::oldCyclicPolyPatch::getCentresAndAnchors
+(
+    const primitivePatch& pp,
+    const faceList& half0Faces,
+    const faceList& half1Faces,
+
+    pointField& ppPoints,
+    pointField& half0Ctrs,
+    pointField& half1Ctrs,
+    pointField& anchors0,
+    scalarField& tols
+) const
+{
+    // Get geometric data on both halves.
+    half0Ctrs = calcFaceCentres(half0Faces, pp.points());
+    anchors0 = getAnchorPoints(half0Faces, pp.points());
+    half1Ctrs = calcFaceCentres(half1Faces, pp.points());
+
+    switch (transform_)
+    {
+        case ROTATIONAL:
+        {
+            label face0 = getConsistentRotationFace(half0Ctrs);
+            label face1 = getConsistentRotationFace(half1Ctrs);
+
+            vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
+            vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
+            n0 /= mag(n0) + VSMALL;
+            n1 /= mag(n1) + VSMALL;
+
+            if (debug)
+            {
+                Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
+                    << " Specified rotation :"
+                    << " n0:" << n0 << " n1:" << n1 << endl;
+            }
+
+            // Rotation (around origin)
+            const tensor reverseT(rotationTensor(n0, -n1));
+
+            // Rotation
+            forAll(half0Ctrs, faceI)
+            {
+                half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
+                anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
+            }
+
+            ppPoints = Foam::transform(reverseT, pp.points());
+
+            break;
+        }
+        //- Problem: usually specified translation is not accurate enough
+        //- to get proper match so keep automatic determination over here.
+        //case TRANSLATIONAL:
+        //{
+        //    // Transform 0 points.
+        //
+        //    if (debug)
+        //    {
+        //        Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
+        //            << "Specified translation : " << separationVector_ << endl;
+        //    }
+        //
+        //    half0Ctrs += separationVector_;
+        //    anchors0 += separationVector_;
+        //    break;
+        //}
+        default:
+        {
+            // Assumes that cyclic is planar. This is also the initial
+            // condition for patches without faces.
+
+            // Determine the face with max area on both halves. These
+            // two faces are used to determine the transformation tensors
+            label max0I = findMaxArea(pp.points(), half0Faces);
+            vector n0 = half0Faces[max0I].normal(pp.points());
+            n0 /= mag(n0) + VSMALL;
+
+            label max1I = findMaxArea(pp.points(), half1Faces);
+            vector n1 = half1Faces[max1I].normal(pp.points());
+            n1 /= mag(n1) + VSMALL;
+
+            if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
+            {
+                if (debug)
+                {
+                    Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
+                        << " Detected rotation :"
+                        << " n0:" << n0 << " n1:" << n1 << endl;
+                }
+
+                // Rotation (around origin)
+                const tensor reverseT(rotationTensor(n0, -n1));
+
+                // Rotation
+                forAll(half0Ctrs, faceI)
+                {
+                    half0Ctrs[faceI] = Foam::transform
+                    (
+                        reverseT,
+                        half0Ctrs[faceI]
+                    );
+                    anchors0[faceI] = Foam::transform
+                    (
+                        reverseT,
+                        anchors0[faceI]
+                    );
+                }
+                ppPoints = Foam::transform(reverseT, pp.points());
+            }
+            else
+            {
+                // Parallel translation. Get average of all used points.
+
+                primitiveFacePatch half0(half0Faces, pp.points());
+                const pointField& half0Pts = half0.localPoints();
+                const point ctr0(sum(half0Pts)/half0Pts.size());
+
+                primitiveFacePatch half1(half1Faces, pp.points());
+                const pointField& half1Pts = half1.localPoints();
+                const point ctr1(sum(half1Pts)/half1Pts.size());
+
+                if (debug)
+                {
+                    Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
+                        << " Detected translation :"
+                        << " n0:" << n0 << " n1:" << n1
+                        << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
+                }
+
+                half0Ctrs += ctr1 - ctr0;
+                anchors0 += ctr1 - ctr0;
+                ppPoints = pp.points() + ctr1 - ctr0;
+            }
+            break;
+        }
+    }
+
+
+    // Calculate typical distance per face
+    tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
+}
+
+
+// Calculates faceMap and rotation. Returns true if all ok.
+bool Foam::oldCyclicPolyPatch::matchAnchors
+(
+    const bool report,
+    const primitivePatch& pp,
+    const labelList& half0ToPatch,
+    const pointField& anchors0,
+
+    const labelList& half1ToPatch,
+    const faceList& half1Faces,
+    const labelList& from1To0,
+
+    const scalarField& tols,
+
+    labelList& faceMap,
+    labelList& rotation
+) const
+{
+    // Set faceMap such that half0 faces get first and corresponding half1
+    // faces last.
+
+    forAll(half0ToPatch, half0FaceI)
+    {
+        // Label in original patch
+        label patchFaceI = half0ToPatch[half0FaceI];
+
+        faceMap[patchFaceI] = half0FaceI;
+
+        // No rotation
+        rotation[patchFaceI] = 0;
+    }
+
+    bool fullMatch = true;
+
+    forAll(from1To0, half1FaceI)
+    {
+        label patchFaceI = half1ToPatch[half1FaceI];
+
+        // This face has to match the corresponding one on half0.
+        label half0FaceI = from1To0[half1FaceI];
+
+        label newFaceI = half0FaceI + pp.size()/2;
+
+        faceMap[patchFaceI] = newFaceI;
+
+        // Rotate patchFaceI such that its f[0] aligns with that of
+        // the corresponding face
+        // (which after shuffling will be at position half0FaceI)
+
+        const point& wantedAnchor = anchors0[half0FaceI];
+
+        rotation[newFaceI] = getRotation
+        (
+            pp.points(),
+            half1Faces[half1FaceI],
+            wantedAnchor,
+            tols[half1FaceI]
+        );
+
+        if (rotation[newFaceI] == -1)
+        {
+            fullMatch = false;
+
+            if (report)
+            {
+                const face& f = half1Faces[half1FaceI];
+                SeriousErrorIn
+                (
+                    "oldCyclicPolyPatch::matchAnchors(..)"
+                )   << "Patch:" << name() << " : "
+                    << "Cannot find point on face " << f
+                    << " with vertices:"
+                    << UIndirectList<point>(pp.points(), f)()
+                    << " that matches point " << wantedAnchor
+                    << " when matching the halves of cyclic patch " << name()
+                    << endl
+                    << "Continuing with incorrect face ordering from now on!"
+                    << endl;
+            }
+        }
+    }
+    return fullMatch;
+}
+
+
+Foam::label Foam::oldCyclicPolyPatch::getConsistentRotationFace
+(
+    const pointField& faceCentres
+) const
+{
+    const scalarField magRadSqr =
+        magSqr((faceCentres - rotationCentre_) ^ rotationAxis_);
+    scalarField axisLen = (faceCentres - rotationCentre_) & rotationAxis_;
+    axisLen = axisLen - min(axisLen);
+    const scalarField magLenSqr = magRadSqr + axisLen*axisLen;
+
+    label rotFace = -1;
+    scalar maxMagLenSqr = -GREAT;
+    scalar maxMagRadSqr = -GREAT;
+    forAll(faceCentres, i)
+    {
+        if (magLenSqr[i] >= maxMagLenSqr)
+        {
+            if (magRadSqr[i] > maxMagRadSqr)
+            {
+                rotFace = i;
+                maxMagLenSqr = magLenSqr[i];
+                maxMagRadSqr = magRadSqr[i];
+            }
+        }
+    }
+
+    if (debug)
+    {
+        Info<< "getConsistentRotationFace(const pointField&)" << nl
+            << "    rotFace = " << rotFace << nl
+            << "    point =  " << faceCentres[rotFace] << endl;
+    }
+
+    return rotFace;
+}
+
+
+// * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
+
+Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
+(
+    const word& name,
+    const label size,
+    const label start,
+    const label index,
+    const polyBoundaryMesh& bm
+)
+:
+    coupledPolyPatch(name, size, start, index, bm),
+    featureCos_(0.9),
+    transform_(UNKNOWN),
+    rotationAxis_(vector::zero),
+    rotationCentre_(point::zero),
+    separationVector_(vector::zero)
+{}
+
+
+Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
+(
+    const word& name,
+    const dictionary& dict,
+    const label index,
+    const polyBoundaryMesh& bm
+)
+:
+    coupledPolyPatch(name, dict, index, bm),
+    featureCos_(0.9),
+    transform_(UNKNOWN),
+    rotationAxis_(vector::zero),
+    rotationCentre_(point::zero),
+    separationVector_(vector::zero)
+{
+    if (dict.found("neighbourPatch"))
+    {
+        FatalIOErrorIn
+        (
+            "oldCyclicPolyPatch::oldCyclicPolyPatch\n"
+            "(\n"
+            "    const word& name,\n"
+            "    const dictionary& dict,\n"
+            "    const label index,\n"
+            "    const polyBoundaryMesh& bm\n"
+            ")",
+            dict
+        )   << "Found \"neighbourPatch\" entry when reading cyclic patch "
+            << name << endl
+            << "Is this mesh already with split cyclics?" << endl
+            << "If so run a newer version that supports it"
+            << ", if not comment out the \"neighbourPatch\" entry and re-run"
+            << exit(FatalIOError);
+    }
+
+    dict.readIfPresent("featureCos", featureCos_);
+
+    if (dict.found("transform"))
+    {
+        transform_ = transformTypeNames.read(dict.lookup("transform"));
+        switch (transform_)
+        {
+            case ROTATIONAL:
+            {
+                dict.lookup("rotationAxis") >> rotationAxis_;
+                dict.lookup("rotationCentre") >> rotationCentre_;
+                break;
+            }
+            case TRANSLATIONAL:
+            {
+                dict.lookup("separationVector") >> separationVector_;
+                break;
+            }
+            default:
+            {
+                // no additional info required
+            }
+        }
+    }
+}
+
+
+Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
+(
+    const oldCyclicPolyPatch& pp,
+    const polyBoundaryMesh& bm
+)
+:
+    coupledPolyPatch(pp, bm),
+    featureCos_(pp.featureCos_),
+    transform_(pp.transform_),
+    rotationAxis_(pp.rotationAxis_),
+    rotationCentre_(pp.rotationCentre_),
+    separationVector_(pp.separationVector_)
+{}
+
+
+Foam::oldCyclicPolyPatch::oldCyclicPolyPatch
+(
+    const oldCyclicPolyPatch& pp,
+    const polyBoundaryMesh& bm,
+    const label index,
+    const label newSize,
+    const label newStart
+)
+:
+    coupledPolyPatch(pp, bm, index, newSize, newStart),
+    featureCos_(pp.featureCos_),
+    transform_(pp.transform_),
+    rotationAxis_(pp.rotationAxis_),
+    rotationCentre_(pp.rotationCentre_),
+    separationVector_(pp.separationVector_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::oldCyclicPolyPatch::~oldCyclicPolyPatch()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::oldCyclicPolyPatch::initGeometry(PstreamBuffers& pBufs)
+{
+    polyPatch::initGeometry(pBufs);
+}
+
+
+void Foam::oldCyclicPolyPatch::calcGeometry
+(
+    const primitivePatch& referPatch,
+    const UList<point>& thisCtrs,
+    const UList<point>& thisAreas,
+    const UList<point>& thisCc,
+    const UList<point>& nbrCtrs,
+    const UList<point>& nbrAreas,
+    const UList<point>& nbrCc
+)
+{}
+
+
+void Foam::oldCyclicPolyPatch::calcGeometry(PstreamBuffers& pBufs)
+{}
+
+
+void Foam::oldCyclicPolyPatch::initMovePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField& p
+)
+{
+    polyPatch::initMovePoints(pBufs, p);
+}
+
+
+void Foam::oldCyclicPolyPatch::movePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField& p
+)
+{
+    polyPatch::movePoints(pBufs, p);
+}
+
+
+void Foam::oldCyclicPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
+{
+    polyPatch::initUpdateMesh(pBufs);
+}
+
+
+void Foam::oldCyclicPolyPatch::updateMesh(PstreamBuffers& pBufs)
+{
+    polyPatch::updateMesh(pBufs);
+}
+
+
+void Foam::oldCyclicPolyPatch::initOrder
+(
+    PstreamBuffers&,
+    const primitivePatch& pp
+) const
+{}
+
+
+//  Return new ordering. Ordering is -faceMap: for every face index
+//  the new face -rotation:for every new face the clockwise shift
+//  of the original face. Return false if nothing changes (faceMap
+//  is identity, rotation is 0)
+bool Foam::oldCyclicPolyPatch::order
+(
+    PstreamBuffers&,
+    const primitivePatch& pp,
+    labelList& faceMap,
+    labelList& rotation
+) const
+{
+    faceMap.setSize(pp.size());
+    faceMap = -1;
+
+    rotation.setSize(pp.size());
+    rotation = 0;
+
+    if (pp.empty())
+    {
+        // No faces, nothing to change.
+        return false;
+    }
+
+    if (pp.size()&1)
+    {
+        FatalErrorIn("oldCyclicPolyPatch::order(..)")
+            << "Size of cyclic " << name() << " should be a multiple of 2"
+            << ". It is " << pp.size() << abort(FatalError);
+    }
+
+    label halfSize = pp.size()/2;
+
+    // Supplied primitivePatch already with new points.
+    // Cyclics are limited to one transformation tensor
+    // currently anyway (i.e. straight plane) so should not be too big a
+    // problem.
+
+
+    // Indices of faces on half0
+    labelList half0ToPatch;
+    // Indices of faces on half1
+    labelList half1ToPatch;
+
+
+    // 1. Test if already correctly oriented by starting from trivial ordering.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    half0ToPatch = identity(halfSize);
+    half1ToPatch = half0ToPatch + halfSize;
+
+    // Get faces
+    faceList half0Faces(UIndirectList<face>(pp, half0ToPatch));
+    faceList half1Faces(UIndirectList<face>(pp, half1ToPatch));
+
+    // Get geometric quantities
+    pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
+    scalarField tols;
+    getCentresAndAnchors
+    (
+        pp,
+        half0Faces,
+        half1Faces,
+
+        ppPoints,
+        half0Ctrs,
+        half1Ctrs,
+        anchors0,
+        tols
+    );
+
+    // Geometric match of face centre vectors
+    labelList from1To0;
+    bool matchedAll = matchPoints
+    (
+        half1Ctrs,
+        half0Ctrs,
+        tols,
+        false,
+        from1To0
+    );
+
+    if (debug)
+    {
+        Pout<< "oldCyclicPolyPatch::order : test if already ordered:"
+            << matchedAll << endl;
+
+        // Dump halves
+        fileName nm0("match1_"+name()+"_half0_faces.obj");
+        Pout<< "oldCyclicPolyPatch::order : Writing half0"
+            << " faces to OBJ file " << nm0 << endl;
+        writeOBJ(nm0, half0Faces, ppPoints);
+
+        fileName nm1("match1_"+name()+"_half1_faces.obj");
+        Pout<< "oldCyclicPolyPatch::order : Writing half1"
+            << " faces to OBJ file " << nm1 << endl;
+        writeOBJ(nm1, half1Faces, pp.points());
+
+        OFstream ccStr
+        (
+            boundaryMesh().mesh().time().path()
+           /"match1_"+ name() + "_faceCentres.obj"
+        );
+        Pout<< "oldCyclicPolyPatch::order : "
+            << "Dumping currently found cyclic match as lines between"
+            << " corresponding face centres to file " << ccStr.name()
+            << endl;
+
+        // Recalculate untransformed face centres
+        //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+        label vertI = 0;
+
+        forAll(half1Ctrs, i)
+        {
+            //if (from1To0[i] != -1)
+            {
+                // Write edge between c1 and c0
+                //const point& c0 = rawHalf0Ctrs[from1To0[i]];
+                //const point& c0 = half0Ctrs[from1To0[i]];
+                const point& c0 = half0Ctrs[i];
+                const point& c1 = half1Ctrs[i];
+                writeOBJ(ccStr, c0, c1, vertI);
+            }
+        }
+    }
+
+
+    // 2. Ordered in pairs (so 0,1 coupled and 2,3 etc.)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (!matchedAll)
+    {
+        label faceI = 0;
+        for (label i = 0; i < halfSize; i++)
+        {
+            half0ToPatch[i] = faceI++;
+            half1ToPatch[i] = faceI++;
+        }
+
+        // And redo all matching
+        half0Faces = UIndirectList<face>(pp, half0ToPatch);
+        half1Faces = UIndirectList<face>(pp, half1ToPatch);
+
+        getCentresAndAnchors
+        (
+            pp,
+            half0Faces,
+            half1Faces,
+
+            ppPoints,
+            half0Ctrs,
+            half1Ctrs,
+            anchors0,
+            tols
+        );
+
+        // Geometric match of face centre vectors
+        matchedAll = matchPoints
+        (
+            half1Ctrs,
+            half0Ctrs,
+            tols,
+            false,
+            from1To0
+        );
+
+        if (debug)
+        {
+            Pout<< "oldCyclicPolyPatch::order : test if pairwise ordered:"
+                << matchedAll << endl;
+            // Dump halves
+            fileName nm0("match2_"+name()+"_half0_faces.obj");
+            Pout<< "oldCyclicPolyPatch::order : Writing half0"
+                << " faces to OBJ file " << nm0 << endl;
+            writeOBJ(nm0, half0Faces, ppPoints);
+
+            fileName nm1("match2_"+name()+"_half1_faces.obj");
+            Pout<< "oldCyclicPolyPatch::order : Writing half1"
+                << " faces to OBJ file " << nm1 << endl;
+            writeOBJ(nm1, half1Faces, pp.points());
+
+            OFstream ccStr
+            (
+                boundaryMesh().mesh().time().path()
+               /"match2_"+name()+"_faceCentres.obj"
+            );
+            Pout<< "oldCyclicPolyPatch::order : "
+                << "Dumping currently found cyclic match as lines between"
+                << " corresponding face centres to file " << ccStr.name()
+                << endl;
+
+            // Recalculate untransformed face centres
+            //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+            label vertI = 0;
+
+            forAll(half1Ctrs, i)
+            {
+                if (from1To0[i] != -1)
+                {
+                    // Write edge between c1 and c0
+                    //const point& c0 = rawHalf0Ctrs[from1To0[i]];
+                    const point& c0 = half0Ctrs[from1To0[i]];
+                    const point& c1 = half1Ctrs[i];
+                    writeOBJ(ccStr, c0, c1, vertI);
+                }
+            }
+        }
+    }
+
+
+    // 3. Baffles(coincident faces) converted into cyclics (e.g. jump)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (!matchedAll)
+    {
+        label baffleI = 0;
+
+        forAll(pp, faceI)
+        {
+            const face& f = pp.localFaces()[faceI];
+            const labelList& pFaces = pp.pointFaces()[f[0]];
+
+            label matchedFaceI = -1;
+
+            forAll(pFaces, i)
+            {
+                label otherFaceI = pFaces[i];
+
+                if (otherFaceI > faceI)
+                {
+                    const face& otherF = pp.localFaces()[otherFaceI];
+
+                    // Note: might pick up two similar oriented faces
+                    //       (but that is illegal anyway)
+                    if (f == otherF)
+                    {
+                        matchedFaceI = otherFaceI;
+                        break;
+                    }
+                }
+            }
+
+            if (matchedFaceI != -1)
+            {
+                half0ToPatch[baffleI] = faceI;
+                half1ToPatch[baffleI] = matchedFaceI;
+                baffleI++;
+            }
+        }
+
+        if (baffleI == halfSize)
+        {
+            // And redo all matching
+            half0Faces = UIndirectList<face>(pp, half0ToPatch);
+            half1Faces = UIndirectList<face>(pp, half1ToPatch);
+
+            getCentresAndAnchors
+            (
+                pp,
+                half0Faces,
+                half1Faces,
+
+                ppPoints,
+                half0Ctrs,
+                half1Ctrs,
+                anchors0,
+                tols
+            );
+
+            // Geometric match of face centre vectors
+            matchedAll = matchPoints
+            (
+                half1Ctrs,
+                half0Ctrs,
+                tols,
+                false,
+                from1To0
+            );
+
+            if (debug)
+            {
+                Pout<< "oldCyclicPolyPatch::order : test if baffles:"
+                    << matchedAll << endl;
+                // Dump halves
+                fileName nm0("match3_"+name()+"_half0_faces.obj");
+                Pout<< "oldCyclicPolyPatch::order : Writing half0"
+                    << " faces to OBJ file " << nm0 << endl;
+                writeOBJ(nm0, half0Faces, ppPoints);
+
+                fileName nm1("match3_"+name()+"_half1_faces.obj");
+                Pout<< "oldCyclicPolyPatch::order : Writing half1"
+                    << " faces to OBJ file " << nm1 << endl;
+                writeOBJ(nm1, half1Faces, pp.points());
+
+                OFstream ccStr
+                (
+                    boundaryMesh().mesh().time().path()
+                   /"match3_"+ name() + "_faceCentres.obj"
+                );
+                Pout<< "oldCyclicPolyPatch::order : "
+                    << "Dumping currently found cyclic match as lines between"
+                    << " corresponding face centres to file " << ccStr.name()
+                    << endl;
+
+                // Recalculate untransformed face centres
+                //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+                label vertI = 0;
+
+                forAll(half1Ctrs, i)
+                {
+                    if (from1To0[i] != -1)
+                    {
+                        // Write edge between c1 and c0
+                        //const point& c0 = rawHalf0Ctrs[from1To0[i]];
+                        const point& c0 = half0Ctrs[from1To0[i]];
+                        const point& c1 = half1Ctrs[i];
+                        writeOBJ(ccStr, c0, c1, vertI);
+                    }
+                }
+            }
+        }
+    }
+
+
+    // 4. Automatic geometric ordering
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    if (!matchedAll)
+    {
+        // Split faces according to feature angle or topology
+        bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
+
+        if (!okSplit)
+        {
+            // Did not split into two equal parts.
+            return false;
+        }
+
+        // And redo all matching
+        half0Faces = UIndirectList<face>(pp, half0ToPatch);
+        half1Faces = UIndirectList<face>(pp, half1ToPatch);
+
+        getCentresAndAnchors
+        (
+            pp,
+            half0Faces,
+            half1Faces,
+
+            ppPoints,
+            half0Ctrs,
+            half1Ctrs,
+            anchors0,
+            tols
+        );
+
+        // Geometric match of face centre vectors
+        matchedAll = matchPoints
+        (
+            half1Ctrs,
+            half0Ctrs,
+            tols,
+            false,
+            from1To0
+        );
+
+        if (debug)
+        {
+            Pout<< "oldCyclicPolyPatch::order : automatic ordering result:"
+                << matchedAll << endl;
+            // Dump halves
+            fileName nm0("match4_"+name()+"_half0_faces.obj");
+            Pout<< "oldCyclicPolyPatch::order : Writing half0"
+                << " faces to OBJ file " << nm0 << endl;
+            writeOBJ(nm0, half0Faces, ppPoints);
+
+            fileName nm1("match4_"+name()+"_half1_faces.obj");
+            Pout<< "oldCyclicPolyPatch::order : Writing half1"
+                << " faces to OBJ file " << nm1 << endl;
+            writeOBJ(nm1, half1Faces, pp.points());
+
+            OFstream ccStr
+            (
+                boundaryMesh().mesh().time().path()
+               /"match4_"+ name() + "_faceCentres.obj"
+            );
+            Pout<< "oldCyclicPolyPatch::order : "
+                << "Dumping currently found cyclic match as lines between"
+                << " corresponding face centres to file " << ccStr.name()
+                << endl;
+
+            // Recalculate untransformed face centres
+            //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+            label vertI = 0;
+
+            forAll(half1Ctrs, i)
+            {
+                if (from1To0[i] != -1)
+                {
+                    // Write edge between c1 and c0
+                    //const point& c0 = rawHalf0Ctrs[from1To0[i]];
+                    const point& c0 = half0Ctrs[from1To0[i]];
+                    const point& c1 = half1Ctrs[i];
+                    writeOBJ(ccStr, c0, c1, vertI);
+                }
+            }
+        }
+    }
+
+
+    if (!matchedAll || debug)
+    {
+        // Dump halves
+        fileName nm0(name()+"_half0_faces.obj");
+        Pout<< "oldCyclicPolyPatch::order : Writing half0"
+            << " faces to OBJ file " << nm0 << endl;
+        writeOBJ(nm0, half0Faces, pp.points());
+
+        fileName nm1(name()+"_half1_faces.obj");
+        Pout<< "oldCyclicPolyPatch::order : Writing half1"
+            << " faces to OBJ file " << nm1 << endl;
+        writeOBJ(nm1, half1Faces, pp.points());
+
+        OFstream ccStr
+        (
+            boundaryMesh().mesh().time().path()
+           /name() + "_faceCentres.obj"
+        );
+        Pout<< "oldCyclicPolyPatch::order : "
+            << "Dumping currently found cyclic match as lines between"
+            << " corresponding face centres to file " << ccStr.name()
+            << endl;
+
+        // Recalculate untransformed face centres
+        //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+        label vertI = 0;
+
+        forAll(half1Ctrs, i)
+        {
+            if (from1To0[i] != -1)
+            {
+                // Write edge between c1 and c0
+                //const point& c0 = rawHalf0Ctrs[from1To0[i]];
+                const point& c0 = half0Ctrs[from1To0[i]];
+                const point& c1 = half1Ctrs[i];
+                writeOBJ(ccStr, c0, c1, vertI);
+            }
+        }
+    }
+
+
+    if (!matchedAll)
+    {
+        SeriousErrorIn
+        (
+            "oldCyclicPolyPatch::order"
+            "(const primitivePatch&, labelList&, labelList&) const"
+        )   << "Patch:" << name() << " : "
+            << "Cannot match vectors to faces on both sides of patch" << endl
+            << "    Perhaps your faces do not match?"
+            << " The obj files written contain the current match." << endl
+            << "    Continuing with incorrect face ordering from now on!"
+            << endl;
+
+            return false;
+    }
+
+
+    // Set faceMap such that half0 faces get first and corresponding half1
+    // faces last.
+    matchAnchors
+    (
+        true,                   // report if anchor matching error
+        pp,
+        half0ToPatch,
+        anchors0,
+        half1ToPatch,
+        half1Faces,
+        from1To0,
+        tols,
+        faceMap,
+        rotation
+    );
+
+    // Return false if no change neccesary, true otherwise.
+
+    forAll(faceMap, faceI)
+    {
+        if (faceMap[faceI] != faceI || rotation[faceI] != 0)
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+void Foam::oldCyclicPolyPatch::write(Ostream& os) const
+{
+    // Replacement of polyPatch::write to write 'cyclic' instead of type():
+    os.writeKeyword("type") << cyclicPolyPatch::typeName
+        << token::END_STATEMENT << nl;
+    patchIdentifier::write(os);
+    os.writeKeyword("nFaces") << size() << token::END_STATEMENT << nl;
+    os.writeKeyword("startFace") << start() << token::END_STATEMENT << nl;
+
+
+    os.writeKeyword("featureCos") << featureCos_ << token::END_STATEMENT << nl;
+    switch (transform_)
+    {
+        case ROTATIONAL:
+        {
+            os.writeKeyword("transform") << transformTypeNames[transform_]
+                << token::END_STATEMENT << nl;
+            os.writeKeyword("rotationAxis") << rotationAxis_
+                << token::END_STATEMENT << nl;
+            os.writeKeyword("rotationCentre") << rotationCentre_
+                << token::END_STATEMENT << nl;
+            break;
+        }
+        case TRANSLATIONAL:
+        {
+            os.writeKeyword("transform") << transformTypeNames[transform_]
+                << token::END_STATEMENT << nl;
+            os.writeKeyword("separationVector") << separationVector_
+                << token::END_STATEMENT << nl;
+            break;
+        }
+        default:
+        {
+            // no additional info to write
+        }
+    }
+
+    WarningIn("oldCyclicPolyPatch::write(Ostream& os) const")
+        << "Please run foamUpgradeCyclics to convert these old-style"
+        << " cyclics into two separate cyclics patches."
+        << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H
new file mode 100644
index 0000000000..f8cfb169c0
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H
@@ -0,0 +1,318 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::oldCyclicPolyPatch
+
+Description
+    'old' style cyclic polyPatch with all faces in single patch. Does ordering
+    but cannot be used to run. Writes 'type cyclic' so foamUpgradeCyclics
+    can be run afterwards.
+    Used to get cyclics from mesh converters that assume cyclics in single
+    patch (e.g. fluent3DMeshToFoam)
+
+SourceFiles
+    oldCyclicPolyPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oldCyclicPolyPatch_H
+#define oldCyclicPolyPatch_H
+
+#include "coupledPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class oldCyclicPolyPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class oldCyclicPolyPatch
+:
+    public coupledPolyPatch
+{
+public:
+
+    enum transformType
+    {
+        UNKNOWN,
+        ROTATIONAL,
+        TRANSLATIONAL
+    };
+    static const NamedEnum<transformType, 3> transformTypeNames;
+
+
+private:
+
+    // Private data
+
+        //- Morph:angle between normals of neighbouring faces.
+        //  Used to split cyclic into halves.
+        scalar featureCos_;
+
+        //- Type of transformation - rotational or translational
+        transformType transform_;
+
+        // For rotation
+
+            //- Axis of rotation for rotational cyclics
+            vector rotationAxis_;
+
+            //- point on axis of rotation for rotational cyclics
+            point rotationCentre_;
+
+        // For translation
+
+            //- Translation vector
+            vector separationVector_;
+
+
+    // Private member functions
+
+        //- Find amongst selected faces the one with the largest area
+        static label findMaxArea(const pointField&, const faceList&);
+
+        void calcTransforms();
+
+        //- Calculate face centres
+        static pointField calcFaceCentres
+        (
+            const UList<face>&,
+            const pointField&
+        );
+
+        //- Get f[0] for all faces
+        static pointField getAnchorPoints
+        (
+            const UList<face>&,
+            const pointField&
+        );
+
+        // Face ordering
+
+            //- Find the two parts of the faces of pp using feature edges.
+            //  Returns true if successfull.
+            bool getGeometricHalves
+            (
+                const primitivePatch&,
+                labelList&,
+                labelList&
+            ) const;
+
+            //- Calculate geometric factors of the two halves.
+            void getCentresAndAnchors
+            (
+                const primitivePatch&,
+                const faceList& half0Faces,
+                const faceList& half1Faces,
+
+                pointField& ppPoints,
+                pointField& half0Ctrs,
+                pointField& half1Ctrs,
+                pointField& anchors0,
+                scalarField& tols
+            ) const;
+
+            //- Given matched faces matches the anchor point. Sets faceMap,
+            //  rotation. Returns true if all matched.
+            bool matchAnchors
+            (
+                const bool report,
+                const primitivePatch&,
+                const labelList&,
+                const pointField&,
+                const labelList&,
+                const faceList&,
+                const labelList&,
+                const scalarField&,
+
+                labelList& faceMap,
+                labelList& rotation
+            ) const;
+
+            //- For rotational cases, try to find a unique face on each side
+            //  of the cyclic.
+            label getConsistentRotationFace
+            (
+                const pointField& faceCentres
+            ) const;
+
+
+protected:
+
+    // Protected Member functions
+
+        //- Initialise the calculation of the patch geometry
+        virtual void initGeometry(PstreamBuffers&);
+
+        //- Calculate the patch geometry
+        virtual void calcGeometry(PstreamBuffers&);
+
+        //- Initialise the patches for moving points
+        virtual void initMovePoints(PstreamBuffers&, const pointField&);
+
+        //- Correct patches after moving points
+        virtual void movePoints(PstreamBuffers&, const pointField&);
+
+        //- Initialise the update of the patch topology
+        virtual void initUpdateMesh(PstreamBuffers&);
+
+        //- Update of the patch topology
+        virtual void updateMesh(PstreamBuffers&);
+
+public:
+
+    //- Runtime type information
+    TypeName("oldCyclic");
+
+
+    // Constructors
+
+        //- Construct from components
+        oldCyclicPolyPatch
+        (
+            const word& name,
+            const label size,
+            const label start,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct from dictionary
+        oldCyclicPolyPatch
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const polyBoundaryMesh& bm
+        );
+
+        //- Construct as copy, resetting the boundary mesh
+        oldCyclicPolyPatch(const oldCyclicPolyPatch&, const polyBoundaryMesh&);
+
+        //- Construct given the original patch and resetting the
+        //  face list and boundary mesh information
+        oldCyclicPolyPatch
+        (
+            const oldCyclicPolyPatch& pp,
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        );
+
+        //- Construct and return a clone, resetting the boundary mesh
+        virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const
+        {
+            return autoPtr<polyPatch>(new oldCyclicPolyPatch(*this, bm));
+        }
+
+        //- Construct and return a clone, resetting the face list
+        //  and boundary mesh
+        virtual autoPtr<polyPatch> clone
+        (
+            const polyBoundaryMesh& bm,
+            const label index,
+            const label newSize,
+            const label newStart
+        ) const
+        {
+            return autoPtr<polyPatch>
+            (
+                new oldCyclicPolyPatch(*this, bm, index, newSize, newStart)
+            );
+        }
+
+
+    // Destructor
+
+        virtual ~oldCyclicPolyPatch();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Does this side own the patch ?
+            virtual bool owner() const
+            {
+                notImplemented("oldCyclicPolyPatch::owner()");
+                return true;
+            }
+
+            //- Transform a patch-based position from other side to this side
+            virtual void transformPosition(pointField& l) const
+            {
+                notImplemented("transformPosition(pointField&)");
+            }
+
+        //- Calculate the patch geometry
+        virtual void calcGeometry
+        (
+            const primitivePatch& referPatch,
+            const UList<point>& thisCtrs,
+            const UList<point>& thisAreas,
+            const UList<point>& thisCc,
+            const UList<point>& nbrCtrs,
+            const UList<point>& nbrAreas,
+            const UList<point>& nbrCc
+        );
+
+        //- Initialize ordering for primitivePatch. Does not
+        //  refer to *this (except for name() and type() etc.)
+        virtual void initOrder
+        (
+            PstreamBuffers&,
+            const primitivePatch&
+        ) const;
+
+        //- Return new ordering for primitivePatch.
+        //  Ordering is -faceMap: for every face
+        //  index of the new face -rotation:for every new face the clockwise
+        //  shift of the original face. Return false if nothing changes
+        //  (faceMap is identity, rotation is 0), true otherwise.
+        virtual bool order
+        (
+            PstreamBuffers&,
+            const primitivePatch&,
+            labelList& faceMap,
+            labelList& rotation
+        ) const;
+
+        //- Write the polyPatch data as a dictionary
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

From 287865823c402d086f1b6615bdb72333c276924c Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 22 Oct 2010 11:51:29 +0100
Subject: [PATCH 13/17] ENH: ptscotchDecomp : added bit pseudo code for writing
 graph file

---
 .../decompose/ptscotchDecomp/ptscotchDecomp.C | 30 ++++++++++++++-----
 1 file changed, 22 insertions(+), 8 deletions(-)

diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
index 98f231cf48..ff262c5645 100644
--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
@@ -180,17 +180,31 @@ Foam::label Foam::ptscotchDecomp::decompose
 //            Info<< "Dumping Scotch graph file to " << str.name() << endl
 //                << "Use this in combination with gpart." << endl;
 //
-//            label version = 0;
+//            // Distributed graph file (.grf)
+//            label version = 2;
 //            str << version << nl;
-//            // Numer of vertices
-//            str << xadj.size()-1 << ' ' << adjncy.size() << nl;
+//            // Number of files
+
+//            // Number of files (procglbnbr)
+//            str << Pstream::nProcs();
+//            // My file number (procloc)
+//            str << ' ' << Pstream::myProcNo() << nl;
+//
+//            // Total number of vertices (vertglbnbr)
+//            str << returnReduce(mesh.nCells(), sumOp<label>());
+//            // Total number of connections (edgeglbnbr)
+//            str << ' ' << returnReduce(xadj[mesh.nCells()], sumOp<label>())
+//                << nl;
+//            // Local number of vertices (vertlocnbr)
+//            str << mesh.nCells();
+//            // Local number of connections (edgelocnbr)
+//            str << ' ' << xadj[mesh.nCells()] << nl;
 //            // Numbering starts from 0
 //            label baseval = 0;
-//            // Has weights?
-//            label hasEdgeWeights = 0;
-//            label hasVertexWeights = 0;
-//            label numericflag = 10*hasEdgeWeights+hasVertexWeights;
-//            str << baseval << ' ' << numericflag << nl;
+//            // Start of my global vertices (procdsptab[proclocnum])
+//            str << ' ' << globalCells.toGlobal(0);
+//            100*hasVertlabels+10*hasEdgeWeights+1*hasVertWeighs
+//            str << ' ' << "0" << nl;
 //            for (label cellI = 0; cellI < xadj.size()-1; cellI++)
 //            {
 //                label start = xadj[cellI];

From 511567be3804135acebdf84f63c315398da27e52 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 22 Oct 2010 12:19:34 +0100
Subject: [PATCH 14/17] STYLE: cyclicPolyPatch : 80 char limit

---
 .../constraint/oldCyclic/oldCyclicPolyPatch.C            | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
index 2dcbfef30e..218debfa5b 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C
@@ -361,7 +361,8 @@ void Foam::oldCyclicPolyPatch::getCentresAndAnchors
         //    if (debug)
         //    {
         //        Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
-        //            << "Specified translation : " << separationVector_ << endl;
+        //            << "Specified translation : " << separationVector_
+        //            << endl;
         //    }
         //
         //    half0Ctrs += separationVector_;
@@ -946,7 +947,6 @@ bool Foam::oldCyclicPolyPatch::order
                 << endl;
 
             // Recalculate untransformed face centres
-            //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
             label vertI = 0;
 
             forAll(half1Ctrs, i)
@@ -954,7 +954,6 @@ bool Foam::oldCyclicPolyPatch::order
                 if (from1To0[i] != -1)
                 {
                     // Write edge between c1 and c0
-                    //const point& c0 = rawHalf0Ctrs[from1To0[i]];
                     const point& c0 = half0Ctrs[from1To0[i]];
                     const point& c1 = half1Ctrs[i];
                     writeOBJ(ccStr, c0, c1, vertI);
@@ -1059,7 +1058,6 @@ bool Foam::oldCyclicPolyPatch::order
                     << endl;
 
                 // Recalculate untransformed face centres
-                //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
                 label vertI = 0;
 
                 forAll(half1Ctrs, i)
@@ -1067,7 +1065,6 @@ bool Foam::oldCyclicPolyPatch::order
                     if (from1To0[i] != -1)
                     {
                         // Write edge between c1 and c0
-                        //const point& c0 = rawHalf0Ctrs[from1To0[i]];
                         const point& c0 = half0Ctrs[from1To0[i]];
                         const point& c1 = half1Ctrs[i];
                         writeOBJ(ccStr, c0, c1, vertI);
@@ -1145,7 +1142,6 @@ bool Foam::oldCyclicPolyPatch::order
                 << endl;
 
             // Recalculate untransformed face centres
-            //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
             label vertI = 0;
 
             forAll(half1Ctrs, i)
@@ -1153,7 +1149,6 @@ bool Foam::oldCyclicPolyPatch::order
                 if (from1To0[i] != -1)
                 {
                     // Write edge between c1 and c0
-                    //const point& c0 = rawHalf0Ctrs[from1To0[i]];
                     const point& c0 = half0Ctrs[from1To0[i]];
                     const point& c1 = half1Ctrs[i];
                     writeOBJ(ccStr, c0, c1, vertI);

From 78ee5c8c6c85fde1ffe9f5197007506826b3336d Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 22 Oct 2010 12:42:32 +0100
Subject: [PATCH 15/17] ENH: libgenericPatchFields for subsetMesh.

---
 .../utilities/mesh/manipulation/subsetMesh/Make/options        | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/applications/utilities/mesh/manipulation/subsetMesh/Make/options b/applications/utilities/mesh/manipulation/subsetMesh/Make/options
index d27c95d033..969020c4af 100644
--- a/applications/utilities/mesh/manipulation/subsetMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/subsetMesh/Make/options
@@ -4,4 +4,5 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools \
+    -lgenericPatchFields

From a67a967d9ac3a2066d06740f2924ccab8c9b7751 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 25 Oct 2010 17:49:04 +0100
Subject: [PATCH 16/17] finiteVolume: Added fvcSmooth which provides field
 smoothing operations: smooth, sweep and spread

---
 src/finiteVolume/Make/files                   |   1 +
 .../finiteVolume/fvc/fvcSmooth/fvcSmooth.C    | 318 ++++++++++++++++++
 .../finiteVolume/fvc/fvcSmooth/fvcSmooth.H    |  73 ++++
 .../finiteVolume/fvc/fvcSmooth/smoothData.H   | 204 +++++++++++
 .../finiteVolume/fvc/fvcSmooth/smoothDataI.H  | 192 +++++++++++
 .../finiteVolume/fvc/fvcSmooth/sweepData.H    | 205 +++++++++++
 .../finiteVolume/fvc/fvcSmooth/sweepDataI.H   | 208 ++++++++++++
 7 files changed, 1201 insertions(+)
 create mode 100644 src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C
 create mode 100644 src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H
 create mode 100644 src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H
 create mode 100644 src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H
 create mode 100644 src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H
 create mode 100644 src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H

diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 21b90fa409..dc325c1ddf 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -337,6 +337,7 @@ $(laplacianSchemes)/laplacianScheme/laplacianSchemes.C
 $(laplacianSchemes)/gaussLaplacianScheme/gaussLaplacianSchemes.C
 
 finiteVolume/fvc/fvcMeshPhi.C
+finiteVolume/fvc/fvcSmooth/fvcSmooth.C
 
 general = cfdTools/general
 $(general)/findRefCell/findRefCell.C
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C
new file mode 100644
index 0000000000..c04d16dcc9
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C
@@ -0,0 +1,318 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvcSmooth.H"
+#include "volFields.H"
+#include "FaceCellWave.H"
+#include "smoothData.H"
+#include "sweepData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::smoothData::maxRatio = 1.0;
+
+void Foam::fvc::smooth
+(
+    volScalarField& field,
+    const scalar coeff
+)
+{
+    const fvMesh& mesh = field.mesh();
+    smoothData::maxRatio = 1 + coeff;
+
+    DynamicList<label> changedFaces(mesh.nFaces()/100 + 100);
+    DynamicList<smoothData> changedFacesInfo(changedFaces.size());
+
+    const unallocLabelList& owner = mesh.owner();
+    const unallocLabelList& neighbour = mesh.neighbour();
+
+    forAll(owner, facei)
+    {
+        const label own = owner[facei];
+        const label nbr = neighbour[facei];
+
+        // Check if owner value much larger than neighbour value or vice versa
+        if (field[own] > smoothData::maxRatio*field[nbr])
+        {
+            changedFaces.append(facei);
+            changedFacesInfo.append(smoothData(field[own]));
+        }
+        else if (field[nbr] > smoothData::maxRatio*field[own])
+        {
+            changedFaces.append(facei);
+            changedFacesInfo.append(smoothData(field[nbr]));
+        }
+    }
+
+    // Insert all faces of coupled patches - FaceCellWave will correct them
+    forAll(mesh.boundaryMesh(), patchi)
+    {
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+        if (patch.coupled())
+        {
+            forAll(patch, patchFacei)
+            {
+                label facei = patch.start() + patchFacei;
+                label own = mesh.faceOwner()[facei];
+
+                changedFaces.append(facei);
+                changedFacesInfo.append(smoothData(field[own]));
+            }
+        }
+    }
+
+    changedFaces.shrink();
+    changedFacesInfo.shrink();
+
+    // Set initial field on cells
+    List<smoothData> cellData(mesh.nCells());
+
+    forAll(field, celli)
+    {
+        cellData[celli] = field[celli];
+    }
+
+    // Set initial field on faces
+    List<smoothData> faceData(mesh.nFaces());
+
+
+    // Propagate information over whole domain
+    FaceCellWave<smoothData > smoothData
+    (
+        mesh,
+        changedFaces,
+        changedFacesInfo,
+        faceData,
+        cellData,
+        mesh.globalData().nTotalCells()    // max iterations
+    );
+
+    forAll(field, celli)
+    {
+        field[celli] = cellData[celli].value();
+    }
+
+    field.correctBoundaryConditions();
+}
+
+
+void Foam::fvc::spread
+(
+    volScalarField& field,
+    const volScalarField& alpha,
+    const label nLayers
+)
+{
+    const fvMesh& mesh = field.mesh();
+    smoothData::maxRatio = 1;
+
+    DynamicList<label> changedFaces(mesh.nFaces()/100 + 100);
+    DynamicList<smoothData> changedFacesInfo(changedFaces.size());
+
+    // Set initial field on cells
+    List<smoothData> cellData(mesh.nCells());
+
+    forAll(field, celli)
+    {
+        cellData[celli] = field[celli];
+    }
+
+    // Set initial field on faces
+    List<smoothData> faceData(mesh.nFaces());
+
+    const unallocLabelList& owner = mesh.owner();
+    const unallocLabelList& neighbour = mesh.neighbour();
+
+    forAll(owner, facei)
+    {
+        const label own = owner[facei];
+        const label nbr = neighbour[facei];
+
+        if
+        (
+            (alpha[own] > 0.01 && alpha[own] < 0.99)
+         || (alpha[nbr] > 0.01 && alpha[nbr] < 0.99)
+        )
+        {
+            if (mag(alpha[own] - alpha[nbr]) > 0.2)
+            {
+                changedFaces.append(facei);
+                changedFacesInfo.append
+                (
+                    smoothData(max(field[own], field[nbr]))
+                );
+            }
+        }
+    }
+
+    // Insert all faces of coupled patches - FaceCellWave will correct them
+    forAll(mesh.boundaryMesh(), patchi)
+    {
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+        if (patch.coupled())
+        {
+            forAll(patch, patchFacei)
+            {
+                label facei = patch.start() + patchFacei;
+                label own = mesh.faceOwner()[facei];
+
+                scalarField alphapn =
+                    alpha.boundaryField()[patchi].patchNeighbourField();
+
+                if
+                (
+                    (alpha[own] > 0.01 && alpha[own] < 0.99)
+                 || (alphapn[patchFacei] > 0.01 && alphapn[patchFacei] < 0.99)
+                )
+                {
+                    if (mag(alpha[own] - alphapn[patchFacei]) > 0.2)
+                    {
+                        changedFaces.append(facei);
+                        changedFacesInfo.append(smoothData(field[own]));
+                    }
+                }
+            }
+        }
+    }
+
+    changedFaces.shrink();
+    changedFacesInfo.shrink();
+
+    // Propagate information over whole domain
+    FaceCellWave<smoothData> smoothData
+    (
+        mesh,
+        faceData,
+        cellData
+    );
+
+    smoothData.setFaceInfo(changedFaces, changedFacesInfo);
+
+    smoothData.iterate(nLayers);
+
+    forAll(field, celli)
+    {
+        field[celli] = cellData[celli].value();
+    }
+
+    field.correctBoundaryConditions();
+}
+
+
+void Foam::fvc::sweep
+(
+    volScalarField& field,
+    const volScalarField& alpha,
+    const label nLayers
+)
+{
+    const fvMesh& mesh = field.mesh();
+
+    DynamicList<label> changedFaces(mesh.nFaces()/100 + 100);
+    DynamicList<sweepData> changedFacesInfo(changedFaces.size());
+
+    // Set initial field on cells
+    List<sweepData> cellData(mesh.nCells());
+
+    // Set initial field on faces
+    List<sweepData> faceData(mesh.nFaces());
+
+    const unallocLabelList& owner = mesh.owner();
+    const unallocLabelList& neighbour = mesh.neighbour();
+    const vectorField& Cf = mesh.faceCentres();
+
+    forAll(owner, facei)
+    {
+        const label own = owner[facei];
+        const label nbr = neighbour[facei];
+
+        if (mag(alpha[own] - alpha[nbr]) > 0.2)
+        {
+            changedFaces.append(facei);
+            changedFacesInfo.append
+            (
+                sweepData(max(field[own], field[nbr]),
+                Cf[facei])
+            );
+        }
+    }
+
+    // Insert all faces of coupled patches - FaceCellWave will correct them
+    forAll(mesh.boundaryMesh(), patchi)
+    {
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+        if (patch.coupled())
+        {
+            forAll(patch, patchFacei)
+            {
+                label facei = patch.start() + patchFacei;
+                label own = mesh.faceOwner()[facei];
+
+                scalarField alphapn =
+                    alpha.boundaryField()[patchi].patchNeighbourField();
+
+                if (mag(alpha[own] - alphapn[patchFacei]) > 0.2)
+                {
+                    changedFaces.append(facei);
+                    changedFacesInfo.append
+                    (
+                        sweepData(field[own], Cf[facei])
+                    );
+                }
+            }
+        }
+    }
+
+    changedFaces.shrink();
+    changedFacesInfo.shrink();
+
+    // Propagate information over whole domain
+    FaceCellWave<sweepData> sweepData
+    (
+        mesh,
+        faceData,
+        cellData
+    );
+
+    sweepData.setFaceInfo(changedFaces, changedFacesInfo);
+
+    sweepData.iterate(nLayers);
+
+    forAll(field, celli)
+    {
+        if (cellData[celli].valid())
+        {
+            field[celli] = max(field[celli], cellData[celli].value());
+        }
+    }
+
+    field.correctBoundaryConditions();
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H
new file mode 100644
index 0000000000..ced147e3dc
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H
@@ -0,0 +1,73 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+function
+    Foam::fvc::smooth
+
+Description
+    Function that uses smoothData to apply spatial smoothing of a
+    volume field using the FaceCellWave algorithm
+
+SourceFiles
+    fvcSmooth.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fvcSmooth_H
+#define fvcSmooth_H
+
+#include "volFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fvc
+{
+    void smooth
+    (
+        volScalarField& field,
+        const scalar coeff
+    );
+
+    void spread
+    (
+        volScalarField& field,
+        const volScalarField& alpha,
+        const label nLayers
+    );
+
+    void sweep
+    (
+        volScalarField& field,
+        const volScalarField& alpha,
+        const label nLayers
+    );
+}
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H
new file mode 100644
index 0000000000..b55c93a52a
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H
@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::smoothData
+
+Description
+    Helper class used by the smoothVolField class
+
+Files
+    smoothData.H
+    smoothDataI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef smoothData_H
+#define smoothData_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class smoothData Declaration
+\*---------------------------------------------------------------------------*/
+
+class smoothData
+{
+    scalar value_;
+
+    // Private Member Functions
+
+        //- Update - gets information from neighbouring face/cell and
+        //  uses this to update itself (if necessary) and return true
+        inline bool update
+        (
+            const smoothData& svf,
+            const scalar scale,
+            const scalar tol
+        );
+
+
+public:
+
+    // Static data members
+
+        //- Field fraction
+        static scalar maxRatio;
+
+
+    // Constructors
+
+        //- Construct null
+        inline smoothData();
+
+        //- Construct from inverse field value
+        inline smoothData(const scalar value);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return value
+            scalar value() const
+            {
+                return value_;
+            }
+
+
+        // Needed by FaceCellWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value
+            inline bool valid() const;
+
+            //- Check for identical geometrical data
+            //  Used for cyclics checking
+            inline bool sameGeometry
+            (
+                const polyMesh&,
+                const smoothData&,
+                const scalar
+            ) const;
+
+            //- Convert any absolute coordinates into relative to
+            //  (patch)face centre
+            inline void leaveDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Reverse of leaveDomain
+            inline void enterDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Apply rotation matrix to any coordinates
+            inline void transform
+            (
+                const polyMesh&,
+                const tensor&
+            );
+
+            //- Influence of neighbouring face
+            inline bool updateCell
+            (
+                const polyMesh&,
+                const label thisCellI,
+                const label neighbourFaceI,
+                const smoothData& svf,
+                const scalar tol
+            );
+
+            //- Influence of neighbouring cell
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const label neighbourCellI,
+                const smoothData& svf,
+                const scalar tol
+            );
+
+            //- Influence of different value on same face
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const smoothData& svf,
+                const scalar tol
+            );
+
+
+        // Member Operators
+
+            inline void operator=(const scalar value);
+
+            // Needed for List IO
+            inline bool operator==(const smoothData&) const;
+
+            inline bool operator!=(const smoothData&) const;
+
+
+        // IOstream Operators
+
+            friend Ostream& operator<<
+            (
+                Ostream& os,
+                const smoothData& svf
+            )
+            {
+                return os  << svf.value_;
+            }
+
+            friend Istream& operator>>(Istream& is, smoothData& svf)
+            {
+                return is  >> svf.value_;
+            }
+};
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "smoothDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H
new file mode 100644
index 0000000000..3faaa694b6
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H
@@ -0,0 +1,192 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline bool Foam::smoothData::update
+(
+    const smoothData::smoothData& svf,
+    const scalar scale,
+    const scalar tol
+)
+{
+    if (!valid() || (value_ < VSMALL))
+    {
+        // My value not set - take over neighbour
+        value_ = svf.value()/scale;
+
+        // Something changed - let caller know
+        return true;
+    }
+    else if (svf.value() > (1 + tol)*scale*value_)
+    {
+        // Neighbour is too big for me - Up my value
+        value_ = svf.value()/scale;
+
+        // Something changed - let caller know
+        return true;
+    }
+    else
+    {
+        // Neighbour is not too big for me or change is too small
+        // Nothing changed
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::smoothData::smoothData()
+:
+    value_(-GREAT)
+{}
+
+
+inline Foam::smoothData::smoothData(const scalar value)
+:
+    value_(value)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::smoothData::valid() const
+{
+    return value_ > -SMALL;
+}
+
+
+inline bool Foam::smoothData::sameGeometry
+(
+    const polyMesh&,
+    const smoothData&,
+    const scalar
+) const
+{
+    return true;
+}
+
+
+inline void Foam::smoothData::leaveDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point&
+)
+{}
+
+
+inline void Foam::smoothData::transform
+(
+    const polyMesh&,
+    const tensor&
+)
+{}
+
+
+inline void Foam::smoothData::enterDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point&
+)
+{}
+
+
+inline bool Foam::smoothData::updateCell
+(
+    const polyMesh&,
+    const label,
+    const label,
+    const smoothData& svf,
+    const scalar tol
+)
+{
+    // Take over info from face if more than deltaRatio larger
+    return update(svf, maxRatio, tol);
+}
+
+
+inline bool Foam::smoothData::updateFace
+(
+    const polyMesh&,
+    const label,
+    const label,
+    const smoothData& svf,
+    const scalar tol
+)
+{
+    // Take over information from cell without any scaling (scale = 1.0)
+    return update(svf, 1.0, tol);
+}
+
+
+// Update this (face) with coupled face information.
+inline bool Foam::smoothData::updateFace
+(
+    const polyMesh&,
+    const label,
+    const smoothData& svf,
+    const scalar tol
+)
+{
+    // Take over information from coupled face without any scaling (scale = 1.0)
+    return update(svf, 1.0, tol);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline void Foam::smoothData::operator=
+(
+    const scalar value
+)
+{
+    value_ = value;
+}
+
+
+inline bool Foam::smoothData::operator==
+(
+    const smoothData& rhs
+) const
+{
+    return value_ == rhs.value();
+}
+
+
+inline bool Foam::smoothData::operator!=
+(
+    const smoothData& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H
new file mode 100644
index 0000000000..ebee81c0f7
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H
@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sweepData
+
+Description
+    Helper class used by fvcSmooth
+
+Files
+    sweepData.H
+    sweepDataI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sweepData_H
+#define sweepData_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class sweepData Declaration
+\*---------------------------------------------------------------------------*/
+
+class sweepData
+{
+    scalar value_;
+    point origin_;
+
+    // Private Member Functions
+
+        //- Update - gets information from neighbouring face/cell and
+        //  uses this to update itself (if necessary) and return true
+        inline bool update
+        (
+            const sweepData& svf,
+            const point& position,
+            const scalar tol
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline sweepData();
+
+        //- Construct from component
+        inline sweepData(const scalar value, const point& origin);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return value
+            scalar value() const
+            {
+                return value_;
+            }
+
+            //- Return origin
+            const point& origin() const
+            {
+                return origin_;
+            }
+
+
+        // Needed by FaceCellWave
+
+            //- Check whether origin has been changed at all or
+            //  still contains original (invalid) value
+            inline bool valid() const;
+
+            //- Check for identical geometrical data
+            //  Used for cyclics checking
+            inline bool sameGeometry
+            (
+                const polyMesh&,
+                const sweepData&,
+                const scalar
+            ) const;
+
+            //- Convert any absolute coordinates into relative to
+            //  (patch)face centre
+            inline void leaveDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Reverse of leaveDomain
+            inline void enterDomain
+            (
+                const polyMesh&,
+                const polyPatch&,
+                const label patchFaceI,
+                const point& faceCentre
+            );
+
+            //- Apply rotation matrix to any coordinates
+            inline void transform
+            (
+                const polyMesh&,
+                const tensor&
+            );
+
+            //- Influence of neighbouring face
+            inline bool updateCell
+            (
+                const polyMesh&,
+                const label thisCellI,
+                const label neighbourFaceI,
+                const sweepData& svf,
+                const scalar tol
+            );
+
+            //- Influence of neighbouring cell
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const label neighbourCellI,
+                const sweepData& svf,
+                const scalar tol
+            );
+
+            //- Influence of different value on same face
+            inline bool updateFace
+            (
+                const polyMesh&,
+                const label thisFaceI,
+                const sweepData& svf,
+                const scalar tol
+            );
+
+
+        // Member Operators
+
+            inline void operator=(const scalar value);
+
+            // Needed for List IO
+            inline bool operator==(const sweepData&) const;
+
+            inline bool operator!=(const sweepData&) const;
+
+
+        // IOstream Operators
+
+            friend Ostream& operator<<
+            (
+                Ostream& os,
+                const sweepData& svf
+            )
+            {
+                return os << svf.value_ << svf.origin_;
+            }
+
+            friend Istream& operator>>(Istream& is, sweepData& svf)
+            {
+                return is >> svf.value_ >> svf.origin_;
+            }
+};
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "sweepDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H
new file mode 100644
index 0000000000..651c5e4e52
--- /dev/null
+++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H
@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "transform.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+inline bool Foam::sweepData::update
+(
+    const sweepData::sweepData& svf,
+    const point& position,
+    const scalar tol
+)
+{
+    if (!valid())
+    {
+        operator=(svf);
+        return true;
+    }
+
+    scalar myDist2 = magSqr(position - origin());
+
+    if (myDist2 < SMALL)
+    {
+        if (svf.value() > value())
+        {
+            operator=(svf);
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+
+    scalar dist2 = magSqr(position - svf.origin());
+
+    if (dist2 < myDist2)
+    {
+        operator=(svf);
+        return true;
+    }
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::sweepData::sweepData()
+:
+    value_(-GREAT),
+    origin_(vector::max)
+{}
+
+
+inline Foam::sweepData::sweepData(const scalar value, const point& origin)
+:
+    value_(value),
+    origin_(origin)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::sweepData::valid() const
+{
+    return value_ > -SMALL;
+}
+
+
+inline bool Foam::sweepData::sameGeometry
+(
+    const polyMesh&,
+    const sweepData&,
+    const scalar
+) const
+{
+    return true;
+}
+
+
+inline void Foam::sweepData::leaveDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point& faceCentre
+)
+{
+    origin_ -= faceCentre;
+}
+
+
+inline void Foam::sweepData::transform
+(
+    const polyMesh&,
+    const tensor& rotTensor
+)
+{
+    origin_ = Foam::transform(rotTensor, origin_);
+}
+
+
+inline void Foam::sweepData::enterDomain
+(
+    const polyMesh&,
+    const polyPatch&,
+    const label,
+    const point& faceCentre
+)
+{
+    // back to absolute form
+    origin_ += faceCentre;
+}
+
+
+inline bool Foam::sweepData::updateCell
+(
+    const polyMesh& mesh,
+    const label thisCellI,
+    const label,
+    const sweepData& svf,
+    const scalar tol
+)
+{
+    return update(svf, mesh.cellCentres()[thisCellI], tol);
+}
+
+
+inline bool Foam::sweepData::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFaceI,
+    const label,
+    const sweepData& svf,
+    const scalar tol
+)
+{
+    return update(svf, mesh.faceCentres()[thisFaceI], tol);
+}
+
+
+// Update this (face) with coupled face information.
+inline bool Foam::sweepData::updateFace
+(
+    const polyMesh& mesh,
+    const label thisFaceI,
+    const sweepData& svf,
+    const scalar tol
+)
+{
+    return update(svf, mesh.faceCentres()[thisFaceI], tol);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+inline void Foam::sweepData::operator=
+(
+    const scalar value
+)
+{
+    value_ = value;
+}
+
+
+inline bool Foam::sweepData::operator==
+(
+    const sweepData& rhs
+) const
+{
+    return origin() == rhs.origin();
+}
+
+
+inline bool Foam::sweepData::operator!=
+(
+    const sweepData& rhs
+) const
+{
+    return !(*this == rhs);
+}
+
+
+// ************************************************************************* //

From 730481aed3e00a66a8a85ecee4643824e42db5fd Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Mon, 25 Oct 2010 17:49:56 +0100
Subject: [PATCH 17/17] LTSInterFoam: Initial version of interFoam supporting
 local time-stepping for acceleration to steady-state

---
 .../solvers/multiphase/interFoam/Allwclean    |   1 +
 .../solvers/multiphase/interFoam/Allwmake     |   1 +
 .../interFoam/LTSInterFoam/LTSInterFoam.C     | 101 +++
 .../multiphase/interFoam/LTSInterFoam/MULES.C |  83 +++
 .../multiphase/interFoam/LTSInterFoam/MULES.H | 136 ++++
 .../interFoam/LTSInterFoam/MULESTemplates.C   | 602 ++++++++++++++++++
 .../interFoam/LTSInterFoam/Make/files         |   4 +
 .../interFoam/LTSInterFoam/Make/options       |  15 +
 .../interFoam/LTSInterFoam/alphaEqn.H         |  36 ++
 .../interFoam/LTSInterFoam/alphaEqnSubCycle.H |  35 +
 .../LTSInterFoam/setInitialrDeltaT.H          |  30 +
 .../interFoam/LTSInterFoam/setrDeltaT.H       | 109 ++++
 12 files changed, 1153 insertions(+)
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
 create mode 100644 applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

diff --git a/applications/solvers/multiphase/interFoam/Allwclean b/applications/solvers/multiphase/interFoam/Allwclean
index 350d4b268b..2c54bfa321 100755
--- a/applications/solvers/multiphase/interFoam/Allwclean
+++ b/applications/solvers/multiphase/interFoam/Allwclean
@@ -6,5 +6,6 @@ wclean
 wclean interDyMFoam
 wclean MRFInterFoam
 wclean porousInterFoam
+wclean LTSInterFoam
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/interFoam/Allwmake b/applications/solvers/multiphase/interFoam/Allwmake
index b25e4440e6..8044426582 100755
--- a/applications/solvers/multiphase/interFoam/Allwmake
+++ b/applications/solvers/multiphase/interFoam/Allwmake
@@ -6,5 +6,6 @@ wmake
 wmake interDyMFoam
 wmake MRFInterFoam
 wmake porousInterFoam
+wmake LTSInterFoam
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
new file mode 100644
index 0000000000..9fcb19a7cc
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C
@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    interFoam
+
+Description
+    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach.
+
+    The momentum and other fluid properties are of the "mixture" and a single
+    momentum equation is solved.
+
+    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
+
+    For a two-fluid approach see twoPhaseEulerFoam.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "MULES.H"
+#include "subCycle.H"
+#include "interfaceProperties.H"
+#include "twoPhaseMixture.H"
+#include "turbulenceModel.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readPISOControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+    #include "setInitialrDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readPISOControls.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "setrDeltaT.H"
+
+        twoPhaseProperties.correct();
+
+        #include "alphaEqnSubCycle.H"
+        turbulence->correct();
+        #include "UEqn.H"
+
+        // --- PISO loop
+        for (int corr=0; corr<nCorr; corr++)
+        {
+            #include "pEqn.H"
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
new file mode 100644
index 0000000000..e281ae575c
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C
@@ -0,0 +1,83 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MULES.H"
+#include "upwind.H"
+#include "uncorrectedSnGrad.H"
+#include "gaussConvectionScheme.H"
+#include "gaussLaplacianScheme.H"
+#include "uncorrectedSnGrad.H"
+#include "surfaceInterpolate.H"
+#include "fvcSurfaceIntegrate.H"
+#include "slicedSurfaceFields.H"
+#include "syncTools.H"
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::MULES::explicitLTSSolve
+(
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    explicitLTSSolve
+    (
+        geometricOneField(),
+        psi,
+        phi,
+        phiPsi,
+        zeroField(), zeroField(),
+        psiMax, psiMin
+    );
+}
+
+
+void Foam::MULES::implicitSolve
+(
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    implicitSolve
+    (
+        geometricOneField(),
+        psi,
+        phi,
+        phiPsi,
+        zeroField(), zeroField(),
+        psiMax, psiMin
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
new file mode 100644
index 0000000000..729b041aa6
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H
@@ -0,0 +1,136 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Global
+    MULES
+
+Description
+    MULES: Multidimensional universal limiter with explicit solution.
+
+    Solve a convective-only transport equation using an explicit universal
+    multi-dimensional limiter.
+
+    Parameters are the variable to solve, the normal convective flux and the
+    actual explicit flux of the variable which is also used to return limited
+    flux used in the bounded-solution.
+
+SourceFiles
+    MULES.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MULES_H
+#define MULES_H
+
+#include "volFields.H"
+#include "surfaceFieldsFwd.H"
+#include "primitiveFieldsFwd.H"
+#include "zeroField.H"
+#include "geometricOneField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace MULES
+{
+
+template<class RhoType, class SpType, class SuType>
+void explicitLTSSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+void explicitLTSSolve
+(
+    volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    surfaceScalarField& phiPsi,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+template<class RhoType, class SpType, class SuType>
+void implicitSolve
+(
+    const RhoType& rho,
+    volScalarField& gamma,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+void implicitSolve
+(
+    volScalarField& gamma,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiCorr,
+    const scalar psiMax,
+    const scalar psiMin
+);
+
+template<class RhoType, class SpType, class SuType>
+void limiter
+(
+    scalarField& allLambda,
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+);
+
+} // End namespace MULES
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "MULESTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
new file mode 100644
index 0000000000..a6d3a3f6b2
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C
@@ -0,0 +1,602 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MULES.H"
+#include "upwind.H"
+#include "uncorrectedSnGrad.H"
+#include "gaussConvectionScheme.H"
+#include "gaussLaplacianScheme.H"
+#include "uncorrectedSnGrad.H"
+#include "surfaceInterpolate.H"
+#include "fvcSurfaceIntegrate.H"
+#include "slicedSurfaceFields.H"
+#include "syncTools.H"
+
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::explicitLTSSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    Info<< "MULES: Solving for " << psi.name() << endl;
+
+    const fvMesh& mesh = psi.mesh();
+    psi.correctBoundaryConditions();
+
+    surfaceScalarField phiBD = upwind<scalar>(psi.mesh(), phi).flux(psi);
+
+    surfaceScalarField& phiCorr = phiPsi;
+    phiCorr -= phiBD;
+
+    scalarField allLambda(mesh.nFaces(), 1.0);
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    limiter
+    (
+        allLambda,
+        rho,
+        psi,
+        phiBD,
+        phiCorr,
+        Sp.field(),
+        Su.field(),
+        psiMax,
+        psiMin,
+        3
+    );
+
+    phiPsi = phiBD + lambda*phiCorr;
+
+    scalarField& psiIf = psi;
+    const scalarField& psi0 = psi.oldTime();
+
+    const volScalarField& rDeltaT =
+        mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
+
+    psiIf = 0.0;
+    fvc::surfaceIntegrate(psiIf, phiPsi);
+
+    if (mesh.moving())
+    {
+        psiIf =
+        (
+            mesh.Vsc0()*rho.oldTime()*psi0*rDeltaT/mesh.Vsc()
+          + Su.field()
+          - psiIf
+        )/(rho*rDeltaT - Sp.field());
+    }
+    else
+    {
+        psiIf =
+        (
+            rho.oldTime()*psi0*rDeltaT
+          + Su.field()
+          - psiIf
+        )/(rho*rDeltaT - Sp.field());
+    }
+
+    psi.correctBoundaryConditions();
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::implicitSolve
+(
+    const RhoType& rho,
+    volScalarField& psi,
+    const surfaceScalarField& phi,
+    surfaceScalarField& phiPsi,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin
+)
+{
+    const fvMesh& mesh = psi.mesh();
+
+    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
+
+    label maxIter
+    (
+        readLabel(MULEScontrols.lookup("maxIter"))
+    );
+
+    label nLimiterIter
+    (
+        readLabel(MULEScontrols.lookup("nLimiterIter"))
+    );
+
+    scalar maxUnboundedness
+    (
+        readScalar(MULEScontrols.lookup("maxUnboundedness"))
+    );
+
+    scalar CoCoeff
+    (
+        readScalar(MULEScontrols.lookup("CoCoeff"))
+    );
+
+    scalarField allCoLambda(mesh.nFaces());
+
+    {
+        surfaceScalarField Cof =
+            mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs()
+           *mag(phi)/mesh.magSf();
+
+        slicedSurfaceScalarField CoLambda
+        (
+            IOobject
+            (
+                "CoLambda",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimless,
+            allCoLambda,
+            false   // Use slices for the couples
+        );
+
+        CoLambda == 1.0/max(CoCoeff*Cof, scalar(1));
+    }
+
+    scalarField allLambda(allCoLambda);
+    //scalarField allLambda(mesh.nFaces(), 1.0);
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    linear<scalar> CDs(mesh);
+    upwind<scalar> UDs(mesh, phi);
+    //fv::uncorrectedSnGrad<scalar> snGrads(mesh);
+
+    fvScalarMatrix psiConvectionDiffusion
+    (
+        fvm::ddt(rho, psi)
+      + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi)
+        //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads)
+        //.fvmLaplacian(Dpsif, psi)
+      - fvm::Sp(Sp, psi)
+      - Su
+    );
+
+    surfaceScalarField phiBD = psiConvectionDiffusion.flux();
+
+    surfaceScalarField& phiCorr = phiPsi;
+    phiCorr -= phiBD;
+
+    for (label i=0; i<maxIter; i++)
+    {
+        if (i != 0 && i < 4)
+        {
+            allLambda = allCoLambda;
+        }
+
+        limiter
+        (
+            allLambda,
+            rho,
+            psi,
+            phiBD,
+            phiCorr,
+            Sp.field(),
+            Su.field(),
+            psiMax,
+            psiMin,
+            nLimiterIter
+        );
+
+        solve
+        (
+            psiConvectionDiffusion + fvc::div(lambda*phiCorr),
+            MULEScontrols
+        );
+
+        scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;
+        scalar minPsi = gMin(psi.internalField());
+
+        scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0));
+
+        if (unboundedness < maxUnboundedness)
+        {
+            break;
+        }
+        else
+        {
+            Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1
+                << " min(" << psi.name() << ") = " << minPsi << endl;
+
+            phiBD = psiConvectionDiffusion.flux();
+
+            /*
+            word gammaScheme("div(phi,gamma)");
+            word gammarScheme("div(phirb,gamma)");
+
+            const surfaceScalarField& phir =
+                mesh.lookupObject<surfaceScalarField>("phir");
+
+            phiCorr =
+                fvc::flux
+                (
+                    phi,
+                    psi,
+                    gammaScheme
+                )
+              + fvc::flux
+                (
+                    -fvc::flux(-phir, scalar(1) - psi, gammarScheme),
+                    psi,
+                    gammarScheme
+                )
+                - phiBD;
+            */
+        }
+    }
+
+    phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr;
+}
+
+
+template<class RhoType, class SpType, class SuType>
+void Foam::MULES::limiter
+(
+    scalarField& allLambda,
+    const RhoType& rho,
+    const volScalarField& psi,
+    const surfaceScalarField& phiBD,
+    const surfaceScalarField& phiCorr,
+    const SpType& Sp,
+    const SuType& Su,
+    const scalar psiMax,
+    const scalar psiMin,
+    const label nLimiterIter
+)
+{
+    const scalarField& psiIf = psi;
+    const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();
+
+    const scalarField& psi0 = psi.oldTime();
+
+    const fvMesh& mesh = psi.mesh();
+
+    const unallocLabelList& owner = mesh.owner();
+    const unallocLabelList& neighb = mesh.neighbour();
+    tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
+    const scalarField& V = tVsc();
+
+    const volScalarField& rDeltaT =
+        mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT");
+
+    const scalarField& phiBDIf = phiBD;
+    const surfaceScalarField::GeometricBoundaryField& phiBDBf =
+        phiBD.boundaryField();
+
+    const scalarField& phiCorrIf = phiCorr;
+    const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
+        phiCorr.boundaryField();
+
+    slicedSurfaceScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh,
+        dimless,
+        allLambda,
+        false   // Use slices for the couples
+    );
+
+    scalarField& lambdaIf = lambda;
+    surfaceScalarField::GeometricBoundaryField& lambdaBf =
+        lambda.boundaryField();
+
+    scalarField psiMaxn(psiIf.size(), psiMin);
+    scalarField psiMinn(psiIf.size(), psiMax);
+
+    scalarField sumPhiBD(psiIf.size(), 0.0);
+
+    scalarField sumPhip(psiIf.size(), VSMALL);
+    scalarField mSumPhim(psiIf.size(), VSMALL);
+
+    forAll(phiCorrIf, facei)
+    {
+        label own = owner[facei];
+        label nei = neighb[facei];
+
+        psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
+        psiMinn[own] = min(psiMinn[own], psiIf[nei]);
+
+        psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
+        psiMinn[nei] = min(psiMinn[nei], psiIf[own]);
+
+        sumPhiBD[own] += phiBDIf[facei];
+        sumPhiBD[nei] -= phiBDIf[facei];
+
+        scalar phiCorrf = phiCorrIf[facei];
+
+        if (phiCorrf > 0.0)
+        {
+            sumPhip[own] += phiCorrf;
+            mSumPhim[nei] += phiCorrf;
+        }
+        else
+        {
+            mSumPhim[own] -= phiCorrf;
+            sumPhip[nei] -= phiCorrf;
+        }
+    }
+
+    forAll(phiCorrBf, patchi)
+    {
+        const fvPatchScalarField& psiPf = psiBf[patchi];
+        const scalarField& phiBDPf = phiBDBf[patchi];
+        const scalarField& phiCorrPf = phiCorrBf[patchi];
+
+        const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+        if (psiPf.coupled())
+        {
+            scalarField psiPNf = psiPf.patchNeighbourField();
+
+            forAll(phiCorrPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]);
+            }
+        }
+        else
+        {
+            forAll(phiCorrPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]);
+                psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]);
+            }
+        }
+
+        forAll(phiCorrPf, pFacei)
+        {
+            label pfCelli = pFaceCells[pFacei];
+
+            sumPhiBD[pfCelli] += phiBDPf[pFacei];
+
+            scalar phiCorrf = phiCorrPf[pFacei];
+
+            if (phiCorrf > 0.0)
+            {
+                sumPhip[pfCelli] += phiCorrf;
+            }
+            else
+            {
+                mSumPhim[pfCelli] -= phiCorrf;
+            }
+        }
+    }
+
+    psiMaxn = min(psiMaxn, psiMax);
+    psiMinn = max(psiMinn, psiMin);
+
+    //scalar smooth = 0.5;
+    //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax);
+    //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin);
+
+    if (mesh.moving())
+    {
+        tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0();
+
+        psiMaxn =
+            V*((rho*rDeltaT - Sp)*psiMaxn - Su)
+          - (V0()*rDeltaT)*rho.oldTime()*psi0
+          + sumPhiBD;
+
+        psiMinn =
+            V*(Su - (rho*rDeltaT - Sp)*psiMinn)
+          + (V0*rDeltaT)*rho.oldTime()*psi0
+          - sumPhiBD;
+    }
+    else
+    {
+        psiMaxn =
+            V*((rho*rDeltaT - Sp)*psiMaxn - (rho.oldTime()*rDeltaT)*psi0 - Su)
+          + sumPhiBD;
+
+        psiMinn =
+            V*((rho*rDeltaT)*psi0 - (rho.oldTime()*rDeltaT - Sp)*psiMinn + Su)
+          - sumPhiBD;
+    }
+
+    scalarField sumlPhip(psiIf.size());
+    scalarField mSumlPhim(psiIf.size());
+
+    for(int j=0; j<nLimiterIter; j++)
+    {
+        sumlPhip = 0.0;
+        mSumlPhim = 0.0;
+
+        forAll(lambdaIf, facei)
+        {
+            label own = owner[facei];
+            label nei = neighb[facei];
+
+            scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei];
+
+            if (lambdaPhiCorrf > 0.0)
+            {
+                sumlPhip[own] += lambdaPhiCorrf;
+                mSumlPhim[nei] += lambdaPhiCorrf;
+            }
+            else
+            {
+                mSumlPhim[own] -= lambdaPhiCorrf;
+                sumlPhip[nei] -= lambdaPhiCorrf;
+            }
+        }
+
+        forAll(lambdaBf, patchi)
+        {
+            scalarField& lambdaPf = lambdaBf[patchi];
+            const scalarField& phiCorrfPf = phiCorrBf[patchi];
+
+            const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+            forAll(lambdaPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei];
+
+                if (lambdaPhiCorrf > 0.0)
+                {
+                    sumlPhip[pfCelli] += lambdaPhiCorrf;
+                }
+                else
+                {
+                    mSumlPhim[pfCelli] -= lambdaPhiCorrf;
+                }
+            }
+        }
+
+        forAll (sumlPhip, celli)
+        {
+            sumlPhip[celli] =
+                max(min
+                (
+                    (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli],
+                    1.0), 0.0
+                );
+
+            mSumlPhim[celli] =
+                max(min
+                (
+                    (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli],
+                    1.0), 0.0
+                );
+        }
+
+        const scalarField& lambdam = sumlPhip;
+        const scalarField& lambdap = mSumlPhim;
+
+        forAll(lambdaIf, facei)
+        {
+            if (phiCorrIf[facei] > 0.0)
+            {
+                lambdaIf[facei] = min
+                (
+                    lambdaIf[facei],
+                    min(lambdap[owner[facei]], lambdam[neighb[facei]])
+                );
+            }
+            else
+            {
+                lambdaIf[facei] = min
+                (
+                    lambdaIf[facei],
+                    min(lambdam[owner[facei]], lambdap[neighb[facei]])
+                );
+            }
+        }
+
+
+        forAll(lambdaBf, patchi)
+        {
+            fvsPatchScalarField& lambdaPf = lambdaBf[patchi];
+            const scalarField& phiCorrfPf = phiCorrBf[patchi];
+
+            const labelList& pFaceCells = mesh.boundary()[patchi].faceCells();
+
+            forAll(lambdaPf, pFacei)
+            {
+                label pfCelli = pFaceCells[pFacei];
+
+                if (phiCorrfPf[pFacei] > 0.0)
+                {
+                    lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]);
+                }
+                else
+                {
+                    lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]);
+                }
+            }
+        }
+
+        syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>());
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
new file mode 100644
index 0000000000..205812e978
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files
@@ -0,0 +1,4 @@
+LTSInterFoam.C
+MULES.C
+
+EXE = $(FOAM_APPBIN)/LTSInterFoam
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options
new file mode 100644
index 0000000000..24349f694e
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options
@@ -0,0 +1,15 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -ltwoPhaseInterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleTurbulenceModel \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
new file mode 100644
index 0000000000..a4cb6e7688
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H
@@ -0,0 +1,36 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phic = mag(phi/mesh.magSf());
+    phic = min(interface.cAlpha()*phic, max(phic));
+    surfaceScalarField phir = phic*interface.nHatf();
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        surfaceScalarField phiAlpha =
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                alpha1,
+                alpharScheme
+            );
+
+        MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0);
+        //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
+
+        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
+    }
+
+    Info<< "Liquid phase volume fraction = "
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
+        << endl;
+}
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
new file mode 100644
index 0000000000..9aae37a8be
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H
@@ -0,0 +1,35 @@
+label nAlphaCorr
+(
+    readLabel(piso.lookup("nAlphaCorr"))
+);
+
+label nAlphaSubCycles
+(
+    readLabel(piso.lookup("nAlphaSubCycles"))
+);
+
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+#       include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+#       include "alphaEqn.H"
+}
+
+interface.correct();
+
+rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
new file mode 100644
index 0000000000..1f2f082fb3
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H
@@ -0,0 +1,30 @@
+scalar maxDeltaT
+(
+    piso.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+);
+
+volScalarField rDeltaT
+(
+    IOobject
+    (
+        "rDeltaT",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
+);
+
+volScalarField rSubDeltaT
+(
+    IOobject
+    (
+        "rSubDeltaT",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT)
+);
diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
new file mode 100644
index 0000000000..f5b8450ddd
--- /dev/null
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
@@ -0,0 +1,109 @@
+{
+    scalar maxCo
+    (
+        piso.lookupOrDefault<scalar>("maxCo", 0.9)
+    );
+
+    scalar maxAlphaCo
+    (
+        piso.lookupOrDefault<scalar>("maxAlphaCo", 0.2)
+    );
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        piso.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    label nAlphaSpreadIter
+    (
+        piso.lookupOrDefault<label>("nAlphaSpreadIter", 1)
+    );
+
+    label nAlphaSweepIter
+    (
+        piso.lookupOrDefault<label>("nAlphaSweepIter", 5)
+    );
+
+    scalar rDeltaTDampingCoeff
+    (
+        piso.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+    );
+
+    scalar maxDeltaT
+    (
+        piso.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    );
+
+    volScalarField rDeltaT0 = rDeltaT;
+
+    // Set the reciprocal time-step using an effective maximum Courant number
+    rDeltaT = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum
+        (
+            mag(rhoPhi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf())
+        )/rho
+    );
+
+    // Limit the time-step further in the region of the interface
+    {
+        surfaceScalarField alphaf = fvc::interpolate(alpha1);
+
+        surfaceScalarField SfUfbyDelta =
+            pos(alphaf - 0.01)*pos(0.99 - alphaf)
+           *mesh.surfaceInterpolation::deltaCoeffs()*mag(phi);
+
+        rDeltaT = max
+        (
+            rDeltaT,
+            fvc::surfaceSum(mag(SfUfbyDelta/(maxAlphaCo*mesh.magSf())))
+        );
+    }
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    if (nAlphaSpreadIter > 0)
+    {
+        fvc::spread(rDeltaT, alpha1, nAlphaSpreadIter);
+    }
+
+    if (nAlphaSweepIter > 0)
+    {
+        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter);
+    }
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        Info<< "Damping rDeltaT" << endl;
+        rDeltaT = rDeltaT0*max(rDeltaT/rDeltaT0, 1.0 - rDeltaTDampingCoeff);
+    }
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    label nAlphaSubCycles
+    (
+        readLabel(piso.lookup("nAlphaSubCycles"))
+    );
+
+    rSubDeltaT = rDeltaT*nAlphaSubCycles;
+}