steadyParticleTracks.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 steadyParticleTracks
29
30Group
31 grpPostProcessingUtilities
32
33Description
34 Generate a legacy VTK file of particle tracks for cases that were
35 computed using a steady-state cloud
36
37 Note:
38 - Case must be re-constructed (if running in parallel) before use
39
40\*---------------------------------------------------------------------------*/
41
42#include "argList.H"
43#include "Cloud.H"
44#include "IOdictionary.H"
45#include "fvMesh.H"
46#include "Time.H"
47#include "timeSelector.H"
48#include "OFstream.H"
49#include "labelPairHashes.H"
50#include "IOField.H"
51#include "IOobjectList.H"
52#include "SortableList.H"
55
56// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
57
58namespace Foam
59{
60
61// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
62
63// Extract list of IOobjects, modifying the input IOobjectList in the
64// process
65IOobjectList preFilterFields
66(
67 IOobjectList& cloudObjects,
68 const wordRes& acceptFields,
70)
71{
72 IOobjectList filteredObjects(cloudObjects.capacity());
74
75 // Selection here is slighly different than usual
76 // - an empty accept filter means "ignore everything"
77
78 if (!acceptFields.empty())
79 {
81
82 const wordList allNames(cloudObjects.sortedNames());
83
84 // Detect missing fields
86 {
87 if
88 (
89 acceptFields[i].isLiteral()
90 && !allNames.found(acceptFields[i])
91 )
92 {
93 missed.append(i);
94 }
95 }
96
97 for (const word& fldName : allNames)
98 {
99 const auto iter = cloudObjects.cfind(fldName);
100 if (!pred(fldName) || !iter.found())
101 {
102 continue; // reject
103 }
104
105 const IOobject& io = *(iter.val());
106
107 if
108 (
109 //OR: fieldTypes::basic.found(io.headerClassName())
116 )
117 {
118 // Transfer from cloudObjects -> filteredObjects
119 filteredObjects.add(cloudObjects.remove(fldName));
120 }
121 }
122 }
123
124 if (missed.size())
125 {
127 << nl
128 << "Cannot find field file matching "
130 }
131
132 return filteredObjects;
133}
134
135
136void readFieldsAndWriteVTK
137(
138 OFstream& os,
139 const List<labelList>& particleMap,
140 const IOobjectList& filteredObjects
141)
142{
143 processFields<label>(os, particleMap, filteredObjects);
144 processFields<scalar>(os, particleMap, filteredObjects);
145 processFields<vector>(os, particleMap, filteredObjects);
146 processFields<sphericalTensor>(os, particleMap, filteredObjects);
147 processFields<symmTensor>(os, particleMap, filteredObjects);
148 processFields<tensor>(os, particleMap, filteredObjects);
149}
150
151} // End namespace Foam
152
153
154using namespace Foam;
155
156// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
157
158int main(int argc, char *argv[])
159{
160 argList::addNote
161 (
162 "Generate a legacy VTK file of particle tracks for cases that were"
163 " computed using a steady-state cloud"
164 );
165
166 argList::noParallel();
167 timeSelector::addOptions();
168 #include "addRegionOption.H"
169 argList::addOption
170 (
171 "dict",
172 "file",
173 "Alternative particleTrackDict dictionary"
174 );
175 argList::addVerboseOption("Additional verbosity");
176
177 #include "setRootCase.H"
178 #include "createTime.H"
179 instantList timeDirs = timeSelector::select0(runTime, args);
180 #include "createNamedMesh.H"
181 #include "createControls.H"
182
183 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
184
185 const fileName vtkPath(runTime.rootPath()/runTime.globalCaseName()/"VTK");
186 mkDir(vtkPath);
187
188 forAll(timeDirs, timeI)
189 {
190 runTime.setTime(timeDirs[timeI], timeI);
191 Info<< "Time = " << runTime.timeName() << endl;
192
193 const fileName vtkTimePath(vtkPath/runTime.timeName());
194 mkDir(vtkTimePath);
195
196 pointField particlePosition;
197 labelList particleToTrack;
198 label nTracks = 0;
199
200 // Transfer particles to (more convenient) list
201 {
202 Info<< " Reading particle positions" << endl;
203
205 Info<< "\n Read " << returnReduce(myCloud.size(), sumOp<label>())
206 << " particles" << endl;
207
208 const label nParticles = myCloud.size();
209
210 particlePosition.resize(nParticles);
211 particleToTrack.resize(nParticles);
212
213 LabelPairMap<label> trackTable;
214
215 label np = 0;
216 for (const passiveParticle& p : myCloud)
217 {
218 const label origId = p.origId();
219 const label origProc = p.origProc();
220 particlePosition[np] = p.position();
221
222 const labelPair key(origProc, origId);
223
224 const auto iter = trackTable.cfind(key);
225
226 if (iter.found())
227 {
228 particleToTrack[np] = *iter;
229 }
230 else
231 {
232 particleToTrack[np] = trackTable.size();
233 trackTable.insert(key, trackTable.size());
234 }
235
236 ++np;
237 }
238
239 nTracks = trackTable.size();
240 }
241
242 if (nTracks == 0)
243 {
244 Info<< "\n No track data" << endl;
245 }
246 else
247 {
248 Info<< "\n Generating " << nTracks << " tracks" << endl;
249
250 // Determine length of each track
251 labelList trackLengths(nTracks, Zero);
252 for (const label tracki : particleToTrack)
253 {
254 ++trackLengths[tracki];
255 }
256
257 // Particle "age" property used to sort the tracks
258 List<SortableList<scalar>> agePerTrack(nTracks);
259 List<List<label>> particleMap(nTracks);
260
261 forAll(trackLengths, i)
262 {
263 const label length = trackLengths[i];
264 agePerTrack[i].setSize(length);
265 particleMap[i].setSize(length);
266 }
267
268 // Store the particle age per track
269 IOobjectList cloudObjects
270 (
271 mesh,
272 runTime.timeName(),
273 cloud::prefix/cloudName
274 );
275
276 // TODO: gather age across all procs
277 {
278 tmp<IOField<scalar>> tage =
279 readParticleField<scalar>("age", cloudObjects);
280
281 const auto& age = tage();
282
283 labelList trackSamples(nTracks, Zero);
284
285 forAll(particleToTrack, i)
286 {
287 const label tracki = particleToTrack[i];
288 const label samplei = trackSamples[tracki];
289 agePerTrack[tracki][samplei] = age[i];
290 particleMap[tracki][samplei] = i;
291 ++trackSamples[tracki];
292 }
293 }
294
295
296 const IOobjectList filteredObjects
297 (
298 preFilterFields(cloudObjects, acceptFields, excludeFields)
299 );
300
301
302 if (Pstream::master())
303 {
304 OFstream os(vtkTimePath/"particleTracks.vtk");
305
306 Info<< "\n Writing particle tracks to " << os.name() << endl;
307
308 label nPoints = sum(trackLengths);
309
310 os << "# vtk DataFile Version 2.0" << nl
311 << "particleTracks" << nl
312 << "ASCII" << nl
313 << "DATASET POLYDATA" << nl
314 << "POINTS " << nPoints << " float" << nl;
315
316 Info<< "\n Writing points" << endl;
317
318 {
319 forAll(agePerTrack, i)
320 {
321 agePerTrack[i].sort();
322
323 const labelList& ids = agePerTrack[i].indices();
324 labelList& particleIds = particleMap[i];
325
326 {
327 // Update addressing
328 List<label> sortedIds(ids);
329 forAll(sortedIds, j)
330 {
331 sortedIds[j] = particleIds[ids[j]];
332 }
333 particleIds = sortedIds;
334 }
335
336 forAll(ids, j)
337 {
338 const label localId = particleIds[j];
339 const point& pos = particlePosition[localId];
340 os << pos.x() << ' ' << pos.y() << ' ' << pos.z()
341 << nl;
342 }
343 }
344 }
345
346
347 // Write track (line) connectivity to file
348
349 Info<< "\n Writing track lines" << endl;
350 os << "\nLINES " << nTracks << ' ' << nPoints + nTracks << nl;
351
352 // Write ids of track points to file
353 {
354 label globalPtI = 0;
355 forAll(particleMap, i)
356 {
357 os << particleMap[i].size() << nl;
358
359 forAll(particleMap[i], j)
360 {
361 os << ' ' << globalPtI++;
362
363 if (((j + 1) % 10 == 0) && (j != 0))
364 {
365 os << nl;
366 }
367 }
368
369 os << nl;
370 }
371 }
372
373
374 const label nFields = filteredObjects.size();
375
376 os << "POINT_DATA " << nPoints << nl
377 << "FIELD attributes " << nFields << nl;
378
379 Info<< "\n Processing fields" << nl << endl;
380
381 readFieldsAndWriteVTK(os, particleMap, filteredObjects);
382 }
383 }
384 Info<< endl;
385 }
386
387 Info<< "End\n" << endl;
388
389 return 0;
390}
391
392
393// ************************************************************************* //
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
A HashTable similar to std::unordered_map.
Definition: HashTable.H:123
const_iterator cfind(const Key &key) const
Find and return an const_iterator set at the hashed entry.
Definition: HashTableI.H:141
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
label capacity() const noexcept
The size of the underlying table.
Definition: HashTableI.H:45
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
A primitive field of type <T> with automated input and output.
Definition: IOField.H:58
List of IOobjects with searching and retrieving facilities.
Definition: IOobjectList.H:59
wordList sortedNames() const
The sorted names of the IOobjects.
Definition: IOobjectList.C:383
bool remove(const IOobject &io)
Remove IOobject from the list.
Definition: IOobjectList.C:259
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & headerClassName() const noexcept
Return name of the class name read from header.
Definition: IOobjectI.H:83
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
Output to file stream, using an OSstream.
Definition: OFstream.H:57
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
bool found(const T &val, label pos=0) const
True if the value if found in the list.
Definition: UListI.H:265
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A class for handling file names.
Definition: fileName.H:76
A Cloud of passive particles.
Copy of base particle.
A class for managing temporary objects.
Definition: tmp.H:65
A List of wordRe with additional matching capabilities.
Definition: wordRes.H:54
A class for handling words, derived from Foam::string.
Definition: word.H:68
volScalarField & p
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
OBJstream os(runTime.globalPath()/outputName)
label nPoints
A HashTable to objects of type <T> with a labelPair key. The hashing is based on labelPair (FixedList...
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
#define WarningInFunction
Report a warning using Foam::Warning.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
Definition: foamGltfBase.H:108
Namespace for OpenFOAM.
dimensionedScalar pos(const dimensionedScalar &ds)
List< word > wordList
A List of words.
Definition: fileName.H:63
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
static const char *const typeName
The type name used in ensight case files.
Functor wrapper of allow/deny lists of wordRe for filtering.
Definition: wordRes.H:174
wordRes excludeFields
wordRes acceptFields
const word cloudName(propsDict.get< word >("cloud"))
mkDir(pdfPath)