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 | Copyright (C) 2011-2016 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
61  labelList order;
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
186  labelList order;
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 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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 sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Various functions to operate on Lists.
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...
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
static const label labelMax
Definition: label.H:62
A bit-packed bool list.
labelList bandCompression(const labelListList &addressing)
Renumbers the addressing to reduce the band of the matrix.
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:224