applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.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 Description
  25     Reads .msh file as written by Gmsh.
  26
  27     Needs surface elements on mesh to be present and aligned with outside faces
  28     of the mesh. I.e. if the mesh is hexes, the outside faces need to be quads
  29
  30     Note: There is something seriously wrong with the ordering written in the
  31     .msh file. Normal operation is to check the ordering and invert prisms
  32     and hexes if found to be wrong way round.
  33     Use the -keepOrientation to keep the raw information.
  34
  35     Note: The code now uses the element (cell,face) physical region id number
  36     to create cell zones and faces zones (similar to
  37     fluentMeshWithInternalFaces).
  38
  39     A use of the cell zone information, is for field initialization with the
  40     "setFields" utility. see the classes:  topoSetSource, zoneToCell.
  41 \*---------------------------------------------------------------------------*/
  42
  43 #include "argList.H"
  44 #include "polyMesh.H"
  45 #include "foamTime.H"
  46 #include "polyMesh.H"
  47 #include "IFstream.H"
  48 #include "cellModeller.H"
  49 #include "repatchPolyTopoChanger.H"
  50 #include "cellSet.H"
  51 #include "faceSet.H"
  52
  53 using namespace Foam;
  54
  55 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  56
  57 // Element type numbers
  58 static label MSHTRI   = 2;
  59 static label MSHQUAD  = 3;
  60 static label MSHTET   = 4;
  61 static label MSHPYR   = 7;
  62 static label MSHPRISM = 6;
  63 static label MSHHEX   = 5;
  64
  65
  66 // Skips till end of section. Returns false if end of file.
  67 bool skipSection(IFstream& inFile)
  68 {
  69     string line;
  70     do
  71     {
  72         inFile.getLine(line);
  73
  74         if (!inFile.good())
  75         {
  76             return false;
  77         }
  78     }
  79     while (line.size() < 4 || line.substr(0, 4) != "$End");
  80
  81     return true;
  82 }
  83
  84
  85 void renumber
  86 (
  87     const Map<label>& mshToFoam,
  88     labelList& labels
  89 )
  90 {
  91     forAll(labels, labelI)
  92     {
  93         labels[labelI] = mshToFoam[labels[labelI]];
  94     }
  95 }
  96
  97
  98 // Find face in pp which uses all vertices in meshF (in mesh point labels)
  99 label findFace(const primitivePatch& pp, const labelList& meshF)
 100 {
 101     const Map<label>& meshPointMap = pp.meshPointMap();
 102
 103     // meshF[0] in pp labels.
 104     if (!meshPointMap.found(meshF[0]))
 105     {
 106         Warning<< "Not using gmsh face " << meshF
 107             << " since zero vertex is not on boundary of polyMesh" << endl;
 108         return -1;
 109     }
 110
 111     // Find faces using first point
 112     const labelList& pFaces = pp.pointFaces()[meshPointMap[meshF[0]]];
 113
 114     // Go through all these faces and check if there is one which uses all of
 115     // meshF vertices (in any order ;-)
 116     forAll(pFaces, i)
 117     {
 118         label faceI = pFaces[i];
 119
 120         const face& f = pp[faceI];
 121
 122         // Count uses of vertices of meshF for f
 123         label nMatched = 0;
 124
 125         forAll(f, fp)
 126         {
 127             if (findIndex(meshF, f[fp]) != -1)
 128             {
 129                 nMatched++;
 130             }
 131         }
 132
 133         if (nMatched == meshF.size())
 134         {
 135             return faceI;
 136         }
 137     }
 138
 139     return -1;
 140 }
 141
 142
 143 // Same but find internal face. Expensive addressing.
 144 label findInternalFace(const primitiveMesh& mesh, const labelList& meshF)
 145 {
 146     const labelList& pFaces = mesh.pointFaces()[meshF[0]];
 147
 148     forAll(pFaces, i)
 149     {
 150         label faceI = pFaces[i];
 151
 152         const face& f = mesh.faces()[faceI];
 153
 154         // Count uses of vertices of meshF for f
 155         label nMatched = 0;
 156
 157         forAll(f, fp)
 158         {
 159             if (findIndex(meshF, f[fp]) != -1)
 160             {
 161                 nMatched++;
 162             }
 163         }
 164
 165         if (nMatched == meshF.size())
 166         {
 167             return faceI;
 168         }
 169     }
 170     return -1;
 171 }
 172
 173
 174 // Determine whether cell is inside-out by checking for any wrong-oriented
 175 // face.
 176 bool correctOrientation(const pointField& points, const cellShape& shape)
 177 {
 178     // Get centre of shape.
 179     point cc(shape.centre(points));
 180
 181     // Get outwards pointing faces.
 182     faceList faces(shape.faces());
 183
 184     forAll(faces, i)
 185     {
 186         const face& f = faces[i];
 187
 188         vector n(f.normal(points));
 189
 190         // Check if vector from any point on face to cc points outwards
 191         if (((points[f[0]] - cc) & n) < 0)
 192         {
 193             // Incorrectly oriented
 194             return false;
 195         }
 196     }
 197
 198     return true;
 199 }
 200
 201
 202 void storeCellInZone
 203 (
 204     const label regPhys,
 205     const label cellI,
 206     Map<label>& physToZone,
 207
 208     labelList& zoneToPhys,
 209     List<DynamicList<label> >& zoneCells
 210 )
 211 {
 212     Map<label>::const_iterator zoneFnd = physToZone.find(regPhys);
 213
 214     if (zoneFnd == physToZone.end())
 215     {
 216         // New region. Allocate zone for it.
 217         label zoneI = zoneCells.size();
 218         zoneCells.setSize(zoneI+1);
 219         zoneToPhys.setSize(zoneI+1);
 220
 221         Info<< "Mapping region " << regPhys << " to Foam cellZone "
 222             << zoneI << endl;
 223         physToZone.insert(regPhys, zoneI);
 224
 225         zoneToPhys[zoneI] = regPhys;
 226         zoneCells[zoneI].append(cellI);
 227     }
 228     else
 229     {
 230         // Existing zone for region
 231         zoneCells[zoneFnd()].append(cellI);
 232     }
 233 }
 234
 235
 236 // Reads points and map
 237 void readPoints(IFstream& inFile, pointField& points, Map<label>& mshToFoam)
 238 {
 239     Info<< "Starting to read points at line " << inFile.lineNumber() << endl;
 240
 241     string line;
 242     inFile.getLine(line);
 243     IStringStream lineStr(line);
 244
 245     label nVerts;
 246     lineStr >> nVerts;
 247
 248     Info<< "Vertices to be read:" << nVerts << endl;
 249
 250     points.setSize(nVerts);
 251     mshToFoam.resize(2*nVerts);
 252
 253     for (label pointI = 0; pointI < nVerts; pointI++)
 254     {
 255         label mshLabel;
 256         scalar xVal, yVal, zVal;
 257
 258         string line;
 259         inFile.getLine(line);
 260         IStringStream lineStr(line);
 261
 262         lineStr >> mshLabel >> xVal >> yVal >> zVal;
 263
 264         point& pt = points[pointI];
 265
 266         pt.x() = xVal;
 267         pt.y() = yVal;
 268         pt.z() = zVal;
 269
 270         mshToFoam.insert(mshLabel, pointI);
 271     }
 272
 273     Info<< "Vertices read:" << mshToFoam.size() << endl;
 274
 275     inFile.getLine(line);
 276     IStringStream tagStr(line);
 277     word tag(tagStr);
 278
 279     if (tag != "$ENDNOD" && tag != "$EndNodes")
 280     {
 281         FatalErrorIn("readPoints(..)")
 282             << "Did not find $ENDNOD tag on line "
 283             << inFile.lineNumber() << exit(FatalError);
 284     }
 285     Info<< endl;
 286 }
 287
 288
 289 // Reads physical names
 290 void readPhysNames(IFstream& inFile, Map<word>& physicalNames)
 291 {
 292     Info<< "Starting to read physical names at line " << inFile.lineNumber()
 293         << endl;
 294
 295     string line;
 296     inFile.getLine(line);
 297     IStringStream lineStr(line);
 298
 299     label nNames;
 300     lineStr >> nNames;
 301
 302     Info<< "Physical names:" << nNames << endl;
 303
 304     physicalNames.resize(nNames);
 305
 306     for (label i = 0; i < nNames; i++)
 307     {
 308         label regionI;
 309         string regionName;
 310
 311         string line;
 312         inFile.getLine(line);
 313         IStringStream lineStr(line);
 314         label nSpaces = lineStr.str().count(' ');
 315
 316         if(nSpaces == 1)
 317         {
 318             lineStr >> regionI >> regionName;
 319
 320             Info<< "    " << regionI << '\t'
 321                 << string::validate<word>(regionName) << endl;
 322         }
 323         else if(nSpaces == 2)
 324         {
 325             // >= Gmsh2.4 physical types has tag in front.
 326             label physType;
 327             lineStr >> physType >> regionI >> regionName;
 328             if (physType == 1)
 329             {
 330                 Info<< "    " << "Line " << regionI << '\t'
 331                     << string::validate<word>(regionName) << endl;
 332             }
 333             else if (physType == 2)
 334             {
 335                 Info<< "    " << "Surface " << regionI << '\t'
 336                     << string::validate<word>(regionName) << endl;
 337             }
 338             else if (physType == 3)
 339             {
 340                 Info<< "    " << "Volume " << regionI << '\t'
 341                     << string::validate<word>(regionName) << endl;
 342             }
 343         }
 344
 345         physicalNames.insert(regionI, string::validate<word>(regionName));
 346     }
 347
 348     inFile.getLine(line);
 349     IStringStream tagStr(line);
 350     word tag(tagStr);
 351
 352     if (tag != "$EndPhysicalNames")
 353     {
 354         FatalErrorIn("readPhysicalNames(..)")
 355             << "Did not find $EndPhysicalNames tag on line "
 356             << inFile.lineNumber() << exit(FatalError);
 357     }
 358     Info<< endl;
 359 }
 360
 361
 362 // Reads cells and patch faces
 363 void readCells
 364 (
 365     const bool version2Format,
 366     const bool keepOrientation,
 367     const pointField& points,
 368     const Map<label>& mshToFoam,
 369     IFstream& inFile,
 370     cellShapeList& cells,
 371
 372     labelList& patchToPhys,
 373     List<DynamicList<face> >& patchFaces,
 374
 375     labelList& zoneToPhys,
 376     List<DynamicList<label> >& zoneCells
 377 )
 378 {
 379     Info<< "Starting to read cells at line " << inFile.lineNumber() << endl;
 380
 381     const cellModel& hex = *(cellModeller::lookup("hex"));
 382     const cellModel& prism = *(cellModeller::lookup("prism"));
 383     const cellModel& pyr = *(cellModeller::lookup("pyr"));
 384     const cellModel& tet = *(cellModeller::lookup("tet"));
 385
 386     face triPoints(3);
 387     face quadPoints(4);
 388     labelList tetPoints(4);
 389     labelList pyrPoints(5);
 390     labelList prismPoints(6);
 391     labelList hexPoints(8);
 392
 393
 394     string line;
 395     inFile.getLine(line);
 396     IStringStream lineStr(line);
 397
 398     label nElems;
 399     lineStr >> nElems;
 400
 401     Info<< "Cells to be read:" << nElems << endl << endl;
 402
 403
 404     // Storage for all cells. Too big. Shrink later
 405     cells.setSize(nElems);
 406
 407     label cellI = 0;
 408     label nTet = 0;
 409     label nPyr = 0;
 410     label nPrism = 0;
 411     label nHex = 0;
 412
 413
 414     // From gmsh physical region to Foam patch
 415     Map<label> physToPatch;
 416
 417     // From gmsh physical region to Foam cellZone
 418     Map<label> physToZone;
 419
 420
 421     for (label elemI = 0; elemI < nElems; elemI++)
 422     {
 423         string line;
 424         inFile.getLine(line);
 425         IStringStream lineStr(line);
 426
 427         label elmNumber, elmType, regPhys;
 428
 429         if (version2Format)
 430         {
 431             lineStr >> elmNumber >> elmType;
 432
 433             label nTags;
 434             lineStr>> nTags;
 435
 436             label regElem, partition;
 437
 438             if (nTags == 3)
 439             {
 440                 lineStr >> regPhys >> regElem >> partition;
 441             }
 442             else if (nTags == 2)
 443             {
 444                 lineStr >> regPhys >> regElem;
 445             }
 446             else
 447             {
 448                 regPhys = 0;
 449                 for (label i = 0; i < nTags; i++)
 450                 {
 451                     label dummy;
 452                     lineStr>> dummy;
 453                 }
 454             }
 455         }
 456         else
 457         {
 458             label regElem, nNodes;
 459             lineStr >> elmNumber >> elmType >> regPhys >> regElem >> nNodes;
 460         }
 461
 462         // regPhys on surface elements is region number.
 463
 464         if (elmType == MSHTRI)
 465         {
 466             lineStr >> triPoints[0] >> triPoints[1] >> triPoints[2];
 467
 468             renumber(mshToFoam, triPoints);
 469
 470             Map<label>::iterator regFnd = physToPatch.find(regPhys);
 471
 472             label patchI = -1;
 473             if (regFnd == physToPatch.end())
 474             {
 475                 // New region. Allocate patch for it.
 476                 patchI = patchFaces.size();
 477
 478                 patchFaces.setSize(patchI + 1);
 479                 patchToPhys.setSize(patchI + 1);
 480
 481                 Info<< "Mapping region " << regPhys << " to Foam patch "
 482                     << patchI << endl;
 483                 physToPatch.insert(regPhys, patchI);
 484                 patchToPhys[patchI] = regPhys;
 485             }
 486             else
 487             {
 488                 // Existing patch for region
 489                 patchI = regFnd();
 490             }
 491
 492             // Add triangle to correct patchFaces.
 493             patchFaces[patchI].append(triPoints);
 494         }
 495         else if (elmType == MSHQUAD)
 496         {
 497             lineStr
 498                 >> quadPoints[0] >> quadPoints[1] >> quadPoints[2]
 499                 >> quadPoints[3];
 500
 501             renumber(mshToFoam, quadPoints);
 502
 503             Map<label>::iterator regFnd = physToPatch.find(regPhys);
 504
 505             label patchI = -1;
 506             if (regFnd == physToPatch.end())
 507             {
 508                 // New region. Allocate patch for it.
 509                 patchI = patchFaces.size();
 510
 511                 patchFaces.setSize(patchI + 1);
 512                 patchToPhys.setSize(patchI + 1);
 513
 514                 Info<< "Mapping region " << regPhys << " to Foam patch "
 515                     << patchI << endl;
 516                 physToPatch.insert(regPhys, patchI);
 517                 patchToPhys[patchI] = regPhys;
 518             }
 519             else
 520             {
 521                 // Existing patch for region
 522                 patchI = regFnd();
 523             }
 524
 525             // Add quad to correct patchFaces.
 526             patchFaces[patchI].append(quadPoints);
 527         }
 528         else if (elmType == MSHTET)
 529         {
 530             storeCellInZone
 531             (
 532                 regPhys,
 533                 cellI,
 534                 physToZone,
 535                 zoneToPhys,
 536                 zoneCells
 537             );
 538
 539             lineStr
 540                 >> tetPoints[0] >> tetPoints[1] >> tetPoints[2]
 541                 >> tetPoints[3];
 542
 543             renumber(mshToFoam, tetPoints);
 544
 545             cells[cellI++] = cellShape(tet, tetPoints);
 546
 547             nTet++;
 548         }
 549         else if (elmType == MSHPYR)
 550         {
 551             storeCellInZone
 552             (
 553                 regPhys,
 554                 cellI,
 555                 physToZone,
 556                 zoneToPhys,
 557                 zoneCells
 558             );
 559
 560             lineStr
 561                 >> pyrPoints[0] >> pyrPoints[1] >> pyrPoints[2]
 562                 >> pyrPoints[3] >> pyrPoints[4];
 563
 564             renumber(mshToFoam, pyrPoints);
 565
 566             cells[cellI++] = cellShape(pyr, pyrPoints);
 567
 568             nPyr++;
 569         }
 570         else if (elmType == MSHPRISM)
 571         {
 572             storeCellInZone
 573             (
 574                 regPhys,
 575                 cellI,
 576                 physToZone,
 577                 zoneToPhys,
 578                 zoneCells
 579             );
 580
 581             lineStr
 582                 >> prismPoints[0] >> prismPoints[1] >> prismPoints[2]
 583                 >> prismPoints[3] >> prismPoints[4] >> prismPoints[5];
 584
 585             renumber(mshToFoam, prismPoints);
 586
 587             cells[cellI] = cellShape(prism, prismPoints);
 588
 589             const cellShape& cell = cells[cellI];
 590
 591             if (!keepOrientation && !correctOrientation(points, cell))
 592             {
 593                 Info<< "Inverting prism " << cellI << endl;
 594                 // Reorder prism.
 595                 prismPoints[0] = cell[0];
 596                 prismPoints[1] = cell[2];
 597                 prismPoints[2] = cell[1];
 598                 prismPoints[3] = cell[3];
 599                 prismPoints[4] = cell[4];
 600                 prismPoints[5] = cell[5];
 601
 602                 cells[cellI] = cellShape(prism, prismPoints);
 603             }
 604
 605             cellI++;
 606
 607             nPrism++;
 608         }
 609         else if (elmType == MSHHEX)
 610         {
 611             storeCellInZone
 612             (
 613                 regPhys,
 614                 cellI,
 615                 physToZone,
 616                 zoneToPhys,
 617                 zoneCells
 618             );
 619
 620             lineStr
 621                 >> hexPoints[0] >> hexPoints[1]
 622                 >> hexPoints[2] >> hexPoints[3]
 623                 >> hexPoints[4] >> hexPoints[5]
 624                 >> hexPoints[6] >> hexPoints[7];
 625
 626             renumber(mshToFoam, hexPoints);
 627
 628             cells[cellI] = cellShape(hex, hexPoints);
 629
 630             const cellShape& cell = cells[cellI];
 631
 632             if (!keepOrientation && !correctOrientation(points, cell))
 633             {
 634                 Info<< "Inverting hex " << cellI << endl;
 635                 // Reorder hex.
 636                 hexPoints[0] = cell[4];
 637                 hexPoints[1] = cell[5];
 638                 hexPoints[2] = cell[6];
 639                 hexPoints[3] = cell[7];
 640                 hexPoints[4] = cell[0];
 641                 hexPoints[5] = cell[1];
 642                 hexPoints[6] = cell[2];
 643                 hexPoints[7] = cell[3];
 644
 645                 cells[cellI] = cellShape(hex, hexPoints);
 646             }
 647
 648             cellI++;
 649
 650             nHex++;
 651         }
 652         else
 653         {
 654             Info<< "Unhandled element " << elmType << " at line "
 655                 << inFile.lineNumber() << endl;
 656         }
 657     }
 658
 659
 660     inFile.getLine(line);
 661     IStringStream tagStr(line);
 662     word tag(tagStr);
 663
 664     if (tag != "$ENDELM" && tag != "$EndElements")
 665     {
 666         FatalErrorIn("readCells(..)")
 667             << "Did not find $ENDELM tag on line "
 668             << inFile.lineNumber() << exit(FatalError);
 669     }
 670
 671
 672     cells.setSize(cellI);
 673
 674     forAll(patchFaces, patchI)
 675     {
 676         patchFaces[patchI].shrink();
 677     }
 678
 679
 680     Info<< "Cells:" << endl
 681     << "    total:" << cells.size() << endl
 682     << "    hex  :" << nHex << endl
 683     << "    prism:" << nPrism << endl
 684     << "    pyr  :" << nPyr << endl
 685     << "    tet  :" << nTet << endl
 686     << endl;
 687
 688     if (cells.size() == 0)
 689     {
 690         FatalErrorIn("readCells(..)")
 691             << "No cells read from file " << inFile.name() << nl
 692             << "Does your file specify any 3D elements (hex=" << MSHHEX
 693             << ", prism=" << MSHPRISM << ", pyramid=" << MSHPYR
 694             << ", tet=" << MSHTET << ")?" << nl
 695             << "Perhaps you have not exported the 3D elements?"
 696             << exit(FatalError);
 697     }
 698
 699     Info<< "CellZones:" << nl
 700         << "Zone\tSize" << endl;
 701
 702     forAll(zoneCells, zoneI)
 703     {
 704         zoneCells[zoneI].shrink();
 705
 706         const labelList& zCells = zoneCells[zoneI];
 707
 708         if (zCells.size())
 709         {
 710             Info<< "    " << zoneI << '\t' << zCells.size() << endl;
 711         }
 712     }
 713     Info<< endl;
 714 }
 715
 716 // Simplified version of the function from createPatch utility.
 717 void removeEmptyPatches(polyMesh& mesh)
 718 {
 719     Info<< "\n";
 720
 721     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 722
 723     DynamicList<polyPatch*> nonEmptyPatches(patches.size());
 724
 725     forAll(patches, idx)
 726     {
 727         const polyPatch& pp = patches[idx];
 728
 729         if (pp.size() > 0)
 730         {
 731             nonEmptyPatches.append
 732             (
 733                 pp.clone
 734                 (
 735                     patches,
 736                     nonEmptyPatches.size(),
 737                     pp.size(),
 738                     pp.start()
 739                 ).ptr()
 740             );
 741         }
 742         else
 743         {
 744             Info<< "Removing empty patch " << pp.name() << endl;
 745         }
 746     }
 747
 748     if (patches.size() != nonEmptyPatches.size())
 749     {
 750         nonEmptyPatches.shrink();
 751         mesh.removeBoundary();
 752         mesh.addPatches(nonEmptyPatches);
 753     }
 754
 755     Info<< "\n";
 756 }
 757
 758
 759 // Main program:
 760
 761 int main(int argc, char *argv[])
 762 {
 763     argList::noParallel();
 764     argList::validArgs.append(".msh file");
 765     argList::validOptions.insert("keepOrientation", "");
 766
 767 #   include "setRootCase.H"
 768 #   include "createTime.H"
 769
 770     fileName mshName(args.additionalArgs()[0]);
 771
 772     bool keepOrientation = args.optionFound("keepOrientation");
 773
 774     // Storage for points
 775     pointField points;
 776     Map<label> mshToFoam;
 777
 778     // Storage for all cells.
 779     cellShapeList cells;
 780
 781     // Map from patch to gmsh physical region
 782     labelList patchToPhys;
 783     // Storage for patch faces.
 784     List<DynamicList<face> > patchFaces(0);
 785
 786     // Map from cellZone to gmsh physical region
 787     labelList zoneToPhys;
 788     // Storage for cell zones.
 789     List<DynamicList<label> > zoneCells(0);
 790
 791     // Name per physical region
 792     Map<word> physicalNames;
 793
 794     // Version 1 or 2 format
 795     bool version2Format = false;
 796
 797
 798     IFstream inFile(mshName);
 799
 800     while (inFile.good())
 801     {
 802         string line;
 803         inFile.getLine(line);
 804         IStringStream lineStr(line);
 805
 806         word tag(lineStr);
 807
 808         if (tag == "$MeshFormat")
 809         {
 810             Info<< "Found $MeshFormat tag; assuming version 2 file format."
 811                 << endl;
 812             version2Format = true;
 813
 814             if (!skipSection(inFile))
 815             {
 816                 break;
 817             }
 818         }
 819         else if (tag == "$PhysicalNames")
 820         {
 821             readPhysNames(inFile, physicalNames);
 822         }
 823         else if (tag == "$NOD" || tag == "$Nodes")
 824         {
 825             readPoints(inFile, points, mshToFoam);
 826         }
 827         else if (tag == "$ELM" || tag == "$Elements")
 828         {
 829             readCells
 830             (
 831                 version2Format,
 832                 keepOrientation,
 833                 points,
 834                 mshToFoam,
 835                 inFile,
 836                 cells,
 837                 patchToPhys,
 838                 patchFaces,
 839                 zoneToPhys,
 840                 zoneCells
 841             );
 842         }
 843         else
 844         {
 845             Info<< "Skipping tag " << tag << " at line "
 846                 << inFile.lineNumber()
 847                 << endl;
 848
 849             if (!skipSection(inFile))
 850             {
 851                 break;
 852             }
 853         }
 854     }
 855
 856
 857     label nValidCellZones = 0;
 858
 859     forAll(zoneCells, zoneI)
 860     {
 861         if (zoneCells[zoneI].size())
 862         {
 863             nValidCellZones++;
 864         }
 865     }
 866
 867
 868     // Problem is that the orientation of the patchFaces does not have to
 869     // be consistent with the outwards orientation of the mesh faces. So
 870     // we have to construct the mesh in two stages:
 871     // 1. define mesh with all boundary faces in one patch
 872     // 2. use the read patchFaces to find the corresponding boundary face
 873     //    and repatch it.
 874
 875
 876     // Create correct number of patches
 877     // (but without any faces in it)
 878     faceListList boundaryFaces(patchFaces.size());
 879
 880     wordList boundaryPatchNames(boundaryFaces.size());
 881
 882     forAll(boundaryPatchNames, patchI)
 883     {
 884         label physReg = patchToPhys[patchI];
 885
 886         Map<word>::const_iterator iter = physicalNames.find(physReg);
 887
 888         if (iter != physicalNames.end())
 889         {
 890             boundaryPatchNames[patchI] = iter();
 891         }
 892         else
 893         {
 894             boundaryPatchNames[patchI] = word("patch") + name(patchI);
 895         }
 896         Info<< "Patch " << patchI << " gets name "
 897             << boundaryPatchNames[patchI] << endl;
 898     }
 899     Info<< endl;
 900
 901     wordList boundaryPatchTypes(boundaryFaces.size(), polyPatch::typeName);
 902     word defaultFacesName = "defaultFaces";
 903     word defaultFacesType = polyPatch::typeName;
 904     wordList boundaryPatchPhysicalTypes
 905     (
 906         boundaryFaces.size(),
 907         polyPatch::typeName
 908     );
 909
 910     polyMesh mesh
 911     (
 912         IOobject
 913         (
 914             polyMesh::defaultRegion,
 915             runTime.constant(),
 916             runTime
 917         ),
 918         xferMove(points),
 919         cells,
 920         boundaryFaces,
 921         boundaryPatchNames,
 922         boundaryPatchTypes,
 923         defaultFacesName,
 924         defaultFacesType,
 925         boundaryPatchPhysicalTypes
 926     );
 927
 928     repatchPolyTopoChanger repatcher(mesh);
 929
 930     // Now use the patchFaces to patch up the outside faces of the mesh.
 931
 932     // Get the patch for all the outside faces (= default patch added as last)
 933     const polyPatch& pp = mesh.boundaryMesh()[mesh.boundaryMesh().size()-1];
 934
 935     // Storage for faceZones.
 936     List<DynamicList<label> > zoneFaces(patchFaces.size());
 937
 938
 939     // Go through all the patchFaces and find corresponding face in pp.
 940     forAll(patchFaces, patchI)
 941     {
 942         const DynamicList<face>& pFaces = patchFaces[patchI];
 943
 944         Info<< "Finding faces of patch " << patchI << endl;
 945
 946         forAll(pFaces, i)
 947         {
 948             const face& f = pFaces[i];
 949
 950             // Find face in pp using all vertices of f.
 951             label patchFaceI = findFace(pp, f);
 952
 953             if (patchFaceI != -1)
 954             {
 955                 label meshFaceI = pp.start() + patchFaceI;
 956
 957                 repatcher.changePatchID(meshFaceI, patchI);
 958             }
 959             else
 960             {
 961                 // Maybe internal face? If so add to faceZone with same index
 962                 // - might be useful.
 963                 label meshFaceI = findInternalFace(mesh, f);
 964
 965                 if (meshFaceI != -1)
 966                 {
 967                     zoneFaces[patchI].append(meshFaceI);
 968                 }
 969                 else
 970                 {
 971                     WarningIn(args.executable())
 972                         << "Could not match gmsh face " << f
 973                         << " to any of the interior or exterior faces"
 974                         << " that share the same 0th point" << endl;
 975                 }
 976             }
 977         }
 978     }
 979     Info<< nl;
 980
 981     // Face zones
 982     label nValidFaceZones = 0;
 983
 984     Info<< "FaceZones:" << nl
 985         << "Zone\tSize" << endl;
 986
 987     forAll(zoneFaces, zoneI)
 988     {
 989         zoneFaces[zoneI].shrink();
 990
 991         const labelList& zFaces = zoneFaces[zoneI];
 992
 993         if (zFaces.size())
 994         {
 995             nValidFaceZones++;
 996
 997             Info<< "    " << zoneI << '\t' << zFaces.size() << endl;
 998         }
 999     }
1000     Info<< endl;
1001
1002
1003     //Get polyMesh to write to constant
1004     runTime.setTime(instant(runTime.constant()), 0);
1005
1006     repatcher.repatch();
1007
1008     List<cellZone*> cz;
1009     List<faceZone*> fz;
1010
1011     // Construct and add the zones. Note that cell ordering does not change
1012     // because of repatch() and neither does internal faces so we can
1013     // use the zoneCells/zoneFaces as is.
1014
1015     if (nValidCellZones > 0)
1016     {
1017         cz.setSize(nValidCellZones);
1018
1019         nValidCellZones = 0;
1020
1021         forAll(zoneCells, zoneI)
1022         {
1023             if (zoneCells[zoneI].size())
1024             {
1025                 label physReg = zoneToPhys[zoneI];
1026
1027                 Map<word>::const_iterator iter = physicalNames.find(physReg);
1028
1029                 word zoneName = "cellZone_" + name(zoneI);
1030                 if (iter != physicalNames.end())
1031                 {
1032                     zoneName = iter();
1033                 }
1034
1035                 Info<< "Writing zone " << zoneI << " to cellZone "
1036                     << zoneName << " and cellSet"
1037                     << endl;
1038
1039                 cellSet cset(mesh, zoneName, labelHashSet(zoneCells[zoneI]));
1040                 cset.write();
1041
1042                 cz[nValidCellZones] = new cellZone
1043                 (
1044                     zoneName,
1045                     zoneCells[zoneI],
1046                     nValidCellZones,
1047                     mesh.cellZones()
1048                 );
1049                 nValidCellZones++;
1050             }
1051         }
1052     }
1053
1054     if (nValidFaceZones > 0)
1055     {
1056         fz.setSize(nValidFaceZones);
1057
1058         nValidFaceZones = 0;
1059
1060         forAll(zoneFaces, zoneI)
1061         {
1062             if (zoneFaces[zoneI].size())
1063             {
1064                 label physReg = zoneToPhys[zoneI];
1065
1066                 Map<word>::const_iterator iter = physicalNames.find(physReg);
1067
1068                 word zoneName = "faceZone_" + name(zoneI);
1069                 if (iter != physicalNames.end())
1070                 {
1071                     zoneName = iter();
1072                 }
1073
1074                 Info<< "Writing zone " << zoneI << " to faceZone "
1075                     << zoneName << " and faceSet"
1076                     << endl;
1077
1078                 faceSet fset(mesh, zoneName, labelHashSet(zoneFaces[zoneI]));
1079                 fset.write();
1080
1081                 fz[nValidFaceZones] = new faceZone
1082                 (
1083                     zoneName,
1084                     zoneFaces[zoneI],
1085                     boolList(zoneFaces[zoneI].size(), true),
1086                     nValidFaceZones,
1087                     mesh.faceZones()
1088                 );
1089                 nValidFaceZones++;
1090             }
1091         }
1092     }
1093
1094     if (cz.size() || fz.size())
1095     {
1096         mesh.addZones(List<pointZone*>(0), fz, cz);
1097     }
1098
1099     removeEmptyPatches(mesh);
1100
1101     mesh.write();
1102
1103     Info<< "End\n" << endl;
1104
1105     return 0;
1106 }
1107
1108
1109 // ************************************************************************* //
1110