bandCompression.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "bandCompression.H"
30#include "bitSet.H"
31#include "CircularBuffer.H"
32#include "CompactListList.H"
33#include "DynamicList.H"
34#include "IOstreams.H"
35
36// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
37
38namespace
39{
40
41// Process connections with the Cuthill-McKee algorithm.
42// The connections are CompactListList<label> or a labelListList.
43template<class ConnectionListListType>
44Foam::labelList cuthill_mckee_algorithm
45(
46 const ConnectionListListType& cellCellAddressing
47)
48{
49 using namespace Foam;
50
51 const label nOldCells(cellCellAddressing.size());
52
53 // Which cells are visited/unvisited
54 bitSet unvisited(nOldCells, true);
55
56 // The new output order
57 labelList newOrder(nOldCells);
58
59
60 // Various work arrays
61 // ~~~~~~~~~~~~~~~~~~~
62
63 // Neighbour cells
64 DynamicList<label> nbrCells;
65
66 // Neighbour ordering
67 DynamicList<label> nbrOrder;
68
69 // Corresponding weights for neighbour cells
70 DynamicList<label> weights;
71
72 // FIFO buffer for renumbering.
73 CircularBuffer<label> queuedCells(1024);
74
75 label cellInOrder = 0;
76
77 while (true)
78 {
79 // For a disconnected region find the lowest connected cell.
80 label currCelli = -1;
81 label minCount = labelMax;
82
83 for (const label celli : unvisited)
84 {
85 const label nbrCount = cellCellAddressing[celli].size();
86
87 if (minCount > nbrCount)
88 {
89 minCount = nbrCount;
90 currCelli = celli;
91 }
92 }
93
94 if (currCelli == -1)
95 {
96 break;
97 }
98
99
100 // Starting from currCelli - walk breadth-first
101
102 queuedCells.append(currCelli);
103
104 // Loop through queuedCells list. Add the first cell into the
105 // cell order if it has not already been visited and ask for its
106 // neighbours. If the neighbour in question has not been visited,
107 // add it to the end of the queuedCells list
108
109 while (!queuedCells.empty())
110 {
111 // Process as FIFO
112 currCelli = queuedCells.first();
113 queuedCells.pop_front();
114
115 if (unvisited.test(currCelli))
116 {
117 // First visit...
118 unvisited.unset(currCelli);
119
120 // Add into cellOrder
121 newOrder[cellInOrder] = currCelli;
122 ++cellInOrder;
123
124 // Add in increasing order of connectivity
125
126 // 1. Count neighbours of unvisited neighbours
127 nbrCells.clear();
128 weights.clear();
129
130 // Find if the neighbours have been visited
131 const auto& neighbours = cellCellAddressing[currCelli];
132
133 for (const label nbr : neighbours)
134 {
135 const label nbrCount = cellCellAddressing[nbr].size();
136
137 if (unvisited.test(nbr))
138 {
139 // Not visited (or removed), add to the list
140 nbrCells.append(nbr);
141 weights.append(nbrCount);
142 }
143 }
144
145 // Resize DynamicList prior to sortedOrder
146 nbrOrder.resize_nocopy(weights.size());
147
148 // 2. Ascending order
149 Foam::sortedOrder(weights, nbrOrder);
150
151 // 3. Add to FIFO in sorted order
152 for (const label nbrIdx : nbrOrder)
153 {
154 queuedCells.append(nbrCells[nbrIdx]);
155 }
156 }
157 }
158 }
159
160 // Now we have new-to-old in newOrder.
161 return newOrder;
162}
163
164} // End anonymous namespace
165
166
167// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
168
170(
171 const labelUList& cellCells,
172 const labelUList& offsets
173)
174{
175 // Protect against zero-sized offset list
176 const label nOldCells = max(0, (offsets.size()-1));
177
178 // Count number of neighbours
179 labelList numNbrs(nOldCells, Zero);
180 for (label celli = 0; celli < nOldCells; ++celli)
181 {
182 const label beg = offsets[celli];
183 const label end = offsets[celli+1];
184
185 for (label idx = beg; idx < end; ++idx)
186 {
187 ++numNbrs[celli];
188 ++numNbrs[cellCells[idx]];
189 }
190 }
191
192
193 // Which cells are visited/unvisited
194 bitSet unvisited(nOldCells, true);
195
196 // The new output order
197 labelList newOrder(nOldCells);
198
199
200 // Various work arrays
201 // ~~~~~~~~~~~~~~~~~~~
202
203 // Neighbour cells
204 DynamicList<label> nbrCells;
205
206 // Neighbour ordering
207 DynamicList<label> nbrOrder;
208
209 // Corresponding weights for neighbour cells
210 DynamicList<label> weights;
211
212 // FIFO buffer for renumbering.
213 CircularBuffer<label> queuedCells(1024);
214
215
216 label cellInOrder = 0;
217
218 while (true)
219 {
220 // Find lowest connected cell that has not been visited yet
221 label currCelli = -1;
222 label minCount = labelMax;
223
224 for (const label celli : unvisited)
225 {
226 const label nbrCount = numNbrs[celli];
227
228 if (minCount > nbrCount)
229 {
230 minCount = nbrCount;
231 currCelli = celli;
232 }
233 }
234
235 if (currCelli == -1)
236 {
237 break;
238 }
239
240
241 // Starting from currCellii - walk breadth-first
242
243 queuedCells.append(currCelli);
244
245 // loop through the nextCell list. Add the first cell into the
246 // cell order if it has not already been visited and ask for its
247 // neighbours. If the neighbour in question has not been visited,
248 // add it to the end of the nextCell list
249
250 // Loop through queuedCells list. Add the first cell into the
251 // cell order if it has not already been visited and ask for its
252 // neighbours. If the neighbour in question has not been visited,
253 // add it to the end of the queuedCells list
254
255 while (!queuedCells.empty())
256 {
257 // Process as FIFO
258 currCelli = queuedCells.first();
259 queuedCells.pop_front();
260
261 if (unvisited.test(currCelli))
262 {
263 // First visit...
264 unvisited.unset(currCelli);
265
266 // Add into cellOrder
267 newOrder[cellInOrder] = currCelli;
268 ++cellInOrder;
269
270 // Add in increasing order of connectivity
271
272 // 1. Count neighbours of unvisited neighbours
273 nbrCells.clear();
274 weights.clear();
275
276 const label beg = offsets[currCelli];
277 const label end = offsets[currCelli+1];
278
279 for (label idx = beg; idx < end; ++idx)
280 {
281 const label nbr = cellCells[idx];
282 const label nbrCount = numNbrs[nbr];
283
284 if (unvisited.test(nbr))
285 {
286 // Not visited (or removed), add to the list
287 nbrCells.append(nbr);
288 weights.append(nbrCount);
289 }
290 }
291
292 // Resize DynamicList prior to sortedOrder
293 nbrOrder.resize_nocopy(weights.size());
294
295 // 2. Ascending order
296 Foam::sortedOrder(weights, nbrOrder);
297
298 // 3. Add to FIFO in sorted order
299 for (const label nbrIdx : nbrOrder)
300 {
301 queuedCells.append(nbrCells[nbrIdx]);
302 }
303 }
304 }
305 }
306
307 // Now we have new-to-old in newOrder.
308
309 return newOrder;
310}
311
312
313// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
314
316(
317 const CompactListList<label>& cellCellAddressing
318)
319{
320 return cuthill_mckee_algorithm(cellCellAddressing);
321}
322
323
325(
326 const labelListList& cellCellAddressing
327)
328{
329 return cuthill_mckee_algorithm(cellCellAddressing);
330}
331
332
333// ************************************************************************* //
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
Y[inertIndex] max(0.0)
The bandCompression function renumbers the addressing such that the band of the matrix is reduced....
A simple list of objects of type <T> that is intended to be used as a circular buffer (eg,...
void pop_front(label n=1)
Shrink by moving the front of the buffer 1 or more times.
T & first()
Access the first element (front). Requires !empty().
bool empty() const noexcept
Empty or exhausted buffer.
void append(const T &val)
Copy append an element to the end of the buffer.
A packed storage unstructured matrix of objects of type <T> using an offset table for access.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:391
void resize_nocopy(const label len)
Definition: DynamicListI.H:363
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
labelList bandCompression(const CompactListList< label > &addressing)
Namespace for OpenFOAM.
constexpr label labelMax
Definition: label.H:61
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.