Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

This commit is contained in:
andy
2010-02-08 17:59:06 +00:00
55 changed files with 290252 additions and 1219 deletions

View File

@ -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"

View File

@ -53,7 +53,7 @@ int main(int argc, char *argv[])
const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
// Test:print shared points
// Test:print (collocated) shared points
{
const labelListList& globalPointSlaves =
globalData.globalPointSlaves();
@ -90,7 +90,7 @@ int main(int argc, char *argv[])
// Test: point to faces addressing
// Test: (collocated) point to faces addressing
{
const labelListList& globalPointBoundaryFaces =
globalData.globalPointBoundaryFaces();
@ -137,7 +137,7 @@ int main(int argc, char *argv[])
// Test:point to cells addressing
// Test:(collocated) point to cells addressing
{
const labelList& boundaryCells = globalData.boundaryCells();
const labelListList& globalPointBoundaryCells =
@ -172,7 +172,7 @@ int main(int argc, char *argv[])
// Test:print shared edges
// Test:print (collocated) shared edges
{
const labelListList& globalEdgeSlaves =
globalData.globalEdgeSlaves();

View File

@ -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();

View File

@ -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];
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();
}

View File

@ -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();

View File

@ -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

View File

@ -1,7 +1,7 @@
// ignore special fields or fields that we don't handle
//
bool variableGood = true;
for (label n1=startTime; n1<endTime && variableGood; ++n1)
for (label n1=0; n1<Times.size() && variableGood; ++n1)
{
// ignore _0 fields
if (fieldName.size() > 2 && fieldName(fieldName.size() - 2, 2) == "_0")

View File

@ -19,7 +19,7 @@ if (Pstream::master())
Info<< "Correcting time values. Adding " << Tcorr << endl;
}
for (int n=startTime; n<endTime; n++)
forAll(Times, n)
{
ensightCaseFile << setw(12) << Times[n].value() + Tcorr << " ";

View File

@ -185,7 +185,6 @@ void writeAllFaceData
const labelList& prims,
const label nPrims,
const Field<Type>& pf,
const labelList& patchProcessors,
OFstream& ensightFile
)
{
@ -199,16 +198,12 @@ void writeAllFaceData
{
writeData(map(pf, prims, cmpt), ensightFile);
forAll(patchProcessors, i)
for (int slave=1; slave<Pstream::nProcs(); slave++)
{
if (patchProcessors[i] != 0)
{
label slave = patchProcessors[i];
IPstream fromSlave(Pstream::scheduled, slave);
scalarField pf(fromSlave);
IPstream fromSlave(Pstream::scheduled, slave);
scalarField pf(fromSlave);
writeData(pf, ensightFile);
}
writeData(pf, ensightFile);
}
}
}
@ -231,7 +226,6 @@ void writeAllFaceDataBinary
const labelList& prims,
const label nPrims,
const Field<Type>& pf,
const labelList& patchProcessors,
std::ofstream& ensightFile
)
{
@ -245,16 +239,12 @@ void writeAllFaceDataBinary
{
writeEnsDataBinary(map(pf, prims, cmpt), ensightFile);
forAll(patchProcessors, i)
for (int slave=1; slave<Pstream::nProcs(); slave++)
{
if (patchProcessors[i] != 0)
{
label slave = patchProcessors[i];
IPstream fromSlave(Pstream::scheduled, slave);
scalarField pf(fromSlave);
IPstream fromSlave(Pstream::scheduled, slave);
scalarField pf(fromSlave);
writeEnsDataBinary(pf, ensightFile);
}
writeEnsDataBinary(pf, ensightFile);
}
}
}
@ -278,7 +268,6 @@ bool writePatchField
const Foam::label ensightPatchI,
const Foam::faceSets& boundaryFaceSet,
const Foam::ensightMesh::nFacePrimitives& nfp,
const Foam::labelList& patchProcessors,
Foam::OFstream& ensightFile
)
{
@ -297,7 +286,6 @@ bool writePatchField
boundaryFaceSet.tris,
nfp.nTris,
pf,
patchProcessors,
ensightFile
);
@ -307,7 +295,6 @@ bool writePatchField
boundaryFaceSet.quads,
nfp.nQuads,
pf,
patchProcessors,
ensightFile
);
@ -317,7 +304,6 @@ bool writePatchField
boundaryFaceSet.polys,
nfp.nPolys,
pf,
patchProcessors,
ensightFile
);
@ -338,7 +324,6 @@ bool writePatchFieldBinary
const Foam::label ensightPatchI,
const Foam::faceSets& boundaryFaceSet,
const Foam::ensightMesh::nFacePrimitives& nfp,
const Foam::labelList& patchProcessors,
std::ofstream& ensightFile
)
{
@ -356,7 +341,6 @@ bool writePatchFieldBinary
boundaryFaceSet.tris,
nfp.nTris,
pf,
patchProcessors,
ensightFile
);
@ -366,7 +350,6 @@ bool writePatchFieldBinary
boundaryFaceSet.quads,
nfp.nQuads,
pf,
patchProcessors,
ensightFile
);
@ -376,7 +359,6 @@ bool writePatchFieldBinary
boundaryFaceSet.polys,
nfp.nPolys,
pf,
patchProcessors,
ensightFile
);
@ -406,7 +388,6 @@ void writePatchField
const List<faceSets>& boundaryFaceSets = eMesh.boundaryFaceSets();
const wordList& allPatchNames = eMesh.allPatchNames();
const List<labelList>& allPatchProcs = eMesh.allPatchProcs();
const HashTable<ensightMesh::nFacePrimitives>&
nPatchPrims = eMesh.nPatchPrims();
@ -425,8 +406,6 @@ void writePatchField
}
const labelList& patchProcessors = allPatchProcs[patchi];
word pfName = patchName + '.' + fieldName;
word timeFile = prepend + itoa(timeIndex);
@ -473,7 +452,6 @@ void writePatchField
ensightPatchI,
boundaryFaceSets[patchi],
nPatchPrims.find(patchName)(),
patchProcessors,
ensightFile
);
}
@ -488,7 +466,6 @@ void writePatchField
ensightPatchI,
nullFaceSets,
nPatchPrims.find(patchName)(),
patchProcessors,
ensightFile
);
}
@ -521,7 +498,6 @@ void ensightFieldAscii
const cellSets& meshCellSets = eMesh.meshCellSets();
const List<faceSets>& boundaryFaceSets = eMesh.boundaryFaceSets();
const wordList& allPatchNames = eMesh.allPatchNames();
const List<labelList>& allPatchProcs = eMesh.allPatchProcs();
const wordHashSet& patchNames = eMesh.patchNames();
const HashTable<ensightMesh::nFacePrimitives>&
nPatchPrims = eMesh.nPatchPrims();
@ -619,50 +595,23 @@ void ensightFieldAscii
forAll(allPatchNames, patchi)
{
const word& patchName = allPatchNames[patchi];
const labelList& patchProcessors = allPatchProcs[patchi];
if (patchNames.empty() || patchNames.found(patchName))
{
if (mesh.boundary()[patchi].size())
{
if
if
(
writePatchField
(
writePatchField
(
vf.boundaryField()[patchi],
patchi,
ensightPatchI,
boundaryFaceSets[patchi],
nPatchPrims.find(patchName)(),
patchProcessors,
ensightFile
)
vf.boundaryField()[patchi],
patchi,
ensightPatchI,
boundaryFaceSets[patchi],
nPatchPrims.find(patchName)(),
ensightFile
)
{
ensightPatchI++;
}
}
else if (Pstream::master())
)
{
faceSets nullFaceSet;
if
(
writePatchField
(
Field<Type>(),
-1,
ensightPatchI,
nullFaceSet,
nPatchPrims.find(patchName)(),
patchProcessors,
ensightFile
)
)
{
ensightPatchI++;
}
ensightPatchI++;
}
}
}
@ -695,7 +644,6 @@ void ensightFieldBinary
const cellSets& meshCellSets = eMesh.meshCellSets();
const List<faceSets>& boundaryFaceSets = eMesh.boundaryFaceSets();
const wordList& allPatchNames = eMesh.allPatchNames();
const List<labelList>& allPatchProcs = eMesh.allPatchProcs();
const wordHashSet& patchNames = eMesh.patchNames();
const HashTable<ensightMesh::nFacePrimitives>&
nPatchPrims = eMesh.nPatchPrims();
@ -819,50 +767,23 @@ void ensightFieldBinary
forAll(allPatchNames, patchi)
{
const word& patchName = allPatchNames[patchi];
const labelList& patchProcessors = allPatchProcs[patchi];
if (patchNames.empty() || patchNames.found(patchName))
{
if (mesh.boundary()[patchi].size())
{
if
if
(
writePatchFieldBinary
(
writePatchFieldBinary
(
vf.boundaryField()[patchi],
patchi,
ensightPatchI,
boundaryFaceSets[patchi],
nPatchPrims.find(patchName)(),
patchProcessors,
ensightFile
)
vf.boundaryField()[patchi],
patchi,
ensightPatchI,
boundaryFaceSets[patchi],
nPatchPrims.find(patchName)(),
ensightFile
)
{
ensightPatchI++;
}
}
else if (Pstream::master())
)
{
faceSets nullFaceSet;
if
(
writePatchFieldBinary
(
Field<Type>(),
-1,
ensightPatchI,
nullFaceSet,
nPatchPrims.find(patchName)(),
patchProcessors,
ensightFile
)
)
{
ensightPatchI++;
}
ensightPatchI++;
}
}
}

View File

@ -42,6 +42,7 @@ SourceFiles
#include "fvMesh.H"
#include "OFstream.H"
#include <fstream>
#include "globalIndex.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -62,14 +63,12 @@ public:
{
public:
label nPoints;
label nTris;
label nQuads;
label nPolys;
nFacePrimitives()
:
nPoints(0),
nTris(0),
nQuads(0),
nPolys(0)
@ -95,8 +94,6 @@ private:
wordList allPatchNames_;
List<labelList> allPatchProcs_;
wordHashSet patchNames_;
HashTable<nFacePrimitives> nPatchPrims_;
@ -119,20 +116,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 +154,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,31 +168,21 @@ private:
const char* key,
const label nPrims,
const cellShapeList& cellShapes,
const labelList& pointOffsets,
OFstream& ensightGeometryFile
) const;
void writeFacePrims
(
const faceList& patchFaces,
const label pointOffset,
OFstream& ensightGeometryFile
) const;
faceList map
(
const faceList& patchFaces,
const labelList& prims
) const;
void writeAllFacePrims
(
const char* key,
const labelList& prims,
const label nPrims,
const faceList& patchFaces,
const labelList& pointOffsets,
const labelList& patchProcessors,
OFstream& ensightGeometryFile
) const;
@ -208,7 +195,6 @@ private:
void writeNSidedPoints
(
const faceList& patchFaces,
const label pointOffset,
OFstream& ensightGeometryFile
) const;
@ -217,8 +203,6 @@ private:
const labelList& prims,
const label nPrims,
const faceList& patchFaces,
const labelList& pointOffsets,
const labelList& patchProcessors,
OFstream& ensightGeometryFile
) const;
@ -227,7 +211,10 @@ private:
const fileName& postProcPath,
const word& prepend,
const label timeIndex,
Ostream& ensightCaseFile
Ostream& ensightCaseFile,
const labelList& pointToGlobal,
const pointField& uniquePoints,
const globalIndex& globalPoints
) const;
void writeBinary
@ -235,13 +222,15 @@ private:
const fileName& postProcPath,
const word& prepend,
const label timeIndex,
Ostream& ensightCaseFile
Ostream& ensightCaseFile,
const labelList& pointToGlobal,
const pointField& uniquePoints,
const globalIndex& globalPoints
) const;
void writePrimsBinary
(
const cellShapeList& cellShapes,
const label pointOffset,
std::ofstream& ensightGeometryFile
) const;
@ -250,7 +239,6 @@ private:
const char* key,
const label nPrims,
const cellShapeList& cellShapes,
const labelList& pointOffsets,
std::ofstream& ensightGeometryFile
) const;
@ -274,13 +262,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;
@ -290,22 +277,18 @@ private:
const labelList& prims,
const label nPrims,
const faceList& patchFaces,
const labelList& pointOffsets,
const labelList& patchProcessors,
std::ofstream& ensightGeometryFile
) const;
void writeFacePrimsBinary
(
const faceList& patchFaces,
const label pointOffset,
std::ofstream& ensightGeometryFile
) const;
void writeNSidedPointsBinary
(
const faceList& patchFaces,
const label pointOffset,
std::ofstream& ensightGeometryFile
) const;
@ -320,8 +303,6 @@ private:
const labelList& prims,
const label nPrims,
const faceList& patchFaces,
const labelList& pointOffsets,
const labelList& patchProcessors,
std::ofstream& ensightGeometryFile
) const;
@ -361,11 +342,6 @@ public:
return allPatchNames_;
}
const List<labelList>& allPatchProcs() const
{
return allPatchProcs_;
}
const wordHashSet& patchNames() const
{
return patchNames_;

View File

@ -93,6 +93,9 @@ bool inFileNameList
int main(int argc, char *argv[])
{
timeSelector::addOptions();
# include "addRegionOption.H"
argList::addBoolOption
(
"ascii",
@ -111,7 +114,6 @@ int main(int argc, char *argv[])
"An empty list suppresses writing the internalMesh."
);
# include "addTimeOptions.H"
# include "setRootCase.H"
// Check options
@ -119,12 +121,7 @@ int main(int argc, char *argv[])
# include "createTime.H"
// get the available time-steps
instantList Times = runTime.times();
# include "checkTimeOptions.H"
runTime.setTime(Times[startTime], startTime);
instantList Times = timeSelector::select0(runTime, args);
# include "createNamedMesh.H"
@ -214,9 +211,9 @@ int main(int argc, char *argv[])
// Identify if lagrangian data exists at each time, and add clouds
// to the 'allCloudNames' hash set
for (label n=startTime; n<endTime; n++)
forAll(Times, timeI)
{
runTime.setTime(Times[n], n);
runTime.setTime(Times[timeI], timeI);
fileNameList cloudDirs = readDir
(
@ -267,9 +264,9 @@ int main(int argc, char *argv[])
// Loop over all times to build list of fields and field types
// for each cloud
for (label n=startTime; n<endTime; n++)
forAll(Times, timeI)
{
runTime.setTime(Times[n], n);
runTime.setTime(Times[timeI], timeI);
IOobjectList cloudObjs
(
@ -296,20 +293,19 @@ int main(int argc, char *argv[])
}
label nTimeSteps = 0;
for (label n=startTime; n<endTime; n++)
forAll(Times, timeIndex)
{
nTimeSteps++;
runTime.setTime(Times[n], n);
label timeIndex = n - startTime;
runTime.setTime(Times[timeIndex], timeIndex);
word timeName = itoa(timeIndex);
word timeFile = prepend + timeName;
Info<< "Translating time = " << runTime.timeName() << nl;
# include "moveMesh.H"
polyMesh::readUpdateState meshState = mesh.readUpdate();
if (timeIndex == 0 || mesh.moving())
if (timeIndex == 0 || (meshState != polyMesh::UNCHANGED))
{
eMesh.write
(

View File

@ -1,28 +0,0 @@
{
IOobject ioPoints
(
"points",
runTime.timeName(),
polyMesh::meshSubDir,
mesh
);
if (ioPoints.headerOk())
{
// Reading new points
pointIOField newPoints
(
IOobject
(
"points",
mesh.time().timeName(),
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
mesh.movePoints(newPoints);
}
}

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

File diff suppressed because it is too large Load Diff

View File

@ -41,6 +41,9 @@ Description
#include "OFstream.H"
#include "meshTools.H"
#include "Random.H"
#include "transform.H"
#include "IOmanip.H"
#include "Pair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -355,6 +358,12 @@ int main(int argc, char *argv[])
);
}
if (m < 0)
{
WarningIn(args.executable() + "::main")
<< "Negative mass detected" << endl;
}
vector eVal = eigenValues(J);
tensor eVec = eigenVectors(J);
@ -380,19 +389,221 @@ int main(int argc, char *argv[])
pertI++;
}
Info<< nl
<< "Density = " << density << nl
<< "Mass = " << m << nl
<< "Centre of mass = " << cM << nl
<< "Inertia tensor around centre of mass = " << J << nl
<< "eigenValues (principal moments) = " << eVal << nl
<< "eigenVectors (principal axes) = "
<< eVec.x() << ' ' << eVec.y() << ' ' << eVec.z()
<< endl;
bool showTransform = true;
if
(
(mag(eVec.x() ^ eVec.y()) > (1.0 - SMALL))
&& (mag(eVec.y() ^ eVec.z()) > (1.0 - SMALL))
&& (mag(eVec.z() ^ eVec.x()) > (1.0 - SMALL))
)
{
// Make the eigenvectors a right handed orthogonal triplet
eVec.z() *= sign((eVec.x() ^ eVec.y()) & eVec.z());
// Finding the most natural transformation. Using Lists
// rather than tensors to allow indexed permutation.
// Cartesian basis vectors - right handed orthogonal triplet
List<vector> cartesian(3);
cartesian[0] = vector(1, 0, 0);
cartesian[1] = vector(0, 1, 0);
cartesian[2] = vector(0, 0, 1);
// Principal axis basis vectors - right handed orthogonal
// triplet
List<vector> principal(3);
principal[0] = eVec.x();
principal[1] = eVec.y();
principal[2] = eVec.z();
scalar maxMagDotProduct = -GREAT;
// Matching axis indices, first: cartesian, second:principal
Pair<label> match(-1, -1);
forAll(cartesian, cI)
{
forAll(principal, pI)
{
scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
if (magDotProduct > maxMagDotProduct)
{
maxMagDotProduct = magDotProduct;
match.first() = cI;
match.second() = pI;
}
}
}
scalar sense = sign
(
cartesian[match.first()] & principal[match.second()]
);
if (sense < 0)
{
// Invert the best match direction and swap the order of
// the other two vectors
List<vector> tPrincipal = principal;
tPrincipal[match.second()] *= -1;
tPrincipal[(match.second() + 1) % 3] =
principal[(match.second() + 2) % 3];
tPrincipal[(match.second() + 2) % 3] =
principal[(match.second() + 1) % 3];
principal = tPrincipal;
vector tEVal = eVal;
tEVal[(match.second() + 1) % 3] = eVal[(match.second() + 2) % 3];
tEVal[(match.second() + 2) % 3] = eVal[(match.second() + 1) % 3];
eVal = tEVal;
}
label permutationDelta = match.second() - match.first();
if (permutationDelta != 0)
{
// Add 3 to the permutationDelta to avoid negative indices
permutationDelta += 3;
List<vector> tPrincipal = principal;
vector tEVal = eVal;
for (label i = 0; i < 3; i++)
{
tPrincipal[i] = principal[(i + permutationDelta) % 3];
tEVal[i] = eVal[(i + permutationDelta) % 3];
}
principal = tPrincipal;
eVal = tEVal;
}
label matchedAlready = match.first();
match =Pair<label>(-1, -1);
maxMagDotProduct = -GREAT;
forAll(cartesian, cI)
{
if (cI == matchedAlready)
{
continue;
}
forAll(principal, pI)
{
if (pI == matchedAlready)
{
continue;
}
scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
if (magDotProduct > maxMagDotProduct)
{
maxMagDotProduct = magDotProduct;
match.first() = cI;
match.second() = pI;
}
}
}
sense = sign
(
cartesian[match.first()] & principal[match.second()]
);
if (sense < 0 || (match.second() - match.first()) != 0)
{
principal[match.second()] *= -1;
List<vector> tPrincipal = principal;
tPrincipal[(matchedAlready + 1) % 3] =
principal[(matchedAlready + 2) % 3]*-sense;
tPrincipal[(matchedAlready + 2) % 3] =
principal[(matchedAlready + 1) % 3]*-sense;
principal = tPrincipal;
vector tEVal = eVal;
tEVal[(matchedAlready + 1) % 3] = eVal[(matchedAlready + 2) % 3];
tEVal[(matchedAlready + 2) % 3] = eVal[(matchedAlready + 1) % 3];
eVal = tEVal;
}
eVec.x() = principal[0];
eVec.y() = principal[1];
eVec.z() = principal[2];
// {
// tensor R = rotationTensor(vector(1, 0, 0), eVec.x());
// R = rotationTensor(R & vector(0, 1, 0), eVec.y()) & R;
// Info<< "R = " << nl << R << endl;
// Info<< "R - eVec.T() " << R - eVec.T() << endl;
// }
}
else
{
WarningIn(args.executable() + "::main")
<< "Non-unique eigenvectors, cannot compute transformation "
<< "from Cartesian axes" << endl;
showTransform = false;
}
Info<< nl << setprecision(10)
<< "Density: " << density << nl
<< "Mass: " << m << nl
<< "Centre of mass: " << cM << nl
<< "Inertia tensor around centre of mass: " << nl << J << nl
<< "eigenValues (principal moments): " << eVal << nl
<< "eigenVectors (principal axes): " << nl
<< eVec.x() << nl << eVec.y() << nl << eVec.z() << endl;
if (showTransform)
{
Info<< "Transform tensor from reference state (Q). " << nl
<< "Rotation tensor required to transform "
"from the body reference frame to the global "
"reference frame, i.e.:" << nl
<< "globalVector = Q & bodyLocalVector"
<< nl << eVec.T()
<< endl;
}
if (calcAroundRefPt)
{
Info << "Inertia tensor relative to " << refPt << " = "
Info << "Inertia tensor relative to " << refPt << ": "
<< applyParallelAxisTheorem(m, cM, J, refPt)
<< endl;
}