src/finiteVolume/fvMesh/extendedStencil/cellToFace/fullStencils/FECCellToFaceStencil.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "FECCellToFaceStencil.H"
  27 #include "syncTools.H"
  28 #include "emptyPolyPatch.H"
  29 #include "dummyTransform.H"
  30
  31 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  32
  33 // Calculates per edge the neighbour data (= edgeCells)
  34 void Foam::FECCellToFaceStencil::calcEdgeBoundaryData
  35 (
  36     const boolList& isValidBFace,
  37     const labelList& boundaryEdges,
  38     EdgeMap<labelList>& neiGlobal
  39 ) const
  40 {
  41     neiGlobal.resize(2*boundaryEdges.size());
  42
  43     labelHashSet edgeGlobals;
  44
  45     forAll(boundaryEdges, i)
  46     {
  47         label edgeI = boundaryEdges[i];
  48
  49         neiGlobal.insert
  50         (
  51             mesh().edges()[edgeI],
  52             calcFaceCells
  53             (
  54                 isValidBFace,
  55                 mesh().edgeFaces(edgeI),
  56                 edgeGlobals
  57             )
  58         );
  59     }
  60
  61     syncTools::syncEdgeMap(mesh(), neiGlobal, unionEqOp(), dummyTransform());
  62 }
  63
  64
  65 // Calculates per face the edge connected data (= cell or boundary in global
  66 // numbering).
  67 void Foam::FECCellToFaceStencil::calcFaceStencil
  68 (
  69     labelListList& faceStencil
  70 ) const
  71 {
  72     const polyBoundaryMesh& patches = mesh().boundaryMesh();
  73     const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
  74     const labelList& own = mesh().faceOwner();
  75     const labelList& nei = mesh().faceNeighbour();
  76
  77
  78
  79     // Determine neighbouring global cell
  80     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  81
  82     labelList neiGlobalCell(nBnd);
  83     forAll(patches, patchI)
  84     {
  85         const polyPatch& pp = patches[patchI];
  86
  87         if (pp.coupled())
  88         {
  89             label faceI = pp.start();
  90
  91             forAll(pp, i)
  92             {
  93                 neiGlobalCell[faceI-mesh().nInternalFaces()] =
  94                     globalNumbering().toGlobal(own[faceI]);
  95                 faceI++;
  96             }
  97         }
  98     }
  99     syncTools::swapBoundaryFaceList(mesh(), neiGlobalCell);
 100
 101
 102
 103     // Determine on coupled edges the edge cells on the other side
 104     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 105
 106     // Calculate edges on coupled patches
 107     labelList boundaryEdges
 108     (
 109         allCoupledFacesPatch()().meshEdges
 110         (
 111             mesh().edges(),
 112             mesh().pointEdges()
 113         )
 114     );
 115
 116     // Mark boundary faces to be included in stencil (i.e. not coupled or empty)
 117     boolList isValidBFace;
 118     validBoundaryFaces(isValidBFace);
 119
 120     // Swap edgeCells for coupled edges. Note: use EdgeMap for now since we've
 121     // got syncTools::syncEdgeMap for those. Should be replaced with Map and
 122     // syncTools functionality to handle those.
 123     EdgeMap<labelList> neiGlobal;
 124     calcEdgeBoundaryData
 125     (
 126         isValidBFace,
 127         boundaryEdges,
 128         neiGlobal
 129     );
 130
 131
 132
 133     // Construct stencil in global numbering
 134     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 135
 136     faceStencil.setSize(mesh().nFaces());
 137
 138     // Do coupled edges first
 139
 140     forAll(boundaryEdges, i)
 141     {
 142         label edgeI = boundaryEdges[i];
 143
 144         const labelList& eGlobals = neiGlobal[mesh().edges()[edgeI]];
 145
 146         // Distribute to all edgeFaces
 147         const labelList& eFaces = mesh().edgeFaces(edgeI);
 148
 149         forAll(eFaces, j)
 150         {
 151             label faceI = eFaces[j];
 152
 153             // Insert eGlobals into faceStencil.
 154             merge(-1, -1, eGlobals, faceStencil[faceI]);
 155         }
 156     }
 157     neiGlobal.clear();
 158
 159
 160     // Do remaining edges by looping over all faces
 161
 162     // Work arrays
 163     DynamicList<label> fEdgesSet;
 164     DynamicList<label> eFacesSet;
 165     labelHashSet faceStencilSet;
 166
 167     for (label faceI = 0; faceI < mesh().nInternalFaces(); faceI++)
 168     {
 169         label globalOwn = globalNumbering().toGlobal(own[faceI]);
 170         label globalNei = globalNumbering().toGlobal(nei[faceI]);
 171
 172         // Convert any existing faceStencil (from coupled edges) into
 173         // set and operate on this.
 174
 175         faceStencilSet.clear();
 176
 177         // Insert all but global owner and neighbour
 178         forAll(faceStencil[faceI], i)
 179         {
 180             label globalI = faceStencil[faceI][i];
 181             if (globalI != globalOwn && globalI != globalNei)
 182             {
 183                 faceStencilSet.insert(globalI);
 184             }
 185         }
 186         faceStencil[faceI].clear();
 187
 188         // Collect all edge connected (internal) cells
 189         const labelList& fEdges = mesh().faceEdges(faceI, fEdgesSet);
 190
 191         forAll(fEdges, i)
 192         {
 193             label edgeI = fEdges[i];
 194
 195             insertFaceCells
 196             (
 197                 globalOwn,
 198                 globalNei,
 199                 isValidBFace,
 200                 mesh().edgeFaces(edgeI, eFacesSet),
 201                 faceStencilSet
 202             );
 203         }
 204
 205         // Extract, guarantee owner first, neighbour second.
 206         faceStencil[faceI].setSize(faceStencilSet.size()+2);
 207         label n = 0;
 208         faceStencil[faceI][n++] = globalOwn;
 209         faceStencil[faceI][n++] = globalNei;
 210         forAllConstIter(labelHashSet, faceStencilSet, iter)
 211         {
 212             if (iter.key() == globalOwn || iter.key() == globalNei)
 213             {
 214                 FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
 215                     << "problem:" << faceStencilSet
 216                     << abort(FatalError);
 217             }
 218             faceStencil[faceI][n++] = iter.key();
 219         }
 220     }
 221     forAll(patches, patchI)
 222     {
 223         const polyPatch& pp = patches[patchI];
 224         label faceI = pp.start();
 225
 226         if (pp.coupled())
 227         {
 228             forAll(pp, i)
 229             {
 230                 label globalOwn = globalNumbering().toGlobal(own[faceI]);
 231                 label globalNei = neiGlobalCell[faceI-mesh().nInternalFaces()];
 232
 233                 // Convert any existing faceStencil (from coupled edges) into
 234                 // set and operate on this.
 235
 236                 faceStencilSet.clear();
 237
 238                 // Insert all but global owner and neighbour
 239                 forAll(faceStencil[faceI], i)
 240                 {
 241                     label globalI = faceStencil[faceI][i];
 242                     if (globalI != globalOwn && globalI != globalNei)
 243                     {
 244                         faceStencilSet.insert(globalI);
 245                     }
 246                 }
 247                 faceStencil[faceI].clear();
 248
 249                 // Collect all edge connected (internal) cells
 250                 const labelList& fEdges = mesh().faceEdges(faceI, fEdgesSet);
 251
 252                 forAll(fEdges, i)
 253                 {
 254                     label edgeI = fEdges[i];
 255
 256                     insertFaceCells
 257                     (
 258                         globalOwn,
 259                         globalNei,
 260                         isValidBFace,
 261                         mesh().edgeFaces(edgeI, eFacesSet),
 262                         faceStencilSet
 263                     );
 264                 }
 265
 266                 // Extract, guarantee owner first, neighbour second.
 267                 faceStencil[faceI].setSize(faceStencilSet.size()+2);
 268                 label n = 0;
 269                 faceStencil[faceI][n++] = globalOwn;
 270                 faceStencil[faceI][n++] = globalNei;
 271                 forAllConstIter(labelHashSet, faceStencilSet, iter)
 272                 {
 273                     if (iter.key() == globalOwn || iter.key() == globalNei)
 274                     {
 275                         FatalErrorIn
 276                         (
 277                             "FECCellToFaceStencil::calcFaceStencil(..)"
 278                         )   << "problem:" << faceStencilSet
 279                             << abort(FatalError);
 280                     }
 281                     faceStencil[faceI][n++] = iter.key();
 282                 }
 283
 284                 if (n != faceStencil[faceI].size())
 285                 {
 286                     FatalErrorIn("problem") << "n:" << n
 287                         << " size:" << faceStencil[faceI].size()
 288                         << abort(FatalError);
 289                 }
 290
 291                 faceI++;
 292             }
 293         }
 294         else if (!isA<emptyPolyPatch>(pp))
 295         {
 296             forAll(pp, i)
 297             {
 298                 label globalOwn = globalNumbering().toGlobal(own[faceI]);
 299
 300                 // Convert any existing faceStencil (from coupled edges) into
 301                 // set and operate on this.
 302
 303                 faceStencilSet.clear();
 304
 305                 // Insert all but global owner and neighbour
 306                 forAll(faceStencil[faceI], i)
 307                 {
 308                     label globalI = faceStencil[faceI][i];
 309                     if (globalI != globalOwn)
 310                     {
 311                         faceStencilSet.insert(globalI);
 312                     }
 313                 }
 314                 faceStencil[faceI].clear();
 315
 316                 // Collect all edge connected (internal) cells
 317                 const labelList& fEdges = mesh().faceEdges(faceI, fEdgesSet);
 318
 319                 forAll(fEdges, i)
 320                 {
 321                     label edgeI = fEdges[i];
 322
 323                     insertFaceCells
 324                     (
 325                         globalOwn,
 326                         -1,
 327                         isValidBFace,
 328                         mesh().edgeFaces(edgeI, eFacesSet),
 329                         faceStencilSet
 330                     );
 331                 }
 332
 333                 // Extract, guarantee owner first, neighbour second.
 334                 faceStencil[faceI].setSize(faceStencilSet.size()+1);
 335                 label n = 0;
 336                 faceStencil[faceI][n++] = globalOwn;
 337                 forAllConstIter(labelHashSet, faceStencilSet, iter)
 338                 {
 339                     if (iter.key() == globalOwn)
 340                     {
 341                         FatalErrorIn
 342                         (
 343                             "FECCellToFaceStencil::calcFaceStencil(..)"
 344                         )   << "problem:" << faceStencilSet
 345                             << abort(FatalError);
 346                     }
 347                     faceStencil[faceI][n++] = iter.key();
 348                 }
 349
 350                 faceI++;
 351             }
 352         }
 353     }
 354
 355
 356     for (label faceI = 0; faceI < mesh().nInternalFaces(); faceI++)
 357     {
 358         label globalOwn = globalNumbering().toGlobal(own[faceI]);
 359         if (faceStencil[faceI][0] != globalOwn)
 360         {
 361             FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
 362                 << "problem:" << faceStencil[faceI]
 363                 << " globalOwn:" << globalOwn
 364                 << abort(FatalError);
 365         }
 366         label globalNei = globalNumbering().toGlobal(nei[faceI]);
 367         if (faceStencil[faceI][1] != globalNei)
 368         {
 369             FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
 370                 << "problem:" << faceStencil[faceI]
 371                 << " globalNei:" << globalNei
 372                 << abort(FatalError);
 373         }
 374     }
 375
 376
 377     forAll(patches, patchI)
 378     {
 379         const polyPatch& pp = patches[patchI];
 380
 381         if (pp.coupled())
 382         {
 383             forAll(pp, i)
 384             {
 385                 label faceI = pp.start()+i;
 386
 387                 label globalOwn = globalNumbering().toGlobal(own[faceI]);
 388                 if (faceStencil[faceI][0] != globalOwn)
 389                 {
 390                     FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
 391                         << "problem:" << faceStencil[faceI]
 392                         << " globalOwn:" << globalOwn
 393                         << abort(FatalError);
 394                 }
 395                 label globalNei = neiGlobalCell[faceI-mesh().nInternalFaces()];
 396                 if (faceStencil[faceI][1] != globalNei)
 397                 {
 398                     FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
 399                         << "problem:" << faceStencil[faceI]
 400                         << " globalNei:" << globalNei
 401                         << abort(FatalError);
 402                 }
 403             }
 404         }
 405         else if (!isA<emptyPolyPatch>(pp))
 406         {
 407             forAll(pp, i)
 408             {
 409                 label faceI = pp.start()+i;
 410
 411                 label globalOwn = globalNumbering().toGlobal(own[faceI]);
 412                 if (faceStencil[faceI][0] != globalOwn)
 413                 {
 414                     FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
 415                         << "problem:" << faceStencil[faceI]
 416                         << " globalOwn:" << globalOwn
 417                         << abort(FatalError);
 418                 }
 419             }
 420         }
 421     }
 422 }
 423
 424
 425 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 426
 427 Foam::FECCellToFaceStencil::FECCellToFaceStencil(const polyMesh& mesh)
 428 :
 429     cellToFaceStencil(mesh)
 430 {
 431     // Calculate per face the (edge) connected cells (in global numbering)
 432     labelListList faceStencil;
 433     calcFaceStencil(faceStencil);
 434
 435     // Transfer to *this
 436     transfer(faceStencil);
 437 }
 438
 439
 440 // ************************************************************************* //