lduAddressing.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "lduAddressing.H"
27 #include "demandDrivenData.H"
28 #include "scalarField.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 labelUList& 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 
176 {
177  deleteDemandDrivenData(losortPtr_);
178  deleteDemandDrivenData(ownerStartPtr_);
179  deleteDemandDrivenData(losortStartPtr_);
180 }
181 
182 
183 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
184 
186 {
187  if (!losortPtr_)
188  {
189  calcLosort();
190  }
191 
192  return *losortPtr_;
193 }
194 
195 
197 {
198  if (!ownerStartPtr_)
199  {
200  calcOwnerStart();
201  }
202 
203  return *ownerStartPtr_;
204 }
205 
206 
208 {
209  if (!losortStartPtr_)
210  {
211  calcLosortStart();
212  }
213 
214  return *losortStartPtr_;
215 }
216 
217 
218 // Return edge index given owner and neighbour label
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 labelUList& 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
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 
253 {
254  const labelUList& owner = lowerAddr();
255  const labelUList& neighbour = upperAddr();
256 
257  labelList cellBandwidth(size(), 0);
258 
259  forAll(neighbour, faceI)
260  {
261  label own = owner[faceI];
262  label nei = neighbour[faceI];
263 
264  // Note: mag not necessary for correct (upper-triangular) ordering.
265  label diff = nei-own;
266  cellBandwidth[nei] = max(cellBandwidth[nei], diff);
267  }
268 
269  label bandwidth = max(cellBandwidth);
270 
271  // Do not use field algebra because of conversion label to scalar
272  scalar profile = 0.0;
273  forAll(cellBandwidth, cellI)
274  {
275  profile += 1.0*cellBandwidth[cellI];
276  }
277 
278  return Tuple2<label, scalar>(bandwidth, profile);
279 }
280 
281 
282 // ************************************************************************* //
label size() const
Return number of equations.
virtual const labelUList & upperAddr() const =0
Return upper addressing.
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:47
UList< label > labelUList
Definition: UList.H:63
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
label triIndex(const label a, const label b) const
Return off-diagonal index given owner and neighbour label.
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:407
const labelUList & ownerStartAddr() const
Return owner start addressing.
void deleteDemandDrivenData(DataPtr &dataPtr)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const labelUList & losortAddr() const
Return losort addressing.
Tuple2< label, scalar > band() const
Calculate bandwidth and profile of addressing.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
Template functions to aid in the implementation of demand driven data.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const labelUList & losortStartAddr() const
Return losort start addressing.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
virtual ~lduAddressing()
Destructor.