FreeFOAM The Cross-Platform CFD Toolkit
pairGAMGAgglomerate.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) 1991-2010 OpenCFD Ltd.
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 \*---------------------------------------------------------------------------*/
25 
26 #include "pairGAMGAgglomeration.H"
27 #include <OpenFOAM/lduAddressing.H>
28 
29 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
30 
32 (
33  label& nCoarseCells,
34  const lduAddressing& fineMatrixAddressing,
35  const scalarField& faceWeights
36 )
37 {
38  const label nFineCells = fineMatrixAddressing.size();
39 
40  const unallocLabelList& upperAddr = fineMatrixAddressing.upperAddr();
41  const unallocLabelList& lowerAddr = fineMatrixAddressing.lowerAddr();
42 
43  // For each cell calculate faces
44  labelList cellFaces(upperAddr.size() + lowerAddr.size());
45  labelList cellFaceOffsets(nFineCells + 1);
46 
47  // memory management
48  {
49  labelList nNbrs(nFineCells, 0);
50 
51  forAll (upperAddr, facei)
52  {
53  nNbrs[upperAddr[facei]]++;
54  }
55 
56  forAll (lowerAddr, facei)
57  {
58  nNbrs[lowerAddr[facei]]++;
59  }
60 
61  cellFaceOffsets[0] = 0;
62  forAll (nNbrs, celli)
63  {
64  cellFaceOffsets[celli+1] = cellFaceOffsets[celli] + nNbrs[celli];
65  }
66 
67  // reset the whole list to use as counter
68  nNbrs = 0;
69 
70  forAll (upperAddr, facei)
71  {
72  cellFaces
73  [
74  cellFaceOffsets[upperAddr[facei]] + nNbrs[upperAddr[facei]]
75  ] = facei;
76 
77  nNbrs[upperAddr[facei]]++;
78  }
79 
80  forAll (lowerAddr, facei)
81  {
82  cellFaces
83  [
84  cellFaceOffsets[lowerAddr[facei]] + nNbrs[lowerAddr[facei]]
85  ] = facei;
86 
87  nNbrs[lowerAddr[facei]]++;
88  }
89  }
90 
91 
92  // go through the faces and create clusters
93 
94  tmp<labelField> tcoarseCellMap(new labelField(nFineCells, -1));
95  labelField& coarseCellMap = tcoarseCellMap();
96 
97  nCoarseCells = 0;
98 
99  for (label celli=0; celli<nFineCells; celli++)
100  {
101  if (coarseCellMap[celli] < 0)
102  {
103  label matchFaceNo = -1;
104  scalar maxFaceWeight = -GREAT;
105 
106  // check faces to find ungrouped neighbour with largest face weight
107  for
108  (
109  label faceOs=cellFaceOffsets[celli];
110  faceOs<cellFaceOffsets[celli+1];
111  faceOs++
112  )
113  {
114  label facei = cellFaces[faceOs];
115 
116  // I don't know whether the current cell is owner or neighbour.
117  // Therefore I'll check both sides
118  if
119  (
120  coarseCellMap[upperAddr[facei]] < 0
121  && coarseCellMap[lowerAddr[facei]] < 0
122  && faceWeights[facei] > maxFaceWeight
123  )
124  {
125  // Match found. Pick up all the necessary data
126  matchFaceNo = facei;
127  maxFaceWeight = faceWeights[facei];
128  }
129  }
130 
131  if (matchFaceNo >= 0)
132  {
133  // Make a new group
134  coarseCellMap[upperAddr[matchFaceNo]] = nCoarseCells;
135  coarseCellMap[lowerAddr[matchFaceNo]] = nCoarseCells;
136  nCoarseCells++;
137  }
138  else
139  {
140  // No match. Find the best neighbouring cluster and
141  // put the cell there
142  label clusterMatchFaceNo = -1;
143  scalar clusterMaxFaceCoeff = -GREAT;
144 
145  for
146  (
147  label faceOs=cellFaceOffsets[celli];
148  faceOs<cellFaceOffsets[celli+1];
149  faceOs++
150  )
151  {
152  label facei = cellFaces[faceOs];
153 
154  if (faceWeights[facei] > clusterMaxFaceCoeff)
155  {
156  clusterMatchFaceNo = facei;
157  clusterMaxFaceCoeff = faceWeights[facei];
158  }
159  }
160 
161  if (clusterMatchFaceNo >= 0)
162  {
163  // Add the cell to the best cluster
164  coarseCellMap[celli] = max
165  (
166  coarseCellMap[upperAddr[clusterMatchFaceNo]],
167  coarseCellMap[lowerAddr[clusterMatchFaceNo]]
168  );
169  }
170  }
171  }
172  }
173 
174 
175  // Check that all cells are part of clusters,
176  // if not create single-cell "clusters" for each
177  for (label celli=0; celli<nFineCells; celli++)
178  {
179  if (coarseCellMap[celli] < 0)
180  {
181  coarseCellMap[celli] = nCoarseCells;
182  nCoarseCells++;
183  }
184  }
185 
186  // Reverese the map ordering to improve the next level of agglomeration
187  // (doesn't always help and is sometimes detrimental)
188  nCoarseCells--;
189 
190  forAll (coarseCellMap, celli)
191  {
192  coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
193  }
194 
195  nCoarseCells++;
196 
197  return tcoarseCellMap;
198 }
199 
200 
202 (
203  const lduMesh& mesh,
204  const scalarField& faceWeights
205 )
206 {
207  // Get the finest-level interfaces from the mesh
208  interfaceLevels_.set
209  (
210  0,
211  new lduInterfacePtrsList(mesh.interfaces())
212  );
213 
214  // Start geometric agglomeration from the given faceWeights
215  scalarField* faceWeightsPtr = const_cast<scalarField*>(&faceWeights);
216 
217  // Agglomerate until the required number of cells in the coarsest level
218  // is reached
219 
220  label nPairLevels = 0;
221  label nCreatedLevels = 0;
222 
223  while (nCreatedLevels < maxLevels_ - 1)
224  {
225  label nCoarseCells = -1;
226 
227  tmp<labelField> finalAgglomPtr = agglomerate
228  (
229  nCoarseCells,
230  meshLevel(nCreatedLevels).lduAddr(),
231  *faceWeightsPtr
232  );
233 
234  if (continueAgglomerating(nCoarseCells))
235  {
236  nCells_[nCreatedLevels] = nCoarseCells;
237  restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
238  }
239  else
240  {
241  break;
242  }
243 
244  agglomerateLduAddressing(nCreatedLevels);
245 
246  // Agglomerate the faceWeights field for the next level
247  {
248  scalarField* aggFaceWeightsPtr
249  (
250  new scalarField
251  (
252  meshLevels_[nCreatedLevels].upperAddr().size(),
253  0.0
254  )
255  );
256 
257  restrictFaceField
258  (
259  *aggFaceWeightsPtr,
260  *faceWeightsPtr,
261  nCreatedLevels
262  );
263 
264  if (nCreatedLevels)
265  {
266  delete faceWeightsPtr;
267  }
268 
269  faceWeightsPtr = aggFaceWeightsPtr;
270  }
271 
272  if (nPairLevels % mergeLevels_)
273  {
274  combineLevels(nCreatedLevels);
275  }
276  else
277  {
278  nCreatedLevels++;
279  }
280 
281  nPairLevels++;
282  }
283 
284  // Shrink the storage of the levels to those created
285  compactLevels(nCreatedLevels);
286 
287  // Delete temporary geometry storage
288  if (nCreatedLevels)
289  {
290  delete faceWeightsPtr;
291  }
292 }
293 
294 
295 // ************************ vim: set sw=4 sts=4 et: ************************ //