src/foam/matrices/lduMatrix/lduAddressing/lduAddressing.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 "lduAddressing.H"
  27 #include "demandDrivenData.H"
  28
  29 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  30
  31 void Foam::lduAddressing::calcLosort() const
  32 {
  33     if (losortPtr_)
  34     {
  35         FatalErrorIn("lduAddressing::calcLosort() const")
  36             << "losort already calculated"
  37             << abort(FatalError);
  38     }
  39
  40     // Scan the neighbour list to find out how many times the cell
  41     // appears as a neighbour of the face. Done this way to avoid guessing
  42     // and resizing list
  43     labelList nNbrOfFace(size(), 0);
  44
  45     const unallocLabelList& nbr = upperAddr();
  46
  47     forAll (nbr, nbrI)
  48     {
  49         nNbrOfFace[nbr[nbrI]]++;
  50     }
  51
  52     // Create temporary neighbour addressing
  53     labelListList cellNbrFaces(size());
  54
  55     forAll (cellNbrFaces, cellI)
  56     {
  57         cellNbrFaces[cellI].setSize(nNbrOfFace[cellI]);
  58     }
  59
  60     // Reset the list of number of neighbours to zero
  61     nNbrOfFace = 0;
  62
  63     // Scatter the neighbour faces
  64     forAll (nbr, nbrI)
  65     {
  66         cellNbrFaces[nbr[nbrI]][nNbrOfFace[nbr[nbrI]]] = nbrI;
  67
  68         nNbrOfFace[nbr[nbrI]]++;
  69     }
  70
  71     // Gather the neighbours into the losort array
  72     losortPtr_ = new labelList(nbr.size(), -1);
  73
  74     labelList& lst = *losortPtr_;
  75
  76     // Set counter for losort
  77     label lstI = 0;
  78
  79     forAll (cellNbrFaces, cellI)
  80     {
  81         const labelList& curNbr = cellNbrFaces[cellI];
  82
  83         forAll (curNbr, curNbrI)
  84         {
  85             lst[lstI] = curNbr[curNbrI];
  86             lstI++;
  87         }
  88     }
  89 }
  90
  91
  92 void Foam::lduAddressing::calcOwnerStart() const
  93 {
  94     if (ownerStartPtr_)
  95     {
  96         FatalErrorIn("lduAddressing::calcOwnerStart() const")
  97             << "owner start already calculated"
  98             << abort(FatalError);
  99     }
 100
 101     const labelList& own = lowerAddr();
 102
 103     ownerStartPtr_ = new labelList(size() + 1, own.size());
 104
 105     labelList& ownStart = *ownerStartPtr_;
 106
 107     // Set up first lookup by hand
 108     ownStart[0] = 0;
 109     label nOwnStart = 0;
 110     label i = 1;
 111
 112     forAll (own, faceI)
 113     {
 114         label curOwn = own[faceI];
 115
 116         if (curOwn > nOwnStart)
 117         {
 118             while (i <= curOwn)
 119             {
 120                 ownStart[i++] = faceI;
 121             }
 122
 123             nOwnStart = curOwn;
 124         }
 125     }
 126 }
 127
 128
 129 void Foam::lduAddressing::calcLosortStart() const
 130 {
 131     if (losortStartPtr_)
 132     {
 133         FatalErrorIn("lduAddressing::calcLosortStart() const")
 134             << "losort start already calculated"
 135             << abort(FatalError);
 136     }
 137
 138     losortStartPtr_ = new labelList(size() + 1, 0);
 139
 140     labelList& lsrtStart = *losortStartPtr_;
 141
 142     const labelList& nbr = upperAddr();
 143
 144     const labelList& lsrt = losortAddr();
 145
 146     // Set up first lookup by hand
 147     lsrtStart[0] = 0;
 148     label nLsrtStart = 0;
 149     label i = 0;
 150
 151     forAll (lsrt, faceI)
 152     {
 153         // Get neighbour
 154         const label curNbr = nbr[lsrt[faceI]];
 155
 156         if (curNbr > nLsrtStart)
 157         {
 158             while (i <= curNbr)
 159             {
 160                 lsrtStart[i++] = faceI;
 161             }
 162
 163             nLsrtStart = curNbr;
 164         }
 165     }
 166
 167     // Set up last lookup by hand
 168     lsrtStart[size()] = nbr.size();
 169 }
 170
 171
 172 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 173
 174 Foam::lduAddressing::~lduAddressing()
 175 {
 176     deleteDemandDrivenData(losortPtr_);
 177     deleteDemandDrivenData(ownerStartPtr_);
 178     deleteDemandDrivenData(losortStartPtr_);
 179 }
 180
 181
 182 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 183
 184 const Foam::unallocLabelList& Foam::lduAddressing::losortAddr() const
 185 {
 186     if (!losortPtr_)
 187     {
 188         calcLosort();
 189     }
 190
 191     return *losortPtr_;
 192 }
 193
 194
 195 const Foam::unallocLabelList& Foam::lduAddressing::ownerStartAddr() const
 196 {
 197     if (!ownerStartPtr_)
 198     {
 199         calcOwnerStart();
 200     }
 201
 202     return *ownerStartPtr_;
 203 }
 204
 205
 206 const Foam::unallocLabelList& Foam::lduAddressing::losortStartAddr() const
 207 {
 208     if (!losortStartPtr_)
 209     {
 210         calcLosortStart();
 211     }
 212
 213     return *losortStartPtr_;
 214 }
 215
 216
 217 // Return edge index given owner and neighbour label
 218 Foam::label Foam::lduAddressing::triIndex(const label a, const label b) const
 219 {
 220     label own = min(a, b);
 221
 222     label nbr = max(a, b);
 223
 224     label startLabel = ownerStartAddr()[own];
 225
 226     label endLabel = ownerStartAddr()[own + 1];
 227
 228     const unallocLabelList& neighbour = upperAddr();
 229
 230     for (label i = startLabel; i < endLabel; i++)
 231     {
 232         if (neighbour[i] == nbr)
 233         {
 234             return i;
 235         }
 236     }
 237
 238     // If neighbour has not been found, something has gone seriously
 239     // wrong with the addressing mechanism
 240     FatalErrorIn
 241     (
 242         "lduAddressing::triIndex(const label owner, const label nbr) const"
 243     )   << "neighbour " << nbr << " not found for owner " << own << ". "
 244         << "Problem with addressing"
 245         << abort(FatalError);
 246
 247     return -1;
 248 }
 249
 250
 251 // ************************************************************************* //