interpolatedFaces.H
Go to the documentation of this file.
1// Interpolation used
2interpolationCellPoint<vector> UInterpolator(HbyA);
3
4// Determine faces on outside of interpolated cells
5bitSet isOwnerInterpolatedFace(mesh.nInternalFaces());
6bitSet isNeiInterpolatedFace(mesh.nInternalFaces());
7
8// Determine donor cells
9labelListList donorCell(mesh.nInternalFaces());
10
11scalarListList weightCellCells(mesh.nInternalFaces());
12
13// Interpolated HbyA faces
14vectorField UIntFaces(mesh.nInternalFaces(), Zero);
15
16// Determine receptor neighbour cells
17labelList receptorNeigCell(mesh.nInternalFaces(), -1);
18
19{
20 const cellCellStencilObject& overlap = Stencil::New(mesh);
21 const labelList& cellTypes = overlap.cellTypes();
22 const labelIOList& zoneID = overlap.zoneID();
23
24 label nZones = gMax(zoneID)+1;
25 PtrList<fvMeshSubset> meshParts(nZones);
26 labelList nCellsPerZone(nZones, Zero);
27
28 // A mesh subset for each zone
30 {
31 meshParts.set
32 (
33 zonei,
34 // Select cells where the zoneID == zonei
35 new fvMeshSubset(mesh, zonei, zoneID)
36 );
37 }
38
39 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
40 {
41 label ownType = cellTypes[mesh.faceOwner()[faceI]];
42 label neiType = cellTypes[mesh.faceNeighbour()[faceI]];
43 if
44 (
45 ownType == cellCellStencil::INTERPOLATED
46 && neiType == cellCellStencil::CALCULATED
47 )
48 {
49 isOwnerInterpolatedFace.set(faceI);
50
51 const vector& fc = mesh.faceCentres()[faceI];
52
53 for (label zoneI = 0; zoneI < nZones; zoneI++)
54 {
55 if (zoneI != zoneID[mesh.faceOwner()[faceI]])
56 {
57 const fvMesh& partMesh = meshParts[zoneI].subMesh();
58 const labelList& cellMap = meshParts[zoneI].cellMap();
59 label cellI = partMesh.findCell(fc);
60
61 if (cellI != -1)
62 {
63 // Determine weights
64 labelList stencil(partMesh.cellCells()[cellI]);
65
66 stencil.append(cellI);
67
68 label st = stencil.size();
69
70 donorCell[faceI].setSize(st);
71
72 weightCellCells[faceI].setSize(st);
73
74 scalarField weights(st);
75
76 forAll(stencil, i)
77 {
78 scalar d = mag
79 (
80 partMesh.cellCentres()[stencil[i]]
81 - fc
82 );
83 weights[i] = 1.0/d;
84 donorCell[faceI][i] = cellMap[stencil[i]];
85 }
86 weights /= sum(weights);
87
88 weightCellCells[faceI] = weights;
89
90 forAll(stencil, i)
91 {
92 UIntFaces[faceI] +=
93 weightCellCells[faceI][i]
94 *UInterpolator.interpolate
95 (
96 fc,
97 donorCell[faceI][i]
98 );
99 }
100
101 break;
102 }
103 }
104 }
105
106 receptorNeigCell[faceI] = mesh.faceNeighbour()[faceI];
107 }
108 else if
109 (
110 ownType == cellCellStencil::CALCULATED
111 && neiType == cellCellStencil::INTERPOLATED
112 )
113 {
114 isNeiInterpolatedFace.set(faceI);
115
116 const vector& fc = mesh.faceCentres()[faceI];
117 for (label zoneI = 0; zoneI < nZones; zoneI++)
118 {
119 if (zoneI != zoneID[mesh.faceNeighbour()[faceI]])
120 {
121 const fvMesh& partMesh = meshParts[zoneI].subMesh();
122 const labelList& cellMap = meshParts[zoneI].cellMap();
123 label cellI = partMesh.findCell(fc);
124
125 if (cellI != -1)
126 {
127 // Determine weights
128 labelList stencil(partMesh.cellCells()[cellI]);
129
130 stencil.append(cellI);
131
132 label st = stencil.size();
133
134 donorCell[faceI].setSize(st);
135
136 weightCellCells[faceI].setSize(st);
137
138 scalarField weights(st);
139
140 forAll(stencil, i)
141 {
142 scalar d = mag
143 (
144 partMesh.cellCentres()[stencil[i]]
145 - fc
146 );
147 weights[i] = 1.0/d;
148 donorCell[faceI][i] = cellMap[stencil[i]];
149 }
150 weights /= sum(weights);
151
152 weightCellCells[faceI] = weights;
153
154 forAll(stencil, i)
155 {
156 UIntFaces[faceI] +=
157 weightCellCells[faceI][i]
158 *UInterpolator.interpolate
159 (
160 fc,
161 donorCell[faceI][i]
162 );
163 }
164
165 break;
166 }
167 }
168 }
169
170 receptorNeigCell[faceI] = mesh.faceOwner()[faceI];
171 }
172 }
173}
174
175// contravariant U
176vectorField U1Contrav(mesh.nInternalFaces(), Zero);
177
178surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf());
179
181{
182 label cellId = -1;
183 if (isNeiInterpolatedFace.test(faceI))
184 {
185 cellId = mesh.faceNeighbour()[faceI];
186 }
187 else if (isOwnerInterpolatedFace.test(faceI))
188 {
189 cellId = mesh.faceOwner()[faceI];
190 }
191
192 if (cellId != -1)
193 {
194 const vector& n = faceNormals[faceI];
195 vector n1(Zero);
196
197 // 2-D cases
198 if (mesh.nSolutionD() == 2)
199 {
200 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
201 {
202 if (mesh.geometricD()[cmpt] == -1)
203 {
204 switch (cmpt)
205 {
206 case vector::X:
207 {
208 n1 = vector(0, n.z(), -n.y());
209 break;
210 }
211
212 case vector::Y:
213 {
214 n1 = vector(n.z(), 0, -n.x());
215 break;
216 }
217
218 case vector::Z:
219 {
220 n1 = vector(n.y(), -n.x(), 0);
221 break;
222 }
223 }
224 }
225 }
226 }
227 else if (mesh.nSolutionD() == 3)
228 {
229 //Determine which is the primary direction
230 if (mag(n.x()) > mag(n.y()) && mag(n.x()) > mag(n.z()))
231 {
232 n1 = vector(n.y(), -n.x(), 0);
233 }
234 else if (mag(n.y()) > mag(n.z()))
235 {
236 n1 = vector(0, n.z(), -n.y());
237 }
238 else
239 {
240 n1 = vector(-n.z(), 0, n.x());
241 }
242 }
243 n1.normalise();
244
245 const vector n2 = normalised(n ^ n1);
246
247 tensor rot =
248 tensor
249 (
250 n.x() ,n.y(), n.z(),
251 n1.x() ,n1.y(), n1.z(),
252 n2.x() ,n2.y(), n2.z()
253 );
254
255// tensor rot =
256// tensor
257// (
258// n & x ,n & y, n & z,
259// n1 & x ,n1 & y, n1 & z,
260// n2 & x ,n2 & y, n2 & z
261// );
262
263 U1Contrav[faceI].x() =
264 2*transform(rot, UIntFaces[faceI]).x()
265 - transform(rot, HbyA[receptorNeigCell[faceI]]).x();
266
267 U1Contrav[faceI].y() = transform(rot, HbyA[cellId]).y();
268
269 U1Contrav[faceI].z() = transform(rot, HbyA[cellId]).z();
270
271 HbyA[cellId] = transform(inv(rot), U1Contrav[faceI]);
272 }
273}
label n
static constexpr direction nComponents
Number of components in bool is 1.
Definition: bool.H:98
Tensor of scalars, i.e. Tensor<scalar>.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
HbyA
Definition: pcEqn.H:74
dynamicFvMesh & mesh
const cellCellStencilObject & overlap
Definition: correctPhi.H:57
label nZones
scalarListList weightCellCells(mesh.nInternalFaces())
interpolationCellPoint< vector > UInterpolator(HbyA)
labelListList donorCell(mesh.nInternalFaces())
surfaceVectorField faceNormals(mesh.Sf()/mesh.magSf())
bitSet isNeiInterpolatedFace(mesh.nInternalFaces())
labelList receptorNeigCell(mesh.nInternalFaces(), -1)
labelList nCellsPerZone(nZones, Zero)
vectorField UIntFaces(mesh.nInternalFaces(), Zero)
bitSet isOwnerInterpolatedFace(mesh.nInternalFaces())
PtrList< fvMeshSubset > meshParts(nZones)
const labelList & cellTypes
vectorField U1Contrav(mesh.nInternalFaces(), Zero)
const labelIOList & zoneID
label cellId
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333