Merge branch 'feature-MPPIC-dynamicMesh' into 'develop'

ENH: MPPIC dynamic mesh

See merge request Development/openfoam!406

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C

@@ -42,8 +42,8 @@ Description
 #include "CorrectPhi.H"
 
 #ifdef MPPIC
-    #include "basicKinematicMPPICCloud.H"
+    #include "basicKinematicCloud.H"
-    #define basicKinematicTypeCloud basicKinematicMPPICCloud
+    #define basicKinematicTypeCloud basicKinematicCloud
 #else
     #include "basicKinematicCollidingCloud.H"
     #define basicKinematicTypeCloud basicKinematicCollidingCloud
@@ -111,7 +111,6 @@ int main(int argc, char *argv[])
         continuousPhaseTransport.correct();
         muc = rhoc*continuousPhaseTransport.nu();
 
-        Info<< "Evolving " << kinematicCloud.name() << endl;
         kinematicCloud.evolve();
 
         // Update continuous phase volume fraction field

applications/solvers/lagrangian/DPMFoam/DPMFoam.C

@@ -43,8 +43,8 @@ Description
 #include "pimpleControl.H"
 
 #ifdef MPPIC
-    #include "basicKinematicMPPICCloud.H"
+    #include "basicKinematicCloud.H"
-    #define basicKinematicTypeCloud basicKinematicMPPICCloud
+    #define basicKinematicTypeCloud basicKinematicCloud
 #else
     #include "basicKinematicCollidingCloud.H"
     #define basicKinematicTypeCloud basicKinematicCollidingCloud
@@ -118,6 +118,12 @@ int main(int argc, char *argv[])
         cloudVolSUSu.correctBoundaryConditions();
         cloudSU.source() = Zero;
 
+//         cloudVolSUSu.primitiveFieldRef() =
+//             (cloudSU.diag()*Uc() - cloudSU.source())/mesh.V();
+//         cloudVolSUSu.correctBoundaryConditions();
+//         cloudSU.source() = cloudSU.diag()*Uc();
+
+
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -7,6 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
     Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -309,7 +310,9 @@ Foam::polyMesh::polyMesh(const IOobject& io, const bool doInit)
     moving_(false),
     topoChanging_(false),
     curMotionTimeIndex_(time().timeIndex()),
-    oldPointsPtr_(nullptr)
+    oldPointsPtr_(nullptr),
+    oldCellCentresPtr_(nullptr),
+    storeOldCellCentres_(false)
 {
     if (!owner_.headerClassName().empty())
     {
@@ -513,7 +516,9 @@ Foam::polyMesh::polyMesh
     moving_(false),
     topoChanging_(false),
     curMotionTimeIndex_(time().timeIndex()),
-    oldPointsPtr_(nullptr)
+    oldPointsPtr_(nullptr),
+    oldCellCentresPtr_(nullptr),
+    storeOldCellCentres_(false)
 {
     // Note: changed that the constructors where values can be supplied
     //       (points, faces, owner/neighbour) use the readOpt. All others
@@ -669,7 +674,9 @@ Foam::polyMesh::polyMesh
     moving_(false),
     topoChanging_(false),
     curMotionTimeIndex_(time().timeIndex()),
-    oldPointsPtr_(nullptr)
+    oldPointsPtr_(nullptr),
+    oldCellCentresPtr_(nullptr),
+    storeOldCellCentres_(false)
 {
     // Note: probably needs io.readOpt() for points/faces/cells etc so
     //       we can run with READ_IF_PRESENT. See constructor above.
@@ -1120,6 +1127,11 @@ const Foam::labelList& Foam::polyMesh::faceNeighbour() const
 
 const Foam::pointField& Foam::polyMesh::oldPoints() const
 {
+    if (!moving_)
+    {
+        return points_;
+    }
+
     if (!oldPointsPtr_)
     {
         if (debug)
@@ -1135,6 +1147,24 @@ const Foam::pointField& Foam::polyMesh::oldPoints() const
 }
 
 
+const Foam::pointField& Foam::polyMesh::oldCellCentres() const
+{
+    storeOldCellCentres_ = true;
+
+    if (!moving_)
+    {
+        return cellCentres();
+    }
+
+    if (!oldCellCentresPtr_)
+    {
+        oldCellCentresPtr_.reset(new pointField(cellCentres()));
+    }
+
+    return *oldCellCentresPtr_;
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 (
     const pointField& newPoints
@@ -1166,6 +1196,12 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
                 << " index " << time().timeIndex() << endl;
         }
 
+        if (storeOldCellCentres_)
+        {
+            oldCellCentresPtr_.clear();
+            oldCellCentresPtr_.reset(new pointField(cellCentres()));
+        }
+
         // Mesh motion in the new time step
         oldPointsPtr_.clear();
         oldPointsPtr_.reset(new pointField(points_));
@@ -1261,6 +1297,7 @@ void Foam::polyMesh::resetMotion() const
 {
     curMotionTimeIndex_ = 0;
     oldPointsPtr_.clear();
+    oldCellCentresPtr_.clear();
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyMesh.H

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -80,10 +81,9 @@ class polyMesh
     public objectRegistry,
     public primitiveMesh
 {
-
 public:
 
-    // Public data types
+    // Public Data
 
         //- Enumeration defining the state of the mesh after a read update.
         //  Used for post-processing applications, where the mesh
@@ -114,7 +114,7 @@ public:
 
 private:
 
-    // Permanent data
+    // Private Data
 
         // Primitive mesh data
 
@@ -189,6 +189,12 @@ private:
             //- Old points (for the last mesh motion)
             mutable autoPtr<pointField> oldPointsPtr_;
 
+            //- Old cell centres (for the last mesh motion)
+            mutable autoPtr<pointField> oldCellCentresPtr_;
+
+            //- Whether or not to store the old cell centres
+            mutable bool storeOldCellCentres_;
+
 
     // Private Member Functions
 
@@ -432,6 +438,9 @@ public:
             //- Return old points for mesh motion
             virtual const pointField& oldPoints() const;
 
+            //- Return old points for mesh motion
+            virtual const pointField& oldCellCentres() const;
+
             //- Return boundary mesh
             const polyBoundaryMesh& boundaryMesh() const
             {

src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C

@@ -7,6 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
     Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -576,7 +577,9 @@ Foam::polyMesh::polyMesh
     moving_(false),
     topoChanging_(false),
     curMotionTimeIndex_(time().timeIndex()),
-    oldPointsPtr_(nullptr)
+    oldPointsPtr_(nullptr),
+    oldCellCentresPtr_(nullptr),
+    storeOldCellCentres_(false)
 {
     DebugInfo
         << "Constructing polyMesh from cell and boundary shapes." << endl;
@@ -856,7 +859,9 @@ Foam::polyMesh::polyMesh
     moving_(false),
     topoChanging_(false),
     curMotionTimeIndex_(time().timeIndex()),
-    oldPointsPtr_(nullptr)
+    oldPointsPtr_(nullptr),
+    oldCellCentresPtr_(nullptr),
+    storeOldCellCentres_(false)
 {
     DebugInfo
         << "Constructing polyMesh from cell and boundary shapes." << endl;

src/OpenFOAM/meshes/polyMesh/polyMeshUpdate.C

@@ -7,6 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
     Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -119,6 +120,30 @@ void Foam::polyMesh::updateMesh(const mapPolyMesh& mpm)
         }
     }
 
+    if (oldCellCentresPtr_)
+    {
+        // Make a copy of the original cell-centres
+        pointField oldMotionCellCentres = oldCellCentresPtr_();
+
+        pointField& newMotionCellCentres = oldCellCentresPtr_();
+
+        // Resize the list to new size
+        newMotionCellCentres.setSize(cellCentres().size());
+
+        // Map the list
+        newMotionCellCentres.map(oldMotionCellCentres, mpm.cellMap());
+
+        // Any points created out-of-nothing get set to the current coordinate
+        // for lack of anything better.
+        forAll(mpm.cellMap(), newCelli)
+        {
+            if (mpm.cellMap()[newCelli] == -1)
+            {
+                newMotionCellCentres[newCelli] = cellCentres()[newCelli];
+            }
+        }
+    }
+
     meshObject::updateMesh<polyMesh>(*this, mpm);
     meshObject::updateMesh<pointMesh>(*this, mpm);
 

src/lagrangian/basic/Cloud/Cloud.C

@@ -7,6 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
     Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -85,6 +86,7 @@ Foam::Cloud<ParticleType>::Cloud
     geometryType_(cloud::geometryType::COORDINATES)
 {
     checkPatches();
+    polyMesh_.oldCellCentres();
 
     // Ask for the tetBasePtIs to trigger all processors to build
     // them, otherwise, if some processors have no particles then
@@ -174,7 +176,7 @@ void Foam::Cloud<ParticleType>::move
     // Initialise the stepFraction moved for the particles
     forAllIters(*this, pIter)
     {
-        pIter().stepFraction() = 0;
+        pIter().reset();
     }
 
     // List of lists of particles to be transferred for all of the
@@ -360,6 +362,7 @@ void Foam::Cloud<ParticleType>::autoMap(const mapPolyMesh& mapper)
     // them, otherwise, if some processors have no particles then
     // there is a comms mismatch.
     polyMesh_.tetBasePtIs();
+    polyMesh_.oldCellCentres();
 
     const vectorField& positions = globalPositionsPtr_();
 

src/lagrangian/basic/Cloud/CloudIO.C

@@ -7,6 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
     Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -177,6 +178,9 @@ Foam::Cloud<ParticleType>::Cloud
 {
     checkPatches();
 
+    polyMesh_.tetBasePtIs();
+    polyMesh_.oldCellCentres();
+
     initCloud(checkClass);
 }
 

src/lagrangian/basic/particle/particle.C

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,7 +41,7 @@ namespace Foam
     defineTypeNameAndDebug(particle, 0);
 }
 
-const Foam::scalar Foam::particle::negativeSpaceDisplacementFactor = 1.01;
+const Foam::label Foam::particle::maxNBehind_ = 10;
 
 Foam::label Foam::particle::particleCount_ = 0;
 
@@ -365,7 +366,6 @@ void Foam::particle::changeCell()
     const label ownerCellI = mesh_.faceOwner()[tetFacei_];
     const bool isOwner = celli_ == ownerCellI;
     celli_ = isOwner ? mesh_.faceNeighbour()[tetFacei_] : ownerCellI;
-
     // Reflect to account for the change of triangle orientation in the new cell
     reflect();
 }
@@ -414,78 +414,84 @@ void Foam::particle::locate
     const string& boundaryMsg
 )
 {
+    if (debug)
+    {
+        Info << "Particle " << origId() << nl << FUNCTION_NAME << nl << endl;
+    }
+
     celli_ = celli;
 
     // Find the cell, if it has not been given
-    if (celli_ < 0)
+    if (celli < 0)
     {
         celli_ = mesh_.cellTree().findInside(position);
     }
-    if (celli_ < 0)
+    if (celli < 0)
     {
         FatalErrorInFunction
             << "Cell not found for particle position " << position << "."
             << exit(FatalError);
     }
 
-    // Put the particle at (almost) the cell centre and in a random tet.
+    const vector displacement = position - mesh_.cellCentres()[celli_];
-    // Note perturbing the cell centre to make sure we find at least one
+
-    // tet containing it. With start point exactly at the cell centre very
+    // Loop all cell tets to find the one containing the position. Track
-    // occasionally it would not get found in any of the tets
+    // through each tet from the cell centre. If a tet contains the position
-    coordinates_ = barycentric(1-3*SMALL, SMALL, SMALL, SMALL);
+    // then the track will end with a single trackToTri.
-    tetFacei_ = mesh_.cells()[celli_][0];
+    const class cell& c = mesh_.cells()[celli_];
-    tetPti_ = 1;
+    scalar minF = VGREAT;
+    label minTetFacei = -1, minTetPti = -1;
+    forAll(c, cellTetFacei)
+    {
+        const class face& f = mesh_.faces()[c[cellTetFacei]];
+        for (label tetPti = 1; tetPti < f.size() - 1; ++tetPti)
+        {
+            coordinates_ = barycentric(1, 0, 0, 0);
+            tetFacei_ = c[cellTetFacei];
+            tetPti_ = tetPti;
             facei_ = -1;
 
-    // Track to the injection point
+            label tetTriI = -1;
-    track(position - this->position(), 0);
+            const scalar f = trackToTri(displacement, 0, tetTriI);
 
+            if (tetTriI == -1)
+            {
+                return;
+            }
+
+            if (f < minF)
+            {
+                minF = f;
+                minTetFacei = tetFacei_;
+                minTetPti = tetPti_;
+            }
+        }
+    }
+
+    // The particle must be (hopefully only slightly) outside the cell. Track
+    // into the tet which got the furthest.
+    coordinates_ = barycentric(1, 0, 0, 0);
+    tetFacei_ = minTetFacei;
+    tetPti_ = minTetPti;
+    facei_ = -1;
+    track(displacement, 0);
     if (!onFace())
     {
         return;
     }
 
-    // We hit a boundary ...
+    // If we are here then we hit a boundary
     if (boundaryFail)
     {
-        FatalErrorInFunction << boundaryMsg
+        FatalErrorInFunction << boundaryMsg << exit(FatalError);
-            << " when tracking from centre " << mesh_.cellCentres()[celli_]
-            << " of cell " << celli_ << " to position " << position
-            << exit(FatalError);
     }
     else
     {
-        // Re-do the track, but this time do the bit tangential to the
-        // direction/patch first. This gets us as close as possible to the
-        // original path/position.
-
-        if (direction == nullptr)
-        {
-            const polyPatch& p = mesh_.boundaryMesh()[patch()];
-            direction = &p.faceNormals()[p.whichFace(facei_)];
-        }
-
-        const vector n = *direction/mag(*direction);
-        const vector s = position - mesh_.cellCentres()[celli_];
-        const vector sN = (s & n)*n;
-        const vector sT = s - sN;
-
-        coordinates_ = barycentric(1, 0, 0, 0);
-        tetFacei_ = mesh_.cells()[celli_][0];
-        tetPti_ = 1;
-        facei_ = -1;
-
-        track(sT, 0);
-        track(sN, 0);
-
         static label nWarnings = 0;
         static const label maxNWarnings = 100;
-        if (nWarnings < maxNWarnings)
+        if ((nWarnings < maxNWarnings) && boundaryFail)
         {
-            WarningInFunction << boundaryMsg.c_str()
+            WarningInFunction << boundaryMsg.c_str() << endl;
-                << " when tracking from centre " << mesh_.cellCentres()[celli_]
-                << " of cell " << celli_ << " to position " << position
-                << endl;
             ++ nWarnings;
         }
         if (nWarnings == maxNWarnings)
@@ -496,6 +502,7 @@ void Foam::particle::locate
             ++ nWarnings;
         }
     }
+
 }
 
 
@@ -516,7 +523,9 @@ Foam::particle::particle
     tetFacei_(tetFacei),
     tetPti_(tetPti),
     facei_(-1),
-    stepFraction_(0.0),
+    stepFraction_(1.0),
+    behind_(0.0),
+    nBehind_(0),
     origProc_(Pstream::myProcNo()),
     origId_(getNewParticleID())
 {}
@@ -535,7 +544,9 @@ Foam::particle::particle
     tetFacei_(-1),
     tetPti_(-1),
     facei_(-1),
-    stepFraction_(0.0),
+    stepFraction_(1.0),
+    behind_(0.0),
+    nBehind_(0),
     origProc_(Pstream::myProcNo()),
     origId_(getNewParticleID())
 {
@@ -566,7 +577,9 @@ Foam::particle::particle
     tetFacei_(tetFacei),
     tetPti_(tetPti),
     facei_(-1),
-    stepFraction_(0.0),
+    stepFraction_(1.0),
+    behind_(0.0),
+    nBehind_(0),
     origProc_(Pstream::myProcNo()),
     origId_(getNewParticleID())
 {
@@ -593,6 +606,8 @@ Foam::particle::particle(const particle& p)
     tetPti_(p.tetPti_),
     facei_(p.facei_),
     stepFraction_(p.stepFraction_),
+    behind_(p.behind_),
+    nBehind_(p.nBehind_),
     origProc_(p.origProc_),
     origId_(p.origId_)
 {}
@@ -607,6 +622,8 @@ Foam::particle::particle(const particle& p, const polyMesh& mesh)
     tetPti_(p.tetPti_),
     facei_(p.facei_),
     stepFraction_(p.stepFraction_),
+    behind_(p.behind_),
+    nBehind_(p.nBehind_),
     origProc_(p.origProc_),
     origId_(p.origId_)
 {}
@@ -645,7 +662,8 @@ Foam::scalar Foam::particle::trackToFace
 
     facei_ = -1;
 
-    while (true)
+    // Loop the tets in the current cell
+    while (nBehind_ < maxNBehind_)
     {
         f *= trackToTri(f*displacement, f*fraction, tetTriI);
 
@@ -666,6 +684,25 @@ Foam::scalar Foam::particle::trackToFace
             changeTet(tetTriI);
         }
     }
+
+    // Warn if stuck, and incorrectly advance the step fraction to completion
+    static label stuckID = -1, stuckProc = -1;
+    if (origId_ != stuckID && origProc_ != stuckProc)
+    {
+        WarningInFunction
+            << "Particle #" << origId_ << " got stuck at " << position()
+            << endl;
+    }
+
+    stuckID = origId_;
+    stuckProc = origProc_;
+
+    stepFraction_ += f*fraction;
+
+    behind_ = 0;
+    nBehind_ = 0;
+
+    return 0;
 }
 
 
@@ -702,11 +739,8 @@ Foam::scalar Foam::particle::trackToStationaryTri
             << "Start local coordinates = " << y0 << endl;
     }
 
-    // Get the factor by which the displacement is increased
-    const scalar f = detA >= 0 ? 1 : negativeSpaceDisplacementFactor;
-
     // Calculate the local tracking displacement
-    barycentric Tx1(f*x1 & T);
+    barycentric Tx1(x1 & T);
 
     if (debug)
     {
@@ -718,7 +752,7 @@ Foam::scalar Foam::particle::trackToStationaryTri
     scalar muH = std::isnormal(detA) && detA <= 0 ? VGREAT : 1/detA;
     for (label i = 0; i < 4; ++ i)
     {
-        if (std::isnormal(Tx1[i]) && Tx1[i] < 0)
+        if (Tx1[i] < - detA*SMALL)
         {
             scalar mu = - y0[i]/Tx1[i];
 
@@ -776,6 +810,30 @@ Foam::scalar Foam::particle::trackToStationaryTri
     // Set the proportion of the track that has been completed
     stepFraction_ += fraction*muH*detA;
 
+    if (debug)
+    {
+        Info << "Step Fraction : " << stepFraction_ << endl;
+        Info << "fraction*muH*detA : " << fraction*muH*detA << endl;
+        Info << "static muH : " << muH << endl;
+        Info << "origId() : " << origId() << endl;
+    }
+
+    // Accumulate displacement behind
+    if (detA <= 0 || nBehind_ > 0)
+    {
+        behind_ += muH*detA*mag(displacement);
+
+        if (behind_ > 0)
+        {
+            behind_ = 0;
+            nBehind_ = 0;
+        }
+        else
+        {
+            ++ nBehind_;
+        }
+    }
+
     return iH != -1 ? 1 - muH*detA : 0;
 }
 
@@ -803,12 +861,20 @@ Foam::scalar Foam::particle::trackToMovingTri
     FixedList<barycentricTensor, 3> T;
     movingTetReverseTransform(fraction, centre, detA, T);
 
-    // Get the factor by which the displacement is increased
+    if (debug)
-    const scalar f = detA[0] >= 0 ? 1 : negativeSpaceDisplacementFactor;
+    {
+        Pair<vector> o, b, v1, v2;
+        movingTetGeometry(fraction, o, b, v1, v2);
+        Info<< "Tet points o=" << o[0] << ", b=" << b[0]
+            << ", v1=" << v1[0] << ", v2=" << v2[0] << endl
+            << "Tet determinant = " << detA[0] << endl
+            << "Start local coordinates = " << y0[0] << endl;
+    }
+
 
     // Get the relative global position
     const vector x0Rel = x0 - centre[0];
-    const vector x1Rel = f*x1 - centre[1];
+    const vector x1Rel = x1 - centre[1];
 
     // Form the determinant and hit equations
     cubicEqn detAEqn(sqr(detA[0])*detA[3], detA[0]*detA[2], detA[1], 1);
@@ -826,6 +892,16 @@ Foam::scalar Foam::particle::trackToMovingTri
         hitEqn[i] = cubicEqn(hitEqnA[i], hitEqnB[i], hitEqnC[i], hitEqnD[i]);
     }
 
+    if (debug)
+    {
+        for (label i = 0; i < 4; ++ i)
+        {
+            Info<< (i ? "             " : "Hit equation ") << i << " = "
+                << hitEqn[i] << endl;
+        }
+        Info<< " DetA equation = " << detA << endl;
+    }
+
     // Calculate the hit fraction
     label iH = -1;
     scalar muH = std::isnormal(detA[0]) && detA[0] <= 0 ? VGREAT : 1/detA[0];
@@ -835,8 +911,33 @@ Foam::scalar Foam::particle::trackToMovingTri
 
         for (label j = 0; j < 3; ++j)
         {
-            if (mu.type(j) == roots::real && hitEqn[i].derivative(mu[j]) < 0)
+            if
+            (
+                mu.type(j) == roots::real
+             && hitEqn[i].derivative(mu[j]) < - detA[0]*SMALL
+            )
             {
+                if (debug)
+                {
+                    const barycentric yH
+                    (
+                        hitEqn[0].value(mu[j]),
+                        hitEqn[1].value(mu[j]),
+                        hitEqn[2].value(mu[j]),
+                        hitEqn[3].value(mu[j])
+                    );
+                    const scalar detAH = detAEqn.value(mu[j]);
+
+                    Info<< "Hit on tet face " << i << " at local coordinate "
+                        << (std::isnormal(detAH) ? name(yH/detAH) : "???")
+                        << ", " << mu[j]*detA[0]*100 << "% of the "
+                        << "way along the track" << endl;
+
+                    Info<< "derivative : " << hitEqn[i].derivative(mu[j]) << nl
+                        << " coord " << j << " mu[j]: " << mu[j] << nl
+                        << " hitEq " << i << endl;
+                }
+
                 if (0 <= mu[j] && mu[j] < muH)
                 {
                     iH = i;
@@ -887,8 +988,26 @@ Foam::scalar Foam::particle::trackToMovingTri
     coordinates_ = yH;
     tetTriI = iH;
 
+    scalar advance = muH*detA[0];
+
     // Set the proportion of the track that has been completed
-    stepFraction_ += fraction*muH*detA[0];
+    stepFraction_ += fraction*advance;
+
+    // Accumulate displacement behind
+    if (detA[0] <= 0 || nBehind_ > 0)
+    {
+        behind_ += muH*detA[0]*mag(displacement);
+
+        if (behind_ > 0)
+        {
+            behind_ = 0;
+            nBehind_ = 0;
+        }
+        else
+        {
+            ++ nBehind_;
+        }
+    }
 
     if (debug)
     {
@@ -900,10 +1019,15 @@ Foam::scalar Foam::particle::trackToMovingTri
         {
             Pout<< "Track hit no tet faces" << endl;
         }
-        Pout<< "End local coordinates = " << yH << endl
+//         Pout<< "End local coordinates = " << yH << endl
-            << "End global coordinates = " << position() << endl;
+//             << "End global coordinates = " << position() << endl
+//             << "Tracking displacement = " << position() - x0 << endl
+//             << muH*detA[0]*100 << "% of the step from " << stepFraction_
+//             << " to " << stepFraction_ + fraction << " completed" << endl
+//             << endl;
     }
 
+
     return iH != -1 ? 1 - muH*detA[0] : 0;
 }
 
@@ -915,7 +1039,7 @@ Foam::scalar Foam::particle::trackToTri
     label& tetTriI
 )
 {
-    if (mesh_.moving())
+    if ((mesh_.moving() && (stepFraction_ != 1 || fraction != 0)))
     {
         return trackToMovingTri(displacement, fraction, tetTriI);
     }
@@ -1049,7 +1173,7 @@ void Foam::particle::correctAfterInteractionListReferral(const label celli)
     // so this approximate topology is good enough. By using the nearby cell we
     // minimise the error associated with the incorrect topology.
     coordinates_ = barycentric(1, 0, 0, 0);
-    if (mesh_.moving())
+    if (mesh_.moving() && stepFraction_ != 1)
     {
         Pair<vector> centre;
         FixedList<scalar, 4> detA;

src/lagrangian/basic/particle/particle.H

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -78,7 +79,7 @@ class particle
 :
     public IDLList<particle>::link
 {
-    // Private member data
+    // Private Data
 
         //- Size in bytes of the position data
         static const std::size_t sizeofPosition;
@@ -86,15 +87,9 @@ class particle
         //- Size in bytes of the fields
         static const std::size_t sizeofFields;
 
-        //- The factor by which the displacement is increased when passing
+        //- The value of nBehind_ at which tracking is abandoned. See the
-        //  through negative space. This should be slightly bigger than one.
+        //  description of nBehind_.
-        //  This is needed as a straight trajectory can form a closed loop
+        static const label maxNBehind_;
-        //  through regions of overlapping positive and negative space, leading
-        //  to a track which never ends. By increasing the rate of displacement
-        //  through negative regions, the change in track fraction over this
-        //  part of the loop no longer exactly cancels the change over the
-        //  positive part, and the track therefore terminates.
-        static const scalar negativeSpaceDisplacementFactor;
 
 
 public:
@@ -103,7 +98,7 @@ public:
     {
     public:
 
-        // Public data
+        // Public Data
 
             //- Flag to switch processor
             bool switchProcessor;
@@ -135,7 +130,7 @@ public:
 
 private:
 
-    // Private data
+    // Private Data
 
         //- Reference to the polyMesh database
         const polyMesh& mesh_;
@@ -147,12 +142,12 @@ private:
         label celli_;
 
         //- Index of the face that owns the decomposed tet that the
-        //  particle is in
+        //- particle is in
         label tetFacei_;
 
         //- Index of the point on the face that defines the decomposed
-        //  tet that the particle is in.  Relative to the face base
+        //- tet that the particle is in.  Relative to the face base
-        //  point.
+        //- point.
         label tetPti_;
 
         //- Face index if the particle is on a face otherwise -1
@@ -161,6 +156,18 @@ private:
         //- Fraction of time-step completed
         scalar stepFraction_;
 
+        //- The distance behind the maximum distance reached so far
+        scalar behind_;
+
+        //- The number of tracks carried out that ended in a distance behind the
+        //- maximum distance reached so far. Once this reaches maxNBehind_,
+        //  tracking is abandoned for the current step. This is needed because
+        //  when tetrahedra are inverted a straight trajectory can form a closed
+        //  loop through regions of overlapping positive and negative space.
+        //  Without this break clause, such loops can result in a valid track
+        //  which never ends.
+        label nBehind_;
+
         //- Originating processor id
         label origProc_;
 
@@ -331,7 +338,7 @@ protected:
 
 public:
 
-    // Static data members
+    // Static Data Members
 
         //- Runtime type information
         TypeName("particle");
@@ -526,6 +533,9 @@ public:
             //- Return current particle position
             inline vector position() const;
 
+            //- Reset particle data
+            inline void reset();
+
 
     // Track
 
@@ -664,6 +674,7 @@ public:
         void relocate(const point& position, const label celli = -1);
 
 
+
     // I-O
 
         //- Write the name representation to stream

src/lagrangian/basic/particle/particleI.H

@@ -6,6 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2018 OpenFOAM Foundation
+    Copyright (C) 2011-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -68,6 +70,7 @@ inline void Foam::particle::movingTetGeometry
 ) const
 {
     const triFace triIs(currentTetIndices().faceTriIs(mesh_));
+
     const pointField& ptsOld = mesh_.oldPoints();
     const pointField& ptsNew = mesh_.points();
 
@@ -75,8 +78,11 @@ inline void Foam::particle::movingTetGeometry
     // we need to put a mesh_.oldCellCentres() method in for this to work. The
     // values obtained from the mesh and those obtained from the cell do not
     // necessarily match. See mantis #1993.
-    const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
+    //const vector ccOld = mesh_.cells()[celli_].centre(ptsOld, mesh_.faces());
-    const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
+    //const vector ccNew = mesh_.cells()[celli_].centre(ptsNew, mesh_.faces());
+
+    const vector ccOld = mesh_.oldCellCentres()[celli_];
+    const vector ccNew = mesh_.cellCentres()[celli_];
 
     // Old and new points and cell centres are not sub-cycled. If we are sub-
     // cycling, then we have to account for the timestep change here by
@@ -265,7 +271,7 @@ inline Foam::tetIndices Foam::particle::currentTetIndices() const
 
 inline Foam::barycentricTensor Foam::particle::currentTetTransform() const
 {
-    if (mesh_.moving())
+    if (mesh_.moving() && stepFraction_ != 1)
     {
         return movingTetTransform(0)[0];
     }
@@ -312,6 +318,14 @@ inline Foam::vector Foam::particle::position() const
 }
 
 
+inline void Foam::particle::reset()
+{
+    stepFraction_ = 0;
+    nBehind_ = 0;
+    behind_ = 0;
+}
+
+
 void Foam::particle::patchData(vector& n, vector& U) const
 {
     if (!onBoundaryFace())
@@ -321,7 +335,7 @@ void Foam::particle::patchData(vector& n, vector& U) const
             << exit(FatalError);
     }
 
-    if (mesh_.moving())
+    if ((mesh_.moving() && stepFraction_ != 1))
     {
         Pair<vector> centre, base, vertex1, vertex2;
         movingTetGeometry(1, centre, base, vertex1, vertex2);

src/lagrangian/basic/particle/particleTemplates.C

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2020 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -293,7 +294,7 @@ void Foam::particle::writeObjects(const CloudType& c, objectRegistry& obr)
 template<class TrackCloudType>
 void Foam::particle::hitFace
 (
-    const vector& direction,
+    const vector& displacement,
     TrackCloudType& cloud,
     trackingData& td
 )
@@ -337,11 +338,11 @@ void Foam::particle::hitFace
             }
             else if (isA<cyclicACMIPolyPatch>(patch))
             {
-                p.hitCyclicACMIPatch(cloud, ttd, direction);
+                p.hitCyclicACMIPatch(cloud, ttd, displacement);
             }
             else if (isA<cyclicAMIPolyPatch>(patch))
             {
-                p.hitCyclicAMIPatch(cloud, ttd, direction);
+                p.hitCyclicAMIPatch(cloud, ttd, displacement);
             }
             else if (isA<processorPolyPatch>(patch))
             {
@@ -459,7 +460,7 @@ void Foam::particle::hitCyclicAMIPatch
 (
     TrackCloudType&,
     trackingData& td,
-    const vector& direction
+    const vector& displacement
 )
 {
     vector pos = position();
@@ -468,7 +469,7 @@ void Foam::particle::hitCyclicAMIPatch
         static_cast<const cyclicAMIPolyPatch&>(mesh_.boundaryMesh()[patch()]);
     const cyclicAMIPolyPatch& receiveCpp = cpp.neighbPatch();
     const label sendFacei = cpp.whichFace(facei_);
-    const label receiveFacei = cpp.pointFace(sendFacei, direction, pos);
+    const label receiveFacei = cpp.pointFace(sendFacei, displacement, pos);
 
     if (receiveFacei < 0)
     {
@@ -485,12 +486,19 @@ void Foam::particle::hitCyclicAMIPatch
     facei_ = tetFacei_ = receiveFacei + receiveCpp.start();
 
     // Locate the particle on the receiving side
-    vector directionT = direction;
+    vector displacementT = displacement;
-    cpp.reverseTransformDirection(directionT, sendFacei);
+    cpp.reverseTransformDirection(displacementT, sendFacei);
+
+    // NOTE: The ray used to find the hit location accross the AMI might not
+    // be consistent in the displacement direction. Therefore a particle can
+    // be looping accross AMI patches indefinitely. Advancing the particle
+    // trajectory inside the cell is a possible solution.
+    const vector dispDir = cpp.fraction()*displacementT;
+    stepFraction_ += cpp.fraction();
     locate
     (
-        pos,
+        pos + dispDir,
-        &directionT,
+        &displacementT,
         mesh_.faceOwner()[facei_],
         false,
         "Particle crossed between " + cyclicAMIPolyPatch::typeName +
@@ -523,6 +531,26 @@ void Foam::particle::hitCyclicAMIPatch
         );
         transformProperties(-s);
     }
+
+    //if (onBoundaryFace())
+    {
+//         vector receiveNormal, receiveDisplacement;
+//         patchData(receiveNormal, receiveDisplacement);
+//
+//         if (((displacementT - fraction*receiveDisplacement)&receiveNormal) > 0)
+//         {
+//             td.keepParticle = false;
+//             WarningInFunction
+//                 << "Particle transfer from " << cyclicAMIPolyPatch::typeName
+//                 << " patches " << cpp.name() << " to " << receiveCpp.name()
+//                 << " failed at position " << pos << " and with displacement "
+//                 << (displacementT - fraction*receiveDisplacement) << nl
+//                 << "    The displacement points into both the source and "
+//                 << "receiving faces, so the tracking cannot proceed" << nl
+//                 << "    The particle has been removed" << nl << endl;
+//             return;
+//         }
+    }
 }
 
 
@@ -531,7 +559,7 @@ void Foam::particle::hitCyclicACMIPatch
 (
     TrackCloudType& cloud,
     trackingData& td,
-    const vector& direction
+    const vector& displacement
 )
 {
     const cyclicACMIPolyPatch& cpp =
@@ -551,20 +579,20 @@ void Foam::particle::hitCyclicACMIPatch
     if (!couple && !nonOverlap)
     {
         vector pos = position();
-        couple = cpp.pointFace(localFacei, direction, pos) >= 0;
+        couple = cpp.pointFace(localFacei, displacement, pos) >= 0;
         nonOverlap = !couple;
     }
 
     if (couple)
     {
-        hitCyclicAMIPatch(cloud, td, direction);
+        hitCyclicAMIPatch(cloud, td, displacement);
     }
     else
     {
         // Move to the face associated with the non-overlap patch and redo the
         // face interaction.
         tetFacei_ = facei_ = cpp.nonOverlapPatch().start() + localFacei;
-        hitFace(direction, cloud, td);
+        hitFace(displacement, cloud, td);
     }
 }
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -38,6 +38,12 @@ License
 #include "SurfaceFilmModel.H"
 #include "profiling.H"
 
+#include "PackingModel.H"
+#include "ParticleStressModel.H"
+#include "DampingModel.H"
+#include "IsotropyModel.H"
+#include "TimeScaleModel.H"
+
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
 template<class CloudType>
@@ -79,6 +85,33 @@ void Foam::KinematicCloud<CloudType>::setModels()
         ).ptr()
     );
 
+    packingModel_.reset
+    (
+        PackingModel<KinematicCloud<CloudType>>::New
+        (
+            subModelProperties_,
+            *this
+        ).ptr()
+    );
+
+    dampingModel_.reset
+    (
+        DampingModel<KinematicCloud<CloudType>>::New
+        (
+            subModelProperties_,
+            *this
+        ).ptr()
+    );
+
+    isotropyModel_.reset
+    (
+        IsotropyModel<KinematicCloud<CloudType>>::New
+        (
+            subModelProperties_,
+            *this
+        ).ptr()
+    );
+
     UIntegrator_.reset
     (
         integrationScheme::New
@@ -210,7 +243,6 @@ void Foam::KinematicCloud<CloudType>::evolveCloud
 
         injectors_.inject(cloud, td);
 
-
         // Assume that motion will update the cellOccupancy as necessary
         // before it is required.
         cloud.motion(cloud, td);
@@ -264,6 +296,15 @@ void Foam::KinematicCloud<CloudType>::postEvolve
             true
         );
     }
+
+    if (this->dampingModel().active())
+    {
+        this->dampingModel().cacheFields(false);
+    }
+    if (this->packingModel().active())
+    {
+        this->packingModel().cacheFields(false);
+    }
 }
 
 
@@ -285,6 +326,10 @@ void Foam::KinematicCloud<CloudType>::cloudReset(KinematicCloud<CloudType>& c)
     stochasticCollisionModel_.reset(c.stochasticCollisionModel_.ptr());
     surfaceFilmModel_.reset(c.surfaceFilmModel_.ptr());
 
+    packingModel_.reset(c.packingModel_.ptr());
+    dampingModel_.reset(c.dampingModel_.ptr());
+    isotropyModel_.reset(c.isotropyModel_.ptr());
+
     UIntegrator_.reset(c.UIntegrator_.ptr());
 }
 
@@ -375,6 +420,11 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     patchInteractionModel_(nullptr),
     stochasticCollisionModel_(nullptr),
     surfaceFilmModel_(nullptr),
+
+    packingModel_(nullptr),
+    dampingModel_(nullptr),
+    isotropyModel_(nullptr),
+
     UIntegrator_(nullptr),
     UTrans_
     (
@@ -458,6 +508,11 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     patchInteractionModel_(c.patchInteractionModel_->clone()),
     stochasticCollisionModel_(c.stochasticCollisionModel_->clone()),
     surfaceFilmModel_(c.surfaceFilmModel_->clone()),
+
+    packingModel_(c.packingModel_->clone()),
+    dampingModel_(c.dampingModel_->clone()),
+    isotropyModel_(c.isotropyModel_->clone()),
+
     UIntegrator_(c.UIntegrator_->clone()),
     UTrans_
     (
@@ -549,6 +604,11 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     patchInteractionModel_(nullptr),
     stochasticCollisionModel_(nullptr),
     surfaceFilmModel_(nullptr),
+
+    packingModel_(nullptr),
+    dampingModel_(nullptr),
+    isotropyModel_(nullptr),
+
     UIntegrator_(nullptr),
     UTrans_(nullptr),
     UCoeff_(nullptr)
@@ -687,11 +747,26 @@ void Foam::KinematicCloud<CloudType>::preEvolve
 
     this->dispersion().cacheFields(true);
     forces_.cacheFields(true);
-    updateCellOccupancy();
 
     pAmbient_ = constProps_.dict().template
         getOrDefault<scalar>("pAmbient", pAmbient_);
 
+    if (this->dampingModel().active() || this->packingModel().active())
+    {
+         const_cast<typename parcelType::trackingData&>(td).updateAverages(*this);
+    }
+
+    if (this->dampingModel().active())
+    {
+        this->dampingModel().cacheFields(true);
+    }
+    if (this->packingModel().active())
+    {
+        this->packingModel().cacheFields(true);
+    }
+
+    updateCellOccupancy();
+
     functions_.preEvolve(td);
 }
 
@@ -702,7 +777,6 @@ void Foam::KinematicCloud<CloudType>::evolve()
     if (solution_.canEvolve())
     {
         typename parcelType::trackingData td(*this);
-
         solve(*this, td);
     }
 }
@@ -719,6 +793,12 @@ void Foam::KinematicCloud<CloudType>::motion
     td.part() = parcelType::trackingData::tpLinearTrack;
     CloudType::move(cloud, td, solution_.trackTime());
 
+    if (isotropyModel_->active())
+    {
+        td.updateAverages(cloud);
+        isotropyModel_->calculate();
+    }
+
     updateCellOccupancy();
 }
 
@@ -749,15 +829,15 @@ void Foam::KinematicCloud<CloudType>::patchData
         // just inside the domain rather than that of the wall itself.
         if (U_.boundaryField()[patchi].fixesValue())
         {
-            const vector Uw1 = U_.boundaryField()[patchi][patchFacei];
+            const vector Uw1(U_.boundaryField()[patchi][patchFacei]);
             const vector& Uw0 =
                 U_.oldTime().boundaryField()[patchi][patchFacei];
 
             const scalar f = p.currentTimeFraction();
 
-            const vector Uw = Uw0 + f*(Uw1 - Uw0);
+            const vector Uw(Uw0 + f*(Uw1 - Uw0));
 
-            const tensor nnw = nw*nw;
+            const tensor nnw(nw*nw);
 
             Up = (nnw & Up) + Uw - (nnw & Uw);
         }
@@ -813,6 +893,48 @@ void Foam::KinematicCloud<CloudType>::info()
     injectors_.info(Info);
     this->surfaceFilm().info(Info);
     this->patchInteraction().info(Info);
+
+    if (this->packingModel().active())
+    {
+        tmp<volScalarField> alpha = this->theta();
+
+        if (this->db().time().writeTime())
+        {
+            alpha().write();
+        }
+
+        const scalar alphaMin = gMin(alpha().primitiveField());
+        const scalar alphaMax = gMax(alpha().primitiveField());
+
+        Info<< "    Min cell volume fraction        = " << alphaMin << endl;
+        Info<< "    Max cell volume fraction        = " << alphaMax << endl;
+
+        if (alphaMax < SMALL)
+        {
+            return;
+        }
+
+        scalar nMin = GREAT;
+
+        forAll(this->mesh().cells(), celli)
+        {
+            const label n = this->cellOccupancy()[celli].size();
+
+            if (n > 0)
+            {
+                const scalar nPack = n*alphaMax/alpha()[celli];
+
+                if (nPack < nMin)
+                {
+                    nMin = nPack;
+                }
+            }
+        }
+
+        reduce(nMin, minOp<scalar>());
+
+        Info<< "    Min dense number of parcels     = " << nMin << endl;
+    }
 }
 
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -95,6 +95,15 @@ class SurfaceFilmModel;
 template<class CloudType>
 class StochasticCollisionModel;
 
+template<class CloudType>
+class PackingModel;
+
+template<class CloudType>
+class DampingModel;
+
+template<class CloudType>
+class IsotropyModel;
+
 
 /*---------------------------------------------------------------------------*\
                        Class KinematicCloud Declaration
@@ -129,7 +138,7 @@ public:
 
 private:
 
-    // Private data
+    // Private Data
 
         //- Cloud copy pointer
         autoPtr<KinematicCloud<CloudType>> cloudCopyPtr_;
@@ -146,7 +155,7 @@ private:
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- References to the mesh and time databases
         const fvMesh& mesh_;
@@ -225,6 +234,18 @@ protected:
             autoPtr<SurfaceFilmModel<KinematicCloud<CloudType>>>
                 surfaceFilmModel_;
 
+            //- Packing model
+            autoPtr<PackingModel<KinematicCloud<CloudType>>>
+                packingModel_;
+
+            //- Damping model
+            autoPtr<DampingModel<KinematicCloud<CloudType>>>
+                dampingModel_;
+
+            //- Exchange model
+            autoPtr<IsotropyModel<KinematicCloud<CloudType>>>
+                isotropyModel_;
+
 
         // Reference to the particle integration schemes
 
@@ -464,6 +485,31 @@ public:
                     surfaceFilm();
 
 
+                //- Return const access to the packing model
+                inline const PackingModel<KinematicCloud<CloudType>>&
+                    packingModel() const;
+
+                //- Return a reference to the packing model
+                inline PackingModel<KinematicCloud<CloudType>>&
+                    packingModel();
+
+                //- Return const access to the damping model
+                inline const DampingModel<KinematicCloud<CloudType>>&
+                    dampingModel() const;
+
+                //- Return a reference to the damping model
+                inline DampingModel<KinematicCloud<CloudType>>&
+                    dampingModel();
+
+                //- Return const access to the isotropy model
+                inline const IsotropyModel<KinematicCloud<CloudType>>&
+                    isotropyModel() const;
+
+                //- Return a reference to the isotropy model
+                inline IsotropyModel<KinematicCloud<CloudType>>&
+                    isotropyModel();
+
+
             // Integration schemes
 
                 //-Return reference to velocity integration
@@ -590,7 +636,10 @@ public:
             void scaleSources();
 
             //- Pre-evolve
-            void preEvolve(const typename parcelType::trackingData& td);
+            void preEvolve
+            (
+                const typename parcelType::trackingData& td
+            );
 
             //- Evolve the cloud
             void evolve();

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -255,6 +255,54 @@ Foam::KinematicCloud<CloudType>::surfaceFilm()
 }
 
 
+template<class CloudType>
+inline const Foam::PackingModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::packingModel() const
+{
+    return *packingModel_;
+}
+
+
+template<class CloudType>
+inline Foam::PackingModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::packingModel()
+{
+    return *packingModel_;
+}
+
+
+template<class CloudType>
+inline const Foam::DampingModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::dampingModel() const
+{
+    return *dampingModel_;
+}
+
+
+template<class CloudType>
+inline Foam::DampingModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::dampingModel()
+{
+    return *dampingModel_;
+}
+
+
+template<class CloudType>
+inline const Foam::IsotropyModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::isotropyModel() const
+{
+    return *isotropyModel_;
+}
+
+
+template<class CloudType>
+inline Foam::IsotropyModel<Foam::KinematicCloud<CloudType>>&
+Foam::KinematicCloud<CloudType>::isotropyModel()
+{
+    return *isotropyModel_;
+}
+
+
 template<class CloudType>
 inline const Foam::integrationScheme&
 Foam::KinematicCloud<CloudType>::UIntegrator() const

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -90,6 +91,28 @@ void Foam::KinematicParcel<ParcelType>::calcDispersion
 }
 
 
+template<class ParcelType>
+template<class TrackCloudType>
+void Foam::KinematicParcel<ParcelType>::calcUCorrection
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar dt
+)
+{
+    typename TrackCloudType::parcelType& p =
+        static_cast<typename TrackCloudType::parcelType&>(*this);
+
+    this->UCorrect_ = Zero;
+
+    this->UCorrect_ =
+        cloud.dampingModel().velocityCorrection(p, dt);
+
+    this->UCorrect_ +=
+        cloud.packingModel().velocityCorrection(p, dt);
+}
+
+
 template<class ParcelType>
 template<class TrackCloudType>
 void Foam::KinematicParcel<ParcelType>::cellValueSourceCorrection
@@ -141,6 +164,7 @@ void Foam::KinematicParcel<ParcelType>::calc
     this->U_ =
         calcVelocity(cloud, td, dt, Re, td.muc(), mass0, Su, dUTrans, Spu);
 
+    this->U_ += this->UCorrect_;
 
     // Accumulate carrier phase source terms
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -249,7 +273,8 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(p.rho_),
     age_(p.age_),
     tTurb_(p.tTurb_),
-    UTurb_(p.UTurb_)
+    UTurb_(p.UTurb_),
+    UCorrect_(p.UCorrect_)
 {}
 
 
@@ -270,7 +295,8 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(p.rho_),
     age_(p.age_),
     tTurb_(p.tTurb_),
-    UTurb_(p.UTurb_)
+    UTurb_(p.UTurb_),
+    UCorrect_(p.UCorrect_)
 {}
 
 
@@ -348,6 +374,8 @@ bool Foam::KinematicParcel<ParcelType>::move
                 p.cellValueSourceCorrection(cloud, ttd, dt);
             }
 
+            p.calcUCorrection(cloud, ttd, dt);
+
             p.calc(cloud, ttd, dt);
         }
 
@@ -357,7 +385,6 @@ bool Foam::KinematicParcel<ParcelType>::move
         {
             cloud.functions().postFace(p, ttd.keepParticle);
         }
-
         cloud.functions().postMove(p, dt, start, ttd.keepParticle);
 
         if (p.active() && p.onFace() && ttd.keepParticle)

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,8 +55,8 @@ SourceFiles
 #include "autoPtr.H"
 #include "interpolation.H"
 #include "demandDrivenEntry.H"
-
+#include "labelFieldIOField.H"
-// #include "ParticleForceList.H" // TODO
+#include "vectorFieldIOField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,6 +66,9 @@ namespace Foam
 template<class ParcelType>
 class KinematicParcel;
 
+template<class Type>
+class AveragingMethod;
+
 // Forward declaration of friend functions
 
 template<class ParcelType>
@@ -84,11 +87,12 @@ class KinematicParcel
 :
     public ParcelType
 {
-    // Private data
+    // Private Data
 
         //- Number of particle tracking attempts before we assume that it stalls
         static label maxTrackAttempts;
 
+
 public:
 
     //- Size in bytes of the fields
@@ -100,7 +104,7 @@ public:
     {
     protected:
 
-        // Protected data
+        // Protected Data
 
             //- Constant properties dictionary
             const dictionary dict_;
@@ -108,7 +112,7 @@ public:
 
     private:
 
-        // Private data
+        // Private Data
 
             //- Parcel type id - used for post-processing to flag the type
             //- of parcels issued by this cloud
@@ -138,7 +142,7 @@ public:
             constantProperties(const dictionary& parentDict);
 
 
-        // Member functions
+        // Member Functions
 
             //- Return const access to the constant properties dictionary
             inline const dictionary& dict() const;
@@ -173,7 +177,7 @@ public:
 
     private:
 
-        // Private data
+        // Private Data
 
             // Interpolators for continuous phase fields
 
@@ -199,6 +203,30 @@ public:
                 scalar muc_;
 
 
+            // MPPIC Averages
+
+                //- Volume average
+                autoPtr<AveragingMethod<scalar>> volumeAverage_;
+
+                //- Radius average [ volume^(1/3) ]
+                autoPtr<AveragingMethod<scalar>> radiusAverage_;
+
+                //- Density average
+                autoPtr<AveragingMethod<scalar>> rhoAverage_;
+
+                //- Velocity average
+                autoPtr<AveragingMethod<vector>> uAverage_;
+
+                //- Magnitude velocity squared average
+                autoPtr<AveragingMethod<scalar>> uSqrAverage_;
+
+                //- Frequency average
+                autoPtr<AveragingMethod<scalar>> frequencyAverage_;
+
+                //- Mass average
+                autoPtr<AveragingMethod<scalar>> massAverage_;
+
+
             //- Local gravitational or other body-force acceleration
             const vector& g_;
 
@@ -220,7 +248,7 @@ public:
             );
 
 
-        // Member functions
+        // Member Functions
 
             //- Return const access to the interpolator for continuous
             //- phase density field
@@ -260,12 +288,17 @@ public:
 
             //- Return access to the part of the tracking operation taking place
             inline trackPart& part();
+
+            //- Update the MPPIC averages
+            template<class TrackCloudType>
+            inline void updateAverages(const TrackCloudType& cloud);
+
     };
 
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         // Parcel properties
 
@@ -300,6 +333,9 @@ protected:
             //- Turbulent velocity fluctuation [m/s]
             vector UTurb_;
 
+            //- Velocity correction due to collisions MPPIC [m/s]
+            vector UCorrect_;
+
 
     // Protected Member Functions
 
@@ -340,6 +376,7 @@ public:
           + " age"
           + " tTurb"
           + " (UTurbx UTurby UTurbz)"
+          + " (UCorrectx UCorrecty UCorrectz)"
         );
 
 
@@ -465,6 +502,9 @@ public:
             //- Return const access to turbulent velocity fluctuation
             inline const vector& UTurb() const;
 
+            //- Return const access to correction velocity
+            inline const vector& UCorrect() const;
+
 
         // Edit
 
@@ -498,6 +538,9 @@ public:
             //- Return access to turbulent velocity fluctuation
             inline vector& UTurb();
 
+            //- Return access to correction velocity
+            inline vector& UCorrect();
+
 
         // Helper functions
 
@@ -611,6 +654,15 @@ public:
                 const scalar dt
             );
 
+            //- Correct U following MP-PIC sub-models
+            template<class TrackCloudType>
+            void calcUCorrection
+            (
+                TrackCloudType& cloud,
+                trackingData& td,
+                const scalar dt
+            );
+
 
         // Tracking
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -304,6 +304,20 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::U()
 }
 
 
+template<class ParcelType>
+inline const Foam::vector& Foam::KinematicParcel<ParcelType>::UCorrect() const
+{
+    return UCorrect_;
+}
+
+
+template<class ParcelType>
+inline Foam::vector& Foam::KinematicParcel<ParcelType>::UCorrect()
+{
+    return UCorrect_;
+}
+
+
 template<class ParcelType>
 inline Foam::scalar& Foam::KinematicParcel<ParcelType>::rho()
 {

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2019 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -67,7 +67,8 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     rho_(0.0),
     age_(0.0),
     tTurb_(0.0),
-    UTurb_(Zero)
+    UTurb_(Zero),
+    UCorrect_(Zero)
 {
     if (readFields)
     {
@@ -82,7 +83,8 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
                 >> rho_
                 >> age_
                 >> tTurb_
-                >> UTurb_;
+                >> UTurb_
+                >> UCorrect_;
         }
         else if (!is.checkLabelSize<>() || !is.checkScalarSize<>())
         {
@@ -100,6 +102,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
             readRawScalar(is, &age_);
             readRawScalar(is, &tTurb_);
             readRawScalar(is, UTurb_.data(), vector::nComponents);
+            readRawScalar(is, UCorrect_.data(), vector::nComponents);
 
             is.endRawRead();
         }
@@ -191,6 +194,13 @@ void Foam::KinematicParcel<ParcelType>::readFields(CloudType& c)
     );
     c.checkFieldIOobject(c, UTurb);
 
+    IOField<vector> UCorrect
+    (
+        c.fieldIOobject("UCorrect", IOobject::MUST_READ),
+        valid
+    );
+    c.checkFieldIOobject(c, UCorrect);
+
     label i = 0;
 
     for (KinematicParcel<ParcelType>& p : c)
@@ -205,6 +215,7 @@ void Foam::KinematicParcel<ParcelType>::readFields(CloudType& c)
         p.age_ = age[i];
         p.tTurb_ = tTurb[i];
         p.UTurb_ = UTurb[i];
+        p.UCorrect_ = UCorrect[i];
 
         ++i;
     }
@@ -234,6 +245,7 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const CloudType& c)
     IOField<scalar> age(c.fieldIOobject("age", IOobject::NO_READ), np);
     IOField<scalar> tTurb(c.fieldIOobject("tTurb", IOobject::NO_READ), np);
     IOField<vector> UTurb(c.fieldIOobject("UTurb", IOobject::NO_READ), np);
+    IOField<vector> UCorrect(c.fieldIOobject("UCorrect", IOobject::NO_READ), np);
 
     label i = 0;
 
@@ -249,6 +261,7 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const CloudType& c)
         age[i] = p.age();
         tTurb[i] = p.tTurb();
         UTurb[i] = p.UTurb();
+        UCorrect[i] = p.UCorrect();
 
         ++i;
     }
@@ -263,6 +276,7 @@ void Foam::KinematicParcel<ParcelType>::writeFields(const CloudType& c)
     age.write(valid);
     tTurb.write(valid);
     UTurb.write(valid);
+    UCorrect.write(valid);
 }
 
 
@@ -291,6 +305,7 @@ void Foam::KinematicParcel<ParcelType>::writeProperties
     writeProp("age", age_);
     writeProp("tTurb", tTurb_);
     writeProp("UTurb", UTurb_);
+    writeProp("UCorrect", UCorrect_);
 
     #undef writeProp
 }
@@ -318,6 +333,7 @@ void Foam::KinematicParcel<ParcelType>::readObjects
     const auto& age = cloud::lookupIOField<scalar>("age", obr);
     const auto& tTurb = cloud::lookupIOField<scalar>("tTurb", obr);
     const auto& UTurb = cloud::lookupIOField<vector>("UTurb", obr);
+    const auto& UCorrect = cloud::lookupIOField<vector>("UCorrect", obr);
 
     label i = 0;
 
@@ -333,6 +349,7 @@ void Foam::KinematicParcel<ParcelType>::readObjects
         p.age_ = age[i];
         p.tTurb_ = tTurb[i];
         p.UTurb_ = UTurb[i];
+        p.UCorrect_ = UCorrect[i];
 
         ++i;
     }
@@ -361,6 +378,7 @@ void Foam::KinematicParcel<ParcelType>::writeObjects
     auto& age = cloud::createIOField<scalar>("age", np, obr);
     auto& tTurb = cloud::createIOField<scalar>("tTurb", np, obr);
     auto&& UTurb = cloud::createIOField<vector>("UTurb", np, obr);
+    auto&& UCorrect = cloud::createIOField<vector>("UCorrect", np, obr);
 
     label i = 0;
 
@@ -376,6 +394,7 @@ void Foam::KinematicParcel<ParcelType>::writeObjects
         age[i] = p.age();
         tTurb[i] = p.tTurb();
         UTurb[i] = p.UTurb();
+        UCorrect[i] = p.UCorrect();
 
         ++i;
     }
@@ -403,7 +422,8 @@ Foam::Ostream& Foam::operator<<
             << token::SPACE << p.rho()
             << token::SPACE << p.age()
             << token::SPACE << p.tTurb()
-            << token::SPACE << p.UTurb();
+            << token::SPACE << p.UTurb()
+            << token::SPACE << p.UCorrect();
     }
     else
     {

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelTrackingDataI.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,6 +26,10 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "AveragingMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 template<class ParcelType>
 template<class TrackCloudType>
 inline Foam::KinematicParcel<ParcelType>::trackingData::trackingData
@@ -61,6 +66,106 @@ inline Foam::KinematicParcel<ParcelType>::trackingData::trackingData
     rhoc_(Zero),
     Uc_(Zero),
     muc_(Zero),
+
+    volumeAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":volumeAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    radiusAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":radiusAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    rhoAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":rhoAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    uAverage_
+    (
+        AveragingMethod<vector>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":uAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    uSqrAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":uSqrAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    frequencyAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":frequencyAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+    massAverage_
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":massAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    ),
+
     g_(cloud.g().value()),
     part_(part)
 {}
@@ -158,4 +263,107 @@ Foam::KinematicParcel<ParcelType>::trackingData::part()
 }
 
 
+template<class ParcelType>
+template<class TrackCloudType>
+inline void Foam::KinematicParcel<ParcelType>::trackingData::
+updateAverages
+(
+    const TrackCloudType& cloud
+)
+{
+    // zero the sums
+    volumeAverage_() = 0;
+    radiusAverage_() = 0;
+    rhoAverage_() = 0;
+    uAverage_() = Zero;
+    uSqrAverage_() = 0;
+    frequencyAverage_() = 0;
+    massAverage_() = 0;
+
+    // temporary weights
+    autoPtr<AveragingMethod<scalar>> weightAveragePtr
+    (
+        AveragingMethod<scalar>::New
+        (
+            IOobject
+            (
+                cloud.name() + ":weightAverage",
+                cloud.db().time().timeName(),
+                cloud.mesh()
+            ),
+            cloud.solution().dict(),
+            cloud.mesh()
+        )
+    );
+    AveragingMethod<scalar>& weightAverage = weightAveragePtr();
+
+    // averaging sums
+    for (const typename TrackCloudType::parcelType& p : cloud)
+    {
+        const tetIndices tetIs = p.currentTetIndices();
+
+        const scalar m = p.nParticle()*p.mass();
+
+        volumeAverage_->add(p.coordinates(), tetIs, p.nParticle()*p.volume());
+        rhoAverage_->add(p.coordinates(), tetIs, m*p.rho());
+        uAverage_->add(p.coordinates(), tetIs, m*p.U());
+        massAverage_->add(p.coordinates(), tetIs, m);
+    }
+    volumeAverage_->average();
+    massAverage_->average();
+    rhoAverage_->average(*massAverage_);
+    uAverage_->average(*massAverage_);
+
+    // squared velocity deviation
+    for (const typename TrackCloudType::parcelType& p : cloud)
+    {
+        const tetIndices tetIs = p.currentTetIndices();
+
+        const vector u = uAverage_->interpolate(p.coordinates(), tetIs);
+
+        uSqrAverage_->add
+        (
+            p.coordinates(),
+            tetIs,
+            p.nParticle()*p.mass()*magSqr(p.U() - u)
+        );
+    }
+    uSqrAverage_->average(*massAverage_);
+
+    // sauter mean radius
+    radiusAverage_() = volumeAverage_();
+    weightAverage = 0;
+    for (const typename TrackCloudType::parcelType& p : cloud)
+    {
+        const tetIndices tetIs = p.currentTetIndices();
+
+        weightAverage.add
+        (
+            p.coordinates(),
+            tetIs,
+            p.nParticle()*pow(p.volume(), 2.0/3.0)
+        );
+    }
+    weightAverage.average();
+    radiusAverage_->average(weightAverage);
+
+    // collision frequency
+    weightAverage = 0;
+    for (const typename TrackCloudType::parcelType& p : cloud)
+    {
+        const tetIndices tetIs = p.currentTetIndices();
+
+        const scalar a = volumeAverage_->interpolate(p.coordinates(), tetIs);
+        const scalar r = radiusAverage_->interpolate(p.coordinates(), tetIs);
+        const vector u = uAverage_->interpolate(p.coordinates(), tetIs);
+
+        const scalar f = 0.75*a/pow3(r)*sqr(0.5*p.d() + r)*mag(p.U() - u);
+
+        frequencyAverage_->add(p.coordinates(), tetIs, p.nParticle()*f*f);
+
+        weightAverage.add(p.coordinates(), tetIs, p.nParticle()*f);
+    }
+    frequencyAverage_->average(weightAverage);
+}
+
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicHeterogeneousReactingParcel/makeBasicHeterogeneousReactingParcelSubmodels.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,6 +45,11 @@ License
 #include "makeReactingMultiphaseParcelCompositionModels.H"
 #include "makeReactingParcelPhaseChangeModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicHeterogeneousReactingCloud);
@@ -68,4 +73,9 @@ makeHeterogeneousReactingParcelHeterogeneousReactingModels
     basicHeterogeneousReactingCloud
 );
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicHeterogeneousReactingCloud);
+makeMPPICParcelIsotropyModels(basicHeterogeneousReactingCloud);
+makeMPPICParcelPackingModels(basicHeterogeneousReactingCloud);
+
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicKinematicCollidingParcel/makeBasicKinematicCollidingParcelSubmodels.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,6 +39,11 @@ License
 #include "makeParcelStochasticCollisionModels.H"
 #include "makeParcelSurfaceFilmModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicKinematicCollidingCloud);
@@ -51,5 +57,9 @@ makeParcelPatchInteractionModels(basicKinematicCollidingCloud);
 makeParcelStochasticCollisionModels(basicKinematicCollidingCloud);
 makeParcelSurfaceFilmModels(basicKinematicCollidingCloud);
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicKinematicCollidingCloud);
+makeMPPICParcelIsotropyModels(basicKinematicCollidingCloud);
+makeMPPICParcelPackingModels(basicKinematicCollidingCloud);
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicKinematicMPPICParcel/makeBasicKinematicMPPICParcelSubmodels.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,9 +57,15 @@ makeParcelStochasticCollisionModels(basicKinematicMPPICCloud);
 makeParcelSurfaceFilmModels(basicKinematicMPPICCloud);
 
 // MPPIC sub-models
+// WIP: These models are defined in Kinematic and MPPIC clouds temporarily
 makeMPPICParcelDampingModels(basicKinematicMPPICCloud);
+makeMPPICCloudParcelDampingModels(basicKinematicMPPICCloud);
+
 makeMPPICParcelIsotropyModels(basicKinematicMPPICCloud);
+makeMPPICCloudParcelIsotropyModels(basicKinematicMPPICCloud);
+
 makeMPPICParcelPackingModels(basicKinematicMPPICCloud);
+makeMPPICCloudParcelPackingModels(basicKinematicMPPICCloud);
 
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelSubmodels.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,6 +38,11 @@ License
 #include "makeParcelStochasticCollisionModels.H"
 #include "makeParcelSurfaceFilmModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicKinematicCloud);
@@ -49,5 +55,10 @@ makeParcelPatchInteractionModels(basicKinematicCloud);
 makeParcelStochasticCollisionModels(basicKinematicCloud);
 makeParcelSurfaceFilmModels(basicKinematicCloud);
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicKinematicCloud);
+makeMPPICParcelIsotropyModels(basicKinematicCloud);
+makeMPPICParcelPackingModels(basicKinematicCloud);
+
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/makeBasicReactingMultiphaseParcelSubmodels.C

@@ -49,6 +49,11 @@ License
 #include "makeReactingMultiphaseParcelDevolatilisationModels.H"
 #include "makeReactingMultiphaseParcelSurfaceReactionModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicReactingMultiphaseCloud);
@@ -84,5 +89,9 @@ makeReactingMultiphaseParcelSurfaceReactionModels
     basicReactingMultiphaseCloud
 );
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicReactingMultiphaseCloud);
+makeMPPICParcelIsotropyModels(basicReactingMultiphaseCloud);
+makeMPPICParcelPackingModels(basicReactingMultiphaseCloud);
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicReactingParcel/makeBasicReactingParcelSubmodels.C

@@ -45,6 +45,11 @@ License
 #include "makeReactingParcelCompositionModels.H"
 #include "makeReactingParcelPhaseChangeModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicReactingCloud);
@@ -64,5 +69,10 @@ makeParcelHeatTransferModels(basicReactingCloud);
 makeReactingParcelCompositionModels(basicReactingCloud);
 makeReactingParcelPhaseChangeModels(basicReactingCloud);
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicReactingCloud);
+makeMPPICParcelIsotropyModels(basicReactingCloud);
+makeMPPICParcelPackingModels(basicReactingCloud);
+
 
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/makeBasicThermoParcelSubmodels.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,6 +41,11 @@ License
 // Thermodynamic
 #include "makeParcelHeatTransferModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicThermoCloud);
@@ -56,4 +62,10 @@ makeParcelSurfaceFilmModels(basicThermoCloud);
 makeParcelHeatTransferModels(basicThermoCloud);
 
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicThermoCloud);
+makeMPPICParcelIsotropyModels(basicThermoCloud);
+makeMPPICParcelPackingModels(basicThermoCloud);
+
+
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/include/makeMPPICParcelDampingModels.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,6 +45,13 @@ License
     makeDampingModelType(Relaxation, CloudType);
 
 
+#define makeMPPICCloudParcelDampingModels(CloudType)                           \
+                                                                               \
+    makeDampingModelMPPIC(CloudType);                                          \
+                                                                               \
+    makeDampingModelTypeMPPIC(NoDamping, CloudType);                           \
+    makeDampingModelTypeMPPIC(Relaxation, CloudType);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/lagrangian/intermediate/parcels/include/makeMPPICParcelIsotropyModels.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,6 +44,13 @@ License
     makeIsotropyModelType(Stochastic, CloudType);
 
 
+#define makeMPPICCloudParcelIsotropyModels(CloudType)                          \
+                                                                               \
+    makeIsotropyModelMPPIC(CloudType);                                         \
+                                                                               \
+    makeIsotropyModelTypeMPPIC(NoIsotropy, CloudType);                         \
+    makeIsotropyModelTypeMPPIC(Stochastic, CloudType);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/lagrangian/intermediate/parcels/include/makeMPPICParcelPackingModels.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -44,6 +45,15 @@ License
     makePackingModelType(Explicit, CloudType);                                 \
     makePackingModelType(Implicit, CloudType);
 
+
+#define makeMPPICCloudParcelPackingModels(CloudType)                           \
+                                                                               \
+    makePackingModelMPPIC(CloudType);                                          \
+                                                                               \
+    makePackingModelTypeMPPIC(NoPacking, CloudType);                           \
+    makePackingModelTypeMPPIC(Explicit, CloudType);                            \
+    makePackingModelTypeMPPIC(Implicit, CloudType);
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C

@@ -294,6 +294,18 @@ bool Foam::LocalInteraction<CloudType>::correct
                 // Calculate motion relative to patch velocity
                 U -= Up;
 
+                if (mag(Up) > 0 && mag(U) < this->Urmax())
+                {
+                    WarningInFunction
+                        << "Particle U the same as patch "
+                        << "    The particle has been removed" << nl << endl;
+
+                    keepParticle = false;
+                    p.active(false);
+                    U = Zero;
+                    break;
+                }
+
                 scalar Un = U & nw;
                 vector Ut = U - Un*nw;
 

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C

@@ -137,7 +137,8 @@ Foam::PatchInteractionModel<CloudType>::PatchInteractionModel
     functionObjects::writeFile(owner, this->localPath(), typeName, false),
     UName_("unknown_U"),
     escapedParcels_(0),
-    escapedMass_(0.0)
+    escapedMass_(0.0),
+    Urmax_(1e-4)
 {}
 
 
@@ -160,7 +161,8 @@ Foam::PatchInteractionModel<CloudType>::PatchInteractionModel
     ),
     UName_(this->coeffDict().template getOrDefault<word>("U", "U")),
     escapedParcels_(0),
-    escapedMass_(0.0)
+    escapedMass_(0.0),
+    Urmax_(this->coeffDict().template getOrDefault<scalar>("UrMax", 0))
 {}
 
 
@@ -174,7 +176,8 @@ Foam::PatchInteractionModel<CloudType>::PatchInteractionModel
     functionObjects::writeFile(pim),
     UName_(pim.UName_),
     escapedParcels_(pim.escapedParcels_),
-    escapedMass_(pim.escapedMass_)
+    escapedMass_(pim.escapedMass_),
+    Urmax_(pim.Urmax_)
 {}
 
 
@@ -187,6 +190,13 @@ const Foam::word& Foam::PatchInteractionModel<CloudType>::UName() const
 }
 
 
+template<class CloudType>
+const Foam::scalar& Foam::PatchInteractionModel<CloudType>::Urmax() const
+{
+    return Urmax_;
+}
+
+
 template<class CloudType>
 void Foam::PatchInteractionModel<CloudType>::addToEscapedParcels
 (

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.H

@@ -99,6 +99,9 @@ protected:
             //- Mass of parcels escaped
             scalar escapedMass_;
 
+            //- Maximum relative U with patch for particle to be removed
+            scalar Urmax_;
+
 
     // Protected Member Functions
 
@@ -162,6 +165,8 @@ public:
         //- Return name of velocity field
         const word& UName() const;
 
+        //- Return Urmax
+        const scalar& Urmax() const;
 
     // Member Functions
 

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C

@@ -211,6 +211,18 @@ bool Foam::StandardWallInteraction<CloudType>::correct
                 // Calculate motion relative to patch velocity
                 U -= Up;
 
+                if (mag(Up) > 0 && mag(U) < this->Urmax())
+                {
+                    WarningInFunction
+                        << "Particle U the same as patch "
+                        << "    The particle has been removed" << nl << endl;
+
+                    keepParticle = false;
+                    p.active(false);
+                    U = Zero;
+                    break;
+                }
+
                 scalar Un = U & nw;
                 vector Ut = U - Un*nw;
 
@@ -224,6 +236,7 @@ bool Foam::StandardWallInteraction<CloudType>::correct
                 // Return velocity to global space
                 U += Up;
 
+
                 break;
             }
             default:

src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2017 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -87,9 +87,10 @@ Foam::AveragingMethod<Type>::New
     const fvMesh& mesh
 )
 {
-    const word modelType(dict.get<word>(typeName));
+    const word modelType
-
+    (
-    //Info<< "Selecting averaging method " << modelType << endl;
+        dict.template getOrDefault<word>(typeName, "basic")
+    );
 
     auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 

src/lagrangian/intermediate/submodels/MPPIC/DampingModels/DampingModel/DampingModel.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,7 +83,10 @@ Foam::DampingModel<CloudType>::New
     CloudType& owner
 )
 {
-    const word modelType(dict.get<word>(typeName));
+    const word modelType
+    (
+        dict.template getOrDefault<word>(typeName, "none")
+    );
 
     Info<< "Selecting damping model " << modelType << endl;
 

src/lagrangian/intermediate/submodels/MPPIC/DampingModels/DampingModel/DampingModel.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -142,6 +143,39 @@ public:
 
 #define makeDampingModel(CloudType)                                            \
                                                                                \
+    typedef Foam::CloudType::kinematicCloudType kinematicCloudType;            \
+    defineNamedTemplateTypeNameAndDebug                                        \
+    (                                                                          \
+        Foam::DampingModel<kinematicCloudType>,                                \
+        0                                                                      \
+    );                                                                         \
+                                                                               \
+    namespace Foam                                                             \
+    {                                                                          \
+        defineTemplateRunTimeSelectionTable                                    \
+        (                                                                      \
+            DampingModel<kinematicCloudType>,                                  \
+            dictionary                                                         \
+        );                                                                     \
+    }
+
+
+#define makeDampingModelType(SS, CloudType)                                    \
+                                                                               \
+    typedef Foam::CloudType::kinematicCloudType kinematicCloudType;            \
+    defineNamedTemplateTypeNameAndDebug                                        \
+        (Foam::DampingModels::SS<kinematicCloudType>, 0);                      \
+                                                                               \
+    Foam::DampingModel<kinematicCloudType>::                                   \
+        adddictionaryConstructorToTable                                        \
+        <Foam::DampingModels::SS<kinematicCloudType>>                          \
+            add##SS##CloudType##kinematicCloudType##ConstructorToTable_;
+
+
+// These are used to defined models in MPPIC cloud specifically
+
+#define makeDampingModelMPPIC(CloudType)                                       \
+                                                                               \
     typedef Foam::CloudType::MPPICCloudType MPPICCloudType;                    \
     defineNamedTemplateTypeNameAndDebug                                        \
     (                                                                          \
@@ -159,7 +193,7 @@ public:
     }
 
 
-#define makeDampingModelType(SS, CloudType)                                    \
+#define makeDampingModelTypeMPPIC(SS, CloudType)                               \
                                                                                \
     typedef Foam::CloudType::MPPICCloudType MPPICCloudType;                    \
     defineNamedTemplateTypeNameAndDebug                                        \

src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/IsotropyModel/IsotropyModel.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -86,7 +86,10 @@ Foam::IsotropyModel<CloudType>::New
     CloudType& owner
 )
 {
-    const word modelType(dict.get<word>(typeName));
+    const word modelType
+    (
+        dict.template getOrDefault<word>(typeName, "none")
+    );
 
     Info<< "Selecting isotropy model " << modelType << endl;
 

src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/IsotropyModel/IsotropyModel.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -136,14 +137,44 @@ public:
 
 #define makeIsotropyModel(CloudType)                                           \
                                                                                \
-    typedef Foam::CloudType::MPPICCloudType MPPICCloudType;                    \
+    typedef Foam::CloudType::kinematicCloudType kinematicCloudType;            \
     defineNamedTemplateTypeNameAndDebug                                        \
     (                                                                          \
-        Foam::IsotropyModel<MPPICCloudType>,                                   \
+        Foam::IsotropyModel<kinematicCloudType>,                               \
         0                                                                      \
     );                                                                         \
                                                                                \
     namespace Foam                                                             \
+    {                                                                          \
+        defineTemplateRunTimeSelectionTable                                    \
+        (                                                                      \
+            IsotropyModel<kinematicCloudType>,                                 \
+            dictionary                                                         \
+        );                                                                     \
+    }
+
+
+#define makeIsotropyModelType(SS, CloudType)                                   \
+                                                                               \
+    typedef Foam::CloudType::kinematicCloudType kinematicCloudType;            \
+    defineNamedTemplateTypeNameAndDebug                                        \
+        (Foam::IsotropyModels::SS<kinematicCloudType>, 0);                     \
+                                                                               \
+    Foam::IsotropyModel<kinematicCloudType>::                                  \
+        adddictionaryConstructorToTable                                        \
+        <Foam::IsotropyModels::SS<kinematicCloudType>>                         \
+            add##SS##CloudType##kinematicCloudType##ConstructorToTable_;
+
+
+// These are used to defined models in MPPIC cloud specifically
+
+#define makeIsotropyModelMPPIC(CloudType)                                  \
+                                                                               \
+    typedef Foam::CloudType::MPPICCloudType MPPICCloudType;                    \
+    defineNamedTemplateTypeNameAndDebug                                        \
+        (Foam::IsotropyModel<MPPICCloudType>, 0);                              \
+                                                                               \
+    namespace Foam                                                             \
     {                                                                          \
         defineTemplateRunTimeSelectionTable                                    \
         (                                                                      \
@@ -153,7 +184,7 @@ public:
     }
 
 
-#define makeIsotropyModelType(SS, CloudType)                                   \
+#define makeIsotropyModelTypeMPPIC(SS, CloudType)                              \
                                                                                \
     typedef Foam::CloudType::MPPICCloudType MPPICCloudType;                    \
     defineNamedTemplateTypeNameAndDebug                                        \
@@ -164,7 +195,6 @@ public:
         <Foam::IsotropyModels::SS<MPPICCloudType>>                             \
             add##SS##CloudType##MPPICCloudType##ConstructorToTable_;
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository

src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2017 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -150,8 +151,10 @@ Foam::vector Foam::PackingModels::Explicit<CloudType>::velocityCorrection
     // interpolated quantities
     const scalar alpha =
         this->volumeAverage_->interpolate(p.coordinates(), tetIs);
+
     const vector alphaGrad =
         this->volumeAverage_->interpolateGrad(p.coordinates(), tetIs);
+
     const vector uMean =
         this->uAverage_->interpolate(p.coordinates(), tetIs);
 
@@ -175,7 +178,7 @@ Foam::vector Foam::PackingModels::Explicit<CloudType>::velocityCorrection
     // correction velocity
     if ((uRelative & alphaGrad) > 0)
     {
-        dU = - deltaT*tauGrad/(p.rho()*alpha/* + deltaT*F.Sp()*/);
+        dU = - deltaT*tauGrad/(p.rho()*(alpha + SMALL)/* + deltaT*F.Sp()*/);
     }
 
     // apply the velocity limiters

src/lagrangian/intermediate/submodels/MPPIC/PackingModels/PackingModel/PackingModel.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -85,8 +85,10 @@ Foam::PackingModel<CloudType>::New
     CloudType& owner
 )
 {
-    const word modelType(dict.get<word>(typeName));
+    const word modelType
-
+    (
+        dict.template getOrDefault<word>(typeName, "none")
+    );
     Info<< "Selecting packing model " << modelType << endl;
 
     auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);

src/lagrangian/intermediate/submodels/MPPIC/PackingModels/PackingModel/PackingModel.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -147,6 +148,39 @@ public:
 
 #define makePackingModel(CloudType)                                            \
                                                                                \
+    typedef Foam::CloudType::kinematicCloudType kinematicCloudType;            \
+    defineNamedTemplateTypeNameAndDebug                                        \
+    (                                                                          \
+        Foam::PackingModel<kinematicCloudType>,                                \
+        0                                                                      \
+    );                                                                         \
+                                                                               \
+    namespace Foam                                                             \
+    {                                                                          \
+        defineTemplateRunTimeSelectionTable                                    \
+        (                                                                      \
+            PackingModel<kinematicCloudType>,                                  \
+            dictionary                                                         \
+        );                                                                     \
+    }
+
+
+#define makePackingModelType(SS, CloudType)                                    \
+                                                                               \
+    typedef Foam::CloudType::kinematicCloudType kinematicCloudType;            \
+    defineNamedTemplateTypeNameAndDebug                                        \
+        (Foam::PackingModels::SS<kinematicCloudType>, 0);                      \
+                                                                               \
+    Foam::PackingModel<kinematicCloudType>::                                   \
+        adddictionaryConstructorToTable                                        \
+        <Foam::PackingModels::SS<kinematicCloudType>>                          \
+            add##SS##CloudType##kinematicCloudType##ConstructorToTable_;
+
+
+// These are used to defined models in MPPIC cloud specifically
+
+#define makePackingModelMPPIC(CloudType)                                   \
+                                                                               \
     typedef Foam::CloudType::MPPICCloudType MPPICCloudType;                    \
     defineNamedTemplateTypeNameAndDebug                                        \
     (                                                                          \
@@ -164,7 +198,7 @@ public:
     }
 
 
-#define makePackingModelType(SS, CloudType)                                    \
+#define makePackingModelTypeMPPIC(SS, CloudType)                               \
                                                                                \
     typedef Foam::CloudType::MPPICCloudType MPPICCloudType;                    \
     defineNamedTemplateTypeNameAndDebug                                        \
@@ -175,7 +209,6 @@ public:
         <Foam::PackingModels::SS<MPPICCloudType>>                              \
             add##SS##CloudType##MPPICCloudType##ConstructorToTable_;
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository

src/lagrangian/spray/parcels/derived/basicSprayParcel/makeBasicSprayParcelSubmodels.C

@@ -52,6 +52,11 @@ License
 #include "makeSprayParcelAtomizationModels.H"
 #include "makeSprayParcelBreakupModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicSprayCloud);
@@ -78,5 +83,9 @@ makeParticleForceModelType(DistortedSphereDragForce, basicSprayCloud);
 makeSprayParcelAtomizationModels(basicSprayCloud);
 makeSprayParcelBreakupModels(basicSprayCloud);
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicSprayCloud);
+makeMPPICParcelIsotropyModels(basicSprayCloud);
+makeMPPICParcelPackingModels(basicSprayCloud);
 
 // ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -1069,15 +1069,17 @@ const
     {
         const face& f = srcPatch[srcFacei];
 
-        pointHit ray = f.ray(tgtPoint, n, srcPoints);
+        pointHit ray =
+            f.ray(tgtPoint, n, srcPoints, intersection::algorithm::VISIBLE);
 
         if (ray.hit())
         {
-            // tgtPoint = ray.rawPoint();
+            tgtPoint = ray.rawPoint();
             return srcFacei;
         }
         else if (ray.distance() < nearest.distance())
         {
+
             nearest = ray;
             nearestFacei = srcFacei;
         }
@@ -1085,7 +1087,7 @@ const
 
     if (nearest.hit() || nearest.eligibleMiss())
     {
-        // tgtPoint = nearest.rawPoint();
+        tgtPoint = nearest.rawPoint();
         return nearestFacei;
     }
 
@@ -1116,23 +1118,26 @@ const
     {
         const face& f = tgtPatch[tgtFacei];
 
-        pointHit ray = f.ray(srcPoint, n, tgtPoints);
+        pointHit ray =
+            f.ray(srcPoint, n, tgtPoints, intersection::algorithm::VISIBLE);
 
-        if (ray.hit() || ray.eligibleMiss())
+        if (ray.hit())
         {
-            // srcPoint = ray.rawPoint();
+            srcPoint = ray.rawPoint();
             return tgtFacei;
         }
-        else if (ray.distance() < nearest.distance())
+        const pointHit near = f.nearestPoint(srcPoint, tgtPoints);
+
+        if (near.distance() < nearest.distance())
         {
-            nearest = ray;
+            nearest = near;
             nearestFacei = tgtFacei;
         }
     }
-
     if (nearest.hit() || nearest.eligibleMiss())
     {
-        // srcPoint = nearest.rawPoint();
+
+        srcPoint = nearest.rawPoint();
         return nearestFacei;
     }
 

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C

@@ -568,6 +568,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     coupledPolyPatch(name, size, start, index, bm, patchType, transform),
     nbrPatchName_(word::null),
     nbrPatchID_(-1),
+    fraction_(Zero),
     rotationAxis_(Zero),
     rotationCentre_(Zero),
     rotationAngleDefined_(false),
@@ -603,6 +604,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     nbrPatchName_(dict.getOrDefault<word>("neighbourPatch", word::null)),
     coupleGroup_(dict),
     nbrPatchID_(-1),
+    fraction_(dict.getOrDefault<scalar>("fraction", Zero)),
     rotationAxis_(Zero),
     rotationCentre_(Zero),
     rotationAngleDefined_(false),
@@ -707,6 +709,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     nbrPatchName_(pp.nbrPatchName_),
     coupleGroup_(pp.coupleGroup_),
     nbrPatchID_(-1),
+    fraction_(pp.fraction_),
     rotationAxis_(pp.rotationAxis_),
     rotationCentre_(pp.rotationCentre_),
     rotationAngleDefined_(pp.rotationAngleDefined_),
@@ -742,6 +745,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     nbrPatchName_(nbrPatchName),
     coupleGroup_(pp.coupleGroup_),
     nbrPatchID_(-1),
+    fraction_(pp.fraction_),
     rotationAxis_(pp.rotationAxis_),
     rotationCentre_(pp.rotationCentre_),
     rotationAngleDefined_(pp.rotationAngleDefined_),
@@ -784,6 +788,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     nbrPatchName_(pp.nbrPatchName_),
     coupleGroup_(pp.coupleGroup_),
     nbrPatchID_(-1),
+    fraction_(pp.fraction_),
     rotationAxis_(pp.rotationAxis_),
     rotationCentre_(pp.rotationCentre_),
     rotationAngleDefined_(pp.rotationAngleDefined_),
@@ -1153,6 +1158,8 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
         os.writeEntry("tgtSize", tgtFaceIDs_.size());
         os.writeEntry("moveFaceCentres", moveFaceCentres_);
     }
+
+    os.writeEntryIfDifferent<scalar>("fraction", Zero, fraction_);
 }
 
 

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H

@@ -97,6 +97,9 @@ protected:
         //- Index of other half
         mutable label nbrPatchID_;
 
+        //- Particle displacement fraction accross AMI
+        const scalar fraction_;
+
 
         // Transformations
 
@@ -357,6 +360,9 @@ public:
             //- Neighbour patch ID
             virtual label neighbPatchID() const;
 
+            //- Particle fraction increase between AMI pathces
+            inline scalar fraction() const;
+
             //- Does this side own the patch?
             virtual bool owner() const;
 

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H

@@ -58,6 +58,13 @@ inline const Foam::word& Foam::cyclicAMIPolyPatch::neighbPatchName() const
     return nbrPatchName_;
 }
 
+
+inline Foam::scalar Foam::cyclicAMIPolyPatch::fraction() const
+{
+    return fraction_;
+}
+
+
 inline const Foam::scalarListList& Foam::cyclicAMIPolyPatch::weights() const
 {
     if (owner())

tutorials/lagrangian/MPPICDyMFoam/column/0/U.air

@@ -0,0 +1,40 @@
+/*--------------------------------*- 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      binary;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    sides
+    {
+        type            symmetry;
+    }
+
+    top
+    {
+        type            noSlip;
+    }
+
+    bottom
+    {
+        type            noSlip;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/0/p

@@ -0,0 +1,42 @@
+/*--------------------------------*- 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       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    sides
+    {
+        type            symmetry;
+    }
+
+    top
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+
+    bottom
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/constant/dynamicMeshDict

@@ -0,0 +1,27 @@
+/*--------------------------------*- 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;
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicMotionSolverFvMesh;
+
+motionSolver    solidBody;
+
+solidBodyMotionFunction  oscillatingLinearMotion;
+
+amplitude       (3 0 0);
+omega           0.1; // rad/s  (.5 rps)
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/constant/g

@@ -0,0 +1,21 @@
+/*--------------------------------*- 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       uniformDimensionedVectorField;
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 0 -9.81);
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/constant/kinematicCloudProperties

@@ -0,0 +1,187 @@
+/*--------------------------------*- 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;
+    object      kinematicCloudProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solution
+{
+    active          true;
+    coupled         true;
+    transient       yes;
+    cellValueSourceCorrection off;
+
+    interpolationSchemes
+    {
+        rho.air         cell;
+        U.air           cellPoint;
+        mu.air          cell;
+    }
+
+    averagingMethod dual;
+
+    integrationSchemes
+    {
+        U               Euler;
+    }
+
+    sourceTerms
+    {
+        schemes
+        {
+            U     semiImplicit 1;
+        }
+    }
+}
+
+constantProperties
+{
+    rho0            2526;
+    alphaMax        0.9;
+}
+
+subModels
+{
+    particleForces
+    {
+        ErgunWenYuDrag
+        {
+            alphac alpha.air;
+        }
+    }
+
+    injectionModels
+    {
+        model1
+        {
+            type            manualInjection;
+            massTotal       0;
+            parcelBasisType fixed;
+            nParticle       750;
+            SOI             0;
+            positionsFile   "kinematicCloudPositions";
+            U0              (0 0 0);
+            sizeDistribution
+            {
+                type        fixedValue;
+                fixedValueDistribution
+                {
+                    value   0.0003;
+                }
+            }
+        }
+    }
+
+    dispersionModel none;
+
+    patchInteractionModel localInteraction;
+
+    localInteractionCoeffs
+    {
+        patches
+        (
+            sides
+            {
+                type rebound;
+                e    1;
+                mu   0.09;
+            }
+            top
+            {
+                type rebound;
+                e    1;
+                mu   0.09;
+            }
+            bottom
+            {
+                type rebound;
+                e    1;
+                mu   0.09;
+            }
+        );
+    }
+
+    heatTransferModel none;
+
+    surfaceFilmModel none;
+
+    packingModel implicit;
+
+    explicitCoeffs
+    {
+        particleStressModel
+        {
+            type HarrisCrighton;
+            alphaPacked 0.6;
+            pSolid 10.0;
+            beta 2.0;
+            eps 1.0e-7;
+        }
+        correctionLimitingMethod
+        {
+            type absolute;
+            e 0.9;
+        }
+    }
+
+    implicitCoeffs
+    {
+        alphaMin 0.0001;
+        rhoMin 1.0;
+        applyLimiting true;
+        applyGravity true;
+        particleStressModel
+        {
+            type HarrisCrighton;
+            alphaPacked 0.6;
+            pSolid 10.0;
+            beta 2.0;
+            eps 1.0e-2;
+        }
+    }
+
+    dampingModel none;//relaxation;
+
+    relaxationCoeffs
+    {
+        timeScaleModel
+        {
+            type nonEquilibrium;
+            alphaPacked 0.6;
+            e 0.9;
+        }
+    }
+
+    isotropyModel stochastic;
+
+    stochasticCoeffs
+    {
+        timeScaleModel
+        {
+            type isotropic;
+            alphaPacked 0.6;
+            e 0.9;
+        }
+    }
+
+    stochasticCollisionModel none;
+
+    radiation off;
+}
+
+
+cloudFunctions
+{}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/constant/transportProperties

@@ -0,0 +1,25 @@
+/*--------------------------------*- 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;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+continuousPhase air;
+
+rho.air         1.2;
+
+transportModel  Newtonian;
+nu              1e-05;
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/constant/turbulenceProperties.air

@@ -0,0 +1,20 @@
+/*--------------------------------*- 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;
+    object      turbulenceProperties.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/system/blockMeshDict

@@ -0,0 +1,74 @@
+/*--------------------------------*- 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;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   0.001;
+
+vertices
+(
+    (-9.2 -9.2   0)
+    ( 9.2 -9.2   0)
+    ( 9.2  9.2   0)
+    (-9.2  9.2   0)
+    (-9.2 -9.2 300)
+    ( 9.2 -9.2 300)
+    ( 9.2  9.2 300)
+    (-9.2  9.2 300)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (1 1 40) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    sides
+    {
+        type symmetry;
+        faces
+        (
+            (0 1 5 4)
+            (1 2 6 5)
+            (2 3 7 6)
+            (3 0 4 7)
+        );
+    }
+
+    top
+    {
+        type wall;
+        faces
+        (
+            (4 5 6 7)
+        );
+    }
+
+    bottom
+    {
+        type wall;
+        faces
+        (
+            (0 1 2 3)
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/system/controlDict

@@ -0,0 +1,48 @@
+/*--------------------------------*- 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;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     MPPICFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         1;
+
+deltaT          2e-4;
+
+writeControl    runTime;
+
+writeInterval   0.01;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  8;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/system/fvSchemes

@@ -0,0 +1,58 @@
+/*--------------------------------*- 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;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+d2dt2Schemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(alphaPhi.air,U.air)    Gauss linearUpwindV unlimited;
+    div(((alpha.air*nuEff.air)*dev2(T(grad(U.air))))) Gauss linear;
+
+    div(phiGByA,kinematicCloud:alpha) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/system/fvSolution

@@ -0,0 +1,89 @@
+/*--------------------------------*- 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;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          GAMG;
+        tolerance       1e-6;
+        relTol          0.01;
+        smoother        GaussSeidel;
+    }
+
+    pFinal
+    {
+        $p;
+        relTol          0;
+    }
+
+    pcorrFinal
+    {
+        $p;
+        relTol          0;
+    }
+
+    U.air
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    U.airFinal
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    kinematicCloud:alpha
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nOuterCorrectors 1;
+    nCorrectors     2;
+    momentumPredictor yes;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+    checkMeshCourantNo  no;
+    correctPhi          no;
+}
+
+relaxationFactors
+{
+    fields
+    {
+        ".*"            1;
+    }
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/0/T

@@ -0,0 +1,48 @@
+/*--------------------------------*- 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       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 350;
+
+boundaryField
+{
+    AMI1a
+    {
+        type            cyclicAMI;
+        value           $internalField;
+    }
+
+    AMI1b
+    {
+        type            cyclicAMI;
+        value           $internalField;
+    }
+
+    ".*"
+    {
+        type            zeroGradient;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/0/U.air

@@ -0,0 +1,55 @@
+/*--------------------------------*- 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       volVectorField;
+    object      U.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    AMI1a
+    {
+        type            cyclicAMI;
+        value           uniform (0 0 0);
+    }
+
+    AMI1b
+    {
+        type            cyclicAMI;
+        value           uniform (0 0 0);
+    }
+
+    rotor
+    {
+
+        type            movingWallVelocity;
+        value           uniform (0 0 0);
+    }
+
+    ".*"
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/0/p

@@ -0,0 +1,54 @@
+/*--------------------------------*- 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       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    AMI1a
+    {
+        type            cyclicAMI;
+        value           $internalField;
+    }
+
+    AMI1b
+    {
+        type            cyclicAMI;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    ".*"
+    {
+        type            zeroGradient;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/constant/dynamicMeshDict

@@ -0,0 +1,34 @@
+/*--------------------------------*- 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;
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicMotionSolverFvMesh;
+
+motionSolverLibs (fvMotionSolvers);
+
+motionSolver    solidBody;
+
+cellZone        rotor;
+
+solidBodyMotionFunction  rotatingMotion;
+rotatingMotionCoeffs
+{
+    origin        (0.225 0.15 0);
+    axis          (0 0 1);
+    omega         5; // rad/s
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/constant/g

@@ -0,0 +1,21 @@
+/*--------------------------------*- 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       uniformDimensionedVectorField;
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 -9.81 0 );
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/constant/kinematicCloudProperties

@@ -0,0 +1,170 @@
+/*--------------------------------*- 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;
+    object      particleProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solution
+{
+    active          true;
+    coupled         false;
+    transient       yes;
+    cellValueSourceCorrection off;
+    maxCo           0.3;
+
+    averagingMethod dual;
+
+    sourceTerms
+    {
+        schemes
+        {
+            U       semiImplicit 1;
+        }
+    }
+
+    interpolationSchemes
+    {
+        rho             cell;
+        U.air           cellPoint;
+        mu.air          cell;
+        rho.air         cell;
+    }
+
+    integrationSchemes
+    {
+        U               Euler;
+    }
+}
+
+constantProperties
+{
+    rho0            1000;
+    alphaMax        0.9;
+}
+
+subModels
+{
+    particleForces
+    {
+        gravity;
+
+        ErgunWenYuDrag
+        {
+           alphac alpha.air;
+        }
+    }
+
+    injectionModels
+    {
+        model1
+        {
+            type            patchInjection;
+            parcelBasisType fixed;
+            patch           inlet;
+            U0              (0.3 0 0);
+            nParticle       1;
+            parcelsPerSecond  20000;
+
+            sizeDistribution
+	    {
+	        type uniform;
+                uniformDistribution
+                {
+                    minValue        5e-03;
+                    maxValue        5e-03;
+                }
+	    }
+
+            flowRateProfile constant 1;
+            massTotal       2000000;
+            SOI             0;
+            duration        2;
+        }
+    }
+
+    dispersionModel none;
+
+    patchInteractionModel none;
+
+    heatTransferModel none;
+
+    surfaceFilmModel none;
+
+    collisionModel none;
+
+    patchInteractionModel standardWallInteraction;
+
+    standardWallInteractionCoeffs
+    {
+        type        rebound;    //  stick, escape
+        e           0.6;        //  elasticity coeff
+        mu          0.09;       //  tan coeff
+        UrMax       1e-4;       //  relative U of particle after collision
+                                //  bellow which the particle is considered
+                                //  at the same U as the patch and deleted
+    }
+
+    stochasticCollisionModel none;
+
+    packingModel explicit;
+
+    explicitCoeffs
+    {
+        particleStressModel
+        {
+            type        HarrisCrighton;
+            alphaPacked 0.6;
+            pSolid      8.0;
+            beta        2;
+            eps         1.0e-7;
+        }
+        correctionLimitingMethod
+        {
+            type        absolute;
+            e           0.7;
+        }
+    }
+
+    dampingModel none;
+
+    relaxationCoeffs
+    {
+        timeScaleModel
+        {
+            type        isotropic;
+            alphaPacked 0.6;
+            e           0.88;
+        }
+    }
+
+    isotropyModel none;
+
+    stochasticCoeffs
+    {
+        timeScaleModel
+        {
+            type isotropic;
+            alphaPacked 0.6;
+            e 0.8;
+        }
+    }
+
+    radiation off;
+}
+
+
+cloudFunctions
+{}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/constant/transportProperties

@@ -0,0 +1,25 @@
+/*--------------------------------*- 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;
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+continuousPhase air;
+
+rho.air         1.2;
+
+transportModel  Newtonian;
+nu              1e-05;
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/constant/turbulenceProperties.air

@@ -0,0 +1,20 @@
+/*--------------------------------*- 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;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/system/blockMeshDict

@@ -0,0 +1,421 @@
+/*--------------------------------*- 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;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale    1;
+
+vertices
+(
+    (0 0 0)         //p0
+    (0 0.1 0)       //p1
+    (0 0.2 0)       //p2
+    (0 0.3 0)       //p3
+
+    (0.075 0 0)     //p4
+    (0.116747 0.0875 0) //p5
+    (0.116747 0.2125 0) //p6
+    (0.075 0.3 0)   //p7
+
+    (0.225 -0.2 0)  //p8
+    (0.225 0 0)     //p9
+    (0.225 0.025 0) //p10
+    (0.225 0.275 0) //p11
+    (0.225 0.3 0)   //p12
+    (0.225 0.5 0)   //p13
+
+    (0.375 -0.2 0)  //p14
+    (0.375 0 0)     //p15
+    (0.333253 0.0875 0) //p16
+    (0.333253 0.2125 0) //p17
+    (0.375 0.3 0)   //p18
+    (0.375 0.5 0)   //p19
+
+    (1.225 -0.2 0)      //p20
+    (1.225 0 0)         //p21
+    (1.225 0.0875 0)    //p22
+    (1.225 0.2125 0)    //p23
+    (1.225 0.3 0)       //p24
+    (1.225 0.5 0)       //p25
+
+ //domain back
+    (0 0 0.10) //p26
+    (0 0.1 0.10) //p27
+    (0 0.2 0.10) //p28
+    (0 0.3 0.10) //p29
+
+    (0.075 0 0.10) //p30
+    (0.116747 0.0875 0.10) //p31
+    (0.116747 0.2125 0.10) //p32
+    (0.075 0.3 0.10) //p33
+
+    (0.225 -0.2 0.10) //p34
+    (0.225 0 0.10) //p35
+    (0.225 0.025 0.10) //p36
+    (0.225 0.275 0.10) //p37
+    (0.225 0.3 0.10) //p38
+    (0.225 0.5 0.10) //p39
+
+    (0.375 -0.2 0.10) //p40
+    (0.375 0 0.10) //p41
+    (0.333253 0.0875 0.10) //p42
+    (0.333253 0.2125 0.10) //p43
+    (0.375 0.3 0.10) //p44
+    (0.375 0.5 0.10) //p45
+
+    (1.225 -0.2 0.10) //p46
+    (1.225 0 0.10) //p47
+    (1.225 0.0875 0.10) //p48
+    (1.225 0.2125 0.10) //p49
+    (1.225 0.3 0.10) //p50
+    (1.225 0.5 0.10) //p51
+
+//rotorfront
+    (0.192524 0.13125 0)  //p52
+    (0.121077 0.09 0) //p53
+    (0.116747 0.0875 0) //p54
+
+    (0.192524 0.16875 0) //p55
+    (0.121077 0.21 0) //p56
+    (0.116747 0.2125 0) //p57
+    (0.192524 0.16875 0) //p58
+    (0.121077 0.21 0) //p59
+    (0.116747 0.2125 0) //p60
+
+    (0.225 0.1875 0) //61
+    (0.225 0.27 0) //62
+    (0.225 0.275 0) //63
+    (0.225 0.1875 0) //64
+    (0.225 0.27 0) //65
+    (0.225 0.275 0) //66
+
+    (0.257476 0.16875 0) //67
+    (0.328923 0.21 0) //68
+    (0.333253 0.2125 0) //69
+    (0.257476 0.16875 0) //70
+    (0.328923 0.21 0) //71
+    (0.333253 0.2125 0) //72
+
+    (0.257476 0.13125 0) //73
+    (0.328923 0.09 0) //74
+    (0.333253 0.0875 0) //75
+    (0.257476 0.13125 0) //76
+    (0.328923 0.09 0) //77
+    (0.333253 0.0875 0) //78
+
+    (0.225 0.1125 0) //79
+    (0.225 0.03 0) //80
+    (0.225 0.025 0) //81
+    (0.225 0.1125 0) //82
+    (0.225 0.03 0) //83
+    (0.225 0.025 0) //84
+
+    (0.192524 0.13125 0)  //p85
+    (0.121077 0.09 0) //p86
+    (0.116747 0.0875 0) //p87
+
+
+//rotorback
+    (0.192524 0.13125 0.1)  //p88
+    (0.121077 0.09 0.1) //p89
+    (0.116747 0.0875 0.1) //p90
+
+    (0.192524 0.16875 0.1) //p91
+    (0.121077 0.21 0.1) //p92
+    (0.116747 0.2125 0.1) //p93
+    (0.192524 0.16875 0.1) //p94
+    (0.121077 0.21 0.1) //p95
+    (0.116747 0.2125 0.1) //p96
+
+    (0.225 0.1875 0.1) //97
+    (0.225 0.27 0.1) //98
+    (0.225 0.275 0.1) //99
+    (0.225 0.1875 0.1) //100
+    (0.225 0.27 0.1) //101
+    (0.225 0.275 0.1) //102
+
+    (0.257476 0.16875 0.1) //103
+    (0.328923 0.21 0.1) //104
+    (0.333253 0.2125 0.1) //105
+    (0.257476 0.16875 0.1) //106
+    (0.328923 0.21 0.1) //107
+    (0.333253 0.2125 0.1) //108
+
+    (0.257476 0.13125 0.1) //109
+    (0.328923 0.09 0.1) //110
+    (0.333253 0.0875 0.1) //111
+    (0.257476 0.13125 0.1) //112
+    (0.328923 0.09 0.1) //113
+    (0.333253 0.0875 0.1) //114
+
+    (0.225 0.1125 0.1) //115
+    (0.225 0.03 0.1) //116
+    (0.225 0.025 0.1) //117
+    (0.225 0.1125 0.1) //118
+    (0.225 0.03 0.1) //119
+    (0.225 0.025 0.1) //120
+
+    (0.192524 0.13125 0.1)  //p121
+    (0.121077 0.09 0.1) //p122
+    (0.116747 0.0875 0.1) //p123
+);
+
+
+blocks
+(
+//domain
+    hex (0 4 5 1 26 30 31 27) domain (5 5 1) simpleGrading (1 1 1) //b0
+    hex (1 5 6 2 27 31 32 28) domain (5 5 1) simpleGrading (1 1 1) //b1
+    hex (2 6 7 3 28 32 33 29) domain (5 5 1) simpleGrading (1 1 1) //b2
+
+    hex (4 9 10 5 30 35 36 31) domain (5 5 1) simpleGrading (1 1 1) //b3
+    hex (6 11 12 7 32 37 38 33) domain (5 5 1) simpleGrading (1 1 1) //b4
+
+    hex (8 14 15 9 34 40 41 35) domain (5 10 1) simpleGrading (1 1 1) //b5
+    hex (9 15 16 10 35 41 42 36) domain (5 5 1) simpleGrading (1 1 1) //b6
+    hex (11 17 18 12 37 43 44 38) domain (5 5 1) simpleGrading (1 1 1) //b7
+    hex (12 18 19 13 38 44 45 39) domain (5 10 1) simpleGrading (1 1 1) //b8
+
+    hex (14 20 21 15 40 46 47 41) domain (20 10 1) simpleGrading (1 1 1)//b9
+    hex (15 21 22 16 41 47 48 42) domain (20 5 1) simpleGrading (1 1 1)//b10
+    hex (16 22 23 17 42 48 49 43) domain (20 5 1) simpleGrading (1 1 1)//b11
+    hex (17 23 24 18 43 49 50 44) domain (20 5 1) simpleGrading (1 1 1)//b12
+    hex (18 24 25 19 44 50 51 45) domain (20 10 1) simpleGrading (1 1 1)//b13
+
+//rotor
+    hex (52 55 56 53 88 91 92 89) rotor (5 5 1) simpleGrading (1 1 1) //b14
+    hex (58 61 62 59 94 97 98 95) rotor (5 5 1) simpleGrading (1 1 1) //b15
+    hex (64 67 68 65 100 103 104 101) rotor (5 5 1) simpleGrading (1 1 1) //b16
+    hex (70 73 74 71 106 109 110 107) rotor (5 5 1) simpleGrading (1 1 1) //b17
+    hex (76 79 80 77 112 115 116 113) rotor (5 5 1) simpleGrading (1 1 1) //b18
+    hex (82 85 86 83 118 121 122 119) rotor (5 5 1) simpleGrading (1 1 1) //b19
+
+    hex (53 56 57 54 89 92 93 90) rotor (5 1 1) simpleGrading (1 1 1) //b20
+    hex (59 62 63 60 95 98 99 96) rotor (5 1 1) simpleGrading (1 1 1) //b21
+    hex (65 68 69 66 101 104 105 102) rotor (5 1 1) simpleGrading (1 1 1) //b22
+    hex (71 74 75 72 107 110 111 108) rotor (5 1 1) simpleGrading (1 1 1) //b23
+    hex (77 80 81 78 113 116 117 114) rotor (5 1 1) simpleGrading (1 1 1) //b24
+    hex (83 86 87 84 119 122 123 120) rotor (5 1 1) simpleGrading (1 1 1) //b25
+);
+
+edges
+(
+///// outer AMI
+    arc 10 5 (0.1625 0.041747 0)
+    arc 5 6 (0.1 0.15 0)
+    arc 6 11 (0.1625 0.258253 0)
+    arc 11 17 (0.2875 0.258253 0)
+    arc 17 16 (0.35 0.15 0)
+    arc 16 10 (0.2875 0.041747 0)
+
+    arc 36 31 (0.1625 0.041747 0.1)
+    arc 31 32 (0.1 0.15 0.1)
+    arc 32 37 (0.1625 0.258253 0.1)
+    arc 37 43 (0.2875 0.258253 0.1)
+    arc 43 42 (0.35 0.15 0.1)
+    arc 42 36 (0.2875 0.041747 0.1)
+
+//inner AMI
+    arc 84 87 (0.1625 0.041747 0)
+    arc 54 57 (0.1 0.15 0)
+    arc 60 63 (0.1625 0.258253 0)
+    arc 66 69 (0.2875 0.258253 0)
+    arc 72 75 (0.35 0.15 0)
+    arc 78 81 (0.2875 0.041747 0)
+
+    arc 120 123 (0.1625 0.041747 0.1)
+    arc 90 93 (0.1 0.15 0.1)
+    arc 96 99 (0.1625 0.258253 0.1)
+    arc 102 105 (0.2875 0.258253 0.1)
+    arc 108 111 (0.35 0.15 0.1)
+    arc 114 117 (0.2875 0.041747 0.1)
+
+//outer rotor
+
+    arc 83 86 (0.165 0.046077 0)
+    arc 53 56 (0.105 0.15 0)
+    arc 59 62 (0.165 0.253923 0)
+    arc 65 68 (0.285 0.253923 0)
+    arc 71 74 (0.345 0.15 0)
+    arc 77 80 (0.285 0.046077 0)
+
+    arc 119 122 (0.165 0.046077 0.1)
+    arc 89 92 (0.105 0.15 0.1)
+    arc 95 98 (0.165 0.253923 0.1)
+    arc 101 104 (0.285 0.253923 0.1)
+    arc 107 110 (0.345 0.15 0.1)
+    arc 113 116 (0.285 0.046077 0.1)
+
+//inner rotor
+
+    arc 82 85 (0.20625 0.117524 0)
+    arc 52 55 (0.1875 0.15 0)
+    arc 58 61 (0.20625 0.182476 0)
+    arc 64 67 (0.24375 0.182476 0)
+    arc 70 73 (0.2625 0.15 0)
+    arc 76 79 (0.24375 0.117524 0)
+
+    arc 118 121 (0.20625 0.117524 0.1)
+    arc 88 91 (0.1875 0.15 0.1)
+    arc 94 97 (0.20625 0.182476 0.1)
+    arc 100 103 (0.24375 0.182476 0.1)
+    arc 106 109 (0.2625 0.15 0.1)
+    arc 112 115 (0.24375 0.117524 0.1)
+);
+
+boundary
+(
+    walls
+    {
+        type wall;
+        faces
+        (
+            (0 4 30 26)
+            (4 9 35 30)
+            (9 8 34 35)
+            (3 7 33 29)
+            (7 12 38 33)
+            (12 13 39 38)
+        );
+    }
+
+    floor
+    {
+        type wall;
+        faces
+        (
+            (8 14 40 34)
+            (14 20 46 40)
+        );
+    }
+
+    roof
+    {
+        type wall;
+        faces
+        (
+            (13 19 45 39)
+            (19 25 51 45)
+        );
+    }
+
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 1 27 26)
+            (1 2 28 27)
+            (2 3 29 28)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (20 21 47 46)
+            (21 22 48 47)
+            (22 23 49 48)
+            (23 24 50 49)
+            (24 25 51 50)
+        );
+    }
+
+    AMI1a
+    {
+        type            cyclicAMI;
+        fraction        0.05;
+        matchTolerance  0.0001;
+        transform       noOrdering;
+        neighbourPatch  AMI1b;
+        faces
+        (
+            (10 5 31 36)
+            (5 6 32 31)
+            (6 11 37 32)
+            (11 17 43 37)
+            (17 16 42 43)
+            (16 10 36 42)
+        );
+    }
+
+    AMI1b
+    {
+        type            cyclicAMI;
+        fraction        0.05;
+        matchTolerance  0.0001;
+        transform       noOrdering;
+        neighbourPatch  AMI1a;
+        faces
+        (
+
+                (84 87 123 120)
+                (54 57 93 90)
+                (60 63 99 96)
+                (66 69 105 102)
+                (72 75 111 108)
+                (78 81 117 114)
+        );
+    }
+
+    rotor
+    {
+        type wall;
+        faces
+        (
+            (82 83 119 118)
+            (83 84 120 119)
+            (79 80 116 115)
+            (80 81 117 116)
+
+            (85 86 122 121)
+            (86 87 123 122)
+            (52 53 89 88)
+            (53 54 90 89)
+
+            (55 56 92 91)
+            (56 57 93 92)
+            (58 59 95 94)
+            (59 60 96 95)
+
+            (61 62 98 97)
+            (62 63 99 98)
+            (64 65 101 100)
+            (65 66 102 101)
+
+            (67 68 104 103)
+            (68 69 105 104)
+            (70 71 107 106)
+            (71 72 108 107)
+
+            (76 77 113 112)
+            (77 78 114 113)
+            (73 74 110 109)
+            (74 75 111 110)
+
+            (82 85 121 118)
+            (52 55 91 88)
+            (58 61 97 94)
+            (64 67 103 100)
+            (70 73 109 106)
+            (76 79 115 112)
+
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/system/controlDict

@@ -0,0 +1,53 @@
+/*--------------------------------*- 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;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DebugSwitches
+{
+    particle 0;
+}
+
+application     MPPICDyMFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2;
+
+deltaT          1e-3;
+
+writeControl    runTime;
+
+writeInterval   0.01;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable no;
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/system/fvSchemes

@@ -0,0 +1,58 @@
+/*--------------------------------*- 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;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+d2dt2Schemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+
+    div(alphaPhi.air,U.air)    Gauss linearUpwindV limited 1;
+    div(((alpha.air*nuEff.air)*dev2(T(grad(U.air))))) Gauss linear;
+
+    div(phiGByA,kinematicCloud:alpha) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/denseRotor2DAMI/system/fvSolution

@@ -0,0 +1,84 @@
+/*--------------------------------*- 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;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "(p|pcorr)"
+    {
+        solver          GAMG;
+        tolerance       1e-7;
+        relTol          0.01;
+        smoother        GaussSeidel;
+        maxIter         20;
+    }
+
+    "(pFinal|pcorrFinal)"
+    {
+        $p;
+        relTol          0;
+    }
+
+    U.air
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    U.airFinal
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    kinematicCloud:alpha
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-06;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nOuterCorrectors 1;
+    nCorrectors     3;
+    momentumPredictor no;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+    checkMeshCourantNo  no;
+    correctPhi          no;
+}
+
+relaxationFactors
+{
+    fields
+    {
+        ".*"            1;
+    }
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/0.orig/T

@@ -0,0 +1,48 @@
+/*--------------------------------*- 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       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 350;
+
+boundaryField
+{
+    AMI1a
+    {
+        type            cyclicAMI;
+        value           $internalField;
+    }
+
+    AMI1b
+    {
+        type            cyclicAMI;
+        value           $internalField;
+    }
+
+    ".*"
+    {
+        type            zeroGradient;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/0.orig/U

@@ -0,0 +1,54 @@
+/*--------------------------------*- 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       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    AMI1a
+    {
+        type            cyclicAMI;
+        value           uniform (0 0 0);
+    }
+
+    AMI1b
+    {
+        type            cyclicAMI;
+        value           uniform (0 0 0);
+    }
+
+    rotor
+    {
+        type            movingWallVelocity;
+        value           uniform (0 0 0);
+    }
+
+    ".*"
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/0.orig/p

@@ -0,0 +1,48 @@
+/*--------------------------------*- 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       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    AMI1a
+    {
+        type            cyclicAMI;
+        value           $internalField;
+    }
+
+    AMI1b
+    {
+        type            cyclicAMI;
+        value           $internalField;
+    }
+
+    ".*"
+    {
+        type            zeroGradient;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/Allclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase0
+
+#------------------------------------------------------------------------------

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/Allrun

@@ -0,0 +1,12 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+restore0Dir
+
+runApplication blockMesh
+
+runApplication $(getApplication)
+
+#------------------------------------------------------------------------------

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/constant/dynamicMeshDict

@@ -0,0 +1,34 @@
+/*--------------------------------*- 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;
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicMotionSolverFvMesh;
+
+motionSolverLibs (fvMotionSolvers);
+
+motionSolver    solidBody;
+
+cellZone        rotor;
+
+solidBodyMotionFunction  rotatingMotion;
+rotatingMotionCoeffs
+{
+    origin        (0.225 0.15 0);
+    axis          (0 0 1);
+    omega         5; // rad/s
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/constant/g

@@ -0,0 +1,21 @@
+/*--------------------------------*- 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       uniformDimensionedVectorField;
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           ( 0 -9.81 0 );
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/constant/kinematicCloudProperties

@@ -0,0 +1,119 @@
+/*--------------------------------*- 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;
+    object      particleProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solution
+{
+    active          true;
+    coupled         false;
+    transient       yes;
+    cellValueSourceCorrection off;
+    maxCo           0.3;
+
+    sourceTerms
+    {
+        schemes
+        {
+        }
+    }
+
+    interpolationSchemes
+    {
+        rho             cell;
+        U               cellPoint;
+        thermo:mu       cell;
+    }
+
+    integrationSchemes
+    {
+        U               Euler;
+    }
+}
+
+constantProperties
+{
+
+    rho0            1000;
+    constantVolume  true;
+
+}
+
+subModels
+{
+    particleForces
+    {
+        //sphereDrag;
+        //gravity;
+    }
+
+    injectionModels
+    {
+         model1
+        {
+            type            patchInjection;
+            parcelBasisType fixed;
+            patch           inlet;
+            U0              (0.1 0 0);
+            nParticle       1;
+            parcelsPerSecond  20000000;
+
+            sizeDistribution
+            {
+                type uniform;
+                uniformDistribution
+                {
+                    minValue        50e-06;
+                    maxValue        50e-06;
+                }
+            }
+
+            flowRateProfile constant 1;
+            massTotal       2000000;
+            SOI             0;
+            duration        0.00001;
+        }
+    }
+
+    dispersionModel none;
+
+    patchInteractionModel none;
+
+    heatTransferModel none;
+
+    surfaceFilmModel none;
+
+    collisionModel none;
+
+    patchInteractionModel standardWallInteraction;
+
+    standardWallInteractionCoeffs
+    {
+        type        rebound;    // stick, escape
+        e           0.9;        // optional - elasticity coeff
+        mu          0.09;       // optional - restitution coeff
+        UrMax       1e-4;
+    }
+
+    stochasticCollisionModel none;
+
+    radiation off;
+}
+
+
+cloudFunctions
+{}
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/constant/thermophysicalProperties

@@ -0,0 +1,47 @@
+/*--------------------------------*- 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;
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            hePsiThermo;
+    mixture         pureMixture;
+    transport       const;
+    thermo          eConst;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    specie
+    {
+        molWeight   28.9;
+    }
+    thermodynamics
+    {
+        Cv          712;
+        Hf          0;
+    }
+    transport
+    {
+        mu          1.8e-05;
+        Pr          0.7;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/constant/turbulenceProperties

@@ -0,0 +1,20 @@
+/*--------------------------------*- 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;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  laminar;
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/system/blockMeshDict

@@ -0,0 +1,420 @@
+/*--------------------------------*- 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;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale    1;
+
+vertices
+(
+    (0 0 0)         //p0
+    (0 0.1 0)       //p1
+    (0 0.2 0)       //p2
+    (0 0.3 0)       //p3
+
+    (0.075 0 0)     //p4
+    (0.116747 0.0875 0) //p5
+    (0.116747 0.2125 0) //p6
+    (0.075 0.3 0)   //p7
+
+    (0.225 -0.2 0)  //p8
+    (0.225 0 0)     //p9
+    (0.225 0.025 0) //p10
+    (0.225 0.275 0) //p11
+    (0.225 0.3 0)   //p12
+    (0.225 0.5 0)   //p13
+
+    (0.375 -0.2 0)  //p14
+    (0.375 0 0)     //p15
+    (0.333253 0.0875 0) //p16
+    (0.333253 0.2125 0) //p17
+    (0.375 0.3 0)   //p18
+    (0.375 0.5 0)   //p19
+
+    (1.225 -0.2 0)      //p20
+    (1.225 0 0)         //p21
+    (1.225 0.0875 0)    //p22
+    (1.225 0.2125 0)    //p23
+    (1.225 0.3 0)       //p24
+    (1.225 0.5 0)       //p25
+
+ //domain back
+    (0 0 0.10) //p26
+    (0 0.1 0.10) //p27
+    (0 0.2 0.10) //p28
+    (0 0.3 0.10) //p29
+
+    (0.075 0 0.10) //p30
+    (0.116747 0.0875 0.10) //p31
+    (0.116747 0.2125 0.10) //p32
+    (0.075 0.3 0.10) //p33
+
+    (0.225 -0.2 0.10) //p34
+    (0.225 0 0.10) //p35
+    (0.225 0.025 0.10) //p36
+    (0.225 0.275 0.10) //p37
+    (0.225 0.3 0.10) //p38
+    (0.225 0.5 0.10) //p39
+
+    (0.375 -0.2 0.10) //p40
+    (0.375 0 0.10) //p41
+    (0.333253 0.0875 0.10) //p42
+    (0.333253 0.2125 0.10) //p43
+    (0.375 0.3 0.10) //p44
+    (0.375 0.5 0.10) //p45
+
+    (1.225 -0.2 0.10) //p46
+    (1.225 0 0.10) //p47
+    (1.225 0.0875 0.10) //p48
+    (1.225 0.2125 0.10) //p49
+    (1.225 0.3 0.10) //p50
+    (1.225 0.5 0.10) //p51
+
+//rotorfront
+
+    (0.192524 0.13125 0)  //p52
+    (0.121077 0.09 0) //p53
+    (0.116747 0.0875 0) //p54
+
+    (0.192524 0.16875 0) //p55
+    (0.121077 0.21 0) //p56
+    (0.116747 0.2125 0) //p57
+    (0.192524 0.16875 0) //p58
+    (0.121077 0.21 0) //p59
+    (0.116747 0.2125 0) //p60
+
+    (0.225 0.1875 0) //61
+    (0.225 0.27 0) //62
+    (0.225 0.275 0) //63
+    (0.225 0.1875 0) //64
+    (0.225 0.27 0) //65
+    (0.225 0.275 0) //66
+
+    (0.257476 0.16875 0) //67
+    (0.328923 0.21 0) //68
+    (0.333253 0.2125 0) //69
+    (0.257476 0.16875 0) //70
+    (0.328923 0.21 0) //71
+    (0.333253 0.2125 0) //72
+
+    (0.257476 0.13125 0) //73
+    (0.328923 0.09 0) //74
+    (0.333253 0.0875 0) //75
+    (0.257476 0.13125 0) //76
+    (0.328923 0.09 0) //77
+    (0.333253 0.0875 0) //78
+
+    (0.225 0.1125 0) //79
+    (0.225 0.03 0) //80
+    (0.225 0.025 0) //81
+    (0.225 0.1125 0) //82
+    (0.225 0.03 0) //83
+    (0.225 0.025 0) //84
+
+    (0.192524 0.13125 0)  //p85
+    (0.121077 0.09 0) //p86
+    (0.116747 0.0875 0) //p87
+
+
+//rotorback
+
+    (0.192524 0.13125 0.1)  //p88
+    (0.121077 0.09 0.1) //p89
+    (0.116747 0.0875 0.1) //p90
+
+    (0.192524 0.16875 0.1) //p91
+    (0.121077 0.21 0.1) //p92
+    (0.116747 0.2125 0.1) //p93
+    (0.192524 0.16875 0.1) //p94
+    (0.121077 0.21 0.1) //p95
+    (0.116747 0.2125 0.1) //p96
+
+    (0.225 0.1875 0.1) //97
+    (0.225 0.27 0.1) //98
+    (0.225 0.275 0.1) //99
+    (0.225 0.1875 0.1) //100
+    (0.225 0.27 0.1) //101
+    (0.225 0.275 0.1) //102
+
+    (0.257476 0.16875 0.1) //103
+    (0.328923 0.21 0.1) //104
+    (0.333253 0.2125 0.1) //105
+    (0.257476 0.16875 0.1) //106
+    (0.328923 0.21 0.1) //107
+    (0.333253 0.2125 0.1) //108
+
+    (0.257476 0.13125 0.1) //109
+    (0.328923 0.09 0.1) //110
+    (0.333253 0.0875 0.1) //111
+    (0.257476 0.13125 0.1) //112
+    (0.328923 0.09 0.1) //113
+    (0.333253 0.0875 0.1) //114
+
+    (0.225 0.1125 0.1) //115
+    (0.225 0.03 0.1) //116
+    (0.225 0.025 0.1) //117
+    (0.225 0.1125 0.1) //118
+    (0.225 0.03 0.1) //119
+    (0.225 0.025 0.1) //120
+
+    (0.192524 0.13125 0.1)  //p121
+    (0.121077 0.09 0.1) //p122
+    (0.116747 0.0875 0.1) //p123
+);
+
+blocks
+(
+//domain
+    hex (0 4 5 1 26 30 31 27) domain (5 5 1) simpleGrading (1 1 1) //b0
+    hex (1 5 6 2 27 31 32 28) domain (5 5 1) simpleGrading (1 1 1) //b1
+    hex (2 6 7 3 28 32 33 29) domain (5 5 1) simpleGrading (1 1 1) //b2
+
+    hex (4 9 10 5 30 35 36 31) domain (5 5 1) simpleGrading (1 1 1) //b3
+    hex (6 11 12 7 32 37 38 33) domain (5 5 1) simpleGrading (1 1 1) //b4
+
+    hex (8 14 15 9 34 40 41 35) domain (5 10 1) simpleGrading (1 1 1) //b5
+    hex (9 15 16 10 35 41 42 36) domain (5 5 1) simpleGrading (1 1 1) //b6
+    hex (11 17 18 12 37 43 44 38) domain (5 5 1) simpleGrading (1 1 1) //b7
+    hex (12 18 19 13 38 44 45 39) domain (5 10 1) simpleGrading (1 1 1) //b8
+
+    hex (14 20 21 15 40 46 47 41) domain (20 10 1) simpleGrading (1 1 1)//b9
+    hex (15 21 22 16 41 47 48 42) domain (20 5 1) simpleGrading (1 1 1)//b10
+    hex (16 22 23 17 42 48 49 43) domain (20 5 1) simpleGrading (1 1 1)//b11
+    hex (17 23 24 18 43 49 50 44) domain (20 5 1) simpleGrading (1 1 1)//b12
+    hex (18 24 25 19 44 50 51 45) domain (20 10 1) simpleGrading (1 1 1)//b13
+
+//rotor
+    hex (52 55 56 53 88 91 92 89) rotor (5 5 1) simpleGrading (1 1 1) //b14
+    hex (58 61 62 59 94 97 98 95) rotor (5 5 1) simpleGrading (1 1 1) //b15
+    hex (64 67 68 65 100 103 104 101) rotor (5 5 1) simpleGrading (1 1 1) //b16
+    hex (70 73 74 71 106 109 110 107) rotor (5 5 1) simpleGrading (1 1 1) //b17
+    hex (76 79 80 77 112 115 116 113) rotor (5 5 1) simpleGrading (1 1 1) //b18
+    hex (82 85 86 83 118 121 122 119) rotor (5 5 1) simpleGrading (1 1 1) //b19
+
+    hex (53 56 57 54 89 92 93 90) rotor (5 1 1) simpleGrading (1 1 1) //b20
+    hex (59 62 63 60 95 98 99 96) rotor (5 1 1) simpleGrading (1 1 1) //b21
+    hex (65 68 69 66 101 104 105 102) rotor (5 1 1) simpleGrading (1 1 1) //b22
+    hex (71 74 75 72 107 110 111 108) rotor (5 1 1) simpleGrading (1 1 1) //b23
+    hex (77 80 81 78 113 116 117 114) rotor (5 1 1) simpleGrading (1 1 1) //b24
+    hex (83 86 87 84 119 122 123 120) rotor (5 1 1) simpleGrading (1 1 1) //b25
+);
+
+edges
+(
+///// outer AMI
+    arc 10 5 (0.1625 0.041747 0)
+    arc 5 6 (0.1 0.15 0)
+    arc 6 11 (0.1625 0.258253 0)
+    arc 11 17 (0.2875 0.258253 0)
+    arc 17 16 (0.35 0.15 0)
+    arc 16 10 (0.2875 0.041747 0)
+
+    arc 36 31 (0.1625 0.041747 0.1)
+    arc 31 32 (0.1 0.15 0.1)
+    arc 32 37 (0.1625 0.258253 0.1)
+    arc 37 43 (0.2875 0.258253 0.1)
+    arc 43 42 (0.35 0.15 0.1)
+    arc 42 36 (0.2875 0.041747 0.1)
+
+//inner AMI
+    arc 84 87 (0.1625 0.041747 0)
+    arc 54 57 (0.1 0.15 0)
+    arc 60 63 (0.1625 0.258253 0)
+    arc 66 69 (0.2875 0.258253 0)
+    arc 72 75 (0.35 0.15 0)
+    arc 78 81 (0.2875 0.041747 0)
+
+    arc 120 123 (0.1625 0.041747 0.1)
+    arc 90 93 (0.1 0.15 0.1)
+    arc 96 99 (0.1625 0.258253 0.1)
+    arc 102 105 (0.2875 0.258253 0.1)
+    arc 108 111 (0.35 0.15 0.1)
+    arc 114 117 (0.2875 0.041747 0.1)
+
+//outer rotor
+
+    arc 83 86 (0.165 0.046077 0)
+    arc 53 56 (0.105 0.15 0)
+    arc 59 62 (0.165 0.253923 0)
+    arc 65 68 (0.285 0.253923 0)
+    arc 71 74 (0.345 0.15 0)
+    arc 77 80 (0.285 0.046077 0)
+
+    arc 119 122 (0.165 0.046077 0.1)
+    arc 89 92 (0.105 0.15 0.1)
+    arc 95 98 (0.165 0.253923 0.1)
+    arc 101 104 (0.285 0.253923 0.1)
+    arc 107 110 (0.345 0.15 0.1)
+    arc 113 116 (0.285 0.046077 0.1)
+
+//inner rotor
+
+    arc 82 85 (0.20625 0.117524 0)
+    arc 52 55 (0.1875 0.15 0)
+    arc 58 61 (0.20625 0.182476 0)
+    arc 64 67 (0.24375 0.182476 0)
+    arc 70 73 (0.2625 0.15 0)
+    arc 76 79 (0.24375 0.117524 0)
+
+    arc 118 121 (0.20625 0.117524 0.1)
+    arc 88 91 (0.1875 0.15 0.1)
+    arc 94 97 (0.20625 0.182476 0.1)
+    arc 100 103 (0.24375 0.182476 0.1)
+    arc 106 109 (0.2625 0.15 0.1)
+    arc 112 115 (0.24375 0.117524 0.1)
+);
+
+boundary
+(
+    walls
+    {
+        type wall;
+        faces
+        (
+            (0 4 30 26)
+            (4 9 35 30)
+            (9 8 34 35)
+            (3 7 33 29)
+            (7 12 38 33)
+            (12 13 39 38)
+        );
+    }
+
+    floor
+    {
+        type wall;
+        faces
+        (
+            (8 14 40 34)
+            (14 20 46 40)
+        );
+    }
+
+    roof
+    {
+        type wall;
+        faces
+        (
+            (13 19 45 39)
+            (19 25 51 45)
+        );
+    }
+
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 1 27 26)
+            (1 2 28 27)
+            (2 3 29 28)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (20 21 47 46)
+            (21 22 48 47)
+            (22 23 49 48)
+            (23 24 50 49)
+            (24 25 51 50)
+        );
+    }
+
+    AMI1a
+    {
+        type            cyclicAMI;
+        matchTolerance  0.0001;
+        transform       noOrdering;
+        neighbourPatch  AMI1b;
+        faces
+        (
+            (10 5 31 36)
+            (5 6 32 31)
+            (6 11 37 32)
+            (11 17 43 37)
+            (17 16 42 43)
+            (16 10 36 42)
+        );
+    }
+
+    AMI1b
+    {
+        type            cyclicAMI;
+        matchTolerance  0.0001;
+        transform       noOrdering;
+        neighbourPatch  AMI1a;
+        faces
+        (
+
+                (84 87 123 120)
+                (54 57 93 90)
+                (60 63 99 96)
+                (66 69 105 102)
+                (72 75 111 108)
+                (78 81 117 114)
+        );
+    }
+
+    rotor
+    {
+        type wall;
+        faces
+        (
+            (82 83 119 118)
+            (83 84 120 119)
+            (79 80 116 115)
+            (80 81 117 116)
+
+            (85 86 122 121)
+            (86 87 123 122)
+            (52 53 89 88)
+            (53 54 90 89)
+
+            (55 56 92 91)
+            (56 57 93 92)
+            (58 59 95 94)
+            (59 60 96 95)
+
+            (61 62 98 97)
+            (62 63 99 98)
+            (64 65 101 100)
+            (65 66 102 101)
+
+            (67 68 104 103)
+            (68 69 105 104)
+            (70 71 107 106)
+            (71 72 108 107)
+
+            (76 77 113 112)
+            (77 78 114 113)
+            (73 74 110 109)
+            (74 75 111 110)
+
+            (82 85 121 118)
+            (52 55 91 88)
+            (58 61 97 94)
+            (64 67 103 100)
+            (70 73 109 106)
+            (76 79 115 112)
+
+        );
+    }
+);
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/system/controlDict

@@ -0,0 +1,53 @@
+/*--------------------------------*- 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;
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+DebugSwitches
+{
+    particle 0;
+}
+
+application     uncoupledKinematicParcelDyMFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         2;
+
+deltaT          1e-3;
+
+writeControl    runTime;
+
+writeInterval   0.01;
+
+purgeWrite      0;
+
+writeFormat     binary;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable yes;
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/system/fvSchemes

@@ -0,0 +1,43 @@
+/*--------------------------------*- 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;
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         none;
+}
+
+divSchemes
+{
+    default         none;
+}
+
+laplacianSchemes
+{
+    default         none;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/uncoupledKinematicParcelDyMFoam/rotor2DAMI/system/fvSolution

@@ -0,0 +1,21 @@
+/*--------------------------------*- 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;
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{}
+
+
+// ************************************************************************* //