diff --git a/ThirdParty/.gitignore b/ThirdParty/.gitignore
new file mode 100644
index 000000000..947ff5a8a
--- /dev/null
+++ b/ThirdParty/.gitignore
@@ -0,0 +1,19 @@
+# git-ls-files --others --exclude-from=.git/info/exclude
+# Lines that start with '#' are comments.
+
+# editor and misc backup files - anywhere
+*~
+.*~
+*.bak
+*.bak[0-9][0-9]
+*.orig
+*.orig[0-9][0-9]
+\#*\#
+
+# ignore the build
+rpmBuild/BUILD
+rpmBuild/RPMS
+rpmBuild/rpmDB
+rpmBuild/tmp
+
+# end-of-file
diff --git a/ThirdParty/rpmBuild/SPECS/ParaView-3.12.0.spec b/ThirdParty/rpmBuild/SPECS/ParaView-3.12.0.spec
index b562a58c4..504d7cb7d 100644
--- a/ThirdParty/rpmBuild/SPECS/ParaView-3.12.0.spec
+++ b/ThirdParty/rpmBuild/SPECS/ParaView-3.12.0.spec
@@ -77,8 +77,8 @@ URL:            http://www.paraview.org/files/v3.12/
 Source: 		%url/%{name}-%{version}.tar.gz
 Prefix: 		%{_prefix}
 Group: 			Development/Tools
-Patch0:                 ParaView-3.12.0.patch_darwin
-Patch1:                 paraview-gcc47.patch
+Patch0:         ParaView-3.12.0.patch_darwin
+Patch1:         paraview-gcc47.patch
 
 %define _installPrefix  %{_prefix}/packages/%{name}-%{version}/platforms/%{_WM_OPTIONS}
 
diff --git a/etc/bashrc b/etc/bashrc
index b4250770d..689e347b0 100644
--- a/etc/bashrc
+++ b/etc/bashrc
@@ -251,22 +251,22 @@ Darwin)
     then
         # Use Mac-Ports-Compiler instead of Apple-gcc-4.2
 	case $WM_COMPILER in
-	    Gcc43)
-		export WM_CC='gcc-mp-4.3'
-		export WM_CXX='g++-mp-4.3'
-		;;
-	    Gcc44)
-		export WM_CC='gcc-mp-4.4'
-		export WM_CXX='g++-mp-4.4'
-		;;
-	    Gcc45)
-		export WM_CC='gcc-mp-4.5'
-		export WM_CXX='g++-mp-4.5'
-		;;
-	    *)
-		export WM_COMPILER=
-		;;
-	esac
+        Gcc43)
+        export WM_CC='gcc-mp-4.3'
+        export WM_CXX='g++-mp-4.3'
+        ;;
+        Gcc44)
+        export WM_CC='gcc-mp-4.4'
+        export WM_CXX='g++-mp-4.4'
+        ;;
+        Gcc45)
+        export WM_CC='gcc-mp-4.5'
+        export WM_CXX='g++-mp-4.5'
+        ;;
+        *)
+        export WM_COMPILER=
+        ;;
+    esac
     fi
     ;;
 
diff --git a/etc/cshrc b/etc/cshrc
index e39a7139c..973c34fac 100644
--- a/etc/cshrc
+++ b/etc/cshrc
@@ -79,9 +79,9 @@ alias _foamSource 'if ($?FOAM_VERBOSE && $?prompt) echo "Sourcing: \!*"; source
 # Add in preset user or site preferences:
 set foamPrefs=`$WM_PROJECT_DIR/bin/foamEtcFile prefs.csh`
 if ( $status == 0 ) then
-    if ($foamPrefs) then
+    if ( $?foamPrefs ) then
         _foamSource $foamPrefs
-     endif
+    endif
 endif
 unset foamPrefs
 
@@ -110,7 +110,7 @@ if ( ! $?WM_OSTYPE ) setenv WM_OSTYPE POSIX
 # Compiler: set to Gcc, Gcc43 or Icc (for Intel's icc)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 if ( ! $?WM_COMPILER ) setenv WM_COMPILER Gcc
-setenv WM_COMPILER Gcc
+#setenv WM_COMPILER Gcc
 #setenv WM_COMPILER Icc
 
 setenv WM_COMPILER_ARCH
@@ -129,7 +129,7 @@ if ( ! $?WM_PRECISION_OPTION ) setenv WM_PRECISION_OPTION DP
 
 # WM_COMPILE_OPTION = Opt | Debug | Prof
 if ( ! $?WM_COMPILE_OPTION ) setenv WM_COMPILE_OPTION Opt
-setenv WM_COMPILE_OPTION Opt
+#setenv WM_COMPILE_OPTION Opt
 #setenv WM_COMPILE_OPTION Debug
 
 # WM_MPLIB = | OPENMPI | MPICH | MPICH-GM | HPMPI | GAMMA | MPI | QSMPI
@@ -247,20 +247,20 @@ case Darwin:
     # Use Mac-Ports-Compiler instead of Apple-gcc-4.2
     if ( ! $?compilerInstall ) setenv compilerInstall System
     if ( $compilerInstall == "System" ) then 
-	switch ($WM_COMPILER)
-	case Gcc43:
-	    setenv WM_CC 'gcc-mp-4.3'
-	    setenv WM_CXX 'g++-mp-4.3'
-	    breaksw;
-	case Gcc44:
-	    setenv WM_CC 'gcc-mp-4.4'
-	    setenv WM_CXX 'g++-mp-4.4'
-	    breaksw;
-	case Gcc45:
-	    setenv WM_CC 'gcc-mp-4.5'
-	    setenv WM_CXX 'g++-mp-4.5'
-	    breaksw;
-	endsw
+    switch ($WM_COMPILER)
+    case Gcc43:
+        setenv WM_CC 'gcc-mp-4.3'
+        setenv WM_CXX 'g++-mp-4.3'
+        breaksw;
+    case Gcc44:
+        setenv WM_CC 'gcc-mp-4.4'
+        setenv WM_CXX 'g++-mp-4.4'
+        breaksw;
+    case Gcc45:
+        setenv WM_CC 'gcc-mp-4.5'
+        setenv WM_CXX 'g++-mp-4.5'
+        breaksw;
+    endsw
         #setenv WM_COMPILER ""
     endif
 
diff --git a/etc/prefs.csh-EXAMPLE b/etc/prefs.csh-EXAMPLE
index 816907aeb..9599a5195 100644
--- a/etc/prefs.csh-EXAMPLE
+++ b/etc/prefs.csh-EXAMPLE
@@ -140,26 +140,13 @@
 #setenv PARAVIEW_DIR path_to_system_installed_paraview
 #setenv PARAVIEW_BIN_DIR $PARAVIEW_DIR/bin
 
-# System installed bison
-#setenv BISON_SYSTEM 1
-
-# System installed flex
-#setenv FLEX_SYSTEM 1
-
-# System installed m4
-#setenv M4_SYSTEM 1
-
-# System installed PyFoam
-#setenv PYFOAM_SYSTEM 1
-
-# System installed swak4Foam
-#setenv SWAK4FOAM_SYSTEM 1
 
 # Specify ParaView version
 # ~~~~~~~~~~~~~~~~~~~~~~~~
 #setenv ParaView_VERSION git        # eg, cvs/git version
 #setenv ParaView_MAJOR 3.7
 
+
 # System identifier for the OpenFOAM CDash test harness on openfoam-extend
 #
 # By default, your system FQN/hostname will be used as the system identifier
diff --git a/etc/prefs.sh-EXAMPLE b/etc/prefs.sh-EXAMPLE
index 44800c433..990207083 100644
--- a/etc/prefs.sh-EXAMPLE
+++ b/etc/prefs.sh-EXAMPLE
@@ -151,17 +151,12 @@ compilerInstall=System
 # System installed m4
 #export M4_SYSTEM=1
 
-# System installed PyFoam
-#export PYFOAM_SYSTEM=1
-
-# System installed swak4Foam
-#export SWAK4FOAM_SYSTEM=1
-
 # Specify ParaView version
 # ~~~~~~~~~~~~~~~~~~~~~~~~
 #export ParaView_VERSION=git        # eg, cvs/git version
 #export ParaView_MAJOR=3.7
 
+
 # System identifier for the OpenFOAM CDash test harness on openfoam-extend
 #
 # By default, your system FQN/hostname will be used as the system identifier
diff --git a/etc/settings.csh b/etc/settings.csh
index 879ba8c2c..3e97c86be 100644
--- a/etc/settings.csh
+++ b/etc/settings.csh
@@ -450,21 +450,23 @@ endif
 
 # QT
 # ~~~~~
-if ( $?QT_SYSTEM == 0 && -e "$WM_THIRD_PARTY_DIR"/packages/qt-everywhere-opensource-src-4.7.0/platforms/$WM_OPTIONS )then
-    _foamSource $WM_THIRD_PARTY_DIR/packages/qt-everywhere-opensource-src-4.7.0/platforms/$WM_OPTIONS/etc/qt-everywhere-opensource-src-4.7.0.csh
+if ( $?QT_SYSTEM == 0 && -e "$WM_THIRD_PARTY_DIR"/packages/qt-everywhere-opensource-src-4.7.4/platforms/$WM_OPTIONS )then
+    _foamSource $WM_THIRD_PARTY_DIR/packages/qt-everywhere-opensource-src-4.7.4/platforms/$WM_OPTIONS/etc/qt-everywhere-opensource-src-4.7.4.csh
+#if ( $?QT_SYSTEM == 0 && -e "$WM_THIRD_PARTY_DIR"/packages/qt-everywhere-opensource-src-4.7.0/platforms/$WM_OPTIONS )then
+#    _foamSource $WM_THIRD_PARTY_DIR/packages/qt-everywhere-opensource-src-4.7.0/platforms/$WM_OPTIONS/etc/qt-everywhere-opensource-src-4.7.0.csh
 endif
 
 # PARAVIEW
 # ~~~~~
 if ( $?PARAVIEW_SYSTEM == 0 && -e "$WM_THIRD_PARTY_DIR"/packages/ParaView-3.12.0/platforms/$WM_OPTIONS ) then
     _foamSource $WM_THIRD_PARTY_DIR/packages/ParaView-3.12.0/platforms/$WM_OPTIONS/etc/ParaView-3.12.0.csh
+
 #if ( $?PARAVIEW_SYSTEM == 0 && -e "$WM_THIRD_PARTY_DIR"/packages/ParaView-3.8.1/platforms/$WM_OPTIONS ) then
 #    _foamSource $WM_THIRD_PARTY_DIR/packages/ParaView-3.8.1/platforms/$WM_OPTIONS/etc/ParaView-3.8.1.csh
-endif
 
-#if ( $?PARAVIEW_SYSTEM == 0 && -e "$WM_THIRD_PARTY_DIR"/packages/ParaView-3.14.1 ) then
+#if ( $?PARAVIEW_SYSTEM == 0 && -e "$WM_THIRD_PARTY_DIR"/packages/ParaView-3.14.1/platforms/$WM_OPTIONS ) then
 #    _foamSource $WM_THIRD_PARTY_DIR/packages/ParaView-3.14.1/platforms/$WM_OPTIONS/etc/ParaView-3.14.1.csh
-#endif
+endif
 
 if ( $WM_ARCH == "darwinIntel" || $WM_ARCH == "darwinIntel64" ) then
     setenv DYLD_LIBRARY_PATH ${LD_LIBRARY_PATH}
diff --git a/etc/settings.sh b/etc/settings.sh
index 5b12f2e5c..c62b7b56a 100644
--- a/etc/settings.sh
+++ b/etc/settings.sh
@@ -487,13 +487,10 @@ export MPI_BUFFER_SIZE
 
 # Load cmake
 # ~~~~~~~~~~
-<<<<<<< Temporary merge branch 1
 [ -z "$CMAKE_SYSTEM" ] && [ -e $WM_THIRD_PARTY_DIR/packages/cmake-2.8.8/platforms/$WM_OPTIONS ] && {
     _foamSource $WM_THIRD_PARTY_DIR/packages/cmake-2.8.8/platforms/$WM_OPTIONS/etc/cmake-2.8.8.sh
-=======
-[ -z "$CMAKE_SYSTEM" ] && [ -e $WM_THIRD_PARTY_DIR/packages/cmake-2.8.6/platforms/$WM_OPTIONS ] && {
-    _foamSource $WM_THIRD_PARTY_DIR/packages/cmake-2.8.6/platforms/$WM_OPTIONS/etc/cmake-2.8.6.sh
->>>>>>> Temporary merge branch 2
+#[ -z "$CMAKE_SYSTEM" ] && [ -e $WM_THIRD_PARTY_DIR/packages/cmake-2.8.6/platforms/$WM_OPTIONS ] && {
+#    _foamSource $WM_THIRD_PARTY_DIR/packages/cmake-2.8.6/platforms/$WM_OPTIONS/etc/cmake-2.8.6.sh
 }
 [ "$FOAM_VERBOSE" -a "$PS1" ] && echo "    CMAKE_DIR is initialized to: $CMAKE_DIR"
 
@@ -521,8 +518,8 @@ export MPI_BUFFER_SIZE
 
 # Load zoltan
 # ~~~~~~~~~~
-[ -z "$ZOLTAN_SYSTEM" ] && [ -e $WM_THIRD_PARTY_DIR/packages/zoltan-3.5/platforms/$WM_OPTIONS ] && {
-    _foamSource $WM_THIRD_PARTY_DIR/packages/zoltan-3.5/platforms/$WM_OPTIONS/etc/zoltan-3.5.sh
+[ -z "$ZOLTAN_SYSTEM" ] && [ -e $WM_THIRD_PARTY_DIR/packages/zoltan_3.5 ] && {
+    _foamSource $WM_THIRD_PARTY_DIR/packages/zoltan_3.5/platforms/$WM_OPTIONS/etc/zoltan_3.5.sh
 }
 [ "$FOAM_VERBOSE" -a "$PS1" ] && echo "    ZOLTAN_DIR is initialized to: $ZOLTAN_DIR"
 
@@ -554,6 +551,8 @@ export MPI_BUFFER_SIZE
 # ~~~~~~~~~~~~~
 [ -z "$PARAVIEW_SYSTEM" ] && [ -e $WM_THIRD_PARTY_DIR/packages/ParaView-3.12.0/platforms/$WM_OPTIONS ] && {
     _foamSource $WM_THIRD_PARTY_DIR/packages/ParaView-3.12.0/platforms/$WM_OPTIONS/etc/ParaView-3.12.0.sh
+#[ -z "$PARAVIEW_SYSTEM" ] && [ -e $WM_THIRD_PARTY_DIR/packages/ParaView-3.8.1/platforms/$WM_OPTIONS ] && {
+#    _foamSource $WM_THIRD_PARTY_DIR/packages/ParaView-3.8.1/platforms/$WM_OPTIONS/etc/ParaView-3.8.1.sh
 }
 [ "$FOAM_VERBOSE" -a "$PS1" ] && echo "    PARAVIEW_DIR is initialized to: $PARAVIEW_DIR"
 
diff --git a/src/OpenFOAM/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C b/src/OpenFOAM/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C
index c6511e42e..45481d726 100644
--- a/src/OpenFOAM/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C
+++ b/src/OpenFOAM/fields/PointPatchFields/constraint/processor/ProcessorPointPatchField.C
@@ -48,7 +48,8 @@ template
     class Type
 >
 void ProcessorPointPatchField
-<PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::resizeBuf
+<PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::
+resizeBuf
 (
     List<char>& buf, 
     const label size
@@ -961,7 +962,8 @@ initInterfaceMatrixUpdate
     const lduMatrix& m,
     const scalarField& coeffs,
     const direction,
-    const Pstream::commsTypes commsType
+    const Pstream::commsTypes commsType,
+    const bool switchToLhs
 ) const
 {
     tmp<scalarField> tlocalMult(new scalarField(this->size(), 0));
@@ -985,95 +987,191 @@ initInterfaceMatrixUpdate
     // use the counter.  
     label coeffI = 0;
 
-    // Owner side
-    // ~~~~~~~~~~
+    if (switchToLhs)
     {
-        const labelList& cutOwn = procPatch_.cutEdgeOwnerIndices();
-        const labelList& cutOwnStart = procPatch_.cutEdgeOwnerStart();
-
-        forAll (mp, pointI)
+        // Owner side
+        // ~~~~~~~~~~
         {
-            label ownIndex = cutOwnStart[pointI];
-            label endOwn = cutOwnStart[pointI + 1];
+            const labelList& cutOwn = procPatch_.cutEdgeOwnerIndices();
+            const labelList& cutOwnStart = procPatch_.cutEdgeOwnerStart();
 
-            for (; ownIndex < endOwn; ownIndex++)
+            forAll (mp, pointI)
             {
-                localMult[pointI] +=
-                    coeffs[coeffI]*psiInternal[U[cutOwn[ownIndex]]];
+                label ownIndex = cutOwnStart[pointI];
+                label endOwn = cutOwnStart[pointI + 1];
 
-                // Multiply the internal side as well using the cut mask
-                result[U[cutOwn[ownIndex]]] +=
-                    cutMask[coeffI]*coeffs[coeffI]*psiInternal[mp[pointI]];
+                for (; ownIndex < endOwn; ownIndex++)
+                {
+                    localMult[pointI] +=
+                        coeffs[coeffI]*psiInternal[U[cutOwn[ownIndex]]];
 
+                    // Multiply the internal side as well using the cut mask
+                    result[U[cutOwn[ownIndex]]] -=
+                        cutMask[coeffI]*coeffs[coeffI]*psiInternal[mp[pointI]];
+
+                    coeffI++;
+                }
+            }
+        }
+
+        // Neighbour side
+        // ~~~~~~~~~~~~~~
+        {
+            const labelList& cutNei = procPatch_.cutEdgeNeighbourIndices();
+            const labelList& cutNeiStart = procPatch_.cutEdgeNeighbourStart();
+
+            forAll (mp, pointI)
+            {
+                label neiIndex = cutNeiStart[pointI];
+                label endNei = cutNeiStart[pointI + 1];
+
+                for (; neiIndex < endNei; neiIndex++)
+                {
+                    localMult[pointI] +=
+                        coeffs[coeffI]*psiInternal[L[cutNei[neiIndex]]];
+
+                    // Multiply the internal side as well using the cut mask
+                    result[L[cutNei[neiIndex]]] -=
+                        cutMask[coeffI]*coeffs[coeffI]*psiInternal[mp[pointI]];
+
+                    coeffI++;
+                }
+            }
+        }
+
+        // Doubly cut coefficients
+        // ~~~~~~~~~~~~~~~~~~~~~~~
+
+        // There exists a possibility of having an internal edge for a
+        // point on the processor patch which is in fact connected to
+        // another point of the same patch.  This particular nastiness
+        // introduces a deformation in the solution because the edge is
+        // either multiplied twice or not at all.  For this purpose, the
+        // offending edges need to be separated out and multiplied
+        // appropriately.  This will only happen for cell tetrahedral
+        // decomposition and is generally nasty.  
+        // No need for cut mask here
+        {
+            const labelList& doubleCut = procPatch_.doubleCutEdgeIndices();
+
+            const labelList& doubleCutOwner = procPatch_.doubleCutOwner();
+            const labelList& doubleCutNeighbour =
+                procPatch_.doubleCutNeighbour();
+
+            forAll (doubleCut, edgeI)
+            {
+                // Owner side
+                localMult[doubleCutOwner[edgeI]] +=
+                    coeffs[coeffI]*psiInternal[U[doubleCut[edgeI]]];
+                coeffI++;
+
+                // Neighbour side
+                localMult[doubleCutNeighbour[edgeI]] +=
+                    coeffs[coeffI]*psiInternal[L[doubleCut[edgeI]]];
                 coeffI++;
             }
         }
-    }
 
-    // Neighbour side
-    // ~~~~~~~~~~~~~~
-    {
-        const labelList& cutNei = procPatch_.cutEdgeNeighbourIndices();
-        const labelList& cutNeiStart = procPatch_.cutEdgeNeighbourStart();
+        // Add the local multiplication to this side as well
 
         forAll (mp, pointI)
         {
-            label neiIndex = cutNeiStart[pointI];
-            label endNei = cutNeiStart[pointI + 1];
+            result[mp[pointI]] -= localMult[pointI];
+        }
+    }
+    else
+    {
+        // Owner side
+        // ~~~~~~~~~~
+        {
+            const labelList& cutOwn = procPatch_.cutEdgeOwnerIndices();
+            const labelList& cutOwnStart = procPatch_.cutEdgeOwnerStart();
 
-            for (; neiIndex < endNei; neiIndex++)
+            forAll (mp, pointI)
             {
-                localMult[pointI] +=
-                    coeffs[coeffI]*psiInternal[L[cutNei[neiIndex]]];
+                label ownIndex = cutOwnStart[pointI];
+                label endOwn = cutOwnStart[pointI + 1];
 
-                // Multiply the internal side as well using the cut mask
-                result[L[cutNei[neiIndex]]] +=
-                    cutMask[coeffI]*coeffs[coeffI]*psiInternal[mp[pointI]];
+                for (; ownIndex < endOwn; ownIndex++)
+                {
+                    localMult[pointI] +=
+                        coeffs[coeffI]*psiInternal[U[cutOwn[ownIndex]]];
 
+                    // Multiply the internal side as well using the cut mask
+                    result[U[cutOwn[ownIndex]]] +=
+                        cutMask[coeffI]*coeffs[coeffI]*psiInternal[mp[pointI]];
+
+                    coeffI++;
+                }
+            }
+        }
+
+        // Neighbour side
+        // ~~~~~~~~~~~~~~
+        {
+            const labelList& cutNei = procPatch_.cutEdgeNeighbourIndices();
+            const labelList& cutNeiStart = procPatch_.cutEdgeNeighbourStart();
+
+            forAll (mp, pointI)
+            {
+                label neiIndex = cutNeiStart[pointI];
+                label endNei = cutNeiStart[pointI + 1];
+
+                for (; neiIndex < endNei; neiIndex++)
+                {
+                    localMult[pointI] +=
+                        coeffs[coeffI]*psiInternal[L[cutNei[neiIndex]]];
+
+                    // Multiply the internal side as well using the cut mask
+                    result[L[cutNei[neiIndex]]] +=
+                        cutMask[coeffI]*coeffs[coeffI]*psiInternal[mp[pointI]];
+
+                    coeffI++;
+                }
+            }
+        }
+
+        // Doubly cut coefficients
+        // ~~~~~~~~~~~~~~~~~~~~~~~
+
+        // There exists a possibility of having an internal edge for a
+        // point on the processor patch which is in fact connected to
+        // another point of the same patch.  This particular nastiness
+        // introduces a deformation in the solution because the edge is
+        // either multiplied twice or not at all.  For this purpose, the
+        // offending edges need to be separated out and multiplied
+        // appropriately.  This will only happen for cell tetrahedral
+        // decomposition and is generally nasty.  
+        // No need for cut mask here
+        {
+            const labelList& doubleCut = procPatch_.doubleCutEdgeIndices();
+
+            const labelList& doubleCutOwner = procPatch_.doubleCutOwner();
+            const labelList& doubleCutNeighbour =
+                procPatch_.doubleCutNeighbour();
+
+            forAll (doubleCut, edgeI)
+            {
+                // Owner side
+                localMult[doubleCutOwner[edgeI]] +=
+                    coeffs[coeffI]*psiInternal[U[doubleCut[edgeI]]];
+                coeffI++;
+
+                // Neighbour side
+                localMult[doubleCutNeighbour[edgeI]] +=
+                    coeffs[coeffI]*psiInternal[L[doubleCut[edgeI]]];
                 coeffI++;
             }
         }
-    }
 
-    // Doubly cut coefficients
-    // ~~~~~~~~~~~~~~~~~~~~~~~
+        // Add the local multiplication to this side as well
 
-    // There exists a possibility of having an internal edge for a
-    // point on the processor patch which is in fact connected to
-    // another point of the same patch.  This particular nastiness
-    // introduces a deformation in the solution because the edge is
-    // either multiplied twice or not at all.  For this purpose, the
-    // offending edges need to be separated out and multiplied
-    // appropriately.  This will only happen for cell tetrahedral
-    // decomposition and is generally nasty.  
-    // No need for cut mask here
-    {
-        const labelList& doubleCut = procPatch_.doubleCutEdgeIndices();
-
-        const labelList& doubleCutOwner = procPatch_.doubleCutOwner();
-        const labelList& doubleCutNeighbour = procPatch_.doubleCutNeighbour();
-
-        forAll (doubleCut, edgeI)
+        forAll (mp, pointI)
         {
-            // Owner side
-            localMult[doubleCutOwner[edgeI]] +=
-                coeffs[coeffI]*psiInternal[U[doubleCut[edgeI]]];
-            coeffI++;
-
-            // Neighbour side
-            localMult[doubleCutNeighbour[edgeI]] +=
-                coeffs[coeffI]*psiInternal[L[doubleCut[edgeI]]];
-            coeffI++;
+            result[mp[pointI]] += localMult[pointI];
         }
     }
 
-    // Add the local multiplication to this side as well
-
-    forAll (mp, pointI)
-    {
-        result[mp[pointI]] += localMult[pointI];
-    }
-
     // Send the localMult
     sendField(tlocalMult, commsType);
 }
@@ -1099,9 +1197,13 @@ updateInterfaceMatrix
     const lduMatrix&,
     const scalarField&,
     const direction,
-    const Pstream::commsTypes commsType
+    const Pstream::commsTypes commsType,
+    const bool switchToLhs
 ) const
 {
+    // Switch to lhs handled in init
+    // HJ, 22/May/2013
+
     // Get the neighbour side multiplication
     tmp<scalarField> tneiMult = receivePointField<scalar>(commsType);
     this->addToInternalField(result, tneiMult());
diff --git a/src/OpenFOAM/fields/PointPatchFields/derived/global/GlobalPointPatchField.C b/src/OpenFOAM/fields/PointPatchFields/derived/global/GlobalPointPatchField.C
index ccedca4c1..e5c5e99f5 100644
--- a/src/OpenFOAM/fields/PointPatchFields/derived/global/GlobalPointPatchField.C
+++ b/src/OpenFOAM/fields/PointPatchFields/derived/global/GlobalPointPatchField.C
@@ -1024,7 +1024,8 @@ void GlobalPointPatchField
     const lduMatrix& m,
     const scalarField& coeffs,
     const direction,
-    const Pstream::commsTypes commsType
+    const Pstream::commsTypes commsType,
+    const bool switchToLhs
 ) const
 {
     tmp<scalarField> tlocalMult(new scalarField(this->size(), 0));
@@ -1047,125 +1048,261 @@ void GlobalPointPatchField
     label coeffI = 0;
     scalarField sumOffDiag(this->size(), 0);
 
-    // Owner side
-    // ~~~~~~~~~~
+    if (switchToLhs)
     {
-        const labelList& cutOwn = globalPointPatch_.cutEdgeOwnerIndices();
-        const labelList& cutOwnStart = globalPointPatch_.cutEdgeOwnerStart();
-
-        forAll (mp, pointI)
+        // Owner side
+        // ~~~~~~~~~~
         {
-            label ownIndex = cutOwnStart[pointI];
-            label endOwn = cutOwnStart[pointI + 1];
+            const labelList& cutOwn = globalPointPatch_.cutEdgeOwnerIndices();
+            const labelList& cutOwnStart =
+                globalPointPatch_.cutEdgeOwnerStart();
 
-            for (; ownIndex < endOwn; ownIndex++)
+            forAll (mp, pointI)
             {
-                localMult[pointI] +=
-                    cutMask[coeffI]*coeffs[coeffI]
-                    *psiInternal[U[cutOwn[ownIndex]]];
+                label ownIndex = cutOwnStart[pointI];
+                label endOwn = cutOwnStart[pointI + 1];
 
-                sumOffDiag[pointI] += cutMask[coeffI]*coeffs[coeffI];
+                for (; ownIndex < endOwn; ownIndex++)
+                {
+                    localMult[pointI] +=
+                        cutMask[coeffI]*coeffs[coeffI]
+                        *psiInternal[U[cutOwn[ownIndex]]];
 
-                // Multiply the internal side as well
-                result[U[cutOwn[ownIndex]]] +=
-                    coeffs[coeffI]*psiInternal[mp[pointI]];
+                    sumOffDiag[pointI] += cutMask[coeffI]*coeffs[coeffI];
+
+                    // Multiply the internal side as well
+                    result[U[cutOwn[ownIndex]]] -=
+                        coeffs[coeffI]*psiInternal[mp[pointI]];
+
+                    coeffI++;
+                }
+            }
+        }
+
+        // Neighbour side
+        // ~~~~~~~~~~~~~~
+        {
+            const labelList& cutNei =
+                globalPointPatch_.cutEdgeNeighbourIndices();
+            const labelList& cutNeiStart =
+                globalPointPatch_.cutEdgeNeighbourStart();
+
+            forAll (mp, pointI)
+            {
+                label neiIndex = cutNeiStart[pointI];
+                label endNei = cutNeiStart[pointI + 1];
+
+                for (; neiIndex < endNei; neiIndex++)
+                {
+                    localMult[pointI] +=
+                        cutMask[coeffI]*coeffs[coeffI]
+                        *psiInternal[L[cutNei[neiIndex]]];
+
+                    sumOffDiag[pointI] += cutMask[coeffI]*coeffs[coeffI];
+
+                    // Multiply the internal side as well
+                    result[L[cutNei[neiIndex]]] -=
+                        coeffs[coeffI]*psiInternal[mp[pointI]];
+
+                    coeffI++;
+                }
+            }
+        }
+
+        // Doubly cut coefficients
+        // ~~~~~~~~~~~~~~~~~~~~~~~
+
+        // There exists a possibility of having an internal edge for a
+        // point on the processor patch which is in fact connected to
+        // another point of the same patch.  This particular nastiness
+        // introduces a deformation in the solution because the edge is
+        // either multiplied twice or not at all.  For this purpose, the
+        // offending edges need to be separated out and multiplied
+        // appropriately.
+        {
+            const labelList& doubleCut =
+                globalPointPatch_.doubleCutEdgeIndices();
+
+            const labelList& doubleCutOwner =
+                globalPointPatch_.doubleCutOwner();
+
+            const labelList& doubleCutNeighbour =
+                globalPointPatch_.doubleCutNeighbour();
+
+            forAll (doubleCut, edgeI)
+            {
+                // Owner side
+                localMult[doubleCutOwner[edgeI]] +=
+                    cutMask[coeffI]*coeffs[coeffI]*
+                    psiInternal[U[doubleCut[edgeI]]];
+
+                sumOffDiag[doubleCutOwner[edgeI]] +=
+                    cutMask[coeffI]*coeffs[coeffI];
+
+                coeffI++;
+
+                // Neighbour side
+                localMult[doubleCutNeighbour[edgeI]] +=
+                    cutMask[coeffI]*coeffs[coeffI]*
+                    psiInternal[L[doubleCut[edgeI]]];
+
+                sumOffDiag[doubleCutNeighbour[edgeI]] +=
+                    cutMask[coeffI]*coeffs[coeffI];
 
                 coeffI++;
             }
         }
-    }
 
-    // Neighbour side
-    // ~~~~~~~~~~~~~~
-    {
-        const labelList& cutNei = globalPointPatch_.cutEdgeNeighbourIndices();
-        const labelList& cutNeiStart =
-            globalPointPatch_.cutEdgeNeighbourStart();
+        // Reduce/extract the result and enforce over all processors
 
-        forAll (mp, pointI)
+        // Requires global sync points to flush buffers before gather-scatter
+        // communications.  Reconsider.  HJ, 29/Mar/2011
+        if (Pstream::defaultCommsType == Pstream::nonBlocking)
         {
-            label neiIndex = cutNeiStart[pointI];
-            label endNei = cutNeiStart[pointI + 1];
+            IPstream::waitRequests();
+            OPstream::waitRequests();
+        }
 
-            for (; neiIndex < endNei; neiIndex++)
+        tmp<Field<scalar> > trpf =
+            reduceExtractPoint<scalar>(localMult);
+
+        Field<scalar>& rpf = trpf();
+
+        // Get addressing
+        const labelList& addr = globalPointPatch_.meshPoints();
+
+        forAll (addr, i)
+        {
+            result[addr[i]] -= rpf[i];
+        }
+    }
+    else
+    {
+        // Owner side
+        // ~~~~~~~~~~
+        {
+            const labelList& cutOwn = globalPointPatch_.cutEdgeOwnerIndices();
+            const labelList& cutOwnStart =
+                globalPointPatch_.cutEdgeOwnerStart();
+
+            forAll (mp, pointI)
             {
-                localMult[pointI] +=
-                    cutMask[coeffI]*coeffs[coeffI]
-                    *psiInternal[L[cutNei[neiIndex]]];
+                label ownIndex = cutOwnStart[pointI];
+                label endOwn = cutOwnStart[pointI + 1];
 
-                sumOffDiag[pointI] += cutMask[coeffI]*coeffs[coeffI];
+                for (; ownIndex < endOwn; ownIndex++)
+                {
+                    localMult[pointI] +=
+                        cutMask[coeffI]*coeffs[coeffI]
+                        *psiInternal[U[cutOwn[ownIndex]]];
 
-                // Multiply the internal side as well
-                result[L[cutNei[neiIndex]]] +=
-                    coeffs[coeffI]*psiInternal[mp[pointI]];
+                    sumOffDiag[pointI] += cutMask[coeffI]*coeffs[coeffI];
+
+                    // Multiply the internal side as well
+                    result[U[cutOwn[ownIndex]]] +=
+                        coeffs[coeffI]*psiInternal[mp[pointI]];
+
+                    coeffI++;
+                }
+            }
+        }
+
+        // Neighbour side
+        // ~~~~~~~~~~~~~~
+        {
+            const labelList& cutNei =
+                globalPointPatch_.cutEdgeNeighbourIndices();
+            const labelList& cutNeiStart =
+                globalPointPatch_.cutEdgeNeighbourStart();
+
+            forAll (mp, pointI)
+            {
+                label neiIndex = cutNeiStart[pointI];
+                label endNei = cutNeiStart[pointI + 1];
+
+                for (; neiIndex < endNei; neiIndex++)
+                {
+                    localMult[pointI] +=
+                        cutMask[coeffI]*coeffs[coeffI]
+                        *psiInternal[L[cutNei[neiIndex]]];
+
+                    sumOffDiag[pointI] += cutMask[coeffI]*coeffs[coeffI];
+
+                    // Multiply the internal side as well
+                    result[L[cutNei[neiIndex]]] +=
+                        coeffs[coeffI]*psiInternal[mp[pointI]];
+
+                    coeffI++;
+                }
+            }
+        }
+
+        // Doubly cut coefficients
+        // ~~~~~~~~~~~~~~~~~~~~~~~
+
+        // There exists a possibility of having an internal edge for a
+        // point on the processor patch which is in fact connected to
+        // another point of the same patch.  This particular nastiness
+        // introduces a deformation in the solution because the edge is
+        // either multiplied twice or not at all.  For this purpose, the
+        // offending edges need to be separated out and multiplied
+        // appropriately.
+        {
+            const labelList& doubleCut =
+                globalPointPatch_.doubleCutEdgeIndices();
+
+            const labelList& doubleCutOwner =
+                globalPointPatch_.doubleCutOwner();
+
+            const labelList& doubleCutNeighbour =
+                globalPointPatch_.doubleCutNeighbour();
+
+            forAll (doubleCut, edgeI)
+            {
+                // Owner side
+                localMult[doubleCutOwner[edgeI]] +=
+                    cutMask[coeffI]*coeffs[coeffI]*
+                    psiInternal[U[doubleCut[edgeI]]];
+
+                sumOffDiag[doubleCutOwner[edgeI]] +=
+                    cutMask[coeffI]*coeffs[coeffI];
+
+                coeffI++;
+
+                // Neighbour side
+                localMult[doubleCutNeighbour[edgeI]] +=
+                    cutMask[coeffI]*coeffs[coeffI]*
+                    psiInternal[L[doubleCut[edgeI]]];
+
+                sumOffDiag[doubleCutNeighbour[edgeI]] +=
+                    cutMask[coeffI]*coeffs[coeffI];
 
                 coeffI++;
             }
         }
-    }
 
-    // Doubly cut coefficients
-    // ~~~~~~~~~~~~~~~~~~~~~~~
+        // Reduce/extract the result and enforce over all processors
 
-    // There exists a possibility of having an internal edge for a
-    // point on the processor patch which is in fact connected to
-    // another point of the same patch.  This particular nastiness
-    // introduces a deformation in the solution because the edge is
-    // either multiplied twice or not at all.  For this purpose, the
-    // offending edges need to be separated out and multiplied
-    // appropriately.
-    {
-        const labelList& doubleCut = globalPointPatch_.doubleCutEdgeIndices();
-
-        const labelList& doubleCutOwner = globalPointPatch_.doubleCutOwner();
-        const labelList& doubleCutNeighbour =
-            globalPointPatch_.doubleCutNeighbour();
-
-        forAll (doubleCut, edgeI)
+        // Requires global sync points to flush buffers before gather-scatter
+        // communications.  Reconsider.  HJ, 29/Mar/2011
+        if (Pstream::defaultCommsType == Pstream::nonBlocking)
         {
-            // Owner side
-            localMult[doubleCutOwner[edgeI]] +=
-                cutMask[coeffI]*coeffs[coeffI]*
-                psiInternal[U[doubleCut[edgeI]]];
-
-            sumOffDiag[doubleCutOwner[edgeI]] +=
-                cutMask[coeffI]*coeffs[coeffI];
-
-            coeffI++;
-
-            // Neighbour side
-            localMult[doubleCutNeighbour[edgeI]] +=
-                cutMask[coeffI]*coeffs[coeffI]*
-                psiInternal[L[doubleCut[edgeI]]];
-
-            sumOffDiag[doubleCutNeighbour[edgeI]] +=
-                cutMask[coeffI]*coeffs[coeffI];
-
-            coeffI++;
+            IPstream::waitRequests();
+            OPstream::waitRequests();
         }
-    }
 
-    // Reduce/extract the result and enforce over all processors
+        tmp<Field<scalar> > trpf =
+            reduceExtractPoint<scalar>(localMult);
 
-    // Requires global sync points to flush buffers before gather-scatter
-    // communications.  Reconsider.  HJ, 29/Mar/2011
-    if (Pstream::defaultCommsType == Pstream::nonBlocking)
-    {
-        IPstream::waitRequests();
-        OPstream::waitRequests();
-    }
+        Field<scalar>& rpf = trpf();
 
-    tmp<Field<scalar> > trpf =
-        reduceExtractPoint<scalar>(localMult);
+        // Get addressing
+        const labelList& addr = globalPointPatch_.meshPoints();
 
-    Field<scalar>& rpf = trpf();
-
-    // Get addressing
-    const labelList& addr = globalPointPatch_.meshPoints();
-
-    forAll (addr, i)
-    {
-        result[addr[i]] += rpf[i];
+        forAll (addr, i)
+        {
+            result[addr[i]] += rpf[i];
+        }
     }
 }
 
diff --git a/src/OpenFOAM/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolate.C b/src/OpenFOAM/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolate.C
index 4dd912548..9b8d3b3e8 100644
--- a/src/OpenFOAM/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolate.C
+++ b/src/OpenFOAM/interpolations/MixingPlaneInterpolation/MixingPlaneInterpolate.C
@@ -33,10 +33,6 @@ Contributor
 
 \*---------------------------------------------------------------------------*/
 
-#include <sys/time.h>
-#include <sstream>
-#include <iomanip>
-
 namespace Foam
 {
 
@@ -44,17 +40,13 @@ namespace Foam
 
 template<class MasterPatch, class SlavePatch>
 template<class Type>
-void
-MixingPlaneInterpolation<MasterPatch, SlavePatch>::
-interpolate
+void MixingPlaneInterpolation<MasterPatch, SlavePatch>::toProfile
 (
     const Field<Type>& srcF,
     const labelListList& srcAddr,
     const scalarListList& srcWeights,
-    const labelListList& dstAddr,
-    const scalarListList& dstWeights,
-    Field<Type>& dstResultF
-) const
+    Field<Type>& profileBandValues
+)
 {
     // The src to profile transfer is done using weighted averaging
     // evaluation of srcF.
@@ -76,121 +68,133 @@ interpolate
     //
     //         average(phi)  ==  sum(w*phi)
 
-
-    int nbrProfileBands = interpolationProfile_.size() - 1;
-
-    List<Type> profileBandValues(nbrProfileBands, pTraits<Type>::zero);
-    scalarField srcScalingValues(nbrProfileBands, 0.0);
-
     forAll (srcAddr, bandI)
     {
         forAll (srcAddr[bandI], faceI)
         {
             profileBandValues[bandI] +=
                 srcF[srcAddr[bandI][faceI]]*srcWeights[bandI][faceI];
-
-            // NB: The next operation should be computed only
-            // once... and should sum up to 1.0 Let's keep this
-            // operation for now, until the mixingPlane interface is
-            // fully debugged (MB, 07/2010)
-            srcScalingValues[bandI] += srcWeights[bandI][faceI];
-
-            if (debug <= -200)
-            {
-                Info << "bande: " << bandI
-                    << "  src valeur: " << srcF[srcAddr[bandI][faceI]]
-                    << "  src weight: " << srcWeights[bandI][faceI] << endl;
-            }
         }
     }
+}
 
-    // We don't need to divide the profileBandValues by the
-    // srcScalingValues because the srcScalingValues are identically
-    // equal to 1.0, thanks to the conservativeness of the GGI
-    // weighting factors
-    //profileBandValues = profileBandValues/srcScalingValues;
 
-    // profileBandValues are now the circumferentially averaged values
-
-    // The profile to dst transfer is done by simply distributing the
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void MixingPlaneInterpolation<MasterPatch, SlavePatch>::fromProfile
+(
+    const Field<Type>& profileBandValues,
+    const labelListList& dstAddr,
+    const scalarListList& dstWeights,
+    Field<Type>& dstResultF
+)
+{
+    // The profile to dst transfer is done by distributing the
     // average value accordingly to the dst weighting factors
     forAll (dstAddr, faceI)
     {
         const labelList& curAddr = dstAddr[faceI];
         const scalarList& curW = dstWeights[faceI];
 
+        dstResultF[faceI] = pTraits<Type>::zero;
+
         forAll (curAddr, bandI)
         {
             dstResultF[faceI] += profileBandValues[curAddr[bandI]]*curW[bandI];
-
-            if (debug <= -200)
-            {
-                Info<< "bande: " << dstAddr[faceI][bandI]
-                    << "  dst valeur: " << dstResultF[faceI]
-                    << "  dst weight: " << dstWeights[faceI][bandI] << endl;
-            }
         }
     }
+}
 
-    if (debug <= -500)
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void MixingPlaneInterpolation<MasterPatch, SlavePatch>::maskedFromProfile
+(
+    const Field<Type>& profileBandValues,
+    const labelListList& dstAddr,
+    const scalarListList& dstWeights,
+    Field<Type>& dstResultF,
+    const labelList& mask
+)
+{
+    // The profile to dst transfer is done by distributing the
+    // average value accordingly to the dst weighting factors
+    forAll (mask, maskI)
     {
-        error::printStack(Info);
+        // Pick the masked face
+        const label faceI = mask[maskI];
 
-        Info<< "srcF      :  " << srcF << nl
-            << "srcAddr   :  " << srcAddr << nl
-            << "srcWeights:  " << srcWeights << nl
-            << "profileBandValues:  " << profileBandValues << nl
-            << "srcScalingValues:  " << srcScalingValues << nl
-            << "dstAddr   :  " << dstAddr << nl
-            << "dstWeights:  " << dstWeights << nl
-            << "dstResultF:  " << dstResultF << nl
-            << "srcScalingValues: " << srcScalingValues << endl;
-    }
+        const labelList& curAddr = dstAddr[faceI];
+        const scalarList& curW = dstWeights[faceI];
 
-    if (debug <= -999)
-    {
-        fileName traceFileDir("./mixingPlaneTraceFiles");
+        dstResultF[maskI] = pTraits<Type>::zero;
 
-        if (!exists(traceFileDir))
+        forAll (curAddr, bandI)
         {
-            mkDir(traceFileDir);
+            dstResultF[maskI] += profileBandValues[curAddr[bandI]]*curW[bandI];
         }
-
-        struct timeval tod;
-        gettimeofday(&tod, NULL);
-
-        //struct timespec tp;
-        //clock_gettime(CLOCK_MONOTONIC,  &tp);
-        std::ostringstream osBuffer;
-
-        osBuffer
-            << Foam::name(tod.tv_sec)
-            << "." << std::setfill('0') << std::setw(6)
-            <<  Foam::name(tod.tv_usec);
-
-        fileName traceFileName(traceFileDir/"profileValues_" + osBuffer.str());
-        OFstream dumpFileSrc(traceFileName + "_orig");
-        OFstream dumpFileDst(traceFileName + "_interpolated");
-        OFstream dumpFileProfile(traceFileName + "_profile");
-
-        //Foam::error::printStack(Info);
-
-        InfoIn
-        (
-            "MixingPlaneInterpolation::interpolate"
-        )   << "Dumping src profiles to: " << traceFileName + "_orig" << nl
-            << "Dumping interpolated profiles to : "
-            << traceFileName + "_interpolated" << nl
-            << "Dumping profile values to: "
-            << traceFileName + "_profile" << endl;
-
-        dumpFileDst << dstResultF << endl;
-        dumpFileSrc << srcF << endl;
-        dumpFileProfile << profileBandValues << endl;
     }
 }
 
 
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void MixingPlaneInterpolation<MasterPatch, SlavePatch>::maskedTransform
+(
+    Field<Type>& transField,
+    const tensorField& t,
+    const Field<Type>& inField,
+    const labelList& mask
+)
+{
+    // The profile to dst transfer is done by distributing the
+    // average value accordingly to the dst weighting factors
+    forAll (mask, maskI)
+    {
+        // Pick the masked face
+        const label faceI = mask[maskI];
+
+        transField[maskI] = transform(t[faceI], inField[maskI]);
+    }
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void MixingPlaneInterpolation<MasterPatch, SlavePatch>::interpolate
+(
+    const Field<Type>& srcF,
+    const labelListList& srcAddr,
+    const scalarListList& srcWeights,
+    const labelListList& dstAddr,
+    const scalarListList& dstWeights,
+    Field<Type>& dstResultF
+) const
+{
+    Field<Type> profileBandValues(nProfileBands(), pTraits<Type>::zero);
+
+    // Interpolate from patch to profile
+    toProfile
+    (
+        srcF,
+        srcAddr,
+        srcWeights,
+        profileBandValues
+    );
+
+    // profileBandValues are now the circumferentially averaged values
+
+    // Collect from profile to patch
+    fromProfile
+    (
+        profileBandValues,
+        dstAddr,
+        dstWeights,
+        dstResultF
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class MasterPatch, class SlavePatch>
@@ -206,7 +210,7 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToSlave
         FatalErrorIn
         (
             "MixingPlaneInterpolation::masterToSlave("
-            "const Field<Type> ff)"
+            "const Field<Type> ff) const"
         )   << "given field does not correspond to patch. Patch size: "
             << masterPatch_.size() << " field size: " << patchFF.size()
             << abort(FatalError);
@@ -216,16 +220,6 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToSlave
     // MB: We need this back
     Field<Type> profileFF = transform(masterPatchToProfileT(), patchFF);
 
-    if (debug > 1)
-    {
-        Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToSlave: "
-            << "patchFF: " << patchFF << endl
-            << "profileFF: " << profileFF << endl
-            << "masterPatchToProfileT(): " << masterPatchToProfileT() << endl
-            << "slaveProfileToPatchT():  " << slaveProfileToPatchT() << endl
-            << endl;
-    }
-
     // Do interpolation
     tmp<Field<Type> > tresult
     (
@@ -243,14 +237,14 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToSlave
         profileFF,                     // Master data in 'profile space'
         masterPatchToProfileAddr(),    // From master: compute the average
         masterPatchToProfileWeights(),
-        slaveProfileToPatchAddr(),     // To slave we distribute the average from
+        slaveProfileToPatchAddr(),     // To slave distribute average from
         slaveProfileToPatchWeights(),  // profile to patch
         result
     );
 
     // Apply transform to bring the slave field back from 'profile space'
     // to 'patch space'
-    transform(result, slaveProfileToPatchT(), result); // MB: We need this back
+    transform(result, slaveProfileToPatchT(), result);
 
     return tresult;
 }
@@ -283,7 +277,7 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToMaster
         FatalErrorIn
         (
             "MixingPlaneInterpolation::slaveToMaster("
-            "const Field<Type> ff)"
+            "const Field<Type> ff) const"
         )   << "given field does not correspond to patch. Patch size: "
             << slavePatch_.size() << " field size: " << patchFF.size()
             << abort(FatalError);
@@ -292,16 +286,6 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToMaster
     // Move slave data from 'patch space' to 'profile space'
      Field<Type> profileFF = transform(slavePatchToProfileT(), patchFF);
 
-    if (debug > 1)
-    {
-        Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToMaster: "
-            << "patchFF: " << patchFF << endl
-            << "profileFF: " << profileFF << endl
-            << "slavePatchToProfileT(): " << slavePatchToProfileT() << endl
-            << "masterProfileToPatchT(): " << masterProfileToPatchT() << endl
-            << endl;
-    }
-
     // Do interpolation
     tmp<Field<Type> > tresult
     (
@@ -345,6 +329,334 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToMaster
     return tint;
 }
 
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+tmp<Field<Type> >
+MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToProfile
+(
+    const Field<Type>& patchFF
+) const
+{
+    if (patchFF.size() != masterPatch_.size())
+    {
+        FatalErrorIn
+        (
+            "MixingPlaneInterpolation::masterToProfile("
+            "const Field<Type> ff) const"
+        )   << "given field does not correspond to patch. Patch size: "
+            << masterPatch_.size() << " field size: " << patchFF.size()
+            << abort(FatalError);
+    }
+
+    // Move master data from 'patch space' to 'profile space'
+    Field<Type> profileFF = transform(masterPatchToProfileT(), patchFF);
+
+    // Do interpolation
+    tmp<Field<Type> > tresult
+    (
+        new Field<Type>
+        (
+            nProfileBands(),
+            pTraits<Type>::zero
+        )
+    );
+
+    Field<Type>& result = tresult();
+        
+    toProfile
+    (
+        profileFF,                      // Master data in 'profile space'
+        masterPatchToProfileAddr(),     // From master: compute the average
+        masterPatchToProfileWeights(),
+        result
+    );
+        
+    return tresult;
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+tmp<Field<Type> >
+MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToProfile
+(
+    const tmp<Field<Type> >& tff
+) const
+{
+    tmp<Field<Type> > tint = masterToProfile(tff());
+    tff.clear();
+    return tint;
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+tmp<Field<Type> >
+MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToProfile
+(
+    const Field<Type>& patchFF
+) const
+{
+    if (patchFF.size() != slavePatch_.size())
+    {
+        FatalErrorIn
+        (
+            "MixingPlaneInterpolation::slaveToProfile("
+            "const Field<Type> ff) const"
+        )   << "given field does not correspond to patch. Patch size: "
+            << slavePatch_.size() << " field size: " << patchFF.size()
+            << abort(FatalError);
+    }
+
+    // Move slave data from 'patch space' to 'profile space'
+    Field<Type> profileFF = transform(slavePatchToProfileT(), patchFF);
+
+    // Do interpolation
+    tmp<Field<Type> > tresult
+    (
+        new Field<Type>
+        (
+            nProfileBands(),
+            pTraits<Type>::zero
+        )
+    );
+
+    Field<Type>& result = tresult();
+
+    toProfile
+    (
+        profileFF,                     // Slave data in 'profile space'
+        slavePatchToProfileAddr(),     // From slave: compute the average
+        slavePatchToProfileWeights(),
+        result
+    );
+        
+    return tresult;
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+tmp<Field<Type> >
+MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToProfile
+(
+    const tmp<Field<Type> >& tff
+) const
+{
+    tmp<Field<Type> > tint = slaveToProfile(tff());
+    tff.clear();
+    return tint;
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+tmp<Field<Type> >
+MixingPlaneInterpolation<MasterPatch, SlavePatch>::profileToMaster
+(
+    const Field<Type>& profileFF
+) const
+{
+    if (profileFF.size() != nProfileBands())
+    {
+        FatalErrorIn
+        (
+            "MixingPlaneInterpolation::profileToMaster("
+            "const Field<Type> ff) const"
+        )   << "given field does not correspond to profile.  Profile size: "
+            << nProfileBands() << " field size: " << profileFF.size()
+            << abort(FatalError);
+    }
+
+    // Do interpolation
+    tmp<Field<Type> > tresult
+    (
+        new Field<Type>
+        (
+            masterPatch_.size(),
+            pTraits<Type>::zero
+        )
+    );
+
+    Field<Type>& result = tresult();
+        
+    fromProfile
+    (
+        profileFF,                      // Master data in 'profile space'
+        masterProfileToPatchAddr(),     // To master: distribute the average
+        masterProfileToPatchWeights(),
+        result
+    );
+        
+    // Apply transform to bring the master field back from 'profile space'
+    // to 'patch space'
+    transform(result, masterProfileToPatchT(), result);
+
+    return tresult;
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+tmp<Field<Type> >
+MixingPlaneInterpolation<MasterPatch, SlavePatch>::profileToMaster
+(
+    const tmp<Field<Type> >& tff
+) const
+{
+    tmp<Field<Type> > tint = profileToMaster(tff());
+    tff.clear();
+    return tint;
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void MixingPlaneInterpolation<MasterPatch, SlavePatch>::maskedProfileToMaster
+(
+    const Field<Type>& profileFF,
+    Field<Type>& result,
+    const labelList& mask
+) const
+{
+    if (profileFF.size() != nProfileBands() || result.size() != mask.size())
+    {
+        FatalErrorIn
+        (
+            "bvoid MixingPlaneInterpolation<MasterPatch, SlavePatch>::"
+            "maskedProfileToMaster\n"
+            "(\n"
+            "    const Field<Type>& profileFF,\n"
+            "    Field<Type>& result,\n"
+            "    const labelList& mask\n"
+            ") const"
+        )   << "given field does not correspond to profile.  Profile size: "
+            << nProfileBands() << " field size: " << profileFF.size()
+            << " result size: " << result.size()
+            << " mask size: " << mask.size()
+            << abort(FatalError);
+    }
+
+    // Do interpolation
+    maskedFromProfile
+    (
+        profileFF,                      // Master data in 'profile space'
+        masterProfileToPatchAddr(),     // To master: distribute the average
+        masterProfileToPatchWeights(),
+        result,
+        mask
+    );
+        
+    // Apply transform to bring the master field back from 'profile space'
+    // to 'patch space'
+    maskedTransform(result, masterProfileToPatchT(), result, mask);
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+tmp<Field<Type> >
+MixingPlaneInterpolation<MasterPatch, SlavePatch>::profileToSlave
+(
+    const Field<Type>& profileFF
+) const
+{
+    if (profileFF.size() != nProfileBands())
+    {
+        FatalErrorIn
+        (
+            "MixingPlaneInterpolation::profileToSlave("
+            "const Field<Type> ff) const"
+        )   << "given field does not correspond to profile.  Profile size: "
+            << nProfileBands() << " field size: " << profileFF.size()
+            << abort(FatalError);
+    }
+
+    // Do interpolation
+    tmp<Field<Type> > tresult
+    (
+        new Field<Type>
+        (
+            slavePatch_.size(),
+            pTraits<Type>::zero
+        )
+    );
+
+    Field<Type>& result = tresult();
+
+    fromProfile
+    (
+        profileFF,                     // Slave data in 'profile space'
+        slaveProfileToPatchAddr(),     // To slave distribute average from
+        slaveProfileToPatchWeights(),  // profile to patch
+        result
+    );
+        
+    // Apply transform to bring the slave field back from 'profile space'
+    // to 'patch space'
+    transform(result, slaveProfileToPatchT(), result);
+
+    return tresult;
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+tmp<Field<Type> >
+MixingPlaneInterpolation<MasterPatch, SlavePatch>::profileToSlave
+(
+    const tmp<Field<Type> >& tff
+) const
+{
+    tmp<Field<Type> > tint = profileToSlave(tff());
+    tff.clear();
+    return tint;
+}
+
+
+template<class MasterPatch, class SlavePatch>
+template<class Type>
+void MixingPlaneInterpolation<MasterPatch, SlavePatch>::maskedProfileToSlave
+(
+    const Field<Type>& profileFF,
+    Field<Type>& result,
+    const labelList& mask
+) const
+{
+    if (profileFF.size() != nProfileBands() || result.size() != mask.size())
+    {
+        FatalErrorIn
+        (
+            "void MixingPlaneInterpolation<MasterPatch, SlavePatch>::"
+            "maskedProfileToSlave\n"
+            "(\n"
+            "    const Field<Type>& profileFF,\n"
+            "    Field<Type>& result,\n"
+            "    const labelList& mask\n"
+            ") const"
+        )   << "given field does not correspond to profile.  Profile size: "
+            << nProfileBands() << " field size: " << profileFF.size()
+            << " result size: " << result.size()
+            << " mask size: " << mask.size()
+            << abort(FatalError);
+    }
+
+    maskedFromProfile
+    (
+        profileFF,                     // Slave data in 'profile space'
+        slaveProfileToPatchAddr(),     // To slave distribute average from
+        slaveProfileToPatchWeights(),  // profile to patch
+        result,
+        mask
+    );
+        
+    // Apply transform to bring the slave field back from 'profile space'
+    // to 'patch space'
+    maskedTransform(result, slaveProfileToPatchT(), result, mask);
+}
+
+
 template<class MasterPatch, class SlavePatch>
 template<class Type>
 tmp<Field<Type> >
@@ -367,15 +679,6 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToMaster
     // Move master data from 'patch space' to 'profile space'
     Field<Type> profileFF = transform(masterPatchToProfileT(), patchFF);
 
-    if (debug > 1)
-    {
-        Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToMaster: "
-            << "patchFF: " << patchFF << endl
-            << "profileFF: " << profileFF << endl
-            << "masterPatchToProfileT(): " << masterPatchToProfileT() << endl
-            << endl;
-    }
-
     // Do interpolation
     tmp<Field<Type> > tresult
     (
@@ -393,7 +696,7 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToMaster
         profileFF,                      // Master data in 'profile space'
         masterPatchToProfileAddr(),     // From master: compute the average
         masterPatchToProfileWeights(),
-        masterProfileToPatchAddr(),    // To master: distribute the average
+        masterProfileToPatchAddr(),     // To master: distribute the average
         masterProfileToPatchWeights(),
         result
     );
@@ -405,6 +708,7 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToMaster
     return tresult;
 }
 
+
 template<class MasterPatch, class SlavePatch>
 template<class Type>
 tmp<Field<Type> >
@@ -418,6 +722,7 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::masterToMaster
     return tint;
 }
 
+
 template<class MasterPatch, class SlavePatch>
 template<class Type>
 tmp<Field<Type> >
@@ -440,15 +745,6 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToSlave
     // Move slave data from 'patch space' to 'profile space'
     Field<Type> profileFF = transform(slavePatchToProfileT(), patchFF);
 
-    if (debug > 1)
-    {
-        Info << "MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToSlave: "
-            << "patchFF: " << patchFF << endl
-            << "profileFF: " << profileFF << endl
-            << "slavePatchToProfileT(): " << slavePatchToProfileT() << endl
-            << endl;
-    }
-
     // Do interpolation
     tmp<Field<Type> > tresult
     (
@@ -462,14 +758,14 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToSlave
     Field<Type>& result = tresult();
 
     interpolate
-        (
-            profileFF,                     // Slave data in 'profile space'
-            slavePatchToProfileAddr(),     // From slave: compute the average
-            slavePatchToProfileWeights(),
-            slaveProfileToPatchAddr(),     // To slave: distribute the average
-            slaveProfileToPatchWeights(),
-            result
-        );
+    (
+        profileFF,                     // Slave data in 'profile space'
+        slavePatchToProfileAddr(),     // From slave: compute the average
+        slavePatchToProfileWeights(),
+        slaveProfileToPatchAddr(),     // To slave: distribute the average
+        slaveProfileToPatchWeights(),
+        result
+    );
         
     // Apply transform to bring the slave field back from 'profile space'
     // to 'patch space'
@@ -478,6 +774,7 @@ MixingPlaneInterpolation<MasterPatch, SlavePatch>::slaveToSlave
     return tresult;
 }
 
+
 template<class MasterPatch, class SlavePatch>
 template<class Type>
 tmp<Field<Type> >
diff --git a/src/OpenFOAM/interpolations/RBFInterpolation/RBFInterpolation.H b/src/OpenFOAM/interpolations/RBFInterpolation/RBFInterpolation.H
index 592a9a769..16e3f1579 100644
--- a/src/OpenFOAM/interpolations/RBFInterpolation/RBFInterpolation.H
+++ b/src/OpenFOAM/interpolations/RBFInterpolation/RBFInterpolation.H
@@ -29,8 +29,8 @@ Description
     Radial basis function interpolation class
 
 Description
-    Interpolation class which uses Radial Basis Functions to interpolate the
-    fluid displacements for given boundary displacements.
+    Interpolation class which uses Radial Basis Functions to interpolate
+    field from given data points to arbitrary set of points.
 
 	The coefficient vectors, alpha and beta are determined by solving
     the system:
@@ -38,25 +38,17 @@ Description
 	| db | = | Mbb Pb | | alpha |
 	| 0  |   | Pb  0  | |  beta |
 
-	where db are the given boundary displacements,
+	where db are the given field values at data carrier points.
 	Mbb the boundary RBF correlation matrix (NbxNb), containing RBF evaluations
     at the boundary nodes, and Pb some linear polynomial matrix (Nbx4).
 
-	Those coefficients are calculated every timestep, with the current
-    boundary displacements db, with the inverse of Mbb. Using those
-    coefficients, the RBF is evaluated at all fluid points every
-    timestep.
-
-	The efficiency of this method is increased by:
-		1) using control points which is a subset of the moving
-           boundary points. Those control points are selected by
-		   a coarsening function.
-		2) The outer boundary points are neglected since a cutoff function
-           is used toward the outer boundaries.
+    In cases where far field data is not of interest, a cutoff function
+    is used to eliminate unnecessary data points in the far field
 
 Author
     Frank Bos, TU Delft.  All rights reserved.
     Dubravko Matijasevic, FSB Zagreb.
+    Reorganisation by Hrvoje Jasak, Wikki Ltd.
 
 SourceFiles
     RBFInterpolation.C
@@ -90,7 +82,7 @@ class RBFInterpolation
         const vectorField& controlPoints_;
 
         //- Reference to all points
-        const vectorField& allPoints_;
+        const vectorField& dataPoints_;
 
         //- RBF function
         autoPtr<RBFFunction> RBF_;
@@ -136,7 +128,7 @@ public:
         (
             const dictionary& dict,
             const vectorField& controlPoints,
-            const vectorField& allPoints
+            const vectorField& dataPoints
         );
 
         //- Construct as copy
@@ -150,6 +142,18 @@ public:
 
     // Member Functions
 
+        //- Return reference to control points
+        const vectorField& controlPoints() const
+        {
+            return controlPoints_;
+        }
+
+        //- Reference to all points
+        const vectorField& dataPoints() const
+        {
+            return dataPoints_;
+        }
+
         //- Interpolate
         template<class Type>
         tmp<Field<Type> > interpolate(const Field<Type>& ctrlField) const;
diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockAmgSolver/BlockAmgSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockAmgSolver/BlockAmgSolver.C
index 61b6ab479..033c1f87a 100644
--- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockAmgSolver/BlockAmgSolver.C
+++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockAmgSolver/BlockAmgSolver.C
@@ -105,7 +105,7 @@ Foam::BlockAmgSolver<Type>::solve
             solverPerf.finalResidual() = gSum(cmptMag(wA))/norm;
             solverPerf.nIterations()++;
 
-        } while (!stop(solverPerf));
+        } while (!this->stop(solverPerf));
     }
 
     return solverPerf;
diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C
index 90090b51c..2decd7488 100644
--- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C
+++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockBiCGStab/BlockBiCGStabSolver.C
@@ -93,7 +93,7 @@ Foam::BlockBiCGStabSolver<Type>::solve
 
     // Check convergence, solve if not converged
 
-    if (!stop(solverPerf))
+    if (!this->stop(solverPerf))
     {
         scalar rho = this->great_;
         scalar rhoOld = rho;
diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C
index df80e23e0..757e9f269 100644
--- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C
+++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockCG/BlockCGSolver.C
@@ -92,7 +92,7 @@ typename Foam::BlockSolverPerformance<Type> Foam::BlockCGSolver<Type>::solve
 
     // Check convergence, solve if not converged
 
-    if (!stop(solverPerf))
+    if (!this->stop(solverPerf))
     {
         scalar rho = this->great_;
         scalar rhoOld = rho;
diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C
index 923eb4d92..b01a6198c 100644
--- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C
+++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGMRES/BlockGMRESSolver.C
@@ -128,7 +128,7 @@ Foam::BlockGMRESSolver<Type>::solve
 
     // Check convergence, solve if not converged
 
-    if (!stop(solverPerf))
+    if (!this->stop(solverPerf))
     {
         // Create the Hesenberg matrix
         scalarSquareMatrix H(nDirs_, 0);
diff --git a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C
index 31de86dc5..7700b0611 100644
--- a/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C
+++ b/src/OpenFOAM/matrices/blockLduMatrix/BlockLduSolvers/BlockGaussSeidel/BlockGaussSeidelSolver.C
@@ -86,7 +86,7 @@ Foam::BlockGaussSeidelSolver<Type>::solve
 
     // Check convergence, solve if not converged
 
-    if (!stop(solverPerf))
+    if (!this->stop(solverPerf))
     {
         // Iteration loop
 
diff --git a/src/engine/Make/files b/src/engine/Make/files
index 1ebef9aed..5c97e2acb 100644
--- a/src/engine/Make/files
+++ b/src/engine/Make/files
@@ -41,11 +41,6 @@ engineTopoChangerMesh/deformingEngineMesh/deformingEngineMeshInitialize.C
 engineTopoChangerMesh/deformingEngineMesh/deformingEngineMeshMove.C
 engineTopoChangerMesh/deformingEngineMesh/addDeformingEngineMeshZones.C
 
-engineTopoChangerMesh/deformingEngineMesh/deformingEngineMesh.C
-engineTopoChangerMesh/deformingEngineMesh/deformingEngineMeshInitialize.C
-engineTopoChangerMesh/deformingEngineMesh/deformingEngineMeshMove.C
-engineTopoChangerMesh/deformingEngineMesh/addDeformingEngineMeshZones.C
-
 engineTopoChangerMesh/engineValveSliding/addEngineValveSlidingMeshModifiers.C
 engineTopoChangerMesh/engineValveSliding/engineValveSliding.C
 engineTopoChangerMesh/engineValveSliding/engineValveSlidingInitialize.C
diff --git a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
index f7fe1b1fd..7d3c6a65d 100644
--- a/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
+++ b/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C
@@ -140,7 +140,7 @@ gaussGrad<Type>::grad
     GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad();
 
     gGrad.rename("grad(" + vsf.name() + ')');
-    correctBoundaryConditions(vsf, gGrad);
+    this->correctBoundaryConditions(vsf, gGrad);
 
     return tgGrad;
 }
diff --git a/wmake/rules/linux64Gcc/c++ b/wmake/rules/linux64Gcc/c++
index a4f5f54cb..767c0162e 100644
--- a/wmake/rules/linux64Gcc/c++
+++ b/wmake/rules/linux64Gcc/c++
@@ -2,7 +2,7 @@
 
 c++WARN     = -Wall -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor
 
-CC          = ccache g++ -m64
+CC          = g++ -m64
 
 include $(RULES)/c++$(WM_COMPILE_OPTION)