applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "ensightMesh.H"
  27 #include "argList.H"
  28 #include "Time.H"
  29 #include "fvMesh.H"
  30 #include "globalMeshData.H"
  31 #include "PstreamCombineReduceOps.H"
  32 #include "processorPolyPatch.H"
  33 #include "cellModeller.H"
  34 #include "IOmanip.H"
  35 #include "itoa.H"
  36 #include "globalIndex.H"
  37 #include "mapDistribute.H"
  38 #include "stringListOps.H"
  39
  40 #include "ensightBinaryStream.H"
  41 #include "ensightAsciiStream.H"
  42
  43 #include <fstream>
  44
  45 // * * * * * * * * * * * * * Private Functions * * * * * * * * * * * * * * //
  46
  47 void Foam::ensightMesh::correct()
  48 {
  49     patchPartOffset_ = 2;
  50     meshCellSets_ = mesh_.nCells();
  51     boundaryFaceSets_.setSize(mesh_.boundary().size());
  52     allPatchNames_.clear();
  53     patchNames_.clear();
  54     nPatchPrims_ = 0;
  55     faceZoneFaceSets_.setSize(mesh_.faceZones().size());
  56     faceZoneNames_.clear();
  57     nFaceZonePrims_ = 0;
  58     boundaryFaceToBeIncluded_.clear();
  59
  60     const cellShapeList& cellShapes = mesh_.cellShapes();
  61
  62     const cellModel& tet = *(cellModeller::lookup("tet"));
  63     const cellModel& pyr = *(cellModeller::lookup("pyr"));
  64     const cellModel& prism = *(cellModeller::lookup("prism"));
  65     const cellModel& wedge = *(cellModeller::lookup("wedge"));
  66     const cellModel& hex = *(cellModeller::lookup("hex"));
  67
  68     if (!noPatches_)
  69     {
  70         // Patches are output. Check that they're synced.
  71         mesh_.boundaryMesh().checkParallelSync(true);
  72
  73         allPatchNames_ = mesh_.boundaryMesh().names();
  74         if (Pstream::parRun())
  75         {
  76             allPatchNames_.setSize
  77             (
  78                 mesh_.boundary().size()
  79               - mesh_.globalData().processorPatches().size()
  80             );
  81         }
  82
  83         if (patches_)
  84         {
  85             if (patchPatterns_.empty())
  86             {
  87                 forAll(allPatchNames_, nameI)
  88                 {
  89                     patchNames_.insert(allPatchNames_[nameI]);
  90                 }
  91             }
  92             else
  93             {
  94                 // Find patch names which match that requested at command-line
  95                 forAll(allPatchNames_, nameI)
  96                 {
  97                     const word& patchName = allPatchNames_[nameI];
  98                     if (findStrings(patchPatterns_, patchName))
  99                     {
 100                         patchNames_.insert(patchName);
 101                     }
 102                 }
 103             }
 104         }
 105     }
 106
 107     if (patchNames_.size())
 108     {
 109         // no internalMesh
 110         patchPartOffset_ = 1;
 111     }
 112     else
 113     {
 114         // Count the shapes
 115         labelList& tets = meshCellSets_.tets;
 116         labelList& pyrs = meshCellSets_.pyrs;
 117         labelList& prisms = meshCellSets_.prisms;
 118         labelList& wedges = meshCellSets_.wedges;
 119         labelList& hexes = meshCellSets_.hexes;
 120         labelList& polys = meshCellSets_.polys;
 121
 122         label nTets = 0;
 123         label nPyrs = 0;
 124         label nPrisms = 0;
 125         label nWedges = 0;
 126         label nHexes = 0;
 127         label nPolys = 0;
 128
 129         forAll(cellShapes, cellI)
 130         {
 131             const cellShape& cellShape = cellShapes[cellI];
 132             const cellModel& cellModel = cellShape.model();
 133
 134             if (cellModel == tet)
 135             {
 136                 tets[nTets++] = cellI;
 137             }
 138             else if (cellModel == pyr)
 139             {
 140                 pyrs[nPyrs++] = cellI;
 141             }
 142             else if (cellModel == prism)
 143             {
 144                 prisms[nPrisms++] = cellI;
 145             }
 146             else if (cellModel == wedge)
 147             {
 148                 wedges[nWedges++] = cellI;
 149             }
 150             else if (cellModel == hex)
 151             {
 152                 hexes[nHexes++] = cellI;
 153             }
 154             else
 155             {
 156                 polys[nPolys++] = cellI;
 157             }
 158         }
 159
 160         tets.setSize(nTets);
 161         pyrs.setSize(nPyrs);
 162         prisms.setSize(nPrisms);
 163         wedges.setSize(nWedges);
 164         hexes.setSize(nHexes);
 165         polys.setSize(nPolys);
 166
 167         meshCellSets_.nTets = nTets;
 168         reduce(meshCellSets_.nTets, sumOp<label>());
 169
 170         meshCellSets_.nPyrs = nPyrs;
 171         reduce(meshCellSets_.nPyrs, sumOp<label>());
 172
 173         meshCellSets_.nPrisms = nPrisms;
 174         reduce(meshCellSets_.nPrisms, sumOp<label>());
 175
 176         meshCellSets_.nHexesWedges = nWedges+nHexes;
 177         reduce(meshCellSets_.nHexesWedges, sumOp<label>());
 178
 179         meshCellSets_.nPolys = nPolys;
 180         reduce(meshCellSets_.nPolys, sumOp<label>());
 181
 182
 183         // Determine parallel shared points
 184         globalPointsPtr_ = mesh_.globalData().mergePoints
 185         (
 186             pointToGlobal_,
 187             uniquePointMap_
 188         );
 189     }
 190
 191     if (!noPatches_)
 192     {
 193         forAll(mesh_.boundary(), patchi)
 194         {
 195             if (mesh_.boundary()[patchi].size())
 196             {
 197                 const polyPatch& p = mesh_.boundaryMesh()[patchi];
 198
 199                 labelList& tris = boundaryFaceSets_[patchi].tris;
 200                 labelList& quads = boundaryFaceSets_[patchi].quads;
 201                 labelList& polys = boundaryFaceSets_[patchi].polys;
 202
 203                 tris.setSize(p.size());
 204                 quads.setSize(p.size());
 205                 polys.setSize(p.size());
 206
 207                 label nTris = 0;
 208                 label nQuads = 0;
 209                 label nPolys = 0;
 210
 211                 forAll(p, faceI)
 212                 {
 213                     const face& f = p[faceI];
 214
 215                     if (f.size() == 3)
 216                     {
 217                         tris[nTris++] = faceI;
 218                     }
 219                     else if (f.size() == 4)
 220                     {
 221                         quads[nQuads++] = faceI;
 222                     }
 223                     else
 224                     {
 225                         polys[nPolys++] = faceI;
 226                     }
 227                 }
 228
 229                 tris.setSize(nTris);
 230                 quads.setSize(nQuads);
 231                 polys.setSize(nPolys);
 232             }
 233         }
 234     }
 235
 236     forAll(allPatchNames_, patchi)
 237     {
 238         const word& patchName = allPatchNames_[patchi];
 239         nFacePrimitives nfp;
 240
 241         if (patchNames_.empty() || patchNames_.found(patchName))
 242         {
 243             if (mesh_.boundary()[patchi].size())
 244             {
 245                 nfp.nTris   = boundaryFaceSets_[patchi].tris.size();
 246                 nfp.nQuads  = boundaryFaceSets_[patchi].quads.size();
 247                 nfp.nPolys  = boundaryFaceSets_[patchi].polys.size();
 248             }
 249         }
 250
 251         reduce(nfp.nTris, sumOp<label>());
 252         reduce(nfp.nQuads, sumOp<label>());
 253         reduce(nfp.nPolys, sumOp<label>());
 254
 255         nPatchPrims_.insert(patchName, nfp);
 256     }
 257
 258     // faceZones
 259     if (faceZones_)
 260     {
 261         const wordList faceZoneNamesAll = mesh_.faceZones().names();
 262
 263         // Find faceZone names which match that requested at command-line
 264         forAll(faceZoneNamesAll, nameI)
 265         {
 266             const word& zoneName = faceZoneNamesAll[nameI];
 267             if (findStrings(faceZonePatterns_, zoneName))
 268             {
 269                 faceZoneNames_.insert(zoneName);
 270             }
 271         }
 272
 273         // Build list of boundary faces to be exported
 274         boundaryFaceToBeIncluded_.setSize
 275         (
 276             mesh_.nFaces()
 277           - mesh_.nInternalFaces(),
 278             1
 279         );
 280
 281         forAll(mesh_.boundaryMesh(), patchI)
 282         {
 283             const polyPatch& pp = mesh_.boundaryMesh()[patchI];
 284             if
 285             (
 286                 isA<processorPolyPatch>(pp)
 287              && !refCast<const processorPolyPatch>(pp).owner()
 288             )
 289             {
 290                 label bFaceI = pp.start()-mesh_.nInternalFaces();
 291                 forAll(pp, i)
 292                 {
 293                     boundaryFaceToBeIncluded_[bFaceI++] = 0;
 294                 }
 295             }
 296         }
 297
 298         // Count face types in each faceZone
 299         forAll(faceZoneNamesAll, zoneI)
 300         {
 301             //const word& zoneName = faceZoneNamesAll[zoneI];
 302
 303             const faceZone& fz = mesh_.faceZones()[zoneI];
 304
 305             if (fz.size())
 306             {
 307                 labelList& tris = faceZoneFaceSets_[zoneI].tris;
 308                 labelList& quads = faceZoneFaceSets_[zoneI].quads;
 309                 labelList& polys = faceZoneFaceSets_[zoneI].polys;
 310
 311                 tris.setSize(fz.size());
 312                 quads.setSize(fz.size());
 313                 polys.setSize(fz.size());
 314
 315                 label nTris = 0;
 316                 label nQuads = 0;
 317                 label nPolys = 0;
 318
 319                 label faceCounter = 0;
 320
 321                 forAll(fz, i)
 322                 {
 323                     label faceI = fz[i];
 324
 325                     // Avoid counting faces on processor boundaries twice
 326                     if (faceToBeIncluded(faceI))
 327                     {
 328                         const face& f = mesh_.faces()[faceI];
 329
 330                         if (f.size() == 3)
 331                         {
 332                             tris[nTris++] = faceCounter;
 333                         }
 334                         else if (f.size() == 4)
 335                         {
 336                             quads[nQuads++] = faceCounter;
 337                         }
 338                         else
 339                         {
 340                             polys[nPolys++] = faceCounter;
 341                         }
 342
 343                         ++faceCounter;
 344                     }
 345                 }
 346
 347                 tris.setSize(nTris);
 348                 quads.setSize(nQuads);
 349                 polys.setSize(nPolys);
 350             }
 351         }
 352
 353         forAll(faceZoneNamesAll, zoneI)
 354         {
 355             const word& zoneName = faceZoneNamesAll[zoneI];
 356             nFacePrimitives nfp;
 357
 358             if (faceZoneNames_.found(zoneName))
 359             {
 360                 if
 361                 (
 362                     faceZoneFaceSets_[zoneI].tris.size()
 363                  || faceZoneFaceSets_[zoneI].quads.size()
 364                  || faceZoneFaceSets_[zoneI].polys.size()
 365                 )
 366                 {
 367                     nfp.nTris   = faceZoneFaceSets_[zoneI].tris.size();
 368                     nfp.nQuads  = faceZoneFaceSets_[zoneI].quads.size();
 369                     nfp.nPolys  = faceZoneFaceSets_[zoneI].polys.size();
 370                 }
 371             }
 372
 373             reduce(nfp.nTris, sumOp<label>());
 374             reduce(nfp.nQuads, sumOp<label>());
 375             reduce(nfp.nPolys, sumOp<label>());
 376
 377             nFaceZonePrims_.insert(zoneName, nfp);
 378         }
 379     }
 380 }
 381
 382
 383 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 384
 385 Foam::ensightMesh::ensightMesh
 386 (
 387     const fvMesh& mesh,
 388     const bool noPatches,
 389
 390     const bool patches,
 391     const wordReList& patchPatterns,
 392
 393     const bool faceZones,
 394     const wordReList& faceZonePatterns,
 395
 396     const bool binary
 397 )
 398 :
 399     mesh_(mesh),
 400     noPatches_(noPatches),
 401     patches_(patches),
 402     patchPatterns_(patchPatterns),
 403     faceZones_(faceZones),
 404     faceZonePatterns_(faceZonePatterns),
 405     binary_(binary),
 406     meshCellSets_(mesh.nCells())
 407 {
 408     correct();
 409 }
 410
 411
 412 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 413
 414 Foam::ensightMesh::~ensightMesh()
 415 {}
 416
 417
 418 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 419
 420 bool Foam::ensightMesh::faceToBeIncluded(const label faceI) const
 421 {
 422     bool res = false;
 423
 424     if (mesh_.isInternalFace(faceI))
 425     {
 426         res = true;
 427     }
 428     else
 429     {
 430         res = boundaryFaceToBeIncluded_[faceI-mesh_.nInternalFaces()];
 431     }
 432
 433     return res;
 434 }
 435
 436
 437 void Foam::ensightMesh::barrier()
 438 {
 439     label appI = 0;
 440     reduce(appI,maxOp<label>());
 441 }
 442
 443
 444 Foam::cellShapeList Foam::ensightMesh::map
 445 (
 446     const cellShapeList& cellShapes,
 447     const labelList& prims,
 448     const labelList& pointToGlobal
 449 ) const
 450 {
 451     cellShapeList mcsl(prims.size());
 452
 453     forAll(prims, i)
 454     {
 455         mcsl[i] = cellShapes[prims[i]];
 456         inplaceRenumber(pointToGlobal, mcsl[i]);
 457     }
 458
 459     return mcsl;
 460 }
 461
 462
 463 Foam::cellShapeList Foam::ensightMesh::map
 464 (
 465     const cellShapeList& cellShapes,
 466     const labelList& hexes,
 467     const labelList& wedges,
 468     const labelList& pointToGlobal
 469 ) const
 470 {
 471     cellShapeList mcsl(hexes.size() + wedges.size());
 472
 473     forAll(hexes, i)
 474     {
 475         mcsl[i] = cellShapes[hexes[i]];
 476         inplaceRenumber(pointToGlobal, mcsl[i]);
 477     }
 478
 479     label offset = hexes.size();
 480
 481     const cellModel& hex = *(cellModeller::lookup("hex"));
 482     labelList hexLabels(8);
 483
 484     forAll(wedges, i)
 485     {
 486         const cellShape& cellPoints = cellShapes[wedges[i]];
 487
 488         hexLabels[0] = cellPoints[0];
 489         hexLabels[1] = cellPoints[1];
 490         hexLabels[2] = cellPoints[0];
 491         hexLabels[3] = cellPoints[2];
 492         hexLabels[4] = cellPoints[3];
 493         hexLabels[5] = cellPoints[4];
 494         hexLabels[6] = cellPoints[6];
 495         hexLabels[7] = cellPoints[5];
 496
 497         mcsl[i + offset] = cellShape(hex, hexLabels);
 498         inplaceRenumber(pointToGlobal, mcsl[i + offset]);
 499     }
 500
 501     return mcsl;
 502 }
 503
 504
 505 void Foam::ensightMesh::writePrims
 506 (
 507     const cellShapeList& cellShapes,
 508     ensightStream& ensightGeometryFile
 509 ) const
 510 {
 511     // Create a temp int array
 512     if (cellShapes.size())
 513     {
 514         if (ensightGeometryFile.ascii())
 515         {
 516             // Workaround for paraview issue : write one cell per line
 517
 518             forAll(cellShapes, i)
 519             {
 520                 const cellShape& cellPoints = cellShapes[i];
 521
 522                 List<int> temp(cellPoints.size());
 523
 524                 forAll(cellPoints, pointI)
 525                 {
 526                     temp[pointI] = cellPoints[pointI] + 1;
 527                 }
 528                 ensightGeometryFile.write(temp);
 529             }
 530         }
 531         else
 532         {
 533             // All the cellShapes have the same number of elements!
 534             int numIntElem = cellShapes.size()*cellShapes[0].size();
 535             List<int> temp(numIntElem);
 536
 537             int n = 0;
 538
 539             forAll(cellShapes, i)
 540             {
 541                 const cellShape& cellPoints = cellShapes[i];
 542
 543                 forAll(cellPoints, pointI)
 544                 {
 545                     temp[n] = cellPoints[pointI] + 1;
 546                     n++;
 547                 }
 548             }
 549             ensightGeometryFile.write(temp);
 550         }
 551     }
 552 }
 553
 554
 555 void Foam::ensightMesh::writePolysNFaces
 556 (
 557     const labelList& polys,
 558     const cellList& cellFaces,
 559     ensightStream& ensightGeometryFile
 560 ) const
 561 {
 562     forAll(polys, i)
 563     {
 564         ensightGeometryFile.write(cellFaces[polys[i]].size());
 565     }
 566 }
 567
 568
 569 void Foam::ensightMesh::writePolysNPointsPerFace
 570 (
 571     const labelList& polys,
 572     const cellList& cellFaces,
 573     const faceList& faces,
 574     ensightStream& ensightGeometryFile
 575 ) const
 576 {
 577     forAll(polys, i)
 578     {
 579         const labelList& cf = cellFaces[polys[i]];
 580
 581         forAll(cf, faceI)
 582         {
 583             ensightGeometryFile.write(faces[cf[faceI]].size());
 584         }
 585     }
 586 }
 587
 588
 589 void Foam::ensightMesh::writePolysPoints
 590 (
 591     const labelList& polys,
 592     const cellList& cellFaces,
 593     const faceList& faces,
 594     ensightStream& ensightGeometryFile
 595 ) const
 596 {
 597     forAll(polys, i)
 598     {
 599         const labelList& cf = cellFaces[polys[i]];
 600
 601         forAll(cf, faceI)
 602         {
 603             const face& f = faces[cf[faceI]];
 604
 605             List<int> temp(f.size());
 606             forAll(f, pointI)
 607             {
 608                 temp[pointI] = f[pointI] + 1;
 609             }
 610             ensightGeometryFile.write(temp);
 611         }
 612     }
 613 }
 614
 615
 616 void Foam::ensightMesh::writeAllPolys
 617 (
 618     const labelList& pointToGlobal,
 619     ensightStream& ensightGeometryFile
 620 ) const
 621 {
 622     if (meshCellSets_.nPolys)
 623     {
 624         const cellList& cellFaces = mesh_.cells();
 625         // Renumber faces to use global point numbers
 626         faceList faces(mesh_.faces());
 627         forAll(faces, i)
 628         {
 629             inplaceRenumber(pointToGlobal, faces[i]);
 630         }
 631
 632         if (Pstream::master())
 633         {
 634             ensightGeometryFile.write("nfaced");
 635             ensightGeometryFile.write(meshCellSets_.nPolys);
 636         }
 637
 638         // Number of faces for each poly cell
 639
 640         if (Pstream::master())
 641         {
 642             // Master
 643             writePolysNFaces
 644             (
 645                 meshCellSets_.polys,
 646                 cellFaces,
 647                 ensightGeometryFile
 648             );
 649             // Slaves
 650             for (int slave=1; slave<Pstream::nProcs(); slave++)
 651             {
 652                 IPstream fromSlave(Pstream::scheduled, slave);
 653                 labelList polys(fromSlave);
 654                 cellList cellFaces(fromSlave);
 655
 656                 writePolysNFaces
 657                 (
 658                     polys,
 659                     cellFaces,
 660                     ensightGeometryFile
 661                 );
 662             }
 663         }
 664         else
 665         {
 666             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
 667             toMaster<< meshCellSets_.polys << cellFaces;
 668         }
 669
 670
 671         // Number of points for each face of the above list
 672         if (Pstream::master())
 673         {
 674             // Master
 675             writePolysNPointsPerFace
 676             (
 677                 meshCellSets_.polys,
 678                 cellFaces,
 679                 faces,
 680                 ensightGeometryFile
 681             );
 682             // Slaves
 683             for (int slave=1; slave<Pstream::nProcs(); slave++)
 684             {
 685                 IPstream fromSlave(Pstream::scheduled, slave);
 686                 labelList polys(fromSlave);
 687                 cellList cellFaces(fromSlave);
 688                 faceList faces(fromSlave);
 689
 690                 writePolysNPointsPerFace
 691                 (
 692                     polys,
 693                     cellFaces,
 694                     faces,
 695                     ensightGeometryFile
 696                 );
 697             }
 698         }
 699         else
 700         {
 701             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
 702             toMaster<< meshCellSets_.polys << cellFaces << faces;
 703         }
 704
 705
 706         // List of points id for each face of the above list
 707         if (Pstream::master())
 708         {
 709             // Master
 710             writePolysPoints
 711             (
 712                 meshCellSets_.polys,
 713                 cellFaces,
 714                 faces,
 715                 ensightGeometryFile
 716             );
 717             // Slaves
 718             for (int slave=1; slave<Pstream::nProcs(); slave++)
 719             {
 720                 IPstream fromSlave(Pstream::scheduled, slave);
 721                 labelList polys(fromSlave);
 722                 cellList cellFaces(fromSlave);
 723                 faceList faces(fromSlave);
 724
 725                 writePolysPoints
 726                 (
 727                     polys,
 728                     cellFaces,
 729                     faces,
 730                     ensightGeometryFile
 731                 );
 732             }
 733         }
 734         else
 735         {
 736             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
 737             toMaster<< meshCellSets_.polys << cellFaces << faces;
 738         }
 739     }
 740 }
 741
 742
 743 void Foam::ensightMesh::writeAllPrims
 744 (
 745     const char* key,
 746     const label nPrims,
 747     const cellShapeList& cellShapes,
 748     ensightStream& ensightGeometryFile
 749 ) const
 750 {
 751     if (nPrims)
 752     {
 753         if (Pstream::master())
 754         {
 755             ensightGeometryFile.write(key);
 756             ensightGeometryFile.write(nPrims);
 757
 758             writePrims(cellShapes, ensightGeometryFile);
 759
 760             for (int slave=1; slave<Pstream::nProcs(); slave++)
 761             {
 762                 IPstream fromSlave(Pstream::scheduled, slave);
 763                 cellShapeList cellShapes(fromSlave);
 764
 765                 writePrims(cellShapes, ensightGeometryFile);
 766             }
 767         }
 768         else
 769         {
 770             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
 771             toMaster<< cellShapes;
 772         }
 773     }
 774 }
 775
 776
 777 void Foam::ensightMesh::writeFacePrims
 778 (
 779     const faceList& patchFaces,
 780     ensightStream& ensightGeometryFile
 781 ) const
 782 {
 783     forAll(patchFaces, i)
 784     {
 785         const face& patchFace = patchFaces[i];
 786
 787         List<int> temp(patchFace.size());
 788         forAll(patchFace, pointI)
 789         {
 790             temp[pointI] = patchFace[pointI] + 1;
 791         }
 792
 793         ensightGeometryFile.write(temp);
 794     }
 795 }
 796
 797
 798 void Foam::ensightMesh::writeAllFacePrims
 799 (
 800     const char* key,
 801     const labelList& prims,
 802     const label nPrims,
 803     const faceList& patchFaces,
 804     ensightStream& ensightGeometryFile
 805 ) const
 806 {
 807     if (nPrims)
 808     {
 809         if (Pstream::master())
 810         {
 811             ensightGeometryFile.write(key);
 812             ensightGeometryFile.write(nPrims);
 813
 814             writeFacePrims
 815             (
 816                 UIndirectList<face>(patchFaces, prims)(),
 817                 ensightGeometryFile
 818             );
 819
 820             for (int slave=1; slave<Pstream::nProcs(); slave++)
 821             {
 822                 IPstream fromSlave(Pstream::scheduled, slave);
 823                 faceList patchFaces(fromSlave);
 824
 825                 writeFacePrims(patchFaces, ensightGeometryFile);
 826             }
 827         }
 828         else
 829         {
 830             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
 831             toMaster<< UIndirectList<face>(patchFaces, prims);
 832         }
 833     }
 834 }
 835
 836
 837 void Foam::ensightMesh::writeNSidedNPointsPerFace
 838 (
 839     const faceList& patchFaces,
 840     ensightStream& ensightGeometryFile
 841 ) const
 842 {
 843     forAll(patchFaces, i)
 844     {
 845         ensightGeometryFile.write(patchFaces[i].size());
 846     }
 847 }
 848
 849
 850 void Foam::ensightMesh::writeNSidedPoints
 851 (
 852     const faceList& patchFaces,
 853     ensightStream& ensightGeometryFile
 854 ) const
 855 {
 856     writeFacePrims(patchFaces, ensightGeometryFile);
 857 }
 858
 859
 860 void Foam::ensightMesh::writeAllNSided
 861 (
 862     const labelList& prims,
 863     const label nPrims,
 864     const faceList& patchFaces,
 865     ensightStream& ensightGeometryFile
 866 ) const
 867 {
 868     if (nPrims)
 869     {
 870         if (Pstream::master())
 871         {
 872             ensightGeometryFile.write("nsided");
 873             ensightGeometryFile.write(nPrims);
 874         }
 875
 876         // Number of points for each face
 877         if (Pstream::master())
 878         {
 879             writeNSidedNPointsPerFace
 880             (
 881                 UIndirectList<face>(patchFaces, prims)(),
 882                 ensightGeometryFile
 883             );
 884
 885             for (int slave=1; slave<Pstream::nProcs(); slave++)
 886             {
 887                 IPstream fromSlave(Pstream::scheduled, slave);
 888                 faceList patchFaces(fromSlave);
 889
 890                 writeNSidedNPointsPerFace
 891                 (
 892                     patchFaces,
 893                     ensightGeometryFile
 894                 );
 895             }
 896         }
 897         else
 898         {
 899             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
 900             toMaster<< UIndirectList<face>(patchFaces, prims);
 901         }
 902
 903         // List of points id for each face
 904         if (Pstream::master())
 905         {
 906             writeNSidedPoints
 907             (
 908                 UIndirectList<face>(patchFaces, prims)(),
 909                 ensightGeometryFile
 910             );
 911
 912             for (int slave=1; slave<Pstream::nProcs(); slave++)
 913             {
 914                 IPstream fromSlave(Pstream::scheduled, slave);
 915                 faceList patchFaces(fromSlave);
 916
 917                 writeNSidedPoints(patchFaces, ensightGeometryFile);
 918             }
 919         }
 920         else
 921         {
 922             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
 923             toMaster<< UIndirectList<face>(patchFaces, prims);
 924         }
 925     }
 926 }
 927
 928
 929 void Foam::ensightMesh::writeAllInternalPoints
 930 (
 931     const pointField& uniquePoints,
 932     const label nPoints,
 933     ensightStream& ensightGeometryFile
 934 ) const
 935 {
 936     barrier();
 937
 938     if (Pstream::master())
 939     {
 940         ensightGeometryFile.writePartHeader(1);
 941         ensightGeometryFile.write("internalMesh");
 942         ensightGeometryFile.write("coordinates");
 943         ensightGeometryFile.write(nPoints);
 944
 945         for (direction d=0; d<vector::nComponents; d++)
 946         {
 947             ensightGeometryFile.write(uniquePoints.component(d));
 948
 949             for (int slave=1; slave<Pstream::nProcs(); slave++)
 950             {
 951                 IPstream fromSlave(Pstream::scheduled, slave);
 952                 scalarField pointsComponent(fromSlave);
 953                 ensightGeometryFile.write(pointsComponent);
 954             }
 955         }
 956     }
 957     else
 958     {
 959         for (direction d=0; d<vector::nComponents; d++)
 960         {
 961             OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
 962             toMaster<< uniquePoints.component(d);
 963         }
 964     }
 965 }
 966
 967
 968 void Foam::ensightMesh::writeAllPatchPoints
 969 (
 970     const label ensightPatchI,
 971     const word& patchName,
 972     const pointField& uniquePoints,
 973     const label nPoints,
 974     ensightStream& ensightGeometryFile
 975 ) const
 976 {
 977     barrier();
 978
 979     if (Pstream::master())
 980     {
 981         ensightGeometryFile.writePartHeader(ensightPatchI);
 982         ensightGeometryFile.write(patchName.c_str());
 983         ensightGeometryFile.write("coordinates");
 984         ensightGeometryFile.write(nPoints);
 985
 986         for (direction d=0; d<vector::nComponents; d++)
 987         {
 988             ensightGeometryFile.write(uniquePoints.component(d));
 989             for (int slave=1; slave<Pstream::nProcs(); slave++)
 990             {
 991                 IPstream fromSlave(Pstream::scheduled, slave);
 992                 scalarField patchPointsComponent(fromSlave);
 993                 ensightGeometryFile.write(patchPointsComponent);
 994             }
 995         }
 996     }
 997     else
 998     {
 999         for (direction d=0; d<vector::nComponents; d++)
1000         {
1001             OPstream toMaster
1002             (
1003                 Pstream::scheduled,
1004                 Pstream::masterNo()
1005             );
1006             toMaster<< uniquePoints.component(d);
1007         }
1008     }
1009 }
1010
1011
1012 void Foam::ensightMesh::write
1013 (
1014     const fileName& postProcPath,
1015     const word& prepend,
1016     const label timeIndex,
1017     Ostream& ensightCaseFile
1018 ) const
1019 {
1020     const Time& runTime = mesh_.time();
1021     const cellShapeList& cellShapes = mesh_.cellShapes();
1022
1023
1024     word timeFile = prepend;
1025
1026     if (timeIndex == 0)
1027     {
1028         timeFile += "000.";
1029     }
1030     else if (mesh_.moving())
1031     {
1032         timeFile += itoa(timeIndex) + '.';
1033     }
1034
1035     // set the filename of the ensight file
1036     fileName ensightGeometryFileName = timeFile + "mesh";
1037
1038     ensightStream* ensightGeometryFilePtr = NULL;
1039     if (Pstream::master())
1040     {
1041         if (binary_)
1042         {
1043             ensightGeometryFilePtr = new ensightBinaryStream
1044             (
1045                 postProcPath/ensightGeometryFileName,
1046                 runTime
1047             );
1048             ensightGeometryFilePtr->write("C binary");
1049         }
1050         else
1051         {
1052             ensightGeometryFilePtr = new ensightAsciiStream
1053             (
1054                 postProcPath/ensightGeometryFileName,
1055                 runTime
1056             );
1057         }
1058     }
1059
1060     ensightStream& ensightGeometryFile = *ensightGeometryFilePtr;
1061
1062     if (Pstream::master())
1063     {
1064         string desc = string("written by OpenFOAM-") + Foam::FOAMversion;
1065
1066         ensightGeometryFile.write("EnSight Geometry File");
1067         ensightGeometryFile.write(desc.c_str());
1068         ensightGeometryFile.write("node id assign");
1069         ensightGeometryFile.write("element id assign");
1070     }
1071
1072     if (patchNames_.empty())
1073     {
1074         label nPoints = globalPoints().size();
1075
1076         const pointField uniquePoints(mesh_.points(), uniquePointMap_);
1077
1078         writeAllInternalPoints
1079         (
1080             uniquePoints,
1081             nPoints,
1082             ensightGeometryFile
1083         );
1084
1085         writeAllPrims
1086         (
1087             "hexa8",
1088             meshCellSets_.nHexesWedges,
1089             map         // Rewrite cellShapes to global numbering
1090             (
1091                 cellShapes,
1092                 meshCellSets_.hexes,
1093                 meshCellSets_.wedges,
1094                 pointToGlobal_
1095             ),
1096             ensightGeometryFile
1097         );
1098
1099         writeAllPrims
1100         (
1101             "penta6",
1102             meshCellSets_.nPrisms,
1103             map(cellShapes, meshCellSets_.prisms, pointToGlobal_),
1104             ensightGeometryFile
1105         );
1106
1107         writeAllPrims
1108         (
1109             "pyramid5",
1110             meshCellSets_.nPyrs,
1111             map(cellShapes, meshCellSets_.pyrs, pointToGlobal_),
1112             ensightGeometryFile
1113         );
1114
1115         writeAllPrims
1116         (
1117             "tetra4",
1118             meshCellSets_.nTets,
1119             map(cellShapes, meshCellSets_.tets, pointToGlobal_),
1120             ensightGeometryFile
1121         );
1122
1123         writeAllPolys
1124         (
1125             pointToGlobal_,
1126             ensightGeometryFile
1127         );
1128     }
1129
1130
1131     label ensightPatchI = patchPartOffset_;
1132
1133     forAll(allPatchNames_, patchi)
1134     {
1135         const word& patchName = allPatchNames_[patchi];
1136
1137         if (patchNames_.empty() || patchNames_.found(patchName))
1138         {
1139             const nFacePrimitives& nfp = nPatchPrims_[patchName];
1140
1141             if (nfp.nTris || nfp.nQuads || nfp.nPolys)
1142             {
1143                 const polyPatch& p = mesh_.boundaryMesh()[patchi];
1144                 const labelList& tris = boundaryFaceSets_[patchi].tris;
1145                 const labelList& quads = boundaryFaceSets_[patchi].quads;
1146                 const labelList& polys = boundaryFaceSets_[patchi].polys;
1147
1148                 // Renumber the patch points/faces into unique points
1149                 labelList pointToGlobal;
1150                 labelList uniqueMeshPointLabels;
1151                 autoPtr<globalIndex> globalPointsPtr =
1152                 mesh_.globalData().mergePoints
1153                 (
1154                     p.meshPoints(),
1155                     p.meshPointMap(),
1156                     pointToGlobal,
1157                     uniqueMeshPointLabels
1158                 );
1159
1160                 pointField uniquePoints(mesh_.points(), uniqueMeshPointLabels);
1161                 // Renumber the patch faces
1162                 faceList patchFaces(p.localFaces());
1163                 forAll(patchFaces, i)
1164                 {
1165                     inplaceRenumber(pointToGlobal, patchFaces[i]);
1166                 }
1167
1168                 writeAllPatchPoints
1169                 (
1170                     ensightPatchI++,
1171                     patchName,
1172                     uniquePoints,
1173                     globalPointsPtr().size(),
1174                     ensightGeometryFile
1175                 );
1176
1177                 writeAllFacePrims
1178                 (
1179                     "tria3",
1180                     tris,
1181                     nfp.nTris,
1182                     patchFaces,
1183                     ensightGeometryFile
1184                 );
1185
1186                 writeAllFacePrims
1187                 (
1188                     "quad4",
1189                     quads,
1190                     nfp.nQuads,
1191                     patchFaces,
1192                     ensightGeometryFile
1193                 );
1194
1195                 writeAllNSided
1196                 (
1197                     polys,
1198                     nfp.nPolys,
1199                     patchFaces,
1200                     ensightGeometryFile
1201                 );
1202             }
1203         }
1204     }
1205
1206     // write faceZones, if requested
1207     forAllConstIter(wordHashSet, faceZoneNames_, iter)
1208     {
1209         const word& faceZoneName = iter.key();
1210
1211         label faceID = mesh_.faceZones().findZoneID(faceZoneName);
1212
1213         const faceZone& fz = mesh_.faceZones()[faceID];
1214
1215         const nFacePrimitives& nfp = nFaceZonePrims_[faceZoneName];
1216
1217         if (nfp.nTris || nfp.nQuads || nfp.nPolys)
1218         {
1219             const labelList& tris = faceZoneFaceSets_[faceID].tris;
1220             const labelList& quads = faceZoneFaceSets_[faceID].quads;
1221             const labelList& polys = faceZoneFaceSets_[faceID].polys;
1222
1223             // Renumber the faceZone points/faces into unique points
1224             labelList pointToGlobal;
1225             labelList uniqueMeshPointLabels;
1226             autoPtr<globalIndex> globalPointsPtr =
1227             mesh_.globalData().mergePoints
1228             (
1229                 fz().meshPoints(),
1230                 fz().meshPointMap(),
1231                 pointToGlobal,
1232                 uniqueMeshPointLabels
1233             );
1234
1235             pointField uniquePoints(mesh_.points(), uniqueMeshPointLabels);
1236
1237             // Find the list of master faces belonging to the faceZone,
1238             // in loacal numbering
1239             faceList faceZoneFaces(fz().localFaces());
1240
1241             // Count how many master faces belong to the faceZone. Is there
1242             // a better way of doing this?
1243             label nMasterFaces = 0;
1244
1245             forAll(fz, faceI)
1246             {
1247                 if (faceToBeIncluded(fz[faceI]))
1248                 {
1249                     ++nMasterFaces;
1250                 }
1251             }
1252
1253             // Create the faceList for the master faces only and fill it.
1254             faceList faceZoneMasterFaces(nMasterFaces);
1255
1256             label currentFace = 0;
1257
1258             forAll(fz, faceI)
1259             {
1260                 if (faceToBeIncluded(fz[faceI]))
1261                 {
1262                     faceZoneMasterFaces[currentFace] = faceZoneFaces[faceI];
1263                     ++currentFace;
1264                 }
1265             }
1266
1267             // Renumber the faceZone master faces
1268             forAll(faceZoneMasterFaces, i)
1269             {
1270                 inplaceRenumber(pointToGlobal, faceZoneMasterFaces[i]);
1271             }
1272
1273             writeAllPatchPoints
1274             (
1275                 ensightPatchI++,
1276                 faceZoneName,
1277                 uniquePoints,
1278                 globalPointsPtr().size(),
1279                 ensightGeometryFile
1280             );
1281
1282             writeAllFacePrims
1283             (
1284                 "tria3",
1285                 tris,
1286                 nfp.nTris,
1287                 faceZoneMasterFaces,
1288                 ensightGeometryFile
1289             );
1290
1291             writeAllFacePrims
1292             (
1293                 "quad4",
1294                 quads,
1295                 nfp.nQuads,
1296                 faceZoneMasterFaces,
1297                 ensightGeometryFile
1298             );
1299
1300             writeAllNSided
1301             (
1302                 polys,
1303                 nfp.nPolys,
1304                 faceZoneMasterFaces,
1305                 ensightGeometryFile
1306             );
1307         }
1308     }
1309
1310     if (Pstream::master())
1311     {
1312         delete ensightGeometryFilePtr;
1313     }
1314 }
1315
1316
1317 // ************************************************************************* //