FreeFOAM The Cross-Platform CFD Toolkit
extendedCellToFaceStencil.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 
27 #include <OpenFOAM/globalIndex.H>
28 #include <OpenFOAM/syncTools.H>
29 #include <OpenFOAM/SortableList.H>
30 
31 /* * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * */
32 
34 
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
39 (
40  Ostream& os,
41  const labelListList& stencil,
42  const mapDistribute& map
43 )
44 {
45  label sumSize = 0;
46  label nSum = 0;
47  label minSize = labelMax;
48  label maxSize = labelMin;
49 
50  forAll(stencil, i)
51  {
52  const labelList& sCells = stencil[i];
53 
54  if (sCells.size() > 0)
55  {
56  sumSize += sCells.size();
57  nSum++;
58  minSize = min(minSize, sCells.size());
59  maxSize = max(maxSize, sCells.size());
60  }
61  }
62  reduce(sumSize, sumOp<label>());
63  reduce(nSum, sumOp<label>());
64 
65  reduce(minSize, minOp<label>());
66  reduce(maxSize, maxOp<label>());
67 
68  os << "Stencil size :" << nl
69  << " average : " << scalar(sumSize)/nSum << nl
70  << " min : " << minSize << nl
71  << " max : " << maxSize << nl
72  << endl;
73 
74  // Sum all sent data
75  label nSent = 0;
76  label nLocal = 0;
77  forAll(map.subMap(), procI)
78  {
79  if (procI != Pstream::myProcNo())
80  {
81  nSent += map.subMap()[procI].size();
82  }
83  else
84  {
85  nLocal += map.subMap()[procI].size();
86  }
87  }
88 
89  os << "Local data size : " << returnReduce(nLocal, sumOp<label>()) << nl
90  << "Sent data size : " << returnReduce(nSent, sumOp<label>()) << nl
91  << endl;
92 }
93 
94 
97 (
98  const polyMesh& mesh,
99  const globalIndex& globalNumbering,
100  labelListList& faceStencil
101 )
102 {
103  // Convert stencil to schedule
104  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
105 
106  // We now know what information we need from other processors. This needs
107  // to be converted into what information I need to send as well
108  // (mapDistribute)
109 
110 
111  // 1. Construct per processor compact addressing of the global cells
112  // needed. The ones from the local processor are not included since
113  // these are always all needed.
114  List<Map<label> > globalToProc(Pstream::nProcs());
115  {
116  const labelList& procPatchMap = mesh.globalData().procPatchMap();
117  const polyBoundaryMesh& patches = mesh.boundaryMesh();
118 
119  // Presize with (as estimate) size of patch to neighbour.
120  forAll(procPatchMap, procI)
121  {
122  if (procPatchMap[procI] != -1)
123  {
124  globalToProc[procI].resize
125  (
126  patches[procPatchMap[procI]].size()
127  );
128  }
129  }
130 
131  // Collect all (non-local) globalcells/faces needed.
132  forAll(faceStencil, faceI)
133  {
134  const labelList& stencilCells = faceStencil[faceI];
135 
136  forAll(stencilCells, i)
137  {
138  label globalCellI = stencilCells[i];
139  label procI = globalNumbering.whichProcID(stencilCells[i]);
140 
141  if (procI != Pstream::myProcNo())
142  {
143  label nCompact = globalToProc[procI].size();
144  globalToProc[procI].insert(globalCellI, nCompact);
145  }
146  }
147  }
148  // Sort global cells needed (not really necessary)
149  forAll(globalToProc, procI)
150  {
151  if (procI != Pstream::myProcNo())
152  {
153  Map<label>& globalMap = globalToProc[procI];
154 
155  SortableList<label> sorted(globalMap.toc().xfer());
156 
157  forAll(sorted, i)
158  {
159  Map<label>::iterator iter = globalMap.find(sorted[i]);
160  iter() = i;
161  }
162  }
163  }
164 
165 
166  //forAll(globalToProc, procI)
167  //{
168  // Pout<< "From processor:" << procI << " want cells/faces:" << endl;
169  // forAllConstIter(Map<label>, globalToProc[procI], iter)
170  // {
171  // Pout<< " global:" << iter.key()
172  // << " local:" << globalNumbering.toLocal(procI, iter.key())
173  // << endl;
174  // }
175  // Pout<< endl;
176  //}
177  }
178 
179 
180  // 2. The overall compact addressing is
181  // - myProcNo data first (uncompacted)
182  // - all other processors consecutively
183 
184  labelList compactStart(Pstream::nProcs());
185  compactStart[Pstream::myProcNo()] = 0;
186  label nCompact = globalNumbering.localSize();
187  forAll(compactStart, procI)
188  {
189  if (procI != Pstream::myProcNo())
190  {
191  compactStart[procI] = nCompact;
192  nCompact += globalToProc[procI].size();
193  }
194  }
195 
196 
197  // 3. Find out what to receive/send in compact addressing.
198  labelListList recvCompact(Pstream::nProcs());
199  for (label procI = 0; procI < Pstream::nProcs(); procI++)
200  {
201  if (procI != Pstream::myProcNo())
202  {
203  labelList wantedGlobals(globalToProc[procI].size());
204  recvCompact[procI].setSize(globalToProc[procI].size());
205 
206  label i = 0;
207  forAllConstIter(Map<label>, globalToProc[procI], iter)
208  {
209  wantedGlobals[i] = iter.key();
210  recvCompact[procI][i] = compactStart[procI]+iter();
211  i++;
212  }
213 
214  // Send the global cell numbers I need from procI
215  OPstream str(Pstream::blocking, procI);
216  str << wantedGlobals;
217  }
218  else
219  {
220  recvCompact[procI] =
221  compactStart[procI]
222  + identity(globalNumbering.localSize());
223  }
224  }
225  labelListList sendCompact(Pstream::nProcs());
226  for (label procI = 0; procI < Pstream::nProcs(); procI++)
227  {
228  if (procI != Pstream::myProcNo())
229  {
230  // See what neighbour wants to receive (= what I need to send)
231 
232  IPstream str(Pstream::blocking, procI);
233  labelList globalCells(str);
234 
235  labelList& procCompact = sendCompact[procI];
236  procCompact.setSize(globalCells.size());
237 
238  // Convert from globalCells (all on my processor!) into compact
239  // addressing
240  forAll(globalCells, i)
241  {
242  label cellI = globalNumbering.toLocal(globalCells[i]);
243  procCompact[i] = compactStart[Pstream::myProcNo()]+cellI;
244  }
245  }
246  else
247  {
248  sendCompact[procI] = recvCompact[procI];
249  }
250  }
251 
252  // Convert stencil to compact numbering
253  forAll(faceStencil, faceI)
254  {
255  labelList& stencilCells = faceStencil[faceI];
256 
257  forAll(stencilCells, i)
258  {
259  label globalCellI = stencilCells[i];
260  label procI = globalNumbering.whichProcID(globalCellI);
261  if (procI != Pstream::myProcNo())
262  {
263  label localCompact = globalToProc[procI][globalCellI];
264  stencilCells[i] = compactStart[procI]+localCompact;
265  }
266  else
267  {
268  label localCompact = globalNumbering.toLocal(globalCellI);
269  stencilCells[i] = compactStart[procI]+localCompact;
270  }
271 
272  }
273  }
274 
275 
276  // Constuct map for distribution of compact data.
278  (
279  new mapDistribute
280  (
281  nCompact,
282  sendCompact,
283  recvCompact,
284  true // reuse send/recv maps.
285  )
286  );
287 
288  return mapPtr;
289 }
290 
291 
292 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
293 
294 Foam::extendedCellToFaceStencil::extendedCellToFaceStencil(const polyMesh& mesh)
295 :
296  mesh_(mesh)
297 {
298  // Check for transformation - not supported.
299  const polyBoundaryMesh& patches = mesh.boundaryMesh();
300 
301  forAll(patches, patchI)
302  {
303  if (isA<coupledPolyPatch>(patches[patchI]))
304  {
305  const coupledPolyPatch& cpp =
306  refCast<const coupledPolyPatch>(patches[patchI]);
307 
308  if (!cpp.parallel() || cpp.separated())
309  {
311  (
312  "extendedCellToFaceStencil::extendedCellToFaceStencil"
313  "(const polyMesh&)"
314  ) << "Coupled patches with transformations not supported."
315  << endl
316  << "Problematic patch " << cpp.name() << exit(FatalError);
317  }
318  }
319  }
320 }
321 
322 
323 // ************************ vim: set sw=4 sts=4 et: ************************ //