Files
openfoam/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C

1420 lines
41 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 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/>.
InClass
decompositionMethod
\*---------------------------------------------------------------------------*/
#include "decompositionMethod.H"
#include "globalIndex.H"
#include "syncTools.H"
#include "Tuple2.H"
#include "faceSet.H"
#include "regionSplit.H"
#include "localPointRegion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(decompositionMethod, 0);
defineRunTimeSelectionTable(decompositionMethod, dictionary);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::autoPtr<Foam::decompositionMethod> Foam::decompositionMethod::New
(
const dictionary& decompositionDict
)
{
word methodType(decompositionDict.lookup("method"));
if (methodType == "scotch" && Pstream::parRun())
{
methodType = "ptscotch";
}
Info<< "Selecting decompositionMethod " << methodType << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(methodType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn
(
"decompositionMethod::New"
"(const dictionary& decompositionDict)"
) << "Unknown decompositionMethod "
<< methodType << nl << nl
<< "Valid decompositionMethods are : " << endl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<decompositionMethod>(cstrIter()(decompositionDict));
}
Foam::labelList Foam::decompositionMethod::decompose
(
const polyMesh& mesh,
const pointField& points
)
{
scalarField weights(points.size(), 1.0);
return decompose(mesh, points, weights);
}
Foam::labelList Foam::decompositionMethod::decompose
(
const polyMesh& mesh,
const labelList& fineToCoarse,
const pointField& coarsePoints,
const scalarField& coarseWeights
)
{
CompactListList<label> coarseCellCells;
calcCellCells
(
mesh,
fineToCoarse,
coarsePoints.size(),
true, // use global cell labels
coarseCellCells
);
// Decompose based on agglomerated points
labelList coarseDistribution
(
decompose
(
coarseCellCells(),
coarsePoints,
coarseWeights
)
);
// Rework back into decomposition for original mesh_
labelList fineDistribution(fineToCoarse.size());
forAll(fineDistribution, i)
{
fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
}
return fineDistribution;
}
Foam::labelList Foam::decompositionMethod::decompose
(
const polyMesh& mesh,
const labelList& fineToCoarse,
const pointField& coarsePoints
)
{
scalarField cWeights(coarsePoints.size(), 1.0);
return decompose
(
mesh,
fineToCoarse,
coarsePoints,
cWeights
);
}
Foam::labelList Foam::decompositionMethod::decompose
(
const labelListList& globalCellCells,
const pointField& cc
)
{
scalarField cWeights(cc.size(), 1.0);
return decompose(globalCellCells, cc, cWeights);
}
void Foam::decompositionMethod::calcCellCells
(
const polyMesh& mesh,
const labelList& agglom,
const label nCoarse,
const bool parallel,
CompactListList<label>& cellCells
)
{
const labelList& faceOwner = mesh.faceOwner();
const labelList& faceNeighbour = mesh.faceNeighbour();
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Create global cell numbers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
globalIndex globalAgglom
(
nCoarse,
Pstream::msgType(),
Pstream::worldComm,
parallel
);
// Get agglomerate owner on other side of coupled faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
{
label faceI = pp.start();
label bFaceI = pp.start() - mesh.nInternalFaces();
forAll(pp, i)
{
globalNeighbour[bFaceI] = globalAgglom.toGlobal
(
agglom[faceOwner[faceI]]
);
bFaceI++;
faceI++;
}
}
}
// Get the cell on the other side of coupled patches
syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
// Count number of faces (internal + coupled)
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Number of faces per coarse cell
labelList nFacesPerCell(nCoarse, 0);
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
label own = agglom[faceOwner[faceI]];
label nei = agglom[faceNeighbour[faceI]];
nFacesPerCell[own]++;
nFacesPerCell[nei]++;
}
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
{
label faceI = pp.start();
label bFaceI = pp.start()-mesh.nInternalFaces();
forAll(pp, i)
{
label own = agglom[faceOwner[faceI]];
label globalNei = globalNeighbour[bFaceI];
if
(
!globalAgglom.isLocal(globalNei)
|| globalAgglom.toLocal(globalNei) != own
)
{
nFacesPerCell[own]++;
}
faceI++;
bFaceI++;
}
}
}
// Fill in offset and data
// ~~~~~~~~~~~~~~~~~~~~~~~
cellCells.setSize(nFacesPerCell);
nFacesPerCell = 0;
labelList& m = cellCells.m();
const labelList& offsets = cellCells.offsets();
// For internal faces is just offsetted owner and neighbour
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
label own = agglom[faceOwner[faceI]];
label nei = agglom[faceNeighbour[faceI]];
m[offsets[own] + nFacesPerCell[own]++] = globalAgglom.toGlobal(nei);
m[offsets[nei] + nFacesPerCell[nei]++] = globalAgglom.toGlobal(own);
}
// For boundary faces is offsetted coupled neighbour
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
{
label faceI = pp.start();
label bFaceI = pp.start()-mesh.nInternalFaces();
forAll(pp, i)
{
label own = agglom[faceOwner[faceI]];
label globalNei = globalNeighbour[bFaceI];
if
(
!globalAgglom.isLocal(globalNei)
|| globalAgglom.toLocal(globalNei) != own
)
{
m[offsets[own] + nFacesPerCell[own]++] = globalNei;
}
faceI++;
bFaceI++;
}
}
}
// Check for duplicates connections between cells
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Done as postprocessing step since we now have cellCells.
label newIndex = 0;
labelHashSet nbrCells;
if (cellCells.size() == 0)
{
return;
}
label startIndex = cellCells.offsets()[0];
forAll(cellCells, cellI)
{
nbrCells.clear();
nbrCells.insert(globalAgglom.toGlobal(cellI));
label endIndex = cellCells.offsets()[cellI+1];
for (label i = startIndex; i < endIndex; i++)
{
if (nbrCells.insert(cellCells.m()[i]))
{
cellCells.m()[newIndex++] = cellCells.m()[i];
}
}
startIndex = endIndex;
cellCells.offsets()[cellI+1] = newIndex;
}
cellCells.m().setSize(newIndex);
//forAll(cellCells, cellI)
//{
// Pout<< "Original: Coarse cell " << cellI << endl;
// forAll(mesh.cellCells()[cellI], i)
// {
// Pout<< " nbr:" << mesh.cellCells()[cellI][i] << endl;
// }
// Pout<< "Compacted: Coarse cell " << cellI << endl;
// const labelUList cCells = cellCells[cellI];
// forAll(cCells, i)
// {
// Pout<< " nbr:" << cCells[i] << endl;
// }
//}
}
//void Foam::decompositionMethod::calcCellCells
//(
// const polyMesh& mesh,
// const boolList& blockedFace,
// const List<labelPair>& explicitConnections,
// const labelList& agglom,
// const label nCoarse,
// const bool parallel,
// CompactListList<label>& cellCells
//)
//{
// const labelList& faceOwner = mesh.faceOwner();
// const labelList& faceNeighbour = mesh.faceNeighbour();
// const polyBoundaryMesh& patches = mesh.boundaryMesh();
//
//
// // Create global cell numbers
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// globalIndex globalAgglom
// (
// nCoarse,
// Pstream::msgType(),
// Pstream::worldComm,
// parallel
// );
//
//
// // Get agglomerate owner on other side of coupled faces
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
//
// forAll(patches, patchI)
// {
// const polyPatch& pp = patches[patchI];
//
// if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
// {
// label faceI = pp.start();
// label bFaceI = pp.start() - mesh.nInternalFaces();
//
// forAll(pp, i)
// {
// globalNeighbour[bFaceI] = globalAgglom.toGlobal
// (
// agglom[faceOwner[faceI]]
// );
//
// bFaceI++;
// faceI++;
// }
// }
// }
//
// // Get the cell on the other side of coupled patches
// syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
//
//
// // Count number of faces (internal + coupled)
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// // Number of faces per coarse cell
// labelList nFacesPerCell(nCoarse, 0);
//
// // 1. Internal faces
// for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
// {
// if (!blockedFace[faceI])
// {
// label own = agglom[faceOwner[faceI]];
// label nei = agglom[faceNeighbour[faceI]];
//
// nFacesPerCell[own]++;
// nFacesPerCell[nei]++;
// }
// }
//
// // 2. Coupled faces
// forAll(patches, patchI)
// {
// const polyPatch& pp = patches[patchI];
//
// if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
// {
// label faceI = pp.start();
// label bFaceI = pp.start()-mesh.nInternalFaces();
//
// forAll(pp, i)
// {
// if (!blockedFace[faceI])
// {
// label own = agglom[faceOwner[faceI]];
//
// label globalNei = globalNeighbour[bFaceI];
// if
// (
// !globalAgglom.isLocal(globalNei)
// || globalAgglom.toLocal(globalNei) != own
// )
// {
// nFacesPerCell[own]++;
// }
//
// faceI++;
// bFaceI++;
// }
// }
// }
// }
//
// // 3. Explicit connections between non-coupled boundary faces
// forAll(explicitConnections, i)
// {
// const labelPair& baffle = explicitConnections[i];
// label f0 = baffle.first();
// label f1 = baffle.second();
//
// if (!blockedFace[f0] && blockedFace[f1])
// {
// label f0Own = agglom[faceOwner[f0]];
// label f1Own = agglom[faceOwner[f1]];
//
// // Always count the connection between the two owner sides
// if (f0Own != f1Own)
// {
// nFacesPerCell[f0Own]++;
// nFacesPerCell[f1Own]++;
// }
//
// // Add any neighbour side connections
// if (mesh.isInternalFace(f0))
// {
// label f0Nei = agglom[faceNeighbour[f0]];
//
// if (mesh.isInternalFace(f1))
// {
// // Internal faces
// label f1Nei = agglom[faceNeighbour[f1]];
//
// if (f0Own != f1Nei)
// {
// nFacesPerCell[f0Own]++;
// nFacesPerCell[f1Nei]++;
// }
// if (f0Nei != f1Own)
// {
// nFacesPerCell[f0Nei]++;
// nFacesPerCell[f1Own]++;
// }
// if (f0Nei != f1Nei)
// {
// nFacesPerCell[f0Nei]++;
// nFacesPerCell[f1Nei]++;
// }
// }
// else
// {
// // f1 boundary face
// if (f0Nei != f1Own)
// {
// nFacesPerCell[f0Nei]++;
// nFacesPerCell[f1Own]++;
// }
// }
// }
// else
// {
// if (mesh.isInternalFace(f1))
// {
// label f1Nei = agglom[faceNeighbour[f1]];
// if (f0Own != f1Nei)
// {
// nFacesPerCell[f0Own]++;
// nFacesPerCell[f1Nei]++;
// }
// }
// }
// }
// }
//
//
// // Fill in offset and data
// // ~~~~~~~~~~~~~~~~~~~~~~~
//
// cellCells.setSize(nFacesPerCell);
//
// nFacesPerCell = 0;
//
// labelList& m = cellCells.m();
// const labelList& offsets = cellCells.offsets();
//
// // 1. For internal faces is just offsetted owner and neighbour
// for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
// {
// if (!blockedFace[faceI])
// {
// label own = agglom[faceOwner[faceI]];
// label nei = agglom[faceNeighbour[faceI]];
//
// m[offsets[own] + nFacesPerCell[own]++] =
// globalAgglom.toGlobal(nei);
// m[offsets[nei] + nFacesPerCell[nei]++] =
// globalAgglom.toGlobal(own);
// }
// }
//
// // 2. For boundary faces is offsetted coupled neighbour
// forAll(patches, patchI)
// {
// const polyPatch& pp = patches[patchI];
//
// if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
// {
// label faceI = pp.start();
// label bFaceI = pp.start()-mesh.nInternalFaces();
//
// forAll(pp, i)
// {
// if (!blockedFace[faceI])
// {
// label own = agglom[faceOwner[faceI]];
//
// label globalNei = globalNeighbour[bFaceI];
//
// if
// (
// !globalAgglom.isLocal(globalNei)
// || globalAgglom.toLocal(globalNei) != own
// )
// {
// m[offsets[own] + nFacesPerCell[own]++] = globalNei;
// }
//
// faceI++;
// bFaceI++;
// }
// }
// }
// }
//
// // 3. Explicit connections between non-coupled boundary faces
// forAll(explicitConnections, i)
// {
// const labelPair& baffle = explicitConnections[i];
// label f0 = baffle.first();
// label f1 = baffle.second();
//
// if (!blockedFace[f0] && blockedFace[f1])
// {
// label f0Own = agglom[faceOwner[f0]];
// label f1Own = agglom[faceOwner[f1]];
//
// // Always count the connection between the two owner sides
// if (f0Own != f1Own)
// {
// m[offsets[f0Own] + nFacesPerCell[f0Own]++] =
// globalAgglom.toGlobal(f1Own);
// m[offsets[f1Own] + nFacesPerCell[f1Own]++] =
// globalAgglom.toGlobal(f0Own);
// }
//
// // Add any neighbour side connections
// if (mesh.isInternalFace(f0))
// {
// label f0Nei = agglom[faceNeighbour[f0]];
//
// if (mesh.isInternalFace(f1))
// {
// // Internal faces
// label f1Nei = agglom[faceNeighbour[f1]];
//
// if (f0Own != f1Nei)
// {
// m[offsets[f0Own] + nFacesPerCell[f0Own]++] =
// globalAgglom.toGlobal(f1Nei);
// m[offsets[f1Nei] + nFacesPerCell[f1Nei]++] =
// globalAgglom.toGlobal(f1Nei);
// }
// if (f0Nei != f1Own)
// {
// m[offsets[f0Nei] + nFacesPerCell[f0Nei]++] =
// globalAgglom.toGlobal(f1Own);
// m[offsets[f1Own] + nFacesPerCell[f1Own]++] =
// globalAgglom.toGlobal(f0Nei);
// }
// if (f0Nei != f1Nei)
// {
// m[offsets[f0Nei] + nFacesPerCell[f0Nei]++] =
// globalAgglom.toGlobal(f1Nei);
// m[offsets[f1Nei] + nFacesPerCell[f1Nei]++] =
// globalAgglom.toGlobal(f0Nei);
// }
// }
// else
// {
// // f1 boundary face
// if (f0Nei != f1Own)
// {
// m[offsets[f0Nei] + nFacesPerCell[f0Nei]++] =
// globalAgglom.toGlobal(f1Own);
// m[offsets[f1Own] + nFacesPerCell[f1Own]++] =
// globalAgglom.toGlobal(f0Nei);
// }
// }
// }
// else
// {
// if (mesh.isInternalFace(f1))
// {
// label f1Nei = agglom[faceNeighbour[f1]];
// if (f0Own != f1Nei)
// {
// m[offsets[f0Own] + nFacesPerCell[f0Own]++] =
// globalAgglom.toGlobal(f1Nei);
// m[offsets[f1Nei] + nFacesPerCell[f1Nei]++] =
// globalAgglom.toGlobal(f0Own);
// }
// }
// }
// }
// }
//
//
// // Check for duplicates connections between cells
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// // Done as postprocessing step since we now have cellCells.
// label newIndex = 0;
// labelHashSet nbrCells;
//
//
// if (cellCells.size() == 0)
// {
// return;
// }
//
// label startIndex = cellCells.offsets()[0];
//
// forAll(cellCells, cellI)
// {
// nbrCells.clear();
// nbrCells.insert(globalAgglom.toGlobal(cellI));
//
// label endIndex = cellCells.offsets()[cellI+1];
//
// for (label i = startIndex; i < endIndex; i++)
// {
// if (nbrCells.insert(cellCells.m()[i]))
// {
// cellCells.m()[newIndex++] = cellCells.m()[i];
// }
// }
// startIndex = endIndex;
// cellCells.offsets()[cellI+1] = newIndex;
// }
//
// cellCells.m().setSize(newIndex);
//
// //forAll(cellCells, cellI)
// //{
// // Pout<< "Original: Coarse cell " << cellI << endl;
// // forAll(mesh.cellCells()[cellI], i)
// // {
// // Pout<< " nbr:" << mesh.cellCells()[cellI][i] << endl;
// // }
// // Pout<< "Compacted: Coarse cell " << cellI << endl;
// // const labelUList cCells = cellCells[cellI];
// // forAll(cCells, i)
// // {
// // Pout<< " nbr:" << cCells[i] << endl;
// // }
// //}
//}
Foam::labelList Foam::decompositionMethod::decompose
(
const polyMesh& mesh,
const scalarField& cellWeights,
//- Whether owner and neighbour should be on same processor
// (takes priority over explicitConnections)
const boolList& blockedFace,
//- Whether whole sets of faces (and point neighbours) need to be kept
// on single processor
const PtrList<labelList>& specifiedProcessorFaces,
const labelList& specifiedProcessor,
//- Additional connections between boundary faces
const List<labelPair>& explicitConnections
)
{
// Any weights specified?
label nWeights = returnReduce(cellWeights.size(), sumOp<label>());
// Any processor sets?
label nProcSets = 0;
forAll(specifiedProcessorFaces, setI)
{
nProcSets += specifiedProcessorFaces[setI].size();
}
reduce(nProcSets, sumOp<label>());
// Any non-mesh connections?
label nConnections = returnReduce
(
explicitConnections.size(),
sumOp<label>()
);
// Any faces not blocked?
label nUnblocked = 0;
forAll(blockedFace, faceI)
{
if (!blockedFace[faceI])
{
nUnblocked++;
}
}
reduce(nUnblocked, sumOp<label>());
// Either do decomposition on cell centres or on agglomeration
labelList finalDecomp;
if (nProcSets+nConnections+nUnblocked == 0)
{
// No constraints, possibly weights
if (nWeights > 0)
{
finalDecomp = decompose
(
mesh,
mesh.cellCentres(),
cellWeights
);
}
else
{
finalDecomp = decompose(mesh, mesh.cellCentres());
}
}
else
{
if (debug)
{
Info<< "Constrained decomposition:" << endl
<< " faces with same owner and neighbour processor : "
<< nUnblocked << endl
<< " baffle faces with same owner processor : "
<< nConnections << endl
<< " faces all on same processor : "
<< nProcSets << endl << endl;
}
// Determine global regions, separated by blockedFaces
regionSplit globalRegion(mesh, blockedFace, explicitConnections);
// Determine region cell centres
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// This just takes the first cell in the region. Otherwise the problem
// is with cyclics - if we'd average the region centre might be
// somewhere in the middle of the domain which might not be anywhere
// near any of the cells.
pointField regionCentres(globalRegion.nRegions(), point::max);
forAll(globalRegion, cellI)
{
label regionI = globalRegion[cellI];
if (regionCentres[regionI] == point::max)
{
regionCentres[regionI] = mesh.cellCentres()[cellI];
}
}
// Do decomposition on agglomeration
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
scalarField regionWeights(globalRegion.nRegions(), 0);
if (nWeights > 0)
{
forAll(globalRegion, cellI)
{
label regionI = globalRegion[cellI];
regionWeights[regionI] += cellWeights[cellI];
}
}
else
{
forAll(globalRegion, cellI)
{
label regionI = globalRegion[cellI];
regionWeights[regionI] += 1.0;
}
}
finalDecomp = decompose
(
mesh,
globalRegion,
regionCentres,
regionWeights
);
// Implement the explicitConnections since above decompose
// does not know about them
forAll(explicitConnections, i)
{
const labelPair& baffle = explicitConnections[i];
label f0 = baffle.first();
label f1 = baffle.second();
if (!blockedFace[f0] && !blockedFace[f1])
{
// Note: what if internal faces and owner and neighbour on
// different processor? So for now just push owner side
// proc
const label procI = finalDecomp[mesh.faceOwner()[f0]];
finalDecomp[mesh.faceOwner()[f1]] = procI;
if (mesh.isInternalFace(f1))
{
finalDecomp[mesh.faceNeighbour()[f1]] = procI;
}
}
else if (blockedFace[f0] != blockedFace[f1])
{
FatalErrorIn
(
"labelList decompose\n"
"(\n"
" const polyMesh&,\n"
" const scalarField&,\n"
" const boolList&,\n"
" const PtrList<labelList>&,\n"
" const labelList&,\n"
" const List<labelPair>&\n"
")"
) << "On explicit connection between faces " << f0
<< " and " << f1
<< " the two blockedFace status are not equal : "
<< blockedFace[f0] << " and " << blockedFace[f1]
<< exit(FatalError);
}
}
// For specifiedProcessorFaces rework the cellToProc to enforce
// all on one processor since we can't guarantee that the input
// to regionSplit was a single region.
// E.g. faceSet 'a' with the cells split into two regions
// by a notch formed by two walls
//
// \ /
// \ /
// ---a----+-----a-----
//
//
// Note that reworking the cellToProc might make the decomposition
// unbalanced.
forAll(specifiedProcessorFaces, setI)
{
const labelList& set = specifiedProcessorFaces[setI];
label procI = specifiedProcessor[setI];
if (procI == -1)
{
// If no processor specified use the one from the
// 0th element
procI = finalDecomp[mesh.faceOwner()[set[0]]];
}
forAll(set, fI)
{
const face& f = mesh.faces()[set[fI]];
forAll(f, fp)
{
const labelList& pFaces = mesh.pointFaces()[f[fp]];
forAll(pFaces, i)
{
label faceI = pFaces[i];
finalDecomp[mesh.faceOwner()[faceI]] = procI;
if (mesh.isInternalFace(faceI))
{
finalDecomp[mesh.faceNeighbour()[faceI]] = procI;
}
}
}
}
}
}
return finalDecomp;
}
void Foam::decompositionMethod::setConstraints
(
const polyMesh& mesh,
boolList& blockedFace,
PtrList<labelList>& specifiedProcessorFaces,
labelList& specifiedProcessor,
List<labelPair>& explicitConnections
)
{
blockedFace.setSize(mesh.nFaces());
blockedFace = true;
//label nUnblocked = 0;
specifiedProcessorFaces.clear();
explicitConnections.clear();
if (decompositionDict_.found("preservePatches"))
{
wordList pNames(decompositionDict_.lookup("preservePatches"));
Info<< nl
<< "Keeping owner of faces in patches " << pNames
<< " on same processor. This only makes sense for cyclics." << endl;
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(pNames, i)
{
const label patchI = patches.findPatchID(pNames[i]);
if (patchI == -1)
{
FatalErrorIn("decompositionMethod::decompose(const polyMesh&)")
<< "Unknown preservePatch " << pNames[i]
<< endl << "Valid patches are " << patches.names()
<< exit(FatalError);
}
const polyPatch& pp = patches[patchI];
forAll(pp, i)
{
if (blockedFace[pp.start() + i])
{
blockedFace[pp.start() + i] = false;
//nUnblocked++;
}
}
}
}
if (decompositionDict_.found("preserveFaceZones"))
{
wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
Info<< nl
<< "Keeping owner and neighbour of faces in zones " << zNames
<< " on same processor" << endl;
const faceZoneMesh& fZones = mesh.faceZones();
forAll(zNames, i)
{
label zoneI = fZones.findZoneID(zNames[i]);
if (zoneI == -1)
{
FatalErrorIn("decompositionMethod::decompose(const polyMesh&)")
<< "Unknown preserveFaceZone " << zNames[i]
<< endl << "Valid faceZones are " << fZones.names()
<< exit(FatalError);
}
const faceZone& fz = fZones[zoneI];
forAll(fz, i)
{
if (blockedFace[fz[i]])
{
blockedFace[fz[i]] = false;
//nUnblocked++;
}
}
}
}
bool preserveBaffles = decompositionDict_.lookupOrDefault
(
"preserveBaffles",
false
);
if (preserveBaffles)
{
Info<< nl
<< "Keeping owner of faces in baffles "
<< " on same processor." << endl;
// Faces to test: all boundary faces
labelList testFaces
(
identity(mesh.nFaces()-mesh.nInternalFaces())
+ mesh.nInternalFaces()
);
// Find correspondencing baffle face (or -1)
labelList duplicateFace
(
localPointRegion::findDuplicateFaces
(
mesh,
testFaces
)
);
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Convert into list of coupled face pairs (mesh face labels).
explicitConnections.setSize(testFaces.size());
label dupI = 0;
forAll(duplicateFace, i)
{
label otherFaceI = duplicateFace[i];
if (otherFaceI != -1 && i < otherFaceI)
{
label meshFace0 = testFaces[i];
label patch0 = patches.whichPatch(meshFace0);
label meshFace1 = testFaces[otherFaceI];
label patch1 = patches.whichPatch(meshFace1);
// Check for illegal topology. Should normally not happen!
if
(
(patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
|| (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
)
{
FatalErrorIn
(
"decompositionMethod::decompose(const polyMesh&)"
) << "One of two duplicate faces is on"
<< " processorPolyPatch."
<< "This is not allowed." << nl
<< "Face:" << meshFace0
<< " is on patch:" << patches[patch0].name()
<< nl
<< "Face:" << meshFace1
<< " is on patch:" << patches[patch1].name()
<< abort(FatalError);
}
explicitConnections[dupI++] = labelPair(meshFace0, meshFace1);
if (blockedFace[meshFace0])
{
blockedFace[meshFace0] = false;
//nUnblocked++;
}
if (blockedFace[meshFace1])
{
blockedFace[meshFace1] = false;
//nUnblocked++;
}
}
}
explicitConnections.setSize(dupI);
}
if
(
decompositionDict_.found("preservePatches")
|| decompositionDict_.found("preserveFaceZones")
|| preserveBaffles
)
{
syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
//reduce(nUnblocked, sumOp<label>());
}
// Specified processor for group of cells connected to faces
label nProcSets = 0;
if (decompositionDict_.found("singleProcessorFaceSets"))
{
List<Tuple2<word, label> > zNameAndProcs
(
decompositionDict_.lookup("singleProcessorFaceSets")
);
specifiedProcessorFaces.setSize(zNameAndProcs.size());
specifiedProcessor.setSize(zNameAndProcs.size());
forAll(zNameAndProcs, setI)
{
Info<< "Keeping all cells connected to faceSet "
<< zNameAndProcs[setI].first()
<< " on processor " << zNameAndProcs[setI].second() << endl;
// Read faceSet
faceSet fz(mesh, zNameAndProcs[setI].first());
specifiedProcessorFaces.set(setI, new labelList(fz.sortedToc()));
specifiedProcessor[setI] = zNameAndProcs[setI].second();
nProcSets += fz.size();
}
reduce(nProcSets, sumOp<label>());
// Unblock all point connected faces
// 1. Mark all points on specifiedProcessorFaces
boolList procFacePoint(mesh.nPoints(), false);
forAll(specifiedProcessorFaces, setI)
{
const labelList& set = specifiedProcessorFaces[setI];
forAll(set, fI)
{
const face& f = mesh.faces()[set[fI]];
forAll(f, fp)
{
procFacePoint[f[fp]] = true;
}
}
}
syncTools::syncPointList(mesh, procFacePoint, orEqOp<bool>(), false);
// 2. Unblock all faces on procFacePoint
forAll(procFacePoint, pointI)
{
if (procFacePoint[pointI])
{
const labelList& pFaces = mesh.pointFaces()[pointI];
forAll(pFaces, i)
{
blockedFace[pFaces[i]] = false;
}
}
}
syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
}
}
Foam::labelList Foam::decompositionMethod::decompose
(
const polyMesh& mesh,
const scalarField& cellWeights
)
{
//labelHashSet sameProcFaces;
//
//if (decompositionDict_.found("preservePatches"))
//{
// wordList pNames(decompositionDict_.lookup("preservePatches"));
//
// Info<< nl
// << "Keeping owner of faces in patches " << pNames
// << " on same processor. This only makes sense for cyclics."
// << endl;
//
// const polyBoundaryMesh& patches = mesh.boundaryMesh();
//
// forAll(pNames, i)
// {
// const label patchI = patches.findPatchID(pNames[i]);
//
// if (patchI == -1)
// {
// FatalErrorIn
// (
// "decompositionMethod::decompose(const polyMesh&)")
// << "Unknown preservePatch " << pNames[i]
// << endl << "Valid patches are " << patches.names()
// << exit(FatalError);
// }
//
// const polyPatch& pp = patches[patchI];
//
// forAll(pp, i)
// {
// sameProcFaces.insert(pp.start() + i);
// }
// }
//}
//if (decompositionDict_.found("preserveFaceZones"))
//{
// wordList zNames(decompositionDict_.lookup("preserveFaceZones"));
//
// Info<< nl
// << "Keeping owner and neighbour of faces in zones " << zNames
// << " on same processor" << endl;
//
// const faceZoneMesh& fZones = mesh.faceZones();
//
// forAll(zNames, i)
// {
// label zoneI = fZones.findZoneID(zNames[i]);
//
// if (zoneI == -1)
// {
// FatalErrorIn
// ("decompositionMethod::decompose(const polyMesh&)")
// << "Unknown preserveFaceZone " << zNames[i]
// << endl << "Valid faceZones are " << fZones.names()
// << exit(FatalError);
// }
//
// const faceZone& fz = fZones[zoneI];
//
// forAll(fz, i)
// {
// sameProcFaces.insert(fz[i]);
// }
// }
//}
//
//
//// Specified processor for group of cells connected to faces
//
////- Sets of faces to move together
//PtrList<labelList> specifiedProcessorFaces;
////- Destination processor
//labelList specifiedProcessor;
//
//label nProcSets = 0;
//if (decompositionDict_.found("singleProcessorFaceSets"))
//{
// List<Tuple2<word, label> > zNameAndProcs
// (
// decompositionDict_.lookup("singleProcessorFaceSets")
// );
//
// specifiedProcessorFaces.setSize(zNameAndProcs.size());
// specifiedProcessor.setSize(zNameAndProcs.size());
//
// forAll(zNameAndProcs, setI)
// {
// Info<< "Keeping all cells connected to faceSet "
// << zNameAndProcs[setI].first()
// << " on processor " << zNameAndProcs[setI].second() << endl;
//
// // Read faceSet
// faceSet fz(mesh, zNameAndProcs[setI].first());
//
// specifiedProcessorFaces.set(setI, new labelList(fz.sortedToc()));
// specifiedProcessor[setI] = zNameAndProcs[setI].second();
// nProcSets += fz.size();
// }
// reduce(nProcSets, sumOp<label>());
//}
//
//
//label nUnblocked = returnReduce(sameProcFaces.size(), sumOp<label>());
//
//if (nProcSets+nUnblocked > 0)
//{
// Info<< "Constrained decomposition:" << endl
// << " faces with same owner and neighbour processor : "
// << nUnblocked << endl
// << " faces all on same processor : "
// << nProcSets << endl << endl;
//}
//
//
//// Faces where owner and neighbour are not 'connected' (= all except
//// sameProcFaces)
//boolList blockedFace(mesh.nFaces(), true);
//{
// forAllConstIter(labelHashSet, sameProcFaces, iter)
// {
// blockedFace[iter.key()] = false;
// }
// syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
//
// // Add all point connected faces
// boolList procFacePoint(mesh.nPoints(), false);
// forAll(specifiedProcessorFaces, setI)
// {
// const labelList& set = specifiedProcessorFaces[setI];
// forAll(set, fI)
// {
// const face& f = mesh.faces()[set[fI]];
// forAll(f, fp)
// {
// procFacePoint[f[fp]] = true;
// }
// }
// }
// syncTools::syncPointList(mesh, procFacePoint, orEqOp<bool>(), false);
//
// forAll(procFacePoint, pointI)
// {
// if (procFacePoint[pointI])
// {
// const labelList& pFaces = mesh.pointFaces()[pointI];
// forAll(pFaces, i)
// {
// blockedFace[pFaces[i]] = false;
// }
// }
// }
// syncTools::syncFaceList(mesh, blockedFace, andEqOp<bool>());
//}
boolList blockedFace;
PtrList<labelList> specifiedProcessorFaces;
labelList specifiedProcessor;
List<labelPair> explicitConnections;
setConstraints
(
mesh,
blockedFace,
specifiedProcessorFaces,
specifiedProcessor,
explicitConnections
);
// Construct decomposition method and either do decomposition on
// cell centres or on agglomeration
labelList finalDecomp = decompose
(
mesh,
cellWeights, // optional weights
blockedFace, // any cells to be combined
specifiedProcessorFaces,// any whole cluster of cells to be kept
specifiedProcessor,
explicitConnections // baffles
);
return finalDecomp;
}
// ************************************************************************* //