ENH: Adding oldCellCentres field to polyMesh

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -309,7 +309,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 +515,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 +673,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 +1126,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 +1146,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 +1195,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 +1296,7 @@ void Foam::polyMesh::resetMotion() const
 {
     curMotionTimeIndex_ = 0;
     oldPointsPtr_.clear();
+    oldCellCentresPtr_.clear();
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyMesh.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -80,10 +80,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 +113,7 @@ public:
 
 private:
 
-    // Permanent data
+    // Private Data
 
         // Primitive mesh data
 
@@ -189,6 +188,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 +437,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

@@ -576,7 +576,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 +858,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

@@ -119,6 +119,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);