ENH: mapDistribute : added transforms

This commit is contained in:
mattijs
2011-01-07 17:01:17 +00:00
parent 7243704068
commit 7a01f5c12c
20 changed files with 1333 additions and 965 deletions

View File

@ -1,3 +0,0 @@
Test-lduMatrix3.C
EXE = $(FOAM_USER_APPBIN)/Test-lduMatrix

View File

@ -1,4 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = -lfiniteVolume

View File

@ -1,163 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 "argList.H"
#include "Time.H"
#include "fvMesh.H"
#include "volFields.H"
#include "LduMatrix.H"
#include "diagTensorField.H"
#include "TPCG.H"
#include "TPBiCG.H"
#include "NoPreconditioner.H"
using namespace Foam;
typedef Foam::LduMatrix<vector, diagTensor, scalar>
lduVectorMatrix;
defineNamedTemplateTypeNameAndDebug(lduVectorMatrix, 0);
typedef Foam::DiagonalSolver<vector, diagTensor, scalar>
lduVectorDiagonalSolver;
defineNamedTemplateTypeNameAndDebug(lduVectorDiagonalSolver, 0);
template<>
const vector lduVectorMatrix::great_(1e15, 1e15, 1e15);
template<>
const vector lduVectorMatrix::small_(1e-15, 1e-15, 1e-15);
namespace Foam
{
typedef LduMatrix<vector, diagTensor, scalar>::preconditioner
lduVectorPreconditioner;
defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, symMatrix);
defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, asymMatrix);
typedef LduMatrix<vector, diagTensor, scalar>::smoother
lduVectorSmoother;
defineTemplateRunTimeSelectionTable(lduVectorSmoother, symMatrix);
defineTemplateRunTimeSelectionTable(lduVectorSmoother, asymMatrix);
typedef LduMatrix<vector, diagTensor, scalar>::solver
lduVectorSolver;
defineTemplateRunTimeSelectionTable(lduVectorSolver, symMatrix);
defineTemplateRunTimeSelectionTable(lduVectorSolver, asymMatrix);
typedef TPCG<vector, diagTensor, scalar> TPCGVector;
defineNamedTemplateTypeNameAndDebug(TPCGVector, 0);
LduMatrix<vector, diagTensor, scalar>::solver::
addsymMatrixConstructorToTable<TPCGVector>
addTPCGSymMatrixConstructorToTable_;
typedef TPBiCG<vector, diagTensor, scalar> TPBiCGVector;
defineNamedTemplateTypeNameAndDebug(TPBiCGVector, 0);
LduMatrix<vector, diagTensor, scalar>::solver::
addasymMatrixConstructorToTable<TPBiCGVector>
addTPBiCGSymMatrixConstructorToTable_;
typedef NoPreconditioner<vector, diagTensor, scalar> NoPreconditionerVector;
defineNamedTemplateTypeNameAndDebug(NoPreconditionerVector, 0);
LduMatrix<vector, diagTensor, scalar>::preconditioner::
addsymMatrixConstructorToTable<NoPreconditionerVector>
addNoPreconditionerSymMatrixConstructorToTable_;
LduMatrix<vector, diagTensor, scalar>::preconditioner::
addasymMatrixConstructorToTable<NoPreconditionerVector>
addNoPreconditionerAsymMatrixConstructorToTable_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
volVectorField psi
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
lduVectorMatrix testMatrix(mesh);
testMatrix.diag() = 2*pTraits<diagTensor>::one;
testMatrix.source() = pTraits<vector>::one;
testMatrix.upper() = 0.1;
testMatrix.lower() = -0.1;
Info<< testMatrix << endl;
FieldField<Field, scalar> boundaryCoeffs(0);
FieldField<Field, scalar> internalCoeffs(0);
autoPtr<lduVectorMatrix::solver> testMatrixSolver =
lduVectorMatrix::solver::New
(
psi.name(),
testMatrix,
boundaryCoeffs,
internalCoeffs,
psi.boundaryField().interfaces(),
IStringStream
(
"PBiCG"
"{"
" preconditioner none;"
" tolerance (1e-05 1e-05 1e-05);"
" relTol (0 0 0);"
"}"
)()
);
lduVectorMatrix::solverPerformance solverPerf =
testMatrixSolver->solve(psi);
solverPerf.print();
Info<< psi << endl;
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,163 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 "argList.H"
#include "Time.H"
#include "fvMesh.H"
#include "volFields.H"
#include "LduMatrix.H"
#include "tensorField.H"
#include "TPCG.H"
#include "TPBiCG.H"
#include "NoPreconditioner.H"
using namespace Foam;
typedef Foam::LduMatrix<vector, tensor, scalar>
lduVectorMatrix;
defineNamedTemplateTypeNameAndDebug(lduVectorMatrix, 0);
typedef Foam::DiagonalSolver<vector, tensor, scalar>
lduVectorDiagonalSolver;
defineNamedTemplateTypeNameAndDebug(lduVectorDiagonalSolver, 0);
template<>
const vector lduVectorMatrix::great_(1e15, 1e15, 1e15);
template<>
const vector lduVectorMatrix::small_(1e-15, 1e-15, 1e-15);
namespace Foam
{
typedef LduMatrix<vector, tensor, scalar>::preconditioner
lduVectorPreconditioner;
defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, symMatrix);
defineTemplateRunTimeSelectionTable(lduVectorPreconditioner, asymMatrix);
typedef LduMatrix<vector, tensor, scalar>::smoother
lduVectorSmoother;
defineTemplateRunTimeSelectionTable(lduVectorSmoother, symMatrix);
defineTemplateRunTimeSelectionTable(lduVectorSmoother, asymMatrix);
typedef LduMatrix<vector, tensor, scalar>::solver
lduVectorSolver;
defineTemplateRunTimeSelectionTable(lduVectorSolver, symMatrix);
defineTemplateRunTimeSelectionTable(lduVectorSolver, asymMatrix);
typedef TPCG<vector, tensor, scalar> TPCGVector;
defineNamedTemplateTypeNameAndDebug(TPCGVector, 0);
LduMatrix<vector, tensor, scalar>::solver::
addsymMatrixConstructorToTable<TPCGVector>
addTPCGSymMatrixConstructorToTable_;
typedef TPBiCG<vector, tensor, scalar> TPBiCGVector;
defineNamedTemplateTypeNameAndDebug(TPBiCGVector, 0);
LduMatrix<vector, tensor, scalar>::solver::
addasymMatrixConstructorToTable<TPBiCGVector>
addTPBiCGSymMatrixConstructorToTable_;
typedef NoPreconditioner<vector, tensor, scalar> NoPreconditionerVector;
defineNamedTemplateTypeNameAndDebug(NoPreconditionerVector, 0);
LduMatrix<vector, tensor, scalar>::preconditioner::
addsymMatrixConstructorToTable<NoPreconditionerVector>
addNoPreconditionerSymMatrixConstructorToTable_;
LduMatrix<vector, tensor, scalar>::preconditioner::
addasymMatrixConstructorToTable<NoPreconditionerVector>
addNoPreconditionerAsymMatrixConstructorToTable_;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
volVectorField psi
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
lduVectorMatrix testMatrix(mesh);
testMatrix.diag() = 2*I;
testMatrix.source() = pTraits<vector>::one;
testMatrix.upper() = 0.1;
testMatrix.lower() = -0.1;
Info<< testMatrix << endl;
FieldField<Field, scalar> boundaryCoeffs(0);
FieldField<Field, scalar> internalCoeffs(0);
autoPtr<lduVectorMatrix::solver> testMatrixSolver =
lduVectorMatrix::solver::New
(
psi.name(),
testMatrix,
//boundaryCoeffs,
//internalCoeffs,
//psi.boundaryField().interfaces(),
IStringStream
(
"PBiCG"
"{"
" preconditioner none;"
" tolerance (1e-05 1e-05 1e-05);"
" relTol (0 0 0);"
"}"
)()
);
lduVectorMatrix::solverPerformance solverPerf =
testMatrixSolver->solve(psi);
solverPerf.print();
Info<< psi << endl;
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,148 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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/>.
Application
icoFoam
Description
Transient solver for incompressible, laminar flow of Newtonian fluids.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "LduMatrix.H"
#include "diagTensorField.H"
typedef LduMatrix<vector, scalar, scalar> lduVectorMatrix;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "readPISOControls.H"
# include "CourantNo.H"
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);
fvVectorMatrix UEqnp(UEqn == -fvc::grad(p));
lduVectorMatrix U3Eqnp(mesh);
U3Eqnp.diag() = UEqnp.diag();
U3Eqnp.upper() = UEqnp.upper();
U3Eqnp.lower() = UEqnp.lower();
U3Eqnp.source() = UEqnp.source();
UEqnp.addBoundaryDiag(U3Eqnp.diag(), 0);
UEqnp.addBoundarySource(U3Eqnp.source(), false);
autoPtr<lduVectorMatrix::solver> U3EqnpSolver =
lduVectorMatrix::solver::New
(
U.name(),
U3Eqnp,
dictionary
(
IStringStream
(
"{"
" solver PBiCG;"
" preconditioner DILU;"
" tolerance (1e-13 1e-13 1e-13);"
" relTol (0 0 0);"
"}"
)()
)
);
U3EqnpSolver->solve(U).print(Info);
// --- PISO loop
for (int corr=0; corr<nCorr; corr++)
{
volScalarField rAU(1.0/UEqn.A());
U = rAU*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phi);
adjustPhi(phi, U, p);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
# include "continuityErrs.H"
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,56 +0,0 @@
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
dimensionedScalar nu
(
transportProperties.lookup("nu")
);
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -129,7 +129,7 @@ int main(int argc, char *argv[])
mapDistribute map(constructSize, sendMap.xfer(), recvMap.xfer());
// Distribute complexData
mapDistribute::distribute(complexData);
map.distribute(complexData);
Pout<< "complexData:" << complexData << endl;
}

View File

@ -42,6 +42,8 @@ primitives/SphericalTensor2D/sphericalTensor2D/sphericalTensor2D.C
primitives/Vector2D/vector2D/vector2D.C
primitives/complex/complex.C
primitives/globalIndexAndTransform/globalIndexAndTransform.C
primitives/globalIndexAndTransform/vectorTensorTransform/vectorTensorTransform.C
primitives/quaternion/quaternion.C
primitives/septernion/septernion.C

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -27,6 +27,11 @@ License
#include "commSchedule.H"
#include "HashSet.H"
#include "globalIndex.H"
#include "globalIndexAndTransform.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::mapDistribute, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -194,6 +199,309 @@ void Foam::mapDistribute::checkReceivedSize
}
void Foam::mapDistribute::printLayout
(
const label localSize,
List<Map<label> >& compactMap,
Ostream& os
) const
{
os << "Layout:" << endl
<< "local (processor " << Pstream::myProcNo() << "):" << endl
<< " start : 0" << endl
<< " size : " << localSize << endl;
label offset = localSize;
forAll(compactMap, procI)
{
if (procI != Pstream::myProcNo())
{
os << "processor " << procI << ':' << endl
<< " start : " << offset << endl
<< " size : " << compactMap[procI].size() << endl;
offset += compactMap[procI].size();
}
}
forAll(transformElements_, trafoI)
{
os << "transform " << trafoI << ':' << endl
<< " start : " << transformStart_[trafoI] << endl
<< " size : " << transformElements_[trafoI].size() << endl;
}
}
// Construct per processor compact addressing of the global elements
// needed. The ones from the local processor are not included since
// these are always all needed.
void Foam::mapDistribute::calcCompactAddressing
(
const globalIndex& globalNumbering,
const labelList& elements,
List<Map<label> >& compactMap
) const
{
compactMap.setSize(Pstream::nProcs());
// Count all (non-local) elements needed. Just for presizing map.
labelList nNonLocal(Pstream::nProcs(), 0);
forAll(elements, i)
{
label globalIndex = elements[i];
if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
{
label procI = globalNumbering.whichProcID(globalIndex);
nNonLocal[procI]++;
}
}
forAll(compactMap, procI)
{
compactMap[procI].clear();
if (procI != Pstream::myProcNo())
{
compactMap[procI].resize(2*nNonLocal[procI]);
}
}
// Collect all (non-local) elements needed.
forAll(elements, i)
{
label globalIndex = elements[i];
if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
{
label procI = globalNumbering.whichProcID(globalIndex);
label index = globalNumbering.toLocal(procI, globalIndex);
label nCompact = compactMap[procI].size();
compactMap[procI].insert(index, nCompact);
}
}
}
void Foam::mapDistribute::calcCompactAddressing
(
const globalIndex& globalNumbering,
const labelListList& cellCells,
List<Map<label> >& compactMap
) const
{
compactMap.setSize(Pstream::nProcs());
// Count all (non-local) elements needed. Just for presizing map.
labelList nNonLocal(Pstream::nProcs(), 0);
forAll(cellCells, cellI)
{
const labelList& cCells = cellCells[cellI];
forAll(cCells, i)
{
label globalIndex = cCells[i];
if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
{
label procI = globalNumbering.whichProcID(globalIndex);
nNonLocal[procI]++;
}
}
}
forAll(compactMap, procI)
{
compactMap[procI].clear();
if (procI != Pstream::myProcNo())
{
compactMap[procI].resize(2*nNonLocal[procI]);
}
}
// Collect all (non-local) elements needed.
forAll(cellCells, cellI)
{
const labelList& cCells = cellCells[cellI];
forAll(cCells, i)
{
label globalIndex = cCells[i];
if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
{
label procI = globalNumbering.whichProcID(globalIndex);
label index = globalNumbering.toLocal(procI, globalIndex);
label nCompact = compactMap[procI].size();
compactMap[procI].insert(index, nCompact);
}
}
}
}
void Foam::mapDistribute::exchangeAddressing
(
const globalIndex& globalNumbering,
labelList& elements,
List<Map<label> >& compactMap,
labelList& compactStart
)
{
// The overall compact addressing is
// - myProcNo data first (uncompacted)
// - all other processors consecutively
compactStart.setSize(Pstream::nProcs());
compactStart[Pstream::myProcNo()] = 0;
constructSize_ = globalNumbering.localSize();
forAll(compactStart, procI)
{
if (procI != Pstream::myProcNo())
{
compactStart[procI] = constructSize_;
constructSize_ += compactMap[procI].size();
}
}
// Find out what to receive/send in compact addressing.
// What I want to receive is what others have to send
labelListList wantedRemoteElements(Pstream::nProcs());
// Compact addressing for received data
constructMap_.setSize(Pstream::nProcs());
forAll(compactMap, procI)
{
if (procI == Pstream::myProcNo())
{
// All my own elements are used
label nLocal = globalNumbering.localSize();
wantedRemoteElements[procI] = identity(nLocal);
constructMap_[procI] = identity(nLocal);
}
else
{
// Remote elements wanted from processor procI
labelList& remoteElem = wantedRemoteElements[procI];
labelList& localElem = constructMap_[procI];
remoteElem.setSize(compactMap[procI].size());
localElem.setSize(compactMap[procI].size());
label i = 0;
forAllIter(Map<label>, compactMap[procI], iter)
{
const label compactI = compactStart[procI] + iter();
remoteElem[i] = iter.key();
localElem[i] = compactI;
iter() = compactI;
i++;
}
}
}
subMap_.setSize(Pstream::nProcs());
labelListList sendSizes;
Pstream::exchange<labelList, label>
(
wantedRemoteElements,
subMap_,
sendSizes
);
// Renumber elements
forAll(elements, i)
{
elements[i] = renumber(globalNumbering, compactMap, elements[i]);
}
}
void Foam::mapDistribute::exchangeAddressing
(
const globalIndex& globalNumbering,
labelListList& cellCells,
List<Map<label> >& compactMap,
labelList& compactStart
)
{
// The overall compact addressing is
// - myProcNo data first (uncompacted)
// - all other processors consecutively
compactStart.setSize(Pstream::nProcs());
compactStart[Pstream::myProcNo()] = 0;
constructSize_ = globalNumbering.localSize();
forAll(compactStart, procI)
{
if (procI != Pstream::myProcNo())
{
compactStart[procI] = constructSize_;
constructSize_ += compactMap[procI].size();
}
}
// Find out what to receive/send in compact addressing.
// What I want to receive is what others have to send
labelListList wantedRemoteElements(Pstream::nProcs());
// Compact addressing for received data
constructMap_.setSize(Pstream::nProcs());
forAll(compactMap, procI)
{
if (procI == Pstream::myProcNo())
{
// All my own elements are used
label nLocal = globalNumbering.localSize();
wantedRemoteElements[procI] = identity(nLocal);
constructMap_[procI] = identity(nLocal);
}
else
{
// Remote elements wanted from processor procI
labelList& remoteElem = wantedRemoteElements[procI];
labelList& localElem = constructMap_[procI];
remoteElem.setSize(compactMap[procI].size());
localElem.setSize(compactMap[procI].size());
label i = 0;
forAllIter(Map<label>, compactMap[procI], iter)
{
const label compactI = compactStart[procI] + iter();
remoteElem[i] = iter.key();
localElem[i] = compactI;
iter() = compactI;
i++;
}
}
}
subMap_.setSize(Pstream::nProcs());
labelListList sendSizes;
Pstream::exchange<labelList, label>
(
wantedRemoteElements,
subMap_,
sendSizes
);
// Renumber elements
forAll(cellCells, cellI)
{
labelList& cCells = cellCells[cellI];
forAll(cCells, i)
{
cCells[i] = renumber(globalNumbering, compactMap, cCells[i]);
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
//- Construct from components
@ -211,6 +519,25 @@ Foam::mapDistribute::mapDistribute
{}
//- Construct from components
Foam::mapDistribute::mapDistribute
(
const label constructSize,
const Xfer<labelListList>& subMap,
const Xfer<labelListList>& constructMap,
const Xfer<labelListList>& transformElements,
const Xfer<labelList>& transformStart
)
:
constructSize_(constructSize),
subMap_(subMap),
constructMap_(constructMap),
transformElements_(transformElements),
transformStart_(transformStart),
schedulePtr_()
{}
Foam::mapDistribute::mapDistribute
(
const labelList& sendProcs,
@ -295,133 +622,48 @@ Foam::mapDistribute::mapDistribute
constructSize_(0),
schedulePtr_()
{
// 1. Construct per processor compact addressing of the global elements
// needed. The ones from the local processor are not included since
// these are always all needed.
compactMap.setSize(Pstream::nProcs());
{
// Count all (non-local) elements needed. Just for presizing map.
labelList nNonLocal(Pstream::nProcs(), 0);
forAll(elements, i)
{
label globalIndex = elements[i];
if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
{
label procI = globalNumbering.whichProcID(globalIndex);
nNonLocal[procI]++;
}
}
forAll(compactMap, procI)
{
compactMap[procI].clear();
if (procI != Pstream::myProcNo())
{
compactMap[procI].resize(2*nNonLocal[procI]);
}
}
// Collect all (non-local) elements needed.
forAll(elements, i)
{
label globalIndex = elements[i];
if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
{
label procI = globalNumbering.whichProcID(globalIndex);
label index = globalNumbering.toLocal(procI, globalIndex);
label nCompact = compactMap[procI].size();
compactMap[procI].insert(index, nCompact);
}
}
//// Sort remote elements needed (not really necessary)
//forAll(compactMap, procI)
//{
// if (procI != Pstream::myProcNo())
// {
// Map<label>& globalMap = compactMap[procI];
//
// SortableList<label> sorted(globalMap.toc().xfer());
//
// forAll(sorted, i)
// {
// Map<label>::iterator iter = globalMap.find(sorted[i]);
// iter() = i;
// }
// }
//}
}
// 2. The overall compact addressing is
// - myProcNo data first (uncompacted)
// - all other processors consecutively
labelList compactStart(Pstream::nProcs());
compactStart[Pstream::myProcNo()] = 0;
constructSize_ = globalNumbering.localSize();
forAll(compactStart, procI)
{
if (procI != Pstream::myProcNo())
{
compactStart[procI] = constructSize_;
constructSize_ += compactMap[procI].size();
}
}
// 3. Find out what to receive/send in compact addressing.
// What I want to receive is what others have to send
labelListList wantedRemoteElements(Pstream::nProcs());
// Compact addressing for received data
constructMap_.setSize(Pstream::nProcs());
forAll(compactMap, procI)
{
if (procI == Pstream::myProcNo())
{
// All my own elements are used
label nLocal = globalNumbering.localSize();
wantedRemoteElements[procI] = identity(nLocal);
constructMap_[procI] = identity(nLocal);
}
else
{
// Remote elements wanted from processor procI
labelList& remoteElem = wantedRemoteElements[procI];
labelList& localElem = constructMap_[procI];
remoteElem.setSize(compactMap[procI].size());
localElem.setSize(compactMap[procI].size());
label i = 0;
forAllIter(Map<label>, compactMap[procI], iter)
{
const label compactI = compactStart[procI] + iter();
remoteElem[i] = iter.key();
localElem[i] = compactI;
iter() = compactI;
i++;
}
}
}
subMap_.setSize(Pstream::nProcs());
labelListList sendSizes;
Pstream::exchange<labelList, label>
// Construct per processor compact addressing of the global elements
// needed. The ones from the local processor are not included since
// these are always all needed.
calcCompactAddressing
(
wantedRemoteElements,
subMap_,
sendSizes
globalNumbering,
elements,
compactMap
);
// Renumber elements
forAll(elements, i)
//// Sort remote elements needed (not really necessary)
//forAll(compactMap, procI)
//{
// if (procI != Pstream::myProcNo())
// {
// Map<label>& globalMap = compactMap[procI];
//
// SortableList<label> sorted(globalMap.toc().xfer());
//
// forAll(sorted, i)
// {
// Map<label>::iterator iter = globalMap.find(sorted[i]);
// iter() = i;
// }
// }
//}
// Exchange what I need with processor that supplies it. Renumber elements
// into compact numbering
labelList compactStart;
exchangeAddressing
(
globalNumbering,
elements,
compactMap,
compactStart
);
if (debug)
{
elements[i] = renumber(globalNumbering, compactMap, elements[i]);
printLayout(globalNumbering.localSize(), compactMap, Pout);
}
}
@ -436,152 +678,272 @@ Foam::mapDistribute::mapDistribute
constructSize_(0),
schedulePtr_()
{
// 1. Construct per processor compact addressing of the global elements
// needed. The ones from the local processor are not included since
// these are always all needed.
compactMap.setSize(Pstream::nProcs());
{
// Count all (non-local) elements needed. Just for presizing map.
labelList nNonLocal(Pstream::nProcs(), 0);
forAll(cellCells, cellI)
{
const labelList& cCells = cellCells[cellI];
forAll(cCells, i)
{
label globalIndex = cCells[i];
if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
{
label procI = globalNumbering.whichProcID(globalIndex);
nNonLocal[procI]++;
}
}
}
forAll(compactMap, procI)
{
compactMap[procI].clear();
if (procI != Pstream::myProcNo())
{
compactMap[procI].resize(2*nNonLocal[procI]);
}
}
// Collect all (non-local) elements needed.
// Collect all (non-local) elements needed.
forAll(cellCells, cellI)
{
const labelList& cCells = cellCells[cellI];
forAll(cCells, i)
{
label globalIndex = cCells[i];
if (globalIndex != -1 && !globalNumbering.isLocal(globalIndex))
{
label procI = globalNumbering.whichProcID(globalIndex);
label index = globalNumbering.toLocal(procI, globalIndex);
label nCompact = compactMap[procI].size();
compactMap[procI].insert(index, nCompact);
}
}
}
//// Sort remote elements needed (not really necessary)
//forAll(compactMap, procI)
//{
// if (procI != Pstream::myProcNo())
// {
// Map<label>& globalMap = compactMap[procI];
//
// SortableList<label> sorted(globalMap.toc().xfer());
//
// forAll(sorted, i)
// {
// Map<label>::iterator iter = globalMap.find(sorted[i]);
// iter() = i;
// }
// }
//}
}
// 2. The overall compact addressing is
// - myProcNo data first (uncompacted)
// - all other processors consecutively
labelList compactStart(Pstream::nProcs());
compactStart[Pstream::myProcNo()] = 0;
constructSize_ = globalNumbering.localSize();
forAll(compactStart, procI)
{
if (procI != Pstream::myProcNo())
{
compactStart[procI] = constructSize_;
constructSize_ += compactMap[procI].size();
}
}
// 3. Find out what to receive/send in compact addressing.
// What I want to receive is what others have to send
labelListList wantedRemoteElements(Pstream::nProcs());
// Compact addressing for received data
constructMap_.setSize(Pstream::nProcs());
forAll(compactMap, procI)
{
if (procI == Pstream::myProcNo())
{
// All my own elements are used
label nLocal = globalNumbering.localSize();
wantedRemoteElements[procI] = identity(nLocal);
constructMap_[procI] = identity(nLocal);
}
else
{
// Remote elements wanted from processor procI
labelList& remoteElem = wantedRemoteElements[procI];
labelList& localElem = constructMap_[procI];
remoteElem.setSize(compactMap[procI].size());
localElem.setSize(compactMap[procI].size());
label i = 0;
forAllIter(Map<label>, compactMap[procI], iter)
{
const label compactI = compactStart[procI] + iter();
remoteElem[i] = iter.key();
localElem[i] = compactI;
iter() = compactI;
i++;
}
}
}
subMap_.setSize(Pstream::nProcs());
labelListList sendSizes;
Pstream::exchange<labelList, label>
// Construct per processor compact addressing of the global elements
// needed. The ones from the local processor are not included since
// these are always all needed.
calcCompactAddressing
(
wantedRemoteElements,
subMap_,
sendSizes
globalNumbering,
cellCells,
compactMap
);
// Renumber elements
forAll(cellCells, cellI)
{
labelList& cCells = cellCells[cellI];
//// Sort remote elements needed (not really necessary)
//forAll(compactMap, procI)
//{
// if (procI != Pstream::myProcNo())
// {
// Map<label>& globalMap = compactMap[procI];
//
// SortableList<label> sorted(globalMap.toc().xfer());
//
// forAll(sorted, i)
// {
// Map<label>::iterator iter = globalMap.find(sorted[i]);
// iter() = i;
// }
// }
//}
forAll(cCells, i)
// Exchange what I need with processor that supplies it. Renumber elements
// into compact numbering
labelList compactStart;
exchangeAddressing
(
globalNumbering,
cellCells,
compactMap,
compactStart
);
if (debug)
{
printLayout(globalNumbering.localSize(), compactMap, Pout);
}
}
Foam::mapDistribute::mapDistribute
(
const globalIndex& globalNumbering,
labelList& elements,
const globalIndexAndTransform& globalTransforms,
const labelPairList& transformedElements,
labelList& transformedIndices,
List<Map<label> >& compactMap
)
:
constructSize_(0),
schedulePtr_()
{
// Construct per processor compact addressing of the global elements
// needed. The ones from the local processor are not included since
// these are always all needed.
calcCompactAddressing
(
globalNumbering,
elements,
compactMap
);
// Add all (non-local) transformed elements needed.
forAll(transformedElements, i)
{
labelPair elem = transformedElements[i];
label procI = globalIndexAndTransform::processor(elem);
if (procI != Pstream::myProcNo())
{
cCells[i] = renumber(globalNumbering, compactMap, cCells[i]);
label index = globalIndexAndTransform::index(elem);
label nCompact = compactMap[procI].size();
compactMap[procI].insert(index, nCompact);
}
}
// Exchange what I need with processor that supplies it. Renumber elements
// into compact numbering
labelList compactStart;
exchangeAddressing
(
globalNumbering,
elements,
compactMap,
compactStart
);
// Renumber the transformed elements
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Count per transformIndex
label nTrafo = globalTransforms.transformPermutations().size();
labelList nPerTransform(nTrafo, 0);
forAll(transformedElements, i)
{
labelPair elem = transformedElements[i];
label trafoI = globalIndexAndTransform::transformIndex(elem);
nPerTransform[trafoI]++;
}
// Offset per transformIndex
transformStart_.setSize(nTrafo);
transformElements_.setSize(nTrafo);
forAll(transformStart_, trafoI)
{
transformStart_[trafoI] = constructSize_;
constructSize_ += nPerTransform[trafoI];
transformElements_[trafoI].setSize(nPerTransform[trafoI]);
}
// Sort transformed elements into their new slot.
nPerTransform = 0;
transformedIndices.setSize(transformedElements.size());
forAll(transformedElements, i)
{
labelPair elem = transformedElements[i];
label procI = globalIndexAndTransform::processor(elem);
label index = globalIndexAndTransform::index(elem);
label trafoI = globalIndexAndTransform::transformIndex(elem);
// Get compact index for untransformed element
label rawElemI =
(
procI == Pstream::myProcNo()
? index
: compactMap[procI][index]
);
label& n = nPerTransform[trafoI];
// index of element to transform
transformElements_[trafoI][n] = rawElemI;
// destination of transformed element
transformedIndices[i] = transformStart_[trafoI]+n;
n++;
}
if (debug)
{
printLayout(globalNumbering.localSize(), compactMap, Pout);
}
}
Foam::mapDistribute::mapDistribute
(
const globalIndex& globalNumbering,
labelListList& cellCells,
const globalIndexAndTransform& globalTransforms,
const List<labelPairList>& transformedElements,
labelListList& transformedIndices,
List<Map<label> >& compactMap
)
:
constructSize_(0),
schedulePtr_()
{
// Construct per processor compact addressing of the global elements
// needed. The ones from the local processor are not included since
// these are always all needed.
calcCompactAddressing
(
globalNumbering,
cellCells,
compactMap
);
// Add all (non-local) transformed elements needed.
forAll(transformedElements, cellI)
{
const labelPairList& elems = transformedElements[cellI];
forAll(elems, i)
{
label procI = globalIndexAndTransform::processor(elems[i]);
if (procI != Pstream::myProcNo())
{
label index = globalIndexAndTransform::index(elems[i]);
label nCompact = compactMap[procI].size();
compactMap[procI].insert(index, nCompact);
}
}
}
// Exchange what I need with processor that supplies it. Renumber elements
// into compact numbering
labelList compactStart;
exchangeAddressing
(
globalNumbering,
cellCells,
compactMap,
compactStart
);
// Renumber the transformed elements
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Count per transformIndex
label nTrafo = globalTransforms.transformPermutations().size();
labelList nPerTransform(nTrafo, 0);
forAll(transformedElements, cellI)
{
const labelPairList& elems = transformedElements[cellI];
forAll(elems, i)
{
label trafoI = globalIndexAndTransform::transformIndex(elems[i]);
nPerTransform[trafoI]++;
}
}
// Offset per transformIndex
transformStart_.setSize(nTrafo);
transformElements_.setSize(nTrafo);
forAll(transformStart_, trafoI)
{
transformStart_[trafoI] = constructSize_;
constructSize_ += nPerTransform[trafoI];
transformElements_[trafoI].setSize(nPerTransform[trafoI]);
}
// Sort transformed elements into their new slot.
nPerTransform = 0;
transformedIndices.setSize(transformedElements.size());
forAll(transformedElements, cellI)
{
const labelPairList& elems = transformedElements[cellI];
transformedIndices[cellI].setSize(elems.size());
forAll(elems, i)
{
label procI = globalIndexAndTransform::processor(elems[i]);
label index = globalIndexAndTransform::index(elems[i]);
label trafoI = globalIndexAndTransform::transformIndex(elems[i]);
// Get compact index for untransformed element
label rawElemI =
(
procI == Pstream::myProcNo()
? index
: compactMap[procI][index]
);
label& n = nPerTransform[trafoI];
// index of element to transform
transformElements_[trafoI][n] = rawElemI;
// destination of transformed element
transformedIndices[cellI][i] = transformStart_[trafoI]+n;
n++;
}
}
if (debug)
{
printLayout(globalNumbering.localSize(), compactMap, Pout);
}
}
@ -590,6 +952,8 @@ Foam::mapDistribute::mapDistribute(const mapDistribute& map)
constructSize_(map.constructSize_),
subMap_(map.subMap_),
constructMap_(map.constructMap_),
transformElements_(map.transformElements_),
transformStart_(map.transformStart_),
schedulePtr_()
{}
@ -774,6 +1138,8 @@ void Foam::mapDistribute::operator=(const mapDistribute& rhs)
constructSize_ = rhs.constructSize_;
subMap_ = rhs.subMap_;
constructMap_ = rhs.constructMap_;
transformElements_ = rhs.transformElements_;
transformStart_ = rhs.transformStart_;
schedulePtr_.clear();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -51,6 +51,7 @@ Note:
SourceFiles
mapDistribute.C
mapDistributeTemplates.C
\*---------------------------------------------------------------------------*/
@ -71,6 +72,7 @@ namespace Foam
class mapPolyMesh;
class globalIndex;
class PstreamBuffers;
class globalIndexAndTransform;
/*---------------------------------------------------------------------------*\
Class mapDistribute Declaration
@ -89,6 +91,16 @@ class mapDistribute
//- Maps from subsetted data to new reconstructed data
labelListList constructMap_;
// Optional transformation
//- For every globalIndexAndTransform::transformPermutations
// gives the elements that need to be transformed
labelListList transformElements_;
//- Destination in constructMap for transformed elements
labelList transformStart_;
//- Schedule
mutable autoPtr<List<labelPair> > schedulePtr_;
@ -102,8 +114,70 @@ class mapDistribute
const label receivedSize
);
void printLayout
(
const label localSize,
List<Map<label> >& compactMap,
Ostream& os
) const;
void calcCompactAddressing
(
const globalIndex& globalNumbering,
const labelList& elements,
List<Map<label> >& compactMap
) const;
void calcCompactAddressing
(
const globalIndex& globalNumbering,
const labelListList& elements,
List<Map<label> >& compactMap
) const;
void exchangeAddressing
(
const globalIndex& globalNumbering,
labelList& elements,
List<Map<label> >& compactMap,
labelList& compactStart
);
void exchangeAddressing
(
const globalIndex& globalNumbering,
labelListList& elements,
List<Map<label> >& compactMap,
labelList& compactStart
);
//- Helper function: copy transformElements without transformation
template<class T>
void copyElements(List<T>& field) const;
template<class T, class CombineOp>
void applyTransforms
(
const globalIndexAndTransform& globalTransforms,
List<T>& field,
const bool isPosition,
const CombineOp& cop
) const;
template<class T, class CombineOp>
void applyInverseTransforms
(
const globalIndexAndTransform& globalTransforms,
List<T>& field,
const bool isPosition,
const CombineOp& cop
) const;
public:
// Declare name of the class and its debug switch
ClassName("mapDistribute");
// Constructors
//- Construct from components
@ -114,6 +188,16 @@ public:
const Xfer<labelListList>& constructMap
);
//- Construct from components
mapDistribute
(
const label constructSize,
const Xfer<labelListList>& subMap,
const Xfer<labelListList>& constructMap,
const Xfer<labelListList>& transformElements,
const Xfer<labelList>& transformStart
);
//- Construct from reverse addressing: per data item the send
// processor and the receive processor. All processors get same data.
mapDistribute
@ -133,6 +217,20 @@ public:
List<Map<label> >& compactMap
);
//- Construct from list of (possibly) remote elements in globalIndex
// numbering (or -1). Determines compact numbering (see above) and
// distribute map to get data into this ordering and renumbers the
// elements to be in compact numbering.
mapDistribute
(
const globalIndex&,
labelList& untransformedElements,
const globalIndexAndTransform&,
const labelPairList& transformedElements,
labelList& transformedIndices,
List<Map<label> >& compactMap
);
//- Special variant that works with the info sorted into bins
// according to local indices. E.g. think cellCells where
// cellCells[localCellI] is a list of global cells
@ -143,6 +241,17 @@ public:
List<Map<label> >& compactMap
);
//- As above but with transformed elements as well
mapDistribute
(
const globalIndex&,
labelListList& cellCells,
const globalIndexAndTransform& ,
const List<labelPairList>& transformedElements,
labelListList& transformedIndices,
List<Map<label> >& compactMap
);
//- Construct copy
mapDistribute(const mapDistribute&);
@ -245,88 +354,30 @@ public:
//- Distribute data using default commsType.
template<class T>
void distribute(List<T>& fld) const
{
if (Pstream::defaultCommsType == Pstream::nonBlocking)
{
distribute
(
Pstream::nonBlocking,
List<labelPair>(),
constructSize_,
subMap_,
constructMap_,
fld
);
}
else if (Pstream::defaultCommsType == Pstream::scheduled)
{
distribute
(
Pstream::scheduled,
schedule(),
constructSize_,
subMap_,
constructMap_,
fld
);
}
else
{
distribute
(
Pstream::blocking,
List<labelPair>(),
constructSize_,
subMap_,
constructMap_,
fld
);
}
}
void distribute(List<T>& fld) const;
//- Same but with transforms
template<class T>
void distribute
(
const globalIndexAndTransform&,
List<T>& fld,
const bool isPosition
) const;
//- Reverse distribute data using default commsType.
template<class T>
void reverseDistribute(const label constructSize, List<T>& fld)
const
{
if (Pstream::defaultCommsType == Pstream::nonBlocking)
{
distribute
(
Pstream::nonBlocking,
List<labelPair>(),
constructSize,
constructMap_,
subMap_,
fld
);
}
else if (Pstream::defaultCommsType == Pstream::scheduled)
{
distribute
(
Pstream::scheduled,
schedule(),
constructSize,
constructMap_,
subMap_,
fld
);
}
else
{
distribute
(
Pstream::blocking,
List<labelPair>(),
constructSize,
constructMap_,
subMap_,
fld
);
}
}
void reverseDistribute(const label constructSize, List<T>&) const;
//- Same but with transforms
template<class T>
void reverseDistribute
(
const globalIndexAndTransform&,
const label constructSize,
List<T>& fld,
const bool isPosition
) const;
//- Reverse distribute data using default commsType.
// Since constructSize might be larger than supplied size supply
@ -337,52 +388,18 @@ public:
const label constructSize,
const T& nullValue,
List<T>& fld
)
const
{
if (Pstream::defaultCommsType == Pstream::nonBlocking)
{
distribute
(
Pstream::nonBlocking,
List<labelPair>(),
constructSize,
constructMap_,
subMap_,
fld,
eqOp<T>(),
nullValue
);
}
else if (Pstream::defaultCommsType == Pstream::scheduled)
{
distribute
(
Pstream::scheduled,
schedule(),
constructSize,
constructMap_,
subMap_,
fld,
eqOp<T>(),
nullValue
);
}
else
{
distribute
(
Pstream::blocking,
List<labelPair>(),
constructSize,
constructMap_,
subMap_,
fld,
eqOp<T>(),
nullValue
);
}
}
) const;
//- Same but with transforms
template<class T>
void reverseDistribute
(
const globalIndexAndTransform&,
const label constructSize,
const T& nullValue,
List<T>& fld,
const bool isPosition
) const;
//- Do all sends using PstreamBuffers
template<class T>

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,6 +26,7 @@ License
#include "Pstream.H"
#include "PstreamBuffers.H"
#include "PstreamCombineReduceOps.H"
#include "globalIndexAndTransform.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -774,4 +775,318 @@ void Foam::mapDistribute::receive(PstreamBuffers& pBufs, List<T>& field) const
}
// In case of no transform: copy elements
template<class T>
void Foam::mapDistribute::copyElements(List<T>& field) const
{
forAll(transformElements_, trafoI)
{
const labelList& elems = transformElements_[trafoI];
label n = transformStart_[trafoI];
forAll(elems, i)
{
field[n++] = field[elems[i]];
}
}
}
// Calculate transformed elements.
template<class T, class CombineOp>
void Foam::mapDistribute::applyTransforms
(
const globalIndexAndTransform& globalTransforms,
List<T>& field,
const bool isPosition,
const CombineOp& cop
) const
{
const List<vectorTensorTransform>& totalTransform =
globalTransforms.transformPermutations();
forAll(totalTransform, trafoI)
{
const vectorTensorTransform& vt = totalTransform[trafoI];
const labelList& elems = transformElements_[trafoI];
label n = transformStart_[trafoI];
// Could be optimised to avoid memory allocations
if (isPosition)
{
Field<T> transformFld(vt.transformPosition(Field<T>(field, elems)));
forAll(transformFld, i)
{
cop(field[n++], transformFld[i]);
}
}
else
{
Field<T> transformFld(transform(vt.R(), Field<T>(field, elems)));
forAll(transformFld, i)
{
cop(field[n++], transformFld[i]);
}
}
}
}
// Calculate transformed elements.
template<class T, class CombineOp>
void Foam::mapDistribute::applyInverseTransforms
(
const globalIndexAndTransform& globalTransforms,
List<T>& field,
const bool isPosition,
const CombineOp& cop
) const
{
const List<vectorTensorTransform>& totalTransform =
globalTransforms.transformPermutations();
forAll(totalTransform, trafoI)
{
const vectorTensorTransform& vt = totalTransform[trafoI];
const labelList& elems = transformElements_[trafoI];
label n = transformStart_[trafoI];
// Could be optimised to avoid memory allocations
if (isPosition)
{
Field<T> transformFld
(
vt.invTransformPosition
(
SubField<T>(field, n, elems.size())
)
);
forAll(transformFld, i)
{
cop(field[elems[i]], transformFld[i]);
}
}
else
{
Field<T> transformFld
(
transform
(
vt.R().T(),
SubField<T>(field, n, elems)
)
);
forAll(transformFld, i)
{
cop(field[elems[i]], transformFld[i]);
}
}
}
}
//- Distribute data using default commsType.
template<class T>
void Foam::mapDistribute::distribute(List<T>& fld) const
{
if (Pstream::defaultCommsType == Pstream::nonBlocking)
{
distribute
(
Pstream::nonBlocking,
List<labelPair>(),
constructSize_,
subMap_,
constructMap_,
fld
);
}
else if (Pstream::defaultCommsType == Pstream::scheduled)
{
distribute
(
Pstream::scheduled,
schedule(),
constructSize_,
subMap_,
constructMap_,
fld
);
}
else
{
distribute
(
Pstream::blocking,
List<labelPair>(),
constructSize_,
subMap_,
constructMap_,
fld
);
}
}
//- Reverse distribute data using default commsType.
template<class T>
void Foam::mapDistribute::reverseDistribute
(
const label constructSize,
List<T>& fld
) const
{
if (Pstream::defaultCommsType == Pstream::nonBlocking)
{
distribute
(
Pstream::nonBlocking,
List<labelPair>(),
constructSize,
constructMap_,
subMap_,
fld
);
}
else if (Pstream::defaultCommsType == Pstream::scheduled)
{
distribute
(
Pstream::scheduled,
schedule(),
constructSize,
constructMap_,
subMap_,
fld
);
}
else
{
distribute
(
Pstream::blocking,
List<labelPair>(),
constructSize,
constructMap_,
subMap_,
fld
);
}
}
//- Reverse distribute data using default commsType.
// Since constructSize might be larger than supplied size supply
// a nullValue
template<class T>
void Foam::mapDistribute::reverseDistribute
(
const label constructSize,
const T& nullValue,
List<T>& fld
) const
{
if (Pstream::defaultCommsType == Pstream::nonBlocking)
{
distribute
(
Pstream::nonBlocking,
List<labelPair>(),
constructSize,
constructMap_,
subMap_,
fld,
eqOp<T>(),
nullValue
);
}
else if (Pstream::defaultCommsType == Pstream::scheduled)
{
distribute
(
Pstream::scheduled,
schedule(),
constructSize,
constructMap_,
subMap_,
fld,
eqOp<T>(),
nullValue
);
}
else
{
distribute
(
Pstream::blocking,
List<labelPair>(),
constructSize,
constructMap_,
subMap_,
fld,
eqOp<T>(),
nullValue
);
}
}
//- Distribute data using default commsType.
template<class T>
void Foam::mapDistribute::distribute
(
const globalIndexAndTransform& git,
List<T>& fld,
const bool isPosition
) const
{
distribute(fld);
applyTransforms(git, fld, isPosition, eqOp<T>());
}
template<class T>
void Foam::mapDistribute::reverseDistribute
(
const globalIndexAndTransform& git,
const label constructSize,
List<T>& fld,
const bool isPosition
) const
{
fld.setSize(constructSize);
// Fill slots with reverse-transformed data
applyInverseTransforms(git, fld, isPosition, eqOp<T>());
// And send back
reverseDistribute(constructSize, fld);
}
template<class T>
void Foam::mapDistribute::reverseDistribute
(
const globalIndexAndTransform& git,
const label constructSize,
const T& nullValue,
List<T>& fld,
const bool isPosition
) const
{
fld = List<T>(constructSize, nullValue);
// Fill slots with reverse-transformed data
applyInverseTransforms(git, fld, isPosition, eqOp<T>());
// And send back
reverseDistribute(constructSize, fld);
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -324,6 +324,11 @@ void Foam::globalIndexAndTransform::determineTransformPermutations()
transformPermutations_[tPI] = transform;
}
// Encode index for 0 sign
labelList permutationIndices(nIndependentTransforms(), 0);
nullTransformIndex_ = encodeTransformIndex(permutationIndices);
}
@ -463,6 +468,7 @@ Foam::globalIndexAndTransform::globalIndexAndTransform
determinePatchTransformSign();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::globalIndexAndTransform::~globalIndexAndTransform()

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -42,11 +42,9 @@ SourceFiles
#ifndef globalIndexAndTransform_H
#define globalIndexAndTransform_H
#include "Pstream.H"
#include "List.H"
#include "labelPair.H"
#include "vectorTensorTransform.H"
#include "polyMesh.H"
#include "HashSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -82,6 +80,9 @@ class globalIndexAndTransform
// transform.
List<vectorTensorTransform> transformPermutations_;
//- Index of identity transform.
label nullTransformIndex_;
//- Mapping from patch index to which transform it matches (or
// -1 for none) (.first()) and what sign to us for it,
// i.e. +/- 1 (.second()).
@ -93,6 +94,7 @@ class globalIndexAndTransform
//- Number of spaces to reserve for transform encoding
static const label base_;
// Private Member Functions
//- Determine all of the independent basic transforms of the
@ -147,6 +149,16 @@ public:
const List<label>& permutationIndices
) const;
//- Add patch transformation to transformIndex. Return new
// transformIndex. (by default the patch is the sending, not the
// receiving, patch)
inline label addToTransformIndex
(
const label transformIndex,
const label patchI,
const bool isSendingSide = true
) const;
//- Encode index and bare index as components on own processor
inline labelPair encode
(
@ -163,13 +175,16 @@ public:
);
//- Index carried by the object
inline label index(const labelPair& globalIAndTransform);
inline static label index(const labelPair& globalIAndTransform);
//- Which processor does this come from?
inline label processor(const labelPair& globalIAndTransform);
inline static label processor(const labelPair& globalIAndTransform);
//- Transform carried by the object
inline label transformIndex(const labelPair& globalIAndTransform);
inline static label transformIndex
(
const labelPair& globalIAndTransform
);
// Access
@ -183,6 +198,10 @@ public:
inline const List<vectorTensorTransform>&
transformPermutations() const;
//- Return the transformIndex (index in transformPermutations)
// of the identity transform
inline label nullTransformIndex() const;
//- Return access to the per-patch transform-sign pairs
inline const List<Pair<label> >& patchTransformSign() const;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -23,8 +23,9 @@ License
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
#include "polyMesh.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::globalIndexAndTransform::encodeTransformIndex
(
@ -74,6 +75,90 @@ Foam::label Foam::globalIndexAndTransform::encodeTransformIndex
}
Foam::label Foam::globalIndexAndTransform::addToTransformIndex
(
const label transformIndex,
const label patchI,
const bool isSendingSide
) const
{
const Pair<label>& transSign = patchTransformSign_[patchI];
label matchTransI = transSign.first();
// Hardcoded for max 3 transforms only!
if (matchTransI > -1 && matchTransI < 3)
{
label t = transformIndex;
// Decode permutation as 3 integers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Note: FixedList for speed reasons.
FixedList<label, 3> permutation;
permutation[0] = (t%3)-1;
t /= 3;
permutation[1] = (t%3)-1;
t /= 3;
permutation[2] = (t%3)-1;
// Add patch transform
// ~~~~~~~~~~~~~~~~~~~
label sign = transSign.second();
if (!isSendingSide)
{
sign = -sign;
}
// If this transform been found already by a patch?
if (permutation[matchTransI] != 0)
{
// If so, if they have opposite signs, then this is
// considered an error. They are allowed to be the
// same sign, but this only results in a single
// transform.
if (permutation[matchTransI] != sign)
{
FatalErrorIn
(
"Foam::label "
"Foam::globalIndexAndTransform::addToTransformIndex\n"
"(\n"
"const label transformIndex,\n"
"const label patchI\n"
") const\n"
) << "More than one patch accessing the same transform "
<< "but not of the same sign." << endl
<< "patch:" << mesh_.boundaryMesh()[patchI].name()
<< " transform:" << matchTransI << " sign:" << sign
<< " cumulative transforms:" << permutation
<< exit(FatalError);
}
}
else
{
permutation[matchTransI] = sign;
}
// Re-ecode permutation
// ~~~~~~~~~~~~~~~~~~~~
return
(permutation[2]+1)*9
+ (permutation[1]+1)*3
+ (permutation[0]+1);
}
else
{
return transformIndex;
}
}
Foam::labelPair Foam::globalIndexAndTransform::encode
(
const label index,
@ -180,6 +265,12 @@ Foam::globalIndexAndTransform::transformPermutations() const
}
Foam::label Foam::globalIndexAndTransform::nullTransformIndex() const
{
return nullTransformIndex_;
}
const Foam::List<Foam::Pair<Foam::label> >&
Foam::globalIndexAndTransform::patchTransformSign() const
{
@ -399,7 +490,10 @@ Foam::pointField Foam::globalIndexAndTransform::transformPatches
forAll(transIs, tII)
{
transPts[tII] = transformPermutations_[transIs[tII]].transform(pt);
transPts[tII] = transformPermutations_[transIs[tII]].transformPosition
(
pt
);
}
return transPts;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -56,7 +56,7 @@ Foam::vectorTensorTransform::vectorTensorTransform(Istream& is)
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::word Foam::name(const vectorTensorTransform& s)
{
@ -68,6 +68,32 @@ Foam::word Foam::name(const vectorTensorTransform& s)
}
template<>
Foam::tmp<Foam::Field<bool> > Foam::vectorTensorTransform::transform
(
const Field<bool>& fld
) const
{
return fld;
}
template<>
Foam::tmp<Foam::Field<Foam::label> > Foam::vectorTensorTransform::transform
(
const Field<label>& fld
) const
{
return fld;
}
template<>
Foam::tmp<Foam::Field<Foam::scalar> > Foam::vectorTensorTransform::transform
(
const Field<scalar>& fld
) const
{
return fld;
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Istream& Foam::operator>>(Istream& is, vectorTensorTransform& tr)

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,13 +31,13 @@ Description
SourceFiles
vectorTensorTransformI.H
vectorTensorTransform.C
vectorTensorTransformTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef vectorTensorTransform_H
#define vectorTensorTransform_H
#include "vector.H"
#include "tensor.H"
#include "word.H"
#include "contiguous.H"
@ -132,17 +132,21 @@ public:
// Transform
//- Transform the given vector
inline vector transform(const vector& v) const;
//- Transform the given position
inline vector transformPosition(const vector& v) const;
//- Transform the given pointField
inline pointField transform(const pointField& pts) const;
inline pointField transformPosition(const pointField& pts) const;
//- Inverse transform the given vector
inline vector invTransform(const vector& v) const;
//- Inverse transform the given position
inline vector invTransformPosition(const vector& v) const;
//- Inverse transform the given pointField
inline pointField invTransform(const pointField& pts) const;
inline pointField invTransformPosition(const pointField& pts) const;
//- Transform the given field
template<class Type>
tmp<Field<Type> > transform(const Field<Type>&) const;
// Member operators
@ -180,6 +184,15 @@ word name(const vectorTensorTransform&);
template<>
inline bool contiguous<vectorTensorTransform>() {return true;}
//- Template specialisations
template<>
tmp<Field<bool> > vectorTensorTransform::transform(const Field<bool>&) const;
template<>
tmp<Field<label> > vectorTensorTransform::transform(const Field<label>&) const;
template<>
tmp<Field<scalar> > vectorTensorTransform::transform(const Field<scalar>&)
const;
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
@ -234,6 +247,10 @@ inline vectorTensorTransform operator&
#include "vectorTensorTransformI.H"
#ifdef NoRepository
# include "vectorTensorTransformTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2010-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -99,7 +99,7 @@ inline Foam::tensor& Foam::vectorTensorTransform::R()
}
inline Foam::vector Foam::vectorTensorTransform::transform
inline Foam::vector Foam::vectorTensorTransform::transformPosition
(
const vector& v
) const
@ -115,7 +115,7 @@ inline Foam::vector Foam::vectorTensorTransform::transform
}
inline Foam::pointField Foam::vectorTensorTransform::transform
inline Foam::pointField Foam::vectorTensorTransform::transformPosition
(
const pointField& pts
) const
@ -134,7 +134,7 @@ inline Foam::pointField Foam::vectorTensorTransform::transform
}
inline Foam::vector Foam::vectorTensorTransform::invTransform
inline Foam::vector Foam::vectorTensorTransform::invTransformPosition
(
const vector& v
) const
@ -150,7 +150,7 @@ inline Foam::vector Foam::vectorTensorTransform::invTransform
}
inline Foam::pointField Foam::vectorTensorTransform::invTransform
inline Foam::pointField Foam::vectorTensorTransform::invTransformPosition
(
const pointField& pts
) const

View File

@ -0,0 +1,45 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ 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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::vectorTensorTransform::transform
(
const Field<Type>& fld
) const
{
if (hasR_)
{
return R() & fld;
}
else
{
return fld;
}
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2008-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -180,7 +180,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
treeBoundBox tempTransformedBb
(
transform.invTransform(cellBbsToExchange[bbI].points())
transform.invTransformPosition(cellBbsToExchange[bbI].points())
);
treeBoundBox extendedBb
@ -397,7 +397,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
treeBoundBox tempTransformedBb
(
transform.invTransform(wallFaceBbsToExchange[bbI].points())
transform.invTransformPosition(wallFaceBbsToExchange[bbI].points())
);
treeBoundBox extendedBb
@ -535,7 +535,7 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
referredWallFaces_[rWFI] = referredWallFace
(
face(identity(f.size())),
transform.invTransform(f.points(mesh_.points())),
transform.invTransformPosition(f.points(mesh_.points())),
patchI
);
}
@ -700,7 +700,7 @@ void Foam::InteractionLists<ParticleType>::findExtendedProcBbsInRange
treeBoundBox extendedReferredProcBb
(
transform.transform
transform.transformPosition
(
allExtendedProcBbs[procI].points()
)
@ -747,7 +747,7 @@ void Foam::InteractionLists<ParticleType>::findExtendedProcBbsInRange
treeBoundBox extendedReferredProcBb
(
transform.transform
transform.transformPosition
(
allExtendedProcBbs[procI].points()
)
@ -790,7 +790,7 @@ void Foam::InteractionLists<ParticleType>::findExtendedProcBbsInRange
treeBoundBox extendedReferredProcBb
(
transform.transform
transform.transformPosition
(
allExtendedProcBbs[procI].points()
)
@ -957,7 +957,7 @@ void Foam::InteractionLists<ParticleType>::prepareParticleToBeReferred
globalTransforms_.transformIndex(ciat)
);
particle->position() = transform.invTransform(particle->position());
particle->position() = transform.invTransformPosition(particle->position());
particle->transformProperties(-transform.t());

View File

@ -4,8 +4,6 @@ indexedParticle = indexedParticle
$(passiveParticle)/passiveParticleCloud.C
$(indexedParticle)/indexedParticleCloud.C
InteractionLists/globalIndexAndTransform/globalIndexAndTransform.C
InteractionLists/globalIndexAndTransform/vectorTensorTransform/vectorTensorTransform.C
InteractionLists/referredWallFace/referredWallFace.C
LIB = $(FOAM_LIBBIN)/liblagrangian