steadyParticleTracks.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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/>.
 
Application
    steadyParticleTracks
 
Description
    Generates a VTK file of particle tracks for cases that were computed using
    a steady-state cloud
 
    Note:
    - Case must be re-constructed (if running in parallel) before use
 
\*---------------------------------------------------------------------------*/
 
#include "argList.H"
#include "Cloud.H"
#include "IOdictionary.H"
#include "fvMesh.H"
#include "Time.H"
#include "timeSelector.H"
#include "OFstream.H"
#include "passiveParticleCloud.H"
 
#include "SortableList.H"
#include "IOobjectList.H"
#include "PtrList.H"
#include "Field.H"
#include "steadyParticleTracksTemplates.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
using namespace Foam;
 
namespace Foam
{
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
label validateFields
(
    const List<word>& userFields,
    const IOobjectList& cloudObjs
)
{
    List<bool> ok(userFields.size(), false);
 
    forAll(userFields, i)
    {
        ok[i] = ok[i] || fieldOk<label>(cloudObjs, userFields[i]);
        ok[i] = ok[i] || fieldOk<scalar>(cloudObjs, userFields[i]);
        ok[i] = ok[i] || fieldOk<vector>(cloudObjs, userFields[i]);
        ok[i] = ok[i] || fieldOk<sphericalTensor>(cloudObjs, userFields[i]);
        ok[i] = ok[i] || fieldOk<symmTensor>(cloudObjs, userFields[i]);
        ok[i] = ok[i] || fieldOk<tensor>(cloudObjs, userFields[i]);
    }
 
    label nOk = 0;
    forAll(ok, i)
    {
        if (ok[i])
        {
            nOk++;
        }
        else
        {
            Info << "\n*** Warning: user specified field '" << userFields[i]
                 << "' unavailable" << endl;
        }
    }
 
    return nOk;
}
 
 
template<>
void writeVTK(OFstream& os, const label& value)
{
    os  << value;
}
 
 
template<>
void writeVTK(OFstream& os, const scalar& value)
{
    os  << value;
}
 
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
int main(int argc, char *argv[])
{
    argList::noParallel();
    timeSelector::addOptions();
    #include "addRegionOption.H"
    #include "addDictOption.H"
 
    #include "setRootCase.H"
 
    #include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
    #include "createNamedMesh.H"
    #include "createFields.H"
 
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
    fileName vtkPath(runTime.path()/"VTK");
    mkDir(vtkPath);
 
    typedef HashTable<label, labelPair, labelPair::Hash<> > trackTableType;
 
    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;
 
        fileName vtkTimePath(runTime.path()/"VTK"/runTime.timeName());
        mkDir(vtkTimePath);
 
        Info<< "    Reading particle positions" << endl;
 
        PtrList<passiveParticle> particles(0);
 
        // transfer particles to (more convenient) list
        {
            passiveParticleCloud ppc(mesh, cloudName);
            Info<< "\n    Read " << returnReduce(ppc.size(), sumOp<label>())
                << " particles" << endl;
 
            particles.setSize(ppc.size());
 
            label i = 0;
            forAllIter(passiveParticleCloud, ppc, iter)
            {
                particles.set(i++, ppc.remove(&iter()));
            }
 
            // myCloud should now be empty
        }
 
        List<label> particleToTrack(particles.size());
        label nTracks = 0;
 
        {
            trackTableType trackTable;
            forAll(particles, i)
            {
                const label origProc = particles[i].origProc();
                const label origId = particles[i].origId();
 
                const trackTableType::const_iterator& iter =
                    trackTable.find(labelPair(origProc, origId));
 
                if (iter == trackTable.end())
                {
                    particleToTrack[i] = nTracks;
                    trackTable.insert(labelPair(origProc, origId), nTracks);
                    nTracks++;
                }
                else
                {
                    particleToTrack[i] = iter();
                }
            }
        }
 
 
        if (nTracks == 0)
        {
            Info<< "\n    No track data" << endl;
        }
        else
        {
            Info<< "\n    Generating " << nTracks << " tracks" << endl;
 
            // determine length of each track
            labelList trackLengths(nTracks, 0);
            forAll(particleToTrack, i)
            {
                const label trackI = particleToTrack[i];
                trackLengths[trackI]++;
            }
 
            // particle "age" property used to sort the tracks
            List<SortableList<scalar> > agePerTrack(nTracks);
            List<List<label> > particleMap(nTracks);
 
            forAll(trackLengths, i)
            {
                const label length = trackLengths[i];
                agePerTrack[i].setSize(length);
                particleMap[i].setSize(length);
            }
 
            // store the particle age per track
            IOobjectList cloudObjs
            (
                mesh,
                runTime.timeName(),
                cloud::prefix/cloudName
            );
 
            // TODO: gather age across all procs
            {
                tmp<scalarField> tage =
                    readParticleField<scalar>("age", cloudObjs);
                const scalarField& age = tage();
                List<label> trackSamples(nTracks, 0);
                forAll(particleToTrack, i)
                {
                    const label trackI = particleToTrack[i];
                    const label sampleI = trackSamples[trackI];
                    agePerTrack[trackI][sampleI] = age[i];
                    particleMap[trackI][sampleI] = i;
                    trackSamples[trackI]++;
                }
                tage.clear();
            }
 
 
            if (Pstream::master())
            {
                OFstream os(vtkTimePath/"particleTracks.vtk");
 
                Info<< "\n    Writing particle tracks to " << os.name() << endl;
 
                label nPoints = sum(trackLengths);
 
                os  << "# vtk DataFile Version 2.0" << nl
                    << "particleTracks" << nl
                    << "ASCII" << nl
                    << "DATASET POLYDATA" << nl
                    << "POINTS " << nPoints << " float" << nl;
 
                Info<< "\n    Writing points" << endl;
 
                {
                    forAll(agePerTrack, i)
                    {
                        agePerTrack[i].sort();
 
                        const labelList& ids = agePerTrack[i].indices();
                        labelList& particleIds = particleMap[i];
 
                        {
                            // update addressing
                            List<label> sortedIds(ids);
                            forAll(sortedIds, j)
                            {
                                sortedIds[j] = particleIds[ids[j]];
                            }
                            particleIds = sortedIds;
                        }
 
                        forAll(ids, j)
                        {
                            const label localId = particleIds[j];
                            const vector& pos = particles[localId].position();
                            os  << pos.x() << ' ' << pos.y() << ' ' << pos.z()
                                << nl;
                        }
                    }
                }
 
 
                // write track (line) connectivity to file
 
                Info<< "\n    Writing track lines" << endl;
                os  << "\nLINES " << nTracks << ' ' << nPoints + nTracks << nl;
 
                // Write ids of track points to file
                {
                    label globalPtI = 0;
                    forAll(particleMap, i)
                    {
                        os  << particleMap[i].size() << nl;
 
                        forAll(particleMap[i], j)
                        {
                            os  << ' ' << globalPtI++;
 
                            if (((j + 1) % 10 == 0) && (j != 0))
                            {
                                os  << nl;
                            }
                        }
 
                        os  << nl;
                    }
                }
 
 
                const label nFields = validateFields(userFields, cloudObjs);
 
                os  << "POINT_DATA " << nPoints << nl
                    << "FIELD attributes " << nFields << nl;
 
                Info<< "\n    Processing fields" << nl << endl;
 
                processFields<label>(os, particleMap, userFields, cloudObjs);
                processFields<scalar>(os, particleMap, userFields, cloudObjs);
                processFields<vector>(os, particleMap, userFields, cloudObjs);
                processFields<sphericalTensor>
                    (os, particleMap, userFields, cloudObjs);
                processFields<symmTensor>
                    (os, particleMap, userFields, cloudObjs);
                processFields<tensor>(os, particleMap, userFields, cloudObjs);
 
            }
        }
        Info<< endl;
    }
 
    Info<< "End" << nl << endl;
 
    return 0;
}
 
 
// ************************************************************************* //