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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     steadyParticleTracks
  26
  27 Description
  28     Generates a VTK file of particle tracks for cases that were computed using
  29     a steady-state cloud
  30     NOTE: case must be re-constructed (if running in parallel) before use
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "argList.H"
  35 #include "Cloud.H"
  36 #include "IOdictionary.H"
  37 #include "fvMesh.H"
  38 #include "Time.H"
  39 #include "timeSelector.H"
  40 #include "OFstream.H"
  41 #include "passiveParticleCloud.H"
  42
  43 #include "SortableList.H"
  44 #include "IOobjectList.H"
  45 #include "PtrList.H"
  46 #include "Field.H"
  47 #include "steadyParticleTracksTemplates.H"
  48
  49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  50
  51 using namespace Foam;
  52
  53 namespace Foam
  54 {
  55
  56 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  57
  58 label validateFields
  59 (
  60     const List<word>& userFields,
  61     const IOobjectList& cloudObjs
  62 )
  63 {
  64     List<bool> ok(userFields.size(), false);
  65
  66     forAll(userFields, i)
  67     {
  68         ok[i] = ok[i] || fieldOk<label>(cloudObjs, userFields[i]);
  69         ok[i] = ok[i] || fieldOk<scalar>(cloudObjs, userFields[i]);
  70         ok[i] = ok[i] || fieldOk<vector>(cloudObjs, userFields[i]);
  71         ok[i] = ok[i] || fieldOk<sphericalTensor>(cloudObjs, userFields[i]);
  72         ok[i] = ok[i] || fieldOk<symmTensor>(cloudObjs, userFields[i]);
  73         ok[i] = ok[i] || fieldOk<tensor>(cloudObjs, userFields[i]);
  74     }
  75
  76     label nOk = 0;
  77     forAll(ok, i)
  78     {
  79         if (ok[i])
  80         {
  81             nOk++;
  82         }
  83         else
  84         {
  85             Info << "\n*** Warning: user specified field '" << userFields[i]
  86                  << "' unavailable" << endl;
  87         }
  88     }
  89
  90     return nOk;
  91 }
  92
  93
  94 template<>
  95 void writeVTK(OFstream& os, const label& value)
  96 {
  97     os  << value;
  98 }
  99
 100
 101 template<>
 102 void writeVTK(OFstream& os, const scalar& value)
 103 {
 104     os  << value;
 105 }
 106
 107 }
 108
 109 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 110
 111 int main(int argc, char *argv[])
 112 {
 113     argList::noParallel();
 114     timeSelector::addOptions();
 115     #include "addRegionOption.H"
 116     argList::validOptions.insert("dict", "");
 117
 118     #include "setRootCase.H"
 119
 120     #include "createTime.H"
 121     instantList timeDirs = timeSelector::select0(runTime, args);
 122     #include "createNamedMesh.H"
 123     #include "createFields.H"
 124
 125     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 126
 127     fileName vtkPath(runTime.path()/"VTK");
 128     mkDir(vtkPath);
 129
 130     typedef HashTable<label, labelPair, labelPair::Hash<> > trackTableType;
 131
 132     forAll(timeDirs, timeI)
 133     {
 134         runTime.setTime(timeDirs[timeI], timeI);
 135         Info<< "Time = " << runTime.timeName() << endl;
 136
 137         fileName vtkTimePath(runTime.path()/"VTK"/runTime.timeName());
 138         mkDir(vtkTimePath);
 139
 140         Info<< "    Reading particle positions" << endl;
 141
 142         PtrList<passiveParticle> particles(0);
 143
 144         // transfer particles to (more convenient) list
 145         {
 146             passiveParticleCloud ppc(mesh, cloudName);
 147             Info<< "\n    Read " << returnReduce(ppc.size(), sumOp<label>())
 148                 << " particles" << endl;
 149
 150             particles.setSize(ppc.size());
 151
 152             label i = 0;
 153             forAllIter(passiveParticleCloud, ppc, iter)
 154             {
 155                 particles.set(i++, ppc.remove(&iter()));
 156             }
 157
 158             // myCloud should now be empty
 159         }
 160
 161         List<label> particleToTrack(particles.size());
 162         label nTracks = 0;
 163
 164         {
 165             trackTableType trackTable;
 166             forAll(particles, i)
 167             {
 168                 const label origProc = particles[i].origProc();
 169                 const label origId = particles[i].origId();
 170
 171                 const trackTableType::const_iterator& iter =
 172                     trackTable.find(labelPair(origProc, origId));
 173
 174                 if (iter == trackTable.end())
 175                 {
 176                     particleToTrack[i] = nTracks;
 177                     trackTable.insert(labelPair(origProc, origId), nTracks);
 178                     nTracks++;
 179                 }
 180                 else
 181                 {
 182                     particleToTrack[i] = iter();
 183                 }
 184             }
 185         }
 186
 187         if (nTracks == 0)
 188         {
 189             Info<< "\n    No track data" << endl;
 190         }
 191         else
 192         {
 193             Info<< "\n    Generating " << nTracks << " tracks" << endl;
 194
 195             // determine length of each track
 196             labelList trackLengths(nTracks, 0);
 197             forAll(particleToTrack, i)
 198             {
 199                 const label trackI = particleToTrack[i];
 200                 trackLengths[trackI]++;
 201             }
 202
 203             // particle "age" property used to sort the tracks
 204             List<SortableList<scalar> > agePerTrack(nTracks);
 205
 206             forAll(trackLengths, i)
 207             {
 208                 const label length = trackLengths[i];
 209                 agePerTrack[i].setSize(length);
 210             }
 211
 212             // store the particle age per track
 213             IOobjectList cloudObjs
 214             (
 215                 mesh,
 216                 runTime.timeName(),
 217                 cloud::prefix/cloudName
 218             );
 219
 220             // TODO: gather age across all procs
 221             {
 222                 tmp<scalarField> tage =
 223                     readParticleField<scalar>("age", cloudObjs);
 224                 const scalarField& age = tage();
 225                 List<label> trackSamples(nTracks, 0);
 226                 forAll(particleToTrack, i)
 227                 {
 228                     const label trackI = particleToTrack[i];
 229                     const label sampleI = trackSamples[trackI];
 230                     agePerTrack[trackI][sampleI] = age[i];
 231                     trackSamples[trackI]++;
 232                 }
 233                 tage.clear();
 234             }
 235
 236
 237             if (Pstream::master())
 238             {
 239                 OFstream os(vtkTimePath/"particleTracks.vtk");
 240
 241                 Info<< "\n    Writing particle tracks to " << os.name() << endl;
 242
 243                 label nPoints = sum(trackLengths);
 244
 245                 os  << "# vtk DataFile Version 2.0" << nl
 246                     << "particleTracks" << nl
 247                     << "ASCII" << nl
 248                     << "DATASET POLYDATA" << nl
 249                     << "POINTS " << nPoints << " float" << nl;
 250
 251                 Info<< "\n    Writing points" << endl;
 252
 253                 {
 254                     label offset = 0;
 255                     forAll(agePerTrack, i)
 256                     {
 257                         agePerTrack[i].sort();
 258                         const labelList& ids = agePerTrack[i].indices();
 259
 260                         forAll(ids, j)
 261                         {
 262                             const label localId = offset + ids[j];
 263                             const vector& pos = particles[localId].position();
 264                             os  << pos.x() << ' ' << pos.y() << ' ' << pos.z()
 265                                 << nl;
 266                         }
 267
 268                         offset += trackLengths[i];
 269                     }
 270                 }
 271
 272
 273                 // write track (line) connectivity to file
 274
 275                 Info<< "\n    Writing track lines" << endl;
 276                 os  << "\nLINES " << nTracks << ' ' << nPoints + nTracks << nl;
 277
 278                 // Write ids of track points to file
 279                 {
 280                     label globalPtI = 0;
 281                     forAll(agePerTrack, i)
 282                     {
 283                         os  << agePerTrack[i].size() << nl;
 284
 285                         forAll(agePerTrack[i], j)
 286                         {
 287                             os  << ' ' << globalPtI++;
 288                             if (((j + 1) % 10 == 0) && (j != 0))
 289                             {
 290                                 os  << nl;
 291                             }
 292                         }
 293
 294                         os  << nl;
 295                     }
 296                 }
 297
 298
 299                 const label nFields = validateFields(userFields, cloudObjs);
 300
 301                 os  << "POINT_DATA " << nPoints << nl
 302                     << "FIELD attributes " << nFields << nl;
 303
 304                 Info<< "\n    Processing fields" << nl << endl;
 305
 306                 processFields<label>(os, agePerTrack, userFields, cloudObjs);
 307                 processFields<scalar>(os, agePerTrack, userFields, cloudObjs);
 308                 processFields<vector>(os, agePerTrack, userFields, cloudObjs);
 309                 processFields<sphericalTensor>
 310                     (os, agePerTrack, userFields, cloudObjs);
 311                 processFields<symmTensor>
 312                     (os, agePerTrack, userFields, cloudObjs);
 313                 processFields<tensor>(os, agePerTrack, userFields, cloudObjs);
 314
 315             }
 316         }
 317         Info<< endl;
 318     }
 319
 320     Info<< "\ndone" << endl;
 321
 322     return 0;
 323 }
 324
 325
 326 // ************************************************************************* //