FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

polyMeshFromShapeMesh.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Description
00025     Create polyMesh from cell and patch shapes
00026 
00027 \*---------------------------------------------------------------------------*/
00028 
00029 #include "polyMesh.H"
00030 #include <OpenFOAM/Time.H>
00031 #include <OpenFOAM/primitiveMesh.H>
00032 #include <OpenFOAM/DynamicList.H>
00033 
00034 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00035 
00036 Foam::labelListList Foam::polyMesh::cellShapePointCells
00037 (
00038     const cellShapeList& c
00039 ) const
00040 {
00041     List<DynamicList<label, primitiveMesh::cellsPerPoint_> >
00042         pc(points().size());
00043 
00044     // For each cell
00045     forAll(c, i)
00046     {
00047         // For each vertex
00048         const labelList& labels = c[i];
00049 
00050         forAll(labels, j)
00051         {
00052             // Set working point label
00053             label curPoint = labels[j];
00054             DynamicList<label, primitiveMesh::cellsPerPoint_>& curPointCells =
00055                 pc[curPoint];
00056 
00057             // Enter the cell label in the point's cell list
00058             curPointCells.append(i);
00059         }
00060     }
00061 
00062     labelListList pointCellAddr(pc.size());
00063 
00064     forAll (pc, pointI)
00065     {
00066         pointCellAddr[pointI].transfer(pc[pointI]);
00067     }
00068 
00069     return pointCellAddr;
00070 }
00071 
00072 
00073 Foam::labelList Foam::polyMesh::facePatchFaceCells
00074 (
00075     const faceList& patchFaces,
00076     const labelListList& pointCells,
00077     const faceListList& cellsFaceShapes,
00078     const label patchID
00079 ) const
00080 {
00081     register bool found;
00082 
00083     labelList FaceCells(patchFaces.size());
00084 
00085     forAll(patchFaces, fI)
00086     {
00087         found = false;
00088 
00089         const face& curFace = patchFaces[fI];
00090         const labelList& facePoints = patchFaces[fI];
00091 
00092         forAll(facePoints, pointI)
00093         {
00094             const labelList& facePointCells = pointCells[facePoints[pointI]];
00095 
00096             forAll(facePointCells, cellI)
00097             {
00098                 faceList cellFaces = cellsFaceShapes[facePointCells[cellI]];
00099 
00100                 forAll(cellFaces, cellFace)
00101                 {
00102                     if (cellFaces[cellFace] == curFace)
00103                     {
00104                         // Found the cell corresponding to this face
00105                         FaceCells[fI] = facePointCells[cellI];
00106 
00107                         found = true;
00108                     }
00109                     if (found) break;
00110                 }
00111                 if (found) break;
00112             }
00113             if (found) break;
00114         }
00115 
00116         if (!found)
00117         {
00118             FatalErrorIn
00119             (
00120                 "polyMesh::facePatchFaceCells(const faceList& patchFaces,"
00121                 "const labelListList& pointCells,"
00122                 "const faceListList& cellsFaceShapes,"
00123                 "const label patchID)"
00124             )   << "face " << fI << " in patch " << patchID
00125                 << " does not have neighbour cell"
00126                 << " face: " << patchFaces[fI]
00127                 << abort(FatalError);
00128         }
00129     }
00130 
00131     return FaceCells;
00132 }
00133 
00134 
00135 Foam::polyMesh::polyMesh
00136 (
00137     const IOobject& io,
00138     const Xfer<pointField>& points,
00139     const cellShapeList& cellsAsShapes,
00140     const faceListList& boundaryFaces,
00141     const wordList& boundaryPatchNames,
00142     const wordList& boundaryPatchTypes,
00143     const word& defaultBoundaryPatchName,
00144     const word& defaultBoundaryPatchType,
00145     const wordList& boundaryPatchPhysicalTypes,
00146     const bool syncPar
00147 )
00148 :
00149     objectRegistry(io),
00150     primitiveMesh(),
00151     points_
00152     (
00153         IOobject
00154         (
00155             "points",
00156             instance(),
00157             meshSubDir,
00158             *this,
00159             IOobject::NO_READ,
00160             IOobject::AUTO_WRITE
00161         ),
00162         points
00163     ),
00164     faces_
00165     (
00166         IOobject
00167         (
00168             "faces",
00169             instance(),
00170             meshSubDir,
00171             *this,
00172             IOobject::NO_READ,
00173             IOobject::AUTO_WRITE
00174         ),
00175         0
00176     ),
00177     owner_
00178     (
00179         IOobject
00180         (
00181             "owner",
00182             instance(),
00183             meshSubDir,
00184             *this,
00185             IOobject::NO_READ,
00186             IOobject::AUTO_WRITE
00187         ),
00188         0
00189     ),
00190     neighbour_
00191     (
00192         IOobject
00193         (
00194             "neighbour",
00195             instance(),
00196             meshSubDir,
00197             *this,
00198             IOobject::NO_READ,
00199             IOobject::AUTO_WRITE
00200         ),
00201         0
00202     ),
00203     clearedPrimitives_(false),
00204     boundary_
00205     (
00206         IOobject
00207         (
00208             "boundary",
00209             instance(),
00210             meshSubDir,
00211             *this,
00212             IOobject::NO_READ,
00213             IOobject::AUTO_WRITE
00214         ),
00215         *this,
00216         boundaryFaces.size() + 1    // add room for a default patch
00217     ),
00218     bounds_(points_, syncPar),
00219     geometricD_(Vector<label>::zero),
00220     solutionD_(Vector<label>::zero),
00221     pointZones_
00222     (
00223         IOobject
00224         (
00225             "pointZones",
00226             instance(),
00227             meshSubDir,
00228             *this,
00229             IOobject::NO_READ,
00230             IOobject::NO_WRITE
00231         ),
00232         *this,
00233         0
00234     ),
00235     faceZones_
00236     (
00237         IOobject
00238         (
00239             "faceZones",
00240             instance(),
00241             meshSubDir,
00242             *this,
00243             IOobject::NO_READ,
00244             IOobject::NO_WRITE
00245         ),
00246         *this,
00247         0
00248     ),
00249     cellZones_
00250     (
00251         IOobject
00252         (
00253             "cellZones",
00254             instance(),
00255             meshSubDir,
00256             *this,
00257             IOobject::NO_READ,
00258             IOobject::NO_WRITE
00259         ),
00260         *this,
00261         0
00262     ),
00263     globalMeshDataPtr_(NULL),
00264     moving_(false),
00265     curMotionTimeIndex_(time().timeIndex()),
00266     oldPointsPtr_(NULL)
00267 {
00268     if (debug)
00269     {
00270         Info<<"Constructing polyMesh from cell and boundary shapes." << endl;
00271     }
00272 
00273     // Remove all of the old mesh files if they exist
00274     removeFiles(instance());
00275 
00276     // Calculate the faces of all cells
00277     // Initialise maximum possible numer of mesh faces to 0
00278     label maxFaces = 0;
00279 
00280     // Set up a list of face shapes for each cell
00281     faceListList cellsFaceShapes(cellsAsShapes.size());
00282     cellList cells(cellsAsShapes.size());
00283 
00284     forAll(cellsFaceShapes, cellI)
00285     {
00286         cellsFaceShapes[cellI] = cellsAsShapes[cellI].faces();
00287 
00288         cells[cellI].setSize(cellsFaceShapes[cellI].size());
00289 
00290         // Initialise cells to -1 to flag undefined faces
00291         static_cast<labelList&>(cells[cellI]) = -1;
00292 
00293         // Count maximum possible numer of mesh faces
00294         maxFaces += cellsFaceShapes[cellI].size();
00295     }
00296 
00297     // Set size of faces array to maximum possible number of mesh faces
00298     faces_.setSize(maxFaces);
00299 
00300     // Initialise number of faces to 0
00301     label nFaces = 0;
00302 
00303     // set reference to point-cell addressing
00304     labelListList PointCells = cellShapePointCells(cellsAsShapes);
00305 
00306     bool found = false;
00307 
00308     forAll(cells, cellI)
00309     {
00310         // Note:
00311         // Insertion cannot be done in one go as the faces need to be
00312         // added into the list in the increasing order of neighbour
00313         // cells.  Therefore, all neighbours will be detected first
00314         // and then added in the correct order.
00315 
00316         const faceList& curFaces = cellsFaceShapes[cellI];
00317 
00318         // Record the neighbour cell
00319         labelList neiCells(curFaces.size(), -1);
00320 
00321         // Record the face of neighbour cell
00322         labelList faceOfNeiCell(curFaces.size(), -1);
00323 
00324         label nNeighbours = 0;
00325 
00326         // For all faces ...
00327         forAll(curFaces, faceI)
00328         {
00329             // Skip faces that have already been matched
00330             if (cells[cellI][faceI] >= 0) continue;
00331 
00332             found = false;
00333 
00334             const face& curFace = curFaces[faceI];
00335 
00336             // Get the list of labels
00337             const labelList& curPoints = curFace;
00338 
00339             // For all points
00340             forAll(curPoints, pointI)
00341             {
00342                 // dGget the list of cells sharing this point
00343                 const labelList& curNeighbours =
00344                     PointCells[curPoints[pointI]];
00345 
00346                 // For all neighbours
00347                 forAll(curNeighbours, neiI)
00348                 {
00349                     label curNei = curNeighbours[neiI];
00350 
00351                     // Reject neighbours with the lower label
00352                     if (curNei > cellI)
00353                     {
00354                         // Get the list of search faces
00355                         const faceList& searchFaces = cellsFaceShapes[curNei];
00356 
00357                         forAll(searchFaces, neiFaceI)
00358                         {
00359                             if (searchFaces[neiFaceI] == curFace)
00360                             {
00361                                 // Match!!
00362                                 found = true;
00363 
00364                                 // Record the neighbour cell and face
00365                                 neiCells[faceI] = curNei;
00366                                 faceOfNeiCell[faceI] = neiFaceI;
00367                                 nNeighbours++;
00368 
00369                                 break;
00370                             }
00371                         }
00372                         if (found) break;
00373                     }
00374                     if (found) break;
00375                 }
00376                 if (found) break;
00377             } // End of current points
00378         }  // End of current faces
00379 
00380         // Add the faces in the increasing order of neighbours
00381         for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
00382         {
00383             // Find the lowest neighbour which is still valid
00384             label nextNei = -1;
00385             label minNei = cells.size();
00386 
00387             forAll (neiCells, ncI)
00388             {
00389                 if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
00390                 {
00391                     nextNei = ncI;
00392                     minNei = neiCells[ncI];
00393                 }
00394             }
00395 
00396             if (nextNei > -1)
00397             {
00398                 // Add the face to the list of faces
00399                 faces_[nFaces] = curFaces[nextNei];
00400 
00401                 // Set cell-face and cell-neighbour-face to current face label
00402                 cells[cellI][nextNei] = nFaces;
00403                 cells[neiCells[nextNei]][faceOfNeiCell[nextNei]] = nFaces;
00404 
00405                 // Stop the neighbour from being used again
00406                 neiCells[nextNei] = -1;
00407 
00408                 // Increment number of faces counter
00409                 nFaces++;
00410             }
00411             else
00412             {
00413                 FatalErrorIn
00414                 (
00415                     "polyMesh::polyMesh\n"
00416                     "(\n"
00417                     "    const IOobject&,\n"
00418                     "    const Xfer<pointField>&,\n"
00419                     "    const cellShapeList& cellsAsShapes,\n"
00420                     "    const faceListList& boundaryFaces,\n"
00421                     "    const wordList& boundaryPatchTypes,\n"
00422                     "    const wordList& boundaryPatchNames,\n"
00423                     "    const word& defaultBoundaryPatchType\n"
00424                     ")"
00425                 )   << "Error in internal face insertion"
00426                     << abort(FatalError);
00427             }
00428         }
00429     }
00430 
00431     // Do boundary faces
00432 
00433     labelList patchSizes(boundaryFaces.size(), -1);
00434     labelList patchStarts(boundaryFaces.size(), -1);
00435 
00436     forAll (boundaryFaces, patchI)
00437     {
00438         const faceList& patchFaces = boundaryFaces[patchI];
00439 
00440         labelList curPatchFaceCells =
00441             facePatchFaceCells
00442             (
00443                 patchFaces,
00444                 PointCells,
00445                 cellsFaceShapes,
00446                 patchI
00447             );
00448 
00449         // Grab the start label
00450         label curPatchStart = nFaces;
00451 
00452         forAll (patchFaces, faceI)
00453         {
00454             const face& curFace = patchFaces[faceI];
00455 
00456             const label cellInside = curPatchFaceCells[faceI];
00457 
00458             faces_[nFaces] = curFace;
00459 
00460             // get faces of the cell inside
00461             const faceList& facesOfCellInside = cellsFaceShapes[cellInside];
00462 
00463             bool found = false;
00464 
00465             forAll (facesOfCellInside, cellFaceI)
00466             {
00467                 if (facesOfCellInside[cellFaceI] == curFace)
00468                 {
00469                     if (cells[cellInside][cellFaceI] >= 0)
00470                     {
00471                         FatalErrorIn
00472                         (
00473                             "polyMesh::polyMesh\n"
00474                             "(\n"
00475                             "    const IOobject&,\n"
00476                             "    const Xfer<pointField>&,\n"
00477                             "    const cellShapeList& cellsAsShapes,\n"
00478                             "    const faceListList& boundaryFaces,\n"
00479                             "    const wordList& boundaryPatchTypes,\n"
00480                             "    const wordList& boundaryPatchNames,\n"
00481                             "    const word& defaultBoundaryPatchType\n"
00482                             ")"
00483                         )   << "Trying to specify a boundary face " << curFace
00484                             << " on the face on cell " << cellInside
00485                             << " which is either an internal face or already "
00486                             << "belongs to some other patch.  This is face "
00487                             << faceI << " of patch "
00488                             << patchI << " named "
00489                             << boundaryPatchNames[patchI] << "."
00490                             << abort(FatalError);
00491                     }
00492 
00493                     found = true;
00494 
00495                     cells[cellInside][cellFaceI] = nFaces;
00496 
00497                     break;
00498                 }
00499             }
00500 
00501             if (!found)
00502             {
00503                 FatalErrorIn("polyMesh::polyMesh(... construct from shapes...)")
00504                     << "face " << faceI << " of patch " << patchI
00505                     << " does not seem to belong to cell " << cellInside
00506                     << " which, according to the addressing, "
00507                     << "should be next to it."
00508                     << abort(FatalError);
00509             }
00510 
00511             // increment the counter of faces
00512             nFaces++;
00513         }
00514 
00515         patchSizes[patchI] = nFaces - curPatchStart;
00516         patchStarts[patchI] = curPatchStart;
00517     }
00518 
00519     // Grab "non-existing" faces and put them into a default patch
00520 
00521     label defaultPatchStart = nFaces;
00522 
00523     forAll(cells, cellI)
00524     {
00525         labelList& curCellFaces = cells[cellI];
00526 
00527         forAll(curCellFaces, faceI)
00528         {
00529             if (curCellFaces[faceI] == -1) // "non-existent" face
00530             {
00531                 curCellFaces[faceI] = nFaces;
00532                 faces_[nFaces] = cellsFaceShapes[cellI][faceI];
00533 
00534                 nFaces++;
00535             }
00536         }
00537     }
00538 
00539     // Reset the size of the face list
00540     faces_.setSize(nFaces);
00541 
00542     // Warning: Patches can only be added once the face list is
00543     // completed, as they hold a subList of the face list
00544     forAll (boundaryFaces, patchI)
00545     {
00546         // add the patch to the list
00547         boundary_.set
00548         (
00549             patchI,
00550             polyPatch::New
00551             (
00552                 boundaryPatchTypes[patchI],
00553                 boundaryPatchNames[patchI],
00554                 patchSizes[patchI],
00555                 patchStarts[patchI],
00556                 patchI,
00557                 boundary_
00558             )
00559         );
00560 
00561         if
00562         (
00563             boundaryPatchPhysicalTypes.size()
00564          && boundaryPatchPhysicalTypes[patchI].size()
00565         )
00566         {
00567             boundary_[patchI].physicalType() =
00568                 boundaryPatchPhysicalTypes[patchI];
00569         }
00570     }
00571 
00572     label nAllPatches = boundaryFaces.size();
00573 
00574     if (nFaces > defaultPatchStart)
00575     {
00576         WarningIn("polyMesh::polyMesh(... construct from shapes...)")
00577             << "Found " << nFaces - defaultPatchStart
00578             << " undefined faces in mesh; adding to default patch." << endl;
00579 
00580         boundary_.set
00581         (
00582             nAllPatches,
00583             polyPatch::New
00584             (
00585                 defaultBoundaryPatchType,
00586                 defaultBoundaryPatchName,
00587                 nFaces - defaultPatchStart,
00588                 defaultPatchStart,
00589                 boundary_.size() - 1,
00590                 boundary_
00591             )
00592         );
00593 
00594         nAllPatches++;
00595     }
00596 
00597     // Reset the size of the boundary
00598     boundary_.setSize(nAllPatches);
00599 
00600     // Set the primitive mesh
00601     initMesh(cells);
00602 
00603     if (syncPar)
00604     {
00605         // Calculate topology for the patches (processor-processor comms etc.)
00606         boundary_.updateMesh();
00607 
00608         // Calculate the geometry for the patches (transformation tensors etc.)
00609         boundary_.calcGeometry();
00610     }
00611 
00612     if (debug)
00613     {
00614         if (checkMesh())
00615         {
00616             Info << "Mesh OK" << endl;
00617         }
00618     }
00619 }
00620 
00621 
00622 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines