src/tetDecompositionFiniteElement/tetPolyMeshCellDecomp/addParallelPointPatchCellDecomp.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "tetPolyMeshCellDecomp.H"
  30 #include "Time.H"
  31 #include "PstreamCombineReduceOps.H"
  32 #include "processorTetPolyPatchCellDecomp.H"
  33 #include "globalTetPolyPatchCellDecomp.H"
  34 #include "Pstream.H"
  35 #include "labelIOList.H"
  36 #include "globalMeshData.H"
  37
  38 using namespace Foam;
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 // Thing label combination class
  43 template<class Type>
  44 class globalPatchCombine
  45 {
  46 public:
  47
  48     void operator()
  49     (
  50         SLList<Type>& globalObjects,
  51         const SLList<Type>& myObjects
  52     ) const
  53     {
  54         // For each of my points check whether it exists in the global
  55         // points list; if not, add it to the global points
  56
  57         for
  58         (
  59             typename SLList<Type>::const_iterator myObjectsIter =
  60                 myObjects.begin();
  61             myObjectsIter != myObjects.end();
  62             ++myObjectsIter
  63         )
  64         {
  65             const Type& curMyType = myObjectsIter();
  66
  67             bool found = false;
  68
  69             for
  70             (
  71                 typename SLList<Type>::iterator globalObjectsIter =
  72                     globalObjects.begin();
  73                 globalObjectsIter != globalObjects.end();
  74                 ++globalObjectsIter
  75             )
  76             {
  77                 if (globalObjectsIter() == curMyType)
  78                 {
  79                     found = true;
  80                     break;
  81                 }
  82             }
  83
  84             if (!found)
  85             {
  86                 globalObjects.append(curMyType);
  87             }
  88         }
  89     }
  90 };
  91
  92
  93 void Foam::tetPolyMeshCellDecomp::addParallelPointPatch()
  94 {
  95     // Note:
  96     // The processor point patch will be added if processor boundaries
  97     // exist in the case.  If the mesh with processor boundaries is
  98     // not created during a parallel run (e.g. decomposePar), the
  99     // addressing will be dummy.  HJ, 19/Mar/2002
 100
 101     if (mesh_.globalData().parallel())
 102     {
 103         // This cannot be done if I do not know global point/face/cell mapping;
 104         // at the very least I need it on all parallel patches.  Having in mind
 105         // the added expense of reading/writing of this data, the global
 106         // mapping files will still be used.  I can do the globally
 107         // shared points because they carry their own global labels,
 108         // but not for cut edges shared by several processors.
 109         // HJ, 21/Apr/2003
 110
 111         if (debug)
 112         {
 113             Info<< "void tetPolyMeshCellDecomp::addParallelPointPatch() : "
 114                 << "Reading pointProcAddressing." << endl;
 115         }
 116
 117         labelIOList pointProcAddressing
 118         (
 119             IOobject
 120             (
 121                 "pointProcAddressing",
 122                 mesh_.time().findInstance(operator()().meshDir(), "faces"),
 123                 mesh_.meshSubDir,
 124                 mesh_,
 125                 IOobject::MUST_READ,
 126                 IOobject::NO_WRITE
 127             )
 128         );
 129
 130         if (debug)
 131         {
 132             Info<< "void tetPolyMeshCellDecomp::addParallelPointPatch() : "
 133                 << "Reading cellProcAddressing." << endl;
 134         }
 135
 136         labelIOList cellProcAddressing
 137         (
 138             IOobject
 139             (
 140                 "cellProcAddressing",
 141                 mesh_.time().findInstance(operator()().meshDir(), "faces"),
 142                 mesh_.meshSubDir,
 143                 mesh_,
 144                 IOobject::MUST_READ,
 145                 IOobject::NO_WRITE
 146             )
 147         );
 148
 149         const label globalCellOffset = mesh_.globalData().nTotalPoints();
 150
 151         // Add the cell centres to the lookup list
 152         label oldSize = pointProcAddressing.size();
 153         pointProcAddressing.setSize(oldSize + cellProcAddressing.size());
 154
 155         forAll (cellProcAddressing, cellI)
 156         {
 157             pointProcAddressing[oldSize + cellI] =
 158                 cellProcAddressing[cellI] + globalCellOffset;
 159         }
 160
 161         // Get the list of local parallel processor points, their addressing
 162         // and global labels
 163         const labelList& localParPoints =
 164           mesh_.globalData().sharedPointLabels();
 165
 166         // Do parallel edges
 167
 168         // Get the list of local parallel processor edges
 169         const edgeList& localParEdges = parallelEdges();
 170
 171         // Extract global numbers for each of the parallel edges
 172         SLList<edge> globalParEdges;
 173
 174         forAll (localParEdges, edgeI)
 175         {
 176             globalParEdges.append
 177             (
 178                 edge
 179                 (
 180                     pointProcAddressing[localParEdges[edgeI].start()],
 181                     pointProcAddressing[localParEdges[edgeI].end()]
 182                 )
 183             );
 184         }
 185
 186         // Create the local-to-global edge mapping
 187         labelList localEdgeMapping(localParEdges.size(), -1);
 188
 189         if (Pstream::parRun())
 190         {
 191             // Create a global list of edges by reduction from all processors
 192             combineReduce(globalParEdges, globalPatchCombine<edge>());
 193
 194             // Find out which of the parallel edges are local.  For
 195             // easier search indexing make a plain list out ot the
 196             // singly linked list of global edges
 197             edgeList ge(globalParEdges);
 198
 199             forAll (localParEdges, edgeI)
 200             {
 201                 edge curGlobal =
 202                     edge
 203                     (
 204                         pointProcAddressing[localParEdges[edgeI].start()],
 205                         pointProcAddressing[localParEdges[edgeI].end()]
 206                     );
 207
 208                 forAll (ge, geI)
 209                 {
 210                     if (curGlobal == ge[geI])
 211                     {
 212                         localEdgeMapping[edgeI] = geI;
 213                         break;
 214                     }
 215                 }
 216             }
 217         }
 218
 219         // Debug check
 220         if (debug)
 221         {
 222             if (localEdgeMapping.size() > 0)
 223             {
 224                 if (min(localEdgeMapping) < 0)
 225                 {
 226                     FatalErrorIn
 227                     (
 228                         "void Foam::tetPolyMeshCellDecomp::"
 229                         "addParallelPointPatch()"
 230                     )   << "Error in parallel points patch: edges"
 231                         << abort(FatalError);
 232                 }
 233             }
 234         }
 235
 236         // Do parallel cut edges
 237
 238         // Algorithm
 239         // Parallel cut edges come in two flavours: the ones local to
 240         // the processor and the ones shared between several
 241         // processors.  The first lot cause no trouble, but for the
 242         // second we need to make sure that only one processor
 243         // multiplies out the global sum part of the matrix.
 244         // At the same time, the local part of the product needs to be
 245         // done on every processor seprately.  This is catered for by
 246         // two mechanisms:
 247         // 1) Local edges will be calculated here for every patch.
 248         //    The list needs to be complete.
 249         // 2) As the lists are calculated, a global list of cut edges
 250         // is assembled.  The first time the edge is added into the
 251         // global list, it is accepted for global sum multiplication.
 252         // If the same global edge is found, the contribution to the
 253         // global sum is blocked.
 254
 255         // Count the maximum number of edges coming from the patch
 256         label maxEdgesOnPatch = 0;
 257
 258         // For every local point get a list of edges
 259         forAll (localParPoints, pointI)
 260         {
 261             maxEdgesOnPatch += nEdgesForPoint(localParPoints[pointI]);
 262         }
 263
 264         edgeList localCutEdges(maxEdgesOnPatch);
 265         label nCutEdges = 0;
 266
 267         // Go through all the local points and get all the edges coming
 268         // from that point.  Check if the edge is local, if not, it is cut
 269
 270         const unallocLabelList& L = lduAddr().lowerAddr();
 271         const unallocLabelList& U = lduAddr().upperAddr();
 272
 273         forAll (localParPoints, pointI)
 274         {
 275             const label curPoint = localParPoints[pointI];
 276
 277             labelList curEdges = edgesForPoint(curPoint);
 278
 279             forAll (curEdges, edgeI)
 280             {
 281                 edge e(L[curEdges[edgeI]], U[curEdges[edgeI]]);
 282                 label otherEnd = -1;
 283
 284                 if (e.start() == curPoint)
 285                 {
 286                     otherEnd = e.end();
 287                 }
 288                 else
 289                 {
 290                     otherEnd = e.start();
 291                 }
 292
 293                 // If the other end is not a local point, this is a cut edge
 294                 bool found = false;
 295
 296                 forAll (localParPoints, compI)
 297                 {
 298                     if (localParPoints[compI] == otherEnd)
 299                     {
 300                         found = true;
 301                         break;
 302                     }
 303                 }
 304
 305                 if (!found)
 306                 {
 307                     // This is a cut edge. Add it to the list
 308                     localCutEdges[nCutEdges] = e;
 309                     nCutEdges++;
 310                 }
 311             }
 312         }
 313
 314         // Reset the size of the local cut edge list
 315         localCutEdges.setSize(nCutEdges);
 316
 317         // Create a masking array.  For the edges that do not contribute
 318         // to the global product, the mask will be set to zero
 319         scalarField localCutEdgeMask(localCutEdges.size(), 1);
 320
 321         if (Pstream::parRun())
 322         {
 323             // Creating the mask
 324             // Every processor goes through all of its edges and creates the
 325             // global version of the edge. If such an edge does not
 326             // already exist, it is added to the list, if it does, its
 327             // mask is set to zero (somebody is already doing the
 328             // multiplication)
 329             if (Pstream::master())
 330             {
 331                 // Master's mask is always one. Add all edges to the list
 332                 SLList<edge> globalCutEdges;
 333
 334                 forAll (localCutEdges, edgeI)
 335                 {
 336                     globalCutEdges.append
 337                     (
 338                         edge
 339                         (
 340                             pointProcAddressing[localCutEdges[edgeI].start()],
 341                             pointProcAddressing[localCutEdges[edgeI].end()]
 342                         )
 343                     );
 344                 }
 345
 346                 // Send the list to the first slave
 347                 OPstream toFirstSlave(Pstream::blocking, Pstream::firstSlave());
 348                 toFirstSlave << globalCutEdges;
 349             }
 350             else
 351             {
 352                 // Slave.  Read the list from the previous slave, do
 353                 // local processing and send it to the next slave
 354                 for
 355                 (
 356                     int slave = Pstream::firstSlave();
 357                     slave <= Pstream::lastSlave();
 358                     slave++
 359                 )
 360                 {
 361                     if (Pstream::myProcNo() == slave)
 362                     {
 363                         // My turn.
 364                         int receiveFrom = slave - 1;
 365                         int sendTo = slave + 1;
 366
 367                         if (Pstream::myProcNo() == Pstream::firstSlave())
 368                         {
 369                             // Receive from master
 370                             receiveFrom = Pstream::masterNo();
 371                         }
 372
 373                         IPstream rf(Pstream::blocking, receiveFrom);
 374                         SLList<edge> globalCutEdges(rf);
 375
 376                         // Check local cut edges against the list
 377                         forAll (localCutEdges, edgeI)
 378                         {
 379                             edge e
 380                             (
 381                                 pointProcAddressing
 382                                     [localCutEdges[edgeI].start()],
 383                                 pointProcAddressing[localCutEdges[edgeI].end()]
 384                             );
 385
 386                             bool found = false;
 387
 388                             for
 389                             (
 390                                 SLList<edge>::iterator gEdgeIter =
 391                                     globalCutEdges.begin();
 392                                 gEdgeIter != globalCutEdges.end();
 393                                 ++gEdgeIter
 394                             )
 395                             {
 396                                 if (gEdgeIter() == e)
 397                                 {
 398                                     found = true;
 399                                     break;
 400                                 }
 401                             }
 402
 403                             if (found)
 404                             {
 405                                 // Edge already exists.  Set mask to zero
 406                                 localCutEdgeMask[edgeI] = 0;
 407                             }
 408                             else
 409                             {
 410                                 // Edge not found. Add it to the list
 411                                 globalCutEdges.append(e);
 412                             }
 413                         }
 414
 415                         // Re-transmit the list to the next processor
 416                         if (slave < Pstream::lastSlave())
 417                         {
 418                             OPstream passOnEdges(Pstream::blocking, sendTo);
 419                             passOnEdges << globalCutEdges;
 420                         }
 421                     }
 422                 }
 423             }
 424         }
 425
 426         if (debug)
 427         {
 428             Info<< "void tetPolyMeshCellDecomp::addParallelPointPatch() : "
 429                 << "Adding globalTetPolyPatchCellDecomp."
 430                 << "  nGlobalPoints: "
 431                 << mesh_.globalData().nGlobalPoints() << endl;
 432         }
 433
 434         // Add the processor point patch
 435         boundary_.setSize(boundary_.size() + 1);
 436
 437         boundary_.set
 438         (
 439             boundary_.size() - 1,
 440             new globalTetPolyPatchCellDecomp
 441             (
 442                 mesh_.globalData().nGlobalPoints(),
 443                 localParPoints,
 444                 mesh_.globalData().sharedPointAddr(),
 445                 globalParEdges.size(),
 446                 localParEdges,
 447                 localEdgeMapping,
 448                 localCutEdges,
 449                 localCutEdgeMask,
 450                 boundary_,
 451                 boundary_.size() - 1
 452             )
 453         );
 454     }
 455 }
 456
 457
 458 // ************************************************************************* //