1415 lines
42 KiB
C
1415 lines
42 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2013-2017 OpenFOAM Foundation
|
|
Copyright (C) 2019,2022 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
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 "lduPrimitiveMesh.H"
|
|
#include "processorLduInterface.H"
|
|
#include "edgeHashes.H"
|
|
#include "labelPair.H"
|
|
#include "processorGAMGInterface.H"
|
|
#include "cyclicAMIGAMGInterface.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
namespace Foam
|
|
{
|
|
defineTypeNameAndDebug(lduPrimitiveMesh, 0);
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
void Foam::lduPrimitiveMesh::checkUpperTriangular
|
|
(
|
|
const label size,
|
|
const labelUList& l,
|
|
const labelUList& u
|
|
)
|
|
{
|
|
forAll(l, facei)
|
|
{
|
|
if (u[facei] < l[facei])
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Reversed face. Problem at face " << facei
|
|
<< " l:" << l[facei] << " u:" << u[facei]
|
|
<< abort(FatalError);
|
|
}
|
|
if (l[facei] < 0 || u[facei] < 0 || u[facei] >= size)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Illegal cell label. Problem at face " << facei
|
|
<< " l:" << l[facei] << " u:" << u[facei]
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
|
|
for (label facei=1; facei < l.size(); ++facei)
|
|
{
|
|
if (l[facei-1] > l[facei])
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Lower not in incremental cell order."
|
|
<< " Problem at face " << facei
|
|
<< " l:" << l[facei] << " u:" << u[facei]
|
|
<< " previous l:" << l[facei-1]
|
|
<< abort(FatalError);
|
|
}
|
|
else if (l[facei-1] == l[facei])
|
|
{
|
|
// Same cell.
|
|
if (u[facei-1] > u[facei])
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Upper not in incremental cell order."
|
|
<< " Problem at face " << facei
|
|
<< " l:" << l[facei] << " u:" << u[facei]
|
|
<< " previous u:" << u[facei-1]
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
Foam::labelListList Foam::lduPrimitiveMesh::globalCellCells
|
|
(
|
|
const lduMesh& mesh,
|
|
const globalIndex& globalNumbering
|
|
)
|
|
{
|
|
const lduAddressing& addr = mesh.lduAddr();
|
|
lduInterfacePtrsList interfaces = mesh.interfaces();
|
|
|
|
const labelList globalIndices
|
|
(
|
|
identity
|
|
(
|
|
addr.size(),
|
|
globalNumbering.localStart(UPstream::myProcNo(mesh.comm()))
|
|
)
|
|
);
|
|
|
|
// Get the interface cells
|
|
PtrList<labelList> nbrGlobalCells(interfaces.size());
|
|
{
|
|
// Initialise transfer of restrict addressing on the interface
|
|
forAll(interfaces, inti)
|
|
{
|
|
if (interfaces.set(inti))
|
|
{
|
|
interfaces[inti].initInternalFieldTransfer
|
|
(
|
|
Pstream::commsTypes::nonBlocking,
|
|
globalIndices
|
|
);
|
|
}
|
|
}
|
|
|
|
// Wait for all
|
|
UPstream::waitRequests();
|
|
|
|
forAll(interfaces, inti)
|
|
{
|
|
if (interfaces.set(inti))
|
|
{
|
|
nbrGlobalCells.set
|
|
(
|
|
inti,
|
|
new labelList
|
|
(
|
|
interfaces[inti].internalFieldTransfer
|
|
(
|
|
Pstream::commsTypes::nonBlocking,
|
|
globalIndices
|
|
)
|
|
)
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Scan the neighbour list to find out how many times the cell
|
|
// appears as a neighbour of the face. Done this way to avoid guessing
|
|
// and resizing list
|
|
labelList nNbrs(addr.size(), Zero);
|
|
|
|
const labelUList& nbr = addr.upperAddr();
|
|
const labelUList& own = addr.lowerAddr();
|
|
|
|
{
|
|
forAll(nbr, facei)
|
|
{
|
|
nNbrs[nbr[facei]]++;
|
|
nNbrs[own[facei]]++;
|
|
}
|
|
|
|
forAll(interfaces, inti)
|
|
{
|
|
if (interfaces.set(inti))
|
|
{
|
|
for (const label celli : interfaces[inti].faceCells())
|
|
{
|
|
nNbrs[celli]++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Create cell-cells addressing
|
|
labelListList cellCells(addr.size());
|
|
|
|
forAll(cellCells, celli)
|
|
{
|
|
cellCells[celli].setSize(nNbrs[celli], -1);
|
|
}
|
|
|
|
// Reset the list of number of neighbours to zero
|
|
nNbrs = 0;
|
|
|
|
// Scatter the neighbour faces
|
|
forAll(nbr, facei)
|
|
{
|
|
const label c0 = own[facei];
|
|
const label c1 = nbr[facei];
|
|
|
|
cellCells[c0][nNbrs[c0]++] = globalIndices[c1];
|
|
cellCells[c1][nNbrs[c1]++] = globalIndices[c0];
|
|
}
|
|
forAll(interfaces, inti)
|
|
{
|
|
if (interfaces.set(inti))
|
|
{
|
|
const labelUList& faceCells = interfaces[inti].faceCells();
|
|
|
|
forAll(faceCells, facei)
|
|
{
|
|
const label c0 = faceCells[facei];
|
|
cellCells[c0][nNbrs[c0]++] = nbrGlobalCells[inti][facei];
|
|
}
|
|
}
|
|
}
|
|
|
|
return cellCells;
|
|
}
|
|
|
|
|
|
Foam::label Foam::lduPrimitiveMesh::totalSize
|
|
(
|
|
const PtrList<lduPrimitiveMesh>& meshes
|
|
)
|
|
{
|
|
label size = 0;
|
|
|
|
for (const lduPrimitiveMesh& msh : meshes)
|
|
{
|
|
size += msh.lduAddr().size();
|
|
}
|
|
return size;
|
|
}
|
|
|
|
|
|
Foam::labelList Foam::lduPrimitiveMesh::upperTriOrder
|
|
(
|
|
const label nCells,
|
|
const labelUList& lower,
|
|
const labelUList& upper
|
|
)
|
|
{
|
|
labelList nNbrs(nCells, Zero);
|
|
|
|
// Count number of upper neighbours
|
|
forAll(lower, facei)
|
|
{
|
|
if (upper[facei] < lower[facei])
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Problem at face:" << facei
|
|
<< " lower:" << lower[facei]
|
|
<< " upper:" << upper[facei]
|
|
<< exit(FatalError);
|
|
}
|
|
nNbrs[lower[facei]]++;
|
|
}
|
|
|
|
// Construct cell-upper cell addressing
|
|
labelList offsets(nCells+1);
|
|
offsets[0] = 0;
|
|
forAll(nNbrs, celli)
|
|
{
|
|
offsets[celli+1] = offsets[celli]+nNbrs[celli];
|
|
}
|
|
|
|
nNbrs = offsets;
|
|
|
|
labelList cellToFaces(offsets.last());
|
|
forAll(upper, facei)
|
|
{
|
|
const label celli = lower[facei];
|
|
cellToFaces[nNbrs[celli]++] = facei;
|
|
}
|
|
|
|
// Sort
|
|
|
|
labelList oldToNew(lower.size());
|
|
|
|
DynamicList<label> order;
|
|
DynamicList<label> nbr;
|
|
|
|
label newFacei = 0;
|
|
|
|
for (label celli = 0; celli < nCells; ++celli)
|
|
{
|
|
const label startOfCell = offsets[celli];
|
|
const label nNbr = offsets[celli+1] - startOfCell;
|
|
|
|
nbr.resize(nNbr);
|
|
order.resize(nNbr);
|
|
|
|
forAll(nbr, i)
|
|
{
|
|
nbr[i] = upper[cellToFaces[offsets[celli]+i]];
|
|
}
|
|
sortedOrder(nbr, order);
|
|
|
|
for (const label index : order)
|
|
{
|
|
oldToNew[cellToFaces[startOfCell + index]] = newFacei++;
|
|
}
|
|
}
|
|
|
|
return oldToNew;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
Foam::lduPrimitiveMesh::lduPrimitiveMesh
|
|
(
|
|
const label nCells,
|
|
labelList& l,
|
|
labelList& u,
|
|
const label comm,
|
|
bool reuse
|
|
)
|
|
:
|
|
lduAddressing(nCells),
|
|
lowerAddr_(l, reuse),
|
|
upperAddr_(u, reuse),
|
|
comm_(comm)
|
|
{
|
|
if (debug && lowerAddr_.size())
|
|
{
|
|
if (max(lowerAddr_) >= nCells || min(lowerAddr_) < 0)
|
|
{
|
|
FatalErrorInFunction << "Illegal lower addressing."
|
|
<< " nCells:" << nCells
|
|
<< " max(lower):" << max(lowerAddr_)
|
|
<< " min(lower):" << min(lowerAddr_)
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
if (debug && upperAddr_.size())
|
|
{
|
|
if (max(upperAddr_) >= nCells || min(upperAddr_) < 0)
|
|
{
|
|
FatalErrorInFunction << "Illegal upper addressing."
|
|
<< " nCells:" << nCells
|
|
<< " max(upper):" << max(upperAddr_)
|
|
<< " min(upper):" << min(upperAddr_)
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::lduPrimitiveMesh::addInterfaces
|
|
(
|
|
lduInterfacePtrsList& interfaces,
|
|
const lduSchedule& ps
|
|
)
|
|
{
|
|
interfaces_ = interfaces;
|
|
patchSchedule_ = ps;
|
|
|
|
// Create interfaces
|
|
primitiveInterfaces_.setSize(interfaces_.size());
|
|
forAll(interfaces_, i)
|
|
{
|
|
if (interfaces_.set(i))
|
|
{
|
|
primitiveInterfaces_.set(i, &interfaces_[i]);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
Foam::lduPrimitiveMesh::lduPrimitiveMesh
|
|
(
|
|
const label nCells
|
|
)
|
|
:
|
|
lduAddressing(nCells),
|
|
comm_(0)
|
|
{}
|
|
|
|
|
|
Foam::lduPrimitiveMesh::lduPrimitiveMesh
|
|
(
|
|
const label nCells,
|
|
labelList& l,
|
|
labelList& u,
|
|
PtrList<const lduInterface>& primitiveInterfaces,
|
|
const lduSchedule& ps,
|
|
const label comm
|
|
)
|
|
:
|
|
lduAddressing(nCells),
|
|
lowerAddr_(l, true),
|
|
upperAddr_(u, true),
|
|
primitiveInterfaces_(),
|
|
patchSchedule_(ps),
|
|
comm_(comm)
|
|
{
|
|
if (debug && lowerAddr_.size())
|
|
{
|
|
if (max(lowerAddr_) >= nCells || min(lowerAddr_) < 0)
|
|
{
|
|
FatalErrorInFunction << "Illegal lower addressing."
|
|
<< " nCells:" << nCells
|
|
<< " max(lower):" << max(lowerAddr_)
|
|
<< " min(lower):" << min(lowerAddr_)
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
if (debug && upperAddr_.size())
|
|
{
|
|
if (max(upperAddr_) >= nCells || min(upperAddr_) < 0)
|
|
{
|
|
FatalErrorInFunction << "Illegal upper addressing."
|
|
<< " nCells:" << nCells
|
|
<< " max(upper):" << max(upperAddr_)
|
|
<< " min(upper):" << min(upperAddr_)
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
|
|
|
|
primitiveInterfaces_.transfer(primitiveInterfaces);
|
|
|
|
// Create interfaces
|
|
interfaces_.setSize(primitiveInterfaces_.size());
|
|
forAll(primitiveInterfaces_, i)
|
|
{
|
|
if (primitiveInterfaces_.set(i))
|
|
{
|
|
interfaces_.set(i, &primitiveInterfaces_[i]);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
Foam::lduPrimitiveMesh::lduPrimitiveMesh
|
|
(
|
|
const label comm,
|
|
const labelList& procAgglomMap,
|
|
|
|
const labelList& procIDs,
|
|
const lduMesh& myMesh,
|
|
const PtrList<lduPrimitiveMesh>& otherMeshes,
|
|
|
|
labelList& cellOffsets,
|
|
labelList& faceOffsets,
|
|
labelListList& faceMap,
|
|
labelListList& boundaryMap,
|
|
labelListListList& boundaryFaceMap
|
|
)
|
|
:
|
|
lduAddressing(myMesh.lduAddr().size() + totalSize(otherMeshes)),
|
|
comm_(comm)
|
|
{
|
|
const label currentComm = myMesh.comm();
|
|
|
|
forAll(otherMeshes, i)
|
|
{
|
|
if (otherMeshes[i].comm() != currentComm)
|
|
{
|
|
WarningInFunction
|
|
<< "Communicator " << otherMeshes[i].comm()
|
|
<< " at index " << i
|
|
<< " differs from that of predecessor "
|
|
<< currentComm
|
|
<< endl; //exit(FatalError);
|
|
}
|
|
}
|
|
|
|
const label nMeshes = otherMeshes.size()+1;
|
|
|
|
const label myAgglom = procAgglomMap[UPstream::myProcNo(currentComm)];
|
|
|
|
if (lduPrimitiveMesh::debug)
|
|
{
|
|
Pout<< "I am " << UPstream::myProcNo(currentComm)
|
|
<< " agglomerating into " << myAgglom
|
|
<< " as are " << findIndices(procAgglomMap, myAgglom)
|
|
<< endl;
|
|
}
|
|
|
|
|
|
forAll(procIDs, i)
|
|
{
|
|
if (procAgglomMap[procIDs[i]] != procAgglomMap[procIDs[0]])
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Processor " << procIDs[i]
|
|
<< " agglomerates to " << procAgglomMap[procIDs[i]]
|
|
<< " whereas other processors " << procIDs
|
|
<< " agglomerate to "
|
|
<< labelUIndList(procAgglomMap, procIDs)
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
|
|
|
|
// Cells get added in order.
|
|
cellOffsets.setSize(nMeshes+1);
|
|
cellOffsets[0] = 0;
|
|
for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
|
|
{
|
|
const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
|
|
|
|
cellOffsets[procMeshI+1] =
|
|
cellOffsets[procMeshI]
|
|
+ procMesh.lduAddr().size();
|
|
}
|
|
|
|
|
|
// Faces initially get added in order, sorted later
|
|
labelList internalFaceOffsets(nMeshes+1);
|
|
internalFaceOffsets[0] = 0;
|
|
for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
|
|
{
|
|
const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
|
|
|
|
internalFaceOffsets[procMeshI+1] =
|
|
internalFaceOffsets[procMeshI]
|
|
+ procMesh.lduAddr().lowerAddr().size();
|
|
}
|
|
|
|
// Count how faces get added. Interfaces inbetween get merged.
|
|
|
|
// Merged interfaces: map from two coarse processors back to
|
|
// - procMeshes
|
|
// - interface in procMesh
|
|
// (estimate size from number of patches of mesh0)
|
|
EdgeMap<labelPairList> mergedMap(2*myMesh.interfaces().size());
|
|
|
|
// Unmerged interfaces: map from two coarse processors back to
|
|
// - procMeshes
|
|
// - interface in procMesh
|
|
EdgeMap<labelPairList> unmergedMap(mergedMap.size());
|
|
|
|
boundaryMap.setSize(nMeshes);
|
|
boundaryFaceMap.setSize(nMeshes);
|
|
|
|
|
|
label nOtherInterfaces = 0;
|
|
labelList nCoupledFaces(nMeshes, Zero);
|
|
|
|
for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
|
|
{
|
|
const lduInterfacePtrsList interfaces =
|
|
mesh(myMesh, otherMeshes, procMeshI).interfaces();
|
|
|
|
// Initialise all boundaries as merged
|
|
boundaryMap[procMeshI].setSize(interfaces.size(), -1);
|
|
boundaryFaceMap[procMeshI].setSize(interfaces.size());
|
|
|
|
// Get sorted order of processors
|
|
forAll(interfaces, intI)
|
|
{
|
|
if (interfaces.set(intI))
|
|
{
|
|
const lduInterface& ldui = interfaces[intI];
|
|
|
|
if (isA<processorLduInterface>(ldui))
|
|
{
|
|
const processorLduInterface& pldui =
|
|
refCast<const processorLduInterface>(ldui);
|
|
|
|
label agglom0 = procAgglomMap[pldui.myProcNo()];
|
|
label agglom1 = procAgglomMap[pldui.neighbProcNo()];
|
|
|
|
const edge procEdge(agglom0, agglom1);
|
|
|
|
if (agglom0 != myAgglom && agglom1 != myAgglom)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "At mesh from processor " << procIDs[procMeshI]
|
|
<< " have interface " << intI
|
|
<< " with myProcNo:" << pldui.myProcNo()
|
|
<< " with neighbProcNo:" << pldui.neighbProcNo()
|
|
<< exit(FatalError);
|
|
}
|
|
else if (agglom0 == myAgglom && agglom1 == myAgglom)
|
|
{
|
|
// Merged interface
|
|
if (debug)
|
|
{
|
|
Pout<< "merged interface: myProcNo:"
|
|
<< pldui.myProcNo()
|
|
<< " nbr:" << pldui.neighbProcNo()
|
|
<< " size:" << ldui.faceCells().size()
|
|
<< endl;
|
|
}
|
|
|
|
const label nbrProcMeshI =
|
|
procIDs.find(pldui.neighbProcNo());
|
|
|
|
if (procMeshI < nbrProcMeshI)
|
|
{
|
|
// I am 'master' since get lowest numbered cells
|
|
nCoupledFaces[procMeshI] += ldui.faceCells().size();
|
|
}
|
|
|
|
mergedMap(procEdge).append
|
|
(
|
|
labelPair(procMeshI, intI)
|
|
);
|
|
}
|
|
else
|
|
{
|
|
if (debug)
|
|
{
|
|
Pout<< "external interface: myProcNo:"
|
|
<< pldui.myProcNo()
|
|
<< " nbr:" << pldui.neighbProcNo()
|
|
<< " size:" << ldui.faceCells().size()
|
|
<< endl;
|
|
}
|
|
|
|
unmergedMap(procEdge).append
|
|
(
|
|
labelPair(procMeshI, intI)
|
|
);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Still external (non proc) interface
|
|
//FatalErrorInFunction
|
|
WarningInFunction
|
|
<< "At mesh from processor " << procIDs[procMeshI]
|
|
<< " have interface " << intI
|
|
<< " of unhandled type " << interfaces[intI].type()
|
|
//<< exit(FatalError);
|
|
<< endl;
|
|
|
|
++nOtherInterfaces;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "Remaining interfaces:" << endl;
|
|
forAllConstIters(unmergedMap, iter)
|
|
{
|
|
Pout<< " agglom procEdge:" << iter.key() << endl;
|
|
const labelPairList& elems = iter.val();
|
|
forAll(elems, i)
|
|
{
|
|
label procMeshI = elems[i][0];
|
|
label interfacei = elems[i][1];
|
|
const lduInterfacePtrsList interfaces =
|
|
mesh(myMesh, otherMeshes, procMeshI).interfaces();
|
|
|
|
const processorLduInterface& pldui =
|
|
refCast<const processorLduInterface>
|
|
(
|
|
interfaces[interfacei]
|
|
);
|
|
|
|
Pout<< " proc:" << procIDs[procMeshI]
|
|
<< " interfacei:" << interfacei
|
|
<< " between:" << pldui.myProcNo()
|
|
<< " and:" << pldui.neighbProcNo()
|
|
<< endl;
|
|
}
|
|
Pout<< endl;
|
|
}
|
|
}
|
|
if (debug)
|
|
{
|
|
Pout<< "Merged interfaces:" << endl;
|
|
forAllConstIters(mergedMap, iter)
|
|
{
|
|
Pout<< " agglom procEdge:" << iter.key() << endl;
|
|
const labelPairList& elems = iter.val();
|
|
|
|
forAll(elems, i)
|
|
{
|
|
label procMeshI = elems[i][0];
|
|
label interfacei = elems[i][1];
|
|
const lduInterfacePtrsList interfaces =
|
|
mesh(myMesh, otherMeshes, procMeshI).interfaces();
|
|
const processorLduInterface& pldui =
|
|
refCast<const processorLduInterface>
|
|
(
|
|
interfaces[interfacei]
|
|
);
|
|
|
|
Pout<< " proc:" << procIDs[procMeshI]
|
|
<< " interfacei:" << interfacei
|
|
<< " between:" << pldui.myProcNo()
|
|
<< " and:" << pldui.neighbProcNo()
|
|
<< endl;
|
|
}
|
|
Pout<< endl;
|
|
}
|
|
}
|
|
|
|
|
|
// Adapt faceOffsets for internal interfaces
|
|
faceOffsets.setSize(nMeshes+1);
|
|
faceOffsets[0] = 0;
|
|
faceMap.setSize(nMeshes);
|
|
for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
|
|
{
|
|
const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
|
|
label nInternal = procMesh.lduAddr().lowerAddr().size();
|
|
|
|
faceOffsets[procMeshI+1] =
|
|
faceOffsets[procMeshI]
|
|
+ nInternal
|
|
+ nCoupledFaces[procMeshI];
|
|
|
|
labelList& map = faceMap[procMeshI];
|
|
map.setSize(nInternal);
|
|
forAll(map, i)
|
|
{
|
|
map[i] = faceOffsets[procMeshI] + i;
|
|
}
|
|
}
|
|
|
|
|
|
// Combine upper and lower
|
|
lowerAddr_.setSize(faceOffsets.last(), -1);
|
|
upperAddr_.setSize(lowerAddr_.size(), -1);
|
|
|
|
|
|
// Old internal faces and resolved coupled interfaces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
|
|
{
|
|
const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
|
|
|
|
const labelUList& l = procMesh.lduAddr().lowerAddr();
|
|
const labelUList& u = procMesh.lduAddr().upperAddr();
|
|
|
|
// Add internal faces
|
|
label allFacei = faceOffsets[procMeshI];
|
|
|
|
forAll(l, facei)
|
|
{
|
|
lowerAddr_[allFacei] = cellOffsets[procMeshI]+l[facei];
|
|
upperAddr_[allFacei] = cellOffsets[procMeshI]+u[facei];
|
|
allFacei++;
|
|
}
|
|
|
|
|
|
// Add merged interfaces
|
|
const lduInterfacePtrsList interfaces = procMesh.interfaces();
|
|
|
|
forAll(interfaces, intI)
|
|
{
|
|
if (interfaces.set(intI))
|
|
{
|
|
const lduInterface& ldui = interfaces[intI];
|
|
|
|
if (isA<processorLduInterface>(ldui))
|
|
{
|
|
const processorLduInterface& pldui =
|
|
refCast<const processorLduInterface>(ldui);
|
|
|
|
// Look up corresponding interfaces
|
|
label myP = pldui.myProcNo();
|
|
label nbrP = pldui.neighbProcNo();
|
|
label nbrProcMeshI = procIDs.find(nbrP);
|
|
|
|
if (procMeshI < nbrProcMeshI)
|
|
{
|
|
// I am 'master' since my cell numbers will be lower
|
|
// since cells get added in procMeshI order.
|
|
|
|
const label agglom0 = procAgglomMap[myP];
|
|
const label agglom1 = procAgglomMap[nbrP];
|
|
|
|
const auto fnd =
|
|
mergedMap.cfind(edge(agglom0, agglom1));
|
|
|
|
if (fnd.found())
|
|
{
|
|
const labelPairList& elems = fnd();
|
|
|
|
// Find nbrP in elems
|
|
label nbrIntI = -1;
|
|
forAll(elems, i)
|
|
{
|
|
label proci = elems[i][0];
|
|
label interfacei = elems[i][1];
|
|
const lduInterfacePtrsList interfaces =
|
|
mesh
|
|
(
|
|
myMesh,
|
|
otherMeshes,
|
|
proci
|
|
).interfaces();
|
|
const processorLduInterface& pldui =
|
|
refCast<const processorLduInterface>
|
|
(
|
|
interfaces[interfacei]
|
|
);
|
|
|
|
if
|
|
(
|
|
elems[i][0] == nbrProcMeshI
|
|
&& pldui.neighbProcNo() == procIDs[procMeshI]
|
|
)
|
|
{
|
|
nbrIntI = elems[i][1];
|
|
break;
|
|
}
|
|
}
|
|
|
|
|
|
if (nbrIntI == -1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "elems:" << elems << abort(FatalError);
|
|
}
|
|
|
|
|
|
const lduInterfacePtrsList nbrInterfaces = mesh
|
|
(
|
|
myMesh,
|
|
otherMeshes,
|
|
nbrProcMeshI
|
|
).interfaces();
|
|
|
|
|
|
const labelUList& faceCells =
|
|
ldui.faceCells();
|
|
const labelUList& nbrFaceCells =
|
|
nbrInterfaces[nbrIntI].faceCells();
|
|
|
|
if (faceCells.size() != nbrFaceCells.size())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "faceCells:" << faceCells
|
|
<< " nbrFaceCells:" << nbrFaceCells
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
|
|
labelList& bfMap =
|
|
boundaryFaceMap[procMeshI][intI];
|
|
labelList& nbrBfMap =
|
|
boundaryFaceMap[nbrProcMeshI][nbrIntI];
|
|
|
|
bfMap.setSize(faceCells.size());
|
|
nbrBfMap.setSize(faceCells.size());
|
|
|
|
forAll(faceCells, pfI)
|
|
{
|
|
lowerAddr_[allFacei] =
|
|
cellOffsets[procMeshI]+faceCells[pfI];
|
|
bfMap[pfI] = allFacei;
|
|
upperAddr_[allFacei] =
|
|
cellOffsets[nbrProcMeshI]+nbrFaceCells[pfI];
|
|
nbrBfMap[pfI] = (-allFacei-1);
|
|
allFacei++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Sort upper-tri order
|
|
{
|
|
labelList oldToNew
|
|
(
|
|
upperTriOrder
|
|
(
|
|
cellOffsets.last(), //nCells
|
|
lowerAddr_,
|
|
upperAddr_
|
|
)
|
|
);
|
|
|
|
forAll(faceMap, procMeshI)
|
|
{
|
|
labelList& map = faceMap[procMeshI];
|
|
forAll(map, i)
|
|
{
|
|
if (map[i] >= 0)
|
|
{
|
|
map[i] = oldToNew[map[i]];
|
|
}
|
|
else
|
|
{
|
|
label allFacei = -map[i]-1;
|
|
map[i] = -oldToNew[allFacei]-1;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
inplaceReorder(oldToNew, lowerAddr_);
|
|
inplaceReorder(oldToNew, upperAddr_);
|
|
|
|
forAll(boundaryFaceMap, proci)
|
|
{
|
|
const labelList& bMap = boundaryMap[proci];
|
|
forAll(bMap, intI)
|
|
{
|
|
if (bMap[intI] == -1)
|
|
{
|
|
// Merged interface
|
|
labelList& bfMap = boundaryFaceMap[proci][intI];
|
|
|
|
forAll(bfMap, i)
|
|
{
|
|
if (bfMap[i] >= 0)
|
|
{
|
|
bfMap[i] = oldToNew[bfMap[i]];
|
|
}
|
|
else
|
|
{
|
|
label allFacei = -bfMap[i]-1;
|
|
bfMap[i] = (-oldToNew[allFacei]-1);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Kept interfaces
|
|
// ~~~~~~~~~~~~~~~
|
|
|
|
interfaces_.setSize(unmergedMap.size() + nOtherInterfaces);
|
|
primitiveInterfaces_.setSize(interfaces_.size());
|
|
|
|
label allInterfacei = 0;
|
|
|
|
//XXXXXXXX
|
|
// Do global interfaces
|
|
{
|
|
const lduInterfacePtrsList myInterfaces = myMesh.interfaces();
|
|
|
|
forAll(myInterfaces, intI)
|
|
{
|
|
if
|
|
(
|
|
myInterfaces.set(intI)
|
|
&& isA<cyclicAMIGAMGInterface>(myInterfaces[intI])
|
|
)
|
|
{
|
|
const auto& myIntf =
|
|
refCast<const GAMGInterface>(myInterfaces[intI]);
|
|
|
|
// Count number of faces on all processors
|
|
label n = myIntf.faceCells().size();
|
|
label nFine = myIntf.faceRestrictAddressing().size();
|
|
forAll(otherMeshes, otherI)
|
|
{
|
|
const lduInterfacePtrsList otherInterfaces =
|
|
otherMeshes[otherI].interfaces();
|
|
|
|
|
|
Pout<< " adding otherI:" << otherI
|
|
<< " with nFaces:"
|
|
<< otherInterfaces[intI].faceCells().size()
|
|
<< endl;
|
|
|
|
if
|
|
(
|
|
intI >= otherInterfaces.size()
|
|
|| !isA<cyclicAMIGAMGInterface>(otherInterfaces[intI])
|
|
)
|
|
{
|
|
FatalErrorInFunction<< "Other mesh " << otherI
|
|
<< " does not have global interface " << intI
|
|
<< " of type " << myIntf.type()
|
|
<< exit(FatalError);
|
|
}
|
|
const auto& otherIntf =
|
|
refCast<const GAMGInterface>(otherInterfaces[intI]);
|
|
n += otherIntf.faceCells().size();
|
|
nFine += otherIntf.faceRestrictAddressing().size();
|
|
}
|
|
|
|
// Size
|
|
labelField allFaceCells(n);
|
|
labelField allFaceRestrictAddressing(nFine);
|
|
labelField faceOffsets(nMeshes+1, -1);
|
|
|
|
n = 0;
|
|
nFine = 0;
|
|
|
|
// Fill
|
|
lduInterfacePtrsList allInterfaces(nMeshes);
|
|
|
|
for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
|
|
{
|
|
const auto& procMesh = mesh(myMesh, otherMeshes, procMeshI);
|
|
const auto interfaces = procMesh.interfaces();
|
|
|
|
faceOffsets[procMeshI] = n;
|
|
|
|
boundaryMap[procMeshI][intI] = allInterfacei;
|
|
labelList& bfMap = boundaryFaceMap[procMeshI][intI];
|
|
|
|
const labelUList& l = interfaces[intI].faceCells();
|
|
|
|
Pout<< " from procMesh:" << procMeshI
|
|
<< " from intI:" << intI
|
|
<< " with nFaces:" << l.size()
|
|
<< endl;
|
|
|
|
|
|
bfMap.setSize(l.size());
|
|
|
|
forAll(l, facei)
|
|
{
|
|
allFaceCells[n] = cellOffsets[procMeshI]+l[facei];
|
|
//allFaceRestrictAddressing[n] = n;
|
|
bfMap[facei] = n;
|
|
n++;
|
|
}
|
|
const auto& r = refCast<const GAMGInterface>
|
|
(
|
|
interfaces[intI]
|
|
).faceRestrictAddressing();
|
|
SubList<label>(allFaceRestrictAddressing, r.size(), nFine) =
|
|
r;
|
|
nFine += r.size();
|
|
|
|
allInterfaces.set(procMeshI, &interfaces[intI]);
|
|
}
|
|
faceOffsets.last() = n;
|
|
|
|
Pout<< "At interface:" << allInterfacei
|
|
<< " copying settings from global " << intI
|
|
<< " type " << myInterfaces[intI].type()
|
|
<< nl
|
|
<< " allFaceCells:" << flatOutput(allFaceCells) << nl
|
|
<< " allFaceRestr:"
|
|
<< flatOutput(allFaceRestrictAddressing) << nl
|
|
<< " faceOffsets:" << flatOutput(faceOffsets) << nl
|
|
<< endl;
|
|
|
|
|
|
// Construct new interface, remap local and remote contributions
|
|
primitiveInterfaces_.set
|
|
(
|
|
allInterfacei,
|
|
new cyclicAMIGAMGInterface
|
|
(
|
|
allInterfacei,
|
|
interfaces_,
|
|
allFaceCells,
|
|
allFaceRestrictAddressing,
|
|
faceOffsets,
|
|
myInterfaces[intI], // reference patch
|
|
allInterfaces, // all patches
|
|
comm_,
|
|
myAgglom,
|
|
procAgglomMap
|
|
)
|
|
);
|
|
interfaces_.set
|
|
(
|
|
allInterfacei,
|
|
&primitiveInterfaces_[allInterfacei]
|
|
);
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "Created " << interfaces_[allInterfacei].type()
|
|
<< " interface at " << allInterfacei
|
|
<< " comm:" << comm_
|
|
<< " myProcNo:" << myAgglom
|
|
<< " procAgglomMap:" << flatOutput(procAgglomMap)
|
|
<< " nFaces:" << allFaceCells.size()
|
|
<< endl;
|
|
}
|
|
|
|
|
|
allInterfacei++;
|
|
}
|
|
}
|
|
}
|
|
//XXXXXXXX
|
|
|
|
|
|
|
|
forAllConstIters(unmergedMap, iter)
|
|
{
|
|
const labelPairList& elems = iter.val();
|
|
|
|
// Sort processors in increasing order so both sides walk through in
|
|
// same order.
|
|
labelPairList procPairs(elems.size());
|
|
forAll(elems, i)
|
|
{
|
|
const labelPair& elem = elems[i];
|
|
label procMeshI = elem[0];
|
|
label interfacei = elem[1];
|
|
const lduInterfacePtrsList interfaces = mesh
|
|
(
|
|
myMesh,
|
|
otherMeshes,
|
|
procMeshI
|
|
).interfaces();
|
|
|
|
const processorLduInterface& pldui =
|
|
refCast<const processorLduInterface>
|
|
(
|
|
interfaces[interfacei]
|
|
);
|
|
label myProcNo = pldui.myProcNo();
|
|
label nbrProcNo = pldui.neighbProcNo();
|
|
procPairs[i] = labelPair
|
|
(
|
|
min(myProcNo, nbrProcNo),
|
|
max(myProcNo, nbrProcNo)
|
|
);
|
|
}
|
|
|
|
labelList order(sortedOrder(procPairs));
|
|
|
|
// Count
|
|
label n = 0;
|
|
|
|
forAll(order, i)
|
|
{
|
|
const labelPair& elem = elems[order[i]];
|
|
label procMeshI = elem[0];
|
|
label interfacei = elem[1];
|
|
const lduInterfacePtrsList interfaces = mesh
|
|
(
|
|
myMesh,
|
|
otherMeshes,
|
|
procMeshI
|
|
).interfaces();
|
|
|
|
n += interfaces[interfacei].faceCells().size();
|
|
}
|
|
|
|
// Size
|
|
labelField allFaceCells(n);
|
|
labelField allFaceRestrictAddressing(n);
|
|
n = 0;
|
|
|
|
// Fill
|
|
forAll(order, i)
|
|
{
|
|
const labelPair& elem = elems[order[i]];
|
|
label procMeshI = elem[0];
|
|
label interfacei = elem[1];
|
|
const lduInterfacePtrsList interfaces = mesh
|
|
(
|
|
myMesh,
|
|
otherMeshes,
|
|
procMeshI
|
|
).interfaces();
|
|
|
|
boundaryMap[procMeshI][interfacei] = allInterfacei;
|
|
labelList& bfMap = boundaryFaceMap[procMeshI][interfacei];
|
|
|
|
const labelUList& l = interfaces[interfacei].faceCells();
|
|
bfMap.setSize(l.size());
|
|
|
|
forAll(l, facei)
|
|
{
|
|
allFaceCells[n] = cellOffsets[procMeshI]+l[facei];
|
|
allFaceRestrictAddressing[n] = n;
|
|
bfMap[facei] = n;
|
|
n++;
|
|
}
|
|
}
|
|
|
|
|
|
// Find out local and remote processor in new communicator
|
|
|
|
label neighbProcNo = -1;
|
|
|
|
// See what the two processors map onto
|
|
|
|
if (iter.key()[0] == myAgglom)
|
|
{
|
|
if (iter.key()[1] == myAgglom)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "problem procEdge:" << iter.key()
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
neighbProcNo = iter.key()[1];
|
|
}
|
|
else
|
|
{
|
|
if (iter.key()[1] != myAgglom)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "problem procEdge:" << iter.key()
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
neighbProcNo = iter.key()[0];
|
|
}
|
|
|
|
primitiveInterfaces_.set
|
|
(
|
|
allInterfacei,
|
|
new processorGAMGInterface
|
|
(
|
|
allInterfacei,
|
|
interfaces_,
|
|
allFaceCells,
|
|
allFaceRestrictAddressing,
|
|
comm_,
|
|
myAgglom,
|
|
neighbProcNo,
|
|
tensorField(), // forwardT
|
|
Pstream::msgType() // tag
|
|
)
|
|
);
|
|
interfaces_.set(allInterfacei, &primitiveInterfaces_[allInterfacei]);
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "Created " << interfaces_[allInterfacei].type()
|
|
<< " interface at " << allInterfacei
|
|
<< " comm:" << comm_
|
|
<< " myProcNo:" << myAgglom
|
|
<< " neighbProcNo:" << neighbProcNo
|
|
<< " nFaces:" << allFaceCells.size()
|
|
<< endl;
|
|
}
|
|
|
|
|
|
allInterfacei++;
|
|
}
|
|
|
|
|
|
patchSchedule_ = nonBlockingSchedule<processorGAMGInterface>(interfaces_);
|
|
|
|
if (debug)
|
|
{
|
|
checkUpperTriangular(cellOffsets.last(), lowerAddr_, upperAddr_);
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
const Foam::lduMesh& Foam::lduPrimitiveMesh::mesh
|
|
(
|
|
const lduMesh& myMesh,
|
|
const PtrList<lduPrimitiveMesh>& otherMeshes,
|
|
const label meshI
|
|
)
|
|
{
|
|
return (meshI == 0 ? myMesh : otherMeshes[meshI-1]);
|
|
}
|
|
|
|
|
|
void Foam::lduPrimitiveMesh::gather
|
|
(
|
|
const label comm,
|
|
const lduMesh& mesh,
|
|
const labelList& procIDs,
|
|
PtrList<lduPrimitiveMesh>& otherMeshes
|
|
)
|
|
{
|
|
// Force calculation of schedule (since does parallel comms)
|
|
(void)mesh.lduAddr().patchSchedule();
|
|
|
|
if (Pstream::myProcNo(comm) == procIDs[0])
|
|
{
|
|
// Master.
|
|
otherMeshes.setSize(procIDs.size()-1);
|
|
|
|
for (label i = 1; i < procIDs.size(); ++i)
|
|
{
|
|
Pout<< "lduPrimitiveMesh::gather : on master :"
|
|
<< " receiving from proc " << procIDs[i] << endl;
|
|
|
|
IPstream fromProc
|
|
(
|
|
Pstream::commsTypes::scheduled,
|
|
procIDs[i],
|
|
0, // bufSize
|
|
Pstream::msgType(),
|
|
comm
|
|
);
|
|
|
|
label nCells = readLabel(fromProc);
|
|
labelList lowerAddr(fromProc);
|
|
labelList upperAddr(fromProc);
|
|
boolList validInterface(fromProc);
|
|
|
|
|
|
// Construct mesh without interfaces
|
|
otherMeshes.set
|
|
(
|
|
i-1,
|
|
new lduPrimitiveMesh
|
|
(
|
|
nCells,
|
|
lowerAddr,
|
|
upperAddr,
|
|
comm,
|
|
true // reuse
|
|
)
|
|
);
|
|
|
|
|
|
// Construct GAMGInterfaces
|
|
lduInterfacePtrsList newInterfaces(validInterface.size());
|
|
forAll(validInterface, intI)
|
|
{
|
|
if (validInterface[intI])
|
|
{
|
|
word coupleType(fromProc);
|
|
|
|
Pout<< "lduPrimitiveMesh::gather : on master :"
|
|
<< " receiving patch " << intI
|
|
<< " type:" << coupleType << endl;
|
|
|
|
newInterfaces.set
|
|
(
|
|
intI,
|
|
GAMGInterface::New
|
|
(
|
|
coupleType,
|
|
intI,
|
|
otherMeshes[i-1].rawInterfaces(),
|
|
fromProc
|
|
).ptr()
|
|
);
|
|
Pout<< "lduPrimitiveMesh::gather : on master :"
|
|
<< " DONE receiving patch " << intI
|
|
<< " type:" << coupleType << endl;
|
|
}
|
|
}
|
|
|
|
otherMeshes[i-1].addInterfaces
|
|
(
|
|
newInterfaces,
|
|
nonBlockingSchedule<processorGAMGInterface>
|
|
(
|
|
newInterfaces
|
|
)
|
|
);
|
|
|
|
Pout<< "lduPrimitiveMesh::gather : on master :"
|
|
<< " DONE received mesh from proc " << procIDs[i] << endl;
|
|
}
|
|
}
|
|
else if (procIDs.found(Pstream::myProcNo(comm)))
|
|
{
|
|
// Send to master
|
|
|
|
const lduAddressing& addressing = mesh.lduAddr();
|
|
lduInterfacePtrsList interfaces(mesh.interfaces());
|
|
boolList validInterface(interfaces.size());
|
|
forAll(interfaces, intI)
|
|
{
|
|
validInterface[intI] = interfaces.set(intI);
|
|
}
|
|
|
|
OPstream toMaster
|
|
(
|
|
Pstream::commsTypes::scheduled,
|
|
procIDs[0],
|
|
0,
|
|
Pstream::msgType(),
|
|
comm
|
|
);
|
|
|
|
Pout<< "lduPrimitiveMesh::gather : on proc :" << Pstream::myProcNo(comm)
|
|
<< " sending to master"
|
|
<< " lowerAddr:" << flatOutput(addressing.lowerAddr())
|
|
<< " upperAddr:" << flatOutput(addressing.upperAddr())
|
|
<< " validInterface:" << flatOutput(validInterface)
|
|
<< endl;
|
|
|
|
toMaster
|
|
<< addressing.size()
|
|
<< addressing.lowerAddr()
|
|
<< addressing.upperAddr()
|
|
<< validInterface;
|
|
|
|
forAll(interfaces, intI)
|
|
{
|
|
if (interfaces.set(intI))
|
|
{
|
|
const GAMGInterface& interface = refCast<const GAMGInterface>
|
|
(
|
|
interfaces[intI]
|
|
);
|
|
|
|
Pout<< " sending to master patch :" << intI
|
|
<< " of type:" << interface.type() << endl;
|
|
|
|
toMaster << interface.type();
|
|
interface.write(toMaster);
|
|
|
|
Pout<< " DONE sending to master patch :" << intI
|
|
<< " of type:" << interface.type() << endl;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|