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