diff --git a/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C b/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C
index 7f695abffb1b7eef3944c0f5d21999589361f6f7..769f3e79d97b5c0525d92dda090dee05d79ad4c2 100644
--- a/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C
+++ b/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C
@@ -41,11 +41,14 @@ Usage
         Specify the time to search from and apply the transformation
         (default is latest)
 
+    -recentre
+        Recentre using the bounding box centre before other operations
+
     -translate vector
-        Translates the points by the given vector before rotations
+        Translate the points by the given vector before rotations
 
     -rotate (vector vector)
-        Rotates the points from the first vector to the second,
+        Rotate the points from the first vector to the second
 
     -rotate-angle (vector angle)
         Rotate angle degrees about vector axis.
@@ -54,7 +57,7 @@ Usage
      or -rollPitchYaw (rolldegrees pitchdegrees yawdegrees)
 
     -scale scalar|vector
-        Scales the points by the given scalar or vector.
+        Scale the points by the given scalar or vector on output.
 
     The any or all of the three options may be specified and are processed
     in the above order.
@@ -150,6 +153,85 @@ void rotateFields(const argList& args, const Time& runTime, const tensor& T)
 }
 
 
+// Retrieve scaling option
+// - size 0 : no scaling
+// - size 1 : uniform scaling
+// - size 3 : non-uniform scaling
+List<scalar> getScalingOpt(const word& optName, const argList& args)
+{
+    // readListIfPresent handles single or multiple values
+    // - accept 1 or 3 values
+
+    List<scalar> scaling;
+    args.readListIfPresent(optName, scaling);
+
+    if (scaling.size() == 1)
+    {
+        // Uniform scaling
+    }
+    else if (scaling.size() == 3)
+    {
+        // Non-uniform, but may actually be uniform
+        if
+        (
+            equal(scaling[0], scaling[1])
+         && equal(scaling[0], scaling[2])
+        )
+        {
+            scaling.resize(1);
+        }
+    }
+    else if (!scaling.empty())
+    {
+        FatalError
+            << "Incorrect number of components, must be 1 or 3." << nl
+            << "    -" << optName << ' ' << args[optName].c_str() << endl
+            << exit(FatalError);
+    }
+
+    if (scaling.size() == 1 && equal(scaling[0], 1))
+    {
+        // Scale factor 1 == no scaling
+        scaling.clear();
+    }
+
+    for (const scalar scale : scaling)
+    {
+        if (scale <= 0)
+        {
+            FatalError
+                << "Invalid scaling value, must be positive." << nl
+                << "    -" << optName << ' ' << args[optName].c_str() << endl
+                << exit(FatalError);
+        }
+    }
+
+    return scaling;
+}
+
+
+void applyScaling(pointField& points, const List<scalar>& scaling)
+{
+    if (scaling.size() == 1)
+    {
+        Info<< "Scaling points uniformly by " << scaling[0] << nl;
+        points *= scaling[0];
+    }
+    else if (scaling.size() == 3)
+    {
+        Info<< "Scaling points by ("
+            << scaling[0] << ' '
+            << scaling[1] << ' '
+            << scaling[2] << ')' << nl;
+
+        points.replace(vector::X, scaling[0]*points.component(vector::X));
+        points.replace(vector::Y, scaling[1]*points.component(vector::Y));
+        points.replace(vector::Z, scaling[2]*points.component(vector::Z));
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
@@ -165,11 +247,21 @@ int main(int argc, char *argv[])
         "Specify the time to search from and apply the transformation"
         " (default is latest)"
     );
+    argList::addBoolOption
+    (
+        "recentre",
+        "Recentre the bounding box before other operations"
+    );
     argList::addOption
     (
         "translate",
         "vector",
-        "Translate by specified <vector> - eg, '(1 0 0)' before rotations"
+        "Translate by specified <vector> before rotations"
+    );
+    argList::addBoolOption
+    (
+        "auto-origin",
+        "Use bounding box centre as origin for rotations"
     );
     argList::addOption
     (
@@ -214,6 +306,9 @@ int main(int argc, char *argv[])
         "use either '(0.001 0.001 0.001)' or simply '0.001'"
     );
 
+    // Compatibility with surfaceTransformPoints
+    argList::addOptionCompat("scale", {"write-scale", 0});
+
     #include "addRegionOption.H"
     #include "setRootCase.H"
 
@@ -223,6 +318,7 @@ int main(int argc, char *argv[])
     {
         const List<word> operationNames
         {
+            "recentre",
             "translate",
             "rotate",
             "rotate-angle",
@@ -286,16 +382,31 @@ int main(int argc, char *argv[])
         )
     );
 
+
+    // Begin operations
+
     vector v;
+    if (args.found("recentre"))
+    {
+        v = boundBox(points).centre();
+        Info<< "Adjust centre " << v << " -> (0 0 0)" << endl;
+        points -= v;
+    }
+
     if (args.readIfPresent("translate", v))
     {
         Info<< "Translating points by " << v << endl;
-
         points += v;
     }
 
     vector origin;
-    const bool useOrigin = args.readIfPresent("origin", origin);
+    bool useOrigin = args.readIfPresent("origin", origin);
+    if (args.found("auto-origin") && !useOrigin)
+    {
+        useOrigin = true;
+        origin = boundBox(points).centre();
+    }
+
     if (useOrigin)
     {
         Info<< "Set origin for rotations to " << origin << endl;
@@ -328,8 +439,8 @@ int main(int argc, char *argv[])
             args.lookup("rotate-angle")()
         );
 
-        const vector& axis  = rotAxisAngle.first();
-        const scalar& angle = rotAxisAngle.second();
+        const vector& axis = rotAxisAngle.first();
+        const scalar angle = rotAxisAngle.second();
 
         Info<< "Rotating points " << nl
             << "    about " << axis << nl
@@ -380,35 +491,8 @@ int main(int argc, char *argv[])
         }
     }
 
-    List<scalar> scaling;
-    if (args.readListIfPresent("scale", scaling))
-    {
-        // readListIfPresent handles single or multiple values
-
-        if (scaling.size() == 1)
-        {
-            Info<< "Scaling points uniformly by " << scaling[0] << nl;
-            points *= scaling[0];
-        }
-        else if (scaling.size() == 3)
-        {
-            Info<< "Scaling points by ("
-                << scaling[0] << " "
-                << scaling[1] << " "
-                << scaling[2] << ")" << nl;
-
-            points.replace(vector::X, scaling[0]*points.component(vector::X));
-            points.replace(vector::Y, scaling[1]*points.component(vector::Y));
-            points.replace(vector::Z, scaling[2]*points.component(vector::Z));
-        }
-        else
-        {
-            FatalError
-                << "-scale with 1 or 3 components only" << nl
-                << "given: " << args["scale"] << endl
-                << exit(FatalError);
-        }
-    }
+    // Output scaling
+    applyScaling(points, getScalingOpt("scale", args));
 
     if (useOrigin)
     {
diff --git a/applications/utilities/surface/surfaceRefineRedGreen/surfaceRefineRedGreen.C b/applications/utilities/surface/surfaceRefineRedGreen/surfaceRefineRedGreen.C
index b5b6b0cd30a9d37c330181685b4dc7372621480a..b92418e1d43d08c9ce0d29286b041851767c71db 100644
--- a/applications/utilities/surface/surfaceRefineRedGreen/surfaceRefineRedGreen.C
+++ b/applications/utilities/surface/surfaceRefineRedGreen/surfaceRefineRedGreen.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -63,33 +64,52 @@ int main(int argc, char *argv[])
     argList::noParallel();
     argList::addArgument("input", "The input surface file");
     argList::addArgument("output", "The output surface file");
+    argList::addOption
+    (
+        "steps",
+        "N",
+        "Number of refinement steps (default: 1)"
+    );
     argList args(argc, argv);
 
-    const fileName surfFileName = args[1];
-    const fileName outFileName  = args[2];
+    const fileName surfFileName(args[1]);
+    const fileName outFileName(args[2]);
 
     Info<< "Reading surface from " << surfFileName << " ..." << endl;
 
-    triSurface surf1(surfFileName);
+    triSurface surf(surfFileName);
 
-    // Refine
-    triSurface surf2 = triSurfaceTools::redGreenRefine
-    (
-        surf1,
-        identity(surf1.size())  //Hack: refine all
-    );
+    Info<< "Original surface:" << nl
+        << "    triangles     :" << surf.size() << nl
+        << "    vertices(used):" << surf.nPoints() << endl;
+
+
+    const label nsteps =
+        args.getCheckOrDefault<label>("steps", 1, labelMinMax::ge(1));
+
+
+    Info<< "Refining " << nsteps << " times" << flush;
+
+    for (label step = 0; step < nsteps; ++step)
+    {
+        Info<< '.' << flush;
 
-    Info<< "Original surface:" << endl
-        << "    triangles     :" << surf1.size() << endl
-        << "    vertices(used):" << surf1.nPoints() << endl
-        << "Refined surface:" << endl
-        << "    triangles     :" << surf2.size() << endl
-        << "    vertices(used):" << surf2.nPoints() << endl << endl;
+        surf = triSurfaceTools::redGreenRefine
+        (
+            surf,
+            identity(surf.size())  // Refine all
+        );
+    }
+    Info<< endl;
 
+    Info<< "Refined surface:" << nl
+        << "    triangles     :" << surf.size() << nl
+        << "    vertices(used):" << surf.nPoints() << endl;
 
-    Info<< "Writing refined surface to " << outFileName << " ..." << endl;
+    Info<< nl
+        << "Writing refined surface to " << outFileName << " ..." << endl;
 
-    surf2.write(outFileName);
+    surf.write(outFileName);
 
     Info<< "End\n" << endl;
 
diff --git a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
index a3ba6687bcbeca4555f016bd454295cee14a4de1..a2302190ce5295578da0886d30de91e8d53b4f75 100644
--- a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
+++ b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
@@ -93,6 +93,84 @@ static bool hasReadWriteTypes(const word& readType, const word& writeType)
 }
 
 
+// Retrieve scaling option
+// - size 0 : no scaling
+// - size 1 : uniform scaling
+// - size 3 : non-uniform scaling
+List<scalar> getScalingOpt(const word& optName, const argList& args)
+{
+    // readListIfPresent handles single or multiple values
+    // - accept 1 or 3 values
+
+    List<scalar> scaling;
+    args.readListIfPresent(optName, scaling);
+
+    if (scaling.size() == 1)
+    {
+        // Uniform scaling
+    }
+    else if (scaling.size() == 3)
+    {
+        // Non-uniform, but may actually be uniform
+        if
+        (
+            equal(scaling[0], scaling[1])
+         && equal(scaling[0], scaling[2])
+        )
+        {
+            scaling.resize(1);
+        }
+    }
+    else if (!scaling.empty())
+    {
+        FatalError
+            << "Incorrect number of components, must be 1 or 3." << nl
+            << "    -" << optName << ' ' << args[optName].c_str() << endl
+            << exit(FatalError);
+    }
+
+    if (scaling.size() == 1 && equal(scaling[0], 1))
+    {
+        // Scale factor 1 == no scaling
+        scaling.clear();
+    }
+
+    for (const scalar scale : scaling)
+    {
+        if (scale <= 0)
+        {
+            FatalError
+                << "Invalid scaling value, must be positive." << nl
+                << "    -" << optName << ' ' << args[optName].c_str() << endl
+                << exit(FatalError);
+        }
+    }
+
+    return scaling;
+}
+
+
+void applyScaling(pointField& points, const List<scalar>& scaling)
+{
+    if (scaling.size() == 1)
+    {
+        Info<< "Scaling points uniformly by " << scaling[0] << nl;
+        points *= scaling[0];
+    }
+    else if (scaling.size() == 3)
+    {
+        Info<< "Scaling points by ("
+            << scaling[0] << ' '
+            << scaling[1] << ' '
+            << scaling[2] << ')' << nl;
+
+        points.replace(vector::X, scaling[0]*points.component(vector::X));
+        points.replace(vector::Y, scaling[1]*points.component(vector::Y));
+        points.replace(vector::Z, scaling[2]*points.component(vector::Z));
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -106,11 +184,21 @@ int main(int argc, char *argv[])
     argList::noParallel();
     argList::addArgument("input", "The input surface file");
     argList::addArgument("output", "The output surface file");
+    argList::addBoolOption
+    (
+        "recentre",
+        "Recentre the bounding box before other operations"
+    );
     argList::addOption
     (
         "translate",
         "vector",
-        "Translate by specified <vector> - eg, '(1 0 0)' before rotations"
+        "Translate by specified <vector> before rotations"
+    );
+    argList::addBoolOption
+    (
+        "auto-origin",
+        "Use bounding box centre as origin for rotations"
     );
     argList::addOption
     (
@@ -144,10 +232,15 @@ int main(int argc, char *argv[])
     );
     argList::addOption
     (
-        "scale",
+        "read-scale",
+        "scalar | vector",
+        "Uniform or non-uniform input scaling"
+    );
+    argList::addOption
+    (
+        "write-scale",
         "scalar | vector",
-        "Scale by the specified amount - Eg, for uniform [mm] to [m] scaling "
-        "use either '(0.001 0.001 0.001)' or simply '0.001'"
+        "Uniform or non-uniform output scaling"
     );
     argList::addOption
     (
@@ -162,18 +255,23 @@ int main(int argc, char *argv[])
         "Output format (default: use file extension)"
     );
 
+    // Backward compatibility and with transformPoints
+    argList::addOptionCompat("write-scale", {"scale", -2006});
+
     argList args(argc, argv);
 
     // Verify that an operation has been specified
     {
         const List<word> operationNames
         {
+            "recentre",
             "translate",
             "rotate",
             "rotate-angle",
             "rollPitchYaw",
             "yawPitchRoll",
-            "scale"
+            "read-scale",
+            "write-scale"
         };
 
         if (!args.count(operationNames))
@@ -225,16 +323,34 @@ int main(int argc, char *argv[])
 
     pointField points(surf1.points());
 
+
+    // Begin operations
+
+    // Input scaling
+    applyScaling(points, getScalingOpt("read-scale", args));
+
     vector v;
+    if (args.found("recentre"))
+    {
+        v = boundBox(points).centre();
+        Info<< "Adjust centre " << v << " -> (0 0 0)" << endl;
+        points -= v;
+    }
+
     if (args.readIfPresent("translate", v))
     {
         Info<< "Translating points by " << v << endl;
-
         points += v;
     }
 
-    vector origin(Zero);
-    const bool useOrigin = args.readIfPresent("origin", origin);
+    vector origin;
+    bool useOrigin = args.readIfPresent("origin", origin);
+    if (args.found("auto-origin") && !useOrigin)
+    {
+        useOrigin = true;
+        origin = boundBox(points).centre();
+    }
+
     if (useOrigin)
     {
         Info<< "Set origin for rotations to " << origin << endl;
@@ -262,8 +378,8 @@ int main(int argc, char *argv[])
             args.lookup("rotate-angle")()
         );
 
-        const vector& axis  = rotAxisAngle.first();
-        const scalar& angle = rotAxisAngle.second();
+        const vector& axis = rotAxisAngle.first();
+        const scalar angle = rotAxisAngle.second();
 
         Info<< "Rotating points " << nl
             << "    about " << axis << nl
@@ -299,44 +415,8 @@ int main(int argc, char *argv[])
         points = transform(rot, points);
     }
 
-    List<scalar> scaling;
-    if (args.readListIfPresent("scale", scaling))
-    {
-        // readListIfPresent handles single or multiple values
-
-        if
-        (
-            scaling.size() == 1
-         ||
-            (
-                scaling.size() == 3
-             && equal(scaling[0], scaling[1])
-             && equal(scaling[0], scaling[2])
-            )
-        )
-        {
-            Info<< "Scaling points uniformly by " << scaling[0] << nl;
-            points *= scaling[0];
-        }
-        else if (scaling.size() == 3)
-        {
-            Info<< "Scaling points by ("
-                << scaling[0] << ' '
-                << scaling[1] << ' '
-                << scaling[2] << ')' << nl;
-
-            points.replace(vector::X, scaling[0]*points.component(vector::X));
-            points.replace(vector::Y, scaling[1]*points.component(vector::Y));
-            points.replace(vector::Z, scaling[2]*points.component(vector::Z));
-        }
-        else
-        {
-            FatalError
-                << "-scale with 1 or 3 components only" << nl
-                << "given: " << args["scale"] << endl
-                << exit(FatalError);
-        }
-    }
+    // Output scaling
+    applyScaling(points, getScalingOpt("write-scale", args));
 
     if (useOrigin)
     {
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/Allclean b/tutorials/compressible/rhoSimpleFoam/squareBend/Allclean
new file mode 100755
index 0000000000000000000000000000000000000000..48e77c0ce28c718866a419c8a8d85755b121c43e
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/Allclean
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase
+
+rm -f constant/triSurface/inlet_sample*.obj
+rm -f constant/triSurface/oversized_sample.obj
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/Allrun b/tutorials/compressible/rhoSimpleFoam/squareBend/Allrun
new file mode 100755
index 0000000000000000000000000000000000000000..2a7e9ad86bc6b285357c7e06427b66a81b648df0
--- /dev/null
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/Allrun
@@ -0,0 +1,27 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+runApplication blockMesh
+
+# Extract some faces for sampling
+runApplication surfaceMeshExtract \
+    -patches inlet -time 0 \
+    constant/triSurface/inlet_sample0.obj
+
+# Recentre and inflate
+runApplication surfaceTransformPoints \
+    -recentre \
+    -scale 3.5 \
+    constant/triSurface/inlet_sample0.obj \
+    constant/triSurface/inlet_sample1.obj
+
+# Finer mesh
+runApplication surfaceRefineRedGreen \
+    constant/triSurface/inlet_sample1.obj \
+    constant/triSurface/oversized_sample.obj
+
+runApplication $(getApplication)
+
+#------------------------------------------------------------------------------