src/tetDecompositionFiniteElement/tetPolyMeshCellDecomp/tetPolyPatches/constraint/global/calcGlobalTetPolyPatchCellDecompAddr.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 global addressing, needed for vector-matrix
  27     multiply on global boundaries.
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "globalTetPolyPatchCellDecomp.H"
  32 #include "tetPolyBoundaryMeshCellDecomp.H"
  33 #include "tetPolyMeshCellDecomp.H"
  34 #include "boolList.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 void globalTetPolyPatchCellDecomp::calcLocalEdgesIndices() const
  44 {
  45     if (debug)
  46     {
  47         Info<< "labelList globalTetPolyPatchCellDecomp::"
  48             << "calcLocalEdgesIndices() const : "
  49             << "calculating local edge indices"
  50             << endl;
  51     }
  52
  53     // Get reference to the mesh
  54     const tetPolyMeshCellDecomp& mesh = boundaryMesh().mesh();
  55
  56     // Get reference to edges
  57     const edgeList& patchEdges = meshEdges();
  58
  59     localEdgeIndicesPtr_ = new labelList(patchEdges.size(), -1);
  60     labelList& localEdgeInd = *localEdgeIndicesPtr_;
  61
  62     const lduAddressing& lduAddr = mesh.lduAddr();
  63
  64     forAll (patchEdges, edgeI)
  65     {
  66         localEdgeInd[edgeI] =
  67             lduAddr.triIndex
  68             (
  69                 patchEdges[edgeI].start(),
  70                 patchEdges[edgeI].end()
  71             );
  72     }
  73
  74 #   ifdef DEBUGtetFemMatrix
  75     if (localEdgeInd.size() > 0 && min(localEdgeInd) < 0)
  76     {
  77         FatalErrorIn
  78         (
  79             "void globalTetPolyPatchCellDecomp::"
  80             "calcLocalEdgesIndices() const"
  81         )   << "Problem in local edge addressing"
  82             << abort(FatalError);
  83     }
  84 #   endif
  85
  86     if (debug)
  87     {
  88         Info<< "void globalTetPolyPatchCellDecomp::"
  89             << "calcLocalEdgesIndices() const : "
  90             << "finished calculating local edge indices"
  91             << endl;
  92     }
  93 }
  94
  95
  96 void globalTetPolyPatchCellDecomp::calcCutEdgeIndices() const
  97 {
  98     if (debug)
  99     {
 100         Info<< "void globalTetPolyPatchCellDecomp::"
 101             << "calcCutEdgeIndices() const : "
 102             << "calculating cut edge indices"
 103             << endl;
 104     }
 105
 106     if (cutEdgeIndicesPtr_)
 107     {
 108         FatalErrorIn
 109         (
 110             "void globalTetPolyPatchCellDecomp::"
 111             "calcCutEdgesIndices() const"
 112         )   << "addressing already allocated"
 113             << abort(FatalError);
 114     }
 115
 116     // Get reference to the mesh
 117     const tetPolyMeshCellDecomp& mesh = boundaryMesh().mesh();
 118
 119     // Get reference to edges
 120     const edgeList& patchCutEdges = meshCutEdges();
 121
 122     cutEdgeIndicesPtr_ = new labelList(patchCutEdges.size(), -1);
 123     labelList& cutEdgeInd = *cutEdgeIndicesPtr_;
 124
 125     const lduAddressing& lduAddr = mesh.lduAddr();
 126
 127     forAll (patchCutEdges, edgeI)
 128     {
 129         cutEdgeInd[edgeI] =
 130             lduAddr.triIndex
 131             (
 132                 patchCutEdges[edgeI].start(),
 133                 patchCutEdges[edgeI].end()
 134             );
 135     }
 136
 137     if (debug)
 138     {
 139         Info<< "void globalTetPolyPatchCellDecomp::"
 140             << "calcCutEdgeIndices() const : "
 141             << "finished calculating cut edge indices"
 142             << endl;
 143     }
 144 }
 145
 146
 147 void globalTetPolyPatchCellDecomp::calcCutEdgeAddressing() const
 148 {
 149     if
 150     (
 151         cutEdgeOwnerIndicesPtr_
 152      || cutEdgeOwnerStartPtr_
 153      || cutEdgeNeighbourIndicesPtr_
 154      || cutEdgeNeighbourStartPtr_
 155      || doubleCutEdgeIndicesPtr_
 156      || doubleCutOwnerPtr_
 157      || doubleCutNeighbourPtr_
 158      || ownNeiDoubleMaskPtr_
 159     )
 160     {
 161         FatalErrorIn
 162         (
 163             "void globalTetPolyPatchCellDecomp::"
 164             "calcCutEdgeAddressing() "
 165             "const"
 166         )   << "addressing already allocated"
 167             << abort(FatalError);
 168     }
 169
 170
 171     // Make a list over all edges in the mesh.  Mark the ones that are local
 172     // to the patch and then collect the rest.
 173     // For doubly cut edges, mark up the local points
 174
 175     const tetPolyMeshCellDecomp& mesh = boundaryMesh().mesh();
 176
 177     // Get reference to mesh points
 178     const labelList& mp = meshPoints();
 179
 180     // Get reference to local edge indices
 181     const labelList& cutEdges = cutEdgeIndices();
 182
 183     // Get the cut edge mask.  It will be used to create the new mask
 184     const scalarField& cutMask = meshCutEdgeMask();
 185
 186     // Mark up the cut edges
 187
 188     labelList cutIndexLookup(mesh.nEdges(), -1);
 189
 190     forAll (cutEdges, edgeI)
 191     {
 192         cutIndexLookup[cutEdges[edgeI]] = edgeI;
 193     }
 194
 195     labelList localPointLabel(mesh.nPoints(), -1);
 196
 197     forAll (mp, pointI)
 198     {
 199         localPointLabel[mp[pointI]] = pointI;
 200     }
 201
 202     // Get reference to addressing
 203     const lduAddressing& ldu = mesh.lduAddr();
 204
 205     const labelList& globalOwner = ldu.lowerAddr();
 206     const labelList& globalNeighbour = ldu.upperAddr();
 207
 208     // Allocate doubly cut arrays
 209     doubleCutEdgeIndicesPtr_ = new labelList(cutEdges.size(), -1);
 210     labelList& doubleCutEdges = *doubleCutEdgeIndicesPtr_;
 211
 212     doubleCutOwnerPtr_ = new labelList(cutEdges.size(), -1);
 213     labelList& doubleCutOwn = *doubleCutOwnerPtr_;
 214
 215     doubleCutNeighbourPtr_ = new labelList(cutEdges.size(), -1);
 216     labelList& doubleCutNei = *doubleCutNeighbourPtr_;
 217
 218     label nDoubleCut = 0;
 219
 220     // Owner side
 221     //~~~~~~~~~~~
 222
 223     // Allocate the array
 224     cutEdgeOwnerIndicesPtr_ = new labelList(cutEdges.size(), -1);
 225     labelList& own = *cutEdgeOwnerIndicesPtr_;
 226     label nOwn = 0;
 227
 228     cutEdgeOwnerStartPtr_ = new labelList(mp.size() + 1, -1);
 229     labelList& ownStart = *cutEdgeOwnerStartPtr_;
 230
 231     // Go through all the local points and get all the edges coming
 232     // from that point.  Check if the edge has been marked as cut;
 233     // add it to the list of cut edges and grab the mask value that goes
 234     // with it.
 235     forAll (mp, pointI)
 236     {
 237         ownStart[pointI] = nOwn;
 238
 239         const label curPointID = mp[pointI];
 240
 241         // get owner edges indices
 242         const label startFaceOwn = ldu.ownerStartAddr()[curPointID];
 243         const label endFaceOwn = ldu.ownerStartAddr()[curPointID + 1];
 244
 245         for
 246         (
 247             label edgeLabel = startFaceOwn;
 248             edgeLabel < endFaceOwn;
 249             edgeLabel++
 250         )
 251         {
 252             if (cutIndexLookup[edgeLabel] >= 0)
 253             {
 254                 if (localPointLabel[globalNeighbour[edgeLabel]] == -1)
 255                 {
 256                     // Singly cut edge
 257                     own[nOwn] = edgeLabel;
 258                     nOwn++;
 259                 }
 260                 else
 261                 {
 262                     // Doubly cut edge
 263                     doubleCutEdges[nDoubleCut] = edgeLabel;
 264                     doubleCutOwn[nDoubleCut] = pointI;
 265                     doubleCutNei[nDoubleCut] =
 266                         localPointLabel[globalNeighbour[edgeLabel]];
 267
 268                     nDoubleCut++;
 269                 }
 270             }
 271         }
 272     }
 273
 274     // reset the size of owner edges
 275     own.setSize(nOwn);
 276
 277     // set the last start label by hand
 278     ownStart[mp.size()] = nOwn;
 279
 280     // Neighbour side
 281     //~~~~~~~~~~~~~~~
 282
 283     // Allocate the array
 284     cutEdgeNeighbourIndicesPtr_ = new labelList(cutEdges.size(), -1);
 285     labelList& nei = *cutEdgeNeighbourIndicesPtr_;
 286     label nNei = 0;
 287
 288     cutEdgeNeighbourStartPtr_ = new labelList(mp.size() + 1, -1);
 289     labelList& neiStart = *cutEdgeNeighbourStartPtr_;
 290
 291     const unallocLabelList& losort = ldu.losortAddr();
 292
 293     // Go through all the local points and get all the edges coming
 294     // from that point.  Check if the edge has been marked as local;
 295     // if not, add it to the list of cut edges.
 296     forAll (mp, pointI)
 297     {
 298         neiStart[pointI] = nNei;
 299
 300         const label curPointID = mp[pointI];
 301
 302         // get neighbour edges indices
 303         const label startFaceNei = ldu.losortStartAddr()[curPointID];
 304         const label endFaceNei = ldu.losortStartAddr()[curPointID + 1];
 305
 306         for
 307         (
 308             label edgeLabel = startFaceNei;
 309             edgeLabel < endFaceNei;
 310             edgeLabel++
 311         )
 312         {
 313             if (cutIndexLookup[losort[edgeLabel]] >= 0)
 314             {
 315                 if (localPointLabel[globalOwner[losort[edgeLabel]]] == -1)
 316                 {
 317                     // Singly cut edge
 318                     nei[nNei] = losort[edgeLabel];
 319                     nNei++;
 320                 }
 321                 else
 322                 {
 323                     // Doubly cut edge
 324                     doubleCutEdges[nDoubleCut] = losort[edgeLabel];
 325                     doubleCutOwn[nDoubleCut] =
 326                         localPointLabel[globalOwner[losort[edgeLabel]]];
 327                     doubleCutNei[nDoubleCut] = pointI;
 328
 329                     nDoubleCut++;
 330                 }
 331             }
 332         }
 333     }
 334
 335     // reset the size of neighbour edges
 336     nei.setSize(nNei);
 337
 338     // set the last start label by hand
 339     neiStart[mp.size()] = nNei;
 340
 341     // Reset the size of double cut edge data
 342     doubleCutEdges.setSize(nDoubleCut);
 343     doubleCutOwn.setSize(nDoubleCut);
 344     doubleCutNei.setSize(nDoubleCut);
 345
 346     // Allocate the reordered mask that corresponds to the owner-neighbour
 347     // coefficient ordering
 348     ownNeiDoubleMaskPtr_ =
 349         new scalarField
 350         (
 351             own.size() + nei.size()
 352           + doubleCutOwn.size() + doubleCutNei.size()
 353         );
 354     scalarField& ownNeiDoubleMask = *ownNeiDoubleMaskPtr_;
 355
 356     label nMask = 0;
 357
 358     forAll (own, i)
 359     {
 360         ownNeiDoubleMask[nMask] = cutMask[cutIndexLookup[own[i]]];
 361         nMask++;
 362     }
 363
 364     forAll (nei, i)
 365     {
 366         ownNeiDoubleMask[nMask] = cutMask[cutIndexLookup[nei[i]]];
 367         nMask++;
 368     }
 369
 370     forAll (doubleCutOwn, i)
 371     {
 372         ownNeiDoubleMask[nMask] = cutMask[cutIndexLookup[doubleCutOwn[i]]];
 373         nMask++;
 374
 375         ownNeiDoubleMask[nMask] = cutMask[cutIndexLookup[doubleCutNei[i]]];
 376         nMask++;
 377     }
 378 }
 379
 380
 381 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 382
 383 const labelList& globalTetPolyPatchCellDecomp::localEdgeIndices() const
 384 {
 385     if (!localEdgeIndicesPtr_)
 386     {
 387         calcLocalEdgesIndices();
 388     }
 389
 390     return *localEdgeIndicesPtr_;
 391 }
 392
 393
 394 const labelList&
 395 globalTetPolyPatchCellDecomp::cutEdgeIndices() const
 396 {
 397     if (!cutEdgeIndicesPtr_)
 398     {
 399         calcCutEdgeIndices();
 400     }
 401
 402     return *cutEdgeIndicesPtr_;
 403 }
 404
 405
 406 const labelList&
 407 globalTetPolyPatchCellDecomp::cutEdgeOwnerIndices() const
 408 {
 409     if (!cutEdgeOwnerIndicesPtr_)
 410     {
 411         calcCutEdgeAddressing();
 412     }
 413
 414     return *cutEdgeOwnerIndicesPtr_;
 415 }
 416
 417
 418 const labelList&
 419 globalTetPolyPatchCellDecomp::cutEdgeOwnerStart() const
 420 {
 421     if (!cutEdgeOwnerStartPtr_)
 422     {
 423         calcCutEdgeAddressing();
 424     }
 425
 426     return *cutEdgeOwnerStartPtr_;
 427 }
 428
 429
 430 const labelList&
 431 globalTetPolyPatchCellDecomp::cutEdgeNeighbourIndices() const
 432 {
 433     if (!cutEdgeNeighbourIndicesPtr_)
 434     {
 435         calcCutEdgeAddressing();
 436     }
 437
 438     return *cutEdgeNeighbourIndicesPtr_;
 439 }
 440
 441
 442 const labelList&
 443 globalTetPolyPatchCellDecomp::cutEdgeNeighbourStart() const
 444 {
 445     if (!cutEdgeNeighbourStartPtr_)
 446     {
 447         calcCutEdgeAddressing();
 448     }
 449
 450     return *cutEdgeNeighbourStartPtr_;
 451 }
 452
 453
 454 const labelList&
 455 globalTetPolyPatchCellDecomp::doubleCutEdgeIndices() const
 456 {
 457     if (!doubleCutEdgeIndicesPtr_)
 458     {
 459         calcCutEdgeAddressing();
 460     }
 461
 462     return *doubleCutEdgeIndicesPtr_;
 463 }
 464
 465
 466 const labelList& globalTetPolyPatchCellDecomp::doubleCutOwner() const
 467 {
 468     if (!doubleCutOwnerPtr_)
 469     {
 470         calcCutEdgeAddressing();
 471     }
 472
 473     return *doubleCutOwnerPtr_;
 474 }
 475
 476
 477 const labelList&
 478 globalTetPolyPatchCellDecomp::doubleCutNeighbour() const
 479 {
 480     if (!doubleCutNeighbourPtr_)
 481     {
 482         calcCutEdgeAddressing();
 483     }
 484
 485     return *doubleCutNeighbourPtr_;
 486 }
 487
 488
 489 const scalarField&
 490 globalTetPolyPatchCellDecomp::ownNeiDoubleMask() const
 491 {
 492     if (!ownNeiDoubleMaskPtr_)
 493     {
 494         calcCutEdgeAddressing();
 495     }
 496
 497     return *ownNeiDoubleMaskPtr_;
 498 }
 499
 500
 501 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 502
 503 } // End namespace Foam
 504
 505 // ************************************************************************* //