diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.C b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.C
index 2d8cdb776..06f5a7409 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.C
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.C
@@ -996,10 +996,10 @@ scalar dynamicTopoFvMesh::testProximity
     if (twoDMesh_)
     {
         // Obtain the face-normal.
-        meshOps::faceNormal(faces_[index], points_, gNormal);
+        gNormal = faces_[index].normal(points_);
 
         // Obtain the face centre.
-        meshOps::faceCentre(faces_[index], points_, gCentre);
+        gCentre = faces_[index].centre(points_);
 
         // Fetch the edge
         const edge& edgeToCheck = edges_[getTriBoundaryEdge(index)];
@@ -1034,16 +1034,7 @@ scalar dynamicTopoFvMesh::testProximity
             if (neighbour_[eFaces[faceI]] == -1)
             {
                 // Obtain the normal.
-                vector gTmp;
-
-                meshOps::faceNormal
-                (
-                    faces_[eFaces[faceI]],
-                    points_,
-                    gTmp
-                );
-
-                gNormal += gTmp;
+                gNormal += faces_[eFaces[faceI]].normal(points_);
             }
         }
 
@@ -2190,8 +2181,7 @@ const changeMap dynamicTopoFvMesh::identifySliverType
         }
 
         // Obtain the unit normal.
-        vector testNormal;
-        meshOps::faceNormal(testFace, points_, testNormal);
+        vector testNormal = testFace.normal(points_);
 
         testNormal /= (mag(testNormal) + VSMALL);
 
@@ -2213,8 +2203,7 @@ const changeMap dynamicTopoFvMesh::identifySliverType
     }
 
     // Obtain the face-normal.
-    vector refArea;
-    meshOps::faceNormal(tFace, points_, refArea);
+    vector refArea = tFace.normal(points_);
 
     // Normalize it.
     vector n = refArea/mag(refArea);
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.H b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.H
index d9838ffa1..e4256d135 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.H
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMesh.H
@@ -50,6 +50,7 @@ SourceFiles
 
 #include "Switch.H"
 #include "tetMetric.H"
+#include "DynamicField.H"
 #include "threadHandler.H"
 #include "dynamicFvMesh.H"
 
@@ -124,7 +125,14 @@ class dynamicTopoFvMesh
             typedef DynamicList<T, 0, 11, 10> Type;
         };
 
-        resizableList<point>::Type oldPoints_, points_;
+        template <class T>
+        class resizableField
+        {
+        public:
+            typedef DynamicField<T> Type;
+        };
+
+        resizableField<point>::Type oldPoints_, points_;
         resizableList<face>::Type faces_;
         resizableList<label>::Type owner_, neighbour_;
         resizableList<cell>::Type cells_;
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C
index 42c809080..9d8ec2d24 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshCheck.C
@@ -276,12 +276,7 @@ bool dynamicTopoFvMesh::checkBoundingCurve(const label eIndex) const
         if ((fPatch = whichPatch(edgeFaces[faceI])) > -1)
         {
             // Obtain the normal.
-            meshOps::faceNormal
-            (
-                faces_[edgeFaces[faceI]],
-                points_,
-                fNorm[count]
-            );
+            fNorm[count] = faces_[edgeFaces[faceI]].normal(points_);
 
             // Normalize it.
             fNorm[count] /= mag(fNorm[count]) + VSMALL;
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshI.H b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshI.H
index fbdce588e..48996badd 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshI.H
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshI.H
@@ -407,12 +407,7 @@ inline scalar dynamicTopoFvMesh::edgeLengthScale
                 if (neighbour_[eFaces[faceI]] == -1)
                 {
                     // Obtain the normal.
-                    meshOps::faceNormal
-                    (
-                        faces_[eFaces[faceI]],
-                        points_,
-                        fNorm[count]
-                    );
+                    fNorm[count] = faces_[eFaces[faceI]].normal(points_);
 
                     // Normalize it.
                     fNorm[count] /= mag(fNorm[count]);
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshMapping.C b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshMapping.C
index a45855703..94e829d22 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshMapping.C
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/dynamicTopoFvMeshMapping.C
@@ -319,7 +319,7 @@ void dynamicTopoFvMesh::computeParents
     {
         offset = boundaryMesh()[whichPatch(index)].start();
 
-        meshOps::faceCentre(faces_[index], oldPoints_, centre);
+        centre = faces_[index].centre(oldPoints_);
     }
     else
     if (dimension == 3)
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeBisect.C b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeBisect.C
index bdb6e3cf2..44ee279a6 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeBisect.C
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeBisect.C
@@ -1615,26 +1615,9 @@ const changeMap dynamicTopoFvMesh::bisectQuadFace
                 // Fetch face-normals
                 vector tfNorm, f0Norm, f1Norm;
 
-                meshOps::faceNormal
-                (
-                    faces_[newFaceIndex[faceI]],
-                    oldPoints_,
-                    tfNorm
-                );
-
-                meshOps::faceNormal
-                (
-                    faces_[c0BdyIndex[0]],
-                    oldPoints_,
-                    f0Norm
-                );
-
-                meshOps::faceNormal
-                (
-                    faces_[c0BdyIndex[1]],
-                    oldPoints_,
-                    f1Norm
-                );
+                tfNorm = faces_[newFaceIndex[faceI]].normal(oldPoints_);
+                f0Norm = faces_[c0BdyIndex[0]].normal(oldPoints_);
+                f1Norm = faces_[c0BdyIndex[1]].normal(oldPoints_);
 
                 // Tri-face on boundary. Perform normal checks
                 // also, because of empty patches.
@@ -4232,7 +4215,7 @@ scalar dynamicTopoFvMesh::computeTrisectionQuality
     point midPoint;
 
     // Fetch the midPoint
-    meshOps::faceCentre(faces_[fIndex], points_, midPoint);
+    midPoint = faces_[fIndex].centre(points_);
 
     FixedList<label,2> apexPoint(-1);
 
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeCollapse.C b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeCollapse.C
index 191f32821..99c1382b2 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeCollapse.C
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeCollapse.C
@@ -400,20 +400,20 @@ const changeMap dynamicTopoFvMesh::collapseQuadFace
             // Compute face position / normal
             if (c0BdyFace[0].which(original[0]) > -1)
             {
-                meshOps::faceCentre(c0BdyFace[0], oldPoints_, xf[0]);
-                meshOps::faceNormal(c0BdyFace[0], oldPoints_, nf[0]);
+                xf[0] = c0BdyFace[0].centre(oldPoints_);
+                nf[0] = c0BdyFace[0].normal(oldPoints_);
 
-                meshOps::faceCentre(c0BdyFace[1], oldPoints_, xf[1]);
-                meshOps::faceNormal(c0BdyFace[1], oldPoints_, nf[1]);
+                xf[1] = c0BdyFace[1].centre(oldPoints_);
+                nf[1] = c0BdyFace[1].normal(oldPoints_);
             }
             else
             if (c0BdyFace[1].which(original[0]) > -1)
             {
-                meshOps::faceCentre(c0BdyFace[1], oldPoints_, xf[0]);
-                meshOps::faceNormal(c0BdyFace[1], oldPoints_, nf[0]);
+                xf[0] = c0BdyFace[1].centre(oldPoints_);
+                nf[0] = c0BdyFace[1].normal(oldPoints_);
 
-                meshOps::faceCentre(c0BdyFace[0], oldPoints_, xf[1]);
-                meshOps::faceNormal(c0BdyFace[0], oldPoints_, nf[1]);
+                xf[1] = c0BdyFace[0].centre(oldPoints_);
+                nf[1] = c0BdyFace[0].normal(oldPoints_);
             }
             else
             {
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeSwap.C b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeSwap.C
index f8d37a7e5..2c879cf70 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeSwap.C
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/edgeSwap.C
@@ -370,19 +370,8 @@ const changeMap dynamicTopoFvMesh::swapQuadFace
         // Assume that centre-plane passes through the origin.
         vector xf, nf;
 
-        meshOps::faceCentre
-        (
-            triFaceOldPoints[faceI],
-            oldPoints_,
-            xf
-        );
-
-        meshOps::faceNormal
-        (
-            triFaceOldPoints[faceI],
-            oldPoints_,
-            nf
-        );
+        xf = triFaceOldPoints[faceI].centre(oldPoints_);
+        nf = triFaceOldPoints[faceI].normal(oldPoints_);
 
         if ((((xf & n) * n) & nf) < 0.0)
         {
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOps.H b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOps.H
index 5b3c96052..a32f17c89 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOps.H
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOps.H
@@ -67,22 +67,6 @@ class polyMesh;
 
 namespace meshOps
 {
-    // Compute the centre for a given face, using UList
-    inline void faceCentre
-    (
-        const face& faceToCheck,
-        const UList<vector>& points,
-        vector& centre
-    );
-
-    // Compute the normal for a given face, using UList
-    inline void faceNormal
-    (
-        const face& faceToCheck,
-        const UList<vector>& points,
-        vector& normal
-    );
-
     // Method to find the common-edge between two faces.
     inline bool findCommonEdge
     (
@@ -168,16 +152,15 @@ namespace meshOps
     );
 
     // Given a cell index, find the centroid / volume of a cell.
-    template <class T1, class T2>
     inline bool cellCentreAndVolume
     (
         const label cIndex,
-        const UList<Vector<T1> >& points,
+        const vectorField& points,
         const UList<face>& faces,
         const UList<cell>& cells,
         const UList<label>& owner,
-        Vector<T2>& centre,
-        T2& volume,
+        vector& centre,
+        scalar& volume,
         bool checkClosed = false
     );
 
diff --git a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOpsI.H b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOpsI.H
index 95ea37aff..e71ca805d 100644
--- a/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOpsI.H
+++ b/src/dynamicMesh/dynamicFvMesh/dynamicTopoFvMesh/meshOpsI.H
@@ -47,127 +47,6 @@ namespace Foam
 namespace meshOps
 {
 
-// Compute the centroid for a given face, using UList
-inline void faceCentre
-(
-    const face& faceToCheck,
-    const UList<vector>& points,
-    vector& centre
-)
-{
-    // Reset to zero
-    centre = vector::zero;
-
-    vector a(vector::zero);
-    vector b(vector::zero);
-    vector c(vector::zero);
-
-    // If the face is a triangle, do a direct calculation
-    // to avoid round-off error-related problems
-    if (faceToCheck.size() == 3)
-    {
-        a = points[faceToCheck[0]];
-        b = points[faceToCheck[1]];
-        c = points[faceToCheck[2]];
-
-        centre = ((1.0 / 3.0) * (a + b + c));
-
-        return;
-    }
-
-    label nPoints = faceToCheck.size();
-
-    register label pI;
-
-    // Store the centre point in c
-    for (pI = 0; pI < nPoints; pI++)
-    {
-        c += points[faceToCheck[pI]];
-    }
-
-    c /= nPoints;
-
-    scalar ta = 0.0;
-    scalar sumA = 0.0;
-    vector ttc = vector::zero;
-    vector sumAc = vector::zero;
-
-    for (pI = 0; pI < nPoints; pI++)
-    {
-        a = points[faceToCheck[pI]];
-        b = points[faceToCheck[(pI + 1) % nPoints]];
-
-        // Calculate 3*triangle centre
-        ttc = (a + b + c);
-
-        // Calculate 2*triangle area
-        ta = mag((a - c)^(b - c));
-
-        sumA += ta;
-        sumAc += ta*ttc;
-    }
-
-    if (sumA > VSMALL)
-    {
-        centre = (sumAc/(3.0 * sumA));
-    }
-    else
-    {
-        centre = c;
-    }
-}
-
-
-// Compute the normal for a given face, using UList
-inline void faceNormal
-(
-    const face& faceToCheck,
-    const UList<vector>& points,
-    vector& normal
-)
-{
-    // Reset to zero
-    normal = vector::zero;
-
-    vector a(vector::zero);
-    vector b(vector::zero);
-    vector c(vector::zero);
-
-    // If the face is a triangle, do a direct calculation
-    // to avoid round-off error-related problems
-    if (faceToCheck.size() == 3)
-    {
-        a = points[faceToCheck[0]];
-        b = points[faceToCheck[1]];
-        c = points[faceToCheck[2]];
-
-        normal = (0.5 * ((b - a)^(c - a)));
-
-        return;
-    }
-
-    label nPoints = faceToCheck.size();
-
-    register label pI;
-
-    // Store the centre point in c
-    for (pI = 0; pI < nPoints; pI++)
-    {
-        c += points[faceToCheck[pI]];
-    }
-
-    c /= nPoints;
-
-    for (pI = 0; pI < nPoints; pI++)
-    {
-        a = points[faceToCheck[pI]];
-        b = points[faceToCheck[(pI + 1) % nPoints]];
-
-        normal += (0.5 * ((b - a)^(c - a)));
-    }
-}
-
-
 // Utility method to find the common edge between two faces.
 inline bool findCommonEdge
 (
@@ -393,35 +272,27 @@ inline label tetApexPoint
 
 // Given a cell index, find the centroid / volume of a cell.
 //  - If checking is enabled, return whether the cell is closed
-template <class T1, class T2>
 inline bool cellCentreAndVolume
 (
     const label cIndex,
-    const UList<Vector<T1> >& points,
+    const vectorField& points,
     const UList<face>& faces,
     const UList<cell>& cells,
     const UList<label>& owner,
-    Vector<T2>& centre,
-    T2& volume,
+    vector& centre,
+    scalar& volume,
     bool checkClosed
 )
 {
     // Reset inputs
-    volume = pTraits<T2>::zero;
-    centre = Vector<T2>::zero;
+    volume = 0.0;
+    centre = vector::zero;
 
     const cell& cellToCheck = cells[cIndex];
 
     // Average face-centres to get an estimate centroid
-    Vector<T2> cEst(Vector<T2>::zero);
-    Vector<T2> fCentre(Vector<T2>::zero);
-    Vector<T2> fArea(Vector<T2>::zero);
-    Vector<T2> sumClosed(Vector<T2>::zero), sumMagClosed(Vector<T2>::zero);
-
-    const T2 three(3.0), four(4.0);
-    const T2 oneThird = (pTraits<T2>::one / three);
-    const T2 oneFourth = (pTraits<T2>::one / four);
-    const T2 threeFourths = (three / four);
+    vector cEst(vector::zero), fCentre(vector::zero), fArea(vector::zero);
+    vector sumClosed(vector::zero), sumMagClosed(vector::zero);
 
     forAll(cellToCheck, faceI)
     {
@@ -432,17 +303,10 @@ inline bool cellCentreAndVolume
             continue;
         }
 
-        meshOps::faceCentre
-        (
-            checkFace,
-            points,
-            fCentre
-        );
-
-        cEst += fCentre;
+        cEst += checkFace.centre(points);
     }
 
-    cEst /= T2(cellToCheck.size());
+    cEst /= cellToCheck.size();
 
     forAll(cellToCheck, faceI)
     {
@@ -453,19 +317,8 @@ inline bool cellCentreAndVolume
             continue;
         }
 
-        meshOps::faceNormal
-        (
-            checkFace,
-            points,
-            fArea
-        );
-
-        meshOps::faceCentre
-        (
-            checkFace,
-            points,
-            fCentre
-        );
+        fArea = checkFace.normal(points);
+        fCentre = checkFace.centre(points);
 
         // Flip if necessary
         if (owner[cellToCheck[faceI]] != cIndex)
@@ -474,10 +327,10 @@ inline bool cellCentreAndVolume
         }
 
         // Calculate 3*face-pyramid volume
-        T2 pyr3Vol = fArea & (fCentre - cEst);
+        scalar pyr3Vol = fArea & (fCentre - cEst);
 
         // Calculate face-pyramid centre
-        Vector<T2> pc = (threeFourths*fCentre) + (oneFourth*cEst);
+        vector pc = ((3.0 / 4.0) * fCentre) + ((1.0 / 4.0) * cEst);
 
         // Accumulate volume-weighted face-pyramid centre
         centre += pyr3Vol*pc;
@@ -494,12 +347,12 @@ inline bool cellCentreAndVolume
     }
 
     centre /= volume + VSMALL;
-    volume *= oneThird;
+    volume *= (1.0 / 3.0);
 
     if (checkClosed)
     {
         // Check the sum across components
-        T2 closed = pTraits<T2>::zero;
+        scalar closed = 0.0;
 
         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
         {
diff --git a/tutorials/incompressible/icoDyMFoam/movingBlockRBF/RBFMotionFunction/Make/options b/tutorials/incompressible/icoDyMFoam/movingBlockRBF/RBFMotionFunction/Make/options
index aff544616..632ec8263 100644
--- a/tutorials/incompressible/icoDyMFoam/movingBlockRBF/RBFMotionFunction/Make/options
+++ b/tutorials/incompressible/icoDyMFoam/movingBlockRBF/RBFMotionFunction/Make/options
@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/meshMotion/RBFMotionSolver/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude