src/finiteVolume/fvMesh/extendedStencil/cellToFace/fullStencils/cellToFaceStencil.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "cellToFaceStencil.H"
  28 #include "syncTools.H"
  29 #include "SortableList.H"
  30 #include "emptyPolyPatch.H"
  31
  32 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  33
  34 // Merge two list and guarantee global0,global1 are first.
  35 void Foam::cellToFaceStencil::merge
  36 (
  37     const label global0,
  38     const label global1,
  39     const labelList& listA,
  40     labelList& listB
  41 )
  42 {
  43     sort(listB);
  44
  45     // See if global0, global1 already present in listB
  46     label nGlobalInsert = 0;
  47
  48     if (global0 != -1)
  49     {
  50         label index0 = findSortedIndex(listB, global0);
  51         if (index0 == -1)
  52         {
  53             nGlobalInsert++;
  54         }
  55     }
  56
  57     if (global1 != -1)
  58     {
  59         label index1 = findSortedIndex(listB, global1);
  60         if (index1 == -1)
  61         {
  62             nGlobalInsert++;
  63         }
  64     }
  65
  66
  67     // For all in listA see if they are present
  68     label nInsert = 0;
  69
  70     forAll(listA, i)
  71     {
  72         label elem = listA[i];
  73
  74         if (elem != global0 && elem != global1)
  75         {
  76             if (findSortedIndex(listB, elem) == -1)
  77             {
  78                 nInsert++;
  79             }
  80         }
  81     }
  82
  83     // Extend B with nInsert and whether global0,global1 need to be inserted.
  84     labelList result(listB.size() + nGlobalInsert + nInsert);
  85
  86     label resultI = 0;
  87
  88     // Insert global0,1 first
  89     if (global0 != -1)
  90     {
  91         result[resultI++] = global0;
  92     }
  93     if (global1 != -1)
  94     {
  95         result[resultI++] = global1;
  96     }
  97
  98
  99     // Insert listB
 100     forAll(listB, i)
 101     {
 102         label elem = listB[i];
 103
 104         if (elem != global0 && elem != global1)
 105         {
 106             result[resultI++] = elem;
 107         }
 108     }
 109
 110
 111     // Insert listA
 112     forAll(listA, i)
 113     {
 114         label elem = listA[i];
 115
 116         if (elem != global0 && elem != global1)
 117         {
 118             if (findSortedIndex(listB, elem) == -1)
 119             {
 120                 result[resultI++] = elem;
 121             }
 122         }
 123     }
 124
 125     if (resultI != result.size())
 126     {
 127         FatalErrorIn("cellToFaceStencil::merge(..)")
 128             << "problem" << abort(FatalError);
 129     }
 130
 131     listB.transfer(result);
 132 }
 133
 134
 135 // Merge two list and guarantee globalI is first.
 136 void Foam::cellToFaceStencil::merge
 137 (
 138     const label globalI,
 139     const labelList& pGlobals,
 140     labelList& cCells
 141 )
 142 {
 143     labelHashSet set;
 144     forAll(cCells, i)
 145     {
 146         if (cCells[i] != globalI)
 147         {
 148             set.insert(cCells[i]);
 149         }
 150     }
 151
 152     forAll(pGlobals, i)
 153     {
 154         if (pGlobals[i] != globalI)
 155         {
 156             set.insert(pGlobals[i]);
 157         }
 158     }
 159
 160     cCells.setSize(set.size()+1);
 161     label n = 0;
 162     cCells[n++] = globalI;
 163
 164     forAllConstIter(labelHashSet, set, iter)
 165     {
 166         cCells[n++] = iter.key();
 167     }
 168 }
 169
 170
 171 void Foam::cellToFaceStencil::validBoundaryFaces(boolList& isValidBFace) const
 172 {
 173     const polyBoundaryMesh& patches = mesh().boundaryMesh();
 174
 175     isValidBFace.setSize(mesh().nFaces()-mesh().nInternalFaces(), true);
 176
 177     forAll(patches, patchI)
 178     {
 179         const polyPatch& pp = patches[patchI];
 180
 181         if (pp.coupled() || isA<emptyPolyPatch>(pp))
 182         {
 183             label bFaceI = pp.start()-mesh().nInternalFaces();
 184             forAll(pp, i)
 185             {
 186                 isValidBFace[bFaceI++] = false;
 187             }
 188         }
 189     }
 190 }
 191
 192
 193 Foam::autoPtr<Foam::indirectPrimitivePatch>
 194 Foam::cellToFaceStencil::allCoupledFacesPatch() const
 195 {
 196     const polyBoundaryMesh& patches = mesh().boundaryMesh();
 197
 198     label nCoupled = 0;
 199
 200     forAll(patches, patchI)
 201     {
 202         const polyPatch& pp = patches[patchI];
 203
 204         if (pp.coupled())
 205         {
 206             nCoupled += pp.size();
 207         }
 208     }
 209     labelList coupledFaces(nCoupled);
 210     nCoupled = 0;
 211
 212     forAll(patches, patchI)
 213     {
 214         const polyPatch& pp = patches[patchI];
 215
 216         if (pp.coupled())
 217         {
 218             label faceI = pp.start();
 219
 220             forAll(pp, i)
 221             {
 222                 coupledFaces[nCoupled++] = faceI++;
 223             }
 224         }
 225     }
 226
 227     return autoPtr<indirectPrimitivePatch>
 228     (
 229         new indirectPrimitivePatch
 230         (
 231             IndirectList<face>
 232             (
 233                 mesh().faces(),
 234                 coupledFaces
 235             ),
 236             mesh().points()
 237         )
 238     );
 239 }
 240
 241
 242 void Foam::cellToFaceStencil::unionEqOp::operator()
 243 (
 244     labelList& x,
 245     const labelList& y
 246 ) const
 247 {
 248     if (y.size())
 249     {
 250         if (x.empty())
 251         {
 252             x = y;
 253         }
 254         else
 255         {
 256             labelHashSet set(x);
 257             forAll(y, i)
 258             {
 259                 set.insert(y[i]);
 260             }
 261             x = set.toc();
 262         }
 263     }
 264 }
 265
 266
 267 void Foam::cellToFaceStencil::insertFaceCells
 268 (
 269     const label exclude0,
 270     const label exclude1,
 271     const boolList& isValidBFace,
 272     const labelList& faceLabels,
 273     labelHashSet& globals
 274 ) const
 275 {
 276     const labelList& own = mesh().faceOwner();
 277     const labelList& nei = mesh().faceNeighbour();
 278
 279     forAll(faceLabels, i)
 280     {
 281         label faceI = faceLabels[i];
 282
 283         label globalOwn = globalNumbering().toGlobal(own[faceI]);
 284         if (globalOwn != exclude0 && globalOwn != exclude1)
 285         {
 286             globals.insert(globalOwn);
 287         }
 288
 289         if (mesh().isInternalFace(faceI))
 290         {
 291             label globalNei = globalNumbering().toGlobal(nei[faceI]);
 292             if (globalNei != exclude0 && globalNei != exclude1)
 293             {
 294                 globals.insert(globalNei);
 295             }
 296         }
 297         else
 298         {
 299             label bFaceI = faceI-mesh().nInternalFaces();
 300
 301             if (isValidBFace[bFaceI])
 302             {
 303                 label globalI = globalNumbering().toGlobal
 304                 (
 305                     mesh().nCells()
 306                   + bFaceI
 307                 );
 308
 309                 if (globalI != exclude0 && globalI != exclude1)
 310                 {
 311                     globals.insert(globalI);
 312                 }
 313             }
 314         }
 315     }
 316 }
 317
 318
 319 Foam::labelList Foam::cellToFaceStencil::calcFaceCells
 320 (
 321     const boolList& isValidBFace,
 322     const labelList& faceLabels,
 323     labelHashSet& globals
 324 ) const
 325 {
 326     globals.clear();
 327
 328     insertFaceCells
 329     (
 330         -1,
 331         -1,
 332         isValidBFace,
 333         faceLabels,
 334         globals
 335     );
 336
 337     return globals.toc();
 338 }
 339
 340
 341 // Calculates per face a list of global cell/face indices.
 342 void Foam::cellToFaceStencil::calcFaceStencil
 343 (
 344     const labelListList& globalCellCells,
 345     labelListList& faceStencil
 346 ) const
 347 {
 348     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 349     const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
 350     const labelList& own = mesh_.faceOwner();
 351     const labelList& nei = mesh_.faceNeighbour();
 352
 353
 354     // Determine neighbouring global cell Cells
 355     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 356
 357     labelListList neiGlobalCellCells(nBnd);
 358     forAll(patches, patchI)
 359     {
 360         const polyPatch& pp = patches[patchI];
 361
 362         if (pp.coupled())
 363         {
 364             label faceI = pp.start();
 365
 366             forAll(pp, i)
 367             {
 368                 neiGlobalCellCells[faceI-mesh_.nInternalFaces()] =
 369                     globalCellCells[own[faceI]];
 370                 faceI++;
 371             }
 372         }
 373     }
 374     syncTools::swapBoundaryFaceList(mesh_, neiGlobalCellCells, false);
 375
 376
 377
 378     // Construct stencil in global numbering
 379     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 380
 381     faceStencil.setSize(mesh_.nFaces());
 382
 383     labelHashSet faceStencilSet;
 384
 385     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 386     {
 387         faceStencilSet.clear();
 388
 389         const labelList& ownCCells = globalCellCells[own[faceI]];
 390         label globalOwn = ownCCells[0];
 391         // Insert cellCells
 392         forAll(ownCCells, i)
 393         {
 394             faceStencilSet.insert(ownCCells[i]);
 395         }
 396
 397         const labelList& neiCCells = globalCellCells[nei[faceI]];
 398         label globalNei = neiCCells[0];
 399         // Insert cellCells
 400         forAll(neiCCells, i)
 401         {
 402             faceStencilSet.insert(neiCCells[i]);
 403         }
 404
 405         // Guarantee owner first, neighbour second.
 406         faceStencil[faceI].setSize(faceStencilSet.size());
 407         label n = 0;
 408         faceStencil[faceI][n++] = globalOwn;
 409         faceStencil[faceI][n++] = globalNei;
 410         forAllConstIter(labelHashSet, faceStencilSet, iter)
 411         {
 412             if (iter.key() != globalOwn && iter.key() != globalNei)
 413             {
 414                 faceStencil[faceI][n++] = iter.key();
 415             }
 416         }
 417         //Pout<< "internalface:" << faceI << " toc:" << faceStencilSet.toc()
 418         //    << " faceStencil:" << faceStencil[faceI] << endl;
 419     }
 420     forAll(patches, patchI)
 421     {
 422         const polyPatch& pp = patches[patchI];
 423         label faceI = pp.start();
 424
 425         if (pp.coupled())
 426         {
 427             forAll(pp, i)
 428             {
 429                 faceStencilSet.clear();
 430
 431                 const labelList& ownCCells = globalCellCells[own[faceI]];
 432                 label globalOwn = ownCCells[0];
 433                 forAll(ownCCells, i)
 434                 {
 435                     faceStencilSet.insert(ownCCells[i]);
 436                 }
 437
 438                 // And the neighbours of the coupled cell
 439                 const labelList& neiCCells =
 440                     neiGlobalCellCells[faceI-mesh_.nInternalFaces()];
 441                 label globalNei = neiCCells[0];
 442                 forAll(neiCCells, i)
 443                 {
 444                     faceStencilSet.insert(neiCCells[i]);
 445                 }
 446
 447                 // Guarantee owner first, neighbour second.
 448                 faceStencil[faceI].setSize(faceStencilSet.size());
 449                 label n = 0;
 450                 faceStencil[faceI][n++] = globalOwn;
 451                 faceStencil[faceI][n++] = globalNei;
 452                 forAllConstIter(labelHashSet, faceStencilSet, iter)
 453                 {
 454                     if (iter.key() != globalOwn && iter.key() != globalNei)
 455                     {
 456                         faceStencil[faceI][n++] = iter.key();
 457                     }
 458                 }
 459
 460                 //Pout<< "coupledface:" << faceI
 461                 //    << " toc:" << faceStencilSet.toc()
 462                 //    << " faceStencil:" << faceStencil[faceI] << endl;
 463
 464                 faceI++;
 465             }
 466         }
 467         else if (!isA<emptyPolyPatch>(pp))
 468         {
 469             forAll(pp, i)
 470             {
 471                 faceStencilSet.clear();
 472
 473                 const labelList& ownCCells = globalCellCells[own[faceI]];
 474                 label globalOwn = ownCCells[0];
 475                 forAll(ownCCells, i)
 476                 {
 477                     faceStencilSet.insert(ownCCells[i]);
 478                 }
 479
 480                 // Guarantee owner first
 481                 faceStencil[faceI].setSize(faceStencilSet.size());
 482                 label n = 0;
 483                 faceStencil[faceI][n++] = globalOwn;
 484                 forAllConstIter(labelHashSet, faceStencilSet, iter)
 485                 {
 486                     if (iter.key() != globalOwn)
 487                     {
 488                         faceStencil[faceI][n++] = iter.key();
 489                     }
 490                 }
 491
 492                 //Pout<< "boundaryface:" << faceI
 493                 //    << " toc:" << faceStencilSet.toc()
 494                 //    << " faceStencil:" << faceStencil[faceI] << endl;
 495
 496                 faceI++;
 497             }
 498         }
 499     }
 500 }
 501
 502
 503 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 504
 505 Foam::cellToFaceStencil::cellToFaceStencil(const polyMesh& mesh)
 506 :
 507     mesh_(mesh),
 508     globalNumbering_(mesh_.nCells()+mesh_.nFaces()-mesh_.nInternalFaces())
 509 {}
 510
 511
 512 // ************************************************************************* //