FreeFOAM The Cross-Platform CFD Toolkit
GAMGAgglomerateLduAddressing.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 "GAMGAgglomeration.H"
27 #include <OpenFOAM/GAMGInterface.H>
28 
29 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
30 
32 (
33  const label fineLevelIndex
34 )
35 {
36  const lduMesh& fineMesh = meshLevel(fineLevelIndex);
37  const lduAddressing& fineMeshAddr = fineMesh.lduAddr();
38 
39  const unallocLabelList& upperAddr = fineMeshAddr.upperAddr();
40  const unallocLabelList& lowerAddr = fineMeshAddr.lowerAddr();
41 
42  label nFineFaces = upperAddr.size();
43 
44  // Get restriction map for current level
45  const labelField& restrictMap = restrictAddressing(fineLevelIndex);
46 
47  if (min(restrictMap) == -1)
48  {
49  FatalErrorIn("GAMGAgglomeration::agglomerateLduAddressing")
50  << "min(restrictMap) == -1" << exit(FatalError);
51  }
52 
53  if (restrictMap.size() != fineMeshAddr.size())
54  {
56  (
57  "GAMGAgglomeration::agglomerateLduAddressing"
58  "(const label fineLevelIndex)"
59  ) << "restrict map does not correspond to fine level. " << endl
60  << " Sizes: restrictMap: " << restrictMap.size()
61  << " nEqns: " << fineMeshAddr.size()
62  << abort(FatalError);
63  }
64 
65 
66  // Get the number of coarse cells
67  const label nCoarseCells = nCells_[fineLevelIndex];
68 
69  // Storage for coarse cell neighbours and coefficients
70 
71  // Guess initial maximum number of neighbours in coarse cell
72  label maxNnbrs = 10;
73 
74  // Number of faces for each coarse-cell
75  labelList cCellnFaces(nCoarseCells, 0);
76 
77  // Setup initial packed storage for coarse-cell faces
78  labelList cCellFaces(maxNnbrs*nCoarseCells);
79 
80  // Create face-restriction addressing
81  faceRestrictAddressing_.set(fineLevelIndex, new labelList(nFineFaces));
82  labelList& faceRestrictAddr = faceRestrictAddressing_[fineLevelIndex];
83 
84  // Initial neighbour array (not in upper-triangle order)
85  labelList initCoarseNeighb(nFineFaces);
86 
87  // Counter for coarse faces
88  label nCoarseFaces = 0;
89 
90  // Loop through all fine faces
91  forAll (upperAddr, fineFacei)
92  {
93  label rmUpperAddr = restrictMap[upperAddr[fineFacei]];
94  label rmLowerAddr = restrictMap[lowerAddr[fineFacei]];
95 
96  if (rmUpperAddr == rmLowerAddr)
97  {
98  // For each fine face inside of a coarse cell keep the address
99  // of the cell corresponding to the face in the faceRestrictAddr
100  // as a negative index
101  faceRestrictAddr[fineFacei] = -(rmUpperAddr + 1);
102  }
103  else
104  {
105  // this face is a part of a coarse face
106 
107  label cOwn = rmUpperAddr;
108  label cNei = rmLowerAddr;
109 
110  // get coarse owner and neighbour
111  if (rmUpperAddr > rmLowerAddr)
112  {
113  cOwn = rmLowerAddr;
114  cNei = rmUpperAddr;
115  }
116 
117  // check the neighbour to see if this face has already been found
118  label* ccFaces = &cCellFaces[maxNnbrs*cOwn];
119 
120  bool nbrFound = false;
121  label& ccnFaces = cCellnFaces[cOwn];
122 
123  for (int i=0; i<ccnFaces; i++)
124  {
125  if (initCoarseNeighb[ccFaces[i]] == cNei)
126  {
127  nbrFound = true;
128  faceRestrictAddr[fineFacei] = ccFaces[i];
129  break;
130  }
131  }
132 
133  if (!nbrFound)
134  {
135  if (ccnFaces >= maxNnbrs)
136  {
137  label oldMaxNnbrs = maxNnbrs;
138  maxNnbrs *= 2;
139 
140  cCellFaces.setSize(maxNnbrs*nCoarseCells);
141 
142  forAllReverse(cCellnFaces, i)
143  {
144  label* oldCcNbrs = &cCellFaces[oldMaxNnbrs*i];
145  label* newCcNbrs = &cCellFaces[maxNnbrs*i];
146 
147  for (int j=0; j<cCellnFaces[i]; j++)
148  {
149  newCcNbrs[j] = oldCcNbrs[j];
150  }
151  }
152 
153  ccFaces = &cCellFaces[maxNnbrs*cOwn];
154  }
155 
156  ccFaces[ccnFaces] = nCoarseFaces;
157  initCoarseNeighb[nCoarseFaces] = cNei;
158  faceRestrictAddr[fineFacei] = nCoarseFaces;
159  ccnFaces++;
160 
161  // new coarse face created
162  nCoarseFaces++;
163  }
164  }
165  } // end for all fine faces
166 
167 
168  // Renumber into upper-triangular order
169 
170  // All coarse owner-neighbour storage
171  labelList coarseOwner(nCoarseFaces);
172  labelList coarseNeighbour(nCoarseFaces);
173  labelList coarseFaceMap(nCoarseFaces);
174 
175  label coarseFacei = 0;
176 
177  forAll (cCellnFaces, cci)
178  {
179  label* cFaces = &cCellFaces[maxNnbrs*cci];
180  label ccnFaces = cCellnFaces[cci];
181 
182  for (int i=0; i<ccnFaces; i++)
183  {
184  coarseOwner[coarseFacei] = cci;
185  coarseNeighbour[coarseFacei] = initCoarseNeighb[cFaces[i]];
186  coarseFaceMap[cFaces[i]] = coarseFacei;
187  coarseFacei++;
188  }
189  }
190 
191  forAll(faceRestrictAddr, fineFacei)
192  {
193  if (faceRestrictAddr[fineFacei] >= 0)
194  {
195  faceRestrictAddr[fineFacei] =
196  coarseFaceMap[faceRestrictAddr[fineFacei]];
197  }
198  }
199 
200  // Clear the temporary storage for the coarse cell data
201  cCellnFaces.setSize(0);
202  cCellFaces.setSize(0);
203  initCoarseNeighb.setSize(0);
204  coarseFaceMap.setSize(0);
205 
206 
207  // Create coarse-level interfaces
208 
209  // Get reference to fine-level interfaces
210  const lduInterfacePtrsList& fineInterfaces =
211  interfaceLevels_[fineLevelIndex];
212 
213  // Create coarse-level interfaces
214  interfaceLevels_.set
215  (
216  fineLevelIndex + 1,
217  new lduInterfacePtrsList(fineInterfaces.size())
218  );
219 
220  lduInterfacePtrsList& coarseInterfaces =
221  interfaceLevels_[fineLevelIndex + 1];
222 
223  labelListList coarseInterfaceAddr(fineInterfaces.size());
224 
225  // Initialise transfer of restrict addressing on the interface
226  forAll (fineInterfaces, inti)
227  {
228  if (fineInterfaces.set(inti))
229  {
230  fineInterfaces[inti].initInternalFieldTransfer
231  (
232  Pstream::blocking,
233  restrictMap
234  );
235  }
236  }
237 
238  // Add the coarse level
239  forAll (fineInterfaces, inti)
240  {
241  if (fineInterfaces.set(inti))
242  {
243  coarseInterfaces.set
244  (
245  inti,
246  GAMGInterface::New
247  (
248  fineInterfaces[inti],
249  fineInterfaces[inti].interfaceInternalField(restrictMap),
250  fineInterfaces[inti].internalFieldTransfer
251  (
252  Pstream::blocking,
253  restrictMap
254  )
255  ).ptr()
256  );
257 
258  coarseInterfaceAddr[inti] = coarseInterfaces[inti].faceCells();
259  }
260  }
261 
262 
263  // Set the coarse ldu addressing onto the list
264  meshLevels_.set
265  (
266  fineLevelIndex,
267  new lduPrimitiveMesh
268  (
269  nCoarseCells,
270  coarseOwner,
271  coarseNeighbour,
272  coarseInterfaceAddr,
273  coarseInterfaces,
274  fineMeshAddr.patchSchedule(),
275  true
276  )
277  );
278 }
279 
280 
281 // ************************ vim: set sw=4 sts=4 et: ************************ //