FreeFOAM The Cross-Platform CFD Toolkit
extendedUpwindCellToFaceStencil.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 
28 #include <OpenFOAM/syncTools.H>
29 #include <OpenFOAM/SortableList.H>
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 void Foam::extendedUpwindCellToFaceStencil::selectOppositeFaces
34 (
35  const boolList& nonEmptyFace,
36  const scalar minOpposedness,
37  const label faceI,
38  const label cellI,
39  DynamicList<label>& oppositeFaces
40 ) const
41 {
42  const vectorField& areas = mesh_.faceAreas();
43  const labelList& own = mesh_.faceOwner();
44  const cell& cFaces = mesh_.cells()[cellI];
45 
46  SortableList<scalar> opposedness(cFaces.size(), -GREAT);
47 
48  // Pick up all the faces that oppose this one.
49  forAll(cFaces, i)
50  {
51  label otherFaceI = cFaces[i];
52 
53  if (otherFaceI != faceI && nonEmptyFace[otherFaceI])
54  {
55  if ((own[otherFaceI] == cellI) == (own[faceI] == cellI))
56  {
57  opposedness[i] = -(areas[otherFaceI] & areas[faceI]);
58  }
59  else
60  {
61  opposedness[i] = (areas[otherFaceI] & areas[faceI]);
62  }
63  }
64  }
65 
66  label sz = opposedness.size();
67 
68  oppositeFaces.clear();
69 
70  scalar myAreaSqr = magSqr(areas[faceI]);
71 
72  if (myAreaSqr > VSMALL)
73  {
74  forAll(opposedness, i)
75  {
76  opposedness[i] /= myAreaSqr;
77  }
78  // Sort in incrementing order
79  opposedness.sort();
80 
81  // Pick largest no matter what
82  oppositeFaces.append(cFaces[opposedness.indices()[sz-1]]);
83 
84  for (label i = sz-2; i >= 0; --i)
85  {
86  if (opposedness[i] < minOpposedness)
87  {
88  break;
89  }
90  oppositeFaces.append(cFaces[opposedness.indices()[i]]);
91  }
92  }
93  else
94  {
95  // Sort in incrementing order
96  opposedness.sort();
97 
98  // Tiny face. Do what?
99  // Pick largest no matter what
100  oppositeFaces.append(cFaces[opposedness.indices()[sz-1]]);
101  }
102 }
103 
104 
105 void Foam::extendedUpwindCellToFaceStencil::transportStencil
106 (
107  const boolList& nonEmptyFace,
108  const labelListList& faceStencil,
109  const scalar minOpposedness,
110  const label faceI,
111  const label cellI,
112  const bool stencilHasNeighbour,
113 
114  DynamicList<label>& oppositeFaces,
115  labelHashSet& faceStencilSet,
116  labelList& transportedStencil
117 ) const
118 {
119  label globalOwn = faceStencil[faceI][0];
120  label globalNei = -1;
121  if (stencilHasNeighbour && faceStencil[faceI].size() >= 2)
122  {
123  globalNei = faceStencil[faceI][1];
124  }
125 
126 
127  selectOppositeFaces
128  (
129  nonEmptyFace,
130  minOpposedness,
131  faceI,
132  cellI,
133  oppositeFaces
134  );
135 
136  // Collect all stencils of oppositefaces
137  faceStencilSet.clear();
138  forAll(oppositeFaces, i)
139  {
140  const labelList& fStencil = faceStencil[oppositeFaces[i]];
141 
142  forAll(fStencil, j)
143  {
144  label globalI = fStencil[j];
145 
146  if (globalI != globalOwn && globalI != globalNei)
147  {
148  faceStencilSet.insert(globalI);
149  }
150  }
151  }
152 
153  // Add my owner and neighbour first.
154  if (stencilHasNeighbour)
155  {
156  transportedStencil.setSize(faceStencilSet.size()+2);
157  label n = 0;
158  transportedStencil[n++] = globalOwn;
159  transportedStencil[n++] = globalNei;
160 
161  forAllConstIter(labelHashSet, faceStencilSet, iter)
162  {
163  if (iter.key() != globalOwn && iter.key() != globalNei)
164  {
165  transportedStencil[n++] = iter.key();
166  }
167  }
168  if (n != transportedStencil.size())
169  {
171  (
172  "extendedUpwindCellToFaceStencil::transportStencil(..)"
173  ) << "problem:" << faceStencilSet
174  << abort(FatalError);
175  }
176  }
177  else
178  {
179  transportedStencil.setSize(faceStencilSet.size()+1);
180  label n = 0;
181  transportedStencil[n++] = globalOwn;
182 
183  forAllConstIter(labelHashSet, faceStencilSet, iter)
184  {
185  if (iter.key() != globalOwn)
186  {
187  transportedStencil[n++] = iter.key();
188  }
189  }
190  if (n != transportedStencil.size())
191  {
193  (
194  "extendedUpwindCellToFaceStencil::transportStencil(..)"
195  ) << "problem:" << faceStencilSet
196  << abort(FatalError);
197  }
198  }
199 }
200 
201 
202 void Foam::extendedUpwindCellToFaceStencil::transportStencils
203 (
204  const labelListList& faceStencil,
205  const scalar minOpposedness,
206  labelListList& ownStencil,
207  labelListList& neiStencil
208 )
209 {
210  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
211  const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
212  const labelList& own = mesh_.faceOwner();
213  const labelList& nei = mesh_.faceNeighbour();
214 
215  // Work arrays
216  DynamicList<label> oppositeFaces;
217  labelHashSet faceStencilSet;
218 
219 
220  // For quick detection of empty faces
221  boolList nonEmptyFace(mesh_.nFaces(), true);
222  forAll(patches, patchI)
223  {
224  const polyPatch& pp = patches[patchI];
225 
226  if (isA<emptyPolyPatch>(pp))
227  {
228  label faceI = pp.start();
229  forAll(pp, i)
230  {
231  nonEmptyFace[faceI++] = false;
232  }
233  }
234  }
235 
236 
237  // Do the owner side
238  // ~~~~~~~~~~~~~~~~~
239  // stencil is synchronised at entry so no need to swap.
240 
241  ownStencil.setSize(mesh_.nFaces());
242 
243  // Internal faces
244  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
245  {
246  // Get stencil as owner + neighbour + stencil from 'opposite' faces
247  transportStencil
248  (
249  nonEmptyFace,
250  faceStencil,
251  minOpposedness,
252  faceI,
253  own[faceI],
254  true, //stencilHasNeighbour
255  oppositeFaces,
256  faceStencilSet,
257  ownStencil[faceI]
258  );
259  }
260  // Boundary faces
261  forAll(patches, patchI)
262  {
263  const polyPatch& pp = patches[patchI];
264  label faceI = pp.start();
265 
266  if (pp.coupled())
267  {
268  forAll(pp, i)
269  {
270  transportStencil
271  (
272  nonEmptyFace,
273  faceStencil,
274  minOpposedness,
275  faceI,
276  own[faceI],
277  true, //stencilHasNeighbour
278 
279  oppositeFaces,
280  faceStencilSet,
281  ownStencil[faceI]
282  );
283  faceI++;
284  }
285  }
286  else if (!isA<emptyPolyPatch>(pp))
287  {
288  forAll(pp, i)
289  {
290  // faceStencil does not contain neighbour
291  transportStencil
292  (
293  nonEmptyFace,
294  faceStencil,
295  minOpposedness,
296  faceI,
297  own[faceI],
298  false, //stencilHasNeighbour
299 
300  oppositeFaces,
301  faceStencilSet,
302  ownStencil[faceI]
303  );
304  faceI++;
305  }
306  }
307  }
308 
309 
310  // Swap coupled boundary stencil
311  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
312 
313  labelListList neiBndStencil(nBnd);
314  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
315  {
316  neiBndStencil[faceI-mesh_.nInternalFaces()] = ownStencil[faceI];
317  }
318  syncTools::swapBoundaryFaceList(mesh_, neiBndStencil, false);
319 
320 
321  // Do the neighbour side
322  // ~~~~~~~~~~~~~~~~~~~~~
323  // - internal faces : get opposite faces on neighbour side
324  // - boundary faces : empty
325  // - coupled faces : in neiBndStencil
326 
327  neiStencil.setSize(mesh_.nFaces());
328 
329  // Internal faces
330  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
331  {
332  transportStencil
333  (
334  nonEmptyFace,
335  faceStencil,
336  minOpposedness,
337  faceI,
338  nei[faceI],
339  true, //stencilHasNeighbour
340 
341  oppositeFaces,
342  faceStencilSet,
343  neiStencil[faceI]
344  );
345  }
346 
347  // Boundary faces
348  forAll(patches, patchI)
349  {
350  const polyPatch& pp = patches[patchI];
351  label faceI = pp.start();
352 
353  if (pp.coupled())
354  {
355  forAll(pp, i)
356  {
357  neiStencil[faceI].transfer
358  (
359  neiBndStencil[faceI-mesh_.nInternalFaces()]
360  );
361  faceI++;
362  }
363  }
364  else
365  {
366  // Boundary has empty neighbour stencil
367  }
368  }
369 }
370 
371 
372 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
373 
374 Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil
375 (
376  const cellToFaceStencil& stencil,
377  const bool pureUpwind,
378  const scalar minOpposedness
379 )
380 :
381  extendedCellToFaceStencil(stencil.mesh()),
382  pureUpwind_(pureUpwind)
383 {
384  //forAll(stencil, faceI)
385  //{
386  // const labelList& fCells = stencil[faceI];
387  //
388  // Pout<< "Face:" << faceI << " at:" << mesh_.faceCentres()[faceI]
389  // << endl;
390  //
391  // forAll(fCells, i)
392  // {
393  // label globalI = fCells[i];
394  //
395  // if (globalI < mesh_.nCells())
396  // {
397  // Pout<< " cell:" << globalI
398  // << " at:" << mesh_.cellCentres()[globalI] << endl;
399  // }
400  // else
401  // {
402  // label faceI = globalI-mesh_.nCells() + mesh_.nInternalFaces();
403  //
404  // Pout<< " boundary:" << faceI
405  // << " at:" << mesh_.faceCentres()[faceI] << endl;
406  // }
407  // }
408  //}
409  //Pout<< endl << endl;
410 
411 
412  // Transport centred stencil to upwind/downwind face
413  transportStencils
414  (
415  stencil,
416  minOpposedness,
417  ownStencil_,
418  neiStencil_
419  );
420 
421  ownMapPtr_ = calcDistributeMap
422  (
423  stencil.mesh(),
424  stencil.globalNumbering(),
425  ownStencil_
426  );
427 
428  neiMapPtr_ = calcDistributeMap
429  (
430  stencil.mesh(),
431  stencil.globalNumbering(),
432  neiStencil_
433  );
434 
435  // stencil now in compact form
436  if (pureUpwind_)
437  {
438  const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh());
439 
440  List<List<point> > stencilPoints(ownStencil_.size());
441 
442  // Owner stencil
443  // ~~~~~~~~~~~~~
444 
445  collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints);
446 
447  // Mask off all stencil points on wrong side of face
448  forAll(stencilPoints, faceI)
449  {
450  const point& fc = mesh.faceCentres()[faceI];
451  const vector& fArea = mesh.faceAreas()[faceI];
452 
453  const List<point>& points = stencilPoints[faceI];
454  const labelList& stencil = ownStencil_[faceI];
455 
456  DynamicList<label> newStencil(stencil.size());
457  forAll(points, i)
458  {
459  if (((points[i]-fc) & fArea) < 0)
460  {
461  newStencil.append(stencil[i]);
462  }
463  }
464  if (newStencil.size() != stencil.size())
465  {
466  ownStencil_[faceI].transfer(newStencil);
467  }
468  }
469 
470 
471  // Neighbour stencil
472  // ~~~~~~~~~~~~~~~~~
473 
474  collectData(neiMapPtr_(), neiStencil_, mesh.C(), stencilPoints);
475 
476  // Mask off all stencil points on wrong side of face
477  forAll(stencilPoints, faceI)
478  {
479  const point& fc = mesh.faceCentres()[faceI];
480  const vector& fArea = mesh.faceAreas()[faceI];
481 
482  const List<point>& points = stencilPoints[faceI];
483  const labelList& stencil = neiStencil_[faceI];
484 
485  DynamicList<label> newStencil(stencil.size());
486  forAll(points, i)
487  {
488  if (((points[i]-fc) & fArea) > 0)
489  {
490  newStencil.append(stencil[i]);
491  }
492  }
493  if (newStencil.size() != stencil.size())
494  {
495  neiStencil_[faceI].transfer(newStencil);
496  }
497  }
498 
499  // Note: could compact schedule as well. for if cells are not needed
500  // across any boundary anymore. However relatively rare.
501  }
502 }
503 
504 
505 Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil
506 (
507  const cellToFaceStencil& stencil
508 )
509 :
510  extendedCellToFaceStencil(stencil.mesh()),
511  pureUpwind_(true)
512 {
513  // Calculate stencil points with full stencil
514 
515  ownStencil_ = stencil;
516 
517  ownMapPtr_ = calcDistributeMap
518  (
519  stencil.mesh(),
520  stencil.globalNumbering(),
521  ownStencil_
522  );
523 
524  const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh());
525 
526  List<List<point> > stencilPoints(ownStencil_.size());
527  collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints);
528 
529  // Split stencil into owner and neighbour
530  neiStencil_.setSize(ownStencil_.size());
531 
532  forAll(stencilPoints, faceI)
533  {
534  const point& fc = mesh.faceCentres()[faceI];
535  const vector& fArea = mesh.faceAreas()[faceI];
536 
537  const List<point>& points = stencilPoints[faceI];
538  const labelList& stencil = ownStencil_[faceI];
539 
540  DynamicList<label> newOwnStencil(stencil.size());
541  DynamicList<label> newNeiStencil(stencil.size());
542  forAll(points, i)
543  {
544  if (((points[i]-fc) & fArea) > 0)
545  {
546  newNeiStencil.append(stencil[i]);
547  }
548  else
549  {
550  newOwnStencil.append(stencil[i]);
551  }
552  }
553  if (newNeiStencil.size() > 0)
554  {
555  ownStencil_[faceI].transfer(newOwnStencil);
556  neiStencil_[faceI].transfer(newNeiStencil);
557  }
558  }
559 
560  // Should compact schedule. Or have both return the same schedule.
561  neiMapPtr_.reset(new mapDistribute(ownMapPtr_()));
562 }
563 
564 
565 // ************************ vim: set sw=4 sts=4 et: ************************ //