applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Application
  26     foamToVTK
  27
  28 Description
  29     Legacy VTK file format writer.
  30
  31     - handles volScalar, volVector, pointScalar, pointVector, surfaceScalar
  32       fields.
  33     - mesh topo changes.
  34     - both ascii and binary.
  35     - single time step writing.
  36     - write subset only.
  37     - automatic decomposition of cells; polygons on boundary undecomposed since
  38       handled by vtk.
  39
  40 Usage
  41
  42     - foamToVTK [OPTION]
  43
  44
  45     @param -ascii \n
  46     Write VTK data in ASCII format instead of binary.
  47
  48     @param -mesh \<name\>\n
  49     Use a different mesh name (instead of -region)
  50
  51     @param -fields \<fields\>\n
  52     Convert selected fields only. For example,
  53     @verbatim
  54          -fields "( p T U )"
  55     @endverbatim
  56     The quoting is required to avoid shell expansions and to pass the
  57     information as a single argument.
  58
  59     @param -surfaceFields \n
  60     Write surfaceScalarFields (e.g., phi)
  61
  62     @param -cellSet \<name\>\n
  63     @param -faceSet \<name\>\n
  64     @param -pointSet \<name\>\n
  65     Restrict conversion to the cellSet, faceSet or pointSet.
  66
  67     @param -nearCellValue \n
  68     Output cell value on patches instead of patch value itself
  69
  70     @param -noInternal \n
  71     Do not generate file for mesh, only for patches
  72
  73     @param -noPointValues \n
  74     No pointFields
  75
  76     @param -noFaceZones \n
  77     No faceZones
  78
  79     @param -noLinks \n
  80     (in parallel) do not link processor files to master
  81
  82     @param -allPatches \n
  83     Combine all patches into a single file
  84
  85     @param -allWallPatches \n
  86     Combine all wall patches into a single file
  87
  88     @param -excludePatches \<patchNames\>\n
  89     Specify patches to exclude. For example,
  90     @verbatim
  91          -excludePatches "( inlet_1 inlet_2 )"
  92     @endverbatim
  93     The quoting is required to avoid shell expansions and to pass the
  94     information as a single argument.
  95
  96     @param -useTimeName \n
  97     use the time index in the VTK file name instead of the time index
  98
  99 Note
 100     mesh subset is handled by vtkMesh. Slight inconsistency in
 101     interpolation: on the internal field it interpolates the whole volfield
 102     to the whole-mesh pointField and then selects only those values it
 103     needs for the subMesh (using the fvMeshSubset cellMap(), pointMap()
 104     functions). For the patches however it uses the
 105     fvMeshSubset.interpolate function to directly interpolate the
 106     whole-mesh values onto the subset patch.
 107
 108 \*---------------------------------------------------------------------------*/
 109
 110 #include "fvCFD.H"
 111 #include "pointMesh.H"
 112 #include "volPointInterpolation.H"
 113 #include "emptyPolyPatch.H"
 114 #include "labelIOField.H"
 115 #include "scalarIOField.H"
 116 #include "sphericalTensorIOField.H"
 117 #include "symmTensorIOField.H"
 118 #include "tensorIOField.H"
 119 #include "faceZoneMesh.H"
 120 #include "Cloud.H"
 121 #include "passiveParticle.H"
 122
 123 #include "faCFD.H"
 124
 125 #include "vtkMesh.H"
 126 #include "readFields.H"
 127 #include "writeFuns.H"
 128
 129 #include "internalWriter.H"
 130 #include "patchWriter.H"
 131 #include "faMeshWriter.H"
 132 #include "lagrangianWriter.H"
 133
 134 #include "writeFaceSet.H"
 135 #include "writePointSet.H"
 136 #include "writePatchGeom.H"
 137 #include "writeSurfFields.H"
 138
 139
 140 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 141
 142
 143 static const label VTK_TETRA = 10;
 144 static const label VTK_PYRAMID = 14;
 145 static const label VTK_WEDGE = 13;
 146 static const label VTK_HEXAHEDRON = 12;
 147
 148
 149 template<class GeoField>
 150 void print(const char* msg, Ostream& os, const PtrList<GeoField>& flds)
 151 {
 152     if (flds.size())
 153     {
 154         os << msg;
 155         forAll(flds, i)
 156         {
 157             os<< ' ' << flds[i].name();
 158         }
 159         os << endl;
 160     }
 161 }
 162
 163
 164 void print(Ostream& os, const wordList& flds)
 165 {
 166     forAll(flds, i)
 167     {
 168         os<< ' ' << flds[i];
 169     }
 170     os << endl;
 171 }
 172
 173
 174 labelList getSelectedPatches
 175 (
 176     const polyBoundaryMesh& patches,
 177     const HashSet<word>& excludePatches
 178 )
 179 {
 180     DynamicList<label> patchIDs(patches.size());
 181
 182     Info<< "Combining patches:" << endl;
 183
 184     forAll(patches, patchI)
 185     {
 186         const polyPatch& pp = patches[patchI];
 187
 188         if
 189         (
 190             isA<emptyPolyPatch>(pp)
 191             || (Pstream::parRun() && isA<processorPolyPatch>(pp))
 192         )
 193         {
 194             Info<< "    discarding empty/processor patch " << patchI
 195                 << " " << pp.name() << endl;
 196         }
 197         else if (!excludePatches.found(pp.name()))
 198         {
 199             patchIDs.append(patchI);
 200             Info<< "    patch " << patchI << " " << pp.name() << endl;
 201         }
 202     }
 203     return patchIDs.shrink();
 204 }
 205
 206
 207
 208
 209 // Main program:
 210
 211 int main(int argc, char *argv[])
 212 {
 213     timeSelector::addOptions();
 214
 215 #   include "addRegionOption.H"
 216
 217     argList::validOptions.insert("fields", "fields");
 218     argList::validOptions.insert("cellSet", "cellSet name");
 219     argList::validOptions.insert("faceSet", "faceSet name");
 220     argList::validOptions.insert("pointSet", "pointSet name");
 221     argList::validOptions.insert("ascii","");
 222     argList::validOptions.insert("surfaceFields","");
 223     argList::validOptions.insert("nearCellValue","");
 224     argList::validOptions.insert("noInternal","");
 225     argList::validOptions.insert("noPointValues","");
 226     argList::validOptions.insert("allPatches","");
 227     argList::validOptions.insert("allWallPatches","");
 228     argList::validOptions.insert("excludePatches","patches to exclude");
 229     argList::validOptions.insert("noFaceZones","");
 230     argList::validOptions.insert("faMesh","");
 231     argList::validOptions.insert("noLinks","");
 232     argList::validOptions.insert("useTimeName","");
 233
 234 #   include "setRootCase.H"
 235 #   include "createTime.H"
 236
 237     bool doWriteInternal = !args.optionFound("noInternal");
 238     bool doFaceZones = !args.optionFound("noFaceZones");
 239     bool doLinks = !args.optionFound("noLinks");
 240     bool binary = !args.optionFound("ascii");
 241     bool useTimeName = args.optionFound("useTimeName");
 242
 243     if (binary && (sizeof(floatScalar) != 4 || sizeof(label) != 4))
 244     {
 245         FatalErrorIn(args.executable())
 246             << "floatScalar and/or label are not 4 bytes in size" << nl
 247             << "Hence cannot use binary VTK format. Please use -ascii"
 248             << exit(FatalError);
 249     }
 250
 251     bool nearCellValue = args.optionFound("nearCellValue");
 252
 253     if (nearCellValue)
 254     {
 255         Info<< "Using neighbouring cell value instead of patch value"
 256             << nl << endl;
 257     }
 258
 259     bool noPointValues = args.optionFound("noPointValues");
 260
 261     if (noPointValues)
 262     {
 263         Info<< "Outputting cell values only" << nl << endl;
 264     }
 265
 266     bool allPatches = args.optionFound("allPatches");
 267     bool allWallPatches = args.optionFound("allWallPatches");
 268
 269     HashSet<word> excludePatches;
 270     if (args.optionFound("excludePatches"))
 271     {
 272         args.optionLookup("excludePatches")() >> excludePatches;
 273
 274         Info<< "Not including patches " << excludePatches << nl << endl;
 275     }
 276
 277     word cellSetName;
 278     string vtkName;
 279
 280     if (args.optionFound("cellSet"))
 281     {
 282         cellSetName = args.option("cellSet");
 283         vtkName = cellSetName;
 284     }
 285     else if (Pstream::parRun())
 286     {
 287         // Strip off leading casename, leaving just processor_DDD ending.
 288         vtkName = runTime.caseName();
 289
 290         string::size_type i = vtkName.rfind("processor");
 291
 292         if (i != string::npos)
 293         {
 294             vtkName = vtkName.substr(i);
 295         }
 296     }
 297     else
 298     {
 299         vtkName = runTime.caseName();
 300     }
 301
 302
 303     instantList timeDirs = timeSelector::select0(runTime, args);
 304
 305 #   include "createNamedMesh.H"
 306
 307     // VTK/ directory in the case
 308     fileName fvPath(runTime.path()/"VTK");
 309     // Directory of mesh (region0 gets filtered out)
 310     fileName regionPrefix = "";
 311
 312     if (regionName != polyMesh::defaultRegion)
 313     {
 314         fvPath = fvPath/regionName;
 315         regionPrefix = regionName;
 316     }
 317
 318     if (isDir(fvPath))
 319     {
 320         if
 321         (
 322             args.optionFound("time")
 323          || args.optionFound("latestTime")
 324          || cellSetName.size()
 325          || regionName != polyMesh::defaultRegion
 326         )
 327         {
 328             Info<< "Keeping old VTK files in " << fvPath << nl << endl;
 329         }
 330         else
 331         {
 332             Info<< "Deleting old VTK files in " << fvPath << nl << endl;
 333
 334             rmDir(fvPath);
 335         }
 336     }
 337
 338     mkDir(fvPath);
 339
 340
 341     // Mesh wrapper; does subsetting and decomposition
 342     vtkMesh vMesh
 343     (
 344         IOobject
 345         (
 346             regionName,
 347             runTime.timeName(),
 348             runTime,
 349             IOobject::MUST_READ
 350         ),
 351         cellSetName
 352     );
 353
 354     // Scan for all possible lagrangian clouds
 355     HashSet<fileName> allCloudDirs;
 356     forAll(timeDirs, timeI)
 357     {
 358         runTime.setTime(timeDirs[timeI], timeI);
 359         fileNameList cloudDirs
 360         (
 361             readDir
 362             (
 363                 runTime.timePath()/regionPrefix/cloud::prefix,
 364                 fileName::DIRECTORY
 365             )
 366         );
 367         forAll(cloudDirs, i)
 368         {
 369             IOobjectList sprayObjs
 370             (
 371                 mesh,
 372                 runTime.timeName(),
 373                 cloud::prefix/cloudDirs[i]
 374             );
 375
 376             IOobject* positionsPtr = sprayObjs.lookup("positions");
 377
 378             if (positionsPtr)
 379             {
 380                 if (allCloudDirs.insert(cloudDirs[i]))
 381                 {
 382                     Info<< "At time: " << runTime.timeName()
 383                         << " detected cloud directory : " << cloudDirs[i]
 384                         << endl;
 385                 }
 386             }
 387         }
 388     }
 389
 390
 391     forAll(timeDirs, timeI)
 392     {
 393         runTime.setTime(timeDirs[timeI], timeI);
 394
 395         Info<< "Time: " << runTime.timeName() << endl;
 396
 397         word timeDesc = "";
 398         if (useTimeName)
 399         {
 400             timeDesc = runTime.timeName();
 401         }
 402         else
 403         {
 404             timeDesc = name(runTime.timeIndex());
 405         }
 406
 407         // Check for new polyMesh/ and update mesh, fvMeshSubset and cell
 408         // decomposition.
 409         polyMesh::readUpdateState meshState = vMesh.readUpdate();
 410
 411         const fvMesh& mesh = vMesh.mesh();
 412
 413         if
 414         (
 415             meshState == polyMesh::TOPO_CHANGE
 416          || meshState == polyMesh::TOPO_PATCH_CHANGE
 417         )
 418         {
 419             Info<< "    Read new mesh" << nl << endl;
 420         }
 421
 422
 423         // If faceSet: write faceSet only (as polydata)
 424         if (args.optionFound("faceSet"))
 425         {
 426             // Load the faceSet
 427             faceSet set(mesh, args.option("faceSet"));
 428
 429             // Filename as if patch with same name.
 430             mkDir(fvPath/set.name());
 431
 432             fileName patchFileName
 433             (
 434                 fvPath/set.name()/set.name()
 435               + "_"
 436               + timeDesc
 437               + ".vtk"
 438             );
 439
 440             Info<< "    FaceSet   : " << patchFileName << endl;
 441
 442             writeFaceSet(binary, vMesh, set, patchFileName);
 443
 444             continue;
 445         }
 446         // If pointSet: write pointSet only (as polydata)
 447         if (args.optionFound("pointSet"))
 448         {
 449             // Load the pointSet
 450             pointSet set(mesh, args.option("pointSet"));
 451
 452             // Filename as if patch with same name.
 453             mkDir(fvPath/set.name());
 454
 455             fileName patchFileName
 456             (
 457                 fvPath/set.name()/set.name()
 458               + "_"
 459               + timeDesc
 460               + ".vtk"
 461             );
 462
 463             Info<< "    pointSet   : " << patchFileName << endl;
 464
 465             writePointSet(binary, vMesh, set, patchFileName);
 466
 467             continue;
 468         }
 469
 470
 471         // Search for list of objects for this time
 472         IOobjectList objects(mesh, runTime.timeName());
 473
 474         HashSet<word> selectedFields;
 475         if (args.optionFound("fields"))
 476         {
 477             args.optionLookup("fields")() >> selectedFields;
 478         }
 479
 480         // Construct the vol fields (on the original mesh if subsetted)
 481
 482         PtrList<volScalarField> vsf;
 483         readFields(vMesh, vMesh, objects, selectedFields, vsf);
 484         readFields(vMesh, vMesh, objects, selectedFields, vsf);
 485         print("    volScalarFields            :", Info, vsf);
 486
 487         PtrList<volVectorField> vvf;
 488         readFields(vMesh, vMesh, objects, selectedFields, vvf);
 489         print("    volVectorFields            :", Info, vvf);
 490
 491         PtrList<volSphericalTensorField> vSpheretf;
 492         readFields(vMesh, vMesh, objects, selectedFields, vSpheretf);
 493         print("    volSphericalTensorFields   :", Info, vSpheretf);
 494
 495         PtrList<volSymmTensorField> vSymmtf;
 496         readFields(vMesh, vMesh, objects, selectedFields, vSymmtf);
 497         print("    volSymmTensorFields        :", Info, vSymmtf);
 498
 499         PtrList<volTensorField> vtf;
 500         readFields(vMesh, vMesh, objects, selectedFields, vtf);
 501         print("    volTensorFields            :", Info, vtf);
 502
 503         const label nVolFields =
 504                 vsf.size()
 505               + vvf.size()
 506               + vSpheretf.size()
 507               + vSymmtf.size()
 508               + vtf.size();
 509
 510
 511         // Construct pointMesh only if nessecary since constructs edge
 512         // addressing (expensive on polyhedral meshes)
 513         if (noPointValues)
 514         {
 515             Info<< "    pointScalarFields : switched off"
 516                 << " (\"-noPointValues\" option)\n";
 517             Info<< "    pointVectorFields : switched off"
 518                 << " (\"-noPointValues\" option)\n";
 519         }
 520
 521         PtrList<pointScalarField> psf;
 522         PtrList<pointVectorField> pvf;
 523         PtrList<pointSphericalTensorField> pSpheretf;
 524         PtrList<pointSymmTensorField> pSymmtf;
 525         PtrList<pointTensorField> ptf;
 526
 527         if (!noPointValues)
 528         {
 529             readFields
 530             (
 531                 vMesh,
 532                 pointMesh::New(vMesh),
 533                 objects,
 534                 selectedFields,
 535                 psf
 536             );
 537             print("    pointScalarFields          :", Info, psf);
 538
 539             readFields
 540             (
 541                 vMesh,
 542                 pointMesh::New(vMesh),
 543                 objects,
 544                 selectedFields,
 545                 pvf
 546             );
 547             print("    pointVectorFields          :", Info, pvf);
 548
 549             readFields
 550             (
 551                 vMesh,
 552                 pointMesh::New(vMesh),
 553                 objects,
 554                 selectedFields,
 555                 pSpheretf
 556             );
 557             print("    pointSphericalTensorFields :", Info, pSpheretf);
 558
 559             readFields
 560             (
 561                 vMesh,
 562                 pointMesh::New(vMesh),
 563                 objects,
 564                 selectedFields,
 565                 pSymmtf
 566             );
 567             print("    pointSymmTensorFields      :", Info, pSymmtf);
 568
 569             readFields
 570             (
 571                 vMesh,
 572                 pointMesh::New(vMesh),
 573                 objects,
 574                 selectedFields,
 575                 ptf
 576             );
 577             print("    pointTensorFields          :", Info, ptf);
 578         }
 579         Info<< endl;
 580
 581         label nPointFields =
 582             psf.size()
 583           + pvf.size()
 584           + pSpheretf.size()
 585           + pSymmtf.size()
 586           + ptf.size();
 587
 588         if (doWriteInternal)
 589         {
 590             //
 591             // Create file and write header
 592             //
 593             fileName vtkFileName
 594             (
 595                 fvPath/vtkName
 596               + "_"
 597               + timeDesc
 598               + ".vtk"
 599             );
 600
 601             Info<< "    Internal  : " << vtkFileName << endl;
 602
 603             // Write mesh
 604             internalWriter writer(vMesh, binary, vtkFileName);
 605
 606             // VolFields + cellID
 607             writeFuns::writeCellDataHeader
 608             (
 609                 writer.os(),
 610                 vMesh.nFieldCells(),
 611                 1+nVolFields
 612             );
 613
 614             // Write cellID field
 615             writer.writeCellIDs();
 616
 617             // Write volFields
 618             writer.write(vsf);
 619             writer.write(vvf);
 620             writer.write(vSpheretf);
 621             writer.write(vSymmtf);
 622             writer.write(vtf);
 623
 624             if (!noPointValues)
 625             {
 626                 writeFuns::writePointDataHeader
 627                 (
 628                     writer.os(),
 629                     vMesh.nFieldPoints(),
 630                     nVolFields+nPointFields
 631                 );
 632
 633                 // pointFields
 634                 writer.write(psf);
 635                 writer.write(pvf);
 636                 writer.write(pSpheretf);
 637                 writer.write(pSymmtf);
 638                 writer.write(ptf);
 639
 640                 // Interpolated volFields
 641                 volPointInterpolation pInterp(mesh);
 642                 writer.write(pInterp, vsf);
 643                 writer.write(pInterp, vvf);
 644                 writer.write(pInterp, vSpheretf);
 645                 writer.write(pInterp, vSymmtf);
 646                 writer.write(pInterp, vtf);
 647             }
 648         }
 649
 650         //---------------------------------------------------------------------
 651         //
 652         // Write surface fields
 653         //
 654         //---------------------------------------------------------------------
 655
 656         if (args.optionFound("surfaceFields"))
 657         {
 658             PtrList<surfaceScalarField> ssf;
 659             readFields
 660             (
 661                 vMesh,
 662                 vMesh,
 663                 objects,
 664                 selectedFields,
 665                 ssf
 666             );
 667             print("    surfScalarFields  :", Info, ssf);
 668
 669             PtrList<surfaceVectorField> svf;
 670             readFields
 671             (
 672                 vMesh,
 673                 vMesh,
 674                 objects,
 675                 selectedFields,
 676                 svf
 677             );
 678             print("    surfVectorFields  :", Info, svf);
 679
 680             if (ssf.size() + svf.size() > 0)
 681             {
 682                 // Rework the scalar fields into vectorfields.
 683                 label sz = svf.size();
 684
 685                 svf.setSize(sz+ssf.size());
 686
 687                 surfaceVectorField n(mesh.Sf()/mesh.magSf());
 688
 689                 forAll(ssf, i)
 690                 {
 691                     svf.set(sz+i, ssf[i]*n);
 692                     svf[sz+i].rename(ssf[i].name());
 693                 }
 694                 ssf.clear();
 695
 696                 mkDir(fvPath / "surfaceFields");
 697
 698                 fileName surfFileName
 699                 (
 700                     fvPath
 701                    /"surfaceFields"
 702                    /"surfaceFields"
 703                    + "_"
 704                    + timeDesc
 705                    + ".vtk"
 706                 );
 707
 708                 writeSurfFields
 709                 (
 710                     binary,
 711                     vMesh,
 712                     surfFileName,
 713                     svf
 714                 );
 715             }
 716         }
 717
 718
 719         //---------------------------------------------------------------------
 720         //
 721         // Write patches (POLYDATA file, one for each patch)
 722         //
 723         //---------------------------------------------------------------------
 724
 725         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 726
 727         if (allPatches)
 728         {
 729             mkDir(fvPath/"allPatches");
 730
 731             fileName patchFileName;
 732
 733             if (vMesh.useSubMesh())
 734             {
 735                 patchFileName =
 736                     fvPath/"allPatches"/cellSetName
 737                   + "_"
 738                   + timeDesc
 739                   + ".vtk";
 740             }
 741             else
 742             {
 743                 patchFileName =
 744                     fvPath/"allPatches"/"allPatches"
 745                   + "_"
 746                   + timeDesc
 747                   + ".vtk";
 748             }
 749
 750             Info<< "    Combined patches     : " << patchFileName << endl;
 751
 752             patchWriter writer
 753             (
 754                 vMesh,
 755                 binary,
 756                 nearCellValue,
 757                 patchFileName,
 758                 getSelectedPatches(patches, excludePatches)
 759             );
 760
 761             // VolFields + patchID
 762             writeFuns::writeCellDataHeader
 763             (
 764                 writer.os(),
 765                 writer.nFaces(),
 766                 1+nVolFields
 767             );
 768
 769             // Write patchID field
 770             writer.writePatchIDs();
 771
 772             // Write volFields
 773             writer.write(vsf);
 774             writer.write(vvf);
 775             writer.write(vSpheretf);
 776             writer.write(vSymmtf);
 777             writer.write(vtf);
 778
 779             if (!noPointValues)
 780             {
 781                 writeFuns::writePointDataHeader
 782                 (
 783                     writer.os(),
 784                     writer.nPoints(),
 785                     nPointFields
 786                 );
 787
 788                 // Write pointFields
 789                 writer.write(psf);
 790                 writer.write(pvf);
 791                 writer.write(pSpheretf);
 792                 writer.write(pSymmtf);
 793                 writer.write(ptf);
 794
 795                 // no interpolated volFields since I cannot be bothered to
 796                 // create the patchInterpolation for all subpatches.
 797             }
 798         }
 799         else if (allWallPatches)
 800         {
 801             mkDir(fvPath/"allWallPatches");
 802
 803             fileName wallPatchFileName;
 804
 805             if (vMesh.useSubMesh())
 806             {
 807                 wallPatchFileName =
 808                     fvPath/"allWallPatches"/cellSetName
 809                   + "_"
 810                   + timeDesc
 811                   + ".vtk";
 812             }
 813             else
 814             {
 815                 wallPatchFileName =
 816                     fvPath/"allWallPatches"/"allWallPatches"
 817                   + "_"
 818                   + timeDesc
 819                   + ".vtk";
 820             }
 821
 822             // Go through the boundary and add all patches that are not
 823             // of type wall onto the exclude{atches list
 824             forAll (patches, patchI)
 825             {
 826                 if (!isA<wallPolyPatch>(patches[patchI]))
 827                 {
 828                     excludePatches.insert(patches[patchI].name());
 829                 }
 830             }
 831
 832             Info<< "    Combined wall patches     : " << wallPatchFileName
 833                 << endl;
 834
 835             patchWriter writer
 836             (
 837                 vMesh,
 838                 binary,
 839                 nearCellValue,
 840                 wallPatchFileName,
 841                 getSelectedPatches(patches, excludePatches)
 842             );
 843
 844             // VolFields + patchID
 845             writeFuns::writeCellDataHeader
 846             (
 847                 writer.os(),
 848                 writer.nFaces(),
 849                 1+nVolFields
 850             );
 851
 852             // Write patchID field
 853             writer.writePatchIDs();
 854
 855             // Write volFields
 856             writer.write(vsf);
 857             writer.write(vvf);
 858             writer.write(vSpheretf);
 859             writer.write(vSymmtf);
 860             writer.write(vtf);
 861
 862             if (!noPointValues)
 863             {
 864                 writeFuns::writePointDataHeader
 865                 (
 866                     writer.os(),
 867                     writer.nPoints(),
 868                     nPointFields
 869                 );
 870
 871                 // Write pointFields
 872                 writer.write(psf);
 873                 writer.write(pvf);
 874                 writer.write(pSpheretf);
 875                 writer.write(pSymmtf);
 876                 writer.write(ptf);
 877
 878                 // no interpolated volFields since I cannot be bothered to
 879                 // create the patchInterpolation for all subpatches.
 880             }
 881         }
 882         else
 883         {
 884             forAll(patches, patchI)
 885             {
 886                 const polyPatch& pp = patches[patchI];
 887
 888                 if (!excludePatches.found(pp.name()))
 889                 {
 890                     mkDir(fvPath/pp.name());
 891
 892                     fileName patchFileName;
 893
 894                     if (vMesh.useSubMesh())
 895                     {
 896                         patchFileName =
 897                             fvPath/pp.name()/cellSetName
 898                           + "_"
 899                           + timeDesc
 900                           + ".vtk";
 901                     }
 902                     else
 903                     {
 904                         patchFileName =
 905                             fvPath/pp.name()/pp.name()
 906                           + "_"
 907                           + timeDesc
 908                           + ".vtk";
 909                     }
 910
 911                     Info<< "    Patch     : " << patchFileName << endl;
 912
 913                     patchWriter writer
 914                     (
 915                         vMesh,
 916                         binary,
 917                         nearCellValue,
 918                         patchFileName,
 919                         labelList(1, patchI)
 920                     );
 921
 922                     if (!isA<emptyPolyPatch>(pp))
 923                     {
 924                         // VolFields + patchID
 925                         writeFuns::writeCellDataHeader
 926                         (
 927                             writer.os(),
 928                             writer.nFaces(),
 929                             1+nVolFields
 930                         );
 931
 932                         // Write patchID field
 933                         writer.writePatchIDs();
 934
 935                         // Write volFields
 936                         writer.write(vsf);
 937                         writer.write(vvf);
 938                         writer.write(vSpheretf);
 939                         writer.write(vSymmtf);
 940                         writer.write(vtf);
 941
 942                         if (!noPointValues)
 943                         {
 944                             writeFuns::writePointDataHeader
 945                             (
 946                                 writer.os(),
 947                                 writer.nPoints(),
 948                                 nVolFields
 949                               + nPointFields
 950                             );
 951
 952                             // Write pointFields
 953                             writer.write(psf);
 954                             writer.write(pvf);
 955                             writer.write(pSpheretf);
 956                             writer.write(pSymmtf);
 957                             writer.write(ptf);
 958
 959                             PrimitivePatchInterpolation<primitivePatch> pInter
 960                             (
 961                                 pp
 962                             );
 963
 964                             // Write interpolated volFields
 965                             writer.write(pInter, vsf);
 966                             writer.write(pInter, vvf);
 967                             writer.write(pInter, vSpheretf);
 968                             writer.write(pInter, vSymmtf);
 969                             writer.write(pInter, vtf);
 970                         }
 971                     }
 972                 }
 973             }
 974         }
 975
 976         //---------------------------------------------------------------------
 977         //
 978         // Write faceZones (POLYDATA file, one for each zone)
 979         //
 980         //---------------------------------------------------------------------
 981
 982         if (doFaceZones)
 983         {
 984             const faceZoneMesh& zones = mesh.faceZones();
 985
 986             forAll(zones, zoneI)
 987             {
 988                 const faceZone& pp = zones[zoneI];
 989
 990                 mkDir(fvPath/pp.name());
 991
 992                 fileName patchFileName;
 993
 994                 if (vMesh.useSubMesh())
 995                 {
 996                     patchFileName =
 997                         fvPath/pp.name()/cellSetName
 998                       + "_"
 999                       + timeDesc
1000                       + ".vtk";
1001                 }
1002                 else
1003                 {
1004                     patchFileName =
1005                         fvPath/pp.name()/pp.name()
1006                       + "_"
1007                       + timeDesc
1008                       + ".vtk";
1009                 }
1010
1011                 Info<< "    FaceZone  : " << patchFileName << endl;
1012
1013                 std::ofstream str(patchFileName.c_str());
1014
1015                 writeFuns::writeHeader(str, binary, pp.name());
1016                 str << "DATASET POLYDATA" << std::endl;
1017
1018                 writePatchGeom
1019                 (
1020                     binary,
1021                     pp().localFaces(),
1022                     pp().localPoints(),
1023                     str
1024                 );
1025             }
1026         }
1027
1028
1029
1030         //---------------------------------------------------------------------
1031         //
1032         // Write finite area data
1033         //
1034         //---------------------------------------------------------------------
1035
1036         if (args.options().found("faMesh"))
1037         {
1038             mkDir(fvPath/"faMesh");
1039
1040             fileName faFileName =
1041                 fvPath/"faMesh"/"faMesh"
1042               + "_"
1043               + name(timeI)
1044               + ".vtk";
1045
1046             // Create FA mesh
1047             faMesh aMesh(vMesh.mesh());
1048
1049             PtrList<areaScalarField> asf;
1050             readFieldsNoSubset(vMesh, aMesh, objects, selectedFields, asf);
1051             print("    areaScalarFields           :", Info, asf);
1052
1053             PtrList<areaVectorField> avf;
1054             readFieldsNoSubset(vMesh, aMesh, objects, selectedFields, avf);
1055             print("    areaVectorFields           :", Info, avf);
1056
1057             const label nAreaFields =
1058                 asf.size()
1059               + avf.size();
1060
1061             faMeshWriter writer
1062             (
1063                 aMesh,
1064                 binary,
1065                 faFileName
1066             );
1067
1068             writeFuns::writeCellDataHeader
1069             (
1070                 writer.os(),
1071                 aMesh.nFaces(),
1072                 nAreaFields
1073             );
1074
1075             // Write areaFields
1076             writer.write(asf);
1077             writer.write(avf);
1078
1079             if (!noPointValues)
1080             {
1081                 writeFuns::writePointDataHeader
1082                 (
1083                     writer.os(),
1084                     aMesh.nPoints(),
1085                     nAreaFields
1086                 );
1087
1088                 // Interpolated areaFields
1089                 PrimitivePatchInterpolation<indirectPrimitivePatch> pInterp
1090                 (
1091                     aMesh.patch()
1092                 );
1093                 writer.write(pInterp, asf);
1094                 writer.write(pInterp, avf);
1095             }
1096         }
1097
1098         //---------------------------------------------------------------------
1099         //
1100         // Write lagrangian data
1101         //
1102         //---------------------------------------------------------------------
1103
1104         forAllConstIter(HashSet<fileName>, allCloudDirs, iter)
1105         {
1106             const fileName& cloudName = iter.key();
1107
1108             // Always create the cloud directory.
1109             mkDir(fvPath/cloud::prefix/cloudName);
1110
1111             fileName lagrFileName
1112             (
1113                 fvPath/cloud::prefix/cloudName/cloudName
1114               + "_" + timeDesc + ".vtk"
1115             );
1116
1117             Info<< "    Lagrangian: " << lagrFileName << endl;
1118
1119
1120             IOobjectList sprayObjs
1121             (
1122                 mesh,
1123                 runTime.timeName(),
1124                 cloud::prefix/cloudName
1125             );
1126
1127             IOobject* positionsPtr = sprayObjs.lookup("positions");
1128
1129             if (positionsPtr)
1130             {
1131                 wordList labelNames(sprayObjs.names(labelIOField::typeName));
1132                 Info<< "        labels            :";
1133                 print(Info, labelNames);
1134
1135                 wordList scalarNames(sprayObjs.names(scalarIOField::typeName));
1136                 Info<< "        scalars           :";
1137                 print(Info, scalarNames);
1138
1139                 wordList vectorNames(sprayObjs.names(vectorIOField::typeName));
1140                 Info<< "        vectors           :";
1141                 print(Info, vectorNames);
1142
1143                 wordList sphereNames
1144                 (
1145                     sprayObjs.names
1146                     (
1147                         sphericalTensorIOField::typeName
1148                     )
1149                 );
1150                 Info<< "        spherical tensors :";
1151                 print(Info, sphereNames);
1152
1153                 wordList symmNames
1154                 (
1155                     sprayObjs.names
1156                     (
1157                         symmTensorIOField::typeName
1158                     )
1159                 );
1160                 Info<< "        symm tensors      :";
1161                 print(Info, symmNames);
1162
1163                 wordList tensorNames(sprayObjs.names(tensorIOField::typeName));
1164                 Info<< "        tensors           :";
1165                 print(Info, tensorNames);
1166
1167                 lagrangianWriter writer
1168                 (
1169                     vMesh,
1170                     binary,
1171                     lagrFileName,
1172                     cloudName,
1173                     false
1174                 );
1175
1176                 // Write number of fields
1177                 writer.writeParcelHeader
1178                 (
1179                     labelNames.size()
1180                   + scalarNames.size()
1181                   + vectorNames.size()
1182                   + sphereNames.size()
1183                   + symmNames.size()
1184                   + tensorNames.size()
1185                 );
1186
1187                 // Fields
1188                 writer.writeIOField<label>(labelNames);
1189                 writer.writeIOField<scalar>(scalarNames);
1190                 writer.writeIOField<vector>(vectorNames);
1191                 writer.writeIOField<sphericalTensor>(sphereNames);
1192                 writer.writeIOField<symmTensor>(symmNames);
1193                 writer.writeIOField<tensor>(tensorNames);
1194             }
1195             else
1196             {
1197                 lagrangianWriter writer
1198                 (
1199                     vMesh,
1200                     binary,
1201                     lagrFileName,
1202                     cloudName,
1203                     true
1204                 );
1205
1206                 // Write number of fields
1207                 writer.writeParcelHeader(0);
1208             }
1209         }
1210     }
1211
1212
1213     //---------------------------------------------------------------------
1214     //
1215     // Link parallel outputs back to undecomposed case for ease of loading
1216     //
1217     //---------------------------------------------------------------------
1218
1219     if (Pstream::parRun() && doLinks)
1220     {
1221         mkDir(runTime.path()/".."/"VTK");
1222         chDir(runTime.path()/".."/"VTK");
1223
1224         Info<< "Linking all processor files to " << runTime.path()/".."/"VTK"
1225             << endl;
1226
1227         // Get list of vtk files
1228         fileName procVTK
1229         (
1230             fileName("..")
1231            /"processor" + name(Pstream::myProcNo())
1232            /"VTK"
1233         );
1234
1235         fileNameList dirs(readDir(procVTK, fileName::DIRECTORY));
1236         label sz = dirs.size();
1237         dirs.setSize(sz+1);
1238         dirs[sz] = ".";
1239
1240         forAll(dirs, i)
1241         {
1242             fileNameList subFiles(readDir(procVTK/dirs[i], fileName::FILE));
1243
1244             forAll(subFiles, j)
1245             {
1246                 fileName procFile(procVTK/dirs[i]/subFiles[j]);
1247
1248                 if (exists(procFile))
1249                 {
1250                     string cmd
1251                     (
1252                         "ln -s "
1253                       + procFile
1254                       + " "
1255                       + "processor"
1256                       + name(Pstream::myProcNo())
1257                       + "_"
1258                       + procFile.name()
1259                     );
1260                     system(cmd.c_str());
1261                 }
1262             }
1263         }
1264     }
1265
1266     Info<< "End\n" << endl;
1267
1268     return 0;
1269 }
1270
1271
1272 // ************************************************************************* //