applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C

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