All Classes Namespaces Files Functions Variables Typedefs Friends Macros Groups
GaussVolPointBase2D.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
25  grpGaussVolPoint
26 Class
27  Foam::fvsc::GaussVolPoint::GaussVolPoint2D
28 Description
29  This is a method for approximating derivatives of tangents to a face (2D case).
30  They are further used in the calculation of the QGD terms.
31 \*---------------------------------------------------------------------------*/
32 #include "GaussVolPointBase2D.H"
33 #include "fvMesh.H"
34 #include "volFields.H"
35 #include "surfaceFields.H"
36 #include "pointFields.H"
37 #include "coupledFvPatch.H"
38 #include "processorFvPatch.H"
39 #include "processorFvPatchFields.H"
40 #include "emptyFvPatch.H"
41 #include "wedgeFvPatch.H"
42 #include "volPointInterpolation.H"
43 
45 :
46  c1_(mesh.neighbour().size(), 0.0),
47  c2_(mesh.neighbour().size(), 0.0),
48  c3_(mesh.neighbour().size(), 0.0),
49  c4_(mesh.neighbour().size(), 0.0),
50  mv42_(mesh.neighbour().size(), 1.0),
51  mv13_(mesh.neighbour().size(), 1.0),
52  ip3_(mesh.neighbour().size(), -1),
53  ip1_(mesh.neighbour().size(), -1),
54  ic4_(mesh.neighbour().size(), -1),
55  ic2_(mesh.neighbour().size(), -1),
56  e1_(1, 0, 0),
57  e2_(0, 1, 0),
58  ie1_(-1),
59  ie2_(-1),
60  ie3_(-1),
61  ordinaryPatches_(0),
62  ip3e_(mesh.boundary().size()),
63  ip1e_(mesh.boundary().size()),
64  ic4e_(mesh.boundary().size()),
65  c1e_(mesh.boundary().size()),
66  c2e_(mesh.boundary().size()),
67  c3e_(mesh.boundary().size()),
68  c4e_(mesh.boundary().size()),
69  mv42e_(mesh.boundary().size()),
70  mv13e_(mesh.boundary().size())
71 {
72  if (mesh.nGeometricD() == 2)
73  {
74  List<scalar> cosa1(mesh.neighbour().size(), 0.0);
75  List<scalar> cosa2(mesh.neighbour().size(), 0.0);
76  List<scalar> sina1(mesh.neighbour().size(), 0.0);
77  List<scalar> sina2(mesh.neighbour().size(), 0.0);
78  List<scalar> den(mesh.neighbour().size(), 1.0);
79  List<vector> v42(mesh.neighbour().size(), vector::one);
80  List<vector> v13(mesh.neighbour().size(), vector::one);
81 
82  List<List<scalar> > cosa1e(mesh.boundary().size());
83  List<List<scalar> > cosa2e(mesh.boundary().size());
84  List<List<scalar> > sina1e(mesh.boundary().size());
85  List<List<scalar> > sina2e(mesh.boundary().size());
86  List<List<scalar> > dene(mesh.boundary().size());
87 
88  label ip1=-1, ic2=-1, ip3=-1, ic4=-1;
89 
90  forAll(mesh.geometricD(), iDir)
91  {
92  if (mesh.geometricD()[iDir] < 1)
93  {
94  ie3_ = iDir;
95  }
96  }
97  if (ie3_ == 0)
98  {
99  e1_ = vector(0, 1, 0);
100  e2_ = vector(0, 0, 1);
101  ie1_ = 1;
102  ie2_ = 2;
103  ie3_ = 0;
104  }
105  if (ie3_ == 1)
106  {
107  e1_ = vector(1, 0, 0);
108  e2_ = vector(0, 0, 1);
109  ie1_ = 0;
110  ie2_ = 2;
111  ie3_ = 1;
112  }
113  if (ie3_ == 2)
114  {
115  e1_ = vector(1, 0, 0);
116  e2_ = vector(0, 1, 0);
117  ie1_ = 0;
118  ie2_ = 1;
119  ie3_ = 2;
120  }
121 
122  forAll(mesh.neighbour(), iFace) //---> for(iFace=0; iFace < mesh_.neighbour().size(); iFace++)
123  {
124  ic2 = mesh.neighbour()[iFace]; //cell center for point #2
125  ic4 = mesh.owner()[iFace]; //cell center for point #4
126 
127  const face& f = mesh.faces()[iFace];
128  #warning "Raise error if number of points on face is not equal to 4"
129  forAll(f, ip)
130  {
131  if (mesh.points()[f[ip]][ie3_] >= mesh.C()[ic2][ie3_])
132  {
133  ip1 = f[ip];
134  break;
135  }
136  }
137  forAll(f, ip)
138  {
139  if (mesh.points()[f[ip]][ie3_] >= mesh.C()[ic2][ie3_])
140  {
141  if (ip1 != f[ip])
142  {
143  ip3 = f[ip];
144  break;
145  }
146  }
147  }
148 
149  ic2_[iFace] = ic2;
150  ic4_[iFace] = ic4;
151  ip3_[iFace] = ip3;
152  ip1_[iFace] = ip1;
153 
154  v42[iFace] = mesh.C()[ic2] - mesh.C()[ic4];
155  v13[iFace] = mesh.points()[ip3] - mesh.points()[ip1];
156  mv42_[iFace] = mag(v42[iFace]);
157  mv13_[iFace] = mag(v13[iFace]);
158 
159  cosa1[iFace] = (v42[iFace]/mv42_[iFace]) & e1_;
160  cosa2[iFace] = (v13[iFace]/mv13_[iFace]) & e1_;
161  sina1[iFace] = (v42[iFace]/mv42_[iFace]) & e2_;
162  sina2[iFace] = (v13[iFace]/mv13_[iFace]) & e2_;
163 
164  den[iFace] = sina2[iFace]*cosa1[iFace] - sina1[iFace]*cosa2[iFace];
165  c1_[iFace] = sina2[iFace] / den[iFace];
166  c2_[iFace] = sina1[iFace] / den[iFace];
167  c3_[iFace] = cosa1[iFace] / den[iFace];
168  c4_[iFace] = cosa2[iFace] / den[iFace];
169  }
170 
171  //Info << "Creating weights" << endl;
172  forAll(mesh.boundary(), iPatch)
173  {
174  label patchId = iPatch;
175  if (
176  !isA<emptyFvPatch>(mesh.boundary()[patchId])
177  &&
178  !isA<wedgeFvPatch>(mesh.boundary()[patchId])
179  )
180  {
181  if
182  (
183  isA<coupledFvPatch>(mesh.boundary()[patchId])
184  &&
185  !isA<processorFvPatch>(mesh.boundary()[patchId])
186  )
187  {
188  continue;
189  }
190  else
191  {
192  if (isA<processorFvPatch>(mesh.boundary()[patchId]))
193  {
194  processorPatch_.append(true);
195  }
196  else
197  {
198  processorPatch_.append(false);
199  }
200  }
201  ordinaryPatches_.append(patchId);
202  List<vector> v42(mesh.boundary()[patchId].size(), vector::one);
203  List<vector> v13(mesh.boundary()[patchId].size(), vector::one);
204  ic4e_[patchId].resize(v42.size());
205  mv42e_[patchId].resize(v42.size());
206  mv13e_[patchId].resize(v42.size());
207  ip3e_[patchId].resize(v42.size());
208  ip1e_[patchId].resize(v42.size());
209  sina1e[patchId].resize(v42.size());
210  sina2e[patchId].resize(v42.size());
211  cosa1e[patchId].resize(v42.size());
212  cosa2e[patchId].resize(v42.size());
213  dene[patchId].resize(v42.size());
214  c1e_[patchId].resize(v42.size());
215  c2e_[patchId].resize(v42.size());
216  c3e_[patchId].resize(v42.size());
217  c4e_[patchId].resize(v42.size());
218 
219  ic4e_[patchId] = mesh.boundary()[patchId].faceCells();
220 
221  //set v42 vector centers for point 2 if patch is processor
222  if (isA<processorFvPatch>(mesh.boundary()[patchId]))
223  {
224  const processorPolyPatch & procPolyPatch =
225  refCast<const processorFvPatch>(mesh.boundary()[patchId]).procPolyPatch();
226 
227  forAll(v42, iFace)
228  {
229  v42[iFace] = procPolyPatch.neighbFaceCellCentres()[iFace] -
230  mesh.C()[ic4e_[patchId][iFace]];
231  }
232  }
233  else
234  {
235  forAll(v42, iFace)
236  {
237  v42[iFace] = 2.0*(mesh.boundary()[patchId].Cf()[iFace] -
238  mesh.C()[ic4e_[patchId][iFace]]);
239  }
240  }
241 
242  forAll(mesh.boundary()[patchId], iFace)
243  {
244  label ip1=-1, ip3=-1, ic4=ic4e_[patchId][iFace];
245  label globalFaceID = mesh.boundary()[patchId].start() + iFace;
246 
247  const face& f = mesh.faces()[globalFaceID];
248  forAll(f, ip)
249  {
250  if (mesh.points()[f[ip]][ie3_] >= mesh.C()[ic4][ie3_])
251  {
252  ip1 = f[ip];
253  break;
254  }
255  }
256  forAll(f, ip)
257  {
258  if (mesh.points()[f[ip]][ie3_] >= mesh.C()[ic4][ie3_])
259  {
260  if (ip1 != f[ip])
261  {
262  ip3 = f[ip];
263  break;
264  }
265  }
266  }
267 
268  ip3e_[patchId][iFace] = ip3;
269  ip1e_[patchId][iFace] = ip1;
270  v13[iFace] = mesh.points()[ip3] - mesh.points()[ip1];
271 
272  mv42e_[patchId][iFace] = mag(v42[iFace]);
273  mv13e_[patchId][iFace] = mag(v13[iFace]);
274 
275  cosa1e[patchId][iFace] = (v42[iFace]/mv42e_[patchId][iFace]) & e1_;
276  cosa2e[patchId][iFace] = (v13[iFace]/mv13e_[patchId][iFace]) & e1_;
277  sina1e[patchId][iFace] = (v42[iFace]/mv42e_[patchId][iFace]) & e2_;
278  sina2e[patchId][iFace] = (v13[iFace]/mv13e_[patchId][iFace]) & e2_;
279 
280  dene[patchId][iFace] =
281  sina2e[patchId][iFace]*cosa1e[patchId][iFace]
282  -
283  sina1e[patchId][iFace]*cosa2e[patchId][iFace];
284 
285  c1e_[patchId][iFace] = sina2e[patchId][iFace] / dene[patchId][iFace];
286  c2e_[patchId][iFace] = sina1e[patchId][iFace] / dene[patchId][iFace];
287  c3e_[patchId][iFace] = cosa1e[patchId][iFace] / dene[patchId][iFace];
288  c4e_[patchId][iFace] = cosa2e[patchId][iFace] / dene[patchId][iFace];
289  }
290  }
291  }
292  }
293 };
295 
297 {
298 }
299 
300 
301 void Foam::fvsc::GaussVolPointBase2D::faceGrad(const volScalarField& f, surfaceVectorField& gradf)
302 {
303  scalar dfdn = 0.0, dfdt = 0.0;
304 
305  if (f.mesh().nGeometricD() == 2)
306  {
307  pointScalarField pF
308  (
309  volPointInterpolation::New(f.mesh()).interpolate
310  (
311  f
312  )
313  );
314 
315  forAll(f.mesh().neighbour(), iFace) //---> for(iFace=0; iFace < mesh_.neighbour().size(); iFace++)
316  {
317  dfdn = (f[ic2_[iFace]] - f[ic4_[iFace]]) / mv42_[iFace];
318  //
319 
320  // df/dl2 (1-3)
321  dfdt = (pF[ip3_[iFace]] - pF[ip1_[iFace]]) / mv13_[iFace];
322  //
323 
324  gradf[iFace][ie1_] = (dfdn*c1_[iFace] - dfdt*c2_[iFace]);
325  gradf[iFace][ie2_] = (dfdt*c3_[iFace] - dfdn*c4_[iFace]);
326  //gradf[iFace][ie1_] = (dfdn*sina2_[iFace] - dfdt*sina1_[iFace])/den_[iFace];
327  //gradf[iFace][ie2_] = (dfdt*cosa1_[iFace] - dfdn*cosa2_[iFace])/den_[iFace];
328  gradf[iFace][ie3_] = 0.0;
329  }
330 
331  List<List<scalar> > psi2(f.boundaryField().size());
332 
333  forAll(ordinaryPatches_, iPatch)
334  {
335  label patchId = ordinaryPatches_[iPatch];
336  if (processorPatch_[iPatch])
337  {
338  psi2[patchId] = refCast<const processorFvPatchField<scalar> >
339  (f.boundaryField()[patchId]).patchNeighbourField();
340  }
341  else
342  {
343  psi2[patchId] = f.boundaryField()[patchId] +
344  f.boundaryField()[patchId].snGrad()
345  *
346  mv42e_[patchId]*0.5;
347  }
348 
349  forAll(f.boundaryField()[patchId], iFace)
350  {
351  //dfdn
352  dfdn = (psi2[patchId][iFace] - f[ic4e_[patchId][iFace]]) / mv42e_[patchId][iFace];
353 
354  //dfdt
355  dfdt = (pF[ip3e_[patchId][iFace]] - pF[ip1e_[patchId][iFace]]) / mv13e_[patchId][iFace];
356 
357  gradf.boundaryFieldRef()[patchId][iFace][ie1_] = (dfdn*c1e_[patchId][iFace] - dfdt*c2e_[patchId][iFace]);
358  gradf.boundaryFieldRef()[patchId][iFace][ie2_] = (dfdt*c3e_[patchId][iFace] - dfdn*c4e_[patchId][iFace]);
359  gradf.boundaryFieldRef()[patchId][iFace][ie3_] = 0.0;
360  }
361  }
362  }
363  else
364  {
365  #warning "Raise error if called not for 2D case"
366  }
367 };
368 
369 void Foam::fvsc::GaussVolPointBase2D::faceDiv(const volVectorField& f, surfaceScalarField& divf)
370 {
371  scalar df1dn = 0.0, df1dt = 0.0,
372  df2dn = 0.0, df2dt = 0.0;
373 
374  if (f.mesh().nGeometricD() == 2)
375  {
376  pointVectorField pF
377  (
378  volPointInterpolation::New(f.mesh()).interpolate
379  (
380  f
381  )
382  );
383 
384  forAll(f.mesh().neighbour(), iFace) //---> for(iFace=0; iFace < mesh_.neighbour().size(); iFace++)
385  {
386  df1dn = (f[ic2_[iFace]][ie1_] - f[ic4_[iFace]][ie1_]) / mv42_[iFace];
387  df2dn = (f[ic2_[iFace]][ie2_] - f[ic4_[iFace]][ie2_]) / mv42_[iFace];
388  //
389 
390  // df/dl2 (1-3)
391  df1dt = (pF[ip3_[iFace]][ie1_] - pF[ip1_[iFace]][ie1_]) / mv13_[iFace];
392  df2dt = (pF[ip3_[iFace]][ie2_] - pF[ip1_[iFace]][ie2_]) / mv13_[iFace];
393  //
394 
395  divf[iFace] = (df1dn*c1_[iFace] - df1dt*c2_[iFace]) +
396  (df2dt*c3_[iFace] - df2dn*c4_[iFace]);
397  }
398 
399  List<List<vector> > psi2(f.boundaryField().size());
400 
401  forAll(ordinaryPatches_, iPatch)
402  {
403  label patchId = ordinaryPatches_[iPatch];
404  if (processorPatch_[iPatch])
405  {
406  psi2[patchId] = refCast<const processorFvPatchField<vector> >
407  (f.boundaryField()[patchId]).patchNeighbourField();
408  }
409  else
410  {
411  psi2[patchId] = f.boundaryField()[patchId] +
412  f.boundaryField()[patchId].snGrad()
413  *
414  mv42e_[patchId]*0.5;
415  }
416 
417  forAll(f.boundaryField()[patchId], iFace)
418  {
419  //dfdn
420  df1dn = (psi2[patchId][iFace][ie1_] - f[ic4e_[patchId][iFace]][ie1_]) / mv42e_[patchId][iFace];
421  df2dn = (psi2[patchId][iFace][ie2_] - f[ic4e_[patchId][iFace]][ie2_]) / mv42e_[patchId][iFace];
422 
423  //dfdt
424  df1dt = (pF[ip3e_[patchId][iFace]][ie1_] - pF[ip1e_[patchId][iFace]][ie1_]) / mv13e_[patchId][iFace];
425  df2dt = (pF[ip3e_[patchId][iFace]][ie2_] - pF[ip1e_[patchId][iFace]][ie2_]) / mv13e_[patchId][iFace];
426 
427  divf.boundaryFieldRef()[patchId][iFace] =
428  (df1dn*c1e_[patchId][iFace] - df1dt*c2e_[patchId][iFace])
429  +
430  (df2dt*c3e_[patchId][iFace] - df2dn*c4e_[patchId][iFace]);
431  }
432  }
433  }
434  else
435  {
436  #warning "Raise error if called not for 2D case"
437  }
438 };
439 
440 void Foam::fvsc::GaussVolPointBase2D::faceDiv(const volTensorField& f, surfaceVectorField& divf)
441 {
442  scalar df11dn = 0.0, df11dt = 0.0,
443  df21dn = 0.0, df21dt = 0.0,
444  df22dn = 0.0, df22dt = 0.0,
445  df12dn = 0.0, df12dt = 0.0;
446 
447  label i11 = ie1_*3 + ie1_;
448  label i21 = ie2_*3 + ie1_;
449  label i22 = ie2_*3 + ie2_;
450  label i12 = ie1_*3 + ie2_;
451 
452  if (f.mesh().nGeometricD() == 2)
453  {
454  pointTensorField pF
455  (
456  volPointInterpolation::New(f.mesh()).interpolate
457  (
458  f
459  )
460  );
461 
462  forAll(f.mesh().neighbour(), iFace) //---> for(iFace=0; iFace < mesh_.neighbour().size(); iFace++)
463  {
464  df11dn = (f[ic2_[iFace]][i11] - f[ic4_[iFace]][i11]) / mv42_[iFace];
465  df21dn = (f[ic2_[iFace]][i21] - f[ic4_[iFace]][i21]) / mv42_[iFace];
466 
467  // df/dl2 (1-3)
468  df11dt = (pF[ip3_[iFace]][i11] - pF[ip1_[iFace]][i11]) / mv13_[iFace];
469  df21dt = (pF[ip3_[iFace]][i21] - pF[ip1_[iFace]][i21]) / mv13_[iFace];
470 
471  divf[iFace][ie1_] =
472  (df11dn*c1_[iFace] - df11dt*c2_[iFace]) + //d/dx
473  (df21dt*c3_[iFace] - df21dn*c4_[iFace]); //d/dy
474 
475  df22dn = (f[ic2_[iFace]][i22] - f[ic4_[iFace]][i22]) / mv42_[iFace];
476  df12dn = (f[ic2_[iFace]][i12] - f[ic4_[iFace]][i12]) / mv42_[iFace];
477 
478  // df/dl2 (1-3)
479  df22dt = (pF[ip3_[iFace]][i22] - pF[ip1_[iFace]][i22]) / mv13_[iFace];
480  df12dt = (pF[ip3_[iFace]][i12] - pF[ip1_[iFace]][i12]) / mv13_[iFace];
481 
482  divf[iFace][ie2_] =
483  (df12dn*c1_[iFace] - df12dt*c2_[iFace]) + //d/dx
484  (df22dt*c3_[iFace] - df22dn*c4_[iFace]); //d/dy
485  }
486 
487  List<List<tensor> > psi2(f.boundaryField().size());
488 
489  forAll(ordinaryPatches_, iPatch)
490  {
491  label patchId = ordinaryPatches_[iPatch];
492  if (processorPatch_[iPatch])
493  {
494  psi2[patchId] = refCast<const processorFvPatchField<tensor> >
495  (f.boundaryField()[patchId]).patchNeighbourField();
496  }
497  else
498  {
499  psi2[patchId] = f.boundaryField()[patchId] +
500  f.boundaryField()[patchId].snGrad()
501  *
502  mv42e_[patchId]*0.5;
503  }
504 
505  forAll(f.boundaryField()[patchId], iFace)
506  {
507  //dfdn
508  df11dn = (psi2[patchId][iFace][i11] - f[ic4e_[patchId][iFace]][i11]) / mv42e_[patchId][iFace];
509  df21dn = (psi2[patchId][iFace][i21] - f[ic4e_[patchId][iFace]][i21]) / mv42e_[patchId][iFace];
510 
511  //dfdt
512  df11dt = (pF[ip3e_[patchId][iFace]][i11] - pF[ip1e_[patchId][iFace]][i11]) / mv13e_[patchId][iFace];
513  df21dt = (pF[ip3e_[patchId][iFace]][i21] - pF[ip1e_[patchId][iFace]][i21]) / mv13e_[patchId][iFace];
514 
515  divf.boundaryFieldRef()[patchId][iFace][ie1_] =
516  (df11dn*c1e_[patchId][iFace] - df11dt*c2e_[patchId][iFace]) //d/dx
517  +
518  (df21dt*c3e_[patchId][iFace] - df21dn*c4e_[patchId][iFace]); //d/dy
519 
520  //dfdn
521  df12dn = (psi2[patchId][iFace][i12] - f[ic4e_[patchId][iFace]][i12]) / mv42e_[patchId][iFace];
522  df22dn = (psi2[patchId][iFace][i22] - f[ic4e_[patchId][iFace]][i22]) / mv42e_[patchId][iFace];
523 
524  //dfdt
525  df12dt = (pF[ip3e_[patchId][iFace]][i12] - pF[ip1e_[patchId][iFace]][i12]) / mv13e_[patchId][iFace];
526  df22dt = (pF[ip3e_[patchId][iFace]][i22] - pF[ip1e_[patchId][iFace]][i22]) / mv13e_[patchId][iFace];
527 
528  divf.boundaryFieldRef()[patchId][iFace][ie2_] =
529  (df12dn*c1e_[patchId][iFace] - df12dt*c2e_[patchId][iFace]) //d/dx
530  +
531  (df22dt*c3e_[patchId][iFace] - df22dn*c4e_[patchId][iFace]); //d/dy
532  }
533  }
534  }
535  else
536  {
537  #warning "Raise error if called not for 2D case"
538  }
539 }
540 
541 //
542 //END-OF-FILE
543 //
544 
545 
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &vf, const surfaceScalarField &dir, const word &reconFieldName=word::null)
Interpolate field vf according to direction dir.
List< label > ordinaryPatches_
ordinary external boundaries
void faceDiv(const volVectorField &f, surfaceScalarField &divf)
forAll(Y, i)
Definition: QGDYEqn.H:36
GaussVolPointBase2D(const fvMesh &mesh)
void faceGrad(const volScalarField &f, surfaceVectorField &gradf)