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