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