All Classes Namespaces Files Functions Variables Typedefs Friends Macros Groups
leastSquaresStencilOpt.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2019 OpenCFD Ltd.
10  Copyright (C) 2016-2019 ISP RAS (www.ispras.ru) UniCFD Group (www.unicfd.ru)
11 -------------------------------------------------------------------------------
12 License
13  This file is part of QGDsolver library, based on OpenFOAM+.
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22  You should have received a copy of the GNU General Public License
23  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
24 Group grpleastSquaresOpt
25  This group contains common part of QGD solvers.
26 Class
27  Foam::fvsc::leastSquaresOpt::leastSquaresStencilOpt
28 Description
29  Methods for optimal calculating of directional derivative.
30  With parallel realisation.
31 \*---------------------------------------------------------------------------*/
32 
33 
34 #include "leastSquaresStencilOpt.H"
35 #include "polyMesh.H"
36 #include "fvMesh.H"
37 #include "word.H"
38 #include "IOstream.H"
39 #include "Ostream.H"
40 #include "wedgeFvPatch.H"
41 #include "symmetryPlaneFvPatch.H"
42 #include "symmetryFvPatch.H"
43 #include "emptyFvPatch.H"
44 #include <HashTable.H>
45 
46 #include "addToRunTimeSelectionTable.H"
47 
48 namespace Foam
49 {
50 namespace fvsc
51 {
52  defineTypeNameAndDebug(leastSquaresOpt,0);
53  addToRunTimeSelectionTable
54  (
55  fvscStencil,
56  leastSquaresOpt,
57  components
58  );
59 }
60 }
61 
62 // constructors
63 Foam::fvsc::leastSquaresOpt::leastSquaresOpt(const IOobject& io)
64 :
65  fvscStencil(io),
66  leastSquaresBase(this->mesh_)
67 {
68 }
69 
70 
71 Foam::fvsc::leastSquaresOpt::~leastSquaresOpt()
72 {
73 }
74 
75 Foam::tmp<Foam::surfaceTensorField> Foam::fvsc::leastSquaresOpt::Grad(const volVectorField& iVF)
76 {
77  return Grad(iVF, linearInterpolate(iVF));
78 }
79 
80 //- Calculate gradient of volume vector field on the faces.
81 //
82 // \param iVF Internal vector field.
83 // Allowable values: constant reference to the volVectorField.
84 //
85 // \return Gradient of iVF (tensor field) which was computed on the faces of mesh.
86 Foam::tmp<Foam::surfaceTensorField> Foam::fvsc::leastSquaresOpt::Grad(const volVectorField& iVF, const surfaceVectorField& sVF)
87 {
88 
89  tmp<surfaceTensorField> tgradIVF(0*nf_*fvc::snGrad(iVF));
90  surfaceScalarField tField = sVF.component(0)*0;
91  surfaceTensorField& gradIVF = tgradIVF.ref();
92 
93  faceScalarDer(iVF.primitiveField().component(0),sVF.primitiveField().component(0),0,tField);
94  gradIVF.primitiveFieldRef().replace(0,tField);
95 
96  faceScalarDer(iVF.primitiveField().component(0),sVF.primitiveField().component(0),1,tField);
97  gradIVF.primitiveFieldRef().replace(3,tField);
98 
99  faceScalarDer(iVF.primitiveField().component(0),sVF.primitiveField().component(0),2,tField);
100  gradIVF.primitiveFieldRef().replace(6,tField);
101 
102  faceScalarDer(iVF.primitiveField().component(1),sVF.primitiveField().component(1),0,tField);
103  gradIVF.primitiveFieldRef().replace(1,tField);
104 
105  faceScalarDer(iVF.primitiveField().component(1),sVF.primitiveField().component(1),1,tField);
106  gradIVF.primitiveFieldRef().replace(4,tField);
107 
108  faceScalarDer(iVF.primitiveField().component(1),sVF.primitiveField().component(1),2,tField);
109  gradIVF.primitiveFieldRef().replace(7,tField);
110 
111  faceScalarDer(iVF.primitiveField().component(2),sVF.primitiveField().component(2),0,tField);
112  gradIVF.primitiveFieldRef().replace(2,tField);
113 
114  faceScalarDer(iVF.primitiveField().component(2),sVF.primitiveField().component(2),1,tField);
115  gradIVF.primitiveFieldRef().replace(5,tField);
116 
117  faceScalarDer(iVF.primitiveField().component(2),sVF.primitiveField().component(2),2,tField);
118  gradIVF.primitiveFieldRef().replace(8,tField);
119 
120 
121  forAll(mesh_.boundaryMesh(), ipatch)
122  {
123  bool notConstrain = true;
124  const fvPatch& fvp = mesh_.boundary()[ipatch];
125  if
126  (
127  isA<emptyFvPatch>(fvp) ||
128  isA<wedgeFvPatch>(fvp) ||
129  isA<coupledFvPatch>(fvp) ||
130  isA<symmetryFvPatch>(fvp) ||
131  isA<symmetryPlaneFvPatch>(fvp)
132  )
133  {
134  notConstrain = false;
135  }
136 
137  if (notConstrain)
138  {
139  gradIVF.boundaryFieldRef()[ipatch] = nf_.boundaryField()[ipatch]*iVF.boundaryField()[ipatch].snGrad();
140  }
141  }
142 
143  if(!Pstream::parRun())
144  {
145  return tgradIVF;
146  }
147 
148  List<List3<scalar>> procVfValues(nProcPatches_);
149 
150  formVfValues(iVF,procVfValues);
151 
152  faceScalarDer(procVfValues[0],sVF.component(0),0,tField);
153  gradIVF.boundaryFieldRef().replace(0,tField.boundaryField());
154 
155  faceScalarDer(procVfValues[0],sVF.component(0),1,tField);
156  gradIVF.boundaryFieldRef().replace(3,tField.boundaryField());
157 
158  faceScalarDer(procVfValues[0],sVF.component(0),2,tField);
159  gradIVF.boundaryFieldRef().replace(6,tField.boundaryField());
160 
161  faceScalarDer(procVfValues[1],sVF.component(1),0,tField);
162  gradIVF.boundaryFieldRef().replace(1,tField.boundaryField());
163 
164  faceScalarDer(procVfValues[1],sVF.component(1),1,tField);
165  gradIVF.boundaryFieldRef().replace(4,tField.boundaryField());
166 
167  faceScalarDer(procVfValues[1],sVF.component(1),2,tField);
168  gradIVF.boundaryFieldRef().replace(7,tField.boundaryField());
169 
170  faceScalarDer(procVfValues[2],sVF.component(2),0,tField);
171  gradIVF.boundaryFieldRef().replace(2,tField.boundaryField());
172 
173  faceScalarDer(procVfValues[2],sVF.component(2),1,tField);
174  gradIVF.boundaryFieldRef().replace(5,tField.boundaryField());
175 
176  faceScalarDer(procVfValues[2],sVF.component(2),2,tField);
177  gradIVF.boundaryFieldRef().replace(8,tField.boundaryField());
178 
179  return tgradIVF;
180 
181 };
182 
183 //- Calculate divergence of volume vector field on the faces.
184 //
185 // \param iVF Internal vector field.
186 // Allowable values: constant reference to the volVectorField.
187 //
188 // \return Divergence of iVF (scalar field) which was computed on the faces of mesh.
189 Foam::tmp<Foam::surfaceScalarField> Foam::fvsc::leastSquaresOpt::Div(const volVectorField& iVF)
190 {
191 
192  surfaceVectorField sVF = linearInterpolate(iVF);
193  surfaceScalarField tField = sVF.component(0)*0;
194  tmp<surfaceScalarField> tdivIVF(0*nf_&fvc::snGrad(iVF));
195  surfaceScalarField& divIVF = tdivIVF.ref();
196 
197  faceScalarDer(iVF.primitiveField().component(0),sVF.primitiveField().component(0),0,tField);
198  divIVF.primitiveFieldRef() = tField;
199 
200  faceScalarDer(iVF.primitiveField().component(1),sVF.primitiveField().component(1),1,tField);
201  divIVF.primitiveFieldRef() += tField;
202 
203  faceScalarDer(iVF.primitiveField().component(2),sVF.primitiveField().component(2),2,tField);
204  divIVF.primitiveFieldRef() += tField;
205 
206  forAll(mesh_.boundaryMesh(), ipatch)
207  {
208  bool notConstrain = true;
209  const fvPatch& fvp = mesh_.boundary()[ipatch];
210  if
211  (
212  isA<emptyFvPatch>(fvp) ||
213  isA<wedgeFvPatch>(fvp) ||
214  isA<coupledFvPatch>(fvp) ||
215  isA<symmetryFvPatch>(fvp) ||
216  isA<symmetryPlaneFvPatch>(fvp)
217  )
218  {
219  notConstrain = false;
220  }
221 
222  if (notConstrain)
223  {
224  divIVF.boundaryFieldRef()[ipatch] = nf_.boundaryField()[ipatch] & iVF.boundaryField()[ipatch].snGrad();
225  }
226  }
227 
228  if(!Pstream::parRun())
229  {
230  return tdivIVF;
231  }
232 
233  List<List3<scalar>> procVfValues(nProcPatches_);
234  formVfValues(iVF,procVfValues);
235 
236 
237  faceScalarDer(procVfValues[0],sVF.component(0),0,tField);
238  divIVF.boundaryFieldRef() = tField.boundaryField();
239 
240  faceScalarDer(procVfValues[1],sVF.component(1),1,tField);
241  divIVF.boundaryFieldRef() += tField.boundaryField();
242 
243  faceScalarDer(procVfValues[2],sVF.component(2),2,tField);
244  divIVF.boundaryFieldRef() += tField.boundaryField();
245 
246 
247  return tdivIVF;
248 };
249 
250 //- Calculate divergence of volume tensor field on the faces.
251 //
252 // \param iTF Internal tensor field.
253 // Allowable values: constant reference to the volTensorField.
254 //
255 // \return Divergence of iTF (vector field) which was computed on the faces of mesh.
256 Foam::tmp<Foam::surfaceVectorField> Foam::fvsc::leastSquaresOpt::Div(const volTensorField& iTF)
257 {
258 
259  surfaceTensorField sTF = linearInterpolate(iTF);
260  surfaceScalarField tField = sTF.component(0)*0;
261  tmp<surfaceVectorField> tdivITF(0*nf_*fvc::snGrad(iTF.component(0)));
262  surfaceVectorField& divITF = tdivITF.ref();
263  surfaceScalarField divComp = tField;
264 
265  faceScalarDer(iTF.primitiveField().component(0),sTF.primitiveField().component(0),0,tField);
266  divComp = tField;
267  faceScalarDer(iTF.primitiveField().component(1),sTF.primitiveField().component(1),1,tField);
268  divComp += tField;
269  faceScalarDer(iTF.primitiveField().component(2),sTF.primitiveField().component(2),2,tField);
270  divComp += tField;
271  divITF.primitiveFieldRef().replace(0,divComp);
272 
273  faceScalarDer(iTF.primitiveField().component(3),sTF.primitiveField().component(3),0,tField);
274  divComp = tField;
275  faceScalarDer(iTF.primitiveField().component(4),sTF.primitiveField().component(4),1,tField);
276  divComp += tField;
277  faceScalarDer(iTF.primitiveField().component(5),sTF.primitiveField().component(5),2,tField);
278  divComp += tField;
279  divITF.primitiveFieldRef().replace(1,divComp);
280 
281  faceScalarDer(iTF.primitiveField().component(6),sTF.primitiveField().component(6),0,tField);
282  divComp = tField;
283  faceScalarDer(iTF.primitiveField().component(7),sTF.primitiveField().component(7),1,tField);
284  divComp += tField;
285  faceScalarDer(iTF.primitiveField().component(8),sTF.primitiveField().component(8),2,tField);
286  divComp += tField;
287  divITF.primitiveFieldRef().replace(2,divComp);
288 
289  forAll(mesh_.boundaryMesh(), ipatch)
290  {
291  bool notConstrain = true;
292  const fvPatch& fvp = mesh_.boundary()[ipatch];
293  if
294  (
295  isA<emptyFvPatch>(fvp) ||
296  isA<wedgeFvPatch>(fvp) ||
297  isA<coupledFvPatch>(fvp) ||
298  isA<symmetryFvPatch>(fvp) ||
299  isA<symmetryPlaneFvPatch>(fvp)
300  )
301  {
302  notConstrain = false;
303  }
304 
305  if (notConstrain)
306  {
307  divITF.boundaryFieldRef()[ipatch] = nf_.boundaryField()[ipatch]
308  & iTF.boundaryField()[ipatch].snGrad();
309  }
310  }
311 
312 
313  if (!Pstream::parRun())
314  {
315  return tdivITF;
316  }
317 
318 
319  List<List3<scalar>> procVfValues(nProcPatches_);
320  formVfValues(iTF,procVfValues);
321 
322  faceScalarDer(procVfValues[0],sTF.component(0),0,tField);
323  divComp = tField;
324  faceScalarDer(procVfValues[1],sTF.component(1),1,tField);
325  divComp += tField;
326  faceScalarDer(procVfValues[2],sTF.component(2),2,tField);
327  divComp += tField;
328  divITF.boundaryFieldRef().replace(0,divComp.boundaryField());
329 
330  faceScalarDer(procVfValues[3],sTF.component(3),0,tField);
331  divComp = tField;
332  faceScalarDer(procVfValues[4],sTF.component(4),1,tField);
333  divComp += tField;
334  faceScalarDer(procVfValues[5],sTF.component(5),2,tField);
335  divComp += tField;
336  divITF.boundaryFieldRef().replace(1,divComp.boundaryField());
337 
338  faceScalarDer(procVfValues[6],sTF.component(6),0,tField);
339  divComp = tField;
340  faceScalarDer(procVfValues[7],sTF.component(7),1,tField);
341  divComp += tField;
342  faceScalarDer(procVfValues[8],sTF.component(8),2,tField);
343  divComp += tField;
344  divITF.boundaryFieldRef().replace(2,divComp.boundaryField());
345 
346  return tdivITF;
347 }
348 
349 //
350 //END-OF-FILE
351 //
352 
353 
Foam::fvsc::fvscStencil
This is a method for calculation the differential operators without tangential derivatives. They are further used in the calculation of the QGD terms.
Definition: fvscStencil.H:38
formVfValues
void formVfValues(const volScalarField &iF, List3< scalar > &procVfValues)
Definition: extendedFaceStencilFormVfValues.H:30
Foam::fvsc::leastSquaresOpt::Div
tmp< surfaceScalarField > Div(const volVectorField &iVF)
Calculate divergence of volume vector field on the faces.
Definition: leastSquaresStencilOpt.C:188
Foam::fvsc::leastSquaresOpt::Grad
tmp< surfaceVectorField > Grad(const volScalarField &iF)
Calculate gradient of volume scalar function on the faces.
Definition: extendedFaceStencilScalarGradOpt.C:50
Foam::fvsc::addToRunTimeSelectionTable
addToRunTimeSelectionTable(fvscStencil, GaussVolPoint, components)
forAll
forAll(Y, i)
Definition: QGDYEqn.H:36
Foam::fvsc::leastSquaresOpt::leastSquaresOpt
leastSquaresOpt(const IOobject &)
Construct from IOobject. Optional flag for if IOobject is the.
Definition: leastSquaresStencilOpt.C:56
Foam::fvsc::defineTypeNameAndDebug
defineTypeNameAndDebug(fvscStencil, 0)

Generated by   doxygen 1.8.5