fvMeshGeometry.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 \*---------------------------------------------------------------------------*/
25 
26 #include "fvMesh.H"
27 #include "Time.H"
28 #include "volFields.H"
29 #include "surfaceFields.H"
30 #include "slicedVolFields.H"
31 #include "slicedSurfaceFields.H"
32 #include "SubField.H"
33 #include "cyclicFvPatchFields.H"
34 #include "cyclicAMIFvPatchFields.H"
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 void fvMesh::makeSf() const
44 {
45  if (debug)
46  {
47  Info<< "void fvMesh::makeSf() : "
48  << "assembling face areas"
49  << endl;
50  }
51 
52  // It is an error to attempt to recalculate
53  // if the pointer is already set
54  if (SfPtr_)
55  {
56  FatalErrorIn("fvMesh::makeSf()")
57  << "face areas already exist"
58  << abort(FatalError);
59  }
60 
61  SfPtr_ = new slicedSurfaceVectorField
62  (
63  IOobject
64  (
65  "S",
67  meshSubDir,
68  *this,
71  false
72  ),
73  *this,
74  dimArea,
75  faceAreas()
76  );
77 }
78 
79 
80 void fvMesh::makeMagSf() const
81 {
82  if (debug)
83  {
84  Info<< "void fvMesh::makeMagSf() : "
85  << "assembling mag face areas"
86  << endl;
87  }
88 
89  // It is an error to attempt to recalculate
90  // if the pointer is already set
91  if (magSfPtr_)
92  {
93  FatalErrorIn("void fvMesh::makeMagSf()")
94  << "mag face areas already exist"
95  << abort(FatalError);
96  }
97 
98  // Note: Added stabilisation for faces with exactly zero area.
99  // These should be caught on mesh checking but at least this stops
100  // the code from producing Nans.
101  magSfPtr_ = new surfaceScalarField
102  (
103  IOobject
104  (
105  "magSf",
106  pointsInstance(),
107  meshSubDir,
108  *this,
111  false
112  ),
113  mag(Sf()) + dimensionedScalar("vs", dimArea, VSMALL)
114  );
115 }
116 
117 
118 void fvMesh::makeC() const
119 {
120  if (debug)
121  {
122  Info<< "void fvMesh::makeC() : "
123  << "assembling cell centres"
124  << endl;
125  }
126 
127  // It is an error to attempt to recalculate
128  // if the pointer is already set
129  if (CPtr_)
130  {
131  FatalErrorIn("fvMesh::makeC()")
132  << "cell centres already exist"
133  << abort(FatalError);
134  }
135 
136  // Construct as slices. Only preserve processor (not e.g. cyclic)
137 
138  CPtr_ = new slicedVolVectorField
139  (
140  IOobject
141  (
142  "C",
143  pointsInstance(),
144  meshSubDir,
145  *this,
148  false
149  ),
150  *this,
151  dimLength,
152  cellCentres(),
153  faceCentres(),
154  true, //preserveCouples
155  true //preserveProcOnly
156  );
157 }
158 
159 
160 void fvMesh::makeCf() const
161 {
162  if (debug)
163  {
164  Info<< "void fvMesh::makeCf() : "
165  << "assembling face centres"
166  << endl;
167  }
168 
169  // It is an error to attempt to recalculate
170  // if the pointer is already set
171  if (CfPtr_)
172  {
173  FatalErrorIn("fvMesh::makeCf()")
174  << "face centres already exist"
175  << abort(FatalError);
176  }
177 
178  CfPtr_ = new slicedSurfaceVectorField
179  (
180  IOobject
181  (
182  "Cf",
183  pointsInstance(),
184  meshSubDir,
185  *this,
188  false
189  ),
190  *this,
191  dimLength,
192  faceCentres()
193  );
194 }
195 
196 
197 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
198 
200 {
201  if (!VPtr_)
202  {
203  if (debug)
204  {
205  Info<< "fvMesh::V() const: "
206  << "constructing from primitiveMesh::cellVolumes()"
207  << endl;
208  }
209 
211  (
212  IOobject
213  (
214  "V",
215  time().timeName(),
216  *this,
219  false
220  ),
221  *this,
222  dimVolume,
223  cellVolumes()
224  );
225  }
226 
227  return *static_cast<slicedVolScalarField::DimensionedInternalField*>(VPtr_);
228 }
229 
230 
232 {
233  if (!V0Ptr_)
234  {
235  FatalErrorIn("fvMesh::V0() const")
236  << "V0 is not available"
237  << abort(FatalError);
238  }
239 
240  return *V0Ptr_;
241 }
242 
243 
245 {
246  if (!V0Ptr_)
247  {
248  FatalErrorIn("fvMesh::setV0()")
249  << "V0 is not available"
250  << abort(FatalError);
251  }
252 
253  return *V0Ptr_;
254 }
255 
256 
258 {
259  if (!V00Ptr_)
260  {
261  if (debug)
262  {
263  Info<< "fvMesh::V00() const: "
264  << "constructing from V0"
265  << endl;
266  }
267 
269  (
270  IOobject
271  (
272  "V00",
273  time().timeName(),
274  *this,
277  false
278  ),
279  V0()
280  );
281 
282  // If V00 is used then V0 should be stored for restart
283  V0Ptr_->writeOpt() = IOobject::AUTO_WRITE;
284  }
285 
286  return *V00Ptr_;
287 }
288 
289 
291 {
292  if (moving() && time().subCycling())
293  {
294  const TimeState& ts = time();
295  const TimeState& ts0 = time().prevTimeState();
296 
297  scalar tFrac =
298  (
299  ts.value() - (ts0.value() - ts0.deltaTValue())
300  )/ts0.deltaTValue();
301 
302  if (tFrac < (1 - SMALL))
303  {
304  return V0() + tFrac*(V() - V0());
305  }
306  else
307  {
308  return V();
309  }
310  }
311  else
312  {
313  return V();
314  }
315 }
316 
317 
319 {
320  if (moving() && time().subCycling())
321  {
322  const TimeState& ts = time();
323  const TimeState& ts0 = time().prevTimeState();
324 
325  scalar t0Frac =
326  (
327  (ts.value() - ts.deltaTValue())
328  - (ts0.value() - ts0.deltaTValue())
329  )/ts0.deltaTValue();
330 
331  if (t0Frac > SMALL)
332  {
333  return V0() + t0Frac*(V() - V0());
334  }
335  else
336  {
337  return V0();
338  }
339  }
340  else
341  {
342  return V0();
343  }
344 }
345 
346 
348 {
349  if (!SfPtr_)
350  {
351  makeSf();
352  }
353 
354  return *SfPtr_;
355 }
356 
357 
359 {
360  if (!magSfPtr_)
361  {
362  makeMagSf();
363  }
364 
365  return *magSfPtr_;
366 }
367 
368 
369 const volVectorField& fvMesh::C() const
370 {
371  if (!CPtr_)
372  {
373  makeC();
374  }
375 
376  return *CPtr_;
377 }
378 
379 
381 {
382  if (!CfPtr_)
383  {
384  makeCf();
385  }
386 
387  return *CfPtr_;
388 }
389 
390 
392 {
393  if (debug)
394  {
395  Info<< "void fvMesh::delta() : "
396  << "calculating face deltas"
397  << endl;
398  }
399 
401  (
403  (
404  IOobject
405  (
406  "delta",
407  pointsInstance(),
408  meshSubDir,
409  *this,
412  false
413  ),
414  *this,
415  dimLength
416  )
417  );
418  surfaceVectorField& delta = tdelta();
419 
420  const volVectorField& C = this->C();
421  const labelUList& owner = this->owner();
422  const labelUList& neighbour = this->neighbour();
423 
424  forAll(owner, facei)
425  {
426  delta[facei] = C[neighbour[facei]] - C[owner[facei]];
427  }
428 
429  forAll(delta.boundaryField(), patchi)
430  {
431  delta.boundaryField()[patchi] = boundary()[patchi].delta();
432  }
433 
434  return tdelta;
435 }
436 
437 
439 {
440  if (!phiPtr_)
441  {
442  FatalErrorIn("fvMesh::phi()")
443  << "mesh flux field does not exist, is the mesh actually moving?"
444  << abort(FatalError);
445  }
446 
447  // Set zero current time
448  // mesh motion fluxes if the time has been incremented
449  if (phiPtr_->timeIndex() != time().timeIndex())
450  {
451  (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
452  }
453 
454  return *phiPtr_;
455 }
456 
457 
459 {
460  if (!phiPtr_)
461  {
462  FatalErrorIn("fvMesh::setPhi()")
463  << "mesh flux field does not exist, is the mesh actually moving?"
464  << abort(FatalError);
465  }
466 
467  return *phiPtr_;
468 }
469 
470 
471 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
472 
473 } // End namespace Foam
474 
475 // ************************************************************************* //
const vectorField & faceAreas() const
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:309
Foam::surfaceFields.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
dimensioned< scalar > mag(const dimensioned< Type > &)
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
const surfaceVectorField & Sf() const
Return cell face area vectors.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
const DimensionedField< scalar, volMesh > & V() const
Return cell volumes.
const surfaceScalarField & phi() const
Return cell face motion fluxes.
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
const DimensionedField< scalar, volMesh > & V00() const
Return old-old-time cell volumes.
const labelUList & owner() const
Internal face owner.
Definition: fvMesh.H:282
DimensionedField< scalar, volMesh > & setV0()
Return old-time cell volumes.
Graphite solid properties.
Definition: C.H:57
SlicedGeometricField< vector, fvsPatchField, slicedFvsPatchField, surfaceMesh > slicedSurfaceVectorField
The time value with time-stepping information, user-defined remapping, etc.
Definition: TimeState.H:49
const vectorField & cellCentres() const
messageStream Info
const surfaceVectorField & Cf() const
Return face centres as surfaceVectorField.
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.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
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
scalar deltaTValue() const
Return time step value.
Definition: TimeState.H:100
const volVectorField & C() const
Return cell centres as volVectorField.
label timeIndex() const
Return current time index.
Definition: TimeState.C:73
#define forAll(list, i)
Definition: UList.H:421
SlicedGeometricField< vector, fvPatchField, slicedFvPatchField, volMesh > slicedVolVectorField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
label patchi
errorManip< error > abort(error &err)
Definition: errorManip.H:131
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
tmp< DimensionedField< scalar, volMesh > > Vsc() const
Return sub-cycle cell volumes.
error FatalError
IOobject(const word &name, const fileName &instance, const objectRegistry &registry, readOption r=NO_READ, writeOption w=NO_WRITE, bool registerObject=true)
Construct from name, instance, registry, io options.
Definition: IOobject.C:131
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:58
surfaceScalarField & setPhi()
Return cell face motion fluxes.
DimensionedField< Type, GeoMesh > DimensionedInternalField
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
label timeIndex() const
Return the time index of the field.
const scalarField & cellVolumes() const
const labelUList & neighbour() const
Internal face neighbour.
Definition: fvMesh.H:288
The internalField of a SlicedGeometricField.
bool moving() const
Is mesh moving.
Definition: polyMesh.H:493
tmp< surfaceVectorField > delta() const
Return face deltas as surfaceVectorField.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
const DimensionedField< scalar, volMesh > & V0() const
Return old-time cell volumes.
tmp< DimensionedField< scalar, volMesh > > Vsc0() const
Return sub-cycl old-time cell volumes.
const Type & value() const
Return const reference to value.
word timeName
Definition: getTimeIndex.H:3
const TimeState & prevTimeState() const
Return previous TimeState if time is being sub-cycled.
Definition: Time.H:459
const vectorField & faceCentres() const
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:552
static int debug
Debug switch.
Definition: fvSchemes.H:103