diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
index d4492f3eb4..1486db9280 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C
@@ -85,73 +85,11 @@ void Foam::primitiveMesh::makeFaceCentresAndAreas
 
     forAll(fs, facei)
     {
-        const labelList& f = fs[facei];
-        label nPoints = f.size();
-
-        // If the face is a triangle, do a direct calculation for efficiency
-        // and to avoid round-off error-related problems
-        if (nPoints == 3)
-        {
-            fCtrs[facei] = (1.0/3.0)*(p[f[0]] + p[f[1]] + p[f[2]]);
-            fAreas[facei] = 0.5*((p[f[1]] - p[f[0]])^(p[f[2]] - p[f[0]]));
-        }
-
-        // For more complex faces, decompose into triangles
-        else
-        {
-            // Compute an estimate of the centre as the average of the points
-            point pAvg = p[f[0]];
-            for (label pi = 1; pi < nPoints; pi++)
-            {
-                pAvg += p[f[pi]];
-            }
-            pAvg /= nPoints;
-
-            // Compute the face area normal and unit normal by summing up the
-            // normals of the triangles formed by connecting each edge to the
-            // point average.
-            vector sumA = Zero;
-            forAll(f, i)
-            {
-                const vector a =
-                    (p[f[f.fcIndex(i)]] - p[f[i]])^(pAvg - p[f[i]]);
-
-                sumA += a;
-            }
-            const vector sumAHat = normalised(sumA);
-
-            // Compute the area-weighted sum of the triangle centres. Note use
-            // the triangle area projected in the direction of the face normal
-            // as the weight, *not* the triangle area magnitude. Only the
-            // former makes the calculation independent of the initial estimate.
-            scalar sumAn = 0.0;
-            vector sumAnc = Zero;
-            forAll(f, i)
-            {
-                const vector a =
-                    (p[f[f.fcIndex(i)]] - p[f[i]])^(pAvg - p[f[i]]);
-                const vector c = p[f[i]] + p[f[f.fcIndex(i)]] + pAvg;
-
-                const scalar an = a & sumAHat;
-
-                sumAn += an;
-                sumAnc += an*c;
-            }
-
-            // Complete calculating centres and areas. If the face is too small
-            // for the sums to be reliably divided then just set the centre to
-            // the initial estimate.
-            if (sumAn > vSmall)
-            {
-                fCtrs[facei] = (1.0/3.0)*sumAnc/sumAn;
-            }
-            else
-            {
-                fCtrs[facei] = pAvg;
-            }
-            fAreas[facei] = 0.5*sumA;
-        }
+        const Tuple2<vector, point> areaAndCentre =
+            face::areaAndCentre(UIndirectList<point>(p, fs[facei]));
 
+        fCtrs[facei] = areaAndCentre.second();
+        fAreas[facei] = areaAndCentre.first();
         magfAreas[facei] = max(mag(fAreas[facei]), rootVSmall);
     }
 }