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

Generated by   doxygen 1.8.5