applications/utilities/postProcessing/dataConversion/foamToFieldview9/foamToFieldview9.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 Description
  25     Write out the OpenFOAM 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 "Cloud.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::addBoolOption
 180     (
 181         "noWall",
 182         "skip setting wall information"
 183     );
 184     timeSelector::addOptions(true, false);
 185
 186 #   include "addRegionOption.H"
 187
 188 #   include "setRootCase.H"
 189 #   include "createTime.H"
 190
 191     instantList timeDirs = timeSelector::select0(runTime, args);
 192
 193 #   include "createNamedMesh.H"
 194
 195     // Initialize name mapping table
 196     FieldviewNames.insert("alpha", "aalpha");
 197     FieldviewNames.insert("Alpha", "AAlpha");
 198     FieldviewNames.insert("fsmach", "ffsmach");
 199     FieldviewNames.insert("FSMach", "FFSMach");
 200     FieldviewNames.insert("re", "rre");
 201     FieldviewNames.insert("Re", "RRe");
 202     FieldviewNames.insert("time", "ttime");
 203     FieldviewNames.insert("Time", "TTime");
 204     FieldviewNames.insert("pi", "ppi");
 205     FieldviewNames.insert("PI", "PPI");
 206     FieldviewNames.insert("x", "xx");
 207     FieldviewNames.insert("X", "XX");
 208     FieldviewNames.insert("y", "yy");
 209     FieldviewNames.insert("Y", "YY");
 210     FieldviewNames.insert("z", "zz");
 211     FieldviewNames.insert("Z", "ZZ");
 212     FieldviewNames.insert("rcyl", "rrcyl");
 213     FieldviewNames.insert("Rcyl", "RRcyl");
 214     FieldviewNames.insert("theta", "ttheta");
 215     FieldviewNames.insert("Theta", "TTheta");
 216     FieldviewNames.insert("rsphere", "rrsphere");
 217     FieldviewNames.insert("Rsphere", "RRsphere");
 218     FieldviewNames.insert("k", "kk");
 219     FieldviewNames.insert("K", "KK");
 220
 221
 222     // Scan for all available fields, in all timesteps
 223     //     volScalarNames  : all scalar fields
 224     //     volVectorNames  : ,,  vector ,,
 225     //     surfScalarNames : surface fields
 226     //     surfVectorNames : ,,
 227     //     sprayScalarNames: spray fields
 228     //     sprayVectorNames: ,,
 229 #   include "getFieldNames.H"
 230
 231     bool hasLagrangian = false;
 232     if (sprayScalarNames.size() || sprayVectorNames.size())
 233     {
 234         hasLagrangian = true;
 235     }
 236
 237     Info<< "All fields:   Foam/Fieldview" << endl;
 238     Info<< "    volScalar   :";
 239     printNames(volScalarNames, Info);
 240     Info<< endl;
 241     Info<< "    volVector   :";
 242     printNames(volVectorNames, Info);
 243     Info<< endl;
 244     Info<< "    surfScalar  :";
 245     printNames(surfScalarNames, Info);
 246     Info<< endl;
 247     Info<< "    surfVector  :";
 248     printNames(surfVectorNames, Info);
 249     Info<< endl;
 250     Info<< "    sprayScalar :";
 251     printNames(sprayScalarNames, Info);
 252     Info<< endl;
 253     Info<< "    sprayVector :";
 254     printNames(sprayVectorNames, Info);
 255     Info<< endl;
 256
 257
 258     //
 259     // Start writing
 260     //
 261
 262     // make a directory called FieldView in the case
 263     fileName fvPath(runTime.path()/"Fieldview");
 264
 265     if (regionName != polyMesh::defaultRegion)
 266     {
 267         fvPath = fvPath/regionName;
 268     }
 269
 270     if (isDir(fvPath))
 271     {
 272         if (regionName != polyMesh::defaultRegion)
 273         {
 274             Info<< "Keeping old FieldView files in " << fvPath << nl << endl;
 275         }
 276         else
 277         {
 278             Info<< "Deleting old FieldView files in " << fvPath << nl << endl;
 279             rmDir(fvPath);
 280         }
 281     }
 282
 283     mkDir(fvPath);
 284
 285     fileName fvParticleFileName(fvPath/runTime.caseName() + ".fvp");
 286     if (hasLagrangian)
 287     {
 288         Info<< "Opening particle file " << fvParticleFileName << endl;
 289     }
 290     std::ofstream fvParticleFile(fvParticleFileName.c_str());
 291
 292     // Write spray file header
 293     if (hasLagrangian)
 294     {
 295 #       include "writeSprayHeader.H"
 296     }
 297
 298     // Current mesh. Start off from unloaded mesh.
 299     autoPtr<fieldviewTopology> topoPtr(NULL);
 300
 301     label fieldViewTime = 0;
 302
 303     forAll(timeDirs, timeI)
 304     {
 305         runTime.setTime(timeDirs[timeI], timeI);
 306         Info<< "Time: " << runTime.timeName() << endl;
 307
 308         fvMesh::readUpdateState state = mesh.readUpdate();
 309
 310         if
 311         (
 312             timeI == 0
 313          || state == fvMesh::TOPO_CHANGE
 314          || state == fvMesh::TOPO_PATCH_CHANGE
 315         )
 316         {
 317             // Mesh topo changed. Update Fieldview topo.
 318
 319             topoPtr.reset
 320             (
 321                 new fieldviewTopology
 322                 (
 323                     mesh,
 324                     !args.optionFound("noWall")
 325                 )
 326             );
 327
 328             Info<< "    Mesh read:" << endl
 329                 << "        tet   : " << topoPtr().nTet() << endl
 330                 << "        hex   : " << topoPtr().nHex() << endl
 331                 << "        prism : " << topoPtr().nPrism() << endl
 332                 << "        pyr   : " << topoPtr().nPyr() << endl
 333                 << "        poly  : " << topoPtr().nPoly() << endl
 334                 << endl;
 335         }
 336         else if (state == fvMesh::POINTS_MOVED)
 337         {
 338             // points exists for time step, let's read them
 339             Info<< "    Points file detected - updating points" << endl;
 340         }
 341
 342         const fieldviewTopology& topo = topoPtr();
 343
 344
 345         //
 346         // Create file and write header
 347         //
 348
 349         fileName fvFileName
 350         (
 351             fvPath/runTime.caseName() + "_" + Foam::name(timeI) + ".uns"
 352         );
 353
 354         Info<< "    file:" << fvFileName.c_str() << endl;
 355
 356
 357         std::ofstream fvFile(fvFileName.c_str());
 358
 359         //Info<< "Writing header ..." << endl;
 360
 361         // Output the magic number.
 362         writeInt(fvFile, FV_MAGIC);
 363
 364         // Output file header and version number.
 365         writeStr80(fvFile, "FIELDVIEW");
 366
 367         // This version of the FIELDVIEW unstructured file is "3.0".
 368         // This is written as two integers.
 369         writeInt(fvFile, 3);
 370         writeInt(fvFile, 0);
 371
 372
 373         // File type code. Grid and results.
 374         writeInt(fvFile, FV_COMBINED_FILE);
 375
 376         // Reserved field, always zero
 377         writeInt(fvFile, 0);
 378
 379         // Output constants for time, fsmach, alpha and re.
 380         float fBuf[4];
 381         fBuf[0] = runTime.value();
 382         fBuf[1] = 0.0;
 383         fBuf[2] = 0.0;
 384         fBuf[3] = 1.0;
 385         fvFile.write(reinterpret_cast<char*>(fBuf), 4*sizeof(float));
 386
 387
 388         // Output the number of grids
 389         writeInt(fvFile, 1);
 390
 391
 392         //
 393         //  Boundary table
 394         //
 395         //Info<< "Writing boundary table ..." << endl;
 396
 397         // num patches
 398         writeInt(fvFile, mesh.boundary().size());
 399
 400         forAll(mesh.boundary(), patchI)
 401         {
 402             const fvPatch& currPatch = mesh.boundary()[patchI];
 403
 404             writeInt(fvFile, 1);   // data present
 405             writeInt(fvFile, 1);   // normals ok
 406
 407             // name
 408             writeStr80(fvFile, currPatch.name().c_str());
 409         }
 410
 411
 412         //
 413         // Create fields:
 414         //     volFieldPtrs     : List of ptrs to all volScalar/Vector fields
 415         //                        (null if field not present at current time)
 416         //     volFieldNames    : FieldView compatible names of volFields
 417         //     surfFieldPtrs    : same for surfaceScalar/Vector
 418         //     surfFieldNames
 419 #       include "createFields.H"
 420
 421
 422
 423         //
 424         // Write Variables table
 425         //
 426
 427         //Info<< "Writing variables table ..." << endl;
 428
 429         writeInt(fvFile, volFieldNames.size());
 430         forAll(volFieldNames, fieldI)
 431         {
 432             if (volFieldPtrs[fieldI] == NULL)
 433             {
 434                 Info<< "    dummy field for "
 435                     << volFieldNames[fieldI].c_str() << endl;
 436             }
 437
 438             writeStr80(fvFile, volFieldNames[fieldI].c_str());
 439         }
 440
 441         //
 442         // Write Boundary Variables table = vol + surface fields
 443         //
 444
 445         //Info<< "Writing boundary variables table ..." << endl;
 446
 447         writeInt
 448         (
 449             fvFile,
 450             volFieldNames.size() + surfFieldNames.size()
 451         );
 452         forAll(volFieldNames, fieldI)
 453         {
 454             writeStr80(fvFile, volFieldNames[fieldI].c_str());
 455         }
 456         forAll(surfFieldNames, fieldI)
 457         {
 458             if (surfFieldPtrs[fieldI] == NULL)
 459             {
 460                 Info<< "    dummy surface field for "
 461                     << surfFieldNames[fieldI].c_str() << endl;
 462             }
 463
 464             writeStr80(fvFile, surfFieldNames[fieldI].c_str());
 465         }
 466
 467
 468         // Output grid data.
 469
 470         //
 471         // Nodes
 472         //
 473
 474         //Info<< "Writing points ..." << endl;
 475
 476         const pointField& points = mesh.points();
 477         label nPoints = points.size();
 478
 479         writeInt(fvFile, FV_NODES);
 480         writeInt(fvFile, nPoints);
 481
 482         for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
 483         {
 484             floatField fField(nPoints);
 485
 486             for (label pointi = 0; pointi<nPoints; pointi++)
 487             {
 488                 fField[pointi] = float(points[pointi][cmpt]);
 489             }
 490
 491             fvFile.write
 492             (
 493                 reinterpret_cast<char*>(fField.begin()),
 494                 fField.size()*sizeof(float)
 495             );
 496         }
 497
 498         //
 499         // Boundary Faces - regular
 500         //
 501
 502         //Info<< "Writing regular boundary faces ..." << endl;
 503
 504         forAll(mesh.boundary(), patchI)
 505         {
 506             label nQuadFaces = topo.quadFaceLabels()[patchI].size()/4;
 507
 508             if (nQuadFaces != 0)
 509             {
 510                 writeInt(fvFile, FV_FACES);
 511                 writeInt(fvFile, patchI + 1);  // patch number
 512                 writeInt(fvFile, nQuadFaces);  // number of faces in patch
 513                 fvFile.write
 514                 (
 515                     reinterpret_cast<const char*>
 516                         (topo.quadFaceLabels()[patchI].begin()),
 517                     nQuadFaces*4*sizeof(int)
 518                 );
 519             }
 520         }
 521
 522         //
 523         // Boundary Faces - arbitrary polygon
 524         //
 525
 526         //Info<< "Write polygonal boundary faces ..." << endl;
 527
 528         forAll(mesh.boundary(), patchI)
 529         {
 530             if (topo.nPolyFaces()[patchI] > 0)
 531             {
 532                 writeInt(fvFile, FV_ARB_POLY_FACES);
 533                 writeInt(fvFile, patchI + 1);
 534
 535                 // number of arb faces in patch
 536                 writeInt(fvFile, topo.nPolyFaces()[patchI]);
 537
 538                 const polyPatch& patchFaces = mesh.boundary()[patchI].patch();
 539
 540                 forAll(patchFaces, faceI)
 541                 {
 542                     const face& f = patchFaces[faceI];
 543
 544                     if (f.size() > 4)
 545                     {
 546                         writeInt(fvFile, f.size());
 547
 548                         forAll(f, fp)
 549                         {
 550                             writeInt(fvFile, f[fp] + 1);
 551                         }
 552                     }
 553                 }
 554             }
 555         }
 556
 557
 558         //
 559         // Write regular topology
 560         //
 561
 562         //Info<< "Writing regular elements ..." << endl;
 563
 564         writeInt(fvFile, FV_ELEMENTS);
 565         writeInt(fvFile, topo.nTet());
 566         writeInt(fvFile, topo.nHex());
 567         writeInt(fvFile, topo.nPrism());
 568         writeInt(fvFile, topo.nPyr());
 569         fvFile.write
 570         (
 571             reinterpret_cast<const char*>(topo.tetLabels().begin()),
 572             topo.nTet()*(1+4)*sizeof(int)
 573         );
 574         fvFile.write
 575         (
 576             reinterpret_cast<const char*>(topo.hexLabels().begin()),
 577             topo.nHex()*(1+8)*sizeof(int)
 578         );
 579         fvFile.write
 580         (
 581             reinterpret_cast<const char*>(topo.prismLabels().begin()),
 582             topo.nPrism()*(1+6)*sizeof(int)
 583         );
 584         fvFile.write
 585         (
 586             reinterpret_cast<const char*>(topo.pyrLabels().begin()),
 587             topo.nPyr()*(1+5)*sizeof(int)
 588         );
 589
 590
 591         //
 592         // Write arbitrary polyhedra
 593         //
 594
 595         //Info<< "Writing polyhedral elements ..." << endl;
 596
 597
 598         const cellShapeList& cellShapes = mesh.cellShapes();
 599         const cellModel& unknown = *(cellModeller::lookup("unknown"));
 600
 601         if (topo.nPoly() > 0)
 602         {
 603             writeInt(fvFile, FV_ARB_POLY_ELEMENTS);
 604             writeInt(fvFile, topo.nPoly());
 605
 606             forAll(cellShapes, celli)
 607             {
 608                 if (cellShapes[celli].model() == unknown)
 609                 {
 610                     const cell& cll = mesh.cells()[celli];
 611
 612                     // number of faces
 613                     writeInt(fvFile, cll.size());
 614                     // number of vertices used (no cell centre)
 615                     writeInt(fvFile, countVerts(mesh, celli));
 616                     // cell centre node id
 617                     writeInt(fvFile, -1);
 618
 619                     forAll(cll, cellFacei)
 620                     {
 621                         label faceI = cll[cellFacei];
 622
 623                         const face& f = mesh.faces()[faceI];
 624
 625                         // Not a wall for now
 626                         writeInt(fvFile, NOT_A_WALL);
 627
 628                         writeInt(fvFile, f.size());
 629
 630                         if (mesh.faceOwner()[faceI] == celli)
 631                         {
 632                             forAll(f, fp)
 633                             {
 634                                 writeInt(fvFile, f[fp]+1);
 635                             }
 636                         }
 637                         else
 638                         {
 639                             for (label fp = f.size()-1; fp >= 0; fp--)
 640                             {
 641                                 writeInt(fvFile, f[fp]+1);
 642                             }
 643                         }
 644
 645                         // Number of hanging nodes
 646                         writeInt(fvFile, 0);
 647                     }
 648                 }
 649             }
 650         }
 651
 652
 653         //
 654         // Variables data
 655         //
 656
 657         //Info<< "Writing variables data ..." << endl;
 658
 659         volPointInterpolation pInterp(mesh);
 660
 661         writeInt(fvFile, FV_VARIABLES);
 662
 663
 664         forAll(volFieldPtrs, fieldI)
 665         {
 666             if (volFieldPtrs[fieldI] != NULL)
 667             {
 668                 const volScalarField& vField = *volFieldPtrs[fieldI];
 669
 670                 // Interpolate to points
 671                 pointScalarField psf(pInterp.interpolate(vField));
 672
 673                 floatField fField(nPoints);
 674
 675                 for (label pointi = 0; pointi<nPoints; pointi++)
 676                 {
 677                     fField[pointi] = float(psf[pointi]);
 678                 }
 679
 680                 fvFile.write
 681                 (
 682                     reinterpret_cast<char*>(fField.begin()),
 683                     fField.size()*sizeof(float)
 684                 );
 685             }
 686             else
 687             {
 688                 // Create dummy field
 689                 floatField dummyField(nPoints, 0.0);
 690
 691                 fvFile.write
 692                 (
 693                     reinterpret_cast<char*>(dummyField.begin()),
 694                     dummyField.size()*sizeof(float)
 695                 );
 696             }
 697         }
 698
 699
 700         //
 701         // Boundary variables data
 702         //     1. volFields
 703         //     2. surfFields
 704
 705         //Info<< "Writing regular boundary elements data ..." << endl;
 706
 707         writeInt(fvFile, FV_BNDRY_VARS);
 708
 709         forAll(volFieldPtrs, fieldI)
 710         {
 711             if (volFieldPtrs[fieldI] != NULL)
 712             {
 713                 const volScalarField& vsf = *volFieldPtrs[fieldI];
 714
 715                 forAll(mesh.boundary(), patchI)
 716                 {
 717                     writeFaceData
 718                     (
 719                         mesh,
 720                         topo,
 721                         patchI,
 722                         vsf.boundaryField()[patchI],
 723                         false,
 724                         fvFile
 725                     );
 726                 }
 727             }
 728             else
 729             {
 730                 forAll(mesh.boundaryMesh(), patchI)
 731                 {
 732                     // Dummy value.
 733                     floatField fField
 734                     (
 735                         mesh.boundaryMesh()[patchI].size()
 736                       - topo.nPolyFaces()[patchI],
 737                         0.0
 738                     );
 739
 740                     fvFile.write
 741                     (
 742                         reinterpret_cast<char*>(fField.begin()),
 743                         fField.size()*sizeof(float)
 744                     );
 745                 }
 746             }
 747         }
 748
 749         // surfFields
 750         forAll(surfFieldPtrs, fieldI)
 751         {
 752             if (surfFieldPtrs[fieldI] != NULL)
 753             {
 754                 const surfaceScalarField& ssf = *surfFieldPtrs[fieldI];
 755
 756                 forAll(mesh.boundary(), patchI)
 757                 {
 758                     writeFaceData
 759                     (
 760                         mesh,
 761                         topo,
 762                         patchI,
 763                         ssf.boundaryField()[patchI],
 764                         false,
 765                         fvFile
 766                     );
 767                 }
 768             }
 769             else
 770             {
 771                 forAll(mesh.boundaryMesh(), patchI)
 772                 {
 773                     // Dummy value.
 774                     floatField fField
 775                     (
 776                         mesh.boundaryMesh()[patchI].size()
 777                       - topo.nPolyFaces()[patchI],
 778                         0.0
 779                     );
 780
 781                     fvFile.write
 782                     (
 783                         reinterpret_cast<char*>(fField.begin()),
 784                         fField.size()*sizeof(float)
 785                     );
 786                 }
 787             }
 788         }
 789
 790         //
 791         // Polygonal faces boundary data
 792         //     1. volFields
 793         //     2. surfFields
 794
 795         //Info<< "Writing polygonal boundary elements data ..." << endl;
 796
 797         writeInt(fvFile, FV_ARB_POLY_BNDRY_VARS);
 798         forAll(volFieldPtrs, fieldI)
 799         {
 800             if (volFieldPtrs[fieldI] != NULL)
 801             {
 802                 const volScalarField& vsf = *volFieldPtrs[fieldI];
 803
 804                 // All non-empty patches
 805                 forAll(mesh.boundary(), patchI)
 806                 {
 807                     writeFaceData
 808                     (
 809                         mesh,
 810                         topo,
 811                         patchI,
 812                         vsf.boundaryField()[patchI],
 813                         true,
 814                         fvFile
 815                     );
 816                 }
 817             }
 818             else
 819             {
 820                 forAll(mesh.boundary(), patchI)
 821                 {
 822                     // Dummy value.
 823                     floatField fField(topo.nPolyFaces()[patchI], 0.0);
 824
 825                     fvFile.write
 826                     (
 827                         reinterpret_cast<char*>(fField.begin()),
 828                         fField.size()*sizeof(float)
 829                     );
 830                 }
 831             }
 832         }
 833
 834         // surfFields
 835         forAll(surfFieldPtrs, fieldI)
 836         {
 837             if (surfFieldPtrs[fieldI] != NULL)
 838             {
 839                 const surfaceScalarField& ssf = *surfFieldPtrs[fieldI];
 840
 841                 // All non-empty patches
 842                 forAll(mesh.boundary(), patchI)
 843                 {
 844                     writeFaceData
 845                     (
 846                         mesh,
 847                         topo,
 848                         patchI,
 849                         ssf.boundaryField()[patchI],
 850                         true,
 851                         fvFile
 852                     );
 853                 }
 854             }
 855             else
 856             {
 857                 forAll(mesh.boundaryMesh(), patchI)
 858                 {
 859                     // Dummy value.
 860                     floatField fField
 861                     (
 862                         mesh.boundaryMesh()[patchI].size()
 863                       - topo.nPolyFaces()[patchI],
 864                         0.0
 865                     );
 866
 867                     fvFile.write
 868                     (
 869                         reinterpret_cast<char*>(fField.begin()),
 870                         fField.size()*sizeof(float)
 871                     );
 872                 }
 873             }
 874         }
 875
 876
 877         //
 878         // Cleanup volume and surface fields
 879         //
 880         forAll(volFieldPtrs, fieldI)
 881         {
 882             delete volFieldPtrs[fieldI];
 883         }
 884         forAll(surfFieldPtrs, fieldI)
 885         {
 886             delete surfFieldPtrs[fieldI];
 887         }
 888
 889
 890
 891
 892         //
 893         // Spray
 894         //
 895         if (hasLagrangian)
 896         {
 897             // Read/create fields:
 898             //     sprayScalarFieldPtrs: List of ptrs to lagrangian scalfields
 899             //     sprayVectorFieldPtrs:               ,,           vec  ,,
 900 #           include "createSprayFields.H"
 901
 902
 903             // Write time header
 904
 905             // Time index (FieldView: has to start from 1)
 906             writeInt(fvParticleFile, fieldViewTime + 1);
 907
 908             // Time value
 909             writeFloat(fvParticleFile, runTime.value());
 910
 911             // Read particles
 912             Cloud<passiveParticle> parcels(mesh);
 913
 914             // Num particles
 915             writeInt(fvParticleFile, parcels.size());
 916
 917             Info<< "    Writing " << parcels.size() << " particles." << endl;
 918
 919
 920             //
 921             // Output data parcelwise
 922             //
 923
 924             label parcelNo = 0;
 925
 926
 927             for
 928             (
 929                 Cloud<passiveParticle>::iterator elmnt = parcels.begin();
 930                 elmnt != parcels.end();
 931                 ++elmnt, parcelNo++
 932             )
 933             {
 934                 writeInt(fvParticleFile, parcelNo+1);
 935
 936                 writeFloat(fvParticleFile, elmnt().position().x());
 937                 writeFloat(fvParticleFile, elmnt().position().y());
 938                 writeFloat(fvParticleFile, elmnt().position().z());
 939
 940                 forAll(sprayScalarFieldPtrs, fieldI)
 941                 {
 942                     if (sprayScalarFieldPtrs[fieldI] != NULL)
 943                     {
 944                         const IOField<scalar>& sprayField =
 945                             *sprayScalarFieldPtrs[fieldI];
 946                         writeFloat
 947                         (
 948                             fvParticleFile,
 949                             sprayField[parcelNo]
 950                         );
 951                     }
 952                     else
 953                     {
 954                         writeFloat(fvParticleFile, 0.0);
 955                     }
 956                 }
 957                 forAll(sprayVectorFieldPtrs, fieldI)
 958                 {
 959                     if (sprayVectorFieldPtrs[fieldI] != NULL)
 960                     {
 961                         const IOField<vector>& sprayVectorField =
 962                             *sprayVectorFieldPtrs[fieldI];
 963                         const vector& val =
 964                             sprayVectorField[parcelNo];
 965
 966                         writeFloat(fvParticleFile, val.x());
 967                         writeFloat(fvParticleFile, val.y());
 968                         writeFloat(fvParticleFile, val.z());
 969                     }
 970                     else
 971                     {
 972                         writeFloat(fvParticleFile, 0.0);
 973                         writeFloat(fvParticleFile, 0.0);
 974                         writeFloat(fvParticleFile, 0.0);
 975                     }
 976                 }
 977             }
 978
 979             // increment fieldView particle time
 980             fieldViewTime++;
 981
 982
 983             //
 984             // Cleanup spray fields
 985             //
 986             forAll(sprayScalarFieldPtrs, fieldI)
 987             {
 988                 delete sprayScalarFieldPtrs[fieldI];
 989             }
 990             forAll(sprayVectorFieldPtrs, fieldI)
 991             {
 992                 delete sprayVectorFieldPtrs[fieldI];
 993             }
 994
 995         } // end of hasLagrangian
 996     }
 997
 998     if (!hasLagrangian)
 999     {
1000         rm(fvParticleFileName);
1001     }
1002
1003     Info<< "End\n" << endl;
1004
1005     return 0;
1006 }
1007
1008
1009 // ************************************************************************* //