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) 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 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 "Time.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 mesh format
 237 scalar readMeshFormat(IFstream& inFile)
 238 {
 239     Info<< "Starting to read mesh format at line "
 240         << inFile.lineNumber()
 241         << endl;
 242
 243     string line;
 244     inFile.getLine(line);
 245     IStringStream lineStr(line);
 246
 247     scalar version;
 248     label asciiFlag, nBytes;
 249     lineStr >> version >> asciiFlag >> nBytes;
 250
 251     Info<< "Read format version " << version << "  ascii " << asciiFlag << endl;
 252
 253     if (asciiFlag != 0)
 254     {
 255         FatalIOErrorIn("readMeshFormat(IFstream&)", inFile)
 256             << "Can only read ascii msh files."
 257             << exit(FatalIOError);
 258     }
 259
 260     inFile.getLine(line);
 261     IStringStream tagStr(line);
 262     word tag(tagStr);
 263
 264     if (tag != "$EndMeshFormat")
 265     {
 266         FatalIOErrorIn("readMeshFormat(IFstream&)", inFile)
 267             << "Did not find $ENDNOD tag on line "
 268             << inFile.lineNumber() << exit(FatalIOError);
 269     }
 270     Info<< endl;
 271
 272     return version;
 273 }
 274
 275
 276 // Reads points and map
 277 void readPoints(IFstream& inFile, pointField& points, Map<label>& mshToFoam)
 278 {
 279     Info<< "Starting to read points at line " << inFile.lineNumber() << endl;
 280
 281     string line;
 282     inFile.getLine(line);
 283     IStringStream lineStr(line);
 284
 285     label nVerts;
 286     lineStr >> nVerts;
 287
 288     Info<< "Vertices to be read:" << nVerts << endl;
 289
 290     points.setSize(nVerts);
 291     mshToFoam.resize(2*nVerts);
 292
 293     for (label pointI = 0; pointI < nVerts; pointI++)
 294     {
 295         label mshLabel;
 296         scalar xVal, yVal, zVal;
 297
 298         string line;
 299         inFile.getLine(line);
 300         IStringStream lineStr(line);
 301
 302         lineStr >> mshLabel >> xVal >> yVal >> zVal;
 303
 304         point& pt = points[pointI];
 305
 306         pt.x() = xVal;
 307         pt.y() = yVal;
 308         pt.z() = zVal;
 309
 310         mshToFoam.insert(mshLabel, pointI);
 311     }
 312
 313     Info<< "Vertices read:" << mshToFoam.size() << endl;
 314
 315     inFile.getLine(line);
 316     IStringStream tagStr(line);
 317     word tag(tagStr);
 318
 319     if (tag != "$ENDNOD" && tag != "$EndNodes")
 320     {
 321         FatalIOErrorIn("readPoints(..)", inFile)
 322             << "Did not find $ENDNOD tag on line "
 323             << inFile.lineNumber() << exit(FatalIOError);
 324     }
 325     Info<< endl;
 326 }
 327
 328
 329 // Reads physical names
 330 void readPhysNames(IFstream& inFile, Map<word>& physicalNames)
 331 {
 332     Info<< "Starting to read physical names at line " << inFile.lineNumber()
 333         << endl;
 334
 335     string line;
 336     inFile.getLine(line);
 337     IStringStream lineStr(line);
 338
 339     label nNames;
 340     lineStr >> nNames;
 341
 342     Info<< "Physical names:" << nNames << endl;
 343
 344     physicalNames.resize(nNames);
 345
 346     for (label i = 0; i < nNames; i++)
 347     {
 348         label regionI;
 349         string regionName;
 350
 351         string line;
 352         inFile.getLine(line);
 353         IStringStream lineStr(line);
 354         label nSpaces = lineStr.str().count(' ');
 355
 356         if (nSpaces == 1)
 357         {
 358             lineStr >> regionI >> regionName;
 359
 360             Info<< "    " << regionI << '\t'
 361                 << string::validate<word>(regionName) << endl;
 362         }
 363         else if (nSpaces == 2)
 364         {
 365             // >= Gmsh2.4 physical types has tag in front.
 366             label physType;
 367             lineStr >> physType >> regionI >> regionName;
 368             if (physType == 1)
 369             {
 370                 Info<< "    " << "Line " << regionI << '\t'
 371                     << string::validate<word>(regionName) << endl;
 372             }
 373             else if (physType == 2)
 374             {
 375                 Info<< "    " << "Surface " << regionI << '\t'
 376                     << string::validate<word>(regionName) << endl;
 377             }
 378             else if (physType == 3)
 379             {
 380                 Info<< "    " << "Volume " << regionI << '\t'
 381                     << string::validate<word>(regionName) << endl;
 382             }
 383         }
 384
 385         physicalNames.insert(regionI, string::validate<word>(regionName));
 386     }
 387
 388     inFile.getLine(line);
 389     IStringStream tagStr(line);
 390     word tag(tagStr);
 391
 392     if (tag != "$EndPhysicalNames")
 393     {
 394         FatalIOErrorIn("readPhysicalNames(..)", inFile)
 395             << "Did not find $EndPhysicalNames tag on line "
 396             << inFile.lineNumber() << exit(FatalIOError);
 397     }
 398     Info<< endl;
 399 }
 400
 401
 402 // Reads cells and patch faces
 403 void readCells
 404 (
 405     const scalar versionFormat,
 406     const bool keepOrientation,
 407     const pointField& points,
 408     const Map<label>& mshToFoam,
 409     IFstream& inFile,
 410     cellShapeList& cells,
 411
 412     labelList& patchToPhys,
 413     List<DynamicList<face> >& patchFaces,
 414
 415     labelList& zoneToPhys,
 416     List<DynamicList<label> >& zoneCells
 417 )
 418 {
 419     Info<< "Starting to read cells at line " << inFile.lineNumber() << endl;
 420
 421     const cellModel& hex = *(cellModeller::lookup("hex"));
 422     const cellModel& prism = *(cellModeller::lookup("prism"));
 423     const cellModel& pyr = *(cellModeller::lookup("pyr"));
 424     const cellModel& tet = *(cellModeller::lookup("tet"));
 425
 426     face triPoints(3);
 427     face quadPoints(4);
 428     labelList tetPoints(4);
 429     labelList pyrPoints(5);
 430     labelList prismPoints(6);
 431     labelList hexPoints(8);
 432
 433
 434     string line;
 435     inFile.getLine(line);
 436     IStringStream lineStr(line);
 437
 438     label nElems;
 439     lineStr >> nElems;
 440
 441     Info<< "Cells to be read:" << nElems << endl << endl;
 442
 443
 444     // Storage for all cells. Too big. Shrink later
 445     cells.setSize(nElems);
 446
 447     label cellI = 0;
 448     label nTet = 0;
 449     label nPyr = 0;
 450     label nPrism = 0;
 451     label nHex = 0;
 452
 453
 454     // From gmsh physical region to Foam patch
 455     Map<label> physToPatch;
 456
 457     // From gmsh physical region to Foam cellZone
 458     Map<label> physToZone;
 459
 460
 461     for (label elemI = 0; elemI < nElems; elemI++)
 462     {
 463         string line;
 464         inFile.getLine(line);
 465         IStringStream lineStr(line);
 466
 467         label elmNumber, elmType, regPhys;
 468
 469         if (versionFormat >= 2)
 470         {
 471             lineStr >> elmNumber >> elmType;
 472
 473             label nTags;
 474             lineStr>> nTags;
 475
 476             if (nTags > 0)
 477             {
 478                 // Assume the first tag is the physical surface
 479                 lineStr >> regPhys;
 480                 for (label i = 1; i < nTags; i++)
 481                 {
 482                     label dummy;
 483                     lineStr>> dummy;
 484                 }
 485             }
 486         }
 487         else
 488         {
 489             label regElem, nNodes;
 490             lineStr >> elmNumber >> elmType >> regPhys >> regElem >> nNodes;
 491         }
 492
 493         // regPhys on surface elements is region number.
 494
 495         if (elmType == MSHTRI)
 496         {
 497             lineStr >> triPoints[0] >> triPoints[1] >> triPoints[2];
 498
 499             renumber(mshToFoam, triPoints);
 500
 501             Map<label>::iterator regFnd = physToPatch.find(regPhys);
 502
 503             label patchI = -1;
 504             if (regFnd == physToPatch.end())
 505             {
 506                 // New region. Allocate patch for it.
 507                 patchI = patchFaces.size();
 508
 509                 patchFaces.setSize(patchI + 1);
 510                 patchToPhys.setSize(patchI + 1);
 511
 512                 Info<< "Mapping region " << regPhys << " to Foam patch "
 513                     << patchI << endl;
 514                 physToPatch.insert(regPhys, patchI);
 515                 patchToPhys[patchI] = regPhys;
 516             }
 517             else
 518             {
 519                 // Existing patch for region
 520                 patchI = regFnd();
 521             }
 522
 523             // Add triangle to correct patchFaces.
 524             patchFaces[patchI].append(triPoints);
 525         }
 526         else if (elmType == MSHQUAD)
 527         {
 528             lineStr
 529                 >> quadPoints[0] >> quadPoints[1] >> quadPoints[2]
 530                 >> quadPoints[3];
 531
 532             renumber(mshToFoam, quadPoints);
 533
 534             Map<label>::iterator regFnd = physToPatch.find(regPhys);
 535
 536             label patchI = -1;
 537             if (regFnd == physToPatch.end())
 538             {
 539                 // New region. Allocate patch for it.
 540                 patchI = patchFaces.size();
 541
 542                 patchFaces.setSize(patchI + 1);
 543                 patchToPhys.setSize(patchI + 1);
 544
 545                 Info<< "Mapping region " << regPhys << " to Foam patch "
 546                     << patchI << endl;
 547                 physToPatch.insert(regPhys, patchI);
 548                 patchToPhys[patchI] = regPhys;
 549             }
 550             else
 551             {
 552                 // Existing patch for region
 553                 patchI = regFnd();
 554             }
 555
 556             // Add quad to correct patchFaces.
 557             patchFaces[patchI].append(quadPoints);
 558         }
 559         else if (elmType == MSHTET)
 560         {
 561             storeCellInZone
 562             (
 563                 regPhys,
 564                 cellI,
 565                 physToZone,
 566                 zoneToPhys,
 567                 zoneCells
 568             );
 569
 570             lineStr
 571                 >> tetPoints[0] >> tetPoints[1] >> tetPoints[2]
 572                 >> tetPoints[3];
 573
 574             renumber(mshToFoam, tetPoints);
 575
 576             cells[cellI++] = cellShape(tet, tetPoints);
 577
 578             nTet++;
 579         }
 580         else if (elmType == MSHPYR)
 581         {
 582             storeCellInZone
 583             (
 584                 regPhys,
 585                 cellI,
 586                 physToZone,
 587                 zoneToPhys,
 588                 zoneCells
 589             );
 590
 591             lineStr
 592                 >> pyrPoints[0] >> pyrPoints[1] >> pyrPoints[2]
 593                 >> pyrPoints[3] >> pyrPoints[4];
 594
 595             renumber(mshToFoam, pyrPoints);
 596
 597             cells[cellI++] = cellShape(pyr, pyrPoints);
 598
 599             nPyr++;
 600         }
 601         else if (elmType == MSHPRISM)
 602         {
 603             storeCellInZone
 604             (
 605                 regPhys,
 606                 cellI,
 607                 physToZone,
 608                 zoneToPhys,
 609                 zoneCells
 610             );
 611
 612             lineStr
 613                 >> prismPoints[0] >> prismPoints[1] >> prismPoints[2]
 614                 >> prismPoints[3] >> prismPoints[4] >> prismPoints[5];
 615
 616             renumber(mshToFoam, prismPoints);
 617
 618             cells[cellI] = cellShape(prism, prismPoints);
 619
 620             const cellShape& cell = cells[cellI];
 621
 622             if (!keepOrientation && !correctOrientation(points, cell))
 623             {
 624                 Info<< "Inverting prism " << cellI << endl;
 625                 // Reorder prism.
 626                 prismPoints[0] = cell[0];
 627                 prismPoints[1] = cell[2];
 628                 prismPoints[2] = cell[1];
 629                 prismPoints[3] = cell[3];
 630                 prismPoints[4] = cell[4];
 631                 prismPoints[5] = cell[5];
 632
 633                 cells[cellI] = cellShape(prism, prismPoints);
 634             }
 635
 636             cellI++;
 637
 638             nPrism++;
 639         }
 640         else if (elmType == MSHHEX)
 641         {
 642             storeCellInZone
 643             (
 644                 regPhys,
 645                 cellI,
 646                 physToZone,
 647                 zoneToPhys,
 648                 zoneCells
 649             );
 650
 651             lineStr
 652                 >> hexPoints[0] >> hexPoints[1]
 653                 >> hexPoints[2] >> hexPoints[3]
 654                 >> hexPoints[4] >> hexPoints[5]
 655                 >> hexPoints[6] >> hexPoints[7];
 656
 657             renumber(mshToFoam, hexPoints);
 658
 659             cells[cellI] = cellShape(hex, hexPoints);
 660
 661             const cellShape& cell = cells[cellI];
 662
 663             if (!keepOrientation && !correctOrientation(points, cell))
 664             {
 665                 Info<< "Inverting hex " << cellI << endl;
 666                 // Reorder hex.
 667                 hexPoints[0] = cell[4];
 668                 hexPoints[1] = cell[5];
 669                 hexPoints[2] = cell[6];
 670                 hexPoints[3] = cell[7];
 671                 hexPoints[4] = cell[0];
 672                 hexPoints[5] = cell[1];
 673                 hexPoints[6] = cell[2];
 674                 hexPoints[7] = cell[3];
 675
 676                 cells[cellI] = cellShape(hex, hexPoints);
 677             }
 678
 679             cellI++;
 680
 681             nHex++;
 682         }
 683         else
 684         {
 685             Info<< "Unhandled element " << elmType << " at line "
 686                 << inFile.lineNumber() << endl;
 687         }
 688     }
 689
 690
 691     inFile.getLine(line);
 692     IStringStream tagStr(line);
 693     word tag(tagStr);
 694
 695     if (tag != "$ENDELM" && tag != "$EndElements")
 696     {
 697         FatalIOErrorIn("readCells(..)", inFile)
 698             << "Did not find $ENDELM tag on line "
 699             << inFile.lineNumber() << exit(FatalIOError);
 700     }
 701
 702
 703     cells.setSize(cellI);
 704
 705     forAll(patchFaces, patchI)
 706     {
 707         patchFaces[patchI].shrink();
 708     }
 709
 710
 711     Info<< "Cells:" << endl
 712     << "    total:" << cells.size() << endl
 713     << "    hex  :" << nHex << endl
 714     << "    prism:" << nPrism << endl
 715     << "    pyr  :" << nPyr << endl
 716     << "    tet  :" << nTet << endl
 717     << endl;
 718
 719     if (cells.size() == 0)
 720     {
 721         FatalIOErrorIn("readCells(..)", inFile)
 722             << "No cells read from file " << inFile.name() << nl
 723             << "Does your file specify any 3D elements (hex=" << MSHHEX
 724             << ", prism=" << MSHPRISM << ", pyramid=" << MSHPYR
 725             << ", tet=" << MSHTET << ")?" << nl
 726             << "Perhaps you have not exported the 3D elements?"
 727             << exit(FatalIOError);
 728     }
 729
 730     Info<< "CellZones:" << nl
 731         << "Zone\tSize" << endl;
 732
 733     forAll(zoneCells, zoneI)
 734     {
 735         zoneCells[zoneI].shrink();
 736
 737         const labelList& zCells = zoneCells[zoneI];
 738
 739         if (zCells.size())
 740         {
 741             Info<< "    " << zoneI << '\t' << zCells.size() << endl;
 742         }
 743     }
 744     Info<< endl;
 745 }
 746
 747
 748 // Main program:
 749
 750 int main(int argc, char *argv[])
 751 {
 752     argList::noParallel();
 753     argList::validArgs.append(".msh file");
 754     argList::addBoolOption
 755     (
 756         "keepOrientation",
 757         "retain raw orientation for prisms/hexs"
 758     );
 759
 760 #   include "addRegionOption.H"
 761
 762 #   include "setRootCase.H"
 763 #   include "createTime.H"
 764
 765     Foam::word regionName;
 766
 767     if (args.optionReadIfPresent("region", regionName))
 768     {
 769         Foam::Info
 770             << "Creating polyMesh for region " << regionName << endl;
 771     }
 772     else
 773     {
 774         regionName = Foam::polyMesh::defaultRegion;
 775     }
 776
 777     const bool keepOrientation = args.optionFound("keepOrientation");
 778     IFstream inFile(args[1]);
 779
 780     // Storage for points
 781     pointField points;
 782     Map<label> mshToFoam;
 783
 784     // Storage for all cells.
 785     cellShapeList cells;
 786
 787     // Map from patch to gmsh physical region
 788     labelList patchToPhys;
 789     // Storage for patch faces.
 790     List<DynamicList<face> > patchFaces(0);
 791
 792     // Map from cellZone to gmsh physical region
 793     labelList zoneToPhys;
 794     // Storage for cell zones.
 795     List<DynamicList<label> > zoneCells(0);
 796
 797     // Name per physical region
 798     Map<word> physicalNames;
 799
 800     // Version 1 or 2 format
 801     scalar versionFormat = 1;
 802
 803
 804     while (inFile.good())
 805     {
 806         string line;
 807         inFile.getLine(line);
 808         IStringStream lineStr(line);
 809
 810         word tag(lineStr);
 811
 812         if (tag == "$MeshFormat")
 813         {
 814             versionFormat = readMeshFormat(inFile);
 815         }
 816         else if (tag == "$PhysicalNames")
 817         {
 818             readPhysNames(inFile, physicalNames);
 819         }
 820         else if (tag == "$NOD" || tag == "$Nodes")
 821         {
 822             readPoints(inFile, points, mshToFoam);
 823         }
 824         else if (tag == "$ELM" || tag == "$Elements")
 825         {
 826             readCells
 827             (
 828                 versionFormat,
 829                 keepOrientation,
 830                 points,
 831                 mshToFoam,
 832                 inFile,
 833                 cells,
 834                 patchToPhys,
 835                 patchFaces,
 836                 zoneToPhys,
 837                 zoneCells
 838             );
 839         }
 840         else
 841         {
 842             Info<< "Skipping tag " << tag << " at line "
 843                 << inFile.lineNumber()
 844                 << endl;
 845
 846             if (!skipSection(inFile))
 847             {
 848                 break;
 849             }
 850         }
 851     }
 852
 853
 854     label nValidCellZones = 0;
 855
 856     forAll(zoneCells, zoneI)
 857     {
 858         if (zoneCells[zoneI].size())
 859         {
 860             nValidCellZones++;
 861         }
 862     }
 863
 864
 865     // Problem is that the orientation of the patchFaces does not have to
 866     // be consistent with the outwards orientation of the mesh faces. So
 867     // we have to construct the mesh in two stages:
 868     // 1. define mesh with all boundary faces in one patch
 869     // 2. use the read patchFaces to find the corresponding boundary face
 870     //    and repatch it.
 871
 872
 873     // Create correct number of patches
 874     // (but without any faces in it)
 875     faceListList boundaryFaces(patchFaces.size());
 876
 877     wordList boundaryPatchNames(boundaryFaces.size());
 878
 879     forAll(boundaryPatchNames, patchI)
 880     {
 881         label physReg = patchToPhys[patchI];
 882
 883         Map<word>::const_iterator iter = physicalNames.find(physReg);
 884
 885         if (iter != physicalNames.end())
 886         {
 887             boundaryPatchNames[patchI] = iter();
 888         }
 889         else
 890         {
 891             boundaryPatchNames[patchI] = word("patch") + name(patchI);
 892         }
 893         Info<< "Patch " << patchI << " gets name "
 894             << boundaryPatchNames[patchI] << endl;
 895     }
 896     Info<< endl;
 897
 898     wordList boundaryPatchTypes(boundaryFaces.size(), polyPatch::typeName);
 899     word defaultFacesName = "defaultFaces";
 900     word defaultFacesType = polyPatch::typeName;
 901     wordList boundaryPatchPhysicalTypes
 902     (
 903         boundaryFaces.size(),
 904         polyPatch::typeName
 905     );
 906
 907     polyMesh mesh
 908     (
 909         IOobject
 910         (
 911             regionName,
 912             runTime.constant(),
 913             runTime
 914         ),
 915         xferMove(points),
 916         cells,
 917         boundaryFaces,
 918         boundaryPatchNames,
 919         boundaryPatchTypes,
 920         defaultFacesName,
 921         defaultFacesType,
 922         boundaryPatchPhysicalTypes
 923     );
 924
 925     repatchPolyTopoChanger repatcher(mesh);
 926
 927     // Now use the patchFaces to patch up the outside faces of the mesh.
 928
 929     // Get the patch for all the outside faces (= default patch added as last)
 930     const polyPatch& pp = mesh.boundaryMesh().last();
 931
 932     // Storage for faceZones.
 933     List<DynamicList<label> > zoneFaces(patchFaces.size());
 934
 935
 936     // Go through all the patchFaces and find corresponding face in pp.
 937     forAll(patchFaces, patchI)
 938     {
 939         const DynamicList<face>& pFaces = patchFaces[patchI];
 940
 941         Info<< "Finding faces of patch " << patchI << endl;
 942
 943         forAll(pFaces, i)
 944         {
 945             const face& f = pFaces[i];
 946
 947             // Find face in pp using all vertices of f.
 948             label patchFaceI = findFace(pp, f);
 949
 950             if (patchFaceI != -1)
 951             {
 952                 label meshFaceI = pp.start() + patchFaceI;
 953
 954                 repatcher.changePatchID(meshFaceI, patchI);
 955             }
 956             else
 957             {
 958                 // Maybe internal face? If so add to faceZone with same index
 959                 // - might be useful.
 960                 label meshFaceI = findInternalFace(mesh, f);
 961
 962                 if (meshFaceI != -1)
 963                 {
 964                     zoneFaces[patchI].append(meshFaceI);
 965                 }
 966                 else
 967                 {
 968                     WarningIn(args.executable())
 969                         << "Could not match gmsh face " << f
 970                         << " to any of the interior or exterior faces"
 971                         << " that share the same 0th point" << endl;
 972                 }
 973             }
 974         }
 975     }
 976     Info<< nl;
 977
 978     // Face zones
 979     label nValidFaceZones = 0;
 980
 981     Info<< "FaceZones:" << nl
 982         << "Zone\tSize" << endl;
 983
 984     forAll(zoneFaces, zoneI)
 985     {
 986         zoneFaces[zoneI].shrink();
 987
 988         const labelList& zFaces = zoneFaces[zoneI];
 989
 990         if (zFaces.size())
 991         {
 992             nValidFaceZones++;
 993
 994             Info<< "    " << zoneI << '\t' << zFaces.size() << endl;
 995         }
 996     }
 997     Info<< endl;
 998
 999
1000     //Get polyMesh to write to constant
1001
1002     runTime.setTime(instant(runTime.constant()), 0);
1003
1004     repatcher.repatch();
1005
1006     List<cellZone*> cz;
1007     List<faceZone*> fz;
1008
1009     // Construct and add the zones. Note that cell ordering does not change
1010     // because of repatch() and neither does internal faces so we can
1011     // use the zoneCells/zoneFaces as is.
1012
1013     if (nValidCellZones > 0)
1014     {
1015         cz.setSize(nValidCellZones);
1016
1017         nValidCellZones = 0;
1018
1019         forAll(zoneCells, zoneI)
1020         {
1021             if (zoneCells[zoneI].size())
1022             {
1023                 label physReg = zoneToPhys[zoneI];
1024
1025                 Map<word>::const_iterator iter = physicalNames.find(physReg);
1026
1027                 word zoneName = "cellZone_" + name(zoneI);
1028                 if (iter != physicalNames.end())
1029                 {
1030                     zoneName = iter();
1031                 }
1032
1033                 Info<< "Writing zone " << zoneI << " to cellZone "
1034                     << zoneName << " and cellSet"
1035                     << endl;
1036
1037                 cellSet cset(mesh, zoneName, zoneCells[zoneI]);
1038                 cset.write();
1039
1040                 cz[nValidCellZones] = new cellZone
1041                 (
1042                     zoneName,
1043                     zoneCells[zoneI],
1044                     nValidCellZones,
1045                     mesh.cellZones()
1046                 );
1047                 nValidCellZones++;
1048             }
1049         }
1050     }
1051
1052     if (nValidFaceZones > 0)
1053     {
1054         fz.setSize(nValidFaceZones);
1055
1056         nValidFaceZones = 0;
1057
1058         forAll(zoneFaces, zoneI)
1059         {
1060             if (zoneFaces[zoneI].size())
1061             {
1062                 label physReg = zoneToPhys[zoneI];
1063
1064                 Map<word>::const_iterator iter = physicalNames.find(physReg);
1065
1066                 word zoneName = "faceZone_" + name(zoneI);
1067                 if (iter != physicalNames.end())
1068                 {
1069                     zoneName = iter();
1070                 }
1071
1072                 Info<< "Writing zone " << zoneI << " to faceZone "
1073                     << zoneName << " and faceSet"
1074                     << endl;
1075
1076                 faceSet fset(mesh, zoneName, zoneFaces[zoneI]);
1077                 fset.write();
1078
1079                 fz[nValidFaceZones] = new faceZone
1080                 (
1081                     zoneName,
1082                     zoneFaces[zoneI],
1083                     boolList(zoneFaces[zoneI].size(), true),
1084                     nValidFaceZones,
1085                     mesh.faceZones()
1086                 );
1087                 nValidFaceZones++;
1088             }
1089         }
1090     }
1091
1092     if (cz.size() || fz.size())
1093     {
1094         mesh.addZones(List<pointZone*>(0), fz, cz);
1095     }
1096
1097     // Remove empty defaultFaces
1098     label defaultPatchID = mesh.boundaryMesh().findPatchID(defaultFacesName);
1099     if (mesh.boundaryMesh()[defaultPatchID].size() == 0)
1100     {
1101         List<polyPatch*> newPatchPtrList((mesh.boundaryMesh().size() - 1));
1102         label newPatchI = 0;
1103         forAll(mesh.boundaryMesh(), patchI)
1104         {
1105             if (patchI != defaultPatchID)
1106             {
1107                 const polyPatch& patch = mesh.boundaryMesh()[patchI];
1108
1109                 newPatchPtrList[newPatchI] = patch.clone
1110                 (
1111                     mesh.boundaryMesh(),
1112                     newPatchI,
1113                     patch.size(),
1114                     patch.start()
1115                 ).ptr();
1116
1117                 newPatchI++;
1118             }
1119         }
1120         repatcher.changePatches(newPatchPtrList);
1121     }
1122
1123     mesh.write();
1124
1125     Info<< "End\n" << endl;
1126
1127     return 0;
1128 }
1129
1130
1131 // ************************************************************************* //