diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
index 1e412a7cd..c215db2d3 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H
@@ -576,6 +576,15 @@ typename cut::opAddResult<AboveOp, BelowOp>::type triCut
     const BelowOp& belowOp
 );
 
+template<class Point, class AboveOp, class BelowOp>
+Tuple2<typename cut::opAddResult<AboveOp, BelowOp>::type, Point> triCutAndCentre
+(
+    const FixedList<vector, 3>& tri,
+    const FixedList<scalar, 3>& level,
+    const AboveOp& aboveOp,
+    const BelowOp& belowOp
+);
+
 //- As above, but with a plane specifying the location of the cut
 template<class AboveOp, class BelowOp>
 typename cut::opAddResult<AboveOp, BelowOp>::type triCut
@@ -596,6 +605,15 @@ typename cut::opAddResult<AboveOp, BelowOp>::type tetCut
     const BelowOp& belowOp
 );
 
+template<class Point, class AboveOp, class BelowOp>
+Tuple2<typename cut::opAddResult<AboveOp, BelowOp>::type, Point> tetCutAndCentre
+(
+    const FixedList<vector, 4>& tet,
+    const FixedList<scalar, 4>& level,
+    const AboveOp& aboveOp,
+    const BelowOp& belowOp
+);
+
 //- As above, but with a plane specifying the location of the cut
 template<class AboveOp, class BelowOp>
 typename cut::opAddResult<AboveOp, BelowOp>::type tetCut
diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H b/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H
index 9b7ca4794..263b2a352 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cutI.H
@@ -286,6 +286,27 @@ inline typename Op::result triCutTri
     return Op(triCutTri(op, f))(triCutTri(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> triCutTriAndCentre
+(
+    const Op& op,
+    const FixedList<Point, 3>& p,
+    const Pair<scalar>& f
+)
+{
+    FixedList<Point, 3> cutFace = triCutTri(p, f);
+    Point centre = Zero;
+    forAll(cutFace, pi)
+    {
+        centre += cutFace[pi];
+    }
+    centre /= 3;
+    return Tuple2<typename Op::result, Point>
+        (
+            Op(triCutTri(op, f))(cutFace),
+            centre
+        );
+}
 
 //- Apply an operation to a quad. Splits the quad into two tris.
 template<class Op, class OpData, class Point>
@@ -317,6 +338,27 @@ inline typename Op::result triCutQuad
     return quadOp<Op>(triCutQuad(op, f), triCutQuad(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> triCutQuadAndCentre
+(
+    const Op& op,
+    const FixedList<Point, 3>& p,
+    const FixedList<scalar, 2>& f
+)
+{
+    FixedList<Point, 4> cutFace = triCutQuad(p, f);
+    Point centre = Zero;
+    forAll(cutFace, pi)
+    {
+        centre += cutFace[pi];
+    }
+    centre /= 4;
+    return Tuple2<typename Op::result, Point>
+        (
+            quadOp<Op>(triCutQuad(op, f), cutFace),
+            centre
+        );
+}
 
 //- Cut a tet from a tet and apply an operation to the result. The cut is made
 //  along the three edges connected to vertex 0, and the cut locations are given
@@ -332,6 +374,28 @@ inline typename Op::result tetCutTet
     return Op(tetCutTet(op, f))(tetCutTet(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> tetCutTetAndCentre
+(
+    const Op& op,
+    const FixedList<Point, 4>& p,
+    const FixedList<scalar, 3>& f
+)
+{
+    FixedList<Point, 4> cutCell = tetCutTet(p, f);
+    Point centre = Zero;
+    forAll(cutCell, pi)
+    {
+        centre += cutCell[pi];
+    }
+    centre /= 4;
+
+    return Tuple2<typename Op::result, Point>
+        (
+            Op(tetCutTet(op, f))(cutCell),
+            centre
+        );
+}
 
 //- Apply an operation to a prism. Splits the prism into three tets.
 template<class Op, class OpData, class Point>
@@ -364,6 +428,27 @@ inline typename Op::result tetCutPrism0
     return prismOp<Op>(tetCutPrism0(op, f), tetCutPrism0(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> tetCutPrism0AndCentre
+(
+    const Op& op,
+    const FixedList<Point, 4>& p,
+    const FixedList<scalar, 3>& f
+)
+{
+    FixedList<Point, 6> cutCell = tetCutPrism0(p, f);
+    Point centre = Zero;
+    forAll(cutCell, pi)
+    {
+        centre += cutCell[pi];
+    }
+    centre /= 6;
+    return Tuple2<typename Op::result, Point>
+        (
+            prismOp<Op>(tetCutPrism0(op, f), cutCell),
+            centre
+        );
+}
 
 //- Cut a prism from a tet and apply an operation to the result. The cut is made
 //  along four edges, not edges 01 or 23, and the cut locations are given as
@@ -379,6 +464,27 @@ inline typename Op::result tetCutPrism01
     return prismOp<Op>(tetCutPrism01(op, f), tetCutPrism01(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> tetCutPrism01AndCentre
+(
+    const Op& op,
+    const FixedList<Point, 4>& p,
+    const FixedList<scalar, 4>& f
+)
+{
+    FixedList<Point, 6> cutCell = tetCutPrism01(p, f);
+    Point centre = Zero;
+    forAll(cutCell, pi)
+    {
+        centre += cutCell[pi];
+    }
+    centre /= 6;
+    return Tuple2<typename Op::result, Point>
+        (
+            prismOp<Op>(tetCutPrism01(op, f), cutCell),
+            centre
+        );
+}
 
 //- Cut a prism from a tet and apply an operation to the result. The cuts are
 //  the same as for tetCutPrism01. The result is the side connected to edge 23.
@@ -393,6 +499,28 @@ inline typename Op::result tetCutPrism23
     return prismOp<Op>(tetCutPrism23(op, f), tetCutPrism23(p, f));
 }
 
+template<class Op, class Point>
+Foam::Tuple2<typename Op::result, Point> tetCutPrism23AndCentre
+(
+    const Op& op,
+    const FixedList<Point, 4>& p,
+    const FixedList<scalar, 4>& f
+)
+{
+    FixedList<Point, 6> cutCell = tetCutPrism23(p, f);
+    Point centre = Zero;
+    forAll(cutCell, pi)
+    {
+        centre += cutCell[pi];
+    }
+    centre /= 6;
+    return Tuple2<typename Op::result, Point>
+        (
+            prismOp<Op>(tetCutPrism23(op, f), cutCell),
+            centre
+        );
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C b/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C
index 11bff4eba..25d1a1106 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C
+++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C
@@ -97,6 +97,96 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
     }
 }
 
+template<class Point, class AboveOp, class BelowOp>
+Foam::Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+Foam::triCutAndCentre
+(
+    const FixedList<Point, 3>& tri,
+    const FixedList<scalar, 3>& level,
+    const AboveOp& aboveOp,
+    const BelowOp& belowOp
+)
+{
+    // Quick return if all levels are zero or have the same sign
+    if (level[0] == 0 && level[1] == 0 && level[2] == 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp() + belowOp(), Zero);
+    }
+    if (level[0] >= 0 && level[1] >= 0 && level[2] >= 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp(tri) + belowOp(), Zero);
+    }
+    if (level[0] <= 0 && level[1] <= 0 && level[2] <= 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp() + belowOp(tri), triPointRef(tri[0], tri[1], tri[2]).centre());
+    }
+
+    // There will be just one edge without a sign change. Find it, and put it
+    // opposite the first vertex. This may change the sign of the tri.
+    FixedList<label, 3> indices({0, 1, 2});
+    label i;
+    for (i = 0; i < 3; ++ i)
+    {
+        if (level[(i + 1)%3]*level[(i + 2)%3] >= 0)
+        {
+            Swap(indices[0], indices[i]);
+            break;
+        }
+    }
+    if (i == 3)
+    {
+        FatalErrorInFunction
+            << "The number of tri vertices above the level set should always "
+            << "be 1" << exit(FatalError);
+    }
+
+    // Correct the sign
+    if (indices[0] != 0)
+    {
+        Swap(indices[1], indices[2]);
+    }
+
+    // Permute the data
+    const FixedList<Point, 3> p = triReorder(tri, indices);
+    const FixedList<scalar, 3> l = triReorder(level, indices);
+    AboveOp a = triReorder(aboveOp, indices);
+    BelowOp b = triReorder(belowOp, indices);
+
+    // Slice off one corner to form a tri and a quad
+    Pair<scalar> f;
+    for (label i = 0; i < 2; ++ i)
+    {
+        f[i] = l[0]/(l[0] - l[i+1]);
+    }
+    if (l[0] > 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+        (
+            triCutTriAndCentre(a, p, f).first() + triCutQuadAndCentre(b, p, f).first(),
+            (
+                triCutTriAndCentre(a, p, f).first()*triCutTriAndCentre(a, p, f).second()
+                + triCutQuadAndCentre(b, p, f).first()*triCutQuadAndCentre(b, p, f).second()
+            )
+          / (triCutTriAndCentre(a, p, f).first() + triCutQuadAndCentre(b, p, f).first())
+        );
+    }
+    else
+    {
+
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+        (
+            triCutQuadAndCentre(a, p, f).first() + triCutTriAndCentre(b, p, f).first(),
+            (
+                triCutQuadAndCentre(a, p, f).first()*triCutQuadAndCentre(a, p, f).second()
+                + triCutTriAndCentre(b, p, f).first()*triCutTriAndCentre(b, p, f).second()
+            )
+          / (triCutQuadAndCentre(a, p, f).first() + triCutTriAndCentre(b, p, f).first())
+        );
+    }
+}
 
 template<class AboveOp, class BelowOp>
 typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::triCut
@@ -246,6 +336,169 @@ typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut
     return aboveOp() + belowOp();
 }
 
+template<class Point, class AboveOp, class BelowOp>
+Foam::Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+Foam::tetCutAndCentre
+(
+    const FixedList<Point, 4>& tet,
+    const FixedList<scalar, 4>& level,
+    const AboveOp& aboveOp,
+    const BelowOp& belowOp
+)
+{
+    // Get the min and max over all four vertices and quick return if
+    // all levels are zero or have the same sign
+    scalar levelMin = vGreat, levelMax = - vGreat;
+    for (label i = 0; i < 4; ++ i)
+    {
+        levelMin = min(levelMin, level[i]);
+        levelMax = max(levelMax, level[i]);
+    }
+    if (levelMin == 0 && levelMax == 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp() + belowOp(), Zero);
+    }
+    if (levelMin >= 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp(tet) + belowOp(), Zero);
+    }
+    if (levelMax <= 0)
+    {
+        return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (aboveOp() + belowOp(tet), tetPointRef(tet[0], tet[1], tet[2], tet[3]).centre());
+    }
+
+    // Partition the level so that positive values are at the start. This is
+    // like a single iteration of quick-sort, except that the pivot is a hard-
+    // coded zero, rather than an element of the array. This can change the sign
+    // of the tet.
+    FixedList<label, 4> indices({0, 1, 2, 3});
+    bool signChange = false;
+    label i = 0, j = 3;
+    while (true)
+    {
+        while (i < j && level[indices[i]] > 0)
+        {
+            i ++;
+        }
+        while (j > i && level[indices[j]] <= 0)
+        {
+            j --;
+        }
+        if (i == j)
+        {
+            break;
+        }
+        Swap(indices[i], indices[j]);
+        signChange = !signChange;
+    }
+
+    // The number of vertices above the slice
+    label n = i;
+
+    // If there are more positives than negatives then reverse the order so that
+    // the negatives are at the start
+    if (n > 2)
+    {
+        n = 4 - n;
+        for (label i = 0; i < 2; ++ i)
+        {
+            Swap(indices[i], indices[3-i]);
+        }
+    }
+
+    // Correct the sign
+    if (signChange)
+    {
+        Swap(indices[2], indices[3]);
+    }
+
+    // Permute the data
+    const FixedList<Point, 4> p = tetReorder(tet, indices);
+    const FixedList<scalar, 4> l = tetReorder(level, indices);
+    AboveOp a = tetReorder(aboveOp, indices);
+    BelowOp b = tetReorder(belowOp, indices);
+
+    // Calculate the integrals above and below the level set
+    if (n == 1)
+    {
+        // Slice off one corner to form a tet and a prism
+        FixedList<scalar, 3> f;
+        for (label i = 0; i < 3; ++ i)
+        {
+            f[i] = l[0]/(l[0] - l[i+1]);
+        }
+        if (l[0] > 0)
+        {
+            return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (
+                tetCutTetAndCentre(a, p, f).first() + tetCutPrism0AndCentre(b, p, f).first(),
+                (
+                    tetCutTetAndCentre(a, p, f).first()*tetCutTetAndCentre(a, p, f).second()
+                    + tetCutPrism0AndCentre(b, p, f).first()*tetCutPrism0AndCentre(b, p, f).second()
+                )
+              / (tetCutTetAndCentre(a, p, f).first() + tetCutPrism0AndCentre(b, p, f).first())
+            );
+        }
+        else
+        {
+            return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (
+                tetCutPrism0AndCentre(a, p, f).first() + tetCutTetAndCentre(b, p, f).first(),
+                (
+                    tetCutPrism0AndCentre(a, p, f).first()*tetCutPrism0AndCentre(a, p, f).second()
+                    + tetCutTetAndCentre(b, p, f).first()*tetCutTetAndCentre(b, p, f).second()
+                )
+              / (tetCutPrism0AndCentre(a, p, f).first() + tetCutTetAndCentre(b, p, f).first())
+            );
+        }
+    }
+    else if (n == 2)
+    {
+        // Slice off two corners to form two prisms
+        FixedList<scalar, 4> f;
+        for (label i = 0; i < 2; ++ i)
+        {
+            for (label j = 0; j < 2; ++ j)
+            {
+                f[2*i+j] = l[i]/(l[i] - l[j+2]);
+            }
+        }
+        if (l[0] > 0)
+        {
+            return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (
+                tetCutPrism01AndCentre(a, p, f).first() + tetCutPrism23AndCentre(b, p, f).first(),
+                (
+                    tetCutPrism01AndCentre(a, p, f).first()*tetCutPrism01AndCentre(a, p, f).second()
+                    + tetCutPrism23AndCentre(b, p, f).first()*tetCutPrism23AndCentre(b, p, f).second()
+                )
+              / (tetCutPrism01AndCentre(a, p, f).first() + tetCutPrism23AndCentre(b, p, f).first())
+            );
+        }
+        else
+        {
+            return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+            (
+                tetCutPrism23AndCentre(a, p, f).first() + tetCutPrism01AndCentre(b, p, f).first(),
+                (
+                    tetCutPrism23AndCentre(a, p, f).first()*tetCutPrism23AndCentre(a, p, f).second()
+                    + tetCutPrism01AndCentre(b, p, f).first()*tetCutPrism01AndCentre(b, p, f).second()
+                )
+              / (tetCutPrism23AndCentre(a, p, f).first() + tetCutPrism01AndCentre(b, p, f).first())
+            );
+        }
+    }
+
+    FatalErrorInFunction
+        << "The number of tet vertices above the level set should always be "
+        << "either 1 or 2" << exit(FatalError);
+
+    return Tuple2<typename Foam::cut::opAddResult<AboveOp, BelowOp>::type, Point>
+        (aboveOp() + belowOp(), Zero);
+}
 
 template<class AboveOp, class BelowOp>
 typename Foam::cut::opAddResult<AboveOp, BelowOp>::type Foam::tetCut