src/conversion/meshReader/starcd/STARCDMeshReader.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "STARCDMeshReader.H"
  27 #include "cyclicPolyPatch.H"
  28 #include "emptyPolyPatch.H"
  29 #include "wallPolyPatch.H"
  30 #include "symmetryPolyPatch.H"
  31 #include "cellModeller.H"
  32 #include "ListOps.H"
  33 #include "IFstream.H"
  34 #include "IOMap.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 const char* const Foam::meshReaders::STARCD::defaultBoundaryName =
  39     "Default_Boundary_Region";
  40
  41 const char* const Foam::meshReaders::STARCD::defaultSolidBoundaryName =
  42     "Default_Boundary_Solid";
  43
  44 bool Foam::meshReaders::STARCD::keepSolids = false;
  45
  46 const int Foam::meshReaders::STARCD::starToFoamFaceAddr[4][6] =
  47 {
  48     { 4, 5, 2, 3, 0, 1 },     // 11 = pro-STAR hex
  49     { 0, 1, 4, -1, 2, 3 },    // 12 = pro-STAR prism
  50     { 3, -1, 2, -1, 1, 0 },   // 13 = pro-STAR tetra
  51     { 0, -1, 4, 2, 1, 3 }     // 14 = pro-STAR pyramid
  52 };
  53
  54
  55 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  56
  57 void Foam::meshReaders::STARCD::readToNewline(IFstream& is)
  58 {
  59     char ch = '\n';
  60     do
  61     {
  62         (is).get(ch);
  63     }
  64     while ((is) && ch != '\n');
  65 }
  66
  67
  68 bool Foam::meshReaders::STARCD::readHeader(IFstream& is, word fileSignature)
  69 {
  70     if (!is.good())
  71     {
  72         FatalErrorIn("meshReaders::STARCD::readHeader()")
  73             << "cannot read " << fileSignature  << "  " << is.name()
  74             << abort(FatalError);
  75     }
  76
  77     word header;
  78     label majorVersion;
  79
  80     is >> header;
  81     is >> majorVersion;
  82
  83     // skip the rest of the line
  84     readToNewline(is);
  85
  86     // add other checks ...
  87     if (header != fileSignature)
  88     {
  89         Info<< "header mismatch " << fileSignature << "  " << is.name()
  90             << endl;
  91     }
  92
  93     return true;
  94 }
  95
  96
  97 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  98
  99 void Foam::meshReaders::STARCD::readAux(const objectRegistry& registry)
 100 {
 101     boundaryRegion_.readDict(registry);
 102     cellTable_.readDict(registry);
 103 }
 104
 105
 106 // read in the points from the .vrt file
 107 //
 108 /*---------------------------------------------------------------------------*\
 109 Line 1:
 110   PROSTAR_VERTEX [newline]
 111
 112 Line 2:
 113   <version> 0 0 0 0 0 0 0 [newline]
 114
 115 Body:
 116   <vertexId>  <x>  <y>  <z> [newline]
 117
 118 \*---------------------------------------------------------------------------*/
 119 void Foam::meshReaders::STARCD::readPoints
 120 (
 121     const fileName& inputName,
 122     const scalar scaleFactor
 123 )
 124 {
 125     const word fileSignature = "PROSTAR_VERTEX";
 126     label nPoints = 0, maxId = 0;
 127
 128     // Pass 1:
 129     // get # points and maximum vertex label
 130     {
 131         IFstream is(inputName);
 132         readHeader(is, fileSignature);
 133
 134         label lineLabel;
 135         scalar x, y, z;
 136
 137         while ((is >> lineLabel).good())
 138         {
 139             nPoints++;
 140             maxId = max(maxId, lineLabel);
 141             is >> x >> y >> z;
 142         }
 143     }
 144
 145     Info<< "Number of points  = " << nPoints << endl;
 146
 147     // set sizes and reset to invalid values
 148
 149     points_.setSize(nPoints);
 150     mapToFoamPointId_.setSize(maxId+1);
 151
 152     //- original Point number for a given vertex
 153     // might need again in the future
 154     ////     labelList origPointId(nPoints);
 155     ////     origPointId = -1;
 156
 157     mapToFoamPointId_ = -1;
 158
 159     // Pass 2:
 160     // construct pointList and conversion table
 161     // from Star vertex numbers to Foam point labels
 162     if (nPoints > 0)
 163     {
 164         IFstream is(inputName);
 165         readHeader(is, fileSignature);
 166
 167         label lineLabel;
 168
 169         label pointI = 0;
 170         while ((is >> lineLabel).good())
 171         {
 172             is  >> points_[pointI].x()
 173                 >> points_[pointI].y()
 174                 >> points_[pointI].z();
 175
 176             // might need again in the future
 177             ////  origPointId[pointI] = lineLabel;
 178             mapToFoamPointId_[lineLabel] = pointI;
 179             pointI++;
 180         }
 181
 182         if (nPoints > pointI)
 183         {
 184             nPoints = pointI;
 185             points_.setSize(nPoints);
 186             // might need again in the future
 187             //// origPointId.setSize(nPoints);
 188         }
 189
 190         if (scaleFactor > 1.0 + SMALL || scaleFactor < 1.0 - SMALL)
 191         {
 192             points_ *= scaleFactor;
 193         }
 194     }
 195     else
 196     {
 197         FatalErrorIn("meshReaders::STARCD::readPoints()")
 198             << "no points in file " << inputName
 199             << abort(FatalError);
 200     }
 201
 202 }
 203
 204
 205 // read in the cells from the .cel file
 206 //
 207 /*---------------------------------------------------------------------------*\
 208 Line 1:
 209   PROSTAR_CELL [newline]
 210
 211 Line 2:
 212   <version> 0 0 0 0 0 0 0 [newline]
 213
 214 Body:
 215   <cellId>  <shapeId>  <nLabels>  <cellTableId>  <typeId> [newline]
 216   <cellId>  <int1> .. <int8>
 217   <cellId>  <int9> .. <int16>
 218
 219  with shapeId:
 220  *   1 = point
 221  *   2 = line
 222  *   3 = shell
 223  *  11 = hexa
 224  *  12 = prism
 225  *  13 = tetra
 226  *  14 = pyramid
 227  * 255 = polyhedron
 228
 229  with typeId
 230  *   1 = fluid
 231  *   2 = solid
 232  *   3 = baffle
 233  *   4 = shell
 234  *   5 = line
 235  *   6 = point
 236
 237 For primitive cell shapes, the number of vertices will never exceed 8 (hexa)
 238 and corresponds to <nLabels>.
 239 For polyhedral, <nLabels> includes an index table comprising beg/end pairs
 240 for each cell face.
 241
 242 Strictly speaking, we only need the cellModeller for adding boundaries.
 243 \*---------------------------------------------------------------------------*/
 244
 245 void Foam::meshReaders::STARCD::readCells(const fileName& inputName)
 246 {
 247     const word fileSignature = "PROSTAR_CELL";
 248     label nFluids = 0, nSolids = 0, nBaffles = 0, nShells = 0;
 249     label maxId = 0;
 250
 251     bool unknownVertices = false;
 252
 253
 254     // Pass 1:
 255     // count nFluids, nSolids, nBaffle, nShell and maxId
 256     // also see if polyhedral cells were used
 257     {
 258         IFstream is(inputName);
 259         readHeader(is, fileSignature);
 260
 261         label lineLabel, shapeId, nLabels, cellTableId, typeId;
 262
 263         while ((is >> lineLabel).good())
 264         {
 265             label starCellId = lineLabel;
 266             is  >> shapeId
 267                 >> nLabels
 268                 >> cellTableId
 269                 >> typeId;
 270
 271             // skip the rest of the line
 272             readToNewline(is);
 273
 274             // max 8 indices per line
 275             while (nLabels > 0)
 276             {
 277                 readToNewline(is);
 278                 nLabels -= 8;
 279             }
 280
 281             if (typeId == starcdFluidType)
 282             {
 283                 nFluids++;
 284                 maxId = max(maxId, starCellId);
 285
 286                 if (!cellTable_.found(cellTableId))
 287                 {
 288                     cellTable_.setName(cellTableId);
 289                     cellTable_.setMaterial(cellTableId, "fluid");
 290                 }
 291             }
 292             else if (typeId == starcdSolidType)
 293             {
 294                 nSolids++;
 295                 if (keepSolids)
 296                 {
 297                     maxId = max(maxId, starCellId);
 298                 }
 299
 300                 if (!cellTable_.found(cellTableId))
 301                 {
 302                     cellTable_.setName(cellTableId);
 303                     cellTable_.setMaterial(cellTableId, "solid");
 304                 }
 305
 306             }
 307             else if (typeId == starcdBaffleType)
 308             {
 309                 // baffles have no cellTable entry
 310                 nBaffles++;
 311                 maxId = max(maxId, starCellId);
 312             }
 313             else if (typeId == starcdShellType)
 314             {
 315                 nShells++;
 316                 if (!cellTable_.found(cellTableId))
 317                 {
 318                     cellTable_.setName(cellTableId);
 319                     cellTable_.setMaterial(cellTableId, "shell");
 320                 }
 321             }
 322
 323         }
 324     }
 325
 326     Info<< "Number of fluids  = " << nFluids << nl
 327         << "Number of baffles = " << nBaffles << nl;
 328     if (keepSolids)
 329     {
 330         Info<< "Number of solids  = " << nSolids << nl;
 331     }
 332     else
 333     {
 334         Info<< "Ignored   solids  = " << nSolids << nl;
 335     }
 336     Info<< "Ignored   shells  = " << nShells << endl;
 337
 338
 339     label nCells;
 340     if (keepSolids)
 341     {
 342         nCells = nFluids + nSolids;
 343     }
 344     else
 345     {
 346         nCells = nFluids;
 347     }
 348
 349     cellFaces_.setSize(nCells);
 350     cellShapes_.setSize(nCells);
 351     cellTableId_.setSize(nCells);
 352
 353     // information for the interfaces
 354     baffleFaces_.setSize(nBaffles);
 355
 356     // extra space for baffles
 357     origCellId_.setSize(nCells + nBaffles);
 358     mapToFoamCellId_.setSize(maxId+1);
 359     mapToFoamCellId_ = -1;
 360
 361
 362     // avoid undefined shapes for polyhedra
 363     cellShape genericShape(*unknownModel, labelList(0));
 364
 365     // Pass 2:
 366     // construct cellFaces_ and possibly cellShapes_
 367     if (nCells <= 0)
 368     {
 369         FatalErrorIn("meshReaders::STARCD::readCells()")
 370             << "no cells in file " << inputName
 371             << abort(FatalError);
 372     }
 373     else
 374     {
 375         IFstream is(inputName);
 376         readHeader(is, fileSignature);
 377
 378         labelList starLabels(64);
 379         label lineLabel, shapeId, nLabels, cellTableId, typeId;
 380
 381         label cellI = 0;
 382         label baffleI = 0;
 383
 384         while ((is >> lineLabel).good())
 385         {
 386             label starCellId = lineLabel;
 387             is  >> shapeId
 388                 >> nLabels
 389                 >> cellTableId
 390                 >> typeId;
 391
 392             if (nLabels > starLabels.size())
 393             {
 394                 starLabels.setSize(nLabels);
 395             }
 396             starLabels = -1;
 397
 398             // read indices - max 8 per line
 399             for (label i = 0; i < nLabels; ++i)
 400             {
 401                 if ((i % 8) == 0)
 402                 {
 403                     is >> lineLabel;
 404                 }
 405                 is >> starLabels[i];
 406             }
 407
 408             // skip solid cells
 409             if (typeId == starcdSolidType && !keepSolids)
 410             {
 411                 continue;
 412             }
 413
 414             // determine the foam cell shape
 415             const cellModel* curModelPtr = NULL;
 416
 417             // fluid/solid cells
 418             switch (shapeId)
 419             {
 420                 case starcdHex:
 421                     curModelPtr = hexModel;
 422                     break;
 423                 case starcdPrism:
 424                     curModelPtr = prismModel;
 425                     break;
 426                 case starcdTet:
 427                     curModelPtr = tetModel;
 428                     break;
 429                 case starcdPyr:
 430                     curModelPtr = pyrModel;
 431                     break;
 432             }
 433
 434             if (curModelPtr)
 435             {
 436                 // primitive cell - use shapes
 437
 438                 // convert orig vertex Id to point label
 439                 bool isBad = false;
 440                 for (label i=0; i < nLabels; ++i)
 441                 {
 442                     label pointId = mapToFoamPointId_[starLabels[i]];
 443                     if (pointId < 0)
 444                     {
 445                         Info<< "Cells inconsistent with vertex file. "
 446                             << "Star vertex " << starLabels[i]
 447                             << " does not exist" << endl;
 448                         isBad = true;
 449                         unknownVertices = true;
 450                     }
 451                     starLabels[i] = pointId;
 452                 }
 453
 454                 if (isBad)
 455                 {
 456                     continue;
 457                 }
 458
 459                 // record original cell number and lookup
 460                 origCellId_[cellI] = starCellId;
 461                 mapToFoamCellId_[starCellId] = cellI;
 462
 463                 cellTableId_[cellI] = cellTableId;
 464                 cellShapes_[cellI] = cellShape
 465                 (
 466                     *curModelPtr,
 467                     SubList<label>(starLabels, nLabels)
 468                 );
 469
 470                 cellFaces_[cellI] = cellShapes_[cellI].faces();
 471                 cellI++;
 472             }
 473             else if (shapeId == starcdPoly)
 474             {
 475                 // polyhedral cell
 476                 label nFaces = starLabels[0] - 1;
 477
 478                 // convert orig vertex id to point label
 479                 // start with offset (skip the index table)
 480                 bool isBad = false;
 481                 for (label i=starLabels[0]; i < nLabels; ++i)
 482                 {
 483                     label pointId = mapToFoamPointId_[starLabels[i]];
 484                     if (pointId < 0)
 485                     {
 486                         Info<< "Cells inconsistent with vertex file. "
 487                             << "Star vertex " << starLabels[i]
 488                             << " does not exist" << endl;
 489                         isBad = true;
 490                         unknownVertices = true;
 491                     }
 492                     starLabels[i] = pointId;
 493                 }
 494
 495                 if (isBad)
 496                 {
 497                     continue;
 498                 }
 499
 500                 // traverse beg/end indices
 501                 faceList faces(nFaces);
 502                 label faceI = 0;
 503                 for (label i=0; i < nFaces; ++i)
 504                 {
 505                     label beg = starLabels[i];
 506                     label n   = starLabels[i+1] - beg;
 507
 508                     face f
 509                     (
 510                         SubList<label>(starLabels, n, beg)
 511                     );
 512
 513                     f.collapse();
 514
 515                     // valid faces only
 516                     if (f.size() >= 3)
 517                     {
 518                         faces[faceI++] = f;
 519                     }
 520                 }
 521
 522                 if (nFaces > faceI)
 523                 {
 524                     Info<< "star cell " << starCellId << " has "
 525                         << (nFaces - faceI)
 526                         << " empty faces - could cause boundary "
 527                         << "addressing problems"
 528                         << endl;
 529
 530                     nFaces = faceI;
 531                     faces.setSize(nFaces);
 532                 }
 533
 534                 if (nFaces < 4)
 535                 {
 536                     FatalErrorIn("meshReaders::STARCD::readCells()")
 537                         << "star cell " << starCellId << " has " << nFaces
 538                         << abort(FatalError);
 539                 }
 540
 541                 // record original cell number and lookup
 542                 origCellId_[cellI] = starCellId;
 543                 mapToFoamCellId_[starCellId] = cellI;
 544
 545                 cellTableId_[cellI] = cellTableId;
 546                 cellShapes_[cellI]  = genericShape;
 547                 cellFaces_[cellI]   = faces;
 548                 cellI++;
 549             }
 550             else if (typeId == starcdBaffleType)
 551             {
 552                 // baffles
 553
 554                 // convert orig vertex id to point label
 555                 bool isBad = false;
 556                 for (label i=0; i < nLabels; ++i)
 557                 {
 558                     label pointId = mapToFoamPointId_[starLabels[i]];
 559                     if (pointId < 0)
 560                     {
 561                         Info<< "Baffles inconsistent with vertex file. "
 562                             << "Star vertex " << starLabels[i]
 563                             << " does not exist" << endl;
 564                         isBad = true;
 565                         unknownVertices = true;
 566                     }
 567                     starLabels[i] = pointId;
 568                 }
 569
 570                 if (isBad)
 571                 {
 572                     continue;
 573                 }
 574
 575
 576                 face f
 577                 (
 578                     SubList<label>(starLabels, nLabels)
 579                 );
 580
 581                 f.collapse();
 582
 583                 // valid faces only
 584                 if (f.size() >= 3)
 585                 {
 586                     baffleFaces_[baffleI] = f;
 587                     // insert lookup addressing in normal list
 588                     mapToFoamCellId_[starCellId]  = nCells + baffleI;
 589                     origCellId_[nCells + baffleI] = starCellId;
 590                     baffleI++;
 591                 }
 592             }
 593         }
 594
 595         baffleFaces_.setSize(baffleI);
 596     }
 597
 598     if (unknownVertices)
 599     {
 600         FatalErrorIn("meshReaders::STARCD::readCells()")
 601             << "cells with unknown vertices"
 602             << abort(FatalError);
 603     }
 604
 605     // truncate lists
 606
 607 #ifdef DEBUG_READING
 608     Info<< "CELLS READ" << endl;
 609 #endif
 610
 611     // cleanup
 612     mapToFoamPointId_.clear();
 613 }
 614
 615
 616 // read in the boundaries from the .bnd file
 617 //
 618 /*---------------------------------------------------------------------------*\
 619 Line 1:
 620   PROSTAR_BOUNDARY [newline]
 621
 622 Line 2:
 623   <version> 0 0 0 0 0 0 0 [newline]
 624
 625 Body:
 626   <boundId>  <cellId>  <cellFace>  <regionId>  0  <boundaryType> [newline]
 627
 628 where boundaryType is truncated to 4 characters from one of the following:
 629 INLET
 630 PRESSSURE
 631 OUTLET
 632 BAFFLE
 633 etc,
 634 \*---------------------------------------------------------------------------*/
 635
 636 void Foam::meshReaders::STARCD::readBoundary(const fileName& inputName)
 637 {
 638     const word fileSignature = "PROSTAR_BOUNDARY";
 639     label nPatches = 0, nFaces = 0, nBafflePatches = 0, maxId = 0;
 640     label lineLabel, starCellId, cellFaceId, starRegion, configNumber;
 641     word patchType;
 642
 643     labelList mapToFoamPatchId(1000, -1);
 644     labelList nPatchFaces(1000, 0);
 645     labelList origRegion(1000, 0);
 646     patchTypes_.setSize(1000);
 647
 648     // this is what we seem to need
 649     // these MUST correspond to starToFoamFaceAddr
 650     //
 651     Map<label> faceLookupIndex;
 652
 653     faceLookupIndex.insert(hexModel->index(), 0);
 654     faceLookupIndex.insert(prismModel->index(), 1);
 655     faceLookupIndex.insert(tetModel->index(), 2);
 656     faceLookupIndex.insert(pyrModel->index(), 3);
 657
 658     // Pass 1:
 659     // collect
 660     // no. of faces (nFaces), no. of patches (nPatches)
 661     // and for each of these patches the number of faces
 662     // (nPatchFaces[patchLabel])
 663     //
 664     // and a conversion table from Star regions to (Foam) patchLabels
 665     //
 666     // additionally note the no. of baffle patches (nBafflePatches)
 667     // so that we sort these to the end of the patch list
 668     // - this makes it easier to transfer them to an adjacent patch if reqd
 669     {
 670         IFstream is(inputName);
 671
 672         if (is.good())
 673         {
 674             readHeader(is, fileSignature);
 675
 676             while ((is >> lineLabel).good())
 677             {
 678                 nFaces++;
 679                 is  >> starCellId
 680                     >> cellFaceId
 681                     >> starRegion
 682                     >> configNumber
 683                     >> patchType;
 684
 685                 // Build translation table to convert star patch to foam patch
 686                 label patchLabel = mapToFoamPatchId[starRegion];
 687                 if (patchLabel == -1)
 688                 {
 689                     patchLabel = nPatches;
 690                     mapToFoamPatchId[starRegion] = patchLabel;
 691                     origRegion[patchLabel] = starRegion;
 692                     patchTypes_[patchLabel] = patchType;
 693
 694                     maxId = max(maxId, starRegion);
 695
 696                     // should actually be case-insensitive
 697                     if (patchType == "BAFF")
 698                     {
 699                         nBafflePatches++;
 700                     }
 701                     nPatches++;
 702                 }
 703
 704                 nPatchFaces[patchLabel]++;
 705             }
 706
 707             if (nPatches == 0)
 708             {
 709                 Info<< "No boundary faces in file " << inputName << endl;
 710             }
 711         }
 712         else
 713         {
 714             Info<< "Could not read boundary file " << inputName << endl;
 715         }
 716     }
 717
 718     // keep empty patch region in reserve
 719     nPatches++;
 720     Info<< "Number of patches = " << nPatches
 721         << " (including extra for missing)" << endl;
 722
 723     // resize
 724     origRegion.setSize(nPatches);
 725     patchTypes_.setSize(nPatches);
 726     patchNames_.setSize(nPatches);
 727     nPatchFaces.setSize(nPatches);
 728
 729     // add our empty patch
 730     origRegion[nPatches-1] = 0;
 731     nPatchFaces[nPatches-1] = 0;
 732     patchTypes_[nPatches-1] = "none";
 733
 734     // create names
 735     // - use 'Label' entry from "constant/boundaryRegion" dictionary
 736     forAll(patchTypes_, patchI)
 737     {
 738         bool foundName = false, foundType = false;
 739
 740         Map<dictionary>::const_iterator
 741             iter = boundaryRegion_.find(origRegion[patchI]);
 742
 743         if
 744         (
 745             iter != boundaryRegion_.end()
 746         )
 747         {
 748             foundType = iter().readIfPresent
 749             (
 750                 "BoundaryType",
 751                 patchTypes_[patchI]
 752             );
 753
 754             foundName = iter().readIfPresent
 755             (
 756                 "Label",
 757                 patchNames_[patchI]
 758             );
 759         }
 760
 761         // consistent names, in long form and in lowercase
 762         if (!foundType)
 763         {
 764             // transform
 765             forAllIter(string, patchTypes_[patchI], i)
 766             {
 767                 *i = tolower(*i);
 768             }
 769
 770             if (patchTypes_[patchI] == "symp")
 771             {
 772                 patchTypes_[patchI] = "symplane";
 773             }
 774             else if (patchTypes_[patchI] == "cycl")
 775             {
 776                 patchTypes_[patchI] = "cyclic";
 777             }
 778             else if (patchTypes_[patchI] == "baff")
 779             {
 780                 patchTypes_[patchI] = "baffle";
 781             }
 782             else if (patchTypes_[patchI] == "moni")
 783             {
 784                 patchTypes_[patchI] = "monitoring";
 785             }
 786         }
 787
 788         // create a name if needed
 789         if (!foundName)
 790         {
 791             patchNames_[patchI] =
 792                 patchTypes_[patchI] + "_" + name(origRegion[patchI]);
 793         }
 794     }
 795
 796     // enforce name "Default_Boundary_Region"
 797     patchNames_[nPatches-1] = defaultBoundaryName;
 798
 799     // sort according to ascending region numbers, but leave
 800     // Default_Boundary_Region as the final patch
 801     {
 802         labelList sortedIndices;
 803         sortedOrder(SubList<label>(origRegion, nPatches-1), sortedIndices);
 804
 805         labelList oldToNew = identity(nPatches);
 806         forAll(sortedIndices, i)
 807         {
 808             oldToNew[sortedIndices[i]] = i;
 809         }
 810
 811         inplaceReorder(oldToNew, origRegion);
 812         inplaceReorder(oldToNew, patchTypes_);
 813         inplaceReorder(oldToNew, patchNames_);
 814         inplaceReorder(oldToNew, nPatchFaces);
 815     }
 816
 817     // re-sort to have baffles near the end
 818     nBafflePatches = 1;
 819     if (nBafflePatches)
 820     {
 821         labelList oldToNew = identity(nPatches);
 822         label newIndex = 0;
 823         label baffleIndex = (nPatches-1 - nBafflePatches);
 824
 825         for (label i=0; i < oldToNew.size()-1; ++i)
 826         {
 827             if (patchTypes_[i] == "baffle")
 828             {
 829                 oldToNew[i] = baffleIndex++;
 830             }
 831             else
 832             {
 833                 oldToNew[i] = newIndex++;
 834             }
 835         }
 836
 837         inplaceReorder(oldToNew, origRegion);
 838         inplaceReorder(oldToNew, patchTypes_);
 839         inplaceReorder(oldToNew, patchNames_);
 840         inplaceReorder(oldToNew, nPatchFaces);
 841     }
 842
 843     mapToFoamPatchId.setSize(maxId+1, -1);
 844     forAll(origRegion, patchI)
 845     {
 846         mapToFoamPatchId[origRegion[patchI]] = patchI;
 847     }
 848
 849     boundaryIds_.setSize(nPatches);
 850     forAll(boundaryIds_, patchI)
 851     {
 852         boundaryIds_[patchI].setSize(nPatchFaces[patchI]);
 853         nPatchFaces[patchI] = 0;
 854     }
 855
 856
 857     // Pass 2:
 858     //
 859     if (nPatches > 1 && mapToFoamCellId_.size() > 1)
 860     {
 861         IFstream is(inputName);
 862         readHeader(is, fileSignature);
 863
 864         while ((is >> lineLabel).good())
 865         {
 866             is
 867                 >> starCellId
 868                 >> cellFaceId
 869                 >> starRegion
 870                 >> configNumber
 871                 >> patchType;
 872
 873             label patchI = mapToFoamPatchId[starRegion];
 874
 875             // zero-based indexing
 876             cellFaceId--;
 877
 878             label cellId = -1;
 879
 880             // convert to FOAM cell number
 881             if (starCellId < mapToFoamCellId_.size())
 882             {
 883                 cellId = mapToFoamCellId_[starCellId];
 884             }
 885
 886             if (cellId < 0)
 887             {
 888                 Info
 889                     << "Boundaries inconsistent with cell file. "
 890                     << "Star cell " << starCellId << " does not exist"
 891                     << endl;
 892             }
 893             else
 894             {
 895                 // restrict lookup to volume cells (no baffles)
 896                 if (cellId < cellShapes_.size())
 897                 {
 898                     label index = cellShapes_[cellId].model().index();
 899                     if (faceLookupIndex.found(index))
 900                     {
 901                         index = faceLookupIndex[index];
 902                         cellFaceId = starToFoamFaceAddr[index][cellFaceId];
 903                     }
 904                 }
 905                 else
 906                 {
 907                     // we currently use cellId >= nCells to tag baffles,
 908                     // we can also use a negative face number
 909                     cellFaceId = -1;
 910                 }
 911
 912                 boundaryIds_[patchI][nPatchFaces[patchI]] =
 913                     cellFaceIdentifier(cellId, cellFaceId);
 914
 915 #ifdef DEBUG_BOUNDARY
 916                 Info<< "bnd " << cellId << " " << cellFaceId << endl;
 917 #endif
 918                 // increment counter of faces in current patch
 919                 nPatchFaces[patchI]++;
 920             }
 921         }
 922     }
 923
 924     // retain original information in patchPhysicalTypes_ - overwrite latter
 925     patchPhysicalTypes_.setSize(patchTypes_.size());
 926
 927
 928     forAll(boundaryIds_, patchI)
 929     {
 930         // resize - avoid invalid boundaries
 931         if (nPatchFaces[patchI] < boundaryIds_[patchI].size())
 932         {
 933             boundaryIds_[patchI].setSize(nPatchFaces[patchI]);
 934         }
 935
 936         word origType = patchTypes_[patchI];
 937         patchPhysicalTypes_[patchI] = origType;
 938
 939         if (origType == "symplane")
 940         {
 941             patchTypes_[patchI] = symmetryPolyPatch::typeName;
 942             patchPhysicalTypes_[patchI] = patchTypes_[patchI];
 943         }
 944         else if (origType == "wall")
 945         {
 946             patchTypes_[patchI] = wallPolyPatch::typeName;
 947             patchPhysicalTypes_[patchI] = patchTypes_[patchI];
 948         }
 949         else if (origType == "cyclic")
 950         {
 951             // incorrect. should be cyclicPatch but this
 952             // requires info on connected faces.
 953             patchTypes_[patchI] = cyclicPolyPatch::typeName;
 954             patchPhysicalTypes_[patchI] = patchTypes_[patchI];
 955         }
 956         else if (origType == "baffle")
 957         {
 958             // incorrect. tag the patch until we get proper support.
 959             // set physical type to a canonical "baffle"
 960             patchTypes_[patchI] = emptyPolyPatch::typeName;
 961             patchPhysicalTypes_[patchI] = "baffle";
 962         }
 963         else
 964         {
 965             patchTypes_[patchI] = polyPatch::typeName;
 966         }
 967
 968         Info<< "patch " << patchI
 969             << " (region " << origRegion[patchI]
 970             << ": " << origType << ") type: '" << patchTypes_[patchI]
 971             << "' name: " << patchNames_[patchI] << endl;
 972     }
 973
 974     // cleanup
 975     mapToFoamCellId_.clear();
 976     cellShapes_.clear();
 977 }
 978
 979
 980 //
 981 // remove unused points
 982 //
 983 void Foam::meshReaders::STARCD::cullPoints()
 984 {
 985     label nPoints = points_.size();
 986     labelList oldToNew(nPoints, -1);
 987
 988     // loop through cell faces and note which points are being used
 989     forAll(cellFaces_, cellI)
 990     {
 991         const faceList& faces = cellFaces_[cellI];
 992         forAll(faces, i)
 993         {
 994             const labelList& labels = faces[i];
 995             forAll(labels, j)
 996             {
 997                 oldToNew[labels[j]]++;
 998             }
 999         }
1000     }
1001
1002     // the new ordering and the count of unused points
1003     label pointI = 0;
1004     forAll(oldToNew, i)
1005     {
1006         if (oldToNew[i] >= 0)
1007         {
1008             oldToNew[i] = pointI++;
1009         }
1010     }
1011
1012     // report unused points
1013     if (nPoints > pointI)
1014     {
1015         Info<< "Unused    points  = " << (nPoints - pointI) << endl;
1016         nPoints = pointI;
1017
1018         // adjust points and truncate
1019         inplaceReorder(oldToNew, points_);
1020         points_.setSize(nPoints);
1021
1022         // adjust cellFaces - with mesh shapes this might be faster
1023         forAll(cellFaces_, cellI)
1024         {
1025             faceList& faces = cellFaces_[cellI];
1026             forAll(faces, i)
1027             {
1028                 inplaceRenumber(oldToNew, faces[i]);
1029             }
1030         }
1031
1032         // adjust baffles
1033         forAll(baffleFaces_, faceI)
1034         {
1035             inplaceRenumber(oldToNew, baffleFaces_[faceI]);
1036         }
1037     }
1038 }
1039
1040
1041 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1042
1043 bool Foam::meshReaders::STARCD::readGeometry(const scalar scaleFactor)
1044 {
1045     // Info<< "called meshReaders::STARCD::readGeometry" << endl;
1046
1047     readPoints(geometryFile_ + ".vrt", scaleFactor);
1048     readCells(geometryFile_ + ".cel");
1049     cullPoints();
1050     readBoundary(geometryFile_ + ".bnd");
1051
1052     return true;
1053 }
1054
1055
1056 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
1057
1058 Foam::meshReaders::STARCD::STARCD
1059 (
1060     const fileName& prefix,
1061     const objectRegistry& registry,
1062     const scalar scaleFactor
1063 )
1064 :
1065     meshReader(prefix, scaleFactor),
1066     cellShapes_(0),
1067     mapToFoamPointId_(0),
1068     mapToFoamCellId_(0)
1069 {
1070     readAux(registry);
1071 }
1072
1073
1074 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1075
1076 Foam::meshReaders::STARCD::~STARCD()
1077 {}
1078
1079
1080 // ************************************************************************* //