FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 \*---------------------------------------------------------------------------*/
25 
26 #include <finiteVolume/fvMesh.H>
27 #include <finiteVolume/volFields.H>
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  defineTypeNameAndDebug(surfaceInterpolation, 0);
38 }
39 
40 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
41 
43 {
44  deleteDemandDrivenData(weightingFactors_);
45  deleteDemandDrivenData(differenceFactors_);
46  deleteDemandDrivenData(correctionVectors_);
47 }
48 
49 
50 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
51 
53 :
54  fvSchemes(static_cast<const objectRegistry&>(fvm)),
55  fvSolution(static_cast<const objectRegistry&>(fvm)),
56  mesh_(fvm),
57  weightingFactors_(NULL),
58  differenceFactors_(NULL),
59  orthogonal_(false),
60  correctionVectors_(NULL)
61 {}
62 
63 
64 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
65 
67 {
68  clearOut();
69 }
70 
71 
72 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
73 
75 {
76  if (!weightingFactors_)
77  {
78  makeWeights();
79  }
80 
81  return (*weightingFactors_);
82 }
83 
84 
86 {
87  if (!differenceFactors_)
88  {
89  makeDeltaCoeffs();
90  }
91 
92  return (*differenceFactors_);
93 }
94 
95 
97 {
98  if (orthogonal_ == false && !correctionVectors_)
99  {
100  makeCorrectionVectors();
101  }
102 
103  return orthogonal_;
104 }
105 
106 
109 {
110  if (orthogonal())
111  {
112  FatalErrorIn("surfaceInterpolation::correctionVectors()")
113  << "cannot return correctionVectors; mesh is orthogonal"
114  << abort(FatalError);
115  }
116 
117  return (*correctionVectors_);
118 }
119 
120 
122 {
123  deleteDemandDrivenData(weightingFactors_);
124  deleteDemandDrivenData(differenceFactors_);
125 
126  orthogonal_ = false;
127  deleteDemandDrivenData(correctionVectors_);
128 
129  return true;
130 }
131 
132 
133 void Foam::surfaceInterpolation::makeWeights() const
134 {
135  if (debug)
136  {
137  Info<< "surfaceInterpolation::makeWeights() : "
138  << "Constructing weighting factors for face interpolation"
139  << endl;
140  }
141 
142 
143  weightingFactors_ = new surfaceScalarField
144  (
145  IOobject
146  (
147  "weightingFactors",
148  mesh_.pointsInstance(),
149  mesh_
150  ),
151  mesh_,
152  dimless
153  );
154  surfaceScalarField& weightingFactors = *weightingFactors_;
155 
156 
157  // Set local references to mesh data
158  // (note that we should not use fvMesh sliced fields at this point yet
159  // since this causes a loop when generating weighting factors in
160  // coupledFvPatchField evaluation phase)
161  const unallocLabelList& owner = mesh_.owner();
162  const unallocLabelList& neighbour = mesh_.neighbour();
163 
164  const vectorField& Cf = mesh_.faceCentres();
165  const vectorField& C = mesh_.cellCentres();
166  const vectorField& Sf = mesh_.faceAreas();
167 
168  // ... and reference to the internal field of the weighting factors
169  scalarField& w = weightingFactors.internalField();
170 
171  forAll(owner, facei)
172  {
173  // Note: mag in the dot-product.
174  // For all valid meshes, the non-orthogonality will be less that
175  // 90 deg and the dot-product will be positive. For invalid
176  // meshes (d & s <= 0), this will stabilise the calculation
177  // but the result will be poor.
178  scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]]));
179  scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei]));
180  w[facei] = SfdNei/(SfdOwn + SfdNei);
181  }
182 
183  forAll(mesh_.boundary(), patchi)
184  {
185  mesh_.boundary()[patchi].makeWeights
186  (
187  weightingFactors.boundaryField()[patchi]
188  );
189  }
190 
191 
192  if (debug)
193  {
194  Info<< "surfaceInterpolation::makeWeights() : "
195  << "Finished constructing weighting factors for face interpolation"
196  << endl;
197  }
198 }
199 
200 
201 void Foam::surfaceInterpolation::makeDeltaCoeffs() const
202 {
203  if (debug)
204  {
205  Info<< "surfaceInterpolation::makeDeltaCoeffs() : "
206  << "Constructing differencing factors array for face gradient"
207  << endl;
208  }
209 
210  // Force the construction of the weighting factors
211  // needed to make sure deltaCoeffs are calculated for parallel runs.
212  weights();
213 
214  differenceFactors_ = new surfaceScalarField
215  (
216  IOobject
217  (
218  "differenceFactors_",
219  mesh_.pointsInstance(),
220  mesh_
221  ),
222  mesh_,
224  );
225  surfaceScalarField& DeltaCoeffs = *differenceFactors_;
226 
227 
228  // Set local references to mesh data
229  const volVectorField& C = mesh_.C();
230  const unallocLabelList& owner = mesh_.owner();
231  const unallocLabelList& neighbour = mesh_.neighbour();
232  const surfaceVectorField& Sf = mesh_.Sf();
233  const surfaceScalarField& magSf = mesh_.magSf();
234 
235  forAll(owner, facei)
236  {
237  vector delta = C[neighbour[facei]] - C[owner[facei]];
238  vector unitArea = Sf[facei]/magSf[facei];
239 
240  // Standard cell-centre distance form
241  //DeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
242 
243  // Slightly under-relaxed form
244  //DeltaCoeffs[facei] = 1.0/mag(delta);
245 
246  // More under-relaxed form
247  //DeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
248 
249  // Stabilised form for bad meshes
250  DeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
251  }
252 
253  forAll(DeltaCoeffs.boundaryField(), patchi)
254  {
255  mesh_.boundary()[patchi].makeDeltaCoeffs
256  (
257  DeltaCoeffs.boundaryField()[patchi]
258  );
259  }
260 }
261 
262 
263 void Foam::surfaceInterpolation::makeCorrectionVectors() const
264 {
265  if (debug)
266  {
267  Info<< "surfaceInterpolation::makeCorrectionVectors() : "
268  << "Constructing non-orthogonal correction vectors"
269  << endl;
270  }
271 
272  correctionVectors_ = new surfaceVectorField
273  (
274  IOobject
275  (
276  "correctionVectors",
277  mesh_.pointsInstance(),
278  mesh_
279  ),
280  mesh_,
281  dimless
282  );
283  surfaceVectorField& corrVecs = *correctionVectors_;
284 
285  // Set local references to mesh data
286  const volVectorField& C = mesh_.C();
287  const unallocLabelList& owner = mesh_.owner();
288  const unallocLabelList& neighbour = mesh_.neighbour();
289  const surfaceVectorField& Sf = mesh_.Sf();
290  const surfaceScalarField& magSf = mesh_.magSf();
291  const surfaceScalarField& DeltaCoeffs = deltaCoeffs();
292 
293  forAll(owner, facei)
294  {
295  vector unitArea = Sf[facei]/magSf[facei];
296  vector delta = C[neighbour[facei]] - C[owner[facei]];
297 
298  corrVecs[facei] = unitArea - delta*DeltaCoeffs[facei];
299  }
300 
301  // Boundary correction vectors set to zero for boundary patches
302  // and calculated consistently with internal corrections for
303  // coupled patches
304 
305  forAll(corrVecs.boundaryField(), patchi)
306  {
307  fvsPatchVectorField& patchcorrVecs = corrVecs.boundaryField()[patchi];
308 
309  if (!patchcorrVecs.coupled())
310  {
311  patchcorrVecs = vector::zero;
312  }
313  else
314  {
315  const fvsPatchScalarField& patchDeltaCoeffs
316  = DeltaCoeffs.boundaryField()[patchi];
317 
318  const fvPatch& p = patchcorrVecs.patch();
319 
320  vectorField patchDeltas = mesh_.boundary()[patchi].delta();
321 
322  forAll(p, patchFacei)
323  {
324  vector unitArea =
325  Sf.boundaryField()[patchi][patchFacei]
326  /magSf.boundaryField()[patchi][patchFacei];
327 
328  const vector& delta = patchDeltas[patchFacei];
329 
330  patchcorrVecs[patchFacei] =
331  unitArea - delta*patchDeltaCoeffs[patchFacei];
332  }
333  }
334  }
335 
336  scalar MaxNonOrthog = 0.0;
337 
338  // Calculate the non-orthogonality for meshes with 1 face or more
339  if (returnReduce(magSf.size(), sumOp<label>()) > 0)
340  {
341  MaxNonOrthog =
342  asin
343  (
344  min
345  (
346  max(mag(corrVecs)).value(),
347  1.0
348  )
349  )*180.0/mathematicalConstant::pi;
350  }
351 
352  if (debug)
353  {
354  Info<< "surfaceInterpolation::makeCorrectionVectors() : "
355  << "maximum non-orthogonality = " << MaxNonOrthog << " deg."
356  << endl;
357  }
358 
359  //MaxNonOrthog = 0.0;
360 
361  if (MaxNonOrthog < 5)
362  {
363  orthogonal_ = true;
364  deleteDemandDrivenData(correctionVectors_);
365  }
366  else
367  {
368  orthogonal_ = false;
369  }
370 
371  if (debug)
372  {
373  Info<< "surfaceInterpolation::makeCorrectionVectors() : "
374  << "Finished constructing non-orthogonal correction vectors"
375  << endl;
376  }
377 }
378 
379 
380 // ************************ vim: set sw=4 sts=4 et: ************************ //