Commit 2f0692e0 authored by Mark Olesen's avatar Mark Olesen
Browse files

Merge commit 'OpenCFD/master' into olesenm

parents 573bd8ae 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;
......
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -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];