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-2013 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_(NULL),
60  deltaCoeffs_(NULL),
61  nonOrthDeltaCoeffs_(NULL),
62  nonOrthCorrectionVectors_(NULL)
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 
124 // Do what is neccessary if the mesh has moved
126 {
127  deleteDemandDrivenData(weights_);
128  deleteDemandDrivenData(deltaCoeffs_);
129  deleteDemandDrivenData(nonOrthDeltaCoeffs_);
130  deleteDemandDrivenData(nonOrthCorrectionVectors_);
131 
132  return true;
133 }
134 
135 
136 void Foam::surfaceInterpolation::makeWeights() const
137 {
138  if (debug)
139  {
140  Pout<< "surfaceInterpolation::makeWeights() : "
141  << "Constructing weighting factors for face interpolation"
142  << endl;
143  }
144 
145  weights_ = new surfaceScalarField
146  (
147  IOobject
148  (
149  "weights",
150  mesh_.pointsInstance(),
151  mesh_,
154  false // Do not register
155  ),
156  mesh_,
157  dimless
158  );
159  surfaceScalarField& weights = *weights_;
160 
161  // Set local references to mesh data
162  // (note that we should not use fvMesh sliced fields at this point yet
163  // since this causes a loop when generating weighting factors in
164  // coupledFvPatchField evaluation phase)
165  const labelUList& owner = mesh_.owner();
166  const labelUList& neighbour = mesh_.neighbour();
167 
168  const vectorField& Cf = mesh_.faceCentres();
169  const vectorField& C = mesh_.cellCentres();
170  const vectorField& Sf = mesh_.faceAreas();
171 
172  // ... and reference to the internal field of the weighting factors
173  scalarField& w = weights.internalField();
174 
175  forAll(owner, facei)
176  {
177  // Note: mag in the dot-product.
178  // For all valid meshes, the non-orthogonality will be less that
179  // 90 deg and the dot-product will be positive. For invalid
180  // meshes (d & s <= 0), this will stabilise the calculation
181  // but the result will be poor.
182  scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]]));
183  scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei]));
184  w[facei] = SfdNei/(SfdOwn + SfdNei);
185  }
186 
187  forAll(mesh_.boundary(), patchi)
188  {
189  mesh_.boundary()[patchi].makeWeights
190  (
191  weights.boundaryField()[patchi]
192  );
193  }
194 
195  if (debug)
196  {
197  Pout<< "surfaceInterpolation::makeWeights() : "
198  << "Finished constructing weighting factors for face interpolation"
199  << endl;
200  }
201 }
202 
203 
204 void Foam::surfaceInterpolation::makeDeltaCoeffs() const
205 {
206  if (debug)
207  {
208  Pout<< "surfaceInterpolation::makeDeltaCoeffs() : "
209  << "Constructing differencing factors array for face gradient"
210  << endl;
211  }
212 
213  // Force the construction of the weighting factors
214  // needed to make sure deltaCoeffs are calculated for parallel runs.
215  weights();
216 
217  deltaCoeffs_ = new surfaceScalarField
218  (
219  IOobject
220  (
221  "deltaCoeffs",
222  mesh_.pointsInstance(),
223  mesh_,
226  false // Do not register
227  ),
228  mesh_,
230  );
231  surfaceScalarField& DeltaCoeffs = *deltaCoeffs_;
232 
233 
234  // Set local references to mesh data
235  const volVectorField& C = mesh_.C();
236  const labelUList& owner = mesh_.owner();
237  const labelUList& neighbour = mesh_.neighbour();
238 
239  forAll(owner, facei)
240  {
241  DeltaCoeffs[facei] = 1.0/mag(C[neighbour[facei]] - C[owner[facei]]);
242  }
243 
244  forAll(DeltaCoeffs.boundaryField(), patchi)
245  {
246  DeltaCoeffs.boundaryField()[patchi] =
247  1.0/mag(mesh_.boundary()[patchi].delta());
248  }
249 }
250 
251 
252 void Foam::surfaceInterpolation::makeNonOrthDeltaCoeffs() const
253 {
254  if (debug)
255  {
256  Pout<< "surfaceInterpolation::makeNonOrthDeltaCoeffs() : "
257  << "Constructing differencing factors array for face gradient"
258  << endl;
259  }
260 
261  // Force the construction of the weighting factors
262  // needed to make sure deltaCoeffs are calculated for parallel runs.
263  weights();
264 
265  nonOrthDeltaCoeffs_ = new surfaceScalarField
266  (
267  IOobject
268  (
269  "nonOrthDeltaCoeffs",
270  mesh_.pointsInstance(),
271  mesh_,
274  false // Do not register
275  ),
276  mesh_,
278  );
279  surfaceScalarField& nonOrthDeltaCoeffs = *nonOrthDeltaCoeffs_;
280 
281 
282  // Set local references to mesh data
283  const volVectorField& C = mesh_.C();
284  const labelUList& owner = mesh_.owner();
285  const labelUList& neighbour = mesh_.neighbour();
286  const surfaceVectorField& Sf = mesh_.Sf();
287  const surfaceScalarField& magSf = mesh_.magSf();
288 
289  forAll(owner, facei)
290  {
291  vector delta = C[neighbour[facei]] - C[owner[facei]];
292  vector unitArea = Sf[facei]/magSf[facei];
293 
294  // Standard cell-centre distance form
295  //NonOrthDeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
296 
297  // Slightly under-relaxed form
298  //NonOrthDeltaCoeffs[facei] = 1.0/mag(delta);
299 
300  // More under-relaxed form
301  //NonOrthDeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
302 
303  // Stabilised form for bad meshes
304  nonOrthDeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
305  }
306 
307  forAll(nonOrthDeltaCoeffs.boundaryField(), patchi)
308  {
309  vectorField delta(mesh_.boundary()[patchi].delta());
310 
311  nonOrthDeltaCoeffs.boundaryField()[patchi] =
312  1.0/max(mesh_.boundary()[patchi].nf() & delta, 0.05*mag(delta));
313  }
314 }
315 
316 
317 void Foam::surfaceInterpolation::makeNonOrthCorrectionVectors() const
318 {
319  if (debug)
320  {
321  Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
322  << "Constructing non-orthogonal correction vectors"
323  << endl;
324  }
325 
326  nonOrthCorrectionVectors_ = new surfaceVectorField
327  (
328  IOobject
329  (
330  "nonOrthCorrectionVectors",
331  mesh_.pointsInstance(),
332  mesh_,
335  false // Do not register
336  ),
337  mesh_,
338  dimless
339  );
340  surfaceVectorField& corrVecs = *nonOrthCorrectionVectors_;
341 
342  // Set local references to mesh data
343  const volVectorField& C = mesh_.C();
344  const labelUList& owner = mesh_.owner();
345  const labelUList& neighbour = mesh_.neighbour();
346  const surfaceVectorField& Sf = mesh_.Sf();
347  const surfaceScalarField& magSf = mesh_.magSf();
348  const surfaceScalarField& NonOrthDeltaCoeffs = nonOrthDeltaCoeffs();
349 
350  forAll(owner, facei)
351  {
352  vector unitArea = Sf[facei]/magSf[facei];
353  vector delta = C[neighbour[facei]] - C[owner[facei]];
354 
355  corrVecs[facei] = unitArea - delta*NonOrthDeltaCoeffs[facei];
356  }
357 
358  // Boundary correction vectors set to zero for boundary patches
359  // and calculated consistently with internal corrections for
360  // coupled patches
361 
362  forAll(corrVecs.boundaryField(), patchi)
363  {
364  fvsPatchVectorField& patchCorrVecs = corrVecs.boundaryField()[patchi];
365 
366  if (!patchCorrVecs.coupled())
367  {
368  patchCorrVecs = vector::zero;
369  }
370  else
371  {
372  const fvsPatchScalarField& patchNonOrthDeltaCoeffs
373  = NonOrthDeltaCoeffs.boundaryField()[patchi];
374 
375  const fvPatch& p = patchCorrVecs.patch();
376 
377  const vectorField patchDeltas(mesh_.boundary()[patchi].delta());
378 
379  forAll(p, patchFacei)
380  {
381  vector unitArea =
382  Sf.boundaryField()[patchi][patchFacei]
383  /magSf.boundaryField()[patchi][patchFacei];
384 
385  const vector& delta = patchDeltas[patchFacei];
386 
387  patchCorrVecs[patchFacei] =
388  unitArea - delta*patchNonOrthDeltaCoeffs[patchFacei];
389  }
390  }
391  }
392 
393  if (debug)
394  {
395  Pout<< "surfaceInterpolation::makeNonOrthCorrectionVectors() : "
396  << "Finished constructing non-orthogonal correction vectors"
397  << endl;
398  }
399 }
400 
401 
402 // ************************************************************************* //
const vectorField & faceAreas() const
virtual bool coupled() const
Return true if this patch field is coupled.
const surfaceScalarField & weights() const
Return reference to linear difference weighting factors.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
Foam::surfaceFields.
dimensioned< scalar > mag(const dimensioned< Type > &)
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
const surfaceVectorField & Sf() const
Return cell face area vectors.
const fvPatch & patch() const
Return patch.
void deleteDemandDrivenData(DataPtr &dataPtr)
InternalField & internalField()
Return internal field.
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:282
Graphite solid properties.
Definition: C.H:57
const surfaceScalarField & deltaCoeffs() const
Return reference to cell-centre difference coefficients.
const vectorField & cellCentres() const
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:818
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
Namespace for OpenFOAM.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
const volVectorField & C() const
Return cell centres as volVectorField.
volScalarField & p
Definition: createFields.H:51
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
label patchi
scalar delta
Template functions to aid in the implementation of demand driven data.
const surfaceVectorField & nonOrthCorrectionVectors() const
Return reference to non-orthogonality correction vectors.
const surfaceScalarField & nonOrthDeltaCoeffs() const
Return reference to non-orthogonal cell-centre difference.
void clearOut()
Clear all geometry and addressing.
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:60
static const Vector zero
Definition: Vector.H:80
surfaceInterpolation(const fvMesh &)
Construct given an fvMesh.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:61
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:288
An abstract base class with a fat-interface to all derived classes covering all possible ways in whic...
Definition: fvsPatchField.H:65
bool movePoints()
Do what is neccessary if the mesh has moved.
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
const vectorField & faceCentres() const
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:552