diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C index 02e5a06567..487fde704a 100644 --- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C +++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anispulation | ------------------------------------------------------------------------------- License @@ -105,6 +105,44 @@ label getBand(const labelList& owner, const labelList& neighbour) } +// Calculate band of matrix +void getBand +( + const label nCells, + const labelList& owner, + const labelList& neighbour, + label& bandwidth, + scalar& profile, // scalar to avoid overflow + scalar& sumSqrIntersect // scalar to avoid overflow +) +{ + labelList cellBandwidth(nCells, 0); + scalarField nIntersect(nCells, 0.0); + + forAll(neighbour, faceI) + { + label own = owner[faceI]; + label nei = neighbour[faceI]; + + // Note: mag not necessary for correct (upper-triangular) ordering. + label diff = nei-own; + cellBandwidth[nei] = max(cellBandwidth[nei], diff); + } + + forAll(nIntersect, cellI) + { + for (label rowI = cellI-cellBandwidth[cellI]; rowI < cellI; rowI++) + { + nIntersect[rowI]++; + } + } + + bandwidth = max(cellBandwidth); + profile = sum(cellBandwidth); + sumSqrIntersect = sum(Foam::sqr(nIntersect)); +} + + // Determine upper-triangular face order labelList getFaceOrder ( @@ -488,21 +526,12 @@ labelList regionRenumber const fvMesh& subMesh = subsetter.subMesh(); - labelList subReverseCellOrder = method.renumber + labelList subCellOrder = method.renumber ( subMesh, subMesh.cellCentres() ); - labelList subCellOrder - ( - invert - ( - subMesh.nCells(), - subReverseCellOrder - ) - ); - // Restore state UPstream::parRun() = oldParRun; @@ -556,13 +585,46 @@ int main(int argc, char *argv[]) const bool overwrite = args.optionFound("overwrite"); - label band = getBand(mesh.faceOwner(), mesh.faceNeighbour()); + label band; + scalar profile; + scalar sumSqrIntersect; + getBand + ( + mesh.nCells(), + mesh.faceOwner(), + mesh.faceNeighbour(), + band, + profile, + sumSqrIntersect + ); - Info<< "Mesh size: " << returnReduce(mesh.nCells(), sumOp