applications/utilities/mesh/conversion/fluentMeshToFoam/fluentMeshToFoam.L

   1 /*--------------------------------*- C++ -*----------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     fluentMeshToFoam
  26
  27 Description
  28     Converts a Fluent mesh to FOAM format
  29     including multiple region and region boundary handling.
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 %{
  34
  35 #undef yyFlexLexer
  36
  37 /* ------------------------------------------------------------------------- *\
  38    ------ local definitions
  39 \* ------------------------------------------------------------------------- */
  40
  41 #include "argList.H"
  42 #include "objectRegistry.H"
  43 #include "foamTime.H"
  44 #include "IStringStream.H"
  45 #include "polyMesh.H"
  46 #include "emptyPolyPatch.H"
  47 #include "wallPolyPatch.H"
  48 #include "symmetryPolyPatch.H"
  49 #include "preservePatchTypes.H"
  50 #include "cellShape.H"
  51 #include "faceSet.H"
  52 #include "cellSet.H"
  53 #include "meshTools.H"
  54 #include "OFstream.H"
  55 #include "wordIOList.H"
  56
  57 #include "readHexLabel.H"
  58 #include "cellShapeRecognition.H"
  59 #include "repatchPolyTopoChanger.H"
  60
  61
  62 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  63
  64 using namespace Foam;
  65
  66 const scalar convertToMeters = 1.0;
  67
  68 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  69
  70
  71 label dimensionOfGrid = 0;
  72
  73 label nPoints = 0;
  74 label nFaces = 0;
  75 label nCells = 0;
  76
  77 pointField points(0);
  78 SLList<label> pointGroupZoneID;
  79 SLList<label> pointGroupStartIndex;
  80 SLList<label> pointGroupEndIndex;
  81
  82
  83 faceList faces(0);
  84 labelList owner(0);
  85 labelList neighbour(0);
  86 SLList<label> faceGroupZoneID;
  87 SLList<label> faceGroupStartIndex;
  88 SLList<label> faceGroupEndIndex;
  89
  90 labelList fluentCellModelID(0);
  91 SLList<label> cellGroupZoneID;
  92 SLList<label> cellGroupStartIndex;
  93 SLList<label> cellGroupEndIndex;
  94 SLList<label> cellGroupType;
  95
  96 // number of zones adjusted at run-time if necessary
  97 label maxZoneID = 100;
  98 label zoneIDBuffer = 10;
  99
 100 wordList patchTypeIDs(maxZoneID);
 101 wordList patchNameIDs(maxZoneID);
 102
 103 // Dummy yywrap to keep yylex happy at compile time.
 104 // It is called by yylex but is not used as the mechanism to change file.
 105 // See <<EOF>>
 106 #if YY_FLEX_MINOR_VERSION < 6 && YY_FLEX_SUBMINOR_VERSION < 34
 107 extern "C" int yywrap()
 108 #else
 109 int yyFlexLexer::yywrap()
 110 #endif
 111 {
 112     return 1;
 113 }
 114
 115 %}
 116
 117 one_space                  [ \t\f]
 118 space                      {one_space}*
 119 some_space                 {one_space}+
 120 cspace                     ","{space}
 121 spaceNl                    ({space}|\n|\r)*
 122
 123 alpha                      [_[:alpha:]]
 124 digit                      [[:digit:]]
 125 decDigit                   [[:digit:]]
 126 octalDigit                 [0-7]
 127 hexDigit                   [[:xdigit:]]
 128
 129 lbrac                      "("
 130 rbrac                      ")"
 131 quote                      \"
 132 dash                       "-"
 133 dotColonDash               [.:-]
 134
 135 schemeSpecialInitial       [!$%&*/\\:<=>?~_^#.@']
 136 schemeSpecialSubsequent    [.+-]
 137 schemeSymbol               (({some_space}|{alpha}|{quote}|{schemeSpecialInitial})({alpha}|{quote}|{digit}|{schemeSpecialInitial}|{schemeSpecialSubsequent})*)
 138
 139
 140 identifier                 {alpha}({alpha}|{digit})*
 141 integer                    {decDigit}+
 142 label                      [1-9]{decDigit}*
 143 hexLabel                   {hexDigit}+
 144 zeroLabel                  {digit}*
 145 signedInteger              [-+]?{integer}
 146 word                       ({alpha}|{digit}|{dotColonDash})*
 147
 148 exponent_part              [eE][-+]?{digit}+
 149 fractional_constant        [-+]?(({digit}*"."{digit}+)|({digit}+".")|({digit}))
 150
 151 double                     ((({fractional_constant}{exponent_part}?)|({digit}+{exponent_part}))|0)
 152
 153 x                          {double}
 154 y                          {double}
 155 z                          {double}
 156 scalar                     {double}
 157 labelListElement           {space}{zeroLabel}
 158 hexLabelListElement        {space}{hexLabel}
 159 scalarListElement          {space}{double}
 160 schemeSymbolListElement    {space}{schemeSymbol}
 161 labelList                  ({labelListElement}+{space})
 162 hexLabelList               ({hexLabelListElement}+{space})
 163 scalarList                 ({scalarListElement}+{space})
 164 schemeSymbolList           ({schemeSymbolListElement}+{space})
 165
 166 starStar                   ("**")
 167 text                       ({space}({word}*{space})*)
 168
 169 dateDDMMYYYY               ({digit}{digit}"/"{digit}{digit}"/"{digit}{digit}{digit}{digit})
 170 dateDDMonYYYY              ((({digit}{digit}{space})|({digit}{space})){alpha}*{space}{digit}{digit}{digit}{digit})
 171 time                       ({digit}{digit}":"{digit}{digit}":"{digit}{digit})
 172
 173 versionNumber              ({digit}|".")*
 174
 175 comment                    {spaceNl}"(0"{space}
 176 header                     {spaceNl}"(1"{space}
 177 dimension                  {spaceNl}"(2"{space}
 178 point                      {spaceNl}"(10"{space}
 179 fluentFace                 {spaceNl}"(13"{space}
 180 cell                       {spaceNl}"(12"{space}
 181 zoneVariant1               {spaceNl}"(39"{space}
 182 zoneVariant2               {spaceNl}"(45"{space}
 183
 184 unknownPeriodicFace        {spaceNl}"(17"{space}
 185 periodicFace               {spaceNl}"(18"{space}
 186 cellTree                   {spaceNl}"(58"{space}
 187 faceTree                   {spaceNl}"(59"{space}
 188 faceParents                {spaceNl}"(61"{space}
 189
 190 endOfSection               {space}")"{space}
 191
 192
 193
 194  /* ------------------------------------------------------------------------- *\
 195                       -----  Exclusive start states -----
 196  \* ------------------------------------------------------------------------- */
 197
 198 %option stack
 199
 200 %x readComment
 201 %x embeddedCommentState
 202 %x readHeader
 203 %x readDimension
 204 %x readPoint
 205 %x readPointHeader
 206 %x readNumberOfPoints
 207 %x readPointGroupData
 208 %x readPointData
 209 %x readPoints2D
 210 %x readPoints3D
 211 %x fluentFace
 212 %x readFaceHeader
 213 %x readNumberOfFaces
 214 %x readFaceGroupData
 215 %x readFaceData
 216 %x readFacesMixed
 217 %x readFacesUniform
 218 %x cell
 219 %x readCellHeader
 220 %x readNumberOfCells
 221 %x readCellGroupData
 222 %x readCellData
 223 %x readCellsMixed
 224 %x readCellsUniform
 225 %x zone
 226 %x readZoneHeader
 227 %x readZoneGroupData
 228 %x readZoneData
 229 %x readZoneBlock
 230
 231 %x periodicFace
 232 %x cellTree
 233 %x faceTree
 234 %x faceParents
 235
 236 %x unknownBlock
 237 %x embeddedUnknownBlock
 238 %%
 239
 240 %{
 241     // Point data
 242     label pointGroupNumberOfComponents = 3;
 243     label pointI = 0; // index used for reading points
 244
 245     // Face data
 246     label faceGroupElementType = -1;
 247     label faceI = 0;
 248
 249     // Cell data
 250     label cellGroupElementType = -1;
 251     label celli = 0;
 252 %}
 253
 254
 255  /* ------------------------------------------------------------------------ *\
 256                             ------ Start Lexing ------
 257  \* ------------------------------------------------------------------------ */
 258
 259  /*                     ------ Reading control header ------                 */
 260
 261 {comment} {
 262         yy_push_state(readComment);
 263     }
 264
 265
 266 <readComment>{quote}{text}{quote} {
 267     }
 268
 269
 270 <readComment>{spaceNl}{endOfSection} {
 271         yy_pop_state();
 272     }
 273
 274 {header} {
 275         BEGIN(readHeader);
 276     }
 277
 278 <readHeader>{quote}{text}{quote} {
 279         Info<< "Reading header: " << YYText() << endl;
 280     }
 281
 282
 283 {dimension} {
 284         BEGIN(readDimension);
 285     }
 286
 287 <readDimension>{space}{label}{space} {
 288         IStringStream dimOfGridStream(YYText());
 289
 290         dimensionOfGrid = readLabel(dimOfGridStream);
 291
 292         Info<< "Dimension of grid: " << dimensionOfGrid << endl;
 293     }
 294
 295
 296 {point} {
 297         BEGIN(readPointHeader);
 298     }
 299
 300 <readPointHeader>{spaceNl}{lbrac}{space}"0"{space}"1"{space} {
 301         BEGIN(readNumberOfPoints);
 302     }
 303
 304 <readNumberOfPoints>{space}{hexLabel}{space}{labelList} {
 305
 306         IStringStream numberOfPointsStream(YYText());
 307
 308         nPoints = readHexLabel(numberOfPointsStream);
 309
 310         Info<< "Number of points: " << nPoints << endl;
 311
 312         // set the size of the points list
 313         points.setSize(nPoints);
 314
 315         // meaningless type skipped
 316         readLabel(numberOfPointsStream);
 317
 318         // this dimension of grid may be checked against global dimension
 319         if (numberOfPointsStream)
 320         {
 321             // check dimension of grid
 322             readLabel(numberOfPointsStream);
 323         }
 324         else
 325         {
 326             Info<< endl;
 327         }
 328     }
 329
 330 <readPointHeader>{spaceNl}{lbrac} {
 331         BEGIN(readPointGroupData);
 332     }
 333
 334 <readPointGroupData>{space}{hexLabel}{space}{hexLabel}{space}{hexLabel}{labelList} {
 335         IStringStream pointGroupDataStream(YYText());
 336
 337         // read point zone-ID, start and end-label
 338         // the indices will be used for checking later.
 339         pointGroupZoneID.append(readHexLabel(pointGroupDataStream));
 340
 341         pointGroupStartIndex.append(readHexLabel(pointGroupDataStream));
 342
 343         pointGroupEndIndex.append(readHexLabel(pointGroupDataStream));
 344
 345         // point group type skipped
 346         readHexLabel(pointGroupDataStream);
 347
 348         // In FOAM, indices start from zero - adjust
 349         pointI = pointGroupStartIndex.last() - 1;
 350
 351         // reset number of components to default
 352         pointGroupNumberOfComponents = 3;
 353
 354         // read number of components in the vector
 355         if (pointGroupDataStream)
 356         {
 357             pointGroupNumberOfComponents = readLabel(pointGroupDataStream);
 358         }
 359     }
 360
 361 <readNumberOfPoints,readPointGroupData>{endOfSection} {
 362         BEGIN(readPointData);
 363     }
 364
 365 <readPointData>{spaceNl}{lbrac} {
 366
 367         Info<< "Reading points" << endl;
 368
 369         if (pointGroupNumberOfComponents == 2)
 370         {
 371             yy_push_state(readPoints2D);
 372         }
 373         else
 374         {
 375             yy_push_state(readPoints3D);
 376         }
 377     }
 378
 379 <readPoints2D>{spaceNl}{scalarList} {
 380
 381         IStringStream vertexXyzStream(YYText());
 382
 383         // Note: coordinates must be read one at the time.
 384         scalar x = readScalar(vertexXyzStream);
 385         scalar y = readScalar(vertexXyzStream);
 386
 387         points[pointI] = point(x, y, 0);
 388         pointI++;
 389     }
 390
 391 <readPoints3D>{spaceNl}{scalarList} {
 392
 393         IStringStream vertexXyzStream(YYText());
 394
 395         // Note: coordinates must be read one at the time.
 396         scalar x = readScalar(vertexXyzStream);
 397         scalar y = readScalar(vertexXyzStream);
 398         scalar z = readScalar(vertexXyzStream);
 399
 400         points[pointI] = convertToMeters*point(x, y, z);
 401         pointI++;
 402     }
 403
 404 <readPoints2D,readPoints3D>{spaceNl}{endOfSection} {
 405
 406         // check read of points
 407         if (pointI != pointGroupEndIndex.last())
 408         {
 409             Info<< "problem with reading points: "
 410                 << "start index: " << pointGroupStartIndex.last()
 411                 << " end index: " << pointGroupEndIndex.last()
 412                 << " last points read: " << pointI << endl;
 413         }
 414
 415         yy_pop_state();
 416     }
 417
 418
 419 {fluentFace} {
 420         BEGIN(readFaceHeader);
 421     }
 422
 423 <readFaceHeader>{spaceNl}{lbrac}{space}"0"{space}"1"{space} {
 424         BEGIN(readNumberOfFaces);
 425     }
 426
 427 <readNumberOfFaces>{space}{hexLabel}{space}{labelListElement}+ {
 428
 429         IStringStream numberOfFacesStream(YYText());
 430
 431         nFaces = readHexLabel(numberOfFacesStream);
 432
 433         Info<< "number of faces: " << nFaces << endl;
 434
 435         faces.setSize(nFaces);
 436         owner.setSize(nFaces);
 437         neighbour.setSize(nFaces);
 438
 439         // Meaningless type and element type not read
 440     }
 441
 442 <readFaceHeader>{spaceNl}{lbrac} {
 443         BEGIN(readFaceGroupData);
 444     }
 445
 446 <readFaceGroupData>{space}{hexLabel}{space}{hexLabel}{space}{hexLabel}{hexLabelListElement}+ {
 447
 448         IStringStream faceGroupDataStream(YYText());
 449
 450         // read fluentFace zone-ID, start and end-label
 451         faceGroupZoneID.append(readHexLabel(faceGroupDataStream));
 452
 453         // the indices will be used for checking later.
 454         faceGroupStartIndex.append(readHexLabel(faceGroupDataStream));
 455
 456         faceGroupEndIndex.append(readHexLabel(faceGroupDataStream));
 457
 458         // face group type
 459         readHexLabel(faceGroupDataStream);
 460
 461         faceGroupElementType = readHexLabel(faceGroupDataStream);
 462
 463         // In FOAM, indices start from zero - adjust
 464         faceI = faceGroupStartIndex.last() - 1;
 465     }
 466
 467 <readNumberOfFaces,readFaceGroupData>{spaceNl}{endOfSection} {
 468         BEGIN(readFaceData);
 469     }
 470
 471 <readFaceData>{spaceNl}{lbrac} {
 472
 473         if (faceGroupElementType == 0)
 474         {
 475             Info<< "Reading mixed faces" << endl;
 476             yy_push_state(readFacesMixed);
 477         }
 478         else
 479         {
 480             Info<< "Reading uniform faces" << endl;
 481             yy_push_state(readFacesUniform);
 482         }
 483     }
 484
 485 <readFacesMixed>{spaceNl}{hexLabelList} {
 486
 487         IStringStream mixedFaceStream(YYText());
 488
 489         face& curFaceLabels = faces[faceI];
 490
 491         // set size of label list
 492         curFaceLabels.setSize(readLabel(mixedFaceStream));
 493
 494         forAll (curFaceLabels, i)
 495         {
 496             curFaceLabels[i] = readHexLabel(mixedFaceStream) - 1;
 497         }
 498
 499         // read neighbour and owner. Neighbour comes first
 500         neighbour[faceI] = readHexLabel(mixedFaceStream) - 1;
 501         owner[faceI] = readHexLabel(mixedFaceStream) - 1;
 502         faceI++;
 503     }
 504
 505 <readFacesUniform>{spaceNl}{hexLabelList} {
 506
 507         IStringStream mixedFaceStream(YYText());
 508
 509         face& curFaceLabels = faces[faceI];
 510
 511         // set size of label list. This is OK because in Fluent the type
 512         // for edge is 2, for triangle is 3 and for quad is 4
 513         curFaceLabels.setSize(faceGroupElementType);
 514
 515         forAll (curFaceLabels, i)
 516         {
 517             curFaceLabels[i] = readHexLabel(mixedFaceStream) - 1;
 518         }
 519
 520         // read neighbour and owner. Neighbour comes first
 521         neighbour[faceI] = readHexLabel(mixedFaceStream) - 1;
 522         owner[faceI] = readHexLabel(mixedFaceStream) - 1;
 523         faceI++;
 524     }
 525
 526 <readFacesMixed,readFacesUniform>{spaceNl}{endOfSection} {
 527
 528         // check read of fluentFaces
 529         if (faceI != faceGroupEndIndex.last())
 530         {
 531             Info<< "problem with reading fluentFaces: "
 532                 << "start index: " << faceGroupStartIndex.last()
 533                 << " end index: " << faceGroupEndIndex.last()
 534                 << " last fluentFaces read: " << faceI << endl;
 535         }
 536
 537         yy_pop_state();
 538     }
 539
 540
 541 {cell} {
 542         BEGIN(readCellHeader);
 543     }
 544
 545 <readCellHeader>{spaceNl}{lbrac}{space}"0"{space}"1"{space} {
 546         BEGIN(readNumberOfCells);
 547     }
 548
 549 <readNumberOfCells>{space}{hexLabel}{space}{labelListElement}+ {
 550
 551         IStringStream numberOfCellsStream(YYText());
 552
 553         nCells = readHexLabel(numberOfCellsStream);
 554
 555         Info<< "Number of cells: " << nCells << endl;
 556
 557         fluentCellModelID.setSize(nCells);
 558
 559         // Meaningless type and element type not read
 560     }
 561
 562 <readCellHeader>{spaceNl}{lbrac} {
 563         BEGIN(readCellGroupData);
 564     }
 565
 566 <readCellGroupData>{space}{hexLabel}{space}{hexLabel}{space}{hexLabel}{space}{hexLabel} {
 567         // Warning. This entry must be above the next one because of the lexing
 568         // rules. It is introduced to deal with the problem of reading
 569         // non-standard cell definition from Tgrid, which misses the type label.
 570
 571         Info<< "Tgrid syntax problem: " << YYText() << endl;
 572         IStringStream cellGroupDataStream(YYText());
 573
 574         // read cell zone-ID, start and end-label
 575         cellGroupZoneID.append(readHexLabel(cellGroupDataStream));
 576
 577         // the indices will be used for checking later.
 578         cellGroupStartIndex.append(readHexLabel(cellGroupDataStream));
 579
 580         cellGroupEndIndex.append(readHexLabel(cellGroupDataStream));
 581
 582         cellGroupType.append(readHexLabel(cellGroupDataStream));
 583
 584         Info<< "cellGroupZoneID:" << cellGroupZoneID.last()
 585             << endl;
 586         Info<< "cellGroupStartIndex:" << cellGroupStartIndex.last()
 587             << endl;
 588         Info<< "cellGroupEndIndex:" << cellGroupEndIndex.last()
 589             << endl;
 590         Info<< "cellGroupType:" << cellGroupType.last()
 591             << endl;
 592
 593
 594         // Note. Potentially skip cell set if type is zero.
 595         // This entry does not exist in Tgrid files.
 596         if (dimensionOfGrid == 2)
 597         {
 598             // Tgrid creating triangles
 599             cellGroupElementType = 1;
 600         }
 601         else
 602         {
 603             cellGroupElementType = 2;
 604         }
 605
 606         // In FOAM, indices start from zero - adjust
 607         celli = cellGroupStartIndex.last() - 1;
 608
 609         if (cellGroupElementType != 0)
 610         {
 611             label lastIndex = cellGroupEndIndex.last();
 612
 613             for (; celli < lastIndex; celli++)
 614             {
 615                 fluentCellModelID[celli] = cellGroupElementType;
 616             }
 617         }
 618     }
 619
 620 <readCellGroupData>{space}{hexLabel}{space}{hexLabel}{space}{hexLabel}{space}{hexLabel}{space}{hexLabel} {
 621         // Warning. See above
 622
 623         Info<< "Other readCellGroupData: " << YYText() << endl;
 624
 625
 626         IStringStream cellGroupDataStream(YYText());
 627
 628         // read cell zone-ID, start and end-label
 629         cellGroupZoneID.append(readHexLabel(cellGroupDataStream));
 630
 631         // the indices will be used for checking later.
 632         cellGroupStartIndex.append(readHexLabel(cellGroupDataStream));
 633
 634         cellGroupEndIndex.append(readHexLabel(cellGroupDataStream));
 635
 636         cellGroupType.append(readHexLabel(cellGroupDataStream));
 637
 638         // Note. Potentially skip cell set if type is zero.
 639
 640         cellGroupElementType = readHexLabel(cellGroupDataStream);
 641
 642         // In FOAM, indices start from zero - adjust
 643         celli = cellGroupStartIndex.last() - 1;
 644
 645         if (cellGroupElementType != 0)
 646         {
 647             Info<< "Reading uniform cells" << endl;
 648             label lastIndex = cellGroupEndIndex.last();
 649
 650             for (; celli < lastIndex; celli++)
 651             {
 652                 fluentCellModelID[celli] = cellGroupElementType;
 653             }
 654         }
 655     }
 656
 657 <readNumberOfCells,readCellGroupData>{endOfSection} {
 658         BEGIN(readCellData);
 659     }
 660
 661 <readCellData>{spaceNl}{lbrac} {
 662         Info<< "Reading mixed cells" << endl;
 663         yy_push_state(readCellsMixed);
 664     }
 665
 666 <readCellsMixed>{spaceNl}{labelList} {
 667
 668         IStringStream fluentCellModelIDStream(YYText());
 669
 670         label celliD;
 671         while (fluentCellModelIDStream.read(celliD))
 672         {
 673             fluentCellModelID[celli] = celliD;
 674             celli++;
 675         }
 676     }
 677
 678 <readCellsMixed>{spaceNl}{endOfSection} {
 679
 680         // check read of cells
 681         if (celli != cellGroupEndIndex.last())
 682         {
 683             Info<< "problem with reading cells: "
 684                 << "start index: " << cellGroupStartIndex.last()
 685                 << " end index: " << cellGroupEndIndex.last()
 686                 << " last cells read: " << celli << endl;
 687         }
 688
 689         yy_pop_state();
 690     }
 691
 692
 693
 694 {zoneVariant1} {
 695         BEGIN(readZoneHeader);
 696     }
 697
 698 {zoneVariant2} {
 699         BEGIN(readZoneHeader);
 700     }
 701
 702 <readZoneHeader>{spaceNl}{lbrac} {
 703         BEGIN(readZoneGroupData);
 704     }
 705
 706 <readZoneGroupData>{space}{label}{space}{word}{space}{word} {
 707
 708         IStringStream zoneDataStream(YYText());
 709
 710         // cell zone-ID not in hexadecimal!!! Inconsistency
 711         label zoneID(readLabel(zoneDataStream));
 712
 713         if (zoneID > maxZoneID - 1)
 714         {
 715             // resize the container
 716             maxZoneID = zoneID + zoneIDBuffer;
 717
 718             patchTypeIDs.setSize(maxZoneID);
 719             patchNameIDs.setSize(maxZoneID);
 720         }
 721
 722         zoneDataStream >> patchTypeIDs[zoneID];
 723         zoneDataStream >> patchNameIDs[zoneID];
 724
 725         Info<< "Read zone1:" << zoneID
 726             << " name:" << patchNameIDs[zoneID]
 727             << " patchTypeID:" << patchTypeIDs[zoneID]
 728             << endl;
 729     }
 730
 731 <readZoneGroupData>{space}{label}{space}{word}{space}{word}{space}{label} {
 732         // Fluent manual inconsistency, version 6.1.22
 733         IStringStream zoneDataStream(YYText());
 734
 735         // cell zone-ID not in hexadecimal!!! Inconsistency
 736         label zoneID(readLabel(zoneDataStream));
 737
 738         if (zoneID > maxZoneID - 1)
 739         {
 740             // resize the container
 741             maxZoneID = zoneID + zoneIDBuffer;
 742
 743             patchTypeIDs.setSize(maxZoneID);
 744             patchNameIDs.setSize(maxZoneID);
 745         }
 746
 747         zoneDataStream >> patchTypeIDs[zoneID];
 748         zoneDataStream >> patchNameIDs[zoneID];
 749
 750         Info<< "Read zone2:" << zoneID
 751             << " name:" << patchNameIDs[zoneID]
 752             << " patchTypeID:" << patchTypeIDs[zoneID]
 753             << endl;
 754     }
 755
 756 <readZoneGroupData>{endOfSection} {
 757         BEGIN(readZoneData);
 758     }
 759
 760 <readZoneData>{spaceNl}{lbrac} {
 761         Info<< "Reading zone data" << endl;
 762         yy_push_state(readZoneBlock);
 763     }
 764
 765 <readZoneBlock>{spaceNl}{schemeSymbolList} {
 766     }
 767
 768 <readZoneBlock>{lbrac} {
 769         yy_push_state(unknownBlock);
 770     }
 771
 772 <readZoneBlock>{endOfSection} {
 773         yy_pop_state();
 774     }
 775
 776
 777
 778  /*             ------ Reading end of section and others ------               */
 779
 780 <readHeader,readDimension,readPointData,readFaceData,readCellData,readZoneData>{spaceNl}{endOfSection} {
 781         BEGIN(INITIAL);
 782     }
 783
 784  /*    ------ Reading unknown type or non-standard comment ------             */
 785
 786
 787 {lbrac}{label} {
 788         yy_push_state(unknownBlock);
 789     }
 790
 791 <readComment,unknownBlock,embeddedUnknownBlock>{spaceNl}{schemeSymbol} {
 792     }
 793
 794 <readComment,unknownBlock,embeddedUnknownBlock>{spaceNl}{lbrac} {
 795         yy_push_state(embeddedUnknownBlock);
 796
 797     }
 798
 799 <readComment,unknownBlock,embeddedUnknownBlock>{spaceNl}{endOfSection} {
 800         yy_pop_state();
 801     }
 802
 803 <unknownBlock,embeddedUnknownBlock>{spaceNl}{labelList} {
 804     }
 805
 806 <unknownBlock,embeddedUnknownBlock>{spaceNl}{hexLabelList} {
 807     }
 808
 809 <unknownBlock,embeddedUnknownBlock>{spaceNl}{scalarList} {
 810     }
 811
 812 <unknownBlock,embeddedUnknownBlock>{spaceNl}{schemeSymbolList} {
 813     }
 814
 815 <unknownBlock,embeddedUnknownBlock>{spaceNl}{text} {
 816     }
 817
 818
 819  /* ------ Ignore remaining space and \n s.  Any other characters are errors. */
 820
 821 <readPoints2D,readPoints3D>.|\n {
 822     Info<< "ERROR! Do not understand characters: " << YYText() << endl;
 823     }
 824 .|\n {}
 825
 826
 827  /*  ------ On EOF return to previous file, if none exists terminate. ------  */
 828
 829 <<EOF>> {
 830             yyterminate();
 831     }
 832 %%
 833
 834
 835 #include "fileName.H"
 836 #include <fstream>
 837
 838 // Find label of face.
 839 label findFace(const primitiveMesh& mesh, const face& f)
 840 {
 841     // Get faces using zeroth vertex of face.
 842     const labelList& pFaces = mesh.pointFaces()[f[0] ];
 843
 844     forAll(pFaces, i)
 845     {
 846         label faceI = pFaces[i];
 847
 848         if (f == mesh.faces()[faceI])
 849         {
 850             return faceI;
 851         }
 852     }
 853
 854     // Didn't find face. Do what?
 855     FatalErrorIn("findFace(const primitiveMesh&, const face&)")
 856         << "Problem : cannot find a single face in the mesh which uses"
 857         << " vertices " << f << exit(FatalError);
 858
 859     return -1;
 860 }
 861
 862
 863 int main(int argc, char *argv[])
 864 {
 865     argList::noParallel();
 866     argList::validArgs.append("Fluent mesh file");
 867     argList::validOptions.insert("scale", "scale factor");
 868     argList::validOptions.insert("writeSets", "");
 869     argList::validOptions.insert("writeZones", "");
 870
 871     argList args(argc, argv);
 872
 873     if (!args.check())
 874     {
 875         FatalError.exit();
 876     }
 877
 878     scalar scaleFactor = 1.0;
 879     args.optionReadIfPresent("scale", scaleFactor);
 880
 881     bool writeSets = args.optionFound("writeSets");
 882     bool writeZones = args.optionFound("writeZones");
 883
 884 #   include "createTime.H"
 885
 886     fileName fluentFile(args.additionalArgs()[0]);
 887     std::ifstream fluentStream(fluentFile.c_str());
 888
 889     if (!fluentStream)
 890     {
 891         FatalErrorIn("fluentToFoam::main(int argc, char *argv[])")
 892             << args.executable()
 893             << ": file " << fluentFile << " not found"
 894             << exit(FatalError);
 895     }
 896
 897     //HR 22.11.09: These global variables do not get initialised.
 898     // Debug-version fails!
 899     patchTypeIDs.setSize(maxZoneID);
 900     patchNameIDs.setSize(maxZoneID);
 901
 902     yyFlexLexer lexer(&fluentStream);
 903     while (lexer.yylex() != 0)
 904     {}
 905
 906     Info<< "\n\nFINISHED LEXING\n\n\n";
 907     // Lookup table giving number of vertices given a fluent cell type ID
 908     // Currently not used.
 909     // label fluentModelNVertices[7] = {-1, 3, 4, 4, 8, 5, 6};
 910
 911     // Lookup table giving number of vertices given a fluent cell type ID
 912     label fluentModelNFaces[7] = {-1, 3, 4, 4, 6, 5, 5};
 913
 914     // Make a list of cell faces to be filled in for owner and neighbour
 915
 916     labelListList cellFaces(nCells);
 917
 918     labelList nFacesInCell(nCells, 0);
 919
 920     forAll (cellFaces, celli)
 921     {
 922         cellFaces[celli].setSize(fluentModelNFaces[fluentCellModelID[celli] ]);
 923     }
 924
 925     // fill in owner and neighbour
 926
 927     forAll (owner, faceI)
 928     {
 929         if (owner[faceI] > -1)
 930         {
 931             label curCell = owner[faceI];
 932
 933             if (nFacesInCell[curCell] >= cellFaces[curCell].size())
 934             {
 935                 FatalErrorIn(args.executable())
 936                     << "Trying to add " << nFacesInCell[curCell] + 1
 937                     << "th face to a cell with " << cellFaces[curCell].size()
 938                     << " faces.  Face index: " << faceI
 939                     << " Current faces: " << cellFaces[curCell]
 940                     << abort(FatalError);
 941             }
 942
 943             cellFaces[curCell][nFacesInCell[curCell] ] = faceI;
 944
 945             nFacesInCell[curCell]++;
 946         }
 947     }
 948
 949     forAll (neighbour, faceI)
 950     {
 951         if (neighbour[faceI] > -1)
 952         {
 953             label curCell = neighbour[faceI];
 954
 955             if (nFacesInCell[curCell] >= cellFaces[curCell].size())
 956             {
 957                 FatalErrorIn(args.executable())
 958                     << "Trying to add " << nFacesInCell[curCell] + 1
 959                     << "th face to a cell with " << cellFaces[curCell].size()
 960                     << " faces.  Face index: " << faceI
 961                     << " Current faces: " << cellFaces[curCell]
 962                     << abort(FatalError);
 963             }
 964
 965             cellFaces[curCell][nFacesInCell[curCell] ] = faceI;
 966
 967             nFacesInCell[curCell]++;
 968         }
 969     }
 970
 971     // Construct shapes from face lists
 972     cellShapeList cellShapes(nCells);
 973
 974     // Extrude 2-D mesh into 3-D
 975
 976     Info<< "dimension of grid: " << dimensionOfGrid << endl;
 977     faceList frontAndBackFaces;
 978
 979     if (dimensionOfGrid == 2)
 980     {
 981         const scalar extrusionFactor = 0.01;
 982
 983         boundBox box(max(points), min(points));
 984
 985         const scalar zOffset = extrusionFactor*box.mag();
 986
 987         // two-dimensional grid. Extrude in z-direction
 988         Info<< "Grid is 2-D. Extruding in z-direction by: "
 989             << 2*zOffset << endl;
 990
 991         pointField oldPoints = points;
 992
 993         const label pointOffset = oldPoints.size();
 994
 995         points.setSize(2*pointOffset);
 996
 997         label nNewPoints = 0;
 998
 999         // Note: In order for the owner-neighbour rules to be right, the
1000         // points given by Fluent need to represent the FRONT plane of the
1001         // geometry. Therefore, the extrusion will be in -z direction
1002
1003         forAll (oldPoints, pointI)
1004         {
1005             points[nNewPoints] = oldPoints[pointI];
1006
1007             points[nNewPoints].z() = zOffset;
1008
1009             nNewPoints++;
1010         }
1011
1012         forAll (oldPoints, pointI)
1013         {
1014             points[nNewPoints] = oldPoints[pointI];
1015
1016             points[nNewPoints].z() = -zOffset;
1017
1018             nNewPoints++;
1019         }
1020
1021         // 2-D shape recognition
1022         Info<< "Creating shapes for 2-D cells"<< endl;
1023
1024         // Set the number of empty faces
1025         frontAndBackFaces.setSize(2*nCells);
1026
1027         forAll (fluentCellModelID, celli)
1028         {
1029             switch (fluentCellModelID[celli])
1030             {
1031                 case 1:
1032                 {
1033                     cellShapes[celli] =
1034                         cellShape
1035                         (
1036                             extrudedTriangleCellShape
1037                             (
1038                                 celli,
1039                                 cellFaces[celli],
1040                                 faces,
1041                                 owner,
1042                                 neighbour,
1043                                 pointOffset,
1044                                 frontAndBackFaces
1045                             )
1046                         );
1047                 }
1048                 break;
1049
1050                 case 3:
1051                 {
1052                     cellShapes[celli] =
1053                         cellShape
1054                         (
1055                             extrudedQuadCellShape
1056                             (
1057                                 celli,
1058                                 cellFaces[celli],
1059                                 faces,
1060                                 owner,
1061                                 neighbour,
1062                                 pointOffset,
1063                                 frontAndBackFaces
1064                             )
1065                         );
1066                 }
1067                 break;
1068
1069                 default:
1070                 {
1071                     FatalErrorIn(args.executable())
1072                         << "unrecognised  2-D cell shape: "
1073                         << fluentCellModelID[celli]
1074                         << abort(FatalError);
1075                 }
1076             }
1077         }
1078
1079         // Create new faces
1080         forAll (faces, faceI)
1081         {
1082
1083             if (faces[faceI].size() != 2)
1084             {
1085                 FatalErrorIn(args.executable())
1086                     << "fluentMeshToFoam: a 2-D face defined with "
1087                     << faces[faceI].size() << " points." << endl;
1088             }
1089
1090             labelList& newFace = faces[faceI];
1091
1092             newFace.setSize(4);
1093
1094             newFace[2] = newFace[1] + pointOffset;
1095
1096             newFace[3] = newFace[0] + pointOffset;
1097         }
1098
1099         // Create new cells from 2-D shapes
1100     }
1101     else
1102     {
1103         // 3-D shape recognition
1104         Info<< "Creating shapes for 3-D cells"<< endl;
1105         forAll (fluentCellModelID, celli)
1106         {
1107             if
1108             (
1109                 fluentCellModelID[celli] == 2    // tet
1110              || fluentCellModelID[celli] == 4    // hex
1111              || fluentCellModelID[celli] == 5    // pyramid
1112              || fluentCellModelID[celli] == 6    // prism
1113             )
1114             {
1115                 cellShapes[celli] =
1116                     cellShape
1117                     (
1118                         create3DCellShape
1119                         (
1120                             celli,
1121                             cellFaces[celli],
1122                             faces,
1123                             owner,
1124                             neighbour,
1125                             fluentCellModelID[celli]
1126                         )
1127                     );
1128             }
1129             else
1130             {
1131                 FatalErrorIn(args.executable())
1132                     << "unrecognised 3-D cell shape: "
1133                     << fluentCellModelID[celli]
1134                     << abort(FatalError);
1135             }
1136         }
1137     }
1138
1139     // boundary faces are oriented such that the owner is zero and the face
1140     // area vector points into the domain. Turn them round before making patches
1141     // for Foam compatibility
1142
1143     forAll (faces, faceI)
1144     {
1145         if (owner[faceI] == -1)
1146         {
1147             // reverse face
1148             labelList oldFace = faces[faceI];
1149
1150             forAllReverse(oldFace, i)
1151             {
1152                 faces[faceI][oldFace.size() - i - 1] =
1153                     oldFace[i];
1154             }
1155         }
1156     }
1157
1158
1159     //make patchless mesh before analysing boundaries
1160
1161     faceListList patches(0);
1162     wordList patchNames(0);
1163     wordList patchTypes(0);
1164     word defaultFacesName = "defaultFaces";
1165     word defaultFacesType = emptyPolyPatch::typeName;
1166     wordList patchPhysicalTypes(0);
1167
1168     // Scale the points
1169
1170     points *= scaleFactor;
1171
1172     Info<< "Building patch-less mesh..." << flush;
1173
1174     polyMesh pShapeMesh
1175     (
1176         IOobject
1177         (
1178             polyMesh::defaultRegion,
1179             runTime.constant(),
1180             runTime
1181         ),
1182         xferMove(points),
1183         cellShapes,
1184         patches,
1185         patchNames,
1186         patchTypes,
1187         defaultFacesName,
1188         defaultFacesType,
1189         patchPhysicalTypes
1190     );
1191
1192     //dont write mesh yet, otherwise preservePatchTypes will be broken
1193     //and zones wont be written
1194     //checkmesh done after patch addition as well
1195     Info<< "done." << endl;
1196
1197
1198     Info<< endl << "Building boundary and internal patches." << endl;
1199     //adding patches after mesh construction allows topological checks
1200     //on whether a patch is internal or external, something fluent
1201     //doesnt seem to mind
1202
1203     // Make boundary patches
1204
1205     SLList<label>::iterator faceGroupZoneIDIter = faceGroupZoneID.begin();
1206     SLList<label>::iterator faceGroupStartIndexIter =
1207         faceGroupStartIndex.begin();
1208     SLList<label>::iterator faceGroupEndIndexIter = faceGroupEndIndex.begin();
1209
1210     // Note. Not all groups of faces will be boundary patches.
1211     // Take care on construction
1212
1213     //2D needs extra space for frontAndBack faces
1214     if (dimensionOfGrid == 2)
1215     {
1216         patches.setSize(faceGroupZoneID.size()+1);
1217         patchNames.setSize(faceGroupZoneID.size()+1);
1218         patchTypes.setSize(faceGroupZoneID.size()+1);
1219         patchPhysicalTypes.setSize(faceGroupZoneID.size()+1);
1220     }
1221     else
1222     {
1223         patches.setSize(faceGroupZoneID.size());
1224         patchNames.setSize(faceGroupZoneID.size());
1225         patchTypes.setSize(faceGroupZoneID.size());
1226         patchPhysicalTypes.setSize(faceGroupZoneID.size());
1227     }
1228
1229     label nPatches = 0;
1230
1231     // Colate information for all patches (internal and external)
1232
1233     // Create a file listing patch type for each zone
1234
1235     label maxZoneID = 0;
1236
1237     for
1238     (
1239         SLList<label>::iterator fgzIDIter = faceGroupZoneID.begin();
1240         fgzIDIter != faceGroupZoneID.end();
1241         ++fgzIDIter
1242     )
1243     {
1244         maxZoneID = max(maxZoneID, fgzIDIter());
1245     }
1246
1247     Info << "maxZoneID: " << maxZoneID << endl;
1248     wordIOList zoneToPatchName
1249     (
1250         IOobject
1251         (
1252             "zoneToPatchName",
1253             runTime.constant(),
1254             pShapeMesh.meshSubDir,
1255             pShapeMesh,
1256             IOobject::NO_READ,
1257             IOobject::NO_WRITE
1258         ),
1259         wordList(maxZoneID + 1, "unknown")
1260     );
1261
1262
1263     for
1264     (
1265         ;
1266         faceGroupZoneIDIter != faceGroupZoneID.end()
1267      && faceGroupStartIndexIter != faceGroupStartIndex.end()
1268      && faceGroupEndIndexIter != faceGroupEndIndex.end();
1269         ++faceGroupZoneIDIter,
1270         ++faceGroupStartIndexIter,
1271         ++faceGroupEndIndexIter
1272     )
1273     {
1274         // Get face type and name
1275         const word& curPatchType = patchTypeIDs[faceGroupZoneIDIter()];
1276
1277         const word& curPatchName = patchNameIDs[faceGroupZoneIDIter()];
1278
1279         Info<< "Creating patch " << nPatches
1280             << " for zone: " << faceGroupZoneIDIter()
1281             << " start: " << faceGroupStartIndexIter()
1282             << " end: " << faceGroupEndIndexIter()
1283             << " type: " << curPatchType << " name: " << curPatchName << endl;
1284
1285         // Record zone index
1286         zoneToPatchName[faceGroupZoneIDIter()] = curPatchName;
1287
1288         // Make patch labels
1289         label faceLabel = faceGroupStartIndexIter() - 1;
1290
1291         faceList patchFaces(faceGroupEndIndexIter() - faceLabel);
1292
1293         forAll (patchFaces, faceI)
1294         {
1295             if
1296             (
1297                 faces[faceLabel].size() == 3
1298              || faces[faceLabel].size() == 4
1299             )
1300             {
1301                 patchFaces[faceI] = face(faces[faceLabel]);
1302                 faceLabel++;
1303             }
1304             else
1305             {
1306                 FatalErrorIn("fluentToFoam::main(int argc, char *argv[])")
1307                     << "unrecognised face shape with "
1308                     << patchFaces[faceI].size() << " vertices"
1309                     << abort(FatalError);
1310             }
1311         }
1312
1313         // Inlets and outlets
1314         if
1315         (
1316             curPatchType == "pressure"
1317          || curPatchType == "pressure-inlet"
1318          || curPatchType == "inlet-vent"
1319          || curPatchType == "intake-fan"
1320          || curPatchType == "pressure-outlet"
1321          || curPatchType == "exhaust-fan"
1322          || curPatchType == "outlet-vent"
1323          || curPatchType == "pressure-far-field"
1324          || curPatchType == "velocity-inlet"
1325          || curPatchType == "mass-flow-inlet"
1326          || curPatchType == "outflow"
1327          || curPatchType == "interface"
1328          || curPatchType == "periodic"
1329          || curPatchType == "shadow"
1330         )
1331         {
1332             patches[nPatches] = patchFaces;
1333             patchTypes[nPatches] = polyPatch::typeName;
1334             patchNames[nPatches] = curPatchName;
1335
1336             nPatches++;
1337         }
1338         else if (curPatchType == "wall" ) //wall boundaries
1339         {
1340             patches[nPatches] = patchFaces;
1341             patchTypes[nPatches] = wallPolyPatch::typeName;
1342             patchNames[nPatches] = curPatchName;
1343
1344             nPatches++;
1345         }
1346         else if
1347         (
1348             curPatchType == "symmetry"
1349          || curPatchType == "axis"
1350         ) //symmetry planes
1351         {
1352             patches[nPatches] = patchFaces;
1353             patchTypes[nPatches] = symmetryPolyPatch::typeName;
1354             patchNames[nPatches] = curPatchName;
1355
1356             nPatches++;
1357         }
1358         else if
1359         (
1360             curPatchType == "interior"
1361          || curPatchType == "interface"
1362          || curPatchType == "internal"
1363          || curPatchType == "solid"
1364          || curPatchType == "fan"
1365          || curPatchType == "radiator"
1366          || curPatchType == "porous-jump"
1367         ) //interior boundaries - will not be added as patches
1368         {
1369             patches[nPatches] = patchFaces;
1370             patchTypes[nPatches] = "internal";
1371             patchNames[nPatches] = curPatchName;
1372
1373             nPatches++;
1374         }
1375         else if
1376         (
1377             curPatchType == ""
1378         ) // Unnamed face regions default to interior patches
1379         {
1380             Info<< "Patch " << faceGroupZoneIDIter()
1381                 << ": Faces are defined but "
1382                 << "not created as a zone." << endl
1383                 << "Null specification is only valid for internal faces."
1384                 << endl;
1385
1386             patches[nPatches] = patchFaces;
1387             patchTypes[nPatches] = "internal";
1388             patchNames[nPatches] = curPatchName;
1389
1390             nPatches++;
1391         }
1392         else //unknown face regions are not handled
1393         {
1394             FatalErrorIn(args.executable())
1395                 << "fluent patch type " << curPatchType << " not recognised."
1396                 << abort(FatalError);
1397         }
1398     }
1399
1400     //add front and back boundaries for 2D meshes
1401     if (dimensionOfGrid == 2)
1402     {
1403         Info<< "Creating patch for front and back planes" << endl << endl;
1404
1405         patches[nPatches] = frontAndBackFaces;
1406         patchTypes[nPatches] = emptyPolyPatch::typeName;
1407         patchNames[nPatches] = "frontAndBackPlanes";
1408
1409         nPatches++;
1410     }
1411
1412     //Now have all patch information,
1413     //check whether each patch is internal or external
1414     //and add boundaries to mesh
1415     //also write face sets of all patches
1416     patches.setSize(nPatches);
1417     patchTypes.setSize(nPatches);
1418     patchNames.setSize(nPatches);
1419
1420
1421     //old polyBoundary
1422     const polyBoundaryMesh& oPatches = pShapeMesh.boundaryMesh();
1423     // new patches.
1424     DynamicList<polyPatch*> newPatches(nPatches);
1425
1426     // For every boundary face the old patch.
1427     labelList facePatchID(pShapeMesh.nFaces()-pShapeMesh.nInternalFaces(), -1);
1428     label cMeshFace = pShapeMesh.nInternalFaces();
1429     label nBoundaries = 0;
1430
1431
1432     forAll(patches, patchI)
1433     {
1434         const faceList& bFaces = patches[patchI];
1435
1436         label sz = bFaces.size();
1437         labelList meshFaces(sz,-1);
1438
1439
1440         //make face set and write (seperate from rest for clarity)
1441         //internal and external Fluent boundaries
1442         {
1443             faceSet pFaceSet(pShapeMesh, patchNames[patchI], sz);
1444
1445             forAll (bFaces, j)
1446             {
1447                 const face& f = bFaces[j];
1448                 label cMeshFace = findFace(pShapeMesh, f);
1449                 meshFaces[j] = cMeshFace;
1450                 pFaceSet.insert(cMeshFace);
1451             }
1452             if (writeSets)
1453             {
1454                 Info<< "Writing patch " << patchNames[patchI]
1455                     << " of size " << sz << " to faceSet." << endl;
1456                 pFaceSet.write();
1457             }
1458         }
1459
1460
1461         //check if patch is internal
1462         //also check internal/external-ness of first patch face
1463         //internal faces cannot become foam boundaries
1464         //if a face is defined as internal but is actually external
1465         //it will be put in a default wall boundary
1466         //internal boundaries are simply ignored
1467
1468         if
1469         (
1470             patchTypes[patchI] != "internal"
1471          && !pShapeMesh.isInternalFace(meshFaces[0])
1472         )
1473         {
1474             // First face is external and has valid non-internal type
1475
1476             //check all faces for externalness just to be sure
1477             //and mark patch number to global list
1478             forAll(meshFaces, i)
1479             {
1480                 label faceI = meshFaces[i];
1481
1482                 if (pShapeMesh.isInternalFace(faceI))
1483                 {
1484                     FatalErrorIn(args.executable())
1485                         << "Face " << faceI << " on new patch "
1486                         << patchNames[patchI]
1487                         << " is not an external face of the mesh." << endl
1488                         << exit(FatalError);
1489                 }
1490
1491                 if(facePatchID[faceI - pShapeMesh.nInternalFaces()]!= -1)
1492                 {
1493                     FatalErrorIn(args.executable())
1494                         << "Face " << faceI << " on new patch "
1495                         << patchNames[patchI]
1496                         << " has already been marked for repatching to"
1497                         << " patch "
1498                         << facePatchID[faceI - pShapeMesh.nInternalFaces()]
1499                         << exit(FatalError);
1500                 }
1501                 facePatchID[faceI - pShapeMesh.nInternalFaces()] = nBoundaries;
1502             }
1503
1504             // Add to boundary patch
1505
1506             Info<< "Adding new patch " << patchNames[patchI]
1507                 << " of type " << patchTypes[patchI]
1508                 << " as patch " << nBoundaries << endl;
1509
1510             // Add patch to new patch list
1511             newPatches.append
1512             (
1513                 polyPatch::New
1514                 (
1515                     patchTypes[patchI],
1516                     patchNames[patchI],
1517                     sz,
1518                     cMeshFace,
1519                     nBoundaries,
1520                     oPatches
1521                 ).ptr()
1522             );
1523             nBoundaries++;
1524             cMeshFace += sz;
1525         }
1526         else
1527         {
1528             Info<< "Patch " << patchNames[patchI]
1529                  << " is internal to the mesh "
1530                  << " and is not being added to the boundary."
1531                  << endl;
1532
1533         }
1534     }
1535
1536     // Check for any remaining boundary faces
1537     // and add them to a default wall patch
1538     // this routine should generally not be invoked
1539     {
1540         DynamicList<label> defaultBoundaryFaces(facePatchID.size());
1541         forAll (facePatchID, idI)
1542         {
1543             if(facePatchID[idI] == -1)
1544             {
1545                 defaultBoundaryFaces.append(idI);
1546                 facePatchID[idI] = nBoundaries;
1547             }
1548         }
1549         defaultBoundaryFaces.shrink();
1550
1551         if (defaultBoundaryFaces.size() != 0)
1552         {
1553             Warning << " fluent mesh has " << defaultBoundaryFaces.size()
1554                     << " undefined boundary faces." << endl
1555                     << " Adding undefined faces to new patch `default_wall`"
1556                     << endl;
1557
1558             // Add patch to new patch list
1559
1560             newPatches.append
1561             (
1562                 polyPatch::New
1563                 (
1564                     wallPolyPatch::typeName,
1565                     "default_wall",
1566                     defaultBoundaryFaces.size(),
1567                     cMeshFace,
1568                     nBoundaries,
1569                     oPatches
1570                 ).ptr()
1571             );
1572             nBoundaries++;
1573             cMeshFace += defaultBoundaryFaces.size();
1574         }
1575     }
1576
1577     newPatches.shrink();
1578
1579     // Use facePatchIDs map to reorder boundary faces into compact regions
1580
1581     repatchPolyTopoChanger repatcher(pShapeMesh);
1582
1583     // Add new list of patches
1584     repatcher.changePatches(newPatches);
1585
1586     // Change patch ids
1587     forAll(facePatchID, idI)
1588     {
1589         label faceI = idI + pShapeMesh.nInternalFaces();
1590
1591         repatcher.changePatchID(faceI, facePatchID[idI]);
1592     }
1593     repatcher.repatch();
1594
1595     preservePatchTypes
1596     (
1597         runTime,
1598         runTime.constant(),
1599         polyMesh::defaultRegion,
1600         patchNames,
1601         patchTypes,
1602         defaultFacesName,
1603         defaultFacesType,
1604         patchPhysicalTypes
1605     );
1606
1607     // Set the precision of the points data to 10
1608     IOstream::defaultPrecision(10);
1609
1610
1611     // Zones
1612     // will write out cell zones and internal faces for those zones
1613     // note: zone boundary faces are not added to face zones
1614     // the names of boundaries bordering on cell zones are written to
1615     // a list containing the boundary name and cellzone it borders on
1616     // interior boundaries are handled via faceSets
1617     // cell zones will only be written if there is more than one
1618
1619     if (writeZones && cellGroupZoneID.size() > 1)
1620     {
1621         Info<< "Adding Zones" << endl;
1622         List<pointZone*> pz(0);
1623
1624         label nrCellZones = cellGroupZoneID.size();
1625         List<cellZone*> cz(nrCellZones);
1626
1627         // Make face zones for cell zones
1628         List<faceZone*> fz(nrCellZones);
1629
1630         // List of patch names and the cellZone(s) they border
1631         // this is just an info file to make MRF easier to setup
1632         List<DynamicList<word> > boundaryZones
1633         (
1634             pShapeMesh.boundaryMesh().size()
1635         );
1636
1637         const polyBoundaryMesh& bPatches = pShapeMesh.boundaryMesh();
1638         forAll(bPatches, pI)
1639         {
1640             boundaryZones[pI].append(bPatches[pI].name());
1641         }
1642
1643         label cnt = 0;
1644         SLList<label>::iterator cg = cellGroupZoneID.begin();
1645         SLList<label>::iterator start = cellGroupStartIndex.begin();
1646         SLList<label>::iterator end = cellGroupEndIndex.begin();
1647
1648         for (; cg != cellGroupZoneID.end(); ++cg, ++start, ++end)
1649         {
1650             const word& name = patchNameIDs[cg()];
1651             const word& type = patchTypeIDs[cg()];
1652
1653             Info<< "Writing cell zone: " << name
1654                 << " of type " << type << " starting at " << start() - 1
1655                 << " ending at " << end() - 1 << " to cellSet." << endl;
1656
1657             labelList cls(end() - start() + 1);
1658
1659             // Mark zone cells, used for finding faces
1660             boolList zoneCell(pShapeMesh.nCells(), false);
1661
1662             // Shift cell indices by 1
1663             label nr = 0;
1664             for (label celli = (start() - 1); celli < end(); celli++)
1665             {
1666                 cls[nr] = celli;
1667                 zoneCell[celli] = true;
1668                 nr++;
1669             }
1670
1671             cz[cnt] = new cellZone
1672             (
1673                 name,
1674                 cls,
1675                 cnt,
1676                 pShapeMesh.cellZones()
1677             );
1678
1679             DynamicList<label> zoneFaces(pShapeMesh.nFaces());
1680             forAll (pShapeMesh.faceNeighbour(), faceI)
1681             {
1682                 label nei = pShapeMesh.faceNeighbour()[faceI];
1683                 label own = pShapeMesh.faceOwner()[faceI];
1684                 if(nei != -1)
1685                 {
1686                     if(zoneCell[nei] && zoneCell[own])
1687                     {
1688                         zoneFaces.append(faceI);
1689                     }
1690                 }
1691             }
1692             zoneFaces.shrink();
1693
1694             fz[cnt] = new faceZone
1695             (
1696                 name,
1697                 zoneFaces,
1698                 boolList(zoneFaces.size(), false),
1699                 cnt,
1700                 pShapeMesh.faceZones()
1701             );
1702
1703             // Add cell zones to patch zone list
1704             forAll (bPatches, pI)
1705             {
1706                 const labelList& faceCells = bPatches[pI].faceCells();
1707                 forAll(faceCells, fcI)
1708                 {
1709                     if(zoneCell[faceCells[fcI] ])
1710                     {
1711                         boundaryZones[pI].append(name);
1712                         break;
1713                     }
1714                 }
1715             }
1716
1717             cnt++;
1718         }
1719
1720         pShapeMesh.addZones(pz, fz, cz);
1721
1722         forAll(bPatches, pI)
1723         {
1724             boundaryZones[pI].shrink();
1725         }
1726
1727         fileName bczf
1728         (
1729             runTime.path()/runTime.constant()
1730            /"polyMesh"/"boundaryAdjacentCellZones"
1731         );
1732
1733         OFstream boundaryCellZonesFile(bczf);
1734         forAll (boundaryZones, bzI)
1735         {
1736             forAll(boundaryZones[bzI], bzII)
1737             {
1738                 boundaryCellZonesFile << boundaryZones[bzI][bzII] << " ";
1739             }
1740
1741             boundaryCellZonesFile << endl;
1742         }
1743     }
1744
1745     Info<< endl << "Writing mesh..." << flush;
1746
1747     Info<< " to " << pShapeMesh.instance()/pShapeMesh.meshDir()
1748         << "  " << flush;
1749
1750     pShapeMesh.setInstance(pShapeMesh.instance());
1751     pShapeMesh.write();
1752
1753     // Write zone to patch name
1754     zoneToPatchName.write();
1755     Info<< "done." << endl << endl;
1756
1757     // Write cellSets for Fluent regions
1758     // allows easy post-processing
1759     // set and zone functionality will be integrated some time
1760     // soon negating the need for double output
1761     if (writeSets)
1762     {
1763         if (cellGroupZoneID.size() > 1)
1764         {
1765             Info<< "Writing cell sets" << endl;
1766
1767             SLList<label>::iterator cg = cellGroupZoneID.begin();
1768             SLList<label>::iterator start = cellGroupStartIndex.begin();
1769             SLList<label>::iterator end = cellGroupEndIndex.begin();
1770
1771             // Note: cellGroupXXX are all Fluent indices (starting at 1)
1772             // so offset before using.
1773
1774             for (; cg != cellGroupZoneID.end(); ++cg, ++start, ++end)
1775             {
1776                 const word& name = patchNameIDs[cg()];
1777                 const word& type = patchTypeIDs[cg()];
1778
1779                 Info<< "Writing cell set: " << name
1780                     << " of type " << type << " starting at " << start() - 1
1781                     << " ending at " << end() - 1 << " to cellSet." << endl;
1782
1783                 cellSet internal(pShapeMesh, name, end() - start());
1784
1785                 // shift cell indizes by 1
1786                 for (label celli = start() - 1; celli <= end() - 1; celli++)
1787                 {
1788                     internal.insert(celli);
1789                 }
1790
1791                 internal.write();
1792             }
1793         }
1794         else
1795         {
1796             Info<< "Only one cell group: no set written\n";
1797         }
1798     }
1799
1800     Info<< "\nEnd\n" << endl;
1801     return 0;
1802 }
1803
1804
1805 /* --------------------------------------------------------------------------*\
1806     ------ End of fluentMeshToFoam.L
1807 \* --------------------------------------------------------------------------*/