All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
extendedUpwindCellToFaceStencil.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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  << "problem:" << faceStencilSet
173  << abort(FatalError);
174  }
175  }
176  else
177  {
178  transportedStencil.setSize(faceStencilSet.size()+1);
179  label n = 0;
180  transportedStencil[n++] = globalOwn;
181 
182  forAllConstIter(labelHashSet, faceStencilSet, iter)
183  {
184  if (iter.key() != globalOwn)
185  {
186  transportedStencil[n++] = iter.key();
187  }
188  }
189  if (n != transportedStencil.size())
190  {
192  << "problem:" << faceStencilSet
193  << abort(FatalError);
194  }
195  }
196 }
197 
198 
199 void Foam::extendedUpwindCellToFaceStencil::transportStencils
200 (
201  const labelListList& faceStencil,
202  const scalar minOpposedness,
203  labelListList& ownStencil,
204  labelListList& neiStencil
205 )
206 {
207  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
208  const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
209  const labelList& own = mesh_.faceOwner();
210  const labelList& nei = mesh_.faceNeighbour();
211 
212  // Work arrays
213  DynamicList<label> oppositeFaces;
214  labelHashSet faceStencilSet;
215 
216 
217  // For quick detection of empty faces
218  boolList nonEmptyFace(mesh_.nFaces(), true);
219  forAll(patches, patchi)
220  {
221  const polyPatch& pp = patches[patchi];
222 
223  if (isA<emptyPolyPatch>(pp))
224  {
225  label facei = pp.start();
226  forAll(pp, i)
227  {
228  nonEmptyFace[facei++] = false;
229  }
230  }
231  }
232 
233 
234  // Do the owner side
235  // ~~~~~~~~~~~~~~~~~
236  // stencil is synchronised at entry so no need to swap.
237 
238  ownStencil.setSize(mesh_.nFaces());
239 
240  // Internal faces
241  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
242  {
243  // Get stencil as owner + neighbour + stencil from 'opposite' faces
244  transportStencil
245  (
246  nonEmptyFace,
247  faceStencil,
248  minOpposedness,
249  facei,
250  own[facei],
251  true, // stencilHasNeighbour
252  oppositeFaces,
253  faceStencilSet,
254  ownStencil[facei]
255  );
256  }
257  // Boundary faces
258  forAll(patches, patchi)
259  {
260  const polyPatch& pp = patches[patchi];
261  label facei = pp.start();
262 
263  if (pp.coupled())
264  {
265  forAll(pp, i)
266  {
267  transportStencil
268  (
269  nonEmptyFace,
270  faceStencil,
271  minOpposedness,
272  facei,
273  own[facei],
274  true, // stencilHasNeighbour
275 
276  oppositeFaces,
277  faceStencilSet,
278  ownStencil[facei]
279  );
280  facei++;
281  }
282  }
283  else if (!isA<emptyPolyPatch>(pp))
284  {
285  forAll(pp, i)
286  {
287  // faceStencil does not contain neighbour
288  transportStencil
289  (
290  nonEmptyFace,
291  faceStencil,
292  minOpposedness,
293  facei,
294  own[facei],
295  false, // stencilHasNeighbour
296 
297  oppositeFaces,
298  faceStencilSet,
299  ownStencil[facei]
300  );
301  facei++;
302  }
303  }
304  }
305 
306 
307  // Swap coupled boundary stencil
308  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
309 
310  labelListList neiBndStencil(nBnd);
311  for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
312  {
313  neiBndStencil[facei-mesh_.nInternalFaces()] = ownStencil[facei];
314  }
315  // syncTools::swapBoundaryFaceList(mesh_, neiBndStencil);
317  (
318  mesh_,
319  neiBndStencil,
320  eqOp<labelList>(),
321  dummyTransform()
322  );
323 
324 
325 
326  // Do the neighbour side
327  // ~~~~~~~~~~~~~~~~~~~~~
328  // - internal faces : get opposite faces on neighbour side
329  // - boundary faces : empty
330  // - coupled faces : in neiBndStencil
331 
332  neiStencil.setSize(mesh_.nFaces());
333 
334  // Internal faces
335  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
336  {
337  transportStencil
338  (
339  nonEmptyFace,
340  faceStencil,
341  minOpposedness,
342  facei,
343  nei[facei],
344  true, // stencilHasNeighbour
345 
346  oppositeFaces,
347  faceStencilSet,
348  neiStencil[facei]
349  );
350  }
351 
352  // Boundary faces
353  forAll(patches, patchi)
354  {
355  const polyPatch& pp = patches[patchi];
356  label facei = pp.start();
357 
358  if (pp.coupled())
359  {
360  forAll(pp, i)
361  {
362  neiStencil[facei].transfer
363  (
364  neiBndStencil[facei-mesh_.nInternalFaces()]
365  );
366  facei++;
367  }
368  }
369  else
370  {
371  // Boundary has empty neighbour stencil
372  }
373  }
374 }
375 
376 
377 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
378 
380 (
381  const cellToFaceStencil& stencil,
382  const bool pureUpwind,
383  const scalar minOpposedness
384 )
385 :
386  extendedCellToFaceStencil(stencil.mesh()),
387  pureUpwind_(pureUpwind)
388 {
389  // forAll(stencil, facei)
390  //{
391  // const labelList& fCells = stencil[facei];
392  //
393  // Pout<< "Face:" << facei << " at:" << mesh_.faceCentres()[facei]
394  // << endl;
395  //
396  // forAll(fCells, i)
397  // {
398  // label globalI = fCells[i];
399  //
400  // if (globalI < mesh_.nCells())
401  // {
402  // Pout<< " cell:" << globalI
403  // << " at:" << mesh_.cellCentres()[globalI] << endl;
404  // }
405  // else
406  // {
407  // label facei = globalI-mesh_.nCells() + mesh_.nInternalFaces();
408  //
409  // Pout<< " boundary:" << facei
410  // << " at:" << mesh_.faceCentres()[facei] << endl;
411  // }
412  // }
413  //}
414  // Pout<< endl << endl;
415 
416 
417  // Transport centred stencil to upwind/downwind face
418  transportStencils
419  (
420  stencil,
421  minOpposedness,
422  ownStencil_,
423  neiStencil_
424  );
425 
426  {
427  List<Map<label>> compactMap(Pstream::nProcs());
428  ownMapPtr_.reset
429  (
430  new mapDistribute
431  (
432  stencil.globalNumbering(),
433  ownStencil_,
434  compactMap
435  )
436  );
437  }
438 
439  {
440 
441  List<Map<label>> compactMap(Pstream::nProcs());
442  neiMapPtr_.reset
443  (
444  new mapDistribute
445  (
446  stencil.globalNumbering(),
447  neiStencil_,
448  compactMap
449  )
450  );
451  }
452 
453  // stencil now in compact form
454  if (pureUpwind_)
455  {
456  const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh());
457 
458  List<List<point>> stencilPoints(ownStencil_.size());
459 
460  // Owner stencil
461  // ~~~~~~~~~~~~~
462 
463  collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints);
464 
465  // Mask off all stencil points on wrong side of face
466  forAll(stencilPoints, facei)
467  {
468  const point& fc = mesh.faceCentres()[facei];
469  const vector& fArea = mesh.faceAreas()[facei];
470 
471  const List<point>& points = stencilPoints[facei];
472  const labelList& stencil = ownStencil_[facei];
473 
474  DynamicList<label> newStencil(stencil.size());
475  forAll(points, i)
476  {
477  if (((points[i]-fc) & fArea) < 0)
478  {
479  newStencil.append(stencil[i]);
480  }
481  }
482  if (newStencil.size() != stencil.size())
483  {
484  ownStencil_[facei].transfer(newStencil);
485  }
486  }
487 
488 
489  // Neighbour stencil
490  // ~~~~~~~~~~~~~~~~~
491 
492  collectData(neiMapPtr_(), neiStencil_, mesh.C(), stencilPoints);
493 
494  // Mask off all stencil points on wrong side of face
495  forAll(stencilPoints, facei)
496  {
497  const point& fc = mesh.faceCentres()[facei];
498  const vector& fArea = mesh.faceAreas()[facei];
499 
500  const List<point>& points = stencilPoints[facei];
501  const labelList& stencil = neiStencil_[facei];
502 
503  DynamicList<label> newStencil(stencil.size());
504  forAll(points, i)
505  {
506  if (((points[i]-fc) & fArea) > 0)
507  {
508  newStencil.append(stencil[i]);
509  }
510  }
511  if (newStencil.size() != stencil.size())
512  {
513  neiStencil_[facei].transfer(newStencil);
514  }
515  }
516 
517  // Note: could compact schedule as well. for if cells are not needed
518  // across any boundary anymore. However relatively rare.
519  }
520 }
521 
522 
524 (
525  const cellToFaceStencil& stencil
526 )
527 :
528  extendedCellToFaceStencil(stencil.mesh()),
529  pureUpwind_(true)
530 {
531  // Calculate stencil points with full stencil
532 
533  ownStencil_ = stencil;
534 
535  {
536  List<Map<label>> compactMap(Pstream::nProcs());
537  ownMapPtr_.reset
538  (
539  new mapDistribute
540  (
541  stencil.globalNumbering(),
542  ownStencil_,
543  compactMap
544  )
545  );
546  }
547 
548  const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh());
549 
550  List<List<point>> stencilPoints(ownStencil_.size());
551  collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints);
552 
553  // Split stencil into owner and neighbour
554  neiStencil_.setSize(ownStencil_.size());
555 
556  forAll(stencilPoints, facei)
557  {
558  const point& fc = mesh.faceCentres()[facei];
559  const vector& fArea = mesh.faceAreas()[facei];
560 
561  const List<point>& points = stencilPoints[facei];
562  const labelList& stencil = ownStencil_[facei];
563 
564  DynamicList<label> newOwnStencil(stencil.size());
565  DynamicList<label> newNeiStencil(stencil.size());
566  forAll(points, i)
567  {
568  if (((points[i]-fc) & fArea) > 0)
569  {
570  newNeiStencil.append(stencil[i]);
571  }
572  else
573  {
574  newOwnStencil.append(stencil[i]);
575  }
576  }
577  if (newNeiStencil.size() > 0)
578  {
579  ownStencil_[facei].transfer(newOwnStencil);
580  neiStencil_[facei].transfer(newNeiStencil);
581  }
582  }
583 
584  // Should compact schedule. Or have both return the same schedule.
585  neiMapPtr_.reset(new mapDistribute(ownMapPtr_()));
586 }
587 
588 
589 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
static void syncBoundaryFaceList(const polyMesh &, UList< T > &, const CombineOp &cop, const TransformOp &top, const bool parRun=Pstream::parRun())
Synchronise values on boundary faces only.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1175
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
label nFaces() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const cellList & cells() const
Base class for extended cell-to-face stencils (face values from neighbouring cells) ...
const globalIndex & globalNumbering() const
Global numbering for cells and boundary faces.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
Dummy transform to be used with syncTools.
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
extendedUpwindCellToFaceStencil(const cellToFaceStencil &, const bool pureUpwind, const scalar minOpposedness)
Construct from mesh and uncompacted centred face stencil.
dynamicFvMesh & mesh
const pointField & points
extendedCellToFaceStencil(const polyMesh &)
Construct from mesh.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
List< label > labelList
A List of labels.
Definition: labelList.H:56
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
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.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const polyMesh & mesh() const
const vectorField & faceCentres() const
void setSize(const label)
Reset size of List.
Definition: List.C:281
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
label patchi
Class containing processor-to-processor mapping information.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
const vectorField & faceAreas() const
Field< vector > vectorField
Specialisation of Field<T> for vector.
const volVectorField & C() const
Return cell centres as volVectorField.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342