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