All Classes Namespaces Files Functions Variables Typedefs Friends Macros Groups
GaussVolPointBase3D.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::GaussVolPoint3D
28 Description
29  This is a method for approximating derivatives of tangents to a face (3D case).
30  They are further used in the calculation of the QGD terms.
31 \*---------------------------------------------------------------------------*/
32 #include "GaussVolPointBase3D.H"
33 #include "fvMesh.H"
34 #include "volFields.H"
35 #include "surfaceFields.H"
36 #include "fvcSnGrad.H"
37 #include "volPointInterpolation.H"
38 #include "processorFvPatch.H"
39 #include "processorFvPatchFields.H"
40 
41 Foam::fvsc::GaussVolPointBase3D::GaussVolPointBase3D(const fvMesh& mesh)
42 :
43  volPoint_
44  (
45  volPointInterpolation::New(mesh)
46  ),
47  nfRef_(mesh.thisDb().lookupObject<surfaceVectorField>("nf")),
48  bgfid_(mesh.boundary().size()),
49  processorPatch_(mesh.boundary().size(), false),
50  qf_(0),
51  aqx_(0),
52  aqy_(0),
53  aqz_(0),
54  vq_(0),
55  bqf_(mesh.boundary().size()),
56  baqx_(mesh.boundary().size()),
57  baqy_(mesh.boundary().size()),
58  baqz_(mesh.boundary().size()),
59  bvq_(mesh.boundary().size()),
60  tf_(0),
61  atx_(0),
62  aty_(0),
63  atz_(0),
64  vt_(0),
65  btf_(mesh.boundary().size()),
66  batx_(mesh.boundary().size()),
67  baty_(mesh.boundary().size()),
68  batz_(mesh.boundary().size()),
69  bvt_(mesh.boundary().size()),
70  bmvON_(mesh.boundary().size()),
71  of_(0),
72  bof_(mesh.boundary().size())
73 {
74  forAll(bgfid_, iPatch)
75  {
76  if (isA<processorFvPatch>(mesh.boundary()[iPatch]))
77  {
78  processorPatch_[iPatch] = true;
79  }
80  const fvPatch& fvp = mesh.boundary()[iPatch];
81  bgfid_[iPatch].resize(fvp.size());
82  forAll(fvp, i)
83  {
84  bgfid_[iPatch][i] = mesh.boundary()[iPatch].start() + i;
85  }
86  }
87 
88  //sort quad, tri faces and other
89  const faceList& faces = mesh.faces();
90 
91  forAll(faces, i)
92  {
93  if (mesh.isInternalFace(i))
94  {
95  if (faces[i].size() == 3)
96  {
97  tf_.append(i);
98  }
99  else if (faces[i].size() == 4)
100  {
101  qf_.append(i);
102  }
103  else
104  {
105  of_.append(i);
106  }
107  }
108  }
109  forAll(bgfid_, iPatch)
110  {
111  const fvPatch& fvp = mesh.boundary()[iPatch];
112  forAll(fvp, i)
113  {
114  if (faces[bgfid_[iPatch][i]].size() == 3)
115  {
116  btf_[iPatch].append(i);
117  }
118  else if (faces[bgfid_[iPatch][i]].size() == 4)
119  {
120  bqf_[iPatch].append(i);
121  }
122  else
123  {
124  bof_[iPatch].append(i);
125  }
126  }
127  }
128 
129  forAll(bgfid_, iPatch)
130  {
131  vectorField vO(mesh.boundary()[iPatch].size(), vector::zero);
132  vectorField vN(mesh.boundary()[iPatch].size(), vector::zero);
133 
134  vO = mesh.boundary()[iPatch].Cn();
135  if (processorPatch_[iPatch])
136  {
137  vN = refCast<const processorFvPatch>(mesh.boundary()[iPatch]).
138  procPolyPatch().neighbFaceCellCentres();
139  }
140  else
141  {
142  vN = vO + 2.0*
143  (
144  mesh.boundary()[iPatch].Cf()
145  -
146  vO
147  );
148  }
149 
150  bmvON_[iPatch] = mag
151  (
152  vO - vN
153  );
154  }
155 
156  triCalcWeights(mesh);
157 
158  quaCalcWeights(mesh);
159 };
160 
161 void Foam::fvsc::GaussVolPointBase3D::triCalcWeights(const fvMesh& mesh)
162 {
163  const pointField& points = mesh.points();
164  const faceList& faces = mesh.faces();
165  label facei = -1;
166  atx_.resize(tf_.size());
167  aty_.resize(tf_.size());
168  atz_.resize(tf_.size());
169  vt_.resize(tf_.size());
170  label own = -1;
171  label nei = -1;
172  label p1 = -1, p2 = -1, p3 = -1;
173  const scalar OneBySix = (1.0 / 6.0);
174 
175  forAll(tf_, i)
176  {
177  facei = tf_[i];
178  own = mesh.owner()[facei];
179  nei = mesh.neighbour()[facei];
180  p1 = faces[facei][0];
181  p2 = faces[facei][1];
182  p3 = faces[facei][2];
183  //p4 - is nei
184  //p5 - is own
185 
186  vt_[i] =
187  (
188  (points[p2] - points[p1]) ^ (points[p3] - points[p1])
189  ) & (mesh.C()[own] - mesh.C()[nei]);
190  vt_[i] *= OneBySix;
191 
192  /* Coefficients for X */
193  atx_[i].resize(5);
194  atx_[i][0] = OneBySix*((mesh.C()[own].z() - mesh.C()[nei].z())*(points[p2].y() - points[p3].y()) +
195  (mesh.C()[nei].y() - mesh.C()[own].y())*(points[p2].z() - points[p3].z()));
196  atx_[i][1] = OneBySix*((mesh.C()[nei].y() - mesh.C()[own].y())*(points[p3].z() - points[p1].z()) +
197  (mesh.C()[own].z() - mesh.C()[nei].z())*(points[p3].y() - points[p1].y()));
198  atx_[i][2] = OneBySix*((mesh.C()[nei].y() - mesh.C()[own].y())*(points[p1].z() - points[p2].z()) +
199  (mesh.C()[own].z() - mesh.C()[nei].z())*(points[p1].y() - points[p2].y()));
200  atx_[i][3] = OneBySix*(points[p1].z()*(points[p2].y() - points[p3].y()) +
201  points[p2].z()*(points[p3].y() - points[p1].y()) +
202  points[p3].z()*(points[p1].y() - points[p2].y()));
203  atx_[i][4] = -atx_[i][3];
204 
205  /* Coefficients for Y */
206  aty_[i].resize(5);
207  aty_[i][0] = OneBySix*((mesh.C()[own].x() - mesh.C()[nei].x())*(points[p2].z() - points[p3].z()) +
208  (mesh.C()[nei].z() - mesh.C()[own].z())*(points[p2].x() - points[p3].x()));
209  aty_[i][1] = OneBySix*((mesh.C()[nei].z() - mesh.C()[own].z())*(points[p3].x() - points[p1].x()) +
210  (mesh.C()[own].x() - mesh.C()[nei].x())*(points[p3].z() - points[p1].z()));
211  aty_[i][2] = OneBySix*((mesh.C()[nei].z() - mesh.C()[own].z())*(points[p1].x() - points[p2].x()) +
212  (mesh.C()[own].x() - mesh.C()[nei].x())*(points[p1].z() - points[p2].z()));
213  aty_[i][3] = OneBySix*(points[p1].x()*(points[p2].z() - points[p3].z()) +
214  points[p2].x()*(points[p3].z() - points[p1].z()) +
215  points[p3].x()*(points[p1].z() - points[p2].z()));
216  aty_[i][4] = -aty_[i][3];
217 
218  /* Coefficients for Z */
219  atz_[i].resize(5);
220  atz_[i][0] = OneBySix*((mesh.C()[own].y() - mesh.C()[nei].y())*(points[p2].x() - points[p3].x()) +
221  (mesh.C()[nei].x() - mesh.C()[own].x())*(points[p2].y() - points[p3].y()));
222  atz_[i][1] = OneBySix*((mesh.C()[nei].x() - mesh.C()[own].x())*(points[p3].y() - points[p1].y()) +
223  (mesh.C()[own].y() - mesh.C()[nei].y())*(points[p3].x() - points[p1].x()));
224  atz_[i][2] = OneBySix*((mesh.C()[nei].x() - mesh.C()[own].x())*(points[p1].y() - points[p2].y()) +
225  (mesh.C()[own].y() - mesh.C()[nei].y())*(points[p1].x() - points[p2].x()));
226  atz_[i][3] = OneBySix*(points[p1].y()*(points[p2].x() - points[p3].x()) +
227  points[p2].y()*(points[p3].x() - points[p1].x()) +
228  points[p3].y()*(points[p1].x() - points[p2].x()));
229  atz_[i][4] = -atz_[i][3];
230  }
231 
232  forAll(btf_, iPatch)
233  {
234  batx_[iPatch].resize(btf_[iPatch].size());
235  baty_[iPatch].resize(btf_[iPatch].size());
236  batz_[iPatch].resize(btf_[iPatch].size());
237  bvt_[iPatch].resize(btf_[iPatch].size());
238 
239  vectorField v5(mesh.boundary()[iPatch].size(), vector::zero);
240  vectorField v4(mesh.boundary()[iPatch].size(), vector::zero);
241 
242  v5 = mesh.boundary()[iPatch].Cn();
243  if (processorPatch_[iPatch])
244  {
245  v4 = refCast<const processorFvPatch>(mesh.boundary()[iPatch]).
246  procPolyPatch().neighbFaceCellCentres();
247  }
248  else
249  {
250  v4 = v5 + 2.0*
251  (
252  mesh.boundary()[iPatch].Cf()
253  -
254  v5
255  );
256  }
257 
258  label gFaceId = -1;
259 
260  forAll(btf_[iPatch], k)
261  {
262  facei = btf_[iPatch][k];
263  gFaceId = bgfid_[iPatch][facei];
264 
265  p1 = faces[gFaceId][0];
266  p2 = faces[gFaceId][1];
267  p3 = faces[gFaceId][2];
268 
269  //p4 - is nei and stored in v4
270  //p5 - is own and stored in v5
271 
272  bvt_[iPatch][k] =
273  ((points[p2] - points[p1]) ^ (points[p3] - points[p1]))
274  & (v5[facei] - v4[facei]);
275  //bvt_[iPatch][k] *= TwoBySix;
276  bvt_[iPatch][k] *= OneBySix;
277 
278  /* Coefficients for X */
279  batx_[iPatch][k].resize(5);
280  batx_[iPatch][k][0] = OneBySix*((v5[facei].z() - v4[facei].z())*(points[p2].y() - points[p3].y()) +
281  (v4[facei].y() - v5[facei].y())*(points[p2].z() - points[p3].z()));
282  batx_[iPatch][k][1] = OneBySix*((v4[facei].y() - v5[facei].y())*(points[p3].z() - points[p1].z()) +
283  (v5[facei].z() - v4[facei].z())*(points[p3].y() - points[p1].y()));
284  batx_[iPatch][k][2] = OneBySix*((v4[facei].y() - v5[facei].y())*(points[p1].z() - points[p2].z()) +
285  (v5[facei].z() - v4[facei].z())*(points[p1].y() - points[p2].y()));
286  batx_[iPatch][k][3] = OneBySix*(points[p1].z()*(points[p2].y() - points[p3].y()) +
287  points[p2].z()*(points[p3].y() - points[p1].y()) +
288  points[p3].z()*(points[p1].y() - points[p2].y()));
289  batx_[iPatch][k][4] = -batx_[iPatch][k][3];
290 
291  /* Coefficients for Y */
292  baty_[iPatch][k].resize(5);
293  baty_[iPatch][k][0] = OneBySix*((v5[facei].x() - v4[facei].x())*(points[p2].z() - points[p3].z()) +
294  (v4[facei].z() - v5[facei].z())*(points[p2].x() - points[p3].x()));
295  baty_[iPatch][k][1] = OneBySix*((v4[facei].z() - v5[facei].z())*(points[p3].x() - points[p1].x()) +
296  (v5[facei].x() - v4[facei].x())*(points[p3].z() - points[p1].z()));
297  baty_[iPatch][k][2] = OneBySix*((v4[facei].z() - v5[facei].z())*(points[p1].x() - points[p2].x()) +
298  (v5[facei].x() - v4[facei].x())*(points[p1].z() - points[p2].z()));
299  baty_[iPatch][k][3] = OneBySix*(points[p1].x()*(points[p2].z() - points[p3].z()) +
300  points[p2].x()*(points[p3].z() - points[p1].z()) +
301  points[p3].x()*(points[p1].z() - points[p2].z()));
302  baty_[iPatch][k][4] = -baty_[iPatch][k][3];
303 
304  /* Coefficients for Z */
305  batz_[iPatch][k].resize(5);
306  batz_[iPatch][k][0] = OneBySix*((v5[facei].y() - v4[facei].y())*(points[p2].x() - points[p3].x()) +
307  (v4[facei].x() - v5[facei].x())*(points[p2].y() - points[p3].y()));
308  batz_[iPatch][k][1] = OneBySix*((v4[facei].x() - v5[facei].x())*(points[p3].y() - points[p1].y()) +
309  (v5[facei].y() - v4[facei].y())*(points[p3].x() - points[p1].x()));
310  batz_[iPatch][k][2] = OneBySix*((v4[facei].x() - v5[facei].x())*(points[p1].y() - points[p2].y()) +
311  (v5[facei].y() - v4[facei].y())*(points[p1].x() - points[p2].x()));
312  batz_[iPatch][k][3] = OneBySix*(points[p1].y()*(points[p2].x() - points[p3].x()) +
313  points[p2].y()*(points[p3].x() - points[p1].x()) +
314  points[p3].y()*(points[p1].x() - points[p2].x()));
315  batz_[iPatch][k][4] = -batz_[iPatch][k][3];
316  }
317  }
318 }
319 
320 void Foam::fvsc::GaussVolPointBase3D::quaCalcWeights(const fvMesh& mesh)
321 {
322  const pointField& points = mesh.points();
323  const faceList& faces = mesh.faces();
324  const scalar OneBySix = (1.0 / 6.0);
325  label facei = -1;
326  aqx_.resize(qf_.size());
327  aqy_.resize(qf_.size());
328  aqz_.resize(qf_.size());
329  vq_.resize(qf_.size());
330  label own = -1;
331  label nei = -1;
332  label p4 = -1, p1 = -1, p2 = -1, p3 = -1;
333 
334  forAll(qf_, i)
335  {
336  facei = qf_[i];
337  own = mesh.owner()[facei];
338  nei = mesh.neighbour()[facei];
339  p1 = faces[facei][0];
340  p2 = faces[facei][1];
341  p3 = faces[facei][2];
342  p4 = faces[facei][3];
343  //p5 - is nei
344  //p6 - is own
345 
346  vq_[i] = (points[p3] - points[p1]) &
347  (
348  (points[p4] - points[p2]) ^ (mesh.C()[own] - mesh.C()[nei])
349  );
350  vq_[i] *= OneBySix;
351 
352  /* coefficient for X */
353  aqx_[i].resize(6);
354  aqx_[i][0] = OneBySix*((mesh.C()[nei].y() - mesh.C()[own].y())*(points[p2].z() - points[p4].z())
355  - (mesh.C()[nei].z() - mesh.C()[own].z())*(points[p2].y() - points[p4].y()));
356  aqx_[i][1] = OneBySix*((mesh.C()[nei].y() - mesh.C()[own].y())*(points[p3].z() - points[p1].z())
357  - (mesh.C()[nei].z() - mesh.C()[own].z())*(points[p3].y() - points[p1].y()));
358  aqx_[i][5] = OneBySix*((points[p1].y() - points[p3].y())*(points[p2].z() - points[p4].z())
359  - (points[p1].z() - points[p3].z())*(points[p2].y() - points[p4].y()));
360 
361  aqx_[i][2] = -aqx_[i][0];
362  aqx_[i][3] = -aqx_[i][1];
363  aqx_[i][4] = -aqx_[i][5];
364 
365  /* coefficient for Y */
366  aqy_[i].resize(6);
367  aqy_[i][0] = OneBySix*((mesh.C()[nei].z() - mesh.C()[own].z())*(points[p2].x() - points[p4].x())
368  - (mesh.C()[nei].x() - mesh.C()[own].x())*(points[p2].z() - points[p4].z()));
369  aqy_[i][1] = OneBySix*((mesh.C()[nei].z() - mesh.C()[own].z())*(points[p3].x() - points[p1].x())
370  - (mesh.C()[nei].x() - mesh.C()[own].x())*(points[p3].z() - points[p1].z()));
371  aqy_[i][5] = OneBySix*((points[p1].z() - points[p3].z())*(points[p2].x() - points[p4].x())
372  - (points[p1].x() - points[p3].x())*(points[p2].z() - points[p4].z()));
373 
374  aqy_[i][2] = -aqy_[i][0];
375  aqy_[i][3] = -aqy_[i][1];
376  aqy_[i][4] = -aqy_[i][5];
377 
378  /* coefficient for Z */
379  aqz_[i].resize(6);
380  aqz_[i][0] = OneBySix*((mesh.C()[nei].x() - mesh.C()[own].x())*(points[p2].y() - points[p4].y())
381  - (mesh.C()[nei].y() - mesh.C()[own].y())*(points[p2].x() - points[p4].x()));
382  aqz_[i][1] = OneBySix*((mesh.C()[nei].x() - mesh.C()[own].x())*(points[p3].y() - points[p1].y())
383  - (mesh.C()[nei].y() - mesh.C()[own].y())*(points[p3].x() - points[p1].x()));
384  aqz_[i][5] = OneBySix*((points[p1].x() - points[p3].x())*(points[p2].y() - points[p4].y())
385  - (points[p1].y() - points[p3].y())*(points[p2].x() - points[p4].x()));
386 
387  aqz_[i][2] = -aqz_[i][0];
388  aqz_[i][3] = -aqz_[i][1];
389  aqz_[i][4] = -aqz_[i][5];
390  }
391  forAll(bqf_, iPatch)
392  {
393  baqx_[iPatch].resize(bqf_[iPatch].size());
394  baqy_[iPatch].resize(bqf_[iPatch].size());
395  baqz_[iPatch].resize(bqf_[iPatch].size());
396  bvq_[iPatch].resize(bqf_[iPatch].size());
397 
398  vectorField v6(mesh.boundary()[iPatch].size(), vector::zero);
399  vectorField v5(mesh.boundary()[iPatch].size(), vector::zero);
400 
401  v6 = mesh.boundary()[iPatch].Cn();
402  if (processorPatch_[iPatch])
403  {
404  v5 = refCast<const processorFvPatch>(mesh.boundary()[iPatch]).
405  procPolyPatch().neighbFaceCellCentres();
406  }
407  else
408  {
409  v5 = v6 + 2.0*
410  (
411  mesh.boundary()[iPatch].Cf()
412  -
413  v6
414  );
415  }
416 
417  label gFaceId = -1;
418  forAll(bqf_[iPatch], k)
419  {
420  facei = bqf_[iPatch][k];
421  gFaceId = bgfid_[iPatch][facei];
422 
423  p1 = faces[gFaceId][0];
424  p2 = faces[gFaceId][1];
425  p3 = faces[gFaceId][2];
426  p4 = faces[gFaceId][3];
427  //p5 - is nei and stored in v5
428  //p6 - is own and stored in v6
429 
430  bvq_[iPatch][k] = (points[p3] - points[p1]) &
431  (
432  (points[p4] - points[p2]) ^ (v6[facei] - v5[facei])
433  );
434  bvq_[iPatch][k] *= OneBySix;
435 
436  /* coefficient for X */
437  baqx_[iPatch][k].resize(6);
438  baqx_[iPatch][k][0] = OneBySix*((v5[facei].y() - v6[facei].y())*(points[p2].z() - points[p4].z())
439  - (v5[facei].z() - v6[facei].z())*(points[p2].y() - points[p4].y()));
440  baqx_[iPatch][k][1] = OneBySix*((v5[facei].y() - v6[facei].y())*(points[p3].z() - points[p1].z())
441  - (v5[facei].z() - v6[facei].z())*(points[p3].y() - points[p1].y()));
442  baqx_[iPatch][k][5] = OneBySix*((points[p1].y() - points[p3].y())*(points[p2].z() - points[p4].z())
443  - (points[p1].z() - points[p3].z())*(points[p2].y() - points[p4].y()));
444 
445  baqx_[iPatch][k][2] = -baqx_[iPatch][k][0];
446  baqx_[iPatch][k][3] = -baqx_[iPatch][k][1];
447  baqx_[iPatch][k][4] = -baqx_[iPatch][k][5];
448 
449  /* coefficient for Y */
450  baqy_[iPatch][k].resize(6);
451  baqy_[iPatch][k][0] = OneBySix*((v5[facei].z() - v6[facei].z())*(points[p2].x() - points[p4].x())
452  - (v5[facei].x() - v6[facei].x())*(points[p2].z() - points[p4].z()));
453  baqy_[iPatch][k][1] = OneBySix*((v5[facei].z() - v6[facei].z())*(points[p3].x() - points[p1].x())
454  - (v5[facei].x() - v6[facei].x())*(points[p3].z() - points[p1].z()));
455  baqy_[iPatch][k][5] = OneBySix*((points[p1].z() - points[p3].z())*(points[p2].x() - points[p4].x())
456  - (points[p1].x() - points[p3].x())*(points[p2].z() - points[p4].z()));
457 
458  baqy_[iPatch][k][2] = -baqy_[iPatch][k][0];
459  baqy_[iPatch][k][3] = -baqy_[iPatch][k][1];
460  baqy_[iPatch][k][4] = -baqy_[iPatch][k][5];
461 
462  /* coefficient for Z */
463  baqz_[iPatch][k].resize(6);
464  baqz_[iPatch][k][0] = OneBySix*((v5[facei].x() - v6[facei].x())*(points[p2].y() - points[p4].y())
465  - (v5[facei].y() - v6[facei].y())*(points[p2].x() - points[p4].x()));
466  baqz_[iPatch][k][1] = OneBySix*((v5[facei].x() - v6[facei].x())*(points[p3].y() - points[p1].y())
467  - (v5[facei].y() - v6[facei].y())*(points[p3].x() - points[p1].x()));
468  baqz_[iPatch][k][5] = OneBySix*((points[p1].x() - points[p3].x())*(points[p2].y() - points[p4].y())
469  - (points[p1].y() - points[p3].y())*(points[p2].x() - points[p4].x()));
470 
471  baqz_[iPatch][k][2] = -baqz_[iPatch][k][0];
472  baqz_[iPatch][k][3] = -baqz_[iPatch][k][1];
473  baqz_[iPatch][k][4] = -baqz_[iPatch][k][5];
474  }
475  }
476 }
477 
478 Foam::fvsc::GaussVolPointBase3D::~GaussVolPointBase3D()
479 {
480 }
481 
482 #define VEC_CMPT(V,CMPT)\
483  V[CMPT]
484 
485 #define SCA_CMPT(V,CMPT)\
486  V
487 
488 #define dfdxif(vf,pf,dfdxfield,fi,vi,ai,icmpt,ocmpt,iop,oop) \
489 { \
490  label celll = -1; \
491  label facei = -1; \
492  const label iown = (ai.size() > 0) ? ai[0].size() - 1 : 0; \
493  const label inei = (ai.size() > 0) ? iown - 1 : 0; \
494  scalar dfdxface = 0.0; \
495  forAll(fi, i) \
496  { \
497  facei = fi[i]; \
498  celll = vf.mesh().neighbour()[facei]; \
499  dfdxface = \
500  iop(vf.primitiveField()[celll],icmpt) * ai[i][inei]; \
501  celll = vf.mesh().owner()[facei]; \
502  dfdxface += \
503  iop(vf.primitiveField()[celll],icmpt) * ai[i][iown]; \
504  \
505  forAll(faces[facei], k) \
506  { \
507  dfdxface += \
508  iop(pf[faces[facei][k]],icmpt) * ai[i][k]; \
509  } \
510  oop(dfdxfield.primitiveFieldRef()[facei],ocmpt) \
511  += (dfdxface / vi[i]); \
512  } \
513 }
514 
515 #define dfdxbf(vf,pf,patchi,dfdxfield,bfi,bvfi,bai,icmpt,ocmpt,iop,oop) \
516 { \
517  label ifacei = -1; \
518  const label iown = (bai[patchi].size() > 0) ? bai[patchi][0].size() - 1 : 0; \
519  const label inei = (bai[patchi].size() > 0) ? iown - 1 : 0; \
520  label gFaceId = -1; \
521  scalar dfdxface = 0.0; \
522  forAll(bfi[patchi], i) \
523  { \
524  ifacei = bfi[patchi][i]; \
525  gFaceId = bgfid_[patchi][ifacei]; \
526  dfdxface = \
527  iop(psin[patchi][ifacei],icmpt) * bai[patchi][i][inei]; \
528  dfdxface += \
529  iop(psio[ifacei],icmpt) * bai[patchi][i][iown]; \
530  forAll(faces[gFaceId], k) \
531  { \
532  dfdxface+= \
533  iop(pf[faces[gFaceId][k]],icmpt) * \
534  bai[patchi][i][k]; \
535  } \
536  oop(dfdxfield.boundaryFieldRef()[patchi][ifacei],ocmpt) += \
537  dfdxface / bvfi[patchi][i]; \
538  } \
539 }
540 
541 
542 
543 void Foam::fvsc::GaussVolPointBase3D::calcDivfIF
544 (
545  const volVectorField& vf,
546  const pointVectorField& pf,
547  const faceList& faces,
548  surfaceScalarField& divf,
549  const surfaceScalarField& dfdn
550 )
551 {
552  //calculate at quad faces
553  dfdxif(vf,pf,divf,qf_,vq_,aqx_,0,0,VEC_CMPT,SCA_CMPT) //X
554  dfdxif(vf,pf,divf,qf_,vq_,aqy_,1,0,VEC_CMPT,SCA_CMPT) //Y
555  dfdxif(vf,pf,divf,qf_,vq_,aqz_,2,0,VEC_CMPT,SCA_CMPT) //Z
556 
557  //calculate at triangular faces
558  dfdxif(vf,pf,divf,tf_,vt_,atx_,0,0,VEC_CMPT,SCA_CMPT) //X
559  dfdxif(vf,pf,divf,tf_,vt_,aty_,1,0,VEC_CMPT,SCA_CMPT) //Y
560  dfdxif(vf,pf,divf,tf_,vt_,atz_,2,0,VEC_CMPT,SCA_CMPT) //Z
561 
562  //other faces
563  {
564  label facei = -1;
565  forAll(of_, i)
566  {
567  facei = of_[i];
568  divf.primitiveFieldRef()[facei] =
569  dfdn.primitiveField()[facei];
570  }
571  }
572 }
573 
574 void Foam::fvsc::GaussVolPointBase3D::calcDivfBF
575 (
576  const volVectorField& vf,
577  const pointVectorField& pf,
578  const faceList& faces,
579  surfaceScalarField& divf,
580  const surfaceScalarField& dfdn
581 )
582 {
583  List<List<vector> > psin (vf.boundaryField().size());
584  forAll(psin, iPatch)
585  {
586  if (processorPatch_[iPatch])
587  {
588  psin[iPatch] = refCast<const processorFvPatchField<vector> >
589  (vf.boundaryField()[iPatch]).patchNeighbourField();
590  }
591  else
592  {
593  psin[iPatch] = vf.boundaryField()[iPatch] +
594  vf.boundaryField()[iPatch].snGrad()
595  *
596  bmvON_[iPatch]*0.5;
597  }
598 
599  vectorField psio (vf.boundaryField()[iPatch].patchInternalField());
600 
601  // quad faces
602  dfdxbf(vf,pf,iPatch,divf,bqf_,bvq_,baqx_,0,0,VEC_CMPT,SCA_CMPT) //X
603  dfdxbf(vf,pf,iPatch,divf,bqf_,bvq_,baqy_,1,0,VEC_CMPT,SCA_CMPT) //Y
604  dfdxbf(vf,pf,iPatch,divf,bqf_,bvq_,baqz_,2,0,VEC_CMPT,SCA_CMPT) //Z
605 
606  // tri faces
607  dfdxbf(vf,pf,iPatch,divf,btf_,bvt_,batx_,0,0,VEC_CMPT,SCA_CMPT) //X
608  dfdxbf(vf,pf,iPatch,divf,btf_,bvt_,baty_,1,0,VEC_CMPT,SCA_CMPT) //Y
609  dfdxbf(vf,pf,iPatch,divf,btf_,bvt_,batz_,2,0,VEC_CMPT,SCA_CMPT) //Z
610 
611  //for other faces - apply surface normal derivative
612  {
613  label facei = -1;
614  forAll(bof_[iPatch], i)
615  {
616  facei = bof_[iPatch][i];
617  divf.boundaryFieldRef()[iPatch][facei] =
618  dfdn.boundaryField()[iPatch][facei];
619  }
620  }
621  }
622 }
623 
624 void Foam::fvsc::GaussVolPointBase3D::calcDivfIF
625 (
626  const volTensorField& tf,
627  const pointTensorField& pf,
628  const faceList& faces,
629  surfaceVectorField& divf,
630  const surfaceVectorField& dfdn
631 )
632 {
633  //calculate at quandrangle faces
634  //X
635  dfdxif(tf,pf,divf,qf_,vq_,aqx_,0,0,VEC_CMPT,VEC_CMPT) // dT_xx / dx
636  dfdxif(tf,pf,divf,qf_,vq_,aqy_,3,0,VEC_CMPT,VEC_CMPT) // dT_yx / dy
637  dfdxif(tf,pf,divf,qf_,vq_,aqz_,6,0,VEC_CMPT,VEC_CMPT) // dT_zx / dz
638  //Y
639  dfdxif(tf,pf,divf,qf_,vq_,aqx_,1,1,VEC_CMPT,VEC_CMPT) // dT_xy / dx
640  dfdxif(tf,pf,divf,qf_,vq_,aqy_,4,1,VEC_CMPT,VEC_CMPT) // dT_yy / dy
641  dfdxif(tf,pf,divf,qf_,vq_,aqz_,7,1,VEC_CMPT,VEC_CMPT) // dT_zy / dz
642  //Z
643  dfdxif(tf,pf,divf,qf_,vq_,aqx_,2,2,VEC_CMPT,VEC_CMPT) // dT_xz / dx
644  dfdxif(tf,pf,divf,qf_,vq_,aqy_,5,2,VEC_CMPT,VEC_CMPT) // dT_yz / dy
645  dfdxif(tf,pf,divf,qf_,vq_,aqz_,8,2,VEC_CMPT,VEC_CMPT) // dT_zz / dz
646 
647  //calculate at triangular faces
648  //X
649  dfdxif(tf,pf,divf,tf_,vt_,atx_,0,0,VEC_CMPT,VEC_CMPT) // dT_xx / dx
650  dfdxif(tf,pf,divf,tf_,vt_,aty_,3,0,VEC_CMPT,VEC_CMPT) // dT_yx / dy
651  dfdxif(tf,pf,divf,tf_,vt_,atz_,6,0,VEC_CMPT,VEC_CMPT) // dT_zx / dz
652  //Y
653  dfdxif(tf,pf,divf,tf_,vt_,atx_,1,1,VEC_CMPT,VEC_CMPT) // dT_xy / dx
654  dfdxif(tf,pf,divf,tf_,vt_,aty_,4,1,VEC_CMPT,VEC_CMPT) // dT_yy / dy
655  dfdxif(tf,pf,divf,tf_,vt_,atz_,7,1,VEC_CMPT,VEC_CMPT) // dT_zy / dz
656  //Z
657  dfdxif(tf,pf,divf,tf_,vt_,atx_,2,2,VEC_CMPT,VEC_CMPT) // dT_xz / dx
658  dfdxif(tf,pf,divf,tf_,vt_,aty_,5,2,VEC_CMPT,VEC_CMPT) // dT_yz / dy
659  dfdxif(tf,pf,divf,tf_,vt_,atz_,8,2,VEC_CMPT,VEC_CMPT) // dT_zz / dz
660 
661  //
662  //other faces
663  {
664  label facei = -1;
665  forAll(of_, i)
666  {
667  facei = of_[i];
668  divf.primitiveFieldRef()[facei] =
669  dfdn.primitiveField()[facei];
670  }
671  }
672 }
673 
674 void Foam::fvsc::GaussVolPointBase3D::calcDivfBF
675 (
676  const volTensorField& tf,
677  const pointTensorField& pf,
678  const faceList& faces,
679  surfaceVectorField& divf,
680  const surfaceVectorField& dfdn
681 )
682 {
683  List<List<tensor> > psin (tf.boundaryField().size());
684  forAll(psin, iPatch)
685  {
686  if (processorPatch_[iPatch])
687  {
688  psin[iPatch] = refCast<const processorFvPatchField<tensor> >
689  (tf.boundaryField()[iPatch]).patchNeighbourField();
690  }
691  else
692  {
693  psin[iPatch] = tf.boundaryField()[iPatch] +
694  tf.boundaryField()[iPatch].snGrad()
695  *
696  bmvON_[iPatch]*0.5;
697  }
698 
699  tensorField psio (tf.boundaryField()[iPatch].patchInternalField());
700 
701  // quad faces
702  dfdxbf(tf,pf,iPatch,divf,bqf_,bvq_,baqx_,0,0,VEC_CMPT,VEC_CMPT) //d/dx
703  dfdxbf(tf,pf,iPatch,divf,bqf_,bvq_,baqy_,3,0,VEC_CMPT,VEC_CMPT) //d/dy
704  dfdxbf(tf,pf,iPatch,divf,bqf_,bvq_,baqz_,6,0,VEC_CMPT,VEC_CMPT) //d/dz
705 
706  dfdxbf(tf,pf,iPatch,divf,bqf_,bvq_,baqx_,1,1,VEC_CMPT,VEC_CMPT) //d/dx
707  dfdxbf(tf,pf,iPatch,divf,bqf_,bvq_,baqy_,4,1,VEC_CMPT,VEC_CMPT) //d/dy
708  dfdxbf(tf,pf,iPatch,divf,bqf_,bvq_,baqz_,7,1,VEC_CMPT,VEC_CMPT) //d/dz
709 
710  dfdxbf(tf,pf,iPatch,divf,bqf_,bvq_,baqx_,2,2,VEC_CMPT,VEC_CMPT) //d/dx
711  dfdxbf(tf,pf,iPatch,divf,bqf_,bvq_,baqy_,5,2,VEC_CMPT,VEC_CMPT) //d/dy
712  dfdxbf(tf,pf,iPatch,divf,bqf_,bvq_,baqz_,8,2,VEC_CMPT,VEC_CMPT) //d/dz
713 
714  // tri faces
715  dfdxbf(tf,pf,iPatch,divf,btf_,bvt_,batx_,0,0,VEC_CMPT,VEC_CMPT) //d/dx
716  dfdxbf(tf,pf,iPatch,divf,btf_,bvt_,baty_,3,0,VEC_CMPT,VEC_CMPT) //d/dy
717  dfdxbf(tf,pf,iPatch,divf,btf_,bvt_,batz_,6,0,VEC_CMPT,VEC_CMPT) //d/dz
718 
719  dfdxbf(tf,pf,iPatch,divf,btf_,bvt_,batx_,1,1,VEC_CMPT,VEC_CMPT) //d/dx
720  dfdxbf(tf,pf,iPatch,divf,btf_,bvt_,baty_,4,1,VEC_CMPT,VEC_CMPT) //d/dy
721  dfdxbf(tf,pf,iPatch,divf,btf_,bvt_,batz_,7,1,VEC_CMPT,VEC_CMPT) //d/dz
722 
723  dfdxbf(tf,pf,iPatch,divf,btf_,bvt_,batx_,2,2,VEC_CMPT,VEC_CMPT) //d/dx
724  dfdxbf(tf,pf,iPatch,divf,btf_,bvt_,baty_,5,2,VEC_CMPT,VEC_CMPT) //d/dy
725  dfdxbf(tf,pf,iPatch,divf,btf_,bvt_,batz_,8,2,VEC_CMPT,VEC_CMPT) //d/dz
726 
727  //for other faces - apply surface normal derivative
728  {
729  label facei = -1;
730  forAll(bof_[iPatch], i)
731  {
732  facei = bof_[iPatch][i];
733  divf.boundaryFieldRef()[iPatch][facei] =
734  dfdn.boundaryField()[iPatch][facei];
735  }
736  }
737  }
738 }
739 
740 void Foam::fvsc::GaussVolPointBase3D::calcGradfIF
741 (
742  const volScalarField& sf,
743  const pointScalarField& pf,
744  const faceList& faces,
745  surfaceVectorField& gradf,
746  const surfaceVectorField& dfdn
747 )
748 {
749  //quad faces
750  dfdxif(sf,pf,gradf,qf_,vq_,aqx_,0,0,SCA_CMPT,VEC_CMPT) //X
751  dfdxif(sf,pf,gradf,qf_,vq_,aqy_,0,1,SCA_CMPT,VEC_CMPT) //Y
752  dfdxif(sf,pf,gradf,qf_,vq_,aqz_,0,2,SCA_CMPT,VEC_CMPT) //Z
753 
754  //triangular faces
755  dfdxif(sf,pf,gradf,tf_,vt_,atx_,0,0,SCA_CMPT,VEC_CMPT) //X
756  dfdxif(sf,pf,gradf,tf_,vt_,aty_,0,1,SCA_CMPT,VEC_CMPT) //Y
757  dfdxif(sf,pf,gradf,tf_,vt_,atz_,0,2,SCA_CMPT,VEC_CMPT) //Z
758 
759  //other faces
760  {
761  label facei = -1;
762  forAll(of_, i)
763  {
764  facei = of_[i];
765  gradf.primitiveFieldRef()[facei] =
766  dfdn.primitiveField()[facei];
767  }
768  }
769 }
770 
771 void Foam::fvsc::GaussVolPointBase3D::calcGradfBF
772 (
773  const volScalarField& sf,
774  const pointScalarField& pf,
775  const faceList& faces,
776  surfaceVectorField& gradf,
777  const surfaceVectorField& dfdn
778 )
779 {
780  List<List<scalar> > psin (sf.boundaryField().size());
781  forAll(psin, iPatch)
782  {
783  if (processorPatch_[iPatch])
784  {
785  psin[iPatch] = refCast<const processorFvPatchField<scalar> >
786  (sf.boundaryField()[iPatch]).patchNeighbourField();
787  }
788  else
789  {
790  psin[iPatch] = sf.boundaryField()[iPatch] +
791  sf.boundaryField()[iPatch].snGrad()
792  *
793  bmvON_[iPatch]*0.5;
794  }
795 
796  scalarField psio (sf.boundaryField()[iPatch].patchInternalField());
797 
798  //quad faces
799  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqx_,0,0,SCA_CMPT,VEC_CMPT) //X
800  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqy_,0,1,SCA_CMPT,VEC_CMPT) //Y
801  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqz_,0,2,SCA_CMPT,VEC_CMPT) //Z
802 
803  //tri faces
804  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,batx_,0,0,SCA_CMPT,VEC_CMPT) //X
805  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,baty_,0,1,SCA_CMPT,VEC_CMPT) //Y
806  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,batz_,0,2,SCA_CMPT,VEC_CMPT) //Z
807 
808  //for other faces - apply surface normal derivative
809  {
810  label facei = -1;
811  forAll(bof_[iPatch], i)
812  {
813  facei = bof_[iPatch][i];
814  gradf.boundaryFieldRef()[iPatch][facei] =
815  dfdn.boundaryField()[iPatch][facei];
816  }
817  }
818  }
819 }
820 
821 void Foam::fvsc::GaussVolPointBase3D::calcGradfIF
822 (
823  const volVectorField& vf,
824  const pointVectorField& pf,
825  const faceList& faces,
826  surfaceTensorField& gradf,
827  const surfaceTensorField& dfdn
828 )
829 {
830  //quad faces
831  dfdxif(vf,pf,gradf,qf_,vq_,aqx_,0,0,VEC_CMPT,VEC_CMPT) //X
832  dfdxif(vf,pf,gradf,qf_,vq_,aqx_,1,1,VEC_CMPT,VEC_CMPT) //Y
833  dfdxif(vf,pf,gradf,qf_,vq_,aqx_,2,2,VEC_CMPT,VEC_CMPT) //Z
834 
835  dfdxif(vf,pf,gradf,qf_,vq_,aqy_,0,3,VEC_CMPT,VEC_CMPT) //X
836  dfdxif(vf,pf,gradf,qf_,vq_,aqy_,1,4,VEC_CMPT,VEC_CMPT) //Y
837  dfdxif(vf,pf,gradf,qf_,vq_,aqy_,2,5,VEC_CMPT,VEC_CMPT) //Z
838 
839  dfdxif(vf,pf,gradf,qf_,vq_,aqz_,0,6,VEC_CMPT,VEC_CMPT) //X
840  dfdxif(vf,pf,gradf,qf_,vq_,aqz_,1,7,VEC_CMPT,VEC_CMPT) //Y
841  dfdxif(vf,pf,gradf,qf_,vq_,aqz_,2,8,VEC_CMPT,VEC_CMPT) //Z
842 
843  //triangular faces
844  dfdxif(vf,pf,gradf,tf_,vt_,atx_,0,0,VEC_CMPT,VEC_CMPT) // X
845  dfdxif(vf,pf,gradf,tf_,vt_,aty_,1,1,VEC_CMPT,VEC_CMPT) // Y
846  dfdxif(vf,pf,gradf,tf_,vt_,atz_,2,2,VEC_CMPT,VEC_CMPT) // Z
847 
848  dfdxif(vf,pf,gradf,tf_,vt_,atx_,0,3,VEC_CMPT,VEC_CMPT) // X
849  dfdxif(vf,pf,gradf,tf_,vt_,aty_,1,4,VEC_CMPT,VEC_CMPT) // Y
850  dfdxif(vf,pf,gradf,tf_,vt_,atz_,2,5,VEC_CMPT,VEC_CMPT) // Z
851 
852  dfdxif(vf,pf,gradf,tf_,vt_,atx_,0,6,VEC_CMPT,VEC_CMPT) // X
853  dfdxif(vf,pf,gradf,tf_,vt_,aty_,1,7,VEC_CMPT,VEC_CMPT) // Y
854  dfdxif(vf,pf,gradf,tf_,vt_,atz_,2,8,VEC_CMPT,VEC_CMPT) // Z
855 
856  //other faces
857  {
858  label facei = -1;
859  forAll(of_, i)
860  {
861  facei = of_[i];
862  gradf.primitiveFieldRef()[facei] =
863  dfdn.primitiveField()[facei];
864  }
865  }
866 }
867 
868 void Foam::fvsc::GaussVolPointBase3D::calcGradfBF
869 (
870  const volVectorField& sf,
871  const pointVectorField& pf,
872  const faceList& faces,
873  surfaceTensorField& gradf,
874  const surfaceTensorField& dfdn
875 )
876 {
877  List<List<vector> > psin (sf.boundaryField().size());
878  forAll(psin, iPatch)
879  {
880  if (processorPatch_[iPatch])
881  {
882  psin[iPatch] = refCast<const processorFvPatchField<vector> >
883  (sf.boundaryField()[iPatch]).patchNeighbourField();
884  }
885  else
886  {
887  psin[iPatch] = sf.boundaryField()[iPatch] +
888  sf.boundaryField()[iPatch].snGrad()
889  *
890  bmvON_[iPatch]*0.5;
891  }
892 
893  vectorField psio (sf.boundaryField()[iPatch].patchInternalField());
894 
895  //quad faces
896  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqx_,0,0,VEC_CMPT,VEC_CMPT) //X
897  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqx_,1,1,VEC_CMPT,VEC_CMPT) //Y
898  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqx_,2,2,VEC_CMPT,VEC_CMPT) //Z
899 
900  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqy_,0,3,VEC_CMPT,VEC_CMPT) //X
901  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqy_,1,4,VEC_CMPT,VEC_CMPT) //Y
902  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqy_,2,5,VEC_CMPT,VEC_CMPT) //Z
903 
904  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqz_,0,6,VEC_CMPT,VEC_CMPT) //X
905  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqz_,1,7,VEC_CMPT,VEC_CMPT) //Y
906  dfdxbf(sf,pf,iPatch,gradf,bqf_,bvq_,baqz_,2,8,VEC_CMPT,VEC_CMPT) //Z
907 
908  //tri faces
909  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,batx_,0,0,VEC_CMPT,VEC_CMPT) //X
910  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,batx_,1,1,VEC_CMPT,VEC_CMPT) //Y
911  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,batx_,2,2,VEC_CMPT,VEC_CMPT) //Z
912 
913  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,baty_,0,3,VEC_CMPT,VEC_CMPT) //X
914  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,baty_,1,4,VEC_CMPT,VEC_CMPT) //Y
915  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,baty_,2,5,VEC_CMPT,VEC_CMPT) //Z
916 
917  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,batz_,0,6,VEC_CMPT,VEC_CMPT) //X
918  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,batz_,1,7,VEC_CMPT,VEC_CMPT) //Y
919  dfdxbf(sf,pf,iPatch,gradf,btf_,bvt_,batz_,2,8,VEC_CMPT,VEC_CMPT) //Z
920 
921  //for other faces - apply surface normal derivative
922  {
923  label facei = -1;
924  forAll(bof_[iPatch], i)
925  {
926  facei = bof_[iPatch][i];
927  gradf.boundaryFieldRef()[iPatch][facei] =
928  dfdn.boundaryField()[iPatch][facei];
929  }
930  }
931  }
932 }
933 
934 void Foam::fvsc::GaussVolPointBase3D::faceGrad(const volScalarField& sf, surfaceVectorField& gradf)
935 {
936  pointScalarField pF
937  (
938  volPoint_.interpolate
939  (
940  sf
941  )
942  );
943  //
944  const surfaceVectorField& nf = nfRef_();
945  surfaceVectorField dfdn
946  (
947  nf * fvc::snGrad(sf)
948  );
949  //
950  const faceList& faces = sf.mesh().faces();
951  /*
952  *
953  * Calculate grad at internal faces
954  *
955  */
956  calcGradfIF(sf, pF, faces, gradf, dfdn);
957  /*
958  *
959  * Calculate grad at boundary faces
960  *
961  */
962  calcGradfBF(sf, pF, faces, gradf, dfdn);
963 };
964 
965 void Foam::fvsc::GaussVolPointBase3D::faceGrad(const volVectorField& vf, surfaceTensorField& gradf)
966 {
967  pointVectorField pF
968  (
969  volPoint_.interpolate
970  (
971  vf
972  )
973  );
974  //
975  const surfaceVectorField& nf = nfRef_();
976  surfaceTensorField dfdn
977  (
978  nf * fvc::snGrad(vf)
979  );
980  //
981  const faceList& faces = vf.mesh().faces();
982  /*
983  *
984  * Calculate grad at internal faces
985  *
986  */
987  calcGradfIF(vf, pF, faces, gradf, dfdn);
988  /*
989  *
990  * Calculate grad at boundary faces
991  *
992  */
993  calcGradfBF(vf, pF, faces, gradf, dfdn);
994 };
995 
996 void Foam::fvsc::GaussVolPointBase3D::faceDiv(const volVectorField& vf, surfaceScalarField& divf)
997 {
998  pointVectorField pF
999  (
1000  volPoint_.interpolate
1001  (
1002  vf
1003  )
1004  );
1005 
1006  surfaceScalarField divfn (nfRef_() & fvc::snGrad(vf));
1007  const faceList& faces = vf.mesh().faces();
1008  /*
1009  *
1010  * Calculate grad at internal faces
1011  *
1012  */
1013  calcDivfIF(vf, pF, faces, divf, divfn);
1014  /*
1015  *
1016  * Calculate grad at boundary faces
1017  *
1018  */
1019  calcDivfBF(vf, pF, faces, divf, divfn);
1020 };
1021 
1022 void Foam::fvsc::GaussVolPointBase3D::faceDiv(const volTensorField& tf, surfaceVectorField& divf)
1023 {
1024  pointTensorField pF
1025  (
1026  volPoint_.interpolate
1027  (
1028  tf
1029  )
1030  );
1031 
1032  surfaceVectorField divfn (nfRef_() & fvc::snGrad(tf));
1033  const faceList& faces = tf.mesh().faces();
1034  /*
1035  *
1036  * Calculate grad at internal faces
1037  *
1038  */
1039  calcDivfIF(tf, pF, faces, divf, divfn);
1040  /*
1041  *
1042  * Calculate grad at boundary faces
1043  *
1044  */
1045  calcDivfBF(tf, pF, faces, divf, divfn);
1046 }
1047 
1048 //
1049 //END-OF-FILE
1050 //
1051 
1052 
Foam::fvsc::GaussVolPointBase3D::faceDiv
void faceDiv(const volVectorField &f, surfaceScalarField &divf)
Definition: GaussVolPointBase3D.C:989
pf
pf
Definition: updateFields.H:21
Foam::fvsc::GaussVolPointBase3D::faceGrad
void faceGrad(const volScalarField &f, surfaceVectorField &gradf)
Definition: GaussVolPointBase3D.C:927
VEC_CMPT
#define VEC_CMPT(V, CMPT)
Definition: GaussVolPointBase3D.C:475
Foam::fvsc::GaussVolPointBase3D::triCalcWeights
void triCalcWeights(const fvMesh &m)
Calculate weights for triangles.
Definition: GaussVolPointBase3D.C:154
Foam::faceList
List< face > faceList
Definition: GaussVolPointBase3D.H:42
dfdxif
#define dfdxif(vf, pf, dfdxfield, fi, vi, ai, icmpt, ocmpt, iop, oop)
Definition: GaussVolPointBase3D.C:481
Foam::fvsc::GaussVolPointBase3D::calcDivfBF
void calcDivfBF(const volVectorField &sf, const pointVectorField &pf, const faceList &faces, surfaceScalarField &divf, const surfaceScalarField &dfdn)
Definition: GaussVolPointBase3D.C:568
forAll
forAll(Y, i)
Definition: QGDYEqn.H:36
Foam::fvsc::GaussVolPointBase3D::calcGradfBF
void calcGradfBF(const volScalarField &sf, const pointScalarField &pf, const faceList &faces, surfaceVectorField &gradf, const surfaceVectorField &dfdn)
Definition: GaussVolPointBase3D.C:765
Foam::fvsc::GaussVolPointBase3D::GaussVolPointBase3D
GaussVolPointBase3D(const fvMesh &mesh)
Definition: GaussVolPointBase3D.C:34
Foam::fvsc::GaussVolPointBase3D::calcGradfIF
void calcGradfIF(const volScalarField &sf, const pointScalarField &pf, const faceList &faces, surfaceVectorField &gradf, const surfaceVectorField &dfdn)
Definition: GaussVolPointBase3D.C:734
SCA_CMPT
#define SCA_CMPT(V, CMPT)
Definition: GaussVolPointBase3D.C:478
Foam::fvsc::GaussVolPointBase3D::quaCalcWeights
void quaCalcWeights(const fvMesh &m)
Calcualte weights for quads.
Definition: GaussVolPointBase3D.C:313
dfdxbf
#define dfdxbf(vf, pf, patchi, dfdxfield, bfi, bvfi, bai, icmpt, ocmpt, iop, oop)
Definition: GaussVolPointBase3D.C:508
Foam::fvsc::GaussVolPointBase3D::calcDivfIF
void calcDivfIF(const volVectorField &sf, const pointVectorField &pf, const faceList &faces, surfaceScalarField &divf, const surfaceScalarField &dfdn)
Definition: GaussVolPointBase3D.C:537

Generated by   doxygen 1.8.5