surfaceInterpolation.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 | Copyright (C) 2011-2016 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 {
39  defineTypeNameAndDebug(surfaceInterpolation, 0);
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  scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]]));
182  scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei]));
183  w[facei] = SfdNei/(SfdOwn + SfdNei);
184  }
185 
186  surfaceScalarField::Boundary& wBf =
187  weights.boundaryFieldRef();
188 
189  forAll(mesh_.boundary(), patchi)
190  {
191  mesh_.boundary()[patchi].makeWeights(wBf[patchi]);
192  }
193 
194  if (debug)
195  {
196  Pout<< "surfaceInterpolation::makeWeights() : "
197  << "Finished constructing weighting factors for face interpolation"
198  << endl;
199  }
200 }
201 
202 
203 void Foam::surfaceInterpolation::makeDeltaCoeffs() const
204 {
205  if (debug)
206  {
207  Pout<< "surfaceInterpolation::makeDeltaCoeffs() : "
208  << "Constructing differencing factors array for face gradient"
209  << endl;
210  }
211 
212  // Force the construction of the weighting factors
213  // needed to make sure deltaCoeffs are calculated for parallel runs.
214  weights();
215 
216  deltaCoeffs_ = new surfaceScalarField
217  (
218  IOobject
219  (
220  "deltaCoeffs",
221  mesh_.pointsInstance(),
222  mesh_,
225  false // Do not register
226  ),
227  mesh_,
229  );
230  surfaceScalarField& deltaCoeffs = *deltaCoeffs_;
231 
232 
233  // Set local references to mesh data
234  const volVectorField& C = mesh_.C();
235  const labelUList& owner = mesh_.owner();
236  const labelUList& neighbour = mesh_.neighbour();
237 
238  forAll(owner, facei)
239  {
240  deltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
241  }
242 
243  surfaceScalarField::Boundary& deltaCoeffsBf =
244  deltaCoeffs.boundaryFieldRef();
245 
246  forAll(deltaCoeffsBf, patchi)
247  {
248  deltaCoeffsBf[patchi] = 1.0/mag(mesh_.boundary()[patchi].delta());
249  }
250 }
251 
252 
253 void Foam::surfaceInterpolation::makeNonOrthDeltaCoeffs() const
254 {
255  if (debug)
256  {
257  Pout<< "surfaceInterpolation::makeNonOrthDeltaCoeffs() : "
258  << "Constructing differencing factors array for face gradient"
259  << endl;
260  }
261 
262  // Force the construction of the weighting factors
263  // needed to make sure deltaCoeffs are calculated for parallel runs.
264  weights();
265 
266  nonOrthDeltaCoeffs_ = new surfaceScalarField
267  (
268  IOobject
269  (
270  "nonOrthDeltaCoeffs",
271  mesh_.pointsInstance(),
272  mesh_,
275  false // Do not register
276  ),
277  mesh_,
279  );
280  surfaceScalarField& nonOrthDeltaCoeffs = *nonOrthDeltaCoeffs_;
281 
282 
283  // Set local references to mesh data
284  const volVectorField& C = mesh_.C();
285  const labelUList& owner = mesh_.owner();
286  const labelUList& neighbour = mesh_.neighbour();
287  const surfaceVectorField& Sf = mesh_.Sf();
288  const surfaceScalarField& magSf = mesh_.magSf();
289 
290  forAll(owner, facei)
291  {
292  vector delta = C[neighbour[facei]] - C[owner[facei]];
293  vector unitArea = Sf[facei]/magSf[facei];
294 
295  // Standard cell-centre distance form
296  //NonOrthDeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
297 
298  // Slightly under-relaxed form
299  //NonOrthDeltaCoeffs[facei] = 1.0/mag(delta);
300 
301  // More under-relaxed form
302  //NonOrthDeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
303 
304  // Stabilised form for bad meshes
305  nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
306  }
307 
308  surfaceScalarField::Boundary& nonOrthDeltaCoeffsBf =
309  nonOrthDeltaCoeffs.boundaryFieldRef();
310 
311  forAll(nonOrthDeltaCoeffsBf, patchi)
312  {
313  vectorField delta(mesh_.boundary()[patchi].delta());
314 
315  nonOrthDeltaCoeffsBf[patchi] =
316  1.0/max(mesh_.boundary()[patchi].nf() & delta, 0.05*mag(delta));
317  }
318 }
319 
320 
321 void Foam::surfaceInterpolation::makeNonOrthCorrectionVectors() const
322 {
323  if (debug)
324  {
325  Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
326  << "Constructing non-orthogonal correction vectors"
327  << endl;
328  }
329 
330  nonOrthCorrectionVectors_ = new surfaceVectorField
331  (
332  IOobject
333  (
334  "nonOrthCorrectionVectors",
335  mesh_.pointsInstance(),
336  mesh_,
339  false // Do not register
340  ),
341  mesh_,
342  dimless
343  );
344  surfaceVectorField& corrVecs = *nonOrthCorrectionVectors_;
345 
346  // Set local references to mesh data
347  const volVectorField& C = mesh_.C();
348  const labelUList& owner = mesh_.owner();
349  const labelUList& neighbour = mesh_.neighbour();
350  const surfaceVectorField& Sf = mesh_.Sf();
351  const surfaceScalarField& magSf = mesh_.magSf();
352  const surfaceScalarField& NonOrthDeltaCoeffs = nonOrthDeltaCoeffs();
353 
354  forAll(owner, facei)
355  {
356  vector unitArea = Sf[facei]/magSf[facei];
357  vector delta = C[neighbour[facei]] - C[owner[facei]];
358 
359  corrVecs[facei] = unitArea - delta*NonOrthDeltaCoeffs[facei];
360  }
361 
362  // Boundary correction vectors set to zero for boundary patches
363  // and calculated consistently with internal corrections for
364  // coupled patches
365 
366  surfaceVectorField::Boundary& corrVecsBf =
367  corrVecs.boundaryFieldRef();
368 
369  forAll(corrVecsBf, patchi)
370  {
371  fvsPatchVectorField& patchCorrVecs = corrVecsBf[patchi];
372 
373  if (!patchCorrVecs.coupled())
374  {
375  patchCorrVecs = Zero;
376  }
377  else
378  {
379  const fvsPatchScalarField& patchNonOrthDeltaCoeffs
380  = NonOrthDeltaCoeffs.boundaryField()[patchi];
381 
382  const fvPatch& p = patchCorrVecs.patch();
383 
384  const vectorField patchDeltas(mesh_.boundary()[patchi].delta());
385 
386  forAll(p, patchFacei)
387  {
388  vector unitArea =
389  Sf.boundaryField()[patchi][patchFacei]
390  /magSf.boundaryField()[patchi][patchFacei];
391 
392  const vector& delta = patchDeltas[patchFacei];
393 
394  patchCorrVecs[patchFacei] =
395  unitArea - delta*patchNonOrthDeltaCoeffs[patchFacei];
396  }
397  }
398  }
399 
400  if (debug)
401  {
402  Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
403  << "Finished constructing non-orthogonal correction vectors"
404  << endl;
405  }
406 }
407 
408 
409 // ************************************************************************* //
Foam::surfaceFields.
scalar delta
Graphite solid properties.
Definition: C.H:48
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
const surfaceVectorField & Sf() const
Return cell face area vectors.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const Boundary & boundaryField() const
Return const-reference to the boundary field.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:61
bool movePoints()
Do what is neccessary if the mesh has moved.
const surfaceVectorField & nonOrthCorrectionVectors() const
Return reference to non-orthogonality correction vectors.
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:288
const surfaceScalarField & nonOrthDeltaCoeffs() const
Return reference to non-orthogonal cell-centre difference.
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:795
const surfaceScalarField & weights() const
Return reference to linear difference weighting factors.
virtual bool coupled() const
Return true if this patch field is coupled.
surfaceInterpolation(const fvMesh &)
Construct given an fvMesh.
static const zero Zero
Definition: zero.H:91
const vectorField & cellCentres() const
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:61
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
defineTypeNameAndDebug(combustionModel, 0)
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:282
const fvPatch & patch() const
Return patch.
const vectorField & faceCentres() const
Template functions to aid in the implementation of demand driven data.
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
void clearOut()
Clear all geometry and addressing.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
const vectorField & faceAreas() const
dimensioned< scalar > mag(const dimensioned< Type > &)
const volVectorField & C() const
Return cell centres as volVectorField.
volScalarField & p
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
void deleteDemandDrivenData(DataPtr &dataPtr)
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:65
const surfaceScalarField & deltaCoeffs() const
Return reference to cell-centre difference coefficients.
Namespace for OpenFOAM.
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:545