src/decompositionMethods/scotchDecomp/engineScotchDecomp/engineScotchDecomp.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "engineScotchDecomp.H"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 namespace Foam
  32 {
  33     defineTypeNameAndDebug(engineScotchDecomp, 0);
  34
  35     addToRunTimeSelectionTable
  36     (
  37         scotchDecomp,
  38         engineScotchDecomp,
  39         dictionaryMesh
  40     );
  41 }
  42
  43 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  44
  45 Foam::engineScotchDecomp::engineScotchDecomp
  46 (
  47     const dictionary& decompositionDict,
  48     const polyMesh& mesh
  49 )
  50 :
  51     scotchDecomp(decompositionDict, mesh),
  52     dict_(decompositionDict.subDict(typeName + "Coeffs")),
  53     slidingPatchPairs_(dict_.lookup("slidingPatchPairs")),
  54     expandSliding_(dict_.lookup("expandSliding"))
  55 {}
  56
  57
  58 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  59
  60 Foam::labelList Foam::engineScotchDecomp::decompose
  61 (
  62     const pointField& points,
  63     const scalarField& pointWeights
  64 )
  65 {
  66     if (points.size() != mesh().nCells())
  67     {
  68         FatalErrorIn
  69         (
  70             "engineScotchDecomp::decompose\n"
  71             "(\n"
  72             "    const pointField&,\n"
  73             "    const scalarField&\n"
  74             ")"
  75         )   << "Can use this decomposition method only for the whole mesh"
  76             << endl
  77             << "and supply one coordinate (cellCentre) for every cell." << endl
  78             << "The number of coordinates " << points.size() << endl
  79             << "The number of cells in the mesh " << mesh().nCells()
  80             << exit(FatalError);
  81     }
  82
  83     // Create clustering to coarse level
  84
  85     labelList fineToCoarse(mesh().nCells(), -1);
  86
  87     // Mask all cells in the cylinder layering region with true
  88     // Used in two-pass decomposition later
  89     boolList cylinderMask(mesh().nCells(), false);
  90     label nClusters = 0;
  91
  92     // Locate piston patch and cluster up colums, using opposite faces
  93     label pistonPatchID = mesh().boundaryMesh().findPatchID("piston");
  94
  95     // Go through the sliding pairs and mark clustering
  96     forAll (slidingPatchPairs_, pairI)
  97     {
  98         // Locate patch and cluster up colums, using opposite faces
  99         label firstPatchID = mesh().boundaryMesh().findPatchID
 100         (
 101             slidingPatchPairs_[pairI].first()
 102         );
 103
 104         label secondPatchID = mesh().boundaryMesh().findPatchID
 105         (
 106             slidingPatchPairs_[pairI].second()
 107         );
 108
 109         if (firstPatchID == -1 || secondPatchID == -1)
 110         {
 111             FatalErrorIn
 112             (
 113                 "labelList engineScotchDecomp::decompose\n"
 114                 "(\n"
 115                 "    const pointField& points,\n"
 116                 "    const scalarField& pointWeights\n"
 117                 ")"
 118             )   << "Cannot find sliding patch pair "
 119                 << slidingPatchPairs_[pairI]
 120                 << abort(FatalError);
 121         }
 122
 123         // Put all cells next to the patch into a cluster
 124         if (expandSliding_)
 125         {
 126             // Use point-cell addressing from patch points
 127             const labelListList& pc = mesh().pointCells();
 128
 129             // First side
 130             const labelList& mp1 =
 131                 mesh().boundaryMesh()[firstPatchID].meshPoints();
 132
 133             forAll (mp1, pointI)
 134             {
 135                 const labelList& curCells = pc[mp1[pointI]];
 136
 137                 forAll (curCells, cellI)
 138                 {
 139                     fineToCoarse[curCells[cellI]] = nClusters;
 140                     cylinderMask[curCells[cellI]] = true;
 141                 }
 142             }
 143
 144             // Second side
 145             {
 146                 const labelList& mp2 =
 147                     mesh().boundaryMesh()[secondPatchID].meshPoints();
 148
 149                 forAll (mp2, pointI)
 150                 {
 151                     const labelList& curCells = pc[mp2[pointI]];
 152
 153                     forAll (curCells, cellI)
 154                     {
 155                         fineToCoarse[curCells[cellI]] = nClusters;
 156                         cylinderMask[curCells[cellI]] = true;
 157                     }
 158                 }
 159             }
 160         }
 161         else
 162         {
 163             // First side
 164             {
 165                 const labelList& fc1 =
 166                     mesh().boundaryMesh()[firstPatchID].faceCells();
 167
 168                 forAll (fc1, fcI)
 169                 {
 170                     fineToCoarse[fc1[fcI]] = nClusters;
 171                     cylinderMask[fc1[fcI]] = true;
 172                 }
 173             }
 174
 175             // Second side
 176             {
 177                 const labelList& fc2 =
 178                     mesh().boundaryMesh()[secondPatchID].faceCells();
 179
 180                 forAll (fc2, fcI)
 181                 {
 182                     fineToCoarse[fc2[fcI]] = nClusters;
 183                     cylinderMask[fc2[fcI]] = true;
 184                 }
 185             }
 186         }
 187
 188         nClusters++;
 189     }
 190
 191     if (pistonPatchID > -1)
 192     {
 193         // Found piston patch
 194         const faceList& f = mesh().allFaces();
 195         const cellList& c = mesh().cells();
 196         const labelList& owner = mesh().faceOwner();
 197         const labelList& neighbour = mesh().faceNeighbour();
 198
 199         const labelList& faceCells =
 200             mesh().boundaryMesh()[pistonPatchID].faceCells();
 201
 202         forAll (faceCells, faceI)
 203         {
 204             // Get face index
 205             label curFaceNo = faceI
 206                 + mesh().boundaryMesh()[pistonPatchID].start();
 207
 208             // Get cell index
 209             label curCellNo = faceCells[faceI];
 210
 211             // Mark cell to cluster
 212             if (fineToCoarse[curCellNo] < 0)
 213             {
 214                 // New cluster
 215                 fineToCoarse[curCellNo] = nClusters;
 216                 cylinderMask[curCellNo] = true;
 217
 218                 for(;;)
 219                 {
 220                     // Attempt to find the next face and cell
 221                     curFaceNo = c[curCellNo].opposingFaceLabel(curFaceNo, f);
 222
 223                     if (curFaceNo > -1)
 224                     {
 225                         // Face found, try for a cell
 226                         if (curFaceNo < mesh().nInternalFaces())
 227                         {
 228                             if (owner[curFaceNo] == curCellNo)
 229                             {
 230                                 curCellNo = neighbour[curFaceNo];
 231                             }
 232                             else if (neighbour[curFaceNo] == curCellNo)
 233                             {
 234                                 curCellNo = owner[curFaceNo];
 235                             }
 236                             else
 237                             {
 238                                 // Error in layering.  Should never happen
 239                                 break;
 240                             }
 241
 242                             // Mark cell to cluster
 243                             fineToCoarse[curCellNo] = nClusters;
 244                             cylinderMask[curCellNo] = true;
 245                         }
 246                         else
 247                         {
 248                             // Hit boundary face
 249                             break;
 250                         }
 251                     }
 252                     else
 253                     {
 254                         // Cannot find opposing face: out of prismatic region
 255                         break;
 256                     }
 257                 }
 258
 259                 // Next cluster
 260                 nClusters++;
 261             }
 262         }
 263     }
 264
 265     // Count cylinder cells from mask
 266     label nCylinderCells = 0;
 267
 268     forAll (cylinderMask, cellI)
 269     {
 270         if (cylinderMask[cellI]) nCylinderCells++;
 271     }
 272
 273     label nStaticCells = mesh().nCells() - nCylinderCells;
 274     label nCylClusters = nClusters;
 275
 276     // Visit all unmarked cells and put them into single clusters
 277     forAll (fineToCoarse, cellI)
 278     {
 279         if (fineToCoarse[cellI] == -1)
 280         {
 281             fineToCoarse[cellI] = nClusters;
 282             nClusters++;
 283         }
 284     }
 285
 286     label nStaticClusters = nClusters - nCylClusters;
 287
 288     Info<< "Number of cells: " << mesh().nCells()
 289         << ", in cylinder + sliding: " << nCylinderCells
 290         << ", in static part: " << nStaticCells << nl
 291         << "Number of cylinder clusters " << nCylClusters
 292         << ", static clusters " << nStaticClusters
 293         << ", total clusters " << nClusters << endl;
 294
 295     // Mark-up complete.  Create point centres and weights for all clusters
 296     vectorField clusterCentres(nClusters, vector::zero);
 297
 298     // Stabilise cluster volumes in case a cluster ends up empty
 299     // It is possible to have empty clusters without connectivity
 300     scalarField clusterVols(nClusters, SMALL);
 301     scalarField clusterWeights(nClusters, 0);
 302
 303     const vectorField& centres = mesh().cellCentres();
 304     const scalarField& vols = mesh().cellVolumes();
 305
 306     forAll (fineToCoarse, cellI)
 307     {
 308         const label& curCoarse = fineToCoarse[cellI];
 309
 310         clusterCentres[curCoarse] += centres[cellI]*vols[cellI];
 311         clusterVols[curCoarse] += vols[cellI];
 312         clusterWeights[curCoarse] += 1;
 313     }
 314
 315     clusterCentres /= clusterVols;
 316
 317     // Execute decomposition, first on cylinder layering zone, then on the rest
 318
 319     // Collect cell-cells in cylinder layering zone and the rest
 320
 321     // Create and cellCells hash lookup on two pieces
 322
 323     // Note: cell clusters come first and they will be done
 324     // on a shortened list.  Static clusters need to be renumbered by
 325     // throwing away the first part of the list
 326
 327     List<labelHashSet> cylCellCellsHash(nCylClusters);
 328     List<labelHashSet> staticCellCellsHash(nStaticClusters);
 329
 330     const labelListList cellCells = mesh().cellCells();
 331
 332     forAll (cellCells, cellI)
 333     {
 334         const labelList& curCC = cellCells[cellI];
 335
 336         label curCluster = fineToCoarse[cellI];
 337
 338         if (cylinderMask[cellI])
 339         {
 340             labelHashSet& curCylAddr = cylCellCellsHash[curCluster];
 341
 342             // Collect neighbour cluster addressing
 343             forAll (curCC, neiI)
 344             {
 345                 // Add neighbour if marked
 346                 if (cylinderMask[curCC[neiI]])
 347                 {
 348                     label nbrCluster = fineToCoarse[curCC[neiI]];
 349
 350                     if (nbrCluster != curCluster)
 351                     {
 352                         if (!curCylAddr.found(nbrCluster))
 353                         {
 354                             curCylAddr.insert(nbrCluster);
 355                         }
 356                     }
 357                 }
 358             }
 359         }
 360         else
 361         {
 362             // Offset index
 363             curCluster -= nCylClusters;
 364
 365             labelHashSet& curStaticAddr = staticCellCellsHash[curCluster];
 366
 367             forAll (curCC, neiI)
 368             {
 369                 // Add neighbour if marked
 370                 if (!cylinderMask[curCC[neiI]])
 371                 {
 372                     label nbrCluster = fineToCoarse[curCC[neiI]]
 373                         - nCylClusters;
 374
 375                     if (nbrCluster != curCluster)
 376                     {
 377                         if (!curStaticAddr.found(nbrCluster))
 378                         {
 379                             curStaticAddr.insert(nbrCluster);
 380                         }
 381                     }
 382                 }
 383             }
 384         }
 385     }
 386
 387     // Pack cellCells on the cylinder
 388     labelListList cylCellCells(nCylClusters);
 389
 390     forAll (cylCellCellsHash, clusterI)
 391     {
 392         cylCellCells[clusterI] = cylCellCellsHash[clusterI].toc();
 393     }
 394
 395     // Decompose cylinder: size of list equals the number of clusters
 396     // in the cylinder region
 397     vectorField clusterCentresCyl = clusterCentres;
 398     scalarField clusterWeightsCyl = clusterWeights;
 399
 400     clusterCentresCyl.setSize(nCylClusters);
 401     clusterWeightsCyl.setSize(nCylClusters);
 402
 403     labelList cylDecomp = scotchDecomp::decompose
 404     (
 405         cylCellCells,
 406         clusterCentresCyl,
 407         clusterWeightsCyl
 408     );
 409
 410     // Decompose static: size of list equals the number of clusters
 411
 412     labelListList staticCellCells(nStaticClusters);
 413
 414     forAll (staticCellCellsHash, clusterI)
 415     {
 416         staticCellCells[clusterI] = staticCellCellsHash[clusterI].toc();
 417     }
 418
 419     vectorField clusterCentresStatic(nStaticClusters);
 420     scalarField clusterWeightsStatic(nStaticClusters);
 421
 422     forAll (clusterCentresStatic, i)
 423     {
 424         clusterCentresStatic[i] = clusterCentres[nCylClusters + i];
 425         clusterWeightsStatic[i] = clusterWeights[nCylClusters + i];
 426     }
 427
 428     labelList staticDecomp = scotchDecomp::decompose
 429     (
 430         staticCellCells,
 431         clusterCentresStatic,
 432         clusterWeightsStatic
 433     );
 434
 435     // Reconstruct final decomposition
 436
 437     labelList finalDecomp(mesh().nCells(), -1);
 438
 439     forAll (cylinderMask, cellI)
 440     {
 441         if (cylinderMask[cellI])
 442         {
 443             // Cylinder cell
 444             finalDecomp[cellI] = cylDecomp[fineToCoarse[cellI]];
 445         }
 446         else
 447         {
 448             // Static cell
 449             finalDecomp[cellI] =
 450                 staticDecomp[fineToCoarse[cellI] - nCylClusters];
 451         }
 452     }
 453
 454     return finalDecomp;
 455 }
 456
 457
 458 // ************************************************************************* //