FreeFOAM The Cross-Platform CFD Toolkit
FECCellToFaceStencil.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 "FECCellToFaceStencil.H"
27 #include <OpenFOAM/syncTools.H>
29 //#include "meshTools.H"
30 //#include "OFstream.H"
31 //#include "Time.H"
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 // Calculates per edge the neighbour data (= edgeCells)
36 void Foam::FECCellToFaceStencil::calcEdgeBoundaryData
37 (
38  const boolList& isValidBFace,
39  const labelList& boundaryEdges,
40  EdgeMap<labelList>& neiGlobal
41 ) const
42 {
43  neiGlobal.resize(2*boundaryEdges.size());
44 
45  labelHashSet edgeGlobals;
46 
47  forAll(boundaryEdges, i)
48  {
49  label edgeI = boundaryEdges[i];
50 
51  neiGlobal.insert
52  (
53  mesh().edges()[edgeI],
55  (
56  isValidBFace,
57  mesh().edgeFaces(edgeI),
58  edgeGlobals
59  )
60  );
61  }
62 
64  (
65  mesh(),
66  neiGlobal,
67  unionEqOp(),
68  false // apply separation
69  );
70 }
71 
72 
73 // Calculates per face the edge connected data (= cell or boundary in global
74 // numbering).
75 void Foam::FECCellToFaceStencil::calcFaceStencil
76 (
77  labelListList& faceStencil
78 ) const
79 {
80  const polyBoundaryMesh& patches = mesh().boundaryMesh();
81  const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
82  const labelList& own = mesh().faceOwner();
83  const labelList& nei = mesh().faceNeighbour();
84 
85 
86 
87  // Determine neighbouring global cell
88  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
89 
90  labelList neiGlobalCell(nBnd);
91  forAll(patches, patchI)
92  {
93  const polyPatch& pp = patches[patchI];
94 
95  if (pp.coupled())
96  {
97  label faceI = pp.start();
98 
99  forAll(pp, i)
100  {
101  neiGlobalCell[faceI-mesh().nInternalFaces()] =
102  globalNumbering().toGlobal(own[faceI]);
103  faceI++;
104  }
105  }
106  }
107  syncTools::swapBoundaryFaceList(mesh(), neiGlobalCell, false);
108 
109 
110 
111  // Determine on coupled edges the edge cells on the other side
112  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
113 
114  // Calculate edges on coupled patches
115  labelList boundaryEdges
116  (
117  allCoupledFacesPatch()().meshEdges
118  (
119  mesh().edges(),
120  mesh().pointEdges()
121  )
122  );
123 
124  // Mark boundary faces to be included in stencil (i.e. not coupled or empty)
125  boolList isValidBFace;
126  validBoundaryFaces(isValidBFace);
127 
128  // Swap edgeCells for coupled edges. Note: use EdgeMap for now since we've
129  // got syncTools::syncEdgeMap for those. Should be replaced with Map and
130  // syncTools functionality to handle those.
131  EdgeMap<labelList> neiGlobal;
132  calcEdgeBoundaryData
133  (
134  isValidBFace,
135  boundaryEdges,
136  neiGlobal
137  );
138 
139 
140 
141  // Construct stencil in global numbering
142  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
143 
144  faceStencil.setSize(mesh().nFaces());
145 
146  // Do coupled edges first
147 
148  forAll(boundaryEdges, i)
149  {
150  label edgeI = boundaryEdges[i];
151 
152  const labelList& eGlobals = neiGlobal[mesh().edges()[edgeI]];
153 
154  // Distribute to all edgeFaces
155  const labelList& eFaces = mesh().edgeFaces(edgeI);
156 
157  forAll(eFaces, j)
158  {
159  label faceI = eFaces[j];
160 
161  // Insert eGlobals into faceStencil.
162  merge(-1, -1, eGlobals, faceStencil[faceI]);
163  }
164  }
165  neiGlobal.clear();
166 
167 
168  // Do remaining edges by looping over all faces
169 
170  // Work arrays
171  DynamicList<label> fEdgesSet;
172  DynamicList<label> eFacesSet;
173  labelHashSet faceStencilSet;
174 
175  for (label faceI = 0; faceI < mesh().nInternalFaces(); faceI++)
176  {
177  label globalOwn = globalNumbering().toGlobal(own[faceI]);
178  label globalNei = globalNumbering().toGlobal(nei[faceI]);
179 
180  // Convert any existing faceStencil (from coupled edges) into
181  // set and operate on this.
182 
183  faceStencilSet.clear();
184 
185  // Insert all but global owner and neighbour
186  forAll(faceStencil[faceI], i)
187  {
188  label globalI = faceStencil[faceI][i];
189  if (globalI != globalOwn && globalI != globalNei)
190  {
191  faceStencilSet.insert(globalI);
192  }
193  }
194  faceStencil[faceI].clear();
195 
196  // Collect all edge connected (internal) cells
197  const labelList& fEdges = mesh().faceEdges(faceI, fEdgesSet);
198 
199  forAll(fEdges, i)
200  {
201  label edgeI = fEdges[i];
202 
203  insertFaceCells
204  (
205  globalOwn,
206  globalNei,
207  isValidBFace,
208  mesh().edgeFaces(edgeI, eFacesSet),
209  faceStencilSet
210  );
211  }
212 
213  // Extract, guarantee owner first, neighbour second.
214  faceStencil[faceI].setSize(faceStencilSet.size()+2);
215  label n = 0;
216  faceStencil[faceI][n++] = globalOwn;
217  faceStencil[faceI][n++] = globalNei;
218  forAllConstIter(labelHashSet, faceStencilSet, iter)
219  {
220  if (iter.key() == globalOwn || iter.key() == globalNei)
221  {
222  FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
223  << "problem:" << faceStencilSet
224  << abort(FatalError);
225  }
226  faceStencil[faceI][n++] = iter.key();
227  }
228  }
229  forAll(patches, patchI)
230  {
231  const polyPatch& pp = patches[patchI];
232  label faceI = pp.start();
233 
234  if (pp.coupled())
235  {
236  forAll(pp, i)
237  {
238  label globalOwn = globalNumbering().toGlobal(own[faceI]);
239  label globalNei = neiGlobalCell[faceI-mesh().nInternalFaces()];
240 
241  // Convert any existing faceStencil (from coupled edges) into
242  // set and operate on this.
243 
244  faceStencilSet.clear();
245 
246  // Insert all but global owner and neighbour
247  forAll(faceStencil[faceI], i)
248  {
249  label globalI = faceStencil[faceI][i];
250  if (globalI != globalOwn && globalI != globalNei)
251  {
252  faceStencilSet.insert(globalI);
253  }
254  }
255  faceStencil[faceI].clear();
256 
257  // Collect all edge connected (internal) cells
258  const labelList& fEdges = mesh().faceEdges(faceI, fEdgesSet);
259 
260  forAll(fEdges, i)
261  {
262  label edgeI = fEdges[i];
263 
264  insertFaceCells
265  (
266  globalOwn,
267  globalNei,
268  isValidBFace,
269  mesh().edgeFaces(edgeI, eFacesSet),
270  faceStencilSet
271  );
272  }
273 
274  // Extract, guarantee owner first, neighbour second.
275  faceStencil[faceI].setSize(faceStencilSet.size()+2);
276  label n = 0;
277  faceStencil[faceI][n++] = globalOwn;
278  faceStencil[faceI][n++] = globalNei;
279  forAllConstIter(labelHashSet, faceStencilSet, iter)
280  {
281  if (iter.key() == globalOwn || iter.key() == globalNei)
282  {
284  (
285  "FECCellToFaceStencil::calcFaceStencil(..)"
286  ) << "problem:" << faceStencilSet
287  << abort(FatalError);
288  }
289  faceStencil[faceI][n++] = iter.key();
290  }
291 
292  if (n != faceStencil[faceI].size())
293  {
294  FatalErrorIn("problem") << "n:" << n
295  << " size:" << faceStencil[faceI].size()
296  << abort(FatalError);
297  }
298 
299  faceI++;
300  }
301  }
302  else if (!isA<emptyPolyPatch>(pp))
303  {
304  forAll(pp, i)
305  {
306  label globalOwn = globalNumbering().toGlobal(own[faceI]);
307 
308  // Convert any existing faceStencil (from coupled edges) into
309  // set and operate on this.
310 
311  faceStencilSet.clear();
312 
313  // Insert all but global owner and neighbour
314  forAll(faceStencil[faceI], i)
315  {
316  label globalI = faceStencil[faceI][i];
317  if (globalI != globalOwn)
318  {
319  faceStencilSet.insert(globalI);
320  }
321  }
322  faceStencil[faceI].clear();
323 
324  // Collect all edge connected (internal) cells
325  const labelList& fEdges = mesh().faceEdges(faceI, fEdgesSet);
326 
327  forAll(fEdges, i)
328  {
329  label edgeI = fEdges[i];
330 
331  insertFaceCells
332  (
333  globalOwn,
334  -1,
335  isValidBFace,
336  mesh().edgeFaces(edgeI, eFacesSet),
337  faceStencilSet
338  );
339  }
340 
341  // Extract, guarantee owner first, neighbour second.
342  faceStencil[faceI].setSize(faceStencilSet.size()+1);
343  label n = 0;
344  faceStencil[faceI][n++] = globalOwn;
345  forAllConstIter(labelHashSet, faceStencilSet, iter)
346  {
347  if (iter.key() == globalOwn)
348  {
350  (
351  "FECCellToFaceStencil::calcFaceStencil(..)"
352  ) << "problem:" << faceStencilSet
353  << abort(FatalError);
354  }
355  faceStencil[faceI][n++] = iter.key();
356  }
357 
358  faceI++;
359  }
360  }
361  }
362 
363 
364  for (label faceI = 0; faceI < mesh().nInternalFaces(); faceI++)
365  {
366  label globalOwn = globalNumbering().toGlobal(own[faceI]);
367  if (faceStencil[faceI][0] != globalOwn)
368  {
369  FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
370  << "problem:" << faceStencil[faceI]
371  << " globalOwn:" << globalOwn
372  << abort(FatalError);
373  }
374  label globalNei = globalNumbering().toGlobal(nei[faceI]);
375  if (faceStencil[faceI][1] != globalNei)
376  {
377  FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
378  << "problem:" << faceStencil[faceI]
379  << " globalNei:" << globalNei
380  << abort(FatalError);
381  }
382  }
383 
384 
385  forAll(patches, patchI)
386  {
387  const polyPatch& pp = patches[patchI];
388 
389  if (pp.coupled())
390  {
391  forAll(pp, i)
392  {
393  label faceI = pp.start()+i;
394 
395  label globalOwn = globalNumbering().toGlobal(own[faceI]);
396  if (faceStencil[faceI][0] != globalOwn)
397  {
398  FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
399  << "problem:" << faceStencil[faceI]
400  << " globalOwn:" << globalOwn
401  << abort(FatalError);
402  }
403  label globalNei = neiGlobalCell[faceI-mesh().nInternalFaces()];
404  if (faceStencil[faceI][1] != globalNei)
405  {
406  FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
407  << "problem:" << faceStencil[faceI]
408  << " globalNei:" << globalNei
409  << abort(FatalError);
410  }
411  }
412  }
413  else if (!isA<emptyPolyPatch>(pp))
414  {
415  forAll(pp, i)
416  {
417  label faceI = pp.start()+i;
418 
419  label globalOwn = globalNumbering().toGlobal(own[faceI]);
420  if (faceStencil[faceI][0] != globalOwn)
421  {
422  FatalErrorIn("FECCellToFaceStencil::calcFaceStencil(..)")
423  << "problem:" << faceStencil[faceI]
424  << " globalOwn:" << globalOwn
425  << abort(FatalError);
426  }
427  }
428  }
429  }
430 }
431 
432 
433 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
434 
436 :
437  cellToFaceStencil(mesh)
438 {
439  // Calculate per face the (edge) connected cells (in global numbering)
440  labelListList faceStencil;
441  calcFaceStencil(faceStencil);
442 
443  // Transfer to *this
444  transfer(faceStencil);
445 }
446 
447 
448 // ************************ vim: set sw=4 sts=4 et: ************************ //