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

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