34 #include "volFields.H"
35 #include "surfaceFields.H"
36 #include "pointFields.H"
37 #include "coupledFvPatch.H"
38 #include "processorFvPatch.H"
39 #include "processorFvPatchFields.H"
40 #include "emptyFvPatch.H"
41 #include "wedgeFvPatch.H"
42 #include "volPointInterpolation.H"
46 c1_(mesh.neighbour().size(), 0.0),
47 c2_(mesh.neighbour().size(), 0.0),
48 c3_(mesh.neighbour().size(), 0.0),
49 c4_(mesh.neighbour().size(), 0.0),
50 mv42_(mesh.neighbour().size(), 1.0),
51 mv13_(mesh.neighbour().size(), 1.0),
52 ip3_(mesh.neighbour().size(), -1),
53 ip1_(mesh.neighbour().size(), -1),
54 ic4_(mesh.neighbour().size(), -1),
55 ic2_(mesh.neighbour().size(), -1),
62 ip3e_(mesh.boundary().size()),
63 ip1e_(mesh.boundary().size()),
64 ic4e_(mesh.boundary().size()),
65 c1e_(mesh.boundary().size()),
66 c2e_(mesh.boundary().size()),
67 c3e_(mesh.boundary().size()),
68 c4e_(mesh.boundary().size()),
69 mv42e_(mesh.boundary().size()),
70 mv13e_(mesh.boundary().size())
72 if (mesh.nGeometricD() == 2)
74 List<scalar> cosa1(mesh.neighbour().size(), 0.0);
75 List<scalar> cosa2(mesh.neighbour().size(), 0.0);
76 List<scalar> sina1(mesh.neighbour().size(), 0.0);
77 List<scalar> sina2(mesh.neighbour().size(), 0.0);
78 List<scalar> den(mesh.neighbour().size(), 1.0);
79 List<vector> v42(mesh.neighbour().size(), vector::one);
80 List<vector> v13(mesh.neighbour().size(), vector::one);
82 List<List<scalar> > cosa1e(mesh.boundary().size());
83 List<List<scalar> > cosa2e(mesh.boundary().size());
84 List<List<scalar> > sina1e(mesh.boundary().size());
85 List<List<scalar> > sina2e(mesh.boundary().size());
86 List<List<scalar> > dene(mesh.boundary().size());
88 label ip1=-1, ic2=-1, ip3=-1, ic4=-1;
90 forAll(mesh.geometricD(), iDir)
92 if (mesh.geometricD()[iDir] < 1)
99 e1_ = vector(0, 1, 0);
100 e2_ = vector(0, 0, 1);
107 e1_ = vector(1, 0, 0);
108 e2_ = vector(0, 0, 1);
115 e1_ = vector(1, 0, 0);
116 e2_ = vector(0, 1, 0);
122 forAll(mesh.neighbour(), iFace)
124 ic2 = mesh.neighbour()[iFace];
125 ic4 = mesh.owner()[iFace];
127 const face& f = mesh.faces()[iFace];
128 #warning "Raise error if number of points on face is not equal to 4"
131 if (mesh.points()[f[ip]][
ie3_] >= mesh.C()[ic2][
ie3_])
139 if (mesh.points()[f[ip]][
ie3_] >= mesh.C()[ic2][
ie3_])
154 v42[iFace] = mesh.C()[ic2] - mesh.C()[ic4];
155 v13[iFace] = mesh.points()[ip3] - mesh.points()[ip1];
156 mv42_[iFace] = mag(v42[iFace]);
157 mv13_[iFace] = mag(v13[iFace]);
159 cosa1[iFace] = (v42[iFace]/
mv42_[iFace]) &
e1_;
160 cosa2[iFace] = (v13[iFace]/
mv13_[iFace]) &
e1_;
161 sina1[iFace] = (v42[iFace]/
mv42_[iFace]) &
e2_;
162 sina2[iFace] = (v13[iFace]/
mv13_[iFace]) &
e2_;
164 den[iFace] = sina2[iFace]*cosa1[iFace] - sina1[iFace]*cosa2[iFace];
165 c1_[iFace] = sina2[iFace] / den[iFace];
166 c2_[iFace] = sina1[iFace] / den[iFace];
167 c3_[iFace] = cosa1[iFace] / den[iFace];
168 c4_[iFace] = cosa2[iFace] / den[iFace];
172 forAll(mesh.boundary(), iPatch)
174 label patchId = iPatch;
176 !isA<emptyFvPatch>(mesh.boundary()[patchId])
178 !isA<wedgeFvPatch>(mesh.boundary()[patchId])
183 isA<coupledFvPatch>(mesh.boundary()[patchId])
185 !isA<processorFvPatch>(mesh.boundary()[patchId])
192 if (isA<processorFvPatch>(mesh.boundary()[patchId]))
202 List<vector> v42(mesh.boundary()[patchId].size(), vector::one);
203 List<vector> v13(mesh.boundary()[patchId].size(), vector::one);
204 ic4e_[patchId].resize(v42.size());
205 mv42e_[patchId].resize(v42.size());
206 mv13e_[patchId].resize(v42.size());
207 ip3e_[patchId].resize(v42.size());
208 ip1e_[patchId].resize(v42.size());
209 sina1e[patchId].resize(v42.size());
210 sina2e[patchId].resize(v42.size());
211 cosa1e[patchId].resize(v42.size());
212 cosa2e[patchId].resize(v42.size());
213 dene[patchId].resize(v42.size());
214 c1e_[patchId].resize(v42.size());
215 c2e_[patchId].resize(v42.size());
216 c3e_[patchId].resize(v42.size());
217 c4e_[patchId].resize(v42.size());
219 ic4e_[patchId] = mesh.boundary()[patchId].faceCells();
222 if (isA<processorFvPatch>(mesh.boundary()[patchId]))
224 const processorPolyPatch & procPolyPatch =
225 refCast<const processorFvPatch>(mesh.boundary()[patchId]).procPolyPatch();
229 v42[iFace] = procPolyPatch.neighbFaceCellCentres()[iFace] -
230 mesh.C()[
ic4e_[patchId][iFace]];
237 v42[iFace] = 2.0*(mesh.boundary()[patchId].Cf()[iFace] -
238 mesh.C()[
ic4e_[patchId][iFace]]);
242 forAll(mesh.boundary()[patchId], iFace)
244 label ip1=-1, ip3=-1, ic4=
ic4e_[patchId][iFace];
245 label globalFaceID = mesh.boundary()[patchId].start() + iFace;
247 const face& f = mesh.faces()[globalFaceID];
250 if (mesh.points()[f[ip]][
ie3_] >= mesh.C()[ic4][
ie3_])
258 if (mesh.points()[f[ip]][
ie3_] >= mesh.C()[ic4][
ie3_])
268 ip3e_[patchId][iFace] = ip3;
269 ip1e_[patchId][iFace] = ip1;
270 v13[iFace] = mesh.points()[ip3] - mesh.points()[ip1];
272 mv42e_[patchId][iFace] = mag(v42[iFace]);
273 mv13e_[patchId][iFace] = mag(v13[iFace]);
275 cosa1e[patchId][iFace] = (v42[iFace]/
mv42e_[patchId][iFace]) &
e1_;
276 cosa2e[patchId][iFace] = (v13[iFace]/
mv13e_[patchId][iFace]) &
e1_;
277 sina1e[patchId][iFace] = (v42[iFace]/
mv42e_[patchId][iFace]) &
e2_;
278 sina2e[patchId][iFace] = (v13[iFace]/
mv13e_[patchId][iFace]) &
e2_;
280 dene[patchId][iFace] =
281 sina2e[patchId][iFace]*cosa1e[patchId][iFace]
283 sina1e[patchId][iFace]*cosa2e[patchId][iFace];
285 c1e_[patchId][iFace] = sina2e[patchId][iFace] / dene[patchId][iFace];
286 c2e_[patchId][iFace] = sina1e[patchId][iFace] / dene[patchId][iFace];
287 c3e_[patchId][iFace] = cosa1e[patchId][iFace] / dene[patchId][iFace];
288 c4e_[patchId][iFace] = cosa2e[patchId][iFace] / dene[patchId][iFace];
303 scalar dfdn = 0.0, dfdt = 0.0;
305 if (f.mesh().nGeometricD() == 2)
315 forAll(f.mesh().neighbour(), iFace)
317 dfdn = (f[ic2_[iFace]] - f[ic4_[iFace]]) / mv42_[iFace];
321 dfdt = (pF[ip3_[iFace]] - pF[ip1_[iFace]]) / mv13_[iFace];
324 gradf[iFace][ie1_] = (dfdn*c1_[iFace] - dfdt*c2_[iFace]);
325 gradf[iFace][ie2_] = (dfdt*c3_[iFace] - dfdn*c4_[iFace]);
328 gradf[iFace][ie3_] = 0.0;
331 List<List<scalar> > psi2(f.boundaryField().size());
333 forAll(ordinaryPatches_, iPatch)
335 label patchId = ordinaryPatches_[iPatch];
336 if (processorPatch_[iPatch])
338 psi2[patchId] = refCast<const processorFvPatchField<scalar> >
339 (f.boundaryField()[patchId]).patchNeighbourField();
343 psi2[patchId] = f.boundaryField()[patchId] +
344 f.boundaryField()[patchId].snGrad()
349 forAll(f.boundaryField()[patchId], iFace)
352 dfdn = (psi2[patchId][iFace] - f[ic4e_[patchId][iFace]]) / mv42e_[patchId][iFace];
355 dfdt = (pF[ip3e_[patchId][iFace]] - pF[ip1e_[patchId][iFace]]) / mv13e_[patchId][iFace];
357 gradf.boundaryFieldRef()[patchId][iFace][ie1_] = (dfdn*c1e_[patchId][iFace] - dfdt*c2e_[patchId][iFace]);
358 gradf.boundaryFieldRef()[patchId][iFace][ie2_] = (dfdt*c3e_[patchId][iFace] - dfdn*c4e_[patchId][iFace]);
359 gradf.boundaryFieldRef()[patchId][iFace][ie3_] = 0.0;
365 #warning "Raise error if called not for 2D case"
371 scalar df1dn = 0.0, df1dt = 0.0,
372 df2dn = 0.0, df2dt = 0.0;
374 if (f.mesh().nGeometricD() == 2)
384 forAll(f.mesh().neighbour(), iFace)
386 df1dn = (f[ic2_[iFace]][ie1_] - f[ic4_[iFace]][ie1_]) / mv42_[iFace];
387 df2dn = (f[ic2_[iFace]][ie2_] - f[ic4_[iFace]][ie2_]) / mv42_[iFace];
391 df1dt = (pF[ip3_[iFace]][ie1_] - pF[ip1_[iFace]][ie1_]) / mv13_[iFace];
392 df2dt = (pF[ip3_[iFace]][ie2_] - pF[ip1_[iFace]][ie2_]) / mv13_[iFace];
395 divf[iFace] = (df1dn*c1_[iFace] - df1dt*c2_[iFace]) +
396 (df2dt*c3_[iFace] - df2dn*c4_[iFace]);
399 List<List<vector> > psi2(f.boundaryField().size());
401 forAll(ordinaryPatches_, iPatch)
403 label patchId = ordinaryPatches_[iPatch];
404 if (processorPatch_[iPatch])
406 psi2[patchId] = refCast<const processorFvPatchField<vector> >
407 (f.boundaryField()[patchId]).patchNeighbourField();
411 psi2[patchId] = f.boundaryField()[patchId] +
412 f.boundaryField()[patchId].snGrad()
417 forAll(f.boundaryField()[patchId], iFace)
420 df1dn = (psi2[patchId][iFace][ie1_] - f[ic4e_[patchId][iFace]][ie1_]) / mv42e_[patchId][iFace];
421 df2dn = (psi2[patchId][iFace][ie2_] - f[ic4e_[patchId][iFace]][ie2_]) / mv42e_[patchId][iFace];
424 df1dt = (pF[ip3e_[patchId][iFace]][ie1_] - pF[ip1e_[patchId][iFace]][ie1_]) / mv13e_[patchId][iFace];
425 df2dt = (pF[ip3e_[patchId][iFace]][ie2_] - pF[ip1e_[patchId][iFace]][ie2_]) / mv13e_[patchId][iFace];
427 divf.boundaryFieldRef()[patchId][iFace] =
428 (df1dn*c1e_[patchId][iFace] - df1dt*c2e_[patchId][iFace])
430 (df2dt*c3e_[patchId][iFace] - df2dn*c4e_[patchId][iFace]);
436 #warning "Raise error if called not for 2D case"
442 scalar df11dn = 0.0, df11dt = 0.0,
443 df21dn = 0.0, df21dt = 0.0,
444 df22dn = 0.0, df22dt = 0.0,
445 df12dn = 0.0, df12dt = 0.0;
447 label i11 = ie1_*3 + ie1_;
448 label i21 = ie2_*3 + ie1_;
449 label i22 = ie2_*3 + ie2_;
450 label i12 = ie1_*3 + ie2_;
452 if (f.mesh().nGeometricD() == 2)
462 forAll(f.mesh().neighbour(), iFace)
464 df11dn = (f[ic2_[iFace]][i11] - f[ic4_[iFace]][i11]) / mv42_[iFace];
465 df21dn = (f[ic2_[iFace]][i21] - f[ic4_[iFace]][i21]) / mv42_[iFace];
468 df11dt = (pF[ip3_[iFace]][i11] - pF[ip1_[iFace]][i11]) / mv13_[iFace];
469 df21dt = (pF[ip3_[iFace]][i21] - pF[ip1_[iFace]][i21]) / mv13_[iFace];
472 (df11dn*c1_[iFace] - df11dt*c2_[iFace]) +
473 (df21dt*c3_[iFace] - df21dn*c4_[iFace]);
475 df22dn = (f[ic2_[iFace]][i22] - f[ic4_[iFace]][i22]) / mv42_[iFace];
476 df12dn = (f[ic2_[iFace]][i12] - f[ic4_[iFace]][i12]) / mv42_[iFace];
479 df22dt = (pF[ip3_[iFace]][i22] - pF[ip1_[iFace]][i22]) / mv13_[iFace];
480 df12dt = (pF[ip3_[iFace]][i12] - pF[ip1_[iFace]][i12]) / mv13_[iFace];
483 (df12dn*c1_[iFace] - df12dt*c2_[iFace]) +
484 (df22dt*c3_[iFace] - df22dn*c4_[iFace]);
487 List<List<tensor> > psi2(f.boundaryField().size());
489 forAll(ordinaryPatches_, iPatch)
491 label patchId = ordinaryPatches_[iPatch];
492 if (processorPatch_[iPatch])
494 psi2[patchId] = refCast<const processorFvPatchField<tensor> >
495 (f.boundaryField()[patchId]).patchNeighbourField();
499 psi2[patchId] = f.boundaryField()[patchId] +
500 f.boundaryField()[patchId].snGrad()
505 forAll(f.boundaryField()[patchId], iFace)
508 df11dn = (psi2[patchId][iFace][i11] - f[ic4e_[patchId][iFace]][i11]) / mv42e_[patchId][iFace];
509 df21dn = (psi2[patchId][iFace][i21] - f[ic4e_[patchId][iFace]][i21]) / mv42e_[patchId][iFace];
512 df11dt = (pF[ip3e_[patchId][iFace]][i11] - pF[ip1e_[patchId][iFace]][i11]) / mv13e_[patchId][iFace];
513 df21dt = (pF[ip3e_[patchId][iFace]][i21] - pF[ip1e_[patchId][iFace]][i21]) / mv13e_[patchId][iFace];
515 divf.boundaryFieldRef()[patchId][iFace][ie1_] =
516 (df11dn*c1e_[patchId][iFace] - df11dt*c2e_[patchId][iFace])
518 (df21dt*c3e_[patchId][iFace] - df21dn*c4e_[patchId][iFace]);
521 df12dn = (psi2[patchId][iFace][i12] - f[ic4e_[patchId][iFace]][i12]) / mv42e_[patchId][iFace];
522 df22dn = (psi2[patchId][iFace][i22] - f[ic4e_[patchId][iFace]][i22]) / mv42e_[patchId][iFace];
525 df12dt = (pF[ip3e_[patchId][iFace]][i12] - pF[ip1e_[patchId][iFace]][i12]) / mv13e_[patchId][iFace];
526 df22dt = (pF[ip3e_[patchId][iFace]][i22] - pF[ip1e_[patchId][iFace]][i22]) / mv13e_[patchId][iFace];
528 divf.boundaryFieldRef()[patchId][iFace][ie2_] =
529 (df12dn*c1e_[patchId][iFace] - df12dt*c2e_[patchId][iFace])
531 (df22dt*c3e_[patchId][iFace] - df22dn*c4e_[patchId][iFace]);
537 #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_
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_