src/OpenFOAM/meshes/polyMesh/polyMeshFromShapeMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     Create polyMesh from cell and patch shapes
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "polyMesh.H"
  31 #include <OpenFOAM/Time.H>
  32 #include <OpenFOAM/primitiveMesh.H>
  33 #include <OpenFOAM/DynamicList.H>
  34
  35 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  36
  37 Foam::labelListList Foam::polyMesh::cellShapePointCells
  38 (
  39     const cellShapeList& c
  40 ) const
  41 {
  42     List<DynamicList<label, primitiveMesh::cellsPerPoint_> >
  43         pc(points().size());
  44
  45     // For each cell
  46     forAll(c, i)
  47     {
  48         // For each vertex
  49         const labelList& labels = c[i];
  50
  51         forAll(labels, j)
  52         {
  53             // Set working point label
  54             label curPoint = labels[j];
  55             DynamicList<label, primitiveMesh::cellsPerPoint_>& curPointCells =
  56                 pc[curPoint];
  57
  58             // Enter the cell label in the point's cell list
  59             curPointCells.append(i);
  60         }
  61     }
  62
  63     labelListList pointCellAddr(pc.size());
  64
  65     forAll (pc, pointI)
  66     {
  67         pointCellAddr[pointI].transfer(pc[pointI]);
  68     }
  69
  70     return pointCellAddr;
  71 }
  72
  73
  74 Foam::labelList Foam::polyMesh::facePatchFaceCells
  75 (
  76     const faceList& patchFaces,
  77     const labelListList& pointCells,
  78     const faceListList& cellsFaceShapes,
  79     const label patchID
  80 ) const
  81 {
  82     register bool found;
  83
  84     labelList FaceCells(patchFaces.size());
  85
  86     forAll(patchFaces, fI)
  87     {
  88         found = false;
  89
  90         const face& curFace = patchFaces[fI];
  91         const labelList& facePoints = patchFaces[fI];
  92
  93         forAll(facePoints, pointI)
  94         {
  95             const labelList& facePointCells = pointCells[facePoints[pointI]];
  96
  97             forAll(facePointCells, cellI)
  98             {
  99                 faceList cellFaces = cellsFaceShapes[facePointCells[cellI]];
 100
 101                 forAll(cellFaces, cellFace)
 102                 {
 103                     if (cellFaces[cellFace] == curFace)
 104                     {
 105                         // Found the cell corresponding to this face
 106                         FaceCells[fI] = facePointCells[cellI];
 107
 108                         found = true;
 109                     }
 110                     if (found) break;
 111                 }
 112                 if (found) break;
 113             }
 114             if (found) break;
 115         }
 116
 117         if (!found)
 118         {
 119             FatalErrorIn
 120             (
 121                 "polyMesh::facePatchFaceCells(const faceList& patchFaces,"
 122                 "const labelListList& pointCells,"
 123                 "const faceListList& cellsFaceShapes,"
 124                 "const label patchID)"
 125             )   << "face " << fI << " in patch " << patchID
 126                 << " does not have neighbour cell"
 127                 << " face: " << patchFaces[fI]
 128                 << abort(FatalError);
 129         }
 130     }
 131
 132     return FaceCells;
 133 }
 134
 135
 136 Foam::polyMesh::polyMesh
 137 (
 138     const IOobject& io,
 139     const Xfer<pointField>& points,
 140     const cellShapeList& cellsAsShapes,
 141     const faceListList& boundaryFaces,
 142     const wordList& boundaryPatchNames,
 143     const wordList& boundaryPatchTypes,
 144     const word& defaultBoundaryPatchName,
 145     const word& defaultBoundaryPatchType,
 146     const wordList& boundaryPatchPhysicalTypes,
 147     const bool syncPar
 148 )
 149 :
 150     objectRegistry(io),
 151     primitiveMesh(),
 152     points_
 153     (
 154         IOobject
 155         (
 156             "points",
 157             instance(),
 158             meshSubDir,
 159             *this,
 160             IOobject::NO_READ,
 161             IOobject::AUTO_WRITE
 162         ),
 163         points
 164     ),
 165     faces_
 166     (
 167         IOobject
 168         (
 169             "faces",
 170             instance(),
 171             meshSubDir,
 172             *this,
 173             IOobject::NO_READ,
 174             IOobject::AUTO_WRITE
 175         ),
 176         0
 177     ),
 178     owner_
 179     (
 180         IOobject
 181         (
 182             "owner",
 183             instance(),
 184             meshSubDir,
 185             *this,
 186             IOobject::NO_READ,
 187             IOobject::AUTO_WRITE
 188         ),
 189         0
 190     ),
 191     neighbour_
 192     (
 193         IOobject
 194         (
 195             "neighbour",
 196             instance(),
 197             meshSubDir,
 198             *this,
 199             IOobject::NO_READ,
 200             IOobject::AUTO_WRITE
 201         ),
 202         0
 203     ),
 204     clearedPrimitives_(false),
 205     boundary_
 206     (
 207         IOobject
 208         (
 209             "boundary",
 210             instance(),
 211             meshSubDir,
 212             *this,
 213             IOobject::NO_READ,
 214             IOobject::AUTO_WRITE
 215         ),
 216         *this,
 217         boundaryFaces.size() + 1    // add room for a default patch
 218     ),
 219     bounds_(points_, syncPar),
 220     geometricD_(Vector<label>::zero),
 221     solutionD_(Vector<label>::zero),
 222     pointZones_
 223     (
 224         IOobject
 225         (
 226             "pointZones",
 227             instance(),
 228             meshSubDir,
 229             *this,
 230             IOobject::NO_READ,
 231             IOobject::NO_WRITE
 232         ),
 233         *this,
 234         0
 235     ),
 236     faceZones_
 237     (
 238         IOobject
 239         (
 240             "faceZones",
 241             instance(),
 242             meshSubDir,
 243             *this,
 244             IOobject::NO_READ,
 245             IOobject::NO_WRITE
 246         ),
 247         *this,
 248         0
 249     ),
 250     cellZones_
 251     (
 252         IOobject
 253         (
 254             "cellZones",
 255             instance(),
 256             meshSubDir,
 257             *this,
 258             IOobject::NO_READ,
 259             IOobject::NO_WRITE
 260         ),
 261         *this,
 262         0
 263     ),
 264     globalMeshDataPtr_(NULL),
 265     moving_(false),
 266     curMotionTimeIndex_(time().timeIndex()),
 267     oldPointsPtr_(NULL)
 268 {
 269     if (debug)
 270     {
 271         Info<<"Constructing polyMesh from cell and boundary shapes." << endl;
 272     }
 273
 274     // Remove all of the old mesh files if they exist
 275     removeFiles(instance());
 276
 277     // Calculate the faces of all cells
 278     // Initialise maximum possible numer of mesh faces to 0
 279     label maxFaces = 0;
 280
 281     // Set up a list of face shapes for each cell
 282     faceListList cellsFaceShapes(cellsAsShapes.size());
 283     cellList cells(cellsAsShapes.size());
 284
 285     forAll(cellsFaceShapes, cellI)
 286     {
 287         cellsFaceShapes[cellI] = cellsAsShapes[cellI].faces();
 288
 289         cells[cellI].setSize(cellsFaceShapes[cellI].size());
 290
 291         // Initialise cells to -1 to flag undefined faces
 292         static_cast<labelList&>(cells[cellI]) = -1;
 293
 294         // Count maximum possible numer of mesh faces
 295         maxFaces += cellsFaceShapes[cellI].size();
 296     }
 297
 298     // Set size of faces array to maximum possible number of mesh faces
 299     faces_.setSize(maxFaces);
 300
 301     // Initialise number of faces to 0
 302     label nFaces = 0;
 303
 304     // set reference to point-cell addressing
 305     labelListList PointCells = cellShapePointCells(cellsAsShapes);
 306
 307     bool found = false;
 308
 309     forAll(cells, cellI)
 310     {
 311         // Note:
 312         // Insertion cannot be done in one go as the faces need to be
 313         // added into the list in the increasing order of neighbour
 314         // cells.  Therefore, all neighbours will be detected first
 315         // and then added in the correct order.
 316
 317         const faceList& curFaces = cellsFaceShapes[cellI];
 318
 319         // Record the neighbour cell
 320         labelList neiCells(curFaces.size(), -1);
 321
 322         // Record the face of neighbour cell
 323         labelList faceOfNeiCell(curFaces.size(), -1);
 324
 325         label nNeighbours = 0;
 326
 327         // For all faces ...
 328         forAll(curFaces, faceI)
 329         {
 330             // Skip faces that have already been matched
 331             if (cells[cellI][faceI] >= 0) continue;
 332
 333             found = false;
 334
 335             const face& curFace = curFaces[faceI];
 336
 337             // Get the list of labels
 338             const labelList& curPoints = curFace;
 339
 340             // For all points
 341             forAll(curPoints, pointI)
 342             {
 343                 // dGget the list of cells sharing this point
 344                 const labelList& curNeighbours =
 345                     PointCells[curPoints[pointI]];
 346
 347                 // For all neighbours
 348                 forAll(curNeighbours, neiI)
 349                 {
 350                     label curNei = curNeighbours[neiI];
 351
 352                     // Reject neighbours with the lower label
 353                     if (curNei > cellI)
 354                     {
 355                         // Get the list of search faces
 356                         const faceList& searchFaces = cellsFaceShapes[curNei];
 357
 358                         forAll(searchFaces, neiFaceI)
 359                         {
 360                             if (searchFaces[neiFaceI] == curFace)
 361                             {
 362                                 // Match!!
 363                                 found = true;
 364
 365                                 // Record the neighbour cell and face
 366                                 neiCells[faceI] = curNei;
 367                                 faceOfNeiCell[faceI] = neiFaceI;
 368                                 nNeighbours++;
 369
 370                                 break;
 371                             }
 372                         }
 373                         if (found) break;
 374                     }
 375                     if (found) break;
 376                 }
 377                 if (found) break;
 378             } // End of current points
 379         }  // End of current faces
 380
 381         // Add the faces in the increasing order of neighbours
 382         for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
 383         {
 384             // Find the lowest neighbour which is still valid
 385             label nextNei = -1;
 386             label minNei = cells.size();
 387
 388             forAll (neiCells, ncI)
 389             {
 390                 if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
 391                 {
 392                     nextNei = ncI;
 393                     minNei = neiCells[ncI];
 394                 }
 395             }
 396
 397             if (nextNei > -1)
 398             {
 399                 // Add the face to the list of faces
 400                 faces_[nFaces] = curFaces[nextNei];
 401
 402                 // Set cell-face and cell-neighbour-face to current face label
 403                 cells[cellI][nextNei] = nFaces;
 404                 cells[neiCells[nextNei]][faceOfNeiCell[nextNei]] = nFaces;
 405
 406                 // Stop the neighbour from being used again
 407                 neiCells[nextNei] = -1;
 408
 409                 // Increment number of faces counter
 410                 nFaces++;
 411             }
 412             else
 413             {
 414                 FatalErrorIn
 415                 (
 416                     "polyMesh::polyMesh\n"
 417                     "(\n"
 418                     "    const IOobject&,\n"
 419                     "    const Xfer<pointField>&,\n"
 420                     "    const cellShapeList& cellsAsShapes,\n"
 421                     "    const faceListList& boundaryFaces,\n"
 422                     "    const wordList& boundaryPatchTypes,\n"
 423                     "    const wordList& boundaryPatchNames,\n"
 424                     "    const word& defaultBoundaryPatchType\n"
 425                     ")"
 426                 )   << "Error in internal face insertion"
 427                     << abort(FatalError);
 428             }
 429         }
 430     }
 431
 432     // Do boundary faces
 433
 434     labelList patchSizes(boundaryFaces.size(), -1);
 435     labelList patchStarts(boundaryFaces.size(), -1);
 436
 437     forAll (boundaryFaces, patchI)
 438     {
 439         const faceList& patchFaces = boundaryFaces[patchI];
 440
 441         labelList curPatchFaceCells =
 442             facePatchFaceCells
 443             (
 444                 patchFaces,
 445                 PointCells,
 446                 cellsFaceShapes,
 447                 patchI
 448             );
 449
 450         // Grab the start label
 451         label curPatchStart = nFaces;
 452
 453         forAll (patchFaces, faceI)
 454         {
 455             const face& curFace = patchFaces[faceI];
 456
 457             const label cellInside = curPatchFaceCells[faceI];
 458
 459             faces_[nFaces] = curFace;
 460
 461             // get faces of the cell inside
 462             const faceList& facesOfCellInside = cellsFaceShapes[cellInside];
 463
 464             bool found = false;
 465
 466             forAll (facesOfCellInside, cellFaceI)
 467             {
 468                 if (facesOfCellInside[cellFaceI] == curFace)
 469                 {
 470                     if (cells[cellInside][cellFaceI] >= 0)
 471                     {
 472                         FatalErrorIn
 473                         (
 474                             "polyMesh::polyMesh\n"
 475                             "(\n"
 476                             "    const IOobject&,\n"
 477                             "    const Xfer<pointField>&,\n"
 478                             "    const cellShapeList& cellsAsShapes,\n"
 479                             "    const faceListList& boundaryFaces,\n"
 480                             "    const wordList& boundaryPatchTypes,\n"
 481                             "    const wordList& boundaryPatchNames,\n"
 482                             "    const word& defaultBoundaryPatchType\n"
 483                             ")"
 484                         )   << "Trying to specify a boundary face " << curFace
 485                             << " on the face on cell " << cellInside
 486                             << " which is either an internal face or already "
 487                             << "belongs to some other patch.  This is face "
 488                             << faceI << " of patch "
 489                             << patchI << " named "
 490                             << boundaryPatchNames[patchI] << "."
 491                             << abort(FatalError);
 492                     }
 493
 494                     found = true;
 495
 496                     cells[cellInside][cellFaceI] = nFaces;
 497
 498                     break;
 499                 }
 500             }
 501
 502             if (!found)
 503             {
 504                 FatalErrorIn("polyMesh::polyMesh(... construct from shapes...)")
 505                     << "face " << faceI << " of patch " << patchI
 506                     << " does not seem to belong to cell " << cellInside
 507                     << " which, according to the addressing, "
 508                     << "should be next to it."
 509                     << abort(FatalError);
 510             }
 511
 512             // increment the counter of faces
 513             nFaces++;
 514         }
 515
 516         patchSizes[patchI] = nFaces - curPatchStart;
 517         patchStarts[patchI] = curPatchStart;
 518     }
 519
 520     // Grab "non-existing" faces and put them into a default patch
 521
 522     label defaultPatchStart = nFaces;
 523
 524     forAll(cells, cellI)
 525     {
 526         labelList& curCellFaces = cells[cellI];
 527
 528         forAll(curCellFaces, faceI)
 529         {
 530             if (curCellFaces[faceI] == -1) // "non-existent" face
 531             {
 532                 curCellFaces[faceI] = nFaces;
 533                 faces_[nFaces] = cellsFaceShapes[cellI][faceI];
 534
 535                 nFaces++;
 536             }
 537         }
 538     }
 539
 540     // Reset the size of the face list
 541     faces_.setSize(nFaces);
 542
 543     // Warning: Patches can only be added once the face list is
 544     // completed, as they hold a subList of the face list
 545     forAll (boundaryFaces, patchI)
 546     {
 547         // add the patch to the list
 548         boundary_.set
 549         (
 550             patchI,
 551             polyPatch::New
 552             (
 553                 boundaryPatchTypes[patchI],
 554                 boundaryPatchNames[patchI],
 555                 patchSizes[patchI],
 556                 patchStarts[patchI],
 557                 patchI,
 558                 boundary_
 559             )
 560         );
 561
 562         if
 563         (
 564             boundaryPatchPhysicalTypes.size()
 565          && boundaryPatchPhysicalTypes[patchI].size()
 566         )
 567         {
 568             boundary_[patchI].physicalType() =
 569                 boundaryPatchPhysicalTypes[patchI];
 570         }
 571     }
 572
 573     label nAllPatches = boundaryFaces.size();
 574
 575     if (nFaces > defaultPatchStart)
 576     {
 577         WarningIn("polyMesh::polyMesh(... construct from shapes...)")
 578             << "Found " << nFaces - defaultPatchStart
 579             << " undefined faces in mesh; adding to default patch." << endl;
 580
 581         boundary_.set
 582         (
 583             nAllPatches,
 584             polyPatch::New
 585             (
 586                 defaultBoundaryPatchType,
 587                 defaultBoundaryPatchName,
 588                 nFaces - defaultPatchStart,
 589                 defaultPatchStart,
 590                 boundary_.size() - 1,
 591                 boundary_
 592             )
 593         );
 594
 595         nAllPatches++;
 596     }
 597
 598     // Reset the size of the boundary
 599     boundary_.setSize(nAllPatches);
 600
 601     // Set the primitive mesh
 602     initMesh(cells);
 603
 604     if (syncPar)
 605     {
 606         // Calculate topology for the patches (processor-processor comms etc.)
 607         boundary_.updateMesh();
 608
 609         // Calculate the geometry for the patches (transformation tensors etc.)
 610         boundary_.calcGeometry();
 611     }
 612
 613     if (debug)
 614     {
 615         if (checkMesh())
 616         {
 617             Info << "Mesh OK" << endl;
 618         }
 619     }
 620 }
 621
 622
 623 // ************************ vim: set sw=4 sts=4 et: ************************ //