src/foam/meshes/polyMesh/polyMeshFromShapeMesh.C

   1 /*---------------------------------------------------------------------------*\
   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 Description
  25     Create polyMesh from cell and patch shapes
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "polyMesh.H"
  30 #include "foamTime.H"
  31 #include "primitiveMesh.H"
  32 #include "DynamicList.H"
  33
  34 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  35
  36 Foam::labelListList Foam::polyMesh::cellShapePointCells
  37 (
  38     const cellShapeList& c
  39 ) const
  40 {
  41     List<DynamicList<label, primitiveMesh::cellsPerPoint_> >
  42         pc(points().size());
  43
  44     // For each cell
  45     forAll(c, i)
  46     {
  47         // For each vertex
  48         const labelList& labels = c[i];
  49
  50         forAll(labels, j)
  51         {
  52             // Set working point label
  53             label curPoint = labels[j];
  54             DynamicList<label, primitiveMesh::cellsPerPoint_>& curPointCells =
  55                 pc[curPoint];
  56
  57             // Enter the cell label in the point's cell list
  58             curPointCells.append(i);
  59         }
  60     }
  61
  62     labelListList pointCellAddr(pc.size());
  63
  64     forAll (pc, pointI)
  65     {
  66         pointCellAddr[pointI].transfer(pc[pointI]);
  67     }
  68
  69     return pointCellAddr;
  70 }
  71
  72
  73 Foam::labelList Foam::polyMesh::facePatchFaceCells
  74 (
  75     const faceList& patchFaces,
  76     const labelListList& pointCells,
  77     const faceListList& cellsFaceShapes,
  78     const label patchID
  79 ) const
  80 {
  81     register bool found;
  82
  83     labelList FaceCells(patchFaces.size());
  84
  85     forAll(patchFaces, fI)
  86     {
  87         found = false;
  88
  89         const face& curFace = patchFaces[fI];
  90         const labelList& facePoints = patchFaces[fI];
  91
  92         forAll(facePoints, pointI)
  93         {
  94             const labelList& facePointCells = pointCells[facePoints[pointI]];
  95
  96             forAll(facePointCells, cellI)
  97             {
  98                 faceList cellFaces = cellsFaceShapes[facePointCells[cellI]];
  99
 100                 forAll(cellFaces, cellFace)
 101                 {
 102                     if (cellFaces[cellFace] == curFace)
 103                     {
 104                         // Found the cell corresponding to this face
 105                         FaceCells[fI] = facePointCells[cellI];
 106
 107                         found = true;
 108                     }
 109                     if (found) break;
 110                 }
 111                 if (found) break;
 112             }
 113             if (found) break;
 114         }
 115
 116         if (!found)
 117         {
 118             FatalErrorIn
 119             (
 120                 "polyMesh::facePatchFaceCells(const faceList& patchFaces,"
 121                 "const labelListList& pointCells,"
 122                 "const faceListList& cellsFaceShapes,"
 123                 "const label patchID)"
 124             )   << "face " << fI << " in patch " << patchID
 125                 << " does not have neighbour cell"
 126                 << " face: " << patchFaces[fI]
 127                 << abort(FatalError);
 128         }
 129     }
 130
 131     return FaceCells;
 132 }
 133
 134
 135 // Set faces, calculate cells and patchStarts.
 136 void Foam::polyMesh::setTopology
 137 (
 138     const cellShapeList& cellsAsShapes,
 139     const faceListList& boundaryFaces,
 140     const wordList& boundaryPatchNames,
 141     labelList& patchSizes,
 142     labelList& patchStarts,
 143     label& defaultPatchStart,
 144     label& nFaces,
 145     cellList& cells
 146 )
 147 {
 148     // Calculate the faces of all cells
 149     // Initialise maximum possible numer of mesh faces to 0
 150     label maxFaces = 0;
 151
 152     // Set up a list of face shapes for each cell
 153     faceListList cellsFaceShapes(cellsAsShapes.size());
 154     cells.setSize(cellsAsShapes.size());
 155
 156     forAll(cellsFaceShapes, cellI)
 157     {
 158         cellsFaceShapes[cellI] = cellsAsShapes[cellI].faces();
 159
 160         cells[cellI].setSize(cellsFaceShapes[cellI].size());
 161
 162         // Initialise cells to -1 to flag undefined faces
 163         static_cast<labelList&>(cells[cellI]) = -1;
 164
 165         // Count maximum possible numer of mesh faces
 166         maxFaces += cellsFaceShapes[cellI].size();
 167     }
 168
 169     // Set size of faces array to maximum possible number of mesh faces
 170     allFaces_.setSize(maxFaces);
 171
 172     // Initialise number of faces to 0
 173     nFaces = 0;
 174
 175     // set reference to point-cell addressing
 176     labelListList PointCells = cellShapePointCells(cellsAsShapes);
 177
 178     bool found = false;
 179
 180     forAll(cells, cellI)
 181     {
 182         // Note:
 183         // Insertion cannot be done in one go as the faces need to be
 184         // added into the list in the increasing order of neighbour
 185         // cells.  Therefore, all neighbours will be detected first
 186         // and then added in the correct order.
 187
 188         const faceList& curFaces = cellsFaceShapes[cellI];
 189
 190         // Record the neighbour cell
 191         labelList neiCells(curFaces.size(), -1);
 192
 193         // Record the face of neighbour cell
 194         labelList faceOfNeiCell(curFaces.size(), -1);
 195
 196         label nNeighbours = 0;
 197
 198         // For all faces ...
 199         forAll(curFaces, faceI)
 200         {
 201             // Skip faces that have already been matched
 202             if (cells[cellI][faceI] >= 0) continue;
 203
 204             found = false;
 205
 206             const face& curFace = curFaces[faceI];
 207
 208             // Get the list of labels
 209             const labelList& curPoints = curFace;
 210
 211             // For all points
 212             forAll(curPoints, pointI)
 213             {
 214                 // Get the list of cells sharing this point
 215                 const labelList& curNeighbours =
 216                     PointCells[curPoints[pointI]];
 217
 218                 // For all neighbours
 219                 forAll(curNeighbours, neiI)
 220                 {
 221                     label curNei = curNeighbours[neiI];
 222
 223                     // Reject neighbours with the lower label
 224                     if (curNei > cellI)
 225                     {
 226                         // Get the list of search faces
 227                         const faceList& searchFaces = cellsFaceShapes[curNei];
 228
 229                         forAll(searchFaces, neiFaceI)
 230                         {
 231                             if (searchFaces[neiFaceI] == curFace)
 232                             {
 233                                 // Match!!
 234                                 found = true;
 235
 236                                 // Record the neighbour cell and face
 237                                 neiCells[faceI] = curNei;
 238                                 faceOfNeiCell[faceI] = neiFaceI;
 239                                 nNeighbours++;
 240
 241                                 break;
 242                             }
 243                         }
 244                         if (found) break;
 245                     }
 246 //                     if (found) break; HJ, not needed.  Never true
 247                 }
 248                 if (found) break;
 249             } // End of current points
 250         }  // End of current faces
 251
 252         // Add the faces in the increasing order of neighbours
 253         for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
 254         {
 255             // Find the lowest neighbour which is still valid
 256             label nextNei = -1;
 257             label minNei = cells.size();
 258
 259             forAll (neiCells, ncI)
 260             {
 261                 if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
 262                 {
 263                     nextNei = ncI;
 264                     minNei = neiCells[ncI];
 265                 }
 266             }
 267
 268             if (nextNei > -1)
 269             {
 270                 // Add the face to the list of faces
 271                allFaces_[nFaces] = curFaces[nextNei];
 272
 273                 // Set cell-face and cell-neighbour-face to current face label
 274                 cells[cellI][nextNei] = nFaces;
 275                 cells[neiCells[nextNei]][faceOfNeiCell[nextNei]] = nFaces;
 276
 277                 // Stop the neighbour from being used again
 278                 neiCells[nextNei] = -1;
 279
 280                 // Increment number of faces counter
 281                 nFaces++;
 282             }
 283             else
 284             {
 285                 FatalErrorIn
 286                 (
 287                     "polyMesh::setTopology\n"
 288                     "(\n"
 289                     "    const cellShapeList& cellsAsShapes,\n"
 290                     "    const faceListList& boundaryFaces,\n"
 291                     "    const wordList& boundaryPatchNames,\n"
 292                     "    labelList& patchSizes,\n"
 293                     "    labelList& patchStarts,\n"
 294                     "    label& defaultPatchStart,\n"
 295                     "    label& nFaces,\n"
 296                     "    cellList& cells\n"
 297                     ")"
 298                 )   << "Error in internal face insertion"
 299                     << abort(FatalError);
 300             }
 301         }
 302     }
 303
 304     // Do boundary faces
 305
 306     patchSizes.setSize(boundaryFaces.size(), -1);
 307     patchStarts.setSize(boundaryFaces.size(), -1);
 308
 309     forAll (boundaryFaces, patchI)
 310     {
 311         const faceList& patchFaces = boundaryFaces[patchI];
 312
 313         labelList curPatchFaceCells =
 314             facePatchFaceCells
 315             (
 316                 patchFaces,
 317                 PointCells,
 318                 cellsFaceShapes,
 319                 patchI
 320             );
 321
 322         // Grab the start label
 323         label curPatchStart = nFaces;
 324
 325         forAll (patchFaces, faceI)
 326         {
 327             const face& curFace = patchFaces[faceI];
 328
 329             const label cellInside = curPatchFaceCells[faceI];
 330
 331             allFaces_[nFaces] = curFace;
 332
 333             // get faces of the cell inside
 334             const faceList& facesOfCellInside = cellsFaceShapes[cellInside];
 335
 336             bool found = false;
 337
 338             forAll (facesOfCellInside, cellFaceI)
 339             {
 340                 if (facesOfCellInside[cellFaceI] == curFace)
 341                 {
 342                     if (cells[cellInside][cellFaceI] >= 0)
 343                     {
 344                         FatalErrorIn
 345                         (
 346                             "polyMesh::setTopology\n"
 347                             "(\n"
 348                             "    const cellShapeList& cellsAsShapes,\n"
 349                             "    const faceListList& boundaryFaces,\n"
 350                             "    const wordList& boundaryPatchNames,\n"
 351                             "    labelList& patchSizes,\n"
 352                             "    labelList& patchStarts,\n"
 353                             "    label& defaultPatchStart,\n"
 354                             "    label& nFaces,\n"
 355                             "    cellList& cells\n"
 356                             ")"
 357                         )   << "Trying to specify a boundary face " << curFace
 358                             << " on the face on cell " << cellInside
 359                             << " which is either an internal face or already "
 360                             << "belongs to some other patch.  This is face "
 361                             << faceI << " of patch "
 362                             << patchI << " named "
 363                             << boundaryPatchNames[patchI] << "."
 364                             << abort(FatalError);
 365                     }
 366
 367                     found = true;
 368
 369                     cells[cellInside][cellFaceI] = nFaces;
 370
 371                     break;
 372                 }
 373             }
 374
 375             if (!found)
 376             {
 377                 FatalErrorIn("polyMesh::polyMesh(... construct from shapes...)")
 378                     << "face " << faceI << " of patch " << patchI
 379                     << " does not seem to belong to cell " << cellInside
 380                     << " which, according to the addressing, "
 381                     << "should be next to it."
 382                     << abort(FatalError);
 383             }
 384
 385             // increment the counter of faces
 386             nFaces++;
 387         }
 388
 389         patchSizes[patchI] = nFaces - curPatchStart;
 390         patchStarts[patchI] = curPatchStart;
 391     }
 392
 393     // Grab "non-existing" faces and put them into a default patch
 394
 395     defaultPatchStart = nFaces;
 396
 397     forAll(cells, cellI)
 398     {
 399         labelList& curCellFaces = cells[cellI];
 400
 401         forAll(curCellFaces, faceI)
 402         {
 403             if (curCellFaces[faceI] == -1) // "non-existent" face
 404             {
 405                 curCellFaces[faceI] = nFaces;
 406                 allFaces_[nFaces] = cellsFaceShapes[cellI][faceI];
 407
 408                 nFaces++;
 409             }
 410         }
 411     }
 412
 413     // Reset the size of the face list
 414     allFaces_.setSize(nFaces);
 415     faces_.reset(allFaces_, nFaces);
 416
 417     return ;
 418 }
 419
 420
 421 Foam::polyMesh::polyMesh
 422 (
 423     const IOobject& io,
 424     const Xfer<pointField>& points,
 425     const cellShapeList& cellsAsShapes,
 426     const faceListList& boundaryFaces,
 427     const wordList& boundaryPatchNames,
 428     const wordList& boundaryPatchTypes,
 429     const word& defaultBoundaryPatchName,
 430     const word& defaultBoundaryPatchType,
 431     const wordList& boundaryPatchPhysicalTypes,
 432     const bool syncPar
 433 )
 434 :
 435     objectRegistry(io),
 436     primitiveMesh(),
 437     allPoints_
 438     (
 439         IOobject
 440         (
 441             "points",
 442             instance(),
 443             meshSubDir,
 444             *this,
 445             IOobject::NO_READ,
 446             IOobject::AUTO_WRITE
 447         ),
 448         points
 449     ),
 450     // To be re-sliced later.  HJ, 19/oct/2008
 451     points_(allPoints_, allPoints_.size()),
 452     allFaces_
 453     (
 454         IOobject
 455         (
 456             "faces",
 457             instance(),
 458             meshSubDir,
 459             *this,
 460             IOobject::NO_READ,
 461             IOobject::AUTO_WRITE
 462         ),
 463         0
 464     ),
 465     faces_(allFaces_, allFaces_.size()),
 466     owner_
 467     (
 468         IOobject
 469         (
 470             "owner",
 471             instance(),
 472             meshSubDir,
 473             *this,
 474             IOobject::NO_READ,
 475             IOobject::AUTO_WRITE
 476         ),
 477         0
 478     ),
 479     neighbour_
 480     (
 481         IOobject
 482         (
 483             "neighbour",
 484             instance(),
 485             meshSubDir,
 486             *this,
 487             IOobject::NO_READ,
 488             IOobject::AUTO_WRITE
 489         ),
 490         0
 491     ),
 492     clearedPrimitives_(false),
 493     boundary_
 494     (
 495         IOobject
 496         (
 497             "boundary",
 498             instance(),
 499             meshSubDir,
 500             *this,
 501             IOobject::NO_READ,
 502             IOobject::AUTO_WRITE
 503         ),
 504         *this,
 505         boundaryFaces.size() + 1    // add room for a default patch
 506     ),
 507     bounds_(allPoints_, syncPar),
 508     geometricD_(Vector<label>::zero),
 509     solutionD_(Vector<label>::zero),
 510     pointZones_
 511     (
 512         IOobject
 513         (
 514             "pointZones",
 515             instance(),
 516             meshSubDir,
 517             *this,
 518             IOobject::NO_READ,
 519             IOobject::NO_WRITE
 520         ),
 521         *this,
 522         0
 523     ),
 524     faceZones_
 525     (
 526         IOobject
 527         (
 528             "faceZones",
 529             instance(),
 530             meshSubDir,
 531             *this,
 532             IOobject::NO_READ,
 533             IOobject::NO_WRITE
 534         ),
 535         *this,
 536         0
 537     ),
 538     cellZones_
 539     (
 540         IOobject
 541         (
 542             "cellZones",
 543             instance(),
 544             meshSubDir,
 545             *this,
 546             IOobject::NO_READ,
 547             IOobject::NO_WRITE
 548         ),
 549         *this,
 550         0
 551     ),
 552     globalMeshDataPtr_(NULL),
 553     moving_(false),
 554     changing_(false),
 555     curMotionTimeIndex_(time().timeIndex()),
 556     oldAllPointsPtr_(NULL),
 557     oldPointsPtr_(NULL)
 558 {
 559     if (debug)
 560     {
 561         Info<<"Constructing polyMesh from cell and boundary shapes." << endl;
 562     }
 563
 564     // Remove all of the old mesh files if they exist
 565     removeFiles(instance());
 566
 567     // Calculate faces and cells
 568     labelList patchSizes;
 569     labelList patchStarts;
 570     label defaultPatchStart;
 571     label nFaces;
 572     cellList cells;
 573     setTopology
 574     (
 575         cellsAsShapes,
 576         boundaryFaces,
 577         boundaryPatchNames,
 578         patchSizes,
 579         patchStarts,
 580         defaultPatchStart,
 581         nFaces,
 582         cells
 583     );
 584
 585     // Warning: Patches can only be added once the face list is
 586     // completed, as they hold a subList of the face list
 587     forAll (boundaryFaces, patchI)
 588     {
 589         // add the patch to the list
 590         boundary_.set
 591         (
 592             patchI,
 593             polyPatch::New
 594             (
 595                 boundaryPatchTypes[patchI],
 596                 boundaryPatchNames[patchI],
 597                 patchSizes[patchI],
 598                 patchStarts[patchI],
 599                 patchI,
 600                 boundary_
 601             )
 602         );
 603
 604         if
 605         (
 606             boundaryPatchPhysicalTypes.size()
 607          && boundaryPatchPhysicalTypes[patchI].size()
 608         )
 609         {
 610             boundary_[patchI].physicalType() =
 611                 boundaryPatchPhysicalTypes[patchI];
 612         }
 613     }
 614
 615     label nAllPatches = boundaryFaces.size();
 616
 617
 618     label nDefaultFaces = nFaces - defaultPatchStart;
 619     if (syncPar)
 620     {
 621         reduce(nDefaultFaces, sumOp<label>());
 622     }
 623
 624     if (nDefaultFaces > 0)
 625     {
 626         WarningIn("polyMesh::polyMesh(... construct from shapes...)")
 627             << "Found " << nDefaultFaces
 628             << " undefined faces in mesh; adding to default patch." << endl;
 629
 630         // Check if there already exists a defaultFaces patch as last patch
 631         // and reuse it.
 632         label patchI = findIndex(boundaryPatchNames, defaultBoundaryPatchName);
 633
 634         if (patchI != -1)
 635         {
 636             if (patchI != boundaryFaces.size()-1 || boundary_[patchI].size())
 637             {
 638                 FatalErrorIn("polyMesh::polyMesh(... construct from shapes...)")
 639                     << "Default patch " << boundary_[patchI].name()
 640                     << " already has faces in it or is not"
 641                     << " last in list of patches." << exit(FatalError);
 642             }
 643
 644             WarningIn("polyMesh::polyMesh(... construct from shapes...)")
 645                 << "Reusing existing patch " << patchI
 646                 << " for undefined faces." << endl;
 647
 648             boundary_.set
 649             (
 650                 patchI,
 651                 polyPatch::New
 652                 (
 653                     boundary_[patchI].type(),
 654                     boundary_[patchI].name(),
 655                     nFaces - defaultPatchStart,
 656                     defaultPatchStart,
 657                     patchI,
 658                     boundary_
 659                 )
 660             );
 661         }
 662         else
 663         {
 664             boundary_.set
 665             (
 666                 nAllPatches,
 667                 polyPatch::New
 668                 (
 669                     defaultBoundaryPatchType,
 670                     defaultBoundaryPatchName,
 671                     nFaces - defaultPatchStart,
 672                     defaultPatchStart,
 673                     boundary_.size() - 1,
 674                     boundary_
 675                 )
 676             );
 677
 678             nAllPatches++;
 679         }
 680     }
 681
 682     // Reset the size of the boundary
 683     boundary_.setSize(nAllPatches);
 684
 685     // Set the primitive mesh
 686     initMesh(cells);
 687
 688     if (syncPar)
 689     {
 690         // Calculate topology for the patches (processor-processor comms etc.)
 691         boundary_.updateMesh();
 692
 693         // Calculate the geometry for the patches (transformation tensors etc.)
 694         boundary_.calcGeometry();
 695     }
 696
 697     if (debug)
 698     {
 699         if (checkMesh())
 700         {
 701             Info << "Mesh OK" << endl;
 702         }
 703     }
 704 }
 705
 706
 707 Foam::polyMesh::polyMesh
 708 (
 709     const IOobject& io,
 710     const Xfer<pointField>& points,
 711     const cellShapeList& cellsAsShapes,
 712     const faceListList& boundaryFaces,
 713     const wordList& boundaryPatchNames,
 714     const PtrList<dictionary>& boundaryDicts,
 715     const word& defaultBoundaryPatchName,
 716     const word& defaultBoundaryPatchType,
 717     const bool syncPar
 718 )
 719 :
 720     objectRegistry(io),
 721     primitiveMesh(),
 722     allPoints_
 723     (
 724         IOobject
 725         (
 726             "points",
 727             instance(),
 728             meshSubDir,
 729             *this,
 730             IOobject::NO_READ,
 731             IOobject::AUTO_WRITE
 732         ),
 733         points
 734     ),
 735     // To be re-sliced later.  HJ, 19/oct/2008
 736     points_(allPoints_, allPoints_.size()),
 737     allFaces_
 738     (
 739         IOobject
 740         (
 741             "faces",
 742             instance(),
 743             meshSubDir,
 744             *this,
 745             IOobject::NO_READ,
 746             IOobject::AUTO_WRITE
 747         ),
 748         0
 749     ),
 750     faces_(allFaces_, allFaces_.size()),
 751     owner_
 752     (
 753         IOobject
 754         (
 755             "owner",
 756             instance(),
 757             meshSubDir,
 758             *this,
 759             IOobject::NO_READ,
 760             IOobject::AUTO_WRITE
 761         ),
 762         0
 763     ),
 764     neighbour_
 765     (
 766         IOobject
 767         (
 768             "neighbour",
 769             instance(),
 770             meshSubDir,
 771             *this,
 772             IOobject::NO_READ,
 773             IOobject::AUTO_WRITE
 774         ),
 775         0
 776     ),
 777     clearedPrimitives_(false),
 778     boundary_
 779     (
 780         IOobject
 781         (
 782             "boundary",
 783             instance(),
 784             meshSubDir,
 785             *this,
 786             IOobject::NO_READ,
 787             IOobject::AUTO_WRITE
 788         ),
 789         *this,
 790         boundaryFaces.size() + 1    // add room for a default patch
 791     ),
 792     bounds_(allPoints_, syncPar),
 793     geometricD_(Vector<label>::zero),
 794     solutionD_(Vector<label>::zero),
 795     pointZones_
 796     (
 797         IOobject
 798         (
 799             "pointZones",
 800             instance(),
 801             meshSubDir,
 802             *this,
 803             IOobject::NO_READ,
 804             IOobject::NO_WRITE
 805         ),
 806         *this,
 807         0
 808     ),
 809     faceZones_
 810     (
 811         IOobject
 812         (
 813             "faceZones",
 814             instance(),
 815             meshSubDir,
 816             *this,
 817             IOobject::NO_READ,
 818             IOobject::NO_WRITE
 819         ),
 820         *this,
 821         0
 822     ),
 823     cellZones_
 824     (
 825         IOobject
 826         (
 827             "cellZones",
 828             instance(),
 829             meshSubDir,
 830             *this,
 831             IOobject::NO_READ,
 832             IOobject::NO_WRITE
 833         ),
 834         *this,
 835         0
 836     ),
 837     globalMeshDataPtr_(NULL),
 838     moving_(false),
 839     changing_(false),
 840     curMotionTimeIndex_(time().timeIndex()),
 841     oldAllPointsPtr_(NULL),
 842     oldPointsPtr_(NULL)
 843 {
 844     if (debug)
 845     {
 846         Info<<"Constructing polyMesh from cell and boundary shapes." << endl;
 847     }
 848
 849     // Remove all of the old mesh files if they exist
 850     removeFiles(instance());
 851
 852     // Calculate faces and cells
 853     labelList patchSizes;
 854     labelList patchStarts;
 855     label defaultPatchStart;
 856     label nFaces;
 857     cellList cells;
 858     setTopology
 859     (
 860         cellsAsShapes,
 861         boundaryFaces,
 862         boundaryPatchNames,
 863         patchSizes,
 864         patchStarts,
 865         defaultPatchStart,
 866         nFaces,
 867         cells
 868     );
 869
 870     // Warning: Patches can only be added once the face list is
 871     // completed, as they hold a subList of the face list
 872     forAll (boundaryDicts, patchI)
 873     {
 874         dictionary patchDict(boundaryDicts[patchI]);
 875
 876         patchDict.set("nFaces", patchSizes[patchI]);
 877         patchDict.set("startFace", patchStarts[patchI]);
 878
 879         // add the patch to the list
 880         boundary_.set
 881         (
 882             patchI,
 883             polyPatch::New
 884             (
 885                 boundaryPatchNames[patchI],
 886                 patchDict,
 887                 patchI,
 888                 boundary_
 889             )
 890         );
 891     }
 892
 893     label nAllPatches = boundaryFaces.size();
 894
 895     label nDefaultFaces = nFaces - defaultPatchStart;
 896     if (syncPar)
 897     {
 898         reduce(nDefaultFaces, sumOp<label>());
 899     }
 900
 901     if (nDefaultFaces > 0)
 902     {
 903         WarningIn("polyMesh::polyMesh(... construct from shapes...)")
 904             << "Found " << nDefaultFaces
 905             << " undefined faces in mesh; adding to default patch." << endl;
 906
 907         // Check if there already exists a defaultFaces patch as last patch
 908         // and reuse it.
 909         label patchI = findIndex(boundaryPatchNames, defaultBoundaryPatchName);
 910
 911         if (patchI != -1)
 912         {
 913             if (patchI != boundaryFaces.size()-1 || boundary_[patchI].size())
 914             {
 915                 FatalErrorIn("polyMesh::polyMesh(... construct from shapes...)")
 916                     << "Default patch " << boundary_[patchI].name()
 917                     << " already has faces in it or is not"
 918                     << " last in list of patches." << exit(FatalError);
 919             }
 920
 921             WarningIn("polyMesh::polyMesh(... construct from shapes...)")
 922                 << "Reusing existing patch " << patchI
 923                 << " for undefined faces." << endl;
 924
 925             boundary_.set
 926             (
 927                 patchI,
 928                 polyPatch::New
 929                 (
 930                     boundary_[patchI].type(),
 931                     boundary_[patchI].name(),
 932                     nFaces - defaultPatchStart,
 933                     defaultPatchStart,
 934                     patchI,
 935                     boundary_
 936                 )
 937             );
 938         }
 939         else
 940         {
 941             boundary_.set
 942             (
 943                 nAllPatches,
 944                 polyPatch::New
 945                 (
 946                     defaultBoundaryPatchType,
 947                     defaultBoundaryPatchName,
 948                     nFaces - defaultPatchStart,
 949                     defaultPatchStart,
 950                     boundary_.size() - 1,
 951                     boundary_
 952                 )
 953             );
 954
 955             nAllPatches++;
 956         }
 957     }
 958
 959     // Reset the size of the boundary
 960     boundary_.setSize(nAllPatches);
 961
 962     // Set the primitive mesh
 963     initMesh(cells);
 964
 965     if (syncPar)
 966     {
 967         // Calculate topology for the patches (processor-processor comms etc.)
 968         boundary_.updateMesh();
 969
 970         // Calculate the geometry for the patches (transformation tensors etc.)
 971         boundary_.calcGeometry();
 972     }
 973
 974     if (debug)
 975     {
 976         if (checkMesh())
 977         {
 978             Info << "Mesh OK" << endl;
 979         }
 980     }
 981 }
 982
 983
 984 // ************************************************************************* //