Commit 2defba00 authored by Andrew Heather's avatar Andrew Heather
Browse files

ENH: Lagrangian - provided backwards compatibility for cases using the

old "positions" file form

The change to barycentric-based tracking changed the contents of the
cloud "positions" file to a new format comprising the barycentric
co-ordinates and other cell position-based info.  This broke
backwards compatibility, providing no option to restart old cases
(v1706 and earlier), and caused difficulties for dependent code, e.g.
for post-processing utilities that could only infer the contents only
after reading.

The barycentric position info is now written to a file called
"coordinates" with provision to restart old cases for which only the
"positions" file is available. Related utilities, e.g. for parallel
running and data conversion have been updated to be able to support both
file types.

To write the "positions" file by default, use set the following option
in the InfoSwitches section of the controlDict:

    writeLagrangianPositions 1;
parent 87e3da80
......@@ -478,6 +478,7 @@ int main(int argc, char *argv[])
if
(
name == "positions"
|| name == "coordinates"
|| name == "origProcId"
|| name == "origId"
)
......
......@@ -745,12 +745,12 @@ int main(int argc, char *argv[])
false
);
IOobject* positionsPtr = sprayObjs.lookup
(
word("positions")
);
// Note: looking up "positions" for backwards compatibility
IOobject* positionsPtr =
sprayObjs.lookup(word("positions"));
IOobject* coordsPtr = sprayObjs.lookup(word("coordinates"));
if (positionsPtr)
if (positionsPtr || coordsPtr)
{
// Read lagrangian particles
// ~~~~~~~~~~~~~~~~~~~~~~~~~
......
......@@ -586,8 +586,10 @@ int main(int argc, char *argv[])
IOobject* positionsPtr =
sprayObjs.lookup(word("positions"));
IOobject* coordsPtr =
sprayObjs.lookup(word("coordinates"));
if (positionsPtr)
if (coordsPtr || positionsPtr)
{
cloudObjects.insert(cloudDirs[i], sprayObjs);
}
......
......@@ -96,9 +96,14 @@ void Foam::parLagrangianRedistributor::findClouds
cloud::prefix/localCloudDirs[i]
);
if (sprayObjs.lookup(word("positions")))
if
(
sprayObjs.lookup(word("coordinates"))
|| sprayObjs.lookup(word("positions"))
)
{
// One of the objects is positions so must be valid cloud
// One of the objects is coordinates/positions so must be valid
// cloud
label cloudI = findIndex(cloudNames, localCloudDirs[i]);
......@@ -107,7 +112,7 @@ void Foam::parLagrangianRedistributor::findClouds
forAllConstIter(IOobjectList, sprayObjs, iter)
{
const word& name = iter.key();
if (name != "positions")
if (name != "coordinates" && name != "positions")
{
objectNames[cloudI][objectI++] = name;
}
......
......@@ -35,8 +35,8 @@ if (timeDirs.size() && !noLagrangian)
cloudPrefix/cloudName
);
// Clouds always have "positions"
if (cloudObjs.found("positions"))
// Clouds always have "positions" (v1706 and lower) or "coordinates"
if (cloudObjs.found("positions") || cloudObjs.found("coordinates"))
{
// Save the cloud fields on a per cloud basis
auto fieldsPerCloud = cloudFields(cloudName);
......@@ -52,9 +52,10 @@ if (timeDirs.size() && !noLagrangian)
}
}
// Prune out "positions" again since it gets treated specially
// Prune out geometry again since it gets treated specially
forAllIters(cloudFields, cloudIter)
{
cloudIter().erase("coordinates");
cloudIter().erase("positions");
}
......@@ -74,4 +75,5 @@ if (cloudNames.size())
cloudNames.writeList(Info) << endl;
}
// ************************************************************************* //
......@@ -65,7 +65,7 @@ void Foam::ensightSerialCloud::writePositions
forAllConstIter(Cloud<passiveParticle>, cloudPtr(), elmnt)
{
const vector& p = elmnt().position();
const vector p(elmnt().position());
os.write(p.x());
os.write(p.y());
......@@ -79,7 +79,7 @@ void Foam::ensightSerialCloud::writePositions
label parcelId = 0;
forAllConstIter(Cloud<passiveParticle>, cloudPtr(), elmnt)
{
const vector& p = elmnt().position();
const vector p(elmnt().position());
os.write(++parcelId, 8); // unusual width
os.write(p.x());
......
......@@ -60,16 +60,16 @@ if (timeDirs.size())
cloudPrefix/cloudName
);
bool hasPositions = false;
bool hasCoordinates = false;
forAllConstIter(IOobjectList, objs, fieldIter)
{
const IOobject obj = *fieldIter();
const word& fieldName = obj.name();
const word& fieldType = obj.headerClassName();
if (fieldName == "positions")
if (fieldName == "positions" || fieldName == "coordinates")
{
hasPositions = true;
hasCoordinates = true;
}
else if (cloudFieldTypes.found(fieldType))
{
......@@ -79,7 +79,7 @@ if (timeDirs.size())
}
// drop this cloud if it has no positions or is otherwise empty
if (!hasPositions || cloudIter().empty())
if (!hasCoordinates || cloudIter().empty())
{
Info<< "removing cloud " << cloudName << endl;
cloudFields.erase(cloudIter);
......
......@@ -381,8 +381,13 @@ int main(int argc, char *argv[])
cloudPrefix/cloudName
);
// Check that the positions field is present for this time
if (!cloudObjs.found("positions"))
// Check that the positions/coordinates field is present for this
// time
if
(
!cloudObjs.found("positions")
|| !cloudObjs.found("coordinates")
)
{
continue;
}
......
......@@ -62,7 +62,7 @@ for (label i=0; i < nTypes; i++)
if (fieldTypes[i] == cloud::prefix)
{
IOobject lagrangianHeader
IOobject positionsHeader
(
"positions",
runTime.timeName(),
......@@ -71,9 +71,22 @@ for (label i=0; i < nTypes; i++)
IOobject::NO_READ
);
IOobject coordinatesHeader
(
"coordinates",
runTime.timeName(),
cloud::prefix,
mesh,
IOobject::NO_READ
);
if
(
lagrangianHeader.typeHeaderOk<IOPosition<Cloud<passiveParticle>>>
positionsHeader.typeHeaderOk<IOPosition<Cloud<passiveParticle>>>
(
false
)
|| coordinatesHeader.typeHeaderOk<IOPosition<Cloud<passiveParticle>>>
(
false
)
......
......@@ -1168,8 +1168,9 @@ int main(int argc, char *argv[])
);
IOobject* positionsPtr = sprayObjs.lookup("positions");
IOobject* coordinatesPtr = sprayObjs.lookup("coordinates");
if (positionsPtr)
if (positionsPtr || coordinatesPtr)
{
mkDir(fvPath/cloud::prefix/cloudDirs[cloudI]);
......
......@@ -23,10 +23,8 @@ if (timeDirs.size() && !noLagrangian)
fileName::DIRECTORY
);
forAll(cloudDirs, cloudI)
for (const word& cloudName : cloudDirs)
{
const word& cloudName = cloudDirs[cloudI];
IOobjectList cloudObjs
(
mesh,
......@@ -34,8 +32,8 @@ if (timeDirs.size() && !noLagrangian)
cloudPrefix/cloudName
);
// clouds always require "positions"
if (cloudObjs.found("positions"))
// Clouds always require "positions"/"coordinates"
if (cloudObjs.found("positions") || cloudObjs.found("coordinates"))
{
if (allCloudDirs.insert(cloudName))
{
......@@ -60,9 +58,9 @@ if (cloudNames.size())
{
// complete the echo information
Info<< "(";
forAll(cloudNames, cloudNo)
for (const word& cloudName : cloudNames)
{
Info<< ' ' << cloudNames[cloudNo];
Info<< ' ' << cloudName;
}
Info<< " ) " << endl;
}
......
......@@ -1344,7 +1344,7 @@ int main(int argc, char *argv[])
cloud::prefix/cloudName
);
if (sprayObjs.found("positions"))
if (sprayObjs.found("positions") || sprayObjs.found("coordinates"))
{
wordList labelNames(sprayObjs.names(labelIOField::typeName));
Info<< " labels :";
......
......@@ -64,7 +64,8 @@ vtkSmartPointer<vtkPolyData> Foam::vtkPVFoam::lagrangianVTKMesh
);
IOobject* positionsPtr = sprayObjs.lookup(word("positions"));
if (positionsPtr)
IOobject* coordinatesPtr = sprayObjs.lookup(word("coordinates"));
if (positionsPtr || coordinatesPtr)
{
Cloud<passiveParticle> parcels(mesh, cloudName, false);
......
......@@ -160,7 +160,7 @@ int USERD_set_filenames
}
}
IOobject sprayHeader
IOobject positionsHeader
(
"positions",
runTime.timeName(),
......@@ -171,7 +171,22 @@ int USERD_set_filenames
false
);
if (sprayHeader.typeHeaderOk<Cloud<passiveParticle>>(false))
IOobject coordinatesHeader
(
"coordinates",
runTime.timeName(),
cloud::prefix,
runTime,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
);
if
(
positionsHeader.typeHeaderOk<Cloud<passiveParticle>>(false)
|| coordinatesHeader.typeHeaderOk<Cloud<passiveParticle>>(false)
)
{
Info<< "[Found lagrangian]" << endl;
......@@ -181,10 +196,8 @@ int USERD_set_filenames
IOobjectList objects(*meshPtr, runTime.timeName(), cloud::prefix);
lagrangianScalarNames =
objects.names(sprayScalarFieldName);
lagrangianVectorNames =
objects.names(sprayVectorFieldName);
lagrangianScalarNames = objects.names(sprayScalarFieldName);
lagrangianVectorNames = objects.names(sprayVectorFieldName);
isSpray[fieldNames.size()] = true;
......
......@@ -120,8 +120,9 @@ void mapLagrangian(const meshToMesh0& meshToMesh0Interp)
);
IOobject* positionsPtr = objects.lookup("positions");
IOobject* coordinatesPtr = objects.lookup("coordinates");
if (positionsPtr)
if (positionsPtr || coordinatesPtr)
{
Info<< nl << " processing cloud " << cloudDirs[cloudI] << endl;
......
......@@ -111,7 +111,10 @@ void mapLagrangian(const meshToMesh& interp)
bool foundPositions =
returnReduce(objects.found("positions"), orOp<bool>());;
if (foundPositions)
bool foundCoordinates =
returnReduce(objects.found("coordinates"), orOp<bool>());;
if (foundPositions || foundCoordinates)
{
Info<< nl << " processing cloud " << cloudDirs[cloudI] << endl;
......
......@@ -44,6 +44,9 @@ InfoSwitches
writeDictionaries 0;
writeOptionalEntries 0;
// Write lagrangian "positions" file in v1706 format (at earlier)
writeLagrangianPositions 0;
// Report list of slaves/pids used (parallel)
writeSlaves 1;
......
......@@ -63,7 +63,7 @@ Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
const point& tetBasePt = pPts[f[faceBasePtI]];
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
for (label tetPtI = 1; tetPtI < f.size() - 1; ++tetPtI)
{
label facePtI = (tetPtI + faceBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
......@@ -158,7 +158,7 @@ Foam::label Foam::polyMeshTetDecomposition::findBasePoint
const point& tetBasePt = pPts[f[faceBasePtI]];
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
for (label tetPtI = 1; tetPtI < f.size() - 1; ++tetPtI)
{
label facePtI = (tetPtI + faceBasePtI) % f.size();
label otherFacePtI = f.fcIndex(facePtI);
......@@ -219,17 +219,16 @@ Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts
label nInternalFaces = mesh.nInternalFaces();
for (label fI = 0; fI < nInternalFaces; fI++)
for (label fI = 0; fI < nInternalFaces; ++fI)
{
tetBasePtIs[fI] = findSharedBasePoint(mesh, fI, tol, report);
}
pointField neighbourCellCentres(mesh.nFaces() - nInternalFaces);
for(label facei = nInternalFaces; facei < mesh.nFaces(); facei++)
for (label facei = nInternalFaces; facei < mesh.nFaces(); ++facei)
{
neighbourCellCentres[facei - nInternalFaces] =
pC[pOwner[facei]];
neighbourCellCentres[facei - nInternalFaces] = pC[pOwner[facei]];
}
syncTools::swapBoundaryFacePositions(mesh, neighbourCellCentres);
......@@ -250,8 +249,7 @@ Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts
fI++, bFI++
)
{
label patchi =
mesh.boundaryMesh().patchID()[bFI];
label patchi = mesh.boundaryMesh().patchID()[bFI];
if (patches[patchi].coupled())
{
......@@ -381,7 +379,7 @@ bool Foam::polyMeshTetDecomposition::checkFaceTets
// Calculate coupled cell centre
pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
{
neiCc[facei - mesh.nInternalFaces()] = cc[own[facei]];
}
......@@ -526,21 +524,21 @@ bool Foam::polyMeshTetDecomposition::checkFaceTets
Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
(
const polyMesh& mesh,
label fI,
label cI
label facei,
label celli
)
{
const faceList& pFaces = mesh.faces();
const face& f = pFaces[fI];
const face& f = pFaces[facei];
label nTets = f.size() - 2;
List<tetIndices> faceTets(nTets);
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI ++)
for (label tetPti = 1; tetPti < f.size() - 1; ++tetPti)
{
faceTets[tetPtI - 1] = tetIndices(cI, fI, tetPtI);
faceTets[tetPti - 1] = tetIndices(celli, facei, tetPti);
}
return faceTets;
......@@ -550,28 +548,26 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
(
const polyMesh& mesh,
label cI
label celli
)
{
const faceList& pFaces = mesh.faces();
const cellList& pCells = mesh.cells();
const cell& thisCell = pCells[cI];
const cell& thisCell = pCells[celli];
label nTets = 0;
forAll(thisCell, cFI)
for (const label facei : thisCell)
{
nTets += pFaces[thisCell[cFI]].size() - 2;
nTets += pFaces[facei].size() - 2;
}
DynamicList<tetIndices> cellTets(nTets);
forAll(thisCell, cFI)
for (const label facei : thisCell)
{
label fI = thisCell[cFI];
cellTets.append(faceTetIndices(mesh, fI, cI));
cellTets.append(faceTetIndices(mesh, facei, celli));
}
return cellTets;
......@@ -581,27 +577,26 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
Foam::tetIndices Foam::polyMeshTetDecomposition::findTet
(
const polyMesh& mesh,
label cI,
label celli,
const point& pt
)
{
const faceList& pFaces = mesh.faces();
const cellList& pCells = mesh.cells();
const cell& thisCell = pCells[cI];
const cell& thisCell = pCells[celli];
tetIndices tetContainingPt;
forAll(thisCell, cFI)
for (const label facei : thisCell)
{
label fI = thisCell[cFI];
const face& f = pFaces[fI];
const face& f = pFaces[facei];
for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
for (label tetPti = 1; tetPti < f.size() - 1; ++tetPti)
{
// Get tetIndices of face triangle
tetIndices faceTetIs(cI, fI, tetPtI);
tetIndices faceTetIs(celli, facei, tetPti);
// Check if inside
if (faceTetIs.tet(mesh).inside(pt))
......
......@@ -65,10 +65,11 @@ Foam::findCellParticle::findCellParticle
(
const polyMesh& mesh,
Istream& is,
bool readFields
bool readFields,
bool newFormat
)
:
particle(mesh, is, readFields)
particle(mesh, is, readFields, newFormat)
{
if (readFields)
{
......
......@@ -146,7 +146,8 @@ public:
(
const polyMesh& mesh,
Istream& is,
bool readFields = true
bool readFields = true,
bool newFormat = true
);
//- Construct and return a clone
......
......@@ -480,10 +480,11 @@ Foam::wallBoundedParticle::wallBoundedParticle
(
const polyMesh& mesh,
Istream& is,
bool readFields
bool readFields,
bool newFormat
)
:
particle(mesh, is, readFields)
particle(mesh, is, readFields, newFormat)
{
if (readFields)
{
......
......@@ -249,7 +249,8 @@ public:
(
const polyMesh& c,
Istream& is,
bool readFields = true
bool readFields = true,
bool newFormat = true
);
//- Construct copy
......
......@@ -161,10 +161,11 @@ Foam::wallBoundedStreamLineParticle::wallBoundedStreamLineParticle
(
const polyMesh& mesh,
Istream& is,
bool readFields
bool readFields,
bool newFormat
)
:
wallBoundedParticle(mesh, is, readFields)