GIT: Initial state after latest Foundation merge

This commit is contained in:
Andrew Heather
2016-09-20 14:49:08 +01:00
4571 changed files with 115696 additions and 74609 deletions

View File

@ -0,0 +1,987 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015-2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "streamLineBase.H"
#include "fvMesh.H"
#include "ReadFields.H"
#include "sampledSet.H"
#include "globalIndex.H"
#include "mapDistribute.H"
#include "interpolationCellPoint.H"
#include "wallPolyPatch.H"
#include "meshSearchMeshObject.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(streamLineBase, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::autoPtr<Foam::indirectPrimitivePatch>
Foam::streamLineBase::wallPatch() const
{
const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
const polyBoundaryMesh& patches = mesh.boundaryMesh();
label nFaces = 0;
forAll(patches, patchI)
{
//if (!polyPatch::constraintType(patches[patchI].type()))
if (isA<wallPolyPatch>(patches[patchI]))
{
nFaces += patches[patchI].size();
}
}
labelList addressing(nFaces);
nFaces = 0;
forAll(patches, patchI)
{
//if (!polyPatch::constraintType(patches[patchI].type()))
if (isA<wallPolyPatch>(patches[patchI]))
{
const polyPatch& pp = patches[patchI];
forAll(pp, i)
{
addressing[nFaces++] = pp.start()+i;
}
}
}
return autoPtr<indirectPrimitivePatch>
(
new indirectPrimitivePatch
(
IndirectList<face>
(
mesh.faces(),
addressing
),
mesh.points()
)
);
}
void Foam::streamLineBase::initInterpolations
(
const label nSeeds,
label& UIndex,
PtrList<volScalarField>& vsFlds,
PtrList<interpolation<scalar>>& vsInterp,
PtrList<volVectorField>& vvFlds,
PtrList<interpolation<vector>>& vvInterp
)
{
const Time& runTime = obr_.time();
const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
// Read or lookup fields
if (loadFromFiles_)
{
IOobjectList allObjects(mesh, runTime.timeName());
IOobjectList objects(2*fields_.size());
forAll(fields_, i)
{
objects.add(*allObjects[fields_[i]]);
}
ReadFields(mesh, objects, vsFlds);
vsInterp.setSize(vsFlds.size());
forAll(vsFlds, i)
{
vsInterp.set
(
i,
interpolation<scalar>::New
(
interpolationScheme_,
vsFlds[i]
)
);
}
ReadFields(mesh, objects, vvFlds);
vvInterp.setSize(vvFlds.size());
forAll(vvFlds, i)
{
vvInterp.set
(
i,
interpolation<vector>::New
(
interpolationScheme_,
vvFlds[i]
)
);
}
}
else
{
label nScalar = 0;
label nVector = 0;
forAll(fields_, i)
{
if (mesh.foundObject<volScalarField>(fields_[i]))
{
nScalar++;
}
else if (mesh.foundObject<volVectorField>(fields_[i]))
{
nVector++;
}
else
{
FatalErrorInFunction
<< "Cannot find field " << fields_[i] << nl
<< "Valid scalar fields are:"
<< mesh.names(volScalarField::typeName) << nl
<< "Valid vector fields are:"
<< mesh.names(volVectorField::typeName)
<< exit(FatalError);
}
}
vsInterp.setSize(nScalar);
nScalar = 0;
vvInterp.setSize(nVector);
nVector = 0;
forAll(fields_, i)
{
if (mesh.foundObject<volScalarField>(fields_[i]))
{
const volScalarField& f = mesh.lookupObject<volScalarField>
(
fields_[i]
);
vsInterp.set
(
nScalar++,
interpolation<scalar>::New
(
interpolationScheme_,
f
)
);
}
else if (mesh.foundObject<volVectorField>(fields_[i]))
{
const volVectorField& f = mesh.lookupObject<volVectorField>
(
fields_[i]
);
if (f.name() == UName_)
{
UIndex = nVector;
}
vvInterp.set
(
nVector++,
interpolation<vector>::New
(
interpolationScheme_,
f
)
);
}
}
}
// Store the names
scalarNames_.setSize(vsInterp.size());
forAll(vsInterp, i)
{
scalarNames_[i] = vsInterp[i].psi().name();
}
vectorNames_.setSize(vvInterp.size());
forAll(vvInterp, i)
{
vectorNames_[i] = vvInterp[i].psi().name();
}
// Check that we know the index of U in the interpolators.
if (UIndex == -1)
{
FatalErrorInFunction
<< "Cannot find field to move particles with : " << UName_ << nl
<< "This field has to be present in the sampled fields " << fields_
<< " and in the objectRegistry."
<< exit(FatalError);
}
// Sampled data
// ~~~~~~~~~~~~
// Size to maximum expected sizes.
allTracks_.clear();
allTracks_.setCapacity(nSeeds);
allScalars_.setSize(vsInterp.size());
forAll(allScalars_, i)
{
allScalars_[i].clear();
allScalars_[i].setCapacity(nSeeds);
}
allVectors_.setSize(vvInterp.size());
forAll(allVectors_, i)
{
allVectors_[i].clear();
allVectors_[i].setCapacity(nSeeds);
}
}
void Foam::streamLineBase::storePoint
(
const label trackI,
const scalar w,
const label leftI,
const label rightI,
DynamicList<point>& newTrack,
DynamicList<scalarList>& newScalars,
DynamicList<vectorList>& newVectors
) const
{
label sz = newTrack.size();
const List<point>& track = allTracks_[trackI];
newTrack.append((1.0-w)*track[leftI] + w*track[rightI]);
// Scalars
{
newScalars.append(scalarList(allScalars_.size()));
scalarList& newVals = newScalars[sz];
forAll(allScalars_, scalarI)
{
const scalarList& trackVals = allScalars_[scalarI][trackI];
newVals[scalarI] = (1.0-w)*trackVals[leftI] + w*trackVals[rightI];
}
}
// Vectors
{
newVectors.append(vectorList(allVectors_.size()));
vectorList& newVals = newVectors[sz];
forAll(allVectors_, vectorI)
{
const vectorList& trackVals = allVectors_[vectorI][trackI];
newVals[vectorI] = (1.0-w)*trackVals[leftI] + w*trackVals[rightI];
}
}
}
// Can split a track into multiple tracks
void Foam::streamLineBase::trimToBox
(
const treeBoundBox& bb,
const label trackI,
PtrList<DynamicList<point>>& newTracks,
PtrList<DynamicList<scalarList>>& newScalars,
PtrList<DynamicList<vectorList>>& newVectors
) const
{
const List<point>& track = allTracks_[trackI];
if (track.size())
{
for
(
label segmentI = 1;
segmentI < track.size();
segmentI++
)
{
const point& startPt = track[segmentI-1];
const point& endPt = track[segmentI];
const vector d(endPt-startPt);
scalar magD = mag(d);
if (magD > ROOTVSMALL)
{
if (bb.contains(startPt))
{
// Store 1.0*track[segmentI-1]+0*track[segmentI]
storePoint
(
trackI,
0.0,
segmentI-1,
segmentI,
newTracks.last(),
newScalars.last(),
newVectors.last()
);
if (!bb.contains(endPt))
{
point clipPt;
if (bb.intersects(endPt, startPt, clipPt))
{
// End of track. Store point and interpolated
// values
storePoint
(
trackI,
mag(clipPt-startPt)/magD,
segmentI-1,
segmentI,
newTracks.last(),
newScalars.last(),
newVectors.last()
);
newTracks.last().shrink();
newScalars.last().shrink();
newVectors.last().shrink();
}
}
}
else
{
// startPt outside box. New track. Get starting point
point clipPt;
if (bb.intersects(startPt, endPt, clipPt))
{
// New track
newTracks.append
(
new DynamicList<point>(track.size()/10)
);
newScalars.append
(
new DynamicList<scalarList>(track.size()/10)
);
newVectors.append
(
new DynamicList<vectorList>(track.size()/10)
);
// Store point and interpolated values
storePoint
(
trackI,
mag(clipPt-startPt)/magD,
segmentI-1,
segmentI,
newTracks.last(),
newScalars.last(),
newVectors.last()
);
if (!bb.contains(endPt))
{
bb.intersects
(
endPt,
point(clipPt),
clipPt
);
// Store point and interpolated values
storePoint
(
trackI,
mag(clipPt-startPt)/magD,
segmentI-1,
segmentI,
newTracks.last(),
newScalars.last(),
newVectors.last()
);
newTracks.last().shrink();
newScalars.last().shrink();
newVectors.last().shrink();
}
}
}
}
}
// Last point
if (bb.contains(track.last()))
{
storePoint
(
trackI,
1.0,
track.size()-2,
track.size()-1,
newTracks.last(),
newScalars.last(),
newVectors.last()
);
}
}
}
void Foam::streamLineBase::trimToBox(const treeBoundBox& bb)
{
// Storage for new tracks. Per track, per sample the coordinate (newTracks)
// or values for all the sampled fields (newScalars, newVectors)
PtrList<DynamicList<point>> newTracks;
PtrList<DynamicList<scalarList>> newScalars;
PtrList<DynamicList<vectorList>> newVectors;
forAll(allTracks_, trackI)
{
const List<point>& track = allTracks_[trackI];
if (track.size())
{
// New track. Assume it consists of the whole track
newTracks.append(new DynamicList<point>(track.size()));
newScalars.append(new DynamicList<scalarList>(track.size()));
newVectors.append(new DynamicList<vectorList>(track.size()));
// Trim, split and append to newTracks
trimToBox(bb, trackI, newTracks, newScalars, newVectors);
}
}
// Transfer newTracks to allTracks_
allTracks_.setSize(newTracks.size());
forAll(allTracks_, trackI)
{
allTracks_[trackI].transfer(newTracks[trackI]);
}
// Replace track scalars
forAll(allScalars_, scalarI)
{
DynamicList<scalarList>& fieldVals = allScalars_[scalarI];
fieldVals.setSize(newTracks.size());
forAll(fieldVals, trackI)
{
scalarList& trackVals = allScalars_[scalarI][trackI];
trackVals.setSize(newScalars[trackI].size());
forAll(trackVals, sampleI)
{
trackVals[sampleI] = newScalars[trackI][sampleI][scalarI];
}
}
}
// Replace track vectors
forAll(allVectors_, vectorI)
{
DynamicList<vectorList>& fieldVals = allVectors_[vectorI];
fieldVals.setSize(newTracks.size());
forAll(fieldVals, trackI)
{
vectorList& trackVals = allVectors_[vectorI][trackI];
trackVals.setSize(newVectors[trackI].size());
forAll(trackVals, sampleI)
{
trackVals[sampleI] = newVectors[trackI][sampleI][vectorI];
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::streamLineBase::streamLineBase
(
const word& name,
const objectRegistry& obr,
const dictionary& dict,
const bool loadFromFiles
)
:
functionObjectState(obr, name),
dict_(dict),
obr_(obr),
loadFromFiles_(loadFromFiles),
log_(true)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::streamLineBase::~streamLineBase()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::streamLineBase::read(const dictionary& dict)
{
if (active_)
{
log_.readIfPresent("log", dict);
if (log_) Info<< type() << " " << name_ << ":" << nl;
dict.lookup("fields") >> fields_;
if (dict.found("UName"))
{
dict.lookup("UName") >> UName_;
}
else
{
UName_ = "U";
if (dict.found("U"))
{
IOWarningInFunction(dict)
<< "Using deprecated entry \"U\"."
<< " Please use \"UName\" instead."
<< endl;
dict.lookup("U") >> UName_;
}
}
if (findIndex(fields_, UName_) == -1)
{
FatalIOErrorInFunction(dict)
<< "Velocity field for tracking " << UName_
<< " should be present in the list of fields " << fields_
<< exit(FatalIOError);
}
dict.lookup("trackForward") >> trackForward_;
dict.lookup("lifeTime") >> lifeTime_;
if (lifeTime_ < 1)
{
FatalErrorInFunction
<< "Illegal value " << lifeTime_ << " for lifeTime"
<< exit(FatalError);
}
trackLength_ = VGREAT;
if (dict.found("trackLength"))
{
dict.lookup("trackLength") >> trackLength_;
if (log_)
{
Info<< type() << " : fixed track length specified : "
<< trackLength_ << nl << endl;
}
}
bounds_ = boundBox::greatBox;
if (dict.readIfPresent("bounds", bounds_))
{
if (log_) Info<< " clipping all segments to " << bounds_ << nl << endl;
}
interpolationScheme_ = dict.lookupOrDefault
(
"interpolationScheme",
interpolationCellPoint<scalar>::typeName
);
//if (log_) Info<< " using interpolation " << interpolationScheme_
// << endl;
cloudName_ = dict.lookupOrDefault<word>("cloudName", type());
dict.lookup("seedSampleSet") >> seedSet_;
const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
const dictionary& coeffsDict = dict.subDict(seedSet_ + "Coeffs");
sampledSetPtr_ = sampledSet::New
(
seedSet_,
mesh,
meshSearchMeshObject::New(mesh),
coeffsDict
);
coeffsDict.lookup("axis") >> sampledSetAxis_;
scalarFormatterPtr_ = writer<scalar>::New(dict.lookup("setFormat"));
vectorFormatterPtr_ = writer<vector>::New(dict.lookup("setFormat"));
}
}
void Foam::streamLineBase::execute()
{}
void Foam::streamLineBase::end()
{}
void Foam::streamLineBase::timeSet()
{}
void Foam::streamLineBase::write()
{
if (active_)
{
if (log_) Info<< type() << " " << name_ << " output:" << nl;
const Time& runTime = obr_.time();
const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
// Do all injection and tracking
track();
if (Pstream::parRun())
{
// Append slave tracks to master ones
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
globalIndex globalTrackIDs(allTracks_.size());
// Construct a distribution map to pull all to the master.
labelListList sendMap(Pstream::nProcs());
labelListList recvMap(Pstream::nProcs());
if (Pstream::master())
{
// Master: receive all. My own first, then consecutive
// processors.
label trackI = 0;
forAll(recvMap, procI)
{
labelList& fromProc = recvMap[procI];
fromProc.setSize(globalTrackIDs.localSize(procI));
forAll(fromProc, i)
{
fromProc[i] = trackI++;
}
}
}
labelList& toMaster = sendMap[0];
toMaster.setSize(globalTrackIDs.localSize());
forAll(toMaster, i)
{
toMaster[i] = i;
}
const mapDistribute distMap
(
globalTrackIDs.size(),
sendMap.xfer(),
recvMap.xfer()
);
// Distribute the track positions. Note: use scheduled comms
// to prevent buffering.
allTracks_.shrink();
mapDistributeBase::distribute
(
Pstream::scheduled,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(),
false,
distMap.constructMap(),
false,
allTracks_,
flipOp()
);
allTracks_.setCapacity(allTracks_.size());
// Distribute the scalars
forAll(allScalars_, scalarI)
{
allScalars_[scalarI].shrink();
mapDistributeBase::distribute
(
Pstream::scheduled,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(),
false,
distMap.constructMap(),
false,
allScalars_[scalarI],
flipOp()
);
allScalars_[scalarI].setCapacity(allScalars_[scalarI].size());
}
// Distribute the vectors
forAll(allVectors_, vectorI)
{
allVectors_[vectorI].shrink();
mapDistributeBase::distribute
(
Pstream::scheduled,
distMap.schedule(),
distMap.constructSize(),
distMap.subMap(),
false,
distMap.constructMap(),
false,
allVectors_[vectorI],
flipOp()
);
allVectors_[vectorI].setCapacity(allVectors_[vectorI].size());
}
}
// Note: filenames scattered below since used in global call
fileName scalarVtkFile;
fileName vectorVtkFile;
if (Pstream::master())
{
if (bounds_ != boundBox::greatBox)
{
// Clip to bounding box
trimToBox(treeBoundBox(bounds_));
}
label nTracks = 0;
label n = 0;
forAll(allTracks_, trackI)
{
if (allTracks_[trackI].size())
{
nTracks++;
n += allTracks_[trackI].size();
}
}
if (log_)
{
Info<< " Tracks:" << nTracks << nl
<< " Total samples:" << n
<< endl;
}
// Massage into form suitable for writers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Make output directory
fileName vtkPath
(
Pstream::parRun()
? runTime.path()/".."/"postProcessing"/"sets"/name()
: runTime.path()/"postProcessing"/"sets"/name()
);
if (mesh.name() != fvMesh::defaultRegion)
{
vtkPath = vtkPath/mesh.name();
}
vtkPath = vtkPath/mesh.time().timeName();
mkDir(vtkPath);
// Convert track positions (and compact out empty tracks)
PtrList<coordSet> tracks(nTracks);
nTracks = 0;
labelList oldToNewTrack(allTracks_.size(), -1);
forAll(allTracks_, trackI)
{
if (allTracks_[trackI].size())
{
tracks.set
(
nTracks,
new coordSet
(
"track" + Foam::name(nTracks),
sampledSetAxis_ //"xyz"
)
);
oldToNewTrack[trackI] = nTracks;
tracks[nTracks].transfer(allTracks_[trackI]);
nTracks++;
}
}
// Convert scalar values
if (allScalars_.size() > 0)
{
List<List<scalarField>> scalarValues(allScalars_.size());
forAll(allScalars_, scalarI)
{
DynamicList<scalarList>& allTrackVals =
allScalars_[scalarI];
scalarValues[scalarI].setSize(nTracks);
forAll(allTrackVals, trackI)
{
scalarList& vals = allTrackVals[trackI];
if (vals.size())
{
label newTrackI = oldToNewTrack[trackI];
scalarValues[scalarI][newTrackI].transfer(vals);
}
}
}
scalarVtkFile = fileName
(
vtkPath
/ scalarFormatterPtr_().getFileName
(
tracks[0],
scalarNames_
)
);
if (log_) Info
<< " Writing data to " << scalarVtkFile.path() << endl;
scalarFormatterPtr_().write
(
true, // writeTracks
tracks,
scalarNames_,
scalarValues,
OFstream(scalarVtkFile)()
);
}
// Convert vector values
if (allVectors_.size() > 0)
{
List<List<vectorField>> vectorValues(allVectors_.size());
forAll(allVectors_, vectorI)
{
DynamicList<vectorList>& allTrackVals =
allVectors_[vectorI];
vectorValues[vectorI].setSize(nTracks);
forAll(allTrackVals, trackI)
{
vectorList& vals = allTrackVals[trackI];
if (vals.size())
{
label newTrackI = oldToNewTrack[trackI];
vectorValues[vectorI][newTrackI].transfer(vals);
}
}
}
vectorVtkFile = fileName
(
vtkPath
/ vectorFormatterPtr_().getFileName
(
tracks[0],
vectorNames_
)
);
//if (log_) Info<< " Writing vector data to "
// << vectorVtkFile << endl;
vectorFormatterPtr_().write
(
true, // writeTracks
tracks,
vectorNames_,
vectorValues,
OFstream(vectorVtkFile)()
);
}
}
// fileNames are generated on the master but setProperty needs to
// be across all procs
Pstream::scatter(scalarVtkFile);
forAll(scalarNames_, nameI)
{
dictionary propsDict;
propsDict.add("file", scalarVtkFile);
const word& fieldName = scalarNames_[nameI];
setProperty(fieldName, propsDict);
}
Pstream::scatter(vectorVtkFile);
forAll(vectorNames_, nameI)
{
dictionary propsDict;
propsDict.add("file", vectorVtkFile);
const word& fieldName = vectorNames_[nameI];
setProperty(fieldName, propsDict);
}
}
}
void Foam::streamLineBase::updateMesh(const mapPolyMesh&)
{
read(dict_);
}
void Foam::streamLineBase::movePoints(const polyMesh&)
{
// Moving mesh affects the search tree
read(dict_);
}
// ************************************************************************* //