FreeFOAM The Cross-Platform CFD Toolkit
leastSquaresGrad.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 "leastSquaresGrad.H"
27 #include "leastSquaresVectors.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  typedef typename outerProduct<vector, Type>::type GradType;
61 
62  const fvMesh& mesh = vsf.mesh();
63 
64  tmp<GeometricField<GradType, fvPatchField, volMesh> > tlsGrad
65  (
66  new GeometricField<GradType, fvPatchField, volMesh>
67  (
68  IOobject
69  (
70  "grad("+vsf.name()+')',
71  vsf.instance(),
72  mesh,
75  ),
76  mesh,
77  dimensioned<GradType>
78  (
79  "zero",
80  vsf.dimensions()/dimLength,
81  pTraits<GradType>::zero
82  ),
83  zeroGradientFvPatchField<GradType>::typeName
84  )
85  );
86  GeometricField<GradType, fvPatchField, volMesh>& lsGrad = tlsGrad();
87 
88  // Get reference to least square vectors
89  const leastSquaresVectors& lsv = leastSquaresVectors::New(mesh);
90 
91  const surfaceVectorField& ownLs = lsv.pVectors();
92  const surfaceVectorField& neiLs = lsv.nVectors();
93 
94  const unallocLabelList& own = mesh.owner();
95  const unallocLabelList& nei = mesh.neighbour();
96 
97  forAll(own, facei)
98  {
99  register label ownFaceI = own[facei];
100  register label neiFaceI = nei[facei];
101 
102  Type deltaVsf = vsf[neiFaceI] - vsf[ownFaceI];
103 
104  lsGrad[ownFaceI] += ownLs[facei]*deltaVsf;
105  lsGrad[neiFaceI] -= neiLs[facei]*deltaVsf;
106  }
107 
108  // Boundary faces
109  forAll(vsf.boundaryField(), patchi)
110  {
111  const fvsPatchVectorField& patchOwnLs = ownLs.boundaryField()[patchi];
112 
113  const unallocLabelList& faceCells =
114  lsGrad.boundaryField()[patchi].patch().faceCells();
115 
116  if (vsf.boundaryField()[patchi].coupled())
117  {
118  Field<Type> neiVsf =
119  vsf.boundaryField()[patchi].patchNeighbourField();
120 
121  forAll(neiVsf, patchFaceI)
122  {
123  lsGrad[faceCells[patchFaceI]] +=
124  patchOwnLs[patchFaceI]
125  *(neiVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
126  }
127  }
128  else
129  {
130  const fvPatchField<Type>& patchVsf = vsf.boundaryField()[patchi];
131 
132  forAll(patchVsf, patchFaceI)
133  {
134  lsGrad[faceCells[patchFaceI]] +=
135  patchOwnLs[patchFaceI]
136  *(patchVsf[patchFaceI] - vsf[faceCells[patchFaceI]]);
137  }
138  }
139  }
140 
141 
142  lsGrad.correctBoundaryConditions();
144 
145  return tlsGrad;
146 }
147 
148 
149 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
150 
151 } // End namespace fv
152 
153 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
154 
155 } // End namespace Foam
156 
157 // ************************ vim: set sw=4 sts=4 et: ************************ //