diff --git a/applications/test/IndirectList/IndirectListTest.C b/applications/test/IndirectList/IndirectListTest.C
index 242ff7dbce..365c58a383 100644
--- a/applications/test/IndirectList/IndirectListTest.C
+++ b/applications/test/IndirectList/IndirectListTest.C
@@ -43,6 +43,10 @@ int main(int argc, char *argv[])
         completeList[i] = 0.1*i;
     }
 
+    Info<< "raw : " << completeList << nl
+        << endl;
+
+
     List<label> addresses(5);
     addresses[0] = 1;
     addresses[1] = 0;
@@ -52,14 +56,23 @@ int main(int argc, char *argv[])
 
     IndirectList<double> idl(completeList, addresses);
 
-    forAll(idl, i)
-    {
-        Info<< idl[i] << token::SPACE;
-    }
+    Info<< "addr: " << idl.addressing() << nl
+        << "list: " << idl() << nl
+        << endl;
 
-    Info<< endl;
+    addresses[4] = 1;
+    addresses[3] = 0;
+    addresses[2] = 7;
+    addresses[1] = 8;
+    addresses[0] = 5;
 
-    Info << "\nEnd\n" << endl;
+    idl.resetAddressing(addresses.xfer());
+
+    Info<< "addr: " << idl.addressing() << nl
+        << "list: " << idl() << nl
+        << endl;
+
+    Info << "End\n" << endl;
 
     return 0;
 }
diff --git a/applications/test/tokenizeTest/Make/files b/applications/test/tokenizeTest/Make/files
new file mode 100644
index 0000000000..c2db390007
--- /dev/null
+++ b/applications/test/tokenizeTest/Make/files
@@ -0,0 +1,3 @@
+tokenizeTest.C
+
+EXE = $(FOAM_USER_APPBIN)/tokenizeTest
diff --git a/applications/test/tokenizeTest/Make/options b/applications/test/tokenizeTest/Make/options
new file mode 100644
index 0000000000..42da6e3418
--- /dev/null
+++ b/applications/test/tokenizeTest/Make/options
@@ -0,0 +1,7 @@
+/*
+  check for consistent behaviour with non-optimized code
+  */
+
+EXE_INC = \
+    -DFULLDEBUG -g -O0
+
diff --git a/applications/test/tokenizeTest/tokenizeTest.C b/applications/test/tokenizeTest/tokenizeTest.C
new file mode 100644
index 0000000000..deb675ee00
--- /dev/null
+++ b/applications/test/tokenizeTest/tokenizeTest.C
@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+
+Description
+    Test the tokenizing of various things
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "IOobject.H"
+#include "IOstreams.H"
+#include "IFstream.H"
+#include "IStringStream.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    argList::validArgs.insert("string .. stringN");
+    argList::validOptions.insert("file", "name");
+
+    argList args(argc, argv, false, true);
+
+    forAll(args.additionalArgs(), argI)
+    {
+        const string& rawArg = args.additionalArgs()[argI];
+        Info<< "input string: " << rawArg << nl;
+
+        IStringStream is(rawArg);
+        
+        while (is.good())
+        {
+            token tok(is);
+            Info<< "token: " << tok.info() << endl;
+        }
+        
+        Info<< nl;
+        IOobject::writeDivider(Info);
+    }  
+    
+    
+    if (args.optionFound("file"))
+    {
+        IFstream is(args.option("file"));
+        
+        Info<< "tokenizing file: " << args.option("file") << nl;
+
+        while (is.good())
+        {
+            token tok(is);
+            Info<< "token: " << tok.info() << endl;
+        }
+        
+        Info<< nl;
+        IOobject::writeDivider(Info);
+    }
+
+    return 0;
+}
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/conversion/Optional/ccm26ToFoam/ccm26ToFoam.C b/applications/utilities/mesh/conversion/Optional/ccm26ToFoam/ccm26ToFoam.C
index 47ec7f9137..f4bf559009 100644
--- a/applications/utilities/mesh/conversion/Optional/ccm26ToFoam/ccm26ToFoam.C
+++ b/applications/utilities/mesh/conversion/Optional/ccm26ToFoam/ccm26ToFoam.C
@@ -908,10 +908,10 @@ int main(int argc, char *argv[])
             runTime.constant(),
             runTime
         ),
-        foamPoints,
-        foamFaces,
-        foamOwner,
-        foamNeighbour
+        xferMove<pointField>(foamPoints),
+        xferMove<faceList>(foamFaces),
+        xferCopy<labelList>(foamOwner),
+        xferMove<labelList>(foamNeighbour)
     );
 
     // Create patches. Use patch types to determine what Foam types to generate.
diff --git a/src/OpenFOAM/containers/Lists/BiIndirectList/BiIndirectList.H b/src/OpenFOAM/containers/Lists/BiIndirectList/BiIndirectList.H
index a7827ec6cb..5afd48d0f3 100644
--- a/src/OpenFOAM/containers/Lists/BiIndirectList/BiIndirectList.H
+++ b/src/OpenFOAM/containers/Lists/BiIndirectList/BiIndirectList.H
@@ -69,6 +69,14 @@ public:
             const UList<label>&
         );
 
+        //- Construct given the complete list and by transferring addressing
+        inline BiIndirectList
+        (
+            const UList<T>& posList,
+            const UList<T>& negList,
+            const Xfer< List<label> >&
+        );
+
 
     // Member Functions
 
@@ -80,12 +88,17 @@ public:
             inline const UList<T>& posList() const;
             inline const UList<T>& negList() const;
             inline const List<label>& addressing() const;
-            inline List<label>& addressing();
 
             //- Calculate index given whether index is into posList or negList
             inline static label posIndex(const label);
             inline static label negIndex(const label);
 
+        // Edit
+
+            //- Reset addressing
+            inline void resetAddressing(const UList<label>&);
+            inline void resetAddressing(const Xfer< List<label> >&);
+
 
         // Member Operators
 
diff --git a/src/OpenFOAM/containers/Lists/BiIndirectList/BiIndirectListI.H b/src/OpenFOAM/containers/Lists/BiIndirectList/BiIndirectListI.H
index b30769d107..e6f33f709e 100644
--- a/src/OpenFOAM/containers/Lists/BiIndirectList/BiIndirectListI.H
+++ b/src/OpenFOAM/containers/Lists/BiIndirectList/BiIndirectListI.H
@@ -40,6 +40,20 @@ inline Foam::BiIndirectList<T>::BiIndirectList
 {}
 
 
+template<class T>
+inline Foam::BiIndirectList<T>::BiIndirectList
+(
+    const UList<T>& posList,
+    const UList<T>& negList,
+    const Xfer< List<label> >& addr
+)
+:
+    posList_(const_cast<UList<T>&>(posList)),
+    negList_(const_cast<UList<T>&>(negList)),
+    addressing_(addr)
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class T>
@@ -79,9 +93,22 @@ inline const Foam::List<Foam::label>& Foam::BiIndirectList<T>::addressing()
 
 
 template<class T>
-inline Foam::List<Foam::label>& Foam::BiIndirectList<T>::addressing()
+inline void Foam::BiIndirectList<T>::resetAddressing
+(
+    const UList<label>& addr
+)
 {
-    return addressing_;
+    addressing_ = addr;
+}
+
+
+template<class T>
+inline void Foam::BiIndirectList<T>::resetAddressing
+(
+    const Xfer< List<label> >& addr
+)
+{
+    addressing_.transfer(addr());
 }
 
 
diff --git a/src/OpenFOAM/containers/Lists/IndirectList/IndirectList.H b/src/OpenFOAM/containers/Lists/IndirectList/IndirectList.H
index 776d7b65db..cfdc03d7dc 100644
--- a/src/OpenFOAM/containers/Lists/IndirectList/IndirectList.H
+++ b/src/OpenFOAM/containers/Lists/IndirectList/IndirectList.H
@@ -63,6 +63,8 @@ public:
         //- Construct given the complete list and the addressing array
         inline IndirectList(const UList<T>&, const UList<label>&);
 
+        //- Construct given the complete list and by transferring addressing
+        inline IndirectList(const UList<T>&, const Xfer< List<label> >&);
 
     // Member Functions
 
@@ -77,7 +79,8 @@ public:
         // Edit
 
             //- Reset addressing
-            void resetAddressing(const UList<label>&);
+            inline void resetAddressing(const UList<label>&);
+            inline void resetAddressing(const Xfer< List<label> >&);
 
 
         // Member Operators
diff --git a/src/OpenFOAM/containers/Lists/IndirectList/IndirectListI.H b/src/OpenFOAM/containers/Lists/IndirectList/IndirectListI.H
index 0ae5311d97..a1e2bf9a54 100644
--- a/src/OpenFOAM/containers/Lists/IndirectList/IndirectListI.H
+++ b/src/OpenFOAM/containers/Lists/IndirectList/IndirectListI.H
@@ -38,6 +38,18 @@ inline Foam::IndirectList<T>::IndirectList
 {}
 
 
+template<class T>
+inline Foam::IndirectList<T>::IndirectList
+(
+    const UList<T>& completeList,
+    const Xfer< List<label> >& addr
+)
+:
+    completeList_(const_cast<UList<T>&>(completeList)),
+    addressing_(addr)
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class T>
@@ -78,6 +90,16 @@ inline void Foam::IndirectList<T>::resetAddressing
 }
 
 
+template<class T>
+inline void Foam::IndirectList<T>::resetAddressing
+(
+    const Xfer< List<label> >& addr
+)
+{
+    addressing_.transfer(addr());
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class T>
diff --git a/src/OpenFOAM/db/IOstreams/Sstreams/ISstream.C b/src/OpenFOAM/db/IOstreams/Sstreams/ISstream.C
index ba7635dfcd..21057c3f39 100644
--- a/src/OpenFOAM/db/IOstreams/Sstreams/ISstream.C
+++ b/src/OpenFOAM/db/IOstreams/Sstreams/ISstream.C
@@ -176,84 +176,33 @@ Foam::Istream& Foam::ISstream::read(token& t)
         //
         // ideally match the equivalent of this regular expression
         //
-        //    /^[-+]?([0-9]+\.?[0-9]*|\.[0-9]+)([Ee][-+]?[0-9]+)?$/
+        //    /[-+]?([0-9]+\.?[0-9]*|\.[0-9]+)([Ee][-+]?[0-9]+)?/
         //
         case '-' :
         case '.' :
         case '0' : case '1' : case '2' : case '3' : case '4' :
         case '5' : case '6' : case '7' : case '8' : case '9' :
         {
-            // has a digit
-            bool hasDigit = isdigit(c);
-
-            // has contents that cannot be label
-            bool notLabel = (c == '.');
-
-            // has contents that cannot be scalar
-            bool notScalar = false;
+            bool asLabel = (c != '.');
 
             unsigned int nChar = 0;
             charBuffer[nChar++] = c;
 
-            // the location of the last '[Ee]' exponent
-            unsigned int exponent = 0;
-
-            while (is_.get(c))
+            // get everything that could reasonable look like a number
+            while
+            (
+                is_.get(c)
+             && (
+                    isdigit(c)
+                 || c == '+'
+                 || c == '-'
+                 || c == '.'
+                 || c == 'E'
+                 || c == 'e'
+                )
+            )
             {
-                if (isdigit(c))
-                {
-                    hasDigit = true;
-                }
-                else if (isalpha(c))
-                {
-                    notLabel = true;
-
-                    if (c == 'E' || c == 'e')
-                    {
-                        if (exponent || !hasDigit)
-                        {
-                            // mantissa had no digits,
-                            // or already saw '[Ee]' before
-                            notScalar = true;
-                        }
-
-                        // remember this location
-                        exponent = nChar;
-                    }
-                    else
-                    {
-                        notScalar = true;
-                    }
-                }
-                else if (c == '+' || c == '-')
-                {
-                    notLabel = true;
-
-                    // only allowed once in exponent
-                    if (!exponent || exponent+1 != nChar)
-                    {
-                        notScalar = true;
-                    }
-                    else
-                    {
-                        // require some digits again
-                        hasDigit = false;
-                    }
-                }
-                else if (c == '.')
-                {
-                    // notLabel means we already saw '.' or '[Ee]' before
-                    // cannot have '.' again
-                    if (notLabel)
-                    {
-                        notScalar = true;
-                    }
-                    notLabel = true;
-                }
-                else
-                {
-                    break;
-                }
+                asLabel = asLabel && isdigit(c);
 
                 charBuffer[nChar++] = c;
                 if (nChar >= sizeof(charBuffer))
@@ -265,13 +214,7 @@ Foam::Istream& Foam::ISstream::read(token& t)
             }
             charBuffer[nChar] = '\0';
 
-            if (!hasDigit)
-            {
-                notLabel = notScalar = true;
-            }
-
             setState(is_.rdstate());
-
             if (!is_.bad())
             {
                 is_.putback(c);
@@ -281,32 +224,31 @@ Foam::Istream& Foam::ISstream::read(token& t)
                     // a single '-' is punctuation
                     t = token::punctuationToken(token::SUBTRACT);
                 }
-                else if (notScalar)
-                {
-                    // not label or scalar: must be a word
-                    t = new word(charBuffer);
-                }
-                else if (notLabel)
-                {
-                    // not label: must be a scalar
-                    t = scalar(atof(charBuffer));
-                }
-                else if (hasDigit)
-                {
-                    // has digit: treat as a label
-                    long lt = atol(charBuffer);
-                    t = label(lt);
-
-                    // return as a scalar if doesn't fit in a label
-                    if (t.labelToken() != lt)
-                    {
-                        t = scalar(atof(charBuffer));
-                    }
-                }
                 else
                 {
-                    // some else: must be a word
-                    t = new word(charBuffer);
+                    char *endptr;
+
+                    if (asLabel)
+                    {
+                        long longval = strtol(charBuffer, &endptr, 10);
+                        t = label(longval);
+
+                        // return as a scalar if doesn't fit in a label
+                        if (t.labelToken() != longval)
+                        {
+                            t = scalar(strtod(charBuffer, &endptr));
+                        }
+                    }
+                    else
+                    {
+                        t = scalar(strtod(charBuffer, &endptr));
+                    }
+
+                    // nothing converted (bad format), or trailing junk
+                    if (endptr == charBuffer || *endptr != '\0')
+                    {
+                        t.setBad();
+                    }
                 }
             }
             else
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 3534519ca7..29841d3ec8 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -626,7 +626,8 @@ bool Foam::primitiveMesh::checkFaceSkewness
         vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
 
         // Skewness vector
-        vector sv = Cpf - ((fAreas[faceI] & Cpf)/(fAreas[faceI] & d))*d;
+        vector sv =
+            Cpf - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + SMALL))*d;
         vector svHat = sv/(mag(sv) + VSMALL);
 
         // Normalisation distance calculated as the approximate distance
diff --git a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C
index 01b1a56d03..c38e0494ae 100644
--- a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -25,8 +25,6 @@ License
 InClass
     decompositionMethod
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "decompositionMethod.H"
@@ -104,14 +102,26 @@ Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
 }
 
 
+Foam::labelList Foam::decompositionMethod::decompose
+(
+    const pointField& points
+)
+{
+    scalarField weights(0);
+
+    return decompose(points, weights);
+}
+
+
 Foam::labelList Foam::decompositionMethod::decompose
 (
     const labelList& fineToCoarse,
-    const pointField& coarsePoints
+    const pointField& coarsePoints,
+    const scalarField& coarseWeights
 )
 {
     // Decompose based on agglomerated points
-    labelList coarseDistribution(decompose(coarsePoints));
+    labelList coarseDistribution(decompose(coarsePoints, coarseWeights));
 
     // Rework back into decomposition for original mesh_
     labelList fineDistribution(fineToCoarse.size());
@@ -125,6 +135,18 @@ Foam::labelList Foam::decompositionMethod::decompose
 }
 
 
+Foam::labelList Foam::decompositionMethod::decompose
+(
+    const labelList& fineToCoarse,
+    const pointField& coarsePoints
+)
+{
+    scalarField coarseWeights(0);
+
+    return decompose(fineToCoarse, coarsePoints, coarseWeights);
+}
+
+
 void Foam::decompositionMethod::calcCellCells
 (
     const polyMesh& mesh,
@@ -170,4 +192,16 @@ void Foam::decompositionMethod::calcCellCells
 }
 
 
+Foam::labelList Foam::decompositionMethod::decompose
+(
+    const labelListList& globalCellCells,
+    const pointField& cc
+)
+{
+    scalarField cWeights(0);
+
+    return decompose(globalCellCells, cc, cWeights);
+}
+
+
 // ************************************************************************* //
diff --git a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.H b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.H
index 7055f24d7b..7f66e9f50d 100644
--- a/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.H
+++ b/src/decompositionMethods/decompositionMethods/decompositionMethod/decompositionMethod.H
@@ -150,20 +150,37 @@ public:
 
         //- Return for every coordinate the wanted processor number. Use the
         //  mesh connectivity (if needed)
-        virtual labelList decompose(const pointField&) = 0;
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        ) = 0;
+
+        //- Like decompose but with uniform weights on the points
+        virtual labelList decompose(const pointField&);
+
 
         //- Return for every coordinate the wanted processor number. Gets
         //  passed agglomeration map (from fine to coarse cells) and coarse cell
         //  location. Can be overridden by decomposers that provide this
         //  functionality natively. Coarse cells are local to the processor
         //  (if in parallel). If you want to have coarse cells spanning
-        //  processors use the next function below instead.
+        //  processors use the globalCellCells instead.
+        virtual labelList decompose
+        (
+            const labelList& cellToRegion,
+            const pointField& regionPoints,
+            const scalarField& regionWeights
+        );
+
+        //- Like decompose but with uniform weights on the regions
         virtual labelList decompose
         (
             const labelList& cellToRegion,
             const pointField& regionPoints
         );
 
+
         //- Return for every coordinate the wanted processor number. Explicitly
         //  provided connectivity - does not use mesh_.
         //  The connectivity is equal to mesh.cellCells() except for
@@ -174,9 +191,17 @@ public:
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc
+            const pointField& cc,
+            const scalarField& cWeights
         ) = 0;
 
+        //- Like decompose but with uniform weights on the cells
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc
+        );
+
 };
 
 
diff --git a/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C b/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
index 72d39d215c..bde1181c09 100644
--- a/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
@@ -142,6 +142,38 @@ Foam::label Foam::hierarchGeomDecomp::findLower
 }
 
 
+// Create a mapping between the index and the weighted size.
+// For convenience, sortedWeightedSize is one size bigger than current. This
+// avoids extra tests.
+void Foam::hierarchGeomDecomp::calculateSortedWeightedSizes
+(
+    const labelList& current,
+    const labelList& indices,
+    const scalarField& weights,
+    const label globalCurrentSize,
+
+    scalarField& sortedWeightedSizes
+)
+{
+    // Evaluate cumulative weights.
+    sortedWeightedSizes[0] = 0;
+    forAll(current, i)
+    {
+        label pointI = current[indices[i]];
+        sortedWeightedSizes[i + 1] = sortedWeightedSizes[i] + weights[pointI];
+    }
+    // Non-dimensionalise and multiply by size.
+    scalar globalCurrentLength = returnReduce
+    (
+        sortedWeightedSizes[current.size()], 
+        sumOp<scalar>()
+    );
+    // Normalise weights by global sum of weights and multiply through
+    // by global size.
+    sortedWeightedSizes *= (globalCurrentSize/globalCurrentLength);
+}
+
+
 // Find position in values so between minIndex and this position there
 // are wantedSize elements.
 void Foam::hierarchGeomDecomp::findBinary
@@ -170,7 +202,6 @@ void Foam::hierarchGeomDecomp::findBinary
     label midPrev = -1;
     label highPrev = -1;
 
-    //while (low <= high)
     while (true)
     {
         label size = returnReduce(mid-minIndex, sumOp<label>());
@@ -204,7 +235,98 @@ void Foam::hierarchGeomDecomp::findBinary
         mid = findLower(values, midValue, low, high);
 
         // Safeguard if same as previous.
-        bool hasNotChanged = (mid == midPrev) && (low == lowPrev) && (high == highPrev);
+        bool hasNotChanged =
+            (mid == midPrev)
+         && (low == lowPrev)
+         && (high == highPrev);
+
+        if (returnReduce(hasNotChanged, andOp<bool>()))
+        {
+            WarningIn("hierarchGeomDecomp::findBinary(..)")
+                << "unable to find desired deomposition split, making do!" 
+                << endl;
+            break;
+        }
+
+        midPrev = mid;
+        lowPrev = low;
+        highPrev = high;
+    }
+}
+
+
+// Find position in values so between minIndex and this position there
+// are wantedSize elements.
+void Foam::hierarchGeomDecomp::findBinary
+(
+    const label sizeTol,
+    const List<scalar>& sortedWeightedSizes,
+    const List<scalar>& values,
+    const label minIndex,       // index of previous value
+    const scalar minValue,      // value at minIndex
+    const scalar maxValue,      // global max of values
+    const scalar wantedSize,    // wanted size
+
+    label& mid,                 // index where size of bin is
+                                // wantedSize (to within sizeTol)
+    scalar& midValue            // value at mid
+)
+{
+    label low = minIndex;
+    scalar lowValue = minValue;
+
+    scalar highValue = maxValue;
+    // (one beyond) index of highValue
+    label high = values.size();
+
+    // Safeguards to avoid infinite loop.
+    label lowPrev = -1;
+    label midPrev = -1;
+    label highPrev = -1;
+
+    while (true)
+    {
+        label weightedSize = returnReduce
+        (
+            sortedWeightedSizes[mid] - sortedWeightedSizes[minIndex],
+            sumOp<label>()
+        );
+
+        if (debug)
+        {
+            Pout<< "low:" << low << " lowValue:" << lowValue
+                << " high:" << high << " highValue:" << highValue
+                << " mid:" << mid << " midValue:" << midValue << nl
+                << "globalSize:" << weightedSize 
+                << " wantedSize:" << wantedSize
+                << " sizeTol:" << sizeTol << endl;
+        }
+
+        if (wantedSize < weightedSize - sizeTol)
+        {
+            high = mid;
+            highValue = midValue;
+        }
+        else if (wantedSize > weightedSize + sizeTol)
+        {
+            low = mid;
+            lowValue = midValue;
+        }
+        else
+        {
+            break;
+        }
+
+        // Update mid, midValue
+        midValue = 0.5*(lowValue+highValue);
+        mid = findLower(values, midValue, low, high);
+
+        // Safeguard if same as previous.
+        bool hasNotChanged =
+            (mid == midPrev)
+         && (low == lowPrev)
+         && (high == highPrev);
+
         if (returnReduce(hasNotChanged, andOp<bool>()))
         {
             WarningIn("hierarchGeomDecomp::findBinary(..)")
@@ -394,6 +516,190 @@ void Foam::hierarchGeomDecomp::sortComponent
 }
 
 
+// Sort points into bins according to one component. Recurses to next component.
+void Foam::hierarchGeomDecomp::sortComponent
+(
+    const label sizeTol,
+    const scalarField& weights,
+    const pointField& points,
+    const labelList& current,       // slice of points to decompose
+    const direction componentIndex, // index in decompOrder_
+    const label mult,               // multiplication factor for finalDecomp
+    labelList& finalDecomp
+)
+{
+    // Current component
+    label compI = decompOrder_[componentIndex];
+
+    if (debug)
+    {
+        Pout<< "sortComponent : Sorting slice of size " << current.size()
+            << " in component " << compI << endl;
+    }
+
+    // Storage for sorted component compI
+    SortableList<scalar> sortedCoord(current.size());
+
+    forAll(current, i)
+    {
+        label pointI = current[i];
+
+        sortedCoord[i] = points[pointI][compI];
+    }
+    sortedCoord.sort();
+
+    label globalCurrentSize = returnReduce(current.size(), sumOp<label>());
+
+    // Now evaluate local cumulative weights, based on the sorting.
+    // Make one bigger than the nodes.
+    scalarField sortedWeightedSizes(current.size()+1, 0);
+    calculateSortedWeightedSizes
+    (
+        current,
+        sortedCoord.indices(),
+        weights,
+        globalCurrentSize,
+        sortedWeightedSizes
+    );
+
+    scalar minCoord = returnReduce
+    (
+        (
+            sortedCoord.size()
+          ? sortedCoord[0]
+          : GREAT
+        ),
+        minOp<scalar>()
+    );
+
+    scalar maxCoord = returnReduce
+    (
+        (
+            sortedCoord.size()
+          ? sortedCoord[sortedCoord.size()-1]
+          : -GREAT
+        ),
+        maxOp<scalar>()
+    );
+
+    if (debug)
+    {   
+        Pout<< "sortComponent : minCoord:" << minCoord
+            << " maxCoord:" << maxCoord << endl;
+    }
+
+    // starting index (in sortedCoord) of bin (= local)
+    label leftIndex = 0;
+    // starting value of bin (= global since coordinate)
+    scalar leftCoord = minCoord;
+
+    // Sort bins of size n
+    for (label bin = 0; bin < n_[compI]; bin++)
+    {
+        // Now we need to determine the size of the bin (dx). This is
+        // determined by the 'pivot' values - everything to the left of this
+        // value goes in the current bin, everything to the right into the next
+        // bins.
+
+        // Local number of elements
+        label localSize = -1;     // offset from leftOffset
+
+        // Value at right of bin (leftIndex+localSize-1)
+        scalar rightCoord = -GREAT;
+        
+        if (bin == n_[compI]-1)
+        {
+            // Last bin. Copy all.
+            localSize = current.size()-leftIndex;
+            rightCoord = maxCoord;                  // note: not used anymore
+        }
+        else
+        {
+            // For the current bin (starting at leftCoord) we want a rightCoord
+            // such that the sum of all weighted sizes are 
+            // globalCurrentLength/n_[compI].
+            // We have to iterate to obtain this.
+
+            label rightIndex = current.size();
+            rightCoord = maxCoord;
+
+            // Calculate rightIndex/rightCoord to have wanted size
+            findBinary
+            (
+                sizeTol,
+                sortedWeightedSizes,
+                sortedCoord,
+                leftIndex,
+                leftCoord,
+                maxCoord,
+                globalCurrentSize/n_[compI],  // wanted size
+
+                rightIndex,
+                rightCoord
+            );
+            localSize = rightIndex - leftIndex;
+        }
+
+        if (debug)
+        {
+            Pout<< "For component " << compI << ", bin " << bin
+                << " copying" << nl
+                << "from " << leftCoord << " at local index "
+                << leftIndex << nl
+                << "to " << rightCoord << " localSize:"
+                << localSize << nl
+                << endl;
+        }
+
+
+        // Copy localSize elements starting from leftIndex.
+        labelList slice(localSize);
+
+        forAll(slice, i)
+        {
+            label pointI = current[sortedCoord.indices()[leftIndex+i]];
+
+            // Mark point into correct bin
+            finalDecomp[pointI] += bin*mult;
+
+            // And collect for next sorting action
+            slice[i] = pointI;
+        }
+
+        // Sort slice in next component
+        if (componentIndex < 2)
+        {
+            string oldPrefix;
+            if (debug)
+            {
+                oldPrefix = Pout.prefix();
+                Pout.prefix() = "  " + oldPrefix;
+            }
+
+            sortComponent
+            (
+                sizeTol,
+                weights,
+                points,
+                slice,
+                componentIndex+1,
+                mult*n_[compI],     // Multiplier to apply to decomposition.
+                finalDecomp
+            );
+
+            if (debug)
+            {
+                Pout.prefix() = oldPrefix;
+            }
+        }
+
+        // Step to next bin.
+        leftIndex += localSize;
+        leftCoord = rightCoord;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::hierarchGeomDecomp::hierarchGeomDecomp
@@ -464,4 +770,47 @@ Foam::labelList Foam::hierarchGeomDecomp::decompose
 }
 
 
+Foam::labelList Foam::hierarchGeomDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& weights
+)
+{
+    // construct a list for the final result
+    labelList finalDecomp(points.size(), 0);
+
+    // Start off with every point sorted onto itself.
+    labelList slice(points.size());
+    forAll(slice, i)
+    {
+        slice[i] = i;
+    }
+
+    pointField rotatedPoints = rotDelta_ & points;
+
+
+    // Calculate tolerance of cell distribution. For large cases finding
+    // distibution to the cell exact would cause too many iterations so allow
+    // some slack.
+    label allSize = points.size();
+    reduce(allSize, sumOp<label>());
+
+    const label sizeTol = max(1, label(1E-3*allSize/nProcessors_));
+
+    // Sort recursive
+    sortComponent
+    (
+        sizeTol,
+        weights,
+        rotatedPoints,
+        slice,
+        0,              // Sort first component in decompOrder.
+        1,              // Offset for different x bins.
+        finalDecomp
+    );
+
+    return finalDecomp;
+}
+
+
 // ************************************************************************* //
diff --git a/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H b/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
index c2ed565211..932ca5dcaa 100644
--- a/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
+++ b/src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.H
@@ -80,6 +80,17 @@ class hierarchGeomDecomp
         //- Convert ordering string ("xyz") into list of components.
         void setDecompOrder();
 
+        //- Evaluates the weighted sizes for each sorted point.
+        static void calculateSortedWeightedSizes
+        (
+            const labelList& current,
+            const labelList& indices,
+            const scalarField& weights,
+            const label globalCurrentSize,
+
+            scalarField& sortedWeightedSizes
+        );
+
         //- Find index of t in list inbetween indices left and right
         static label findLower
         (
@@ -104,6 +115,22 @@ class hierarchGeomDecomp
             scalar& midValue            // value at mid
         );
 
+        //- Find midValue (at local index mid) such that the number of
+        //  elements between mid and leftIndex are (globally summed) the
+        //  wantedSize. Binary search.
+        static void findBinary
+        (
+            const label sizeTol,        // size difference considered acceptible
+            const List<scalar>& sortedWeightedSizes,
+            const List<scalar>&,
+            const label leftIndex,      // index of previous value
+            const scalar leftValue,     // value at leftIndex
+            const scalar maxValue,      // global max of values
+            const scalar wantedSize,    // wanted size
+            label& mid,                 // index where size of bin is wantedSize
+            scalar& midValue            // value at mid
+        );
+
         //- Recursively sort in x,y,z (or rather acc. to decompOrder_)
         void sortComponent
         (
@@ -115,6 +142,19 @@ class hierarchGeomDecomp
             labelList& finalDecomp          // overall decomposition
         );
 
+        //- Recursively sort in x,y,z (or rather acc. to decompOrder_)
+        //- using weighted points.
+        void sortComponent
+        (
+            const label sizeTol,
+            const scalarField& weights,
+            const pointField&,
+            const labelList& slice,         // slice of points to decompose
+            const direction componentIndex, // index in decompOrder_
+            const label prevMult,           // multiplication factor
+            labelList& finalDecomp          // overall decomposition
+        );
+
 
         //- Disallow default bitwise copy construct and assignment
         void operator=(const hierarchGeomDecomp&);
@@ -156,6 +196,14 @@ public:
 
         //- Return for every coordinate the wanted processor number. Use the
         //  mesh connectivity (if needed)
+        virtual labelList decompose
+        (
+            const pointField&,
+            const scalarField& weights
+        );
+
+        //- Without weights. Code for weighted decomposition is a bit complex
+        //  so kept separate for now.
         virtual labelList decompose(const pointField&);
 
         //- Return for every coordinate the wanted processor number. Explicitly
@@ -165,6 +213,16 @@ public:
         //    from 0 at processor0 and then incrementing all through the
         //    processors)
         //  - the connections are across coupled patches
+        virtual labelList decompose
+        (
+            const labelListList& globalCellCells,
+            const pointField& cc,
+            const scalarField& cWeights
+        )
+        {
+            return decompose(cc, cWeights);
+        }
+
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
diff --git a/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.C b/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.C
index fd2ff9562d..b852e1cf8d 100644
--- a/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.C
@@ -82,7 +82,11 @@ Foam::manualDecomp::manualDecomp
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::labelList Foam::manualDecomp::decompose(const pointField& points)
+Foam::labelList Foam::manualDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& pointWeights
+)
 {
     const polyMesh& mesh = *meshPtr_;
 
@@ -103,8 +107,10 @@ Foam::labelList Foam::manualDecomp::decompose(const pointField& points)
 
     if (finalDecomp.size() != points.size())
     {
-        FatalErrorIn("manualDecomp::decompose(const pointField&)")
-            << "Size of decomposition list does not correspond "
+        FatalErrorIn
+        (
+            "manualDecomp::decompose(const pointField&, const scalarField&)"
+        )   << "Size of decomposition list does not correspond "
             << "to the number of points.  Size: "
             << finalDecomp.size() << " Number of points: "
             << points.size()
@@ -115,8 +121,10 @@ Foam::labelList Foam::manualDecomp::decompose(const pointField& points)
 
     if (min(finalDecomp) < 0 || max(finalDecomp) > nProcessors_ - 1)
     {
-        FatalErrorIn("labelList manualDecomp::decompose(const pointField&)")
-            << "According to the decomposition, cells assigned to "
+        FatalErrorIn
+        (
+            "manualDecomp::decompose(const pointField&, const scalarField&)"
+        )   << "According to the decomposition, cells assigned to "
             << "impossible processor numbers.  Min processor = "
             << min(finalDecomp) << " Max processor = " << max(finalDecomp)
             << ".\n" << "Manual decomposition data read from file "
diff --git a/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.H b/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.H
index ed5ad4e8c0..0829f3e023 100644
--- a/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.H
+++ b/src/decompositionMethods/decompositionMethods/manualDecomp/manualDecomp.H
@@ -99,7 +99,11 @@ public:
 
         //- Return for every coordinate the wanted processor number. Use the
         //  mesh connectivity (if needed)
-        virtual labelList decompose(const pointField&);
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        );
 
         //- Return for every coordinate the wanted processor number. Explicitly
         //  provided connectivity - does not use mesh_.
@@ -111,10 +115,11 @@ public:
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc
+            const pointField& cc,
+            const scalarField& cWeights
         )
         {
-            return decompose(cc);
+            return decompose(cc, cWeights);
         }
 };
 
diff --git a/src/decompositionMethods/decompositionMethods/metisDecomp/metisDecomp.C b/src/decompositionMethods/decompositionMethods/metisDecomp/metisDecomp.C
index b4a643562f..3b2c5fb58c 100644
--- a/src/decompositionMethods/decompositionMethods/metisDecomp/metisDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/metisDecomp/metisDecomp.C
@@ -60,6 +60,7 @@ Foam::label Foam::metisDecomp::decompose
 (
     const List<int>& adjncy,
     const List<int>& xadj,
+    const scalarField& cWeights,
 
     List<int>& finalDecomp
 )
@@ -72,6 +73,8 @@ Foam::label Foam::metisDecomp::decompose
     // k-way: multi-level k-way
     word method("k-way");
 
+    int numCells = xadj.size()-1;
+
     // decomposition options. 0 = use defaults
     List<int> options(5, 0);
 
@@ -85,6 +88,30 @@ Foam::label Foam::metisDecomp::decompose
     // face weights (so on the edges of the dual)
     List<int> faceWeights;
 
+
+    // Check for externally provided cellweights and if so initialise weights
+    scalar maxWeights = gMax(cWeights);
+    if (cWeights.size() > 0)
+    {
+        if (cWeights.size() != numCells)
+        {
+            FatalErrorIn
+            (
+                "metisDecomp::decompose"
+                "(const pointField&, const scalarField&)"
+            )   << "Number of cell weights " << cWeights.size()
+                << " does not equal number of cells " << numCells
+                << exit(FatalError);
+        }
+        // Convert to integers.
+        cellWeights.setSize(cWeights.size());
+        forAll(cellWeights, i)
+        {
+            cellWeights[i] = int(1000000*cWeights[i]/maxWeights);
+        }
+    }
+
+
     // Check for user supplied weights and decomp options
     if (decompositionDict_.found("metisCoeffs"))
     {
@@ -161,52 +188,8 @@ Foam::label Foam::metisDecomp::decompose
                     << exit(FatalError);
             }
         }
-
-        //- faceWeights disabled. Only makes sense for cellCells from mesh.
-        //if (metisCoeffs.readIfPresent("faceWeightsFile", weightsFile))
-        //{
-        //    Info<< "metisDecomp : Using face-based weights." << endl;
-        //
-        //    IOList<int> weights
-        //    (
-        //        IOobject
-        //        (
-        //            weightsFile,
-        //            mesh_.time().timeName(),
-        //            mesh_,
-        //            IOobject::MUST_READ,
-        //            IOobject::AUTO_WRITE
-        //        )
-        //    );
-        //
-        //    if (weights.size() != adjncy.size()/2)
-        //    {
-        //        FatalErrorIn("metisDecomp::decompose(const pointField&)")
-        //            << "Number of face weights " << weights.size()
-        //            << " does not equal number of internal faces "
-        //            << adjncy.size()/2
-        //            << exit(FatalError);
-        //    }
-        //
-        //    // Assume symmetric weights. Keep same ordering as adjncy.
-        //    faceWeights.setSize(adjncy.size());
-        //
-        //    labelList nFacesPerCell(mesh_.nCells(), 0);
-        //
-        //    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-        //    {
-        //        label w = weights[faceI];
-        //
-        //        label own = mesh_.faceOwner()[faceI];
-        //        label nei = mesh_.faceNeighbour()[faceI];
-        //
-        //        faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
-        //        faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
-        //    }
-        //}
     }
 
-    int numCells = xadj.size()-1;
     int nProcs = nProcessors_;
 
     // output: cell -> processor addressing
@@ -327,12 +310,18 @@ Foam::metisDecomp::metisDecomp
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
+Foam::labelList Foam::metisDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& pointWeights
+)
 {
     if (points.size() != mesh_.nCells())
     {
-        FatalErrorIn("metisDecomp::decompose(const pointField&)")
-            << "Can use this decomposition method only for the whole mesh"
+        FatalErrorIn
+        (
+            "metisDecomp::decompose(const pointField&,const scalarField&)"
+        )   << "Can use this decomposition method only for the whole mesh"
             << endl
             << "and supply one coordinate (cellCentre) for every cell." << endl
             << "The number of coordinates " << points.size() << endl
@@ -340,19 +329,49 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
             << exit(FatalError);
     }
 
+    List<int> adjncy;
+    List<int> xadj;
+    calcMetisCSR
+    (
+        mesh_,
+        adjncy,
+        xadj
+    );
+
+    // Decompose using default weights
+    List<int> finalDecomp;
+    decompose(adjncy, xadj, pointWeights, finalDecomp);
+
+    // Copy back to labelList
+    labelList decomp(finalDecomp.size());
+    forAll(decomp, i)
+    {
+        decomp[i] = finalDecomp[i];
+    }
+    return decomp;
+}
+
+
+void Foam::metisDecomp::calcMetisCSR
+(
+    const polyMesh& mesh,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
     // Make Metis CSR (Compressed Storage Format) storage
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
 
-    List<int> xadj(mesh_.nCells()+1);
+    xadj.setSize(mesh.nCells()+1);
 
     // Initialise the number of internal faces of the cells to twice the
     // number of internal faces
-    label nInternalFaces = 2*mesh_.nInternalFaces();
+    label nInternalFaces = 2*mesh.nInternalFaces();
 
     // Check the boundary for coupled patches and add to the number of
     // internal faces
-    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
     forAll(pbm, patchi)
     {
@@ -364,17 +383,17 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
 
     // Create the adjncy array the size of the total number of internal and
     // coupled faces
-    List<int> adjncy(nInternalFaces);
+    adjncy.setSize(nInternalFaces);
 
     // Fill in xadj
     // ~~~~~~~~~~~~
     label freeAdj = 0;
 
-    for (label cellI = 0; cellI < mesh_.nCells(); cellI++)
+    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
     {
         xadj[cellI] = freeAdj;
 
-        const labelList& cFaces = mesh_.cells()[cellI];
+        const labelList& cFaces = mesh.cells()[cellI];
 
         forAll(cFaces, i)
         {
@@ -382,7 +401,7 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
 
             if
             (
-                mesh_.isInternalFace(faceI)
+                mesh.isInternalFace(faceI)
              || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
             )
             {
@@ -390,19 +409,19 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
             }
         }
     }
-    xadj[mesh_.nCells()] = freeAdj;
+    xadj[mesh.nCells()] = freeAdj;
 
 
     // Fill in adjncy
     // ~~~~~~~~~~~~~~
 
-    labelList nFacesPerCell(mesh_.nCells(), 0);
+    labelList nFacesPerCell(mesh.nCells(), 0);
 
     // Internal faces
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
     {
-        label own = mesh_.faceOwner()[faceI];
-        label nei = mesh_.faceNeighbour()[faceI];
+        label own = mesh.faceOwner()[faceI];
+        label nei = mesh.faceNeighbour()[faceI];
 
         adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
         adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
@@ -427,18 +446,6 @@ Foam::labelList Foam::metisDecomp::decompose(const pointField& points)
             }
         }
     }
-
-    // Decompose using default weights
-    List<int> finalDecomp;
-    decompose(adjncy, xadj, finalDecomp);
-
-    // Copy back to labelList
-    labelList decomp(finalDecomp.size());
-    forAll(decomp, i)
-    {
-        decomp[i] = finalDecomp[i];
-    }
-    return decomp;
 }
 
 
@@ -487,14 +494,16 @@ void Foam::metisDecomp::calcMetisCSR
 Foam::labelList Foam::metisDecomp::decompose
 (
     const labelList& agglom,
-    const pointField& agglomPoints
+    const pointField& agglomPoints,
+    const scalarField& agglomWeights
 )
 {
     if (agglom.size() != mesh_.nCells())
     {
         FatalErrorIn
         (
-            "parMetisDecomp::decompose(const labelList&, const pointField&)"
+            "metisDecomp::decompose"
+            "(const labelList&, const pointField&, const scalarField&)"
         )   << "Size of cell-to-coarse map " << agglom.size()
             << " differs from number of cells in mesh " << mesh_.nCells()
             << exit(FatalError);
@@ -521,7 +530,7 @@ Foam::labelList Foam::metisDecomp::decompose
 
     // Decompose using default weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, finalDecomp);
+    decompose(adjncy, xadj, agglomWeights, finalDecomp);
 
 
     // Rework back into decomposition for original mesh_
@@ -539,14 +548,16 @@ Foam::labelList Foam::metisDecomp::decompose
 Foam::labelList Foam::metisDecomp::decompose
 (
     const labelListList& globalCellCells,
-    const pointField& cellCentres
+    const pointField& cellCentres,
+    const scalarField& cellWeights
 )
 {
     if (cellCentres.size() != globalCellCells.size())
     {
         FatalErrorIn
         (
-            "metisDecomp::decompose(const pointField&, const labelListList&)"
+            "metisDecomp::decompose"
+            "(const pointField&, const labelListList&, const scalarField&)"
         )   << "Inconsistent number of cells (" << globalCellCells.size()
             << ") and number of cell centres (" << cellCentres.size()
             << ")." << exit(FatalError);
@@ -564,7 +575,7 @@ Foam::labelList Foam::metisDecomp::decompose
 
     // Decompose using default weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, finalDecomp);
+    decompose(adjncy, xadj, cellWeights, finalDecomp);
 
     // Copy back to labelList
     labelList decomp(finalDecomp.size());
diff --git a/src/decompositionMethods/decompositionMethods/metisDecomp/metisDecomp.H b/src/decompositionMethods/decompositionMethods/metisDecomp/metisDecomp.H
index a91427a8e9..6a38e25c8a 100644
--- a/src/decompositionMethods/decompositionMethods/metisDecomp/metisDecomp.H
+++ b/src/decompositionMethods/decompositionMethods/metisDecomp/metisDecomp.H
@@ -60,6 +60,7 @@ class metisDecomp
         (
             const List<int>& adjncy,
             const List<int>& xadj,
+            const scalarField& cellWeights,
             List<int>& finalDecomp
         );
 
@@ -100,7 +101,11 @@ public:
 
         //- Return for every coordinate the wanted processor number. Use the
         //  mesh connectivity (if needed)
-        virtual labelList decompose(const pointField&);
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        );
 
         //- Return for every coordinate the wanted processor number. Gets
         //  passed agglomeration map (from fine to coarse cells) and coarse cell
@@ -109,7 +114,8 @@ public:
         virtual labelList decompose
         (
             const labelList& agglom,
-            const pointField&
+            const pointField& regionPoints,
+            const scalarField& regionWeights
         );
 
         //- Return for every coordinate the wanted processor number. Explicitly
@@ -122,10 +128,21 @@ public:
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc
+            const pointField& cc,
+            const scalarField& cWeights
         );
 
-        //- Helper to convert cellcells into Metis storage
+        //- Helper to convert connectivity (supplied as owner,neighbour) into
+        //  Metis storage
+        static void calcMetisCSR
+        (
+            const polyMesh& mesh,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        //- Helper to convert connectivity (supplied as cellcells) into
+        //  Metis storage
         static void calcMetisCSR
         (
             const labelListList& globalCellCells,
diff --git a/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C b/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C
index 994a824192..61a9483837 100644
--- a/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.C
@@ -65,6 +65,7 @@ License
 #include "Time.H"
 #include "cyclicPolyPatch.H"
 #include "OFstream.H"
+#include "metisDecomp.H"
 
 extern "C"
 {
@@ -125,6 +126,7 @@ Foam::label Foam::scotchDecomp::decompose
 (
     const List<int>& adjncy,
     const List<int>& xadj,
+    const scalarField& cWeights,
 
     List<int>& finalDecomp
 )
@@ -160,6 +162,33 @@ Foam::label Foam::scotchDecomp::decompose
     // Graph
     // ~~~~~
 
+    List<int> velotab;
+
+
+    // Check for externally provided cellweights and if so initialise weights
+    scalar maxWeights = gMax(cWeights);
+    if (cWeights.size() > 0)
+    {
+        if (cWeights.size() != xadj.size()-1)
+        {
+            FatalErrorIn
+            (
+                "parMetisDecomp::decompose"
+                "(const pointField&, const scalarField&)"
+            )   << "Number of cell weights " << cWeights.size()
+                << " does not equal number of cells " << xadj.size()-1
+                << exit(FatalError);
+        }
+        // Convert to integers.
+        velotab.setSize(cWeights.size());
+        forAll(velotab, i)
+        {
+            velotab[i] = int(1000000*cWeights[i]/maxWeights);
+        }
+    }
+
+
+
     SCOTCH_Graph grafdat;
     check(SCOTCH_graphInit(&grafdat), "SCOTCH_graphInit");
     check
@@ -171,7 +200,7 @@ Foam::label Foam::scotchDecomp::decompose
             xadj.size()-1,          // vertnbr, nCells
             xadj.begin(),           // verttab, start index per cell into adjncy
             &xadj[1],               // vendtab, end index  ,,
-            NULL,                   // velotab, vertex weights
+            velotab.begin(),        // velotab, vertex weights
             NULL,                   // vlbltab
             adjncy.size(),          // edgenbr, number of arcs
             adjncy.begin(),         // edgetab
@@ -340,7 +369,11 @@ Foam::scotchDecomp::scotchDecomp
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::labelList Foam::scotchDecomp::decompose(const pointField& points)
+Foam::labelList Foam::scotchDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& pointWeights
+)
 {
     if (points.size() != mesh_.nCells())
     {
@@ -356,94 +389,18 @@ Foam::labelList Foam::scotchDecomp::decompose(const pointField& points)
     // Make Metis CSR (Compressed Storage Format) storage
     //   adjncy      : contains neighbours (= edges in graph)
     //   xadj(celli) : start of information in adjncy for celli
-
-    List<int> xadj(mesh_.nCells()+1);
-
-    // Initialise the number of internal faces of the cells to twice the
-    // number of internal faces
-    label nInternalFaces = 2*mesh_.nInternalFaces();
-
-    // Check the boundary for coupled patches and add to the number of
-    // internal faces
-    const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
-
-    forAll(pbm, patchi)
-    {
-        if (isA<cyclicPolyPatch>(pbm[patchi]))
-        {
-            nInternalFaces += pbm[patchi].size();
-        }
-    }
-
-    // Create the adjncy array the size of the total number of internal and
-    // coupled faces
-    List<int> adjncy(nInternalFaces);
-
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-    label freeAdj = 0;
-
-    for (label cellI = 0; cellI < mesh_.nCells(); cellI++)
-    {
-        xadj[cellI] = freeAdj;
-
-        const labelList& cFaces = mesh_.cells()[cellI];
-
-        forAll(cFaces, i)
-        {
-            label faceI = cFaces[i];
-
-            if
-            (
-                mesh_.isInternalFace(faceI)
-             || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
-            )
-            {
-                freeAdj++;
-            }
-        }
-    }
-    xadj[mesh_.nCells()] = freeAdj;
-
-
-    // Fill in adjncy
-    // ~~~~~~~~~~~~~~
-
-    labelList nFacesPerCell(mesh_.nCells(), 0);
-
-    // Internal faces
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-    {
-        label own = mesh_.faceOwner()[faceI];
-        label nei = mesh_.faceNeighbour()[faceI];
-
-        adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-        adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
-    }
-
-    // Coupled faces. Only cyclics done.
-    forAll(pbm, patchi)
-    {
-        if (isA<cyclicPolyPatch>(pbm[patchi]))
-        {
-            const unallocLabelList& faceCells = pbm[patchi].faceCells();
-
-            label sizeby2 = faceCells.size()/2;
-
-            for (label facei=0; facei<sizeby2; facei++)
-            {
-                label own = faceCells[facei];
-                label nei = faceCells[facei + sizeby2];
-
-                adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-                adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
-            }
-        }
-    }
+    List<int> adjncy;
+    List<int> xadj;
+    metisDecomp::calcMetisCSR
+    (
+        mesh_,
+        adjncy,
+        xadj
+    );
 
     // Decompose using default weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, finalDecomp);
+    decompose(adjncy, xadj, pointWeights, finalDecomp);
 
     // Copy back to labelList
     labelList decomp(finalDecomp.size());
@@ -455,52 +412,11 @@ Foam::labelList Foam::scotchDecomp::decompose(const pointField& points)
 }
 
 
-// From cell-cell connections to Metis format (like CompactListList)
-void Foam::scotchDecomp::calcMetisCSR
-(
-    const labelListList& cellCells,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    // Count number of internal faces
-    label nConnections = 0;
-
-    forAll(cellCells, coarseI)
-    {
-        nConnections += cellCells[coarseI].size();
-    }
-
-    // Create the adjncy array as twice the size of the total number of
-    // internal faces
-    adjncy.setSize(nConnections);
-
-    xadj.setSize(cellCells.size()+1);
-
-
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-    label freeAdj = 0;
-
-    forAll(cellCells, coarseI)
-    {
-        xadj[coarseI] = freeAdj;
-
-        const labelList& cCells = cellCells[coarseI];
-
-        forAll(cCells, i)
-        {
-            adjncy[freeAdj++] = cCells[i];
-        }
-    }
-    xadj[cellCells.size()] = freeAdj;
-}
-
-
 Foam::labelList Foam::scotchDecomp::decompose
 (
     const labelList& agglom,
-    const pointField& agglomPoints
+    const pointField& agglomPoints,
+    const scalarField& pointWeights
 )
 {
     if (agglom.size() != mesh_.nCells())
@@ -529,13 +445,12 @@ Foam::labelList Foam::scotchDecomp::decompose
             cellCells
         );
 
-        calcMetisCSR(cellCells, adjncy, xadj);
+        metisDecomp::calcMetisCSR(cellCells, adjncy, xadj);
     }
 
-    // Decompose using default weights
+    // Decompose using weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, finalDecomp);
-
+    decompose(adjncy, xadj, pointWeights, finalDecomp);
 
     // Rework back into decomposition for original mesh_
     labelList fineDistribution(agglom.size());
@@ -552,7 +467,8 @@ Foam::labelList Foam::scotchDecomp::decompose
 Foam::labelList Foam::scotchDecomp::decompose
 (
     const labelListList& globalCellCells,
-    const pointField& cellCentres
+    const pointField& cellCentres,
+    const scalarField& cWeights
 )
 {
     if (cellCentres.size() != globalCellCells.size())
@@ -572,12 +488,11 @@ Foam::labelList Foam::scotchDecomp::decompose
 
     List<int> adjncy;
     List<int> xadj;
-    calcMetisCSR(globalCellCells, adjncy, xadj);
+    metisDecomp::calcMetisCSR(globalCellCells, adjncy, xadj);
 
-
-    // Decompose using default weights
+    // Decompose using weights
     List<int> finalDecomp;
-    decompose(adjncy, xadj, finalDecomp);
+    decompose(adjncy, xadj, cWeights, finalDecomp);
 
     // Copy back to labelList
     labelList decomp(finalDecomp.size());
diff --git a/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.H b/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.H
index e470f84f98..54efb1bafa 100644
--- a/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.H
+++ b/src/decompositionMethods/decompositionMethods/scotchDecomp/scotchDecomp.H
@@ -63,6 +63,7 @@ class scotchDecomp
         (
             const List<int>& adjncy,
             const List<int>& xadj,
+            const scalarField& cWeights,
             List<int>& finalDecomp
         );
 
@@ -103,7 +104,11 @@ public:
 
         //- Return for every coordinate the wanted processor number. Use the
         //  mesh connectivity (if needed)
-        virtual labelList decompose(const pointField&);
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        );
 
         //- Return for every coordinate the wanted processor number. Gets
         //  passed agglomeration map (from fine to coarse cells) and coarse cell
@@ -112,7 +117,8 @@ public:
         virtual labelList decompose
         (
             const labelList& agglom,
-            const pointField&
+            const pointField& regionPoints,
+            const scalarField& regionWeights
         );
 
         //- Return for every coordinate the wanted processor number. Explicitly
@@ -125,15 +131,8 @@ public:
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc
-        );
-
-        //- Helper to convert cellcells into compact row storage
-        static void calcMetisCSR
-        (
-            const labelListList& globalCellCells,
-            List<int>& adjncy,
-            List<int>& xadj
+            const pointField& cc,
+            const scalarField& cWeights
         );
 
 };
diff --git a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
index 57cf17e55f..50d5842fb8 100644
--- a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
+++ b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.C
@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "simpleGeomDecomp.H"
@@ -93,6 +91,49 @@ void Foam::simpleGeomDecomp::assignToProcessorGroup
 }
 
 
+void Foam::simpleGeomDecomp::assignToProcessorGroup
+(
+    labelList& processorGroup,
+    const label nProcGroup,
+    const labelList& indices,
+    const scalarField& weights,
+    const scalar summedWeights
+)
+{
+    // This routine gets the sorted points.
+    // Easiest to explain with an example.
+    // E.g. 400 points, sum of weights : 513.
+    // Now with number of divisions in this direction (nProcGroup) : 4
+    // gives the split at 513/4 = 128
+    // So summed weight from 0..128 goes into bin 0,
+    //     ,,              128..256 goes into bin 1
+    //   etc.
+    // Finally any remaining ones go into the last bin (3).
+
+    const scalar jump = summedWeights/nProcGroup;
+    const label nProcGroupM1 = nProcGroup - 1;
+    scalar sumWeights = 0;
+    label ind = 0;
+    label j = 0;
+
+    // assign cells to all except last group.
+    for (j=0; j<nProcGroupM1; j++)
+    {
+        const scalar limit = jump*scalar(j + 1);
+        while (sumWeights < limit)
+        {
+            sumWeights += weights[indices[ind]];
+            processorGroup[ind++] = j;
+        }
+    }
+    // Ensure last included.
+    while (ind < processorGroup.size())
+    {
+       processorGroup[ind++] = nProcGroupM1;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::simpleGeomDecomp::simpleGeomDecomp(const dictionary& decompositionDict)
@@ -121,7 +162,10 @@ Foam::labelList Foam::simpleGeomDecomp::decompose(const pointField& points)
     labelList processorGroups(points.size());
 
     labelList pointIndices(points.size());
-    forAll(pointIndices, i) {pointIndices[i] = i;}
+    forAll(pointIndices, i)
+    {
+        pointIndices[i] = i;
+    }
 
     pointField rotatedPoints = rotDelta_ & points;
 
@@ -175,9 +219,101 @@ Foam::labelList Foam::simpleGeomDecomp::decompose(const pointField& points)
         finalDecomp[pointIndices[i]] += n_.x()*n_.y()*processorGroups[i];
     }
 
-
     return finalDecomp;
 }
 
 
+Foam::labelList Foam::simpleGeomDecomp::decompose
+(
+    const pointField& points,
+    const scalarField& weights
+)
+{
+    // construct a list for the final result
+    labelList finalDecomp(points.size());
+
+    labelList processorGroups(points.size());
+
+    labelList pointIndices(points.size());
+    forAll(pointIndices, i)
+    {
+        pointIndices[i] = i;
+    }
+
+    pointField rotatedPoints = rotDelta_ & points;
+
+    // and one to take the processor group id's. For each direction.
+    // we assign the processors to groups of processors labelled
+    // 0..nX to give a banded structure on the mesh. Then we
+    // construct the actual processor number by treating this as
+    // the units part of the processor number.
+    sort
+    (
+        pointIndices,
+        UList<scalar>::less(rotatedPoints.component(vector::X))
+    );
+    
+    const scalar summedWeights = sum(weights);
+    assignToProcessorGroup
+    (
+        processorGroups, 
+        n_.x(), 
+        pointIndices, 
+        weights, 
+        summedWeights
+    );
+
+    forAll(points, i)
+    {
+        finalDecomp[pointIndices[i]] = processorGroups[i];
+    }
+
+
+    // now do the same thing in the Y direction. These processor group
+    // numbers add multiples of nX to the proc. number (columns)
+    sort
+    (
+        pointIndices,
+        UList<scalar>::less(rotatedPoints.component(vector::Y))
+    );
+
+    assignToProcessorGroup
+    (
+        processorGroups,
+        n_.y(), 
+        pointIndices, 
+        weights, 
+        summedWeights
+    );
+
+    forAll(points, i)
+    {
+        finalDecomp[pointIndices[i]] += n_.x()*processorGroups[i];
+    }
+
+
+    // finally in the Z direction. Now we add multiples of nX*nY to give
+    // layers
+    sort
+    (
+        pointIndices,
+        UList<scalar>::less(rotatedPoints.component(vector::Z))
+    );
+
+    assignToProcessorGroup
+    (
+        processorGroups, 
+        n_.z(), 
+        pointIndices, 
+        weights, 
+        summedWeights
+    );
+
+    forAll(points, i)
+    {
+        finalDecomp[pointIndices[i]] += n_.x()*n_.y()*processorGroups[i];
+    }
+
+    return finalDecomp;
+} 
 // ************************************************************************* //
diff --git a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
index e1f312b5f3..595db50c2c 100644
--- a/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
+++ b/src/decompositionMethods/decompositionMethods/simpleGeomDecomp/simpleGeomDecomp.H
@@ -52,6 +52,15 @@ class simpleGeomDecomp
 
         void assignToProcessorGroup(labelList& processorGroup, const label);
 
+        void assignToProcessorGroup
+        (
+            labelList& processorGroup, 
+            const label nProcGroup,
+            const labelList& indices,
+            const scalarField& weights,
+            const scalar summedWeights
+        );
+
         //- Disallow default bitwise copy construct and assignment
         void operator=(const simpleGeomDecomp&);
         simpleGeomDecomp(const simpleGeomDecomp&);
@@ -86,21 +95,31 @@ public:
 
         virtual bool parallelAware() const
         {
-            // simple decomp does not attempt to do anything across proc
+            // simpleDecomp does not attempt to do anything across proc
             // boundaries
             return false;
         }
 
-        virtual labelList decompose(const pointField&);
+        virtual labelList decompose
+        (
+            const pointField& points
+        );
+
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        );
 
         //- Explicitly provided connectivity
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc
+            const pointField& cc,
+            const scalarField& cWeights
         )
         {
-            return decompose(cc);
+            return decompose(cc, cWeights);
         }
 };
 
diff --git a/src/decompositionMethods/parMetisDecomp/parMetisDecomp.C b/src/decompositionMethods/parMetisDecomp/parMetisDecomp.C
index c6c6e061e9..c4dde21d65 100644
--- a/src/decompositionMethods/parMetisDecomp/parMetisDecomp.C
+++ b/src/decompositionMethods/parMetisDecomp/parMetisDecomp.C
@@ -24,9 +24,9 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "syncTools.H"
 #include "parMetisDecomp.H"
 #include "metisDecomp.H"
+#include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 #include "floatScalar.H"
 #include "polyMesh.H"
@@ -370,15 +370,22 @@ Foam::parMetisDecomp::parMetisDecomp
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
+Foam::labelList Foam::parMetisDecomp::decompose
+(
+    const pointField& cc,
+    const scalarField& cWeights
+)
 {
-    if (points.size() != mesh_.nCells())
+    if (cc.size() != mesh_.nCells())
     {
-        FatalErrorIn("parMetisDecomp::decompose(const pointField&)")
-            << "Can use this decomposition method only for the whole mesh"
+        FatalErrorIn
+        (
+            "parMetisDecomp::decompose"
+            "(const pointField&, const scalarField&)"
+        )   << "Can use this decomposition method only for the whole mesh"
             << endl
             << "and supply one coordinate (cellCentre) for every cell." << endl
-            << "The number of coordinates " << points.size() << endl
+            << "The number of coordinates " << cc.size() << endl
             << "The number of cells in the mesh " << mesh_.nCells()
             << exit(FatalError);
     }
@@ -386,156 +393,24 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
     // For running sequential ...
     if (Pstream::nProcs() <= 1)
     {
-        return metisDecomp(decompositionDict_, mesh_).decompose(points);
-    }
-
-    // Create global cell numbers
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Get number of cells on all processors
-    List<int> nLocalCells(Pstream::nProcs());
-    nLocalCells[Pstream::myProcNo()] = mesh_.nCells();
-    Pstream::gatherList(nLocalCells);
-    Pstream::scatterList(nLocalCells);
-
-    // Get cell offsets.
-    List<int> cellOffsets(Pstream::nProcs()+1);
-    int nGlobalCells = 0;
-    forAll(nLocalCells, procI)
-    {
-        cellOffsets[procI] = nGlobalCells;
-        nGlobalCells += nLocalCells[procI];
-    }
-    cellOffsets[Pstream::nProcs()] = nGlobalCells;
-
-    int myOffset = cellOffsets[Pstream::myProcNo()];
-
-
-    //
-    // Make Metis Distributed CSR (Compressed Storage Format) storage
-    //   adjncy      : contains cellCells (= edges in graph)
-    //   xadj(celli) : start of information in adjncy for celli
-    //
-
-
-
-    const labelList& faceOwner = mesh_.faceOwner();
-    const labelList& faceNeighbour = mesh_.faceNeighbour();
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-
-    // Get renumbered owner on other side of coupled faces
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    List<int> globalNeighbour(mesh_.nFaces()-mesh_.nInternalFaces());
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-            label bFaceI = pp.start() - mesh_.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                globalNeighbour[bFaceI++] = faceOwner[faceI++] + myOffset;
-            }
-        }
-    }
-
-    // Get the cell on the other side of coupled patches
-    syncTools::swapBoundaryFaceList(mesh_, globalNeighbour, false);
-
-
-    // Count number of faces (internal + coupled)
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    // Number of faces per cell
-    List<int> nFacesPerCell(mesh_.nCells(), 0);
-
-    // Number of coupled faces
-    label nCoupledFaces = 0;
-
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-    {
-        nFacesPerCell[faceOwner[faceI]]++;
-        nFacesPerCell[faceNeighbour[faceI]]++;
-    }
-    // Handle coupled faces
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-
-            forAll(pp, i)
-            {
-                nCoupledFaces++;
-                nFacesPerCell[faceOwner[faceI++]]++;
-            }
-        }
+        return metisDecomp(decompositionDict_, mesh_).decompose
+        (
+            cc,
+            cWeights
+        );
     }
 
 
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-
-    Field<int> xadj(mesh_.nCells()+1, -1);
-
-    int freeAdj = 0;
-
-    for (label cellI = 0; cellI < mesh_.nCells(); cellI++)
-    {
-        xadj[cellI] = freeAdj;
-
-        freeAdj += nFacesPerCell[cellI];
-    }
-    xadj[mesh_.nCells()] = freeAdj;
-
-
-
-    // Fill in adjncy
-    // ~~~~~~~~~~~~~~
-
-    Field<int> adjncy(2*mesh_.nInternalFaces() + nCoupledFaces, -1);
-
-    nFacesPerCell = 0;
-
-    // For internal faces is just offsetted owner and neighbour
-    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-    {
-        label own = faceOwner[faceI];
-        label nei = faceNeighbour[faceI];
-
-        adjncy[xadj[own] + nFacesPerCell[own]++] = nei + myOffset;
-        adjncy[xadj[nei] + nFacesPerCell[nei]++] = own + myOffset;
-    }
-    // For boundary faces is offsetted coupled neighbour
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
-        {
-            label faceI = pp.start();
-            label bFaceI = pp.start()-mesh_.nInternalFaces();
-
-            forAll(pp, i)
-            {
-                label own = faceOwner[faceI];
-                adjncy[xadj[own] + nFacesPerCell[own]++] =
-                    globalNeighbour[bFaceI];
-
-                faceI++;
-                bFaceI++;
-            }
-        }
-    }
-
+    // Connections
+    Field<int> adjncy;
+    // Offsets into adjncy
+    Field<int> xadj;
+    calcMetisDistributedCSR
+    (
+        mesh_,
+        adjncy,
+        xadj
+    );
 
 
     // decomposition options. 0 = use defaults
@@ -550,6 +425,30 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
     // face weights (so on the edges of the dual)
     Field<int> faceWeights;
 
+
+    // Check for externally provided cellweights and if so initialise weights
+    scalar maxWeights = gMax(cWeights);
+    if (cWeights.size() > 0)
+    {
+        if (cWeights.size() != mesh_.nCells())
+        {
+            FatalErrorIn
+            (
+                "metisDecomp::decompose"
+                "(const pointField&, const scalarField&)"
+            )   << "Number of cell weights " << cWeights.size()
+                << " does not equal number of cells " << mesh_.nCells()
+                << exit(FatalError);
+        }
+        // Convert to integers.
+        cellWeights.setSize(cWeights.size());
+        forAll(cellWeights, i)
+        {
+            cellWeights[i] = int(1000000*cWeights[i]/maxWeights);
+        }
+    }
+
+
     // Check for user supplied weights and decomp options
     if (decompositionDict_.found("metisCoeffs"))
     {
@@ -577,8 +476,11 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
 
             if (cellWeights.size() != mesh_.nCells())
             {
-                FatalErrorIn("parMetisDecomp::decompose(const pointField&)")
-                    << "Number of cell weights " << cellWeights.size()
+                FatalErrorIn
+                (
+                    "parMetisDecomp::decompose"
+                    "(const pointField&, const scalarField&)"
+                )   << "Number of cell weights " << cellWeights.size()
                     << " read from " << cellIOWeights.objectPath()
                     << " does not equal number of cells " << mesh_.nCells()
                     << exit(FatalError);
@@ -604,30 +506,35 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
 
             if (weights.size() != mesh_.nFaces())
             {
-                FatalErrorIn("parMetisDecomp::decompose(const pointField&)")
-                    << "Number of face weights " << weights.size()
+                FatalErrorIn
+                (
+                    "parMetisDecomp::decompose"
+                    "(const pointField&, const scalarField&)"
+                )   << "Number of face weights " << weights.size()
                     << " does not equal number of internal and boundary faces "
                     << mesh_.nFaces()
                     << exit(FatalError);
             }
 
-            faceWeights.setSize(2*mesh_.nInternalFaces()+nCoupledFaces);
+            faceWeights.setSize(adjncy.size());
 
             // Assume symmetric weights. Keep same ordering as adjncy.
-            nFacesPerCell = 0;
+            List<int> nFacesPerCell(mesh_.nCells(), 0);
 
             // Handle internal faces
             for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
             {
                 label w = weights[faceI];
 
-                label own = faceOwner[faceI];
-                label nei = faceNeighbour[faceI];
+                label own = mesh_.faceOwner()[faceI];
+                label nei = mesh_.faceNeighbour()[faceI];
 
                 faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
                 faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
             }
             // Coupled boundary faces
+            const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
             forAll(patches, patchI)
             {
                 const polyPatch& pp = patches[patchI];
@@ -639,8 +546,8 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
                     forAll(pp, i)
                     {
                         label w = weights[faceI];
-                        label own = faceOwner[faceI];
-                        adjncy[xadj[own] + nFacesPerCell[own]++] = w;
+                        label own = mesh_.faceOwner()[faceI];
+                        faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
                         faceI++;
                     }
                 }
@@ -654,8 +561,11 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
 
             if (options.size() != 3)
             {
-                FatalErrorIn("parMetisDecomp::decompose(const pointField&)")
-                    << "Number of options " << options.size()
+                FatalErrorIn
+                (
+                    "parMetisDecomp::decompose"
+                    "(const pointField&, const scalarField&)"
+                )   << "Number of options " << options.size()
                     << " should be three." << exit(FatalError);
             }
         }
@@ -668,7 +578,7 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
     (
         xadj,
         adjncy,
-        points,
+        cc,
         cellWeights,
         faceWeights,
         options,
@@ -689,7 +599,8 @@ Foam::labelList Foam::parMetisDecomp::decompose(const pointField& points)
 Foam::labelList Foam::parMetisDecomp::decompose
 (
     const labelList& cellToRegion,
-    const pointField& regionPoints
+    const pointField& regionPoints,
+    const scalarField& regionWeights
 )
 {
     const labelList& faceOwner = mesh_.faceOwner();
@@ -798,7 +709,15 @@ Foam::labelList Foam::parMetisDecomp::decompose
         }
     }
 
-    labelList regionDecomp(decompose(globalRegionRegions, regionPoints));
+    labelList regionDecomp
+    (
+        decompose
+        (
+            globalRegionRegions,
+            regionPoints,
+            regionWeights
+        )
+    );
 
     // Rework back into decomposition for original mesh_
     labelList cellDistribution(cellToRegion.size());
@@ -814,24 +733,26 @@ Foam::labelList Foam::parMetisDecomp::decompose
 Foam::labelList Foam::parMetisDecomp::decompose
 (
     const labelListList& globalCellCells,
-    const pointField& cellCentres
+    const pointField& cellCentres,
+    const scalarField& cWeights
 )
 {
     if (cellCentres.size() != globalCellCells.size())
     {
         FatalErrorIn
         (
-            "parMetisDecomp::decompose(const labelListList&, const pointField&)"
+            "parMetisDecomp::decompose(const labelListList&"
+            ", const pointField&, const scalarField&)"
         )   << "Inconsistent number of cells (" << globalCellCells.size()
             << ") and number of cell centres (" << cellCentres.size()
-            << ")." << exit(FatalError);
+            << ") or weights (" << cWeights.size() << ")." << exit(FatalError);
     }
 
     // For running sequential ...
     if (Pstream::nProcs() <= 1)
     {
         return metisDecomp(decompositionDict_, mesh_)
-            .decompose(globalCellCells, cellCentres);
+            .decompose(globalCellCells, cellCentres, cWeights);
     }
 
 
@@ -855,106 +776,35 @@ Foam::labelList Foam::parMetisDecomp::decompose
     // face weights (so on the edges of the dual)
     Field<int> faceWeights;
 
+
+    // Check for externally provided cellweights and if so initialise weights
+    scalar maxWeights = gMax(cWeights);
+    if (cWeights.size() > 0)
+    {
+        if (cWeights.size() != globalCellCells.size())
+        {
+            FatalErrorIn
+            (
+                "parMetisDecomp::decompose(const labelListList&"
+                ", const pointField&, const scalarField&)"
+            )   << "Number of cell weights " << cWeights.size()
+                << " does not equal number of cells " << globalCellCells.size()
+                << exit(FatalError);
+        }
+        // Convert to integers.
+        cellWeights.setSize(cWeights.size());
+        forAll(cellWeights, i)
+        {
+            cellWeights[i] = int(1000000*cWeights[i]/maxWeights);
+        }
+    }
+
+
     // Check for user supplied weights and decomp options
     if (decompositionDict_.found("metisCoeffs"))
     {
         const dictionary& metisCoeffs =
             decompositionDict_.subDict("metisCoeffs");
-        word weightsFile;
-
-        if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
-        {
-            Info<< "parMetisDecomp : Using cell-based weights read from "
-                << weightsFile << endl;
-
-            labelIOField cellIOWeights
-            (
-                IOobject
-                (
-                    weightsFile,
-                    mesh_.time().timeName(),
-                    mesh_,
-                    IOobject::MUST_READ,
-                    IOobject::AUTO_WRITE
-                )
-            );
-            cellWeights.transfer(cellIOWeights);
-
-            if (cellWeights.size() != cellCentres.size())
-            {
-                FatalErrorIn
-                (
-                    "parMetisDecomp::decompose"
-                    "(const labelListList&, const pointField&)"
-                )   << "Number of cell weights " << cellWeights.size()
-                    << " read from " << cellIOWeights.objectPath()
-                    << " does not equal number of cells " << cellCentres.size()
-                    << exit(FatalError);
-            }
-        }
-
-        //- faceWeights disabled. Only makes sense for cellCells from mesh.
-        //if (metisCoeffs.readIfPresent("faceWeightsFile", weightsFile))
-        //{
-        //    Info<< "parMetisDecomp : Using face-based weights read from "
-        //        << weightsFile << endl;
-        //
-        //    labelIOField weights
-        //    (
-        //        IOobject
-        //        (
-        //            weightsFile,
-        //            mesh_.time().timeName(),
-        //            mesh_,
-        //            IOobject::MUST_READ,
-        //            IOobject::AUTO_WRITE
-        //        )
-        //    );
-        //
-        //    if (weights.size() != mesh_.nFaces())
-        //    {
-        //        FatalErrorIn("parMetisDecomp::decompose(const pointField&)")
-        //            << "Number of face weights " << weights.size()
-        //            << " does not equal number of internal and boundary faces "
-        //            << mesh_.nFaces()
-        //            << exit(FatalError);
-        //    }
-        //
-        //    faceWeights.setSize(2*mesh_.nInternalFaces()+nCoupledFaces);
-        //
-        //    // Assume symmetric weights. Keep same ordering as adjncy.
-        //    nFacesPerCell = 0;
-        //
-        //    // Handle internal faces
-        //    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-        //    {
-        //        label w = weights[faceI];
-        //
-        //        label own = faceOwner[faceI];
-        //        label nei = faceNeighbour[faceI];
-        //
-        //        faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
-        //        faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
-        //    }
-        //    // Coupled boundary faces
-        //    forAll(patches, patchI)
-        //    {
-        //        const polyPatch& pp = patches[patchI];
-        //
-        //        if (pp.coupled())
-        //        {
-        //            label faceI = pp.start();
-        //
-        //            forAll(pp, i)
-        //            {
-        //                label w = weights[faceI];
-        //                label own = faceOwner[faceI];
-        //                adjncy[xadj[own] + nFacesPerCell[own]++] = w;
-        //                faceI++;
-        //            }
-        //        }
-        //    }
-        //}
 
         if (metisCoeffs.readIfPresent("options", options))
         {
@@ -965,8 +815,8 @@ Foam::labelList Foam::parMetisDecomp::decompose
             {
                 FatalErrorIn
                 (
-                    "parMetisDecomp::decompose"
-                    "(const labelListList&, const pointField&)"
+                    "parMetisDecomp::decompose(const labelListList&"
+                    ", const pointField&, const scalarField&)"
                 )   << "Number of options " << options.size()
                     << " should be three." << exit(FatalError);
             }
@@ -998,4 +848,146 @@ Foam::labelList Foam::parMetisDecomp::decompose
 }
 
 
+void Foam::parMetisDecomp::calcMetisDistributedCSR
+(
+    const polyMesh& mesh,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
+    // Create global cell numbers
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    globalIndex globalCells(mesh.nCells());
+
+
+    //
+    // Make Metis Distributed CSR (Compressed Storage Format) storage
+    //   adjncy      : contains cellCells (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+    //
+
+
+    const labelList& faceOwner = mesh.faceOwner();
+    const labelList& faceNeighbour = mesh.faceNeighbour();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+
+    // Get renumbered owner on other side of coupled faces
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    List<int> globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            label bFaceI = pp.start() - mesh.nInternalFaces();
+
+            forAll(pp, i)
+            {
+                globalNeighbour[bFaceI++] = globalCells.toGlobal
+                (
+                    faceOwner[faceI++]
+                );
+            }
+        }
+    }
+
+    // Get the cell on the other side of coupled patches
+    syncTools::swapBoundaryFaceList(mesh, globalNeighbour, false);
+
+
+    // Count number of faces (internal + coupled)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Number of faces per cell
+    List<int> nFacesPerCell(mesh.nCells(), 0);
+
+    // Number of coupled faces
+    label nCoupledFaces = 0;
+
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        nFacesPerCell[faceOwner[faceI]]++;
+        nFacesPerCell[faceNeighbour[faceI]]++;
+    }
+    // Handle coupled faces
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+
+            forAll(pp, i)
+            {
+                nCoupledFaces++;
+                nFacesPerCell[faceOwner[faceI++]]++;
+            }
+        }
+    }
+
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+
+    xadj.setSize(mesh.nCells()+1);
+
+    int freeAdj = 0;
+
+    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
+    {
+        xadj[cellI] = freeAdj;
+
+        freeAdj += nFacesPerCell[cellI];
+    }
+    xadj[mesh.nCells()] = freeAdj;
+
+
+
+    // Fill in adjncy
+    // ~~~~~~~~~~~~~~
+
+    adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces);
+
+    nFacesPerCell = 0;
+
+    // For internal faces is just offsetted owner and neighbour
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        label own = faceOwner[faceI];
+        label nei = faceNeighbour[faceI];
+
+        adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
+        adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
+    }
+    // For boundary faces is offsetted coupled neighbour
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            label faceI = pp.start();
+            label bFaceI = pp.start()-mesh.nInternalFaces();
+
+            forAll(pp, i)
+            {
+                label own = faceOwner[faceI];
+                adjncy[xadj[own] + nFacesPerCell[own]++] =
+                    globalNeighbour[bFaceI];
+
+                faceI++;
+                bFaceI++;
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/decompositionMethods/parMetisDecomp/parMetisDecomp.H b/src/decompositionMethods/parMetisDecomp/parMetisDecomp.H
index 13aa57277c..dda8173557 100644
--- a/src/decompositionMethods/parMetisDecomp/parMetisDecomp.H
+++ b/src/decompositionMethods/parMetisDecomp/parMetisDecomp.H
@@ -112,7 +112,11 @@ public:
 
         //- Return for every coordinate the wanted processor number. Use the
         //  mesh connectivity (if needed)
-        virtual labelList decompose(const pointField&);
+        virtual labelList decompose
+        (
+            const pointField& points,
+            const scalarField& pointWeights
+        );
 
         //- Return for every coordinate the wanted processor number. Gets
         //  passed agglomeration map (from fine to coarse cells) and coarse cell
@@ -120,8 +124,9 @@ public:
         //  functionality natively.
         virtual labelList decompose
         (
-            const labelList& agglom,
-            const pointField&
+            const labelList& cellToRegion,
+            const pointField& regionPoints,
+            const scalarField& regionWeights
         );
 
         //- Return for every coordinate the wanted processor number. Explicitly
@@ -134,7 +139,16 @@ public:
         virtual labelList decompose
         (
             const labelListList& globalCellCells,
-            const pointField& cc
+            const pointField& cc,
+            const scalarField& cWeights
+        );
+
+        //- Helper to convert mesh connectivity into distributed CSR
+        static void calcMetisDistributedCSR
+        (
+            const polyMesh&,
+            List<int>& adjncy,
+            List<int>& xadj
         );
 };
 
diff --git a/src/postProcessing/functionObjects/Allwmake b/src/postProcessing/functionObjects/Allwmake
index 6c642dfe7b..b068fe51cb 100755
--- a/src/postProcessing/functionObjects/Allwmake
+++ b/src/postProcessing/functionObjects/Allwmake
@@ -7,5 +7,6 @@ wmake libso forces
 wmake libso IO
 wmake libso utilities
 wmake libso systemCall
+wmake libso zones
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C
index ec77782df4..31e511377e 100644
--- a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C
+++ b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMax.C
@@ -197,13 +197,13 @@ void Foam::fieldMinMax::calcMinMaxFields<Foam::scalar>
 {
     if (obr_.foundObject<volScalarField>(fieldName))
     {
+        const volScalarField& field =
+            obr_.lookupObject<volScalarField>(fieldName);
+        scalar minValue = min(field).value();
+        scalar maxValue = max(field).value();
+
         if (Pstream::master())
         {
-            const volScalarField& field =
-                obr_.lookupObject<volScalarField>(fieldName);
-            scalar minValue = min(field).value();
-            scalar maxValue = max(field).value();
-
             fieldMinMaxFilePtr_() << obr_.time().value() << tab
                 << fieldName << tab << minValue << tab << maxValue << endl;
 
diff --git a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
index 25daf62dd6..66f212f4c3 100644
--- a/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldMinMax/fieldMinMaxTemplates.C
@@ -38,16 +38,16 @@ void Foam::fieldMinMax::calcMinMaxFields(const word& fieldName)
 
     if (obr_.foundObject<fieldType>(fieldName))
     {
-        if (Pstream::master())
+        const fieldType& field = obr_.lookupObject<fieldType>(fieldName);
+        switch (mode_)
         {
-            const fieldType& field = obr_.lookupObject<fieldType>(fieldName);
-            switch (mode_)
+            case mdMag:
             {
-                case mdMag:
-                {
-                    scalar minValue = min(mag(field)).value();
-                    scalar maxValue = max(mag(field)).value();
+                scalar minValue = min(mag(field)).value();
+                scalar maxValue = max(mag(field)).value();
 
+                if (Pstream::master())
+                {
                     fieldMinMaxFilePtr_() << obr_.time().value() << tab
                         << fieldName << tab << minValue << tab << maxValue
                         << endl;
@@ -61,13 +61,16 @@ void Foam::fieldMinMax::calcMinMaxFields(const word& fieldName)
                             << maxValue << nl
                             << endl;
                     }
-                    break;
                 }
-                case mdCmpt:
-                {
-                    Type minValue = min(field).value();
-                    Type maxValue = max(field).value();
+                break;
+            }
+            case mdCmpt:
+            {
+                Type minValue = min(field).value();
+                Type maxValue = max(field).value();
 
+                if (Pstream::master())
+                {
                     fieldMinMaxFilePtr_() << obr_.time().value() << tab
                         << fieldName << tab << minValue << tab << maxValue
                         << endl;
@@ -81,17 +84,17 @@ void Foam::fieldMinMax::calcMinMaxFields(const word& fieldName)
                             << maxValue << nl
                             << endl;
                     }
-                    break;
-                }
-                default:
-                {
-                    FatalErrorIn
-                    (
-                        "Foam::fieldMinMax::calcMinMaxFields"
-                        "(const word& fieldName)"
-                    )<< "Unknown min/max mode: " << modeTypeNames_[mode_]
-                     << exit(FatalError);
                 }
+                break;
+            }
+            default:
+            {
+                FatalErrorIn
+                (
+                    "Foam::fieldMinMax::calcMinMaxFields"
+                    "(const word& fieldName)"
+                )<< "Unknown min/max mode: " << modeTypeNames_[mode_]
+                 << exit(FatalError);
             }
         }
     }
diff --git a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.C b/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.C
index 0b8314667c..e048920899 100644
--- a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.C
+++ b/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.C
@@ -191,7 +191,7 @@ void Foam::faceZonesIntegration::write()
     if (active_)
     {
         makeFile();
-        scalar dm = 0.0;
+
         forAll(fItems_, fieldI)
         {
             const word& fieldName = fItems_[fieldI];
@@ -201,8 +201,30 @@ void Foam::faceZonesIntegration::write()
 
             const fvMesh& mesh = fD.mesh();
 
+            // 1. integrate over all face zones
+
+            scalarField integralVals(faceZonesSet_.size());
+
+            forAll(faceZonesSet_, setI)
+            {
+                const word name = faceZonesSet_[setI];
+
+                label zoneID = mesh.faceZones().findZoneID(name);
+
+                const faceZone& fz = mesh.faceZones()[zoneID];
+
+                integralVals[setI] = returnReduce
+                (
+                    calcFaceZonesIntegral(fD, fz),
+                    sumOp<scalar>()
+                );
+            }
+
+
             unsigned int w = IOstream::defaultPrecision() + 7;
 
+            // 2. Write only on master
+
             if
             (
                 Pstream::master()
@@ -213,26 +235,14 @@ void Foam::faceZonesIntegration::write()
 
                 os  << obr_.time().value();
 
-                const faceZoneMesh& faceZoneList =  mesh.faceZones();
-
-                forAll(faceZonesSet_, zoneI)
+                forAll(integralVals, setI)
                 {
-                    const word name = faceZonesSet_[zoneI];
-
-                    label zoneID = faceZoneList.findZoneID(name);
-
-                    const faceZone& fz = mesh.faceZones()[zoneID];
-
-                    dm = calcFaceZonesIntegral(fD, fz);
-
-                    reduce(dm, sumOp<scalar>());
-
-                    os  << ' ' << setw(w) << dm;
+                    os  << ' ' << setw(w) << integralVals[setI];
 
                     if (log_)
                     {
                         Info<< "faceZonesIntegration output:" << nl
-                            << "    Integration" << dm << endl;
+                            << "    Integration" << integralVals[setI] << endl;
                     }
                 }
 
@@ -258,7 +268,7 @@ Foam::scalar Foam::faceZonesIntegration::calcFaceZonesIntegral
 
         if (mesh.isInternalFace(faceI))
         {
-            if (fz.flipMap()[faceI])
+            if (fz.flipMap()[i])
             {
                 dm -= fD[faceI];
             }
@@ -275,7 +285,7 @@ Foam::scalar Foam::faceZonesIntegration::calcFaceZonesIntegral
             {
                 if (refCast<const processorPolyPatch>(pp).owner())
                 {
-                    if (fz.flipMap()[faceI])
+                    if (fz.flipMap()[i])
                     {
                         dm -= fD.boundaryField()[patchI][pp.whichFace(faceI)];
                     }
@@ -290,7 +300,7 @@ Foam::scalar Foam::faceZonesIntegration::calcFaceZonesIntegral
                 label patchFaceI = faceI - pp.start();
                 if (patchFaceI < pp.size()/2)
                 {
-                    if (fz.flipMap()[patchFaceI])
+                    if (fz.flipMap()[i])
                     {
                         dm -= fD.boundaryField()[patchI][patchFaceI];
                     }
@@ -303,7 +313,7 @@ Foam::scalar Foam::faceZonesIntegration::calcFaceZonesIntegral
             else if (!isA<emptyPolyPatch>(pp))
             {
                 label patchFaceI = faceI - pp.start();
-                if (fz.flipMap()[patchFaceI])
+                if (fz.flipMap()[i])
                 {
                     dm -= fD.boundaryField()[patchI][patchFaceI];
                 }
diff --git a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.H b/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.H
index 5fc957bba7..87659ac340 100644
--- a/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.H
+++ b/src/postProcessing/functionObjects/zones/faceZoneIntegration/faceZonesIntegration.H
@@ -73,17 +73,14 @@ protected:
 
         const objectRegistry& obr_;
 
-        //- on/off switch
-        bool active_;
-
-        //- Switch to send output to Info as well as to file
-        Switch log_;
-
-        //- Current open files
-        HashPtrTable<OFstream> faceZonesIntegrationFilePtr_;
-
         // Read from dictionary
 
+            //- on/off switch
+            bool active_;
+
+            //- Switch to send output to Info as well as to file
+            Switch log_;
+
             //- faceZones to integrate over
             wordList faceZonesSet_;
 
@@ -91,6 +88,11 @@ protected:
             wordList fItems_;
 
 
+        //- Current open files
+        HashPtrTable<OFstream> faceZonesIntegrationFilePtr_;
+
+
+
     // Private Member Functions
 
         //- If the integration file has not been created create it