applications/utilities/postProcessing/dataConversion/foamToFieldview9/foamToFieldview9.C

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