src/tetDecompositionFiniteElement/tetPolyMeshCellDecomp/tetPolyPatches/constraint/processor/calcProcessorTetPolyPatchCellDecompAddr.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     Calculate cut edge processor addressing, needed for vector-matrix
  27     multiply on processor boundaries.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "processorTetPolyPatchCellDecomp.H"
  32 #include "tetPolyBoundaryMeshCellDecomp.H"
  33 #include "tetPolyMeshCellDecomp.H"
  34 #include "boolList.H"
  35 #include "globalTetPolyPatchCellDecomp.H"
  36 #include "primitiveFacePatch.H"
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42
  43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  44
  45 void processorTetPolyPatchCellDecomp::calcCutEdgeIndices() const
  46 {
  47     if (debug)
  48     {
  49         Info<< "void processorTetPolyPatchCellDecomp::"
  50             << "calcCutEdgeIndices() const : "
  51             << "calculating cut edge indices"
  52             << endl;
  53     }
  54
  55     if (cutEdgeIndicesPtr_)
  56     {
  57         FatalErrorIn
  58         (
  59             "void processorTetPolyPatchCellDecomp::calcCutEdgesIndices() const"
  60         )   << "addressing already allocated"
  61             << abort(FatalError);
  62     }
  63
  64     // Make a list over all edges in the mesh.  Mark the ones that are local
  65     // to the patch and then collect the rest
  66
  67     boolList isLocal(boundaryMesh().mesh().nEdges(), false);
  68
  69     // Get reference to local edge indices
  70     const labelList& localEdges = localEdgeIndices();
  71
  72     forAll (localEdges, edgeI)
  73     {
  74         isLocal[localEdges[edgeI]] = true;
  75     }
  76
  77     // Get reference to parallel shared edges
  78     const labelList& sharedEdges =
  79         boundaryMesh().globalPatch().localEdgeIndices();
  80
  81     forAll (sharedEdges, sharedI)
  82     {
  83         isLocal[sharedEdges[sharedI]] = true;
  84     }
  85
  86     const labelList& sharedCutEdges =
  87         boundaryMesh().globalPatch().cutEdgeIndices();
  88
  89     forAll (sharedCutEdges, sharedCutI)
  90     {
  91         isLocal[sharedCutEdges[sharedCutI]] = true;
  92     }
  93
  94     // Count the maximum number of edges coming from the patch
  95     label maxEdgesOnPatch = 0;
  96
  97     const tetPolyMeshCellDecomp& mesh = boundaryMesh().mesh();
  98
  99     const labelList& mp = meshPoints();
 100
 101     forAll (mp, pointI)
 102     {
 103         maxEdgesOnPatch += mesh.nEdgesForPoint(mp[pointI]);
 104     }
 105
 106     // Allocate the array
 107     cutEdgeIndicesPtr_ = new labelList(maxEdgesOnPatch, -1);
 108     labelList& cutEdgeInd = *cutEdgeIndicesPtr_;
 109
 110     label nCutEdgeInd = 0;
 111
 112     // Go through all the local points and get all the edges coming
 113     // from that point.  Check if the edge has been marked as local;
 114     // if not, add it to the list of cut edges.
 115     forAll (mp, pointI)
 116     {
 117         labelList curEdges = mesh.edgesForPoint(mp[pointI]);
 118
 119         forAll (curEdges, edgeI)
 120         {
 121             if (!isLocal[curEdges[edgeI]])
 122             {
 123                 cutEdgeInd[nCutEdgeInd] = curEdges[edgeI];
 124                 nCutEdgeInd++;
 125             }
 126         }
 127     }
 128
 129     // Reset the size of the edge list
 130     cutEdgeInd.setSize(nCutEdgeInd);
 131
 132     if (debug)
 133     {
 134         Info<< "void processorTetPolyPatchCellDecomp::"
 135             << "calcCutEdgeIndices() const : "
 136             << "finished calculating cut edge indices"
 137             << endl;
 138     }
 139 }
 140
 141
 142 void processorTetPolyPatchCellDecomp::calcCutEdgeAddressing() const
 143 {
 144     if
 145     (
 146         cutEdgeOwnerIndicesPtr_
 147      || cutEdgeOwnerStartPtr_
 148      || cutEdgeNeighbourIndicesPtr_
 149      || cutEdgeNeighbourStartPtr_
 150      || doubleCutEdgeIndicesPtr_
 151      || doubleCutOwnerPtr_
 152      || doubleCutNeighbourPtr_
 153     )
 154     {
 155         FatalErrorIn
 156         (
 157             "void processorTetPolyPatchCellDecomp::calcCutEdgeAddressing() "
 158             "const"
 159         )   << "addressing already allocated"
 160             << abort(FatalError);
 161     }
 162
 163     // Make a list over all edges in the mesh.  Mark the ones that are local
 164     // to the patch and then collect the rest.
 165     // For doubly cut edges, mark up the local points
 166
 167     const tetPolyMeshCellDecomp& mesh = boundaryMesh().mesh();
 168
 169     // Get reference to local edge indices
 170     const labelList& localEdges = localEdgeIndices();
 171
 172     // Get reference to mesh point addressing for the patch
 173     const labelList& mp = meshPoints();
 174
 175     // Mark up the local edges
 176
 177     boolList isLocal(mesh.nEdges(), false);
 178
 179     forAll (localEdges, edgeI)
 180     {
 181         isLocal[localEdges[edgeI]] = true;
 182     }
 183
 184     // Get reference to parallel shared edges
 185     const labelList& sharedEdges =
 186         boundaryMesh().globalPatch().localEdgeIndices();
 187
 188     forAll (sharedEdges, sharedI)
 189     {
 190         isLocal[sharedEdges[sharedI]] = true;
 191     }
 192
 193     const labelList& sharedCutEdges =
 194         boundaryMesh().globalPatch().cutEdgeIndices();
 195
 196     forAll (sharedCutEdges, sharedCutI)
 197     {
 198         isLocal[sharedCutEdges[sharedCutI]] = true;
 199     }
 200
 201     labelList localPointLabel(mesh.nPoints(), -1);
 202
 203     forAll (mp, pointI)
 204     {
 205         localPointLabel[mp[pointI]] = pointI;
 206     }
 207
 208     // Count the maximum number of edges coming from the patch
 209     label maxEdgesOnPatch = 0;
 210
 211     forAll (mp, pointI)
 212     {
 213         maxEdgesOnPatch += mesh.nEdgesForPoint(mp[pointI]);
 214     }
 215
 216     // Get reference to addressing
 217     const lduAddressing& ldu = mesh.lduAddr();
 218
 219     // Allocate doubly cut arrays
 220     doubleCutEdgeIndicesPtr_ = new labelList(maxEdgesOnPatch, -1);
 221     labelList& doubleCutEdges = *doubleCutEdgeIndicesPtr_;
 222
 223     doubleCutOwnerPtr_ = new labelList(maxEdgesOnPatch, -1);
 224     labelList& doubleCutOwn = *doubleCutOwnerPtr_;
 225
 226     doubleCutNeighbourPtr_ = new labelList(maxEdgesOnPatch, -1);
 227     labelList& doubleCutNei = *doubleCutNeighbourPtr_;
 228
 229     label nDoubleCut = 0;
 230
 231     const labelList& globalOwner = ldu.lowerAddr();
 232     const labelList& globalNeighbour = ldu.upperAddr();
 233
 234     // Owner side
 235     //~~~~~~~~~~~
 236
 237     // Allocate the array
 238     cutEdgeOwnerIndicesPtr_ = new labelList(maxEdgesOnPatch, -1);
 239     labelList& own = *cutEdgeOwnerIndicesPtr_;
 240     label nOwn = 0;
 241
 242     cutEdgeOwnerStartPtr_ = new labelList(meshPoints().size() + 1, -1);
 243     labelList& ownStart = *cutEdgeOwnerStartPtr_;
 244
 245     // Go through all the local points and get all the edges coming
 246     // from that point.  Check if the edge has been marked as local;
 247     // if not, add it to the list of cut edges.
 248     forAll (mp, pointI)
 249     {
 250         ownStart[pointI] = nOwn;
 251
 252         const label curPointID = mp[pointI];
 253
 254         // get owner edges indices
 255         const label startFaceOwn = ldu.ownerStartAddr()[curPointID];
 256         const label endFaceOwn = ldu.ownerStartAddr()[curPointID + 1];
 257
 258         for
 259         (
 260             label edgeLabel = startFaceOwn;
 261             edgeLabel < endFaceOwn;
 262             edgeLabel++
 263         )
 264         {
 265             if (!isLocal[edgeLabel])
 266             {
 267                 isLocal[edgeLabel] = true;
 268
 269                 if (localPointLabel[globalNeighbour[edgeLabel]] == -1)
 270                 {
 271                     // Singly cut edge
 272                     own[nOwn] = edgeLabel;
 273                     nOwn++;
 274                 }
 275                 else
 276                 {
 277                     // Doubly cut edge
 278                     doubleCutEdges[nDoubleCut] = edgeLabel;
 279                     doubleCutOwn[nDoubleCut] = pointI;
 280                     doubleCutNei[nDoubleCut] =
 281                         localPointLabel[globalNeighbour[edgeLabel]];
 282
 283                     nDoubleCut++;
 284                 }
 285             }
 286         }
 287     }
 288
 289     // reset the size of owner edges
 290     own.setSize(nOwn);
 291
 292     // set the last start label by hand
 293     ownStart[meshPoints().size()] = nOwn;
 294
 295     // Neighbour side
 296     //~~~~~~~~~~~~~~~
 297
 298     // Allocate the array
 299     cutEdgeNeighbourIndicesPtr_ = new labelList(maxEdgesOnPatch, -1);
 300     labelList& nei = *cutEdgeNeighbourIndicesPtr_;
 301     label nNei = 0;
 302
 303     cutEdgeNeighbourStartPtr_ = new labelList(meshPoints().size() + 1, -1);
 304     labelList& neiStart = *cutEdgeNeighbourStartPtr_;
 305
 306     const unallocLabelList& losort = ldu.losortAddr();
 307
 308     // Go through all the local points and get all the edges coming
 309     // from that point.  Check if the edge has been marked as local;
 310     // if not, add it to the list of cut edges.
 311     forAll (mp, pointI)
 312     {
 313         neiStart[pointI] = nNei;
 314
 315         const label curPointID = mp[pointI];
 316
 317         // get neighbour edges indices
 318         const label startFaceNei = ldu.losortStartAddr()[curPointID];
 319         const label endFaceNei = ldu.losortStartAddr()[curPointID + 1];
 320
 321         for
 322         (
 323             label edgeLabel = startFaceNei;
 324             edgeLabel < endFaceNei;
 325             edgeLabel++
 326         )
 327         {
 328             if (!isLocal[losort[edgeLabel]])
 329             {
 330                 isLocal[losort[edgeLabel]] = true;
 331
 332                 if (localPointLabel[globalOwner[losort[edgeLabel]]] == -1)
 333                 {
 334                     // Singly cut edge
 335                     nei[nNei] = losort[edgeLabel];
 336                     nNei++;
 337                 }
 338                 else
 339                 {
 340                     // Doubly cut edge
 341                     doubleCutEdges[nDoubleCut] = losort[edgeLabel];
 342                     doubleCutOwn[nDoubleCut] =
 343                         localPointLabel[globalOwner[losort[edgeLabel]]];
 344                     doubleCutNei[nDoubleCut] = pointI;
 345
 346                     nDoubleCut++;
 347                 }
 348             }
 349         }
 350     }
 351
 352     // reset the size of neighbour edges
 353     nei.setSize(nNei);
 354
 355     // set the last start label by hand
 356     neiStart[meshPoints().size()] = nNei;
 357
 358     // Reset the size of double cut edge data
 359     doubleCutEdges.setSize(nDoubleCut);
 360     doubleCutOwn.setSize(nDoubleCut);
 361     doubleCutNei.setSize(nDoubleCut);
 362 }
 363
 364
 365 void processorTetPolyPatchCellDecomp::calcOwnNeiDoubleMask() const
 366 {
 367     if (ownNeiDoubleMaskPtr_)
 368     {
 369         FatalErrorIn
 370         (
 371             "void processorTetPolyPatchCellDecomp::calcOwnNeiDoubleMask() "
 372             "const"
 373         )   << "ownNeiDoubleMaskPtr_ already allocated"
 374             << abort(FatalError);
 375     }
 376
 377     // Algorithm:
 378     // Make a point lookup field which marks all points belonging to
 379     // Other processor patches.
 380     // Go through all the cut edges and cheeck the points of the edge.
 381     // If the point belongs to another parallel patch, set the mask to zero
 382
 383     const tetPolyMeshCellDecomp& mesh = boundaryMesh().mesh();
 384
 385     boolList otherProcPoint(mesh.nPoints(), false);
 386
 387     forAll (boundaryMesh(), patchI)
 388     {
 389         if
 390         (
 391             isA<processorTetPolyPatchCellDecomp>(boundaryMesh()[patchI])
 392         )
 393         {
 394             // Get the local points and mark up the list
 395             const labelList& mp = boundaryMesh()[patchI].meshPoints();
 396
 397             forAll (mp, i)
 398             {
 399                 otherProcPoint[mp[i]] = true;
 400             }
 401         }
 402     }
 403
 404     // Go through all cut edges and set the mask
 405
 406     // Get matrix addressing
 407     const labelList& mp = meshPoints();
 408
 409     const unallocLabelList& L = mesh.lduAddr().lowerAddr();
 410     const unallocLabelList& U = mesh.lduAddr().upperAddr();
 411
 412     const labelList& cutOwn = cutEdgeOwnerIndices();
 413     const labelList& cutNei = cutEdgeNeighbourIndices();
 414     const labelList& doubleCut = doubleCutEdgeIndices();
 415
 416     ownNeiDoubleMaskPtr_ =
 417         new scalarField
 418         (
 419             cutOwn.size() + cutNei.size() + doubleCut.size(),
 420             1.0
 421         );
 422     scalarField& mask = *ownNeiDoubleMaskPtr_;
 423
 424     label coeffI = 0;
 425
 426     // Owner side
 427     // ~~~~~~~~~~
 428     {
 429         const labelList& cutOwnStart = cutEdgeOwnerStart();
 430
 431         forAll (mp, pointI)
 432         {
 433             label ownIndex = cutOwnStart[pointI];
 434             label endOwn = cutOwnStart[pointI + 1];
 435
 436             for (; ownIndex < endOwn; ownIndex++)
 437             {
 438                 if(otherProcPoint[U[cutOwn[ownIndex]]])
 439                 {
 440                     mask[coeffI] = 0;
 441                 }
 442                 coeffI++;
 443             }
 444         }
 445     }
 446
 447     // Neighbour side
 448     // ~~~~~~~~~~~~~~
 449     {
 450         const labelList& cutNeiStart = cutEdgeNeighbourStart();
 451
 452         forAll (mp, pointI)
 453         {
 454             label neiIndex = cutNeiStart[pointI];
 455             label endNei = cutNeiStart[pointI + 1];
 456
 457             for (; neiIndex < endNei; neiIndex++)
 458             {
 459                 if (otherProcPoint[L[cutNei[neiIndex]]])
 460                 {
 461                     mask[coeffI] = 0;
 462                 }
 463
 464                 coeffI++;
 465             }
 466         }
 467     }
 468
 469     // Doubly cut coefficients
 470     // ~~~~~~~~~~~~~~~~~~~~~~~
 471     // Not needed.  HJ, 18/Apr/2002
 472 }
 473
 474
 475 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
 476
 477 const labelList& processorTetPolyPatchCellDecomp::localEdgeIndices() const
 478 {
 479     if (!localEdgeIndicesPtr_)
 480     {
 481         if (isMaster())
 482         {
 483             localEdgeIndicesPtr_ =
 484                 new labelList(calcProcLocalEdgesIndices(procPolyPatch_));
 485         }
 486         else
 487         {
 488             // Slave side.  Create the reversed patch and pick up its points
 489             // so that the order is correct
 490             const polyPatch& pp = patch();
 491
 492             faceList masterFaces(pp.size());
 493
 494             forAll (pp, faceI)
 495             {
 496                 masterFaces[faceI] = pp[faceI].reverseFace();
 497             }
 498
 499             // To adhere to the interface, the primitive patch is created
 500             // from master faces locally and lost after use
 501             localEdgeIndicesPtr_ =
 502                 new labelList
 503                 (
 504                     calcProcLocalEdgesIndices
 505                     (
 506                         primitiveFacePatch
 507                         (
 508                             masterFaces,
 509                             pp.points()
 510                         )
 511                     )
 512                 );
 513         }
 514     }
 515
 516     return *localEdgeIndicesPtr_;
 517 }
 518
 519
 520 const labelList& processorTetPolyPatchCellDecomp::cutEdgeIndices() const
 521 {
 522     if (!cutEdgeIndicesPtr_)
 523     {
 524         calcCutEdgeIndices();
 525     }
 526
 527     return *cutEdgeIndicesPtr_;
 528 }
 529
 530
 531 const labelList& processorTetPolyPatchCellDecomp::cutEdgeOwnerIndices() const
 532 {
 533     if (!cutEdgeOwnerIndicesPtr_)
 534     {
 535         calcCutEdgeAddressing();
 536     }
 537
 538     return *cutEdgeOwnerIndicesPtr_;
 539 }
 540
 541
 542 const labelList& processorTetPolyPatchCellDecomp::cutEdgeOwnerStart() const
 543 {
 544     if (!cutEdgeOwnerStartPtr_)
 545     {
 546         calcCutEdgeAddressing();
 547     }
 548
 549     return *cutEdgeOwnerStartPtr_;
 550 }
 551
 552
 553 const labelList&
 554 processorTetPolyPatchCellDecomp::cutEdgeNeighbourIndices() const
 555 {
 556     if (!cutEdgeNeighbourIndicesPtr_)
 557     {
 558         calcCutEdgeAddressing();
 559     }
 560
 561     return *cutEdgeNeighbourIndicesPtr_;
 562 }
 563
 564
 565 const labelList& processorTetPolyPatchCellDecomp::cutEdgeNeighbourStart() const
 566 {
 567     if (!cutEdgeNeighbourStartPtr_)
 568     {
 569         calcCutEdgeAddressing();
 570     }
 571
 572     return *cutEdgeNeighbourStartPtr_;
 573 }
 574
 575
 576 const labelList& processorTetPolyPatchCellDecomp::doubleCutEdgeIndices() const
 577 {
 578     if (!doubleCutEdgeIndicesPtr_)
 579     {
 580         calcCutEdgeAddressing();
 581     }
 582
 583     return *doubleCutEdgeIndicesPtr_;
 584 }
 585
 586
 587 const labelList& processorTetPolyPatchCellDecomp::doubleCutOwner() const
 588 {
 589     if (!doubleCutOwnerPtr_)
 590     {
 591         calcCutEdgeAddressing();
 592     }
 593
 594     return *doubleCutOwnerPtr_;
 595 }
 596
 597
 598 const labelList& processorTetPolyPatchCellDecomp::doubleCutNeighbour() const
 599 {
 600     if (!doubleCutNeighbourPtr_)
 601     {
 602         calcCutEdgeAddressing();
 603     }
 604
 605     return *doubleCutNeighbourPtr_;
 606 }
 607
 608
 609 const scalarField& processorTetPolyPatchCellDecomp::ownNeiDoubleMask() const
 610 {
 611     if (!ownNeiDoubleMaskPtr_)
 612     {
 613         calcOwnNeiDoubleMask();
 614     }
 615
 616     return *ownNeiDoubleMaskPtr_;
 617 }
 618
 619
 620 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 621
 622 } // End namespace Foam
 623
 624 // ************************************************************************* //