FreeFOAM The Cross-Platform CFD Toolkit
particleTracks.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd.
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 Usage
32  - particleTracks [OPTION]
33 
34  @param -region <name>\n
35  Only apply to named mesh region.
36 
37  @param -noZero \n
38  Ignore timestep 0.
39 
40  @param -constant \n
41  Include the constant directory.
42 
43  @param -time <time>\n
44  Apply only to specific time.
45 
46  @param -latestTime \n
47  Only apply to latest time step.
48 
49  @param -case <dir> \n
50  Specify the case directory
51 
52  @param -parallel \n
53  Run the case in parallel
54 
55  @param -help \n
56  Display short usage message
57 
58  @param -doc \n
59  Display Doxygen documentation page
60 
61  @param -srcDoc \n
62  Display source code
63 
64 \*---------------------------------------------------------------------------*/
65 
66 #include <OpenFOAM/argList.H>
67 #include <lagrangian/Cloud.H>
68 #include <OpenFOAM/IOdictionary.H>
69 #include <finiteVolume/fvMesh.H>
70 #include <OpenFOAM/Time.H>
71 #include <OpenFOAM/timeSelector.H>
72 #include <OpenFOAM/OFstream.H>
74 
75 using namespace Foam;
76 
77 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
78 
79 int main(int argc, char *argv[])
80 {
83 
84  #include <OpenFOAM/setRootCase.H>
85 
86  #include <OpenFOAM/createTime.H>
89  #include "createFields.H"
90 
91  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
92 
93  fileName vtkPath(runTime.path()/"VTK");
94  mkDir(vtkPath);
95 
96  Info<< "Scanning times to determine track data for cloud " << cloudName
97  << nl << endl;
98 
99  labelList maxIds(Pstream::nProcs(), -1);
100  forAll(timeDirs, timeI)
101  {
102  runTime.setTime(timeDirs[timeI], timeI);
103  Info<< "Time = " << runTime.timeName() << endl;
104 
105  Info<< " Reading particle positions" << endl;
107  Info<< " Read " << returnReduce(myCloud.size(), sumOp<label>())
108  << " particles" << endl;
109 
110  forAllConstIter(passiveParticleCloud, myCloud, iter)
111  {
112  label origId = iter().origId();
113  label origProc = iter().origProc();
114 
115  maxIds[origProc] = max(maxIds[origProc], origId);
116  }
117  }
120 
121  labelList numIds = maxIds + 1;
122 
123  Info<< nl << "Particle statistics:" << endl;
124  forAll(maxIds, procI)
125  {
126  Info<< " Found " << numIds[procI] << " particles originating"
127  << " from processor " << procI << endl;
128  }
129  Info<< nl << endl;
130 
131 
132  // calc starting ids for particles on each processor
133  List<label> startIds(numIds.size(), 0);
134  for (label i = 0; i < numIds.size()-1; i++)
135  {
136  startIds[i+1] += startIds[i] + numIds[i];
137  }
138  label nParticle = startIds[startIds.size()-1] + numIds[startIds.size()-1];
139 
140 
141 
142  // number of tracks to generate
143  label nTracks = nParticle/sampleFrequency;
144 
145  // storage for all particle tracks
146  List<DynamicList<vector> > allTracks(nTracks);
147 
148  Info<< "\nGenerating " << nTracks << " particle tracks for cloud "
149  << cloudName << nl << endl;
150 
151  forAll(timeDirs, timeI)
152  {
153  runTime.setTime(timeDirs[timeI], timeI);
154  Info<< "Time = " << runTime.timeName() << endl;
155 
156  List<pointField> allPositions(Pstream::nProcs());
157  List<labelField> allOrigIds(Pstream::nProcs());
158  List<labelField> allOrigProcs(Pstream::nProcs());
159 
160  // Read particles. Will be size 0 if no particles.
161  Info<< " Reading particle positions" << endl;
163 
164  // collect the track data on all processors that have positions
165  allPositions[Pstream::myProcNo()].setSize
166  (
167  myCloud.size(),
169  );
170  allOrigIds[Pstream::myProcNo()].setSize(myCloud.size(), 0);
171  allOrigProcs[Pstream::myProcNo()].setSize(myCloud.size(), 0);
172 
173  label i = 0;
174  forAllConstIter(passiveParticleCloud, myCloud, iter)
175  {
176  allPositions[Pstream::myProcNo()][i] = iter().position();
177  allOrigIds[Pstream::myProcNo()][i] = iter().origId();
178  allOrigProcs[Pstream::myProcNo()][i] = iter().origProc();
179  i++;
180  }
181 
182  // collect the track data on the master processor
183  Pstream::gatherList(allPositions);
184  Pstream::gatherList(allOrigIds);
185  Pstream::gatherList(allOrigProcs);
186 
187  Info<< " Constructing tracks" << nl << endl;
188  if (Pstream::master())
189  {
190  forAll(allPositions, procI)
191  {
192  forAll(allPositions[procI], i)
193  {
194  label globalId =
195  startIds[allOrigProcs[procI][i]]
196  + allOrigIds[procI][i];
197 
198  if (globalId % sampleFrequency == 0)
199  {
200  label trackId = globalId/sampleFrequency;
201  if (allTracks[trackId].size() < maxPositions)
202  {
203  allTracks[trackId].append
204  (
205  allPositions[procI][i]
206  );
207  }
208  }
209  }
210  }
211  }
212  }
213 
214  if (Pstream::master())
215  {
216  OFstream vtkTracks(vtkPath/"particleTracks.vtk");
217 
218  Info<< "\nWriting particle tracks to " << vtkTracks.name()
219  << nl << endl;
220 
221  // Total number of points in tracks + 1 per track
222  label nPoints = 0;
223  forAll(allTracks, trackI)
224  {
225  nPoints += allTracks[trackI].size();
226  }
227 
228  vtkTracks
229  << "# vtk DataFile Version 2.0" << nl
230  << "particleTracks" << nl
231  << "ASCII" << nl
232  << "DATASET POLYDATA" << nl
233  << "POINTS " << nPoints << " float" << nl;
234 
235  // Write track points to file
236  forAll(allTracks, trackI)
237  {
238  forAll(allTracks[trackI], i)
239  {
240  const vector& pt = allTracks[trackI][i];
241  vtkTracks << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
242  }
243  }
244 
245  // write track (line) connectivity to file
246  vtkTracks << "LINES " << nTracks << ' ' << nPoints+nTracks << nl;
247 
248  // Write ids of track points to file
249  label globalPtI = 0;
250  forAll(allTracks, trackI)
251  {
252  vtkTracks << allTracks[trackI].size();
253 
254  forAll(allTracks[trackI], i)
255  {
256  vtkTracks << ' ' << globalPtI;
257  globalPtI++;
258  }
259 
260  vtkTracks << nl;
261  }
262 
263  Info<< "end" << endl;
264  }
265 
266  return 0;
267 }
268 
269 
270 // ************************ vim: set sw=4 sts=4 et: ************************ //