applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.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     particleTracks
  26
  27 Description
  28     Generates a VTK file of particle tracks for cases that were computed using
  29     a tracked-parcel-type cloud
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "argList.H"
  34 #include "Cloud.H"
  35 #include "IOdictionary.H"
  36 #include "fvMesh.H"
  37 #include "Time.H"
  38 #include "timeSelector.H"
  39 #include "OFstream.H"
  40 #include "passiveParticleCloud.H"
  41
  42 using namespace Foam;
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 int main(int argc, char *argv[])
  47 {
  48     timeSelector::addOptions();
  49     #include "addRegionOption.H"
  50
  51     #include "setRootCase.H"
  52
  53     #include "createTime.H"
  54     instantList timeDirs = timeSelector::select0(runTime, args);
  55     #include "createNamedMesh.H"
  56     #include "createFields.H"
  57
  58     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  59
  60     fileName vtkPath(runTime.path()/"VTK");
  61     mkDir(vtkPath);
  62
  63     Info<< "Scanning times to determine track data for cloud " << cloudName
  64         << nl << endl;
  65
  66     labelList maxIds(Pstream::nProcs(), -1);
  67     forAll(timeDirs, timeI)
  68     {
  69         runTime.setTime(timeDirs[timeI], timeI);
  70         Info<< "Time = " << runTime.timeName() << endl;
  71
  72         Info<< "    Reading particle positions" << endl;
  73         passiveParticleCloud myCloud(mesh, cloudName);
  74         Info<< "    Read " << returnReduce(myCloud.size(), sumOp<label>())
  75             << " particles" << endl;
  76
  77         forAllConstIter(passiveParticleCloud, myCloud, iter)
  78         {
  79             label origId = iter().origId();
  80             label origProc = iter().origProc();
  81
  82             maxIds[origProc] = max(maxIds[origProc], origId);
  83         }
  84     }
  85     Pstream::listCombineGather(maxIds, maxEqOp<label>());
  86     Pstream::listCombineScatter(maxIds);
  87
  88     labelList numIds = maxIds + 1;
  89
  90     Info<< nl << "Particle statistics:" << endl;
  91     forAll(maxIds, procI)
  92     {
  93         Info<< "    Found " << numIds[procI] << " particles originating"
  94             << " from processor " << procI << endl;
  95     }
  96     Info<< nl << endl;
  97
  98
  99     // calc starting ids for particles on each processor
 100     List<label> startIds(numIds.size(), 0);
 101     for (label i = 0; i < numIds.size()-1; i++)
 102     {
 103         startIds[i+1] += startIds[i] + numIds[i];
 104     }
 105     label nParticle = startIds.last() + numIds[startIds.size()-1];
 106
 107
 108
 109     // number of tracks to generate
 110     label nTracks = nParticle/sampleFrequency;
 111
 112     // storage for all particle tracks
 113     List<DynamicList<vector> > allTracks(nTracks);
 114
 115     Info<< "\nGenerating " << nTracks << " particle tracks for cloud "
 116         << cloudName << nl << endl;
 117
 118     forAll(timeDirs, timeI)
 119     {
 120         runTime.setTime(timeDirs[timeI], timeI);
 121         Info<< "Time = " << runTime.timeName() << endl;
 122
 123         List<pointField> allPositions(Pstream::nProcs());
 124         List<labelField> allOrigIds(Pstream::nProcs());
 125         List<labelField> allOrigProcs(Pstream::nProcs());
 126
 127         // Read particles. Will be size 0 if no particles.
 128         Info<< "    Reading particle positions" << endl;
 129         passiveParticleCloud myCloud(mesh, cloudName);
 130
 131         // collect the track data on all processors that have positions
 132         allPositions[Pstream::myProcNo()].setSize
 133         (
 134             myCloud.size(),
 135             point::zero
 136         );
 137         allOrigIds[Pstream::myProcNo()].setSize(myCloud.size(), 0);
 138         allOrigProcs[Pstream::myProcNo()].setSize(myCloud.size(), 0);
 139
 140         label i = 0;
 141         forAllConstIter(passiveParticleCloud, myCloud, iter)
 142         {
 143             allPositions[Pstream::myProcNo()][i] = iter().position();
 144             allOrigIds[Pstream::myProcNo()][i] = iter().origId();
 145             allOrigProcs[Pstream::myProcNo()][i] = iter().origProc();
 146             i++;
 147         }
 148
 149         // collect the track data on the master processor
 150         Pstream::gatherList(allPositions);
 151         Pstream::gatherList(allOrigIds);
 152         Pstream::gatherList(allOrigProcs);
 153
 154         Info<< "    Constructing tracks" << nl << endl;
 155         if (Pstream::master())
 156         {
 157             forAll(allPositions, procI)
 158             {
 159                 forAll(allPositions[procI], i)
 160                 {
 161                     label globalId =
 162                         startIds[allOrigProcs[procI][i]]
 163                       + allOrigIds[procI][i];
 164
 165                     if (globalId % sampleFrequency == 0)
 166                     {
 167                         label trackId = globalId/sampleFrequency;
 168                         if (allTracks[trackId].size() < maxPositions)
 169                         {
 170                             allTracks[trackId].append
 171                             (
 172                                 allPositions[procI][i]
 173                             );
 174                         }
 175                     }
 176                 }
 177             }
 178         }
 179     }
 180
 181     if (Pstream::master())
 182     {
 183         OFstream vtkTracks(vtkPath/"particleTracks.vtk");
 184
 185         Info<< "\nWriting particle tracks to " << vtkTracks.name()
 186             << nl << endl;
 187
 188         // Total number of points in tracks + 1 per track
 189         label nPoints = 0;
 190         forAll(allTracks, trackI)
 191         {
 192             nPoints += allTracks[trackI].size();
 193         }
 194
 195         vtkTracks
 196             << "# vtk DataFile Version 2.0" << nl
 197             << "particleTracks" << nl
 198             << "ASCII" << nl
 199             << "DATASET POLYDATA" << nl
 200             << "POINTS " << nPoints << " float" << nl;
 201
 202         // Write track points to file
 203         forAll(allTracks, trackI)
 204         {
 205             forAll(allTracks[trackI], i)
 206             {
 207                 const vector& pt = allTracks[trackI][i];
 208                 vtkTracks << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
 209             }
 210         }
 211
 212         // write track (line) connectivity to file
 213         vtkTracks << "LINES " << nTracks << ' ' << nPoints+nTracks << nl;
 214
 215         // Write ids of track points to file
 216         label globalPtI = 0;
 217         forAll(allTracks, trackI)
 218         {
 219             vtkTracks << allTracks[trackI].size();
 220
 221             forAll(allTracks[trackI], i)
 222             {
 223                 vtkTracks << ' ' << globalPtI;
 224                 globalPtI++;
 225             }
 226
 227             vtkTracks << nl;
 228         }
 229
 230         Info<< "end" << endl;
 231     }
 232
 233     return 0;
 234 }
 235
 236
 237 // ************************************************************************* //