diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C index 6a8ae4a947..666293c558 100644 --- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C +++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C @@ -555,7 +555,8 @@ autoPtr reorderMesh labelListList(0), // cellZoneMap, pointField(0), // preMotionPoints, patchStarts, // oldPatchStarts, - oldPatchNMeshPoints // oldPatchNMeshPoints + oldPatchNMeshPoints, // oldPatchNMeshPoints + autoPtr() // oldCellVolumes ) ); } diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/cellMapper/cellMapper.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/cellMapper/cellMapper.C index eeb077ed01..84dfacf47c 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/cellMapper/cellMapper.C +++ b/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& 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,32 +165,74 @@ 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 - //w[cellI] = scalarList(mo.size(), 1.0/mo.size()); + if (mpm_.hasOldCellVolumes()) + { + // Volume weighted - //- volume based - w[cellI].setSize(mo.size()); + const scalarField& V = mpm_.oldCellVolumes(); - if (mo.size()) + if (V.size() != sizeBeforeMapping()) { - scalar sumV = 0; - forAll(mo, ci) + 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(); + + w[cellI].setSize(mo.size()); + + if (mo.size()) { - w[cellI][ci] = V[mo[ci]]; - sumV += V[mo[ci]]; - } - forAll(mo, ci) - { - w[cellI][ci] /= sumV; + scalar sumV = 0; + forAll(mo, ci) + { + 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()); + } + } // Do mapped faces. Note that can already be set from cellsFromCells diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.C index faf2922955..85e69b02dc 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.C +++ b/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& 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& 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) { diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.H index 2bad5824a1..e21ff4c736 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapPolyMesh.H +++ b/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 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& oldCellVolumesPtr ); //- Construct from components and optionally reuse storage @@ -344,6 +348,7 @@ public: pointField& preMotionPoints, labelList& oldPatchStarts, labelList& oldPatchNMeshPoints, + autoPtr& 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_(); + } + }; diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C index 23e39f203d..bb53ef1aa0 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C @@ -3145,6 +3145,7 @@ Foam::autoPtr Foam::polyTopoChange::changeMesh const label nOldPoints(mesh.nPoints()); const label nOldFaces(mesh.nFaces()); const label nOldCells(mesh.nCells()); + autoPtr oldCellVolumes(new scalarField(mesh.cellVolumes())); // Change the mesh @@ -3323,6 +3324,9 @@ Foam::autoPtr Foam::polyTopoChange::changeMesh newPoints, // if empty signals no inflation. oldPatchStarts, oldPatchNMeshPoints, + + oldCellVolumes, + true // steal storage. ) ); @@ -3408,6 +3412,7 @@ Foam::autoPtr Foam::polyTopoChange::makeMesh const label nOldPoints(mesh.nPoints()); const label nOldFaces(mesh.nFaces()); const label nOldCells(mesh.nCells()); + autoPtr oldCellVolumes(new scalarField(mesh.cellVolumes())); // Create the mesh @@ -3617,6 +3622,7 @@ Foam::autoPtr Foam::polyTopoChange::makeMesh newPoints, // if empty signals no inflation. oldPatchStarts, oldPatchNMeshPoints, + oldCellVolumes, true // steal storage. ) ); diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C index 6713435149..26c3b33423 100644 --- a/src/finiteVolume/fvMesh/fvMesh.C +++ b/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 + ( + 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,8 +613,31 @@ 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::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_) - { - *V00Ptr_ = *V0Ptr_; - } - - if (V0Ptr_) - { - *V0Ptr_ = V(); - } - else - { - V0Ptr_ = new DimensionedField - ( - IOobject - ( - "V0", - time().timeName(), - *this, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - V() - ); - } - - curTimeIndex_ = time().timeIndex(); + storeOldVol(V()); } - 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 diff --git a/src/finiteVolume/fvMesh/fvMesh.H b/src/finiteVolume/fvMesh/fvMesh.H index cf3da4b115..c6c234c495 100644 --- a/src/finiteVolume/fvMesh/fvMesh.H +++ b/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