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