applications/utilities/parallelProcessing/decomposePar/domainDecompositionMesh.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 InClass
  25     domainDecomposition
  26
  27 Description
  28     Private member of domainDecomposition.
  29     Decomposes the mesh into bits
  30
  31 \*---------------------------------------------------------------------------*/
  32
  33 #include "domainDecomposition.H"
  34 #include "IOstreams.H"
  35 #include "SLPtrList.H"
  36 #include "boolList.H"
  37 #include "primitiveMesh.H"
  38 #include "cyclicPolyPatch.H"
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 void Foam::domainDecomposition::append(labelList& lst, const label elem)
  43 {
  44     label sz = lst.size();
  45     lst.setSize(sz+1);
  46     lst[sz] = elem;
  47 }
  48
  49
  50 void Foam::domainDecomposition::addInterProcFace
  51 (
  52     const label facei,
  53     const label ownerProc,
  54     const label nbrProc,
  55
  56     List<Map<label> >& nbrToInterPatch,
  57     List<DynamicList<DynamicList<label> > >& interPatchFaces
  58 ) const
  59 {
  60     Map<label>::iterator patchIter = nbrToInterPatch[ownerProc].find(nbrProc);
  61
  62     // Introduce turning index only for internal faces (are duplicated).
  63     label ownerIndex = facei+1;
  64     label nbrIndex = -(facei+1);
  65
  66     if (patchIter != nbrToInterPatch[ownerProc].end())
  67     {
  68         // Existing interproc patch. Add to both sides.
  69         label toNbrProcPatchI = patchIter();
  70         interPatchFaces[ownerProc][toNbrProcPatchI].append(ownerIndex);
  71
  72         if (isInternalFace(facei))
  73         {
  74             label toOwnerProcPatchI = nbrToInterPatch[nbrProc][ownerProc];
  75             interPatchFaces[nbrProc][toOwnerProcPatchI].append(nbrIndex);
  76         }
  77     }
  78     else
  79     {
  80         // Create new interproc patches.
  81         label toNbrProcPatchI = nbrToInterPatch[ownerProc].size();
  82         nbrToInterPatch[ownerProc].insert(nbrProc, toNbrProcPatchI);
  83         DynamicList<label> oneFace;
  84         oneFace.append(ownerIndex);
  85         interPatchFaces[ownerProc].append(oneFace);
  86
  87         if (isInternalFace(facei))
  88         {
  89             label toOwnerProcPatchI = nbrToInterPatch[nbrProc].size();
  90             nbrToInterPatch[nbrProc].insert(ownerProc, toOwnerProcPatchI);
  91             oneFace.clear();
  92             oneFace.append(nbrIndex);
  93             interPatchFaces[nbrProc].append(oneFace);
  94         }
  95     }
  96 }
  97
  98
  99 void Foam::domainDecomposition::decomposeMesh()
 100 {
 101     // Decide which cell goes to which processor
 102     distributeCells();
 103
 104     // Distribute the cells according to the given processor label
 105
 106     // calculate the addressing information for the original mesh
 107     Info<< "\nCalculating original mesh data" << endl;
 108
 109     // set references to the original mesh
 110     const polyBoundaryMesh& patches = boundaryMesh();
 111     const faceList& fcs = faces();
 112     const labelList& owner = faceOwner();
 113     const labelList& neighbour = faceNeighbour();
 114
 115     // loop through the list of processor labels for the cell and add the
 116     // cell shape to the list of cells for the appropriate processor
 117
 118     Info<< "\nDistributing cells to processors" << endl;
 119
 120     // Cells per processor
 121     procCellAddressing_ = invertOneToMany(nProcs_, cellToProc_);
 122
 123     Info<< "\nDistributing faces to processors" << endl;
 124
 125     // Loop through all internal faces and decide which processor they belong to
 126     // First visit all internal faces. If cells at both sides belong to the
 127     // same processor, the face is an internal face. If they are different,
 128     // it belongs to both processors.
 129
 130     procFaceAddressing_.setSize(nProcs_);
 131
 132     // Internal faces
 133     forAll(neighbour, facei)
 134     {
 135         if (cellToProc_[owner[facei]] == cellToProc_[neighbour[facei]])
 136         {
 137             // Face internal to processor. Notice no turning index.
 138             procFaceAddressing_[cellToProc_[owner[facei]]].append(facei+1);
 139         }
 140     }
 141
 142     // for all processors, set the size of start index and patch size
 143     // lists to the number of patches in the mesh
 144     forAll(procPatchSize_, procI)
 145     {
 146         procPatchSize_[procI].setSize(patches.size());
 147         procPatchStartIndex_[procI].setSize(patches.size());
 148     }
 149
 150     forAll(patches, patchi)
 151     {
 152         // Reset size and start index for all processors
 153         forAll(procPatchSize_, procI)
 154         {
 155             procPatchSize_[procI][patchi] = 0;
 156             procPatchStartIndex_[procI][patchi] =
 157                 procFaceAddressing_[procI].size();
 158         }
 159
 160         const label patchStart = patches[patchi].start();
 161
 162         if (!isA<cyclicPolyPatch>(patches[patchi]))
 163         {
 164             // Normal patch. Add faces to processor where the cell
 165             // next to the face lives
 166
 167             const labelUList& patchFaceCells =
 168                 patches[patchi].faceCells();
 169
 170             forAll(patchFaceCells, facei)
 171             {
 172                 const label curProc = cellToProc_[patchFaceCells[facei]];
 173
 174                 // add the face without turning index
 175                 procFaceAddressing_[curProc].append(patchStart+facei+1);
 176
 177                 // increment the number of faces for this patch
 178                 procPatchSize_[curProc][patchi]++;
 179             }
 180         }
 181         else
 182         {
 183             const cyclicPolyPatch& pp = refCast<const cyclicPolyPatch>
 184             (
 185                 patches[patchi]
 186             );
 187             // cyclic: check opposite side on this processor
 188             const labelUList& patchFaceCells = pp.faceCells();
 189
 190             const labelUList& nbrPatchFaceCells =
 191                 pp.neighbPatch().faceCells();
 192
 193             forAll(patchFaceCells, facei)
 194             {
 195                 const label curProc = cellToProc_[patchFaceCells[facei]];
 196                 const label nbrProc = cellToProc_[nbrPatchFaceCells[facei]];
 197                 if (curProc == nbrProc)
 198                 {
 199                     // add the face without turning index
 200                     procFaceAddressing_[curProc].append(patchStart+facei+1);
 201                     // increment the number of faces for this patch
 202                     procPatchSize_[curProc][patchi]++;
 203                 }
 204             }
 205         }
 206     }
 207
 208
 209     // Done internal bits of the new mesh and the ordinary patches.
 210
 211
 212     // Per processor, from neighbour processor to the interprocessorpatch that
 213     // communicates with that neighbour.
 214     List<Map<label> > procNbrToInterPatch(nProcs_);
 215     // Per processor the faces per interprocessorpatch.
 216     List<DynamicList<DynamicList<label> > > interPatchFaces(nProcs_);
 217
 218     // Processor boundaries from internal faces
 219     forAll(neighbour, facei)
 220     {
 221         label ownerProc = cellToProc_[owner[facei]];
 222         label nbrProc = cellToProc_[neighbour[facei]];
 223
 224         if (ownerProc != nbrProc)
 225         {
 226             // inter - processor patch face found.
 227             addInterProcFace
 228             (
 229                 facei,
 230                 ownerProc,
 231                 nbrProc,
 232
 233                 procNbrToInterPatch,
 234                 interPatchFaces
 235             );
 236         }
 237     }
 238
 239     // Add the proper processor faces to the sub information. For faces
 240     // originating from internal faces this is always -1.
 241     List<labelListList> subPatchIDs(nProcs_);
 242     List<labelListList> subPatchStarts(nProcs_);
 243     forAll(interPatchFaces, procI)
 244     {
 245         label nInterfaces = interPatchFaces[procI].size();
 246
 247         subPatchIDs[procI].setSize(nInterfaces, labelList(1, -1));
 248         subPatchStarts[procI].setSize(nInterfaces, labelList(1, 0));
 249     }
 250
 251     // Processor boundaries from split cyclics
 252     forAll(patches, patchi)
 253     {
 254         if (isA<cyclicPolyPatch>(patches[patchi]))
 255         {
 256             const cyclicPolyPatch& pp = refCast<const cyclicPolyPatch>
 257             (
 258                 patches[patchi]
 259             );
 260
 261             // cyclic: check opposite side on this processor
 262             const labelUList& patchFaceCells = pp.faceCells();
 263             const labelUList& nbrPatchFaceCells =
 264                 pp.neighbPatch().faceCells();
 265
 266             // Store old sizes. Used to detect which inter-proc patches
 267             // have been added to.
 268             labelListList oldInterfaceSizes(nProcs_);
 269             forAll(oldInterfaceSizes, procI)
 270             {
 271                 labelList& curOldSizes = oldInterfaceSizes[procI];
 272
 273                 curOldSizes.setSize(interPatchFaces[procI].size());
 274                 forAll(curOldSizes, interI)
 275                 {
 276                     curOldSizes[interI] =
 277                         interPatchFaces[procI][interI].size();
 278                 }
 279             }
 280
 281             // Add faces with different owner and neighbour processors
 282             forAll(patchFaceCells, facei)
 283             {
 284                 const label ownerProc = cellToProc_[patchFaceCells[facei]];
 285                 const label nbrProc = cellToProc_[nbrPatchFaceCells[facei]];
 286                 if (ownerProc != nbrProc)
 287                 {
 288                     // inter - processor patch face found.
 289                     addInterProcFace
 290                     (
 291                         pp.start()+facei,
 292                         ownerProc,
 293                         nbrProc,
 294                         procNbrToInterPatch,
 295                         interPatchFaces
 296                     );
 297                 }
 298             }
 299
 300             // 1. Check if any faces added to existing interfaces
 301             forAll(oldInterfaceSizes, procI)
 302             {
 303                 const labelList& curOldSizes = oldInterfaceSizes[procI];
 304
 305                 forAll(curOldSizes, interI)
 306                 {
 307                     label oldSz = curOldSizes[interI];
 308                     if (interPatchFaces[procI][interI].size() > oldSz)
 309                     {
 310                         // Added faces to this interface. Add an entry
 311                         append(subPatchIDs[procI][interI], patchi);
 312                         append(subPatchStarts[procI][interI], oldSz);
 313                     }
 314                 }
 315             }
 316
 317             // 2. Any new interfaces
 318             forAll(subPatchIDs, procI)
 319             {
 320                 label nIntfcs = interPatchFaces[procI].size();
 321                 subPatchIDs[procI].setSize(nIntfcs, labelList(1, patchi));
 322                 subPatchStarts[procI].setSize(nIntfcs, labelList(1, 0));
 323             }
 324         }
 325     }
 326
 327
 328     // Shrink processor patch face addressing
 329     forAll(interPatchFaces, procI)
 330     {
 331         DynamicList<DynamicList<label> >& curInterPatchFaces =
 332             interPatchFaces[procI];
 333
 334         forAll(curInterPatchFaces, i)
 335         {
 336             curInterPatchFaces[i].shrink();
 337         }
 338         curInterPatchFaces.shrink();
 339     }
 340
 341
 342     // Sort inter-proc patch by neighbour
 343     labelList order;
 344     forAll(procNbrToInterPatch, procI)
 345     {
 346         label nInterfaces = procNbrToInterPatch[procI].size();
 347
 348         procNeighbourProcessors_[procI].setSize(nInterfaces);
 349         procProcessorPatchSize_[procI].setSize(nInterfaces);
 350         procProcessorPatchStartIndex_[procI].setSize(nInterfaces);
 351         procProcessorPatchSubPatchIDs_[procI].setSize(nInterfaces);
 352         procProcessorPatchSubPatchStarts_[procI].setSize(nInterfaces);
 353
 354         //Info<< "Processor " << procI << endl;
 355
 356         // Get sorted neighbour processors
 357         const Map<label>& curNbrToInterPatch = procNbrToInterPatch[procI];
 358         labelList nbrs = curNbrToInterPatch.toc();
 359         sortedOrder(nbrs, order);
 360
 361         DynamicList<DynamicList<label> >& curInterPatchFaces =
 362             interPatchFaces[procI];
 363
 364         forAll(order, i)
 365         {
 366             const label nbrProc = nbrs[i];
 367             const label interPatch = curNbrToInterPatch[nbrProc];
 368
 369             procNeighbourProcessors_[procI][i] =
 370                 nbrProc;
 371             procProcessorPatchSize_[procI][i] =
 372                 curInterPatchFaces[interPatch].size();
 373             procProcessorPatchStartIndex_[procI][i] =
 374                 procFaceAddressing_[procI].size();
 375
 376             // Add size as last element to substarts and transfer
 377             append
 378             (
 379                 subPatchStarts[procI][interPatch],
 380                 curInterPatchFaces[interPatch].size()
 381             );
 382             procProcessorPatchSubPatchIDs_[procI][i].transfer
 383             (
 384                 subPatchIDs[procI][interPatch]
 385             );
 386             procProcessorPatchSubPatchStarts_[procI][i].transfer
 387             (
 388                 subPatchStarts[procI][interPatch]
 389             );
 390
 391             //Info<< "    nbr:" << nbrProc << endl;
 392             //Info<< "    interpatch:" << interPatch << endl;
 393             //Info<< "    size:" << procProcessorPatchSize_[procI][i] << endl;
 394             //Info<< "    start:" << procProcessorPatchStartIndex_[procI][i]
 395             //    << endl;
 396             //Info<< "    subPatches:"
 397             //    << procProcessorPatchSubPatchIDs_[procI][i]
 398             //    << endl;
 399             //Info<< "    subStarts:"
 400             //    << procProcessorPatchSubPatchStarts_[procI][i] << endl;
 401
 402             // And add all the face labels for interPatch
 403             DynamicList<label>& interPatchFaces =
 404                 curInterPatchFaces[interPatch];
 405
 406             forAll(interPatchFaces, j)
 407             {
 408                 procFaceAddressing_[procI].append(interPatchFaces[j]);
 409             }
 410             interPatchFaces.clearStorage();
 411         }
 412         curInterPatchFaces.clearStorage();
 413         procFaceAddressing_[procI].shrink();
 414     }
 415
 416
 417 ////XXXXXXX
 418 //// Print a bit
 419 //    forAll(procPatchStartIndex_, procI)
 420 //    {
 421 //        Info<< "Processor:" << procI << endl;
 422 //
 423 //        Info<< "    total faces:" << procFaceAddressing_[procI].size()
 424 //            << endl;
 425 //
 426 //        const labelList& curProcPatchStartIndex = procPatchStartIndex_[procI];
 427 //
 428 //        forAll(curProcPatchStartIndex, patchI)
 429 //        {
 430 //            Info<< "    patch:" << patchI
 431 //                << "\tstart:" << curProcPatchStartIndex[patchI]
 432 //                << "\tsize:" << procPatchSize_[procI][patchI]
 433 //                << endl;
 434 //        }
 435 //    }
 436 //    Info<< endl;
 437 //
 438 //    forAll(procNeighbourProcessors_, procI)
 439 //    {
 440 //        Info<< "Processor " << procI << endl;
 441 //
 442 //        forAll(procNeighbourProcessors_[procI], i)
 443 //        {
 444 //            Info<< "    nbr:" << procNeighbourProcessors_[procI][i] << endl;
 445 //            Info<< "    size:" << procProcessorPatchSize_[procI][i] << endl;
 446 //            Info<< "    start:" << procProcessorPatchStartIndex_[procI][i]
 447 //                << endl;
 448 //        }
 449 //    }
 450 //    Info<< endl;
 451 //
 452 //    forAll(procFaceAddressing_, procI)
 453 //    {
 454 //        Info<< "Processor:" << procI << endl;
 455 //
 456 //        Info<< "    faces:" << procFaceAddressing_[procI] << endl;
 457 //    }
 458
 459
 460
 461     Info<< "\nDistributing points to processors" << endl;
 462     // For every processor, loop through the list of faces for the processor.
 463     // For every face, loop through the list of points and mark the point as
 464     // used for the processor. Collect the list of used points for the
 465     // processor.
 466
 467     forAll(procPointAddressing_, procI)
 468     {
 469         boolList pointLabels(nPoints(), false);
 470
 471         // Get reference to list of used faces
 472         const labelList& procFaceLabels = procFaceAddressing_[procI];
 473
 474         forAll(procFaceLabels, facei)
 475         {
 476             // Because of the turning index, some labels may be negative
 477             const labelList& facePoints = fcs[mag(procFaceLabels[facei]) - 1];
 478
 479             forAll(facePoints, pointi)
 480             {
 481                 // Mark the point as used
 482                 pointLabels[facePoints[pointi]] = true;
 483             }
 484         }
 485
 486         // Collect the used points
 487         labelList& procPointLabels = procPointAddressing_[procI];
 488
 489         procPointLabels.setSize(pointLabels.size());
 490
 491         label nUsedPoints = 0;
 492
 493         forAll(pointLabels, pointi)
 494         {
 495             if (pointLabels[pointi])
 496             {
 497                 procPointLabels[nUsedPoints] = pointi;
 498
 499                 nUsedPoints++;
 500             }
 501         }
 502
 503         // Reset the size of used points
 504         procPointLabels.setSize(nUsedPoints);
 505     }
 506 }
 507
 508 // ************************************************************************* //