globalIndexAndTransform: Support any number of transforms but no more than 3 per point

Patch contributed by Mattijs Janssens
Resolves bug-report http://bugs.openfoam.org/view.php?id=1815
This commit is contained in:
Henry Weller
2016-09-22 14:10:45 +01:00
parent 40f8709488
commit c0f841e4f7
8 changed files with 384 additions and 370 deletions

View File

@ -24,19 +24,14 @@ License
\*---------------------------------------------------------------------------*/
#include "globalMeshData.H"
#include "Time.H"
#include "Pstream.H"
#include "PstreamCombineReduceOps.H"
#include "processorPolyPatch.H"
#include "demandDrivenData.H"
#include "globalPoints.H"
#include "polyMesh.H"
#include "mapDistribute.H"
#include "labelIOList.H"
#include "PackedList.H"
#include "mergePoints.H"
#include "matchPoints.H"
#include "OFstream.H"
#include "globalIndexAndTransform.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -577,7 +572,7 @@ void Foam::globalMeshData::calcPointConnectivity
labelPairList myData(globalPointSlavesMap().constructSize());
forAll(slaves, pointi)
{
myData[pointi] = globalIndexAndTransform::encode
myData[pointi] = transforms.encode
(
Pstream::myProcNo(),
pointi,
@ -624,9 +619,9 @@ void Foam::globalMeshData::calcPointConnectivity
);
// Add transform to connectivity
const labelPair& n = myData[pTransformSlaves[i]];
label proci = globalIndexAndTransform::processor(n);
label index = globalIndexAndTransform::index(n);
pConnectivity[connI++] = globalIndexAndTransform::encode
label proci = transforms.processor(n);
label index = transforms.index(n);
pConnectivity[connI++] = transforms.encode
(
proci,
index,
@ -678,6 +673,8 @@ void Foam::globalMeshData::calcGlobalPointEdges
const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
const labelListList& slaves = globalPointSlaves();
const labelListList& transformedSlaves = globalPointTransformedSlaves();
const globalIndexAndTransform& transforms = globalTransforms();
// Create local version
globalPointEdges.setSize(globalPointSlavesMap().constructSize());
@ -697,11 +694,11 @@ void Foam::globalMeshData::calcGlobalPointEdges
forAll(pEdges, i)
{
label otherPointi = edges[pEdges[i]].otherVertex(pointi);
globalPPoints[i] = globalIndexAndTransform::encode
globalPPoints[i] = transforms.encode
(
Pstream::myProcNo(),
otherPointi,
globalTransforms().nullTransformIndex()
transforms.nullTransformIndex()
);
}
}
@ -790,9 +787,9 @@ void Foam::globalMeshData::calcGlobalPointEdges
{
// Add transform to connectivity
const labelPair& n = otherData[j];
label proci = globalIndexAndTransform::processor(n);
label index = globalIndexAndTransform::index(n);
globalPPoints[sz++] = globalIndexAndTransform::encode
label proci = transforms.processor(n);
label index = transforms.index(n);
globalPPoints[sz++] = transforms.encode
(
proci,
index,
@ -834,16 +831,18 @@ Foam::label Foam::globalMeshData::findTransform
const label localPoint
) const
{
const label remoteProci = globalIndexAndTransform::processor(remotePoint);
const label remoteIndex = globalIndexAndTransform::index(remotePoint);
const globalIndexAndTransform& transforms = globalTransforms();
const label remoteProci = transforms.processor(remotePoint);
const label remoteIndex = transforms.index(remotePoint);
label remoteTransformI = -1;
label localTransformI = -1;
forAll(info, i)
{
label proci = globalIndexAndTransform::processor(info[i]);
label pointi = globalIndexAndTransform::index(info[i]);
label transformI = globalIndexAndTransform::transformIndex(info[i]);
label proci = transforms.processor(info[i]);
label pointi = transforms.index(info[i]);
label transformI = transforms.transformIndex(info[i]);
if (proci == Pstream::myProcNo() && pointi == localPoint)
{
@ -875,7 +874,7 @@ Foam::label Foam::globalMeshData::findTransform
<< abort(FatalError);
}
return globalTransforms().subtractTransformIndex
return transforms.subtractTransformIndex
(
remoteTransformI,
localTransformI
@ -893,6 +892,7 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
const edgeList& edges = coupledPatch().edges();
const globalIndex& globalEdgeNumbers = globalEdgeNumbering();
const globalIndexAndTransform& transforms = globalTransforms();
// The whole problem with deducting edge-connectivity from
@ -941,11 +941,11 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
// Append myself.
eEdges.append
(
globalIndexAndTransform::encode
transforms.encode
(
Pstream::myProcNo(),
edgeI,
globalTransforms().nullTransformIndex()
transforms.nullTransformIndex()
)
);
@ -986,7 +986,7 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
label proci = globalEdgeNumbers.whichProcID(pEdges0[i]);
eEdges.append
(
globalIndexAndTransform::encode
transforms.encode
(
proci,
globalEdgeNumbers.toLocal(proci, pEdges0[i]),
@ -999,7 +999,11 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
}
allEdgeConnectivity[edgeI].transfer(eEdges);
sort(allEdgeConnectivity[edgeI], globalIndexAndTransform::less());
sort
(
allEdgeConnectivity[edgeI],
globalIndexAndTransform::less(transforms)
);
}
// We now have - in allEdgeConnectivity - a list of edges which are shared
@ -1020,10 +1024,10 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
if
(
(
globalIndexAndTransform::processor(masterInfo)
transforms.processor(masterInfo)
== Pstream::myProcNo()
)
&& (globalIndexAndTransform::index(masterInfo) == edgeI)
&& (transforms.index(masterInfo) == edgeI)
)
{
// Sort into transformed and untransformed
@ -1039,14 +1043,14 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
for (label i = 1; i < edgeInfo.size(); i++)
{
const labelPair& info = edgeInfo[i];
label proci = globalIndexAndTransform::processor(info);
label index = globalIndexAndTransform::index(info);
label transform = globalIndexAndTransform::transformIndex
label proci = transforms.processor(info);
label index = transforms.index(info);
label transform = transforms.transformIndex
(
info
);
if (transform == globalTransforms().nullTransformIndex())
if (transform == transforms.nullTransformIndex())
{
eEdges[nonTransformI++] = globalEdgeNumbers.toGlobal
(
@ -1078,7 +1082,7 @@ void Foam::globalMeshData::calcGlobalEdgeSlaves() const
globalEdgeNumbers,
globalEdgeSlaves,
globalTransforms(),
transforms,
transformedEdges,
globalEdgeTransformedSlavesPtr_(),
@ -1351,6 +1355,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
const labelListList& pointSlaves = globalPointSlaves();
const labelListList& pointTransformSlaves =
globalPointTransformedSlaves();
const globalIndexAndTransform& transforms = globalTransforms();
// Any faces coming in through transformation
@ -1432,7 +1437,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
label proci = globalIndices.whichProcID(slave);
label facei = globalIndices.toLocal(proci, slave);
myBFaces[n++] = globalIndexAndTransform::encode
myBFaces[n++] = transforms.encode
(
proci,
facei,
@ -1466,7 +1471,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
globalIndices,
globalPointBoundaryFaces,
globalTransforms(),
transforms,
transformedFaces,
globalPointTransformedBoundaryFacesPtr_(),
@ -1581,6 +1586,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
const labelListList& pointSlaves = globalPointSlaves();
const labelListList& pointTransformSlaves =
globalPointTransformedSlaves();
const globalIndexAndTransform& transforms = globalTransforms();
List<labelPairList> transformedCells(pointSlaves.size());
@ -1660,7 +1666,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
{
label proci = globalIndices.whichProcID(slave);
label celli = globalIndices.toLocal(proci, slave);
myBCells[n++] = globalIndexAndTransform::encode
myBCells[n++] = transforms.encode
(
proci,
celli,
@ -1693,7 +1699,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
globalIndices,
globalPointBoundaryCells,
globalTransforms(),
transforms,
transformedCells,
globalPointTransformedBoundaryCellsPtr_(),
@ -1765,12 +1771,12 @@ Foam::globalMeshData::globalMeshData(const polyMesh& mesh)
processorPatchIndices_(0),
processorPatchNeighbours_(0),
nGlobalPoints_(-1),
sharedPointLabelsPtr_(nullptr),
sharedPointAddrPtr_(nullptr),
sharedPointGlobalLabelsPtr_(nullptr),
sharedPointLabelsPtr_(NULL),
sharedPointAddrPtr_(NULL),
sharedPointGlobalLabelsPtr_(NULL),
nGlobalEdges_(-1),
sharedEdgeLabelsPtr_(nullptr),
sharedEdgeAddrPtr_(nullptr)
sharedEdgeLabelsPtr_(NULL),
sharedEdgeAddrPtr_(NULL)
{
updateMesh();
}

View File

@ -72,7 +72,6 @@ See also
mapDistribute
globalIndexAndTransform
SourceFiles
globalMeshData.C
globalMeshDataTemplates.C

View File

@ -64,15 +64,15 @@ Foam::label Foam::globalPoints::findSamePoint
const labelPair& info
) const
{
const label proci = globalIndexAndTransform::processor(info);
const label index = globalIndexAndTransform::index(info);
const label proci = globalTransforms_.processor(info);
const label index = globalTransforms_.index(info);
forAll(allInfo, i)
{
if
(
globalIndexAndTransform::processor(allInfo[i]) == proci
&& globalIndexAndTransform::index(allInfo[i]) == index
globalTransforms_.processor(allInfo[i]) == proci
&& globalTransforms_.index(allInfo[i]) == index
)
{
return i;
@ -98,21 +98,21 @@ Foam::labelPairList Foam::globalPoints::addSendTransform
forAll(info, i)
{
//Pout<< " adding send transform to" << nl
// << " proc:" << globalIndexAndTransform::processor(info[i])
// << " proc:" << globalTransforms_.processor(info[i])
// << nl
// << " index:" << globalIndexAndTransform::index(info[i]) << nl
// << " index:" << globalTransforms_.index(info[i]) << nl
// << " trafo:"
// << globalTransforms_.decodeTransformIndex
// (globalIndexAndTransform::transformIndex(info[i]))
// (globalTransforms_.transformIndex(info[i]))
// << endl;
sendInfo[i] = globalIndexAndTransform::encode
sendInfo[i] = globalTransforms_.encode
(
globalIndexAndTransform::processor(info[i]),
globalIndexAndTransform::index(info[i]),
globalTransforms_.processor(info[i]),
globalTransforms_.index(info[i]),
globalTransforms_.addToTransformIndex
(
globalIndexAndTransform::transformIndex(info[i]),
globalTransforms_.transformIndex(info[i]),
patchi,
true, // patchi is sending side
tol // tolerance for comparison
@ -204,11 +204,11 @@ bool Foam::globalPoints::mergeInfo
}
else
{
label myTransform = globalIndexAndTransform::transformIndex
label myTransform = globalTransforms_.transformIndex
(
myInfo[index]
);
label nbrTransform = globalIndexAndTransform::transformIndex
label nbrTransform = globalTransforms_.transformIndex
(
nbrInfo[i]
);
@ -294,7 +294,7 @@ bool Foam::globalPoints::mergeInfo
labelPairList knownInfo
(
1,
globalIndexAndTransform::encode
globalTransforms_.encode
(
Pstream::myProcNo(),
localPointi,
@ -356,9 +356,9 @@ void Foam::globalPoints::printProcPoint
const labelPair& pointInfo
) const
{
label proci = globalIndexAndTransform::processor(pointInfo);
label index = globalIndexAndTransform::index(pointInfo);
label trafoI = globalIndexAndTransform::transformIndex(pointInfo);
label proci = globalTransforms_.processor(pointInfo);
label index = globalTransforms_.index(pointInfo);
label trafoI = globalTransforms_.transformIndex(pointInfo);
Pout<< " proc:" << proci;
Pout<< " localpoint:";
@ -421,7 +421,7 @@ void Foam::globalPoints::initOwnPoints
labelPairList knownInfo
(
1,
globalIndexAndTransform::encode
globalTransforms_.encode
(
Pstream::myProcNo(),
localPointi,
@ -457,7 +457,7 @@ void Foam::globalPoints::initOwnPoints
labelPairList knownInfo
(
1,
globalIndexAndTransform::encode
globalTransforms_.encode
(
Pstream::myProcNo(),
localPointi,
@ -750,8 +750,8 @@ void Foam::globalPoints::remove
// is in it. This would be an ordinary connection and can be
// handled by normal face-face connectivity.
label proc0 = globalIndexAndTransform::processor(pointInfo[0]);
label proc1 = globalIndexAndTransform::processor(pointInfo[1]);
label proc0 = globalTransforms_.processor(pointInfo[0]);
label proc1 = globalTransforms_.processor(pointInfo[1]);
if
(
@ -759,14 +759,14 @@ void Foam::globalPoints::remove
proc0 == Pstream::myProcNo()
&& directNeighbours.found
(
globalIndexAndTransform::index(pointInfo[0])
globalTransforms_.index(pointInfo[0])
)
)
|| (
proc1 == Pstream::myProcNo()
&& directNeighbours.found
(
globalIndexAndTransform::index(pointInfo[1])
globalTransforms_.index(pointInfo[1])
)
)
)
@ -776,14 +776,14 @@ void Foam::globalPoints::remove
{
//Pout<< "Removing direct neighbour:"
// << mesh_.points()
// [globalIndexAndTransform::index(pointInfo[0])]
// [globalTransforms_.index(pointInfo[0])]
// << endl;
}
else if (proc1 == Pstream::myProcNo())
{
//Pout<< "Removing direct neighbour:"
// << mesh_.points()
// [globalIndexAndTransform::index(pointInfo[1])]
// [globalTransforms_.index(pointInfo[1])]
// << endl;
}
}
@ -812,11 +812,11 @@ void Foam::globalPoints::remove
// So this meshPoint will have info of size one only.
if
(
globalIndexAndTransform::processor(pointInfo[0])
globalTransforms_.processor(pointInfo[0])
!= Pstream::myProcNo()
|| !directNeighbours.found
(
globalIndexAndTransform::index(pointInfo[0])
globalTransforms_.index(pointInfo[0])
)
)
{
@ -995,7 +995,7 @@ void Foam::globalPoints::calculateSharedPoints
forAllConstIter(Map<label>, meshToProcPoint_, iter)
{
labelPairList& pointInfo = procPoints_[iter()];
sort(pointInfo, globalIndexAndTransform::less());
sort(pointInfo, globalIndexAndTransform::less(globalTransforms_));
}
@ -1017,10 +1017,10 @@ void Foam::globalPoints::calculateSharedPoints
if
(
(
globalIndexAndTransform::processor(masterInfo)
globalTransforms_.processor(masterInfo)
== Pstream::myProcNo()
)
&& (globalIndexAndTransform::index(masterInfo) == iter.key())
&& (globalTransforms_.index(masterInfo) == iter.key())
)
{
labelList& pPoints = pointPoints_[iter.key()];
@ -1035,9 +1035,9 @@ void Foam::globalPoints::calculateSharedPoints
for (label i = 1; i < pointInfo.size(); i++)
{
const labelPair& info = pointInfo[i];
label proci = globalIndexAndTransform::processor(info);
label index = globalIndexAndTransform::index(info);
label transform = globalIndexAndTransform::transformIndex
label proci = globalTransforms_.processor(info);
label index = globalTransforms_.index(info);
label transform = globalTransforms_.transformIndex
(
info
);

View File

@ -270,10 +270,10 @@ Foam::mapDistribute::mapDistribute
forAll(transformedElements, i)
{
labelPair elem = transformedElements[i];
label proci = globalIndexAndTransform::processor(elem);
label proci = globalTransforms.processor(elem);
if (proci != Pstream::myProcNo())
{
label index = globalIndexAndTransform::index(elem);
label index = globalTransforms.index(elem);
label nCompact = compactMap[proci].size();
compactMap[proci].insert(index, nCompact);
}
@ -301,7 +301,7 @@ Foam::mapDistribute::mapDistribute
forAll(transformedElements, i)
{
labelPair elem = transformedElements[i];
label trafoI = globalIndexAndTransform::transformIndex(elem);
label trafoI = globalTransforms.transformIndex(elem);
nPerTransform[trafoI]++;
}
// Offset per transformIndex
@ -321,9 +321,9 @@ Foam::mapDistribute::mapDistribute
forAll(transformedElements, i)
{
labelPair elem = transformedElements[i];
label proci = globalIndexAndTransform::processor(elem);
label index = globalIndexAndTransform::index(elem);
label trafoI = globalIndexAndTransform::transformIndex(elem);
label proci = globalTransforms.processor(elem);
label index = globalTransforms.index(elem);
label trafoI = globalTransforms.transformIndex(elem);
// Get compact index for untransformed element
label rawElemI =
@ -378,10 +378,10 @@ Foam::mapDistribute::mapDistribute
forAll(elems, i)
{
label proci = globalIndexAndTransform::processor(elems[i]);
label proci = globalTransforms.processor(elems[i]);
if (proci != Pstream::myProcNo())
{
label index = globalIndexAndTransform::index(elems[i]);
label index = globalTransforms.index(elems[i]);
label nCompact = compactMap[proci].size();
compactMap[proci].insert(index, nCompact);
}
@ -413,7 +413,7 @@ Foam::mapDistribute::mapDistribute
forAll(elems, i)
{
label trafoI = globalIndexAndTransform::transformIndex(elems[i]);
label trafoI = globalTransforms.transformIndex(elems[i]);
nPerTransform[trafoI]++;
}
}
@ -438,9 +438,9 @@ Foam::mapDistribute::mapDistribute
forAll(elems, i)
{
label proci = globalIndexAndTransform::processor(elems[i]);
label index = globalIndexAndTransform::index(elems[i]);
label trafoI = globalIndexAndTransform::transformIndex(elems[i]);
label proci = globalTransforms.processor(elems[i]);
label index = globalTransforms.index(elems[i]);
label trafoI = globalTransforms.transformIndex(elems[i]);
// Get compact index for untransformed element
label rawElemI =

View File

@ -461,7 +461,7 @@ void Foam::mapDistributeBase::distribute
{
// Set up sends to neighbours
List<List<T>> sendFields(Pstream::nProcs());
List<List<T > > sendFields(Pstream::nProcs());
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
@ -495,7 +495,7 @@ void Foam::mapDistributeBase::distribute
// Set up receives from neighbours
List<List<T>> recvFields(Pstream::nProcs());
List<List<T > > recvFields(Pstream::nProcs());
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
@ -938,7 +938,7 @@ void Foam::mapDistributeBase::distribute
{
// Set up sends to neighbours
List<List<T>> sendFields(Pstream::nProcs());
List<List<T > > sendFields(Pstream::nProcs());
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
@ -972,7 +972,7 @@ void Foam::mapDistributeBase::distribute
// Set up receives from neighbours
List<List<T>> recvFields(Pstream::nProcs());
List<List<T > > recvFields(Pstream::nProcs());
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{

View File

@ -25,13 +25,14 @@ License
#include "globalIndexAndTransform.H"
#include "cyclicPolyPatch.H"
#include "DynamicField.H"
#include "globalMeshData.H"
// * * * * * * * * * * * * Private Static Data Members * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(globalIndexAndTransform, 0);
const label globalIndexAndTransform::base_ = 32;
defineTypeNameAndDebug(globalIndexAndTransform, 0);
}
@ -127,10 +128,8 @@ void Foam::globalIndexAndTransform::determineTransforms()
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
transforms_ = List<vectorTensorTransform>(6);
scalarField maxTol(6);
label nextTrans = 0;
DynamicList<vectorTensorTransform> localTransforms;
DynamicField<scalar> localTols;
label dummyMatch = -1;
@ -170,7 +169,7 @@ void Foam::globalIndexAndTransform::determineTransforms()
(
matchTransform
(
transforms_,
localTransforms,
dummyMatch,
transform,
cpp.matchTolerance(),
@ -178,15 +177,8 @@ void Foam::globalIndexAndTransform::determineTransforms()
) == 0
)
{
if (nextTrans == 6)
{
FatalErrorInFunction
<< "More than six unsigned transforms"
<< " detected:" << nl << transforms_
<< exit(FatalError);
}
transforms_[nextTrans] = transform;
maxTol[nextTrans++] = cpp.matchTolerance();
localTransforms.append(transform);
localTols.append(cpp.matchTolerance());
}
}
}
@ -207,7 +199,7 @@ void Foam::globalIndexAndTransform::determineTransforms()
(
matchTransform
(
transforms_,
localTransforms,
dummyMatch,
transform,
cpp.matchTolerance(),
@ -215,15 +207,8 @@ void Foam::globalIndexAndTransform::determineTransforms()
) == 0
)
{
if (nextTrans == 6)
{
FatalErrorInFunction
<< "More than six unsigned transforms"
<< " detected:" << nl << transforms_
<< exit(FatalError);
}
transforms_[nextTrans] = transform;
maxTol[nextTrans++] = cpp.matchTolerance();
localTransforms.append(transform);
localTols.append(cpp.matchTolerance());
}
}
}
@ -233,21 +218,18 @@ void Foam::globalIndexAndTransform::determineTransforms()
// Collect transforms on master
List<List<vectorTensorTransform>> allTransforms(Pstream::nProcs());
allTransforms[Pstream::myProcNo()] = transforms_;
allTransforms[Pstream::myProcNo()] = localTransforms;
Pstream::gatherList(allTransforms);
// Collect matching tolerance on master
List<scalarField> allTols(Pstream::nProcs());
allTols[Pstream::myProcNo()] = maxTol;
allTols[Pstream::myProcNo()] = localTols;
Pstream::gatherList(allTols);
if (Pstream::master())
{
transforms_ = List<vectorTensorTransform>(3);
label nextTrans = 0;
localTransforms.clear();
forAll(allTransforms, proci)
{
@ -264,45 +246,23 @@ void Foam::globalIndexAndTransform::determineTransforms()
(
matchTransform
(
transforms_,
localTransforms,
dummyMatch,
transform,
allTols[proci][pSVI],
true
) == 0
) == 0
)
{
transforms_[nextTrans++] = transform;
}
if (nextTrans > 3)
{
FatalErrorInFunction
<< "More than three independent basic "
<< "transforms detected:" << nl
<< allTransforms
<< transforms_
<< exit(FatalError);
localTransforms.append(transform);
}
}
}
}
transforms_.setSize(nextTrans);
}
transforms_.transfer(localTransforms);
Pstream::scatter(transforms_);
if (transforms_.size() > 3)
{
WarningInFunction
<< "More than three independent basic "
<< "transforms detected:" << nl
<< transforms_ << nl
<< "This is not a space filling tiling and will probably"
<< " give problems for e.g. lagrangian tracking or interpolation"
<< endl;
}
}
@ -351,16 +311,12 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
patchTransformSign_.setSize(patches.size(), Pair<label>(-1, 0));
label matchTransI = -1;
patchTransformSign_.setSize(patches.size(), labelPair(-1, 0));
forAll(patches, patchi)
{
const polyPatch& pp = patches[patchi];
// Pout<< nl << patchi << " " << pp.name() << endl;
// Note: special check for unordered cyclics. These are in fact
// transform bcs and should probably be split off.
if
@ -375,15 +331,12 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
)
)
{
const coupledPolyPatch& cpp =
refCast<const coupledPolyPatch>(pp);
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
if (cpp.separated())
{
const vectorField& sepVecs = cpp.separation();
// Pout<< "sepVecs " << sepVecs << endl;
// This loop is implicitly expecting only a single
// value for separation()
forAll(sepVecs, sVI)
@ -394,6 +347,7 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
{
vectorTensorTransform t(sepVec);
label matchTransI;
label sign = matchTransform
(
transforms_,
@ -402,22 +356,8 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
cpp.matchTolerance(),
true
);
// Pout<< sign << " " << matchTransI << endl;
// List<label> permutation(transforms_.size(), 0);
// permutation[matchTransI] = sign;
// Pout<< encodeTransformIndex(permutation) << nl
// << transformPermutations_
// [
// encodeTransformIndex(permutation)
// ]
// << endl;
patchTransformSign_[patchi] =
Pair<label>(matchTransI, sign);
labelPair(matchTransI, sign);
}
}
@ -426,8 +366,6 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
{
const tensorField& transTensors = cpp.reverseT();
// Pout<< "transTensors " << transTensors << endl;
// This loop is implicitly expecting only a single
// value for reverseT()
forAll(transTensors, tTI)
@ -438,6 +376,7 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
{
vectorTensorTransform t(transT);
label matchTransI;
label sign = matchTransform
(
transforms_,
@ -447,37 +386,65 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
true
);
// Pout<< sign << " " << matchTransI << endl;
// List<label> permutation(transforms_.size(), 0);
// permutation[matchTransI] = sign;
// Pout<< encodeTransformIndex(permutation) << nl
// << transformPermutations_
// [
// encodeTransformIndex(permutation)
// ]
// << endl;
patchTransformSign_[patchi] =
Pair<label>(matchTransI, sign);
labelPair(matchTransI, sign);
}
}
}
}
}
}
// Pout<< patchTransformSign_ << endl;
bool Foam::globalIndexAndTransform::uniqueTransform
(
const point& pt,
labelPairList& trafos,
const label patchi,
const labelPair& patchTrafo
) const
{
if (findIndex(trafos, patchTrafo) == -1)
{
// New transform. Check if already have 3
if (trafos.size() == 3)
{
if (patchi > -1)
{
WarningInFunction
<< "Point " << pt
<< " is on patch " << mesh_.boundaryMesh()[patchi].name();
}
else
{
WarningInFunction
<< "Point " << pt << " is on a coupled patch";
}
Warning
<< " with transformation " << patchTrafo
<< " but also on 3 other patches with transforms "
<< trafos << nl
<< "This is not a space filling tiling and might"
<< " indicate a setup problem and give problems"
<< " for e.g. lagrangian tracking or interpolation" << endl;
// Already warned so no need to extend more
trafos.clear();
return false;
}
return true;
}
else
{
return false;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::globalIndexAndTransform::globalIndexAndTransform
(
const polyMesh& mesh
)
Foam::globalIndexAndTransform::globalIndexAndTransform(const polyMesh& mesh)
:
mesh_(mesh),
transforms_(),
@ -546,13 +513,102 @@ Foam::globalIndexAndTransform::globalIndexAndTransform
Info<< "nullTransformIndex:" << nullTransformIndex() << endl
<< endl;
}
if (transforms_.size() > 0)
{
// Check that the transforms are space filling : any point
// can only use up to three transforms
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// 1. Collect transform&sign per point and do local check
List<labelPairList> pointToTrafos(mesh_.nPoints());
forAll(patches, patchi)
{
const polyPatch& pp = patches[patchi];
const labelPair& transSign = patchTransformSign_[patchi];
if (transSign.first() > -1)
{
const labelList& mp = pp.meshPoints();
forAll(mp, i)
{
labelPairList& trafos = pointToTrafos[mp[i]];
bool newTransform = uniqueTransform
(
mesh_.points()[mp[i]],
trafos,
patchi,
transSign
);
if (newTransform)
{
trafos.append(transSign);
}
}
}
}
// Synchronise across collocated (= untransformed) points
// TBD: there is a big problem in that globalMeshData uses
// globalIndexAndTransform. Triggers recursion.
if (false)
{
const globalMeshData& gmd = mesh_.globalData();
const indirectPrimitivePatch& cpp = gmd.coupledPatch();
const labelList& meshPoints = cpp.meshPoints();
const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
const labelListList& slaves = gmd.globalCoPointSlaves();
List<labelPairList> elems(slavesMap.constructSize());
forAll(meshPoints, i)
{
elems[i] = pointToTrafos[meshPoints[i]];
}
// Pull slave data onto master. No need to update transformed slots.
slavesMap.distribute(elems, false);
// Combine master data with slave data
forAll(slaves, i)
{
labelPairList& trafos = elems[i];
const labelList& slavePoints = slaves[i];
// Combine master with untransformed slave data
forAll(slavePoints, j)
{
const labelPairList& slaveTrafos = elems[slavePoints[j]];
forAll(slaveTrafos, slaveI)
{
bool newTransform = uniqueTransform
(
mesh_.points()[meshPoints[i]],
trafos,
-1,
slaveTrafos[slaveI]
);
if (newTransform)
{
trafos.append(slaveTrafos[slaveI]);
}
}
}
}
}
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::globalIndexAndTransform::~globalIndexAndTransform()
{}
// ************************************************************************* //

View File

@ -28,7 +28,9 @@ Description
Determination and storage of the possible independent transforms
introduced by coupledPolyPatches, as well as all of the possible
permutations of these transforms generated by the presence of
multiple coupledPolyPatches, i.e. more than one cyclic boundary.
multiple coupledPolyPatches, i.e. more than one cyclic boundary. Note that
any given point can be on maximum 3 transforms only (and these transforms
have to be perpendicular)
Also provides global index encoding and decoding for entity
(i.e. cell) index, processor index and transform index (0 or
@ -71,8 +73,15 @@ public:
// - transform
class less
{
const globalIndexAndTransform& gi_;
public:
less(const globalIndexAndTransform& gi)
:
gi_(gi)
{}
inline bool operator()(const labelPair&, const labelPair&) const;
};
@ -87,10 +96,8 @@ private:
//- The possible independent (non-permuted) transforms of the
// geometry, i.e. for a geometry with two cyclics, this
// stores the two transforms, not the eight permutations.
// There may not be more than three transforms in the range
// of coupledPolyPatch geometries (separated XOR
// non-parallel) and symmetries (cuboid periodicity only)
// supported.
// Any point can not be on more than three transforms but overall
// the mesh can have more than three.
List<vectorTensorTransform> transforms_;
//- The permutations of the transforms, stored for lookup
@ -105,13 +112,7 @@ private:
//- Mapping from patch index to which transform it matches (or
// -1 for none) (.first()) and what sign to use for it,
// i.e. +/- 1 (.second()).
List<Pair<label>> patchTransformSign_;
// Private static data
//- Number of spaces to reserve for transform encoding
static const label base_;
List<labelPair> patchTransformSign_;
// Private Member Functions
@ -124,7 +125,7 @@ private:
void determineTransformPermutations();
//- Determine which patch uses which transform (if any) and which
//- Sign to use
// sign to use
void determinePatchTransformSign();
//- Test a list of reference transforms to see if the test
@ -139,16 +140,14 @@ private:
bool checkBothSigns
) const;
//- Encode transform index. Hardcoded to 3 independent transforms max.
inline label encodeTransformIndex
//- Return true if transform is not yet present in trafos. Issues
// warning if too many transforms
bool uniqueTransform
(
const FixedList<Foam::label, 3>& permutationIndices
) const;
//- Decode transform index. Hardcoded to 3 independent transforms max.
inline FixedList<label, 3> decodeTransformIndex
(
const label transformIndex
const point& pt,
labelPairList& trafos,
const label patchi,
const labelPair& patchTrafo
) const;
//- Disallow default bitwise copy construct
@ -160,10 +159,6 @@ private:
public:
//- Declare friendship with the entry class for IO
friend class globalPoints;
// Declare name of the class and its debug switch
ClassName("globalIndexAndTransform");
@ -174,10 +169,6 @@ public:
globalIndexAndTransform(const polyMesh& mesh);
//- Destructor
~globalIndexAndTransform();
// Member Functions
//- Generate a transform index from the permutation indices of
@ -185,9 +176,12 @@ public:
// only be -1, 0 or +1.
inline label encodeTransformIndex
(
const List<label>& permutationIndices
const labelList& permutationIndices
) const;
//- Decode transform index
inline labelList decodeTransformIndex(const label transformIndex) const;
//- Add patch transformation to transformIndex. Return new
// transformIndex. (by default the patch is the sending, not the
// receiving, patch)
@ -221,74 +215,71 @@ public:
) const;
//- Encode index and bare index as components on own processor
inline static labelPair encode
inline labelPair encode
(
const label index,
const label transformIndex
);
) const;
//- Encode index and bare index as components on given processor
inline static labelPair encode
inline labelPair encode
(
const label proci,
const label index,
const label transformIndex
);
) const;
//- Index carried by the object
inline static label index(const labelPair& globalIAndTransform);
inline label index(const labelPair& globalIAndTransform) const;
//- Which processor does this come from?
inline static label processor(const labelPair& globalIAndTransform);
inline label processor(const labelPair& globalIAndTransform) const;
//- Transform carried by the object
inline static label transformIndex
inline label transformIndex(const labelPair& globalIAndTransform) const;
//- Return the number of independent transforms
inline label nIndependentTransforms() const;
//- Return access to the stored independent transforms
inline const List<vectorTensorTransform>& transforms() const;
//- Return access to the permuted transforms
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<labelPair>& patchTransformSign() const;
//- Access the overall (permuted) transform corresponding
// to the transformIndex
inline const vectorTensorTransform& transform
(
const labelPair& globalIAndTransform
);
label transformIndex
) const;
// Access
//- Access the all of the indices of the transform
// permutations corresponding the transforms of the
// listed patch indices. This only allows a maximum of three
// transformations (since routine is used to transform points and
// any given point can only be on 3 or less transforms)
inline labelList transformIndicesForPatches
(
const labelHashSet& patchIs
) const;
//- Return the number of independent transforms
inline label nIndependentTransforms() const;
//- Return access to the stored independent transforms
inline const List<vectorTensorTransform>& transforms() const;
//- Return access to the permuted transforms
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;
//- Access the overall (permuted) transform corresponding
// to the transformIndex
inline const vectorTensorTransform& transform
(
label transformIndex
) const;
//- Access the all of the indices of the transform
// permutations corresponding the transforms of the
// listed patch indices
inline labelList transformIndicesForPatches
(
const labelHashSet& patchIs
) const;
//- Apply all of the transform permutations
// corresponding the transforms of the listed patch
// indices to the supplied point
inline pointField transformPatches
(
const labelHashSet& patchIs,
const point& pt
) const;
//- Apply all of the transform permutations
// corresponding the transforms of the listed patch
// indices to the supplied point
inline pointField transformPatches
(
const labelHashSet& patchIs,
const point& pt
) const;
};

View File

@ -33,8 +33,8 @@ bool Foam::globalIndexAndTransform::less::operator()
const labelPair& b
) const
{
label procA = globalIndexAndTransform::processor(a);
label procB = globalIndexAndTransform::processor(b);
label procA = gi_.processor(a);
label procB = gi_.processor(b);
if (procA < procB)
{
@ -47,8 +47,8 @@ bool Foam::globalIndexAndTransform::less::operator()
else
{
// Equal proc.
label indexA = globalIndexAndTransform::index(a);
label indexB = globalIndexAndTransform::index(b);
label indexA = gi_.index(a);
label indexB = gi_.index(b);
if (indexA < indexB)
{
@ -61,8 +61,8 @@ bool Foam::globalIndexAndTransform::less::operator()
else
{
// Equal index
label transformA = globalIndexAndTransform::transformIndex(a);
label transformB = globalIndexAndTransform::transformIndex(b);
label transformA = gi_.transformIndex(a);
label transformB = gi_.transformIndex(b);
return transformA < transformB;
}
@ -72,7 +72,7 @@ bool Foam::globalIndexAndTransform::less::operator()
Foam::label Foam::globalIndexAndTransform::encodeTransformIndex
(
const List<label>& permutationIndices
const labelList& permutationIndices
) const
{
if (permutationIndices.size() != transforms_.size())
@ -106,68 +106,20 @@ Foam::label Foam::globalIndexAndTransform::encodeTransformIndex
}
Foam::label Foam::globalIndexAndTransform::encodeTransformIndex
(
const FixedList<Foam::label, 3>& permutation
) const
{
if (nIndependentTransforms() == 0)
{
return 0;
}
if (nIndependentTransforms() == 1)
{
return permutation[0]+1;
}
else if (nIndependentTransforms() == 2)
{
return (permutation[1]+1)*3 + (permutation[0]+1);
}
else
{
return
(permutation[2]+1)*9
+ (permutation[1]+1)*3
+ (permutation[0]+1);
}
}
Foam::FixedList<Foam::label, 3>
Foam::globalIndexAndTransform::decodeTransformIndex
Foam::labelList Foam::globalIndexAndTransform::decodeTransformIndex
(
const label transformIndex
) const
{
FixedList<label, 3> permutation(label(0));
labelList permutation(transforms_.size(), 0);
label t = transformIndex;
if (nIndependentTransforms() > 0)
forAll(permutation, i)
{
permutation[0] = (t%3)-1;
if (nIndependentTransforms() > 1)
{
t /= 3;
permutation[1] = (t%3)-1;
if (nIndependentTransforms() > 2)
{
t /= 3;
permutation[2] = (t%3)-1;
}
}
permutation[i] = (t%3)-1;
t /= 3;
}
#ifdef FULLDEBUG
t /= 3;
if (t != 0)
{
FatalErrorInFunction
<< "transformIndex : " << transformIndex
<< " has more than 3 fields."
<< abort(FatalError);
}
#endif
return permutation;
}
@ -180,15 +132,28 @@ Foam::label Foam::globalIndexAndTransform::addToTransformIndex
const scalar tol
) const
{
const Pair<label>& transSign = patchTransformSign_[patchi];
const labelPair& transSign = patchTransformSign_[patchi];
label matchTransI = transSign.first();
// Hardcoded for max 3 transforms only!
if (matchTransI > -1 && matchTransI < 3)
if (matchTransI >= transforms_.size())
{
FixedList<label, 3> permutation = decodeTransformIndex(transformIndex);
FatalErrorInFunction
<< "patch:" << mesh_.boundaryMesh()[patchi].name()
<< " transform:" << matchTransI
<< " out of possible transforms:" << transforms_
<< exit(FatalError);
return labelMin;
}
else if (matchTransI == -1)
{
// No additional transformation for this patch
return transformIndex;
}
else
{
// Decode current set of transforms
labelList permutation(decodeTransformIndex(transformIndex));
// Add patch transform
@ -267,10 +232,6 @@ Foam::label Foam::globalIndexAndTransform::addToTransformIndex
return encodeTransformIndex(permutation);
}
else
{
return transformIndex;
}
}
@ -287,7 +248,7 @@ Foam::label Foam::globalIndexAndTransform::minimumTransformIndex
// Count number of transforms
FixedList<label, 3> permutation0 = decodeTransformIndex(transformIndex0);
labelList permutation0(decodeTransformIndex(transformIndex0));
label n0 = 0;
forAll(permutation0, i)
{
@ -297,7 +258,7 @@ Foam::label Foam::globalIndexAndTransform::minimumTransformIndex
}
}
FixedList<label, 3> permutation1 = decodeTransformIndex(transformIndex1);
labelList permutation1(decodeTransformIndex(transformIndex1));
label n1 = 0;
forAll(permutation1, i)
{
@ -324,8 +285,8 @@ Foam::label Foam::globalIndexAndTransform::subtractTransformIndex
const label transformIndex1
) const
{
FixedList<label, 3> permutation0 = decodeTransformIndex(transformIndex0);
FixedList<label, 3> permutation1 = decodeTransformIndex(transformIndex1);
labelList permutation0(decodeTransformIndex(transformIndex0));
labelList permutation1(decodeTransformIndex(transformIndex1));
forAll(permutation0, i)
{
@ -340,7 +301,7 @@ Foam::labelPair Foam::globalIndexAndTransform::encode
(
const label index,
const label transformIndex
)
) const
{
return encode(Pstream::myProcNo(), index, transformIndex);
}
@ -351,21 +312,22 @@ Foam::labelPair Foam::globalIndexAndTransform::encode
const label proci,
const label index,
const label transformIndex
)
) const
{
if (transformIndex < 0 || transformIndex >= base_)
if (transformIndex < 0 || transformIndex >= transformPermutations_.size())
{
FatalErrorInFunction
<< "TransformIndex " << transformIndex
<< " is outside allowed range of 0 to "
<< base_ - 1
<< transformPermutations_.size() - 1
<< abort(FatalError);
}
if (proci > labelMax/base_)
if (proci > labelMax/transformPermutations_.size())
{
FatalErrorInFunction
<< "Overflow : encoding processor " << proci << " in base " << base_
<< "Overflow : encoding processor " << proci
<< " in base " << transformPermutations_.size()
<< " exceeds capability of label (" << labelMax
<< "). Please recompile with larger datatype for label."
<< exit(FatalError);
@ -374,7 +336,7 @@ Foam::labelPair Foam::globalIndexAndTransform::encode
return labelPair
(
index,
transformIndex + proci*base_
transformIndex + proci*transformPermutations_.size()
);
}
@ -382,7 +344,7 @@ Foam::labelPair Foam::globalIndexAndTransform::encode
Foam::label Foam::globalIndexAndTransform::index
(
const labelPair& globalIAndTransform
)
) const
{
return globalIAndTransform.first();
}
@ -391,18 +353,18 @@ Foam::label Foam::globalIndexAndTransform::index
Foam::label Foam::globalIndexAndTransform::processor
(
const labelPair& globalIAndTransform
)
) const
{
return globalIAndTransform.second()/base_;
return globalIAndTransform.second()/transformPermutations_.size();
}
Foam::label Foam::globalIndexAndTransform::transformIndex
(
const labelPair& globalIAndTransform
)
) const
{
return globalIAndTransform.second() % base_;
return globalIAndTransform.second()%transformPermutations_.size();
}
@ -432,7 +394,7 @@ Foam::label Foam::globalIndexAndTransform::nullTransformIndex() const
}
const Foam::List<Foam::Pair<Foam::label>>&
const Foam::labelPairList&
Foam::globalIndexAndTransform::patchTransformSign() const
{
return patchTransformSign_;
@ -453,7 +415,7 @@ Foam::labelList Foam::globalIndexAndTransform::transformIndicesForPatches
const labelHashSet& patchis
) const
{
List<label> permutation(transforms_.size(), 0);
labelList permutation(transforms_.size(), 0);
labelList selectedTransformIs(0);
@ -466,7 +428,7 @@ Foam::labelList Foam::globalIndexAndTransform::transformIndicesForPatches
{
label patchi = iter.key();
const Pair<label>& transSign = patchTransformSign_[patchi];
const labelPair& transSign = patchTransformSign_[patchi];
label matchTransI = transSign.first();
@ -520,7 +482,7 @@ Foam::labelList Foam::globalIndexAndTransform::transformIndicesForPatches
}
case 2:
{
List<label> tempPermutation = permutation;
labelList tempPermutation = permutation;
label a = 0;
label b = 1;
@ -565,7 +527,7 @@ Foam::labelList Foam::globalIndexAndTransform::transformIndicesForPatches
}
case 3:
{
List<label> tempPermutation = permutation;
labelList tempPermutation = permutation;
tempPermutation[0] = 0;
tempPermutation[1] = 0;