diff --git a/applications/utilities/postProcessing/lumped/lumpedPointForces/lumpedPointForces.C b/applications/utilities/postProcessing/lumped/lumpedPointForces/lumpedPointForces.C
index 8e7423557a..0ba2b9a0b3 100644
--- a/applications/utilities/postProcessing/lumped/lumpedPointForces/lumpedPointForces.C
+++ b/applications/utilities/postProcessing/lumped/lumpedPointForces/lumpedPointForces.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,8 +27,8 @@ Application
     lumpedPointForces
 
 Description
-    Extract force/moment information from existing calculations based
-    on the segmentation of the pressure integration zones.
+    Extract force/moment information from simulation results that
+    use the lumped points movement description.
 
 \*---------------------------------------------------------------------------*/
 
@@ -37,7 +37,7 @@ Description
 #include "timeSelector.H"
 #include "volFields.H"
 #include "IOobjectList.H"
-
+#include "foamVtkSeriesWriter.H"
 #include "lumpedPointTools.H"
 #include "lumpedPointIOMovement.H"
 
@@ -83,8 +83,8 @@ int main(int argc, char *argv[])
 {
     argList::addNote
     (
-        "Extract force/moment information from existing calculations based"
-        " on the lumpedPoints pressure zones."
+        "Extract force/moment information from simulation results that"
+        " use the lumped points movement description."
     );
 
     argList::addBoolOption
@@ -102,31 +102,33 @@ int main(int argc, char *argv[])
     const bool withVTK = args.found("vtk");
 
     #include "createTime.H"
+
     instantList timeDirs = timeSelector::select0(runTime, args);
 
     #include "createNamedMesh.H"
 
-    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New
-    (
-        runTime
-    );
+    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New(mesh);
 
-    if (!movement.valid())
+    if (!movement)
     {
-        Info<< "no valid movement given" << endl;
+        Info<< "No valid movement found" << endl;
         return 1;
     }
 
-    const labelList patchLst = lumpedPointTools::lumpedPointPatchList(mesh);
-    if (patchLst.empty())
+    const label nPatches = lumpedPointTools::setPatchControls(mesh);
+    if (!nPatches)
     {
-        Info<< "no patch list found" << endl;
+        Info<< "No point patches with lumped movement found" << endl;
         return 2;
     }
 
-    movement().setMapping(mesh, patchLst, lumpedPointTools::points0Field(mesh));
+    Info<<"Lumped point patch controls set on " << nPatches
+        << " patches" << nl;
 
+
+    vtk::seriesWriter forceSeries;
     List<vector> forces, moments;
+
     forAll(timeDirs, timei)
     {
         runTime.setTime(timeDirs[timei], timei);
@@ -164,11 +166,21 @@ int main(int argc, char *argv[])
                     forces,
                     moments
                 );
+
+                forceSeries.append(runTime.timeIndex(), outputName);
             }
         }
     }
 
-    Info<< "End\n" << endl;
+
+    // Create file series
+
+    if (forceSeries.size())
+    {
+        forceSeries.write("forces.vtp");
+    }
+
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }
diff --git a/applications/utilities/postProcessing/lumped/lumpedPointMovement/lumpedPointMovement.C b/applications/utilities/postProcessing/lumped/lumpedPointMovement/lumpedPointMovement.C
index a6ca67f851..1925d5d569 100644
--- a/applications/utilities/postProcessing/lumped/lumpedPointMovement/lumpedPointMovement.C
+++ b/applications/utilities/postProcessing/lumped/lumpedPointMovement/lumpedPointMovement.C
@@ -40,15 +40,50 @@ Description
 
 #include "argList.H"
 #include "Time.H"
-#include "timeSelector.H"
-
 #include "OFstream.H"
-
+#include "foamVtkSeriesWriter.H"
 #include "lumpedPointTools.H"
 #include "lumpedPointIOMovement.H"
 
 using namespace Foam;
 
+
+inline List<lumpedPointStateTuple> getResponseTable
+(
+    const fileName& file,
+    const lumpedPointState& state0
+)
+{
+    return lumpedPointTools::lumpedPointStates
+    (
+        file,
+        state0.rotationOrder(),
+        state0.degrees()
+    );
+}
+
+
+void echoTableLimits
+(
+    const List<lumpedPointStateTuple>& tbl,
+    const label span,
+    const label maxOut
+)
+{
+    Info<< "Using response table with " << tbl.size() << " entries" << nl;
+
+    if (span)
+    {
+        Info<< "Increment input by " << span << nl;
+    }
+
+    if (maxOut)
+    {
+        Info<< "Stopping after " << maxOut << " outputs" << nl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -59,7 +94,6 @@ int main(int argc, char *argv[])
         " for diagnostic purposes."
     );
 
-    argList::noParallel();
     argList::noFunctionObjects();  // Never use function objects
     argList::addOption
     (
@@ -71,7 +105,7 @@ int main(int argc, char *argv[])
     (
         "span",
         "N",
-        "Increment each input by factor N (default: 1)"
+        "Increment each input by N (default: 1)"
     );
     argList::addOption
     (
@@ -79,6 +113,17 @@ int main(int argc, char *argv[])
         "factor",
         "Relaxation/scaling factor for movement (default: 1)"
     );
+    argList::addOption
+    (
+        "visual-length",
+        "len",
+        "Visualization length for planes (visualized as triangles)"
+    );
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Test movement without a mesh"
+    );
     argList::addBoolOption
     (
         "removeLock",
@@ -96,45 +141,70 @@ int main(int argc, char *argv[])
     const label maxOut = Foam::max(0, args.getOrDefault<label>("max", 0));
     const label span   = Foam::max(1, args.getOrDefault<label>("span", 1));
 
-    const scalar relax = args.getOrDefault<scalar>("scale", 1);
-
+    // Control parameters
+    const bool dryrun = args.found("dry-run");
     const bool slave = args.found("slave");
     const bool removeLock = args.found("removeLock");
 
-    #include "createTime.H"
+    const scalar relax = args.getOrDefault<scalar>("scale", 1);
 
-    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New
-    (
-        runTime
-    );
+    args.readIfPresent("visual-length", lumpedPointState::visLength);
 
-    if (!movement.valid())
-    {
-        Info<< "no valid movement given" << endl;
-        return 1;
-    }
+    const fileName responseFile(args[1]);
 
-    List<lumpedPointStateTuple> responseTable =
-        lumpedPointTools::lumpedPointStates(args[1]);
-
-    Info<< "Using response table with " << responseTable.size()
-        << " entries" << endl;
-
-    Info << "Increment input by " << span << nl;
-
-    if (maxOut)
-    {
-        Info<< "Stopping after " << maxOut << " outputs" << endl;
-    }
+    // ----------------------------------------------------------------------
+    // Slave mode
+    // ----------------------------------------------------------------------
 
     if (slave)
     {
         Info<< "Running as slave responder" << endl;
 
-        externalFileCoupler& coupler = movement().coupler();
+        if (Pstream::parRun())
+        {
+            FatalErrorInFunction
+                << "Running as slave responder is not permitted in parallel"
+                << nl
+                << exit(FatalError);
+        }
 
-        label count = 0;
-        for (label index = 0; index < responseTable.size(); index += span)
+        #include "createTime.H"
+
+        // Create movement without a mesh
+        autoPtr<lumpedPointIOMovement> movementPtr =
+            lumpedPointIOMovement::New(runTime);
+
+        if (!movementPtr)
+        {
+            Info<< "No valid movement found" << endl;
+            return 1;
+        }
+        auto& movement = *movementPtr;
+
+        // Reference state0
+        const lumpedPointState& state0 = movement.state0();
+
+        List<lumpedPointStateTuple> responseTable =
+            getResponseTable(responseFile, state0);
+
+        echoTableLimits(responseTable, span, maxOut);
+
+        if (dryrun)
+        {
+            Info<< "dry-run: response table with " << responseTable.size()
+                << " entries" << nl
+                << "\nEnd\n" << endl;
+            return 0;
+        }
+
+        externalFileCoupler& coupler = movement.coupler();
+
+        for
+        (
+            label timei = 0, outputCount = 0;
+            timei < responseTable.size();
+            timei += span
+        )
         {
             Info<< args.executable() << ": waiting for master" << endl;
 
@@ -146,14 +216,14 @@ int main(int argc, char *argv[])
                 break;
             }
 
-            lumpedPointState state = responseTable[index].second();
-            state.relax(relax, movement().state0());
+            lumpedPointState state = responseTable[timei].second();
+            state.relax(relax, state0);
 
             // Generate input for OpenFOAM
-            OFstream os(coupler.resolveFile(movement().inputName()));
+            OFstream os(coupler.resolveFile(movement.inputName()));
             if
             (
-                movement().inputFormat()
+                movement.inputFormat()
              == lumpedPointState::inputFormatType::PLAIN
             )
             {
@@ -161,101 +231,236 @@ int main(int argc, char *argv[])
             }
             else
             {
-                os.writeEntry("time", responseTable[index].first());
+                os.writeEntry("time", responseTable[timei].first());
                 state.writeDict(os);
             }
 
             Info<< args.executable()
-                << ": updated to state " << index
+                << ": updated to state " << timei
                 << " - switch to master"
                 << endl;
 
             // Let OpenFOAM know that it can continue
             coupler.useMaster();
 
-            if (maxOut && ++count >= maxOut)
+            ++outputCount;
+
+            if (maxOut && outputCount >= maxOut)
             {
                 Info<< args.executable()
-                    <<": stopping after " << maxOut << " outputs" << endl;
+                    << ": stopping after " << maxOut << " outputs" << endl;
                 break;
             }
         }
 
         if (removeLock)
         {
-            Info<< args.executable() <<": removing lock file" << endl;
+            Info<< args.executable() << ": removing lock file" << endl;
             coupler.useSlave();  // This removes the lock-file
         }
+
+        Info<< args.executable() << ": finishing" << nl;
+
+        Info<< "\nEnd\n" << endl;
+        return 0;
     }
-    else
+
+
+    // ----------------------------------------------------------------------
+    // dry-run
+    // ----------------------------------------------------------------------
+
+    if (dryrun)
     {
-        runTime.setTime(instant(0, runTime.constant()), 0);
+        Info<< "dry-run: creating states only" << nl;
 
-        #include "createNamedPolyMesh.H"
+        #include "createTime.H"
 
-        const labelList patchLst = lumpedPointTools::lumpedPointPatchList(mesh);
-        if (patchLst.empty())
+        // Create movement without a mesh
+        autoPtr<lumpedPointIOMovement> movementPtr =
+            lumpedPointIOMovement::New(runTime);
+
+        if (!movementPtr)
         {
-            Info<< "no patch list found" << endl;
-            return 2;
+            Info<< "No valid movement found" << endl;
+            return 1;
         }
+        auto& movement = *movementPtr;
 
-        pointIOField points0 = lumpedPointTools::points0Field(mesh);
-        movement().setBoundBox(mesh, patchLst, points0);
+        // Reference state0
+        const lumpedPointState& state0 = movement.state0();
 
-        label index = 0;
+        List<lumpedPointStateTuple> responseTable =
+            getResponseTable(responseFile, state0);
 
-        // Initial geometry
-        movement().writeVTP("geom_init.vtp", mesh, patchLst, points0);
+        echoTableLimits(responseTable, span, maxOut);
 
-        forAll(responseTable, i)
+
+        vtk::seriesWriter stateSeries;
+
+        for
+        (
+            label timei = 0, outputCount = 0;
+            timei < responseTable.size();
+            timei += span
+        )
         {
-            const bool output = ((i % span) == 0);
-            lumpedPointState state = responseTable[i].second();
-            state.relax(relax, movement().state0());
+            lumpedPointState state = responseTable[timei].second();
 
-            if (output)
-            {
-                Info<<"output [" << i << "/"
-                    << responseTable.size() << "]" << endl;
-            }
-            else
-            {
-                continue;
-            }
+            state += movement.origin();
+            movement.scalePoints(state);
+            state.relax(relax, state0);
+
+            Info<< "output [" << timei << '/' << responseTable.size() << ']';
 
             // State/response = what comes back from FEM
             {
-                const word outputName = word::printf("state_%06d.vtp", index);
+                const word outputName =
+                    word::printf("state_%06d.vtp", outputCount);
 
-                Info<<"    " << outputName << endl;
+                Info<< "  " << outputName;
 
-                state.writeVTP(outputName, movement().axis());
+                movement.writeStateVTP(state, outputName);
+                stateSeries.append(outputCount, outputName);
             }
 
+            Info<< endl;
+
+            ++outputCount;
+
+            if (maxOut && outputCount >= maxOut)
             {
-                const word outputName = word::printf("geom_%06d.vtp", index);
-
-                Info<<"    " << outputName << endl;
-
-                movement().writeVTP(outputName, state, mesh, patchLst, points0);
+                Info<< "Max output " << maxOut << " ... stopping" << endl;
+                break;
             }
+        }
 
-            {
-                ++index;
+        // Write file series
 
-                bool canOutput = !maxOut || (index <= maxOut);
-                if (!canOutput)
-                {
-                    Info<<"stopping output after "
-                        << maxOut << " outputs" << endl;
-                    break;
-                }
-            }
+        if (stateSeries.size())
+        {
+            Info<< nl << "write state.vtp.series" << nl;
+            stateSeries.write("state.vtp");
+        }
+
+        Info<< "\nEnd\n" << endl;
+        return 0;
+    }
+
+
+    // ----------------------------------------------------------------------
+    // test patch movement
+    // ----------------------------------------------------------------------
+
+    #include "createTime.H"
+
+    runTime.setTime(instant(runTime.constant()), 0);
+
+    #include "createNamedMesh.H"
+
+    // Create movement with mesh
+    autoPtr<lumpedPointIOMovement> movementPtr =
+        lumpedPointIOMovement::New(mesh);
+
+    if (!movementPtr)
+    {
+        Info<< "No valid movement found" << endl;
+        return 1;
+    }
+    auto& movement = *movementPtr;
+
+    // Reference state0
+    const lumpedPointState& state0 = movement.state0();
+
+    List<lumpedPointStateTuple> responseTable =
+        getResponseTable(responseFile, state0);
+
+    echoTableLimits(responseTable, span, maxOut);
+
+    pointIOField points0(lumpedPointTools::points0Field(mesh));
+
+    const label nPatches = lumpedPointTools::setPatchControls(mesh, points0);
+    if (!nPatches)
+    {
+        Info<< "No point patches with lumped movement found" << endl;
+        return 2;
+    }
+
+    Info<< "Lumped point patch controls set on "
+        << nPatches << " patches" << nl;
+
+    lumpedPointTools::setInterpolators(mesh, points0);
+
+
+    // Output vtk file series
+    vtk::seriesWriter stateSeries;
+    vtk::seriesWriter geomSeries;
+
+    // Initial geometry
+    movement.writeVTP("geom_init.vtp", state0, mesh, points0);
+
+    lumpedPointTools::setInterpolators(mesh);
+
+    for
+    (
+        label timei = 0, outputCount = 0;
+        timei < responseTable.size();
+        timei += span
+    )
+    {
+        lumpedPointState state = responseTable[timei].second();
+
+        state += movement.origin();
+        movement.scalePoints(state);
+        state.relax(relax, state0);
+
+        Info<< "output [" << timei << '/' << responseTable.size() << ']';
+
+        // State/response = what comes back from FEM
+        {
+            const word outputName =
+                word::printf("state_%06d.vtp", outputCount);
+
+            Info<< "  " << outputName;
+
+            movement.writeStateVTP(state, outputName);
+            stateSeries.append(outputCount, outputName);
+        }
+
+        {
+            const word outputName =
+                word::printf("geom_%06d.vtp", outputCount);
+
+            Info<< "  " << outputName;
+
+            movement.writeVTP(outputName, state, mesh, points0);
+            geomSeries.append(outputCount, outputName);
+        }
+
+        Info<< endl;
+
+        ++outputCount;
+
+        if (maxOut && outputCount >= maxOut)
+        {
+            Info<< "Max output " << maxOut << " ... stopping" << endl;
+            break;
         }
     }
 
-    Info<< args.executable() << ": finishing" << nl;
+
+    // Write file series
+
+    if (geomSeries.size())
+    {
+        Info<< nl << "write geom.vtp.series" << nl;
+        geomSeries.write("geom.vtp");
+    }
+    if (stateSeries.size())
+    {
+        Info<< nl << "write state.vtp.series" << nl;
+        stateSeries.write("state.vtp");
+    }
 
     Info<< "\nEnd\n" << endl;
 
diff --git a/applications/utilities/postProcessing/lumped/lumpedPointZones/Make/options b/applications/utilities/postProcessing/lumped/lumpedPointZones/Make/options
index 5821a660d0..df8ec7d4f0 100644
--- a/applications/utilities/postProcessing/lumped/lumpedPointZones/Make/options
+++ b/applications/utilities/postProcessing/lumped/lumpedPointZones/Make/options
@@ -2,9 +2,12 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/lumpedPointMotion/lnInclude
+    -I$(LIB_SRC)/lumpedPointMotion/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfileFormats \
     -lmeshTools \
-    -llumpedPointMotion
+    -llumpedPointMotion \
+    -ldynamicMesh
diff --git a/applications/utilities/postProcessing/lumped/lumpedPointZones/lumpedPointZones.C b/applications/utilities/postProcessing/lumped/lumpedPointZones/lumpedPointZones.C
index 86793e3054..d0edb92330 100644
--- a/applications/utilities/postProcessing/lumped/lumpedPointZones/lumpedPointZones.C
+++ b/applications/utilities/postProcessing/lumped/lumpedPointZones/lumpedPointZones.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +39,7 @@ Description
 
 #include "lumpedPointTools.H"
 #include "lumpedPointIOMovement.H"
+#include "fvMesh.H"
 
 using namespace Foam;
 
@@ -52,9 +53,26 @@ int main(int argc, char *argv[])
         " pressure integration zones used by lumpedPoint BC."
     );
 
-    argList::noParallel();    // The VTP writer is not yet in parallel
     argList::noFunctionObjects();  // Never use function objects
 
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Test initial lumped points state without a mesh"
+    );
+    argList::addOption
+    (
+        "visual-length",
+        "len",
+        "Visualization length for planes (visualized as triangles)"
+    );
+
+    argList::addBoolOption
+    (
+        "no-interpolate",
+        "Suppress calculation/display of point interpolators"
+    );
+
     argList::addBoolOption
     (
         "verbose",
@@ -62,47 +80,89 @@ int main(int argc, char *argv[])
     );
 
     #include "addRegionOption.H"
+
     #include "setRootCase.H"
 
+    const bool noInterpolate = args.found("no-interpolate");
+
+    const bool dryrun = args.found("dry-run");
+
     // const bool verbose = args.found("verbose");
 
+    args.readIfPresent("visual-length", lumpedPointState::visLength);
+
     #include "createTime.H"
 
-    runTime.setTime(instant(0, runTime.constant()), 0);
+    if (dryrun)
+    {
+        // Create without a mesh
+        autoPtr<lumpedPointIOMovement> movement =
+            lumpedPointIOMovement::New(runTime);
 
-    #include "createNamedPolyMesh.H"
+        if (!movement.valid())
+        {
+            Info<< "No valid movement found" << endl;
+            return 1;
+        }
 
-    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New
-    (
-        runTime
-    );
+        const word outputName("state.vtp");
+
+        Info<< "dry-run: writing " << outputName << nl;
+
+        movement().writeStateVTP(movement().state0(), outputName);
+
+        Info<< "\nEnd\n" << endl;
+
+        return 0;
+    }
+
+
+    runTime.setTime(instant(runTime.constant()), 0);
+
+    #include "createNamedMesh.H"
+
+    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New(mesh);
 
     if (!movement.valid())
     {
-        Info<< "no valid movement found" << endl;
+        Info<< "No valid movement found" << endl;
         return 1;
     }
 
-    const labelList patchLst = lumpedPointTools::lumpedPointPatchList(mesh);
-    if (patchLst.empty())
+    // Initial positions/rotations
+    movement().writeStateVTP("state.vtp");
+
+    pointIOField points0(lumpedPointTools::points0Field(mesh));
+
+    const label nPatches = lumpedPointTools::setPatchControls(mesh, points0);
+    if (!nPatches)
     {
-        Info<< "no patch list found" << endl;
+        Info<< "No point patches with lumped movement found" << endl;
         return 2;
     }
 
-    pointIOField points0 = lumpedPointTools::points0Field(mesh);
-    movement().setMapping(mesh, patchLst, points0);
+    Info<<"Lumped point patch controls set on "
+        << nPatches << " patches" << nl;
 
-    // Initial geometry, but with zone colouring
-    movement().writeZonesVTP("lumpedPointZones.vtp", mesh, points0);
+    Info<<"Areas per point: " << flatOutput(movement().areas(mesh)) << nl;
 
-    // Initial positions/rotations
-    movement().writeStateVTP("initialState.vtp");
+    if (noInterpolate)
+    {
+        // Initial geometry, with zones
+        movement().writeZonesVTP("lumpedPointZones.vtp", mesh, points0);
+    }
+    else
+    {
+        lumpedPointTools::setInterpolators(mesh, points0);
+
+        // Initial geometry, with zones and interpolations
+        movement().writeVTP("lumpedPointZones.vtp", mesh, points0);
+    }
 
     Info<< nl
+        << "wrote 'state.vtp' (reference state)" << nl
         << "wrote 'lumpedPointZones.vtp'" << nl
-        << "wrote 'initialState.vtp'" << nl
-        << "End\n" << endl;
+        << "\nEnd\n" << endl;
 
     return 0;
 }
diff --git a/src/lumpedPointMotion/Make/files b/src/lumpedPointMotion/Make/files
index 3125d35a6f..7548fe1e6b 100644
--- a/src/lumpedPointMotion/Make/files
+++ b/src/lumpedPointMotion/Make/files
@@ -1,3 +1,5 @@
+controller/lumpedPointController.C
+
 state/lumpedPointState.C
 state/lumpedPointStateWriter.C
 
diff --git a/src/lumpedPointMotion/controller/lumpedPointController.C b/src/lumpedPointMotion/controller/lumpedPointController.C
new file mode 100644
index 0000000000..8b7bc592d3
--- /dev/null
+++ b/src/lumpedPointMotion/controller/lumpedPointController.C
@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "lumpedPointController.H"
+#include "dictionary.H"
+#include "labelField.H"
+#include "Map.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::lumpedPointController::lumpedPointController() noexcept
+:
+    pointLabels_()
+{}
+
+
+Foam::lumpedPointController::lumpedPointController
+(
+    const labelUList& pointLabels
+)
+:
+    pointLabels_(pointLabels)
+{}
+
+
+Foam::lumpedPointController::lumpedPointController
+(
+    labelList&& pointLabels
+)
+:
+    pointLabels_(std::move(pointLabels))
+{}
+
+
+Foam::lumpedPointController::lumpedPointController
+(
+    const dictionary& dict
+)
+:
+    pointLabels_(dict.get<labelList>("pointLabels"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::lumpedPointController::remapPointLabels
+(
+    const label nPoints,
+    const Map<label>& originalIds
+)
+{
+    if (originalIds.size())
+    {
+        for (label& pointi : pointLabels_)
+        {
+            pointi = originalIds[pointi];
+        }
+    }
+
+    if (min(pointLabels_) < 0 || max(pointLabels_) >= nPoints)
+    {
+        FatalErrorInFunction
+            << "Point id out-of-range: " << flatOutput(pointLabels_) << nl
+            << exit(FatalError);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lumpedPointMotion/controller/lumpedPointController.H b/src/lumpedPointMotion/controller/lumpedPointController.H
index 7cec1b4541..f9aae10623 100644
--- a/src/lumpedPointMotion/controller/lumpedPointController.H
+++ b/src/lumpedPointMotion/controller/lumpedPointController.H
@@ -1 +1,139 @@
-#warning File removed - left for old dependency check only
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::lumpedPointController
+
+Description
+    Simple connectivity of point labels to specify a controller for lumped
+    point movement.
+
+    \heading Dictionary parameters
+    \table
+        Property    | Description                           | Required | Default
+        pointLabels | List of point labels                  | yes |
+    \endtable
+
+Note
+    If the calling program itself specified a point-label mapping
+    (eg, original ids from FEA), these can be used initially and remapped.
+
+SourceFiles
+    lumpedPointController.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lumpedPointController_H
+#define lumpedPointController_H
+
+#include "List.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class dictionary;
+template<class T> class Map;
+
+/*---------------------------------------------------------------------------*\
+                    Class lumpedPointController Declaration
+\*---------------------------------------------------------------------------*/
+
+class lumpedPointController
+{
+    // Private Data
+
+        //- The lumped points associated with the controller
+        labelList pointLabels_;
+
+
+public:
+
+    // Constructors
+
+        //- Default construct
+        lumpedPointController() noexcept;
+
+        //- Copy construct from point ids
+        explicit lumpedPointController(const labelUList& pointLabels);
+
+        //- Move construct from point ids
+        explicit lumpedPointController(labelList&& pointLabels);
+
+        //- Construct from dictionary
+        explicit lumpedPointController(const dictionary& dict);
+
+
+        //- Factory method
+        static autoPtr<lumpedPointController> New(const dictionary& dict)
+        {
+            return autoPtr<lumpedPointController>::New(dict);
+        }
+
+
+    //- Destructor
+    ~lumpedPointController() = default;
+
+
+    // Member Functions
+
+        //- No controller points specified?
+        bool empty() const
+        {
+            return pointLabels_.empty();
+        }
+
+        //- Number of controller points specified
+        label size() const
+        {
+            return pointLabels_.size();
+        }
+
+        //- The controller points
+        const labelList& pointLabels() const
+        {
+            return pointLabels_;
+        }
+
+        //- Check point mapping (the count) or remap the point labels
+        void remapPointLabels
+        (
+            const label nPoints,
+            const Map<label>& originalIds
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lumpedPointMotion/interpolation/lumpedPointInterpolator.H b/src/lumpedPointMotion/interpolation/lumpedPointInterpolator.H
index 7cec1b4541..5da13742fc 100644
--- a/src/lumpedPointMotion/interpolation/lumpedPointInterpolator.H
+++ b/src/lumpedPointMotion/interpolation/lumpedPointInterpolator.H
@@ -1 +1,160 @@
-#warning File removed - left for old dependency check only
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::lumpedPointInterpolator
+
+Description
+    A simple linear interpolator between two locations, which are
+    referenced by index.
+
+    When the interpolator is built for based on searching for nearest and
+    next-nearest points, the interpolation will typically cover a (0-0.5) range
+    rather than a (0-1) range. The (0.5-1) range implies a flip in the
+    relationship of nearest vs. next-nearest.
+
+SourceFiles
+    lumpedPointInterpolatorI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lumpedPointInterpolator_H
+#define lumpedPointInterpolator_H
+
+#include "Pair.H"
+#include "triFace.H"
+#include "scalarList.H"
+#include "barycentric2D.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class lumpedPointInterpolator Declaration
+\*---------------------------------------------------------------------------*/
+
+class lumpedPointInterpolator
+{
+    // Private Data
+
+        //- The id of the nearest point
+        label nearest_;
+
+        //- The id of the neighbour point(s)
+        label next1_, next2_;
+
+        //- The interpolation weight for the neighbour point(s)
+        scalar weight1_, weight2_;
+
+
+public:
+
+    // Constructors
+
+        //- Default construct, with zero weighting and invalid ids
+        inline lumpedPointInterpolator();
+
+        //- Construct with nearest id
+        explicit inline lumpedPointInterpolator(const label id);
+
+
+    //- Destructor
+    ~lumpedPointInterpolator() = default;
+
+
+    // Member Functions
+
+    // Access
+
+        //- Valid if there is an associated nearest point
+        inline bool valid() const;
+
+        //- The nearest control point, or -1 if invalid
+        inline label nearest() const;
+
+        //- The first neighbour control point - identical to next1()
+        inline label next() const;
+
+        //- The first neighbour control point, or -1 if invalid
+        inline label next1() const;
+
+        //- The second neighbour control point, or -1 if invalid
+        inline label next2() const;
+
+        //- The weighting for the nearest point
+        inline scalar weight0() const;
+
+        //- The weighting for the first neighbour point,
+        //- this also corresponds to the logical location (interval 0-1)
+        inline scalar weight1() const;
+
+        //- The weighting for the second neighbour point,
+        //- this also corresponds to the logical location (interval 0-1)
+        inline scalar weight2() const;
+
+
+    // Edit
+
+        //- Assign the nearest point, clearing any neighbour
+        inline void nearest(const label id);
+
+        //- Assign single neighbour control point and its weight
+        inline void next(const label id, const scalar weight);
+
+        //- Assign the neighbour control point and its weight
+        inline void next
+        (
+            const label id,
+            const scalar weight,
+            const label position    //!< Use 0 or 1 for first/second
+        );
+
+        //- Assign all control points and their weights
+        //  The triFace points [0,1,2] correspond to [nearest,next1,next2],
+        //  respectively
+        inline void set(const triFace& ids, const barycentric2D& weights);
+
+
+    // Evalulate
+
+        //- Linear interpolated value between nearest and next locations
+        template<class T>
+        inline T interpolate(const UList<T>& input) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "lumpedPointInterpolatorI.H"
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lumpedPointMotion/interpolation/lumpedPointInterpolatorI.H b/src/lumpedPointMotion/interpolation/lumpedPointInterpolatorI.H
new file mode 100644
index 0000000000..8556694aa3
--- /dev/null
+++ b/src/lumpedPointMotion/interpolation/lumpedPointInterpolatorI.H
@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::lumpedPointInterpolator::lumpedPointInterpolator()
+:
+    lumpedPointInterpolator(-1)
+{}
+
+
+inline Foam::lumpedPointInterpolator::lumpedPointInterpolator(const label id)
+:
+    nearest_(id),
+    next1_(id),
+    next2_(id),
+    weight1_(Zero),
+    weight2_(Zero)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::lumpedPointInterpolator::valid() const
+{
+    return nearest_ != -1;
+}
+
+
+inline Foam::label Foam::lumpedPointInterpolator::nearest() const
+{
+    return nearest_;
+}
+
+
+inline Foam::label Foam::lumpedPointInterpolator::next() const
+{
+    return next1_;
+}
+
+
+inline Foam::label Foam::lumpedPointInterpolator::next1() const
+{
+    return next1_;
+}
+
+
+inline Foam::label Foam::lumpedPointInterpolator::next2() const
+{
+    return next2_;
+}
+
+
+inline Foam::scalar Foam::lumpedPointInterpolator::weight0() const
+{
+    return scalar(1) - weight1_ - weight2_;
+}
+
+
+inline Foam::scalar Foam::lumpedPointInterpolator::weight1() const
+{
+    return weight1_;
+}
+
+
+inline Foam::scalar Foam::lumpedPointInterpolator::weight2() const
+{
+    return weight2_;
+}
+
+
+inline void Foam::lumpedPointInterpolator::nearest(const label id)
+{
+    nearest_ = id;
+    next(id, Zero);
+}
+
+
+inline void Foam::lumpedPointInterpolator::next
+(
+    const label id,
+    scalar weight
+)
+{
+    if (weight < scalar(0))
+    {
+        weight = 0;
+    }
+    else if (weight > scalar(1))
+    {
+        weight = 1;
+    }
+
+    next1_ = id;
+    next2_ = id;
+    weight1_ = weight;
+    weight2_ = Zero;
+}
+
+
+inline void Foam::lumpedPointInterpolator::next
+(
+    const label id,
+    scalar weight,
+    const label position
+)
+{
+    if (weight < scalar(0))
+    {
+        weight = 0;
+    }
+    else if (weight > scalar(1))
+    {
+        weight = 1;
+    }
+
+    if (!position)
+    {
+        next1_ = id;
+        weight1_ = weight;
+    }
+    else
+    {
+        next2_ = id;
+        weight2_ = weight;
+    }
+}
+
+
+inline void Foam::lumpedPointInterpolator::set
+(
+    const triFace& ids,
+    const barycentric2D& weights
+)
+{
+    nearest_ = ids[0];
+    next1_ = ids[1];
+    next2_ = ids[2];
+
+    weight1_ = weights[1];
+    weight2_ = weights[2];
+}
+
+
+template<class T>
+inline T Foam::lumpedPointInterpolator::interpolate(const UList<T>& input) const
+{
+    if (nearest_ == -1)
+    {
+        return Zero;
+    }
+    else if (next1_ == -1 || next1_ == nearest_)
+    {
+        return input[nearest_];
+    }
+
+    return
+    (
+        input[nearest_] * (1-weight1_-weight2_)
+      + input[next1_] * (weight1_)
+      + input[next2_] * (weight2_)
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/lumpedPointMotion/lumpedPointMovementDict b/src/lumpedPointMotion/lumpedPointMovementDict
index 4b21551430..5feadb91e7 100644
--- a/src/lumpedPointMotion/lumpedPointMovementDict
+++ b/src/lumpedPointMotion/lumpedPointMovementDict
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v1912                                 |
+|  \\    /   O peration     | Version:  v2006                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -15,36 +15,59 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Reference axis for the locations
-axis            (0 0 1);
+//- The initial locations of the lumped points
+points
+(
 
-// Locations of the lumped points
-locations       11(0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5);
+);
 
-// Division for pressure forces (0-1)
-division        0.5;
+//- The FEA point labels for the points (optional)
+// pointLabels
+// (
+//
+// );
 
-//- If present, the offset of patch points compared to the locations
-//  Otherwise determined from the bounding box
-// centre          (0 0 0);
+//- The motion controllers.
+//  The pointLabels are an adjacency list of point connectivity.
+//
+//  Controllers are joined by common point labels between lists
+//  (or internally - eg, ring structure) when they appear in the corresponding
+//  point displacement patch
+controllers
+{
+    control1
+    {
+        pointLabels (0 1 2 3 4);
+    }
+    control2
+    {
+        pointLabels (1 5 6);
+    }
+}
 
-//- The interpolation scheme
-interpolationScheme linear;
 
-//- Relaxation/scaling factor when updating positions
-relax           1.0;
+//- Input offset to shift ALL lumped points when reading (default: 0 0 0)
+// origin  (0 0 0);
 
+//- The Euler rotation order (default: zxz)
+// rotationOrder zxz;
+
+//- Input rotations in degrees instead of radians (default: false)
+// degrees false;
+
+//- Relaxation/scaling factor when updating positions (default: 1)
+// relax   1.0;
 
 forces
 {
-    //- The pressure name (default: p)
-    p           p;
+    //- Pressure name (default: p)
+    p       p;
 
     //- Reference pressure [Pa] (default: 0)
-    pRef        0;
+    pRef    0;
 
     //- Reference density for incompressible calculations (default: 1)
-    rhoRef      1;
+    rhoRef  1;
 }
 
 
@@ -60,6 +83,9 @@ communication
 
     initByExternal  false;
 
+    // Coupling frequency in time-steps (default: 1-steps)
+    // calcFrequency   1;
+
     // Input file of positions/rotation, written by external application
     inputName       positions.in;
 
@@ -70,26 +96,28 @@ communication
     logName         movement.log;
 
     inputFormat     dictionary;
+
     outputFormat    dictionary;
 
-    // Scaling applied to values read from 'inputName'
+    // Scaling applied to values read from 'inputName',
+    // also applies to points (above)
     scaleInput
     {
         //- Length multiplier (to metres). Eg 0.001 for [mm] -> [m]
-        length      1;
+        length  1;
     }
 
     // Scaling applied to values written to 'outputName'
     scaleOutput
     {
         //- Length multiplier (from metres). Eg 1000 for [m] -> [mm]
-        length      1;
+        length  1;
 
         //- Force units multiplier (from Pa)
-        force       1;
+        force   1;
 
         //- Moment units multiplier (from N.m)
-        moment      1;
+        moment  1;
     }
 }
 
diff --git a/src/lumpedPointMotion/movement/lumpedPointIOMovement.C b/src/lumpedPointMotion/movement/lumpedPointIOMovement.C
index 4de29c33a1..786d987394 100644
--- a/src/lumpedPointMotion/movement/lumpedPointIOMovement.C
+++ b/src/lumpedPointMotion/movement/lumpedPointIOMovement.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,16 +32,16 @@ License
 
 namespace Foam
 {
-    defineTypeNameAndDebug(lumpedPointIOMovement, 0);
+    defineTypeName(lumpedPointIOMovement);
 }
 
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-const Foam::lumpedPointIOMovement*
-Foam::lumpedPointIOMovement::lookupInRegistry(const objectRegistry& obr)
+Foam::lumpedPointIOMovement*
+Foam::lumpedPointIOMovement::getMovementObject(const objectRegistry& obr)
 {
-    return obr.findObject<lumpedPointIOMovement>
+    return obr.getObjectPtr<lumpedPointIOMovement>
     (
         lumpedPointMovement::canonicalName
     );
@@ -120,7 +120,7 @@ bool Foam::lumpedPointIOMovement::writeData(Ostream& os) const
 }
 
 
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const lumpedPointIOMovement& obj)
 {
diff --git a/src/lumpedPointMotion/movement/lumpedPointIOMovement.H b/src/lumpedPointMotion/movement/lumpedPointIOMovement.H
index d05f4d6ab0..7945b591f3 100644
--- a/src/lumpedPointMotion/movement/lumpedPointIOMovement.H
+++ b/src/lumpedPointMotion/movement/lumpedPointIOMovement.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,11 +46,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
-class lumpedPointIOMovement;
-
-Ostream& operator<<(Ostream& os, const lumpedPointIOMovement& obj);
-
 /*---------------------------------------------------------------------------*\
                     Class lumpedPointIOMovement Declaration
 \*---------------------------------------------------------------------------*/
@@ -68,22 +63,23 @@ class lumpedPointIOMovement
         //- No copy assignment
         void operator=(const lumpedPointIOMovement&) = delete;
 
+
 public:
 
     //- Runtime type information
-    TypeName("lumpedPointMovement");
+    TypeNameNoDebug("lumpedPointMovement");
 
 
     // Static Member Functions
 
-        //- Lookup pointer to object in the object-registry,
-        //  return nullptr if found.
-        static const lumpedPointIOMovement* lookupInRegistry
+        //- Find the movement object or nullptr if not found
+        static lumpedPointIOMovement* getMovementObject
         (
             const objectRegistry& obr
         );
 
-        //- Create a new object in the registry by reading system dictionary
+        //- Create a movement object in the registry by
+        //- reading system dictionary
         static autoPtr<lumpedPointIOMovement> New
         (
             const objectRegistry& obr,
@@ -94,36 +90,27 @@ public:
     // Constructors
 
         //- Construct from IOobject, optionally with some owner information
-        explicit lumpedPointIOMovement
-        (
-            const IOobject& io,
-            label ownerId = -1
-        );
+        explicit lumpedPointIOMovement(const IOobject& io, label ownerId = -1);
 
 
     //- Destructor
-    ~lumpedPointIOMovement() = default;
+    virtual ~lumpedPointIOMovement() = default;
 
 
     // Member Functions
 
-        //- readData member function used by regIOobject
+        //- The readData member function used by regIOobject
         bool readData(Istream& is);
 
-        //- writeData member function required by regIOobject
+        //- The writeData member function required by regIOobject
         bool writeData(Ostream& os) const;
-
-
-    // IOstream Operators
-
-        friend Ostream& operator<<
-        (
-            Ostream& os,
-            const lumpedPointIOMovement& obj
-        );
 };
 
 
+// IOstream Operators
+Ostream& operator<<(Ostream& os, const lumpedPointIOMovement& obj);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/lumpedPointMotion/movement/lumpedPointMovement.C b/src/lumpedPointMotion/movement/lumpedPointMovement.C
index c5f485b946..6bf3f63d3f 100644
--- a/src/lumpedPointMotion/movement/lumpedPointMovement.C
+++ b/src/lumpedPointMotion/movement/lumpedPointMovement.C
@@ -27,18 +27,30 @@ License
 
 #include "lumpedPointMovement.H"
 #include "lumpedPointIOMovement.H"
-#include "demandDrivenData.H"
-#include "linearInterpolationWeights.H"
 #include "Fstream.H"
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "PtrMap.H"
-#include "faceZoneMesh.H"
+#include "triFace.H"
+#include "labelPair.H"
+#include "indexedOctree.H"
+#include "treeDataPoint.H"
+#include "pointIndexHit.H"
+#include "pointPatch.H"
 #include "PstreamReduceOps.H"
 
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
+int Foam::lumpedPointMovement::debug
+(
+    ::Foam::debug::debugSwitch("lumpedPointMovement", 0)
+);
+
+
+const Foam::word
+Foam::lumpedPointMovement::canonicalName("lumpedPointMovement");
+
+
 const Foam::Enum
 <
     Foam::lumpedPointMovement::outputFormatType
@@ -62,9 +74,6 @@ Foam::lumpedPointMovement::scalingNames
 });
 
 
-const Foam::word
-Foam::lumpedPointMovement::canonicalName("lumpedPointMovement");
-
 
 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
@@ -73,21 +82,9 @@ namespace Foam
 
 //! \cond fileScope
 //- Space-separated vector value (ASCII)
-static inline Ostream& putPlain(Ostream& os, const vector& val)
+static inline Ostream& putPlain(Ostream& os, const vector& v)
 {
-    os  << val.x() << ' ' << val.y() << ' ' << val.z();
-    return os;
-}
-
-
-//! \cond fileScope
-//- Space-separated vector value (ASCII)
-static inline Ostream& putTime(Ostream& os, const Time& t)
-{
-    os  <<"Time index=" << t.timeIndex()
-        << " value=" << t.timeOutputValue();
-
-    return os;
+    return os << v.x() << ' ' << v.y() << ' ' << v.z();
 }
 
 
@@ -119,61 +116,18 @@ static void writeList(Ostream& os, const string& header, const UList<T>& list)
 } // End namespace Foam
 
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::lumpedPointMovement::calcThresholds() const
-{
-    thresholdPtr_ = new scalarField(locations_);
-
-    scalarField& thresh = *thresholdPtr_;
-
-    for (label i=1; i < thresh.size(); ++i)
-    {
-        thresh[i-1] =
-        (
-            locations_[i-1]
-          + (division_ * (locations_[i] - locations_[i-1]))
-        );
-    }
-}
-
-
-Foam::label Foam::lumpedPointMovement::threshold(scalar pos) const
-{
-    const scalarField& thresh = this->thresholds();
-    const label n = thresh.size();
-
-    // Clamp above and below
-    //
-    // ? may need special treatment to clamp values below the first threshold ?
-    // ? special treatment when 'axis' is negative ?
-
-    for (label i=0; i < n; ++i)
-    {
-        if (pos < thresh[i])
-        {
-            return i;
-        }
-    }
-
-    // everything above goes into the last zone too
-    return n-1;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::lumpedPointMovement::lumpedPointMovement()
 :
-    axis_(0,0,1),
-    locations_(),
-    division_(0),
+    origin_(Zero),
+    state0_(),
+    state_(),
+    originalIds_(),
+    controllers_(),
+    patchControls_(),
     relax_(1),
-    interpolationScheme_(linearInterpolationWeights::typeName),
     ownerId_(-1),
-    boundBox_(),
-    centre_(),
-    autoCentre_(true),
     forcesDict_(),
     coupler_(),
     inputName_("positions.in"),
@@ -181,12 +135,10 @@ Foam::lumpedPointMovement::lumpedPointMovement()
     logName_("movement.log"),
     inputFormat_(lumpedPointState::inputFormatType::DICTIONARY),
     outputFormat_(outputFormatType::DICTIONARY),
-    scaleInput_(-1.0),
-    scaleOutput_(-1.0),
-    state0_(),
-    state_(),
-    thresholdPtr_(0),
-    interpolatorPtr_()
+    scaleInput_(-1),
+    scaleOutput_(-1),
+    calcFrequency_(1),
+    lastTrigger_(-1)
 {}
 
 
@@ -196,107 +148,139 @@ Foam::lumpedPointMovement::lumpedPointMovement
     label ownerId
 )
 :
-    axis_(0,0,1),
-    locations_(),
-    division_(0),
+    origin_(Zero),
+    state0_(),
+    state_(),
+    originalIds_(),
+    controllers_(),
+    patchControls_(),
     relax_(1),
-    interpolationScheme_
-    (
-        dict.getOrDefault<word>
-        (
-            "interpolationScheme",
-            linearInterpolationWeights::typeName
-        )
-    ),
     ownerId_(ownerId),
-    boundBox_(),
-    centre_(),
-    autoCentre_(true),
     forcesDict_(),
     coupler_(),
     inputName_("positions.in"),
     outputName_("forces.out"),
     logName_("movement.log"),
-    scaleInput_(-1.0),
-    scaleOutput_(-1.0),
-    state0_(),
-    state_(),
-    thresholdPtr_(0),
-    interpolatorPtr_()
+    scaleInput_(-1),
+    scaleOutput_(-1),
+    calcFrequency_(1),
+    lastTrigger_(-1)
 {
     readDict(dict);
 }
 
 
-// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::lumpedPointMovement::~lumpedPointMovement()
+bool Foam::lumpedPointMovement::couplingPending(const label timeIndex) const
 {
-    deleteDemandDrivenData(thresholdPtr_);
+    return (timeIndex >= lastTrigger_ + calcFrequency_);
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+void Foam::lumpedPointMovement::couplingCompleted(const label timeIndex) const
+{
+    lastTrigger_ = timeIndex;
+}
+
 
 void Foam::lumpedPointMovement::readDict(const dictionary& dict)
 {
-    // assume the worst
-    deleteDemandDrivenData(thresholdPtr_);
+    origin_ = dict.getOrDefault<point>("origin", Zero);
 
-    dict.readEntry("axis", axis_);
+    // Initial point locations (zero rotation)
+    auto tpoints0 = tmp<pointField>::New();
+    auto& points0 = tpoints0.ref();
 
-    division_ = 0;
-    if (dict.readIfPresent("division", division_))
+    dict.readEntry("points", points0);
+    points0 += origin_;
+
+    originalIds_.clear();
+    controllers_.clear();
+    patchControls_.clear();
+
+
+    // The FEA ids for the points (optional)
+    Map<label> pointIdMap;
+
+    if (dict.readIfPresent("pointLabels", originalIds_) && originalIds_.size())
     {
-        if (division_ < 0)
+        if (originalIds_.size() != points0.size())
         {
-            division_ = 0;
+            FatalIOErrorInFunction(dict)
+                << "Incorrect number of pointLabels. Had "
+                << originalIds_.size() << " for " << points0.size() << " points"
+                << nl << endl
+                << exit(FatalIOError);
         }
-        else if (division_ > 1)
+
+        labelHashSet duplicates;
+
+        label pointi = 0;
+
+        for (const label id : originalIds_)
         {
-            division_ = 1;
+            if (!pointIdMap.insert(id, pointi))
+            {
+                duplicates.set(id);
+            }
+
+            ++pointi;
+        }
+
+        if (!duplicates.empty())
+        {
+            FatalIOErrorInFunction(dict)
+                << "Found duplicate point ids "
+                << flatOutput(duplicates.sortedToc()) << nl << endl
+                << exit(FatalIOError);
         }
     }
 
-    dict.readIfPresent("relax", relax_);
-
-    dict.readEntry("locations", locations_);
-
-    if (dict.readIfPresent("interpolationScheme", interpolationScheme_))
+    const dictionary* dictptr = dict.findDict("controllers");
+    if (dictptr)
     {
-        interpolatorPtr_.clear();
+        controllers_ = HashPtrTable<lumpedPointController>(*dictptr);
+
+        // Verify the input
+        forAllIters(controllers_, iter)
+        {
+            (*iter)->remapPointLabels(points0.size(), pointIdMap);
+        }
+    }
+    else
+    {
+        // Add single global controller
+        // Warning?
+
+        controllers_.clear();
     }
 
-    autoCentre_ = !dict.readIfPresent("centre", centre_);
-
-    // Initial state from initial locations, with zero rotation
-    tmp<vectorField> pts = locations_ * axis_;
-    state0_ = lumpedPointState(pts);
-    state_  = state0_;
+    relax_ = dict.getOrDefault<scalar>("relax", 1);
+    scalarMinMax::zero_one().inplaceClip(relax_);
 
     forcesDict_.merge(dict.subOrEmptyDict("forces"));
 
     const dictionary& commDict = dict.subDict("communication");
     coupler_.readDict(commDict);
 
-    // TODO: calcFrequency_  = dict.getOrDefault("calcFrequency", 1);
+    calcFrequency_ = commDict.getOrDefault<label>("calcFrequency", 1);
+
+    // Leave trigger intact
 
     commDict.readEntry("inputName", inputName_);
     commDict.readEntry("outputName", outputName_);
     commDict.readIfPresent("logName", logName_);
 
-    inputFormat_ =
-        lumpedPointState::formatNames.get("inputFormat", commDict);
-
-    outputFormat_ =
-        lumpedPointMovement::formatNames.get("outputFormat", commDict);
+    inputFormat_ = lumpedPointState::formatNames.get("inputFormat", commDict);
+    outputFormat_ = formatNames.get("outputFormat", commDict);
 
     scaleInput_  = -1;
     scaleOutput_ = -1;
 
     const dictionary* scaleDict = nullptr;
 
-    if ((scaleDict = commDict.findDict("scaleInput")))
+    if ((scaleDict = commDict.findDict("scaleInput")) != nullptr)
     {
         for (int i=0; i < scaleInput_.size(); ++i)
         {
@@ -314,7 +298,7 @@ void Foam::lumpedPointMovement::readDict(const dictionary& dict)
         }
     }
 
-    if ((scaleDict = commDict.findDict("scaleOutput")))
+    if ((scaleDict = commDict.findDict("scaleOutput")) != nullptr)
     {
         for (int i=0; i < scaleOutput_.size(); ++i)
         {
@@ -331,128 +315,578 @@ void Foam::lumpedPointMovement::readDict(const dictionary& dict)
             }
         }
     }
+
+    state0_ = lumpedPointState
+    (
+        tpoints0,
+        quaternion::eulerOrderNames.getOrDefault
+        (
+            "rotationOrder",
+            dict,
+            quaternion::eulerOrder::ZXZ
+        ),
+        dict.getOrDefault("degrees", false)
+    );
+
+    state0_.scalePoints(scaleInput_[scalingType::LENGTH]);
+
+    state_ = state0_;
 }
 
 
-void Foam::lumpedPointMovement::setBoundBox
+bool Foam::lumpedPointMovement::hasPatchControl
 (
-    const polyMesh& mesh,
-    const labelUList& patchIds,
-    const pointField& points0
-)
+    const polyPatch& pp
+) const
 {
-    boundBox_ = boundBox::invertedBox;
+    return hasPatchControl(pp.index());
+}
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    for (const label patchId : patchIds)
+
+bool Foam::lumpedPointMovement::hasInterpolator
+(
+    const pointPatch& fpatch
+) const
+{
+    return hasInterpolator(fpatch.index());
+}
+
+
+void Foam::lumpedPointMovement::checkPatchControl
+(
+    const polyPatch& pp
+) const
+{
+    const auto ctrlIter = patchControls_.cfind(pp.index());
+
+    if (!ctrlIter.good())
     {
-        boundBox_.add(points0, patches[patchId].meshPoints());
+        FatalErrorInFunction
+            << "No controllers for patch " << pp.name()
+            << exit(FatalError);
     }
 
-    boundBox_.reduce();
+    const patchControl& ctrl = *ctrlIter;
 
-    if (autoCentre_)
+    for (const word& ctrlName : ctrl.names_)
     {
-        centre_ = boundBox_.centre();
-        centre_ -= (centre_ & axis_) * axis_;
-        if (lumpedPointIOMovement::debug)
+        const auto iter = controllers_.cfind(ctrlName);
+
+        if (!iter.good())
         {
-            Pout<<"autoCentre on " << centre_
-                << " boundBox: " << boundBox_ << endl;
-        }
-    }
-    else
-    {
-        // User-specified centre
-        if (lumpedPointIOMovement::debug)
-        {
-            Pout<<"centre on " << centre_
-                << " boundBox: " << boundBox_ << endl;
+            FatalErrorInFunction
+                << "No controller: " << ctrlName << nl
+                << " For patch " << pp.name()
+                << exit(FatalError);
         }
     }
 }
 
 
-void Foam::lumpedPointMovement::setMapping
+void Foam::lumpedPointMovement::setPatchControl
 (
-    const polyMesh& mesh,
-    const labelUList& patchIds,
+    const polyPatch& pp,
+    const wordList& ctrlNames,
     const pointField& points0
 )
 {
-    setBoundBox(mesh, patchIds, points0);
+    // Reference mass centres
+    const pointField& lumpedCentres0 = state0().points();
 
-    faceZones_.clear();
+    const label patchIndex = pp.index();
 
-    // Local storage location (boundary faces only)
-    const label firstFace = mesh.nInternalFaces();
-    labelList faceToZoneID(mesh.nFaces() - firstFace, -1);
+    // Info<<"Add patch control for patch " << patchIndex << " "
+    //     << pp.name() << nl;
 
-    // Number of faces per zone
-    labelList nFaces(thresholds().size(), Zero);
+    patchControl& ctrl = patchControls_(patchIndex);
+    ctrl.names_ = ctrlNames;
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    labelList& faceToPoint = ctrl.faceToPoint_;
+    faceToPoint.resize(pp.size(), -1);
 
-    for (const label patchId : patchIds)
+    checkPatchControl(pp);
+
+    const faceList& faces = pp.boundaryMesh().mesh().faces();
+
+    // Subset of points to search (if specified)
+    labelHashSet subsetPointIds;
+
+    for (const word& ctrlName : ctrl.names_)
     {
-        const polyPatch& pp = patches[patchId];
+        const auto iter = controllers_.cfind(ctrlName);
 
-        // Classify faces into respective zones according to their centres
-        label meshFaceI = pp.start();
-
-        forAll(pp, i)
+        if (!iter.good())
         {
-            const face& f = mesh.faces()[meshFaceI];
-            label zoneId = this->threshold(f.centre(points0));
-
-            // localize storage location
-            faceToZoneID[meshFaceI - firstFace] = zoneId;
-            ++nFaces[zoneId];
-
-            ++meshFaceI;
+            FatalErrorInFunction
+                << "No controller: " << ctrlName << nl
+                << exit(FatalError);
         }
+
+        const labelList& pointLabels = (*iter)->pointLabels();
+
+        subsetPointIds.insert(pointLabels);
     }
 
-    if (lumpedPointIOMovement::debug)
+    if (ctrl.names_.size() && subsetPointIds.empty())
     {
-        Pout<<"faces per zone:" << nFaces << endl;
+        FatalErrorInFunction
+            << "Controllers specified, but without any points" << nl
+            << exit(FatalError);
     }
 
-    faceZones_.setSize(nFaces.size());
 
-    label nZonedFaces = 0;
-    forAll(nFaces, zonei)
+    treeDataPoint treePoints
+    (
+        lumpedCentres0,
+        subsetPointIds.sortedToc(),
+        subsetPointIds.size()
+    );
+
+
+    treeBoundBox bb(lumpedCentres0);
+    bb.inflate(0.01);
+
+    indexedOctree<treeDataPoint> ppTree
+    (
+        treePoints,
+        bb,     // overall search domain
+        8,      // maxLevel
+        10,     // leafsize
+        3.0     // duplicity
+    );
+
+    const scalar searchDistSqr(sqr(GREAT));
+
+    const label patchStart = pp.start();
+    forAll(pp, patchFacei)
     {
-        label n = nFaces[zonei];
-        labelList& addr = faceZones_[zonei];
-        addr.setSize(n);
+        const point fc(faces[patchStart + patchFacei].centre(points0));
 
-        n = 0;
-        forAll(faceToZoneID, faceI)
+        // Store the original pointId, not subset id
+        faceToPoint[patchFacei] =
+            treePoints.pointLabel
+            (
+                ppTree.findNearest(fc, searchDistSqr).index()
+            );
+    }
+
+    if (debug)
+    {
+        Pout<<"Added face mapping for patch: " << patchIndex << endl;
+    }
+}
+
+
+void Foam::lumpedPointMovement::setInterpolator
+(
+    const pointPatch& fpatch,
+    const pointField& points0
+)
+{
+    // Reference mass centres
+    const pointField& lumpedCentres0 = state0().points();
+
+    const label patchIndex = fpatch.index();
+
+    patchControl& ctrl = patchControls_(patchIndex);
+
+    List<lumpedPointInterpolator>& interpList = ctrl.interp_;
+    interpList.clear();
+
+    // The connectivity, adjacency list
+    Map<labelHashSet> adjacency;
+
+    // Subset of points to search (if specified)
+    labelHashSet subsetPointIds;
+
+    for (const word& ctrlName : ctrl.names_)
+    {
+        const auto iter = controllers_.cfind(ctrlName);
+
+        if (!iter.good())
         {
-            if (faceToZoneID[faceI] == zonei)
+            FatalErrorInFunction
+                << "No controller: " << ctrlName << nl
+                << exit(FatalError);
+        }
+
+        const labelList& pointLabels = (*iter)->pointLabels();
+
+        subsetPointIds.insert(pointLabels);
+
+        // Adjacency lists
+        forAll(pointLabels, i)
+        {
+            const label curr = pointLabels[i];
+
+            if (i)
             {
-                // from local storage to global mesh face
-                label meshFaceI = faceI + firstFace;
+                const label prev = pointLabels[i-1];
+                adjacency(prev).insert(curr);
+                adjacency(curr).insert(prev);
+            }
+            else if (!adjacency.found(curr))
+            {
+                adjacency(curr).clear();
+            }
+        }
+    }
 
-                addr[n] = meshFaceI;
-                ++n;
+    if (ctrl.names_.empty())
+    {
+        // Adjacency lists
+        const label len = state0().size();
+
+        for (label i=0; i < len; ++i)
+        {
+            const label curr = i;
+
+            if (i)
+            {
+                const label prev = i-1;
+                adjacency(prev).insert(curr);
+                adjacency(curr).insert(prev);
+            }
+            else if (!adjacency.found(curr))
+            {
+                adjacency(curr).clear();
+            }
+        }
+    }
+
+    if (ctrl.names_.size() && subsetPointIds.empty())
+    {
+        FatalErrorInFunction
+            << "Controllers specified, but without any points" << nl
+            << exit(FatalError);
+    }
+
+
+    // Pairs defining connecting points as triangles
+    Map<labelPairList> adjacencyPairs;
+
+    barycentric2D bary;
+
+    {
+        // Pairs for the ends
+        DynamicList<labelPair> pairs;
+
+        // Mag sin(angle) around the triangle point 0,
+        // used to sort the generated triangles according to the acute angle
+        DynamicList<scalar> acuteAngles;
+
+        forAllConstIters(adjacency, iter)
+        {
+            const label nearest = iter.key();
+
+            labelList neighbours(iter.val().sortedToc());
+
+            const label len = neighbours.size();
+
+            pairs.clear();
+            acuteAngles.clear();
+
+            const point& nearPt = lumpedCentres0[nearest];
+
+            for (label j=1; j < len; ++j)
+            {
+                for (label i=0; i < j; ++i)
+                {
+                    labelPair neiPair(neighbours[i], neighbours[j]);
+
+                    triPointRef tri
+                    (
+                        nearPt,
+                        lumpedCentres0[neiPair.first()],
+                        lumpedCentres0[neiPair.second()]
+                    );
+
+                    // Require non-degenerate triangles
+                    if (tri.pointToBarycentric(tri.a(), bary) > SMALL)
+                    {
+                        // Triangle OK
+                        pairs.append(neiPair);
+
+                        vector ab(normalised(tri.b() - tri.a()));
+                        vector ac(normalised(tri.c() - tri.a()));
+
+                        // Angle between neighbouring edges
+                        // Use negative cosine to map [0-180] -> [-1 .. +1]
+
+                        acuteAngles.append(-(ab & ac));
+                    }
+                }
+            }
+
+            if (pairs.size() > 1)
+            {
+                // Sort by acute angle
+                labelList order(sortedOrder(acuteAngles));
+                inplaceReorder(order, pairs);
+            }
+
+            adjacencyPairs.insert(nearest, pairs);
+        }
+    }
+
+    if (debug & 2)
+    {
+        Info<< "Adjacency table for patch: " << fpatch.name() << nl;
+
+        for (const label own : adjacency.sortedToc())
+        {
+            Info<< own << " =>";
+            for (const label nei : adjacency[own].sortedToc())
+            {
+                Info<< ' ' << nei;
+            }
+
+            if (originalIds_.size())
+            {
+                Info<< "  # " << originalIds_[own] << " =>";
+                for (const label nei : adjacency[own].sortedToc())
+                {
+                    Info<< ' ' << originalIds_[nei];
+                }
+            }
+
+            Info<< "  # tri " << flatOutput(adjacencyPairs[own]);
+            Info<< nl;
+        }
+    }
+
+    treeDataPoint treePoints
+    (
+        lumpedCentres0,
+        subsetPointIds.sortedToc(),
+        subsetPointIds.size()
+    );
+
+    treeBoundBox bb(lumpedCentres0);
+    bb.inflate(0.01);
+
+    indexedOctree<treeDataPoint> ppTree
+    (
+        treePoints,
+        bb,     // overall search domain
+        8,      // maxLevel
+        10,     // leafsize
+        3.0     // duplicity
+    );
+
+
+    // Searching
+
+    const scalar searchDistSqr(sqr(GREAT));
+
+    const labelList& meshPoints = fpatch.meshPoints();
+
+    interpList.resize(meshPoints.size());
+
+    DynamicList<scalar> unsortedNeiWeightDist;
+    DynamicList<label>  unsortedNeighbours;
+
+    forAll(meshPoints, pointi)
+    {
+        const point& ptOnMesh = points0[meshPoints[pointi]];
+
+        // Nearest (original) point id
+        const label nearest =
+            treePoints.pointLabel
+            (
+                ppTree.findNearest(ptOnMesh, searchDistSqr).index()
+            );
+
+        interpList[pointi].nearest(nearest);
+
+        if (nearest == -1)
+        {
+            // Should not really happen
+            continue;
+        }
+
+        // Have the nearest lumped point, now find the next nearest
+        // but check that the direction is also correct.
+
+        // OK: within the 0-1 bounds
+        // 1+----+0
+        //      .
+        //     .
+        //    +pt
+
+        const point& nearPt = lumpedCentres0[nearest];
+
+        const linePointRef toMeshPt(nearPt, ptOnMesh);
+
+        const labelPairList& adjacentPairs = adjacencyPairs[nearest];
+
+        unsortedNeighbours = adjacency[nearest].toc();
+        unsortedNeiWeightDist.resize(unsortedNeighbours.size());
+
+        forAll(unsortedNeighbours, nbri)
+        {
+            unsortedNeiWeightDist[nbri] =
+                magSqr(ptOnMesh - lumpedCentres0[unsortedNeighbours[nbri]]);
+        }
+
+        // Sort by distance
+        labelList distOrder(sortedOrder(unsortedNeiWeightDist));
+
+        label ngood = 0;
+
+        // Recalculate distance as weighting
+        for (const label nbri : distOrder)
+        {
+            const label nextPointi = unsortedNeighbours[nbri];
+
+            const point& nextPt = lumpedCentres0[nextPointi];
+
+            const linePointRef toNextPt(nearPt, nextPt);
+
+            const scalar weight =
+                (toMeshPt.vec() & toNextPt.unitVec()) / toNextPt.mag();
+
+            if (weight < 0)
+            {
+                // Reject: wrong direction or other bad value
+                continue;
+            }
+            else
+            {
+                // Store weight
+                unsortedNeiWeightDist[nbri] = weight;
+
+                // Retain good weight
+                distOrder[ngood] = nbri;
+                ++ngood;
             }
         }
 
-        addr.setSize(n);
-        nZonedFaces += n;
+        distOrder.resize(ngood);
 
-        if (lumpedPointIOMovement::debug)
+        if (distOrder.size() < 1)
         {
-            Pout<< "Created faceZone[" << zonei << "] "
-                << returnReduce(n, sumOp<label>()) << " faces, "
-                << thresholds()[zonei] << " threshold" << nl;
+            continue;
+        }
+
+        UIndirectList<label>  neighbours(unsortedNeighbours, distOrder);
+        UIndirectList<scalar> neiWeight(unsortedNeiWeightDist, distOrder);
+
+        bool useFirst = true;
+
+        if (neighbours.size() > 1 && adjacentPairs.size())
+        {
+            // Check for best two neighbours
+
+            bitSet neiPointid;
+            neiPointid.set(neighbours);
+
+            for (const labelPair& ends : adjacentPairs)
+            {
+                label nei1 = ends.first();
+                label nei2 = ends.second();
+
+                if (!neiPointid.test(nei1) || !neiPointid.test(nei2))
+                {
+                    // Reject, invalid combination for this point location
+                    continue;
+                }
+                else if (neighbours.find(nei2) < neighbours.find(nei1))
+                {
+                    // Order by distance, which is not really needed,
+                    // but helps with diagnostics
+                    Swap(nei1, nei2);
+                }
+
+
+                triFace triF(nearest, nei1, nei2);
+
+                if
+                (
+                    triF.tri(lumpedCentres0).pointToBarycentric(ptOnMesh, bary)
+                  > SMALL
+                 && !bary.outside()
+                )
+                {
+                    // Use barycentric weights
+                    interpList[pointi].set(triF, bary);
+
+                    useFirst = false;
+                    break;
+                }
+            }
+        }
+
+        if (useFirst)
+        {
+            // Use geometrically closest neighbour
+            interpList[pointi].next(neighbours.first(), neiWeight.first());
         }
     }
 }
 
 
+Foam::List<Foam::scalar> Foam::lumpedPointMovement::areas
+(
+    const polyMesh& pmesh
+) const
+{
+    const label nLumpedPoints = state0().size();
+
+    List<scalar> zoneAreas(nLumpedPoints, Zero);
+
+    if (patchControls_.empty())
+    {
+        WarningInFunction
+            << "Attempted to calculate areas without setMapping()"
+            << nl << endl;
+        return zoneAreas;
+    }
+
+    const polyBoundaryMesh& patches = pmesh.boundaryMesh();
+
+    // fvMesh and has pressure field
+    if (isA<fvMesh>(pmesh))
+    {
+        const fvMesh& mesh = dynamicCast<const fvMesh>(pmesh);
+
+        // Face areas (on patches)
+        const surfaceVectorField::Boundary& patchSf =
+            mesh.Sf().boundaryField();
+
+        forAllConstIters(patchControls_, iter)
+        {
+            const label patchIndex = iter.key();
+            const patchControl& ctrl = iter.val();
+
+            const labelList& faceToPoint = ctrl.faceToPoint_;
+
+            const polyPatch& pp = patches[patchIndex];
+
+            forAll(pp, patchFacei)
+            {
+                // Force from the patch-face into sum
+                const label pointIndex = faceToPoint[patchFacei];
+
+                if (pointIndex < 0)
+                {
+                    // Unmapped, for whatever reason?
+                    continue;
+                }
+
+                // Force from the patch-face into sum
+                zoneAreas[pointIndex] += mag(patchSf[patchIndex][patchFacei]);
+            }
+        }
+    }
+
+    Pstream::listCombineGather(zoneAreas, plusEqOp<scalar>());
+    Pstream::listCombineScatter(zoneAreas);
+
+    return zoneAreas;
+}
+
+
 bool Foam::lumpedPointMovement::forcesAndMoments
 (
     const polyMesh& pmesh,
@@ -460,17 +894,19 @@ bool Foam::lumpedPointMovement::forcesAndMoments
     List<vector>& moments
 ) const
 {
-    forces.setSize(faceZones_.size());
-    moments.setSize(faceZones_.size());
+    const label nLumpedPoints = state0().size();
 
-    if (faceZones_.empty())
+    forces.resize(nLumpedPoints);
+    moments.resize(nLumpedPoints);
+
+    if (patchControls_.empty())
     {
         WarningInFunction
             << "Attempted to calculate forces without setMapping()"
             << nl << endl;
 
-        forces.setSize(thresholds().size(), Zero);
-        moments.setSize(thresholds().size(), Zero);
+        forces.resize(nLumpedPoints, Zero);
+        moments.resize(nLumpedPoints, Zero);
         return false;
     }
 
@@ -478,7 +914,7 @@ bool Foam::lumpedPointMovement::forcesAndMoments
     forces = Zero;
     moments = Zero;
 
-    // The mass centres
+    // Current mass centres
     const pointField& lumpedCentres = state().points();
 
     const polyBoundaryMesh& patches = pmesh.boundaryMesh();
@@ -499,17 +935,14 @@ bool Foam::lumpedPointMovement::forcesAndMoments
         const fvMesh& mesh = dynamicCast<const fvMesh>(pmesh);
         const volScalarField& p = *pPtr;
 
-        // face centres (on patches)
-        const surfaceVectorField::Boundary& patchCf =
-            mesh.Cf().boundaryField();
+        // Face centres (on patches)
+        const surfaceVectorField::Boundary& patchCf = mesh.Cf().boundaryField();
 
-        // face areas (on patches)
-        const surfaceVectorField::Boundary& patchSf =
-            mesh.Sf().boundaryField();
+        // Face areas (on patches)
+        const surfaceVectorField::Boundary& patchSf = mesh.Sf().boundaryField();
 
         // Pressure (on patches)
-        const volScalarField::Boundary& patchPress =
-            p.boundaryField();
+        const volScalarField::Boundary& patchPress = p.boundaryField();
 
         // rhoRef if the pressure field is dynamic, i.e. p/rho otherwise 1
         rhoRef = (p.dimensions() == dimPressure ? 1.0 : rhoRef);
@@ -517,62 +950,63 @@ bool Foam::lumpedPointMovement::forcesAndMoments
         // Scale pRef by density for incompressible simulations
         pRef /= rhoRef;
 
-        forAll(faceZones_, zonei)
+        forAllConstIters(patchControls_, iter)
         {
-            const labelList& zn = faceZones_[zonei];
+            const label patchIndex = iter.key();
+            const patchControl& ctrl = iter.val();
 
-            forAll(zn, localFaceI)
+            const labelList& faceToPoint = ctrl.faceToPoint_;
+
+            if (!forceOnPatches.found(patchIndex))
             {
-                const label meshFaceI = zn[localFaceI];
+                // Patch faces are +ve outwards,
+                // so the forces (exerted by fluid on solid)
+                // already have the correct sign
+                forceOnPatches.set
+                (
+                    patchIndex,
+                    (
+                        rhoRef
+                      * patchSf[patchIndex] * (patchPress[patchIndex] - pRef)
+                    ).ptr()
+                );
+            }
 
-                // locate which patch this face is on,
-                // and its offset within the patch
-                const label patchI = patches.whichPatch(meshFaceI);
-                if (patchI == -1)
+            const vectorField& forceOnPatch = *forceOnPatches[patchIndex];
+
+            const polyPatch& pp = patches[patchIndex];
+
+            forAll(pp, patchFacei)
+            {
+                // Force from the patch-face into sum
+                const label pointIndex = faceToPoint[patchFacei];
+
+                if (pointIndex < 0)
                 {
-                    // Should not happen - could also warn
+                    // Unmapped, for whatever reason?
                     continue;
                 }
 
-                const label patchFaceI = meshFaceI - patches[patchI].start();
-
-                if (!forceOnPatches.found(patchI))
-                {
-                    // Patch faces are +ve outwards,
-                    // so the forces (exerted by fluid on solid)
-                    // already have the correct sign
-                    forceOnPatches.set
-                    (
-                        patchI,
-                        (
-                            rhoRef
-                          * patchSf[patchI] * (patchPress[patchI] - pRef)
-                        ).ptr()
-                    );
-                }
-
-                const vectorField& forceOnPatch = *forceOnPatches[patchI];
-
                 // Force from the patch-face into sum
-                forces[zonei] += forceOnPatch[patchFaceI];
+                forces[pointIndex] += forceOnPatch[patchFacei];
 
-                // effective torque arm:
+                // Effective torque arm:
                 // - translated into the lumped-points coordinate system
                 //   prior to determining the distance
-                const vector lever =
+                const vector lever
                 (
-                    patchCf[patchI][patchFaceI] - centre_
-                  - lumpedCentres[zonei]
+                    patchCf[patchIndex][patchFacei]
+                  - lumpedCentres[pointIndex]
                 );
 
                 // Moment from the patch-face into sum
-                moments[zonei] += lever ^ forceOnPatch[patchFaceI];
+                moments[pointIndex] += lever ^ forceOnPatch[patchFacei];
             }
         }
     }
     else
     {
-        Info<<"no pressure" << endl;
+        Info<<"No pressure field" << endl;
     }
 
     Pstream::listCombineGather(forces, plusEqOp<vector>());
@@ -585,90 +1019,105 @@ bool Foam::lumpedPointMovement::forcesAndMoments
 }
 
 
-// \cond fileScope
-template<class T>
-static inline T interpolatedValue
-(
-    const Foam::UList<T> input,
-    const Foam::labelUList& indices,
-    const Foam::scalarField& weights
-)
-{
-    T result = weights[0] * input[indices[0]];
-    for (Foam::label i=1; i < indices.size(); ++i)
-    {
-        result += weights[i] * input[indices[i]];
-    }
-
-    return result;
-}
-// \endcond
-
-
 Foam::tmp<Foam::pointField>
-Foam::lumpedPointMovement::displacePoints
+Foam::lumpedPointMovement::pointsDisplacement
 (
-    const pointField& points0,
-    const labelList& pointLabels
+    const pointPatch& fpatch,
+    const pointField& points0
 ) const
 {
-    return displacePoints(state(), points0, pointLabels);
+    return pointsDisplacement(state(), fpatch, points0);
 }
 
 
 Foam::tmp<Foam::pointField>
-Foam::lumpedPointMovement::displacePoints
+Foam::lumpedPointMovement::pointsDisplacement
 (
     const lumpedPointState& state,
-    const pointField& points0,
-    const labelList& pointLabels
+    const pointPatch& fpatch,
+    const pointField& points0
 ) const
 {
-    // storage for the interpolation indices/weights
-    labelList   indices;
-    scalarField weights;
-    const interpolationWeights& interp = interpolator();
+    const label patchIndex = fpatch.index();
 
-    // The mass centres
+    // Reference mass centres
+    const pointField& lumpedCentres0 = state0().points();
+
+    // Current mass centres
     const pointField& lumpedCentres = state.points();
 
-    // local-to-global transformation tensor
+    // The local-to-global transformation tensor
     const tensorField& localToGlobal = state.rotations();
 
-    // could also verify the sizes (state vs original)
+    const labelList& meshPoints = fpatch.meshPoints();
 
-    tmp<pointField> tdisp(new pointField(pointLabels.size()));
-    pointField& disp = tdisp.ref();
+    // Could also verify the sizes (state vs original)
 
-    forAll(pointLabels, ptI)
+    auto tdisp = tmp<pointField>::New(fpatch.size());
+    auto& disp = tdisp.ref();
+
+    // The interpolator
+    const List<lumpedPointInterpolator>& interpList
+        = patchControls_[patchIndex].interp_;
+
+    forAll(meshPoints, pointi)
     {
-        const point& p0 = points0[pointLabels[ptI]];
+        const lumpedPointInterpolator& interp = interpList[pointi];
 
-        // Interpolation factor based on the axis component
-        // of the original point (location)
-        scalar where = p0 & axis_;
+        const point& p0 = points0[meshPoints[pointi]];
 
-        interp.valueWeights(where, indices, weights);
+        const vector origin0 = interp.interpolate(lumpedCentres0);
+        const vector origin = interp.interpolate(lumpedCentres);
+        const tensor rotTensor = interp.interpolate(localToGlobal);
 
-        vector origin    = interpolatedValue(lumpedCentres, indices, weights);
-        tensor rotTensor = interpolatedValue(localToGlobal, indices, weights);
+        disp[pointi] = (rotTensor & (p0 - origin0)) + origin - p0;
+    }
 
-        if (indices.size() == 1)
-        {
-            // out-of-bounds, use corresponding end-point
-            where = locations_[indices[0]];
-        }
+    return tdisp;
+}
 
-        // local point:
-        // - translated into the lumped-points coordinate system
-        // - relative to the interpolation (where) location
-        const vector local = (p0 - (where * axis_)) - centre_;
 
-        // local-to-global rotation and translation
-        // - translate back out of the lumped-points coordinate system
-        //
-        // -> subtract p0 to return displacements
-        disp[ptI] = (rotTensor & local) + origin + centre_ - p0;
+Foam::tmp<Foam::pointField>
+Foam::lumpedPointMovement::pointsPosition
+(
+    const lumpedPointState& state,
+    const pointPatch& fpatch,
+    const pointField& points0
+) const
+{
+    const label patchIndex = fpatch.index();
+
+    // Reference mass centres
+    const pointField& lumpedCentres0 = state0().points();
+
+    // Current mass centres
+    const pointField& lumpedCentres = state.points();
+
+    // The local-to-global transformation tensor
+    const tensorField& localToGlobal = state.rotations();
+
+    const labelList& meshPoints = fpatch.meshPoints();
+
+    // Could also verify the sizes (state vs original)
+
+    auto tdisp = tmp<pointField>::New(fpatch.size());
+    auto& disp = tdisp.ref();
+
+    // The interpolator
+    const List<lumpedPointInterpolator>& interpList =
+        patchControls_[patchIndex].interp_;
+
+    forAll(meshPoints, pointi)
+    {
+        const lumpedPointInterpolator& interp = interpList[pointi];
+
+        const point& p0 = points0[meshPoints[pointi]];
+
+        const vector origin0 = interp.interpolate(lumpedCentres0);
+        const vector origin = interp.interpolate(lumpedCentres);
+        const tensor rotTensor = interp.interpolate(localToGlobal);
+
+        disp[pointi] = (rotTensor & (p0 - origin0)) + origin;
     }
 
     return tdisp;
@@ -677,10 +1126,8 @@ Foam::lumpedPointMovement::displacePoints
 
 void Foam::lumpedPointMovement::writeDict(Ostream& os) const
 {
-    os.writeEntry("axis", axis_);
-    os.writeEntry("locations", locations_);
-    os.writeEntry("division", division_);
-    // os.writeEntry("interpolationScheme", interpolationScheme_);
+    // os.writeEntry("axis", axis_);
+    // os.writeEntry("locations", locations_);
 }
 
 
@@ -688,13 +1135,15 @@ bool Foam::lumpedPointMovement::readState()
 {
     lumpedPointState prev = state_;
 
-    bool status = state_.readData
+    const bool status = state_.readData
     (
         inputFormat_,
-        coupler().resolveFile(inputName_)
+        coupler().resolveFile(inputName_),
+        state0().rotationOrder(),
+        state0().degrees()
     );
 
-    state_.scalePoints(scaleInput_[scalingType::LENGTH]);
+    scalePoints(state_);
 
     state_.relax(relax_, prev);
 
@@ -718,8 +1167,10 @@ bool Foam::lumpedPointMovement::writeData
         os  <<"########" << nl;
         if (timeinfo)
         {
-            os <<"# ";
-            putTime(os, *timeinfo) << nl;
+            const Time& t = *timeinfo;
+
+            os  <<"# Time index=" << t.timeIndex() << nl
+                <<"# Time value=" << t.timeOutputValue() << nl;
         }
         os  <<"# size=" << this->size() << nl
             <<"# columns (points) (forces)";
@@ -781,11 +1232,11 @@ bool Foam::lumpedPointMovement::writeData
 
         os <<"########" << nl;
 
-        forAll(locations_, i)
+        forAll(state0().points(), i)
         {
-            const vector pos = scaleLength * (locations_[i] * axis_);
+            const point& pt = state0().points()[i];
 
-            putPlain(os, pos) << ' ';
+            putPlain(os, scaleLength * pt) << ' ';
 
             if (i < forces.size())
             {
@@ -823,12 +1274,14 @@ bool Foam::lumpedPointMovement::writeData
         os  <<"////////" << nl;
         if (timeinfo)
         {
-            os <<"// ";
-            putTime(os, *timeinfo) << nl;
+            const Time& t = *timeinfo;
+
+            os  <<"// Time index=" << t.timeIndex() << nl;
+            os.writeEntry("time", t.timeOutputValue());
         }
         os  << nl;
 
-        writeList(os, "points", (locations_*axis_)());
+        writeList(os, "points", state0().points());
         writeList(os, "forces", forces);
 
         if (writeMoments)
@@ -881,20 +1334,4 @@ bool Foam::lumpedPointMovement::writeData
 }
 
 
-const Foam::interpolationWeights&
-Foam::lumpedPointMovement::interpolator() const
-{
-    if (!interpolatorPtr_.valid())
-    {
-        interpolatorPtr_ = interpolationWeights::New
-        (
-            interpolationScheme_,
-            locations_
-        );
-    }
-
-    return *interpolatorPtr_;
-}
-
-
 // ************************************************************************* //
diff --git a/src/lumpedPointMotion/movement/lumpedPointMovement.H b/src/lumpedPointMotion/movement/lumpedPointMovement.H
index f5cfd71c78..adbc403f05 100644
--- a/src/lumpedPointMotion/movement/lumpedPointMovement.H
+++ b/src/lumpedPointMotion/movement/lumpedPointMovement.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,9 +27,9 @@ Class
     Foam::lumpedPointMovement
 
 Description
-    The movement \a driver that describes initial point locations, the
-    segmentation for pressure integration, the current state of the
-    points/rotations, and forwarding to the externalFileCoupler
+    The movement \a driver that describes initial point locations,
+    the current state of the points/rotations,
+    and forwarding to the externalFileCoupler
     communication coordinator.
 
     The lumpedPointIOMovement class is simply a registered version of
@@ -41,14 +41,15 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                           | Required | Default
-        axis        | Reference axis for the locations      | yes |
-        locations   | List of lumped point locations        | yes |
-        centre      | Offset of patch points to locations   | no | automatic
-        division    | Division (0-1) for pressure forces    | no | 0
-        relax       | Relaxation/scaling for updating positions | no  |
-        interpolationScheme | The interpolation scheme      | yes |
-        forces      | Optional forces dictionary            | no |
-        communication | Required communication dictionary   | yes |
+        points      | Initial locations of lumped points    | yes |
+        pointLabels | The FEA ids for the points            | no  |
+        origin      | Shift offset when reading points      | no  | (0 0 0)
+        rotationOrder | The Euler rotation order            | no  | zxz
+        degrees     | Input rotations in degrees            | no  | false
+        relax       | Relaxation/scaling for updating positions | no  | 1
+        controllers | Motion controllers (dictionary)       | yes |
+        forces      | Force settings (dictionary)           | no  |
+        communication | Communication settings (dictionary) | yes |
     \endtable
 
     \heading Parameters for communication dictionary
@@ -60,7 +61,8 @@ Description
         inputFormat | Input format: dictionary/plain    | yes |
         outputFormat | Output format: dictionary/plain  | yes |
         scaleInput  | Input scaling parameter dictionary | no |
-        scaleOutput | Output scaling parameter dictionary | no |
+        scaleOutput | Output scaling parameter dictionary | no  |
+        calcFrequency | Calculation/coupling frequency    | no  | 1
     \endtable
 
     \heading Parameters for optional communication/scaleInput dictionary
@@ -95,28 +97,25 @@ SourceFiles
 
 #include "dictionary.H"
 #include "scalarList.H"
-#include "scalarField.H"
-#include "pointField.H"
-#include "vectorField.H"
-#include "tensorField.H"
-#include "vector.H"
-#include "interpolationWeights.H"
+#include "primitiveFields.H"
 #include "IOobject.H"
 #include "tmp.H"
-#include "faceZoneMeshFwd.H"
+#include "HashPtrTable.H"
 #include "externalFileCoupler.H"
+#include "lumpedPointController.H"
+#include "lumpedPointInterpolator.H"
 #include "lumpedPointState.H"
-#include "boundBox.H"
 #include "Enum.H"
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class polyMesh;
+class polyPatch;
+class pointPatch;
 class Time;
 
 /*---------------------------------------------------------------------------*\
@@ -127,118 +126,116 @@ class lumpedPointMovement
 {
 public:
 
-    //- Output format types
-    enum class outputFormatType
-    {
-        PLAIN,      //!< "plain" is a simple ASCII format
-        DICTIONARY  //!< "dictionary" is the OpenFOAM dictionary format
-    };
+    // Data Types
 
-    //- Output format types
-    enum scalingType
-    {
-        LENGTH = 0, //!< The "length" scaling
-        FORCE,      //!< The "force" scaling
-        MOMENT      //!< The "moment" scaling
-    };
+        //- Output format types
+        enum class outputFormatType
+        {
+            PLAIN,      //!< "plain" is a simple ASCII format
+            DICTIONARY  //!< "dictionary" is the OpenFOAM dictionary format
+        };
 
-    // Static data
+        //- Output format types
+        enum scalingType
+        {
+            LENGTH = 0, //!< The "length" scaling
+            FORCE,      //!< The "force" scaling
+            MOMENT      //!< The "moment" scaling
+        };
 
-    //- Names for the output format types
-    static const Enum<outputFormatType> formatNames;
 
-    //- Names for the scaling types
-    static const Enum<scalingType> scalingNames;
+    // Static Data
+
+        //- Names for the output format types
+        static const Enum<outputFormatType> formatNames;
+
+        //- Names for the scaling types
+        static const Enum<scalingType> scalingNames;
 
 
 private:
 
-    // Private data
+    // Private Class
 
-        //- Reference axis for the locations
-        vector axis_;
+        //- The controller names and faceCentre mapping for a patch
+        struct patchControl
+        {
+            wordList names_;
+            labelList faceToPoint_;
+            List<lumpedPointInterpolator> interp_;
+        };
 
-        //- The locations of the lumped points within the reference axis
-        //  The interpolator needs scalarField, not scalarList.
-        scalarField locations_;
 
-        //- The division when calculating pressure forces (0-1)
-        scalar division_;
+    // Private Data
+
+        //- The offset for lumped points, used on input.
+        point origin_;
+
+        //- The initial state of positions with zero rotations.
+        lumpedPointState state0_;
+
+        //- The current state (positions, rotations)
+        //  Eg, as response from external application
+        lumpedPointState state_;
+
+        //- The original point ids (optional information)
+        labelList originalIds_;
+
+        //- Connectivity for the controllers
+        HashPtrTable<lumpedPointController> controllers_;
+
+        //- The controller names and faceCentre mapping for patches
+        Map<patchControl> patchControls_;
 
         //- The relaxation factor when moving points (default: 1)
         scalar relax_;
 
-        //- The interpolation type (linear|spline)
-        word interpolationScheme_;
-
-        //- Optional owner information
+        //- Optional owner information (patch owner)
         label ownerId_;
 
-        //- The bounding box (if set)
-        boundBox boundBox_;
-
-        //- The offset of patch points to compared to the locations
-        point centre_;
-
-        //- Calculate centre based on the bounding box
-        bool autoCentre_;
-
         //- Dictionary of controls for force calculation
         dictionary forcesDict_;
 
         //- Communication control
         externalFileCoupler coupler_;
 
-        //- File io
+        //- File IO names
         word inputName_;
         word outputName_;
         word logName_;
 
+        //- The input format for points, rotations
         lumpedPointState::inputFormatType inputFormat_;
+
+        //- The output format for forces, moments
         outputFormatType outputFormat_;
 
-        //- Optional scale factors for input/output files
+        //- Scale factors for input/output files (optional)
         FixedList<scalar, 1> scaleInput_;
         FixedList<scalar, 3> scaleOutput_;
 
+        //- Calculation frequency
+        label calcFrequency_;
 
-    // Demand-driven private data
-
-        //- The initial state (positions, rotations)
-        lumpedPointState state0_;
-
-        //- The current state (positions, rotations)
-        //  Eg, as response from external application
-        lumpedPointState state_;
-
-        //- Threshold locations for pressure forces
-        mutable scalarField* thresholdPtr_;
-
-        //- User-specified interpolator
-        mutable autoPtr<interpolationWeights> interpolatorPtr_;
-
-        //- Pressure zones (only used from the master patch)
-        mutable List<labelList> faceZones_;
+        //- The last timeIndex when coupling was triggered
+        mutable label lastTrigger_;
 
 
     // Private Member Functions
 
-        //- Calculate threshold locations
-        void calcThresholds() const;
-
-        //- Classify the position to be located in one of the threshold zones
-        label threshold(scalar pos) const;
-
-
         //- No copy construct
         lumpedPointMovement(const lumpedPointMovement&) = delete;
 
         //- No copy assignment
         void operator=(const lumpedPointMovement&) = delete;
 
+
 public:
 
-    // Static data members
+    // Static Data Members
+
+        //- Debug switch
+        static int debug;
 
         //- The canonical name ("lumpedPointMovement") for the dictionary
         static const word canonicalName;
@@ -246,15 +243,15 @@ public:
 
     // Constructors
 
-        //- Construct null
+        //- Default construct
         lumpedPointMovement();
 
         //- Construct from dictionary, optionally with some owner information
-        lumpedPointMovement(const dictionary& dict, label ownerId=-1);
+        explicit lumpedPointMovement(const dictionary& dict, label ownerId=-1);
 
 
     //- Destructor
-    virtual ~lumpedPointMovement();
+    virtual ~lumpedPointMovement() = default;
 
 
     // Member Functions
@@ -269,41 +266,36 @@ public:
         //- The number of lumped points (number of locations)
         inline label size() const;
 
-        //- The normalized reference axis
-        inline const vector& axis() const;
-
-        //- Read access to the locations
-        inline const scalarField& locations() const;
-
-        //- The division (0-1) when calculating pressure forces
-        inline scalar division() const;
-
         //- An owner Id, if needed for bookkeeping purposes
         inline label ownerId() const;
 
         //- Change the owner id, if needed for bookkeeping purposes
         inline void ownerId(label id);
 
-
-        //- Threshold locations for pressure forces
-        inline const scalarField& thresholds() const;
-
-        //- Classify the position to be located in one of the threshold zones
-        inline label threshold(const point& position) const;
-
-
         //- Communication control
         inline const externalFileCoupler& coupler() const;
 
         //- Communication control
         inline externalFileCoupler& coupler();
 
+        //- Check if coupling is pending (according to the calcFrequency)
+        bool couplingPending(const label timeIndex) const;
+
+        //- Register that coupling is completed at this calcFrequency
+        void couplingCompleted(const label timeIndex) const;
+
         //- The initial state (positions/rotations)
         inline const lumpedPointState& state0() const;
 
         //- The current state (positions/rotations)
         inline const lumpedPointState& state() const;
 
+        //- The offset for lumped points, used on input.
+        inline const point& origin() const;
+
+        //- Scale the lumped points (on input).
+        inline void scalePoints(lumpedPointState& state) const;
+
         //- The relaxation factor when changing states
         inline scalar relax() const;
 
@@ -325,25 +317,40 @@ public:
         //- The output (forces) file format
         inline lumpedPointMovement::outputFormatType outputFormat() const;
 
+        //- The Euler-angle rotation order
+        inline quaternion::eulerOrder rotationOrder() const;
 
-        //- Define the bounding-box required to enclose the specified patches.
-        //  Calculates the centre as required.
+        //-  Rotation angles in degrees
+        inline bool degrees() const;
+
+        //- Check if patch control exists for specified patch
+        inline bool hasPatchControl(const label patchIndex) const;
+
+        //- Check if patch control exists for specified patch
+        bool hasPatchControl(const polyPatch& pp) const;
+
+        //- Check if patch control exists for specified patch
+        bool hasInterpolator(const pointPatch& fpatch) const;
+
+        //- Check if patch control exists for specified patch
+        void checkPatchControl(const polyPatch& pp) const;
+
+        //- Define pressure-zones mapping for faces in the specified patches.
+        //  The face centres are compared to the controller points,
         //
-        //  \param mesh     The volume mesh reference
-        //  \param patchIds The patch ids to be included in the bounding box.
+        //  \param pp       The patch with a control
+        //  \param ctrlNames The patch ids to be included in the mapping
         //  \param points0  The initial mesh points, prior to movement
-        void setBoundBox
+        void setPatchControl
         (
-            const polyMesh& mesh,
-            const labelUList& patchIds,
+            const polyPatch& pp,
+            const wordList& ctrlNames,
             const pointField& points0
         );
 
 
-        //- Define the pressure-zones mapping for faces in the specified
-        //  patches.
-        //  The face centres are compared to the threshold positions,
-        //  which are determined by locations along the defined axis.
+        //- Define pressure-zones mapping for faces in the specified patches.
+        //  The face centres are compared to the controller points,
         //
         //  \param mesh     The volume mesh reference
         //  \param patchIds The patch ids to be included in the mapping
@@ -356,13 +363,21 @@ public:
         );
 
 
+        //- Check if patch control exists for specified patch
+        void setInterpolator
+        (
+            const pointPatch& fpatch,
+            const pointField& points0
+        );
+
         //- True if the pressure-zones mapping has already been performed
         inline bool hasMapping() const;
 
-        //- Return the pressure-zones mapping with the associated
-        //  patch face ids.
-        inline const List<labelList>& zones() const;
+        //- Check if patch interpolator exists for specified patch
+        inline bool hasInterpolator(const label patchIndex) const;
 
+        //- The areas for each pressure-zone.
+        List<scalar> areas(const polyMesh& pmesh) const;
 
         //- The forces and moments acting on each pressure-zone.
         //  The zones must be previously defined via setMapping.
@@ -375,24 +390,28 @@ public:
 
 
         //- Displace points according to the current state
-        tmp<pointField> displacePoints
+        tmp<pointField> pointsDisplacement
         (
-            const pointField& points0,
-            const labelList& pointLabels
+            const pointPatch& fpatch,
+            const pointField& points0
         ) const;
 
-
         //- Displace points according to specified state
-        tmp<pointField> displacePoints
+        tmp<pointField> pointsDisplacement
         (
             const lumpedPointState& state,
-            const pointField& points0,
-            const labelList& pointLabels
+            const pointPatch& fpatch,
+            const pointField& points0
         ) const;
 
+        //- The points absolute position according to specified state
+        tmp<pointField> pointsPosition
+        (
+            const lumpedPointState& state,
+            const pointPatch& fpatch,
+            const pointField& points0
+        ) const;
 
-        //- Interpolation weights
-        const interpolationWeights& interpolator() const;
 
         //- Write axis, locations, division as a dictionary
         void writeDict(Ostream& os) const;
@@ -418,6 +437,13 @@ public:
         //- Read state from file, applying relaxation as requested
         bool readState();
 
+        //- Write state as VTK PolyData format.
+        void writeStateVTP
+        (
+            const lumpedPointState& state,
+            const fileName& file
+        ) const;
+
         //- Write state as VTK PolyData format.
         void writeStateVTP(const fileName& file) const;
 
@@ -438,14 +464,12 @@ public:
             const pointField& points0
         ) const;
 
-
         //- Write displaced geometry according to the current state,
         //  write as VTK PolyData format.
         void writeVTP
         (
             const fileName& file,
             const polyMesh& mesh,
-            const labelUList& patchIds,
             const pointField& points0
         ) const;
 
@@ -456,25 +480,8 @@ public:
             const fileName& file,
             const lumpedPointState& state,
             const polyMesh& mesh,
-            const labelUList& patchLst,
             const pointField& points0
         ) const;
-
-
-    // Housekeeping
-
-        //- Deprecated(2018-08) compatibility method
-        //  \deprecated(2018-08) - use IOobject::selectIO directly
-        static IOobject selectIO
-        (
-            const IOobject& io,
-            const fileName& altFile,
-            const word& ioName = ""
-        )
-        {
-            return IOobject::selectIO(io, altFile, ioName);
-        }
-
 };
 
 
diff --git a/src/lumpedPointMotion/movement/lumpedPointMovementI.H b/src/lumpedPointMotion/movement/lumpedPointMovementI.H
index ead2000b02..e1f8bc7abf 100644
--- a/src/lumpedPointMotion/movement/lumpedPointMovementI.H
+++ b/src/lumpedPointMotion/movement/lumpedPointMovementI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,31 +27,13 @@ License
 
 inline bool Foam::lumpedPointMovement::empty() const
 {
-    return locations_.empty();
+    return state0_.empty();
 }
 
 
 inline Foam::label Foam::lumpedPointMovement::size() const
 {
-    return locations_.size();
-}
-
-
-inline const Foam::vector& Foam::lumpedPointMovement::axis() const
-{
-    return axis_;
-}
-
-
-inline const Foam::scalarField& Foam::lumpedPointMovement::locations() const
-{
-    return locations_;
-}
-
-
-inline Foam::scalar Foam::lumpedPointMovement::division() const
-{
-    return division_;
+    return state0_.size();
 }
 
 
@@ -67,21 +49,19 @@ inline void Foam::lumpedPointMovement::ownerId(label id)
 }
 
 
-inline const Foam::scalarField& Foam::lumpedPointMovement::thresholds() const
+inline bool
+Foam::lumpedPointMovement::hasPatchControl(const label patchIndex) const
 {
-    if (!thresholdPtr_)
-    {
-        calcThresholds();
-    }
-
-    return *thresholdPtr_;
+    return patchControls_.found(patchIndex);
 }
 
 
-inline Foam::label
-Foam::lumpedPointMovement::threshold(const point& position) const
+inline bool
+Foam::lumpedPointMovement::hasInterpolator(const label patchIndex) const
 {
-    return threshold(position & axis_);
+    const auto iter = patchControls_.cfind(patchIndex);
+
+    return (iter.good() && iter().interp_.size());
 }
 
 
@@ -98,7 +78,6 @@ inline Foam::externalFileCoupler& Foam::lumpedPointMovement::coupler()
 }
 
 
-//- The initial state (positions/rotations)
 inline const Foam::lumpedPointState& Foam::lumpedPointMovement::state0() const
 {
     return state0_;
@@ -111,6 +90,21 @@ inline const Foam::lumpedPointState& Foam::lumpedPointMovement::state() const
 }
 
 
+inline const Foam::point& Foam::lumpedPointMovement::origin() const
+{
+    return origin_;
+}
+
+
+inline void Foam::lumpedPointMovement::scalePoints
+(
+    lumpedPointState& state
+) const
+{
+    state.scalePoints(scaleInput_[scalingType::LENGTH]);
+}
+
+
 inline Foam::scalar Foam::lumpedPointMovement::relax() const
 {
     return relax_;
@@ -155,16 +149,22 @@ Foam::lumpedPointMovement::outputFormat() const
 }
 
 
-inline bool Foam::lumpedPointMovement::hasMapping() const
+inline Foam::quaternion::eulerOrder
+Foam::lumpedPointMovement::rotationOrder() const
 {
-    return !faceZones_.empty();
+    return state0().rotationOrder();
 }
 
 
-inline const Foam::List<Foam::labelList>&
-Foam::lumpedPointMovement::zones() const
+inline bool Foam::lumpedPointMovement::degrees() const
 {
-    return faceZones_;
+    return state0().degrees();
+}
+
+
+inline bool Foam::lumpedPointMovement::hasMapping() const
+{
+    return !patchControls_.empty();
 }
 
 
diff --git a/src/lumpedPointMotion/movement/lumpedPointMovementWriter.C b/src/lumpedPointMotion/movement/lumpedPointMovementWriter.C
index 02e96a1b96..84cd47ada6 100644
--- a/src/lumpedPointMotion/movement/lumpedPointMovementWriter.C
+++ b/src/lumpedPointMotion/movement/lumpedPointMovementWriter.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,16 +27,54 @@ License
 
 #include "lumpedPointMovement.H"
 #include "polyMesh.H"
+#include "pointMesh.H"
 #include "OFstream.H"
-#include "uindirectPrimitivePatch.H"
-
 #include "foamVtkOutput.H"
+#include "foamVtkSurfaceWriter.H"
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+void Foam::lumpedPointMovement::writeStateVTP
+(
+    const lumpedPointState& state,
+    const fileName& file
+) const
+{
+    if (!Pstream::master())
+    {
+        // No extra information available from slaves, write on master only.
+        return;
+    }
+
+    labelListList lines;
+
+    label nLines = controllers_.size();
+
+    if (nLines)
+    {
+        lines.resize(nLines);
+        nLines = 0;
+
+        for (const word& ctrlName : controllers_.sortedToc())
+        {
+            lines[nLines] = controllers_[ctrlName]->pointLabels();
+            ++nLines;
+        }
+    }
+    else
+    {
+        // Default - global with all points as a single line
+        lines.resize(1);
+        lines.first() = identity(state.size());
+    }
+
+    state.writeVTP(file, lines, originalIds_);
+}
+
+
 void Foam::lumpedPointMovement::writeStateVTP(const fileName& file) const
 {
-    state().writeVTP(file, axis());
+    writeStateVTP(state(), file);
 }
 
 
@@ -47,6 +85,12 @@ void Foam::lumpedPointMovement::writeForcesAndMomentsVTP
     const UList<vector>& moments
 ) const
 {
+    if (!Pstream::master())
+    {
+        // Force, moments already reduced
+        return;
+    }
+
     OFstream fos(file);
     std::ostream& os = fos.stdStream();
 
@@ -60,7 +104,7 @@ void Foam::lumpedPointMovement::writeForcesAndMomentsVTP
         .beginVTKFile<vtk::fileTag::POLY_DATA>();
 
     //
-    // The 'spine' of lumped mass points
+    // The 'backbone' of lumped mass points
     //
     const label nPoints = state().points().size();
 
@@ -176,126 +220,38 @@ void Foam::lumpedPointMovement::writeZonesVTP
     const pointField& points0
 ) const
 {
-    OFstream fos(file);
-    std::ostream& os = fos.stdStream();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const labelList patchIds(patchControls_.sortedToc());
 
-    autoPtr<vtk::formatter> format = vtk::newFormatter
+    vtk::surfaceWriter writer
     (
-        os,
-        vtk::formatType::INLINE_ASCII
+        pointField::null(),
+        faceList::null(),
+        vtk::formatType::INLINE_ASCII,
+        file
     );
 
-    format().xmlHeader()
-        .beginVTKFile<vtk::fileTag::POLY_DATA>();
-
-    forAll(faceZones_, zoneI)
+    for (const label patchi : patchIds)
     {
-        uindirectPrimitivePatch pp
+        const labelList& faceToPoint = patchControls_[patchi].faceToPoint_;
+
+        primitivePatch pp
         (
-            UIndirectList<face>(mesh.faces(), faceZones_[zoneI]),
+            SubList<face>(mesh.faces(), patches[patchi].range()),
             points0
         );
 
-        format()
-            .tag
-            (
-                vtk::fileTag::PIECE,
-                vtk::fileAttr::NUMBER_OF_POINTS, pp.nPoints(),
-                vtk::fileAttr::NUMBER_OF_POLYS,  pp.size()
-            );
+        writer.piece(pp.localPoints(), pp.localFaces());
 
-        // 'points'
-        {
-            const uint64_t payLoad = vtk::sizeofData<float, 3>(pp.nPoints());
+        writer.writeGeometry();
 
-            format()
-                .tag(vtk::fileTag::POINTS)
-                .beginDataArray<float, 3>(vtk::dataArrayAttr::POINTS);
+        writer.beginCellData(2);
 
-            format().writeSize(payLoad);
-            vtk::writeList(format(), pp.localPoints());
-            format().flush();
+        writer.writeUniform("patchId", patchi);
+        writer.write("lumpedId", faceToPoint);
 
-            format()
-                .endDataArray()
-                .endTag(vtk::fileTag::POINTS);
-        }
-
-        // <Polys>
-        format().tag(vtk::fileTag::POLYS);
-
-        //
-        // 'connectivity'
-        //
-        {
-            label nVerts = 0;
-            for (const face& f : pp)
-            {
-                nVerts += f.size();
-            }
-
-            const uint64_t payLoad = vtk::sizeofData<label>(nVerts);
-
-            format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
-            format().writeSize(payLoad);
-
-            for (const face& f : pp.localFaces())
-            {
-                vtk::writeList(format(), f);
-            }
-            format().flush();
-
-            format().endDataArray();
-        }
-
-        //
-        // 'offsets'  (connectivity offsets)
-        //
-        {
-            const uint64_t payLoad = vtk::sizeofData<label>(pp.size());
-
-            format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
-            format().writeSize(payLoad);
-
-            label off = 0;
-            forAll(pp, facei)
-            {
-                const face& f = pp[facei];
-                off += f.size();
-
-                format().write(off);
-            }
-            format().flush();
-
-            format().endDataArray();
-        }
-
-        format().endTag(vtk::fileTag::POLYS);
-
-
-        format().beginCellData();
-
-        // zone Id
-        {
-            const uint64_t payLoad = vtk::sizeofData<label>(pp.size());
-
-            format().beginDataArray<label>("zoneId");
-            format().writeSize(payLoad);
-
-            vtk::write(format(), zoneI, pp.size());
-
-            format().flush();
-
-            format().endDataArray();
-        }
-
-        format().endCellData();
-
-        format().endPiece();
+        writer.endCellData();
     }
-
-    format().endTag(vtk::fileTag::POLY_DATA)
-        .endVTKFile();
 }
 
 
@@ -303,11 +259,10 @@ void Foam::lumpedPointMovement::writeVTP
 (
     const fileName& file,
     const polyMesh& mesh,
-    const labelUList& patchIds,
     const pointField& points0
 ) const
 {
-    writeVTP(file, state(), mesh, patchIds, points0);
+    writeVTP(file, state(), mesh, points0);
 }
 
 
@@ -316,119 +271,100 @@ void Foam::lumpedPointMovement::writeVTP
     const fileName& file,
     const lumpedPointState& state,
     const polyMesh& mesh,
-    const labelUList& patchIds,
     const pointField& points0
 ) const
 {
-    const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const labelList patchIds(patchControls_.sortedToc());
 
-    OFstream fos(file);
-    std::ostream& os = fos.stdStream();
+    pointMesh ptMesh(mesh);
 
-    autoPtr<vtk::formatter> format = vtk::newFormatter
+    vtk::surfaceWriter writer
     (
-        os,
-        vtk::formatType::INLINE_ASCII
+        pointField::null(),
+        faceList::null(),
+        vtk::formatType::INLINE_ASCII,
+        file
     );
 
-    format().xmlHeader()
-        .beginVTKFile<vtk::fileTag::POLY_DATA>();
-
-    for (const label patchId : patchIds)
+    for (const label patchi : patchIds)
     {
-        const polyPatch& pp = boundaryMesh[patchId];
+        const polyPatch& pp = patches[patchi];
 
-        format()
-            .tag
-            (
-                vtk::fileTag::PIECE,
-                vtk::fileAttr::NUMBER_OF_POINTS, pp.nPoints(),
-                vtk::fileAttr::NUMBER_OF_POLYS,  pp.size()
-            );
+        const pointPatch& ptPatch = ptMesh.boundary()[patchi];
 
-        // 'points'
+        // Current position (not displacement)
+        tmp<pointField> tpts = pointsPosition(state, ptPatch, points0);
+
+        writer.piece(tpts(), pp.localFaces());
+
+        writer.writeGeometry();
+
+        // Face mapping
+        const labelList& faceToPoint = patchControls_[patchi].faceToPoint_;
+
+        writer.beginCellData(2);
+
+        writer.writeUniform("patchId", patchi);
+        writer.write("lumpedId", faceToPoint);
+
+        writer.endCellData();
+
+        // The interpolator
+        const List<lumpedPointInterpolator>& interpList
+            = patchControls_[patchi].interp_;
+
+        writer.beginPointData(3);
+
+        // Nearest, Next
         {
-            const uint64_t payLoad = vtk::sizeofData<float, 3>(pp.nPoints());
+            labelList intData(interpList.size());
 
-            format()
-                .tag(vtk::fileTag::POINTS)
-                .beginDataArray<float, 3>(vtk::dataArrayAttr::POINTS);
+            forAll(interpList, i)
+            {
+                intData[i] = interpList[i].nearest();
+            }
+            writer.write("nearest", intData);
 
-            // Could be more efficient, but not often needed
-            tmp<pointField> tpts = displacePoints
-            (
-                state,
-                points0,
-                pp.meshPoints()
-            ) + pointField(points0, pp.meshPoints());
+            forAll(interpList, i)
+            {
+                intData[i] = interpList[i].next1();
+            }
+            writer.write("next1", intData);
 
-            const pointField& pts = tpts();
 
-            format().writeSize(payLoad);
-            vtk::writeList(format(), pts);
-            format().flush();
-
-            format()
-                .endDataArray()
-                .endTag(vtk::fileTag::POINTS);
+            forAll(interpList, i)
+            {
+                intData[i] = interpList[i].next2();
+            }
+            writer.write("next2", intData);
         }
 
-        // <Polys>
-        format().tag(vtk::fileTag::POLYS);
-
-
-        //
-        // 'connectivity'
-        //
+        // Weights
         {
-            label nVerts = 0;
-            for (const face& f : pp)
+            scalarList floatData(interpList.size());
+
+            forAll(interpList, i)
             {
-                nVerts += f.size();
+                floatData[i] = interpList[i].weight0();
             }
+            writer.write("weight", floatData);
 
-            const uint64_t payLoad = vtk::sizeofData<label>(nVerts);
-
-            format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
-            format().writeSize(payLoad);
-
-            for (const face& f : pp.localFaces())
+            forAll(interpList, i)
             {
-                vtk::writeList(format(), f);
+                floatData[i] = interpList[i].weight1();
             }
-            format().flush();
+            writer.write("weight1", floatData);
 
-            format().endDataArray();
+            forAll(interpList, i)
+            {
+                floatData[i] = interpList[i].weight2();
+            }
+            writer.write("weight2", floatData);
         }
 
-        //
-        // 'offsets'  (connectivity offsets)
-        //
-        {
-            const uint64_t payLoad = vtk::sizeofData<label>(pp.size());
-
-            format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
-            format().writeSize(payLoad);
-
-            label off = 0;
-            for (const face& f : pp)
-            {
-                off += f.size();
-
-                format().write(off);
-            }
-            format().flush();
-
-            format().endDataArray();
-        }
-
-        format().endTag(vtk::fileTag::POLYS);
-
-        format().endPiece();
+        writer.endPointData();
     }
-
-    format().endTag(vtk::fileTag::POLY_DATA)
-        .endVTKFile();
 }
 
 
diff --git a/src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.C b/src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.C
index e467b73fdd..a81a50e2a0 100644
--- a/src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.C
+++ b/src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,74 @@ namespace Foam
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+Foam::label Foam::lumpedPointDisplacementPointPatchVectorField::setPatchControls
+(
+    const pointVectorField& pvf,
+    const pointField& points0
+)
+{
+    label count = 0;
+
+    const pointVectorField::Boundary& bf = pvf.boundaryField();
+    const polyBoundaryMesh& patches = pvf.mesh().mesh().boundaryMesh();
+
+    forAll(bf, patchi)
+    {
+        // Patch of this type
+        const auto* p = isA<patchType>(bf[patchi]);
+
+        if (p)
+        {
+            // Patch controls (mapping) for calculating forces/moments
+            const_cast<lumpedPointMovement&>(p->movement())
+                .setPatchControl
+                (
+                    patches[patchi],
+                    p->controllers(),
+                    points0
+                );
+
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+Foam::label Foam::lumpedPointDisplacementPointPatchVectorField::setInterpolators
+(
+    const pointVectorField& pvf,
+    const pointField& points0
+)
+{
+    label count = 0;
+
+    const pointVectorField::Boundary& bf = pvf.boundaryField();
+
+    forAll(bf, patchi)
+    {
+        // Patch of this type
+        const auto* p = isA<patchType>(bf[patchi]);
+
+        if (p)
+        {
+            // Patch controls (mapping) for calculating forces/moments
+            const_cast<lumpedPointMovement&>(p->movement())
+                .setInterpolator
+                (
+                    p->patch(),
+                    points0
+                );
+
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
 Foam::labelList
 Foam::lumpedPointDisplacementPointPatchVectorField::patchIds
 (
@@ -61,7 +129,7 @@ Foam::lumpedPointDisplacementPointPatchVectorField::patchIds
     DynamicList<label> patchLst(bf.size());
     forAll(bf, patchi)
     {
-        // All patches of this type
+        // Patch of this type
         if (isA<patchType>(bf[patchi]))
         {
             patchLst.append(patchi);
@@ -80,11 +148,33 @@ Foam::lumpedPointDisplacementPointPatchVectorField::points0() const
 {
     const objectRegistry& obr = this->patch().boundaryMesh().mesh().db();
 
-    // Obtain starting locations from the motionSolver
-    return obr.lookupObject<displacementMotionSolver>
-    (
-        "dynamicMeshDict"
-    ).points0();
+    // Obtain starting locations from the motionSolver (when possible)
+    const auto* solver =
+        obr.cfindObject<displacementMotionSolver>("dynamicMeshDict");
+
+    if (solver)
+    {
+        if (points0Ptr_)
+        {
+            points0Ptr_.reset(nullptr);
+        }
+        return solver->points0();
+    }
+    else if (!points0Ptr_)
+    {
+        points0Ptr_.reset
+        (
+            new pointIOField
+            (
+                points0MotionSolver::points0IO
+                (
+                    this->patch().boundaryMesh().mesh().mesh()
+                )
+            )
+        );
+    }
+
+    return *points0Ptr_;
 }
 
 
@@ -92,22 +182,19 @@ const Foam::lumpedPointMovement&
 Foam::lumpedPointDisplacementPointPatchVectorField::movement() const
 {
     const objectRegistry& obr = this->patch().boundaryMesh().mesh().db();
-    const lumpedPointIOMovement* ptr =
-        lumpedPointIOMovement::lookupInRegistry(obr);
+
+    lumpedPointIOMovement* ptr =
+        lumpedPointIOMovement::getMovementObject(obr);
 
     if (ptr)
     {
         return *ptr; // Already exists
     }
 
-    // create and register with this patch as the owner
-    autoPtr<lumpedPointIOMovement> obj = lumpedPointIOMovement::New
-    (
-        obr,
-        this->patch().index()
-    );
+    // Create and register with this patch as the owner
+    ptr = lumpedPointIOMovement::New(obr, this->patch().index()).ptr();
 
-    return objectRegistry::store(obj);
+    return objectRegistry::store(ptr);
 }
 
 
@@ -120,7 +207,8 @@ lumpedPointDisplacementPointPatchVectorField
     const DimensionedField<vector, pointMesh>& iF
 )
 :
-    fixedValuePointPatchField<vector>(p, iF)
+    fixedValuePointPatchField<vector>(p, iF),
+    controllers_()
 {}
 
 
@@ -132,8 +220,20 @@ lumpedPointDisplacementPointPatchVectorField
     const dictionary& dict
 )
 :
-    fixedValuePointPatchField<vector>(p, iF, dict)
-{}
+    fixedValuePointPatchField<vector>(p, iF, dict),
+    controllers_()
+{
+    dict.readIfPresent("controllers", controllers_);
+
+    if (controllers_.empty())
+    {
+        WarningInFunction
+            << "No controllers specified, using all lumped points for patch: "
+            << this->patch().name() << nl << nl;
+    }
+
+    // controllers_ : check? remove duplicates?
+}
 
 
 Foam::lumpedPointDisplacementPointPatchVectorField::
@@ -145,7 +245,8 @@ lumpedPointDisplacementPointPatchVectorField
     const pointPatchFieldMapper& mapper
 )
 :
-    fixedValuePointPatchField<vector>(pf, p, iF, mapper)
+    fixedValuePointPatchField<vector>(pf, p, iF, mapper),
+    controllers_(pf.controllers_)
 {}
 
 
@@ -156,7 +257,8 @@ lumpedPointDisplacementPointPatchVectorField
     const DimensionedField<vector, pointMesh>& iF
 )
 :
-    fixedValuePointPatchField<vector>(pf, iF)
+    fixedValuePointPatchField<vector>(pf, iF),
+    controllers_(pf.controllers_)
 {}
 
 
@@ -166,16 +268,17 @@ Foam::lumpedPointDisplacementPointPatchVectorField::
 ~lumpedPointDisplacementPointPatchVectorField()
 {
     // de-register movement if in use and managed by this patch
-    const lumpedPointIOMovement* ptr = lumpedPointIOMovement::lookupInRegistry
-    (
-        this->patch().boundaryMesh().mesh().db()
-    );
+    lumpedPointIOMovement* ptr =
+        lumpedPointIOMovement::getMovementObject
+        (
+            this->patch().boundaryMesh().mesh().db()
+        );
 
     if (ptr && ptr->ownerId() == this->patch().index())
     {
         movement().coupler().shutdown();
 
-        const_cast<lumpedPointIOMovement*>(ptr)->checkOut();
+        ptr->checkOut();
     }
 }
 
@@ -189,58 +292,68 @@ void Foam::lumpedPointDisplacementPointPatchVectorField::updateCoeffs()
         return;
     }
 
+    const label timeIndex = this->db().time().timeIndex();
+
     enum Time::stopAtControls action = Time::stopAtControls::saUnknown;
 
-    const bool masterPatch = (movement().ownerId() == this->patch().index());
-    if (masterPatch)
+    if (movement().ownerId() == this->patch().index())
     {
-        if (lumpedPointIOMovement::debug)
+        // The ownerId is always the lowest patch number,
+        // thus will always be triggered first
+
+        if (lumpedPointMovement::debug)
         {
             Pout<<"masterPatch: " << this->patch().index() << endl;
         }
 
         const polyMesh& mesh = this->patch().boundaryMesh().mesh().mesh();
 
-        // need face 'zones' for calculating forces
-        // likely need bounding box for the movement
-        // -> do both now if required
+        // Mapping for calculating forces
+        // as well as face point interpolation
         if (!movement().hasMapping())
         {
-            const_cast<lumpedPointMovement&>(movement()).setMapping
+            // Add mapping for calculating forces/moments
+            setPatchControls
             (
-                mesh,
-                // All patches of this type
-                patchIds
+                static_cast<const pointVectorField&>
                 (
-                    static_cast<const pointVectorField&>
-                    (
-                        this->internalField()
-                    )
+                    this->internalField()
                 ),
                 this->points0()
             );
         }
 
+        int triggered = 0;
 
         if
         (
-            movement().coupler().initialized()
-         || !movement().coupler().slaveFirst()
+            movement().coupler().slaveFirst()
+         && !movement().coupler().initialized()
         )
+        {
+            // Master does nothing yet, but slave will
+            triggered = -1;
+        }
+        else if (movement().couplingPending(timeIndex))
+        {
+            // Trigger is pending, or coupling not yet not initialized
+            triggered = 1;
+        }
+
+        if (triggered > 0)
         {
             // Synchronized for all processes
             List<vector> forces, moments;
             movement().forcesAndMoments(mesh, forces, moments);
 
-            if (lumpedPointIOMovement::debug)
+            if (lumpedPointMovement::debug)
             {
                 Pout<<"gatherForces: " << forces << " called from patch "
                     << this->patch().index() << endl;
 
                 if (Pstream::master())
                 {
-                    Pout<<"output forces to file: "
-                        << movement().locations() << " called from patch "
+                    Pout<<"output forces to file: called from patch "
                         << this->patch().index() << nl
                         <<"# " << forces.size() << " force entries" << nl
                         <<"# fx fy fz" << nl
@@ -252,29 +365,44 @@ void Foam::lumpedPointDisplacementPointPatchVectorField::updateCoeffs()
 
             if (Pstream::master())
             {
-                movement().writeData(forces, moments, &(db().time()));
+                movement().writeData(forces, moments, &(this->db().time()));
 
                 // Signal external source to execute
                 movement().coupler().useSlave();
             }
         }
 
-        // Wait for slave to provide data (includes MPI barrier)
-        // and catch any abort information sent from slave
-        action = movement().coupler().waitForSlave();
+        if (triggered)
+        {
+            // Wait for slave to provide data (includes MPI barrier)
+            // and catch any abort information sent from slave
+            action = movement().coupler().waitForSlave();
 
-        // Read data passed back from external source - includes MPI barrier
-        const_cast<lumpedPointMovement&>(movement()).readState();
+            const_cast<lumpedPointMovement&>(movement()).readState();
+
+            movement().couplingCompleted(timeIndex);
+        }
     }
 
-    tmp<pointField> tdisp = movement().displacePoints
-    (
-        this->points0(),
-        this->patch().meshPoints()
-    );
+    if (!movement().hasInterpolator(this->patch()))
+    {
+        const_cast<lumpedPointMovement&>(movement()).setInterpolator
+        (
+            this->patch(),
+            this->points0()
+        );
+    }
+
+    tmp<pointField> tdisp =
+        movement().pointsDisplacement
+        (
+            this->patch(),
+            this->points0()
+        );
 
     this->operator==(tdisp);
 
+
     fixedValuePointPatchField<vector>::updateCoeffs();
 
     // Process any abort information sent from slave
@@ -293,6 +421,12 @@ void Foam::lumpedPointDisplacementPointPatchVectorField::write(Ostream& os)
 const
 {
     pointPatchField<vector>::write(os);
+
+    if (controllers_.size())
+    {
+        os.writeEntry("controllers", controllers_);
+    }
+
     writeEntry("value", os);
 }
 
diff --git a/src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.H b/src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.H
index 5098bb84e0..e5ae9c8ad7 100644
--- a/src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.H
+++ b/src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,6 +42,7 @@ Description
     {
         type    lumpedPointDisplacement;
         value   uniform (0 0 0);
+        controllers ( controllerName1 controllerName2 );
     }
     \endverbatim
 
@@ -67,7 +68,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class interpolationWeights;
 
 /*---------------------------------------------------------------------------*\
@@ -78,7 +79,14 @@ class lumpedPointDisplacementPointPatchVectorField
 :
     public fixedValuePointPatchField<vector>
 {
-    // Private data
+    // Private Data
+
+        //- Names of the movement controller(s) in use
+        wordList controllers_;
+
+        //- Backup method for getting "points0" without a motion solver
+        mutable autoPtr<pointIOField> points0Ptr_;
+
 
     //- Convenience typedefs
     typedef lumpedPointDisplacementPointPatchVectorField patchType;
@@ -167,7 +175,27 @@ public:
     virtual ~lumpedPointDisplacementPointPatchVectorField();
 
 
-    // Member functions
+    // Member Functions
+
+        //- The controller names for this patch
+        const wordList& controllers() const
+        {
+            return controllers_;
+        }
+
+        //- Set all patch controls for patches of this type
+        static label setPatchControls
+        (
+            const pointVectorField& pvf,
+            const pointField& points0
+        );
+
+        //- Set all patch controls for patches of this type
+        static label setInterpolators
+        (
+            const pointVectorField& pvf,
+            const pointField& points0
+        );
 
         //- The ids for all patches of this type
         static labelList patchIds(const pointVectorField& pvf);
diff --git a/src/lumpedPointMotion/state/lumpedPointState.C b/src/lumpedPointMotion/state/lumpedPointState.C
index 2beb0630a7..83a847664b 100644
--- a/src/lumpedPointMotion/state/lumpedPointState.C
+++ b/src/lumpedPointMotion/state/lumpedPointState.C
@@ -26,7 +26,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "lumpedPointState.H"
-#include "demandDrivenData.H"
 #include "unitConversion.H"
 #include "EulerCoordinateRotation.H"
 #include "IFstream.H"
@@ -44,6 +43,11 @@ Foam::lumpedPointState::formatNames
 });
 
 
+Foam::scalar Foam::lumpedPointState::visLength = 0.1;
+
+bool Foam::lumpedPointState::visUnused = true;
+
+
 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
 //! \cond fileScope
@@ -69,7 +73,7 @@ static Foam::string getLineNoComment
 
 void Foam::lumpedPointState::calcRotations() const
 {
-    rotationPtr_ = new tensorField(angles_.size());
+    rotationPtr_.reset(new tensorField(angles_.size()));
 
     auto rotIter = rotationPtr_->begin();
 
@@ -83,20 +87,26 @@ void Foam::lumpedPointState::calcRotations() const
 }
 
 
-void Foam::lumpedPointState::readDict(const dictionary& dict)
+void Foam::lumpedPointState::readDict
+(
+    const dictionary& dict,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
 {
     dict.readEntry("points", points_);
     dict.readEntry("angles", angles_);
     order_ =
         quaternion::eulerOrderNames.getOrDefault
         (
-            "order",
+            "rotationOrder",
             dict,
-            quaternion::eulerOrder::ZXZ
+            rotOrder
         );
-    degrees_ = dict.getOrDefault("degrees", false);
 
-    deleteDemandDrivenData(rotationPtr_);
+    degrees_ = dict.getOrDefault("degrees", degrees);
+
+    rotationPtr_.reset(nullptr);
 }
 
 
@@ -122,43 +132,82 @@ Foam::lumpedPointState::lumpedPointState(const lumpedPointState& rhs)
 {}
 
 
-Foam::lumpedPointState::lumpedPointState(const pointField& pts)
+Foam::lumpedPointState::lumpedPointState
+(
+    const pointField& pts,
+    const vectorField& ang,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
 :
     points_(pts),
-    angles_(points_.size(), Zero),
-    order_(quaternion::eulerOrder::ZXZ),
-    degrees_(false),
-    rotationPtr_(nullptr)
-{}
-
-
-Foam::lumpedPointState::lumpedPointState(tmp<pointField>& pts)
-:
-    points_(pts),
-    angles_(points_.size(), Zero),
-    order_(quaternion::eulerOrder::ZXZ),
-    degrees_(false),
-    rotationPtr_(nullptr)
-{}
-
-
-Foam::lumpedPointState::lumpedPointState(const dictionary& dict)
-:
-    points_(),
-    angles_(),
-    order_(quaternion::eulerOrder::ZXZ),
-    degrees_(false),
+    angles_(ang),
+    order_(rotOrder),
+    degrees_(degrees),
     rotationPtr_(nullptr)
 {
-    readDict(dict);
+    if (points_.size() != angles_.size())
+    {
+        #ifdef FULLDEBUG
+        FatalErrorInFunction
+            << "Have " << points_.size() << " points but "
+            << angles_.size() << " angles" << nl
+            exit(FatalError);
+        #else
+        WarningInFunction
+            << "Have " << points_.size() << " points but "
+            << angles_.size() << " angles, resizing angles to match" << nl;
+        #endif
+
+        angles_.resize(points_.size(), Zero);
+    }
 }
 
 
-// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+Foam::lumpedPointState::lumpedPointState
+(
+    const pointField& pts,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
+:
+    points_(pts),
+    angles_(points_.size(), Zero),
+    order_(rotOrder),
+    degrees_(degrees),
+    rotationPtr_(nullptr)
+{}
 
-Foam::lumpedPointState::~lumpedPointState()
+
+Foam::lumpedPointState::lumpedPointState
+(
+    tmp<pointField>& pts,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
+:
+    points_(pts),
+    angles_(points_.size(), Zero),
+    order_(rotOrder),
+    degrees_(false),
+    rotationPtr_(nullptr)
+{}
+
+
+Foam::lumpedPointState::lumpedPointState
+(
+    const dictionary& dict,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
+:
+    points_(),
+    angles_(),
+    order_(rotOrder),
+    degrees_(degrees),
+    rotationPtr_(nullptr)
 {
-    deleteDemandDrivenData(rotationPtr_);
+    readDict(dict, rotOrder, degrees);
 }
 
 
@@ -171,7 +220,16 @@ void Foam::lumpedPointState::operator=(const lumpedPointState& rhs)
     order_   = rhs.order_;
     degrees_ = rhs.degrees_;
 
-    deleteDemandDrivenData(rotationPtr_);
+    rotationPtr_.reset(nullptr);
+}
+
+
+void Foam::lumpedPointState::operator+=(const point& origin)
+{
+    for (point& p : points_)
+    {
+        p += origin;
+    }
 }
 
 
@@ -209,11 +267,16 @@ void Foam::lumpedPointState::relax
 
     angles_ = convert*prev.angles_ + alpha*(angles_ - convert*prev.angles_);
 
-    deleteDemandDrivenData(rotationPtr_);
+    rotationPtr_.reset(nullptr);
 }
 
 
-bool Foam::lumpedPointState::readPlain(Istream& is)
+bool Foam::lumpedPointState::readPlain
+(
+    Istream& is,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
 {
     // Assume generic input stream so we can do line-based input
     ISstream& iss = dynamic_cast<ISstream&>(is);
@@ -226,8 +289,8 @@ bool Foam::lumpedPointState::readPlain(Istream& is)
         isstr >> count;
     }
 
-    points_.setSize(count);
-    angles_.setSize(count);
+    points_.resize(count);
+    angles_.resize(count);
 
     count = 0;
     forAll(points_, i)
@@ -242,21 +305,26 @@ bool Foam::lumpedPointState::readPlain(Istream& is)
         ++count;
     }
 
-    points_.setSize(count);
-    angles_.setSize(count);
+    points_.resize(count);
+    angles_.resize(count);
     order_ = quaternion::eulerOrder::ZXZ;
     degrees_ = false;
 
-    deleteDemandDrivenData(rotationPtr_);
+    rotationPtr_.reset(nullptr);
 
     return count;
 }
 
 
-bool Foam::lumpedPointState::readData(Istream& is)
+bool Foam::lumpedPointState::readData
+(
+    Istream& is,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
 {
     dictionary dict(is);
-    readDict(dict);
+    readDict(dict, rotOrder, degrees);
 
     return points_.size();
 }
@@ -292,9 +360,9 @@ void Foam::lumpedPointState::writePlain(Ostream& os) const
 
     forAll(points_, i)
     {
-        const vector& pt = points_[i];
+        const vector& p = points_[i];
 
-        os  << pt.x() << ' ' << pt.y() << ' ' << pt.z();
+        os  << p.x() << ' ' << p.y() << ' ' << p.z();
 
         if (i < angles_.size())
         {
@@ -312,8 +380,10 @@ void Foam::lumpedPointState::writePlain(Ostream& os) const
 
 bool Foam::lumpedPointState::readData
 (
-    const inputFormatType& fmt,
-    const fileName& file
+    const inputFormatType fmt,
+    const fileName& file,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
 )
 {
     bool ok = false;
@@ -323,11 +393,11 @@ bool Foam::lumpedPointState::readData
 
         if (fmt == inputFormatType::PLAIN)
         {
-            ok = this->readPlain(is);
+            ok = this->readPlain(is, rotOrder, degrees);
         }
         else
         {
-            ok = this->readData(is);
+            ok = this->readData(is, rotOrder, degrees);
         }
     }
 
@@ -375,7 +445,7 @@ bool Foam::lumpedPointState::readData
             toBelow << points_ << angles_ << degrees_;
         }
 
-        deleteDemandDrivenData(rotationPtr_);
+        rotationPtr_.reset(nullptr);
 
         // MPI barrier
         Pstream::scatter(ok);
diff --git a/src/lumpedPointMotion/state/lumpedPointState.H b/src/lumpedPointMotion/state/lumpedPointState.H
index 392c97a2e6..2871bf240a 100644
--- a/src/lumpedPointMotion/state/lumpedPointState.H
+++ b/src/lumpedPointMotion/state/lumpedPointState.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,7 +37,7 @@ Description
         Property    | Description                           | Required | Default
         points      | List of points                        | yes |
         angles      | List of Euler rotation angles         | yes |
-        order       | The Euler-angle rotation order        | no  | zxz
+        rotationOrder | The Euler-angle rotation order      | no  | zxz
         degrees     | Rotation angles in degrees            | no  | false
     \endtable
 
@@ -77,6 +77,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class Istream;
 class Ostream;
 
@@ -88,17 +89,20 @@ class lumpedPointState
 {
 public:
 
-    //- Input format types
-    enum class inputFormatType
-    {
-        PLAIN,      //!< "plain" is a simple ASCII format
-        DICTIONARY  //!< "dictionary" is the OpenFOAM dictionary format
-    };
+    // Data Types
 
-    // Static data
+        //- Input format types
+        enum class inputFormatType
+        {
+            PLAIN,      //!< "plain" is a simple ASCII format
+            DICTIONARY  //!< "dictionary" is the OpenFOAM dictionary format
+        };
 
-    //- Names for the input format types
-    static const Enum<inputFormatType> formatNames;
+
+    // Static Data
+
+        //- Names for the input format types
+        static const Enum<inputFormatType> formatNames;
 
 
 private:
@@ -118,37 +122,75 @@ private:
         bool degrees_;
 
         //- Tensor rotation of lumped points
-        mutable tensorField* rotationPtr_;
+        mutable unique_ptr<tensorField> rotationPtr_;
 
 
     // Private Member Functions
 
         void calcRotations() const;
 
-        void readDict(const dictionary& dict);
+        void readDict
+        (
+            const dictionary& dict,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
 public:
 
+    // Public Data
+
+        //- Enable/disable visualization of unused points
+        static bool visUnused;
+
+        //- The length for visualization triangles
+        static scalar visLength;
+
+
     // Constructors
 
-        //- Construct null
+        //- Default construct
         lumpedPointState();
 
-        //- Copy constructor
+        //- Copy construct
         lumpedPointState(const lumpedPointState& rhs);
 
-        //- Construct from points with zero-rotation
-        lumpedPointState(const pointField& pts);
+        //- Copy construct from points and angles
+        lumpedPointState
+        (
+            const pointField& pts,
+            const vectorField& ang,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
+
+        //- Copy construct from points with zero-rotation
+        explicit lumpedPointState
+        (
+            const pointField& pts,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
         //- Construct from points with zero-rotation
-        lumpedPointState(tmp<pointField>& pts);
+        explicit lumpedPointState
+        (
+            tmp<pointField>& pts,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
         //- Construct from dictionary
-        lumpedPointState(const dictionary& dict);
+        explicit lumpedPointState
+        (
+            const dictionary& dict,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
 
     //- Destructor
-    virtual ~lumpedPointState();
+    virtual ~lumpedPointState() = default;
 
 
     // Member Functions
@@ -175,6 +217,12 @@ public:
         //  Zero and negative values are ignored.
         void scalePoints(const scalar scaleFactor);
 
+        //- The Euler-angle rotation order
+        inline quaternion::eulerOrder rotationOrder() const;
+
+        //-  Rotation angles in degrees
+        inline bool degrees() const;
+
         //- Relax the state
         //  alpha = 1 : no relaxation
         //  alpha < 1 : relaxation
@@ -182,7 +230,12 @@ public:
         void relax(const scalar alpha, const lumpedPointState& prev);
 
         //- Read input as dictionary content
-        bool readData(Istream& is);
+        bool readData
+        (
+            Istream& is,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
         //- Output as dictionary content
         bool writeData(Ostream& os) const;
@@ -191,24 +244,42 @@ public:
         void writeDict(Ostream& os) const;
 
         //- Read input as plain content
-        bool readPlain(Istream& is);
+        bool readPlain
+        (
+            Istream& is,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
         //- Output as plain content
         void writePlain(Ostream& os) const;
 
         //- Read input from file (master process only) using specified format
-        bool readData(const inputFormatType& fmt, const fileName& file);
+        bool readData
+        (
+            const inputFormatType fmt,
+            const fileName& file,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
-        //- Output as VTK file for debugging/visualization
-        //  The points are joined as lines, the rotation is visualized
-        //  by planes, write as VTK PolyData format.
-        void writeVTP(const fileName& outputFile, const vector& axis) const;
+        //- Output points/rotations as VTK file for debugging/visualization
+        //  The points are written as vertices, rotation as a triangle
+        void writeVTP
+        (
+            const fileName& outputFile,
+            const labelListList& lines = labelListList(),
+            const labelList& pointIds = labelList::null()
+        ) const;
 
 
     // Member Operators
 
         //- Assignment operator
         void operator=(const lumpedPointState& rhs);
+
+        //- Shift points by specified origin
+        void operator+=(const point& origin);
 };
 
 
diff --git a/src/lumpedPointMotion/state/lumpedPointStateI.H b/src/lumpedPointMotion/state/lumpedPointStateI.H
index bf3784c393..ef3cc4a0a0 100644
--- a/src/lumpedPointMotion/state/lumpedPointStateI.H
+++ b/src/lumpedPointMotion/state/lumpedPointStateI.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -66,4 +66,17 @@ inline const Foam::tensorField& Foam::lumpedPointState::rotations() const
 }
 
 
+inline Foam::quaternion::eulerOrder
+Foam::lumpedPointState::rotationOrder() const
+{
+    return order_;
+}
+
+
+inline bool Foam::lumpedPointState::degrees() const
+{
+    return degrees_;
+}
+
+
 // ************************************************************************* //
diff --git a/src/lumpedPointMotion/state/lumpedPointStateWriter.C b/src/lumpedPointMotion/state/lumpedPointStateWriter.C
index 5e49598748..a44d864b9e 100644
--- a/src/lumpedPointMotion/state/lumpedPointStateWriter.C
+++ b/src/lumpedPointMotion/state/lumpedPointStateWriter.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,30 +28,23 @@ License
 #include "lumpedPointState.H"
 #include "OFstream.H"
 #include "sliceRange.H"
-#include "axesRotation.H"
-#include "coordinateSystem.H"
 #include "foamVtkOutput.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-// file-local
-const static Foam::FixedList<Foam::point, 4> standardCorners
-{
-    {-0.1, -0.1, 0},
-    {+0.1, -0.1, 0},
-    {+0.1, +0.1, 0},
-    {-0.1, +0.1, 0}
-};
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::lumpedPointState::writeVTP
 (
     const fileName& outputFile,
-    const vector& axis
+    const labelListList& lines,
+    const labelList& pointIds
 ) const
 {
+    if (!Pstream::master())
+    {
+        // No extra information available from slaves, write on master only.
+        return;
+    }
+
     // local-to-global transformation tensors
     const tensorField& localToGlobal = rotations();
 
@@ -68,28 +61,52 @@ void Foam::lumpedPointState::writeVTP
         .beginVTKFile<vtk::fileTag::POLY_DATA>();
 
     //
-    // The 'spine' of lumped mass points
+    // Lumped mass points and connections,
+    // with triangles to visualize location/rotation
     //
     {
+        const label nPoints = 3*points_.size();  // 3 points per triangle
+        const label nPolys  = points_.size();
+
         format()
             .tag
             (
                 vtk::fileTag::PIECE,
-                vtk::fileAttr::NUMBER_OF_POINTS, points_.size(),
+                vtk::fileAttr::NUMBER_OF_POINTS, nPoints,
                 vtk::fileAttr::NUMBER_OF_VERTS,  points_.size(),
-                vtk::fileAttr::NUMBER_OF_LINES,  1
+                vtk::fileAttr::NUMBER_OF_LINES,  lines.size(),
+                vtk::fileAttr::NUMBER_OF_POLYS,  nPolys
             );
 
         // 'points'
         {
-            const uint64_t payLoad = vtk::sizeofData<float, 3>(points_.size());
+            const uint64_t payLoad = vtk::sizeofData<float, 3>(nPoints);
 
             format()
                 .tag(vtk::fileTag::POINTS)
                 .beginDataArray<float, 3>(vtk::dataArrayAttr::POINTS);
 
             format().writeSize(payLoad);
+
+            // The lumped points first
             vtk::writeList(format(), points_);
+
+            // Other points (for the triangles) next
+            forAll(points_, posi)
+            {
+                const point& origin = points_[posi];
+                const tensor& rotTensor =
+                (
+                    posi < localToGlobal.size()
+                  ? localToGlobal[posi]
+                  : pTraits<tensor>::I
+                );
+
+                // Local-to-global rotation and translation
+                vtk::write(format(), 2*visLength*rotTensor.cx() + origin);
+                vtk::write(format(), 1*visLength*rotTensor.cy() + origin);
+            }
+
             format().flush();
 
             format()
@@ -140,16 +157,25 @@ void Foam::lumpedPointState::writeVTP
         // <Lines>
         format().tag(vtk::fileTag::LINES);
 
+        label nLinePoints = 0;
+        for (const labelList& linePoints : lines)
+        {
+            nLinePoints += linePoints.size();
+        }
+
         //
         // 'connectivity'
         //
         {
-            const uint64_t payLoad = vtk::sizeofData<label>(points_.size());
+            const uint64_t payLoad = vtk::sizeofData<label>(nLinePoints);
 
             format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
             format().writeSize(payLoad);
 
-            vtk::writeIdentity(format(), points_.size());
+            for (const labelList& linePoints : lines)
+            {
+                vtk::writeList(format(), linePoints);
+            }
 
             format().flush();
 
@@ -158,102 +184,52 @@ void Foam::lumpedPointState::writeVTP
 
         //
         // 'offsets'  (connectivity offsets)
-        // = single line
+        // = N lines
         //
         {
-            const uint64_t payLoad = vtk::sizeofData<label>(1);
+            const uint64_t payLoad = vtk::sizeofData<label>(lines.size());
 
             format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
             format().writeSize(payLoad);
 
-            format().write(points_.size());
+            nLinePoints = 0;
+            for (const labelList& linePoints : lines)
+            {
+                nLinePoints += linePoints.size();
+                format().write(nLinePoints);
+            }
+
             format().flush();
 
             format().endDataArray();
         }
 
         format().endTag(vtk::fileTag::LINES);
-        format().endPiece();
-    }
+        // </Lines>
 
-    // Standard corners in local axis
-    FixedList<point, 4> corners;
-
-    {
-        coordinateRotations::axes orient(axis);
-        coordinateSystem cornerTransform(orient);
-
-        forAll(standardCorners, corni)
-        {
-            corners[corni] = cornerTransform.transform(standardCorners[corni]);
-        }
-    }
-
-    //
-    // Planes to visualize location/rotation
-    //
-    {
-        const label nPoints = 4*points_.size();  // 4 points per quad
-        const label nPolys  = points_.size();
-
-        format()
-            .tag
-            (
-                vtk::fileTag::PIECE,
-                vtk::fileAttr::NUMBER_OF_POINTS, nPoints,
-                vtk::fileAttr::NUMBER_OF_POLYS,  nPolys
-            );
-
-        // 'points'
-        {
-            const uint64_t payLoad = vtk::sizeofData<float, 3>(nPoints);
-
-            format()
-                .tag(vtk::fileTag::POINTS)
-                .beginDataArray<float, 3>(vtk::dataArrayAttr::POINTS);
-
-            format().writeSize(payLoad);
-
-            forAll(points_, posI)
-            {
-                const point& origin = points_[posI];
-                const tensor& rotTensor =
-                (
-                    posI < localToGlobal.size()
-                  ? localToGlobal[posI]
-                  : pTraits<tensor>::I
-                );
-
-
-                for (const point& cornerPt : corners)
-                {
-                    // Local-to-global rotation and translation
-                    const point pt = (rotTensor & cornerPt) + origin;
-
-                    vtk::write(format(), pt);
-                }
-            }
-
-            format().flush();
-
-            format()
-                .endDataArray()
-                .endTag(vtk::fileTag::POINTS);
-        }
 
         // <Polys>
         format().tag(vtk::fileTag::POLYS);
 
         //
-        // 'connectivity' - 4 points (ie, quad)
+        // 'connectivity' - 3 points (ie, tri)
+        // origins appear first, followed by a point pair for each triangle
+        // Eg,
+        // - tri 0: (0 N   N+1)
+        // - tri 1: (1 N+2 N+3)
         //
         {
-            const uint64_t payLoad = vtk::sizeofData<label>(4*nPolys);
+            const uint64_t payLoad = vtk::sizeofData<label>(3*nPolys);
 
             format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
             format().writeSize(payLoad);
 
-            vtk::writeIdentity(format(), 4*nPolys);
+            for (label pointi=0, nei=nPolys; pointi < nPolys; ++pointi)
+            {
+                format().write(pointi);
+                format().write(nei); ++nei;
+                format().write(nei); ++nei;
+            }
 
             format().flush();
 
@@ -262,7 +238,7 @@ void Foam::lumpedPointState::writeVTP
 
         //
         // 'offsets'  (connectivity offsets)
-        // = single quad
+        // = single tri
         //
         {
             const uint64_t payLoad = vtk::sizeofData<label>(nPolys);
@@ -270,7 +246,7 @@ void Foam::lumpedPointState::writeVTP
             format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
             format().writeSize(payLoad);
 
-            for (const label off : sliceRange(4, nPolys, 4))
+            for (const label off : sliceRange(3, nPolys, 3))
             {
                 format().write(off);
             }
@@ -280,17 +256,27 @@ void Foam::lumpedPointState::writeVTP
         }
 
         format().endTag(vtk::fileTag::POLYS);
+        // </Polys>
 
-#if 0
+
+        // CELL_DATA
         format().beginCellData();
 
-        // zone Id
+        // point id
         {
-            const uint64_t payLoad = vtk::sizeofData<label>(nPolys);
+            const uint64_t payLoad =
+                vtk::sizeofData<label>(points_.size() + lines.size());
 
-            format().beginDataArray<label>("zoneId");
+            format().beginDataArray<label>("pointId");
             format().writeSize(payLoad);
 
+            // <Verts>
+            vtk::writeIdentity(format(), points_.size());
+
+            // <Lines>
+            vtk::write(format(), label(-1), lines.size());
+
+            // <Poly>
             vtk::writeIdentity(format(), nPolys);
 
             format().flush();
@@ -298,16 +284,99 @@ void Foam::lumpedPointState::writeVTP
             format().endDataArray();
         }
 
+        // original id
+        if (pointIds.size() == points_.size())
+        {
+            const uint64_t payLoad =
+                vtk::sizeofData<label>(points_.size() + lines.size());
+
+            format().beginDataArray<label>("originalId");
+            format().writeSize(payLoad);
+
+            // <Verts>
+            vtk::writeList(format(), pointIds);
+
+            // <Lines>
+            vtk::write(format(), label(-1), lines.size());
+
+            // <Poly>
+            vtk::writeList(format(), pointIds);
+
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        // line id
+        {
+            const uint64_t payLoad =
+                vtk::sizeofData<label>(points_.size() + lines.size());
+
+            format().beginDataArray<label>("lineId");
+            format().writeSize(payLoad);
+
+            // <Verts>
+            vtk::write(format(), label(-1), points_.size());
+
+            // <Lines>
+            vtk::writeIdentity(format(), lines.size());
+
+            // <Poly>
+            vtk::write(format(), label(-1), nPolys);
+
+            format().flush();
+
+            format().endDataArray();
+        }
+
         format().endCellData();
-#endif
+
+
+        // POINT_DATA
+        format().beginPointData();
+
+        // point id
+        {
+            const uint64_t payLoad = vtk::sizeofData<label>(nPoints);
+
+            format().beginDataArray<label>("pointId");
+            format().writeSize(payLoad);
+
+            // The lumped points first
+            vtk::writeIdentity(format(), points_.size());
+
+            // Tag other (triangle) points as -1
+            vtk::write(format(), label(-1), 2*points_.size());
+
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        // original id
+        if (pointIds.size() == points_.size())
+        {
+            const uint64_t payLoad = vtk::sizeofData<label>(nPoints);
+
+            format().beginDataArray<label>("originalId");
+            format().writeSize(payLoad);
+
+            // The lumped points first
+            vtk::writeList(format(), pointIds);
+
+            // Tag other (triangle) points as -1
+            vtk::write(format(), label(-1), 2*points_.size());
+
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        format().endPointData();
 
         format().endPiece();
     }
 
-    // Finally
-    // could add a 'ghost' level above to visualize extrapolated values
-    // draw as two triangles to distinguish from real levels ...
-
     format().endTag(vtk::fileTag::POLY_DATA)
         .endVTKFile();
 }
diff --git a/src/lumpedPointMotion/tools/lumpedPointTools.C b/src/lumpedPointMotion/tools/lumpedPointTools.C
index b589d45c32..d58b0c4420 100644
--- a/src/lumpedPointMotion/tools/lumpedPointTools.C
+++ b/src/lumpedPointMotion/tools/lumpedPointTools.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,14 +28,13 @@ License
 #include "IFstream.H"
 #include "IOobjectList.H"
 #include "volFields.H"
-
+#include "points0MotionSolver.H"
 #include "lumpedPointDisplacementPointPatchVectorField.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    // file-scope
     template<class GeoFieldType>
     static autoPtr<GeoFieldType> loadPointField
     (
@@ -72,21 +71,35 @@ namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 Foam::List<Foam::lumpedPointStateTuple>
-Foam::lumpedPointTools::lumpedPointStates(Istream& is)
+Foam::lumpedPointTools::lumpedPointStates
+(
+    const dictionary& dict,
+    quaternion::eulerOrder rotOrder,
+    bool degrees
+)
 {
-    dictionary contents(is);
-    List<dictionary> entries(contents.lookup("response"));
+    quaternion::eulerOrderNames.readIfPresent
+    (
+        "rotationOrder",
+        dict,
+        rotOrder
+    );
+
+    dict.readIfPresent("degrees", degrees);
+
+    Info<<"Reading states\n";
+    List<dictionary> entries(dict.lookup("response"));
 
     DynamicList<Tuple2<scalar, lumpedPointState>> states(entries.size());
 
-    for (const dictionary& dict : entries)
+    for (const dictionary& subDict : entries)
     {
         states.append
         (
             lumpedPointStateTuple
             (
-                dict.get<scalar>("time"),
-                lumpedPointState(dict)
+                subDict.get<scalar>("time"),
+                lumpedPointState(subDict)
             )
         );
     }
@@ -96,35 +109,38 @@ Foam::lumpedPointTools::lumpedPointStates(Istream& is)
 
 
 Foam::List<Foam::lumpedPointStateTuple>
-Foam::lumpedPointTools::lumpedPointStates(const fileName& file)
+Foam::lumpedPointTools::lumpedPointStates
+(
+    Istream& is,
+    quaternion::eulerOrder rotOrder,
+    bool degrees
+)
+{
+    dictionary dict(is);
+    return lumpedPointStates(dict, rotOrder, degrees);
+}
+
+
+Foam::List<Foam::lumpedPointStateTuple>
+Foam::lumpedPointTools::lumpedPointStates
+(
+    const fileName& file,
+    quaternion::eulerOrder rotOrder,
+    bool degrees
+)
 {
     IFstream is(file);
-    return lumpedPointStates(is);
+    return lumpedPointStates(is, rotOrder, degrees);
 }
 
 
 Foam::pointIOField
 Foam::lumpedPointTools::points0Field(const polyMesh& mesh)
 {
-    pointIOField pts
-    (
-        IOobject
-        (
-            "points",
-            mesh.time().constant(),
-            polyMesh::meshSubDir,
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE,
-            false // Do not re-register
-        )
-    );
-
-    return pts;
+    return pointIOField(points0MotionSolver::points0IO(mesh));
 }
 
 
-
 Foam::labelList
 Foam::lumpedPointTools::lumpedPointPatchList(const pointVectorField& pvf)
 {
@@ -141,19 +157,127 @@ Foam::labelList Foam::lumpedPointTools::lumpedPointPatchList
 
     pointMesh pMesh(mesh);
 
-    autoPtr<pointVectorField> displacePtr = loadPointField<pointVectorField>
-    (
-        pMesh,
-        objects0.findObject("pointDisplacement")
-    );
+    autoPtr<pointVectorField> displacePtr =
+        loadPointField<pointVectorField>
+        (
+            pMesh,
+            objects0.findObject("pointDisplacement")
+        );
 
-    if (!displacePtr.valid())
+    if (!displacePtr)
     {
-        Info<< "no valid pointDisplacement" << endl;
+        Info<< "No valid pointDisplacement" << endl;
         return labelList();
     }
 
-    return lumpedPointPatchList(displacePtr());
+    return lumpedPointPatchList(*displacePtr);
+}
+
+
+Foam::label
+Foam::lumpedPointTools::setPatchControls
+(
+    const pointVectorField& pvf,
+    const pointField& points0
+)
+{
+    return
+        lumpedPointDisplacementPointPatchVectorField::setPatchControls
+        (
+            pvf,
+            points0
+        );
+}
+
+
+Foam::label Foam::lumpedPointTools::setPatchControls
+(
+    const fvMesh& mesh,
+    const pointField& points0
+)
+{
+    IOobjectList objects0(mesh, "0");
+
+    pointMesh pMesh(mesh);
+
+    autoPtr<pointVectorField> displacePtr =
+        loadPointField<pointVectorField>
+        (
+            pMesh,
+            objects0.findObject("pointDisplacement")
+        );
+
+    if (!displacePtr)
+    {
+        Info<< "No valid pointDisplacement" << endl;
+        return 0;
+    }
+
+    return setPatchControls(*displacePtr, points0);
+}
+
+
+Foam::label Foam::lumpedPointTools::setPatchControls
+(
+    const fvMesh& mesh
+)
+{
+    pointIOField points0(points0Field(mesh));
+
+    return setPatchControls(mesh, points0);
+}
+
+
+Foam::label Foam::lumpedPointTools::setInterpolators
+(
+    const pointVectorField& pvf,
+    const pointField& points0
+)
+{
+    return
+        lumpedPointDisplacementPointPatchVectorField::setInterpolators
+        (
+            pvf,
+            points0
+        );
+}
+
+
+Foam::label Foam::lumpedPointTools::setInterpolators
+(
+    const fvMesh& mesh,
+    const pointField& points0
+)
+{
+    IOobjectList objects0(mesh, "0");
+
+    pointMesh pMesh(mesh);
+
+    autoPtr<pointVectorField> displacePtr =
+        loadPointField<pointVectorField>
+        (
+            pMesh,
+            objects0.findObject("pointDisplacement")
+        );
+
+    if (!displacePtr)
+    {
+        Info<< "No valid pointDisplacement" << endl;
+        return 0;
+    }
+
+    return setInterpolators(*displacePtr, points0);
+}
+
+
+Foam::label Foam::lumpedPointTools::setInterpolators
+(
+    const fvMesh& mesh
+)
+{
+    pointIOField points0(points0Field(mesh));
+
+    return setInterpolators(mesh, points0);
 }
 
 
diff --git a/src/lumpedPointMotion/tools/lumpedPointTools.H b/src/lumpedPointMotion/tools/lumpedPointTools.H
index 52b693a09f..f39d00c15c 100644
--- a/src/lumpedPointMotion/tools/lumpedPointTools.H
+++ b/src/lumpedPointMotion/tools/lumpedPointTools.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,11 +39,10 @@ SourceFiles
 #define lumpedPointTools_H
 
 #include "labelList.H"
-#include "polyMesh.H"
+#include "fvMesh.H"
 #include "pointMesh.H"
 #include "pointFields.H"
 #include "Tuple2.H"
-
 #include "lumpedPointState.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -51,6 +50,7 @@ SourceFiles
 namespace Foam
 {
 
+// Typedefs
 typedef Tuple2<scalar, lumpedPointState> lumpedPointStateTuple;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -58,23 +58,63 @@ typedef Tuple2<scalar, lumpedPointState> lumpedPointStateTuple;
 namespace lumpedPointTools
 {
 
+//- Load a list of states from a dictionary
+List<lumpedPointStateTuple> lumpedPointStates
+(
+    const dictionary& dict,
+    quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+    bool degrees = false
+);
+
 //- Load a list of states from an Istream
-List<lumpedPointStateTuple> lumpedPointStates(Istream& is);
+List<lumpedPointStateTuple> lumpedPointStates
+(
+    Istream& is,
+    quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+    bool degrees = false
+);
 
 //- Load a list of states from a file
-List<lumpedPointStateTuple> lumpedPointStates(const fileName& file);
-
+List<lumpedPointStateTuple> lumpedPointStates
+(
+    const fileName& file,
+    quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+    bool degrees = false
+);
 
 //- Return the 0 or constant points field
 pointIOField points0Field(const polyMesh& mesh);
 
-//- Return the patch-ids associated a "lumpedPointDisplacement" type
+//- Return the patch-ids associated with a "lumpedPointDisplacement" type
 labelList lumpedPointPatchList(const pointVectorField& pvf);
 
 //- Get the "pointDisplacement" at time 0 and use that to determine which
-//  patches have a "lumpedPointDisplacement" type
+//- patches have a "lumpedPointDisplacement" type
 labelList lumpedPointPatchList(const polyMesh& mesh);
 
+//- Return the patch-ids associated with a "lumpedPointDisplacement" type
+label setPatchControls(const pointVectorField& pvf, const pointField& points0);
+
+//- Get the "pointDisplacement" at time 0 and use that to determine which
+//- patches have a "lumpedPointDisplacement" type
+label setPatchControls(const fvMesh& mesh, const pointField& points0);
+
+//- Get the "pointDisplacement" at time 0 and use that to determine which
+//- patches have a "lumpedPointDisplacement" type
+label setPatchControls(const fvMesh& mesh);
+
+
+//- Return the patch-ids associated with a "lumpedPointDisplacement" type
+label setInterpolators(const pointVectorField& pvf, const pointField& points0);
+
+//- Get the "pointDisplacement" at time 0 and use that to determine which
+//- patches have a "lumpedPointDisplacement" type
+label setInterpolators(const fvMesh& mesh, const pointField& points0);
+
+//- Get the "pointDisplacement" at time 0 and use that to determine which
+//- patches have a "lumpedPointDisplacement" type
+label setInterpolators(const fvMesh& mesh);
+
 } // End namespace lumpedPointTools
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/incompressible/lumpedPointMotion/building/Allclean b/tutorials/incompressible/lumpedPointMotion/building/Allclean
index 970e3bfab8..e560111894 100755
--- a/tutorials/incompressible/lumpedPointMotion/building/Allclean
+++ b/tutorials/incompressible/lumpedPointMotion/building/Allclean
@@ -1,5 +1,8 @@
 #!/bin/sh
-cd "${0%/*}" || exit    # Run from this directory
+cd "${0%/*}" || exit                                # Run from this directory
+#------------------------------------------------------------------------------
+
+(cd code && ./Allclean)
 
 (cd steady && ./Allclean)
 
diff --git a/tutorials/incompressible/lumpedPointMotion/building/Allrun b/tutorials/incompressible/lumpedPointMotion/building/Allrun
index 3132b323ff..db3fdc34c5 100755
--- a/tutorials/incompressible/lumpedPointMotion/building/Allrun
+++ b/tutorials/incompressible/lumpedPointMotion/building/Allrun
@@ -3,94 +3,45 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
-# 1) First run steady-state to establish a good initial field.
-# 2) Copy the latest state-state results for the transient case,
+# 1) Run steady-state to establish a good initial field
+# 2) Copy state-state results -> transient case,
 #    but need to copy the pointDisplacement from the 0/ directory
 #    since it will not have been used for the steady-state case
 # 3) Relocate this initial solution to coincide with the first deltaT
 #    to avoid overwriting the 0/ directory at all.
 
-#
-# copyParallelPointDisplacement caseDir timeName
-#
-# Copy pointDisplacement from caseDir/0/ to caseDir/timeName/
-#
-copyParallelPointDisplacement()
-{
-    local src=$1
-    local dstTime=$2
-    local file=pointDisplacement
-
-    [ -d "$src" ] || {
-        echo "Error: no directory: $src"
-        return 1
-    }
-
-    # Copy select directories
-    echo "    copy processor '$file' from 0/ -> $dstTime"
-    if [ -n "$dstTime" ]
-    then
-    (
-        cd "$src" || exit
-
-        for proc in processor*
-        do
-            [ -d "$proc/0" -a -d "$proc/$dstTime" ] && \
-                cp $proc/0/$file $proc/$dstTime/$file
-        done
-    )
-    else
-        echo "    no destination time"
-    fi
-
-    # Restart from latestTime
-    foamDictionary $src/system/controlDict \
-        -entry startFrom -set latestTime
-
-    deltaT=$(foamDictionary $src/system/controlDict -entry deltaT -value)
-    latestTime=$(foamListTimes -case $src -noZero -latestTime -processor)
-
-    # Restart using steady results as first deltaT interval
-    echo "deltaT=$deltaT latestTime=$latestTime"
-    if [ -n "$latestTime" -a "$deltaT" != "$latestTime" ]
-    then
-    (
-        cd "$src" || exit
-
-        for proc in processor*
-        do
-            if [ -d "$proc/$latestTime" -a ! -d "$proc/$deltaT" ]
-            then
-                mv $proc/$latestTime $proc/$deltaT
-                rm -rf "$proc/$deltaT/uniform"
-            fi
-        done
-    )
-    fi
-
-    return 0
-}
-
-
 # Do steady-state case
 (cd steady && foamRunTutorials)
 
 if notTest "$@"
 then
-    latestTime=$(\cd steady && foamListTimes -noZero -latestTime -processor)
+    if canCompile
+    then
+        (cd code && wmake)
+    else
+        exit 0
+    fi
 
-    # Clone the steady-state case to transient
-    cloneParallelCase steady transient 0 $latestTime
+    . files/RunFunctions
 
-    copyParallelPointDisplacement transient $latestTime
+    caseName="transient"
+
+    latestTime=$(foamListTimes -case steady -noZero -latestTime -processor)
+
+    # Clone steady-state case to transient
+    cloneParallelCase steady "$caseName" 0 "$latestTime"
+
+    copyParallelPointDisplacement "$caseName" "$latestTime"
 
     # Adjust application (from simpleFoam -> pimpleFoam)
-    foamDictionary transient/system/controlDict \
+    foamDictionary "$caseName"/system/controlDict \
         -entry application -set pimpleFoam
 
-    # Do the transient case
-    \cp files/Allrun.transient transient/Allrun
-    (\cd transient && foamRunTutorials)
+    # Copy/link support files
+    linkFiles files "$caseName"
+
+    # Run
+    "$caseName/Allrun.$caseName" $*
 fi
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/lumpedPointMotion/building/Allrun.move b/tutorials/incompressible/lumpedPointMotion/building/Allrun.move
index 6c0c6a4f08..727311b2ed 100755
--- a/tutorials/incompressible/lumpedPointMotion/building/Allrun.move
+++ b/tutorials/incompressible/lumpedPointMotion/building/Allrun.move
@@ -4,62 +4,32 @@ cd "${0%/*}" || exit                                # Run from this directory
 #------------------------------------------------------------------------------
 
 # 1) Run meshing
-# 2) Reconstruct
-# 3) Test input zones and movement
+# 2) Test input zones and movement
 
-#
-# linkParallelCase srcDir dstDir
-#
-linkParallelCase()
-{
-    local src=$1
-    local dst=$2
-    shift 2
-
-    if [ -e "$dst" ]
-    then
-        echo "Case already linked: remove case directory $dst prior to linking"
-        return 1
-    elif [ ! -d "$src" ]
-    then
-        echo "Error: no directory to link: $src"
-        return 1
-    fi
-
-    echo "Linking $dst parallel case from $src"
-    mkdir $dst
-
-    # Copy system - may wish to change things
-    for i in system 0
-    do
-        echo "    copy $i/"
-        ( cd $dst && cp -r ../$src/$i . )
-    done
-
-    echo "    link constant/"
-    ( cd $dst && ln -sf ../$src/constant . )
-
-    echo "    link processor*/ with $# times: $@"
-    for proc in $(cd $src && \ls -d processor*)
-    do
-        ( cd $dst && ln -sf ../$src/$proc . )
-    done
-
-    return 0
-}
-
-
-# Do steady-state case
+# Meshing
 (cd steady && ./Allrun.pre)
 
-if notTest "$@"
+if notTest
 then
-    # Copy/link the steady-state case to movement
-    linkParallelCase steady movement
+    if canCompile
+    then
+        (cd code && wmake)
+    else
+        exit 0
+    fi
 
-    # Test movement
-    \cp files/Allrun.movement movement/Allrun
-    (cd movement && foamRunTutorials)
+    . files/RunFunctions
+
+    caseName="movement"
+
+    # Copy/link the steady-state case to movement
+    linkParallelCase steady "$caseName"
+
+    # Copy/link support files
+    linkFiles files "$caseName"
+
+    # Run
+    "$caseName/Allrun.$caseName" $*
 fi
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/lumpedPointMotion/building/code/Allclean b/tutorials/incompressible/lumpedPointMotion/building/code/Allclean
new file mode 100755
index 0000000000..95f32868e3
--- /dev/null
+++ b/tutorials/incompressible/lumpedPointMotion/building/code/Allclean
@@ -0,0 +1,13 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+#------------------------------------------------------------------------------
+
+wclean
+
+# Remove executable
+rm -f building-motion
+
+# Remove known output/debug files
+rm -f *.txt *.vtp *.vtp.series
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/lumpedPointMotion/building/code/Make/files b/tutorials/incompressible/lumpedPointMotion/building/code/Make/files
new file mode 100644
index 0000000000..fc78c4c5fd
--- /dev/null
+++ b/tutorials/incompressible/lumpedPointMotion/building/code/Make/files
@@ -0,0 +1,3 @@
+building-motion.C
+
+EXE = $(PWD)/building-motion
diff --git a/tutorials/incompressible/lumpedPointMotion/building/code/Make/options b/tutorials/incompressible/lumpedPointMotion/building/code/Make/options
new file mode 100644
index 0000000000..6211233098
--- /dev/null
+++ b/tutorials/incompressible/lumpedPointMotion/building/code/Make/options
@@ -0,0 +1,13 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/lumpedPointMotion/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfileFormats \
+    -lmeshTools \
+    -ldynamicMesh \
+    -llumpedPointMotion
diff --git a/tutorials/incompressible/lumpedPointMotion/building/code/building-motion.C b/tutorials/incompressible/lumpedPointMotion/building/code/building-motion.C
new file mode 100644
index 0000000000..14de04b9ed
--- /dev/null
+++ b/tutorials/incompressible/lumpedPointMotion/building/code/building-motion.C
@@ -0,0 +1,682 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    building-motion
+
+Description
+    Forced oscillation code for fluid-structure interface check.
+    Generates position and rotation angle of each node.
+    The top displacements of the target building calculated as sin() function.
+    Horizontal displacements of each node interpolated to the vertical
+    direction according to cantilever beam theory.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "Fstream.H"
+#include "unitConversion.H"
+#include "foamVtkSeriesWriter.H"
+#include "lumpedPointTools.H"
+#include "lumpedPointState.H"
+#include "lumpedPointIOMovement.H"
+
+using namespace Foam;
+
+//- Oscillator generator
+class position_generator
+{
+    // Private Member Functions
+
+        //- Calculate position/rotation at given time
+        lumpedPointState calc(scalar currTime) const
+        {
+            // Point positions
+            pointField points_(nDivisions+1, Zero);
+
+            // Point rotations
+            vectorField angles_(nDivisions+1, Zero);
+
+            // Set node heights (z)
+            forAll(points_, divi)
+            {
+                points_[divi].z() = (height * divi)/scalar(nDivisions);
+            }
+
+            const vector sines
+            (
+                Foam::sin(2*constant::mathematical::pi * currTime/period.x()),
+                Foam::sin(2*constant::mathematical::pi * currTime/period.y()),
+                Foam::sin(2*constant::mathematical::pi * currTime/period.z())
+            );
+
+            for (label divi = 1; divi <= nDivisions; ++divi)
+            {
+                const scalar zpos = points_[divi].z();
+                const scalar height1 = (height - zpos);
+
+                const scalar pos_factor =
+                (
+                    1.0/3.0 / pow4(height)
+                  * (
+                        3*pow4(height)
+                      - 4*pow3(height)*(height1)
+                      + pow4(height1)
+                    )
+                );
+
+                const scalar ang_factor =
+                (
+                    1.0/3.0 / pow4(height)
+                  * (
+                        4*pow3(height)
+                      - 4*pow3(height1)
+                    )
+                );
+
+                vector here
+                (
+                    (amplitude.x() * sines.x() * pos_factor),
+                    (amplitude.y() * sines.y() * pos_factor),
+                    zpos // Z position is invariant
+                );
+
+                vector rot
+                (
+                    Foam::atan(amplitude.x() * sines.x() * ang_factor),
+                    Foam::atan(amplitude.y() * sines.y() * ang_factor),
+                    Foam::atan(amplitude.z() * sines.z() * ang_factor)
+                );
+
+                // Assign
+                points_[divi] = here;
+
+                // The x<->y swap is intentional
+                angles_[divi] = vector{rot.y(), rot.x(), rot.z()};
+            }
+
+            // Return as lumpedPoint state
+            return lumpedPointState{points_, angles_};
+        }
+
+
+public:
+
+    // Control parameters
+
+        // The number of oscillating nodes
+        label nDivisions = 10;
+
+        // Height of target building [m]
+        scalar height = 0.5;
+
+        // Proper period (sec)
+        vector period = vector{1, 0.5, 1};
+
+        // Amplitude
+        vector amplitude = vector{0.03, 0.05, 0.3};
+
+
+    // Constructors
+
+        //- Default construct
+        position_generator() = default;
+
+
+    // Member Functions
+
+        //- Calculate position/rotation at given time
+        lumpedPointState state(const scalar currTime) const
+        {
+            return calc(currTime);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Forced oscillation code for fluid-structure interface check."
+        " Generates position and rotation angle of each node."
+        " The top displacements of the target building calculated as sin()"
+        " function."
+        " Horizontal displacements of each node interpolated to the vertical"
+        " direction according to cantilever beam theory."
+     );
+
+    argList::noBanner();
+    argList::noParallel();
+
+    // Geometric
+    argList::addOption
+    (
+        "nodes",
+        "N",
+        "The number of oscillating nodes (default: 10)"
+    );
+    argList::addOption
+    (
+        "height",
+        "value",
+        "Height of target building (m)"
+    );
+    argList::addOption
+    (
+        "period",
+        "(time time time)",
+        "The proper period (sec)"
+    );
+    argList::addOption
+    (
+        "amplitude",
+        "(value value value)",
+        "The amplitude"
+    );
+
+    // Time control
+    argList::addOption
+    (
+        "time",
+        "value",
+        "The time to use"
+    );
+    argList::addOption
+    (
+        "deltaT",
+        "value",
+        "The time increment for multiple time loops"
+    );
+    argList::addOption
+    (
+        "nTimes",
+        "value",
+        "The number of time loops"
+    );
+
+    // Query, output
+    argList::addBoolOption
+    (
+        "query",
+        "Report values only and exit"
+    );
+    argList::addOption
+    (
+        "output",
+        "file",
+        "write to file, with header"
+    );
+
+    argList::addOption
+    (
+        "scale",
+        "factor",
+        "Scaling factor for movement (default: 1)"
+    );
+    argList::addOption
+    (
+        "visual-length",
+        "len",
+        "Visualization length for planes (visualized as triangles)"
+    );
+
+    // Run controls
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Test movement without a mesh"
+    );
+    argList::addBoolOption
+    (
+        "removeLock",
+        "Remove lock-file on termination of slave"
+    );
+    argList::addBoolOption
+    (
+        "slave",
+        "Invoke as a slave responder for testing"
+    );
+    #include "setRootCase.H"
+
+
+    // The oscillator
+    position_generator gen;
+    args.readIfPresent("nodes", gen.nDivisions);
+    args.readIfPresent("height", gen.height);
+    args.readIfPresent("period", gen.period);
+    args.readIfPresent("amplitude", gen.amplitude);
+
+
+    // Control parameters
+    const bool dryrun = args.found("dry-run");
+    const bool slave = args.found("slave");
+    const bool removeLock = args.found("removeLock");
+
+    const bool optQuery = args.found("query");
+    const fileName outputFile(args.getOrDefault<fileName>("output", ""));
+
+    const scalar relax = args.getOrDefault<scalar>("scale", 1);
+
+    args.readIfPresent("visual-length", lumpedPointState::visLength);
+
+    // Time parameters
+    scalar currTime = args.getOrDefault<scalar>("time", 0);
+    const scalar deltaT = args.getOrDefault("deltaT", 0.001);
+    const label nTimes = args.getOrDefault<label>("nTimes", 1);
+
+    // Loop handling for slave
+    const bool infiniteLoop = slave && !args.found("nTimes");
+
+
+    // ----------------------------------------------------------------------
+    // Slave mode
+    // ----------------------------------------------------------------------
+
+    if (slave)
+    {
+        Info<< "Running as slave responder" << endl;
+
+        if (Pstream::parRun())
+        {
+            FatalErrorInFunction
+                << "Running as slave responder is not permitted in parallel"
+                << nl
+                << exit(FatalError);
+        }
+
+        #include "createTime.H"
+
+        // Create movement without a mesh
+        autoPtr<lumpedPointIOMovement> movementPtr =
+            lumpedPointIOMovement::New(runTime);
+
+        if (!movementPtr)
+        {
+            Info<< "No valid movement found" << endl;
+            return 1;
+        }
+        auto& movement = *movementPtr;
+
+        // Reference state0
+        const lumpedPointState& state0 = movement.state0();
+
+        externalFileCoupler& coupler = movement.coupler();
+
+        for (label timei = 0; infiniteLoop || (timei < nTimes); ++timei)
+        {
+            Info<< args.executable() << ": waiting for master" << endl;
+
+            // Wait for master, but stop if status=done was seen
+            if (!coupler.waitForMaster())
+            {
+                Info<< args.executable()
+                    << ": stopping status=done was detected" << endl;
+                break;
+            }
+
+            scalar timeValue = currTime;
+
+            if (infiniteLoop)
+            {
+                // Get output file
+                IFstream is(coupler.resolveFile(movement.outputName()));
+
+                dictionary dict;
+                is >> dict;
+
+                timeValue = dict.get<scalar>("time");
+            }
+
+            lumpedPointState state(gen.state(timeValue));
+            state.relax(relax, state0);
+
+            // Generate input for OpenFOAM
+            {
+                OFstream os(coupler.resolveFile(movement.inputName()));
+                if
+                (
+                    movement.inputFormat()
+                 == lumpedPointState::inputFormatType::PLAIN
+                )
+                {
+                    state.writePlain(os);
+                }
+                else
+                {
+                    os.writeEntry("time", timeValue);
+                    state.writeDict(os);
+                }
+            }
+
+            Info<< args.executable()
+                << ": updating state " << timei
+                << " - switch to master"
+                << endl;
+
+            // Let OpenFOAM know that it can continue
+            coupler.useMaster();
+
+            currTime += deltaT;
+        }
+
+        if (removeLock)
+        {
+            Info<< args.executable() << ": removing lock file" << endl;
+            coupler.useSlave();  // This removes the lock-file
+        }
+
+        Info<< args.executable() << ": finishing" << nl;
+
+        Info<< "\nEnd\n" << endl;
+        return 0;
+    }
+
+
+    // ----------------------------------------------------------------------
+    // dry-run
+    // ----------------------------------------------------------------------
+
+    if (dryrun)
+    {
+        Info<< "dry-run: creating states only" << nl;
+
+        autoPtr<Time> runTimePtr;
+        autoPtr<lumpedPointIOMovement> movementPtr;
+
+        const bool throwingIOError = FatalIOError.throwExceptions();
+        const bool throwingError = FatalError.throwExceptions();
+
+        try
+        {
+            Info<< "Create time" << flush;
+
+            runTimePtr = Time::New(args);
+
+            // Create movement without a mesh
+            movementPtr = lumpedPointIOMovement::New(*runTimePtr);
+        }
+        catch (...)
+        {
+            Info<< " ... failed (optional for dry-run)";
+        }
+        Info<< nl << endl;
+
+        FatalError.throwExceptions(throwingError);
+        FatalIOError.throwExceptions(throwingIOError);
+
+        if (!movementPtr)
+        {
+            Info<< "No time, run without movement information\n" << endl;
+        }
+
+        const lumpedPointState state0(gen.state(0));
+
+        vtk::seriesWriter stateSeries;
+
+        for
+        (
+            label timei = 0, outputCount = 0;
+            timei < nTimes;
+            ++timei
+        )
+        {
+            lumpedPointState state(gen.state(currTime));
+            state.relax(relax, state0);
+
+            Info<< "output [" << timei << '/' << nTimes << ']';
+
+            // State/response = what comes back from FEM
+            {
+                const word outputName =
+                    word::printf("state_%06d.vtp", outputCount);
+
+                Info<< "  " << outputName;
+
+                if (movementPtr)
+                {
+                    movementPtr->writeStateVTP(state, outputName);
+                }
+                else
+                {
+                    state.writeVTP(outputName);
+                }
+                stateSeries.append(outputCount, outputName);
+            }
+
+            Info<< endl;
+
+            ++outputCount;
+            currTime += deltaT;
+        }
+
+        // Write file series
+        if (stateSeries.size())
+        {
+            Info<< nl << "write state.vtp.series" << nl;
+            stateSeries.write("state.vtp");
+        }
+
+        Info<< "\nEnd\n" << endl;
+        return 0;
+    }
+
+
+    // ----------------------------------------------------------------------
+    // Report values or generate a file
+    // ----------------------------------------------------------------------
+
+    if (optQuery || !outputFile.empty())
+    {
+        autoPtr<OFstream> osPtr;
+
+        if (!outputFile.empty())
+        {
+            osPtr.reset(new OFstream(outputFile));
+            auto& os = *osPtr;
+
+            os.precision(8);
+
+            // One file with everything, output using OpenFOAM syntax
+
+            IOobject::writeBanner(os)
+                << "FoamFile\n{\n"
+                << "    version     " << os.version() << ";\n"
+                << "    format      " << os.format() << ";\n"
+                << "    class       " << "dictionary" << ";\n"
+                << "    object      " << "response" << ";\n"
+                << "}\n";
+
+            IOobject::writeDivider(os) << nl;
+
+            os << "// angles are Euler angles z-x-z (intrinsic)" << nl;
+            os.writeEntry("degrees", "false");
+            os << nl;
+            os << "response" << nl;
+            os << '(' << nl;
+        }
+        else
+        {
+            Info().precision(8);
+        }
+
+
+        for (label timei = 0; timei < nTimes; ++timei)
+        {
+            lumpedPointState state(gen.state(currTime));
+
+            if (osPtr)
+            {
+                // Report position/angle
+                auto& os = *osPtr;
+
+                os.beginBlock();
+
+                os.writeEntry("time", currTime);
+
+                state.writeDict(os);
+
+                os.endBlock();
+            }
+            else
+            {
+                // Report position/angle
+                auto& os = Info();
+
+                os.writeEntry("time", currTime);
+                state.writeDict(os);
+            }
+
+            currTime += deltaT;
+        }
+
+
+        if (osPtr)
+        {
+            auto& os = *osPtr;
+
+            os << ')' << token::END_STATEMENT << nl;
+
+            IOobject::writeEndDivider(os);
+
+            Info<< "\nEnd\n" << endl;
+        }
+
+        return 0;
+    }
+
+
+    // ----------------------------------------------------------------------
+    // test patch movement
+    // ----------------------------------------------------------------------
+
+    #include "createTime.H"
+
+    runTime.setTime(instant(runTime.constant()), 0);
+
+    #include "createNamedMesh.H"
+
+    // Create movement with mesh
+    autoPtr<lumpedPointIOMovement> movementPtr =
+        lumpedPointIOMovement::New(mesh);
+
+    if (!movementPtr)
+    {
+        Info<< "No valid movement found" << endl;
+        return 1;
+    }
+    auto& movement = *movementPtr;
+
+    // Reference state0
+    const lumpedPointState& state0 = movement.state0();
+
+    pointIOField points0(lumpedPointTools::points0Field(mesh));
+
+    const label nPatches = lumpedPointTools::setPatchControls(mesh, points0);
+    if (!nPatches)
+    {
+        Info<< "No point patches with lumped movement found" << endl;
+        return 2;
+    }
+
+    Info<< "Lumped point patch controls set on "
+        << nPatches << " patches" << nl;
+
+    lumpedPointTools::setInterpolators(mesh, points0);
+
+
+    // Output vtk file series
+    vtk::seriesWriter stateSeries;
+    vtk::seriesWriter geomSeries;
+
+    // Initial geometry
+    movement.writeVTP("geom_init.vtp", state0, mesh, points0);
+
+    lumpedPointTools::setInterpolators(mesh);
+
+    for
+    (
+        label timei = 0, outputCount = 0;
+        timei < nTimes;
+        ++timei
+    )
+    {
+        lumpedPointState state(gen.state(currTime));
+
+        state += movement.origin();
+        movement.scalePoints(state);
+        state.relax(relax, state0);
+
+        Info<< "output [" << timei << '/' << nTimes << ']';
+
+        // State/response = what comes back from FEM
+        {
+            const word outputName =
+                word::printf("state_%06d.vtp", outputCount);
+
+            Info<< "  " << outputName;
+
+            movement.writeStateVTP(state, outputName);
+            stateSeries.append(outputCount, outputName);
+        }
+
+        {
+            const word outputName =
+                word::printf("geom_%06d.vtp", outputCount);
+
+            Info<< "  " << outputName;
+
+            movement.writeVTP(outputName, state, mesh, points0);
+            geomSeries.append(outputCount, outputName);
+        }
+
+        Info<< endl;
+
+        ++outputCount;
+        currTime += deltaT;
+    }
+
+
+    // Write file series
+
+    if (geomSeries.size())
+    {
+        Info<< nl << "write geom.vtp.series" << nl;
+        geomSeries.write("geom.vtp");
+    }
+    if (stateSeries.size())
+    {
+        Info<< nl << "write state.vtp.series" << nl;
+        stateSeries.write("state.vtp");
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.moveMesh b/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.moveMesh
new file mode 100755
index 0000000000..268ba7aa20
--- /dev/null
+++ b/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.moveMesh
@@ -0,0 +1,34 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+# Cleanup old junk that may prevent things from starting
+rm -f comms/OpenFOAM.lock
+
+# If OpenFOAM stops prematurely, trigger the external solver to stop
+trap '[ -e comms/OpenFOAM.lock ] && echo "status=done" > comms/OpenFOAM.lock' EXIT TERM INT
+
+# Simulated external solver.
+# Using -scale 0.01 since input movement table included visual scaling
+
+if false
+then
+    # Create response file
+    runApplication -overwrite \
+        ../code/building-motion -deltaT 0.001 -nTimes 5001 -output response.txt
+
+    # Use response file for states
+    runApplication -overwrite \
+        lumpedPointMovement -scale 0.01 -removeLock -slave response.txt &
+else
+
+    # Generate states on demand
+    runApplication -overwrite \
+        ../code/building-motion -scale 0.01 -removeLock -slave &
+fi
+
+# Run moveMesh with deltaT corresponding to dynamicMeshDict updateInterval
+runParallel moveMesh -deltaT 0.001
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.movement b/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.movement
index c6133ada33..7cf6286b09 100755
--- a/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.movement
+++ b/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.movement
@@ -12,7 +12,22 @@ runApplication reconstructParMesh -constant -withZero -time 0
 runApplication lumpedPointZones
 
 # Simulated external solver
-# Using -scale=1 to see the excessively large movements
-runApplication lumpedPointMovement -span 25 -scale 1 ../files/response.txt
+# Use -scale=1 to see the excessively large movements
+
+if false
+then
+    # Create response file
+    runApplication -overwrite \
+        ../code/building-motion -deltaT 0.001 -nTimes 5001 -output response.txt
+
+    # Use response file for states
+    runApplication -overwrite \
+        lumpedPointMovement -span 25 -scale 1 response.txt
+else
+
+    # Generate states on demand
+    runApplication -overwrite \
+        ../code/building-motion -scale 1 -deltaT 0.025 -nTimes 201
+fi
 
 #------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.transient b/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.transient
index 5138ae30ca..998f759d90 100755
--- a/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.transient
+++ b/tutorials/incompressible/lumpedPointMotion/building/files/Allrun.transient
@@ -10,9 +10,23 @@ rm -f comms/OpenFOAM.lock
 trap '[ -e comms/OpenFOAM.lock ] && echo "status=done" > comms/OpenFOAM.lock' EXIT TERM INT
 
 # Simulated external solver.
-# Using -scale since the input movement table is excessively large
-runApplication -overwrite \
-    lumpedPointMovement -scale 0.01 -removeLock -slave ../files/response.txt &
+# Using -scale 0.01 since input movement table included visual scaling
+
+if false
+then
+    # Create response file
+    runApplication -overwrite \
+        ../code/building-motion -deltaT 0.001 -nTimes 5001 -output response.txt
+
+    # Use response file for states
+    runApplication -overwrite \
+        lumpedPointMovement -scale 0.01 -removeLock -slave response.txt &
+else
+
+    # Generate states on demand
+    runApplication -overwrite \
+        ../code/building-motion -scale 0.01 -removeLock -slave &
+fi
 
 # Run OpenFOAM
 runParallel $(getApplication)
diff --git a/tutorials/incompressible/lumpedPointMotion/building/files/RunFunctions b/tutorials/incompressible/lumpedPointMotion/building/files/RunFunctions
new file mode 100644
index 0000000000..9cd3880c93
--- /dev/null
+++ b/tutorials/incompressible/lumpedPointMotion/building/files/RunFunctions
@@ -0,0 +1,147 @@
+#---------------------------------*- sh -*-------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | www.openfoam.com
+#    \\/     M anipulation  |
+#------------------------------------------------------------------------------
+#     Copyright (C) 2020 OpenCFD Ltd.
+#------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+#
+# Script
+#     RunFunctions
+#
+# Description
+#     Additional functions for copy/linking FSI cases
+#
+#------------------------------------------------------------------------------
+
+#
+# copyParallelPointDisplacement caseDir timeName
+#
+# Copy pointDisplacement from caseDir/0/ to caseDir/timeName/
+#
+copyParallelPointDisplacement()
+{
+    local src="$1"
+    local dstTime="$2"
+    local file=pointDisplacement
+
+    [ -d "$src" ] || {
+        echo "Error: no directory: $src"
+        return 1
+    }
+
+    # Copy select directories
+    echo "    copy processor '$file' from 0/ -> $dstTime"
+    if [ -n "$dstTime" ]
+    then
+    (
+        cd "$src" || exit
+
+        for proc in processor*
+        do
+            if [ -d "$proc/0" ] && [ -d "$proc/$dstTime" ]
+            then
+                cp "$proc/0/$file" "$proc/$dstTime/$file"
+                if [ -d "$proc/0/include" ]
+                then
+                    cp -r "$proc/0/include" "$proc/$dstTime"
+                fi
+            fi
+        done
+    )
+    else
+        echo "    no destination time"
+    fi
+
+    # Restart from latestTime
+    foamDictionary "$src"/system/controlDict \
+        -entry startFrom -set latestTime
+
+    deltaT=$(foamDictionary "$src"/system/controlDict -entry deltaT -value)
+    latestTime=$(foamListTimes -case "$src" -noZero -latestTime -processor)
+
+    # Restart using steady results as first deltaT interval
+    echo "deltaT=$deltaT latestTime=$latestTime"
+    if [ -n "$latestTime" ] && [ "$deltaT" != "$latestTime" ]
+    then
+    (
+        cd "$src" || exit
+
+        for proc in processor*
+        do
+            if [ -d "$proc/$latestTime" ] && [ ! -d "$proc/$deltaT" ]
+            then
+                mv "$proc/$latestTime" "$proc/$deltaT"
+                rm -rf "$proc/$deltaT/uniform"
+            fi
+        done
+    )
+    fi
+
+    return 0
+}
+
+#
+# linkParallelCase srcDir dstDir
+#
+linkParallelCase()
+{
+    local src="$1"
+    local dst="$2"
+    shift 2
+
+    if [ -e "$dst" ]
+    then
+        echo "Case already linked: remove case directory $dst prior to linking"
+        return 1
+    elif [ ! -d "$src" ]
+    then
+        echo "Error: no directory to link: $src"
+        return 1
+    fi
+
+    echo "Linking $dst parallel case from $src"
+    mkdir -p "$dst"
+
+    # Copy system - may wish to change things
+    for i in system 0
+    do
+        echo "    copy $i/"
+        ( cd "$dst" && cp -r "../$src/$i" . )
+    done
+
+    echo "    link constant/"
+    ( cd "$dst" && ln -sf "../$src/constant" . )
+
+    echo "    link processor*/ with $# times: $@"
+    for proc in $(cd "$src" && \ls -d processor*)
+    do
+        ( cd "$dst" && ln -sf "../$src/$proc" . )
+    done
+
+    return 0
+}
+
+#
+# linkFiles srcDir dstDir
+#
+linkFiles()
+{
+    local src="$1"
+    local dst="$2"
+    shift
+
+    echo "Linking $dst control files from $src"
+    mkdir -p "$dst"
+
+    ( cd "$dst" && ln -sf ../"$src"/* . )
+
+    return 0
+}
+
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/lumpedPointMotion/building/files/response.txt.gz b/tutorials/incompressible/lumpedPointMotion/building/files/response.txt.gz
deleted file mode 100644
index 30b6604355..0000000000
Binary files a/tutorials/incompressible/lumpedPointMotion/building/files/response.txt.gz and /dev/null differ
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/0.orig/pointDisplacement b/tutorials/incompressible/lumpedPointMotion/building/steady/0.orig/pointDisplacement
index 931c8e7918..05e9f62c96 100644
--- a/tutorials/incompressible/lumpedPointMotion/building/steady/0.orig/pointDisplacement
+++ b/tutorials/incompressible/lumpedPointMotion/building/steady/0.orig/pointDisplacement
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v1912                                 |
+|  \\    /   O peration     | Version:  v2006                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -67,6 +67,7 @@ boundaryField
     {
         type            lumpedPointDisplacement;
         value           uniform (0 0 0);
+        controllers     ( vertical );
     }
 }
 
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/Allclean b/tutorials/incompressible/lumpedPointMotion/building/steady/Allclean
index 4b409fd7cd..143cbed889 100755
--- a/tutorials/incompressible/lumpedPointMotion/building/steady/Allclean
+++ b/tutorials/incompressible/lumpedPointMotion/building/steady/Allclean
@@ -8,4 +8,7 @@ cleanCase0
 # Clean up copied/derived files
 rm -rf constant/triSurface
 
+# Clean up debug/setup files
+rm -f *.txt *.vtp *.vtp.series
+
 #------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun b/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun
index 0ef141508a..57b865a1e7 100755
--- a/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun
+++ b/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun
@@ -21,5 +21,4 @@ else
 
 fi
 
-
 #------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun.init b/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun.init
new file mode 100755
index 0000000000..d3e302def2
--- /dev/null
+++ b/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun.init
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+# Copy geometry from resources directory
+
+mkdir -p constant/triSurface/
+cp "$FOAM_TUTORIALS"/resources/geometry/building_wtc2.obj constant/triSurface/
+
+#------------------------------------------------------------------------------
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun.pre b/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun.pre
index 26f4dff042..60a531539e 100755
--- a/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun.pre
+++ b/tutorials/incompressible/lumpedPointMotion/building/steady/Allrun.pre
@@ -3,9 +3,8 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
-# Copy building from resources directory
-mkdir -p constant/triSurface/
-cp $FOAM_TUTORIALS/resources/geometry/building_wtc2.obj constant/triSurface/
+# Get geometry and other resources
+./Allrun.init
 
 # runApplication surfaceFeatureExtract
 runApplication blockMesh
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/constant/dynamicMeshDict b/tutorials/incompressible/lumpedPointMotion/building/steady/constant/dynamicMeshDict
index c306e18ca9..4d57bcf202 100644
--- a/tutorials/incompressible/lumpedPointMotion/building/steady/constant/dynamicMeshDict
+++ b/tutorials/incompressible/lumpedPointMotion/building/steady/constant/dynamicMeshDict
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v1912                                 |
+|  \\    /   O peration     | Version:  v2006                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -13,9 +13,13 @@ FoamFile
     location    "system";
     object      dynamicMeshDict;
 }
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Less frequent mesh motion
+updateControl   runTime;
+
+updateInterval  0.001;
+
 dynamicFvMesh   dynamicMotionSolverFvMesh;
 
 motionSolverLibs (fvMotionSolvers);
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointControllers b/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointControllers
new file mode 100644
index 0000000000..49f7f73bf9
--- /dev/null
+++ b/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointControllers
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "locations and connectivity";
+    object      dictionary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Locations of the lumped points
+points
+(
+    (0  0  0.00)
+    (0  0  0.05)
+    (0  0  0.10)
+    (0  0  0.15)
+    (0  0  0.20)
+    (0  0  0.25)
+    (0  0  0.30)
+    (0  0  0.35)
+    (0  0  0.40)
+    (0  0  0.45)
+    (0  0  0.50)
+);
+
+//- Connectivity for motion controllers
+controllers
+{
+    vertical
+    {
+        pointLabels (0 1 2 3 4 5 6 7 8 9 10);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointMovement b/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointMovement
index c3725c8fbf..e968b1da09 100644
--- a/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointMovement
+++ b/tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointMovement
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v1912                                 |
+|  \\    /   O peration     | Version:  v2006                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -15,36 +15,24 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Reference axis for the locations
-axis            (0 0 1);
+#include "lumpedPointControllers"
 
-// Locations of the lumped points
-locations       11(0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5);
-
-// Division for pressure forces (0-1)
-division        0.5;
-
-//- If present, the offset of patch points compared to the locations
-//  Otherwise determined from the bounding box
-// centre          (0 0 0);
-
-//- The interpolation scheme
-interpolationScheme linear;
+//- Offset/shift for lumped points
+origin  (0.3 0.15 0);
 
 //- Relaxation/scaling factor when updating positions
-relax           1.0;
-
+relax   1.0;
 
 forces
 {
     //- The pressure name (default: p)
-    p           p;
+    p       p;
 
     //- Reference pressure [Pa] (default: 0)
-    pRef        0;
+    pRef    0;
 
     //- Reference density for incompressible calculations (default: 1)
-    rhoRef      1;
+    rhoRef  1;
 }
 
 
@@ -56,6 +44,9 @@ communication
 
     waitInterval    1;
 
+    // Coupling frequency (default: 1)
+    //calcFrequency   1;
+
     timeOut         100;
 
     initByExternal  false;
@@ -78,20 +69,20 @@ communication
     scaleInput
     {
         //- Length multiplier (to metres). Eg 0.001 for [mm] -> [m]
-        length      1;
+        length  1;
     }
 
     // Scaling applied to values written to 'outputName'
     scaleOutput
     {
         //- Length multiplier (from metres). Eg 1000 for [m] -> [mm]
-        length      1;
+        length  1;
 
         //- Force units multiplier (from Pa)
-        force       1;
+        force   1;
 
         //- Moment units multiplier (from N.m)
-        moment      1;
+        moment  1;
     }
 }