lduAddressing.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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  {
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  {
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  {
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 
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 labelUList& 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
241  << "neighbour " << nbr << " not found for owner " << own << ". "
242  << "Problem with addressing"
243  << abort(FatalError);
244 
245  return -1;
246 }
247 
248 
250 {
251  const labelUList& owner = lowerAddr();
252  const labelUList& neighbour = upperAddr();
253 
254  labelList cellBandwidth(size(), 0);
255 
256  forAll(neighbour, facei)
257  {
258  label own = owner[facei];
259  label nei = neighbour[facei];
260 
261  // Note: mag not necessary for correct (upper-triangular) ordering.
262  label diff = nei-own;
263  cellBandwidth[nei] = max(cellBandwidth[nei], diff);
264  }
265 
266  label bandwidth = max(cellBandwidth);
267 
268  // Do not use field algebra because of conversion label to scalar
269  scalar profile = 0.0;
270  forAll(cellBandwidth, celli)
271  {
272  profile += 1.0*cellBandwidth[celli];
273  }
274 
275  return Tuple2<label, scalar>(bandwidth, profile);
276 }
277 
278 
279 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A 2-tuple for storing two objects of different types.
Definition: Tuple2.H:66
const labelUList & ownerStartAddr() const
Return owner start addressing.
virtual ~lduAddressing()
Destructor.
const labelUList & losortStartAddr() const
Return losort start addressing.
Tuple2< label, scalar > band() const
Calculate bandwidth and profile of addressing.
label size() const
Return number of equations.
virtual const labelUList & upperAddr() const =0
Return upper addressing.
const labelUList & losortAddr() const
Return losort addressing.
label triIndex(const label a, const label b) const
Return off-diagonal index given owner and neighbour label.
Template functions to aid in the implementation of demand driven data.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
volScalarField & b
Definition: createFields.H:25
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
void deleteDemandDrivenData(DataType *&dataPtr)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:409
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
error FatalError
UList< label > labelUList
Definition: UList.H:65