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