FreeFOAM The Cross-Platform CFD Toolkit
faceMDLimitedGrads.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 "faceMDLimitedGrad.H"
28 #include <finiteVolume/gaussGrad.H>
29 #include <finiteVolume/fvMesh.H>
30 #include <finiteVolume/volMesh.H>
32 #include <finiteVolume/volFields.H>
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 
42 namespace fv
43 {
44 
45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46 
48 
49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50 
51 // FaceLimited scalar gradient
52 template<>
54 (
55  const volScalarField& vsf
56 ) const
57 {
58  const fvMesh& mesh = vsf.mesh();
59 
60  tmp<volVectorField> tGrad = basicGradScheme_().grad(vsf);
61 
62  if (k_ < SMALL)
63  {
64  return tGrad;
65  }
66 
67  volVectorField& g = tGrad();
68 
69  const unallocLabelList& owner = mesh.owner();
70  const unallocLabelList& neighbour = mesh.neighbour();
71 
72  const volVectorField& C = mesh.C();
73  const surfaceVectorField& Cf = mesh.Cf();
74 
75  scalar rk = (1.0/k_ - 1.0);
76 
77  forAll(owner, facei)
78  {
79  label own = owner[facei];
80  label nei = neighbour[facei];
81 
82  scalar vsfOwn = vsf[own];
83  scalar vsfNei = vsf[nei];
84 
85  scalar maxFace = max(vsfOwn, vsfNei);
86  scalar minFace = min(vsfOwn, vsfNei);
87 
88  if (k_ < 1.0)
89  {
90  scalar maxMinFace = rk*(maxFace - minFace);
91  maxFace += maxMinFace;
92  minFace -= maxMinFace;
93  }
94 
95  // owner side
97  (
98  g[own],
99  maxFace - vsfOwn,
100  minFace - vsfOwn,
101  Cf[facei] - C[own]
102  );
103 
104  // neighbour side
106  (
107  g[nei],
108  maxFace - vsfNei,
109  minFace - vsfNei,
110  Cf[facei] - C[nei]
111  );
112  }
113 
114  const volScalarField::GeometricBoundaryField& bsf = vsf.boundaryField();
115 
116  forAll(bsf, patchi)
117  {
118  const fvPatchScalarField& psf = bsf[patchi];
119 
120  const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
121  const vectorField& pCf = Cf.boundaryField()[patchi];
122 
123  if (psf.coupled())
124  {
125  scalarField psfNei = psf.patchNeighbourField();
126 
127  forAll(pOwner, pFacei)
128  {
129  label own = pOwner[pFacei];
130 
131  scalar vsfOwn = vsf[own];
132  scalar vsfNei = psfNei[pFacei];
133 
134  scalar maxFace = max(vsfOwn, vsfNei);
135  scalar minFace = min(vsfOwn, vsfNei);
136 
137  if (k_ < 1.0)
138  {
139  scalar maxMinFace = rk*(maxFace - minFace);
140  maxFace += maxMinFace;
141  minFace -= maxMinFace;
142  }
143 
145  (
146  g[own],
147  maxFace - vsfOwn,
148  minFace - vsfOwn,
149  pCf[pFacei] - C[own]
150  );
151  }
152  }
153  else if (psf.fixesValue())
154  {
155  forAll(pOwner, pFacei)
156  {
157  label own = pOwner[pFacei];
158 
159  scalar vsfOwn = vsf[own];
160  scalar vsfNei = psf[pFacei];
161 
162  scalar maxFace = max(vsfOwn, vsfNei);
163  scalar minFace = min(vsfOwn, vsfNei);
164 
165  if (k_ < 1.0)
166  {
167  scalar maxMinFace = rk*(maxFace - minFace);
168  maxFace += maxMinFace;
169  minFace -= maxMinFace;
170  }
171 
173  (
174  g[own],
175  maxFace - vsfOwn,
176  minFace - vsfOwn,
177  pCf[pFacei] - C[own]
178  );
179  }
180  }
181  }
182 
185 
186  return tGrad;
187 }
188 
189 
190 template<>
192 (
193  const volVectorField& vvf
194 ) const
195 {
196  const fvMesh& mesh = vvf.mesh();
197 
198  tmp<volTensorField> tGrad = basicGradScheme_().grad(vvf);
199 
200  if (k_ < SMALL)
201  {
202  return tGrad;
203  }
204 
205  volTensorField& g = tGrad();
206 
207  const unallocLabelList& owner = mesh.owner();
208  const unallocLabelList& neighbour = mesh.neighbour();
209 
210  const volVectorField& C = mesh.C();
211  const surfaceVectorField& Cf = mesh.Cf();
212 
213  scalar rk = (1.0/k_ - 1.0);
214 
215  forAll(owner, facei)
216  {
217  label own = owner[facei];
218  label nei = neighbour[facei];
219 
220  vector vvfOwn = vvf[own];
221  vector vvfNei = vvf[nei];
222 
223  vector maxFace = max(vvfOwn, vvfNei);
224  vector minFace = min(vvfOwn, vvfNei);
225 
226  if (k_ < 1.0)
227  {
228  vector maxMinFace = rk*(maxFace - minFace);
229  maxFace += maxMinFace;
230  minFace -= maxMinFace;
231  }
232 
233  // owner side
235  (
236  g[own],
237  maxFace - vvfOwn,
238  minFace - vvfOwn,
239  Cf[facei] - C[nei]
240  );
241 
242 
243  // neighbour side
245  (
246  g[nei],
247  maxFace - vvfNei,
248  minFace - vvfNei,
249  Cf[facei] - C[nei]
250  );
251  }
252 
253 
254  const volVectorField::GeometricBoundaryField& bvf = vvf.boundaryField();
255 
256  forAll(bvf, patchi)
257  {
258  const fvPatchVectorField& psf = bvf[patchi];
259 
260  const unallocLabelList& pOwner = mesh.boundary()[patchi].faceCells();
261  const vectorField& pCf = Cf.boundaryField()[patchi];
262 
263  if (psf.coupled())
264  {
265  vectorField psfNei = psf.patchNeighbourField();
266 
267  forAll(pOwner, pFacei)
268  {
269  label own = pOwner[pFacei];
270 
271  vector vvfOwn = vvf[own];
272  vector vvfNei = psfNei[pFacei];
273 
274  vector maxFace = max(vvfOwn, vvfNei);
275  vector minFace = min(vvfOwn, vvfNei);
276 
277  if (k_ < 1.0)
278  {
279  vector maxMinFace = rk*(maxFace - minFace);
280  maxFace += maxMinFace;
281  minFace -= maxMinFace;
282  }
283 
285  (
286  g[own],
287  maxFace - vvfOwn, minFace - vvfOwn,
288  pCf[pFacei] - C[own]
289  );
290  }
291  }
292  else if (psf.fixesValue())
293  {
294  forAll(pOwner, pFacei)
295  {
296  label own = pOwner[pFacei];
297 
298  vector vvfOwn = vvf[own];
299  vector vvfNei = psf[pFacei];
300 
301  vector maxFace = max(vvfOwn, vvfNei);
302  vector minFace = min(vvfOwn, vvfNei);
303 
304  if (k_ < 1.0)
305  {
306  vector maxMinFace = rk*(maxFace - minFace);
307  maxFace += maxMinFace;
308  minFace -= maxMinFace;
309  }
310 
312  (
313  g[own],
314  maxFace - vvfOwn,
315  minFace - vvfOwn,
316  pCf[pFacei] - C[own]
317  );
318  }
319  }
320  }
321 
324 
325  return tGrad;
326 }
327 
328 
329 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
330 
331 } // End namespace fv
332 
333 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
334 
335 } // End namespace Foam
336 
337 // ************************ vim: set sw=4 sts=4 et: ************************ //