FreeFOAM The Cross-Platform CFD Toolkit
fourthGrad.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 "fourthGrad.H"
28 #include <finiteVolume/gaussGrad.H>
29 #include <finiteVolume/fvMesh.H>
30 #include <finiteVolume/volMesh.H>
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 
42 namespace fv
43 {
44 
45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
46 
47 template<class Type>
48 tmp
49 <
50  GeometricField
51  <
52  typename outerProduct<vector, Type>::type, fvPatchField, volMesh
53  >
54 >
56 (
57  const GeometricField<Type, fvPatchField, volMesh>& vsf
58 ) const
59 {
60  // The fourth-order gradient is calculated in two passes. First,
61  // the standard least-square gradient is assembled. Then, the
62  // fourth-order correction is added to the second-order accurate
63  // gradient to complete the accuracy.
64 
65  typedef typename outerProduct<vector, Type>::type GradType;
66 
67  const fvMesh& mesh = vsf.mesh();
68 
69  // Assemble the second-order least-square gradient
70  // Calculate the second-order least-square gradient
71  tmp<GeometricField<GradType, fvPatchField, volMesh> > tsecondfGrad
72  = leastSquaresGrad<Type>(mesh).grad(vsf);
73  const GeometricField<GradType, fvPatchField, volMesh>& secondfGrad =
74  tsecondfGrad();
75 
76  tmp<GeometricField<GradType, fvPatchField, volMesh> > tfGrad
77  (
78  new GeometricField<GradType, fvPatchField, volMesh>
79  (
80  IOobject
81  (
82  "grad("+vsf.name()+')',
83  vsf.instance(),
84  mesh,
87  ),
88  secondfGrad
89  )
90  );
91  GeometricField<GradType, fvPatchField, volMesh>& fGrad = tfGrad();
92 
93  const vectorField& C = mesh.C();
94 
95  const surfaceScalarField& lambda = mesh.weights();
96 
97  // Get reference to least square vectors
98  const leastSquaresVectors& lsv = leastSquaresVectors::New(mesh);
99  const surfaceVectorField& ownLs = lsv.pVectors();
100  const surfaceVectorField& neiLs = lsv.nVectors();
101 
102  // owner/neighbour addressing
103  const unallocLabelList& own = mesh.owner();
104  const unallocLabelList& nei = mesh.neighbour();
105 
106  // Assemble the fourth-order gradient
107 
108  // Internal faces
109  forAll(own, facei)
110  {
111  Type dDotGradDelta = 0.5*
112  (
113  (C[nei[facei]] - C[own[facei]])
114  & (secondfGrad[nei[facei]] - secondfGrad[own[facei]])
115  );
116 
117  fGrad[own[facei]] -= lambda[facei]*ownLs[facei]*dDotGradDelta;
118  fGrad[nei[facei]] -= (1.0 - lambda[facei])*neiLs[facei]*dDotGradDelta;
119  }
120 
121  // Boundary faces
122  forAll(vsf.boundaryField(), patchi)
123  {
124  if (secondfGrad.boundaryField()[patchi].coupled())
125  {
126  const fvsPatchVectorField& patchOwnLs =
127  ownLs.boundaryField()[patchi];
128 
129  const scalarField& lambdap = lambda.boundaryField()[patchi];
130 
131  // Build the d-vectors
132  vectorField pd =
133  mesh.Sf().boundaryField()[patchi]
134  /(
135  mesh.magSf().boundaryField()[patchi]
136  *mesh.deltaCoeffs().boundaryField()[patchi]
137  );
138 
139  if (!mesh.orthogonal())
140  {
141  pd -= mesh.correctionVectors().boundaryField()[patchi]
142  /mesh.deltaCoeffs().boundaryField()[patchi];
143  }
144 
145  const unallocLabelList& faceCells =
146  fGrad.boundaryField()[patchi].patch().faceCells();
147 
148  Field<GradType> neighbourSecondfGrad =
149  secondfGrad.boundaryField()[patchi].patchNeighbourField();
150 
151  forAll(faceCells, patchFaceI)
152  {
153  fGrad[faceCells[patchFaceI]] -=
154  0.5*lambdap[patchFaceI]*patchOwnLs[patchFaceI]
155  *(
156  pd[patchFaceI]
157  & (
158  neighbourSecondfGrad[patchFaceI]
159  - secondfGrad[faceCells[patchFaceI]]
160  )
161  );
162  }
163  }
164  }
165 
166  fGrad.correctBoundaryConditions();
168 
169  return tfGrad;
170 }
171 
172 
173 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
174 
175 } // End namespace fv
176 
177 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
178 
179 } // End namespace Foam
180 
181 // ************************ vim: set sw=4 sts=4 et: ************************ //