surfaceInterpolation.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Description
25  Cell to face interpolation scheme. Included in fvMesh.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "fvMesh.H"
30 #include "volFields.H"
31 #include "surfaceFields.H"
32 #include "demandDrivenData.H"
33 #include "coupledFvPatch.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
40 }
41 
42 
43 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
44 
46 {
47  deleteDemandDrivenData(weights_);
48  deleteDemandDrivenData(deltaCoeffs_);
49  deleteDemandDrivenData(nonOrthDeltaCoeffs_);
50  deleteDemandDrivenData(nonOrthCorrectionVectors_);
51 }
52 
53 
54 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
55 
57 :
58  mesh_(fvm),
59  weights_(nullptr),
60  deltaCoeffs_(nullptr),
61  nonOrthDeltaCoeffs_(nullptr),
62  nonOrthCorrectionVectors_(nullptr)
63 {}
64 
65 
66 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
67 
69 {
70  clearOut();
71 }
72 
73 
74 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
75 
78 {
79  if (!weights_)
80  {
81  makeWeights();
82  }
83 
84  return (*weights_);
85 }
86 
87 
90 {
91  if (!deltaCoeffs_)
92  {
93  makeDeltaCoeffs();
94  }
95 
96  return (*deltaCoeffs_);
97 }
98 
99 
102 {
103  if (!nonOrthDeltaCoeffs_)
104  {
105  makeNonOrthDeltaCoeffs();
106  }
107 
108  return (*nonOrthDeltaCoeffs_);
109 }
110 
111 
114 {
115  if (!nonOrthCorrectionVectors_)
116  {
117  makeNonOrthCorrectionVectors();
118  }
119 
120  return (*nonOrthCorrectionVectors_);
121 }
122 
123 
125 {
126  deleteDemandDrivenData(weights_);
127  deleteDemandDrivenData(deltaCoeffs_);
128  deleteDemandDrivenData(nonOrthDeltaCoeffs_);
129  deleteDemandDrivenData(nonOrthCorrectionVectors_);
130 
131  return true;
132 }
133 
134 
135 void Foam::surfaceInterpolation::makeWeights() const
136 {
137  if (debug)
138  {
139  Pout<< "surfaceInterpolation::makeWeights() : "
140  << "Constructing weighting factors for face interpolation"
141  << endl;
142  }
143 
144  weights_ = new surfaceScalarField
145  (
146  IOobject
147  (
148  "weights",
149  mesh_.pointsInstance(),
150  mesh_,
153  false // Do not register
154  ),
155  mesh_,
156  dimless
157  );
158  surfaceScalarField& weights = *weights_;
159 
160  // Set local references to mesh data
161  // (note that we should not use fvMesh sliced fields at this point yet
162  // since this causes a loop when generating weighting factors in
163  // coupledFvPatchField evaluation phase)
164  const labelUList& owner = mesh_.owner();
165  const labelUList& neighbour = mesh_.neighbour();
166 
167  const vectorField& Cf = mesh_.faceCentres();
168  const vectorField& C = mesh_.cellCentres();
169  const vectorField& Sf = mesh_.faceAreas();
170 
171  // ... and reference to the internal field of the weighting factors
172  scalarField& w = weights.primitiveFieldRef();
173 
174  forAll(owner, facei)
175  {
176  // Note: mag in the dot-product.
177  // For all valid meshes, the non-orthogonality will be less that
178  // 90 deg and the dot-product will be positive. For invalid
179  // meshes (d & s <= 0), this will stabilise the calculation
180  // but the result will be poor.
181  const scalar SfdOwn = mag(Sf[facei]&(Cf[facei] - C[owner[facei]]));
182  const scalar SfdNei = mag(Sf[facei]&(C[neighbour[facei]] - Cf[facei]));
183  const scalar SfdOwnNei = SfdOwn + SfdNei;
184 
185  if (SfdNei/vGreat < SfdOwnNei)
186  {
187  w[facei] = SfdNei/SfdOwnNei;
188  }
189  else
190  {
191  const scalar dOwn = mag(Cf[facei] - C[owner[facei]]);
192  const scalar dNei = mag(C[neighbour[facei]] - Cf[facei]);
193  const scalar dOwnNei = dOwn + dNei;
194 
195  w[facei] = dNei/dOwnNei;
196  }
197  }
198 
200  weights.boundaryFieldRef();
201 
202  forAll(mesh_.boundary(), patchi)
203  {
204  mesh_.boundary()[patchi].makeWeights(wBf[patchi]);
205  }
206 
207  if (debug)
208  {
209  Pout<< "surfaceInterpolation::makeWeights() : "
210  << "Finished constructing weighting factors for face interpolation"
211  << endl;
212  }
213 }
214 
215 
216 void Foam::surfaceInterpolation::makeDeltaCoeffs() const
217 {
218  if (debug)
219  {
220  Pout<< "surfaceInterpolation::makeDeltaCoeffs() : "
221  << "Constructing interpolation factors array for face gradient"
222  << endl;
223  }
224 
225  // Force the construction of the weighting factors
226  // needed to make sure deltaCoeffs are calculated for parallel runs.
227  weights();
228 
229  deltaCoeffs_ = new surfaceScalarField
230  (
231  IOobject
232  (
233  "deltaCoeffs",
234  mesh_.pointsInstance(),
235  mesh_,
238  false // Do not register
239  ),
240  mesh_,
242  );
243  surfaceScalarField& deltaCoeffs = *deltaCoeffs_;
244 
245 
246  // Set local references to mesh data
247  const volVectorField& C = mesh_.C();
248  const labelUList& owner = mesh_.owner();
249  const labelUList& neighbour = mesh_.neighbour();
250 
251  forAll(owner, facei)
252  {
253  deltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
254  }
255 
256  surfaceScalarField::Boundary& deltaCoeffsBf =
257  deltaCoeffs.boundaryFieldRef();
258 
259  forAll(deltaCoeffsBf, patchi)
260  {
261  deltaCoeffsBf[patchi] = 1.0/mag(mesh_.boundary()[patchi].delta());
262  }
263 }
264 
265 
266 void Foam::surfaceInterpolation::makeNonOrthDeltaCoeffs() const
267 {
268  if (debug)
269  {
270  Pout<< "surfaceInterpolation::makeNonOrthDeltaCoeffs() : "
271  << "Constructing interpolation factors array for face gradient"
272  << endl;
273  }
274 
275  // Force the construction of the weighting factors
276  // needed to make sure deltaCoeffs are calculated for parallel runs.
277  weights();
278 
279  nonOrthDeltaCoeffs_ = new surfaceScalarField
280  (
281  IOobject
282  (
283  "nonOrthDeltaCoeffs",
284  mesh_.pointsInstance(),
285  mesh_,
288  false // Do not register
289  ),
290  mesh_,
292  );
293  surfaceScalarField& nonOrthDeltaCoeffs = *nonOrthDeltaCoeffs_;
294 
295 
296  // Set local references to mesh data
297  const volVectorField& C = mesh_.C();
298  const labelUList& owner = mesh_.owner();
299  const labelUList& neighbour = mesh_.neighbour();
300  const surfaceVectorField& Sf = mesh_.Sf();
301  const surfaceScalarField& magSf = mesh_.magSf();
302 
303  forAll(owner, facei)
304  {
305  vector delta = C[neighbour[facei]] - C[owner[facei]];
306  vector unitArea = Sf[facei]/magSf[facei];
307 
308  // Standard cell-centre distance form
309  // NonOrthDeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
310 
311  // Slightly under-relaxed form
312  // NonOrthDeltaCoeffs[facei] = 1.0/mag(delta);
313 
314  // More under-relaxed form
315  // NonOrthDeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + vSmall);
316 
317  // Stabilised form for bad meshes
318  nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
319  }
320 
321  surfaceScalarField::Boundary& nonOrthDeltaCoeffsBf =
322  nonOrthDeltaCoeffs.boundaryFieldRef();
323 
324  forAll(nonOrthDeltaCoeffsBf, patchi)
325  {
326  vectorField delta(mesh_.boundary()[patchi].delta());
327 
328  nonOrthDeltaCoeffsBf[patchi] =
329  1.0/max(mesh_.boundary()[patchi].nf() & delta, 0.05*mag(delta));
330  }
331 }
332 
333 
334 void Foam::surfaceInterpolation::makeNonOrthCorrectionVectors() const
335 {
336  if (debug)
337  {
338  Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
339  << "Constructing non-orthogonal correction vectors"
340  << endl;
341  }
342 
343  nonOrthCorrectionVectors_ = new surfaceVectorField
344  (
345  IOobject
346  (
347  "nonOrthCorrectionVectors",
348  mesh_.pointsInstance(),
349  mesh_,
352  false // Do not register
353  ),
354  mesh_,
355  dimless
356  );
357  surfaceVectorField& corrVecs = *nonOrthCorrectionVectors_;
358 
359  // Set local references to mesh data
360  const volVectorField& C = mesh_.C();
361  const labelUList& owner = mesh_.owner();
362  const labelUList& neighbour = mesh_.neighbour();
363  const surfaceVectorField& Sf = mesh_.Sf();
364  const surfaceScalarField& magSf = mesh_.magSf();
365  const surfaceScalarField& NonOrthDeltaCoeffs = nonOrthDeltaCoeffs();
366 
367  forAll(owner, facei)
368  {
369  vector unitArea = Sf[facei]/magSf[facei];
370  vector delta = C[neighbour[facei]] - C[owner[facei]];
371 
372  corrVecs[facei] = unitArea - delta*NonOrthDeltaCoeffs[facei];
373  }
374 
375  // Boundary correction vectors set to zero for boundary patches
376  // and calculated consistently with internal corrections for
377  // coupled patches
378 
379  surfaceVectorField::Boundary& corrVecsBf =
380  corrVecs.boundaryFieldRef();
381 
382  forAll(corrVecsBf, patchi)
383  {
384  fvsPatchVectorField& patchCorrVecs = corrVecsBf[patchi];
385 
386  if (!patchCorrVecs.coupled())
387  {
388  patchCorrVecs = Zero;
389  }
390  else
391  {
392  const fvsPatchScalarField& patchNonOrthDeltaCoeffs
393  = NonOrthDeltaCoeffs.boundaryField()[patchi];
394 
395  const fvPatch& p = patchCorrVecs.patch();
396 
397  const vectorField patchDeltas(mesh_.boundary()[patchi].delta());
398 
399  forAll(p, patchFacei)
400  {
401  vector unitArea =
402  Sf.boundaryField()[patchi][patchFacei]
403  /magSf.boundaryField()[patchi][patchFacei];
404 
405  const vector& delta = patchDeltas[patchFacei];
406 
407  patchCorrVecs[patchFacei] =
408  unitArea - delta*patchNonOrthDeltaCoeffs[patchFacei];
409  }
410  }
411  }
412 
413  if (debug)
414  {
415  Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
416  << "Finished constructing non-orthogonal correction vectors"
417  << endl;
418  }
419 }
420 
421 
422 // ************************************************************************* //
static const Foam::dimensionedScalar C("C", Foam::dimTemperature, 234.5)
scalar delta
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Generic GeometricField class.
GeometricBoundaryField< Type, PatchField, GeoMesh > Boundary
Type of the boundary field.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:101
Cell to surface interpolation scheme. Included in fvMesh.
bool movePoints()
Do what is necessary if the mesh has moved.
const surfaceVectorField & nonOrthCorrectionVectors() const
Return reference to non-orthogonality correction vectors.
const surfaceScalarField & weights() const
Return reference to linear difference weighting factors.
const surfaceScalarField & deltaCoeffs() const
Return reference to cell-centre difference coefficients.
surfaceInterpolation(const fvMesh &)
Construct given an fvMesh.
const surfaceScalarField & nonOrthDeltaCoeffs() const
Return reference to non-orthogonal cell-centre difference.
void clearOut()
Clear all geometry and addressing.
Template functions to aid in the implementation of demand driven data.
label patchi
Namespace for OpenFOAM.
static const zero Zero
Definition: zero.H:97
VolField< vector > volVectorField
Definition: volFieldsFwd.H:62
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const dimensionSet dimless
SurfaceField< scalar > surfaceScalarField
const dimensionSet dimLength
void deleteDemandDrivenData(DataPtr &dataPtr)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Field< vector > vectorField
Specialisation of Field<T> for vector.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
SurfaceField< vector > surfaceVectorField
UList< label > labelUList
Definition: UList.H:65
fvsPatchField< vector > fvsPatchVectorField
fvsPatchField< scalar > fvsPatchScalarField
volScalarField & p
Foam::surfaceFields.