FreeFOAM The Cross-Platform CFD Toolkit
foamToFieldview9.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) 1991-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  foamToFieldview9
26 
27 Description
28  Write out the OpenFOAM mesh in Version 3.0 Fieldview-UNS format (binary).
29 
30  See Fieldview Release 9 Reference Manual - Appendix D
31  (Unstructured Data Format)
32  Borrows various from uns/write_binary_uns.c from FieldView dist.
33 
34 Usage
35 
36  - foamToFieldview9 [OPTIONS]
37 
38  @param -noWall \n
39  Do not export the walls.
40 
41  @param -noZero \n
42  Ignore timestep 0.
43 
44  @param -constant \n
45  Include the constant directory.
46 
47  @param -time <time>\n
48  Apply only to specific time.
49 
50  @param -latestTime \n
51  Only apply to latest time step.
52 
53  @param -case <dir>\n
54  Case directory.
55 
56  @param -help \n
57  Display help message.
58 
59  @param -doc \n
60  Display Doxygen API documentation page for this application.
61 
62  @param -srcDoc \n
63  Display Doxygen source documentation page for this application.
64 
65 \*---------------------------------------------------------------------------*/
66 
67 #include <OpenFOAM/argList.H>
68 #include <OpenFOAM/timeSelector.H>
69 #include <finiteVolume/volFields.H>
71 #include <OpenFOAM/pointFields.H>
72 #include <OpenFOAM/scalarIOField.H>
76 
77 #include <lagrangian/Cloud.H>
79 
80 #include <OpenFOAM/IOobjectList.H>
81 #include <OpenFOAM/boolList.H>
82 #include <OpenFOAM/stringList.H>
83 #include <OpenFOAM/cellModeller.H>
84 
85 #include <OpenFOAM/floatScalar.H>
86 #include "calcFaceAddressing.H"
87 #include "writeFunctions.H"
88 #include "fieldviewTopology.H"
89 
90 #include <fstream>
91 
92 #include "fv_reader_tags.h"
93 
94 extern "C"
95 {
96  unsigned int fv_encode_elem_header(int elem_type, int wall_info[]);
97 }
98 
99 using namespace Foam;
100 
101 typedef Field<floatScalar> floatField;
102 
103 
104 static HashTable<word> FieldviewNames;
105 
106 
107 static word getFieldViewName(const word& foamName)
108 {
109  if (FieldviewNames.found(foamName))
110  {
111  return FieldviewNames[foamName];
112  }
113  else
114  {
115  return foamName;
116  }
117 }
118 
119 
120 static void printNames(const wordList& names, Ostream& os)
121 {
122  forAll(names, fieldI)
123  {
124  Info<< " " << names[fieldI] << '/' << getFieldViewName(names[fieldI]);
125  }
126 }
127 
128 
129 // Count number of vertices used by celli
130 static label countVerts(const primitiveMesh& mesh, const label celli)
131 {
132  const cell& cll = mesh.cells()[celli];
133 
134  // Count number of vertices used
135  labelHashSet usedVerts(10*cll.size());
136 
137  forAll(cll, cellFacei)
138  {
139  const face& f = mesh.faces()[cll[cellFacei]];
140 
141  forAll(f, fp)
142  {
143  if (!usedVerts.found(f[fp]))
144  {
145  usedVerts.insert(f[fp]);
146  }
147  }
148  }
149  return usedVerts.toc().size();
150 }
151 
152 
153 static void writeFaceData
154 (
155  const polyMesh& mesh,
156  const fieldviewTopology& topo,
157  const label patchI,
158  const scalarField& patchField,
159  const bool writePolyFaces,
160  std::ofstream& fvFile
161 )
162 {
163  const polyPatch& pp = mesh.boundaryMesh()[patchI];
164 
165  // Start off with dummy value.
166  if (writePolyFaces)
167  {
168  floatField fField(topo.nPolyFaces()[patchI], 0.0);
169 
170  // Fill selected faces with field values
171  label polyFaceI = 0;
172  forAll(patchField, faceI)
173  {
174  if (pp[faceI].size() > 4)
175  {
176  fField[polyFaceI++] = float(patchField[faceI]);
177  }
178  }
179 
180  fvFile.write
181  (
182  reinterpret_cast<char*>(fField.begin()), fField.size()*sizeof(float)
183  );
184  }
185  else
186  {
187  floatField fField(pp.size() - topo.nPolyFaces()[patchI], 0.0);
188 
189  // Fill selected faces with field values
190  label quadFaceI = 0;
191  forAll(patchField, faceI)
192  {
193  if (pp[faceI].size() <= 4)
194  {
195  fField[quadFaceI++] = float(patchField[faceI]);
196  }
197  }
198 
199  fvFile.write
200  (
201  reinterpret_cast<char*>(fField.begin()), fField.size()*sizeof(float)
202  );
203  }
204 }
205 
206 
207 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
208 // Main program:
209 
210 int main(int argc, char *argv[])
211 {
213  argList::validOptions.insert("noWall", "");
214  timeSelector::addOptions(true, false);
215 
216 # include <OpenFOAM/setRootCase.H>
217 # include <OpenFOAM/createTime.H>
218 
220 
221 # include <OpenFOAM/createMesh.H>
222 
223  // Initialize name mapping table
224  FieldviewNames.insert("alpha", "aalpha");
225  FieldviewNames.insert("Alpha", "AAlpha");
226  FieldviewNames.insert("fsmach", "ffsmach");
227  FieldviewNames.insert("FSMach", "FFSMach");
228  FieldviewNames.insert("re", "rre");
229  FieldviewNames.insert("Re", "RRe");
230  FieldviewNames.insert("time", "ttime");
231  FieldviewNames.insert("Time", "TTime");
232  FieldviewNames.insert("pi", "ppi");
233  FieldviewNames.insert("PI", "PPI");
234  FieldviewNames.insert("x", "xx");
235  FieldviewNames.insert("X", "XX");
236  FieldviewNames.insert("y", "yy");
237  FieldviewNames.insert("Y", "YY");
238  FieldviewNames.insert("z", "zz");
239  FieldviewNames.insert("Z", "ZZ");
240  FieldviewNames.insert("rcyl", "rrcyl");
241  FieldviewNames.insert("Rcyl", "RRcyl");
242  FieldviewNames.insert("theta", "ttheta");
243  FieldviewNames.insert("Theta", "TTheta");
244  FieldviewNames.insert("rsphere", "rrsphere");
245  FieldviewNames.insert("Rsphere", "RRsphere");
246  FieldviewNames.insert("k", "kk");
247  FieldviewNames.insert("Kcond", "KKcond");
248 
249 
250  // Scan for all available fields, in all timesteps
251  // volScalarNames : all scalar fields
252  // volVectorNames : ,, vector ,,
253  // surfScalarNames : surface fields
254  // surfVectorNames : ,,
255  // sprayScalarNames: spray fields
256  // sprayVectorNames: ,,
257 # include "getFieldNames.H"
258 
259  bool hasLagrangian = false;
261  {
262  hasLagrangian = true;
263  }
264 
265  Info<< "All fields: Foam/Fieldview" << endl;
266  Info<< " volScalar :";
267  printNames(volScalarNames, Info);
268  Info<< endl;
269  Info<< " volVector :";
270  printNames(volVectorNames, Info);
271  Info<< endl;
272  Info<< " surfScalar :";
273  printNames(surfScalarNames, Info);
274  Info<< endl;
275  Info<< " surfVector :";
276  printNames(surfVectorNames, Info);
277  Info<< endl;
278  Info<< " sprayScalar :";
279  printNames(sprayScalarNames, Info);
280  Info<< endl;
281  Info<< " sprayVector :";
282  printNames(sprayVectorNames, Info);
283  Info<< endl;
284 
285 
286  //
287  // Start writing
288  //
289 
290  // make a directory called FieldView in the case
291  fileName fvPath(runTime.path()/"Fieldview");
292 
293  if (isDir(fvPath))
294  {
295  rmDir(fvPath);
296  }
297 
298  mkDir(fvPath);
299 
300  fileName fvParticleFileName(fvPath/runTime.caseName() + ".fvp");
301  if (hasLagrangian)
302  {
303  Info<< "Opening particle file " << fvParticleFileName << endl;
304  }
305  std::ofstream fvParticleFile(fvParticleFileName.c_str());
306 
307  // Write spray file header
308  if (hasLagrangian)
309  {
310 # include "writeSprayHeader.H"
311  }
312 
313  // Current mesh. Start off from unloaded mesh.
314  autoPtr<fieldviewTopology> topoPtr(NULL);
315 
316  label fieldViewTime = 0;
317 
318  forAll(timeDirs, timeI)
319  {
320  runTime.setTime(timeDirs[timeI], timeI);
321  Info<< "Time: " << runTime.timeName() << endl;
322 
323  fvMesh::readUpdateState state = mesh.readUpdate();
324 
325  if
326  (
327  timeI == 0
328  || state == fvMesh::TOPO_CHANGE
329  || state == fvMesh::TOPO_PATCH_CHANGE
330  )
331  {
332  // Mesh topo changed. Update Fieldview topo.
333 
334  topoPtr.reset
335  (
337  (
338  mesh,
339  !args.optionFound("noWall")
340  )
341  );
342 
343  Info<< " Mesh read:" << endl
344  << " tet : " << topoPtr().nTet() << endl
345  << " hex : " << topoPtr().nHex() << endl
346  << " prism : " << topoPtr().nPrism() << endl
347  << " pyr : " << topoPtr().nPyr() << endl
348  << " poly : " << topoPtr().nPoly() << endl
349  << endl;
350  }
351  else if (state == fvMesh::POINTS_MOVED)
352  {
353  // points exists for time step, let's read them
354  Info<< " Points file detected - updating points" << endl;
355  }
356 
357  const fieldviewTopology& topo = topoPtr();
358 
359 
360  //
361  // Create file and write header
362  //
363 
364  fileName fvFileName
365  (
366  fvPath/runTime.caseName() + "_" + Foam::name(timeI) + ".uns"
367  );
368 
369  Info<< " file:" << fvFileName.c_str() << endl;
370 
371 
372  std::ofstream fvFile(fvFileName.c_str());
373 
374  //Info<< "Writing header ..." << endl;
375 
376  // Output the magic number.
377  writeInt(fvFile, FV_MAGIC);
378 
379  // Output file header and version number.
380  writeStr80(fvFile, "FIELDVIEW");
381 
382  // This version of the FIELDVIEW unstructured file is "3.0".
383  // This is written as two integers.
384  writeInt(fvFile, 3);
385  writeInt(fvFile, 0);
386 
387 
388  // File type code. Grid and results.
389  writeInt(fvFile, FV_COMBINED_FILE);
390 
391  // Reserved field, always zero
392  writeInt(fvFile, 0);
393 
394  // Output constants for time, fsmach, alpha and re.
395  float fBuf[4];
396  fBuf[0] = runTime.value();
397  fBuf[1] = 0.0;
398  fBuf[2] = 0.0;
399  fBuf[3] = 1.0;
400  fvFile.write(reinterpret_cast<char*>(fBuf), 4*sizeof(float));
401 
402 
403  // Output the number of grids
404  writeInt(fvFile, 1);
405 
406 
407  //
408  // Boundary table
409  //
410  //Info<< "Writing boundary table ..." << endl;
411 
412  // num patches
413  writeInt(fvFile, mesh.boundary().size());
414 
415  forAll (mesh.boundary(), patchI)
416  {
417  const fvPatch& currPatch = mesh.boundary()[patchI];
418 
419  writeInt(fvFile, 1); // data present
420  writeInt(fvFile, 1); // normals ok
421 
422  // name
423  writeStr80(fvFile, currPatch.name().c_str());
424  }
425 
426 
427  //
428  // Create fields:
429  // volFieldPtrs : List of ptrs to all volScalar/Vector fields
430  // (null if field not present at current time)
431  // volFieldNames : FieldView compatible names of volFields
432  // surfFieldPtrs : same for surfaceScalar/Vector
433  // surfFieldNames
434 # include "createFields.H"
435 
436 
437 
438  //
439  // Write Variables table
440  //
441 
442  //Info<< "Writing variables table ..." << endl;
443 
444  writeInt(fvFile, volFieldNames.size());
445  forAll(volFieldNames, fieldI)
446  {
447  if (volFieldPtrs[fieldI] == NULL)
448  {
449  Info<< " dummy field for "
450  << volFieldNames[fieldI].c_str() << endl;
451  }
452 
453  writeStr80(fvFile, volFieldNames[fieldI].c_str());
454  }
455 
456  //
457  // Write Boundary Variables table = vol + surface fields
458  //
459 
460  //Info<< "Writing boundary variables table ..." << endl;
461 
462  writeInt
463  (
464  fvFile,
466  );
467  forAll(volFieldNames, fieldI)
468  {
469  writeStr80(fvFile, volFieldNames[fieldI].c_str());
470  }
471  forAll(surfFieldNames, fieldI)
472  {
473  if (surfFieldPtrs[fieldI] == NULL)
474  {
475  Info<< " dummy surface field for "
476  << surfFieldNames[fieldI].c_str() << endl;
477  }
478 
479  writeStr80(fvFile, surfFieldNames[fieldI].c_str());
480  }
481 
482 
483  // Output grid data.
484 
485  //
486  // Nodes
487  //
488 
489  //Info<< "Writing points ..." << endl;
490 
491  const pointField& points = mesh.points();
492  label nPoints = points.size();
493 
494  writeInt(fvFile, FV_NODES);
495  writeInt(fvFile, nPoints);
496 
497  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
498  {
499  floatField fField(nPoints);
500 
501  for (label pointi = 0; pointi<nPoints; pointi++)
502  {
503  fField[pointi] = float(points[pointi][cmpt]);
504  }
505 
506  fvFile.write
507  (
508  reinterpret_cast<char*>(fField.begin()),
509  fField.size()*sizeof(float)
510  );
511  }
512 
513  //
514  // Boundary Faces - regular
515  //
516 
517  //Info<< "Writing regular boundary faces ..." << endl;
518 
519  forAll(mesh.boundary(), patchI)
520  {
521  label nQuadFaces = topo.quadFaceLabels()[patchI].size()/4;
522 
523  if (nQuadFaces != 0)
524  {
525  writeInt(fvFile, FV_FACES);
526  writeInt(fvFile, patchI + 1); // patch number
527  writeInt(fvFile, nQuadFaces); // number of faces in patch
528  fvFile.write
529  (
530  reinterpret_cast<const char*>
531  (topo.quadFaceLabels()[patchI].begin()),
532  nQuadFaces*4*sizeof(int)
533  );
534  }
535  }
536 
537  //
538  // Boundary Faces - arbitrary polygon
539  //
540 
541  //Info<< "Write polygonal boundary faces ..." << endl;
542 
543  forAll(mesh.boundary(), patchI)
544  {
545  if (topo.nPolyFaces()[patchI] > 0)
546  {
547  writeInt(fvFile, FV_ARB_POLY_FACES);
548  writeInt(fvFile, patchI + 1);
549 
550  // number of arb faces in patch
551  writeInt(fvFile, topo.nPolyFaces()[patchI]);
552 
553  const polyPatch& patchFaces = mesh.boundary()[patchI].patch();
554 
555  forAll(patchFaces, faceI)
556  {
557  const face& f = patchFaces[faceI];
558 
559  if (f.size() > 4)
560  {
561  writeInt(fvFile, f.size());
562 
563  forAll(f, fp)
564  {
565  writeInt(fvFile, f[fp] + 1);
566  }
567  }
568  }
569  }
570  }
571 
572 
573  //
574  // Write regular topology
575  //
576 
577  //Info<< "Writing regular elements ..." << endl;
578 
579  writeInt(fvFile, FV_ELEMENTS);
580  writeInt(fvFile, topo.nTet());
581  writeInt(fvFile, topo.nHex());
582  writeInt(fvFile, topo.nPrism());
583  writeInt(fvFile, topo.nPyr());
584  fvFile.write
585  (
586  reinterpret_cast<const char*>(topo.tetLabels().begin()),
587  topo.nTet()*(1+4)*sizeof(int)
588  );
589  fvFile.write
590  (
591  reinterpret_cast<const char*>(topo.hexLabels().begin()),
592  topo.nHex()*(1+8)*sizeof(int)
593  );
594  fvFile.write
595  (
596  reinterpret_cast<const char*>(topo.prismLabels().begin()),
597  topo.nPrism()*(1+6)*sizeof(int)
598  );
599  fvFile.write
600  (
601  reinterpret_cast<const char*>(topo.pyrLabels().begin()),
602  topo.nPyr()*(1+5)*sizeof(int)
603  );
604 
605 
606  //
607  // Write arbitrary polyhedra
608  //
609 
610  //Info<< "Writing polyhedral elements ..." << endl;
611 
612 
613  const cellShapeList& cellShapes = mesh.cellShapes();
614  const cellModel& unknown = *(cellModeller::lookup("unknown"));
615 
616  if (topo.nPoly() > 0)
617  {
618  writeInt(fvFile, FV_ARB_POLY_ELEMENTS);
619  writeInt(fvFile, topo.nPoly());
620 
621  forAll(cellShapes, celli)
622  {
623  if (cellShapes[celli].model() == unknown)
624  {
625  const cell& cll = mesh.cells()[celli];
626 
627  // number of faces
628  writeInt(fvFile, cll.size());
629  // number of vertices used (no cell centre)
630  writeInt(fvFile, countVerts(mesh, celli));
631  // cell centre node id
632  writeInt(fvFile, -1);
633 
634  forAll(cll, cellFacei)
635  {
636  label faceI = cll[cellFacei];
637 
638  const face& f = mesh.faces()[faceI];
639 
640  // Not a wall for now
641  writeInt(fvFile, NOT_A_WALL);
642 
643  writeInt(fvFile, f.size());
644 
645  if (mesh.faceOwner()[faceI] == celli)
646  {
647  forAll(f, fp)
648  {
649  writeInt(fvFile, f[fp]+1);
650  }
651  }
652  else
653  {
654  for (label fp = f.size()-1; fp >= 0; fp--)
655  {
656  writeInt(fvFile, f[fp]+1);
657  }
658  }
659 
660  // Number of hanging nodes
661  writeInt(fvFile, 0);
662  }
663  }
664  }
665  }
666 
667 
668  //
669  // Variables data
670  //
671 
672  //Info<< "Writing variables data ..." << endl;
673 
674  volPointInterpolation pInterp(mesh);
675 
676  writeInt(fvFile, FV_VARIABLES);
677 
678 
679  forAll(volFieldPtrs, fieldI)
680  {
681  if (volFieldPtrs[fieldI] != NULL)
682  {
683  const volScalarField& vField = *volFieldPtrs[fieldI];
684 
685  // Interpolate to points
686  pointScalarField psf(pInterp.interpolate(vField));
687 
688  floatField fField(nPoints);
689 
690  for (label pointi = 0; pointi<nPoints; pointi++)
691  {
692  fField[pointi] = float(psf[pointi]);
693  }
694 
695  fvFile.write
696  (
697  reinterpret_cast<char*>(fField.begin()),
698  fField.size()*sizeof(float)
699  );
700  }
701  else
702  {
703  // Create dummy field
704  floatField dummyField(nPoints, 0.0);
705 
706  fvFile.write
707  (
708  reinterpret_cast<char*>(dummyField.begin()),
709  dummyField.size()*sizeof(float)
710  );
711  }
712  }
713 
714 
715  //
716  // Boundary variables data
717  // 1. volFields
718  // 2. surfFields
719 
720  //Info<< "Writing regular boundary elements data ..." << endl;
721 
722  writeInt(fvFile, FV_BNDRY_VARS);
723 
724  forAll(volFieldPtrs, fieldI)
725  {
726  if (volFieldPtrs[fieldI] != NULL)
727  {
728  const volScalarField& vsf = *volFieldPtrs[fieldI];
729 
730  forAll (mesh.boundary(), patchI)
731  {
732  writeFaceData
733  (
734  mesh,
735  topo,
736  patchI,
737  vsf.boundaryField()[patchI],
738  false,
739  fvFile
740  );
741  }
742  }
743  else
744  {
745  forAll (mesh.boundaryMesh(), patchI)
746  {
747  // Dummy value.
748  floatField fField
749  (
750  mesh.boundaryMesh()[patchI].size()
751  - topo.nPolyFaces()[patchI],
752  0.0
753  );
754 
755  fvFile.write
756  (
757  reinterpret_cast<char*>(fField.begin()),
758  fField.size()*sizeof(float)
759  );
760  }
761  }
762  }
763 
764  // surfFields
765  forAll(surfFieldPtrs, fieldI)
766  {
767  if (surfFieldPtrs[fieldI] != NULL)
768  {
769  const surfaceScalarField& ssf = *surfFieldPtrs[fieldI];
770 
771  forAll (mesh.boundary(), patchI)
772  {
773  writeFaceData
774  (
775  mesh,
776  topo,
777  patchI,
778  ssf.boundaryField()[patchI],
779  false,
780  fvFile
781  );
782  }
783  }
784  else
785  {
786  forAll (mesh.boundaryMesh(), patchI)
787  {
788  // Dummy value.
789  floatField fField
790  (
791  mesh.boundaryMesh()[patchI].size()
792  - topo.nPolyFaces()[patchI],
793  0.0
794  );
795 
796  fvFile.write
797  (
798  reinterpret_cast<char*>(fField.begin()),
799  fField.size()*sizeof(float)
800  );
801  }
802  }
803  }
804 
805  //
806  // Polygonal faces boundary data
807  // 1. volFields
808  // 2. surfFields
809 
810  //Info<< "Writing polygonal boundary elements data ..." << endl;
811 
812  writeInt(fvFile, FV_ARB_POLY_BNDRY_VARS);
813  forAll(volFieldPtrs, fieldI)
814  {
815  if (volFieldPtrs[fieldI] != NULL)
816  {
817  const volScalarField& vsf = *volFieldPtrs[fieldI];
818 
819  // All non-empty patches
820  forAll (mesh.boundary(), patchI)
821  {
822  writeFaceData
823  (
824  mesh,
825  topo,
826  patchI,
827  vsf.boundaryField()[patchI],
828  true,
829  fvFile
830  );
831  }
832  }
833  else
834  {
835  forAll (mesh.boundary(), patchI)
836  {
837  // Dummy value.
838  floatField fField(topo.nPolyFaces()[patchI], 0.0);
839 
840  fvFile.write
841  (
842  reinterpret_cast<char*>(fField.begin()),
843  fField.size()*sizeof(float)
844  );
845  }
846  }
847  }
848 
849  // surfFields
850  forAll(surfFieldPtrs, fieldI)
851  {
852  if (surfFieldPtrs[fieldI] != NULL)
853  {
854  const surfaceScalarField& ssf = *surfFieldPtrs[fieldI];
855 
856  // All non-empty patches
857  forAll(mesh.boundary(), patchI)
858  {
859  writeFaceData
860  (
861  mesh,
862  topo,
863  patchI,
864  ssf.boundaryField()[patchI],
865  true,
866  fvFile
867  );
868  }
869  }
870  else
871  {
872  forAll (mesh.boundaryMesh(), patchI)
873  {
874  // Dummy value.
875  floatField fField
876  (
877  mesh.boundaryMesh()[patchI].size()
878  - topo.nPolyFaces()[patchI],
879  0.0
880  );
881 
882  fvFile.write
883  (
884  reinterpret_cast<char*>(fField.begin()),
885  fField.size()*sizeof(float)
886  );
887  }
888  }
889  }
890 
891 
892  //
893  // Cleanup volume and surface fields
894  //
895  forAll(volFieldPtrs, fieldI)
896  {
897  delete volFieldPtrs[fieldI];
898  }
899  forAll(surfFieldPtrs, fieldI)
900  {
901  delete surfFieldPtrs[fieldI];
902  }
903 
904 
905 
906 
907  //
908  // Spray
909  //
910  if (hasLagrangian)
911  {
912  // Read/create fields:
913  // sprayScalarFieldPtrs: List of ptrs to lagrangian scalfields
914  // sprayVectorFieldPtrs: ,, vec ,,
915 # include "createSprayFields.H"
916 
917 
918  // Write time header
919 
920  // Time index (FieldView: has to start from 1)
921  writeInt(fvParticleFile, fieldViewTime + 1);
922 
923  // Time value
924  writeFloat(fvParticleFile, runTime.value());
925 
926  // Read particles
927  Cloud<passiveParticle> parcels(mesh);
928 
929  // Num particles
930  writeInt(fvParticleFile, parcels.size());
931 
932  Info<< " Writing " << parcels.size() << " particles." << endl;
933 
934 
935  //
936  // Output data parcelwise
937  //
938 
939  label parcelNo = 0;
940 
941 
942  for
943  (
944  Cloud<passiveParticle>::iterator elmnt = parcels.begin();
945  elmnt != parcels.end();
946  ++elmnt, parcelNo++
947  )
948  {
949  writeInt(fvParticleFile, parcelNo+1);
950 
951  writeFloat(fvParticleFile, elmnt().position().x());
952  writeFloat(fvParticleFile, elmnt().position().y());
953  writeFloat(fvParticleFile, elmnt().position().z());
954 
956  {
957  if (sprayScalarFieldPtrs[fieldI] != NULL)
958  {
959  const IOField<scalar>& sprayField =
960  *sprayScalarFieldPtrs[fieldI];
961  writeFloat
962  (
963  fvParticleFile,
964  sprayField[parcelNo]
965  );
966  }
967  else
968  {
969  writeFloat(fvParticleFile, 0.0);
970  }
971  }
973  {
974  if (sprayVectorFieldPtrs[fieldI] != NULL)
975  {
976  const IOField<vector>& sprayVectorField =
977  *sprayVectorFieldPtrs[fieldI];
978  const vector& val =
979  sprayVectorField[parcelNo];
980 
981  writeFloat(fvParticleFile, val.x());
982  writeFloat(fvParticleFile, val.y());
983  writeFloat(fvParticleFile, val.z());
984  }
985  else
986  {
987  writeFloat(fvParticleFile, 0.0);
988  writeFloat(fvParticleFile, 0.0);
989  writeFloat(fvParticleFile, 0.0);
990  }
991  }
992  }
993 
994  // increment fieldView particle time
995  fieldViewTime++;
996 
997 
998  //
999  // Cleanup spray fields
1000  //
1001  forAll(sprayScalarFieldPtrs, fieldI)
1002  {
1003  delete sprayScalarFieldPtrs[fieldI];
1004  }
1005  forAll(sprayVectorFieldPtrs, fieldI)
1006  {
1007  delete sprayVectorFieldPtrs[fieldI];
1008  }
1009 
1010  } // end of hasLagrangian
1011  }
1012 
1013  if (!hasLagrangian)
1014  {
1015  rm(fvParticleFileName);
1016  }
1017 
1018  Info << "End\n" << endl;
1019 
1020  return 0;
1021 }
1022 
1023 
1024 // ************************ vim: set sw=4 sts=4 et: ************************ //