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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
29#include "lduAddressing.H"
30
31// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
32
34(
35 const lduMesh& mesh,
36 const scalarField& faceWeights
37)
38{
39 // Start geometric agglomeration from the given faceWeights
40 scalarField* faceWeightsPtr = const_cast<scalarField*>(&faceWeights);
41
42 // Agglomerate until the required number of cells in the coarsest level
43 // is reached
44
45 label nPairLevels = 0;
46 label nCreatedLevels = 0;
47
48 while (nCreatedLevels < maxLevels_ - 1)
49 {
50 label nCoarseCells = -1;
51
52 tmp<labelField> finalAgglomPtr = agglomerate
53 (
54 nCoarseCells,
55 meshLevel(nCreatedLevels).lduAddr(),
56 *faceWeightsPtr
57 );
58
59 if (continueAgglomerating(finalAgglomPtr().size(), nCoarseCells))
60 {
61 nCells_[nCreatedLevels] = nCoarseCells;
62 restrictAddressing_.set(nCreatedLevels, finalAgglomPtr);
63 }
64 else
65 {
66 break;
67 }
68
69 agglomerateLduAddressing(nCreatedLevels);
70
71 // Agglomerate the faceWeights field for the next level
72 {
73 scalarField* aggFaceWeightsPtr
74 (
75 new scalarField
76 (
77 meshLevels_[nCreatedLevels].upperAddr().size(),
78 0.0
79 )
80 );
81
83 (
84 *aggFaceWeightsPtr,
85 *faceWeightsPtr,
86 nCreatedLevels
87 );
88
89 if (nCreatedLevels)
90 {
91 delete faceWeightsPtr;
92 }
93
94 faceWeightsPtr = aggFaceWeightsPtr;
95 }
96
97 if (nPairLevels % mergeLevels_)
98 {
99 combineLevels(nCreatedLevels);
100 }
101 else
102 {
103 nCreatedLevels++;
104 }
105
106 nPairLevels++;
107 }
108
109 // Shrink the storage of the levels to those created
110 compactLevels(nCreatedLevels);
111
112 // Delete temporary geometry storage
113 if (nCreatedLevels)
114 {
115 delete faceWeightsPtr;
116 }
117}
118
119
120// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
121
123(
124 label& nCoarseCells,
125 const lduAddressing& fineMatrixAddressing,
126 const scalarField& faceWeights
127)
128{
129 const label nFineCells = fineMatrixAddressing.size();
130
131 const labelUList& upperAddr = fineMatrixAddressing.upperAddr();
132 const labelUList& lowerAddr = fineMatrixAddressing.lowerAddr();
133
134 // For each cell calculate faces
135 labelList cellFaces(upperAddr.size() + lowerAddr.size());
136 labelList cellFaceOffsets(nFineCells + 1);
137
138 // memory management
139 {
140 labelList nNbrs(nFineCells, Zero);
141
142 forAll(upperAddr, facei)
143 {
144 nNbrs[upperAddr[facei]]++;
145 }
146
147 forAll(lowerAddr, facei)
148 {
149 nNbrs[lowerAddr[facei]]++;
150 }
151
152 cellFaceOffsets[0] = 0;
153 forAll(nNbrs, celli)
154 {
155 cellFaceOffsets[celli+1] = cellFaceOffsets[celli] + nNbrs[celli];
156 }
157
158 // reset the whole list to use as counter
159 nNbrs = 0;
160
161 forAll(upperAddr, facei)
162 {
163 cellFaces
164 [
165 cellFaceOffsets[upperAddr[facei]] + nNbrs[upperAddr[facei]]
166 ] = facei;
167
168 nNbrs[upperAddr[facei]]++;
169 }
170
171 forAll(lowerAddr, facei)
172 {
173 cellFaces
174 [
175 cellFaceOffsets[lowerAddr[facei]] + nNbrs[lowerAddr[facei]]
176 ] = facei;
177
178 nNbrs[lowerAddr[facei]]++;
179 }
180 }
181
182
183 // go through the faces and create clusters
184
185 tmp<labelField> tcoarseCellMap(new labelField(nFineCells, -1));
186 labelField& coarseCellMap = tcoarseCellMap.ref();
187
188 nCoarseCells = 0;
189 label celli;
190
191 for (label cellfi=0; cellfi<nFineCells; cellfi++)
192 {
193 // Change cell ordering depending on direction for this level
194 celli = forward_ ? cellfi : nFineCells - cellfi - 1;
195
196 if (coarseCellMap[celli] < 0)
197 {
198 label matchFaceNo = -1;
199 scalar maxFaceWeight = -GREAT;
200
201 // check faces to find ungrouped neighbour with largest face weight
202 for
203 (
204 label faceOs=cellFaceOffsets[celli];
205 faceOs<cellFaceOffsets[celli+1];
206 faceOs++
207 )
208 {
209 label facei = cellFaces[faceOs];
210
211 // I don't know whether the current cell is owner or neighbour.
212 // Therefore I'll check both sides
213 if
214 (
215 coarseCellMap[upperAddr[facei]] < 0
216 && coarseCellMap[lowerAddr[facei]] < 0
217 && faceWeights[facei] > maxFaceWeight
218 )
219 {
220 // Match found. Pick up all the necessary data
221 matchFaceNo = facei;
222 maxFaceWeight = faceWeights[facei];
223 }
224 }
225
226 if (matchFaceNo >= 0)
227 {
228 // Make a new group
229 coarseCellMap[upperAddr[matchFaceNo]] = nCoarseCells;
230 coarseCellMap[lowerAddr[matchFaceNo]] = nCoarseCells;
231 nCoarseCells++;
232 }
233 else
234 {
235 // No match. Find the best neighbouring cluster and
236 // put the cell there
237 label clusterMatchFaceNo = -1;
238 scalar clusterMaxFaceCoeff = -GREAT;
239
240 for
241 (
242 label faceOs=cellFaceOffsets[celli];
243 faceOs<cellFaceOffsets[celli+1];
244 faceOs++
245 )
246 {
247 label facei = cellFaces[faceOs];
248
249 if (faceWeights[facei] > clusterMaxFaceCoeff)
250 {
251 clusterMatchFaceNo = facei;
252 clusterMaxFaceCoeff = faceWeights[facei];
253 }
254 }
255
256 if (clusterMatchFaceNo >= 0)
257 {
258 // Add the cell to the best cluster
259 coarseCellMap[celli] = max
260 (
261 coarseCellMap[upperAddr[clusterMatchFaceNo]],
262 coarseCellMap[lowerAddr[clusterMatchFaceNo]]
263 );
264 }
265 }
266 }
267 }
268
269 // Check that all cells are part of clusters,
270 // if not create single-cell "clusters" for each
271 for (label cellfi=0; cellfi<nFineCells; cellfi++)
272 {
273 // Change cell ordering depending on direction for this level
274 celli = forward_ ? cellfi : nFineCells - cellfi - 1;
275
276 if (coarseCellMap[celli] < 0)
277 {
278 coarseCellMap[celli] = nCoarseCells;
279 nCoarseCells++;
280 }
281 }
282
283 if (!forward_)
284 {
285 nCoarseCells--;
286
287 forAll(coarseCellMap, celli)
288 {
289 coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
290 }
291
292 nCoarseCells++;
293 }
294
295 // Reverse the map ordering for the next level
296 // to improve the next level of agglomeration
297 forward_ = !forward_;
298
299 return tcoarseCellMap;
300}
301
302
303// ************************************************************************* //
void agglomerateLduAddressing(const label fineLevelIndex)
Assemble coarse mesh addressing.
const label maxLevels_
Max number of levels.
void compactLevels(const label nCreatedLevels)
Shrink the number of levels to that specified.
PtrList< lduPrimitiveMesh > meshLevels_
Hierarchy of mesh addressing.
bool continueAgglomerating(const label nCells, const label nCoarseCells) const
Check the need for further agglomeration.
labelList nCells_
The number of cells in each level.
void restrictFaceField(Field< Type > &cf, const Field< Type > &ff, const label fineLevelIndex) const
Restrict (integrate by summation) face field.
void combineLevels(const label curLevel)
Combine a level with the previous one.
PtrList< labelField > restrictAddressing_
Cell restriction addressing array.
const lduMesh & meshLevel(const label leveli) const
Return LDU mesh of given level.
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
The class contains the addressing required by the lduMatrix: upper, lower and losort.
virtual const labelUList & upperAddr() const =0
Return upper addressing.
label size() const
Return number of equations.
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
Abstract base class for meshes which provide LDU addressing for the construction of lduMatrix and LDU...
Definition: lduMesh.H:63
void agglomerate(const lduMesh &mesh, const scalarField &faceWeights)
Agglomerate all levels starting from the given face weights.
A class for managing temporary objects.
Definition: tmp.H:65
T & ref() const
Definition: tmpI.H:227
dynamicFvMesh & mesh
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Field< label > labelField
Specialisation of Field<T> for label.
Definition: labelField.H:52
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333