leastSquaresVectors.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) 2011-2016 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 
26 #include "leastSquaresVectors.H"
27 #include "volFields.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineTypeNameAndDebug(leastSquaresVectors, 0);
34 }
35 
36 
37 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
38 
40 :
41  MeshObject<fvMesh, Foam::MoveableMeshObject, leastSquaresVectors>(mesh),
42  pVectors_
43  (
44  IOobject
45  (
46  "LeastSquaresP",
47  mesh_.pointsInstance(),
48  mesh_,
49  IOobject::NO_READ,
50  IOobject::NO_WRITE,
51  false
52  ),
53  mesh_,
55  ),
56  nVectors_
57  (
58  IOobject
59  (
60  "LeastSquaresN",
61  mesh_.pointsInstance(),
62  mesh_,
63  IOobject::NO_READ,
64  IOobject::NO_WRITE,
65  false
66  ),
67  mesh_,
69  )
70 {
71  calcLeastSquaresVectors();
72 }
73 
74 
75 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
76 
78 {}
79 
80 
81 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
82 
83 void Foam::leastSquaresVectors::calcLeastSquaresVectors()
84 {
85  if (debug)
86  {
87  InfoInFunction << "Calculating least square gradient vectors" << endl;
88  }
89 
90  const fvMesh& mesh = mesh_;
91 
92  // Set local references to mesh data
93  const labelUList& owner = mesh_.owner();
94  const labelUList& neighbour = mesh_.neighbour();
95 
96  const volVectorField& C = mesh.C();
97  const surfaceScalarField& w = mesh.weights();
98  const surfaceScalarField& magSf = mesh.magSf();
99 
100 
101  // Set up temporary storage for the dd tensor (before inversion)
103 
104  forAll(owner, facei)
105  {
106  label own = owner[facei];
107  label nei = neighbour[facei];
108 
109  vector d = C[nei] - C[own];
110  symmTensor wdd = (magSf[facei]/magSqr(d))*sqr(d);
111 
112  dd[own] += (1 - w[facei])*wdd;
113  dd[nei] += w[facei]*wdd;
114  }
115 
116 
117  surfaceVectorField::Boundary& pVectorsBf =
118  pVectors_.boundaryFieldRef();
119 
120  forAll(pVectorsBf, patchi)
121  {
122  const fvsPatchScalarField& pw = w.boundaryField()[patchi];
123  const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
124 
125  const fvPatch& p = pw.patch();
126  const labelUList& faceCells = p.patch().faceCells();
127 
128  // Build the d-vectors
129  vectorField pd(p.delta());
130 
131  if (pw.coupled())
132  {
133  forAll(pd, patchFacei)
134  {
135  const vector& d = pd[patchFacei];
136 
137  dd[faceCells[patchFacei]] +=
138  ((1 - pw[patchFacei])*pMagSf[patchFacei]/magSqr(d))*sqr(d);
139  }
140  }
141  else
142  {
143  forAll(pd, patchFacei)
144  {
145  const vector& d = pd[patchFacei];
146 
147  dd[faceCells[patchFacei]] +=
148  (pMagSf[patchFacei]/magSqr(d))*sqr(d);
149  }
150  }
151  }
152 
153 
154  // Invert the dd tensor
155  const symmTensorField invDd(inv(dd));
156 
157 
158  // Revisit all faces and calculate the pVectors_ and nVectors_ vectors
159  forAll(owner, facei)
160  {
161  label own = owner[facei];
162  label nei = neighbour[facei];
163 
164  vector d = C[nei] - C[own];
165  scalar magSfByMagSqrd = magSf[facei]/magSqr(d);
166 
167  pVectors_[facei] = (1 - w[facei])*magSfByMagSqrd*(invDd[own] & d);
168  nVectors_[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
169  }
170 
171  forAll(pVectorsBf, patchi)
172  {
173  fvsPatchVectorField& patchLsP = pVectorsBf[patchi];
174 
175  const fvsPatchScalarField& pw = w.boundaryField()[patchi];
176  const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
177 
178  const fvPatch& p = pw.patch();
179  const labelUList& faceCells = p.faceCells();
180 
181  // Build the d-vectors
182  vectorField pd(p.delta());
183 
184  if (pw.coupled())
185  {
186  forAll(pd, patchFacei)
187  {
188  const vector& d = pd[patchFacei];
189 
190  patchLsP[patchFacei] =
191  ((1 - pw[patchFacei])*pMagSf[patchFacei]/magSqr(d))
192  *(invDd[faceCells[patchFacei]] & d);
193  }
194  }
195  else
196  {
197  forAll(pd, patchFacei)
198  {
199  const vector& d = pd[patchFacei];
200 
201  patchLsP[patchFacei] =
202  pMagSf[patchFacei]*(1.0/magSqr(d))
203  *(invDd[faceCells[patchFacei]] & d);
204  }
205  }
206  }
207 
208  if (debug)
209  {
211  << "Finished calculating least square gradient vectors" << endl;
212  }
213 }
214 
215 
217 {
218  calcLeastSquaresVectors();
219  return true;
220 }
221 
222 
223 // ************************************************************************* //
fvsPatchField< vector > fvsPatchVectorField
fvsPatchField< scalar > fvsPatchScalarField
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
Field< symmTensor > symmTensorField
Specialisation of Field<T> for symmTensor.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
leastSquaresVectors(const fvMesh &)
Construct given an fvMesh.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
UList< label > labelUList
Definition: UList.H:64
virtual bool movePoints()
Delete the least square vectors when the mesh moves.
label nCells() const
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:288
SymmTensor< scalar > symmTensor
SymmTensor of scalars.
Definition: symmTensor.H:48
static const zero Zero
Definition: zero.H:91
dimensioned< scalar > magSqr(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
label patchi
virtual ~leastSquaresVectors()
Destructor.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
Field< vector > vectorField
Specialisation of Field<T> for vector.
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:282
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Namespace for OpenFOAM.
#define InfoInFunction
Report an information message using Foam::Info.