src/conversion/meshReader/createPolyCells.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 cellPolys
  27     - use pointCells when searching for connectivity
  28     - initialize the cell connectivity with '-1'
  29     - find both cell faces corresponding to the baffles and mark them
  30       to prevent a connection
  31     - standard connectivity checks
  32
  33     - added baffle support
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "meshReader.H"
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 void Foam::meshReader::createPolyCells()
  42 {
  43     // loop through all cell faces and create connectivity. This will produce
  44     // a global face list and will describe all cells as lists of face labels
  45
  46     const faceListList& cFaces = cellFaces();
  47
  48     // count the maximum number of faces and set the size of the cellPolys_
  49     cellPolys_.setSize(cFaces.size());
  50
  51     label maxFaces = 0;
  52
  53     forAll(cellPolys_, cellI)
  54     {
  55         cellPolys_[cellI].setSize(cFaces[cellI].size(), -1);
  56
  57         maxFaces += cFaces[cellI].size();
  58     }
  59
  60     Info<< "Maximum possible number of faces in mesh: " << maxFaces << endl;
  61
  62     meshFaces_.setSize(maxFaces);
  63
  64     // set reference to point-cell addressing
  65     const labelListList& ptCells = pointCells();
  66
  67     // size the baffle lists and initialize to -1
  68     baffleIds_.setSize(baffleFaces_.size());
  69     forAll(baffleIds_, baffleI)
  70     {
  71         baffleIds_[baffleI].setSize(2);
  72     }
  73
  74     // block off baffles first
  75     //
  76     // To prevent internal faces, we'll mark the cell faces
  77     // with negative cell ids (offset by nCells).
  78     // eg,
  79     //    cellI = -(nCells + baffleI)
  80     //
  81     // To distinguish these from the normal '-1' marker, we require
  82     //    cellI = -(nCells + baffleI) < -1
  83     //
  84     // This condition is met provided that nCells > 1.
  85     // ie., baffles require at least 2 volume cells
  86
  87     label baffleOffset = cFaces.size();
  88     forAll(baffleFaces_, baffleI)
  89     {
  90         label cellI = -(baffleOffset + baffleI);
  91         const face& curFace = baffleFaces_[baffleI];
  92
  93         // get the list of labels
  94         const labelList& curPoints = curFace;
  95
  96         // a baffle is a single face - only need to match one face
  97         // get the list of cells sharing this point
  98         const labelList& curNeighbours = ptCells[curPoints[0]];
  99
 100         label nNeighbours = 0;
 101
 102         // For all neighbours
 103         forAll(curNeighbours, neiI)
 104         {
 105             label curNei = curNeighbours[neiI];
 106
 107             // get the list of search faces
 108             const faceList& searchFaces = cFaces[curNei];
 109
 110             forAll(searchFaces, neiFaceI)
 111             {
 112                 int cmp = face::compare(curFace, searchFaces[neiFaceI]);
 113
 114                 if (cmp)
 115                 {
 116                     // maintain baffle orientation
 117                     // side0: baffle normal same as attached face
 118                     // side1: baffle normal opposite from attached face
 119                     //
 120                     label side = 0;
 121                     if (cmp < 0)
 122                     {
 123                         side = 1;
 124                     }
 125
 126 #ifdef DEBUG_FACE_ORDERING
 127                     Info<< "cmp " << cmp << " matched " << curFace
 128                         << " with " << searchFaces[neiFaceI]
 129                         << endl;
 130
 131
 132                     Info<< "match " << baffleI
 133                         << " (" << origCellId_[baffleOffset+baffleI] << ")"
 134                         << " side " << side
 135                         << " against cell " << curNei
 136                         << " face " << neiFaceI
 137                         << " curFace " << curFace[1]
 138                         << " neiFace " << searchFaces[neiFaceI][1]
 139                         << endl;
 140 #endif
 141
 142                     if (baffleIds_[baffleI][side].unused())
 143                     {
 144                         baffleIds_[baffleI][side] = cellFaceIdentifier
 145                         (
 146                             curNei,
 147                             neiFaceI
 148                         );
 149
 150                         nNeighbours++;
 151                     }
 152                     else
 153                     {
 154                         Info<< "multiple matches for side " << side
 155                             << " of baffle " << baffleI
 156                             << " (original cell "
 157                             << origCellId_[baffleOffset+baffleI] << ")"
 158                             << endl;
 159                     }
 160                     break;
 161                 }
 162             }
 163             if (nNeighbours >= 2) break;
 164         }
 165
 166         if (nNeighbours == 2)
 167         {
 168             for (label side = 0; side < nNeighbours; ++side)
 169             {
 170                 label neiCell = baffleIds_[baffleI][side].cell;
 171                 label neiFace = baffleIds_[baffleI][side].face;
 172
 173                 if (baffleIds_[baffleI][side].used())
 174                 {
 175                     cellPolys_[neiCell][neiFace] = cellI;
 176                 }
 177             }
 178         }
 179         else
 180         {
 181             Info<< "drop baffle " << baffleI
 182                 << " (original cell "
 183                 << origCellId_[baffleOffset+baffleI] << ")"
 184                 << " with " << nNeighbours << " neighbours" << endl;
 185
 186             baffleFaces_[baffleI].clear();
 187             baffleIds_[baffleI].clear();
 188         }
 189     }
 190
 191 #ifdef DEBUG_CELLPOLY
 192     Info<< "cellPolys_" << cellPolys_ << endl;
 193     Info<< "baffleFaces_" << baffleFaces_ << endl;
 194     Info<< "baffleIds_"   << baffleIds_ << endl;
 195 #endif
 196
 197     bool found = false;
 198
 199     nInternalFaces_ = 0;
 200
 201     forAll(cFaces, cellI)
 202     {
 203         // Note:
 204         // Insertion cannot be done in one go as the faces need to be
 205         // added into the list in the increasing order of neighbour
 206         // cells.  Therefore, all neighbours will be detected first
 207         // and then added in the correct order.
 208
 209         const faceList& curFaces = cFaces[cellI];
 210
 211         // Record the neighbour cell
 212         labelList neiCells(curFaces.size(), -1);
 213
 214         // Record the face of neighbour cell
 215         labelList faceOfNeiCell(curFaces.size(), -1);
 216
 217         label nNeighbours = 0;
 218
 219         // For all faces ...
 220         forAll(curFaces, faceI)
 221         {
 222             // Skip already matched faces or those tagged by baffles
 223             if (cellPolys_[cellI][faceI] != -1) continue;
 224
 225             found = false;
 226
 227             const face& curFace = curFaces[faceI];
 228
 229             // get the list of labels
 230             const labelList& curPoints = curFace;
 231
 232             // For all points
 233             forAll(curPoints, pointI)
 234             {
 235                 // get the list of cells sharing this point
 236                 const labelList& curNeighbours = ptCells[curPoints[pointI]];
 237
 238                 // For all neighbours
 239                 forAll(curNeighbours, neiI)
 240                 {
 241                     label curNei = curNeighbours[neiI];
 242
 243                     // reject neighbours with the lower label. This should
 244                     // also reject current cell.
 245                     if (curNei > cellI)
 246                     {
 247                         // get the list of search faces
 248                         const faceList& searchFaces = cFaces[curNei];
 249
 250                         forAll(searchFaces, neiFaceI)
 251                         {
 252                             if (searchFaces[neiFaceI] == curFace)
 253                             {
 254                                 // Record the neighbour cell and face
 255                                 neiCells[faceI] = curNei;
 256                                 faceOfNeiCell[faceI] = neiFaceI;
 257                                 nNeighbours++;
 258 #ifdef DEBUG_FACE_ORDERING
 259                                 Info<< " cell " << cellI
 260                                     << " face " << faceI
 261                                     << " point " << pointI
 262                                     << " nei " << curNei
 263                                     << " neiFace " << neiFaceI
 264                                     << endl;
 265 #endif
 266                                 found = true;
 267                                 break;
 268                             }
 269                         }
 270                         if (found) break;
 271                     }
 272                     if (found) break;
 273                 }
 274                 if (found) break;
 275             } // End of current points
 276         } // End of current faces
 277
 278         // Add the faces in the increasing order of neighbours
 279         for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
 280         {
 281             // Find the lowest neighbour which is still valid
 282             label nextNei = -1;
 283             label minNei = cellPolys_.size();
 284
 285             forAll(neiCells, ncI)
 286             {
 287                 if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
 288                 {
 289                     nextNei = ncI;
 290                     minNei = neiCells[ncI];
 291                 }
 292             }
 293
 294             if (nextNei > -1)
 295             {
 296                 // Add the face to the list of faces
 297                 meshFaces_[nInternalFaces_] = curFaces[nextNei];
 298
 299                 // Mark for owner
 300                 cellPolys_[cellI][nextNei] = nInternalFaces_;
 301
 302                 // Mark for neighbour
 303                 cellPolys_[neiCells[nextNei]][faceOfNeiCell[nextNei]] =
 304                     nInternalFaces_;
 305
 306                 // Stop the neighbour from being used again
 307                 neiCells[nextNei] = -1;
 308
 309                 // Increment number of faces counter
 310                 nInternalFaces_++;
 311             }
 312             else
 313             {
 314               FatalErrorIn("meshReader::createPolyCells()")
 315                   << "Error in internal face insertion"
 316                   << abort(FatalError);
 317             }
 318         }
 319     }
 320
 321 #ifdef DEBUG_CELLPOLY
 322     Info<< "cellPolys = " << cellPolys_ << endl;
 323 #endif
 324
 325     // don't reset the size of internal faces, because more faces will be
 326     // added in createPolyBoundary()
 327 }
 328
 329
 330 // ************************************************************************* //