ENH: integration of lumpedPointMotion

- This provides a mechanism for moving mesh patches based on external
  input (eg, from an external structures solver). The patch points are
  influenced by the position and rotation of the lumped points.

  BC:  lumpedPointDisplacementPointPatchVectorField

  Controlling mechanisms:
  - externalCoupler
    for coordinating the master/slave

  - lumpedPointMovement
    manages the patch-points motion, but also for extracting forces/moments

  - lumpedPointState
    represents the positions/rotations of the controlling points

  Utils:
  - lumpedPointZones
    diagnostic for visualizing the correspondence between controlling
    points and patch faces

  - lumpedPointMovement
    Test that the patch motion is as desired without invoking moveMesh.
    With the -slave option, return items from a precalculated table
    for the lumpedPointDisplacementPointPatchVectorField BC.

applications/test/externalCoupler/Make/files

@@ -0,0 +1,3 @@
+Test-externalCoupler.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-externalCoupler

applications/test/externalCoupler/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/lumpedPointMotion/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -llumpedPointMotion

applications/test/externalCoupler/Test-externalCoupler.C

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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
+    Test-externalCoupler
+
+Description
+    Test of master/slave communication etc.
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "externalCoupler.H"
+
+using namespace Foam;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    argList::addOption("max", "N", "max number of calls (default: 1000)");
+    argList::addBoolOption("slave", "run as slave");
+
+    #include "setRootCase.H"
+
+    const label maxCount = args.optionLookupOrDefault<label>("max", 1000);
+
+    externalCoupler coupler;
+
+    if (args.optionFound("slave"))
+    {
+        const word role = "slave";
+        Info<< "Running as " << role << " max=" << maxCount << endl;
+
+        for (label count = 0; count < maxCount; ++count)
+        {
+            // Wait for master, but stop if status=done was seen
+
+            Info<< role << ": waiting for master" << endl;
+            if (!coupler.waitForMaster())
+            {
+                Info<< role << ": stopping. status=done was detected" << endl;
+                break;
+            }
+
+            Info<< role << ": switch to master" << endl;
+            coupler.useMaster();
+        }
+    }
+    else
+    {
+        const word role = "master";
+        Info<< "Running as " << role << " max=" << maxCount << endl;
+
+        for (label count = 0; count < maxCount; ++count)
+        {
+            // Wait for slave, but stop if status=done was seen
+
+            Info<< role << ": waiting for slave" << endl;
+            if (!coupler.waitForSlave())
+            {
+                Info<< role << ": stopping. status=done was detected" << endl;
+                break;
+            }
+
+            Info<< role << ": switch to slave" << endl;
+            coupler.useSlave();
+        }
+    }
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/postProcessing/lumped/lumpedPointForces/Make/files

@@ -0,0 +1,3 @@
+lumpedPointForces.C
+
+EXE = $(FOAM_APPBIN)/lumpedPointForces

applications/utilities/postProcessing/lumped/lumpedPointForces/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lumpedPointMotion/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llumpedPointMotion

applications/utilities/postProcessing/lumped/lumpedPointForces/lumpedPointForces.C

@@ -0,0 +1,177 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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
+    lumpedPointForces
+
+Description
+    Extract force/moment information from existing calculations based
+    on the segmentation of the pressure integration zones.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "timeSelector.H"
+#include "volFields.H"
+#include "IOobjectList.H"
+
+#include "lumpedPointTools.H"
+#include "lumpedPointIOMovement.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class GeoFieldType>
+autoPtr<GeoFieldType> loadField
+(
+    const volMesh::Mesh& mesh,
+    const IOobject* io
+)
+{
+    if (io && io->headerClassName() == GeoFieldType::typeName)
+    {
+        Info<< "Reading " << GeoFieldType::typeName
+            << ' ' << io->name() << endl;
+
+        return autoPtr<GeoFieldType>
+        (
+            new GeoFieldType
+            (
+                IOobject
+                (
+                    io->name(),
+                    io->instance(),
+                    io->local(),
+                    io->db(),
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE,
+                    io->registerObject()
+                ),
+                mesh
+            )
+        );
+    }
+
+    return autoPtr<GeoFieldType>();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Extract force/moment information from existing calculations based "
+        "on the lumpedPoints pressure zones."
+    );
+
+    argList::addBoolOption
+    (
+        "vtk",
+        "create visualization files of the forces"
+    );
+
+    timeSelector::addOptions(true, false);
+    argList::noFunctionObjects();  // Never use function objects
+
+    #include "addRegionOption.H"
+    #include "setRootCase.H"
+
+    const bool withVTK = args.optionFound("vtk");
+
+    #include "createTime.H"
+    instantList timeDirs = timeSelector::select0(runTime, args);
+
+    #include "createNamedMesh.H"
+
+    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New
+    (
+        runTime
+    );
+
+    if (!movement.valid())
+    {
+        Info<< "no valid movement given" << endl;
+        return 1;
+    }
+
+    const labelList patchLst = lumpedPointTools::lumpedPointPatchList(mesh);
+    if (patchLst.empty())
+    {
+        Info<< "no patch list found" << endl;
+        return 2;
+    }
+
+    movement().setMapping(mesh, patchLst, lumpedPointTools::points0Field(mesh));
+
+    List<vector> forces, moments;
+    forAll(timeDirs, timei)
+    {
+        runTime.setTime(timeDirs[timei], timei);
+
+        Info<< "Time = " << runTime.timeName() << endl;
+
+        if (mesh.readUpdate())
+        {
+            Info<< "    Read new mesh" << nl;
+        }
+
+        // Search for list of objects for this time
+        IOobjectList objects(mesh, runTime.timeName());
+
+        // Pressure field
+        autoPtr<volScalarField> pField
+            = loadField<volScalarField>(mesh, objects.lookup("p"));
+
+        // The forces per zone
+        if (movement().forcesAndMoments(mesh, forces, moments))
+        {
+            Info<<"forces per zone: " << forces << endl;
+            Info<<"moments per zone: " << moments << endl;
+
+            if (withVTK && Pstream::master())
+            {
+                const word outputName =
+                    Foam::name("forces_%06d.vtp", runTime.timeIndex());
+
+                Info<<"    " << outputName << endl;
+
+                movement().writeForcesAndMomentsVTP
+                (
+                    outputName,
+                    forces,
+                    moments
+                );
+            }
+        }
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/utilities/postProcessing/lumped/lumpedPointMovement/Make/files

@@ -0,0 +1,3 @@
+lumpedPointMovement.C
+
+EXE = $(FOAM_APPBIN)/lumpedPointMovement

applications/utilities/postProcessing/lumped/lumpedPointMovement/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lumpedPointMotion/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llumpedPointMotion

applications/utilities/postProcessing/lumped/lumpedPointMovement/lumpedPointMovement.C

@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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
+    lumpedPointMovement
+
+Description
+    Thus utility can be used to produce VTK files to visualize the response
+    points/rotations and the corresponding movement of the building surfaces.
+    Uses the tabulated responses from the specified file.
+    Optionally, it can also be used to a dummy responder for the
+    externalCoupler logic, which makes it useful as a debugging facility
+    as well demonstrating how an external application could communicate
+    with the lumpedPointMovement point-patch boundary condition.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "timeSelector.H"
+
+#include "OFstream.H"
+
+#include "lumpedPointTools.H"
+#include "lumpedPointIOMovement.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Visualize lumpedPoint movements or provide a slave responder "
+        "for diagnostic purposes."
+    );
+
+    argList::noParallel();
+    argList::noFunctionObjects();  // Never use function objects
+    argList::addOption
+    (
+        "max",
+        "N",
+        "maximum number of outputs"
+    );
+    argList::addOption
+    (
+        "span",
+        "N",
+        "increment each input by factor N (default: 1)"
+    );
+    argList::addOption
+    (
+        "scale",
+        "factor",
+        "relaxation/scaling factor for movement (default: 1)"
+    );
+    argList::addBoolOption
+    (
+        "removeLock",
+        "remove lock-file on termination of slave"
+    );
+    argList::addBoolOption
+    (
+        "slave",
+        "invoke as a slave responder for testing"
+    );
+    argList::validArgs.insert("responseFile");
+
+    #include "setRootCase.H"
+
+    const label maxOut =
+        Foam::max(0, args.optionLookupOrDefault<label>("max", 0));
+
+    const label span =
+        Foam::max(1, args.optionLookupOrDefault<label>("span", 1));
+
+    const scalar relax = args.optionLookupOrDefault<scalar>("scale", 1);
+
+    const bool slave = args.optionFound("slave");
+    const bool removeLock = args.optionFound("removeLock");
+
+    #include "createTime.H"
+
+    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New
+    (
+        runTime
+    );
+
+    if (!movement.valid())
+    {
+        Info<< "no valid movement given" << endl;
+        return 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;
+    }
+
+    if (slave)
+    {
+        Info<< "Running as slave responder" << endl;
+
+        externalCoupler& coupler = movement().coupler();
+
+        label count = 0;
+        for (label index = 0; index < responseTable.size(); index += span)
+        {
+            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;
+            }
+
+            lumpedPointState state = responseTable[index].second();
+            state.relax(relax, movement().state0());
+
+            // Generate input for OpenFOAM
+            OFstream os(coupler.resolveFile(movement().inputName()));
+            if
+            (
+                movement().inputFormat()
+             == lumpedPointState::inputFormatType::PLAIN
+            )
+            {
+                state.writePlain(os);
+            }
+            else
+            {
+                os.writeEntry("time", responseTable[index].first());
+                state.writeDict(os);
+            }
+
+            Info<< args.executable()
+                << ": updated to state " << index
+                << " - switch to master"
+                << endl;
+
+            // Let OpenFOAM know that it can continue
+            coupler.useMaster();
+
+            if (maxOut && ++count >= maxOut)
+            {
+                Info<< args.executable()
+                    <<": stopping after " << maxOut << " outputs" << endl;
+                break;
+            }
+        }
+
+        if (removeLock)
+        {
+            Info<< args.executable() <<": removing lock file" << endl;
+            coupler.useSlave();  // This removes the lock-file
+        }
+    }
+    else
+    {
+        runTime.setTime(instant(0, runTime.constant()), 0);
+
+        #include "createNamedPolyMesh.H"
+
+        const labelList patchLst = lumpedPointTools::lumpedPointPatchList(mesh);
+        if (patchLst.empty())
+        {
+            Info<< "no patch list found" << endl;
+            return 2;
+        }
+
+        pointIOField points0 = lumpedPointTools::points0Field(mesh);
+        movement().setBoundBox(mesh, patchLst, points0);
+
+        label index = 0;
+
+        // Initial geometry
+        movement().writeVTP("geom_init.vtp", mesh, patchLst, points0);
+
+        forAll(responseTable, i)
+        {
+            const bool output = ((i % span) == 0);
+            lumpedPointState state = responseTable[i].second();
+            state.relax(relax, movement().state0());
+
+            if (output)
+            {
+                Info<<"output [" << i << "/"
+                    << responseTable.size() << "]" << endl;
+            }
+            else
+            {
+                continue;
+            }
+
+            // State/response = what comes back from FEM
+            {
+                const word outputName = Foam::name("state_%06d.vtp", index);
+                Info<<"    " << outputName << endl;
+
+                state.writeVTP(outputName, movement().axis());
+            }
+
+            {
+                const word outputName = Foam::name("geom_%06d.vtp", index);
+                Info<<"    " << outputName << endl;
+
+                movement().writeVTP(outputName, state, mesh, patchLst, points0);
+            }
+
+            {
+                ++index;
+
+                bool canOutput = !maxOut || (index <= maxOut);
+                if (!canOutput)
+                {
+                    Info<<"stopping output after "
+                        << maxOut << " outputs" << endl;
+                    break;
+                }
+            }
+        }
+
+    }
+
+    Info<< args.executable() << ": End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/postProcessing/lumped/lumpedPointZones/Make/files

@@ -0,0 +1,3 @@
+lumpedPointZones.C
+
+EXE = $(FOAM_APPBIN)/lumpedPointZones

applications/utilities/postProcessing/lumped/lumpedPointZones/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/lumpedPointMotion/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfiniteVolume \
+    -llumpedPointMotion

applications/utilities/postProcessing/lumped/lumpedPointZones/lumpedPointZones.C

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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
+    lumpedPointZones
+
+Description
+    Produces a VTK PolyData file \c lumpedPointZones.vtp in which the
+    segmentation of the pressure integration zones can be visualized
+    for diagnostic purposes. Does not use external coupling.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "timeSelector.H"
+
+#include "lumpedPointTools.H"
+#include "lumpedPointIOMovement.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Create lumpedPointZones.vtp to verify the segmentation of "
+        "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
+    (
+        "verbose",
+        "increased verbosity"
+    );
+
+    #include "addRegionOption.H"
+    #include "setRootCase.H"
+
+    // const bool verbose = args.optionFound("verbose");
+
+    #include "createTime.H"
+
+    runTime.setTime(instant(0, runTime.constant()), 0);
+
+    #include "createNamedPolyMesh.H"
+
+    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New
+    (
+        runTime
+    );
+
+    if (!movement.valid())
+    {
+        Info<< "no valid movement found" << endl;
+        return 1;
+    }
+
+    const labelList patchLst = lumpedPointTools::lumpedPointPatchList(mesh);
+    if (patchLst.empty())
+    {
+        Info<< "no patch list found" << endl;
+        return 2;
+    }
+
+    pointIOField points0 = lumpedPointTools::points0Field(mesh);
+    movement().setMapping(mesh, patchLst, points0);
+
+    // Initial geometry, but with zone colouring
+    movement().writeZonesVTP("lumpedPointZones.vtp", mesh, points0);
+
+    // Initial positions/rotations
+    movement().writeStateVTP("initialState.vtp");
+
+    Info<< nl
+        << "wrote 'lumpedPointZones.vtp'" << nl
+        << "wrote 'initialState.vtp'" << nl
+        << "End\n" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

src/Allwmake

@@ -76,6 +76,7 @@ wmake $targetType regionCoupled
 
 functionObjects/Allwmake $targetType $*
 
+wmake $targetType lumpedPointMotion
 wmake $targetType sixDoFRigidBodyMotion
 wmake $targetType rigidBodyDynamics
 wmake $targetType rigidBodyMeshMotion

src/lumpedPointMotion/Make/files

@@ -0,0 +1,13 @@
+externalCoupler.C
+
+lumpedPointMovement.C
+lumpedPointMovementWriter.C
+lumpedPointState.C
+lumpedPointStateWriter.C
+lumpedPointIOMovement.C
+
+lumpedPointDisplacementPointPatchVectorField.C
+
+lumpedPointTools.C
+
+LIB = $(FOAM_LIBBIN)/liblumpedPointMotion

src/lumpedPointMotion/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -ldynamicMesh \
+    -lmeshTools

src/lumpedPointMotion/externalCoupler.C

@@ -0,0 +1,309 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "externalCoupler.H"
+#include "Pstream.H"
+#include "PstreamReduceOps.H"
+#include "OSspecific.H"
+#include "IFstream.H"
+#include "OFstream.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(externalCoupler, 0);
+}
+
+Foam::word Foam::externalCoupler::lockName = "OpenFOAM";
+
+
+// file-scope
+// check file (must exist) for "status=done" content
+static bool checkIsDone(const std::string& lck)
+{
+    std::string content;
+    std::ifstream is(lck.c_str());
+    is >> content;
+
+    return (content.find("done") != std::string::npos);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+const Foam::fileName& Foam::externalCoupler::baseDir() const
+{
+    return commsDir_;
+}
+
+
+Foam::fileName Foam::externalCoupler::lockFile() const
+{
+    return resolveFile(lockName + ".lock");
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::externalCoupler::externalCoupler()
+:
+    runState_(NONE),
+    commsDir_("$FOAM_CASE/comms"),
+    waitInterval_(1u),
+    timeOut_(100u),
+    slaveFirst_(false),
+    log(false)
+{
+    commsDir_.expand();
+    commsDir_.clean();
+}
+
+
+Foam::externalCoupler::externalCoupler(const dictionary& dict)
+:
+    externalCoupler()
+{
+    readDict(dict);
+
+    if (Pstream::master())
+    {
+        mkDir(baseDir());
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::externalCoupler::~externalCoupler()
+{
+    shutdown();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::externalCoupler::readDict(const dictionary& dict)
+{
+    // Normally cannot change directory or initialization
+    // if things have already been initialized
+    dict.lookup("commsDir") >> commsDir_;
+    commsDir_.expand();
+    commsDir_.clean();
+    slaveFirst_ = dict.lookupOrDefault<bool>("initByExternal", false);
+
+    waitInterval_ = dict.lookupOrDefault("waitInterval", 1u);
+    if (!waitInterval_)
+    {
+        // Enforce non-zero sleep
+        waitInterval_ = 1u;
+    }
+
+    timeOut_ = dict.lookupOrDefault("timeOut", 100*waitInterval_);
+
+    log = dict.lookupOrDefault<bool>("log", false);
+
+    return true;
+}
+
+
+bool Foam::externalCoupler::initialized() const
+{
+    return runState_ != NONE;
+}
+
+
+bool Foam::externalCoupler::slaveFirst() const
+{
+    return slaveFirst_;
+}
+
+
+void Foam::externalCoupler::useMaster(const bool wait) const
+{
+    const bool wasInit = initialized();
+    runState_ = MASTER;
+
+    if (Pstream::master())
+    {
+        if (!wasInit)
+        {
+            // First time
+            Foam::mkDir(commsDir_);
+        }
+
+        const fileName lck(lockFile());
+
+        // Create lock file - only if it doesn't already exist
+        if (!Foam::isFile(lck))
+        {
+            Log << type() << ": creating lock file" << endl;
+
+            OFstream os(lck);
+            os << "status=openfoam\n";
+            os.flush();
+        }
+    }
+
+    if (wait)
+    {
+        waitForMaster();
+    }
+}
+
+
+void Foam::externalCoupler::useSlave(const bool wait) const
+{
+    const bool wasInit = initialized();
+    runState_ = SLAVE;
+
+    if (Pstream::master())
+    {
+        if (!wasInit)
+        {
+            // First time
+            Foam::mkDir(commsDir_);
+        }
+
+        Log << type() << ": removing lock file" << endl;
+
+        Foam::rm(lockFile());
+    }
+
+    if (wait)
+    {
+        waitForSlave();
+    }
+}
+
+
+bool Foam::externalCoupler::waitForMaster() const
+{
+    if (!initialized())
+    {
+        useMaster(); // was not initialized
+    }
+
+    bool isDone = false;
+    if (Pstream::master())
+    {
+        const fileName lck(lockFile());
+
+        double prevTime = -1;
+        double modTime = 0;
+
+        // Wait until file disappears (modTime == 0)
+        // But also check for status=done content in the file
+        while ((modTime = highResLastModified(lck)) > 0)
+        {
+            if (prevTime < modTime)
+            {
+                prevTime = modTime;
+                isDone = checkIsDone(lck);
+                if (isDone)
+                {
+                    break;
+                }
+            }
+            sleep(waitInterval_);
+        }
+    }
+
+    // MPI barrier
+    Pstream::scatter(isDone);
+
+    return !isDone;
+}
+
+
+bool Foam::externalCoupler::waitForSlave() const
+{
+    if (!initialized())
+    {
+        useSlave(); // was not initialized
+    }
+
+    bool isDone = false;
+    if (Pstream::master())
+    {
+        const fileName lck(lockFile());
+        unsigned totalTime = 0;
+
+        Log << type() << ": waiting for lock file to appear " << lck << endl;
+
+        while (!Foam::isFile(lck))
+        {
+            sleep(waitInterval_);
+
+            if (timeOut_ && (totalTime += waitInterval_) > timeOut_)
+            {
+                FatalErrorInFunction
+                    << "Wait time exceeded timeout of " << timeOut_
+                    << " s" << abort(FatalError);
+            }
+
+            Log << type() << ": wait time = " << totalTime << endl;
+        }
+
+        // But check for status=done content in the file
+        isDone = checkIsDone(lck);
+
+        Log << type() << ": found lock file " << lck << endl;
+    }
+
+    // MPI barrier
+    Pstream::scatter(isDone);
+
+    return !isDone;
+}
+
+
+Foam::fileName Foam::externalCoupler::resolveFile
+(
+    const word& file
+) const
+{
+    return fileName(baseDir()/file);
+}
+
+
+void Foam::externalCoupler::shutdown() const
+{
+    if (Pstream::master() && runState_ == MASTER && Foam::isDir(commsDir_))
+    {
+        const fileName lck(lockFile());
+
+        Log << type() << ": lock file status=done" << endl;
+        OFstream os(lck);
+        os  << "status=done\n";
+        os.flush();
+    }
+
+    runState_ = DONE;   // Avoid re-triggering in destructor
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/externalCoupler.H

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify i
+    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::externalCoupler
+
+Description
+    Encapsulates the logic for coordinating between OpenFOAM and an
+    external application.
+
+    This class provides a simple interface for explicit coupling with an
+    external application using plain text files located in the user-specified
+    communications directory.
+    These files are to be read/written on the master processor only.
+
+    The explicit coupling follows a master/slave model in which OpenFOAM
+    is the 'master' and the external application is the 'slave'.
+    The readiness to exchange information in either direction is handled
+    by a lock file (ie, OpenFOAM.lock).
+    If the lock file is present, the slave (external application) should wait
+    for the master (OpenFOAM) to complete.
+
+    When the master is finished its tasks and has prepared data for
+    the slave, the lock file is removed, instructing the external
+    source to take control of the program execution.
+
+    When the slave has completed its tasks, it will reinstate the lock file.
+
+    Example of the communication specification:
+    \verbatim
+    communication
+    {
+        commsDir        "${FOAM_CASE}/comms";
+        waitInterval    1;
+        timeOut         100;
+        initByExternal  no;
+    }
+    \endverbatim
+
+SourceFiles
+    externalCoupler.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef externalCoupler_H
+#define externalCoupler_H
+
+#include "fileName.H"
+#include "dictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class externalCoupler Declaration
+\*---------------------------------------------------------------------------*/
+
+class externalCoupler
+{
+    //- The run state (ie, who is currently in charge)
+    enum runState
+    {
+        NONE,       //!< Not initialized
+        MASTER,     //!< The master (OpenFOAM) is in charge
+        SLAVE,      //!< The slave (external program) is in charge
+        DONE        //!< Finished
+    };
+
+
+    // Private data
+
+        //- The current run (and initialization) state
+        mutable runState runState_;
+
+        //- Local path to communications directory
+        fileName commsDir_;
+
+        //- Interval time between checking for return data [s]
+        unsigned waitInterval_;
+
+        //- Timeout [s] while waiting for the external application
+        unsigned timeOut_;
+
+        //- Flag to indicate values are initialized by external application
+        bool slaveFirst_;
+
+        //- Local logging/verbosity flag
+        bool log;
+
+
+    // Private Member Functions
+
+        //- Return the file path to the base communications directory
+        const fileName& baseDir() const;
+
+        //- Return the file path to the lock file
+        fileName lockFile() const;
+
+        //- Disallow default bitwise copy construc
+        externalCoupler(const externalCoupler&) = delete;
+
+        //- Disallow default bitwise assignmen
+        void operator=(const externalCoupler&) = delete;
+
+public:
+
+    //- Runtime type information
+    TypeName("externalCoupler");
+
+
+    // Static data members
+
+        //- Name of the lock file
+        static word lockName;
+
+
+    // Constructors
+
+        //- Construct null using standard defaults
+        externalCoupler();
+
+        //- Construct from dictionary
+        externalCoupler(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~externalCoupler();
+
+
+    // Member Functions
+
+      // Access
+
+        //- True if state has been initialized
+        bool initialized() const;
+
+        //- External application provides initial values
+        bool slaveFirst() const;
+
+        //- Create lock file to indicate that OpenFOAM is in charge
+        //  Optionally wait for master as well.
+        void useMaster(const bool wait=false) const;
+
+        //- Wait for indication that OpenFOAM has supplied output.
+        //  This is when the lock file disappears, or it exists but with
+        //  "status=done" content.
+        //  \return False if lock file contains "status=done"
+        bool waitForMaster() const;
+
+        //- Remove lock file to indicate that the external program is in charge
+        //  Optionally wait for slave as well.
+        void useSlave(const bool wait=false) const;
+
+        //- Wait for indication that the external program has supplied input.
+        //  This is when the lock file appears.
+        //  \return False if lock file contains "status=done"
+        bool waitForSlave() const;
+
+        //- Return the file path in the communications directory
+        fileName resolveFile(const word& file) const;
+
+        //- Generate status=done in lock (only when run-state = master)
+        void shutdown() const;
+
+        //- Remove files written by OpenFOAM
+        void removeDirectory() const;
+
+
+      // Edit
+
+        //- Read communication settings from dictionary
+        bool readDict(const dictionary& dict);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointDisplacementPointPatchVectorField.C

@@ -0,0 +1,285 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "lumpedPointDisplacementPointPatchVectorField.H"
+#include "lumpedPointMovement.H"
+#include "lumpedPointIOMovement.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointFields.H"
+#include "surfaceFields.H"
+#include "volFields.H"
+#include "Time.H"
+#include "polyMesh.H"
+#include "displacementMotionSolver.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePointPatchTypeField
+    (
+        pointPatchVectorField,
+        lumpedPointDisplacementPointPatchVectorField
+    );
+}
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::labelList
+Foam::lumpedPointDisplacementPointPatchVectorField::patchIds
+(
+    const pointVectorField& pvf
+)
+{
+    const pointVectorField::Boundary& bf = pvf.boundaryField();
+
+    DynamicList<label> patchLst(bf.size());
+    forAll(bf, patchi)
+    {
+        // All patches of this type
+        if (isA<patchType>(bf[patchi]))
+        {
+            patchLst.append(patchi);
+            // or patchLst.append(bf[patchi].patch().index());
+        }
+    }
+
+    return patchLst.shrink();
+}
+
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+const Foam::pointField&
+Foam::lumpedPointDisplacementPointPatchVectorField::points0() const
+{
+    const objectRegistry& obr = this->patch().boundaryMesh().mesh().db();
+
+    // Obtain starting locations from the motionSolver
+    return obr.lookupObject<displacementMotionSolver>
+    (
+        "dynamicMeshDict"
+    ).points0();
+}
+
+
+const Foam::lumpedPointMovement&
+Foam::lumpedPointDisplacementPointPatchVectorField::movement() const
+{
+    const objectRegistry& obr = this->patch().boundaryMesh().mesh().db();
+    const lumpedPointIOMovement* ptr =
+        lumpedPointIOMovement::lookupInRegistry(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()
+    );
+
+    return objectRegistry::store(obj);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::lumpedPointDisplacementPointPatchVectorField::
+lumpedPointDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    fixedValuePointPatchField<vector>(p, iF)
+{}
+
+
+Foam::lumpedPointDisplacementPointPatchVectorField::
+lumpedPointDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValuePointPatchField<vector>(p, iF, dict)
+{}
+
+
+Foam::lumpedPointDisplacementPointPatchVectorField::
+lumpedPointDisplacementPointPatchVectorField
+(
+    const lumpedPointDisplacementPointPatchVectorField& pf,
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const pointPatchFieldMapper& mapper
+)
+:
+    fixedValuePointPatchField<vector>(pf, p, iF, mapper)
+{}
+
+
+Foam::lumpedPointDisplacementPointPatchVectorField::
+lumpedPointDisplacementPointPatchVectorField
+(
+    const lumpedPointDisplacementPointPatchVectorField& pf,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    fixedValuePointPatchField<vector>(pf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::lumpedPointDisplacementPointPatchVectorField::
+~lumpedPointDisplacementPointPatchVectorField()
+{
+    // de-register movement if in use and managed by this patch
+    const lumpedPointIOMovement* ptr = lumpedPointIOMovement::lookupInRegistry
+    (
+        this->patch().boundaryMesh().mesh().db()
+    );
+
+    if (ptr && ptr->ownerId() == this->patch().index())
+    {
+        movement().coupler().shutdown();
+
+        const_cast<lumpedPointIOMovement*>(ptr)->checkOut();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::lumpedPointDisplacementPointPatchVectorField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const bool masterPatch = (movement().ownerId() == this->patch().index());
+    if (masterPatch)
+    {
+        if (lumpedPointIOMovement::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
+        if (!movement().hasMapping())
+        {
+            const_cast<lumpedPointMovement&>(movement()).setMapping
+            (
+                mesh,
+                // All patches of this type
+                patchIds
+                (
+                    static_cast<const pointVectorField&>
+                    (
+                        this->internalField()
+                    )
+                ),
+                this->points0()
+            );
+        }
+
+
+        if
+        (
+            movement().coupler().initialized()
+         || !movement().coupler().slaveFirst()
+        )
+        {
+            // Synchronized for all processes
+            List<vector> forces, moments;
+            movement().forcesAndMoments(mesh, forces, moments);
+
+            if (lumpedPointIOMovement::debug)
+            {
+                Pout<<"gatherForces: " << forces << " called from patch "
+                    << this->patch().index() << endl;
+
+                if (Pstream::master())
+                {
+                    Pout<<"output forces to file: "
+                        << movement().locations() << " called from patch "
+                        << this->patch().index() << nl
+                        <<"# " << forces.size() << " force entries" << nl
+                        <<"# fx fy fz" << nl
+                        <<"output forces to file: "
+                        << forces << " called from patch "
+                        << this->patch().index() << endl;
+                }
+            }
+
+            if (Pstream::master())
+            {
+                movement().writeData(forces, moments);
+
+                // signal external source to execute
+                movement().coupler().useSlave();
+            }
+        }
+
+        // Wait for slave to provide data - includes MPI barrier
+        movement().coupler().waitForSlave();
+
+        // Read data passed back from external source - includes MPI barrier
+        const_cast<lumpedPointMovement&>(movement()).readState();
+    }
+
+    tmp<pointField> tdisp = movement().displacePoints
+    (
+        this->points0(),
+        this->patch().meshPoints()
+    );
+
+    this->operator==(tdisp);
+
+    fixedValuePointPatchField<vector>::updateCoeffs();
+}
+
+
+void Foam::lumpedPointDisplacementPointPatchVectorField::write(Ostream& os)
+const
+{
+    pointPatchField<vector>::write(os);
+    writeEntry("value", os);
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointDisplacementPointPatchVectorField.H

@@ -0,0 +1,200 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::lumpedPointDisplacementPointPatchVectorField
+
+Description
+
+    Interpolates pre-specified motion with motion specified as
+    pointVectorFields.
+
+    This is the point-patch responsible for managing the force
+    integration on a 'lumped-point' basis, waiting for the external
+    application, reading back the response from the external program
+    and updating the locations of the associated patch points
+    accordingly.
+
+    The internal patch type name is 'lumpedPointDisplacement'.
+
+    \heading Patch usage
+
+    Example:
+    \verbatim
+    walls
+    {
+        type                lumpedPointDisplacement;
+        value               uniform (0 0 0);
+        fieldName           wantedDisplacement;
+        interpolationScheme linear;
+    }
+    \endverbatim
+
+    This will scan the case for \a wantedDisplacement pointVectorFields
+    and interpolate those in time (using \c linear interpolation) to
+    obtain the current displacement.
+    The advantage of specifying displacement in this way is that it
+    automatically works through decomposePar.
+
+SourceFiles
+    lumpedPointDisplacementPointPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lumpedPointDisplacementPointPatchVectorField_H
+#define lumpedPointDisplacementPointPatchVectorField_H
+
+#include "fixedValuePointPatchField.H"
+#include "lumpedPointMovement.H"
+#include "lumpedPointState.H"
+#include "lumpedPointIOMovement.H"
+#include "labelList.H"
+#include "tmp.H"
+#include "pointField.H"
+#include "pointFieldsFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class interpolationWeights;
+
+/*---------------------------------------------------------------------------*\
+        Class lumpedPointDisplacementPointPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class lumpedPointDisplacementPointPatchVectorField
+:
+    public fixedValuePointPatchField<vector>
+{
+    // Private data
+
+    //- Convenience typedefs
+    typedef lumpedPointDisplacementPointPatchVectorField patchType;
+    typedef DimensionedField<vector, pointMesh> fieldType;
+
+    // Private Member Functions
+
+protected:
+
+    //- The starting locations (obtained from the motionSolver).
+    const pointField& points0() const;
+
+    //- The auto-vivifying singleton for movement.
+    const lumpedPointMovement& movement() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("lumpedPointDisplacement");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        lumpedPointDisplacementPointPatchVectorField
+        (
+            const pointPatch& p,
+            const DimensionedField<vector, pointMesh>& iF
+        );
+
+        //- Construct from patch, internal field and dictionary
+        lumpedPointDisplacementPointPatchVectorField
+        (
+            const pointPatch& p,
+            const DimensionedField<vector, pointMesh>& iF,
+            const dictionary& dict
+        );
+
+        //- Construct by mapping given patchField<vector> onto a new patch
+        lumpedPointDisplacementPointPatchVectorField
+        (
+            const lumpedPointDisplacementPointPatchVectorField& pf,
+            const pointPatch& p,
+            const DimensionedField<vector, pointMesh>& iF,
+            const pointPatchFieldMapper& mapper
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<pointPatchField<vector>> clone() const
+        {
+            return autoPtr<pointPatchField<vector>>
+            (
+                new lumpedPointDisplacementPointPatchVectorField
+                (
+                    *this
+                )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        lumpedPointDisplacementPointPatchVectorField
+        (
+            const lumpedPointDisplacementPointPatchVectorField& pf,
+            const DimensionedField<vector, pointMesh>& iF
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual autoPtr<pointPatchField<vector>> clone
+        (
+            const DimensionedField<vector, pointMesh>& iF
+        ) const
+        {
+            return autoPtr<pointPatchField<vector>>
+            (
+                new lumpedPointDisplacementPointPatchVectorField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+    //- Destructor
+    virtual ~lumpedPointDisplacementPointPatchVectorField();
+
+
+    // Member functions
+
+        //- The ids for all patches of this type
+        static labelList patchIds(const pointVectorField& pvf);
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointIOMovement.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "lumpedPointIOMovement.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(lumpedPointIOMovement, 0);
+}
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+const Foam::lumpedPointIOMovement*
+Foam::lumpedPointIOMovement::lookupInRegistry(const objectRegistry& obr)
+{
+    return obr.lookupObjectPtr<lumpedPointIOMovement>
+    (
+        lumpedPointMovement::dictionaryName
+    );
+}
+
+
+Foam::autoPtr<Foam::lumpedPointIOMovement>
+Foam::lumpedPointIOMovement::New
+(
+    const objectRegistry& obr,
+    label ownerId
+)
+{
+    return autoPtr<lumpedPointIOMovement>
+    (
+        new lumpedPointIOMovement
+        (
+            IOobject
+            (
+                lumpedPointMovement::dictionaryName,
+                obr.time().caseSystem(),
+                obr,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                true // register object
+            ),
+            ownerId  // tag this patch as owner too
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::lumpedPointIOMovement::lumpedPointIOMovement
+(
+    const IOobject& io,
+    label ownerId
+)
+:
+    lumpedPointMovement(),
+    regIOobject(io)
+{
+    bool ok =
+    (
+        readOpt() == IOobject::MUST_READ
+     || readOpt() == IOobject::MUST_READ_IF_MODIFIED
+    );
+
+    if (ok)
+    {
+        ok = readData(readStream(typeName));
+        close();
+
+        if (ok)
+        {
+            this->ownerId(ownerId);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::lumpedPointIOMovement::~lumpedPointIOMovement()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+bool Foam::lumpedPointIOMovement::readData(Istream& is)
+{
+    dictionary dict(is);
+
+    readDict(dict);
+
+    return is.check(FUNCTION_NAME);
+}
+
+
+bool Foam::lumpedPointIOMovement::writeData(Ostream& os) const
+{
+    os  << *this;
+    return os.good();
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const lumpedPointIOMovement& obj)
+{
+    obj.writeDict(os);
+
+    os.check(FUNCTION_NAME);
+    return os;
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointIOMovement.H

@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::lumpedPointIOMovement
+
+Description
+    IO-registered version of lumpedPointMovement.
+
+SourceFiles
+    lumpedPointMovement.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lumpedPointIOMovement_H
+#define lumpedPointIOMovement_H
+
+#include "lumpedPointMovement.H"
+#include "regIOobject.H"
+#include "className.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+class lumpedPointIOMovement;
+
+Ostream& operator<<(Ostream& os, const lumpedPointIOMovement& obj);
+
+/*---------------------------------------------------------------------------*\
+                    Class lumpedPointIOMovement Declaration
+\*---------------------------------------------------------------------------*/
+
+class lumpedPointIOMovement
+:
+    public lumpedPointMovement,
+    public regIOobject
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        lumpedPointIOMovement(const lumpedPointIOMovement&) = delete;
+
+        //- Disallow default bitwise assignment
+        void operator=(const lumpedPointIOMovement&) = delete;
+
+public:
+
+    //- Runtime type information
+    TypeName("lumpedPointMovement");
+
+
+    // Static Member Functions
+
+        //- Lookup pointer to object in the object-registry,
+        //  return nullptr if found.
+        static const lumpedPointIOMovement* lookupInRegistry
+        (
+            const objectRegistry& obr
+        );
+
+        //- Create a new object in the registry by reading system dictionary
+        static autoPtr<lumpedPointIOMovement> New
+        (
+            const objectRegistry& obr,
+            label ownerId = -1
+        );
+
+
+    // Constructors
+
+        //- Construct from IOobject, optionally with some owner information
+        explicit lumpedPointIOMovement
+        (
+            const IOobject& io,
+            label ownerId = -1
+        );
+
+
+    //- Destructor
+    ~lumpedPointIOMovement();
+
+
+    // Member Functions
+
+        //- readData member function used by regIOobject
+        bool readData(Istream& is);
+
+        //- writeData member function required by regIOobject
+        bool writeData(Ostream& os) const;
+
+
+    // IOstream Operators
+
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const lumpedPointIOMovement& obj
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointMovement.C

@@ -0,0 +1,790 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "lumpedPointMovement.H"
+#include "lumpedPointIOMovement.H"
+#include "demandDrivenData.H"
+#include "linearInterpolationWeights.H"
+#include "IFstream.H"
+#include "OFstream.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "PtrMap.H"
+#include "faceZoneMesh.H"
+#include "PstreamReduceOps.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::Enum<Foam::lumpedPointMovement::outputFormatType>
+    Foam::lumpedPointMovement::formatNames
+    {
+        { outputFormatType::PLAIN, "plain" },
+        { outputFormatType::DICTIONARY, "dictionary" }
+    };
+
+
+const Foam::word
+Foam::lumpedPointMovement::dictionaryName("lumpedPointMovement");
+
+
+namespace Foam
+{
+//! \cond fileScope
+//- Write list content with size, bracket, content, bracket one-per-line.
+//  This makes for consistent for parsing, regardless of the list length.
+template <class T>
+static void writeList(Ostream& os, const string& header, const UList<T>& lst)
+{
+    // Header string
+    os  << header.c_str() << nl;
+
+    // Write size and start delimiter
+    os  << lst.size() << nl
+        << token::BEGIN_LIST << nl;
+
+    // Write contents
+    forAll(lst, i)
+    {
+        os << lst[i] << nl;
+    }
+
+    // Write end delimiter
+    os << token::END_LIST << token::END_STATEMENT << nl << nl;
+}
+//! \endcond
+
+}
+// namespace Foam
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::IOobject Foam::lumpedPointMovement::selectIO
+(
+    const IOobject& io,
+    const fileName& f
+)
+{
+    return
+    (
+        f.size()
+      ? IOobject        // construct from filePath instead
+        (
+            f,
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            io.registerObject(),
+            io.globalObject()
+        )
+      : io
+    );
+}
+
+
+// * * * * * * * * * * * * * 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),
+    relax_(1),
+    interpolationScheme_(linearInterpolationWeights::typeName),
+    ownerId_(-1),
+    boundBox_(),
+    centre_(),
+    autoCentre_(true),
+    forcesDict_(),
+    coupler_(),
+    inputName_("positions.in"),
+    outputName_("forces.out"),
+    inputFormat_(lumpedPointState::inputFormatType::DICTIONARY),
+    outputFormat_(outputFormatType::DICTIONARY),
+    state0_(),
+    state_(),
+    thresholdPtr_(0),
+    interpolatorPtr_()
+{}
+
+
+Foam::lumpedPointMovement::lumpedPointMovement
+(
+    const dictionary& dict,
+    label ownerId
+)
+:
+    axis_(0,0,1),
+    locations_(),
+    division_(0),
+    relax_(1),
+    interpolationScheme_
+    (
+        dict.lookupOrDefault<word>
+        (
+            "interpolationScheme",
+            linearInterpolationWeights::typeName
+        )
+    ),
+    ownerId_(ownerId),
+    boundBox_(),
+    centre_(),
+    autoCentre_(true),
+    forcesDict_(),
+    coupler_(),
+    state0_(),
+    state_(),
+    thresholdPtr_(0),
+    interpolatorPtr_()
+{
+    readDict(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+
+Foam::lumpedPointMovement::~lumpedPointMovement()
+{
+    deleteDemandDrivenData(thresholdPtr_);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::lumpedPointMovement::readDict(const dictionary& dict)
+{
+    // assume the worst
+    deleteDemandDrivenData(thresholdPtr_);
+
+    dict.lookup("axis") >> axis_;
+
+    division_ = 0;
+    if (dict.readIfPresent("division", division_))
+    {
+        if (division_ < 0)
+        {
+            division_ = 0;
+        }
+        else if (division_ > 1)
+        {
+            division_ = 1;
+        }
+    }
+
+    dict.readIfPresent("relax", relax_);
+
+    dict.lookup("locations") >> locations_;
+
+    if (dict.readIfPresent("interpolationScheme", interpolationScheme_))
+    {
+        interpolatorPtr_.clear();
+    }
+
+    autoCentre_ = !dict.readIfPresent("centre", centre_);
+
+    // Initial state from initial locations, with zero rotation
+    tmp<vectorField> pts = locations_ * axis_;
+    state0_ = lumpedPointState(pts);
+    state_  = state0_;
+
+    forcesDict_.merge(dict.subOrEmptyDict("forces"));
+
+    const dictionary& commDict = dict.subDict("communication");
+    coupler_.readDict(commDict);
+
+    // TODO: calcFrequency_  = dict.lookupOrDefault("calcFrequency", 1);
+
+    commDict.lookup("inputName")  >> inputName_;
+    commDict.lookup("outputName") >> outputName_;
+
+    inputFormat_ = lumpedPointState::formatNames.lookup
+    (
+        "inputFormat",
+        commDict
+    );
+
+    outputFormat_ = lumpedPointMovement::formatNames.lookup
+    (
+        "outputFormat",
+        commDict
+    );
+}
+
+
+void Foam::lumpedPointMovement::setBoundBox
+(
+    const polyMesh& mesh,
+    const labelUList& patchIds,
+    const pointField& points0
+)
+{
+    boundBox_ = boundBox::invertedBox;
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    for (const label patchId : patchIds)
+    {
+        boundBox_.add(points0, patches[patchId].meshPoints());
+    }
+
+    boundBox_.reduce();
+
+    if (autoCentre_)
+    {
+        centre_ = boundBox_.midpoint();
+        centre_ -= (centre_ & axis_) * axis_;
+        if (lumpedPointIOMovement::debug)
+        {
+            Pout<<"autoCentre on " << centre_ << " boundBox: " << boundBox_ << endl;
+        }
+    }
+    else
+    {
+        // User-specified centre
+        if (lumpedPointIOMovement::debug)
+        {
+            Pout<<"centre on " << centre_ << " boundBox: " << boundBox_ << endl;
+        }
+    }
+}
+
+
+void Foam::lumpedPointMovement::setMapping
+(
+    const polyMesh& mesh,
+    const labelUList& patchIds,
+    const pointField& points0
+)
+{
+    setBoundBox(mesh, patchIds, points0);
+
+    faceZones_.clear();
+
+    // Local storage location (boundary faces only)
+    const label firstFace = mesh.nInternalFaces();
+    labelList faceToZoneID(mesh.nFaces() - firstFace, -1);
+
+    // Number of faces per zone
+    labelList nFaces(thresholds().size(), 0);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    for (const label patchId : patchIds)
+    {
+        const polyPatch& pp = patches[patchId];
+
+        // Classify faces into respective zones according to their centres
+        label meshFaceI = pp.start();
+
+        forAll(pp, i)
+        {
+            const face& f = mesh.faces()[meshFaceI];
+            label zoneId = this->threshold(f.centre(points0));
+
+            // localize storage location
+            faceToZoneID[meshFaceI - firstFace] = zoneId;
+            ++nFaces[zoneId];
+
+            ++meshFaceI;
+        }
+    }
+
+    if (lumpedPointIOMovement::debug)
+    {
+        Pout<<"faces per zone:" << nFaces << endl;
+    }
+
+    faceZones_.setSize(nFaces.size());
+
+    label nZonedFaces = 0;
+    forAll(nFaces, zonei)
+    {
+        label n = nFaces[zonei];
+        labelList& addr = faceZones_[zonei];
+        addr.setSize(n);
+
+        n = 0;
+        forAll(faceToZoneID, faceI)
+        {
+            if (faceToZoneID[faceI] == zonei)
+            {
+                // from local storage to global mesh face
+                label meshFaceI = faceI + firstFace;
+
+                addr[n] = meshFaceI;
+                ++n;
+            }
+        }
+
+        addr.setSize(n);
+        nZonedFaces += n;
+
+        if (lumpedPointIOMovement::debug)
+        {
+            Pout<< "Created faceZone[" << zonei << "] "
+                << returnReduce(n, sumOp<label>()) << " faces, "
+                << thresholds()[zonei] << " threshold" << nl;
+        }
+    }
+}
+
+
+bool Foam::lumpedPointMovement::forcesAndMoments
+(
+    const polyMesh& pmesh,
+    List<vector>& forces,
+    List<vector>& moments
+) const
+{
+    forces.setSize(faceZones_.size());
+    moments.setSize(faceZones_.size());
+
+    if (faceZones_.empty())
+    {
+        WarningInFunction
+            << "Attempted to calculate forces without setMapping()"
+            << nl << endl;
+
+        forces.setSize(thresholds().size(), Zero);
+        moments.setSize(thresholds().size(), Zero);
+        return false;
+    }
+
+    // Initialize with zero
+    forces = Zero;
+    moments = Zero;
+
+    // The mass centres
+    const pointField& lumpedCentres = state().points();
+
+    const polyBoundaryMesh& patches = pmesh.boundaryMesh();
+
+    const word pName = forcesDict_.lookupOrDefault<word>("p", "p");
+
+    scalar pRef   = forcesDict_.lookupOrDefault<scalar>("pRef",   0.0);
+    scalar rhoRef = forcesDict_.lookupOrDefault<scalar>("rhoRef", 1.0);
+
+
+    // Calculated force per patch - cache
+    PtrMap<vectorField> forceOnPatches;
+
+    const volScalarField* pPtr = pmesh.lookupObjectPtr<volScalarField>(pName);
+
+    // fvMesh and has pressure field
+    if (isA<fvMesh>(pmesh) && pPtr)
+    {
+        const fvMesh& mesh = dynamicCast<const fvMesh>(pmesh);
+        const volScalarField& p = *pPtr;
+
+        // face centres (on patches)
+        const surfaceVectorField::Boundary& patchCf =
+            mesh.Cf().boundaryField();
+
+        // face areas (on patches)
+        const surfaceVectorField::Boundary& patchSf =
+            mesh.Sf().boundaryField();
+
+        // Pressure (on patches)
+        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);
+
+        // Scale pRef by density for incompressible simulations
+        pRef /= rhoRef;
+
+        forAll(faceZones_, zonei)
+        {
+            const labelList& zn = faceZones_[zonei];
+
+            forAll(zn, localFaceI)
+            {
+                const label meshFaceI = zn[localFaceI];
+
+                // locate which patch this face is on,
+                // and its offset within the patch
+                const label patchI = patches.whichPatch(meshFaceI);
+                if (patchI == -1)
+                {
+                    // Should not happen - could also warn
+                    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];
+
+                // effective torque arm:
+                // - translated into the lumped-points coordinate system
+                //   prior to determining the distance
+                const vector lever =
+                (
+                    patchCf[patchI][patchFaceI] - centre_
+                  - lumpedCentres[zonei]
+                );
+
+                // Moment from the patch-face into sum
+                moments[zonei] += lever ^ forceOnPatch[patchFaceI];
+            }
+        }
+    }
+    else
+    {
+        Info<<"no pressure" << endl;
+    }
+
+    Pstream::listCombineGather(forces, plusEqOp<vector>());
+    Pstream::listCombineScatter(forces);
+
+    Pstream::listCombineGather(moments, plusEqOp<vector>());
+    Pstream::listCombineScatter(moments);
+
+    return true;
+}
+
+
+// \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
+(
+    const pointField& points0,
+    const labelList& pointLabels
+) const
+{
+    return displacePoints(state(), points0, pointLabels);
+}
+
+
+Foam::tmp<Foam::pointField>
+Foam::lumpedPointMovement::displacePoints
+(
+    const lumpedPointState& state,
+    const pointField& points0,
+    const labelList& pointLabels
+) const
+{
+    // storage for the interpolation indices/weights
+    labelList   indices;
+    scalarField weights;
+    const interpolationWeights& interp = interpolator();
+
+    // The mass centres
+    const pointField& lumpedCentres = state.points();
+
+    // local-to-global transformation tensor
+    const tensorField& localToGlobal = state.rotations();
+
+    // could also verify the sizes (state vs original)
+
+    tmp<pointField> tdisp(new pointField(pointLabels.size()));
+    pointField& disp = tdisp.ref();
+
+    forAll(pointLabels, ptI)
+    {
+        const point& p0 = points0[pointLabels[ptI]];
+
+        // Interpolation factor based on the axis component
+        // of the original point (location)
+        scalar where = p0 & axis_;
+
+        interp.valueWeights(where, indices, weights);
+
+        vector origin    = interpolatedValue(lumpedCentres, indices, weights);
+        tensor rotTensor = interpolatedValue(localToGlobal, indices, weights);
+
+        if (indices.size() == 1)
+        {
+            // out-of-bounds, use corresponding end-point
+            where = locations_[indices[0]];
+        }
+
+        // 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;
+    }
+
+    return tdisp;
+}
+
+
+void Foam::lumpedPointMovement::writeDict(Ostream& os) const
+{
+    os.writeEntry("axis", axis_);
+    os.writeEntry("locations", locations_);
+    os.writeEntry("division", division_);
+    // os.writeEntry("interpolationScheme", interpolationScheme_);
+}
+
+
+bool Foam::lumpedPointMovement::readState()
+{
+    lumpedPointState prev = state_;
+
+    bool status = state_.readData
+    (
+        inputFormat_,
+        coupler().resolveFile(inputName_)
+    );
+
+    state_.relax(relax_, prev);
+
+    return status;
+}
+
+
+bool Foam::lumpedPointMovement::writeData
+(
+    const UList<vector>& forces
+) const
+{
+    if (!Pstream::master())
+    {
+        return false;
+    }
+
+    const fileName output(coupler().resolveFile(outputName_));
+    OFstream os(output); // ASCII
+
+    if (outputFormat_ == outputFormatType::PLAIN)
+    {
+        os  <<"# output from OpenFOAM" << nl
+            <<"# N, points, forces" << nl
+            << this->size() << nl;
+
+        const char* zeroVector = "0 0 0";
+
+        forAll(locations_, i)
+        {
+            const vector pos = locations_[i] * axis_;
+
+            os  << pos.x() << ' '
+                << pos.y() << ' '
+                << pos.z() << ' ';
+
+            if (i < forces.size())
+            {
+                os  << forces[i].x() << ' '
+                    << forces[i].y() << ' '
+                    << forces[i].z();
+            }
+            else
+            {
+                os << zeroVector;
+            }
+
+            os  << nl;
+        }
+    }
+    else
+    {
+        // Make it easier for external programs to parse
+        // - exclude the usual OpenFOAM 'FoamFile' header
+        // - ensure lists have consistent format
+
+        os  <<"// output from OpenFOAM" << nl << nl;
+
+        writeList(os, "points", (locations_*axis_)());
+        writeList(os, "forces", forces);
+    }
+
+    return true;
+}
+
+
+bool Foam::lumpedPointMovement::writeData
+(
+    const UList<vector>& forces,
+    const UList<vector>& moments
+) const
+{
+    if (!Pstream::master())
+    {
+        return false;
+    }
+
+    const fileName output(coupler().resolveFile(outputName_));
+    OFstream os(output); // ASCII
+
+    if (outputFormat_ == outputFormatType::PLAIN)
+    {
+        os  <<"# output from OpenFOAM" << nl
+            <<"# N, points, forces, moments" << nl
+            << this->size() << nl;
+
+        const char* zeroVector = "0 0 0";
+
+        forAll(locations_, i)
+        {
+            const vector pos = locations_[i] * axis_;
+
+            os  << pos.x() << ' '
+                << pos.y() << ' '
+                << pos.z() << ' ';
+
+            if (i < forces.size())
+            {
+                os  << forces[i].x() << ' '
+                    << forces[i].y() << ' '
+                    << forces[i].z() << ' ';
+            }
+            else
+            {
+                os << zeroVector << ' ';
+            }
+
+            if (i < moments.size())
+            {
+                os  << moments[i].x() << ' '
+                    << moments[i].y() << ' '
+                    << moments[i].z();
+            }
+            else
+            {
+                os  << zeroVector;
+            }
+            os  << nl;
+        }
+    }
+    else
+    {
+        // Make it easier for external programs to parse
+        // - exclude the usual OpenFOAM 'FoamFile' header
+        // - ensure lists have consistent format
+
+        os  <<"// output from OpenFOAM" << nl << nl;
+
+        writeList(os, "points", (locations_*axis_)());
+        writeList(os, "forces", forces);
+        writeList(os, "moments", moments);
+    }
+
+    return true;
+}
+
+
+const Foam::interpolationWeights&
+Foam::lumpedPointMovement::interpolator() const
+{
+    if (!interpolatorPtr_.valid())
+    {
+        interpolatorPtr_ = interpolationWeights::New
+        (
+            interpolationScheme_,
+            locations_
+        );
+    }
+
+    return interpolatorPtr_();
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointMovement.H

@@ -0,0 +1,399 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::lumpedPointMovement
+
+Description
+    The movement 'driver' that describes initial point locations, the
+    segmentation for pressure integration, the current state of the
+    points/rotations, and forwarding to the communication
+    coordinator. The 'lumpedPointIOMovement' class is simply a
+    registered version of the same.
+
+SourceFiles
+    lumpedPointMovement.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lumpedPointMovement_H
+#define lumpedPointMovement_H
+
+#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 "IOobject.H"
+#include "tmp.H"
+#include "faceZoneMeshFwd.H"
+#include "externalCoupler.H"
+#include "lumpedPointState.H"
+#include "boundBox.H"
+#include "Enum.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+// Forward declarations
+class polyMesh;
+
+/*---------------------------------------------------------------------------*\
+                     Class lumpedPointMovement Declaration
+\*---------------------------------------------------------------------------*/
+
+class lumpedPointMovement
+{
+public:
+
+    //- Output format types
+    enum class outputFormatType
+    {
+        PLAIN,
+        DICTIONARY
+    };
+
+    // Static data
+
+    //- Names for the output format types
+    static const Enum<outputFormatType> formatNames;
+
+
+private:
+
+    // Private data
+
+        //- Reference axis for the locations
+        vector axis_;
+
+        //- 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_;
+
+        //- The relaxation factor when moving points (default: 1)
+        scalar relax_;
+
+        //- The interpolation type (linear|spline)
+        word interpolationScheme_;
+
+        //- Optional owner information
+        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
+        externalCoupler coupler_;
+
+        //- File io
+        word inputName_;
+        word outputName_;
+        lumpedPointState::inputFormatType inputFormat_;
+        outputFormatType outputFormat_;
+
+
+    // 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_;
+
+
+    // 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;
+
+
+        //- Disallow default bitwise copy constructor
+        lumpedPointMovement(const lumpedPointMovement&) = delete;
+
+        //- Disallow default bitwise assignment
+        void operator=(const lumpedPointMovement&) = delete;
+
+public:
+
+    // Static data members
+
+        //- The standard dictionary name
+        static const word dictionaryName;
+
+        //- Helper: return IOobject with optionally absolute path provided
+        static IOobject selectIO(const IOobject&, const fileName&);
+
+
+    // Constructors
+
+        //- Construct null
+        lumpedPointMovement();
+
+        //- Construct from dictionary, optionally with some owner information
+        lumpedPointMovement(const dictionary& dict, label ownerId=-1);
+
+
+    //- Destructor
+    virtual ~lumpedPointMovement();
+
+
+    // Member Functions
+
+        //- Update settings from dictionary
+        void readDict(const dictionary& dict);
+
+
+        //- If no number of lumped points (locations) were specified
+        inline bool empty() const;
+
+        //- 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 externalCoupler& coupler() const;
+
+        //- Communication control
+        inline externalCoupler& coupler();
+
+        //- The initial state (positions/rotations)
+        inline const lumpedPointState& state0() const;
+
+        //- The current state (positions/rotations)
+        inline const lumpedPointState& state() const;
+
+        //- The relaxation factor when changing states
+        inline scalar relax() const;
+
+        //- The relaxation factor when changing states
+        inline scalar& relax();
+
+        //- The input (state) file name
+        inline const word& inputName() const;
+
+        //- The output (forces) file name
+        inline const word& outputName() const;
+
+        //- The input (state) file format
+        inline lumpedPointState::inputFormatType inputFormat() const;
+
+        //- The output (forces) file format
+        inline lumpedPointMovement::outputFormatType outputFormat() const;
+
+
+        //- Define the bounding-box required to enclose the specified patches.
+        //  Calculates the centre as required.
+        //
+        //  \param mesh     The volume mesh reference
+        //  \param patchIds The patch ids to be included in the bounding box.
+        //  \param points0  The initial mesh points, prior to movement
+        void setBoundBox
+        (
+            const polyMesh& mesh,
+            const labelUList& patchIds,
+            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.
+        //
+        //  \param mesh     The volume mesh reference
+        //  \param patchIds The patch ids to be included in the mapping
+        //  \param points0  The initial mesh points, prior to movement
+        void setMapping
+        (
+            const polyMesh& mesh,
+            const labelUList& patchIds,
+            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;
+
+
+        //- The forces and moments acting on each pressure-zone.
+        //  The zones must be previously defined via setMapping.
+        bool forcesAndMoments
+        (
+            const polyMesh& pmesh,
+            List<vector>& forces,
+            List<vector>& moments
+        ) const;
+
+
+        //- Displace points according to the current state
+        tmp<pointField> displacePoints
+        (
+            const pointField& points0,
+            const labelList& pointLabels
+        ) const;
+
+
+        //- Displace points according to specified state
+        tmp<pointField> displacePoints
+        (
+            const lumpedPointState& state,
+            const pointField& points0,
+            const labelList& pointLabels
+        ) const;
+
+
+        //- Interpolation weights
+        const interpolationWeights& interpolator() const;
+
+        //- Write axis, locations, division as a dictionary
+        void writeDict(Ostream& os) const;
+
+
+        //- Write points, forces
+        bool writeData
+        (
+            const UList<vector>& forces
+        ) const;
+
+        //- Write points, forces, moments
+        bool writeData
+        (
+            const UList<vector>& forces,
+            const UList<vector>& moments
+        ) const;
+
+
+        //- Read state from file, applying relaxation as requested
+        bool readState();
+
+        //- Write state as VTK PolyData format.
+        void writeStateVTP(const fileName& file) const;
+
+        //- Write forces on points as VTK PolyData format.
+        void writeForcesAndMomentsVTP
+        (
+            const fileName& file,
+            const UList<vector>& forces,
+            const UList<vector>& moments
+        ) const;
+
+
+        //- Write pressure-zones geometry, write as VTK PolyData format.
+        void writeZonesVTP
+        (
+            const fileName& file,
+            const polyMesh& mesh,
+            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;
+
+        //- Write displaced geometry according to the specified state,
+        //  write as VTK PolyData format.
+        void writeVTP
+        (
+            const fileName& file,
+            const lumpedPointState& state,
+            const polyMesh& mesh,
+            const labelUList& patchLst,
+            const pointField& points0
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "lumpedPointMovementI.H"
+
+#endif
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointMovementI.H

@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+\*---------------------------------------------------------------------------*/
+
+
+inline bool Foam::lumpedPointMovement::empty() const
+{
+    return locations_.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_;
+}
+
+
+inline Foam::label Foam::lumpedPointMovement::ownerId() const
+{
+    return ownerId_;
+}
+
+
+inline void Foam::lumpedPointMovement::ownerId(label id)
+{
+    ownerId_ = id;
+}
+
+
+inline const Foam::scalarField& Foam::lumpedPointMovement::thresholds() const
+{
+    if (!thresholdPtr_)
+    {
+        calcThresholds();
+    }
+
+    return *thresholdPtr_;
+}
+
+
+inline Foam::label
+Foam::lumpedPointMovement::threshold(const point& position) const
+{
+    return threshold(position & axis_);
+}
+
+
+inline const Foam::externalCoupler& Foam::lumpedPointMovement::coupler() const
+{
+    return coupler_;
+}
+
+
+inline Foam::externalCoupler& Foam::lumpedPointMovement::coupler()
+{
+    return coupler_;
+}
+
+
+//- The initial state (positions/rotations)
+inline const Foam::lumpedPointState& Foam::lumpedPointMovement::state0() const
+{
+    return state0_;
+}
+
+
+inline const Foam::lumpedPointState& Foam::lumpedPointMovement::state() const
+{
+    return state_;
+}
+
+
+inline Foam::scalar Foam::lumpedPointMovement::relax() const
+{
+    return relax_;
+}
+
+
+inline Foam::scalar& Foam::lumpedPointMovement::relax()
+{
+    return relax_;
+}
+
+
+inline const Foam::word& Foam::lumpedPointMovement::inputName() const
+{
+    return inputName_;
+}
+
+
+inline const Foam::word& Foam::lumpedPointMovement::outputName() const
+{
+    return outputName_;
+}
+
+
+inline Foam::lumpedPointState::inputFormatType
+Foam::lumpedPointMovement::inputFormat() const
+{
+    return inputFormat_;
+}
+
+
+Foam::lumpedPointMovement::outputFormatType
+Foam::lumpedPointMovement::outputFormat() const
+{
+    return outputFormat_;
+}
+
+
+inline bool Foam::lumpedPointMovement::hasMapping() const
+{
+    return !faceZones_.empty();
+}
+
+
+inline const Foam::List<Foam::labelList>&
+Foam::lumpedPointMovement::zones() const
+{
+    return faceZones_;
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointMovementWriter.C

@@ -0,0 +1,447 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "lumpedPointMovement.H"
+#include "polyMesh.H"
+#include "OFstream.H"
+#include "uindirectPrimitivePatch.H"
+
+#include "foamVtkOutput.H"
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::lumpedPointMovement::writeStateVTP(const fileName& file) const
+{
+    state().writeVTP(file, axis());
+}
+
+
+void Foam::lumpedPointMovement::writeForcesAndMomentsVTP
+(
+    const fileName& file,
+    const UList<vector>& forces,
+    const UList<vector>& moments
+) const
+{
+    OFstream fos(file);
+    std::ostream& os = fos.stdStream();
+
+    autoPtr<vtk::formatter> format = vtk::newFormatter
+    (
+        os,
+        vtk::formatType::INLINE_ASCII
+    );
+
+    format().xmlHeader()
+        .beginVTKFile(vtk::fileTag::POLY_DATA, "0.1");
+
+    //
+    // The 'spine' of lumped mass points
+    //
+    const label nPoints = state().points().size();
+
+    {
+        format()
+            .openTag(vtk::fileTag::PIECE)
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_POINTS, nPoints)
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_VERTS,  nPoints)
+            .closeTag();
+
+        // 'points'
+        {
+            const uint64_t payLoad = (nPoints*3* sizeof(float));
+
+            format()
+                .tag(vtk::fileTag::POINTS)
+                .openDataArray<float, 3>(vtk::dataArrayAttr::POINTS)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            vtk::writeList(format(), state().points());
+            format().flush();
+
+            format()
+                .endDataArray()
+                .endTag(vtk::fileTag::POINTS);
+        }
+
+        // <Verts>
+        format().tag(vtk::fileTag::VERTS);
+
+        //
+        // 'connectivity'
+        //
+        {
+            const uint64_t payLoad = (nPoints*sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            for (label i=0; i<nPoints; ++i)
+            {
+                format().write(i);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        //
+        // 'offsets'  (connectivity offsets)
+        // = linear mapping onto points (with 1 offset)
+        //
+        {
+            const uint64_t payLoad = (nPoints*sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::OFFSETS)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            for (label i=0; i<nPoints; ++i)
+            {
+                format().write(i+1);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        format().endTag(vtk::fileTag::VERTS);
+        // </Verts>
+    }
+
+    format().tag(vtk::fileTag::POINT_DATA);
+
+    // forces
+    if (forces.size() == nPoints)
+    {
+        const uint64_t payLoad = (nPoints * 3 * sizeof(float));
+
+        format().openDataArray<float, 3>("forces")
+            .closeTag();
+
+        format().writeSize(payLoad);
+        vtk::writeList(format(), forces);
+        format().flush();
+
+        format().endDataArray();
+    }
+
+    // moments
+    if (moments.size() == nPoints)
+    {
+        const uint64_t payLoad = (nPoints * 3 * sizeof(float));
+
+        format().openDataArray<float, 3>("moments")
+            .closeTag();
+
+        format().writeSize(payLoad);
+        vtk::writeList(format(), moments);
+        format().flush();
+
+        format().endDataArray();
+    }
+
+    format().endTag(vtk::fileTag::POINT_DATA);
+
+    format().endTag(vtk::fileTag::PIECE);
+    format().endTag(vtk::fileTag::POLY_DATA);
+    format().endVTKFile();
+}
+
+
+void Foam::lumpedPointMovement::writeZonesVTP
+(
+    const fileName& file,
+    const polyMesh& mesh,
+    const pointField& points0
+) const
+{
+    OFstream fos(file);
+    std::ostream& os = fos.stdStream();
+
+    autoPtr<vtk::formatter> format = vtk::newFormatter
+    (
+        os,
+        vtk::formatType::INLINE_ASCII
+    );
+
+    format().xmlHeader()
+        .beginVTKFile(vtk::fileTag::POLY_DATA, "0.1");
+
+    forAll(faceZones_, zoneI)
+    {
+        uindirectPrimitivePatch pp
+        (
+            UIndirectList<face>(mesh.faces(), faceZones_[zoneI]),
+            points0
+        );
+
+        format()
+            .openTag(vtk::fileTag::PIECE)
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_POINTS, pp.nPoints())
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_POLYS,  pp.size())
+            .closeTag();
+
+        // 'points'
+        {
+            const uint64_t payLoad = (pp.nPoints()*3* sizeof(float));
+
+            format()
+                .tag(vtk::fileTag::POINTS)
+                .openDataArray<float, 3>(vtk::dataArrayAttr::POINTS)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            vtk::writeList(format(), pp.localPoints());
+            format().flush();
+
+            format()
+                .endDataArray()
+                .endTag(vtk::fileTag::POINTS);
+        }
+
+        // <Polys>
+        format().tag(vtk::fileTag::POLYS);
+
+        //
+        // 'connectivity'
+        //
+        {
+            uint64_t payLoad = 0;
+            forAll(pp, facei)
+            {
+                const face& f = pp[facei];
+                payLoad += f.size();
+            }
+
+            format().openDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY)
+                .closeTag();
+
+            format().writeSize(payLoad * sizeof(label));
+
+            for (const face& f : pp.localFaces())
+            {
+                vtk::writeList(format(), f);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        //
+        // 'offsets'  (connectivity offsets)
+        //
+        {
+            const uint64_t payLoad = (pp.size() * sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::OFFSETS)
+                .closeTag();
+
+            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().tag(vtk::fileTag::CELL_DATA);
+
+        // zone Id
+        {
+            const uint64_t payLoad = (pp.size() * sizeof(label));
+
+            format().openDataArray<label>("zoneId")
+                .closeTag();
+
+            format().writeSize(payLoad);
+
+            forAll(pp, facei)
+            {
+                format().write(zoneI);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        format().endTag(vtk::fileTag::CELL_DATA);
+
+        format().endTag(vtk::fileTag::PIECE);
+    }
+
+    format().endTag(vtk::fileTag::POLY_DATA);
+    format().endVTKFile();
+}
+
+
+void Foam::lumpedPointMovement::writeVTP
+(
+    const fileName& file,
+    const polyMesh& mesh,
+    const labelUList& patchIds,
+    const pointField& points0
+) const
+{
+    writeVTP(file, state(), mesh, patchIds, points0);
+}
+
+
+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();
+
+    OFstream fos(file);
+    std::ostream& os = fos.stdStream();
+
+    autoPtr<vtk::formatter> format = vtk::newFormatter
+    (
+        os,
+        vtk::formatType::INLINE_ASCII
+    );
+
+    format().xmlHeader()
+        .beginVTKFile(vtk::fileTag::POLY_DATA, "0.1");
+
+    for (const label patchId : patchIds)
+    {
+        const polyPatch& pp = boundaryMesh[patchId];
+
+        format()
+            .openTag(vtk::fileTag::PIECE)
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_POINTS, pp.nPoints())
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_POLYS,  pp.size())
+            .closeTag();
+
+        // 'points'
+        {
+            const uint64_t payLoad = (pp.nPoints()*3* sizeof(float));
+
+            format()
+                .tag(vtk::fileTag::POINTS)
+                .openDataArray<float, 3>(vtk::dataArrayAttr::POINTS)
+                .closeTag();
+
+            // Could be more efficient, but not often needed
+            tmp<pointField> tpts = displacePoints
+            (
+                state,
+                points0,
+                pp.meshPoints()
+            ) + pointField(points0, pp.meshPoints());
+
+            const pointField& pts = tpts();
+
+            format().writeSize(payLoad);
+            vtk::writeList(format(), pts);
+            format().flush();
+
+            format()
+                .endDataArray()
+                .endTag(vtk::fileTag::POINTS);
+        }
+
+        // <Polys>
+        format().tag(vtk::fileTag::POLYS);
+
+
+        //
+        // 'connectivity'
+        //
+        {
+            uint64_t payLoad = 0;
+            for (const face& f : pp)
+            {
+                payLoad += f.size();
+            }
+
+            format().openDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY)
+                .closeTag();
+
+            format().writeSize(payLoad * sizeof(label));
+
+            for (const face& f : pp.localFaces())
+            {
+                vtk::writeList(format(), f);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        //
+        // 'offsets'  (connectivity offsets)
+        //
+        {
+            const uint64_t payLoad = (pp.size() * sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::OFFSETS)
+                .closeTag();
+
+            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().endTag(vtk::fileTag::PIECE);
+    }
+
+    format().endTag(vtk::fileTag::POLY_DATA);
+    format().endVTKFile();
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointState.C

@@ -0,0 +1,362 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "lumpedPointState.H"
+#include "demandDrivenData.H"
+#include "EulerCoordinateRotation.H"
+#include "unitConversion.H"
+
+#include "ISstream.H"
+#include "IFstream.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::Enum<Foam::lumpedPointState::inputFormatType>
+    Foam::lumpedPointState::formatNames
+    {
+        { inputFormatType::PLAIN, "plain" },
+        { inputFormatType::DICTIONARY, "dictionary" }
+    };
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+//! \cond fileScope
+static Foam::string getLineNoComment(Foam::ISstream& is)
+{
+    Foam::string line;
+    do
+    {
+        is.getLine(line);
+    }
+    while ((line.empty() || line[0] == '#') && is.good());
+
+    return line;
+}
+//! \endcond
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::lumpedPointState::calcRotations() const
+{
+    rotationPtr_ = new tensorField(angles_.size());
+    forAll(angles_, itemi)
+    {
+        rotationPtr_->operator[](itemi) = EulerCoordinateRotation
+        (
+            angles_[itemi],
+            degrees_  // true=degrees, false=radians
+        ).R();
+    }
+}
+
+
+void Foam::lumpedPointState::readDict(const dictionary& dict)
+{
+    dict.lookup("points") >> points_;
+    dict.lookup("angles") >> angles_;
+    degrees_ = dict.lookupOrDefault("degrees", false);
+    deleteDemandDrivenData(rotationPtr_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::lumpedPointState::lumpedPointState()
+:
+    points_(0),
+    angles_(0),
+    degrees_(false),
+    rotationPtr_(nullptr)
+{}
+
+
+Foam::lumpedPointState::lumpedPointState(const lumpedPointState& rhs)
+:
+    points_(rhs.points_),
+    angles_(rhs.angles_),
+    degrees_(rhs.degrees_),
+    rotationPtr_(nullptr)
+{}
+
+
+Foam::lumpedPointState::lumpedPointState
+(
+    const pointField& pts
+)
+:
+    points_(pts),
+    angles_(points_.size(), Zero),
+    degrees_(false),
+    rotationPtr_(nullptr)
+{}
+
+
+Foam::lumpedPointState::lumpedPointState
+(
+    tmp<pointField>& pts
+)
+:
+    points_(pts),
+    angles_(points_.size(), Zero),
+    degrees_(false),
+    rotationPtr_(nullptr)
+{}
+
+
+Foam::lumpedPointState::lumpedPointState
+(
+    const dictionary& dict
+)
+:
+    points_(0),
+    angles_(0),
+    degrees_(false),
+    rotationPtr_(nullptr)
+{
+    readDict(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+
+Foam::lumpedPointState::~lumpedPointState()
+{
+    deleteDemandDrivenData(rotationPtr_);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::lumpedPointState::operator=(const lumpedPointState& rhs)
+{
+    points_  = rhs.points_;
+    angles_  = rhs.angles_;
+    degrees_ = rhs.degrees_;
+
+    deleteDemandDrivenData(rotationPtr_);
+}
+
+
+void Foam::lumpedPointState::relax
+(
+    const scalar alpha,
+    const lumpedPointState& prev
+)
+{
+    points_ = prev.points_ + alpha*(points_ - prev.points_);
+
+    scalar convert = 1.0;
+    if (degrees_ != prev.degrees_)
+    {
+        if (prev.degrees_)
+        {
+            // Was degrees, now radians
+            convert = degToRad(1);
+        }
+        else
+        {
+            // Was radians, now degrees
+            convert = radToDeg(1);
+        }
+    }
+
+    angles_ = convert*prev.angles_ + alpha*(angles_ - convert*prev.angles_);
+
+    deleteDemandDrivenData(rotationPtr_);
+}
+
+
+bool Foam::lumpedPointState::readPlain(Istream& is)
+{
+    // Assume generic input stream so we can do line-based input
+    ISstream& iss = dynamic_cast<ISstream&>(is);
+
+    string line = getLineNoComment(iss);
+
+    label count;
+    {
+        IStringStream isstr(line);
+        isstr >> count;
+    }
+
+    points_.setSize(count);
+    angles_.setSize(count);
+
+    count = 0;
+    forAll(points_, i)
+    {
+        iss.getLine(line);
+        IStringStream isstr(line);
+
+        isstr
+            >> points_[count].x() >> points_[count].y() >> points_[count].z()
+            >> angles_[count].x() >> angles_[count].y() >> angles_[count].z();
+
+        ++count;
+    }
+
+    points_.setSize(count);
+    angles_.setSize(count);
+
+    degrees_ = false;
+    deleteDemandDrivenData(rotationPtr_);
+
+    return count;
+}
+
+
+bool Foam::lumpedPointState::readData(Istream& is)
+{
+    dictionary dict(is);
+    readDict(dict);
+
+    return points_.size();
+}
+
+
+bool Foam::lumpedPointState::writeData(Ostream& os) const
+{
+    writeDict(os);
+    return true;
+}
+
+
+void Foam::lumpedPointState::writeDict(Ostream& os) const
+{
+    os.writeEntry("points", points_);
+    os.writeEntry("angles", angles_);
+    if (degrees_)
+    {
+        os.writeKeyword("degrees") << "true;" << nl;
+    }
+}
+
+
+void Foam::lumpedPointState::writePlain(Ostream& os) const
+{
+    os  <<"# input for OpenFOAM\n"
+        <<"# N, points, angles\n"
+        << points_.size() << "\n";
+
+    forAll(points_, i)
+    {
+        const vector& pt = points_[i];
+
+        os  << pt.x() << ' '
+            << pt.y() << ' '
+            << pt.z() << ' ';
+
+        if (i < angles_.size())
+        {
+            os  << angles_[i].x() << ' '
+                << angles_[i].y() << ' '
+                << angles_[i].z() << '\n';
+        }
+        else
+        {
+            os  << "0 0 0\n";
+        }
+    }
+}
+
+
+bool Foam::lumpedPointState::readData
+(
+    const inputFormatType& fmt,
+    const fileName& file
+)
+{
+    bool ok = false;
+    if (Pstream::master())
+    {
+        IFstream is(file);
+
+        if (fmt == inputFormatType::PLAIN)
+        {
+            ok = this->readPlain(is);
+        }
+        else
+        {
+            ok = this->readData(is);
+        }
+    }
+
+    if (Pstream::parRun())
+    {
+        // Scatter master data using communication scheme
+
+        const List<Pstream::commsStruct>& comms =
+        (
+            (Pstream::nProcs() < Pstream::nProcsSimpleSum)
+          ? Pstream::linearCommunication()
+          : Pstream::treeCommunication()
+        );
+
+        // Get my communication order
+        const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];
+
+        // Receive from up
+        if (myComm.above() != -1)
+        {
+            IPstream fromAbove
+            (
+                UPstream::commsTypes::scheduled,
+                myComm.above(),
+                0,
+                Pstream::msgType(),
+                Pstream::worldComm
+            );
+
+            fromAbove >> points_ >> angles_ >> degrees_;
+        }
+
+        // Send to downstairs neighbours
+        forAllReverse(myComm.below(), belowI)
+        {
+            OPstream toBelow
+            (
+                UPstream::commsTypes::scheduled,
+                myComm.below()[belowI],
+                0,
+                Pstream::msgType(),
+                Pstream::worldComm
+            );
+
+            toBelow << points_ << angles_ << degrees_;
+        }
+
+        deleteDemandDrivenData(rotationPtr_);
+
+        // MPI barrier
+        Pstream::scatter(ok);
+    }
+
+    return ok;
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointState.H

@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::lumpedPointState
+
+Description
+    The 'state' of lumped points corresponds to positions and rotations.
+
+    Encapsulates the response from the external application.
+    Entry place for applying relaxation and sub-stepping etc.
+
+SourceFiles
+    lumpedPointState.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lumpedPointState_H
+#define lumpedPointState_H
+
+#include "dictionary.H"
+#include "scalarList.H"
+#include "pointField.H"
+#include "scalarField.H"
+#include "vectorField.H"
+#include "tensorField.H"
+#include "Enum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class Istream;
+class Ostream;
+
+/*---------------------------------------------------------------------------*\
+                      Class lumpedPointState Declaration
+\*---------------------------------------------------------------------------*/
+
+//- Bundling of position/rotation
+class lumpedPointState
+{
+public:
+
+    //- Input format types
+    enum class inputFormatType
+    {
+        PLAIN,
+        DICTIONARY
+    };
+
+    // Static data
+
+    //- Names for the input format types
+    static const Enum<inputFormatType> formatNames;
+
+
+private:
+
+    // Private data
+
+        //- Positions of lumped points
+        pointField  points_;
+
+        //- Orientation of lumped points (as Euler angles)
+        vectorField angles_;
+
+        //- Euler angles in degrees instead radians
+        bool degrees_;
+
+        //- Tensor rotation of lumped points
+        mutable tensorField* rotationPtr_;
+
+
+    // Private Member Functions
+
+        void calcRotations() const;
+
+        void readDict(const dictionary& dict);
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        lumpedPointState();
+
+        //- Copy constructor
+        lumpedPointState(const lumpedPointState& rhs);
+
+        //- Construct from points with zero-rotation
+        lumpedPointState(const pointField& pts);
+
+        //- Construct from points with zero-rotation
+        lumpedPointState(tmp<pointField>& pts);
+
+        //- Construct from dictionary
+        lumpedPointState(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~lumpedPointState();
+
+
+    // Member Functions
+
+        //- Has positions and consistent number of rotations?
+        inline bool valid() const;
+
+        //- If no points were specified
+        inline bool empty() const;
+
+        //- The number of points
+        inline label size() const;
+
+        //- The points corresponding to mass centres
+        inline const pointField& points() const;
+
+        //- The orientation of the points (mass centres)
+        inline const vectorField& angles() const;
+
+        //- The local-to-global transformation for each point
+        inline const tensorField& rotations() const;
+
+        //- Relax the state
+        //  alpha = 1 : no relaxation
+        //  alpha < 1 : relaxation
+        //  alpha = 0 : do nothing
+        void relax(const scalar alpha, const lumpedPointState& prev);
+
+
+        //- Read input as dictionary content
+        bool readData(Istream& is);
+
+        //- Output as dictionary content
+        bool writeData(Ostream& os) const;
+
+        //- Output as dictionary content
+        void writeDict(Ostream& os) const;
+
+        //- Read input as plain content
+        bool readPlain(Istream& is);
+
+        //- 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);
+
+        //- 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;
+
+
+    // Member Operators
+
+        //- Assignment operator
+        void operator=(const lumpedPointState& rhs);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "lumpedPointStateI.H"
+
+#endif
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointStateI.H

@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+inline bool Foam::lumpedPointState::valid() const
+{
+    return points_.size() && points_.size() == angles_.size();
+}
+
+
+inline bool Foam::lumpedPointState::empty() const
+{
+    return points_.empty();
+}
+
+
+inline Foam::label Foam::lumpedPointState::size() const
+{
+    return points_.size();
+}
+
+
+inline const Foam::pointField& Foam::lumpedPointState::points() const
+{
+    return points_;
+}
+
+
+inline const Foam::vectorField& Foam::lumpedPointState::angles() const
+{
+    return angles_;
+}
+
+
+inline const Foam::tensorField& Foam::lumpedPointState::rotations() const
+{
+    if (!rotationPtr_)
+    {
+        calcRotations();
+    }
+
+    return *rotationPtr_;
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointStateWriter.C

@@ -0,0 +1,326 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "lumpedPointState.H"
+#include "OFstream.H"
+#include "axesRotation.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
+{
+    // local-to-global transformation tensors
+    const tensorField& localToGlobal = rotations();
+
+    OFstream fos(outputFile);
+    std::ostream& os = fos.stdStream();
+
+    autoPtr<vtk::formatter> format = vtk::newFormatter
+    (
+        os,
+        vtk::formatType::INLINE_ASCII
+    );
+
+    format().xmlHeader()
+        .beginVTKFile(vtk::fileTag::POLY_DATA, "0.1");
+
+    //
+    // The 'spine' of lumped mass points
+    //
+    {
+        format()
+            .openTag(vtk::fileTag::PIECE)
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_POINTS, points_.size())
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_VERTS,  points_.size())
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_LINES,  1)
+            .closeTag();
+
+        // 'points'
+        {
+            const uint64_t payLoad = (points_.size()*3* sizeof(float));
+
+            format()
+                .tag(vtk::fileTag::POINTS)
+                .openDataArray<float, 3>(vtk::dataArrayAttr::POINTS)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            vtk::writeList(format(), points_);
+            format().flush();
+
+            format()
+                .endDataArray()
+                .endTag(vtk::fileTag::POINTS);
+        }
+
+        // <Verts>
+        format().tag(vtk::fileTag::VERTS);
+
+        //
+        // 'connectivity'
+        //
+        {
+            const uint64_t payLoad = (points_.size()*sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            forAll(points_, i)
+            {
+                format().write(i);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        //
+        // 'offsets'  (connectivity offsets)
+        // = linear mapping onto points (with 1 offset)
+        //
+        {
+            const uint64_t payLoad = (points_.size()*sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::OFFSETS)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            forAll(points_, i)
+            {
+                format().write(i+1);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        format().endTag(vtk::fileTag::VERTS);
+        // </Verts>
+
+
+        // <Lines>
+        format().tag(vtk::fileTag::LINES);
+
+        //
+        // 'connectivity'
+        //
+        {
+            const uint64_t payLoad = (points_.size()*sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            forAll(points_, i)
+            {
+                format().write(i);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        //
+        // 'offsets'  (connectivity offsets)
+        // = single line
+        //
+        {
+            const uint64_t payLoad = (1*sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::OFFSETS)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            format().write(points_.size());
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        format().endTag(vtk::fileTag::LINES);
+        format().endTag(vtk::fileTag::PIECE);
+    }
+
+    // Standard corners in local axis
+    const axesRotation cornerTransform(axis);
+
+    FixedList<point, 4> corners;
+    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()
+            .openTag(vtk::fileTag::PIECE)
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_POINTS, nPoints)
+            .xmlAttr(vtk::fileAttr::NUMBER_OF_POLYS,  nPolys)
+            .closeTag();
+
+        // 'points'
+        {
+            const uint64_t payLoad = (nPoints*3*sizeof(float));
+
+            format()
+                .tag(vtk::fileTag::POINTS)
+                .openDataArray<float, 3>(vtk::dataArrayAttr::POINTS)
+                .closeTag();
+
+            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)
+        //
+        {
+            const uint64_t payLoad = (4*nPolys*sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            for (label i=0; i < 4*nPolys; ++i)
+            {
+                format().write(i);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        //
+        // 'offsets'  (connectivity offsets)
+        // = single quad
+        //
+        {
+            const uint64_t payLoad = (nPolys*sizeof(label));
+
+            format().openDataArray<label>(vtk::dataArrayAttr::OFFSETS)
+                .closeTag();
+
+            format().writeSize(payLoad);
+            for (label i=0; i < nPolys; ++i)
+            {
+                const label off = 4 * (i+1);
+                format().write(off);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        format().endTag(vtk::fileTag::POLYS);
+
+#if 0
+        format().tag(vtk::fileTag::CELL_DATA);
+
+        // zone Id
+        {
+            const uint64_t payLoad = (points_.size()*sizeof(label));
+
+            format().openDataArray<label>("zoneId")
+                .closeTag();
+
+            format().writeSize(payLoad);
+            for (label i=0; i < nPolys; ++i)
+            {
+                format().write(i);
+            }
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        format().endTag(vtk::fileTag::CELL_DATA);
+#endif
+
+        format().endTag(vtk::fileTag::PIECE);
+    }
+
+    // Finally
+    // could add a 'ghost' level above to visualize extrapolated values
+    // draw as two triangles to distingush from real levels ...
+
+    format().endTag(vtk::fileTag::POLY_DATA);
+    format().endVTKFile();
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointTools.C

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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 "lumpedPointTools.H"
+#include "IFstream.H"
+#include "IOobjectList.H"
+#include "volFields.H"
+
+#include "lumpedPointDisplacementPointPatchVectorField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    // file-scope
+    template<class GeoFieldType>
+    static autoPtr<GeoFieldType> loadPointField
+    (
+        const pointMesh::Mesh& mesh,
+        const IOobject* io
+    )
+    {
+        if (io && io->headerClassName() == GeoFieldType::typeName)
+        {
+            Info<< "Reading " << GeoFieldType::typeName
+                << ' ' << io->name() << endl;
+
+            return autoPtr<GeoFieldType>
+            (
+                new GeoFieldType
+                (
+                    IOobject
+                    (
+                        io->name(),
+                        io->instance(),
+                        io->local(),
+                        io->db(),
+                        IOobject::MUST_READ,
+                        IOobject::AUTO_WRITE,
+                        io->registerObject()
+                    ),
+                    mesh
+                )
+            );
+        }
+
+        return autoPtr<GeoFieldType>();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::List<Foam::lumpedPointStateTuple>
+Foam::lumpedPointTools::lumpedPointStates(Istream& is)
+{
+    dictionary contents(is);
+    List<dictionary> entries(contents.lookup("response"));
+
+    DynamicList<Tuple2<scalar, lumpedPointState>> states(entries.size());
+
+    for (const dictionary& dict : entries)
+    {
+        states.append
+        (
+            lumpedPointStateTuple
+            (
+                readScalar(dict.lookup("time")),
+                lumpedPointState(dict)
+            )
+        );
+    }
+
+    return states.shrink();
+}
+
+
+Foam::List<Foam::lumpedPointStateTuple>
+Foam::lumpedPointTools::lumpedPointStates(const fileName& file)
+{
+    IFstream is(file);
+    return lumpedPointStates(is);
+}
+
+
+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;
+}
+
+
+
+Foam::labelList
+Foam::lumpedPointTools::lumpedPointPatchList(const pointVectorField& pvf)
+{
+    return lumpedPointDisplacementPointPatchVectorField::patchIds(pvf);
+}
+
+
+Foam::labelList Foam::lumpedPointTools::lumpedPointPatchList
+(
+    const polyMesh& mesh
+)
+{
+    IOobjectList objects0(mesh, "0");
+
+    pointMesh pMesh(mesh);
+
+    autoPtr<pointVectorField> displacePtr = loadPointField<pointVectorField>
+    (
+        pMesh,
+        objects0.lookup("pointDisplacement")
+    );
+
+    if (!displacePtr.valid())
+    {
+        Info<< "no valid pointDisplacement" << endl;
+        return labelList();
+    }
+
+    return lumpedPointPatchList(displacePtr());
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointTools.H

@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::lumpedPointTools
+
+Description
+    A collection of utility functions for handling IO related to the
+    lumped-mass movement.
+
+SourceFiles
+    lumpedPointTools.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lumpedPointTools_H
+#define lumpedPointTools_H
+
+#include "labelList.H"
+#include "polyMesh.H"
+#include "pointMesh.H"
+#include "pointFields.H"
+#include "Tuple2.H"
+
+#include "lumpedPointState.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+typedef Tuple2<scalar, lumpedPointState> lumpedPointStateTuple;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace lumpedPointTools
+{
+
+//- Load a list of states from an Istream
+List<lumpedPointStateTuple> lumpedPointStates(Istream& is);
+
+//- Load a list of states from a file
+List<lumpedPointStateTuple> lumpedPointStates(const fileName& file);
+
+
+//- Return the 0 or constant points field
+pointIOField points0Field(const polyMesh& mesh);
+
+//- Return the patch-ids associated 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
+labelList lumpedPointPatchList(const polyMesh& mesh);
+
+} // End namespace lumpedPointTools
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //