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"
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),
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())
70 if (mesh.nGeometricD() == 2)
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);
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());
86 label ip1=-1, ic2=-1, ip3=-1, ic4=-1;
88 forAll(mesh.geometricD(), iDir)
90 if (mesh.geometricD()[iDir] < 1)
97 e1_ = vector(0, 1, 0);
98 e2_ = vector(0, 0, 1);
105 e1_ = vector(1, 0, 0);
106 e2_ = vector(0, 0, 1);
113 e1_ = vector(1, 0, 0);
114 e2_ = vector(0, 1, 0);
120 forAll(mesh.neighbour(), iFace)
122 ic2 = mesh.neighbour()[iFace];
123 ic4 = mesh.owner()[iFace];
125 const face& f = mesh.faces()[iFace];
126 #warning "Raise error if number of points on face is not equal to 4"
129 if (mesh.points()[f[ip]][
ie3_] >= mesh.C()[ic2][
ie3_])
137 if (mesh.points()[f[ip]][
ie3_] >= mesh.C()[ic2][
ie3_])
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]);
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_;
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];
170 forAll(mesh.boundary(), iPatch)
172 label patchId = iPatch;
174 !isA<emptyFvPatch>(mesh.boundary()[patchId])
176 !isA<wedgeFvPatch>(mesh.boundary()[patchId])
181 isA<coupledFvPatch>(mesh.boundary()[patchId])
183 !isA<processorFvPatch>(mesh.boundary()[patchId])
190 if (isA<processorFvPatch>(mesh.boundary()[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());
217 ic4e_[patchId] = mesh.boundary()[patchId].faceCells();
220 if (isA<processorFvPatch>(mesh.boundary()[patchId]))
222 const processorPolyPatch & procPolyPatch =
223 refCast<const processorFvPatch>(mesh.boundary()[patchId]).procPolyPatch();
227 v42[iFace] = procPolyPatch.neighbFaceCellCentres()[iFace] -
228 mesh.C()[
ic4e_[patchId][iFace]];
235 v42[iFace] = 2.0*(mesh.boundary()[patchId].Cf()[iFace] -
236 mesh.C()[
ic4e_[patchId][iFace]]);
240 forAll(mesh.boundary()[patchId], iFace)
242 label ip1=-1, ip3=-1, ic4=
ic4e_[patchId][iFace];
243 label globalFaceID = mesh.boundary()[patchId].start() + iFace;
245 const face& f = mesh.faces()[globalFaceID];
248 if (mesh.points()[f[ip]][
ie3_] >= mesh.C()[ic4][
ie3_])
256 if (mesh.points()[f[ip]][
ie3_] >= mesh.C()[ic4][
ie3_])
266 ip3e_[patchId][iFace] = ip3;
267 ip1e_[patchId][iFace] = ip1;
268 v13[iFace] = mesh.points()[ip3] - mesh.points()[ip1];
270 mv42e_[patchId][iFace] = mag(v42[iFace]);
271 mv13e_[patchId][iFace] = mag(v13[iFace]);
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_;
278 dene[patchId][iFace] =
279 sina2e[patchId][iFace]*cosa1e[patchId][iFace]
281 sina1e[patchId][iFace]*cosa2e[patchId][iFace];
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];
301 scalar dfdn = 0.0, dfdt = 0.0;
303 if (f.mesh().nGeometricD() == 2)
313 forAll(f.mesh().neighbour(), iFace)
315 dfdn = (f[ic2_[iFace]] - f[ic4_[iFace]]) / mv42_[iFace];
319 dfdt = (pF[ip3_[iFace]] - pF[ip1_[iFace]]) / mv13_[iFace];
322 gradf[iFace][ie1_] = (dfdn*c1_[iFace] - dfdt*c2_[iFace]);
323 gradf[iFace][ie2_] = (dfdt*c3_[iFace] - dfdn*c4_[iFace]);
326 gradf[iFace][ie3_] = 0.0;
329 List<List<scalar> > psi2(f.boundaryField().size());
331 forAll(ordinaryPatches_, iPatch)
333 label patchId = ordinaryPatches_[iPatch];
334 if (processorPatch_[iPatch])
336 psi2[patchId] = refCast<const processorFvPatchField<scalar> >
337 (f.boundaryField()[patchId]).patchNeighbourField();
341 psi2[patchId] = f.boundaryField()[patchId] +
342 f.boundaryField()[patchId].snGrad()
347 forAll(f.boundaryField()[patchId], iFace)
350 dfdn = (psi2[patchId][iFace] - f[ic4e_[patchId][iFace]]) / mv42e_[patchId][iFace];
353 dfdt = (pF[ip3e_[patchId][iFace]] - pF[ip1e_[patchId][iFace]]) / mv13e_[patchId][iFace];
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;
363 #warning "Raise error if called not for 2D case"
369 scalar df1dn = 0.0, df1dt = 0.0,
370 df2dn = 0.0, df2dt = 0.0;
372 if (f.mesh().nGeometricD() == 2)
382 forAll(f.mesh().neighbour(), iFace)
384 df1dn = (f[ic2_[iFace]][ie1_] - f[ic4_[iFace]][ie1_]) / mv42_[iFace];
385 df2dn = (f[ic2_[iFace]][ie2_] - f[ic4_[iFace]][ie2_]) / mv42_[iFace];
389 df1dt = (pF[ip3_[iFace]][ie1_] - pF[ip1_[iFace]][ie1_]) / mv13_[iFace];
390 df2dt = (pF[ip3_[iFace]][ie2_] - pF[ip1_[iFace]][ie2_]) / mv13_[iFace];
393 divf[iFace] = (df1dn*c1_[iFace] - df1dt*c2_[iFace]) +
394 (df2dt*c3_[iFace] - df2dn*c4_[iFace]);
397 List<List<vector> > psi2(f.boundaryField().size());
399 forAll(ordinaryPatches_, iPatch)
401 label patchId = ordinaryPatches_[iPatch];
402 if (processorPatch_[iPatch])
404 psi2[patchId] = refCast<const processorFvPatchField<vector> >
405 (f.boundaryField()[patchId]).patchNeighbourField();
409 psi2[patchId] = f.boundaryField()[patchId] +
410 f.boundaryField()[patchId].snGrad()
415 forAll(f.boundaryField()[patchId], iFace)
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];
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];
425 divf.boundaryFieldRef()[patchId][iFace] =
426 (df1dn*c1e_[patchId][iFace] - df1dt*c2e_[patchId][iFace])
428 (df2dt*c3e_[patchId][iFace] - df2dn*c4e_[patchId][iFace]);
434 #warning "Raise error if called not for 2D case"
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;
445 label i11 = ie1_*3 + ie1_;
446 label i21 = ie2_*3 + ie1_;
447 label i22 = ie2_*3 + ie2_;
448 label i12 = ie1_*3 + ie2_;
450 if (f.mesh().nGeometricD() == 2)
460 forAll(f.mesh().neighbour(), iFace)
462 df11dn = (f[ic2_[iFace]][i11] - f[ic4_[iFace]][i11]) / mv42_[iFace];
463 df21dn = (f[ic2_[iFace]][i21] - f[ic4_[iFace]][i21]) / mv42_[iFace];
466 df11dt = (pF[ip3_[iFace]][i11] - pF[ip1_[iFace]][i11]) / mv13_[iFace];
467 df21dt = (pF[ip3_[iFace]][i21] - pF[ip1_[iFace]][i21]) / mv13_[iFace];
470 (df11dn*c1_[iFace] - df11dt*c2_[iFace]) +
471 (df21dt*c3_[iFace] - df21dn*c4_[iFace]);
473 df22dn = (f[ic2_[iFace]][i22] - f[ic4_[iFace]][i22]) / mv42_[iFace];
474 df12dn = (f[ic2_[iFace]][i12] - f[ic4_[iFace]][i12]) / mv42_[iFace];
477 df22dt = (pF[ip3_[iFace]][i22] - pF[ip1_[iFace]][i22]) / mv13_[iFace];
478 df12dt = (pF[ip3_[iFace]][i12] - pF[ip1_[iFace]][i12]) / mv13_[iFace];
481 (df12dn*c1_[iFace] - df12dt*c2_[iFace]) +
482 (df22dt*c3_[iFace] - df22dn*c4_[iFace]);
485 List<List<tensor> > psi2(f.boundaryField().size());
487 forAll(ordinaryPatches_, iPatch)
489 label patchId = ordinaryPatches_[iPatch];
490 if (processorPatch_[iPatch])
492 psi2[patchId] = refCast<const processorFvPatchField<tensor> >
493 (f.boundaryField()[patchId]).patchNeighbourField();
497 psi2[patchId] = f.boundaryField()[patchId] +
498 f.boundaryField()[patchId].snGrad()
503 forAll(f.boundaryField()[patchId], iFace)
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];
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];
513 divf.boundaryFieldRef()[patchId][iFace][ie1_] =
514 (df11dn*c1e_[patchId][iFace] - df11dt*c2e_[patchId][iFace])
516 (df21dt*c3e_[patchId][iFace] - df21dn*c4e_[patchId][iFace]);
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];
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];
526 divf.boundaryFieldRef()[patchId][iFace][ie2_] =
527 (df12dn*c1e_[patchId][iFace] - df12dt*c2e_[patchId][iFace])
529 (df22dt*c3e_[patchId][iFace] - df22dn*c4e_[patchId][iFace]);
535 #warning "Raise error if called not for 2D case"
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &vf, const surfaceScalarField &dir, const word &reconFieldName=word::null)
Interpolate field vf according to direction dir.
List< List< scalar > > c4e_
List< List< label > > ip1e_
List< List< scalar > > c3e_
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)
List< List< scalar > > c2e_
List< List< scalar > > mv42e_
List< List< scalar > > mv13e_
List< List< scalar > > c1e_
List< bool > processorPatch_
GaussVolPointBase2D(const fvMesh &mesh)
void faceGrad(const volScalarField &f, surfaceVectorField &gradf)
List< List< label > > ic4e_
List< List< label > > ip3e_