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