zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: renumberMesh: use renumberMethods library

Browse Source
This commit is contained in:
mattijs
2011-12-08 16:28:32 +00:00
parent f1911657d1
commit d5268bf30c
4 changed files with 433 additions and 101 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/utilities/mesh/manipulation/renumberMesh/Make/options
View File

@ -2,6 +2,7 @@ EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/renumberMethods/lnInclude \
-I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude
EXE_LIBS = \
@ -9,4 +10,5 @@ EXE_LIBS = \
-ldynamicMesh \
-lfiniteVolume \
-lgenericPatchFields \
-lrenumberMethods \
-ldecompositionMethods

335
applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
View File

@ -37,16 +37,51 @@ Description
#include "ReadFields.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "bandCompression.H"
#include "faceSet.H"
#include "SortableList.H"
#include "decompositionMethod.H"
#include "fvMeshSubset.H"
#include "renumberMethod.H"
#include "zeroGradientFvPatchFields.H"
using namespace Foam;
// Create named field from labelList for postprocessing
tmp<volScalarField> createScalarField
(
const fvMesh& mesh,
const word& name,
const labelList& elems
)
{
tmp<volScalarField> tfld
(
new volScalarField
(
IOobject
(
name,
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
mesh,
dimensionedScalar("zero", dimless, 0),
zeroGradientFvPatchScalarField::typeName
)
);
volScalarField& fld = tfld();
forAll(fld, cellI)
{
fld[cellI] = elems[cellI];
}
return tfld;
}
// Calculate band of matrix
label getBand(const labelList& owner, const labelList& neighbour)
{
@ -65,56 +100,111 @@ label getBand(const labelList& owner, const labelList& neighbour)
}
// Return new to old cell numbering
labelList regionBandCompression
// Determine upper-triangular face order
labelList getFaceOrder
(
const fvMesh& mesh,
const labelList& cellToRegion
const primitiveMesh& mesh,
const labelList& cellOrder // new to old cell
)
{
Pout<< "Determining cell order:" << endl;
labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
labelList cellOrder(cellToRegion.size());
labelList oldToNewFace(mesh.nFaces(), -1);
label nRegions = max(cellToRegion)+1;
label newFaceI = 0;
labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
labelList nbr;
labelList order;
label cellI = 0;
forAll(regionToCells, regionI)
forAll(cellOrder, newCellI)
{
Pout<< " region " << regionI << " starts at " << cellI << endl;
label oldCellI = cellOrder[newCellI];
// Per region do a reordering.
fvMeshSubset subsetter(mesh);
subsetter.setLargeCellSubset(cellToRegion, regionI);
const fvMesh& subMesh = subsetter.subMesh();
labelList subCellOrder(bandCompression(subMesh.cellCells()));
const cell& cFaces = mesh.cells()[oldCellI];
const labelList& cellMap = subsetter.cellMap();
// Neighbouring cells
nbr.setSize(cFaces.size());
forAll(subCellOrder, i)
forAll(cFaces, i)
{
cellOrder[cellI++] = cellMap[subCellOrder[i]];
label faceI = cFaces[i];
if (mesh.isInternalFace(faceI))
{
// Internal face. Get cell on other side.
label nbrCellI = reverseCellOrder[mesh.faceNeighbour()[faceI]];
if (nbrCellI == newCellI)
{
nbrCellI = reverseCellOrder[mesh.faceOwner()[faceI]];
}
if (newCellI < nbrCellI)
{
// CellI is master
nbr[i] = nbrCellI;
}
else
{
// nbrCell is master. Let it handle this face.
nbr[i] = -1;
}
}
else
{
// External face. Do later.
nbr[i] = -1;
}
}
Pout<< endl;
return cellOrder;
order.setSize(nbr.size());
sortedOrder(nbr, order);
forAll(order, i)
{
label index = order[i];
if (nbr[index] != -1)
{
oldToNewFace[cFaces[index]] = newFaceI++;
}
}
}
// Leave patch faces intact.
for (label faceI = newFaceI; faceI < mesh.nFaces(); faceI++)
{
oldToNewFace[faceI] = faceI;
}
// Check done all faces.
forAll(oldToNewFace, faceI)
{
if (oldToNewFace[faceI] == -1)
{
FatalErrorIn
(
"getFaceOrder"
"(const primitiveMesh&, const labelList&, const labelList&)"
) << "Did not determine new position"
<< " for face " << faceI
<< abort(FatalError);
}
}
return invert(mesh.nFaces(), oldToNewFace);
}
// Determine face order such that inside region faces are sorted
// upper-triangular but inbetween region faces are handled like boundary faces.
labelList regionFaceOrder
labelList getRegionFaceOrder
(
const primitiveMesh& mesh,
const labelList& cellOrder, // new to old cell
const labelList& cellToRegion // old cell to region
)
{
Pout<< "Determining face order:" << endl;
//Pout<< "Determining face order:" << endl;
labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
@ -131,8 +221,8 @@ labelList regionFaceOrder
if (cellToRegion[oldCellI] != prevRegion)
{
prevRegion = cellToRegion[oldCellI];
Pout<< " region " << prevRegion << " internal faces start at "
<< newFaceI << endl;
//Pout<< " region " << prevRegion << " internal faces start at "
// << newFaceI << endl;
}
const cell& cFaces = mesh.cells()[oldCellI];
@ -219,9 +309,9 @@ labelList regionFaceOrder
if (prevKey != key)
{
Pout<< " faces inbetween region " << key/nRegions
<< " and " << key%nRegions
<< " start at " << newFaceI << endl;
//Pout<< " faces inbetween region " << key/nRegions
// << " and " << key%nRegions
// << " start at " << newFaceI << endl;
prevKey = key;
}
@ -243,14 +333,14 @@ labelList regionFaceOrder
{
FatalErrorIn
(
"polyDualMesh::getFaceOrder"
"(const labelList&, const labelList&, const label) const"
"getRegionFaceOrder"
"(const primitveMesh&, const labelList&, const labelList&)"
) << "Did not determine new position"
<< " for face " << faceI
<< abort(FatalError);
}
}
Pout<< endl;
//Pout<< endl;
return invert(mesh.nFaces(), oldToNewFace);
}
@ -365,10 +455,11 @@ autoPtr<mapPolyMesh> reorderMesh
int main(int argc, char *argv[])
{
argList::addBoolOption
argList::addOption
(
"blockOrder",
"order cells into regions (using decomposition)"
"blockSize",
"block size",
"order cells into blocks (using decomposition) before ordering"
);
argList::addBoolOption
(
@ -400,12 +491,23 @@ int main(int argc, char *argv[])
# include "createNamedMesh.H"
const word oldInstance = mesh.pointsInstance();
const bool blockOrder = args.optionFound("blockOrder");
if (blockOrder)
label blockSize = 0;
args.optionReadIfPresent("blockSize", blockSize, 0);
if (blockSize > 0)
{
Info<< "Ordering cells into regions (using decomposition);"
Info<< "Ordering cells into regions of size " << blockSize
<< " (using decomposition);"
<< " ordering faces into region-internal and region-external." << nl
<< endl;
if (blockSize < 0 || blockSize >= mesh.nCells())
{
FatalErrorIn(args.executable())
<< "Block size " << blockSize << " should be positive integer"
<< " and less than the number of cells in the mesh."
<< exit(FatalError);
}
}
const bool orderPoints = args.optionFound("orderPoints");
@ -432,6 +534,24 @@ int main(int argc, char *argv[])
<< returnReduce(band, maxOp<label>()) << nl << endl;
// Construct renumberMethod
IOdictionary renumberDict
(
IOobject
(
"renumberMeshDict",
runTime.system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
autoPtr<renumberMethod> renumberPtr = renumberMethod::New(renumberDict);
Info<< "Selecting renumberMethod " << renumberPtr().type() << endl;
// Read parallel reconstruct maps
labelIOList cellProcAddressing
(
@ -522,13 +642,19 @@ int main(int argc, char *argv[])
ReadFields(mesh, objects, stFlds);
autoPtr<mapPolyMesh> map;
// From renumbering:
// - from new cell/face back to original cell/face
labelList cellOrder;
labelList faceOrder;
if (blockOrder)
if (blockSize > 0)
{
// Renumbering in two phases. Should be done in one so mapping of
// fields is done correctly!
label nBlocks = mesh.nCells() / blockSize;
Info<< "nBlocks = " << nBlocks << endl;
// Read decomposePar dictionary
IOdictionary decomposeDict
(
@ -541,6 +667,8 @@ int main(int argc, char *argv[])
IOobject::NO_WRITE
)
);
decomposeDict.set("numberOfSubdomains", nBlocks);
autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
(
decomposeDict
@ -556,80 +684,99 @@ int main(int argc, char *argv[])
);
// For debugging: write out region
{
volScalarField cellDist
(
IOobject
createScalarField
(
mesh,
"cellDist",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh,
dimensionedScalar("cellDist", dimless, 0),
zeroGradientFvPatchScalarField::typeName
);
forAll(cellToRegion, cellI)
{
cellDist[cellI] = cellToRegion[cellI];
}
cellDist.write();
cellToRegion
)().write();
Info<< nl << "Written decomposition as volScalarField to "
<< cellDist.name() << " for use in postprocessing."
<< "cellDist for use in postprocessing."
<< nl << endl;
// Find point per region
pointField regionPoints(nBlocks, vector::zero);
forAll(cellToRegion, cellI)
{
regionPoints[cellToRegion[cellI]] = mesh.cellCentres()[cellI];
}
// Use block based renumbering.
//labelList cellOrder(bandCompression(mesh.cellCells()));
labelList cellOrder(regionBandCompression(mesh, cellToRegion));
// Detemines old to new cell ordering
labelList reverseCellOrder = renumberPtr().renumber
(
mesh,
cellToRegion,
regionPoints
);
cellOrder = invert(mesh.nCells(), reverseCellOrder);
// Determine new to old face order with new cell numbering
labelList faceOrder
(
regionFaceOrder
faceOrder = getRegionFaceOrder
(
mesh,
cellOrder,
cellToRegion
)
);
if (!overwrite)
{
runTime++;
}
// Change the mesh.
map = reorderMesh(mesh, cellOrder, faceOrder);
}
else
{
// Use built-in renumbering.
// Detemines old to new cell ordering
labelList reverseCellOrder = renumberPtr().renumber
(
mesh,
mesh.cellCentres()
);
cellOrder = invert(mesh.nCells(), reverseCellOrder);
// Determine new to old face order with new cell numbering
faceOrder = getFaceOrder
(
mesh,
cellOrder // new to old cell
);
}
polyTopoChange meshMod(mesh);
if (!overwrite)
{
runTime++;
}
// Change the mesh.
map = meshMod.changeMesh
autoPtr<mapPolyMesh> map = reorderMesh(mesh, cellOrder, faceOrder);
if (orderPoints)
{
polyTopoChange meshMod(mesh);
autoPtr<mapPolyMesh> pointOrderMap = meshMod.changeMesh
(
mesh,
false, //inflate
true, // parallel sync
true, // cell ordering
orderPoints // point ordering
true, //syncParallel
false, //orderCells
orderPoints //orderPoints
);
// Combine point reordering into map.
const_cast<labelList&>(map().pointMap()) = UIndirectList<label>
(
map().pointMap(),
pointOrderMap().pointMap()
)();
inplaceRenumber
(
pointOrderMap().reversePointMap(),
const_cast<labelList&>(map().reversePointMap())
);
}
// Update fields
mesh.updateMesh(map);
@ -767,6 +914,24 @@ int main(int argc, char *argv[])
if (writeMaps)
{
// For debugging: write out region
createScalarField
(
mesh,
"origCellID",
map().cellMap()
)().write();
createScalarField
(
mesh,
"cellID",
identity(mesh.nCells())
)().write();
Info<< nl << "Written current cellID and origCellID as volScalarField"
<< " for use in postprocessing."
<< nl << endl;
labelIOList
(
IOobject

51
applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict Normal file
View File

@ -0,0 +1,51 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: dev |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
note "mesh renumbering dictionary";
object renumberMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
method CuthillMcKee;
//method manual;
//method random;
//method spring;
CuthillMcKeeCoeffs
{
// Reverse CuthillMcKee (RCM) or plain
reverse true;
}
manualCoeffs
{
// In system directory: old to new labelIOList
dataFile "cellMap";
}
springCoeffs
{
// Maximum jump of cell indices. Is fraction of number of cells
maxCo 0.1;
// Limit the amount of movement; the fraction maxCo gets decreased
// with every iteration
freezeFraction 0.9;
// Maximum number of iterations
maxIter 1000;
}
// ************************************************************************* //

114
src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C Normal file
View File

@ -0,0 +1,114 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "CuthillMcKeeRenumber.H"
#include "addToRunTimeSelectionTable.H"
#include "bandCompression.H"
#include "decompositionMethod.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(CuthillMcKeeRenumber, 0);
addToRunTimeSelectionTable
(
renumberMethod,
CuthillMcKeeRenumber,
dictionary
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::CuthillMcKeeRenumber::CuthillMcKeeRenumber(const dictionary& renumberDict)
:
renumberMethod(renumberDict),
reverse_
(
renumberDict.found(typeName + "Coeffs")
? Switch(renumberDict.subDict(typeName + "Coeffs").lookup("reverse"))
: Switch(true)
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::labelList Foam::CuthillMcKeeRenumber::renumber
(
const polyMesh& mesh,
const pointField& points
)
{
CompactListList<label> cellCells;
decompositionMethod::calcCellCells
(
mesh,
identity(mesh.nCells()),
mesh.nCells(),
false, // local only
cellCells
);
labelList orderedToOld = bandCompression(cellCells());
if (reverse_)
{
Pout<< "** REverseing.." << endl;
for (label i = 0; i < orderedToOld.size()/2; i++)
{
Swap(orderedToOld[i], orderedToOld[orderedToOld.size()-i-1]);
}
}
return invert(orderedToOld.size(), orderedToOld);
}
Foam::labelList Foam::CuthillMcKeeRenumber::renumber
(
const labelListList& cellCells,
const pointField& points
)
{
labelList orderedToOld = bandCompression(cellCells);
if (reverse_)
{
Pout<< "** REverseing.." << endl;
for (label i = 0; i < orderedToOld.size()/2; i++)
{
Swap(orderedToOld[i], orderedToOld[orderedToOld.size()-i-1]);
}
}
return invert(orderedToOld.size(), orderedToOld);
}
// ************************************************************************* //