diff --git a/applications/utilities/mesh/manipulation/renumberMesh/Make/options b/applications/utilities/mesh/manipulation/renumberMesh/Make/options
index b5e293282a48fcc11891c2222ad3a6b75ee23d2f..2981cf30bea351bd387501cbaa5c65443ae6b769 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/renumberMesh/Make/options
@@ -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
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 8d9b0c91ebab1b0c35cb032be821c1c711c7292a..4bb6512a35fbaa331ec7753811f43c2330981143 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -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 cellOrder(cellToRegion.size());
+ labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
- label nRegions = max(cellToRegion)+1;
+ labelList oldToNewFace(mesh.nFaces(), -1);
- labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
+ label newFaceI = 0;
- label cellI = 0;
+ labelList nbr;
+ labelList order;
- 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)
+ {
+ 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;
+ }
+ }
+
+ 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)
{
- cellOrder[cellI++] = cellMap[subCellOrder[i]];
+ FatalErrorIn
+ (
+ "getFaceOrder"
+ "(const primitiveMesh&, const labelList&, const labelList&)"
+ ) << "Did not determine new position"
+ << " for face " << faceI
+ << abort(FatalError);
}
}
- Pout<< endl;
- return cellOrder;
+ 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
- (
- "cellDist",
- runTime.timeName(),
- mesh,
- IOobject::NO_READ,
- IOobject::NO_WRITE,
- false
- ),
- mesh,
- dimensionedScalar("cellDist", dimless, 0),
- zeroGradientFvPatchScalarField::typeName
- );
+ createScalarField
+ (
+ mesh,
+ "cellDist",
+ cellToRegion
+ )().write();
- forAll(cellToRegion, cellI)
- {
- cellDist[cellI] = cellToRegion[cellI];
- }
+ Info<< nl << "Written decomposition as volScalarField to "
+ << "cellDist for use in postprocessing."
+ << nl << endl;
- cellDist.write();
- Info<< nl << "Written decomposition as volScalarField to "
- << cellDist.name() << " 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
+ faceOrder = getRegionFaceOrder
(
- regionFaceOrder
- (
- mesh,
- cellOrder,
- cellToRegion
- )
+ mesh,
+ cellOrder,
+ cellToRegion
+ );
+ }
+ else
+ {
+ // Detemines old to new cell ordering
+ labelList reverseCellOrder = renumberPtr().renumber
+ (
+ mesh,
+ mesh.cellCentres()
);
- if (!overwrite)
- {
- runTime++;
- }
+ cellOrder = invert(mesh.nCells(), reverseCellOrder);
- // Change the mesh.
- map = reorderMesh(mesh, cellOrder, faceOrder);
+ // Determine new to old face order with new cell numbering
+ faceOrder = getFaceOrder
+ (
+ mesh,
+ cellOrder // new to old cell
+ );
}
- else
+
+
+ if (!overwrite)
{
- // Use built-in renumbering.
+ runTime++;
+ }
- polyTopoChange meshMod(mesh);
- if (!overwrite)
- {
- runTime++;
- }
+ // Change the mesh.
+ autoPtr<mapPolyMesh> map = reorderMesh(mesh, cellOrder, faceOrder);
- // Change the mesh.
- map = meshMod.changeMesh
+
+ if (orderPoints)
+ {
+ polyTopoChange meshMod(mesh);
+ autoPtr<mapPolyMesh> pointOrderMap = meshMod.changeMesh
(
mesh,
- false, // inflate
- true, // parallel sync
- true, // cell ordering
- orderPoints // point ordering
+ false, //inflate
+ 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
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
new file mode 100644
index 0000000000000000000000000000000000000000..7cd0fb8c5168a8d5eb3c7bcf3b958b2e24173c56
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
@@ -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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C b/src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
new file mode 100644
index 0000000000000000000000000000000000000000..29e9e19683217128fb3888e319fe7cb3b1892f65
--- /dev/null
+++ b/src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
@@ -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);
+}
+
+
+// ************************************************************************* //