src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.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 "extendedUpwindCellToFaceStencil.H"
  27 #include "cellToFaceStencil.H"
  28 #include "syncTools.H"
  29 #include "SortableList.H"
  30 #include "dummyTransform.H"
  31
  32 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  33
  34 void Foam::extendedUpwindCellToFaceStencil::selectOppositeFaces
  35 (
  36     const boolList& nonEmptyFace,
  37     const scalar minOpposedness,
  38     const label faceI,
  39     const label cellI,
  40     DynamicList<label>& oppositeFaces
  41 ) const
  42 {
  43     const vectorField& areas = mesh_.faceAreas();
  44     const labelList& own = mesh_.faceOwner();
  45     const cell& cFaces = mesh_.cells()[cellI];
  46
  47     SortableList<scalar> opposedness(cFaces.size(), -GREAT);
  48
  49     // Pick up all the faces that oppose this one.
  50     forAll(cFaces, i)
  51     {
  52         label otherFaceI = cFaces[i];
  53
  54         if (otherFaceI != faceI && nonEmptyFace[otherFaceI])
  55         {
  56             if ((own[otherFaceI] == cellI) == (own[faceI] == cellI))
  57             {
  58                 opposedness[i] = -(areas[otherFaceI] & areas[faceI]);
  59             }
  60             else
  61             {
  62                 opposedness[i] = (areas[otherFaceI] & areas[faceI]);
  63             }
  64         }
  65     }
  66
  67     label sz = opposedness.size();
  68
  69     oppositeFaces.clear();
  70
  71     scalar myAreaSqr = magSqr(areas[faceI]);
  72
  73     if (myAreaSqr > VSMALL)
  74     {
  75         forAll(opposedness, i)
  76         {
  77             opposedness[i] /= myAreaSqr;
  78         }
  79         // Sort in incrementing order
  80         opposedness.sort();
  81
  82         // Pick largest no matter what
  83         oppositeFaces.append(cFaces[opposedness.indices()[sz-1]]);
  84
  85         for (label i = sz-2; i >= 0; --i)
  86         {
  87             if (opposedness[i] < minOpposedness)
  88             {
  89                 break;
  90             }
  91             oppositeFaces.append(cFaces[opposedness.indices()[i]]);
  92         }
  93     }
  94     else
  95     {
  96         // Sort in incrementing order
  97         opposedness.sort();
  98
  99         // Tiny face. Do what?
 100         // Pick largest no matter what
 101         oppositeFaces.append(cFaces[opposedness.indices()[sz-1]]);
 102     }
 103 }
 104
 105
 106 void Foam::extendedUpwindCellToFaceStencil::transportStencil
 107 (
 108     const boolList& nonEmptyFace,
 109     const labelListList& faceStencil,
 110     const scalar minOpposedness,
 111     const label faceI,
 112     const label cellI,
 113     const bool stencilHasNeighbour,
 114
 115     DynamicList<label>& oppositeFaces,
 116     labelHashSet& faceStencilSet,
 117     labelList& transportedStencil
 118 ) const
 119 {
 120     label globalOwn = faceStencil[faceI][0];
 121     label globalNei = -1;
 122     if (stencilHasNeighbour && faceStencil[faceI].size() >= 2)
 123     {
 124         globalNei = faceStencil[faceI][1];
 125     }
 126
 127
 128     selectOppositeFaces
 129     (
 130         nonEmptyFace,
 131         minOpposedness,
 132         faceI,
 133         cellI,
 134         oppositeFaces
 135     );
 136
 137     // Collect all stencils of oppositefaces
 138     faceStencilSet.clear();
 139     forAll(oppositeFaces, i)
 140     {
 141         const labelList& fStencil = faceStencil[oppositeFaces[i]];
 142
 143         forAll(fStencil, j)
 144         {
 145             label globalI = fStencil[j];
 146
 147             if (globalI != globalOwn && globalI != globalNei)
 148             {
 149                 faceStencilSet.insert(globalI);
 150             }
 151         }
 152     }
 153
 154     // Add my owner and neighbour first.
 155     if (stencilHasNeighbour)
 156     {
 157         transportedStencil.setSize(faceStencilSet.size()+2);
 158         label n = 0;
 159         transportedStencil[n++] = globalOwn;
 160         transportedStencil[n++] = globalNei;
 161
 162         forAllConstIter(labelHashSet, faceStencilSet, iter)
 163         {
 164             if (iter.key() != globalOwn && iter.key() != globalNei)
 165             {
 166                 transportedStencil[n++] = iter.key();
 167             }
 168         }
 169         if (n != transportedStencil.size())
 170         {
 171             FatalErrorIn
 172             (
 173                 "extendedUpwindCellToFaceStencil::transportStencil(..)"
 174             )   << "problem:" << faceStencilSet
 175                 << abort(FatalError);
 176         }
 177     }
 178     else
 179     {
 180         transportedStencil.setSize(faceStencilSet.size()+1);
 181         label n = 0;
 182         transportedStencil[n++] = globalOwn;
 183
 184         forAllConstIter(labelHashSet, faceStencilSet, iter)
 185         {
 186             if (iter.key() != globalOwn)
 187             {
 188                 transportedStencil[n++] = iter.key();
 189             }
 190         }
 191         if (n != transportedStencil.size())
 192         {
 193             FatalErrorIn
 194             (
 195                 "extendedUpwindCellToFaceStencil::transportStencil(..)"
 196             )   << "problem:" << faceStencilSet
 197                 << abort(FatalError);
 198         }
 199     }
 200 }
 201
 202
 203 void Foam::extendedUpwindCellToFaceStencil::transportStencils
 204 (
 205     const labelListList& faceStencil,
 206     const scalar minOpposedness,
 207     labelListList& ownStencil,
 208     labelListList& neiStencil
 209 )
 210 {
 211     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 212     const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
 213     const labelList& own = mesh_.faceOwner();
 214     const labelList& nei = mesh_.faceNeighbour();
 215
 216     // Work arrays
 217     DynamicList<label> oppositeFaces;
 218     labelHashSet faceStencilSet;
 219
 220
 221     // For quick detection of empty faces
 222     boolList nonEmptyFace(mesh_.nFaces(), true);
 223     forAll(patches, patchI)
 224     {
 225         const polyPatch& pp = patches[patchI];
 226
 227         if (isA<emptyPolyPatch>(pp))
 228         {
 229             label faceI = pp.start();
 230             forAll(pp, i)
 231             {
 232                 nonEmptyFace[faceI++] = false;
 233             }
 234         }
 235     }
 236
 237
 238     // Do the owner side
 239     // ~~~~~~~~~~~~~~~~~
 240     // stencil is synchronised at entry so no need to swap.
 241
 242     ownStencil.setSize(mesh_.nFaces());
 243
 244     // Internal faces
 245     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 246     {
 247         // Get stencil as owner + neighbour + stencil from 'opposite' faces
 248         transportStencil
 249         (
 250             nonEmptyFace,
 251             faceStencil,
 252             minOpposedness,
 253             faceI,
 254             own[faceI],
 255             true,                   //stencilHasNeighbour
 256             oppositeFaces,
 257             faceStencilSet,
 258             ownStencil[faceI]
 259         );
 260     }
 261     // Boundary faces
 262     forAll(patches, patchI)
 263     {
 264         const polyPatch& pp = patches[patchI];
 265         label faceI = pp.start();
 266
 267         if (pp.coupled())
 268         {
 269             forAll(pp, i)
 270             {
 271                 transportStencil
 272                 (
 273                     nonEmptyFace,
 274                     faceStencil,
 275                     minOpposedness,
 276                     faceI,
 277                     own[faceI],
 278                     true,                   //stencilHasNeighbour
 279
 280                     oppositeFaces,
 281                     faceStencilSet,
 282                     ownStencil[faceI]
 283                 );
 284                 faceI++;
 285             }
 286         }
 287         else if (!isA<emptyPolyPatch>(pp))
 288         {
 289             forAll(pp, i)
 290             {
 291                 // faceStencil does not contain neighbour
 292                 transportStencil
 293                 (
 294                     nonEmptyFace,
 295                     faceStencil,
 296                     minOpposedness,
 297                     faceI,
 298                     own[faceI],
 299                     false,                  //stencilHasNeighbour
 300
 301                     oppositeFaces,
 302                     faceStencilSet,
 303                     ownStencil[faceI]
 304                 );
 305                 faceI++;
 306             }
 307         }
 308     }
 309
 310
 311     // Swap coupled boundary stencil
 312     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 313
 314     labelListList neiBndStencil(nBnd);
 315     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
 316     {
 317         neiBndStencil[faceI-mesh_.nInternalFaces()] = ownStencil[faceI];
 318     }
 319     //syncTools::swapBoundaryFaceList(mesh_, neiBndStencil);
 320     syncTools::syncBoundaryFaceList
 321     (
 322         mesh_,
 323         neiBndStencil,
 324         eqOp<labelList>(),
 325         dummyTransform()
 326     );
 327
 328
 329
 330     // Do the neighbour side
 331     // ~~~~~~~~~~~~~~~~~~~~~
 332     // - internal faces : get opposite faces on neighbour side
 333     // - boundary faces : empty
 334     // - coupled faces  : in neiBndStencil
 335
 336     neiStencil.setSize(mesh_.nFaces());
 337
 338     // Internal faces
 339     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 340     {
 341         transportStencil
 342         (
 343             nonEmptyFace,
 344             faceStencil,
 345             minOpposedness,
 346             faceI,
 347             nei[faceI],
 348             true,                   //stencilHasNeighbour
 349
 350             oppositeFaces,
 351             faceStencilSet,
 352             neiStencil[faceI]
 353         );
 354     }
 355
 356     // Boundary faces
 357     forAll(patches, patchI)
 358     {
 359         const polyPatch& pp = patches[patchI];
 360         label faceI = pp.start();
 361
 362         if (pp.coupled())
 363         {
 364             forAll(pp, i)
 365             {
 366                 neiStencil[faceI].transfer
 367                 (
 368                     neiBndStencil[faceI-mesh_.nInternalFaces()]
 369                 );
 370                 faceI++;
 371             }
 372         }
 373         else
 374         {
 375             // Boundary has empty neighbour stencil
 376         }
 377     }
 378 }
 379
 380
 381 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 382
 383 Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil
 384 (
 385     const cellToFaceStencil& stencil,
 386     const bool pureUpwind,
 387     const scalar minOpposedness
 388 )
 389 :
 390     extendedCellToFaceStencil(stencil.mesh()),
 391     pureUpwind_(pureUpwind)
 392 {
 393     //forAll(stencil, faceI)
 394     //{
 395     //    const labelList& fCells = stencil[faceI];
 396     //
 397     //    Pout<< "Face:" << faceI << " at:" << mesh_.faceCentres()[faceI]
 398     //        << endl;
 399     //
 400     //    forAll(fCells, i)
 401     //    {
 402     //        label globalI = fCells[i];
 403     //
 404     //        if (globalI < mesh_.nCells())
 405     //        {
 406     //            Pout<< "    cell:" << globalI
 407     //                << " at:" << mesh_.cellCentres()[globalI] << endl;
 408     //        }
 409     //        else
 410     //        {
 411     //            label faceI = globalI-mesh_.nCells() + mesh_.nInternalFaces();
 412     //
 413     //            Pout<< "    boundary:" << faceI
 414     //                << " at:" << mesh_.faceCentres()[faceI] << endl;
 415     //        }
 416     //    }
 417     //}
 418     //Pout<< endl << endl;
 419
 420
 421     // Transport centred stencil to upwind/downwind face
 422     transportStencils
 423     (
 424         stencil,
 425         minOpposedness,
 426         ownStencil_,
 427         neiStencil_
 428     );
 429
 430     {
 431         List<Map<label> > compactMap(Pstream::nProcs());
 432         ownMapPtr_.reset
 433         (
 434             new mapDistribute
 435             (
 436                 stencil.globalNumbering(),
 437                 ownStencil_,
 438                 compactMap
 439             )
 440         );
 441     }
 442
 443     {
 444
 445         List<Map<label> > compactMap(Pstream::nProcs());
 446         neiMapPtr_.reset
 447         (
 448             new mapDistribute
 449             (
 450                 stencil.globalNumbering(),
 451                 neiStencil_,
 452                 compactMap
 453             )
 454         );
 455     }
 456
 457     // stencil now in compact form
 458     if (pureUpwind_)
 459     {
 460         const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh());
 461
 462         List<List<point> > stencilPoints(ownStencil_.size());
 463
 464         // Owner stencil
 465         // ~~~~~~~~~~~~~
 466
 467         collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints);
 468
 469         // Mask off all stencil points on wrong side of face
 470         forAll(stencilPoints, faceI)
 471         {
 472             const point& fc = mesh.faceCentres()[faceI];
 473             const vector& fArea = mesh.faceAreas()[faceI];
 474
 475             const List<point>& points = stencilPoints[faceI];
 476             const labelList& stencil = ownStencil_[faceI];
 477
 478             DynamicList<label> newStencil(stencil.size());
 479             forAll(points, i)
 480             {
 481                 if (((points[i]-fc) & fArea) < 0)
 482                 {
 483                     newStencil.append(stencil[i]);
 484                 }
 485             }
 486             if (newStencil.size() != stencil.size())
 487             {
 488                 ownStencil_[faceI].transfer(newStencil);
 489             }
 490         }
 491
 492
 493         // Neighbour stencil
 494         // ~~~~~~~~~~~~~~~~~
 495
 496         collectData(neiMapPtr_(), neiStencil_, mesh.C(), stencilPoints);
 497
 498         // Mask off all stencil points on wrong side of face
 499         forAll(stencilPoints, faceI)
 500         {
 501             const point& fc = mesh.faceCentres()[faceI];
 502             const vector& fArea = mesh.faceAreas()[faceI];
 503
 504             const List<point>& points = stencilPoints[faceI];
 505             const labelList& stencil = neiStencil_[faceI];
 506
 507             DynamicList<label> newStencil(stencil.size());
 508             forAll(points, i)
 509             {
 510                 if (((points[i]-fc) & fArea) > 0)
 511                 {
 512                     newStencil.append(stencil[i]);
 513                 }
 514             }
 515             if (newStencil.size() != stencil.size())
 516             {
 517                 neiStencil_[faceI].transfer(newStencil);
 518             }
 519         }
 520
 521         // Note: could compact schedule as well. for if cells are not needed
 522         // across any boundary anymore. However relatively rare.
 523     }
 524 }
 525
 526
 527 Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil
 528 (
 529     const cellToFaceStencil& stencil
 530 )
 531 :
 532     extendedCellToFaceStencil(stencil.mesh()),
 533     pureUpwind_(true)
 534 {
 535     // Calculate stencil points with full stencil
 536
 537     ownStencil_ = stencil;
 538
 539     {
 540         List<Map<label> > compactMap(Pstream::nProcs());
 541         ownMapPtr_.reset
 542         (
 543             new mapDistribute
 544             (
 545                 stencil.globalNumbering(),
 546                 ownStencil_,
 547                 compactMap
 548             )
 549         );
 550     }
 551
 552     const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh());
 553
 554     List<List<point> > stencilPoints(ownStencil_.size());
 555     collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints);
 556
 557     // Split stencil into owner and neighbour
 558     neiStencil_.setSize(ownStencil_.size());
 559
 560     forAll(stencilPoints, faceI)
 561     {
 562         const point& fc = mesh.faceCentres()[faceI];
 563         const vector& fArea = mesh.faceAreas()[faceI];
 564
 565         const List<point>& points = stencilPoints[faceI];
 566         const labelList& stencil = ownStencil_[faceI];
 567
 568         DynamicList<label> newOwnStencil(stencil.size());
 569         DynamicList<label> newNeiStencil(stencil.size());
 570         forAll(points, i)
 571         {
 572             if (((points[i]-fc) & fArea) > 0)
 573             {
 574                 newNeiStencil.append(stencil[i]);
 575             }
 576             else
 577             {
 578                 newOwnStencil.append(stencil[i]);
 579             }
 580         }
 581         if (newNeiStencil.size() > 0)
 582         {
 583             ownStencil_[faceI].transfer(newOwnStencil);
 584             neiStencil_[faceI].transfer(newNeiStencil);
 585         }
 586     }
 587
 588     // Should compact schedule. Or have both return the same schedule.
 589     neiMapPtr_.reset(new mapDistribute(ownMapPtr_()));
 590 }
 591
 592
 593 // ************************************************************************* //