ENH: fvMesh: map fields using old volume

applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C

@@ -555,7 +555,8 @@ autoPtr<mapPolyMesh> reorderMesh
             labelListList(0),           // cellZoneMap,
             pointField(0),              // preMotionPoints,
             patchStarts,                // oldPatchStarts,
-            oldPatchNMeshPoints         // oldPatchNMeshPoints
+            oldPatchNMeshPoints,        // oldPatchNMeshPoints
+            autoPtr<scalarField>()      // oldCellVolumes
         )
     );
 }

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/cellMapper/cellMapper.C

@@ -149,20 +149,10 @@ void Foam::cellMapper::calcAddressing() const
             w[cellI] = scalarList(mo.size(), 1.0/mo.size());
         }
 
+        // Volume conservative mapping if possible
+
         const List<objectMap>& cfc = mpm_.cellsFromCellsMap();
 
-        const scalarField& V = mesh_.cellVolumes();
-
-        if (V.size() != sizeBeforeMapping())
-        {
-            FatalErrorIn("void cellMapper::calcAddressing() const")
-                << "cellVolumes size " << V.size()
-                << " is not the old number of cells " << sizeBeforeMapping()
-                << ". Are your cellVolumes already mapped?"
-                << " (new number of cells " << size() << ")"
-                << abort(FatalError);
-        }
-
         forAll(cfc, cfcI)
         {
             // Get addressing
@@ -175,16 +165,36 @@ void Foam::cellMapper::calcAddressing() const
                 FatalErrorIn("void cellMapper::calcAddressing() const")
                     << "Master cell " << cellI
                     << " mapped from cell cells " << mo
-                    << " already destination of mapping." << abort(FatalError);
+                    << " already destination of mapping."
+                    << abort(FatalError);
             }
 
             // Map from masters
             addr[cellI] = mo;
+        }
 
-            //- uniform weights
+        if (mpm_.hasOldCellVolumes())
-            //w[cellI] = scalarList(mo.size(), 1.0/mo.size());
+        {
+            // Volume weighted
+
+            const scalarField& V = mpm_.oldCellVolumes();
+
+            if (V.size() != sizeBeforeMapping())
+            {
+                FatalErrorIn("void cellMapper::calcAddressing() const")
+                    << "cellVolumes size " << V.size()
+                    << " is not the old number of cells " << sizeBeforeMapping()
+                    << ". Are your cellVolumes already mapped?"
+                    << " (new number of cells " << size() << ")"
+                    << abort(FatalError);
+            }
+
+            forAll(cfc, cfcI)
+            {
+                const labelList& mo = cfc[cfcI].masterObjects();
+
+                label cellI = cfc[cfcI].index();
 
-            //- volume based
                 w[cellI].setSize(mo.size());
 
                 if (mo.size())
@@ -195,11 +205,33 @@ void Foam::cellMapper::calcAddressing() const
                         w[cellI][ci] = V[mo[ci]];
                         sumV += V[mo[ci]];
                     }
+                    if (sumV > VSMALL)
+                    {
                         forAll(mo, ci)
                         {
                             w[cellI][ci] /= sumV;
                         }
                     }
+                    else
+                    {
+                        // Exception: zero volume. Use uniform mapping
+                        w[cellI] = scalarList(mo.size(), 1.0/mo.size());
+                    }
+                }
+            }
+        }
+        else
+        {
+            // Uniform weighted
+
+            forAll(cfc, cfcI)
+            {
+                const labelList& mo = cfc[cfcI].masterObjects();
+
+                label cellI = cfc[cfcI].index();
+
+                w[cellI] = scalarList(mo.size(), 1.0/mo.size());
+            }
         }
 
 

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.C

@@ -56,7 +56,8 @@ Foam::mapPolyMesh::mapPolyMesh
     const labelListList& cellZoneMap,
     const pointField& preMotionPoints,
     const labelList& oldPatchStarts,
-    const labelList& oldPatchNMeshPoints
+    const labelList& oldPatchNMeshPoints,
+    const autoPtr<scalarField>& oldCellVolumesPtr
 )
 :
     mesh_(mesh),
@@ -86,7 +87,8 @@ Foam::mapPolyMesh::mapPolyMesh
     preMotionPoints_(preMotionPoints),
     oldPatchSizes_(oldPatchStarts.size()),
     oldPatchStarts_(oldPatchStarts),
-    oldPatchNMeshPoints_(oldPatchNMeshPoints)
+    oldPatchNMeshPoints_(oldPatchNMeshPoints),
+    oldCellVolumesPtr_(oldCellVolumesPtr)
 {
     // Calculate old patch sizes
     for (label patchI = 0; patchI < oldPatchStarts_.size() - 1; patchI++)
@@ -141,6 +143,7 @@ Foam::mapPolyMesh::mapPolyMesh
     pointField& preMotionPoints,
     labelList& oldPatchStarts,
     labelList& oldPatchNMeshPoints,
+    autoPtr<scalarField>& oldCellVolumesPtr,
     const bool reUse
 )
 :
@@ -171,7 +174,8 @@ Foam::mapPolyMesh::mapPolyMesh
     preMotionPoints_(preMotionPoints, reUse),
     oldPatchSizes_(oldPatchStarts.size()),
     oldPatchStarts_(oldPatchStarts, reUse),
-    oldPatchNMeshPoints_(oldPatchNMeshPoints, reUse)
+    oldPatchNMeshPoints_(oldPatchNMeshPoints, reUse),
+    oldCellVolumesPtr_(oldCellVolumesPtr, reUse)
 {
     if (oldPatchStarts_.size() > 0)
     {

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -268,6 +268,9 @@ class mapPolyMesh
         //- List of numbers of mesh points per old patch
         const labelList oldPatchNMeshPoints_;
 
+        //- Optional old cell volumes (for mapping)
+        autoPtr<scalarField> oldCellVolumesPtr_;
+
 
     // Private Member Functions
 
@@ -282,7 +285,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from components. Copy (except for oldCellVolumes).
         mapPolyMesh
         (
             const polyMesh& mesh,
@@ -311,7 +314,8 @@ public:
             const labelListList& cellZoneMap,
             const pointField& preMotionPoints,
             const labelList& oldPatchStarts,
-            const labelList& oldPatchNMeshPoints
+            const labelList& oldPatchNMeshPoints,
+            const autoPtr<scalarField>& oldCellVolumesPtr
         );
 
         //- Construct from components and optionally reuse storage
@@ -344,6 +348,7 @@ public:
             pointField& preMotionPoints,
             labelList& oldPatchStarts,
             labelList& oldPatchNMeshPoints,
+            autoPtr<scalarField>& oldCellVolumesPtr,
             const bool reUse
         );
 
@@ -631,6 +636,20 @@ public:
             {
                 return oldPatchNMeshPoints_;
             }
+
+
+            // Geometric mapping data
+
+                bool hasOldCellVolumes() const
+                {
+                    return oldCellVolumesPtr_.valid();
+                }
+
+                const scalarField& oldCellVolumes() const
+                {
+                    return oldCellVolumesPtr_();
+                }
+
 };
 
 

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C

@@ -3145,6 +3145,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
     const label nOldPoints(mesh.nPoints());
     const label nOldFaces(mesh.nFaces());
     const label nOldCells(mesh.nCells());
+    autoPtr<scalarField> oldCellVolumes(new scalarField(mesh.cellVolumes()));
 
 
     // Change the mesh
@@ -3323,6 +3324,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
             newPoints,          // if empty signals no inflation.
             oldPatchStarts,
             oldPatchNMeshPoints,
+
+            oldCellVolumes,
+
             true                // steal storage.
         )
     );
@@ -3408,6 +3412,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
     const label nOldPoints(mesh.nPoints());
     const label nOldFaces(mesh.nFaces());
     const label nOldCells(mesh.nCells());
+    autoPtr<scalarField> oldCellVolumes(new scalarField(mesh.cellVolumes()));
 
 
     // Create the mesh
@@ -3617,6 +3622,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
             newPoints,          // if empty signals no inflation.
             oldPatchStarts,
             oldPatchNMeshPoints,
+            oldCellVolumes,
             true                // steal storage.
         )
     );

src/finiteVolume/fvMesh/fvMesh.C

@@ -118,6 +118,73 @@ void Foam::fvMesh::clearAddressing()
 }
 
 
+void Foam::fvMesh::storeOldVol(const scalarField& V)
+{
+    if (curTimeIndex_ < time().timeIndex())
+    {
+        if (debug)
+        {
+            Info<< "fvMesh::storeOldVol(const scalarField&) :"
+                << " Storing old time volumes since from time " << curTimeIndex_
+                << " and time now " << time().timeIndex()
+                << " V:" << V.size()
+                << endl;
+        }
+
+
+        if (V00Ptr_ && V0Ptr_)
+        {
+            // Copy V0 into V00 storage
+            *V00Ptr_ = *V0Ptr_;
+        }
+
+        if (V0Ptr_)
+        {
+            // Copy V into V0 storage
+            V0Ptr_->scalarField::operator=(V);
+        }
+        else
+        {
+            // Allocate V0 storage, fill with V
+            V0Ptr_ = new DimensionedField<scalar, volMesh>
+            (
+                IOobject
+                (
+                    "V0",
+                    time().timeName(),
+                    *this,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                *this,
+                dimVolume
+            );
+            scalarField& V0 = *V0Ptr_;
+            // Note: V0 now sized with current mesh, not with (potentially
+            //       different size) V.
+            V0.setSize(V.size());
+            V0 = V;
+        }
+
+        curTimeIndex_ = time().timeIndex();
+
+        if (debug)
+        {
+            Info<< "fvMesh::storeOldVol() :"
+                << " Stored old time volumes V0:" << V0Ptr_->size()
+                << endl;
+            if (V00Ptr_)
+            {
+                Info<< "fvMesh::storeOldVol() :"
+                    << " Stored oldold time volumes V00:" << V00Ptr_->size()
+                    << endl;
+            }
+        }
+    }
+}
+
+
 void Foam::fvMesh::clearOut()
 {
     clearGeom();
@@ -458,6 +525,35 @@ const Foam::lduAddressing& Foam::fvMesh::lduAddr() const
 
 void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
 {
+    if (debug)
+    {
+        Info<< "fvMesh::mapFields :"
+            << " nOldCells:" << meshMap.nOldCells()
+            << " nCells:" << nCells()
+            << " nOldFaces:" << meshMap.nOldFaces()
+            << " nFaces:" << nFaces()
+            << endl;
+    }
+
+
+    // We require geometric properties valid for the old mesh
+    if
+    (
+        meshMap.cellMap().size() != nCells()
+     || meshMap.faceMap().size() != nFaces()
+    )
+    {
+        FatalErrorIn("fvMesh::mapFields(const mapPolyMesh&)")
+            << "mapPolyMesh does not correspond to the old mesh."
+            << " nCells:" << nCells()
+            << " cellMap:" << meshMap.cellMap().size()
+            << " nOldCells:" << meshMap.nOldCells()
+            << " nFaces:" << nFaces()
+            << " faceMap:" << meshMap.faceMap().size()
+            << " nOldFaces:" << meshMap.nOldFaces()
+            << exit(FatalError);
+    }
+
     // Create a mapper
     const fvMeshMapper mapper(*this, meshMap);
 
@@ -517,7 +613,30 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
                 V0[i] = 0.0;
             }
         }
+
+        // Inject volume of merged cells
+        label nMerged = 0;
+        forAll(meshMap.reverseCellMap(), oldCellI)
+        {
+            label index = meshMap.reverseCellMap()[oldCellI];
+
+            if (index < -1)
+            {
+                label cellI = -index-2;
+
+                V0[cellI] += savedV0[oldCellI];
+
+                nMerged++;
             }
+        }
+
+        if (debug)
+        {
+            Info<< "Mapping old time volume V0. Merged "
+                << nMerged << " out of " << nCells() << " cells" << endl;
+        }
+    }
+
 
     // Map the old-old volume. Just map to new cell labels.
     if (V00Ptr_)
@@ -538,6 +657,27 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
                 V00[i] = 0.0;
             }
         }
+
+        // Inject volume of merged cells
+        label nMerged = 0;
+        forAll(meshMap.reverseCellMap(), oldCellI)
+        {
+            label index = meshMap.reverseCellMap()[oldCellI];
+
+            if (index < -1)
+            {
+                label cellI = -index-2;
+
+                V00[cellI] += savedV00[oldCellI];
+                nMerged++;
+            }
+        }
+
+        if (debug)
+        {
+            Info<< "Mapping old time volume V00. Merged "
+                << nMerged << " out of " << nCells() << " cells" << endl;
+        }
     }
 }
 
@@ -545,38 +685,12 @@ void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
 Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
 {
     // Grab old time volumes if the time has been incremented
+    // This will update V0, V00
     if (curTimeIndex_ < time().timeIndex())
     {
-        if (V00Ptr_ && V0Ptr_)
+        storeOldVol(V());
-        {
-            *V00Ptr_ = *V0Ptr_;
     }
 
-        if (V0Ptr_)
-        {
-            *V0Ptr_ = V();
-        }
-        else
-        {
-            V0Ptr_ = new DimensionedField<scalar, volMesh>
-            (
-                IOobject
-                (
-                    "V0",
-                    time().timeName(),
-                    *this,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE,
-                    false
-                ),
-                V()
-            );
-        }
-
-        curTimeIndex_ = time().timeIndex();
-    }
-
-
     if (!phiPtr_)
     {
         // Create mesh motion flux
@@ -649,10 +763,44 @@ void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
     // Update polyMesh. This needs to keep volume existent!
     polyMesh::updateMesh(mpm);
 
+    if (VPtr_)
+    {
+        // Grab old time volumes if the time has been incremented
+        // This will update V0, V00
+        storeOldVol(mpm.oldCellVolumes());
+
+        // Few checks
+        if (VPtr_ && (V().size() != mpm.nOldCells()))
+        {
+            FatalErrorIn("fvMesh::updateMesh(const mapPolyMesh&)")
+                << "V:" << V().size()
+                << " not equal to the number of old cells "
+                << mpm.nOldCells()
+                << exit(FatalError);
+        }
+        if (V0Ptr_ && (V0Ptr_->size() != mpm.nOldCells()))
+        {
+            FatalErrorIn("fvMesh::updateMesh(const mapPolyMesh&)")
+                << "V0:" << V0Ptr_->size()
+                << " not equal to the number of old cells "
+                << mpm.nOldCells()
+                << exit(FatalError);
+        }
+        if (V00Ptr_ && (V00Ptr_->size() != mpm.nOldCells()))
+        {
+            FatalErrorIn("fvMesh::updateMesh(const mapPolyMesh&)")
+                << "V0:" << V00Ptr_->size()
+                << " not equal to the number of old cells "
+                << mpm.nOldCells()
+                << exit(FatalError);
+        }
+    }
+
+
     // Clear the sliced fields
     clearGeomNotOldVol();
 
-    // Map all fields using current (i.e. not yet mapped) volume
+    // Map all fields
     mapFields(mpm);
 
     // Clear the current volume and other geometry factors

src/finiteVolume/fvMesh/fvMesh.H

@@ -143,6 +143,9 @@ class fvMesh
             //- Clear addressing
             void clearAddressing();
 
+            //- Preserve old volume(s)
+            void storeOldVol(const scalarField&);
+
 
        // Make geometric data