ENH: lumped point motion using local linear basic functions (#1341)

- the earlier implementation of externally controlled lumped point
  motion (see merge request !120 and OpenFOAM-v1706 release notes) was
  conceived for the motion of simple structures such as buildings or
  simple beams. The motion controller was simply defined in terms of
  an orientation axis and divisions along that axis.

  To include complex structures, multiple motion controllers are
  defined in terms of support points and connectivity.

  The points can have additional node Ids associated with them, which
  makes it easier to map to/from FEA models.

  OLD system/lumpedPointMovement specification
  --------------------------------------------

      //- Reference axis for the locations
      axis            (0 0 1);

      //- Locations of the lumped points
      locations       (0 0.05 .. 0.5);

  NEW system/lumpedPointMovement specification
  --------------------------------------------

      // Locations of the lumped points
      points
      (
          (0  0  0.00)
          (0  0  0.05)
          ...
          (0  0  0.50)
      );

      //- Connectivity for motion controllers
      controllers
      {
          vertical
          {
              pointLabels (0 1 2 3 4 5 6 7 8 9 10);
          }
      }

  And the controller(s) must be associated with the given
  pointDisplacement patch. Eg,

     somePatch
     {
         type            lumpedPointDisplacement;
         value           uniform (0 0 0);
         controllers     ( vertical );   // <-- NEW
     }

TUT: adjust building motion tutorial

- use new controllor definitions
- replace building response file with executable
- add updateControl in dynamicMeshDict for slowly moving structure

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,8 +27,8 @@ Application
     lumpedPointForces
 
 Description
-    Extract force/moment information from existing calculations based
+    Extract force/moment information from simulation results that
-    on the segmentation of the pressure integration zones.
+    use the lumped points movement description.
 
 \*---------------------------------------------------------------------------*/
 
@@ -37,7 +37,7 @@ Description
 #include "timeSelector.H"
 #include "volFields.H"
 #include "IOobjectList.H"
-
+#include "foamVtkSeriesWriter.H"
 #include "lumpedPointTools.H"
 #include "lumpedPointIOMovement.H"
 
@@ -83,8 +83,8 @@ int main(int argc, char *argv[])
 {
     argList::addNote
     (
-        "Extract force/moment information from existing calculations based"
+        "Extract force/moment information from simulation results that"
-        " on the lumpedPoints pressure zones."
+        " use the lumped points movement description."
     );
 
     argList::addBoolOption
@@ -102,31 +102,33 @@ int main(int argc, char *argv[])
     const bool withVTK = args.found("vtk");
 
     #include "createTime.H"
+
     instantList timeDirs = timeSelector::select0(runTime, args);
 
     #include "createNamedMesh.H"
 
-    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New
+    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New(mesh);
-    (
-        runTime
-    );
 
-    if (!movement.valid())
+    if (!movement)
     {
-        Info<< "no valid movement given" << endl;
+        Info<< "No valid movement found" << endl;
         return 1;
     }
 
-    const labelList patchLst = lumpedPointTools::lumpedPointPatchList(mesh);
+    const label nPatches = lumpedPointTools::setPatchControls(mesh);
-    if (patchLst.empty())
+    if (!nPatches)
     {
-        Info<< "no patch list found" << endl;
+        Info<< "No point patches with lumped movement found" << endl;
         return 2;
     }
 
-    movement().setMapping(mesh, patchLst, lumpedPointTools::points0Field(mesh));
+    Info<<"Lumped point patch controls set on " << nPatches
+        << " patches" << nl;
 
+
+    vtk::seriesWriter forceSeries;
     List<vector> forces, moments;
+
     forAll(timeDirs, timei)
     {
         runTime.setTime(timeDirs[timei], timei);
@@ -164,11 +166,21 @@ int main(int argc, char *argv[])
                     forces,
                     moments
                 );
+
+                forceSeries.append(runTime.timeIndex(), outputName);
             }
         }
     }
 
-    Info<< "End\n" << endl;
+
+    // Create file series
+
+    if (forceSeries.size())
+    {
+        forceSeries.write("forces.vtp");
+    }
+
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }

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

@@ -40,15 +40,50 @@ Description
 
 #include "argList.H"
 #include "Time.H"
-#include "timeSelector.H"
-
 #include "OFstream.H"
-
+#include "foamVtkSeriesWriter.H"
 #include "lumpedPointTools.H"
 #include "lumpedPointIOMovement.H"
 
 using namespace Foam;
 
+
+inline List<lumpedPointStateTuple> getResponseTable
+(
+    const fileName& file,
+    const lumpedPointState& state0
+)
+{
+    return lumpedPointTools::lumpedPointStates
+    (
+        file,
+        state0.rotationOrder(),
+        state0.degrees()
+    );
+}
+
+
+void echoTableLimits
+(
+    const List<lumpedPointStateTuple>& tbl,
+    const label span,
+    const label maxOut
+)
+{
+    Info<< "Using response table with " << tbl.size() << " entries" << nl;
+
+    if (span)
+    {
+        Info<< "Increment input by " << span << nl;
+    }
+
+    if (maxOut)
+    {
+        Info<< "Stopping after " << maxOut << " outputs" << nl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -59,7 +94,6 @@ int main(int argc, char *argv[])
         " for diagnostic purposes."
     );
 
-    argList::noParallel();
     argList::noFunctionObjects();  // Never use function objects
     argList::addOption
     (
@@ -71,7 +105,7 @@ int main(int argc, char *argv[])
     (
         "span",
         "N",
-        "Increment each input by factor N (default: 1)"
+        "Increment each input by N (default: 1)"
     );
     argList::addOption
     (
@@ -79,6 +113,17 @@ int main(int argc, char *argv[])
         "factor",
         "Relaxation/scaling factor for movement (default: 1)"
     );
+    argList::addOption
+    (
+        "visual-length",
+        "len",
+        "Visualization length for planes (visualized as triangles)"
+    );
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Test movement without a mesh"
+    );
     argList::addBoolOption
     (
         "removeLock",
@@ -96,45 +141,70 @@ int main(int argc, char *argv[])
     const label maxOut = Foam::max(0, args.getOrDefault<label>("max", 0));
     const label span   = Foam::max(1, args.getOrDefault<label>("span", 1));
 
-    const scalar relax = args.getOrDefault<scalar>("scale", 1);
+    // Control parameters
-
+    const bool dryrun = args.found("dry-run");
     const bool slave = args.found("slave");
     const bool removeLock = args.found("removeLock");
 
-    #include "createTime.H"
+    const scalar relax = args.getOrDefault<scalar>("scale", 1);
 
-    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New
+    args.readIfPresent("visual-length", lumpedPointState::visLength);
-    (
-        runTime
-    );
 
-    if (!movement.valid())
+    const fileName responseFile(args[1]);
-    {
-        Info<< "no valid movement given" << endl;
-        return 1;
-    }
 
-    List<lumpedPointStateTuple> responseTable =
+    // ----------------------------------------------------------------------
-        lumpedPointTools::lumpedPointStates(args[1]);
+    // Slave mode
-
+    // ----------------------------------------------------------------------
-    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;
 
-        externalFileCoupler& coupler = movement().coupler();
+        if (Pstream::parRun())
+        {
+            FatalErrorInFunction
+                << "Running as slave responder is not permitted in parallel"
+                << nl
+                << exit(FatalError);
+        }
 
-        label count = 0;
+        #include "createTime.H"
-        for (label index = 0; index < responseTable.size(); index += span)
+
+        // Create movement without a mesh
+        autoPtr<lumpedPointIOMovement> movementPtr =
+            lumpedPointIOMovement::New(runTime);
+
+        if (!movementPtr)
+        {
+            Info<< "No valid movement found" << endl;
+            return 1;
+        }
+        auto& movement = *movementPtr;
+
+        // Reference state0
+        const lumpedPointState& state0 = movement.state0();
+
+        List<lumpedPointStateTuple> responseTable =
+            getResponseTable(responseFile, state0);
+
+        echoTableLimits(responseTable, span, maxOut);
+
+        if (dryrun)
+        {
+            Info<< "dry-run: response table with " << responseTable.size()
+                << " entries" << nl
+                << "\nEnd\n" << endl;
+            return 0;
+        }
+
+        externalFileCoupler& coupler = movement.coupler();
+
+        for
+        (
+            label timei = 0, outputCount = 0;
+            timei < responseTable.size();
+            timei += span
+        )
         {
             Info<< args.executable() << ": waiting for master" << endl;
 
@@ -146,14 +216,14 @@ int main(int argc, char *argv[])
                 break;
             }
 
-            lumpedPointState state = responseTable[index].second();
+            lumpedPointState state = responseTable[timei].second();
-            state.relax(relax, movement().state0());
+            state.relax(relax, state0);
 
             // Generate input for OpenFOAM
-            OFstream os(coupler.resolveFile(movement().inputName()));
+            OFstream os(coupler.resolveFile(movement.inputName()));
             if
             (
-                movement().inputFormat()
+                movement.inputFormat()
              == lumpedPointState::inputFormatType::PLAIN
             )
             {
@@ -161,101 +231,236 @@ int main(int argc, char *argv[])
             }
             else
             {
-                os.writeEntry("time", responseTable[index].first());
+                os.writeEntry("time", responseTable[timei].first());
                 state.writeDict(os);
             }
 
             Info<< args.executable()
-                << ": updated to state " << index
+                << ": updated to state " << timei
                 << " - switch to master"
                 << endl;
 
             // Let OpenFOAM know that it can continue
             coupler.useMaster();
 
-            if (maxOut && ++count >= maxOut)
+            ++outputCount;
+
+            if (maxOut && outputCount >= maxOut)
             {
                 Info<< args.executable()
-                    <<": stopping after " << maxOut << " outputs" << endl;
+                    << ": stopping after " << maxOut << " outputs" << endl;
                 break;
             }
         }
 
         if (removeLock)
         {
-            Info<< args.executable() <<": removing lock file" << endl;
+            Info<< args.executable() << ": removing lock file" << endl;
             coupler.useSlave();  // This removes the lock-file
         }
+
+        Info<< args.executable() << ": finishing" << nl;
+
+        Info<< "\nEnd\n" << endl;
+        return 0;
     }
-    else
+
+
+    // ----------------------------------------------------------------------
+    // dry-run
+    // ----------------------------------------------------------------------
+
+    if (dryrun)
     {
-        runTime.setTime(instant(0, runTime.constant()), 0);
+        Info<< "dry-run: creating states only" << nl;
 
-        #include "createNamedPolyMesh.H"
+        #include "createTime.H"
 
-        const labelList patchLst = lumpedPointTools::lumpedPointPatchList(mesh);
+        // Create movement without a mesh
-        if (patchLst.empty())
+        autoPtr<lumpedPointIOMovement> movementPtr =
+            lumpedPointIOMovement::New(runTime);
+
+        if (!movementPtr)
         {
-            Info<< "no patch list found" << endl;
+            Info<< "No valid movement found" << endl;
-            return 2;
+            return 1;
         }
+        auto& movement = *movementPtr;
 
-        pointIOField points0 = lumpedPointTools::points0Field(mesh);
+        // Reference state0
-        movement().setBoundBox(mesh, patchLst, points0);
+        const lumpedPointState& state0 = movement.state0();
 
-        label index = 0;
+        List<lumpedPointStateTuple> responseTable =
+            getResponseTable(responseFile, state0);
 
-        // Initial geometry
+        echoTableLimits(responseTable, span, maxOut);
-        movement().writeVTP("geom_init.vtp", mesh, patchLst, points0);
 
-        forAll(responseTable, i)
+
+        vtk::seriesWriter stateSeries;
+
+        for
+        (
+            label timei = 0, outputCount = 0;
+            timei < responseTable.size();
+            timei += span
+        )
         {
-            const bool output = ((i % span) == 0);
+            lumpedPointState state = responseTable[timei].second();
-            lumpedPointState state = responseTable[i].second();
-            state.relax(relax, movement().state0());
 
-            if (output)
+            state += movement.origin();
-            {
+            movement.scalePoints(state);
-                Info<<"output [" << i << "/"
+            state.relax(relax, state0);
-                    << responseTable.size() << "]" << endl;
+
-            }
+            Info<< "output [" << timei << '/' << responseTable.size() << ']';
-            else
-            {
-                continue;
-            }
 
             // State/response = what comes back from FEM
             {
-                const word outputName = word::printf("state_%06d.vtp", index);
+                const word outputName =
+                    word::printf("state_%06d.vtp", outputCount);
 
-                Info<<"    " << outputName << endl;
+                Info<< "  " << outputName;
 
-                state.writeVTP(outputName, movement().axis());
+                movement.writeStateVTP(state, outputName);
+                stateSeries.append(outputCount, outputName);
             }
 
+            Info<< endl;
+
+            ++outputCount;
+
+            if (maxOut && outputCount >= maxOut)
             {
-                const word outputName = word::printf("geom_%06d.vtp", index);
+                Info<< "Max output " << maxOut << " ... stopping" << endl;
-
+                break;
-                Info<<"    " << outputName << endl;
-
-                movement().writeVTP(outputName, state, mesh, patchLst, points0);
             }
+        }
 
-            {
+        // Write file series
-                ++index;
 
-                bool canOutput = !maxOut || (index <= maxOut);
+        if (stateSeries.size())
-                if (!canOutput)
+        {
-                {
+            Info<< nl << "write state.vtp.series" << nl;
-                    Info<<"stopping output after "
+            stateSeries.write("state.vtp");
-                        << maxOut << " outputs" << endl;
+        }
-                    break;
+
-                }
+        Info<< "\nEnd\n" << endl;
-            }
+        return 0;
+    }
+
+
+    // ----------------------------------------------------------------------
+    // test patch movement
+    // ----------------------------------------------------------------------
+
+    #include "createTime.H"
+
+    runTime.setTime(instant(runTime.constant()), 0);
+
+    #include "createNamedMesh.H"
+
+    // Create movement with mesh
+    autoPtr<lumpedPointIOMovement> movementPtr =
+        lumpedPointIOMovement::New(mesh);
+
+    if (!movementPtr)
+    {
+        Info<< "No valid movement found" << endl;
+        return 1;
+    }
+    auto& movement = *movementPtr;
+
+    // Reference state0
+    const lumpedPointState& state0 = movement.state0();
+
+    List<lumpedPointStateTuple> responseTable =
+        getResponseTable(responseFile, state0);
+
+    echoTableLimits(responseTable, span, maxOut);
+
+    pointIOField points0(lumpedPointTools::points0Field(mesh));
+
+    const label nPatches = lumpedPointTools::setPatchControls(mesh, points0);
+    if (!nPatches)
+    {
+        Info<< "No point patches with lumped movement found" << endl;
+        return 2;
+    }
+
+    Info<< "Lumped point patch controls set on "
+        << nPatches << " patches" << nl;
+
+    lumpedPointTools::setInterpolators(mesh, points0);
+
+
+    // Output vtk file series
+    vtk::seriesWriter stateSeries;
+    vtk::seriesWriter geomSeries;
+
+    // Initial geometry
+    movement.writeVTP("geom_init.vtp", state0, mesh, points0);
+
+    lumpedPointTools::setInterpolators(mesh);
+
+    for
+    (
+        label timei = 0, outputCount = 0;
+        timei < responseTable.size();
+        timei += span
+    )
+    {
+        lumpedPointState state = responseTable[timei].second();
+
+        state += movement.origin();
+        movement.scalePoints(state);
+        state.relax(relax, state0);
+
+        Info<< "output [" << timei << '/' << responseTable.size() << ']';
+
+        // State/response = what comes back from FEM
+        {
+            const word outputName =
+                word::printf("state_%06d.vtp", outputCount);
+
+            Info<< "  " << outputName;
+
+            movement.writeStateVTP(state, outputName);
+            stateSeries.append(outputCount, outputName);
+        }
+
+        {
+            const word outputName =
+                word::printf("geom_%06d.vtp", outputCount);
+
+            Info<< "  " << outputName;
+
+            movement.writeVTP(outputName, state, mesh, points0);
+            geomSeries.append(outputCount, outputName);
+        }
+
+        Info<< endl;
+
+        ++outputCount;
+
+        if (maxOut && outputCount >= maxOut)
+        {
+            Info<< "Max output " << maxOut << " ... stopping" << endl;
+            break;
         }
     }
 
-    Info<< args.executable() << ": finishing" << nl;
+
+    // Write file series
+
+    if (geomSeries.size())
+    {
+        Info<< nl << "write geom.vtp.series" << nl;
+        geomSeries.write("geom.vtp");
+    }
+    if (stateSeries.size())
+    {
+        Info<< nl << "write state.vtp.series" << nl;
+        stateSeries.write("state.vtp");
+    }
 
     Info<< "\nEnd\n" << endl;
 

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

@@ -2,9 +2,12 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/lumpedPointMotion/lnInclude
+    -I$(LIB_SRC)/lumpedPointMotion/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lfileFormats \
     -lmeshTools \
-    -llumpedPointMotion
+    -llumpedPointMotion \
+    -ldynamicMesh

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +39,7 @@ Description
 
 #include "lumpedPointTools.H"
 #include "lumpedPointIOMovement.H"
+#include "fvMesh.H"
 
 using namespace Foam;
 
@@ -52,9 +53,26 @@ int main(int argc, char *argv[])
         " pressure integration zones used by lumpedPoint BC."
     );
 
-    argList::noParallel();    // The VTP writer is not yet in parallel
     argList::noFunctionObjects();  // Never use function objects
 
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Test initial lumped points state without a mesh"
+    );
+    argList::addOption
+    (
+        "visual-length",
+        "len",
+        "Visualization length for planes (visualized as triangles)"
+    );
+
+    argList::addBoolOption
+    (
+        "no-interpolate",
+        "Suppress calculation/display of point interpolators"
+    );
+
     argList::addBoolOption
     (
         "verbose",
@@ -62,47 +80,89 @@ int main(int argc, char *argv[])
     );
 
     #include "addRegionOption.H"
+
     #include "setRootCase.H"
 
+    const bool noInterpolate = args.found("no-interpolate");
+
+    const bool dryrun = args.found("dry-run");
+
     // const bool verbose = args.found("verbose");
 
+    args.readIfPresent("visual-length", lumpedPointState::visLength);
+
     #include "createTime.H"
 
-    runTime.setTime(instant(0, runTime.constant()), 0);
+    if (dryrun)
+    {
+        // Create without a mesh
+        autoPtr<lumpedPointIOMovement> movement =
+            lumpedPointIOMovement::New(runTime);
 
-    #include "createNamedPolyMesh.H"
+        if (!movement.valid())
+        {
+            Info<< "No valid movement found" << endl;
+            return 1;
+        }
 
-    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New
+        const word outputName("state.vtp");
-    (
+
-        runTime
+        Info<< "dry-run: writing " << outputName << nl;
-    );
+
+        movement().writeStateVTP(movement().state0(), outputName);
+
+        Info<< "\nEnd\n" << endl;
+
+        return 0;
+    }
+
+
+    runTime.setTime(instant(runTime.constant()), 0);
+
+    #include "createNamedMesh.H"
+
+    autoPtr<lumpedPointIOMovement> movement = lumpedPointIOMovement::New(mesh);
 
     if (!movement.valid())
     {
-        Info<< "no valid movement found" << endl;
+        Info<< "No valid movement found" << endl;
         return 1;
     }
 
-    const labelList patchLst = lumpedPointTools::lumpedPointPatchList(mesh);
+    // Initial positions/rotations
-    if (patchLst.empty())
+    movement().writeStateVTP("state.vtp");
+
+    pointIOField points0(lumpedPointTools::points0Field(mesh));
+
+    const label nPatches = lumpedPointTools::setPatchControls(mesh, points0);
+    if (!nPatches)
     {
-        Info<< "no patch list found" << endl;
+        Info<< "No point patches with lumped movement found" << endl;
         return 2;
     }
 
-    pointIOField points0 = lumpedPointTools::points0Field(mesh);
+    Info<<"Lumped point patch controls set on "
-    movement().setMapping(mesh, patchLst, points0);
+        << nPatches << " patches" << nl;
 
-    // Initial geometry, but with zone colouring
+    Info<<"Areas per point: " << flatOutput(movement().areas(mesh)) << nl;
-    movement().writeZonesVTP("lumpedPointZones.vtp", mesh, points0);
 
-    // Initial positions/rotations
+    if (noInterpolate)
-    movement().writeStateVTP("initialState.vtp");
+    {
+        // Initial geometry, with zones
+        movement().writeZonesVTP("lumpedPointZones.vtp", mesh, points0);
+    }
+    else
+    {
+        lumpedPointTools::setInterpolators(mesh, points0);
+
+        // Initial geometry, with zones and interpolations
+        movement().writeVTP("lumpedPointZones.vtp", mesh, points0);
+    }
 
     Info<< nl
+        << "wrote 'state.vtp' (reference state)" << nl
         << "wrote 'lumpedPointZones.vtp'" << nl
-        << "wrote 'initialState.vtp'" << nl
+        << "\nEnd\n" << endl;
-        << "End\n" << endl;
 
     return 0;
 }

src/lumpedPointMotion/Make/files

@@ -1,3 +1,5 @@
+controller/lumpedPointController.C
+
 state/lumpedPointState.C
 state/lumpedPointStateWriter.C
 

src/lumpedPointMotion/controller/lumpedPointController.C

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "lumpedPointController.H"
+#include "dictionary.H"
+#include "labelField.H"
+#include "Map.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::lumpedPointController::lumpedPointController() noexcept
+:
+    pointLabels_()
+{}
+
+
+Foam::lumpedPointController::lumpedPointController
+(
+    const labelUList& pointLabels
+)
+:
+    pointLabels_(pointLabels)
+{}
+
+
+Foam::lumpedPointController::lumpedPointController
+(
+    labelList&& pointLabels
+)
+:
+    pointLabels_(std::move(pointLabels))
+{}
+
+
+Foam::lumpedPointController::lumpedPointController
+(
+    const dictionary& dict
+)
+:
+    pointLabels_(dict.get<labelList>("pointLabels"))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::lumpedPointController::remapPointLabels
+(
+    const label nPoints,
+    const Map<label>& originalIds
+)
+{
+    if (originalIds.size())
+    {
+        for (label& pointi : pointLabels_)
+        {
+            pointi = originalIds[pointi];
+        }
+    }
+
+    if (min(pointLabels_) < 0 || max(pointLabels_) >= nPoints)
+    {
+        FatalErrorInFunction
+            << "Point id out-of-range: " << flatOutput(pointLabels_) << nl
+            << exit(FatalError);
+    }
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/controller/lumpedPointController.H

@@ -1 +1,139 @@
-#warning File removed - left for old dependency check only
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::lumpedPointController
+
+Description
+    Simple connectivity of point labels to specify a controller for lumped
+    point movement.
+
+    \heading Dictionary parameters
+    \table
+        Property    | Description                           | Required | Default
+        pointLabels | List of point labels                  | yes |
+    \endtable
+
+Note
+    If the calling program itself specified a point-label mapping
+    (eg, original ids from FEA), these can be used initially and remapped.
+
+SourceFiles
+    lumpedPointController.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lumpedPointController_H
+#define lumpedPointController_H
+
+#include "List.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class dictionary;
+template<class T> class Map;
+
+/*---------------------------------------------------------------------------*\
+                    Class lumpedPointController Declaration
+\*---------------------------------------------------------------------------*/
+
+class lumpedPointController
+{
+    // Private Data
+
+        //- The lumped points associated with the controller
+        labelList pointLabels_;
+
+
+public:
+
+    // Constructors
+
+        //- Default construct
+        lumpedPointController() noexcept;
+
+        //- Copy construct from point ids
+        explicit lumpedPointController(const labelUList& pointLabels);
+
+        //- Move construct from point ids
+        explicit lumpedPointController(labelList&& pointLabels);
+
+        //- Construct from dictionary
+        explicit lumpedPointController(const dictionary& dict);
+
+
+        //- Factory method
+        static autoPtr<lumpedPointController> New(const dictionary& dict)
+        {
+            return autoPtr<lumpedPointController>::New(dict);
+        }
+
+
+    //- Destructor
+    ~lumpedPointController() = default;
+
+
+    // Member Functions
+
+        //- No controller points specified?
+        bool empty() const
+        {
+            return pointLabels_.empty();
+        }
+
+        //- Number of controller points specified
+        label size() const
+        {
+            return pointLabels_.size();
+        }
+
+        //- The controller points
+        const labelList& pointLabels() const
+        {
+            return pointLabels_;
+        }
+
+        //- Check point mapping (the count) or remap the point labels
+        void remapPointLabels
+        (
+            const label nPoints,
+            const Map<label>& originalIds
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lumpedPointMotion/interpolation/lumpedPointInterpolator.H

@@ -1 +1,160 @@
-#warning File removed - left for old dependency check only
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::lumpedPointInterpolator
+
+Description
+    A simple linear interpolator between two locations, which are
+    referenced by index.
+
+    When the interpolator is built for based on searching for nearest and
+    next-nearest points, the interpolation will typically cover a (0-0.5) range
+    rather than a (0-1) range. The (0.5-1) range implies a flip in the
+    relationship of nearest vs. next-nearest.
+
+SourceFiles
+    lumpedPointInterpolatorI.H
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lumpedPointInterpolator_H
+#define lumpedPointInterpolator_H
+
+#include "Pair.H"
+#include "triFace.H"
+#include "scalarList.H"
+#include "barycentric2D.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class lumpedPointInterpolator Declaration
+\*---------------------------------------------------------------------------*/
+
+class lumpedPointInterpolator
+{
+    // Private Data
+
+        //- The id of the nearest point
+        label nearest_;
+
+        //- The id of the neighbour point(s)
+        label next1_, next2_;
+
+        //- The interpolation weight for the neighbour point(s)
+        scalar weight1_, weight2_;
+
+
+public:
+
+    // Constructors
+
+        //- Default construct, with zero weighting and invalid ids
+        inline lumpedPointInterpolator();
+
+        //- Construct with nearest id
+        explicit inline lumpedPointInterpolator(const label id);
+
+
+    //- Destructor
+    ~lumpedPointInterpolator() = default;
+
+
+    // Member Functions
+
+    // Access
+
+        //- Valid if there is an associated nearest point
+        inline bool valid() const;
+
+        //- The nearest control point, or -1 if invalid
+        inline label nearest() const;
+
+        //- The first neighbour control point - identical to next1()
+        inline label next() const;
+
+        //- The first neighbour control point, or -1 if invalid
+        inline label next1() const;
+
+        //- The second neighbour control point, or -1 if invalid
+        inline label next2() const;
+
+        //- The weighting for the nearest point
+        inline scalar weight0() const;
+
+        //- The weighting for the first neighbour point,
+        //- this also corresponds to the logical location (interval 0-1)
+        inline scalar weight1() const;
+
+        //- The weighting for the second neighbour point,
+        //- this also corresponds to the logical location (interval 0-1)
+        inline scalar weight2() const;
+
+
+    // Edit
+
+        //- Assign the nearest point, clearing any neighbour
+        inline void nearest(const label id);
+
+        //- Assign single neighbour control point and its weight
+        inline void next(const label id, const scalar weight);
+
+        //- Assign the neighbour control point and its weight
+        inline void next
+        (
+            const label id,
+            const scalar weight,
+            const label position    //!< Use 0 or 1 for first/second
+        );
+
+        //- Assign all control points and their weights
+        //  The triFace points [0,1,2] correspond to [nearest,next1,next2],
+        //  respectively
+        inline void set(const triFace& ids, const barycentric2D& weights);
+
+
+    // Evalulate
+
+        //- Linear interpolated value between nearest and next locations
+        template<class T>
+        inline T interpolate(const UList<T>& input) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "lumpedPointInterpolatorI.H"
+
+#endif
+
+// ************************************************************************* //

src/lumpedPointMotion/interpolation/lumpedPointInterpolatorI.H

@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+inline Foam::lumpedPointInterpolator::lumpedPointInterpolator()
+:
+    lumpedPointInterpolator(-1)
+{}
+
+
+inline Foam::lumpedPointInterpolator::lumpedPointInterpolator(const label id)
+:
+    nearest_(id),
+    next1_(id),
+    next2_(id),
+    weight1_(Zero),
+    weight2_(Zero)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline bool Foam::lumpedPointInterpolator::valid() const
+{
+    return nearest_ != -1;
+}
+
+
+inline Foam::label Foam::lumpedPointInterpolator::nearest() const
+{
+    return nearest_;
+}
+
+
+inline Foam::label Foam::lumpedPointInterpolator::next() const
+{
+    return next1_;
+}
+
+
+inline Foam::label Foam::lumpedPointInterpolator::next1() const
+{
+    return next1_;
+}
+
+
+inline Foam::label Foam::lumpedPointInterpolator::next2() const
+{
+    return next2_;
+}
+
+
+inline Foam::scalar Foam::lumpedPointInterpolator::weight0() const
+{
+    return scalar(1) - weight1_ - weight2_;
+}
+
+
+inline Foam::scalar Foam::lumpedPointInterpolator::weight1() const
+{
+    return weight1_;
+}
+
+
+inline Foam::scalar Foam::lumpedPointInterpolator::weight2() const
+{
+    return weight2_;
+}
+
+
+inline void Foam::lumpedPointInterpolator::nearest(const label id)
+{
+    nearest_ = id;
+    next(id, Zero);
+}
+
+
+inline void Foam::lumpedPointInterpolator::next
+(
+    const label id,
+    scalar weight
+)
+{
+    if (weight < scalar(0))
+    {
+        weight = 0;
+    }
+    else if (weight > scalar(1))
+    {
+        weight = 1;
+    }
+
+    next1_ = id;
+    next2_ = id;
+    weight1_ = weight;
+    weight2_ = Zero;
+}
+
+
+inline void Foam::lumpedPointInterpolator::next
+(
+    const label id,
+    scalar weight,
+    const label position
+)
+{
+    if (weight < scalar(0))
+    {
+        weight = 0;
+    }
+    else if (weight > scalar(1))
+    {
+        weight = 1;
+    }
+
+    if (!position)
+    {
+        next1_ = id;
+        weight1_ = weight;
+    }
+    else
+    {
+        next2_ = id;
+        weight2_ = weight;
+    }
+}
+
+
+inline void Foam::lumpedPointInterpolator::set
+(
+    const triFace& ids,
+    const barycentric2D& weights
+)
+{
+    nearest_ = ids[0];
+    next1_ = ids[1];
+    next2_ = ids[2];
+
+    weight1_ = weights[1];
+    weight2_ = weights[2];
+}
+
+
+template<class T>
+inline T Foam::lumpedPointInterpolator::interpolate(const UList<T>& input) const
+{
+    if (nearest_ == -1)
+    {
+        return Zero;
+    }
+    else if (next1_ == -1 || next1_ == nearest_)
+    {
+        return input[nearest_];
+    }
+
+    return
+    (
+        input[nearest_] * (1-weight1_-weight2_)
+      + input[next1_] * (weight1_)
+      + input[next2_] * (weight2_)
+    );
+}
+
+
+// ************************************************************************* //

src/lumpedPointMotion/lumpedPointMovementDict

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v1912                                 |
+|  \\    /   O peration     | Version:  v2006                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -15,36 +15,59 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Reference axis for the locations
+//- The initial locations of the lumped points
-axis            (0 0 1);
+points
+(
 
-// Locations of the lumped points
+);
-locations       11(0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5);
 
-// Division for pressure forces (0-1)
+//- The FEA point labels for the points (optional)
-division        0.5;
+// pointLabels
+// (
+//
+// );
 
-//- If present, the offset of patch points compared to the locations
+//- The motion controllers.
-//  Otherwise determined from the bounding box
+//  The pointLabels are an adjacency list of point connectivity.
-// centre          (0 0 0);
+//
+//  Controllers are joined by common point labels between lists
+//  (or internally - eg, ring structure) when they appear in the corresponding
+//  point displacement patch
+controllers
+{
+    control1
+    {
+        pointLabels (0 1 2 3 4);
+    }
+    control2
+    {
+        pointLabels (1 5 6);
+    }
+}
 
-//- The interpolation scheme
-interpolationScheme linear;
 
-//- Relaxation/scaling factor when updating positions
+//- Input offset to shift ALL lumped points when reading (default: 0 0 0)
-relax           1.0;
+// origin  (0 0 0);
 
+//- The Euler rotation order (default: zxz)
+// rotationOrder zxz;
+
+//- Input rotations in degrees instead of radians (default: false)
+// degrees false;
+
+//- Relaxation/scaling factor when updating positions (default: 1)
+// relax   1.0;
 
 forces
 {
-    //- The pressure name (default: p)
+    //- Pressure name (default: p)
-    p           p;
+    p       p;
 
     //- Reference pressure [Pa] (default: 0)
-    pRef        0;
+    pRef    0;
 
     //- Reference density for incompressible calculations (default: 1)
-    rhoRef      1;
+    rhoRef  1;
 }
 
 
@@ -60,6 +83,9 @@ communication
 
     initByExternal  false;
 
+    // Coupling frequency in time-steps (default: 1-steps)
+    // calcFrequency   1;
+
     // Input file of positions/rotation, written by external application
     inputName       positions.in;
 
@@ -70,26 +96,28 @@ communication
     logName         movement.log;
 
     inputFormat     dictionary;
+
     outputFormat    dictionary;
 
-    // Scaling applied to values read from 'inputName'
+    // Scaling applied to values read from 'inputName',
+    // also applies to points (above)
     scaleInput
     {
         //- Length multiplier (to metres). Eg 0.001 for [mm] -> [m]
-        length      1;
+        length  1;
     }
 
     // Scaling applied to values written to 'outputName'
     scaleOutput
     {
         //- Length multiplier (from metres). Eg 1000 for [m] -> [mm]
-        length      1;
+        length  1;
 
         //- Force units multiplier (from Pa)
-        force       1;
+        force   1;
 
         //- Moment units multiplier (from N.m)
-        moment      1;
+        moment  1;
     }
 }
 

src/lumpedPointMotion/movement/lumpedPointIOMovement.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,16 +32,16 @@ License
 
 namespace Foam
 {
-    defineTypeNameAndDebug(lumpedPointIOMovement, 0);
+    defineTypeName(lumpedPointIOMovement);
 }
 
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-const Foam::lumpedPointIOMovement*
+Foam::lumpedPointIOMovement*
-Foam::lumpedPointIOMovement::lookupInRegistry(const objectRegistry& obr)
+Foam::lumpedPointIOMovement::getMovementObject(const objectRegistry& obr)
 {
-    return obr.findObject<lumpedPointIOMovement>
+    return obr.getObjectPtr<lumpedPointIOMovement>
     (
         lumpedPointMovement::canonicalName
     );
@@ -120,7 +120,7 @@ bool Foam::lumpedPointIOMovement::writeData(Ostream& os) const
 }
 
 
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const lumpedPointIOMovement& obj)
 {

src/lumpedPointMotion/movement/lumpedPointIOMovement.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,11 +46,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
-class lumpedPointIOMovement;
-
-Ostream& operator<<(Ostream& os, const lumpedPointIOMovement& obj);
-
 /*---------------------------------------------------------------------------*\
                     Class lumpedPointIOMovement Declaration
 \*---------------------------------------------------------------------------*/
@@ -68,22 +63,23 @@ class lumpedPointIOMovement
         //- No copy assignment
         void operator=(const lumpedPointIOMovement&) = delete;
 
+
 public:
 
     //- Runtime type information
-    TypeName("lumpedPointMovement");
+    TypeNameNoDebug("lumpedPointMovement");
 
 
     // Static Member Functions
 
-        //- Lookup pointer to object in the object-registry,
+        //- Find the movement object or nullptr if not found
-        //  return nullptr if found.
+        static lumpedPointIOMovement* getMovementObject
-        static const lumpedPointIOMovement* lookupInRegistry
         (
             const objectRegistry& obr
         );
 
-        //- Create a new object in the registry by reading system dictionary
+        //- Create a movement object in the registry by
+        //- reading system dictionary
         static autoPtr<lumpedPointIOMovement> New
         (
             const objectRegistry& obr,
@@ -94,36 +90,27 @@ public:
     // Constructors
 
         //- Construct from IOobject, optionally with some owner information
-        explicit lumpedPointIOMovement
+        explicit lumpedPointIOMovement(const IOobject& io, label ownerId = -1);
-        (
-            const IOobject& io,
-            label ownerId = -1
-        );
 
 
     //- Destructor
-    ~lumpedPointIOMovement() = default;
+    virtual ~lumpedPointIOMovement() = default;
 
 
     // Member Functions
 
-        //- readData member function used by regIOobject
+        //- The readData member function used by regIOobject
         bool readData(Istream& is);
 
-        //- writeData member function required by regIOobject
+        //- The writeData member function required by regIOobject
         bool writeData(Ostream& os) const;
-
-
-    // IOstream Operators
-
-        friend Ostream& operator<<
-        (
-            Ostream& os,
-            const lumpedPointIOMovement& obj
-        );
 };
 
 
+// IOstream Operators
+Ostream& operator<<(Ostream& os, const lumpedPointIOMovement& obj);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lumpedPointMotion/movement/lumpedPointMovement.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2018 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,9 +27,9 @@ Class
     Foam::lumpedPointMovement
 
 Description
-    The movement \a driver that describes initial point locations, the
+    The movement \a driver that describes initial point locations,
-    segmentation for pressure integration, the current state of the
+    the current state of the points/rotations,
-    points/rotations, and forwarding to the externalFileCoupler
+    and forwarding to the externalFileCoupler
     communication coordinator.
 
     The lumpedPointIOMovement class is simply a registered version of
@@ -41,14 +41,15 @@ Description
     \heading Dictionary parameters
     \table
         Property    | Description                           | Required | Default
-        axis        | Reference axis for the locations      | yes |
+        points      | Initial locations of lumped points    | yes |
-        locations   | List of lumped point locations        | yes |
+        pointLabels | The FEA ids for the points            | no  |
-        centre      | Offset of patch points to locations   | no | automatic
+        origin      | Shift offset when reading points      | no  | (0 0 0)
-        division    | Division (0-1) for pressure forces    | no | 0
+        rotationOrder | The Euler rotation order            | no  | zxz
-        relax       | Relaxation/scaling for updating positions | no  |
+        degrees     | Input rotations in degrees            | no  | false
-        interpolationScheme | The interpolation scheme      | yes |
+        relax       | Relaxation/scaling for updating positions | no  | 1
-        forces      | Optional forces dictionary            | no |
+        controllers | Motion controllers (dictionary)       | yes |
-        communication | Required communication dictionary   | yes |
+        forces      | Force settings (dictionary)           | no  |
+        communication | Communication settings (dictionary) | yes |
     \endtable
 
     \heading Parameters for communication dictionary
@@ -60,7 +61,8 @@ Description
         inputFormat | Input format: dictionary/plain    | yes |
         outputFormat | Output format: dictionary/plain  | yes |
         scaleInput  | Input scaling parameter dictionary | no |
-        scaleOutput | Output scaling parameter dictionary | no |
+        scaleOutput | Output scaling parameter dictionary | no  |
+        calcFrequency | Calculation/coupling frequency    | no  | 1
     \endtable
 
     \heading Parameters for optional communication/scaleInput dictionary
@@ -95,28 +97,25 @@ SourceFiles
 
 #include "dictionary.H"
 #include "scalarList.H"
-#include "scalarField.H"
+#include "primitiveFields.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 "HashPtrTable.H"
 #include "externalFileCoupler.H"
+#include "lumpedPointController.H"
+#include "lumpedPointInterpolator.H"
 #include "lumpedPointState.H"
-#include "boundBox.H"
 #include "Enum.H"
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class polyMesh;
+class polyPatch;
+class pointPatch;
 class Time;
 
 /*---------------------------------------------------------------------------*\
@@ -127,118 +126,116 @@ class lumpedPointMovement
 {
 public:
 
-    //- Output format types
+    // Data Types
-    enum class outputFormatType
-    {
-        PLAIN,      //!< "plain" is a simple ASCII format
-        DICTIONARY  //!< "dictionary" is the OpenFOAM dictionary format
-    };
 
-    //- Output format types
+        //- Output format types
-    enum scalingType
+        enum class outputFormatType
-    {
+        {
-        LENGTH = 0, //!< The "length" scaling
+            PLAIN,      //!< "plain" is a simple ASCII format
-        FORCE,      //!< The "force" scaling
+            DICTIONARY  //!< "dictionary" is the OpenFOAM dictionary format
-        MOMENT      //!< The "moment" scaling
+        };
-    };
 
-    // Static data
+        //- Output format types
+        enum scalingType
+        {
+            LENGTH = 0, //!< The "length" scaling
+            FORCE,      //!< The "force" scaling
+            MOMENT      //!< The "moment" scaling
+        };
 
-    //- Names for the output format types
-    static const Enum<outputFormatType> formatNames;
 
-    //- Names for the scaling types
+    // Static Data
-    static const Enum<scalingType> scalingNames;
+
+        //- Names for the output format types
+        static const Enum<outputFormatType> formatNames;
+
+        //- Names for the scaling types
+        static const Enum<scalingType> scalingNames;
 
 
 private:
 
-    // Private data
+    // Private Class
 
-        //- Reference axis for the locations
+        //- The controller names and faceCentre mapping for a patch
-        vector axis_;
+        struct patchControl
+        {
+            wordList names_;
+            labelList faceToPoint_;
+            List<lumpedPointInterpolator> interp_;
+        };
 
-        //- The locations of the lumped points within the reference axis
-        //  The interpolator needs scalarField, not scalarList.
-        scalarField locations_;
 
-        //- The division when calculating pressure forces (0-1)
+    // Private Data
-        scalar division_;
+
+        //- The offset for lumped points, used on input.
+        point origin_;
+
+        //- The initial state of positions with zero rotations.
+        lumpedPointState state0_;
+
+        //- The current state (positions, rotations)
+        //  Eg, as response from external application
+        lumpedPointState state_;
+
+        //- The original point ids (optional information)
+        labelList originalIds_;
+
+        //- Connectivity for the controllers
+        HashPtrTable<lumpedPointController> controllers_;
+
+        //- The controller names and faceCentre mapping for patches
+        Map<patchControl> patchControls_;
 
         //- The relaxation factor when moving points (default: 1)
         scalar relax_;
 
-        //- The interpolation type (linear|spline)
+        //- Optional owner information (patch owner)
-        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
         externalFileCoupler coupler_;
 
-        //- File io
+        //- File IO names
         word inputName_;
         word outputName_;
         word logName_;
 
+        //- The input format for points, rotations
         lumpedPointState::inputFormatType inputFormat_;
+
+        //- The output format for forces, moments
         outputFormatType outputFormat_;
 
-        //- Optional scale factors for input/output files
+        //- Scale factors for input/output files (optional)
         FixedList<scalar, 1> scaleInput_;
         FixedList<scalar, 3> scaleOutput_;
 
+        //- Calculation frequency
+        label calcFrequency_;
 
-    // Demand-driven private data
+        //- The last timeIndex when coupling was triggered
-
+        mutable label lastTrigger_;
-        //- 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;
-
-
         //- No copy construct
         lumpedPointMovement(const lumpedPointMovement&) = delete;
 
         //- No copy assignment
         void operator=(const lumpedPointMovement&) = delete;
 
+
 public:
 
-    // Static data members
+    // Static Data Members
+
+        //- Debug switch
+        static int debug;
 
         //- The canonical name ("lumpedPointMovement") for the dictionary
         static const word canonicalName;
@@ -246,15 +243,15 @@ public:
 
     // Constructors
 
-        //- Construct null
+        //- Default construct
         lumpedPointMovement();
 
         //- Construct from dictionary, optionally with some owner information
-        lumpedPointMovement(const dictionary& dict, label ownerId=-1);
+        explicit lumpedPointMovement(const dictionary& dict, label ownerId=-1);
 
 
     //- Destructor
-    virtual ~lumpedPointMovement();
+    virtual ~lumpedPointMovement() = default;
 
 
     // Member Functions
@@ -269,41 +266,36 @@ public:
         //- The number of lumped points (number of locations)
         inline label size() const;
 
-        //- The normalized reference axis
-        inline const vector& axis() const;
-
-        //- Read access to the locations
-        inline const scalarField& locations() const;
-
-        //- The division (0-1) when calculating pressure forces
-        inline scalar division() const;
-
         //- An owner Id, if needed for bookkeeping purposes
         inline label ownerId() const;
 
         //- Change the owner id, if needed for bookkeeping purposes
         inline void ownerId(label id);
 
-
-        //- Threshold locations for pressure forces
-        inline const scalarField& thresholds() const;
-
-        //- Classify the position to be located in one of the threshold zones
-        inline label threshold(const point& position) const;
-
-
         //- Communication control
         inline const externalFileCoupler& coupler() const;
 
         //- Communication control
         inline externalFileCoupler& coupler();
 
+        //- Check if coupling is pending (according to the calcFrequency)
+        bool couplingPending(const label timeIndex) const;
+
+        //- Register that coupling is completed at this calcFrequency
+        void couplingCompleted(const label timeIndex) const;
+
         //- The initial state (positions/rotations)
         inline const lumpedPointState& state0() const;
 
         //- The current state (positions/rotations)
         inline const lumpedPointState& state() const;
 
+        //- The offset for lumped points, used on input.
+        inline const point& origin() const;
+
+        //- Scale the lumped points (on input).
+        inline void scalePoints(lumpedPointState& state) const;
+
         //- The relaxation factor when changing states
         inline scalar relax() const;
 
@@ -325,25 +317,40 @@ public:
         //- The output (forces) file format
         inline lumpedPointMovement::outputFormatType outputFormat() const;
 
+        //- The Euler-angle rotation order
+        inline quaternion::eulerOrder rotationOrder() const;
 
-        //- Define the bounding-box required to enclose the specified patches.
+        //-  Rotation angles in degrees
-        //  Calculates the centre as required.
+        inline bool degrees() const;
+
+        //- Check if patch control exists for specified patch
+        inline bool hasPatchControl(const label patchIndex) const;
+
+        //- Check if patch control exists for specified patch
+        bool hasPatchControl(const polyPatch& pp) const;
+
+        //- Check if patch control exists for specified patch
+        bool hasInterpolator(const pointPatch& fpatch) const;
+
+        //- Check if patch control exists for specified patch
+        void checkPatchControl(const polyPatch& pp) const;
+
+        //- Define pressure-zones mapping for faces in the specified patches.
+        //  The face centres are compared to the controller points,
         //
-        //  \param mesh     The volume mesh reference
+        //  \param pp       The patch with a control
-        //  \param patchIds The patch ids to be included in the bounding box.
+        //  \param ctrlNames The patch ids to be included in the mapping
         //  \param points0  The initial mesh points, prior to movement
-        void setBoundBox
+        void setPatchControl
         (
-            const polyMesh& mesh,
+            const polyPatch& pp,
-            const labelUList& patchIds,
+            const wordList& ctrlNames,
             const pointField& points0
         );
 
 
-        //- Define the pressure-zones mapping for faces in the specified
+        //- Define pressure-zones mapping for faces in the specified patches.
-        //  patches.
+        //  The face centres are compared to the controller points,
-        //  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
@@ -356,13 +363,21 @@ public:
         );
 
 
+        //- Check if patch control exists for specified patch
+        void setInterpolator
+        (
+            const pointPatch& fpatch,
+            const pointField& points0
+        );
+
         //- True if the pressure-zones mapping has already been performed
         inline bool hasMapping() const;
 
-        //- Return the pressure-zones mapping with the associated
+        //- Check if patch interpolator exists for specified patch
-        //  patch face ids.
+        inline bool hasInterpolator(const label patchIndex) const;
-        inline const List<labelList>& zones() const;
 
+        //- The areas for each pressure-zone.
+        List<scalar> areas(const polyMesh& pmesh) const;
 
         //- The forces and moments acting on each pressure-zone.
         //  The zones must be previously defined via setMapping.
@@ -375,24 +390,28 @@ public:
 
 
         //- Displace points according to the current state
-        tmp<pointField> displacePoints
+        tmp<pointField> pointsDisplacement
         (
-            const pointField& points0,
+            const pointPatch& fpatch,
-            const labelList& pointLabels
+            const pointField& points0
         ) const;
 
-
         //- Displace points according to specified state
-        tmp<pointField> displacePoints
+        tmp<pointField> pointsDisplacement
         (
             const lumpedPointState& state,
-            const pointField& points0,
+            const pointPatch& fpatch,
-            const labelList& pointLabels
+            const pointField& points0
         ) const;
 
+        //- The points absolute position according to specified state
+        tmp<pointField> pointsPosition
+        (
+            const lumpedPointState& state,
+            const pointPatch& fpatch,
+            const pointField& points0
+        ) const;
 
-        //- Interpolation weights
-        const interpolationWeights& interpolator() const;
 
         //- Write axis, locations, division as a dictionary
         void writeDict(Ostream& os) const;
@@ -418,6 +437,13 @@ public:
         //- Read state from file, applying relaxation as requested
         bool readState();
 
+        //- Write state as VTK PolyData format.
+        void writeStateVTP
+        (
+            const lumpedPointState& state,
+            const fileName& file
+        ) const;
+
         //- Write state as VTK PolyData format.
         void writeStateVTP(const fileName& file) const;
 
@@ -438,14 +464,12 @@ public:
             const pointField& points0
         ) const;
 
-
         //- Write displaced geometry according to the current state,
         //  write as VTK PolyData format.
         void writeVTP
         (
             const fileName& file,
             const polyMesh& mesh,
-            const labelUList& patchIds,
             const pointField& points0
         ) const;
 
@@ -456,25 +480,8 @@ public:
             const fileName& file,
             const lumpedPointState& state,
             const polyMesh& mesh,
-            const labelUList& patchLst,
             const pointField& points0
         ) const;
-
-
-    // Housekeeping
-
-        //- Deprecated(2018-08) compatibility method
-        //  \deprecated(2018-08) - use IOobject::selectIO directly
-        static IOobject selectIO
-        (
-            const IOobject& io,
-            const fileName& altFile,
-            const word& ioName = ""
-        )
-        {
-            return IOobject::selectIO(io, altFile, ioName);
-        }
-
 };
 
 

src/lumpedPointMotion/movement/lumpedPointMovementI.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,31 +27,13 @@ License
 
 inline bool Foam::lumpedPointMovement::empty() const
 {
-    return locations_.empty();
+    return state0_.empty();
 }
 
 
 inline Foam::label Foam::lumpedPointMovement::size() const
 {
-    return locations_.size();
+    return state0_.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_;
 }
 
 
@@ -67,21 +49,19 @@ inline void Foam::lumpedPointMovement::ownerId(label id)
 }
 
 
-inline const Foam::scalarField& Foam::lumpedPointMovement::thresholds() const
+inline bool
+Foam::lumpedPointMovement::hasPatchControl(const label patchIndex) const
 {
-    if (!thresholdPtr_)
+    return patchControls_.found(patchIndex);
-    {
-        calcThresholds();
-    }
-
-    return *thresholdPtr_;
 }
 
 
-inline Foam::label
+inline bool
-Foam::lumpedPointMovement::threshold(const point& position) const
+Foam::lumpedPointMovement::hasInterpolator(const label patchIndex) const
 {
-    return threshold(position & axis_);
+    const auto iter = patchControls_.cfind(patchIndex);
+
+    return (iter.good() && iter().interp_.size());
 }
 
 
@@ -98,7 +78,6 @@ inline Foam::externalFileCoupler& Foam::lumpedPointMovement::coupler()
 }
 
 
-//- The initial state (positions/rotations)
 inline const Foam::lumpedPointState& Foam::lumpedPointMovement::state0() const
 {
     return state0_;
@@ -111,6 +90,21 @@ inline const Foam::lumpedPointState& Foam::lumpedPointMovement::state() const
 }
 
 
+inline const Foam::point& Foam::lumpedPointMovement::origin() const
+{
+    return origin_;
+}
+
+
+inline void Foam::lumpedPointMovement::scalePoints
+(
+    lumpedPointState& state
+) const
+{
+    state.scalePoints(scaleInput_[scalingType::LENGTH]);
+}
+
+
 inline Foam::scalar Foam::lumpedPointMovement::relax() const
 {
     return relax_;
@@ -155,16 +149,22 @@ Foam::lumpedPointMovement::outputFormat() const
 }
 
 
-inline bool Foam::lumpedPointMovement::hasMapping() const
+inline Foam::quaternion::eulerOrder
+Foam::lumpedPointMovement::rotationOrder() const
 {
-    return !faceZones_.empty();
+    return state0().rotationOrder();
 }
 
 
-inline const Foam::List<Foam::labelList>&
+inline bool Foam::lumpedPointMovement::degrees() const
-Foam::lumpedPointMovement::zones() const
 {
-    return faceZones_;
+    return state0().degrees();
+}
+
+
+inline bool Foam::lumpedPointMovement::hasMapping() const
+{
+    return !patchControls_.empty();
 }
 
 

src/lumpedPointMotion/movement/lumpedPointMovementWriter.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,16 +27,54 @@ License
 
 #include "lumpedPointMovement.H"
 #include "polyMesh.H"
+#include "pointMesh.H"
 #include "OFstream.H"
-#include "uindirectPrimitivePatch.H"
-
 #include "foamVtkOutput.H"
+#include "foamVtkSurfaceWriter.H"
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+void Foam::lumpedPointMovement::writeStateVTP
+(
+    const lumpedPointState& state,
+    const fileName& file
+) const
+{
+    if (!Pstream::master())
+    {
+        // No extra information available from slaves, write on master only.
+        return;
+    }
+
+    labelListList lines;
+
+    label nLines = controllers_.size();
+
+    if (nLines)
+    {
+        lines.resize(nLines);
+        nLines = 0;
+
+        for (const word& ctrlName : controllers_.sortedToc())
+        {
+            lines[nLines] = controllers_[ctrlName]->pointLabels();
+            ++nLines;
+        }
+    }
+    else
+    {
+        // Default - global with all points as a single line
+        lines.resize(1);
+        lines.first() = identity(state.size());
+    }
+
+    state.writeVTP(file, lines, originalIds_);
+}
+
+
 void Foam::lumpedPointMovement::writeStateVTP(const fileName& file) const
 {
-    state().writeVTP(file, axis());
+    writeStateVTP(state(), file);
 }
 
 
@@ -47,6 +85,12 @@ void Foam::lumpedPointMovement::writeForcesAndMomentsVTP
     const UList<vector>& moments
 ) const
 {
+    if (!Pstream::master())
+    {
+        // Force, moments already reduced
+        return;
+    }
+
     OFstream fos(file);
     std::ostream& os = fos.stdStream();
 
@@ -60,7 +104,7 @@ void Foam::lumpedPointMovement::writeForcesAndMomentsVTP
         .beginVTKFile<vtk::fileTag::POLY_DATA>();
 
     //
-    // The 'spine' of lumped mass points
+    // The 'backbone' of lumped mass points
     //
     const label nPoints = state().points().size();
 
@@ -176,126 +220,38 @@ void Foam::lumpedPointMovement::writeZonesVTP
     const pointField& points0
 ) const
 {
-    OFstream fos(file);
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    std::ostream& os = fos.stdStream();
+    const labelList patchIds(patchControls_.sortedToc());
 
-    autoPtr<vtk::formatter> format = vtk::newFormatter
+    vtk::surfaceWriter writer
     (
-        os,
+        pointField::null(),
-        vtk::formatType::INLINE_ASCII
+        faceList::null(),
+        vtk::formatType::INLINE_ASCII,
+        file
     );
 
-    format().xmlHeader()
+    for (const label patchi : patchIds)
-        .beginVTKFile<vtk::fileTag::POLY_DATA>();
-
-    forAll(faceZones_, zoneI)
     {
-        uindirectPrimitivePatch pp
+        const labelList& faceToPoint = patchControls_[patchi].faceToPoint_;
+
+        primitivePatch pp
         (
-            UIndirectList<face>(mesh.faces(), faceZones_[zoneI]),
+            SubList<face>(mesh.faces(), patches[patchi].range()),
             points0
         );
 
-        format()
+        writer.piece(pp.localPoints(), pp.localFaces());
-            .tag
-            (
-                vtk::fileTag::PIECE,
-                vtk::fileAttr::NUMBER_OF_POINTS, pp.nPoints(),
-                vtk::fileAttr::NUMBER_OF_POLYS,  pp.size()
-            );
 
-        // 'points'
+        writer.writeGeometry();
-        {
-            const uint64_t payLoad = vtk::sizeofData<float, 3>(pp.nPoints());
 
-            format()
+        writer.beginCellData(2);
-                .tag(vtk::fileTag::POINTS)
-                .beginDataArray<float, 3>(vtk::dataArrayAttr::POINTS);
 
-            format().writeSize(payLoad);
+        writer.writeUniform("patchId", patchi);
-            vtk::writeList(format(), pp.localPoints());
+        writer.write("lumpedId", faceToPoint);
-            format().flush();
 
-            format()
+        writer.endCellData();
-                .endDataArray()
-                .endTag(vtk::fileTag::POINTS);
-        }
-
-        // <Polys>
-        format().tag(vtk::fileTag::POLYS);
-
-        //
-        // 'connectivity'
-        //
-        {
-            label nVerts = 0;
-            for (const face& f : pp)
-            {
-                nVerts += f.size();
-            }
-
-            const uint64_t payLoad = vtk::sizeofData<label>(nVerts);
-
-            format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
-            format().writeSize(payLoad);
-
-            for (const face& f : pp.localFaces())
-            {
-                vtk::writeList(format(), f);
-            }
-            format().flush();
-
-            format().endDataArray();
-        }
-
-        //
-        // 'offsets'  (connectivity offsets)
-        //
-        {
-            const uint64_t payLoad = vtk::sizeofData<label>(pp.size());
-
-            format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
-            format().writeSize(payLoad);
-
-            label off = 0;
-            forAll(pp, facei)
-            {
-                const face& f = pp[facei];
-                off += f.size();
-
-                format().write(off);
-            }
-            format().flush();
-
-            format().endDataArray();
-        }
-
-        format().endTag(vtk::fileTag::POLYS);
-
-
-        format().beginCellData();
-
-        // zone Id
-        {
-            const uint64_t payLoad = vtk::sizeofData<label>(pp.size());
-
-            format().beginDataArray<label>("zoneId");
-            format().writeSize(payLoad);
-
-            vtk::write(format(), zoneI, pp.size());
-
-            format().flush();
-
-            format().endDataArray();
-        }
-
-        format().endCellData();
-
-        format().endPiece();
     }
-
-    format().endTag(vtk::fileTag::POLY_DATA)
-        .endVTKFile();
 }
 
 
@@ -303,11 +259,10 @@ void Foam::lumpedPointMovement::writeVTP
 (
     const fileName& file,
     const polyMesh& mesh,
-    const labelUList& patchIds,
     const pointField& points0
 ) const
 {
-    writeVTP(file, state(), mesh, patchIds, points0);
+    writeVTP(file, state(), mesh, points0);
 }
 
 
@@ -316,119 +271,100 @@ void Foam::lumpedPointMovement::writeVTP
     const fileName& file,
     const lumpedPointState& state,
     const polyMesh& mesh,
-    const labelUList& patchIds,
     const pointField& points0
 ) const
 {
-    const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+    const labelList patchIds(patchControls_.sortedToc());
 
-    OFstream fos(file);
+    pointMesh ptMesh(mesh);
-    std::ostream& os = fos.stdStream();
 
-    autoPtr<vtk::formatter> format = vtk::newFormatter
+    vtk::surfaceWriter writer
     (
-        os,
+        pointField::null(),
-        vtk::formatType::INLINE_ASCII
+        faceList::null(),
+        vtk::formatType::INLINE_ASCII,
+        file
     );
 
-    format().xmlHeader()
+    for (const label patchi : patchIds)
-        .beginVTKFile<vtk::fileTag::POLY_DATA>();
-
-    for (const label patchId : patchIds)
     {
-        const polyPatch& pp = boundaryMesh[patchId];
+        const polyPatch& pp = patches[patchi];
 
-        format()
+        const pointPatch& ptPatch = ptMesh.boundary()[patchi];
-            .tag
-            (
-                vtk::fileTag::PIECE,
-                vtk::fileAttr::NUMBER_OF_POINTS, pp.nPoints(),
-                vtk::fileAttr::NUMBER_OF_POLYS,  pp.size()
-            );
 
-        // 'points'
+        // Current position (not displacement)
+        tmp<pointField> tpts = pointsPosition(state, ptPatch, points0);
+
+        writer.piece(tpts(), pp.localFaces());
+
+        writer.writeGeometry();
+
+        // Face mapping
+        const labelList& faceToPoint = patchControls_[patchi].faceToPoint_;
+
+        writer.beginCellData(2);
+
+        writer.writeUniform("patchId", patchi);
+        writer.write("lumpedId", faceToPoint);
+
+        writer.endCellData();
+
+        // The interpolator
+        const List<lumpedPointInterpolator>& interpList
+            = patchControls_[patchi].interp_;
+
+        writer.beginPointData(3);
+
+        // Nearest, Next
         {
-            const uint64_t payLoad = vtk::sizeofData<float, 3>(pp.nPoints());
+            labelList intData(interpList.size());
 
-            format()
+            forAll(interpList, i)
-                .tag(vtk::fileTag::POINTS)
+            {
-                .beginDataArray<float, 3>(vtk::dataArrayAttr::POINTS);
+                intData[i] = interpList[i].nearest();
+            }
+            writer.write("nearest", intData);
 
-            // Could be more efficient, but not often needed
+            forAll(interpList, i)
-            tmp<pointField> tpts = displacePoints
+            {
-            (
+                intData[i] = interpList[i].next1();
-                state,
+            }
-                points0,
+            writer.write("next1", intData);
-                pp.meshPoints()
-            ) + pointField(points0, pp.meshPoints());
 
-            const pointField& pts = tpts();
 
-            format().writeSize(payLoad);
+            forAll(interpList, i)
-            vtk::writeList(format(), pts);
+            {
-            format().flush();
+                intData[i] = interpList[i].next2();
-
+            }
-            format()
+            writer.write("next2", intData);
-                .endDataArray()
-                .endTag(vtk::fileTag::POINTS);
         }
 
-        // <Polys>
+        // Weights
-        format().tag(vtk::fileTag::POLYS);
-
-
-        //
-        // 'connectivity'
-        //
         {
-            label nVerts = 0;
+            scalarList floatData(interpList.size());
-            for (const face& f : pp)
+
+            forAll(interpList, i)
             {
-                nVerts += f.size();
+                floatData[i] = interpList[i].weight0();
             }
+            writer.write("weight", floatData);
 
-            const uint64_t payLoad = vtk::sizeofData<label>(nVerts);
+            forAll(interpList, i)
-
-            format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
-            format().writeSize(payLoad);
-
-            for (const face& f : pp.localFaces())
             {
-                vtk::writeList(format(), f);
+                floatData[i] = interpList[i].weight1();
             }
-            format().flush();
+            writer.write("weight1", floatData);
 
-            format().endDataArray();
+            forAll(interpList, i)
+            {
+                floatData[i] = interpList[i].weight2();
+            }
+            writer.write("weight2", floatData);
         }
 
-        //
+        writer.endPointData();
-        // 'offsets'  (connectivity offsets)
-        //
-        {
-            const uint64_t payLoad = vtk::sizeofData<label>(pp.size());
-
-            format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
-            format().writeSize(payLoad);
-
-            label off = 0;
-            for (const face& f : pp)
-            {
-                off += f.size();
-
-                format().write(off);
-            }
-            format().flush();
-
-            format().endDataArray();
-        }
-
-        format().endTag(vtk::fileTag::POLYS);
-
-        format().endPiece();
     }
-
-    format().endTag(vtk::fileTag::POLY_DATA)
-        .endVTKFile();
 }
 
 

src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,6 +50,74 @@ namespace Foam
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
+Foam::label Foam::lumpedPointDisplacementPointPatchVectorField::setPatchControls
+(
+    const pointVectorField& pvf,
+    const pointField& points0
+)
+{
+    label count = 0;
+
+    const pointVectorField::Boundary& bf = pvf.boundaryField();
+    const polyBoundaryMesh& patches = pvf.mesh().mesh().boundaryMesh();
+
+    forAll(bf, patchi)
+    {
+        // Patch of this type
+        const auto* p = isA<patchType>(bf[patchi]);
+
+        if (p)
+        {
+            // Patch controls (mapping) for calculating forces/moments
+            const_cast<lumpedPointMovement&>(p->movement())
+                .setPatchControl
+                (
+                    patches[patchi],
+                    p->controllers(),
+                    points0
+                );
+
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+Foam::label Foam::lumpedPointDisplacementPointPatchVectorField::setInterpolators
+(
+    const pointVectorField& pvf,
+    const pointField& points0
+)
+{
+    label count = 0;
+
+    const pointVectorField::Boundary& bf = pvf.boundaryField();
+
+    forAll(bf, patchi)
+    {
+        // Patch of this type
+        const auto* p = isA<patchType>(bf[patchi]);
+
+        if (p)
+        {
+            // Patch controls (mapping) for calculating forces/moments
+            const_cast<lumpedPointMovement&>(p->movement())
+                .setInterpolator
+                (
+                    p->patch(),
+                    points0
+                );
+
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
 Foam::labelList
 Foam::lumpedPointDisplacementPointPatchVectorField::patchIds
 (
@@ -61,7 +129,7 @@ Foam::lumpedPointDisplacementPointPatchVectorField::patchIds
     DynamicList<label> patchLst(bf.size());
     forAll(bf, patchi)
     {
-        // All patches of this type
+        // Patch of this type
         if (isA<patchType>(bf[patchi]))
         {
             patchLst.append(patchi);
@@ -80,11 +148,33 @@ Foam::lumpedPointDisplacementPointPatchVectorField::points0() const
 {
     const objectRegistry& obr = this->patch().boundaryMesh().mesh().db();
 
-    // Obtain starting locations from the motionSolver
+    // Obtain starting locations from the motionSolver (when possible)
-    return obr.lookupObject<displacementMotionSolver>
+    const auto* solver =
-    (
+        obr.cfindObject<displacementMotionSolver>("dynamicMeshDict");
-        "dynamicMeshDict"
+
-    ).points0();
+    if (solver)
+    {
+        if (points0Ptr_)
+        {
+            points0Ptr_.reset(nullptr);
+        }
+        return solver->points0();
+    }
+    else if (!points0Ptr_)
+    {
+        points0Ptr_.reset
+        (
+            new pointIOField
+            (
+                points0MotionSolver::points0IO
+                (
+                    this->patch().boundaryMesh().mesh().mesh()
+                )
+            )
+        );
+    }
+
+    return *points0Ptr_;
 }
 
 
@@ -92,22 +182,19 @@ const Foam::lumpedPointMovement&
 Foam::lumpedPointDisplacementPointPatchVectorField::movement() const
 {
     const objectRegistry& obr = this->patch().boundaryMesh().mesh().db();
-    const lumpedPointIOMovement* ptr =
+
-        lumpedPointIOMovement::lookupInRegistry(obr);
+    lumpedPointIOMovement* ptr =
+        lumpedPointIOMovement::getMovementObject(obr);
 
     if (ptr)
     {
         return *ptr; // Already exists
     }
 
-    // create and register with this patch as the owner
+    // Create and register with this patch as the owner
-    autoPtr<lumpedPointIOMovement> obj = lumpedPointIOMovement::New
+    ptr = lumpedPointIOMovement::New(obr, this->patch().index()).ptr();
-    (
-        obr,
-        this->patch().index()
-    );
 
-    return objectRegistry::store(obj);
+    return objectRegistry::store(ptr);
 }
 
 
@@ -120,7 +207,8 @@ lumpedPointDisplacementPointPatchVectorField
     const DimensionedField<vector, pointMesh>& iF
 )
 :
-    fixedValuePointPatchField<vector>(p, iF)
+    fixedValuePointPatchField<vector>(p, iF),
+    controllers_()
 {}
 
 
@@ -132,8 +220,20 @@ lumpedPointDisplacementPointPatchVectorField
     const dictionary& dict
 )
 :
-    fixedValuePointPatchField<vector>(p, iF, dict)
+    fixedValuePointPatchField<vector>(p, iF, dict),
-{}
+    controllers_()
+{
+    dict.readIfPresent("controllers", controllers_);
+
+    if (controllers_.empty())
+    {
+        WarningInFunction
+            << "No controllers specified, using all lumped points for patch: "
+            << this->patch().name() << nl << nl;
+    }
+
+    // controllers_ : check? remove duplicates?
+}
 
 
 Foam::lumpedPointDisplacementPointPatchVectorField::
@@ -145,7 +245,8 @@ lumpedPointDisplacementPointPatchVectorField
     const pointPatchFieldMapper& mapper
 )
 :
-    fixedValuePointPatchField<vector>(pf, p, iF, mapper)
+    fixedValuePointPatchField<vector>(pf, p, iF, mapper),
+    controllers_(pf.controllers_)
 {}
 
 
@@ -156,7 +257,8 @@ lumpedPointDisplacementPointPatchVectorField
     const DimensionedField<vector, pointMesh>& iF
 )
 :
-    fixedValuePointPatchField<vector>(pf, iF)
+    fixedValuePointPatchField<vector>(pf, iF),
+    controllers_(pf.controllers_)
 {}
 
 
@@ -166,16 +268,17 @@ Foam::lumpedPointDisplacementPointPatchVectorField::
 ~lumpedPointDisplacementPointPatchVectorField()
 {
     // de-register movement if in use and managed by this patch
-    const lumpedPointIOMovement* ptr = lumpedPointIOMovement::lookupInRegistry
+    lumpedPointIOMovement* ptr =
-    (
+        lumpedPointIOMovement::getMovementObject
-        this->patch().boundaryMesh().mesh().db()
+        (
-    );
+            this->patch().boundaryMesh().mesh().db()
+        );
 
     if (ptr && ptr->ownerId() == this->patch().index())
     {
         movement().coupler().shutdown();
 
-        const_cast<lumpedPointIOMovement*>(ptr)->checkOut();
+        ptr->checkOut();
     }
 }
 
@@ -189,58 +292,68 @@ void Foam::lumpedPointDisplacementPointPatchVectorField::updateCoeffs()
         return;
     }
 
+    const label timeIndex = this->db().time().timeIndex();
+
     enum Time::stopAtControls action = Time::stopAtControls::saUnknown;
 
-    const bool masterPatch = (movement().ownerId() == this->patch().index());
+    if (movement().ownerId() == this->patch().index())
-    if (masterPatch)
     {
-        if (lumpedPointIOMovement::debug)
+        // The ownerId is always the lowest patch number,
+        // thus will always be triggered first
+
+        if (lumpedPointMovement::debug)
         {
             Pout<<"masterPatch: " << this->patch().index() << endl;
         }
 
         const polyMesh& mesh = this->patch().boundaryMesh().mesh().mesh();
 
-        // need face 'zones' for calculating forces
+        // Mapping for calculating forces
-        // likely need bounding box for the movement
+        // as well as face point interpolation
-        // -> do both now if required
         if (!movement().hasMapping())
         {
-            const_cast<lumpedPointMovement&>(movement()).setMapping
+            // Add mapping for calculating forces/moments
+            setPatchControls
             (
-                mesh,
+                static_cast<const pointVectorField&>
-                // All patches of this type
-                patchIds
                 (
-                    static_cast<const pointVectorField&>
+                    this->internalField()
-                    (
-                        this->internalField()
-                    )
                 ),
                 this->points0()
             );
         }
 
+        int triggered = 0;
 
         if
         (
-            movement().coupler().initialized()
+            movement().coupler().slaveFirst()
-         || !movement().coupler().slaveFirst()
+         && !movement().coupler().initialized()
         )
+        {
+            // Master does nothing yet, but slave will
+            triggered = -1;
+        }
+        else if (movement().couplingPending(timeIndex))
+        {
+            // Trigger is pending, or coupling not yet not initialized
+            triggered = 1;
+        }
+
+        if (triggered > 0)
         {
             // Synchronized for all processes
             List<vector> forces, moments;
             movement().forcesAndMoments(mesh, forces, moments);
 
-            if (lumpedPointIOMovement::debug)
+            if (lumpedPointMovement::debug)
             {
                 Pout<<"gatherForces: " << forces << " called from patch "
                     << this->patch().index() << endl;
 
                 if (Pstream::master())
                 {
-                    Pout<<"output forces to file: "
+                    Pout<<"output forces to file: called from patch "
-                        << movement().locations() << " called from patch "
                         << this->patch().index() << nl
                         <<"# " << forces.size() << " force entries" << nl
                         <<"# fx fy fz" << nl
@@ -252,29 +365,44 @@ void Foam::lumpedPointDisplacementPointPatchVectorField::updateCoeffs()
 
             if (Pstream::master())
             {
-                movement().writeData(forces, moments, &(db().time()));
+                movement().writeData(forces, moments, &(this->db().time()));
 
                 // Signal external source to execute
                 movement().coupler().useSlave();
             }
         }
 
-        // Wait for slave to provide data (includes MPI barrier)
+        if (triggered)
-        // and catch any abort information sent from slave
+        {
-        action = movement().coupler().waitForSlave();
+            // Wait for slave to provide data (includes MPI barrier)
+            // and catch any abort information sent from slave
+            action = movement().coupler().waitForSlave();
 
-        // Read data passed back from external source - includes MPI barrier
+            const_cast<lumpedPointMovement&>(movement()).readState();
-        const_cast<lumpedPointMovement&>(movement()).readState();
+
+            movement().couplingCompleted(timeIndex);
+        }
     }
 
-    tmp<pointField> tdisp = movement().displacePoints
+    if (!movement().hasInterpolator(this->patch()))
-    (
+    {
-        this->points0(),
+        const_cast<lumpedPointMovement&>(movement()).setInterpolator
-        this->patch().meshPoints()
+        (
-    );
+            this->patch(),
+            this->points0()
+        );
+    }
+
+    tmp<pointField> tdisp =
+        movement().pointsDisplacement
+        (
+            this->patch(),
+            this->points0()
+        );
 
     this->operator==(tdisp);
 
+
     fixedValuePointPatchField<vector>::updateCoeffs();
 
     // Process any abort information sent from slave
@@ -293,6 +421,12 @@ void Foam::lumpedPointDisplacementPointPatchVectorField::write(Ostream& os)
 const
 {
     pointPatchField<vector>::write(os);
+
+    if (controllers_.size())
+    {
+        os.writeEntry("controllers", controllers_);
+    }
+
     writeEntry("value", os);
 }
 

src/lumpedPointMotion/pointPatchFields/lumpedPointDisplacementPointPatchVectorField.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,6 +42,7 @@ Description
     {
         type    lumpedPointDisplacement;
         value   uniform (0 0 0);
+        controllers ( controllerName1 controllerName2 );
     }
     \endverbatim
 
@@ -67,7 +68,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declarations
+// Forward Declarations
 class interpolationWeights;
 
 /*---------------------------------------------------------------------------*\
@@ -78,7 +79,14 @@ class lumpedPointDisplacementPointPatchVectorField
 :
     public fixedValuePointPatchField<vector>
 {
-    // Private data
+    // Private Data
+
+        //- Names of the movement controller(s) in use
+        wordList controllers_;
+
+        //- Backup method for getting "points0" without a motion solver
+        mutable autoPtr<pointIOField> points0Ptr_;
+
 
     //- Convenience typedefs
     typedef lumpedPointDisplacementPointPatchVectorField patchType;
@@ -167,7 +175,27 @@ public:
     virtual ~lumpedPointDisplacementPointPatchVectorField();
 
 
-    // Member functions
+    // Member Functions
+
+        //- The controller names for this patch
+        const wordList& controllers() const
+        {
+            return controllers_;
+        }
+
+        //- Set all patch controls for patches of this type
+        static label setPatchControls
+        (
+            const pointVectorField& pvf,
+            const pointField& points0
+        );
+
+        //- Set all patch controls for patches of this type
+        static label setInterpolators
+        (
+            const pointVectorField& pvf,
+            const pointField& points0
+        );
 
         //- The ids for all patches of this type
         static labelList patchIds(const pointVectorField& pvf);

src/lumpedPointMotion/state/lumpedPointState.C

@@ -26,7 +26,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "lumpedPointState.H"
-#include "demandDrivenData.H"
 #include "unitConversion.H"
 #include "EulerCoordinateRotation.H"
 #include "IFstream.H"
@@ -44,6 +43,11 @@ Foam::lumpedPointState::formatNames
 });
 
 
+Foam::scalar Foam::lumpedPointState::visLength = 0.1;
+
+bool Foam::lumpedPointState::visUnused = true;
+
+
 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
 //! \cond fileScope
@@ -69,7 +73,7 @@ static Foam::string getLineNoComment
 
 void Foam::lumpedPointState::calcRotations() const
 {
-    rotationPtr_ = new tensorField(angles_.size());
+    rotationPtr_.reset(new tensorField(angles_.size()));
 
     auto rotIter = rotationPtr_->begin();
 
@@ -83,20 +87,26 @@ void Foam::lumpedPointState::calcRotations() const
 }
 
 
-void Foam::lumpedPointState::readDict(const dictionary& dict)
+void Foam::lumpedPointState::readDict
+(
+    const dictionary& dict,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
 {
     dict.readEntry("points", points_);
     dict.readEntry("angles", angles_);
     order_ =
         quaternion::eulerOrderNames.getOrDefault
         (
-            "order",
+            "rotationOrder",
             dict,
-            quaternion::eulerOrder::ZXZ
+            rotOrder
         );
-    degrees_ = dict.getOrDefault("degrees", false);
 
-    deleteDemandDrivenData(rotationPtr_);
+    degrees_ = dict.getOrDefault("degrees", degrees);
+
+    rotationPtr_.reset(nullptr);
 }
 
 
@@ -122,43 +132,82 @@ Foam::lumpedPointState::lumpedPointState(const lumpedPointState& rhs)
 {}
 
 
-Foam::lumpedPointState::lumpedPointState(const pointField& pts)
+Foam::lumpedPointState::lumpedPointState
+(
+    const pointField& pts,
+    const vectorField& ang,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
 :
     points_(pts),
-    angles_(points_.size(), Zero),
+    angles_(ang),
-    order_(quaternion::eulerOrder::ZXZ),
+    order_(rotOrder),
-    degrees_(false),
+    degrees_(degrees),
-    rotationPtr_(nullptr)
-{}
-
-
-Foam::lumpedPointState::lumpedPointState(tmp<pointField>& pts)
-:
-    points_(pts),
-    angles_(points_.size(), Zero),
-    order_(quaternion::eulerOrder::ZXZ),
-    degrees_(false),
-    rotationPtr_(nullptr)
-{}
-
-
-Foam::lumpedPointState::lumpedPointState(const dictionary& dict)
-:
-    points_(),
-    angles_(),
-    order_(quaternion::eulerOrder::ZXZ),
-    degrees_(false),
     rotationPtr_(nullptr)
 {
-    readDict(dict);
+    if (points_.size() != angles_.size())
+    {
+        #ifdef FULLDEBUG
+        FatalErrorInFunction
+            << "Have " << points_.size() << " points but "
+            << angles_.size() << " angles" << nl
+            exit(FatalError);
+        #else
+        WarningInFunction
+            << "Have " << points_.size() << " points but "
+            << angles_.size() << " angles, resizing angles to match" << nl;
+        #endif
+
+        angles_.resize(points_.size(), Zero);
+    }
 }
 
 
-// * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
+Foam::lumpedPointState::lumpedPointState
+(
+    const pointField& pts,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
+:
+    points_(pts),
+    angles_(points_.size(), Zero),
+    order_(rotOrder),
+    degrees_(degrees),
+    rotationPtr_(nullptr)
+{}
 
-Foam::lumpedPointState::~lumpedPointState()
+
+Foam::lumpedPointState::lumpedPointState
+(
+    tmp<pointField>& pts,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
+:
+    points_(pts),
+    angles_(points_.size(), Zero),
+    order_(rotOrder),
+    degrees_(false),
+    rotationPtr_(nullptr)
+{}
+
+
+Foam::lumpedPointState::lumpedPointState
+(
+    const dictionary& dict,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
+:
+    points_(),
+    angles_(),
+    order_(rotOrder),
+    degrees_(degrees),
+    rotationPtr_(nullptr)
 {
-    deleteDemandDrivenData(rotationPtr_);
+    readDict(dict, rotOrder, degrees);
 }
 
 
@@ -171,7 +220,16 @@ void Foam::lumpedPointState::operator=(const lumpedPointState& rhs)
     order_   = rhs.order_;
     degrees_ = rhs.degrees_;
 
-    deleteDemandDrivenData(rotationPtr_);
+    rotationPtr_.reset(nullptr);
+}
+
+
+void Foam::lumpedPointState::operator+=(const point& origin)
+{
+    for (point& p : points_)
+    {
+        p += origin;
+    }
 }
 
 
@@ -209,11 +267,16 @@ void Foam::lumpedPointState::relax
 
     angles_ = convert*prev.angles_ + alpha*(angles_ - convert*prev.angles_);
 
-    deleteDemandDrivenData(rotationPtr_);
+    rotationPtr_.reset(nullptr);
 }
 
 
-bool Foam::lumpedPointState::readPlain(Istream& is)
+bool Foam::lumpedPointState::readPlain
+(
+    Istream& is,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
 {
     // Assume generic input stream so we can do line-based input
     ISstream& iss = dynamic_cast<ISstream&>(is);
@@ -226,8 +289,8 @@ bool Foam::lumpedPointState::readPlain(Istream& is)
         isstr >> count;
     }
 
-    points_.setSize(count);
+    points_.resize(count);
-    angles_.setSize(count);
+    angles_.resize(count);
 
     count = 0;
     forAll(points_, i)
@@ -242,21 +305,26 @@ bool Foam::lumpedPointState::readPlain(Istream& is)
         ++count;
     }
 
-    points_.setSize(count);
+    points_.resize(count);
-    angles_.setSize(count);
+    angles_.resize(count);
     order_ = quaternion::eulerOrder::ZXZ;
     degrees_ = false;
 
-    deleteDemandDrivenData(rotationPtr_);
+    rotationPtr_.reset(nullptr);
 
     return count;
 }
 
 
-bool Foam::lumpedPointState::readData(Istream& is)
+bool Foam::lumpedPointState::readData
+(
+    Istream& is,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
+)
 {
     dictionary dict(is);
-    readDict(dict);
+    readDict(dict, rotOrder, degrees);
 
     return points_.size();
 }
@@ -292,9 +360,9 @@ void Foam::lumpedPointState::writePlain(Ostream& os) const
 
     forAll(points_, i)
     {
-        const vector& pt = points_[i];
+        const vector& p = points_[i];
 
-        os  << pt.x() << ' ' << pt.y() << ' ' << pt.z();
+        os  << p.x() << ' ' << p.y() << ' ' << p.z();
 
         if (i < angles_.size())
         {
@@ -312,8 +380,10 @@ void Foam::lumpedPointState::writePlain(Ostream& os) const
 
 bool Foam::lumpedPointState::readData
 (
-    const inputFormatType& fmt,
+    const inputFormatType fmt,
-    const fileName& file
+    const fileName& file,
+    const quaternion::eulerOrder rotOrder,
+    const bool degrees
 )
 {
     bool ok = false;
@@ -323,11 +393,11 @@ bool Foam::lumpedPointState::readData
 
         if (fmt == inputFormatType::PLAIN)
         {
-            ok = this->readPlain(is);
+            ok = this->readPlain(is, rotOrder, degrees);
         }
         else
         {
-            ok = this->readData(is);
+            ok = this->readData(is, rotOrder, degrees);
         }
     }
 
@@ -375,7 +445,7 @@ bool Foam::lumpedPointState::readData
             toBelow << points_ << angles_ << degrees_;
         }
 
-        deleteDemandDrivenData(rotationPtr_);
+        rotationPtr_.reset(nullptr);
 
         // MPI barrier
         Pstream::scatter(ok);

src/lumpedPointMotion/state/lumpedPointState.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,7 +37,7 @@ Description
         Property    | Description                           | Required | Default
         points      | List of points                        | yes |
         angles      | List of Euler rotation angles         | yes |
-        order       | The Euler-angle rotation order        | no  | zxz
+        rotationOrder | The Euler-angle rotation order      | no  | zxz
         degrees     | Rotation angles in degrees            | no  | false
     \endtable
 
@@ -77,6 +77,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class Istream;
 class Ostream;
 
@@ -88,17 +89,20 @@ class lumpedPointState
 {
 public:
 
-    //- Input format types
+    // Data Types
-    enum class inputFormatType
-    {
-        PLAIN,      //!< "plain" is a simple ASCII format
-        DICTIONARY  //!< "dictionary" is the OpenFOAM dictionary format
-    };
 
-    // Static data
+        //- Input format types
+        enum class inputFormatType
+        {
+            PLAIN,      //!< "plain" is a simple ASCII format
+            DICTIONARY  //!< "dictionary" is the OpenFOAM dictionary format
+        };
 
-    //- Names for the input format types
+
-    static const Enum<inputFormatType> formatNames;
+    // Static Data
+
+        //- Names for the input format types
+        static const Enum<inputFormatType> formatNames;
 
 
 private:
@@ -118,37 +122,75 @@ private:
         bool degrees_;
 
         //- Tensor rotation of lumped points
-        mutable tensorField* rotationPtr_;
+        mutable unique_ptr<tensorField> rotationPtr_;
 
 
     // Private Member Functions
 
         void calcRotations() const;
 
-        void readDict(const dictionary& dict);
+        void readDict
+        (
+            const dictionary& dict,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
 public:
 
+    // Public Data
+
+        //- Enable/disable visualization of unused points
+        static bool visUnused;
+
+        //- The length for visualization triangles
+        static scalar visLength;
+
+
     // Constructors
 
-        //- Construct null
+        //- Default construct
         lumpedPointState();
 
-        //- Copy constructor
+        //- Copy construct
         lumpedPointState(const lumpedPointState& rhs);
 
-        //- Construct from points with zero-rotation
+        //- Copy construct from points and angles
-        lumpedPointState(const pointField& pts);
+        lumpedPointState
+        (
+            const pointField& pts,
+            const vectorField& ang,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
+
+        //- Copy construct from points with zero-rotation
+        explicit lumpedPointState
+        (
+            const pointField& pts,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
         //- Construct from points with zero-rotation
-        lumpedPointState(tmp<pointField>& pts);
+        explicit lumpedPointState
+        (
+            tmp<pointField>& pts,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
         //- Construct from dictionary
-        lumpedPointState(const dictionary& dict);
+        explicit lumpedPointState
+        (
+            const dictionary& dict,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
 
     //- Destructor
-    virtual ~lumpedPointState();
+    virtual ~lumpedPointState() = default;
 
 
     // Member Functions
@@ -175,6 +217,12 @@ public:
         //  Zero and negative values are ignored.
         void scalePoints(const scalar scaleFactor);
 
+        //- The Euler-angle rotation order
+        inline quaternion::eulerOrder rotationOrder() const;
+
+        //-  Rotation angles in degrees
+        inline bool degrees() const;
+
         //- Relax the state
         //  alpha = 1 : no relaxation
         //  alpha < 1 : relaxation
@@ -182,7 +230,12 @@ public:
         void relax(const scalar alpha, const lumpedPointState& prev);
 
         //- Read input as dictionary content
-        bool readData(Istream& is);
+        bool readData
+        (
+            Istream& is,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
         //- Output as dictionary content
         bool writeData(Ostream& os) const;
@@ -191,24 +244,42 @@ public:
         void writeDict(Ostream& os) const;
 
         //- Read input as plain content
-        bool readPlain(Istream& is);
+        bool readPlain
+        (
+            Istream& is,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
         //- Output as plain content
         void writePlain(Ostream& os) const;
 
         //- Read input from file (master process only) using specified format
-        bool readData(const inputFormatType& fmt, const fileName& file);
+        bool readData
+        (
+            const inputFormatType fmt,
+            const fileName& file,
+            const quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+            const bool degrees = false
+        );
 
-        //- Output as VTK file for debugging/visualization
+        //- Output points/rotations as VTK file for debugging/visualization
-        //  The points are joined as lines, the rotation is visualized
+        //  The points are written as vertices, rotation as a triangle
-        //  by planes, write as VTK PolyData format.
+        void writeVTP
-        void writeVTP(const fileName& outputFile, const vector& axis) const;
+        (
+            const fileName& outputFile,
+            const labelListList& lines = labelListList(),
+            const labelList& pointIds = labelList::null()
+        ) const;
 
 
     // Member Operators
 
         //- Assignment operator
         void operator=(const lumpedPointState& rhs);
+
+        //- Shift points by specified origin
+        void operator+=(const point& origin);
 };
 
 

src/lumpedPointMotion/state/lumpedPointStateI.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -66,4 +66,17 @@ inline const Foam::tensorField& Foam::lumpedPointState::rotations() const
 }
 
 
+inline Foam::quaternion::eulerOrder
+Foam::lumpedPointState::rotationOrder() const
+{
+    return order_;
+}
+
+
+inline bool Foam::lumpedPointState::degrees() const
+{
+    return degrees_;
+}
+
+
 // ************************************************************************* //

src/lumpedPointMotion/state/lumpedPointStateWriter.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,30 +28,23 @@ License
 #include "lumpedPointState.H"
 #include "OFstream.H"
 #include "sliceRange.H"
-#include "axesRotation.H"
-#include "coordinateSystem.H"
 #include "foamVtkOutput.H"
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-// file-local
-const static Foam::FixedList<Foam::point, 4> standardCorners
-{
-    {-0.1, -0.1, 0},
-    {+0.1, -0.1, 0},
-    {+0.1, +0.1, 0},
-    {-0.1, +0.1, 0}
-};
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::lumpedPointState::writeVTP
 (
     const fileName& outputFile,
-    const vector& axis
+    const labelListList& lines,
+    const labelList& pointIds
 ) const
 {
+    if (!Pstream::master())
+    {
+        // No extra information available from slaves, write on master only.
+        return;
+    }
+
     // local-to-global transformation tensors
     const tensorField& localToGlobal = rotations();
 
@@ -68,28 +61,52 @@ void Foam::lumpedPointState::writeVTP
         .beginVTKFile<vtk::fileTag::POLY_DATA>();
 
     //
-    // The 'spine' of lumped mass points
+    // Lumped mass points and connections,
+    // with triangles to visualize location/rotation
     //
     {
+        const label nPoints = 3*points_.size();  // 3 points per triangle
+        const label nPolys  = points_.size();
+
         format()
             .tag
             (
                 vtk::fileTag::PIECE,
-                vtk::fileAttr::NUMBER_OF_POINTS, points_.size(),
+                vtk::fileAttr::NUMBER_OF_POINTS, nPoints,
                 vtk::fileAttr::NUMBER_OF_VERTS,  points_.size(),
-                vtk::fileAttr::NUMBER_OF_LINES,  1
+                vtk::fileAttr::NUMBER_OF_LINES,  lines.size(),
+                vtk::fileAttr::NUMBER_OF_POLYS,  nPolys
             );
 
         // 'points'
         {
-            const uint64_t payLoad = vtk::sizeofData<float, 3>(points_.size());
+            const uint64_t payLoad = vtk::sizeofData<float, 3>(nPoints);
 
             format()
                 .tag(vtk::fileTag::POINTS)
                 .beginDataArray<float, 3>(vtk::dataArrayAttr::POINTS);
 
             format().writeSize(payLoad);
+
+            // The lumped points first
             vtk::writeList(format(), points_);
+
+            // Other points (for the triangles) next
+            forAll(points_, posi)
+            {
+                const point& origin = points_[posi];
+                const tensor& rotTensor =
+                (
+                    posi < localToGlobal.size()
+                  ? localToGlobal[posi]
+                  : pTraits<tensor>::I
+                );
+
+                // Local-to-global rotation and translation
+                vtk::write(format(), 2*visLength*rotTensor.cx() + origin);
+                vtk::write(format(), 1*visLength*rotTensor.cy() + origin);
+            }
+
             format().flush();
 
             format()
@@ -140,16 +157,25 @@ void Foam::lumpedPointState::writeVTP
         // <Lines>
         format().tag(vtk::fileTag::LINES);
 
+        label nLinePoints = 0;
+        for (const labelList& linePoints : lines)
+        {
+            nLinePoints += linePoints.size();
+        }
+
         //
         // 'connectivity'
         //
         {
-            const uint64_t payLoad = vtk::sizeofData<label>(points_.size());
+            const uint64_t payLoad = vtk::sizeofData<label>(nLinePoints);
 
             format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
             format().writeSize(payLoad);
 
-            vtk::writeIdentity(format(), points_.size());
+            for (const labelList& linePoints : lines)
+            {
+                vtk::writeList(format(), linePoints);
+            }
 
             format().flush();
 
@@ -158,102 +184,52 @@ void Foam::lumpedPointState::writeVTP
 
         //
         // 'offsets'  (connectivity offsets)
-        // = single line
+        // = N lines
         //
         {
-            const uint64_t payLoad = vtk::sizeofData<label>(1);
+            const uint64_t payLoad = vtk::sizeofData<label>(lines.size());
 
             format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
             format().writeSize(payLoad);
 
-            format().write(points_.size());
+            nLinePoints = 0;
+            for (const labelList& linePoints : lines)
+            {
+                nLinePoints += linePoints.size();
+                format().write(nLinePoints);
+            }
+
             format().flush();
 
             format().endDataArray();
         }
 
         format().endTag(vtk::fileTag::LINES);
-        format().endPiece();
+        // </Lines>
-    }
 
-    // Standard corners in local axis
-    FixedList<point, 4> corners;
-
-    {
-        coordinateRotations::axes orient(axis);
-        coordinateSystem cornerTransform(orient);
-
-        forAll(standardCorners, corni)
-        {
-            corners[corni] = cornerTransform.transform(standardCorners[corni]);
-        }
-    }
-
-    //
-    // Planes to visualize location/rotation
-    //
-    {
-        const label nPoints = 4*points_.size();  // 4 points per quad
-        const label nPolys  = points_.size();
-
-        format()
-            .tag
-            (
-                vtk::fileTag::PIECE,
-                vtk::fileAttr::NUMBER_OF_POINTS, nPoints,
-                vtk::fileAttr::NUMBER_OF_POLYS,  nPolys
-            );
-
-        // 'points'
-        {
-            const uint64_t payLoad = vtk::sizeofData<float, 3>(nPoints);
-
-            format()
-                .tag(vtk::fileTag::POINTS)
-                .beginDataArray<float, 3>(vtk::dataArrayAttr::POINTS);
-
-            format().writeSize(payLoad);
-
-            forAll(points_, posI)
-            {
-                const point& origin = points_[posI];
-                const tensor& rotTensor =
-                (
-                    posI < localToGlobal.size()
-                  ? localToGlobal[posI]
-                  : pTraits<tensor>::I
-                );
-
-
-                for (const point& cornerPt : corners)
-                {
-                    // Local-to-global rotation and translation
-                    const point pt = (rotTensor & cornerPt) + origin;
-
-                    vtk::write(format(), pt);
-                }
-            }
-
-            format().flush();
-
-            format()
-                .endDataArray()
-                .endTag(vtk::fileTag::POINTS);
-        }
 
         // <Polys>
         format().tag(vtk::fileTag::POLYS);
 
         //
-        // 'connectivity' - 4 points (ie, quad)
+        // 'connectivity' - 3 points (ie, tri)
+        // origins appear first, followed by a point pair for each triangle
+        // Eg,
+        // - tri 0: (0 N   N+1)
+        // - tri 1: (1 N+2 N+3)
         //
         {
-            const uint64_t payLoad = vtk::sizeofData<label>(4*nPolys);
+            const uint64_t payLoad = vtk::sizeofData<label>(3*nPolys);
 
             format().beginDataArray<label>(vtk::dataArrayAttr::CONNECTIVITY);
             format().writeSize(payLoad);
 
-            vtk::writeIdentity(format(), 4*nPolys);
+            for (label pointi=0, nei=nPolys; pointi < nPolys; ++pointi)
+            {
+                format().write(pointi);
+                format().write(nei); ++nei;
+                format().write(nei); ++nei;
+            }
 
             format().flush();
 
@@ -262,7 +238,7 @@ void Foam::lumpedPointState::writeVTP
 
         //
         // 'offsets'  (connectivity offsets)
-        // = single quad
+        // = single tri
         //
         {
             const uint64_t payLoad = vtk::sizeofData<label>(nPolys);
@@ -270,7 +246,7 @@ void Foam::lumpedPointState::writeVTP
             format().beginDataArray<label>(vtk::dataArrayAttr::OFFSETS);
             format().writeSize(payLoad);
 
-            for (const label off : sliceRange(4, nPolys, 4))
+            for (const label off : sliceRange(3, nPolys, 3))
             {
                 format().write(off);
             }
@@ -280,17 +256,27 @@ void Foam::lumpedPointState::writeVTP
         }
 
         format().endTag(vtk::fileTag::POLYS);
+        // </Polys>
 
-#if 0
+
+        // CELL_DATA
         format().beginCellData();
 
-        // zone Id
+        // point id
         {
-            const uint64_t payLoad = vtk::sizeofData<label>(nPolys);
+            const uint64_t payLoad =
+                vtk::sizeofData<label>(points_.size() + lines.size());
 
-            format().beginDataArray<label>("zoneId");
+            format().beginDataArray<label>("pointId");
             format().writeSize(payLoad);
 
+            // <Verts>
+            vtk::writeIdentity(format(), points_.size());
+
+            // <Lines>
+            vtk::write(format(), label(-1), lines.size());
+
+            // <Poly>
             vtk::writeIdentity(format(), nPolys);
 
             format().flush();
@@ -298,16 +284,99 @@ void Foam::lumpedPointState::writeVTP
             format().endDataArray();
         }
 
+        // original id
+        if (pointIds.size() == points_.size())
+        {
+            const uint64_t payLoad =
+                vtk::sizeofData<label>(points_.size() + lines.size());
+
+            format().beginDataArray<label>("originalId");
+            format().writeSize(payLoad);
+
+            // <Verts>
+            vtk::writeList(format(), pointIds);
+
+            // <Lines>
+            vtk::write(format(), label(-1), lines.size());
+
+            // <Poly>
+            vtk::writeList(format(), pointIds);
+
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        // line id
+        {
+            const uint64_t payLoad =
+                vtk::sizeofData<label>(points_.size() + lines.size());
+
+            format().beginDataArray<label>("lineId");
+            format().writeSize(payLoad);
+
+            // <Verts>
+            vtk::write(format(), label(-1), points_.size());
+
+            // <Lines>
+            vtk::writeIdentity(format(), lines.size());
+
+            // <Poly>
+            vtk::write(format(), label(-1), nPolys);
+
+            format().flush();
+
+            format().endDataArray();
+        }
+
         format().endCellData();
-#endif
+
+
+        // POINT_DATA
+        format().beginPointData();
+
+        // point id
+        {
+            const uint64_t payLoad = vtk::sizeofData<label>(nPoints);
+
+            format().beginDataArray<label>("pointId");
+            format().writeSize(payLoad);
+
+            // The lumped points first
+            vtk::writeIdentity(format(), points_.size());
+
+            // Tag other (triangle) points as -1
+            vtk::write(format(), label(-1), 2*points_.size());
+
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        // original id
+        if (pointIds.size() == points_.size())
+        {
+            const uint64_t payLoad = vtk::sizeofData<label>(nPoints);
+
+            format().beginDataArray<label>("originalId");
+            format().writeSize(payLoad);
+
+            // The lumped points first
+            vtk::writeList(format(), pointIds);
+
+            // Tag other (triangle) points as -1
+            vtk::write(format(), label(-1), 2*points_.size());
+
+            format().flush();
+
+            format().endDataArray();
+        }
+
+        format().endPointData();
 
         format().endPiece();
     }
 
-    // Finally
-    // could add a 'ghost' level above to visualize extrapolated values
-    // draw as two triangles to distinguish from real levels ...
-
     format().endTag(vtk::fileTag::POLY_DATA)
         .endVTKFile();
 }

src/lumpedPointMotion/tools/lumpedPointTools.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2017 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,14 +28,13 @@ License
 #include "IFstream.H"
 #include "IOobjectList.H"
 #include "volFields.H"
-
+#include "points0MotionSolver.H"
 #include "lumpedPointDisplacementPointPatchVectorField.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    // file-scope
     template<class GeoFieldType>
     static autoPtr<GeoFieldType> loadPointField
     (
@@ -72,21 +71,35 @@ namespace Foam
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 Foam::List<Foam::lumpedPointStateTuple>
-Foam::lumpedPointTools::lumpedPointStates(Istream& is)
+Foam::lumpedPointTools::lumpedPointStates
+(
+    const dictionary& dict,
+    quaternion::eulerOrder rotOrder,
+    bool degrees
+)
 {
-    dictionary contents(is);
+    quaternion::eulerOrderNames.readIfPresent
-    List<dictionary> entries(contents.lookup("response"));
+    (
+        "rotationOrder",
+        dict,
+        rotOrder
+    );
+
+    dict.readIfPresent("degrees", degrees);
+
+    Info<<"Reading states\n";
+    List<dictionary> entries(dict.lookup("response"));
 
     DynamicList<Tuple2<scalar, lumpedPointState>> states(entries.size());
 
-    for (const dictionary& dict : entries)
+    for (const dictionary& subDict : entries)
     {
         states.append
         (
             lumpedPointStateTuple
             (
-                dict.get<scalar>("time"),
+                subDict.get<scalar>("time"),
-                lumpedPointState(dict)
+                lumpedPointState(subDict)
             )
         );
     }
@@ -96,35 +109,38 @@ Foam::lumpedPointTools::lumpedPointStates(Istream& is)
 
 
 Foam::List<Foam::lumpedPointStateTuple>
-Foam::lumpedPointTools::lumpedPointStates(const fileName& file)
+Foam::lumpedPointTools::lumpedPointStates
+(
+    Istream& is,
+    quaternion::eulerOrder rotOrder,
+    bool degrees
+)
+{
+    dictionary dict(is);
+    return lumpedPointStates(dict, rotOrder, degrees);
+}
+
+
+Foam::List<Foam::lumpedPointStateTuple>
+Foam::lumpedPointTools::lumpedPointStates
+(
+    const fileName& file,
+    quaternion::eulerOrder rotOrder,
+    bool degrees
+)
 {
     IFstream is(file);
-    return lumpedPointStates(is);
+    return lumpedPointStates(is, rotOrder, degrees);
 }
 
 
 Foam::pointIOField
 Foam::lumpedPointTools::points0Field(const polyMesh& mesh)
 {
-    pointIOField pts
+    return pointIOField(points0MotionSolver::points0IO(mesh));
-    (
-        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)
 {
@@ -141,19 +157,127 @@ Foam::labelList Foam::lumpedPointTools::lumpedPointPatchList
 
     pointMesh pMesh(mesh);
 
-    autoPtr<pointVectorField> displacePtr = loadPointField<pointVectorField>
+    autoPtr<pointVectorField> displacePtr =
-    (
+        loadPointField<pointVectorField>
-        pMesh,
+        (
-        objects0.findObject("pointDisplacement")
+            pMesh,
-    );
+            objects0.findObject("pointDisplacement")
+        );
 
-    if (!displacePtr.valid())
+    if (!displacePtr)
     {
-        Info<< "no valid pointDisplacement" << endl;
+        Info<< "No valid pointDisplacement" << endl;
         return labelList();
     }
 
-    return lumpedPointPatchList(displacePtr());
+    return lumpedPointPatchList(*displacePtr);
+}
+
+
+Foam::label
+Foam::lumpedPointTools::setPatchControls
+(
+    const pointVectorField& pvf,
+    const pointField& points0
+)
+{
+    return
+        lumpedPointDisplacementPointPatchVectorField::setPatchControls
+        (
+            pvf,
+            points0
+        );
+}
+
+
+Foam::label Foam::lumpedPointTools::setPatchControls
+(
+    const fvMesh& mesh,
+    const pointField& points0
+)
+{
+    IOobjectList objects0(mesh, "0");
+
+    pointMesh pMesh(mesh);
+
+    autoPtr<pointVectorField> displacePtr =
+        loadPointField<pointVectorField>
+        (
+            pMesh,
+            objects0.findObject("pointDisplacement")
+        );
+
+    if (!displacePtr)
+    {
+        Info<< "No valid pointDisplacement" << endl;
+        return 0;
+    }
+
+    return setPatchControls(*displacePtr, points0);
+}
+
+
+Foam::label Foam::lumpedPointTools::setPatchControls
+(
+    const fvMesh& mesh
+)
+{
+    pointIOField points0(points0Field(mesh));
+
+    return setPatchControls(mesh, points0);
+}
+
+
+Foam::label Foam::lumpedPointTools::setInterpolators
+(
+    const pointVectorField& pvf,
+    const pointField& points0
+)
+{
+    return
+        lumpedPointDisplacementPointPatchVectorField::setInterpolators
+        (
+            pvf,
+            points0
+        );
+}
+
+
+Foam::label Foam::lumpedPointTools::setInterpolators
+(
+    const fvMesh& mesh,
+    const pointField& points0
+)
+{
+    IOobjectList objects0(mesh, "0");
+
+    pointMesh pMesh(mesh);
+
+    autoPtr<pointVectorField> displacePtr =
+        loadPointField<pointVectorField>
+        (
+            pMesh,
+            objects0.findObject("pointDisplacement")
+        );
+
+    if (!displacePtr)
+    {
+        Info<< "No valid pointDisplacement" << endl;
+        return 0;
+    }
+
+    return setInterpolators(*displacePtr, points0);
+}
+
+
+Foam::label Foam::lumpedPointTools::setInterpolators
+(
+    const fvMesh& mesh
+)
+{
+    pointIOField points0(points0Field(mesh));
+
+    return setInterpolators(mesh, points0);
 }
 
 

src/lumpedPointMotion/tools/lumpedPointTools.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,11 +39,10 @@ SourceFiles
 #define lumpedPointTools_H
 
 #include "labelList.H"
-#include "polyMesh.H"
+#include "fvMesh.H"
 #include "pointMesh.H"
 #include "pointFields.H"
 #include "Tuple2.H"
-
 #include "lumpedPointState.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -51,6 +50,7 @@ SourceFiles
 namespace Foam
 {
 
+// Typedefs
 typedef Tuple2<scalar, lumpedPointState> lumpedPointStateTuple;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -58,23 +58,63 @@ typedef Tuple2<scalar, lumpedPointState> lumpedPointStateTuple;
 namespace lumpedPointTools
 {
 
+//- Load a list of states from a dictionary
+List<lumpedPointStateTuple> lumpedPointStates
+(
+    const dictionary& dict,
+    quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+    bool degrees = false
+);
+
 //- Load a list of states from an Istream
-List<lumpedPointStateTuple> lumpedPointStates(Istream& is);
+List<lumpedPointStateTuple> lumpedPointStates
+(
+    Istream& is,
+    quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+    bool degrees = false
+);
 
 //- Load a list of states from a file
-List<lumpedPointStateTuple> lumpedPointStates(const fileName& file);
+List<lumpedPointStateTuple> lumpedPointStates
-
+(
+    const fileName& file,
+    quaternion::eulerOrder rotOrder = quaternion::eulerOrder::ZXZ,
+    bool degrees = false
+);
 
 //- Return the 0 or constant points field
 pointIOField points0Field(const polyMesh& mesh);
 
-//- Return the patch-ids associated a "lumpedPointDisplacement" type
+//- Return the patch-ids associated with a "lumpedPointDisplacement" type
 labelList lumpedPointPatchList(const pointVectorField& pvf);
 
 //- Get the "pointDisplacement" at time 0 and use that to determine which
-//  patches have a "lumpedPointDisplacement" type
+//- patches have a "lumpedPointDisplacement" type
 labelList lumpedPointPatchList(const polyMesh& mesh);
 
+//- Return the patch-ids associated with a "lumpedPointDisplacement" type
+label setPatchControls(const pointVectorField& pvf, const pointField& points0);
+
+//- Get the "pointDisplacement" at time 0 and use that to determine which
+//- patches have a "lumpedPointDisplacement" type
+label setPatchControls(const fvMesh& mesh, const pointField& points0);
+
+//- Get the "pointDisplacement" at time 0 and use that to determine which
+//- patches have a "lumpedPointDisplacement" type
+label setPatchControls(const fvMesh& mesh);
+
+
+//- Return the patch-ids associated with a "lumpedPointDisplacement" type
+label setInterpolators(const pointVectorField& pvf, const pointField& points0);
+
+//- Get the "pointDisplacement" at time 0 and use that to determine which
+//- patches have a "lumpedPointDisplacement" type
+label setInterpolators(const fvMesh& mesh, const pointField& points0);
+
+//- Get the "pointDisplacement" at time 0 and use that to determine which
+//- patches have a "lumpedPointDisplacement" type
+label setInterpolators(const fvMesh& mesh);
+
 } // End namespace lumpedPointTools
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

tutorials/incompressible/lumpedPointMotion/building/Allclean

@@ -1,5 +1,8 @@
 #!/bin/sh
-cd "${0%/*}" || exit    # Run from this directory
+cd "${0%/*}" || exit                                # Run from this directory
+#------------------------------------------------------------------------------
+
+(cd code && ./Allclean)
 
 (cd steady && ./Allclean)
 

tutorials/incompressible/lumpedPointMotion/building/Allrun

@@ -3,94 +3,45 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
-# 1) First run steady-state to establish a good initial field.
+# 1) Run steady-state to establish a good initial field
-# 2) Copy the latest state-state results for the transient case,
+# 2) Copy state-state results -> transient case,
 #    but need to copy the pointDisplacement from the 0/ directory
 #    since it will not have been used for the steady-state case
 # 3) Relocate this initial solution to coincide with the first deltaT
 #    to avoid overwriting the 0/ directory at all.
 
-#
-# copyParallelPointDisplacement caseDir timeName
-#
-# Copy pointDisplacement from caseDir/0/ to caseDir/timeName/
-#
-copyParallelPointDisplacement()
-{
-    local src=$1
-    local dstTime=$2
-    local file=pointDisplacement
-
-    [ -d "$src" ] || {
-        echo "Error: no directory: $src"
-        return 1
-    }
-
-    # Copy select directories
-    echo "    copy processor '$file' from 0/ -> $dstTime"
-    if [ -n "$dstTime" ]
-    then
-    (
-        cd "$src" || exit
-
-        for proc in processor*
-        do
-            [ -d "$proc/0" -a -d "$proc/$dstTime" ] && \
-                cp $proc/0/$file $proc/$dstTime/$file
-        done
-    )
-    else
-        echo "    no destination time"
-    fi
-
-    # Restart from latestTime
-    foamDictionary $src/system/controlDict \
-        -entry startFrom -set latestTime
-
-    deltaT=$(foamDictionary $src/system/controlDict -entry deltaT -value)
-    latestTime=$(foamListTimes -case $src -noZero -latestTime -processor)
-
-    # Restart using steady results as first deltaT interval
-    echo "deltaT=$deltaT latestTime=$latestTime"
-    if [ -n "$latestTime" -a "$deltaT" != "$latestTime" ]
-    then
-    (
-        cd "$src" || exit
-
-        for proc in processor*
-        do
-            if [ -d "$proc/$latestTime" -a ! -d "$proc/$deltaT" ]
-            then
-                mv $proc/$latestTime $proc/$deltaT
-                rm -rf "$proc/$deltaT/uniform"
-            fi
-        done
-    )
-    fi
-
-    return 0
-}
-
-
 # Do steady-state case
 (cd steady && foamRunTutorials)
 
 if notTest "$@"
 then
-    latestTime=$(\cd steady && foamListTimes -noZero -latestTime -processor)
+    if canCompile
+    then
+        (cd code && wmake)
+    else
+        exit 0
+    fi
 
-    # Clone the steady-state case to transient
+    . files/RunFunctions
-    cloneParallelCase steady transient 0 $latestTime
 
-    copyParallelPointDisplacement transient $latestTime
+    caseName="transient"
+
+    latestTime=$(foamListTimes -case steady -noZero -latestTime -processor)
+
+    # Clone steady-state case to transient
+    cloneParallelCase steady "$caseName" 0 "$latestTime"
+
+    copyParallelPointDisplacement "$caseName" "$latestTime"
 
     # Adjust application (from simpleFoam -> pimpleFoam)
-    foamDictionary transient/system/controlDict \
+    foamDictionary "$caseName"/system/controlDict \
         -entry application -set pimpleFoam
 
-    # Do the transient case
+    # Copy/link support files
-    \cp files/Allrun.transient transient/Allrun
+    linkFiles files "$caseName"
-    (\cd transient && foamRunTutorials)
+
+    # Run
+    "$caseName/Allrun.$caseName" $*
 fi
 
 #------------------------------------------------------------------------------

tutorials/incompressible/lumpedPointMotion/building/Allrun.move

@@ -4,62 +4,32 @@ cd "${0%/*}" || exit                                # Run from this directory
 #------------------------------------------------------------------------------
 
 # 1) Run meshing
-# 2) Reconstruct
+# 2) Test input zones and movement
-# 3) Test input zones and movement
 
-#
+# Meshing
-# linkParallelCase srcDir dstDir
-#
-linkParallelCase()
-{
-    local src=$1
-    local dst=$2
-    shift 2
-
-    if [ -e "$dst" ]
-    then
-        echo "Case already linked: remove case directory $dst prior to linking"
-        return 1
-    elif [ ! -d "$src" ]
-    then
-        echo "Error: no directory to link: $src"
-        return 1
-    fi
-
-    echo "Linking $dst parallel case from $src"
-    mkdir $dst
-
-    # Copy system - may wish to change things
-    for i in system 0
-    do
-        echo "    copy $i/"
-        ( cd $dst && cp -r ../$src/$i . )
-    done
-
-    echo "    link constant/"
-    ( cd $dst && ln -sf ../$src/constant . )
-
-    echo "    link processor*/ with $# times: $@"
-    for proc in $(cd $src && \ls -d processor*)
-    do
-        ( cd $dst && ln -sf ../$src/$proc . )
-    done
-
-    return 0
-}
-
-
-# Do steady-state case
 (cd steady && ./Allrun.pre)
 
-if notTest "$@"
+if notTest
 then
-    # Copy/link the steady-state case to movement
+    if canCompile
-    linkParallelCase steady movement
+    then
+        (cd code && wmake)
+    else
+        exit 0
+    fi
 
-    # Test movement
+    . files/RunFunctions
-    \cp files/Allrun.movement movement/Allrun
+
-    (cd movement && foamRunTutorials)
+    caseName="movement"
+
+    # Copy/link the steady-state case to movement
+    linkParallelCase steady "$caseName"
+
+    # Copy/link support files
+    linkFiles files "$caseName"
+
+    # Run
+    "$caseName/Allrun.$caseName" $*
 fi
 
 #------------------------------------------------------------------------------

tutorials/incompressible/lumpedPointMotion/building/code/Allclean

@@ -0,0 +1,13 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+#------------------------------------------------------------------------------
+
+wclean
+
+# Remove executable
+rm -f building-motion
+
+# Remove known output/debug files
+rm -f *.txt *.vtp *.vtp.series
+
+#------------------------------------------------------------------------------

tutorials/incompressible/lumpedPointMotion/building/code/Make/files

@@ -0,0 +1,3 @@
+building-motion.C
+
+EXE = $(PWD)/building-motion

tutorials/incompressible/lumpedPointMotion/building/code/Make/options

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

tutorials/incompressible/lumpedPointMotion/building/code/building-motion.C

@@ -0,0 +1,682 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    building-motion
+
+Description
+    Forced oscillation code for fluid-structure interface check.
+    Generates position and rotation angle of each node.
+    The top displacements of the target building calculated as sin() function.
+    Horizontal displacements of each node interpolated to the vertical
+    direction according to cantilever beam theory.
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "Fstream.H"
+#include "unitConversion.H"
+#include "foamVtkSeriesWriter.H"
+#include "lumpedPointTools.H"
+#include "lumpedPointState.H"
+#include "lumpedPointIOMovement.H"
+
+using namespace Foam;
+
+//- Oscillator generator
+class position_generator
+{
+    // Private Member Functions
+
+        //- Calculate position/rotation at given time
+        lumpedPointState calc(scalar currTime) const
+        {
+            // Point positions
+            pointField points_(nDivisions+1, Zero);
+
+            // Point rotations
+            vectorField angles_(nDivisions+1, Zero);
+
+            // Set node heights (z)
+            forAll(points_, divi)
+            {
+                points_[divi].z() = (height * divi)/scalar(nDivisions);
+            }
+
+            const vector sines
+            (
+                Foam::sin(2*constant::mathematical::pi * currTime/period.x()),
+                Foam::sin(2*constant::mathematical::pi * currTime/period.y()),
+                Foam::sin(2*constant::mathematical::pi * currTime/period.z())
+            );
+
+            for (label divi = 1; divi <= nDivisions; ++divi)
+            {
+                const scalar zpos = points_[divi].z();
+                const scalar height1 = (height - zpos);
+
+                const scalar pos_factor =
+                (
+                    1.0/3.0 / pow4(height)
+                  * (
+                        3*pow4(height)
+                      - 4*pow3(height)*(height1)
+                      + pow4(height1)
+                    )
+                );
+
+                const scalar ang_factor =
+                (
+                    1.0/3.0 / pow4(height)
+                  * (
+                        4*pow3(height)
+                      - 4*pow3(height1)
+                    )
+                );
+
+                vector here
+                (
+                    (amplitude.x() * sines.x() * pos_factor),
+                    (amplitude.y() * sines.y() * pos_factor),
+                    zpos // Z position is invariant
+                );
+
+                vector rot
+                (
+                    Foam::atan(amplitude.x() * sines.x() * ang_factor),
+                    Foam::atan(amplitude.y() * sines.y() * ang_factor),
+                    Foam::atan(amplitude.z() * sines.z() * ang_factor)
+                );
+
+                // Assign
+                points_[divi] = here;
+
+                // The x<->y swap is intentional
+                angles_[divi] = vector{rot.y(), rot.x(), rot.z()};
+            }
+
+            // Return as lumpedPoint state
+            return lumpedPointState{points_, angles_};
+        }
+
+
+public:
+
+    // Control parameters
+
+        // The number of oscillating nodes
+        label nDivisions = 10;
+
+        // Height of target building [m]
+        scalar height = 0.5;
+
+        // Proper period (sec)
+        vector period = vector{1, 0.5, 1};
+
+        // Amplitude
+        vector amplitude = vector{0.03, 0.05, 0.3};
+
+
+    // Constructors
+
+        //- Default construct
+        position_generator() = default;
+
+
+    // Member Functions
+
+        //- Calculate position/rotation at given time
+        lumpedPointState state(const scalar currTime) const
+        {
+            return calc(currTime);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Forced oscillation code for fluid-structure interface check."
+        " Generates position and rotation angle of each node."
+        " The top displacements of the target building calculated as sin()"
+        " function."
+        " Horizontal displacements of each node interpolated to the vertical"
+        " direction according to cantilever beam theory."
+     );
+
+    argList::noBanner();
+    argList::noParallel();
+
+    // Geometric
+    argList::addOption
+    (
+        "nodes",
+        "N",
+        "The number of oscillating nodes (default: 10)"
+    );
+    argList::addOption
+    (
+        "height",
+        "value",
+        "Height of target building (m)"
+    );
+    argList::addOption
+    (
+        "period",
+        "(time time time)",
+        "The proper period (sec)"
+    );
+    argList::addOption
+    (
+        "amplitude",
+        "(value value value)",
+        "The amplitude"
+    );
+
+    // Time control
+    argList::addOption
+    (
+        "time",
+        "value",
+        "The time to use"
+    );
+    argList::addOption
+    (
+        "deltaT",
+        "value",
+        "The time increment for multiple time loops"
+    );
+    argList::addOption
+    (
+        "nTimes",
+        "value",
+        "The number of time loops"
+    );
+
+    // Query, output
+    argList::addBoolOption
+    (
+        "query",
+        "Report values only and exit"
+    );
+    argList::addOption
+    (
+        "output",
+        "file",
+        "write to file, with header"
+    );
+
+    argList::addOption
+    (
+        "scale",
+        "factor",
+        "Scaling factor for movement (default: 1)"
+    );
+    argList::addOption
+    (
+        "visual-length",
+        "len",
+        "Visualization length for planes (visualized as triangles)"
+    );
+
+    // Run controls
+    argList::addBoolOption
+    (
+        "dry-run",
+        "Test movement without a mesh"
+    );
+    argList::addBoolOption
+    (
+        "removeLock",
+        "Remove lock-file on termination of slave"
+    );
+    argList::addBoolOption
+    (
+        "slave",
+        "Invoke as a slave responder for testing"
+    );
+    #include "setRootCase.H"
+
+
+    // The oscillator
+    position_generator gen;
+    args.readIfPresent("nodes", gen.nDivisions);
+    args.readIfPresent("height", gen.height);
+    args.readIfPresent("period", gen.period);
+    args.readIfPresent("amplitude", gen.amplitude);
+
+
+    // Control parameters
+    const bool dryrun = args.found("dry-run");
+    const bool slave = args.found("slave");
+    const bool removeLock = args.found("removeLock");
+
+    const bool optQuery = args.found("query");
+    const fileName outputFile(args.getOrDefault<fileName>("output", ""));
+
+    const scalar relax = args.getOrDefault<scalar>("scale", 1);
+
+    args.readIfPresent("visual-length", lumpedPointState::visLength);
+
+    // Time parameters
+    scalar currTime = args.getOrDefault<scalar>("time", 0);
+    const scalar deltaT = args.getOrDefault("deltaT", 0.001);
+    const label nTimes = args.getOrDefault<label>("nTimes", 1);
+
+    // Loop handling for slave
+    const bool infiniteLoop = slave && !args.found("nTimes");
+
+
+    // ----------------------------------------------------------------------
+    // Slave mode
+    // ----------------------------------------------------------------------
+
+    if (slave)
+    {
+        Info<< "Running as slave responder" << endl;
+
+        if (Pstream::parRun())
+        {
+            FatalErrorInFunction
+                << "Running as slave responder is not permitted in parallel"
+                << nl
+                << exit(FatalError);
+        }
+
+        #include "createTime.H"
+
+        // Create movement without a mesh
+        autoPtr<lumpedPointIOMovement> movementPtr =
+            lumpedPointIOMovement::New(runTime);
+
+        if (!movementPtr)
+        {
+            Info<< "No valid movement found" << endl;
+            return 1;
+        }
+        auto& movement = *movementPtr;
+
+        // Reference state0
+        const lumpedPointState& state0 = movement.state0();
+
+        externalFileCoupler& coupler = movement.coupler();
+
+        for (label timei = 0; infiniteLoop || (timei < nTimes); ++timei)
+        {
+            Info<< args.executable() << ": waiting for master" << endl;
+
+            // Wait for master, but stop if status=done was seen
+            if (!coupler.waitForMaster())
+            {
+                Info<< args.executable()
+                    << ": stopping status=done was detected" << endl;
+                break;
+            }
+
+            scalar timeValue = currTime;
+
+            if (infiniteLoop)
+            {
+                // Get output file
+                IFstream is(coupler.resolveFile(movement.outputName()));
+
+                dictionary dict;
+                is >> dict;
+
+                timeValue = dict.get<scalar>("time");
+            }
+
+            lumpedPointState state(gen.state(timeValue));
+            state.relax(relax, state0);
+
+            // Generate input for OpenFOAM
+            {
+                OFstream os(coupler.resolveFile(movement.inputName()));
+                if
+                (
+                    movement.inputFormat()
+                 == lumpedPointState::inputFormatType::PLAIN
+                )
+                {
+                    state.writePlain(os);
+                }
+                else
+                {
+                    os.writeEntry("time", timeValue);
+                    state.writeDict(os);
+                }
+            }
+
+            Info<< args.executable()
+                << ": updating state " << timei
+                << " - switch to master"
+                << endl;
+
+            // Let OpenFOAM know that it can continue
+            coupler.useMaster();
+
+            currTime += deltaT;
+        }
+
+        if (removeLock)
+        {
+            Info<< args.executable() << ": removing lock file" << endl;
+            coupler.useSlave();  // This removes the lock-file
+        }
+
+        Info<< args.executable() << ": finishing" << nl;
+
+        Info<< "\nEnd\n" << endl;
+        return 0;
+    }
+
+
+    // ----------------------------------------------------------------------
+    // dry-run
+    // ----------------------------------------------------------------------
+
+    if (dryrun)
+    {
+        Info<< "dry-run: creating states only" << nl;
+
+        autoPtr<Time> runTimePtr;
+        autoPtr<lumpedPointIOMovement> movementPtr;
+
+        const bool throwingIOError = FatalIOError.throwExceptions();
+        const bool throwingError = FatalError.throwExceptions();
+
+        try
+        {
+            Info<< "Create time" << flush;
+
+            runTimePtr = Time::New(args);
+
+            // Create movement without a mesh
+            movementPtr = lumpedPointIOMovement::New(*runTimePtr);
+        }
+        catch (...)
+        {
+            Info<< " ... failed (optional for dry-run)";
+        }
+        Info<< nl << endl;
+
+        FatalError.throwExceptions(throwingError);
+        FatalIOError.throwExceptions(throwingIOError);
+
+        if (!movementPtr)
+        {
+            Info<< "No time, run without movement information\n" << endl;
+        }
+
+        const lumpedPointState state0(gen.state(0));
+
+        vtk::seriesWriter stateSeries;
+
+        for
+        (
+            label timei = 0, outputCount = 0;
+            timei < nTimes;
+            ++timei
+        )
+        {
+            lumpedPointState state(gen.state(currTime));
+            state.relax(relax, state0);
+
+            Info<< "output [" << timei << '/' << nTimes << ']';
+
+            // State/response = what comes back from FEM
+            {
+                const word outputName =
+                    word::printf("state_%06d.vtp", outputCount);
+
+                Info<< "  " << outputName;
+
+                if (movementPtr)
+                {
+                    movementPtr->writeStateVTP(state, outputName);
+                }
+                else
+                {
+                    state.writeVTP(outputName);
+                }
+                stateSeries.append(outputCount, outputName);
+            }
+
+            Info<< endl;
+
+            ++outputCount;
+            currTime += deltaT;
+        }
+
+        // Write file series
+        if (stateSeries.size())
+        {
+            Info<< nl << "write state.vtp.series" << nl;
+            stateSeries.write("state.vtp");
+        }
+
+        Info<< "\nEnd\n" << endl;
+        return 0;
+    }
+
+
+    // ----------------------------------------------------------------------
+    // Report values or generate a file
+    // ----------------------------------------------------------------------
+
+    if (optQuery || !outputFile.empty())
+    {
+        autoPtr<OFstream> osPtr;
+
+        if (!outputFile.empty())
+        {
+            osPtr.reset(new OFstream(outputFile));
+            auto& os = *osPtr;
+
+            os.precision(8);
+
+            // One file with everything, output using OpenFOAM syntax
+
+            IOobject::writeBanner(os)
+                << "FoamFile\n{\n"
+                << "    version     " << os.version() << ";\n"
+                << "    format      " << os.format() << ";\n"
+                << "    class       " << "dictionary" << ";\n"
+                << "    object      " << "response" << ";\n"
+                << "}\n";
+
+            IOobject::writeDivider(os) << nl;
+
+            os << "// angles are Euler angles z-x-z (intrinsic)" << nl;
+            os.writeEntry("degrees", "false");
+            os << nl;
+            os << "response" << nl;
+            os << '(' << nl;
+        }
+        else
+        {
+            Info().precision(8);
+        }
+
+
+        for (label timei = 0; timei < nTimes; ++timei)
+        {
+            lumpedPointState state(gen.state(currTime));
+
+            if (osPtr)
+            {
+                // Report position/angle
+                auto& os = *osPtr;
+
+                os.beginBlock();
+
+                os.writeEntry("time", currTime);
+
+                state.writeDict(os);
+
+                os.endBlock();
+            }
+            else
+            {
+                // Report position/angle
+                auto& os = Info();
+
+                os.writeEntry("time", currTime);
+                state.writeDict(os);
+            }
+
+            currTime += deltaT;
+        }
+
+
+        if (osPtr)
+        {
+            auto& os = *osPtr;
+
+            os << ')' << token::END_STATEMENT << nl;
+
+            IOobject::writeEndDivider(os);
+
+            Info<< "\nEnd\n" << endl;
+        }
+
+        return 0;
+    }
+
+
+    // ----------------------------------------------------------------------
+    // test patch movement
+    // ----------------------------------------------------------------------
+
+    #include "createTime.H"
+
+    runTime.setTime(instant(runTime.constant()), 0);
+
+    #include "createNamedMesh.H"
+
+    // Create movement with mesh
+    autoPtr<lumpedPointIOMovement> movementPtr =
+        lumpedPointIOMovement::New(mesh);
+
+    if (!movementPtr)
+    {
+        Info<< "No valid movement found" << endl;
+        return 1;
+    }
+    auto& movement = *movementPtr;
+
+    // Reference state0
+    const lumpedPointState& state0 = movement.state0();
+
+    pointIOField points0(lumpedPointTools::points0Field(mesh));
+
+    const label nPatches = lumpedPointTools::setPatchControls(mesh, points0);
+    if (!nPatches)
+    {
+        Info<< "No point patches with lumped movement found" << endl;
+        return 2;
+    }
+
+    Info<< "Lumped point patch controls set on "
+        << nPatches << " patches" << nl;
+
+    lumpedPointTools::setInterpolators(mesh, points0);
+
+
+    // Output vtk file series
+    vtk::seriesWriter stateSeries;
+    vtk::seriesWriter geomSeries;
+
+    // Initial geometry
+    movement.writeVTP("geom_init.vtp", state0, mesh, points0);
+
+    lumpedPointTools::setInterpolators(mesh);
+
+    for
+    (
+        label timei = 0, outputCount = 0;
+        timei < nTimes;
+        ++timei
+    )
+    {
+        lumpedPointState state(gen.state(currTime));
+
+        state += movement.origin();
+        movement.scalePoints(state);
+        state.relax(relax, state0);
+
+        Info<< "output [" << timei << '/' << nTimes << ']';
+
+        // State/response = what comes back from FEM
+        {
+            const word outputName =
+                word::printf("state_%06d.vtp", outputCount);
+
+            Info<< "  " << outputName;
+
+            movement.writeStateVTP(state, outputName);
+            stateSeries.append(outputCount, outputName);
+        }
+
+        {
+            const word outputName =
+                word::printf("geom_%06d.vtp", outputCount);
+
+            Info<< "  " << outputName;
+
+            movement.writeVTP(outputName, state, mesh, points0);
+            geomSeries.append(outputCount, outputName);
+        }
+
+        Info<< endl;
+
+        ++outputCount;
+        currTime += deltaT;
+    }
+
+
+    // Write file series
+
+    if (geomSeries.size())
+    {
+        Info<< nl << "write geom.vtp.series" << nl;
+        geomSeries.write("geom.vtp");
+    }
+    if (stateSeries.size())
+    {
+        Info<< nl << "write state.vtp.series" << nl;
+        stateSeries.write("state.vtp");
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

tutorials/incompressible/lumpedPointMotion/building/files/Allrun.moveMesh

@@ -0,0 +1,34 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+# Cleanup old junk that may prevent things from starting
+rm -f comms/OpenFOAM.lock
+
+# If OpenFOAM stops prematurely, trigger the external solver to stop
+trap '[ -e comms/OpenFOAM.lock ] && echo "status=done" > comms/OpenFOAM.lock' EXIT TERM INT
+
+# Simulated external solver.
+# Using -scale 0.01 since input movement table included visual scaling
+
+if false
+then
+    # Create response file
+    runApplication -overwrite \
+        ../code/building-motion -deltaT 0.001 -nTimes 5001 -output response.txt
+
+    # Use response file for states
+    runApplication -overwrite \
+        lumpedPointMovement -scale 0.01 -removeLock -slave response.txt &
+else
+
+    # Generate states on demand
+    runApplication -overwrite \
+        ../code/building-motion -scale 0.01 -removeLock -slave &
+fi
+
+# Run moveMesh with deltaT corresponding to dynamicMeshDict updateInterval
+runParallel moveMesh -deltaT 0.001
+
+#------------------------------------------------------------------------------

tutorials/incompressible/lumpedPointMotion/building/files/Allrun.movement

@@ -12,7 +12,22 @@ runApplication reconstructParMesh -constant -withZero -time 0
 runApplication lumpedPointZones
 
 # Simulated external solver
-# Using -scale=1 to see the excessively large movements
+# Use -scale=1 to see the excessively large movements
-runApplication lumpedPointMovement -span 25 -scale 1 ../files/response.txt
+
+if false
+then
+    # Create response file
+    runApplication -overwrite \
+        ../code/building-motion -deltaT 0.001 -nTimes 5001 -output response.txt
+
+    # Use response file for states
+    runApplication -overwrite \
+        lumpedPointMovement -span 25 -scale 1 response.txt
+else
+
+    # Generate states on demand
+    runApplication -overwrite \
+        ../code/building-motion -scale 1 -deltaT 0.025 -nTimes 201
+fi
 
 #------------------------------------------------------------------------------

tutorials/incompressible/lumpedPointMotion/building/files/Allrun.transient

@@ -10,9 +10,23 @@ rm -f comms/OpenFOAM.lock
 trap '[ -e comms/OpenFOAM.lock ] && echo "status=done" > comms/OpenFOAM.lock' EXIT TERM INT
 
 # Simulated external solver.
-# Using -scale since the input movement table is excessively large
+# Using -scale 0.01 since input movement table included visual scaling
-runApplication -overwrite \
+
-    lumpedPointMovement -scale 0.01 -removeLock -slave ../files/response.txt &
+if false
+then
+    # Create response file
+    runApplication -overwrite \
+        ../code/building-motion -deltaT 0.001 -nTimes 5001 -output response.txt
+
+    # Use response file for states
+    runApplication -overwrite \
+        lumpedPointMovement -scale 0.01 -removeLock -slave response.txt &
+else
+
+    # Generate states on demand
+    runApplication -overwrite \
+        ../code/building-motion -scale 0.01 -removeLock -slave &
+fi
 
 # Run OpenFOAM
 runParallel $(getApplication)

tutorials/incompressible/lumpedPointMotion/building/files/RunFunctions

@@ -0,0 +1,147 @@
+#---------------------------------*- sh -*-------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | www.openfoam.com
+#    \\/     M anipulation  |
+#------------------------------------------------------------------------------
+#     Copyright (C) 2020 OpenCFD Ltd.
+#------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+#
+# Script
+#     RunFunctions
+#
+# Description
+#     Additional functions for copy/linking FSI cases
+#
+#------------------------------------------------------------------------------
+
+#
+# copyParallelPointDisplacement caseDir timeName
+#
+# Copy pointDisplacement from caseDir/0/ to caseDir/timeName/
+#
+copyParallelPointDisplacement()
+{
+    local src="$1"
+    local dstTime="$2"
+    local file=pointDisplacement
+
+    [ -d "$src" ] || {
+        echo "Error: no directory: $src"
+        return 1
+    }
+
+    # Copy select directories
+    echo "    copy processor '$file' from 0/ -> $dstTime"
+    if [ -n "$dstTime" ]
+    then
+    (
+        cd "$src" || exit
+
+        for proc in processor*
+        do
+            if [ -d "$proc/0" ] && [ -d "$proc/$dstTime" ]
+            then
+                cp "$proc/0/$file" "$proc/$dstTime/$file"
+                if [ -d "$proc/0/include" ]
+                then
+                    cp -r "$proc/0/include" "$proc/$dstTime"
+                fi
+            fi
+        done
+    )
+    else
+        echo "    no destination time"
+    fi
+
+    # Restart from latestTime
+    foamDictionary "$src"/system/controlDict \
+        -entry startFrom -set latestTime
+
+    deltaT=$(foamDictionary "$src"/system/controlDict -entry deltaT -value)
+    latestTime=$(foamListTimes -case "$src" -noZero -latestTime -processor)
+
+    # Restart using steady results as first deltaT interval
+    echo "deltaT=$deltaT latestTime=$latestTime"
+    if [ -n "$latestTime" ] && [ "$deltaT" != "$latestTime" ]
+    then
+    (
+        cd "$src" || exit
+
+        for proc in processor*
+        do
+            if [ -d "$proc/$latestTime" ] && [ ! -d "$proc/$deltaT" ]
+            then
+                mv "$proc/$latestTime" "$proc/$deltaT"
+                rm -rf "$proc/$deltaT/uniform"
+            fi
+        done
+    )
+    fi
+
+    return 0
+}
+
+#
+# linkParallelCase srcDir dstDir
+#
+linkParallelCase()
+{
+    local src="$1"
+    local dst="$2"
+    shift 2
+
+    if [ -e "$dst" ]
+    then
+        echo "Case already linked: remove case directory $dst prior to linking"
+        return 1
+    elif [ ! -d "$src" ]
+    then
+        echo "Error: no directory to link: $src"
+        return 1
+    fi
+
+    echo "Linking $dst parallel case from $src"
+    mkdir -p "$dst"
+
+    # Copy system - may wish to change things
+    for i in system 0
+    do
+        echo "    copy $i/"
+        ( cd "$dst" && cp -r "../$src/$i" . )
+    done
+
+    echo "    link constant/"
+    ( cd "$dst" && ln -sf "../$src/constant" . )
+
+    echo "    link processor*/ with $# times: $@"
+    for proc in $(cd "$src" && \ls -d processor*)
+    do
+        ( cd "$dst" && ln -sf "../$src/$proc" . )
+    done
+
+    return 0
+}
+
+#
+# linkFiles srcDir dstDir
+#
+linkFiles()
+{
+    local src="$1"
+    local dst="$2"
+    shift
+
+    echo "Linking $dst control files from $src"
+    mkdir -p "$dst"
+
+    ( cd "$dst" && ln -sf ../"$src"/* . )
+
+    return 0
+}
+
+
+#------------------------------------------------------------------------------

tutorials/incompressible/lumpedPointMotion/building/steady/0.orig/pointDisplacement

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v1912                                 |
+|  \\    /   O peration     | Version:  v2006                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -67,6 +67,7 @@ boundaryField
     {
         type            lumpedPointDisplacement;
         value           uniform (0 0 0);
+        controllers     ( vertical );
     }
 }
 

tutorials/incompressible/lumpedPointMotion/building/steady/Allclean

@@ -8,4 +8,7 @@ cleanCase0
 # Clean up copied/derived files
 rm -rf constant/triSurface
 
+# Clean up debug/setup files
+rm -f *.txt *.vtp *.vtp.series
+
 #------------------------------------------------------------------------------

tutorials/incompressible/lumpedPointMotion/building/steady/Allrun

@@ -21,5 +21,4 @@ else
 
 fi
 
-
 #------------------------------------------------------------------------------

tutorials/incompressible/lumpedPointMotion/building/steady/Allrun.init

@@ -0,0 +1,11 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+# Copy geometry from resources directory
+
+mkdir -p constant/triSurface/
+cp "$FOAM_TUTORIALS"/resources/geometry/building_wtc2.obj constant/triSurface/
+
+#------------------------------------------------------------------------------

tutorials/incompressible/lumpedPointMotion/building/steady/Allrun.pre

@@ -3,9 +3,8 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
-# Copy building from resources directory
+# Get geometry and other resources
-mkdir -p constant/triSurface/
+./Allrun.init
-cp $FOAM_TUTORIALS/resources/geometry/building_wtc2.obj constant/triSurface/
 
 # runApplication surfaceFeatureExtract
 runApplication blockMesh

tutorials/incompressible/lumpedPointMotion/building/steady/constant/dynamicMeshDict

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v1912                                 |
+|  \\    /   O peration     | Version:  v2006                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -13,9 +13,13 @@ FoamFile
     location    "system";
     object      dynamicMeshDict;
 }
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+// Less frequent mesh motion
+updateControl   runTime;
+
+updateInterval  0.001;
+
 dynamicFvMesh   dynamicMotionSolverFvMesh;
 
 motionSolverLibs (fvMotionSolvers);

tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointControllers

@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "locations and connectivity";
+    object      dictionary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Locations of the lumped points
+points
+(
+    (0  0  0.00)
+    (0  0  0.05)
+    (0  0  0.10)
+    (0  0  0.15)
+    (0  0  0.20)
+    (0  0  0.25)
+    (0  0  0.30)
+    (0  0  0.35)
+    (0  0  0.40)
+    (0  0  0.45)
+    (0  0  0.50)
+);
+
+//- Connectivity for motion controllers
+controllers
+{
+    vertical
+    {
+        pointLabels (0 1 2 3 4 5 6 7 8 9 10);
+    }
+}
+
+
+// ************************************************************************* //

tutorials/incompressible/lumpedPointMotion/building/steady/system/lumpedPointMovement

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v1912                                 |
+|  \\    /   O peration     | Version:  v2006                                 |
 |   \\  /    A nd           | Website:  www.openfoam.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -15,36 +15,24 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Reference axis for the locations
+#include "lumpedPointControllers"
-axis            (0 0 1);
 
-// Locations of the lumped points
+//- Offset/shift for lumped points
-locations       11(0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5);
+origin  (0.3 0.15 0);
-
-// Division for pressure forces (0-1)
-division        0.5;
-
-//- If present, the offset of patch points compared to the locations
-//  Otherwise determined from the bounding box
-// centre          (0 0 0);
-
-//- The interpolation scheme
-interpolationScheme linear;
 
 //- Relaxation/scaling factor when updating positions
-relax           1.0;
+relax   1.0;
-
 
 forces
 {
     //- The pressure name (default: p)
-    p           p;
+    p       p;
 
     //- Reference pressure [Pa] (default: 0)
-    pRef        0;
+    pRef    0;
 
     //- Reference density for incompressible calculations (default: 1)
-    rhoRef      1;
+    rhoRef  1;
 }
 
 
@@ -56,6 +44,9 @@ communication
 
     waitInterval    1;
 
+    // Coupling frequency (default: 1)
+    //calcFrequency   1;
+
     timeOut         100;
 
     initByExternal  false;
@@ -78,20 +69,20 @@ communication
     scaleInput
     {
         //- Length multiplier (to metres). Eg 0.001 for [mm] -> [m]
-        length      1;
+        length  1;
     }
 
     // Scaling applied to values written to 'outputName'
     scaleOutput
     {
         //- Length multiplier (from metres). Eg 1000 for [m] -> [mm]
-        length      1;
+        length  1;
 
         //- Force units multiplier (from Pa)
-        force       1;
+        force   1;
 
         //- Moment units multiplier (from N.m)
-        moment      1;
+        moment  1;
     }
 }