exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Bugfixes and improvements in DLB and AMR. Vuko Vukcevic

Browse source
This commit is contained in:
Hrvoje Jasak 2019-04-09 12:00:24 +01:00
commit 71a25fc3c3
53 changed files with 1892 additions and 830 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/conjugateHeatTransfer/fvPatchFields/chtRcThermalDiffusivity/chtRcThermalDiffusivityFvPatchScalarField.C
View file

@ -261,7 +261,7 @@ Foam::chtRcThermalDiffusivityFvPatchScalarField::calcThermalDiffusivity
const scalarField kOwn = fOwn/(1.0 - p.weights())/mld.magDelta(p.index());
const scalarField kNei = fNei/p.weights()/mld.magDelta(p.index());
tmp<scalarField> kTmp(new scalarField(p.size()));
scalarField& k = kTmp();

794
src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/decomposeMesh.C
View file

File diff suppressed because it is too large Load diff

39
src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/distributeCells.C
View file

@ -184,6 +184,7 @@ void Foam::domainDecomposition::distributeCells()
{
const fvBoundaryMesh& patches = mesh_.boundary();
// Send cellToProc data
forAll (patches, patchI)
{
if (isA<processorFvPatch>(patches[patchI]))
@ -212,7 +213,43 @@ void Foam::domainDecomposition::distributeCells()
cpPatch.internalFieldTransfer
(
Pstream::blocking,
cellToProc_
cellToProc_ // Dummy argument
);
}
}
// Send face cells for correct ordering of faces in parallel load
// balancing when certain processor faces become internal faces on a
// given processor
forAll (patches, patchI)
{
if (isA<processorFvPatch>(patches[patchI]))
// if (patches[patchI].coupled())
{
const lduInterface& cpPatch =
refCast<const lduInterface>(patches[patchI]);
cpPatch.initTransfer
(
Pstream::blocking,
patches[patchI].faceCells()
);
}
}
forAll (patches, patchI)
{
if (isA<processorFvPatch>(patches[patchI]))
// if (patches[patchI].coupled())
{
const lduInterface& cpPatch =
refCast<const lduInterface>(patches[patchI]);
patchNbrFaceCells_[patchI] =
cpPatch.transfer
(
Pstream::blocking,
patches[patchI].faceCells() // Dummy argument
);
}
}

35
src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C
View file

@ -34,6 +34,7 @@ License
#include "OSspecific.H"
#include "Map.H"
#include "DynamicList.H"
#include "globalMeshData.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -53,8 +54,10 @@ Foam::domainDecomposition::domainDecomposition
nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
distributed_(false),
gfIndex_(mesh_),
gpIndex_(mesh_),
cellToProc_(mesh_.nCells()),
patchNbrCellToProc_(mesh_.boundaryMesh().size()),
patchNbrFaceCells_(mesh_.boundaryMesh().size()),
procPointAddressing_(nProcs_),
procFaceAddressing_(nProcs_),
nInternalProcFaces_(nProcs_),
@ -122,8 +125,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
forAll (curFaceLabels, faceI)
{
// Mark the original face as used
// Remember to decrement the index by one (turning index)
// HJ, 5/Dec/2001
// Remember to decrement the index by one (turning index) (see
// procFaceAddressing_ data member comment). HJ, 5/Dec/2001
label curF = mag(curFaceLabels[faceI]) - 1;
faceLookup[curF] = faceI;
@ -174,8 +177,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
// Note: loop over owner, not all faces: sizes are different
forAll (procOwner, faceI)
{
// Remember to decrement the index by one (turning index)
// HJ, 28/Mar/2009
// Remember to decrement the index by one (turning index) (see
// procFaceAddressing_ data member comment). HJ, 5/Dec/2001
label curF = mag(curFaceLabels[faceI]) - 1;
if (curFaceLabels[faceI] >= 0)
@ -193,8 +196,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
// Note: loop over neighbour, not all faces: sizes are different
forAll (procNeighbour, faceI)
{
// Remember to decrement the index by one (turning index)
// HJ, 28/Mar/2009
// Remember to decrement the index by one (turning index) (see
// procFaceAddressing_ data member comment). HJ, 5/Dec/2001
label curF = mag(curFaceLabels[faceI]) - 1;
if (curFaceLabels[faceI] >= 0)
@ -285,6 +288,8 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
forAll (patchGlobalIndex, fI)
{
// Remember to decrement the index by one (turning index) (see
// procFaceAddressing_ data member comment). HJ, 5/Dec/2001
patchGlobalIndex[fI] =
globalIndex[mag(curFaceLabels[curPatchStart + fI]) - 1];
}
@ -503,6 +508,24 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
}
Foam::labelList Foam::domainDecomposition::globalPointIndex
(
const label procI
) const
{
const labelList& gppi = gpIndex_.globalLabel();
const labelList& ppAddr = procPointAddressing_[procI];
labelList globalPointIndex(ppAddr.size());
forAll (globalPointIndex, gpI)
{
globalPointIndex[gpI] = gppi[ppAddr[gpI]];
}
return globalPointIndex;
}
bool Foam::domainDecomposition::writeDecomposition()
{
Info<< "\nConstructing processor meshes" << endl;

34
src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H
View file

@ -43,6 +43,7 @@ SourceFiles
#include "PtrList.H"
#include "point.H"
#include "globalProcFaceIndex.H"
#include "globalProcPointIndex.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -72,6 +73,9 @@ class domainDecomposition
//- Global face index
globalProcFaceIndex gfIndex_;
//- Global face index
globalProcPointIndex gpIndex_;
//- Processor label for each cell
labelList cellToProc_;
@ -79,6 +83,11 @@ class domainDecomposition
// Data is used when running decomposition in parallel
labelListList patchNbrCellToProc_;
//- Face cells for neighbour cell for each processor boundary
// Data is used for correct upper triangular (neighbour list) ordering
// when running decomposition in parallel
labelListList patchNbrFaceCells_;
//- Labels of points for each processor
labelListList procPointAddressing_;
@ -87,9 +96,9 @@ class domainDecomposition
// Only the processor boundary faces are affected: if the sign of the
// index is negative, the processor face is the reverse of the
// original face. In order to do this properly, all face
// indices will be incremented by 1 and the decremented as
// necessary t avoid the problem of face number zero having no
// sign. HJ, 5/Dec/2001
// indices will be incremented by 1 and then decremented as
// necessary to avoid the problem of face number zero having no
// sign. HJ, 5/Dec/2001
labelListList procFaceAddressing_;
//- Number of internal faces for each processor
@ -133,12 +142,28 @@ class domainDecomposition
//- Create cell-to processor list using domain decomposition tools
void distributeCells();
//- Helper function for decomposeMesh. Given two lists of inter processor
// boundaries and inter processor boundary faces, add neighbouring procJ
// to list of the procI. Also adds faceI into the procI list and return
// whether the procJ is already present in the list of visited
// processors (needed for error handling on cyclic patches)
bool addInterProcessorBoundaryData
(
const label& procI,
const label& procJ,
const label faceI,
List<SLList<label> >& interProcBoundaries,
List<SLList<SLList<label> > >& interProcBFaces
) const;
public:
// Declare name of the class and its debug switch
ClassName("domainDecomposition");
// Constructors
//- Construct from mesh and dictionary
@ -205,6 +230,9 @@ public:
const bool createPassiveProcPatches = false
) const;
//- Create and return global point index
labelList globalPointIndex(const label procI) const;
//- Return processor point addressing
const labelListList& procPointAddressing() const
{

8
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/fvFieldReconstructorReconstructFields.C
View file

@ -44,10 +44,10 @@ void Foam::fvFieldReconstructor::reconstructField
) const
{
// Get references to internal and boundary field
Field<Type>& iField = reconField.internalField();
Field<Type>& iField = reconField;
typename GeometricField<Type, fvPatchField, volMesh>::
GeometricBoundaryField& bouField = reconField.boundaryField();
GeometricBoundaryField& bouField = reconField.boundaryFieldNoStoreOldTimes();
forAll (procFields, procI)
{
@ -143,10 +143,10 @@ void Foam::fvFieldReconstructor::reconstructField
) const
{
// Create the internalField
Field<Type>& iField = reconField.internalField();
Field<Type>& iField = reconField;
typename GeometricField<Type, fvsPatchField, surfaceMesh>::
GeometricBoundaryField& bouField = reconField.boundaryField();
GeometricBoundaryField& bouField = reconField.boundaryFieldNoStoreOldTimes();
forAll (procFields, procI)
{

204
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
View file

@ -437,15 +437,13 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
}
}
// Prepare patch reconstruction
HashTable<label, word> patchNameLookup
(
meshes_[firstValidMesh()].boundaryMesh().size()
);
DynamicList<word> patchTypeLookup
(
meshes_[firstValidMesh()].boundaryMesh().size()
);
// Get first boundary mesh size
const label firstBndryMeshSize =
meshes_[firstValidMesh()].boundaryMesh().size();
// Prepare patch reconstruction. Note: default key type is word in HashTable
HashTable<label> patchNameLookup(firstBndryMeshSize);
DynamicList<word> patchTypeLookup(firstBndryMeshSize);
label nReconPatches = 0;
@ -459,35 +457,37 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
{
if (!isA<processorPolyPatch>(procPatches[patchI]))
{
const word& patchName = procPatches[patchI].name();
// Regular patch, get the patch and the name
const polyPatch& patch = procPatches[patchI];
const word& patchName = patch.name();
// Regular patch. Try to find it in a list
if (!patchNameLookup.found(patchName))
// Insert the pair into the hash table
if (patchNameLookup.insert(patchName, nReconPatches))
{
// Patch not found. Add it
patchNameLookup.insert(patchName, nReconPatches);
patchTypeLookup.append(procPatches[patchI].type());
nReconPatches++;
// This is the first time we're inserting this patch,
// add its type into the list and increment the counter
patchTypeLookup.append(patch.type());
++nReconPatches;
}
// else already present in the table
}
}
}
}
// Fill in patch names and types
// Fill in patch names
wordList reconPatchNames(patchNameLookup.size());
wordList reconPatchTypes(patchNameLookup.size());
wordList patchNameToc = patchNameLookup.toc();
forAll (patchNameToc, pnI)
// Loop through all the patch names and collect names and types
forAllConstIter(HashTable<label>, patchNameLookup, iter)
{
const label pnIndex = patchNameLookup.find(patchNameToc[pnI])();
reconPatchNames[pnIndex] = patchNameToc[pnI];
reconPatchTypes[pnIndex] = patchTypeLookup[pnIndex];
reconPatchNames[iter()] = iter.key();
}
// Transfer the contents of the patchTypeLookup dynamic list into the
// ordinary list. Note: the ordering remains the same
const wordList reconPatchTypes(patchTypeLookup.xfer());
// Prepare point, face and patch reconstruction
label nReconPoints = 0;
label nReconFaces = 0;
@ -498,30 +498,31 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
{
if (meshes_.set(procI))
{
// Count total number of points and faces
nReconPoints += meshes_[procI].nPoints();
nReconFaces += meshes_[procI].allFaces().size();
nReconCells += meshes_[procI].nCells();
// Get current mesh
const fvMesh& procMesh = meshes_[procI];
const polyBoundaryMesh& procPatches = meshes_[procI].boundaryMesh();
// Count total number of points and faces
nReconPoints += procMesh.nPoints();
nReconFaces += procMesh.allFaces().size();
nReconCells += procMesh.nCells();
const polyBoundaryMesh& procPatches = procMesh.boundaryMesh();
forAll (procPatches, patchI)
{
if (!isA<processorPolyPatch>(procPatches[patchI]))
{
// Get current patch
const polyPatch& patch = procPatches[patchI];
// Find processor patch index in reconstructed boundary
const label pnIndex = patchNameLookup.find
(
procPatches[patchI].name()
)();
const label pnIndex = patchNameLookup.find(patch.name())();
// Check patch name and type
if
(
procPatches[patchI].name()
!= reconPatchNames[pnIndex]
|| procPatches[patchI].type()
!= reconPatchTypes[pnIndex]
patch.name() != reconPatchNames[pnIndex]
|| patch.type() != reconPatchTypes[pnIndex]
)
{
FatalErrorIn
@ -530,13 +531,12 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
"reconstructMesh(const Time& db)"
) << "Patch name and type does not match "
<< "across processors for patch "
<< procPatches[patchI].name() << " type: "
<< procPatches[patchI].type()
<< patch.name() << " type: " << patch.type()
<< abort(FatalError);
}
// Record number of faces in patch
reconPatchSizes[pnIndex] += procPatches[patchI].size();
reconPatchSizes[pnIndex] += patch.size();
}
}
}
@ -553,10 +553,11 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
forAll (reconPatchFaces, patchI)
{
// Set size of reconstructed patch faces
reconPatchFaces[patchI].setSize(reconPatchSizes[patchI]);
reconPatchOwner[patchI].setSize(reconPatchSizes[patchI]);
reconPatchOwner[patchI] = -1;
// Set size of reconstructed patch face owners with default value of -1
reconPatchOwner[patchI].setSize(reconPatchSizes[patchI], -1);
}
// Size the mapping arrays
@ -652,10 +653,30 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
reconPatchSizes = 0;
// Prepare handling for globally shared points. This is equivalent
// to parallel processor points, but working on a PtrList of meshes
// on the same processor
sharedPoints sharedData(meshes_);
// HR 21.12.18 : A bit of a hack for the moment in order to support the new
// method (using global point addressing) and old method (syncing across
// processor BCs) of constructing the shared points. The old method is still
// needed since global point addressing does not exist when the sharedPoint
// object is constructed from reconstructParMesh.
//
// TODO: Unify the methods by constructing global point addressing from the
// mesh pieces when. Or even better, calculate procPointAddressing directly
// which will simplify the mesh merging immensely.
autoPtr<sharedPoints> sharedDataPtr;
if(globalPointIndex_.size())
{
// Determine globally shared points using global point
// adressing
sharedDataPtr.set(new sharedPoints(meshes_, globalPointIndex_));
}
else
{
// Prepare handling for globally shared points. This is equivalent
// to parallel processor points, but working on a PtrList of meshes
// on the same processor
sharedDataPtr.set(new sharedPoints(meshes_));
}
sharedPoints& sharedData = sharedDataPtr();
// Record global point index for shared points
labelList globalPointMapping(sharedData.nGlobalPoints(), -1);
@ -784,7 +805,11 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
// Dump first valid mesh without checking
{
const label fvmId = firstValidMesh();
Pout<< "Dump mesh " << fvmId << endl;
if (debug)
{
Pout<< "Dump mesh " << fvmId << endl;
}
cellOffset[fvmId] = 0;
@ -874,14 +899,14 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
}
}
// Get sorting order. Note make a copy of indices because
// Get sorting order. Note: make a copy of indices because
// sortable list will be deleted
labelList procVisitOrder =
const labelList procVisitOrder =
SortableList<label>(nrbProcIndex).indices();
forAll (procVisitOrder, pvoI)
{
const label patchI = procVisitOrder[pvoI];
const label& patchI = procVisitOrder[pvoI];
if (isA<processorPolyPatch>(procPatches[patchI]))
{
@ -944,15 +969,14 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
forAll (indices, i)
{
// Location in reconstructed patch where the face
// is inserted
const label insertSlot = indices[i];
// Location in reconstructed patch where the face is
// inserted
const label& insertSlot = indices[i];
// Calculate face index depending on the ordering
const label faceI = patchStart + i;
face newFace = curFaces[faceI];
inplaceRenumber(ppAddr, newFace);
// Insert into correct slot
@ -1015,7 +1039,9 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
cpAddr[cellI] = cellI; // + cellOffset[firstValidMesh()];
}
// Set cell offset for the next mesh
// Set cell offset for the next mesh. Note: if the next mesh is not
// valid, the cell offset is propagated to the next one in the processor
// loop below.
if (cellOffset.size() > firstValidMesh() + 1)
{
cellOffset[firstValidMesh() + 1] =
@ -1025,14 +1051,25 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
// Dump all other meshes, merging the processor boundaries
for (label procI = firstValidMesh() + 1; procI < meshes_.size(); procI++)
{
if (meshes_.set(procI))
if (!meshes_.set(procI))
{
Pout<< "Dump mesh " << procI
<< " cell offset: " << cellOffset[procI]
<< endl;
// Next mesh is not valid: simply propagate cell offset
if (cellOffset.size() > procI + 1)
{
cellOffset[procI + 1] = cellOffset[procI];
}
}
else
{
// Valid mesh, combine it
if (debug)
{
Pout<< "Dump mesh " << procI
<< " cell offset: " << cellOffset[procI]
<< endl;
}
const polyMesh& curMesh = meshes_[procI];
const polyBoundaryMesh& procPatches = curMesh.boundaryMesh();
@ -1170,7 +1207,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
}
} // End of "is neighbour"
} // End of "is processor"
}
} // End for all patches
// Dump unmarked points into the global point list
label nMergedPoints = 0;
@ -1274,7 +1311,7 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
// Get sorting order. Note make a copy of indices because
// sortable list will be deleted
labelList procVisitOrder =
const labelList procVisitOrder =
SortableList<label>(nrbProcIndex).indices();
forAll (procVisitOrder, pvoI)
@ -1347,13 +1384,14 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
faceI++
)
{
label faceInPatch = faceI - patchStart;
const label faceInPatch = faceI - patchStart;
// Calculate master index
label masterIndex = masterProcPatch.start()
+ faceInPatch;
const label masterIndex =
masterProcPatch.start() + faceInPatch;
label masterFp = masterFaceAddr[masterIndex] - 1;
const label masterFp =
masterFaceAddr[masterIndex] - 1;
// Record face-cells for the neighbour
fpAddr[faceI] = -masterFaceAddr[masterIndex];
@ -1470,14 +1508,6 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
cellOffset[procI + 1] = cellOffset[procI] + meshes_[procI].nCells();
}
}
else
{
// No valid mesh. Propagate cell offset
if (cellOffset.size() > procI + 1)
{
cellOffset[procI + 1] = cellOffset[procI];
}
}
}
// Resize the lists
@ -1519,12 +1549,24 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
// Mesh assembly completed
Info<< "Global mesh size (final): " << nl
<< " nPoints = " << reconPoints.size() << nl
<< " nFaces = " << reconFaces.size() << nl
<< " nCells = " << nReconCells << nl
<< " nPatches = " << reconPatchSizes.size() << nl
<< " nPatchFaces = " << reconPatchSizes << endl;
if (!Pstream::parRun())
{
Info<< "Global mesh size (final): " << nl
<< " nPoints = " << reconPoints.size() << nl
<< " nFaces = " << reconFaces.size() << nl
<< " nCells = " << nReconCells << nl
<< " nPatches = " << reconPatchSizes.size() << nl
<< " nPatchFaces = " << reconPatchSizes << endl;
}
else if (debug)
{
Pout<< "Local mesh size (final): " << nl
<< " nPoints = " << reconPoints.size() << nl
<< " nFaces = " << reconFaces.size() << nl
<< " nCells = " << nReconCells << nl
<< " nPatches = " << reconPatchSizes.size() << nl
<< " nPatchFaces = " << reconPatchSizes << endl;
}
// Renumber the face-processor addressing list for all pieces
// now that the number of internal faces is known

7
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C
View file

@ -25,6 +25,11 @@ License
#include "processorMeshesReconstructor.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::processorMeshesReconstructor, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::processorMeshesReconstructor::readMeshes(PtrList<Time>& databases)
@ -71,6 +76,7 @@ Foam::processorMeshesReconstructor::processorMeshesReconstructor
:
meshName_(meshName),
meshes_(),
globalPointIndex_(),
pointProcAddressing_(),
faceProcAddressing_(),
cellProcAddressing_(),
@ -86,6 +92,7 @@ Foam::processorMeshesReconstructor::processorMeshesReconstructor
:
meshName_(meshName),
meshes_(databases.size()),
globalPointIndex_(),
pointProcAddressing_(),
faceProcAddressing_(),
cellProcAddressing_(),

15
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.H
View file

@ -29,8 +29,6 @@ Description
matching processor boundaries and build a single combined mesh by
matching processor patches to each other.
In order to
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved.
@ -72,6 +70,9 @@ class processorMeshesReconstructor
//- List of processor meshes
PtrList<fvMesh> meshes_;
//- global point index per sub-mesh
labelListList globalPointIndex_;
//- List of processor point addressing lists
PtrList<labelIOList> pointProcAddressing_;
@ -121,6 +122,10 @@ class processorMeshesReconstructor
public:
//- Declare name of the class and its debug switch
ClassName("processorMeshesReconstructor");
// Constructors
//- Construct given name. Set meshes later
@ -161,6 +166,12 @@ public:
return meshes_;
}
//- Return to global point index
labelListList& globalPointIndex()
{
return globalPointIndex_;
}
//- Return point-processor addressing arrays
const PtrList<labelIOList>& pointProcAddressing() const;

171
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C
View file

@ -451,9 +451,160 @@ void Foam::sharedPoints::calcSharedPoints()
}
void Foam::sharedPoints::calcSharedPoints
(
const labelListList& globalPointIndex
)
{
typedef HashTable<Pair<label>, label, Hash<label> > markTable;
markTable marker;
List<labelHashSet> procSharedPoints(meshes_.size());
// Mark up points for the first time
forAll (meshes_, meshI)
{
if (meshes_.set(meshI))
{
const labelList& curGlobalPointIndex = globalPointIndex[meshI];
labelHashSet& sharedPoints = procSharedPoints[meshI];
forAll (curGlobalPointIndex, pI)
{
const label gpi = curGlobalPointIndex[pI];
markTable::iterator iter = marker.find(gpi);
if (iter == Map<label>::end())
{
// Record meshI and pI
marker.insert(gpi, Pair<label>(-meshI - 1, pI));
}
else
{
if (iter().first() < 0)
{
if(iter().first() != -meshI - 1)
{
// Shared point detected. Add for both meshes
// Add for first mesh
const label firstMesh = -iter().first() - 1;
const label firstPointI = iter().second();
procSharedPoints[firstMesh].insert(firstPointI);
// Add for current mesh
sharedPoints.insert(pI);
// Count shared points and update bookkeeping
iter().first() = nGlobalPoints_;
nGlobalPoints_++;
}
}
else
{
// Existing shared point. Add for current mesh
sharedPoints.insert(pI);
}
}
}
}
}
forAll (meshes_, meshI)
{
if (meshes_.set(meshI))
{
const labelList& curGlobalPointIndex = globalPointIndex[meshI];
labelList& curSharedPoints = sharedPointLabels_[meshI];
curSharedPoints = procSharedPoints[meshI].toc();
labelList& curSharedAddr = sharedPointAddr_[meshI];
curSharedAddr.setSize(curSharedPoints.size());
forAll(curSharedPoints, i)
{
const label gpi = curSharedPoints[i];
curSharedAddr[i] = marker.find(curGlobalPointIndex[gpi])().first();
}
}
}
// debug
{
pointField gp(nGlobalPoints_);
boolList gpSet(nGlobalPoints_, false);
forAll (meshes_, meshI)
{
if (meshes_.set(meshI))
{
// Get points
const pointField& P = meshes_[meshI].points();
// Get list of local point labels that are globally shared
const labelList& curShared = sharedPointLabels_[meshI];
// Get index in global point list
const labelList& curAddr = sharedPointAddr_[meshI];
// Loop through all local points
forAll (curShared, i)
{
if (!gpSet[curAddr[i]])
{
// Point not set: set it
gp[curAddr[i]] = P[curShared[i]];
gpSet[curAddr[i]] = true;
}
else
{
// Point already set: check location
if (mag(gp[curAddr[i]] - P[curShared[i]]) > SMALL)
{
Info<< "MERGE MISMATCH: mesh" << meshI
<< " point: " << curShared[i]
<< " dist: " << gp[curAddr[i]] << " "
<< P[curShared[i]]
<< endl;
}
}
}
}
}
/*
// Grab marked points
OFstream ppp("points.vtk");
ppp << "# vtk DataFile Version 2.0" << nl
<< "points.vtk" << nl
<< "ASCII" << nl
<< "DATASET POLYDATA" << nl
<< "POINTS " << nGlobalPoints_ << " float" << nl;
Pout << "Global marking " << nGlobalPoints_ << endl;
forAll (gp, i)
{
ppp << float(gp[i].x()) << ' '
<< float(gp[i].y()) << ' '
<< float(gp[i].z())
<< nl;
Pout << gp[i] << endl;
}
*/
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sharedPoints::sharedPoints(const PtrList<fvMesh>& meshes)
Foam::sharedPoints::sharedPoints
(
const PtrList<fvMesh>& meshes
)
:
meshes_(meshes),
sharedPointAddr_(meshes_.size()),
@ -463,11 +614,19 @@ Foam::sharedPoints::sharedPoints(const PtrList<fvMesh>& meshes)
calcSharedPoints();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
// Foam::sharedPoints::~sharedPoints()
// {}
Foam::sharedPoints::sharedPoints
(
const PtrList<fvMesh>& meshes,
const labelListList& globalPointIndex
)
:
meshes_(meshes),
sharedPointAddr_(meshes_.size()),
sharedPointLabels_(meshes_.size()),
nGlobalPoints_(0)
{
calcSharedPoints(globalPointIndex);
}
// ************************************************************************* //

15
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H
View file

@ -117,13 +117,26 @@ class sharedPoints
//- Calculate shared points
void calcSharedPoints();
//- Calculate shared points
void calcSharedPoints(const labelListList& globalPointIndex);
public:
// Constructors
//- Construct from the list of meshes
sharedPoints(const PtrList<fvMesh>& meshes);
sharedPoints
(
const PtrList<fvMesh>& meshes
);
//- Construct from the list of meshes and global point addressing
sharedPoints
(
const PtrList<fvMesh>& meshes,
const labelListList& globalPointIndex
);
//- Destructor

30
src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/polyhedralRefinement/polyhedralRefinement.C
View file

@ -31,6 +31,7 @@ Author
#include "faceSet.H"
#include "pointSet.H"
#include "addToRunTimeSelectionTable.H"
#include "meshTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -2042,7 +2043,16 @@ void Foam::polyhedralRefinement::setCellsToRefine
// layers
for (label i = 0; i < nRefinementBufferLayers_; ++i)
{
extendMarkedCellsAcrossFaces(refineCell);
meshTools::extendMarkedCellsAcrossFaces(mesh_, refineCell);
}
// Remove all cells that would exceed the maximum refinement level
forAll (refineCell, cellI)
{
if (refineCell[cellI] && (cellLevel_[cellI] + 1 > maxRefinementLevel_))
{
refineCell[cellI] = false;
}
}
// Remove all cells that would exceed the maximum refinement level
@ -2086,6 +2096,12 @@ void Foam::polyhedralRefinement::setCellsToRefine
++nIters;
nTotalAddCells += nAddCells;
if (debug)
{
Info<< "Added " << nAddCells << " cells in iteration "
<< nIters << nl;
}
} while (nAddCells > 0);
Info<< "Added " << nTotalAddCells // nTotalAddCells already reduced
@ -2205,7 +2221,7 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
// Note: if there is no dynamic load balancing, points at the boundary
// cannot be split points by definition. However, in dynamic load balancing
// runs, it is possible that a split point end on processor boundary, in
// runs, it is possible that a split point ends on processor boundary, in
// which case we will simply avoid (actually delay) unrefining until this
// becomes internal point again. VV, 4/Jul/2018.
const label nInternalFaces = mesh_.nInternalFaces();
@ -2213,7 +2229,7 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
for (label faceI = nInternalFaces; faceI < nFaces; ++faceI)
{
// Get the face and make sure that the points are unarked
// Get the face and make sure that the points are unmarked
const face& f = meshFaces[faceI];
forAll (f, fpI)
@ -2260,7 +2276,7 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
// unrefinement buffer layers
for (label i = 0; i < nUnrefinementBufferLayers_; ++i)
{
extendMarkedCellsAcrossPoints(protectedCell);
meshTools::extendMarkedCellsAcrossPoints(mesh_, protectedCell);
}
// Loop through all cells and if the cell should be protected, protect all
@ -2368,6 +2384,12 @@ void Foam::polyhedralRefinement::setSplitPointsToUnrefine
++nIters;
nTotalRemCells += nRemCells;
if (debug)
{
Info<< "Removed " << nRemCells << " cells in iteration "
<< nIters << nl;
}
} while (nRemCells > 0);
Info<< "Removed " << nTotalRemCells // nTotalRemCells already reduced

21
src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/prismatic2DRefinement/prismatic2DRefinement.C
View file

@ -33,6 +33,7 @@ Author
#include "emptyPolyPatch.H"
#include "wedgePolyPatch.H"
#include "addToRunTimeSelectionTable.H"
#include "meshTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -2601,7 +2602,16 @@ void Foam::prismatic2DRefinement::setCellsToRefine
// layers
for (label i = 0; i < nRefinementBufferLayers_; ++i)
{
extendMarkedCellsAcrossFaces(refineCell);
meshTools::extendMarkedCellsAcrossFaces(mesh_, refineCell);
}
// Remove all cells that exceed the maximum refinement level
forAll (refineCell, cellI)
{
if (refineCell[cellI] && (cellLevel_[cellI] + 1 > maxRefinementLevel_))
{
refineCell[cellI] = false;
}
}
// Remove all cells that exceed the maximum refinement level
@ -2695,7 +2705,7 @@ void Foam::prismatic2DRefinement::setSplitPointsToUnrefine
// 1. Has pointLevel_ > 0 (obviously),
// 2. A point that has the same pointLevel_ as ALL of the points of its
// edges. In other words, for each point, we will look through all the
// edges of the point. For each edges, we will visit both points and
// edges of the point. For each edge, we will visit both points and
// check point levels. All point levels must be the same for this point
// candidate to be a split point. This is quite useful since there is no
// need to store the refinement history
@ -2779,11 +2789,14 @@ void Foam::prismatic2DRefinement::setSplitPointsToUnrefine
if (isA<processorPolyPatch>(patch))
{
// Get patch start
const label startIndex = patch.start();
// Loop through all the faces
forAll (patch, i)
{
// Get global face index and face
const label faceI = patch.size() + i;
const label faceI = startIndex + i;
const face& f = meshFaces[faceI];
// Make sure that we don't split around point at all points of
@ -2835,7 +2848,7 @@ void Foam::prismatic2DRefinement::setSplitPointsToUnrefine
// unrefinement buffer layers
for (label i = 0; i < nUnrefinementBufferLayers_; ++i)
{
extendMarkedCellsAcrossPoints(protectedCell);
meshTools::extendMarkedCellsAcrossPoints(mesh_, protectedCell);
}
// Loop through all cells and if the cell should be protected, protect all

264
src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.C
View file

@ -66,130 +66,6 @@ void Foam::refinement::setInstance(const fileName& inst) const
}
void Foam::refinement::extendMarkedCellsAcrossFaces
(
boolList& markedCell
) const
{
// Mark all faces for all marked cells
const label nFaces = mesh_.nFaces();
boolList markedFace(nFaces, false);
// Get mesh cells
const cellList& meshCells = mesh_.cells();
// Loop through all cells
forAll (markedCell, cellI)
{
if (markedCell[cellI])
{
// This cell is marked, get its faces
const cell& cFaces = meshCells[cellI];
forAll (cFaces, i)
{
markedFace[cFaces[i]] = true;
}
}
}
// Snyc the face list across processor boundaries
syncTools::syncFaceList(mesh_, markedFace, orEqOp<bool>(), false);
// Get necessary mesh data
const label nInternalFaces = mesh_.nInternalFaces();
const labelList& owner = mesh_.faceOwner();
const labelList& neighbour = mesh_.faceNeighbour();
// Internal faces
for (label faceI = 0; faceI < nInternalFaces; ++faceI)
{
if (markedFace[faceI])
{
// Face is marked, mark both owner and neighbour if the maximum
// refinement level is not exceeded
const label& own = owner[faceI];
const label& nei = neighbour[faceI];
// Mark owner and neighbour cells
markedCell[own] = true;
markedCell[nei] = true;
}
}
// Boundary faces
for (label faceI = nInternalFaces; faceI < nFaces; ++faceI)
{
if (markedFace[faceI])
{
// Face is marked, mark owner if the maximum refinement level is not
// exceeded
const label& own = owner[faceI];
// Mark owner
markedCell[own] = true;
}
}
}
void Foam::refinement::extendMarkedCellsAcrossPoints
(
boolList& markedCell
) const
{
// Mark all points for all marked cells
const label nPoints = mesh_.nPoints();
boolList markedPoint(nPoints, false);
// Get cell points
const labelListList& meshCellPoints = mesh_.cellPoints();
// Loop through all cells
forAll (markedCell, cellI)
{
if (markedCell[cellI])
{
// This cell is marked, get its points
const labelList& cPoints = meshCellPoints[cellI];
forAll (cPoints, i)
{
markedPoint[cPoints[i]] = true;
}
}
}
// Snyc point list across processor boundaries
syncTools::syncPointList
(
mesh_,
markedPoint,
orEqOp<bool>(),
true, // Default value
true // Apply separation for parallel cyclics
);
// Get point cells
const labelListList& meshPointCells = mesh_.pointCells();
// Loop through all points
forAll (markedPoint, pointI)
{
if (markedPoint[pointI])
{
// This point is marked, mark all of its cells
const labelList& pCells = meshPointCells[pointI];
forAll (pCells, i)
{
markedCell[pCells[i]] = true;
}
}
}
}
Foam::label Foam::refinement::addFace
(
polyTopoChange& ref,
@ -797,7 +673,14 @@ Foam::label Foam::refinement::faceConsistentRefinement
// Note: we are using more stringent 1:1 consistency across coupled
// boundaries in order to simplify handling of edge based consistency
// checks for parallel runs
if (neiLevel[i] > curOwnLevel)
// Bugfix related to PLB: Check whether owner is already marked for
// refinement. Will allow 2:1 consistency across certain processor faces
// where we have a new processor boundary. VV, 23/Jan/2019.
if
(
(neiLevel[i] > curOwnLevel)
&& !cellsToRefine[own]
)
{
// Neighbour level is higher than owner level, owner must be
// marked for refinement
@ -949,9 +832,11 @@ Foam::label Foam::refinement::faceConsistentUnrefinement
<< "Try increasing nUnrefinementBufferLayers. "
<< abort(FatalError);
}
cellsToUnrefine[own] = false;
++nRemCells;
else
{
cellsToUnrefine[own] = false;
++nRemCells;
}
}
else if (neiLevel < (ownLevel - 1))
{
@ -977,9 +862,11 @@ Foam::label Foam::refinement::faceConsistentUnrefinement
<< "Try increasing nUnrefinementBufferLayers. "
<< abort(FatalError);
}
cellsToUnrefine[nei] = false;
++nRemCells;
else
{
cellsToUnrefine[nei] = false;
++nRemCells;
}
}
}
@ -1011,7 +898,7 @@ Foam::label Foam::refinement::faceConsistentUnrefinement
// Note: we are using more stringent 1:1 consistency across coupled
// boundaries in order to simplify handling of edge based consistency
// checkes for parallel runs
// checks for parallel runs
if (curOwnLevel < neiLevel[i])
{
// Owner level is smaller than neighbour level, we must not
@ -1046,7 +933,11 @@ Foam::label Foam::refinement::faceConsistentUnrefinement
<< "This is probably because the refinement and "
<< "unrefinement regions are very close." << nl
<< "Try increasing nUnrefinementBufferLayers. "
<< abort(FatalError);
<< nl
<< "Another possibility is that you are running "
<< "with Dynamic Load Balancing, in which case "
<< "this should be fine."
<< endl;
}
}
else
@ -1123,26 +1014,31 @@ Foam::label Foam::refinement::edgeConsistentUnrefinement
// unrefinement
if (!cellsToUnrefine[cellI])
{
FatalErrorIn
(
"label refinement::"
"edgeConsistentUnrefinement"
"(boolList& cellsToUnrefine)"
) << "Cell not marked for unrefinement, indicating a"
<< " previous unnoticed problem with unrefinement."
<< nl
<< "cellI: " << cellI << ", cellJ: " << cellJ
<< nl
<< "Level of cellI: " << cellILevel
<< ", level of cellJ: " << cellJLevel << nl
<< "This is probably because the refinement and "
<< "unrefinement regions are very close." << nl
<< "Try increasing nUnrefinementBufferLayers. "
<< abort(FatalError);
if (debug)
{
WarningIn
(
"label refinement::"
"edgeConsistentUnrefinement"
"(boolList& cellsToUnrefine)"
) << "Cell not marked for unrefinement, indicating a"
<< " previous unnoticed problem with unrefinement."
<< nl
<< "cellI: " << cellI << ", cellJ: " << cellJ
<< nl
<< "Level of cellI: " << cellILevel
<< ", level of cellJ: " << cellJLevel << nl
<< "This is probably because the refinement and "
<< "unrefinement regions are very close." << nl
<< "Try increasing nUnrefinementBufferLayers. "
<< endl;
}
}
else
{
cellsToUnrefine[cellI] = false;
++nRemCells;
}
cellsToUnrefine[cellI] = false;
++nRemCells;
}
else if (cellJLevel < cellILevel - 1)
{
@ -1153,26 +1049,31 @@ Foam::label Foam::refinement::edgeConsistentUnrefinement
// unrefinement
if (!cellsToUnrefine[cellJ])
{
FatalErrorIn
(
"label refinement::"
"edgeConsistentUnrefinement"
"(boolList& cellsToUnrefine)"
) << "Cell not marked for unrefinement, indicating a"
<< " previous unnoticed problem with unrefinement."
<< nl
<< "cellI: " << cellI << ", cellJ: " << cellJ
<< nl
<< "Level of cellI: " << cellILevel
<< ", level of cellJ: " << cellJLevel << nl
<< "This is probably because the refinement and "
<< "unrefinement regions are very close." << nl
<< "Try increasing nUnrefinementBufferLayers. "
<< abort(FatalError);
if (debug)
{
WarningIn
(
"label refinement::"
"edgeConsistentUnrefinement"
"(boolList& cellsToUnrefine)"
) << "Cell not marked for unrefinement, indicating a"
<< " previous unnoticed problem with unrefinement."
<< nl
<< "cellI: " << cellI << ", cellJ: " << cellJ
<< nl
<< "Level of cellI: " << cellILevel
<< ", level of cellJ: " << cellJLevel << nl
<< "This is probably because the refinement and "
<< "unrefinement regions are very close." << nl
<< "Try increasing nUnrefinementBufferLayers. "
<< endl;
}
}
else
{
cellsToUnrefine[cellJ] = false;
++nRemCells;
}
cellsToUnrefine[cellJ] = false;
++nRemCells;
}
}
}
@ -1442,7 +1343,7 @@ void Foam::refinement::setRefinement(polyTopoChange& ref) const
mesh_,
pointLevel_,
maxEqOp<label>(),
label(0), // Null value
0, // Null value
true // Apply separation for parallel cyclics
);
@ -1599,10 +1500,21 @@ void Foam::refinement::updateMesh(const mapPolyMesh& map)
}
}
// Note: new point level is going to be synced at processor boundaries just
// before the next step in setRefinement. Need to investigate why the sync
// is not done properly if I put it here. Something is not updated yet.
// VV, 31/Jan/2018.
// Sync the new point level. Note: here, we assume that the call to
// updateMesh happened after polyBoundaryMesh::updateMesh where the
// processor data is fully rebuilt. If this is not the case, the point
// level will remain unsynced and will cause all kinds of trouble that
// are extremely difficult to spot. See the change in
// polyTopoChanger::changeMesh order of calling polyMesh::updateMesh and
// polyTopoChanger::update. VV, 19/Feb/2019.
syncTools::syncPointList
(
mesh_,
newPointLevel,
maxEqOp<label>(),
0, // Null value
true // Apply separation for parallel cyclics
);
// Transfer the new point level into the data member
pointLevel_.transfer(newPointLevel);

10
src/dynamicMesh/dynamicMesh/polyMeshModifiers/refinement/refinement/refinement.H
View file

@ -159,12 +159,6 @@ protected:
//- Set file instance for cellLevel_ and pointLevel_
void setInstance(const fileName& inst) const;
//- Extend marked cells across faces
void extendMarkedCellsAcrossFaces(boolList& markedCell) const;
//- Extend marked cells across points
void extendMarkedCellsAcrossPoints(boolList& markedCell) const;
// Local topology modification functions (operate on cells/faces)
@ -189,8 +183,8 @@ protected:
const label nei
) const;
//- Modifies existing faceI for either new owner/neighbour or new face
// points. Reverses if nessecary
//- Modifies existing faceI for either new owner/neighbour or new
// face points. Reverses if nessecary
void modifyFace
(
polyTopoChange& ref,

11
src/dynamicMesh/dynamicMesh/polyTopoChange/polyTopoChanger/polyTopoChangerChangeMesh.C
View file

@ -2490,10 +2490,17 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChanger::changeMesh()
topoChangeRequest()()
);
update(topoChangeMap());
// Mesh data needs to be updated before the update of polyMeshModifiers
// because polyMeshModifiers might need all the new polyMesh data (see
// below for further comments). VV, 19/Feb/2019
mesh_.updateMesh(topoChangeMap());
// Bugfix: call to polyTopoChanger::update must happen after
// polyMesh::updateMesh where all the relevant mesh bits for parallel
// comms are updated. First noticed when the syncying of pointLevel in
// refinement::updateMesh was not syncying properly. VV, 19/Feb/2019
update(topoChangeMap());
// Increment the morph index
morphIndex_++;

32
src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
View file

@ -271,15 +271,11 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
(
refinementSelectionPtr_->refinementCellCandidates()()
);
if (debug)
{
Pout<< "Selected " << refCandidates.size()
<< " refinement candidates."
<< endl;
}
Pout<< "Selected " << refCandidates.size()
<< " refinement candidates."
<< endl;
}
else if (debug)
else
{
Pout<< "Skipping refinement for this time-step..." << endl;
}
@ -301,15 +297,11 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
(
refinementSelectionPtr_->unrefinementPointCandidates()()
);
if (debug)
{
Pout<< "Selected " << unrefCandidates.size()
<< " unrefinement candidates."
<< endl;
}
Pout<< "Selected " << unrefCandidates.size()
<< " unrefinement candidates."
<< endl;
}
else if (debug)
else
{
Pout<< "Skipping unrefinement for this time-step..." << endl;
}
@ -352,13 +344,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
// some topo changes
if (sizeCellMap)
{
Info<< "Successfully performed polyhedral refinement. "
Pout<< "Successfully performed polyhedral refinement. "
<< "Changed from " << nOldCells << " to " << sizeCellMap
<< " cells." << endl;
}
else
{
Info<< "Refinement/unrefinement not performed in this time step "
Pout<< "Refinement/unrefinement not performed in this time step "
<< "since no cells were selected." << endl;
}
@ -370,9 +362,7 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
// per time step
curTimeIndex_ = time().timeIndex();
}
Info<< "No refinement/unrefinement" << endl;
Pout<< "No refinement/unrefinement" << endl;
// No refinement/unrefinement at this time step. Return false
return false;
}

113
src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
View file

@ -33,6 +33,7 @@ License
#include "addToRunTimeSelectionTable.H"
#include "cloud.H"
#include "cloudDistribute.H"
#include "meshObjectBase.H"
#include "IFstream.H"
#include "OFstream.H"
@ -107,7 +108,11 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
// Now that each processor has filled in its own part, combine the data
Pstream::gatherList(migratedCells);
Pstream::scatterList(migratedCells);
Info<< "Migrated cells per processor: " << migratedCells << endl;
if (debug)
{
Info<< "Migrated cells per processor: " << migratedCells << endl;
}
// Reading through second index now tells how many cells will arrive
// from which processor
@ -141,10 +146,14 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
// Prepare receiving side
// Create the reconstructor
processorMeshesReconstructor meshRecon("reconstructed");
// HR 21.12.18 : Use empty domainname to avoid auto-created of
// fvSchemes/fvSolution
processorMeshesReconstructor meshRecon("");
PtrList<fvMesh>& procMeshes = meshRecon.meshes();
procMeshes.setSize(meshDecomp.nProcs());
meshRecon.globalPointIndex().setSize(meshDecomp.nProcs());
// Collect local fields for decomposition
clearOut();
@ -245,6 +254,29 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
}
}
// HR 14.12.18: Create dummy database pointing into the non-parallel case
// directory and disable the runTimeModifiable switch. dummyTime is used
// for decomposed, received and reconstructed fvMeshes (ie. before its data
// is moved into *this).
//
// The pointing into the non-parallel case directory is somewhat a hack to
// prevent auto-creation of fvSchemes and fvSolution (they exist in the root
// case). This is potentially dangerous if we write anything, but fvMeshes
// derived from this database are NO_WRITE so we should be fine. Making if
// non-runTimeModifiable prevents registration of fvSolution with the
// fileMonitor (addWatch). Not all processors will necessarily receive a
// mesh and watches will then cause dead-locks!
Time dummyTime
(
time().rootPath(),
time().globalCaseName(),
"system",
"constant",
false
);
const_cast<Switch&>(dummyTime.runTimeModifiable()) = false;
for (label procI = 0; procI < meshDecomp.nProcs(); procI++)
{
// Check if there is a mesh to send
@ -256,8 +288,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
autoPtr<fvMesh> procMeshPtr = meshDecomp.processorMesh
(
procI,
time(),
"processorPart" + Foam::name(procI),
dummyTime,
"", // HR 21.12.18 : Use empty domainname to avoid
// auto-created offvSchemes/fvSolution
true // Create passive processor patches
);
fvMesh& procMesh = procMeshPtr();
@ -289,6 +322,7 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
// Send the mesh and fields to target processor
toProc << procMesh << nl;
toProc << meshDecomp.globalPointIndex(procI) << nl;
// Send fields
sendFields(volScalarFields, fieldDecomposer, toProc);
@ -345,6 +379,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
procMeshPtr
);
meshRecon.globalPointIndex()[procI] =
meshDecomp.globalPointIndex(procI);
// Set local fields
// Note: first index is field index and second index is procI
insertFields
@ -413,7 +450,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
);
}
//HJ Insert clouds missing. HJ, 12/Oct/2018
// HR, 18.11.2018. Distribution of clouds is trivial and is
// treated in Cloud<ParticleType>::split, which is called in
// the constructor of CloudDistribute.
}
}
}
@ -450,9 +489,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
(
IOobject
(
"processorPart" + Foam::name(procI),
time().timeName(),
time(),
"", // "processorPart" + Foam::name(procI),
dummyTime.timeName(),
dummyTime,
IOobject::NO_READ,
IOobject::NO_WRITE
),
@ -461,6 +500,10 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
)
);
// Receive the global points addr
labelList gppi(fromProc);
meshRecon.globalPointIndex()[procI] = gppi;
// Receive fields
// Note: first index is field index and second index is procI
receiveFields
@ -524,34 +567,40 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
}
}
forAll (procMeshes, procI)
if (debug)
{
if (procMeshes.set(procI))
forAll (procMeshes, procI)
{
Pout<< "procMesh " << procI
<< " points " << procMeshes[procI].nPoints()
<< " faces: " << procMeshes[procI].nFaces()
<< " internal: " << procMeshes[procI].nInternalFaces()
<< " cells: " << procMeshes[procI].nCells()
<< " patches: " << procMeshes[procI].boundary().size()
<< endl;
if (procMeshes.set(procI))
{
Pout<< "procMesh " << procI
<< " points " << procMeshes[procI].nPoints()
<< " faces: " << procMeshes[procI].nFaces()
<< " internal: " << procMeshes[procI].nInternalFaces()
<< " cells: " << procMeshes[procI].nCells()
<< " patches: " << procMeshes[procI].boundary().size()
<< endl;
}
}
}
// Create the reconstructed mesh
autoPtr<fvMesh> reconstructedMeshPtr =
meshRecon.reconstructMesh(time());
meshRecon.reconstructMesh(dummyTime);
fvMesh& reconMesh = reconstructedMeshPtr();
Pout<< "Reconstructed mesh stats: "
<< " nCells: " << reconMesh.nCells()
<< " nFaces: " << reconMesh.nFaces()
<< " nIntFaces: " << reconMesh.nInternalFaces()
<< " polyPatches: "
<< reconMesh.boundaryMesh().size()
<< " patches: "
<< reconMesh.boundary().size()
<< endl;
if (debug)
{
Pout<< "Reconstructed mesh stats: "
<< " nCells: " << reconMesh.nCells()
<< " nFaces: " << reconMesh.nFaces()
<< " nIntFaces: " << reconMesh.nInternalFaces()
<< " polyPatches: "
<< reconMesh.boundaryMesh().size()
<< " patches: "
<< reconMesh.boundary().size()
<< endl;
}
// Apply changes to the local mesh:
// - refactor the boundary to match new patches. Note: processor
@ -920,8 +969,14 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
);
}
// Debug: remove? HJ, 22/Oct/2018
// checkMesh(true);
// HR 13.12.18: Update the mesh objects
meshObjectBase::allUpdateTopology<polyMesh>(*this, meshMap);
if (debug)
{
Info<< "Checking reconstructed mesh after load balancing..." << endl;
checkMesh(true);
}
return true;
}

77
src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C
View file

@ -89,7 +89,6 @@ void Foam::topoChangerFvMesh::insertFields
iter
)
{
localFields[fI].set
(
Pstream::myProcNo(),
@ -126,7 +125,7 @@ void Foam::topoChangerFvMesh::receiveFields
receivedFields.setSize(nScalarFields);
fromProc.readBegin("topoChangerFvMeshReceiveFields");
forAll (receivedFields, fI)
{
word fieldName(fromProc);
@ -181,7 +180,10 @@ void Foam::topoChangerFvMesh::rebuildFields
const PtrList<GeoField>& partFields = receivedFields[fieldI];
Pout<< "Rebuilding field " << masterField.name() << endl;
if (debug)
{
Pout<< "Rebuilding field " << masterField.name() << endl;
}
// Check name match. Note: there may be holes
word partName;
@ -263,15 +265,19 @@ void Foam::topoChangerFvMesh::rebuildFields
// Resize internal field
if
(
masterField.internalField().size()
masterField.size()
!= GeoField::GeoMeshType::size(masterField.mesh())
)
{
Pout<< "Resizing internal field: old size = "
<< masterField.internalField().size()
<< " new size = "
<< GeoField::GeoMeshType::size(masterField.mesh())
<< endl;
if (debug)
{
Pout<< "Resizing internal field: old size = "
<< masterField.size()
<< " new size = "
<< GeoField::GeoMeshType::size(masterField.mesh())
<< endl;
}
masterField.setSize
(
GeoField::GeoMeshType::size(masterField.mesh())
@ -280,13 +286,16 @@ void Foam::topoChangerFvMesh::rebuildFields
// Resize boundary (number of patches)
typename GeoField::GeometricBoundaryField& patchFields =
masterField.boundaryField();
masterField.boundaryFieldNoStoreOldTimes();
if (patchFields.size() != masterField.mesh().boundary().size())
{
Pout<< "Resizing boundary field: "
if (debug)
{
Pout<< "Resizing boundary field: "
<< masterField.mesh().boundary().size()
<< endl;
}
patchFields.setSize(masterField.mesh().boundary().size());
}
@ -297,9 +306,12 @@ void Foam::topoChangerFvMesh::rebuildFields
if (meshMap.resetPatchFlag()[patchI])
{
// Create a new constrained patch field
Pout<< "Inserting constrained patch field for patch "
<< masterField.mesh().boundary()[patchI].name()
<< endl;
if (debug)
{
Pout<< "Inserting constrained patch field for patch "
<< masterField.mesh().boundary()[patchI].name()
<< endl;
}
patchFields.set
(
@ -326,12 +338,15 @@ void Foam::topoChangerFvMesh::rebuildFields
)
{
// Resize patch field
Pout<< "Resizing patch field for patch "
<< masterField.mesh().boundary()[patchI].name()
<< " old size: " << patchFields[patchI].size()
<< " new size: "
<< masterField.mesh().boundary()[patchI].size()
<< endl;
if (debug)
{
Pout<< "Resizing patch field for patch "
<< masterField.mesh().boundary()[patchI].name()
<< " old size: " << patchFields[patchI].size()
<< " new size: "
<< masterField.mesh().boundary()[patchI].size()
<< endl;
}
// Reset patch field size
patchFields[patchI].autoMap
@ -350,7 +365,27 @@ void Foam::topoChangerFvMesh::rebuildFields
// Increment field counter
fieldI++;
Pout<< "... done" << endl;
if (debug)
{
Pout<< "... done" << endl;
}
}
// HR 14.12.18: We create new processor boundary faces from internal
// faces. The values on these faces could be initialised by interpolation.
// Instead we choose to fix the values by evaluating the boundaries.
// I tried to execute evaluateCoupled() at the end of
// fvFieldReconstructor::reconstructField, but this fails in a strange way.
forAllConstIter
(
typename HashTable<const GeoField*>,
geoFields,
iter
)
{
GeoField& masterField = const_cast<GeoField&>(*iter());
masterField.boundaryField().evaluateCoupled();
}
}

21
src/finiteVolume/cfdTools/compressible/compressibleCourantNo.H
View file

@ -33,7 +33,12 @@ scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
scalar velMag = 0.0;
if (mesh.nInternalFaces())
// HR 26.06.18: A parallel run has at least two cells and therefore at least
// one internal face in the global mesh. It may be a processor boundary, but
// this is captured by max(mag(phi)).
// Old formulation hangs on parallel cases where one partition is degenerated
// to a single cell.
if (mesh.nInternalFaces() || Pstream::parRun())
{
surfaceScalarField phiOverRho = mag(phi)/fvc::interpolate(rho);
@ -47,6 +52,20 @@ if (mesh.nInternalFaces())
velMag = max(phiOverRho/mesh.magSf()).value();
}
else
{
// Single cell mesh: Co is still defined; use cell formulation
const scalar deltaT = runTime.deltaT().value();
const scalar deltaX = Foam::cbrt(mesh.V()[0]);
velMag = mag(U[0]);
CoNum = velMag*deltaT/deltaX;
meanCoNum = CoNum;
}
Info<< "Courant Number mean: " << meanCoNum
<< " max: " << CoNum

2
src/finiteVolume/cfdTools/incompressible/CourantNo.H
View file

@ -33,7 +33,7 @@ scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
scalar velMag = 0.0;
// HR 26.06.28: A parallel run has at least two cells and therefore at least
// HR 26.06.18: A parallel run has at least two cells and therefore at least
// one internal face in the global mesh. It may be a processor boundary, but
// this is captured by max(mag(phi)).
// Old formulation hangs on parallel cases where one partition is degenerated

8
src/finiteVolume/fields/fvPatchFields/basic/fixedGradient/fixedGradientFvPatchField.C
View file

@ -67,11 +67,11 @@ fixedGradientFvPatchField<Type>::fixedGradientFvPatchField
const dictionary& dict
)
:
fvPatchField<Type>(p, iF, dict),
// HR 15.12.18: Must not call evaluate during construction. Read value
// instead. This is needed for PLB.
fvPatchField<Type>(p, iF, dict, true),
gradient_("gradient", dict, p.size())
{
evaluate();
}
{}
template<class Type>

8
src/finiteVolume/fields/fvPatchFields/basic/mixed/mixedFvPatchField.C
View file

@ -54,13 +54,13 @@ mixedFvPatchField<Type>::mixedFvPatchField
const dictionary& dict
)
:
fvPatchField<Type>(p, iF, dict),
// HR 15.12.18: Must not call evaluate during construction. Read value
// instead. This is needed for PLB.
fvPatchField<Type>(p, iF, dict, true),
refValue_("refValue", dict, p.size()),
refGrad_("refGradient", dict, p.size()),
valueFraction_("valueFraction", dict, p.size())
{
evaluate();
}
{}
template<class Type>

2
src/finiteVolume/fields/fvPatchFields/derived/inletOutlet/inletOutletFvPatchField.C
View file

@ -62,7 +62,7 @@ inletOutletFvPatchField<Type>::inletOutletFvPatchField
{
// Read patch type
this->readPatchType(dict);
this->refValue() = Field<Type>("inletValue", dict, p.size());
if (dict.found("value"))

29
src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
View file

@ -70,7 +70,12 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
"fvSchemes",
obr.time().system(),
obr,
IOobject::READ_IF_PRESENT_IF_MODIFIED, // Allow default dictionary creation
(
obr.readOpt() == IOobject::MUST_READ
|| obr.readOpt() == IOobject::READ_IF_PRESENT
? IOobject::MUST_READ_IF_MODIFIED
: obr.readOpt()
),
IOobject::NO_WRITE
)
),
@ -168,21 +173,17 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
fluxRequired_(), // Do not read flux required option
defaultFluxRequired_(false)
{
if (!headerOk())
// HR 21.12.18 : Writing a default fvSchemes is not useful in PLB. Please
// specify MUST_READ on obr if you need this.
if
(
readOpt() == IOobject::MUST_READ
|| readOpt() == IOobject::MUST_READ_IF_MODIFIED
|| (readOpt() == IOobject::READ_IF_PRESENT && headerOk())
)
{
if (debug)
{
InfoIn
(
"fvSchemes::fvSchemes(const objectRegistry& obr)"
) << "fvSchemes dictionary not found. Creating default."
<< endl;
}
regIOobject::write();
read();
}
read();
}

3
src/finiteVolume/fvMesh/fvMesh.C
View file

@ -424,6 +424,9 @@ void Foam::fvMesh::resetFvPrimitives
// Reset fvPatches HJ, 16/Apr/2018
boundary_.resetFvPatches(resetFvPatchFlag);
// HR 14.12.18: Indicate that the mesh is changing to e.g. nearWallDist
polyMesh::changing(true);
// Clear all mesh data
clearOut();
}

2
src/finiteVolume/fvMesh/wallDist/nearWallDist.C
View file

@ -84,7 +84,7 @@ void Foam::nearWallDist::correct()
{
setSize(mesh_.boundary().size());
}
// Update size of GeometricBoundaryField
forAll (mesh_.boundary(), patchI)
{

1
src/foam/Make/files
View file

@ -447,6 +447,7 @@ $(globalMeshData)/globalMeshData.C
$(globalMeshData)/globalPoints.C
$(globalMeshData)/globalIndex.C
$(globalMeshData)/globalProcFaceIndex.C
$(globalMeshData)/globalProcPointIndex.C
$(polyMesh)/syncTools/syncTools.C

10
src/foam/fields/GeometricFields/GeometricField/GeometricField.C
View file

@ -728,6 +728,16 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryField()
}
// Return reference to GeometricBoundaryField
template<class Type, template<class> class PatchField, class GeoMesh>
typename
Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField&
Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryFieldNoStoreOldTimes()
{
return boundaryField_;
}
// Store old-time field
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTimes() const

4
src/foam/fields/GeometricFields/GeometricField/GeometricField.H
View file

@ -448,6 +448,10 @@ public:
//- Return reference to GeometricBoundaryField for const field
inline const GeometricBoundaryField& boundaryField() const;
//- Return reference to GeometricBoundaryField without storing old
//- times
GeometricBoundaryField& boundaryFieldNoStoreOldTimes();
//- Return the time index of the field
inline label timeIndex() const;

7
src/foam/fields/cloud/cloud.C
View file

@ -64,6 +64,13 @@ Foam::autoPtr<Foam::cloudDistribute> Foam::cloud::cloudDist
}
Foam::labelList Foam::cloud::nParticlesPerCell() const
{
NotImplemented;
return labelList(0);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::cloud::~cloud()

4
src/foam/fields/cloud/cloud.H
View file

@ -101,7 +101,9 @@ public:
//- Return size of the cloud
virtual label size() const;
//- Count and return number of particles per cell
virtual labelList nParticlesPerCell() const;
// Edit
//- Remap the cells of particles corresponding to the

125
src/foam/meshes/meshTools/meshTools.C
View file

@ -27,6 +27,7 @@ License
#include "polyMesh.H"
#include "hexMatcher.H"
#include "faceZone.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -927,4 +928,128 @@ void Foam::meshTools::setFaceInfo
}
void Foam::meshTools::extendMarkedCellsAcrossFaces
(
const polyMesh& mesh,
boolList& markedCell
)
{
// Mark all faces for all marked cells
const label nFaces = mesh.nFaces();
boolList markedFace(nFaces, false);
// Get mesh cells
const cellList& meshCells = mesh.cells();
// Loop through all cells
forAll (markedCell, cellI)
{
if (markedCell[cellI])
{
// This cell is marked, get its faces
const cell& cFaces = meshCells[cellI];
forAll (cFaces, i)
{
markedFace[cFaces[i]] = true;
}
}
}
// Snyc the face list across processor boundaries
syncTools::syncFaceList(mesh, markedFace, orEqOp<bool>(), false);
// Get necessary mesh data
const label nInternalFaces = mesh.nInternalFaces();
const labelList& owner = mesh.faceOwner();
const labelList& neighbour = mesh.faceNeighbour();
// Internal faces
for (label faceI = 0; faceI < nInternalFaces; ++faceI)
{
if (markedFace[faceI])
{
// Face is marked, mark both owner and neighbour
const label& own = owner[faceI];
const label& nei = neighbour[faceI];
// Mark owner and neighbour cells
markedCell[own] = true;
markedCell[nei] = true;
}
}
// Boundary faces
for (label faceI = nInternalFaces; faceI < nFaces; ++faceI)
{
if (markedFace[faceI])
{
// Face is marked, mark owner
const label& own = owner[faceI];
// Mark owner
markedCell[own] = true;
}
}
}
void Foam::meshTools::extendMarkedCellsAcrossPoints
(
const polyMesh& mesh,
boolList& markedCell
)
{
// Mark all points for all marked cells
const label nPoints = mesh.nPoints();
boolList markedPoint(nPoints, false);
// Get cell points
const labelListList& meshCellPoints = mesh.cellPoints();
// Loop through all cells
forAll (markedCell, cellI)
{
if (markedCell[cellI])
{
// This cell is marked, get its points
const labelList& cPoints = meshCellPoints[cellI];
forAll (cPoints, i)
{
markedPoint[cPoints[i]] = true;
}
}
}
// Snyc point list across processor boundaries
syncTools::syncPointList
(
mesh,
markedPoint,
orEqOp<bool>(),
true, // Default value
true // Apply separation for parallel cyclics
);
// Get point cells
const labelListList& meshPointCells = mesh.pointCells();
// Loop through all points
forAll (markedPoint, pointI)
{
if (markedPoint[pointI])
{
// This point is marked, mark all of its cells
const labelList& pCells = meshPointCells[pointI];
forAll (pCells, i)
{
markedCell[pCells[i]] = true;
}
}
}
}
// ************************************************************************* //

19
src/foam/meshes/meshTools/meshTools.H
View file

@ -318,6 +318,25 @@ namespace meshTools
label& zoneFlip
);
// Mark-up of mesh bits. Relocated from refinement polyMeshModifier
//- Extend marked cells across faces given a bool list of already marked
// cells
void extendMarkedCellsAcrossFaces
(
const polyMesh& mesh,
boolList& markedCell
);
//- Extend marked cells across points given a bool list of already
// marked cells
void extendMarkedCellsAcrossPoints
(
const polyMesh& mesh,
boolList& markedCell
);
} // End namespace meshTools

32
src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
View file

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "pointBoundaryMesh.H"
#include "polyMesh.H"
#include "polyBoundaryMesh.H"
#include "facePointPatch.H"
#include "globalPointPatch.H"
@ -105,10 +106,39 @@ void Foam::pointBoundaryMesh::movePoints()
}
void Foam::pointBoundaryMesh::updateMesh()
void Foam::pointBoundaryMesh::updateMesh
(
const polyMesh& pMesh
)
{
const polyBoundaryMesh& pBoundary = pMesh.boundaryMesh();
pointPatchList& patches = *this;
// 21.1.19 HR : Patches need to be recreated in PLB
patches.clear();
patches.setSize(pBoundary.size());
forAll(patches, patchI)
{
patches.set
(
patchI,
facePointPatch::New(pBoundary[patchI], *this).ptr()
);
}
// Add the globalPointPatch
if(pMesh.globalData().nGlobalPoints())
{
patches.setSize(pBoundary.size() + 1);
patches.set
(
patches.size() - 1,
new globalPointPatch(*this, patches.size() - 1)
);
}
forAll(patches, patchi)
{
patches[patchi].initUpdateMesh();

6
src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
View file

@ -46,6 +46,7 @@ namespace Foam
class pointMesh;
class polyBoundaryMesh;
class globalPointPatch;
class polyMesh;
/*---------------------------------------------------------------------------*\
Class pointBoundaryMesh Declaration
@ -105,7 +106,10 @@ public:
void movePoints();
//- Correct polyBoundaryMesh after topology update
void updateMesh();
void updateMesh
(
const polyMesh& pMesh
);
};

5
src/foam/meshes/pointMesh/pointMesh.C
View file

@ -122,7 +122,10 @@ bool Foam::pointMesh::updateMesh(const mapPolyMesh& mpm) const
{
// Casting const-ness to answer the interface of meshObject
// HJ, 30/Aug/2010
const_cast<pointBoundaryMesh&>(boundary_).updateMesh();
const_cast<pointBoundaryMesh&>(boundary_).updateMesh
(
GeoMesh<polyMesh>::mesh_
);
// Map all registered point fields
mapFields(mpm);

8
src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.C
View file

@ -35,7 +35,7 @@ void Foam::globalProcFaceIndex::calcFaceIndex()
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Assing unique face label to all master processor faces
// Assign unique face label to all master processor faces
// Count faces and processor faces per processor
forAll (patches, patchI)
@ -118,7 +118,8 @@ void Foam::globalProcFaceIndex::calcFaceIndex()
OPstream toProc
(
Pstream::nonBlocking,
// HR 12.12.18: nonBlocking fails on PLB of Aachen bomb
Pstream::blocking,
procPatch.neighbProcNo()
);
toProc<< curFaceLabels;
@ -154,7 +155,8 @@ void Foam::globalProcFaceIndex::calcFaceIndex()
// Receive the data from master and insert into the list
IPstream fromProc
(
Pstream::nonBlocking,
// HR 12.12.18: nonBlocking fails on PLB of Aachen bomb
Pstream::blocking,
procPatch.neighbProcNo()
);

6
src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.H
View file

@ -33,7 +33,7 @@ Description
Face offsets counts a number of unique faces on each processor,
excluding slave processor patch faces, which are given master face index.
If needed, global face indes from all faces can be derived from this data
If needed, global face index from all faces can be derived from this data
Currently, faces are ordered with internal faces first, followed by patch
faces in patch order, excluding slave processor patches.
@ -120,9 +120,9 @@ public:
return procFaceOffset_;
}
//- Return lobal face label for all faces of current mesh
//- Return global face label for all faces of current mesh
// Sized to number of live faces in the mesh
const labelList globalLabel() const
const labelList& globalLabel() const
{
return globalLabel_;
}

319
src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.C Normal file
View file

@ -0,0 +1,319 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "globalProcPointIndex.H"
#include "processorPolyPatch.H"
#include "globalMeshData.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::globalProcPointIndex::calcPointIndex()
{
// Count and mark points in the following order
// 1. global points
// 2. processor patches
// 3. regular patches
// 4. anything else ("internal points")
//
// The marking is a follows
// -4 : Not visited yet
// -3 : Found in a regular patch
// -2 : Found in slave processor patch
// -1 : Found in a master processor patch
// 0 - nGlobalPoints-1 : Global point ID
// 1. Mark the global points
const labelList& spl = mesh_.globalData().sharedPointLabels();
const labelList& spa = mesh_.globalData().sharedPointAddr();
forAll(spl, spI)
{
globalLabel_[spl[spI]] = spa[spI];
}
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
labelList procPointCounter(Pstream::nProcs(), 0);
labelList patchPointCounter(patches.size(), 0);
label& nProcPoints = procPointCounter[Pstream::myProcNo()];
// Master processor patches - first pass
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& patchMeshPoints = pp.meshPoints();
if (isA<processorPolyPatch>(pp))
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
label& nPatchPoints = patchPointCounter[patchI];
if (procPatch.master())
{
forAll (patchMeshPoints, mpI)
{
label& pl = globalLabel_[patchMeshPoints[mpI]];
if (pl < 0)
{
nProcPoints++;
nPatchPoints++;
pl = -1;
}
}
}
}
}
// Slave processor patches - first pass
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& patchMeshPoints = pp.meshPoints();
if (isA<processorPolyPatch>(pp))
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
if (!procPatch.master())
{
forAll (patchMeshPoints, mpI)
{
label& pl = globalLabel_[patchMeshPoints[mpI]];
if (pl == -4)
{
pl = -2;
}
}
}
}
}
// Regular patches - first pass
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& patchMeshPoints = pp.meshPoints();
if (!isA<processorPolyPatch>(pp))
{
label& nPatchPoints = patchPointCounter[patchI];
forAll (patchMeshPoints, mpI)
{
label& pl = globalLabel_[patchMeshPoints[mpI]];
if (pl == -4)
{
nProcPoints++;
nPatchPoints++;
pl = -3;
}
}
}
}
// "Internal" points - first pass
label nInternalPoints = 0;
forAll (globalLabel_, pI)
{
if (globalLabel_[pI] == -4)
{
nInternalPoints++;
nProcPoints++;
}
}
// Gather-Scatter counter data
Pstream::gatherList(procPointCounter);
Pstream::scatterList(procPointCounter);
// Convert counters to offsets
procPointOffset_[0] = mesh_.globalData().nGlobalPoints();
for (label procI = 1; procI < procPointOffset_.size(); procI++)
{
procPointOffset_[procI] =
procPointCounter[procI-1] + procPointOffset_[procI-1];
}
patchPointOffset_[0] =
procPointOffset_[Pstream::myProcNo()] + nInternalPoints;
for (label patchI = 1; patchI < patchPointOffset_.size(); patchI++)
{
patchPointOffset_[patchI] =
patchPointCounter[patchI-1] + patchPointOffset_[patchI-1];
}
const pointField& points = mesh_.points();
// Master processor patches - second pass
// Send patch offset to slave and assign global labels to points on
// master processor patches and regular patches
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& patchMeshPoints = pp.meshPoints();
const label patchPO = patchPointOffset_[patchI];
if (isA<processorPolyPatch>(pp))
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
if (procPatch.master())
{
// Master processor patch
// Send the offset to the slave and mark the points through a
// face loop to establish an order that can be untangled on
// slave side
OPstream toProc
(
Pstream::nonBlocking,
procPatch.neighbProcNo()
);
toProc<< patchPO;
pointField pointLocs(patchMeshPoints.size());
label glPointLabelsI = 0;
forAll (pp, faceI)
{
const face& curFace = pp[faceI];
forAll (curFace, fpI)
{
const label pointI = curFace[fpI];
label& pl = globalLabel_[pointI];
if (pl == -1)
{
pl = patchPO + glPointLabelsI;
pointLocs[glPointLabelsI] = points[pointI];
glPointLabelsI++;
}
}
}
}
}
}
// Slave processor and regular patches - second pass
// Receive data on slave and assign global labels to points on
// master processor patches and regular patches
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (isA<processorPolyPatch>(pp))
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
if (!procPatch.master())
{
// Slave processor patch - second pass
// Receive the offset from master and mark the points by taking
// into account that the points in the faces are in reverse
// order
IPstream fromProc
(
Pstream::nonBlocking,
procPatch.neighbProcNo()
);
label masterPatchPO(readLabel(fromProc));
label glPointLabelsI = 0;
forAll (pp, faceI)
{
const face curFace = pp[faceI].reverseFace();
forAll (curFace, fpI)
{
const label pointI = curFace[fpI];
label& pl = globalLabel_[pointI];
if (pl == -2)
{
pl = masterPatchPO + glPointLabelsI;
glPointLabelsI++;
}
}
}
}
}
else
{
// Regular patch - second pass
const labelList& patchMeshPoints = pp.meshPoints();
const label patchPO = patchPointOffset_[patchI];
label glPointLabelsI = 0;
forAll (patchMeshPoints, mpI)
{
label& pl = globalLabel_[patchMeshPoints[mpI]];
if (pl == -3)
{
pl = patchPO + glPointLabelsI;
glPointLabelsI++;
}
}
}
}
// "Internal" points - second pass
nInternalPoints = 0;
forAll (globalLabel_, pI)
{
label& pl = globalLabel_[pI];
if (pl == -4)
{
pl = procPointOffset_[Pstream::myProcNo()] + nInternalPoints;
nInternalPoints++;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::globalProcPointIndex::globalProcPointIndex(const polyMesh& mesh)
:
mesh_(mesh),
procPointOffset_(Pstream::nProcs()+1),
patchPointOffset_(mesh_.boundaryMesh().size()+1),
globalLabel_(mesh_.nPoints(), -4)
{
calcPointIndex();
}
// ************************************************************************* //

135
src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.H Normal file
View file

@ -0,0 +1,135 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::globalProcPointIndex
Description
The class creates a unique global point index for each processor points
(pair) in the mesh. Master and slave processor points carry the same index.
Point offsets counts a number of unique points on each processor,
excluding slave processor patch points and global points.
Currently, faces are ordered with internal faces first, followed by patch
faces in patch order, excluding slave processor patches.
If needed, this can be changed to group processor faces with internal faces
to facilitate parallel I/O. HJ, 4/May/2018
Author
Henrik Rusche, Wikki GmbH
SourceFiles
globalProcPointIndex.C
\*---------------------------------------------------------------------------*/
#ifndef globalProcPointIndex_H
#define globalProcPointIndex_H
#include "polyMesh.H"
#include "labelList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class globalProcPointIndex Declaration
\*---------------------------------------------------------------------------*/
class globalProcPointIndex
{
// Private data
//- Mesh reference
const polyMesh& mesh_;
//- Processor point index offset
labelList procPointOffset_;
//- point index offset for patches in local mesh
labelList patchPointOffset_;
//- Global poins label for all points of current mesh
// Sized to number of live points in the mesh
labelList globalLabel_;
// Private Member Functions
//- Disallow default bitwise copy construct
globalProcPointIndex(const globalProcPointIndex&);
//- Disallow default bitwise assignment
void operator=(const globalProcPointIndex&);
//- Calculate point index
void calcPointIndex();
public:
// Constructors
//- Construct from mesh
globalProcPointIndex(const polyMesh&);
//- Destructor - default
// Member Functions
//- Return point index offset per processor
inline const labelList& procPointOffset() const
{
return procPointOffset_;
}
inline const labelList& patchPointOffset() const
{
return patchPointOffset_;
}
//- Return global point label for all points of current mesh
// Sized to number of live points in the mesh
const labelList globalLabel() const
{
return globalLabel_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C
View file

@ -127,7 +127,7 @@ bool immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::update()
{
hasChanged = loadBalance(refinementDict());
}
// Execute dummy mesh motion for the background mesh
const pointField oldPoints = allPoints();
fvMesh::movePoints(oldPoints);

2
src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C
View file

@ -68,7 +68,7 @@ immersedBoundaryEpsilonWallFunctionFvPatchScalarField
p,
Switch(dict.lookup("setDeadValue")),
readScalar(dict.lookup("deadValue"))
)
)
{}

2
src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C
View file

@ -153,7 +153,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
(
ptf.ibPatch(),
ptf.setDeadValue(),
ptf.deadValue()
ptf.deadValue()
)
{
this->setPatchType(ptf);

2
src/lagrangian/basic/CloudDistributeTemplate/CloudDistributeTemplate.H
View file

@ -47,7 +47,7 @@ class Cloud;
class polyMesh;
class cloud;
/*---------------------------------------------------------------------------*\
Class Cloud Declaration
\*---------------------------------------------------------------------------*/

34
src/lagrangian/basic/CloudTemplate/CloudTemplate.C
View file

@ -506,6 +506,40 @@ void Foam::Cloud<ParticleType>::rebuild
}
template<class ParticleType>
Foam::labelList Foam::Cloud<ParticleType>::nParticlesPerCell() const
{
Pout << "In nParticlesPerCell" << endl;
labelList nppc (pMesh().nCells(), 0);
forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
{
const ParticleType& p = pIter();
const label celli = p.cell();
// Check
if (celli < 0 || celli >= pMesh().nCells())
{
FatalErrorIn
(
"Foam::Cloud<ParticleType>::nParticlesPerCell()"
)
<< "Illegal cell number " << celli
<< " at position " << p.position() << nl
<< "Cell number should be between 0 and "
<< pMesh().nCells()-1 << nl
<< "On this mesh the particle should be in cell "
<< pMesh().findCell(p.position())
<< exit(FatalError);
}
nppc[celli]++;
}
return nppc;
}
template<class ParticleType>
void Foam::Cloud<ParticleType>::writePositions() const
{

3
src/lagrangian/basic/CloudTemplate/CloudTemplate.H
View file

@ -290,6 +290,9 @@ public:
const PtrList<labelIOList>& faceProcAddressing
);
//- Count and return number of particles per cell
virtual labelList nParticlesPerCell() const;
// Read

2
src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
View file

@ -370,7 +370,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
// constrainToMeshCentre is inside of an if-statement
mesh.geometricD();
mesh.bounds();
// Prepare for next time step
label parcelsAdded = 0;
scalar massAdded = 0.0;

10
src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C
View file

@ -95,8 +95,8 @@ void Foam::adaptiveOverlapFringe::suitabilityFractionSlope
forAllConstIter(FIFOStack<iterationData>, iterHist, it)
{
n += (it().iteration() - iterMean)*
(it().suitability() - suitabilityMean);
d += sqr(it().iteration() - iterMean);
(it().suitability() - suitabilityMean);
d += Foam::sqr((it().iteration() - iterMean));
}
reduce(n, sumOp<scalar>());
@ -389,9 +389,9 @@ Foam::adaptiveOverlapFringe::adaptiveOverlapFringe
)
:
oversetFringe(mesh, region, dict),
fringeHolesPtr_(nullptr),
acceptorsPtr_(nullptr),
finalDonorAcceptorsPtr_(nullptr),
fringeHolesPtr_(NULL),
acceptorsPtr_(NULL),
finalDonorAcceptorsPtr_(NULL),
holesZoneName_(dict.lookupOrDefault<word>("holes", word())),
initPatchNames_

17
src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
View file

@ -57,8 +57,21 @@ Foam::basicChemistryModel::basicChemistryModel
chemistry_(lookup("chemistry")),
deltaTChem_
(
mesh.nCells(),
readScalar(lookup("initialChemicalTimeStep"))
IOobject
(
"deltaTChem",
mesh.time().constant(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar
(
"deltaTChem0",
dimTime,
readScalar(lookup("initialChemicalTimeStep"))
)
)
{}

6
src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
View file

@ -38,7 +38,7 @@ SourceFiles
#include "IOdictionary.H"
#include "Switch.H"
#include "scalarField.H"
#include "volFields.H"
#include "volFieldsFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -77,7 +77,9 @@ protected:
Switch chemistry_;
//- Latest estimation of integration step
scalarField deltaTChem_;
// HR 14.12.18: Upgraded to volScalarField in order to enable automatic
// mapping/resizing after a topo change. Also needed for PLB.
volScalarField deltaTChem_;
// Protected member functions

2
src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H
View file

@ -298,7 +298,7 @@ public:
) const;
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct() = 0;
virtual void correct();
//- Read RASProperties dictionary
virtual bool read() = 0;

2
src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
View file

@ -207,7 +207,7 @@ public:
virtual tmp<fvVectorMatrix> divDevReff() const = 0;
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct() = 0;
virtual void correct();
//- Read turbulenceProperties dictionary
virtual bool read() = 0;