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-2023 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
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 
38 namespace
39 {
40 
41 // Process connections with the Cuthill-McKee algorithm.
42 // The connections are CompactListList<label> or a labelListList.
43 template<class ConnectionListListType>
44 Foam::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.push_back(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.front();
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.push_back(nbr);
141  weights.push_back(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.push_back(nbrCells[nbrIdx]);
155  }
156  }
157  }
158  }
159 
160  // Now we have new-to-old in newOrder.
161  return newOrder;
162 }
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.push_back(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.front();
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.push_back(nbr);
288  weights.push_back(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.push_back(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 // ************************************************************************* //
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:116
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:40
void push_back(const T &val)
Append an element at the end of the list.
Definition: ListI.H:227
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
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...
A packed storage of objects of type <T> using an offset table for access.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:201
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:405
void push_back(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:555
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:59
constexpr label labelMax
Definition: label.H:55
void resize_nocopy(const label len)
Alter addressable list size, allocating new space if required without necessarily recovering old cont...
Definition: DynamicListI.H:375
Namespace for OpenFOAM.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127
labelList bandCompression(const CompactListList< label > &addressing)
Renumber (mesh) addressing to reduce the band of the matrix, using the Cuthill-McKee algorithm...