ENH: Added steady particle tracks utility

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/files

+steadyParticleTracks.C
+
+EXE = $(FOAM_APPBIN)/steadyParticleTracks

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/options

+EXE_INC = \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -llagrangian \
+    -lmeshTools \
+    -lfiniteVolume

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/createFields.H

+word dictName(args.optionLookupOrDefault<word>("dict", "particleTrackDict"));
+
+IOdictionary propsDict
+(
+    IOobject
+    (
+        dictName,
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED
+    )
+);
+
+word cloudName(propsDict.lookup("cloudName"));
+
+List<word> userFields(propsDict.lookup("fields"));
\ No newline at end of file

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/particleTrackDict

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      particleTrackDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+cloudName       reactingCloud1Tracks;
+
+fields ( d U T );
+
+// ************************************************************************* //
+

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+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"
+    argList::validOptions.insert("dict", "");
+
+    #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);
+
+            forAll(trackLengths, i)
+            {
+                const label length = trackLengths[i];
+                agePerTrack[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];
+                    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;
+
+                {
+                    label offset = 0;
+                    forAll(agePerTrack, i)
+                    {
+                        agePerTrack[i].sort();
+                        const labelList& ids = agePerTrack[i].indices();
+
+                        forAll(ids, j)
+                        {
+                            const label localId = offset + ids[j];
+                            const vector& pos = particles[localId].position();
+                            os  << pos.x() << ' ' << pos.y() << ' ' << pos.z()
+                                << nl;
+                        }
+
+                        offset += trackLengths[i];
+                    }
+                }
+
+
+                // 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(agePerTrack, i)
+                    {
+                        os  << agePerTrack[i].size() << nl;
+
+                        forAll(agePerTrack[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, agePerTrack, userFields, cloudObjs);
+                processFields<scalar>(os, agePerTrack, userFields, cloudObjs);
+                processFields<vector>(os, agePerTrack, userFields, cloudObjs);
+                processFields<sphericalTensor>
+                    (os, agePerTrack, userFields, cloudObjs);
+                processFields<symmTensor>
+                    (os, agePerTrack, userFields, cloudObjs);
+                processFields<tensor>(os, agePerTrack, userFields, cloudObjs);
+
+            }
+        }
+        Info<< endl;
+    }
+
+    Info<< "\ndone" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "steadyParticleTracksTemplates.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+bool fieldOk(const IOobjectList& cloudObjs, const word& name)
+{
+    IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
+
+    return (objects.lookup(name) != NULL);
+}
+
+
+template<class Type>
+tmp<Field<Type> > readParticleField
+(
+    const word& name,
+    const IOobjectList cloudObjs
+)
+{
+    IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
+
+    const IOobject* obj = objects.lookup(name);
+    if (obj != NULL)
+    {
+        IOField<Type> newField(*obj);
+        return tmp<Field<Type> >(new Field<Type>(newField.xfer()));
+    }
+
+    Info<< "error: cloud field name " << name << " not found" << endl;
+
+    return Field<Type>::null();
+}
+
+
+template<class Type>
+PtrList<List<Type> > readFields
+(
+    PtrList<List<Type> >& values,
+    const List<word>& fields,
+    const IOobjectList& cloudObjs
+)
+{
+    IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
+
+    label fieldI = 0;
+    forAllConstIter(IOobjectList, objects, iter)
+    {
+        const IOobject& obj = *iter();
+        forAll(fields, j)
+        {
+            if (obj.name() == fields[j])
+            {
+                Info<< "        reading field " << obj.name() << endl;
+                IOField<Type> newField(obj);
+                values.set(fieldI++, new List<Type>(newField.xfer()));
+                break;
+            }
+        }
+    }
+
+    return values;
+}
+
+
+template<class Type>
+void writeVTK(OFstream& os, const Type& value)
+{
+    os  << value.component(0);
+    for (label i=1; i<pTraits<Type>::nComponents; i++)
+    {
+        os  << ' ' << value.component(i);
+    }
+}
+
+
+template<class Type>
+void writeVTKFields
+(
+    OFstream& os,
+    const PtrList<List<Type> >& values,
+    const List<SortableList<scalar> >& agePerTrack,
+    const List<word>& fieldNames
+)
+{
+    label step = max(floor(8/pTraits<Type>::nComponents), 1);
+
+    forAll(values, fieldI)
+    {
+        Info<< "        writing field " << fieldNames[fieldI] << endl;
+        os  << nl << fieldNames[fieldI] << ' ' << pTraits<Type>::nComponents << ' '
+            << values[fieldI].size() << " float" << nl;
+        label offset = 0;
+        forAll(agePerTrack, trackI)
+        {
+            const List<label> ids = agePerTrack[trackI].indices() + offset;
+
+            List<Type> data(UIndirectList<Type>(values[fieldI], ids));
+            label nData = data.size() - 1;
+            forAll(data, i)
+            {
+                writeVTK<Type>(os, data[i]);
+                if (((i + 1) % step == 0) || (i == nData))
+                {
+                    os  << nl;
+                }
+                else
+                {
+                    os  << ' ';
+                }
+            }
+            offset += ids.size();
+        }
+    }
+}
+
+
+template<class Type>
+void processFields
+(
+    OFstream& os,
+    const List<SortableList<scalar> >& agePerTrack,
+    const List<word>& userFieldNames,
+    const IOobjectList& cloudObjs
+)
+{
+    IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
+
+    if (objects.size())
+    {
+        DynamicList<word> fieldNames(objects.size());
+        forAll(userFieldNames, i)
+        {
+            IOobject* obj = objects.lookup(userFieldNames[i]);
+            if (obj != NULL)
+            {
+                fieldNames.append(obj->name());
+            }
+        }
+        fieldNames.shrink();
+
+        PtrList<List<Type> > values(fieldNames.size());
+        readFields<Type>(values, fieldNames, cloudObjs);
+
+        writeVTKFields<Type>
+        (
+            os,
+            values,
+            agePerTrack,
+            fieldNames.xfer()
+        );
+    }
+}
+
+} // End namespace Foam
+
+// ************************************************************************* //