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