src/conversion/meshReader/createPolyCells.C

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