Commit ea36c756 authored by graham's avatar graham
Browse files

Merge branch 'master' into sixDofPatch

parents cbf807ca 3eb95ac5
......@@ -55,7 +55,7 @@ int main(int argc, char *argv[])
#include "initContinuityErrs.H"
while (runTime.run())
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
......@@ -81,8 +81,6 @@ int main(int argc, char *argv[])
#include "convergenceCheck.H"
}
runTime++;
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
......
......@@ -46,9 +46,6 @@ Usage
@param -fields \n
Use existing geometry decomposition and convert fields only.
@param -filterPatches \n
Remove empty patches when decomposing the geometry.
@param -force \n
Remove any existing @a processor subdirectories before decomposing the
geometry.
......@@ -99,11 +96,6 @@ int main(int argc, char *argv[])
"use existing geometry decomposition and convert fields only"
);
argList::addBoolOption
(
"filterPatches",
"remove empty patches when decomposing the geometry"
);
argList::addBoolOption
(
"force",
"remove existing processor*/ subdirs before decomposing the geometry"
......@@ -128,7 +120,6 @@ int main(int argc, char *argv[])
bool writeCellDist = args.optionFound("cellDist");
bool copyUniform = args.optionFound("copyUniform");
bool decomposeFieldsOnly = args.optionFound("fields");
bool filterPatches = args.optionFound("filterPatches");
bool forceOverwrite = args.optionFound("force");
bool ifRequiredDecomposition = args.optionFound("ifRequired");
......@@ -249,7 +240,7 @@ int main(int argc, char *argv[])
// Decompose the mesh
if (!decomposeFieldsOnly)
{
mesh.decomposeMesh(filterPatches);
mesh.decomposeMesh();
mesh.writeDecomposition();
......
......@@ -69,6 +69,24 @@ void Foam::domainDecomposition::mark
Foam::domainDecomposition::domainDecomposition(const IOobject& io)
:
fvMesh(io),
facesInstancePointsPtr_
(
pointsInstance() != facesInstance()
? new pointIOField
(
IOobject
(
"points",
facesInstance(),
polyMesh::meshSubDir,
*this,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
)
: NULL
),
decompositionDict_
(
IOobject
......@@ -86,7 +104,6 @@ Foam::domainDecomposition::domainDecomposition(const IOobject& io)
procPointAddressing_(nProcs_),
procFaceAddressing_(nProcs_),
procCellAddressing_(nProcs_),
procBoundaryAddressing_(nProcs_),
procPatchSize_(nProcs_),
procPatchStartIndex_(nProcs_),
procNeighbourProcessors_(nProcs_),
......@@ -263,24 +280,67 @@ bool Foam::domainDecomposition::writeDecomposition()
"system",
"constant"
);
processorDb.setTime(time());
// create the mesh
polyMesh procMesh
(
IOobject
// create the mesh. Two situations:
// - points and faces come from the same time ('instance'). The mesh
// will get constructed in the same instance.
// - points come from a different time (moving mesh cases).
// It will read the points belonging to the faces instance and
// construct the procMesh with it which then gets handled as above.
// (so with 'old' geometry).
// Only at writing time will it additionally write the current
// points.
autoPtr<polyMesh> procMeshPtr;
if (facesInstancePointsPtr_.valid())
{
// Construct mesh from facesInstance.
pointField facesInstancePoints
(
this->polyMesh::name(), // region name of undecomposed mesh
pointsInstance(),
processorDb
),
xferMove(procPoints),
xferMove(procFaces),
xferMove(procCells)
);
facesInstancePointsPtr_(),
curPointLabels
);
procMeshPtr.reset
(
new polyMesh
(
IOobject
(
this->polyMesh::name(), // region of undecomposed mesh
facesInstance(),
processorDb
),
xferMove(facesInstancePoints),
xferMove(procFaces),
xferMove(procCells)
)
);
}
else
{
procMeshPtr.reset
(
new polyMesh
(
IOobject
(
this->polyMesh::name(), // region of undecomposed mesh
facesInstance(),
processorDb
),
xferMove(procPoints),
xferMove(procFaces),
xferMove(procCells)
)
);
}
polyMesh& procMesh = procMeshPtr();
// Create processor boundary patches
const labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
// Create processor boundary patches
const labelList& curPatchSizes = procPatchSize_[procI];
const labelList& curPatchStarts = procPatchStartIndex_[procI];
......@@ -309,8 +369,7 @@ bool Foam::domainDecomposition::writeDecomposition()
{
// Get the face labels consistent with the field mapping
// (reuse the patch field mappers)
const polyPatch& meshPatch =
meshPatches[curBoundaryAddressing[patchi]];
const polyPatch& meshPatch = meshPatches[patchi];
fvFieldDecomposer::patchFieldDecomposer patchMapper
(
......@@ -324,14 +383,13 @@ bool Foam::domainDecomposition::writeDecomposition()
);
// Map existing patches
procPatches[nPatches] =
meshPatches[curBoundaryAddressing[patchi]].clone
(
procMesh.boundaryMesh(),
nPatches,
patchMapper.directAddressing(),
curPatchStarts[patchi]
).ptr();
procPatches[nPatches] = meshPatch.clone
(
procMesh.boundaryMesh(),
nPatches,
patchMapper.directAddressing(),
curPatchStarts[patchi]
).ptr();
nPatches++;
}
......@@ -589,6 +647,26 @@ bool Foam::domainDecomposition::writeDecomposition()
procMesh.write();
// Write points if pointsInstance differing from facesInstance
if (facesInstancePointsPtr_.valid())
{
pointIOField pointsInstancePoints
(
IOobject
(
"points",
pointsInstance(),
polyMesh::meshSubDir,
procMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
xferMove(procPoints)
);
pointsInstancePoints.write();
}
Info<< endl
<< "Processor " << procI << nl
<< " Number of cells = " << procMesh.nCells()
......@@ -678,6 +756,16 @@ bool Foam::domainDecomposition::writeDecomposition()
);
cellProcAddressing.write();
// Write patch map for backwards compatibility.
// (= identity map for original patches, -1 for processor patches)
label nMeshPatches = curPatchSizes.size();
labelList procBoundaryAddressing(identity(nMeshPatches));
procBoundaryAddressing.setSize
(
nMeshPatches+curProcessorPatchSizes.size(),
-1
);
labelIOList boundaryProcAddressing
(
IOobject
......@@ -689,7 +777,7 @@ bool Foam::domainDecomposition::writeDecomposition()
IOobject::NO_READ,
IOobject::NO_WRITE
),
procBoundaryAddressing_[procI]
procBoundaryAddressing
);
boundaryProcAddressing.write();
}
......
......@@ -55,6 +55,9 @@ class domainDecomposition
{
// Private data
//- Optional: points at the facesInstance
autoPtr<pointIOField> facesInstancePointsPtr_;
//- Mesh decomposition control dictionary
IOdictionary decompositionDict_;
......@@ -83,9 +86,6 @@ class domainDecomposition
//- Labels of cells for each processor
labelListList procCellAddressing_;
//- Original patch index for every processor patch
labelListList procBoundaryAddressing_;
//- Sizes for processor mesh patches
// Excludes inter-processor boundaries
labelListList procPatchSize_;
......@@ -149,8 +149,8 @@ public:
return distributed_;
}
//- Decompose mesh. Optionally remove zero-sized patches.
void decomposeMesh(const bool filterEmptyPatches);
//- Decompose mesh.
void decomposeMesh();
//- Write decomposition
bool writeDecomposition();
......
......@@ -40,7 +40,7 @@ Description
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
void Foam::domainDecomposition::decomposeMesh()
{
// Decide which cell goes to which processor
distributeCells();
......@@ -598,64 +598,6 @@ void Foam::domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
}
}
Info<< "\nCalculating processor boundary addressing" << endl;
// For every patch of processor boundary, find the index of the original
// patch. Mis-alignment is caused by the fact that patches with zero size
// are omitted. For processor patches, set index to -1.
// At the same time, filter the procPatchSize_ and procPatchStartIndex_
// lists to exclude zero-size patches
forAll(procPatchSize_, procI)
{
// Make a local copy of old lists
const labelList oldPatchSizes = procPatchSize_[procI];
const labelList oldPatchStarts = procPatchStartIndex_[procI];
labelList& curPatchSizes = procPatchSize_[procI];
labelList& curPatchStarts = procPatchStartIndex_[procI];
const labelList& curProcessorPatchSizes =
procProcessorPatchSize_[procI];
labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
curBoundaryAddressing.setSize
(
oldPatchSizes.size()
+ curProcessorPatchSizes.size()
);
label nPatches = 0;
forAll(oldPatchSizes, patchi)
{
if (!filterEmptyPatches || oldPatchSizes[patchi] > 0)
{
curBoundaryAddressing[nPatches] = patchi;
curPatchSizes[nPatches] = oldPatchSizes[patchi];
curPatchStarts[nPatches] = oldPatchStarts[patchi];
nPatches++;
}
}
// reset to the size of live patches
curPatchSizes.setSize(nPatches);
curPatchStarts.setSize(nPatches);
forAll(curProcessorPatchSizes, procPatchI)
{
curBoundaryAddressing[nPatches] = -1;
nPatches++;
}
curBoundaryAddressing.setSize(nPatches);
}
Info<< "\nDistributing points to processors" << endl;
// For every processor, loop through the list of faces for the processor.
// For every face, loop through the list of points and mark the point as
......
......@@ -42,6 +42,7 @@ SourceFiles
#include "fvMesh.H"
#include "OFstream.H"
#include <fstream>
#include "globalIndex.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -110,6 +111,16 @@ private:
//- Disallow default bitwise assignment
void operator=(const ensightMesh&);
//- Construct map from mesh points to merged points.
// pointToGlobal : from mesh point to global point
// uniquePoints : my set of unique points
globalIndex mergeMeshPoints
(
labelList& pointToGlobal,
pointField& uniquePoints
) const;
void writePoints
(
const scalarField& pointsComponent,
......@@ -119,20 +130,21 @@ private:
cellShapeList map
(
const cellShapeList& cellShapes,
const labelList& prims
const labelList& prims,
const labelList& pointToGlobal
) const;
cellShapeList map
(
const cellShapeList& cellShapes,
const labelList& hexes,
const labelList& wedges
const labelList& wedges,
const labelList& pointToGlobal
) const;
void writePrims
(
const cellShapeList& cellShapes,
const label pointOffset,
OFstream& ensightGeometryFile
) const;
......@@ -156,13 +168,12 @@ private:
const labelList& polys,
const cellList& cellFaces,
const faceList& faces,
const label pointOffset,
OFstream& ensightGeometryFile
) const;
void writeAllPolys
(
const labelList& pointOffsets,
const labelList& pointToGlobal,
OFstream& ensightGeometryFile
) const;
......@@ -171,7 +182,6 @@ private:
const char* key,
const label nPrims,
const cellShapeList& cellShapes,
const labelList& pointOffsets,
OFstream& ensightGeometryFile
) const;
......@@ -182,12 +192,6 @@ private:
OFstream& ensightGeometryFile
) const;
faceList map
(
const faceList& patchFaces,
const labelList& prims
) const;
void writeAllFacePrims
(
const char* key,
......@@ -241,7 +245,6 @@ private:
void writePrimsBinary
(
const cellShapeList& cellShapes,
const label pointOffset,
std::ofstream& ensightGeometryFile
) const;
......@@ -250,7 +253,6 @@ private:
const char* key,
const label nPrims,
const cellShapeList& cellShapes,
const labelList& pointOffsets,
std::ofstream& ensightGeometryFile
) const;
......@@ -274,13 +276,12 @@ private:
const labelList& polys,
const cellList& cellFaces,
const faceList& faces,
const label pointOffset,
std::ofstream& ensightGeometryFile
) const;
void writeAllPolysBinary
(
const labelList& pointOffsets,
const labelList& pointToGlobal,
std::ofstream& ensightGeometryFile
) const;
......
......@@ -676,6 +676,10 @@ Foam::Time& Foam::Time::operator++()
deltaT0_ = deltaTSave_;
deltaTSave_ = deltaT_;
// Save old time name
const word oldTimeName = dimensionedScalar::name();
setTime(value() + deltaT_, timeIndex_ + 1);
if (!subCycling_)
......@@ -685,7 +689,30 @@ Foam::Time& Foam::Time::operator++()
{
setTime(0.0, timeIndex_);
}
}
// Check that new time representation differs from old one
if (dimensionedScalar::name() == oldTimeName)
{
int oldPrecision = precision_;
do
{
precision_++;
setTime(value(), timeIndex());
}
while (precision_ < 100 && dimensionedScalar::name() == oldTimeName);
WarningIn("Time::operator++()")
<< "Increased the timePrecision from " << oldPrecision
<< " to " << precision_
<< " to distinguish between timeNames at time " << value()
<< endl;
}
if (!subCycling_)
{
switch (writeControl_)
{
case wcTimeStep:
......
......@@ -410,33 +410,30 @@ Foam::label Foam::globalMeshData::countCoincidentFaces
}
void Foam::globalMeshData::calcGlobalPointSlaves() const
void Foam::globalMeshData::calcGlobalPointSlaves
(
const globalPoints& globalData,
autoPtr<globalIndex>& globalIndicesPtr,
autoPtr<labelListList>& globalPointSlavesPtr,
autoPtr<mapDistribute>& globalPointSlavesMapPtr
) const
{
if (debug)
{
Pout<< "globalMeshData::calcGlobalPointSlaves() :"
<< " calculating coupled master to slave point addressing."
<< endl;
}
// Calculate connected points for master points
globalPoints globalData(mesh_, coupledPatch(), true);
const Map<label>& meshToProcPoint = globalData.meshToProcPoint();
// Create global numbering for coupled points
globalPointNumberingPtr_.reset
globalIndicesPtr.reset
(
new globalIndex(globalData.globalIndices())
);
const globalIndex& globalIndices = globalPointNumberingPtr_();
const globalIndex& globalIndices = globalIndicesPtr();
// Create master to slave addressing. Empty for slave points.
globalPointSlavesPtr_.reset
globalPointSlavesPtr.reset
(
new labelListList(coupledPatch().nPoints())
);
labelListList& globalPointSlaves = globalPointSlavesPtr_();
labelListList& globalPointSlaves = globalPointSlavesPtr();
const Map<label>& meshToProcPoint = globalData.meshToProcPoint();
forAllConstIter(Map<label>, meshToProcPoint, iter)
{
......@@ -465,7 +462,7 @@ void Foam::globalMeshData::calcGlobalPointSlaves() const
// Changes globalPointSlaves to be indices into compact data
List<Map<label> > compactMap(Pstream::nProcs());
globalPointSlavesMapPtr_.reset
globalPointSlavesMapPtr.reset
(
new mapDistribute
(
......@@ -477,41 +474,50 @@ void Foam::globalMeshData::calcGlobalPointSlaves() const
if (debug)
{
Pout<< "globalMeshData::calcGlobalPointSlaves() :"
Pout<< "globalMeshData::calcGlobalPointSlaves(..) :"
<< " coupled points:" << coupledPatch().nPoints()
<< " additional remote points:"
<< globalPointSlavesMapPtr_().constructSize()
<< globalPointSlavesMapPtr().constructSize()
- coupledPatch().nPoints()
<< endl;
}
}
void Foam::globalMeshData::calcGlobalEdgeSlaves() const
void Foam::globalMeshData::calcGlobalPointSlaves() const
{
if (debug)
{
Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
<< " calculating coupled master to slave edge addressing."
Pout<< "globalMeshData::calcGlobalPointSlaves() :"
<< " calculating coupled master to collocated"
<< " slave point addressing."
<< endl;
}
const labelListList& globalPointSlaves