bandCompression.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-2021 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 Description
25  The function renumbers the addressing such that the band of the
26  matrix is reduced. The algorithm uses a simple search through the
27  neighbour list
28 
29  See http://en.wikipedia.org/wiki/Cuthill-McKee_algorithm
30 
31 \*---------------------------------------------------------------------------*/
32 
33 #include "bandCompression.H"
34 #include "SLList.H"
35 #include "IOstreams.H"
36 #include "DynamicList.H"
37 #include "ListOps.H"
38 #include "PackedBoolList.H"
39 
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 
43 {
44  labelList newOrder(cellCellAddressing.size());
45 
46  // the business bit of the renumbering
47  SLList<label> nextCell;
48 
49  PackedBoolList visited(cellCellAddressing.size());
50 
51  label cellInOrder = 0;
52 
53 
54  // Work arrays. Kept outside of loop to minimise allocations.
55  // - neighbour cells
56  DynamicList<label> nbrs;
57  // - corresponding weights
58  DynamicList<label> weights;
59 
60  // - ordering
62 
63 
64  while (true)
65  {
66  // For a disconnected region find the lowest connected cell.
67 
68  label currentCell = -1;
69  label minWeight = labelMax;
70 
71  forAll(visited, celli)
72  {
73  // find the lowest connected cell that has not been visited yet
74  if (!visited[celli])
75  {
76  if (cellCellAddressing[celli].size() < minWeight)
77  {
78  minWeight = cellCellAddressing[celli].size();
79  currentCell = celli;
80  }
81  }
82  }
83 
84 
85  if (currentCell == -1)
86  {
87  break;
88  }
89 
90 
91  // Starting from currentCell walk breadth-first
92 
93 
94  // use this cell as a start
95  nextCell.append(currentCell);
96 
97  // loop through the nextCell list. Add the first cell into the
98  // cell order if it has not already been visited and ask for its
99  // neighbours. If the neighbour in question has not been visited,
100  // add it to the end of the nextCell list
101 
102  while (nextCell.size())
103  {
104  currentCell = nextCell.removeHead();
105 
106  if (!visited[currentCell])
107  {
108  visited[currentCell] = 1;
109 
110  // add into cellOrder
111  newOrder[cellInOrder] = currentCell;
112  cellInOrder++;
113 
114  // find if the neighbours have been visited
115  const labelList& neighbours = cellCellAddressing[currentCell];
116 
117  // Add in increasing order of connectivity
118 
119  // 1. Count neighbours of unvisited neighbours
120  nbrs.clear();
121  weights.clear();
122 
123  forAll(neighbours, nI)
124  {
125  label nbr = neighbours[nI];
126  if (!visited[nbr])
127  {
128  // not visited, add to the list
129  nbrs.append(nbr);
130  weights.append(cellCellAddressing[nbr].size());
131  }
132  }
133  // 2. Sort in ascending order
134  sortedOrder(weights, order);
135  // 3. Add in sorted order
136  forAll(order, i)
137  {
138  nextCell.append(nbrs[i]);
139  }
140  }
141  }
142  }
143 
144  return newOrder;
145 }
146 
147 
149 (
150  const labelList& cellCells,
151  const labelList& offsets
152 )
153 {
154  // Count number of neighbours
155  labelList numNbrs(offsets.size()-1, 0);
156  forAll(numNbrs, celli)
157  {
158  label start = offsets[celli];
159  label end = offsets[celli+1];
160 
161  for (label facei = start; facei < end; facei++)
162  {
163  numNbrs[celli]++;
164  numNbrs[cellCells[facei]]++;
165  }
166  }
167 
168 
169  labelList newOrder(offsets.size()-1);
170 
171  // the business bit of the renumbering
172  SLList<label> nextCell;
173 
174  PackedBoolList visited(offsets.size()-1);
175 
176  label cellInOrder = 0;
177 
178 
179  // Work arrays. Kept outside of loop to minimise allocations.
180  // - neighbour cells
181  DynamicList<label> nbrs;
182  // - corresponding weights
183  DynamicList<label> weights;
184 
185  // - ordering
187 
188 
189  while (true)
190  {
191  // For a disconnected region find the lowest connected cell.
192 
193  label currentCell = -1;
194  label minWeight = labelMax;
195 
196  forAll(visited, celli)
197  {
198  // find the lowest connected cell that has not been visited yet
199  if (!visited[celli])
200  {
201  if (numNbrs[celli] < minWeight)
202  {
203  minWeight = numNbrs[celli];
204  currentCell = celli;
205  }
206  }
207  }
208 
209 
210  if (currentCell == -1)
211  {
212  break;
213  }
214 
215 
216  // Starting from currentCell walk breadth-first
217 
218 
219  // use this cell as a start
220  nextCell.append(currentCell);
221 
222  // loop through the nextCell list. Add the first cell into the
223  // cell order if it has not already been visited and ask for its
224  // neighbours. If the neighbour in question has not been visited,
225  // add it to the end of the nextCell list
226 
227  while (nextCell.size())
228  {
229  currentCell = nextCell.removeHead();
230 
231  if (!visited[currentCell])
232  {
233  visited[currentCell] = 1;
234 
235  // add into cellOrder
236  newOrder[cellInOrder] = currentCell;
237  cellInOrder++;
238 
239  // Add in increasing order of connectivity
240 
241  // 1. Count neighbours of unvisited neighbours
242  nbrs.clear();
243  weights.clear();
244 
245  label start = offsets[currentCell];
246  label end = offsets[currentCell+1];
247 
248  for (label facei = start; facei < end; facei++)
249  {
250  label nbr = cellCells[facei];
251  if (!visited[nbr])
252  {
253  // not visited, add to the list
254  nbrs.append(nbr);
255  weights.append(numNbrs[nbr]);
256  }
257  }
258  // 2. Sort in ascending order
259  sortedOrder(weights, order);
260  // 3. Add in sorted order
261  forAll(order, i)
262  {
263  nextCell.append(nbrs[i]);
264  }
265  }
266  }
267  }
268 
269  return newOrder;
270 }
271 
272 
273 // ************************************************************************* //
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
Various functions to operate on Lists.
Non-intrusive singly-linked list.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
The bandCompression function renumbers the addressing such that the band of the matrix is reduced....
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
Template class for non-intrusive linked lists.
Definition: LList.H:76
T removeHead()
Remove and return head.
Definition: LList.H:178
void append(const T &a)
Add at tail of list.
Definition: LList.H:172
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
A bit-packed bool list.
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
int order(const scalar s)
labelList bandCompression(const labelListList &addressing)
Renumbers the addressing to reduce the band of the matrix.
static const label labelMax
Definition: label.H:62
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.