All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
leastSquaresVectors.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 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  DemandDrivenMeshObject
42  <
43  fvMesh,
44  MoveableMeshObject,
45  leastSquaresVectors
46  >(mesh),
47  pVectors_
48  (
49  IOobject
50  (
51  "LeastSquaresP",
52  mesh.pointsInstance(),
53  mesh,
54  IOobject::NO_READ,
55  IOobject::NO_WRITE,
56  false
57  ),
58  mesh,
60  ),
61  nVectors_
62  (
63  IOobject
64  (
65  "LeastSquaresN",
66  mesh.pointsInstance(),
67  mesh,
68  IOobject::NO_READ,
69  IOobject::NO_WRITE,
70  false
71  ),
72  mesh,
74  )
75 {
76  calcLeastSquaresVectors();
77 }
78 
79 
80 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
81 
83 {}
84 
85 
86 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
87 
88 void Foam::leastSquaresVectors::calcLeastSquaresVectors()
89 {
90  if (debug)
91  {
92  InfoInFunction << "Calculating least square gradient vectors" << endl;
93  }
94 
95  const fvMesh& mesh = this->mesh();
96 
97  // Set local references to mesh data
98  const labelUList& owner = mesh.owner();
99  const labelUList& neighbour = mesh.neighbour();
100 
101  const volVectorField& C = mesh.C();
102  const surfaceScalarField& w = mesh.weights();
103  const surfaceScalarField& magSf = mesh.magSf();
104 
105 
106  // Set up temporary storage for the dd tensor (before inversion)
107  symmTensorField dd(mesh().nCells(), Zero);
108 
109  forAll(owner, facei)
110  {
111  label own = owner[facei];
112  label nei = neighbour[facei];
113 
114  vector d = C[nei] - C[own];
115  symmTensor wdd = (magSf[facei]/magSqr(d))*sqr(d);
116 
117  dd[own] += (1 - w[facei])*wdd;
118  dd[nei] += w[facei]*wdd;
119  }
120 
121 
122  surfaceVectorField::Boundary& pVectorsBf =
123  pVectors_.boundaryFieldRef();
124 
125  forAll(pVectorsBf, patchi)
126  {
127  const fvsPatchScalarField& pw = w.boundaryField()[patchi];
128  const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
129 
130  const fvPatch& p = pw.patch();
131  const labelUList& faceCells = p.faceCells();
132 
133  // Build the d-vectors
134  vectorField pd(p.delta());
135 
136  if (pw.coupled())
137  {
138  forAll(pd, patchFacei)
139  {
140  const vector& d = pd[patchFacei];
141 
142  dd[faceCells[patchFacei]] +=
143  ((1 - pw[patchFacei])*pMagSf[patchFacei]/magSqr(d))*sqr(d);
144  }
145  }
146  else
147  {
148  forAll(pd, patchFacei)
149  {
150  const vector& d = pd[patchFacei];
151 
152  dd[faceCells[patchFacei]] +=
153  (pMagSf[patchFacei]/magSqr(d))*sqr(d);
154  }
155  }
156  }
157 
158 
159  // Invert the dd tensor
160  const symmTensorField invDd(inv(dd));
161 
162 
163  // Revisit all faces and calculate the pVectors_ and nVectors_ vectors
164  forAll(owner, facei)
165  {
166  label own = owner[facei];
167  label nei = neighbour[facei];
168 
169  vector d = C[nei] - C[own];
170  scalar magSfByMagSqrd = magSf[facei]/magSqr(d);
171 
172  pVectors_[facei] = (1 - w[facei])*magSfByMagSqrd*(invDd[own] & d);
173  nVectors_[facei] = -w[facei]*magSfByMagSqrd*(invDd[nei] & d);
174  }
175 
176  forAll(pVectorsBf, patchi)
177  {
178  fvsPatchVectorField& patchLsP = pVectorsBf[patchi];
179 
180  const fvsPatchScalarField& pw = w.boundaryField()[patchi];
181  const fvsPatchScalarField& pMagSf = magSf.boundaryField()[patchi];
182 
183  const fvPatch& p = pw.patch();
184  const labelUList& faceCells = p.faceCells();
185 
186  // Build the d-vectors
187  vectorField pd(p.delta());
188 
189  if (pw.coupled())
190  {
191  forAll(pd, patchFacei)
192  {
193  const vector& d = pd[patchFacei];
194 
195  patchLsP[patchFacei] =
196  ((1 - pw[patchFacei])*pMagSf[patchFacei]/magSqr(d))
197  *(invDd[faceCells[patchFacei]] & d);
198  }
199  }
200  else
201  {
202  forAll(pd, patchFacei)
203  {
204  const vector& d = pd[patchFacei];
205 
206  patchLsP[patchFacei] =
207  pMagSf[patchFacei]*(1.0/magSqr(d))
208  *(invDd[faceCells[patchFacei]] & d);
209  }
210  }
211  }
212 
213  if (debug)
214  {
216  << "Finished calculating least square gradient vectors" << endl;
217  }
218 }
219 
220 
222 {
223  calcLeastSquaresVectors();
224  return true;
225 }
226 
227 
228 // ************************************************************************* //
static const Foam::dimensionedScalar C("C", Foam::dimTemperature, 234.5)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
GeometricBoundaryField< Type, PatchField, GeoMesh > Boundary
Type of the boundary field.
virtual bool movePoints()
Delete the least square vectors when the mesh moves.
virtual ~leastSquaresVectors()
Destructor.
leastSquaresVectors(const fvMesh &)
Construct given an fvMesh.
label patchi
#define InfoInFunction
Report an information message using Foam::Info.
Namespace for OpenFOAM.
static const zero Zero
Definition: zero.H:97
VolField< vector > volVectorField
Definition: volFieldsFwd.H:65
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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
const dimensionSet dimless
SymmTensor< scalar > symmTensor
SymmTensor of scalars.
Definition: symmTensor.H:48
SurfaceField< scalar > surfaceScalarField
const dimensionSet dimLength
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Field< vector > vectorField
Specialisation of Field<T> for vector.
Field< symmTensor > symmTensorField
Specialisation of Field<T> for symmTensor.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
UList< label > labelUList
Definition: UList.H:65
fvsPatchField< vector > fvsPatchVectorField
dimensioned< scalar > magSqr(const dimensioned< Type > &)
fvsPatchField< scalar > fvsPatchScalarField
volScalarField & p