FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2011 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 "fvMesh.H"
27 #include <OpenFOAM/Time.H>
28 #include <finiteVolume/volFields.H>
32 #include <OpenFOAM/SubField.H>
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void fvMesh::makeSf() const
43 {
44  if (debug)
45  {
46  Info<< "void fvMesh::makeSf() : "
47  << "assembling face areas"
48  << endl;
49  }
50 
51  // It is an error to attempt to recalculate
52  // if the pointer is already set
53  if (SfPtr_)
54  {
55  FatalErrorIn("fvMesh::makeSf()")
56  << "face areas already exist"
57  << abort(FatalError);
58  }
59 
60  SfPtr_ = new slicedSurfaceVectorField
61  (
62  IOobject
63  (
64  "S",
66  meshSubDir,
67  *this,
70  false
71  ),
72  *this,
73  dimArea,
74  faceAreas()
75  );
76 }
77 
78 
79 void fvMesh::makeMagSf() const
80 {
81  if (debug)
82  {
83  Info<< "void fvMesh::makeMagSf() : "
84  << "assembling mag face areas"
85  << endl;
86  }
87 
88  // It is an error to attempt to recalculate
89  // if the pointer is already set
90  if (magSfPtr_)
91  {
92  FatalErrorIn("void fvMesh::makeMagSf()")
93  << "mag face areas already exist"
94  << abort(FatalError);
95  }
96 
97  // Note: Added stabilisation for faces with exactly zero area.
98  // These should be caught on mesh checking but at least this stops
99  // the code from producing Nans.
100  magSfPtr_ = new surfaceScalarField
101  (
102  IOobject
103  (
104  "magSf",
105  pointsInstance(),
106  meshSubDir,
107  *this,
110  false
111  ),
112  mag(Sf()) + dimensionedScalar("vs", dimArea, VSMALL)
113  );
114 }
115 
116 
117 void fvMesh::makeC() const
118 {
119  if (debug)
120  {
121  Info<< "void fvMesh::makeC() : "
122  << "assembling cell centres"
123  << endl;
124  }
125 
126  // It is an error to attempt to recalculate
127  // if the pointer is already set
128  if (CPtr_)
129  {
130  FatalErrorIn("fvMesh::makeC()")
131  << "cell centres already exist"
132  << abort(FatalError);
133  }
134 
135  CPtr_ = new slicedVolVectorField
136  (
137  IOobject
138  (
139  "C",
140  pointsInstance(),
141  meshSubDir,
142  *this,
145  false
146  ),
147  *this,
148  dimLength,
149  cellCentres(),
150  faceCentres()
151  );
152 
153 
154  // Need to correct for cyclics transformation since absolute quantity.
155  // Ok on processor patches since hold opposite cell centre (no
156  // transformation)
157  slicedVolVectorField& C = *CPtr_;
158 
159  forAll(C.boundaryField(), patchi)
160  {
161  if (isA<cyclicFvPatchVectorField>(C.boundaryField()[patchi]))
162  {
163  // Note: cyclic is not slice but proper field
164  C.boundaryField()[patchi] == static_cast<const vectorField&>
165  (
166  static_cast<const List<vector>&>
167  (
168  boundary_[patchi].patchSlice(faceCentres())
169  )
170  );
171  }
172  }
173 }
174 
175 
176 void fvMesh::makeCf() const
177 {
178  if (debug)
179  {
180  Info<< "void fvMesh::makeCf() : "
181  << "assembling face centres"
182  << endl;
183  }
184 
185  // It is an error to attempt to recalculate
186  // if the pointer is already set
187  if (CfPtr_)
188  {
189  FatalErrorIn("fvMesh::makeCf()")
190  << "face centres already exist"
191  << abort(FatalError);
192  }
193 
194  CfPtr_ = new slicedSurfaceVectorField
195  (
196  IOobject
197  (
198  "Cf",
199  pointsInstance(),
200  meshSubDir,
201  *this,
204  false
205  ),
206  *this,
207  dimLength,
208  faceCentres()
209  );
210 }
211 
212 
213 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
214 
216 {
217  if (!VPtr_)
218  {
220  (
221  IOobject
222  (
223  "V",
224  time().timeName(),
225  *this,
228  ),
229  *this,
230  dimVolume,
231  cellVolumes()
232  );
233  }
234 
235  return *static_cast<slicedVolScalarField::DimensionedInternalField*>(VPtr_);
236 }
237 
238 
240 {
241  if (!V0Ptr_)
242  {
243  FatalErrorIn("fvMesh::V0() const")
244  << "V0 is not available"
245  << abort(FatalError);
246  }
247 
248  return *V0Ptr_;
249 }
250 
251 
253 {
254  if (!V0Ptr_)
255  {
256  FatalErrorIn("fvMesh::setV0()")
257  << "V0 is not available"
258  << abort(FatalError);
259  }
260 
261  return *V0Ptr_;
262 }
263 
264 
266 {
267  if (!V00Ptr_)
268  {
270  (
271  IOobject
272  (
273  "V00",
274  time().timeName(),
275  *this,
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 (!phiPtr_)
394  {
395  FatalErrorIn("fvMesh::phi()")
396  << "mesh flux field does not exists, is the mesh actually moving?"
397  << exit(FatalError);
398  }
399 
400  // Set zero current time
401  // mesh motion fluxes if the time has been incremented
402  if (phiPtr_->timeIndex() != time().timeIndex())
403  {
404  (*phiPtr_) = dimensionedScalar("0", dimVolume/dimTime, 0.0);
405  }
406 
407  return *phiPtr_;
408 }
409 
410 
412 {
413  if (!phiPtr_)
414  {
415  FatalErrorIn("fvMesh::setPhi()")
416  << "mesh flux field does not exists, is the mesh actually moving?"
417  << exit(FatalError);
418  }
419 
420  return *phiPtr_;
421 }
422 
423 
424 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
425 
426 } // End namespace Foam
427 
428 // ************************ vim: set sw=4 sts=4 et: ************************ //