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 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     if (cells.size() == 0)
 661     {
 662         FatalErrorIn("readCells(..)")
 663             << "No cells read from file " << inFile.name() << nl
 664             << "Does your file specify any 3D elements (hex=" << MSHHEX
 665             << ", prism=" << MSHPRISM << ", pyramid=" << MSHPYR
 666             << ", tet=" << MSHTET << ")?" << nl
 667             << "Perhaps you have not exported the 3D elements?"
 668             << exit(FatalError);
 669     }
 670
 671     Info<< "CellZones:" << nl
 672         << "Zone\tSize" << endl;
 673
 674     forAll(zoneCells, zoneI)
 675     {
 676         zoneCells[zoneI].shrink();
 677
 678         const labelList& zCells = zoneCells[zoneI];
 679
 680         if (zCells.size())
 681         {
 682             Info<< "    " << zoneI << '\t' << zCells.size() << endl;
 683         }
 684     }
 685     Info<< endl;
 686 }
 687
 688
 689 // Main program:
 690
 691 int main(int argc, char *argv[])
 692 {
 693     argList::noParallel();
 694     argList::validArgs.append(".msh file");
 695     argList::validOptions.insert("keepOrientation", "");
 696
 697 #   include "setRootCase.H"
 698 #   include "createTime.H"
 699
 700     fileName mshName(args.additionalArgs()[0]);
 701
 702     bool keepOrientation = args.optionFound("keepOrientation");
 703
 704     // Storage for points
 705     pointField points;
 706     Map<label> mshToFoam;
 707
 708     // Storage for all cells.
 709     cellShapeList cells;
 710
 711     // Map from patch to gmsh physical region
 712     labelList patchToPhys;
 713     // Storage for patch faces.
 714     List<DynamicList<face> > patchFaces(0);
 715
 716     // Map from cellZone to gmsh physical region
 717     labelList zoneToPhys;
 718     // Storage for cell zones.
 719     List<DynamicList<label> > zoneCells(0);
 720
 721     // Name per physical region
 722     Map<word> physicalNames;
 723
 724     // Version 1 or 2 format
 725     bool version2Format = false;
 726
 727
 728     IFstream inFile(mshName);
 729
 730     while (inFile.good())
 731     {
 732         string line;
 733         inFile.getLine(line);
 734         IStringStream lineStr(line);
 735
 736         word tag(lineStr);
 737
 738         if (tag == "$MeshFormat")
 739         {
 740             Info<< "Found $MeshFormat tag; assuming version 2 file format."
 741                 << endl;
 742             version2Format = true;
 743
 744             if (!skipSection(inFile))
 745             {
 746                 break;
 747             }
 748         }
 749         else if (tag == "$PhysicalNames")
 750         {
 751             readPhysNames(inFile, physicalNames);
 752         }
 753         else if (tag == "$NOD" || tag == "$Nodes")
 754         {
 755             readPoints(inFile, points, mshToFoam);
 756         }
 757         else if (tag == "$ELM" || tag == "$Elements")
 758         {
 759             readCells
 760             (
 761                 version2Format,
 762                 keepOrientation,
 763                 points,
 764                 mshToFoam,
 765                 inFile,
 766                 cells,
 767                 patchToPhys,
 768                 patchFaces,
 769                 zoneToPhys,
 770                 zoneCells
 771             );
 772         }
 773         else
 774         {
 775             Info<< "Skipping tag " << tag << " at line "
 776                 << inFile.lineNumber()
 777                 << endl;
 778
 779             if (!skipSection(inFile))
 780             {
 781                 break;
 782             }
 783         }
 784     }
 785
 786
 787     label nValidCellZones = 0;
 788
 789     forAll(zoneCells, zoneI)
 790     {
 791         if (zoneCells[zoneI].size())
 792         {
 793             nValidCellZones++;
 794         }
 795     }
 796
 797
 798     // Problem is that the orientation of the patchFaces does not have to
 799     // be consistent with the outwards orientation of the mesh faces. So
 800     // we have to construct the mesh in two stages:
 801     // 1. define mesh with all boundary faces in one patch
 802     // 2. use the read patchFaces to find the corresponding boundary face
 803     //    and repatch it.
 804
 805
 806     // Create correct number of patches
 807     // (but without any faces in it)
 808     faceListList boundaryFaces(patchFaces.size());
 809
 810     wordList boundaryPatchNames(boundaryFaces.size());
 811
 812     forAll(boundaryPatchNames, patchI)
 813     {
 814         label physReg = patchToPhys[patchI];
 815
 816         Map<word>::const_iterator iter = physicalNames.find(physReg);
 817
 818         if (iter != physicalNames.end())
 819         {
 820             boundaryPatchNames[patchI] = iter();
 821         }
 822         else
 823         {
 824             boundaryPatchNames[patchI] = word("patch") + name(patchI);
 825         }
 826         Info<< "Patch " << patchI << " gets name "
 827             << boundaryPatchNames[patchI] << endl;
 828     }
 829     Info<< endl;
 830
 831     wordList boundaryPatchTypes(boundaryFaces.size(), polyPatch::typeName);
 832     word defaultFacesName = "defaultFaces";
 833     word defaultFacesType = polyPatch::typeName;
 834     wordList boundaryPatchPhysicalTypes
 835     (
 836         boundaryFaces.size(),
 837         polyPatch::typeName
 838     );
 839
 840     polyMesh mesh
 841     (
 842         IOobject
 843         (
 844             polyMesh::defaultRegion,
 845             runTime.constant(),
 846             runTime
 847         ),
 848         xferMove(points),
 849         cells,
 850         boundaryFaces,
 851         boundaryPatchNames,
 852         boundaryPatchTypes,
 853         defaultFacesName,
 854         defaultFacesType,
 855         boundaryPatchPhysicalTypes
 856     );
 857
 858     repatchPolyTopoChanger repatcher(mesh);
 859
 860     // Now use the patchFaces to patch up the outside faces of the mesh.
 861
 862     // Get the patch for all the outside faces (= default patch added as last)
 863     const polyPatch& pp = mesh.boundaryMesh()[mesh.boundaryMesh().size()-1];
 864
 865     // Storage for faceZones.
 866     List<DynamicList<label> > zoneFaces(patchFaces.size());
 867
 868
 869     // Go through all the patchFaces and find corresponding face in pp.
 870     forAll(patchFaces, patchI)
 871     {
 872         const DynamicList<face>& pFaces = patchFaces[patchI];
 873
 874         Info<< "Finding faces of patch " << patchI << endl;
 875
 876         forAll(pFaces, i)
 877         {
 878             const face& f = pFaces[i];
 879
 880             // Find face in pp using all vertices of f.
 881             label patchFaceI = findFace(pp, f);
 882
 883             if (patchFaceI != -1)
 884             {
 885                 label meshFaceI = pp.start() + patchFaceI;
 886
 887                 repatcher.changePatchID(meshFaceI, patchI);
 888             }
 889             else
 890             {
 891                 // Maybe internal face? If so add to faceZone with same index
 892                 // - might be useful.
 893                 label meshFaceI = findInternalFace(mesh, f);
 894
 895                 if (meshFaceI != -1)
 896                 {
 897                     zoneFaces[patchI].append(meshFaceI);
 898                 }
 899                 else
 900                 {
 901                     WarningIn(args.executable())
 902                         << "Could not match gmsh face " << f
 903                         << " to any of the interior or exterior faces"
 904                         << " that share the same 0th point" << endl;
 905                 }
 906             }
 907         }
 908     }
 909     Info<< nl;
 910
 911     // Face zones
 912     label nValidFaceZones = 0;
 913
 914     Info<< "FaceZones:" << nl
 915         << "Zone\tSize" << endl;
 916
 917     forAll(zoneFaces, zoneI)
 918     {
 919         zoneFaces[zoneI].shrink();
 920
 921         const labelList& zFaces = zoneFaces[zoneI];
 922
 923         if (zFaces.size())
 924         {
 925             nValidFaceZones++;
 926
 927             Info<< "    " << zoneI << '\t' << zFaces.size() << endl;
 928         }
 929     }
 930     Info<< endl;
 931
 932
 933     //Get polyMesh to write to constant
 934     runTime.setTime(instant(runTime.constant()), 0);
 935
 936     repatcher.repatch();
 937
 938     List<cellZone*> cz;
 939     List<faceZone*> fz;
 940
 941     // Construct and add the zones. Note that cell ordering does not change
 942     // because of repatch() and neither does internal faces so we can
 943     // use the zoneCells/zoneFaces as is.
 944
 945     if (nValidCellZones > 0)
 946     {
 947         cz.setSize(nValidCellZones);
 948
 949         nValidCellZones = 0;
 950
 951         forAll(zoneCells, zoneI)
 952         {
 953             if (zoneCells[zoneI].size())
 954             {
 955                 label physReg = zoneToPhys[zoneI];
 956
 957                 Map<word>::const_iterator iter = physicalNames.find(physReg);
 958
 959                 word zoneName = "cellZone_" + name(zoneI);
 960                 if (iter != physicalNames.end())
 961                 {
 962                     zoneName = iter();
 963                 }
 964
 965                 Info<< "Writing zone " << zoneI << " to cellZone "
 966                     << zoneName << " and cellSet"
 967                     << endl;
 968
 969                 cellSet cset(mesh, zoneName, zoneCells[zoneI]);
 970                 cset.write();
 971
 972                 cz[nValidCellZones] = new cellZone
 973                 (
 974                     zoneName,
 975                     zoneCells[zoneI],
 976                     nValidCellZones,
 977                     mesh.cellZones()
 978                 );
 979                 nValidCellZones++;
 980             }
 981         }
 982     }
 983
 984     if (nValidFaceZones > 0)
 985     {
 986         fz.setSize(nValidFaceZones);
 987
 988         nValidFaceZones = 0;
 989
 990         forAll(zoneFaces, zoneI)
 991         {
 992             if (zoneFaces[zoneI].size())
 993             {
 994                 label physReg = zoneToPhys[zoneI];
 995
 996                 Map<word>::const_iterator iter = physicalNames.find(physReg);
 997
 998                 word zoneName = "faceZone_" + name(zoneI);
 999                 if (iter != physicalNames.end())
1000                 {
1001                     zoneName = iter();
1002                 }
1003
1004                 Info<< "Writing zone " << zoneI << " to faceZone "
1005                     << zoneName << " and faceSet"
1006                     << endl;
1007
1008                 faceSet fset(mesh, zoneName, zoneFaces[zoneI]);
1009                 fset.write();
1010
1011                 fz[nValidFaceZones] = new faceZone
1012                 (
1013                     zoneName,
1014                     zoneFaces[zoneI],
1015                     boolList(zoneFaces[zoneI].size(), true),
1016                     nValidFaceZones,
1017                     mesh.faceZones()
1018                 );
1019                 nValidFaceZones++;
1020             }
1021         }
1022     }
1023
1024     if (cz.size() || fz.size())
1025     {
1026         mesh.addZones(List<pointZone*>(0), fz, cz);
1027     }
1028
1029     mesh.write();
1030
1031     Info<< "End\n" << endl;
1032
1033     return 0;
1034 }
1035
1036
1037 // ************************************************************************* //
1038