diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 64011258c6ae2f67ac49b81ed44f5caaa56ac67f..3702cfe8ce81df764b4124e1112847ebe3428566 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -145,6 +145,7 @@ Usage
#include "hexRef8Data.H"
#include "regionProperties.H"
#include "polyMeshTools.H"
+#include "subsetAdjacency.H"
using namespace Foam;
@@ -194,21 +195,56 @@ tmp<volScalarField> createScalarField
}
-// Calculate band of matrix
-label getBand(const labelList& owner, const labelList& neighbour)
+// Calculate band of mesh
+// label getBand(const labelUList& owner, const labelUList& neighbour)
+// {
+// label bandwidth = 0;
+//
+// forAll(neighbour, facei)
+// {
+// const label width = neighbour[facei] - owner[facei];
+//
+// if (bandwidth < width)
+// {
+// bandwidth = width;
+// }
+// }
+// return bandwidth;
+// }
+
+
+// Calculate band and profile of matrix. Profile is scalar to avoid overflow
+Tuple2<label, scalar> getBand
+(
+ const CompactListList<label>& mat
+)
{
- label band = 0;
+ Tuple2<label, scalar> metrics(0, 0);
- forAll(neighbour, facei)
+ auto& bandwidth = metrics.first();
+ auto& profile = metrics.second();
+
+ forAll(mat, celli)
{
- label diff = neighbour[facei] - owner[facei];
+ const auto& neighbours = mat[celli];
+
+ const label nNbr = neighbours.size();
- if (diff > band)
+ if (nNbr)
{
- band = diff;
+ // Max distance
+ const label width = (neighbours[nNbr-1] - celli);
+
+ if (bandwidth < width)
+ {
+ bandwidth = width;
+ }
+
+ profile += scalar(width);
}
}
- return band;
+
+ return metrics;
}
@@ -217,27 +253,35 @@ void getBand
(
const bool calculateIntersect,
const label nCells,
- const labelList& owner,
- const labelList& neighbour,
+ const labelUList& owner,
+ const labelUList& neighbour,
label& bandwidth,
scalar& profile, // scalar to avoid overflow
scalar& sumSqrIntersect // scalar to avoid overflow
)
{
labelList cellBandwidth(nCells, Foam::zero{});
- scalarField nIntersect(nCells, Foam::zero{});
+
+ bandwidth = 0;
forAll(neighbour, facei)
{
- label own = owner[facei];
- label nei = neighbour[facei];
+ const label own = owner[facei];
+ const label nei = neighbour[facei];
// Note: mag not necessary for correct (upper-triangular) ordering.
- label diff = nei-own;
- cellBandwidth[nei] = max(cellBandwidth[nei], diff);
- }
+ const label width = nei - own;
+
+ if (cellBandwidth[nei] < width)
+ {
+ cellBandwidth[nei] = width;
- bandwidth = max(cellBandwidth);
+ if (bandwidth < width)
+ {
+ bandwidth = width;
+ }
+ }
+ }
// Do not use field algebra because of conversion label to scalar
profile = 0;
@@ -246,14 +290,16 @@ void getBand
profile += scalar(width);
}
- sumSqrIntersect = 0.0;
+ sumSqrIntersect = 0;
if (calculateIntersect)
{
+ scalarField nIntersect(nCells, Foam::zero{});
+
forAll(nIntersect, celli)
{
for (label colI = celli-cellBandwidth[celli]; colI <= celli; colI++)
{
- nIntersect[colI] += 1.0;
+ nIntersect[colI] += scalar(1);
}
}
@@ -675,9 +721,8 @@ CompactListList<label> regionRenumber
forAll(regionCellOrder, regioni)
{
- // Info<< " region " << regioni
- // << " starts at " << regionCellOrder.localStart(regioni)
- // << nl;
+ // Info<< " region " << regioni << " starts at "
+ // << regionCellOrder.localStart(regioni) << nl;
// No parallel communication
const bool oldParRun = UPstream::parRun(false);
@@ -703,27 +748,42 @@ CompactListList<label> regionRenumber
{
timer.resetTimeIncrement();
+ // Create adjacency matrix of the full mesh and subset subsequently.
+ // This is more efficient than creating adjacency matrices of
+ // sub-meshes.
+
+ // No parallel communication
+ const bool oldParRun = UPstream::parRun(false);
+
+ // The local connectivity of the full (non-subsetted) mesh
+ CompactListList<label> meshCellCells;
+ globalMeshData::calcCellCells(mesh, meshCellCells);
+ UPstream::parRun(oldParRun); // Restore parallel state
+
+ timings[TimingType::CELL_CELLS] += timer.timeIncrement();
+
+ // For the respective subMesh selections
+ bitSet subsetCells(mesh.nCells());
+
forAll(regionCellOrder, regioni)
{
- // Info<< " region " << regioni
- // << " starts at " << regionCellOrder.localStart(regioni)
- // << nl;
+ // Info<< " region " << regioni << " starts at "
+ // << regionCellOrder.localStart(regioni) << nl;
- // No parallel communication
- const bool oldParRun = UPstream::parRun(false);
+ subsetCells = false;
+ subsetCells.set(regionCellOrder[regioni]);
// Connectivity of local sub-mesh
- CompactListList<label> cellCells;
- labelList cellMap = globalMeshData::calcCellCells
- (
- mesh,
- regionCellOrder[regioni],
- cellCells
- );
+ labelList cellMap;
+ CompactListList<label> subCellCells =
+ subsetAdjacency(subsetCells, meshCellCells, cellMap);
timings[TimingType::CELL_CELLS] += timer.timeIncrement();
- labelList subCellOrder = method.renumber(cellCells);
+ // No parallel communication
+ const bool oldParRun = UPstream::parRun(false);
+
+ labelList subCellOrder = method.renumber(subCellCells);
UPstream::parRun(oldParRun); // Restore parallel state
@@ -835,10 +895,10 @@ int main(int argc, char *argv[])
const bool dryrun = args.dryRun();
const bool readDict = args.found("dict");
- const bool doFrontWidth = args.found("frontWidth");
+ const bool doDecompose = args.found("decompose");
const bool overwrite = args.found("overwrite");
const bool doFields = !args.found("no-fields");
- const bool doDecompose = args.found("decompose");
+ const bool doFrontWidth = args.found("frontWidth") && !doDecompose;
word renumberMethodName;
args.readIfPresent("renumber-method", renumberMethodName);
@@ -846,8 +906,7 @@ int main(int argc, char *argv[])
if (doDecompose && UPstream::parRun())
{
FatalErrorIn(args.executable())
- << "Cannot use -decompose option in parallel"
- << " ... giving up" << nl
+ << "Cannot use -decompose option in parallel ... giving up" << nl
<< exit(FatalError);
}
@@ -908,21 +967,21 @@ int main(int argc, char *argv[])
reduce(band, maxOp<label>());
reduce(profile, sumOp<scalar>());
- reduce(sumSqrIntersect, sumOp<scalar>());
-
- scalar rmsFrontwidth = Foam::sqrt
- (
- sumSqrIntersect/mesh.globalData().nTotalCells()
- );
Info<< "Mesh " << mesh.name()
<< " size: " << mesh.globalData().nTotalCells() << nl
- << "Before renumbering :" << nl
+ << "Before renumbering" << nl
<< " band : " << band << nl
<< " profile : " << profile << nl;
if (doFrontWidth)
{
+ reduce(sumSqrIntersect, sumOp<scalar>());
+ scalar rmsFrontwidth = Foam::sqrt
+ (
+ sumSqrIntersect/mesh.globalData().nTotalCells()
+ );
+
Info<< " rms frontwidth : " << rmsFrontwidth << nl;
}
@@ -1091,10 +1150,7 @@ int main(int argc, char *argv[])
);
- // List of objects read from time directory
- // List of stored objects to clear from mesh
-
- IOobjectList objects;
+ // List of stored objects to clear from mesh (after reading)
DynamicList<regIOobject*> storedObjects;
if (!dryrun && doFields)
@@ -1103,55 +1159,60 @@ int main(int argc, char *argv[])
timer.resetTimeIncrement();
- objects = IOobjectList(mesh, runTime.timeName());
+ IOobjectList objects(mesh, runTime.timeName());
storedObjects.reserve(objects.size());
const predicates::always nameMatcher;
// Read GeometricFields
- #undef ReadFields
- #define ReadFields(FieldType) \
+ #undef doLocalCode
+ #define doLocalCode(FieldType) \
readFields<FieldType>(mesh, objects, nameMatcher, storedObjects);
// Read volume fields
- ReadFields(volScalarField);
- ReadFields(volVectorField);
- ReadFields(volSphericalTensorField);
- ReadFields(volSymmTensorField);
- ReadFields(volTensorField);
+ doLocalCode(volScalarField);
+ doLocalCode(volVectorField);
+ doLocalCode(volSphericalTensorField);
+ doLocalCode(volSymmTensorField);
+ doLocalCode(volTensorField);
// Read internal fields
- ReadFields(volScalarField::Internal);
- ReadFields(volVectorField::Internal);
- ReadFields(volSphericalTensorField::Internal);
- ReadFields(volSymmTensorField::Internal);
- ReadFields(volTensorField::Internal);
+ doLocalCode(volScalarField::Internal);
+ doLocalCode(volVectorField::Internal);
+ doLocalCode(volSphericalTensorField::Internal);
+ doLocalCode(volSymmTensorField::Internal);
+ doLocalCode(volTensorField::Internal);
// Read surface fields
- ReadFields(surfaceScalarField);
- ReadFields(surfaceVectorField);
- ReadFields(surfaceSphericalTensorField);
- ReadFields(surfaceSymmTensorField);
- ReadFields(surfaceTensorField);
+ doLocalCode(surfaceScalarField);
+ doLocalCode(surfaceVectorField);
+ doLocalCode(surfaceSphericalTensorField);
+ doLocalCode(surfaceSymmTensorField);
+ doLocalCode(surfaceTensorField);
// Read point fields
const pointMesh& pMesh = pointMesh::New(mesh);
- #undef ReadPointFields
- #define ReadPointFields(FieldType) \
+ #undef doLocalCode
+ #define doLocalCode(FieldType) \
readFields<FieldType>(pMesh, objects, nameMatcher, storedObjects);
- ReadPointFields(pointScalarField);
- ReadPointFields(pointVectorField);
- ReadPointFields(pointSphericalTensorField);
- ReadPointFields(pointSymmTensorField);
- ReadPointFields(pointTensorField);
+ doLocalCode(pointScalarField);
+ doLocalCode(pointVectorField);
+ doLocalCode(pointSphericalTensorField);
+ doLocalCode(pointSymmTensorField);
+ doLocalCode(pointTensorField);
- #undef ReadFields
- #undef ReadPointFields
+ #undef doLocalCode
timings[TimingType::READ_FIELDS] += timer.timeIncrement();
+
+ // Write loaded fields when mesh.write() is called
+ for (auto* fldptr : storedObjects)
+ {
+ fldptr->writeOpt(IOobject::AUTO_WRITE);
+ }
}
@@ -1224,7 +1285,13 @@ int main(int argc, char *argv[])
CompactListList<label> regionCellOrder =
- regionRenumber(renumberPtr(), mesh, cellToRegion);
+ regionRenumber
+ (
+ renumberPtr(),
+ mesh,
+ cellToRegion,
+ decomposePtr().nDomains()
+ );
cellOrder = regionCellOrder.values();
@@ -1586,19 +1653,25 @@ int main(int argc, char *argv[])
);
reduce(band, maxOp<label>());
reduce(profile, sumOp<scalar>());
- reduce(sumSqrIntersect, sumOp<scalar>());
- scalar rmsFrontwidth = Foam::sqrt
- (
- sumSqrIntersect/mesh.globalData().nTotalCells()
- );
+ Info<< "After renumbering";
+ if (doDecompose)
+ {
+ Info<< " [values are misleading with -decompose option]";
+ }
- Info<< "After renumbering :" << nl
+ Info<< nl
<< " band : " << band << nl
<< " profile : " << profile << nl;
if (doFrontWidth)
{
+ reduce(sumSqrIntersect, sumOp<scalar>());
+ scalar rmsFrontwidth = Foam::sqrt
+ (
+ sumSqrIntersect/mesh.globalData().nTotalCells()
+ );
+
Info<< " rms frontwidth : " << rmsFrontwidth << nl;
}
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/subsetAdjacency.H b/applications/utilities/mesh/manipulation/renumberMesh/subsetAdjacency.H
new file mode 100644
index 0000000000000000000000000000000000000000..0c098616cf2396faa4daa85f53bc0742157e0464
--- /dev/null
+++ b/applications/utilities/mesh/manipulation/renumberMesh/subsetAdjacency.H
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
+
+Description
+ Subsetting of an adjacency matrix (as CompactListList).
+ Can be relocated elsewhere.
+
+\*---------------------------------------------------------------------------*/
+
+#include "CompactListList.H"
+#include "bitSet.H"
+#include "ListOps.H"
+#include "Map.H"
+
+namespace Foam
+{
+
+// Perform a subset of the adjacency matrix
+CompactListList<label> subsetAdjacency
+(
+ const bitSet& select, // could also be labelHashSet
+ const CompactListList<label>& input,
+ labelList& subMap
+)
+{
+ // Corresponds to cellMap etc (the original selection)
+ subMap = select.sortedToc();
+
+ // Ensure that the subMap corresponds to a valid subset
+ {
+ label validSize = 0;
+
+ const label nTotal = input.size();
+
+ forAllReverse(subMap, i)
+ {
+ if (subMap[i] < nTotal)
+ {
+ validSize = i + 1;
+ break;
+ }
+ }
+
+ subMap.resize(validSize);
+ }
+
+
+ // Assumed to be sparse - use Map for reverse lookup
+ const Map<label> reverseMap(invertToMap(subMap));
+
+
+ // Pass 1: determine the selected sub-sizes
+ labelList sizes(subMap.size(), Foam::zero{});
+
+ forAll(subMap, idx)
+ {
+ for (const label nbr : input[subMap[idx]])
+ {
+ if
+ (
+ select.test(nbr)
+ && reverseMap.contains(nbr) // extra consistency (paranoid)
+ )
+ {
+ ++sizes[idx];
+ }
+ }
+ }
+
+
+ CompactListList<label> output(sizes);
+
+ // Reuse sizes as output offset into output.values()
+ sizes = labelList::subList(output.offsets(), output.size());
+ labelList& values = output.values();
+
+
+ // Pass 2: extract sub-adjacent matrix
+
+ label newNbr = -1;
+
+ forAll(subMap, idx)
+ {
+ for (const label nbr : input[subMap[idx]])
+ {
+ if
+ (
+ select.test(nbr)
+ && (newNbr = reverseMap.lookup(nbr, -1)) >= 0
+ )
+ {
+ values[sizes[idx]++] = newNbr;
+ }
+ }
+ }
+
+ return output;
+}
+
+
+// Perform a subset of the adjacency matrix
+CompactListList<label> subsetAdjacency
+(
+ const labelRange& slice,
+ const CompactListList<label>& input,
+ labelList& subMap
+)
+{
+ // Ensure that the selection corresponds to a valid subset
+ const labelRange select = slice.subset0(input.size());
+
+ // Corresponds to cellMap etc (the original selection)
+ subMap = Foam::identity(select);
+
+
+ // Pass 1: determine the selected sub-sizes
+ labelList sizes(subMap.size(), Foam::zero{});
+
+ forAll(subMap, idx)
+ {
+ for (const label nbr : input[subMap[idx]])
+ {
+ if (select.contains(nbr))
+ {
+ ++sizes[idx];
+ }
+ }
+ }
+
+
+ CompactListList<label> output(sizes);
+
+ // Reuse sizes as output offset into output.values()
+ sizes = labelList::subList(output.offsets(), output.size());
+ labelList& values = output.values();
+
+
+ // Pass 2: extract sub-adjacent matrix
+
+ const label localOffset = select.start();
+
+ forAll(subMap, idx)
+ {
+ for (const label nbr : input[subMap[idx]])
+ {
+ if (select.contains(nbr))
+ {
+ values[sizes[idx]++] = nbr - localOffset;
+ }
+ }
+ }
+
+ return output;
+}
+
+} // End namespace Foam
+
+
+// ************************************************************************* //