FreeFOAM The Cross-Platform CFD Toolkit
fvMesh.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 "fvMesh.H"
27 #include <finiteVolume/volFields.H>
31 #include <OpenFOAM/SubField.H>
33 #include "fvMeshLduAddressing.H"
35 #include <OpenFOAM/mapPolyMesh.H>
38 #include <OpenFOAM/mapClouds.H>
39 
48 //#include "quadraticFitSnGradData.H"
50 
51 
60 
62 
63 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
64 
66 
67 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
68 
69 void Foam::fvMesh::clearGeomNotOldVol()
70 {
74  VPtr_ = NULL;
75 
76  deleteDemandDrivenData(SfPtr_);
77  deleteDemandDrivenData(magSfPtr_);
79  deleteDemandDrivenData(CfPtr_);
80 }
81 
82 
83 void Foam::fvMesh::clearGeom()
84 {
85  clearGeomNotOldVol();
86 
87  deleteDemandDrivenData(V0Ptr_);
88  deleteDemandDrivenData(V00Ptr_);
89 
90  // Mesh motion flux cannot be deleted here because the old-time flux
91  // needs to be saved.
92 
93  // Things geometry dependent that are not updated.
101  //quadraticFitSnGradData::Delete(*this);
102 }
103 
104 
105 void Foam::fvMesh::clearAddressing()
106 {
107  deleteDemandDrivenData(lduPtr_);
108 
109  // Hack until proper callbacks. Below are all the fvMesh-MeshObjects.
110 
118  //quadraticFitSnGradData::Delete(*this);
119 
124  // Is this geometry related - cells distorting to upwind direction?
129 
131 }
132 
133 
135 {
136  clearGeom();
138 
139  clearAddressing();
140 
141  // Clear mesh motion flux
142  deleteDemandDrivenData(phiPtr_);
143 
145 }
146 
147 
148 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
149 
150 Foam::fvMesh::fvMesh(const IOobject& io)
151 :
152  polyMesh(io),
153  surfaceInterpolation(*this),
154  boundary_(*this, boundaryMesh()),
155  lduPtr_(NULL),
156  curTimeIndex_(time().timeIndex()),
157  VPtr_(NULL),
158  V0Ptr_(NULL),
159  V00Ptr_(NULL),
160  SfPtr_(NULL),
161  magSfPtr_(NULL),
162  CPtr_(NULL),
163  CfPtr_(NULL),
164  phiPtr_(NULL)
165 {
166  if (debug)
167  {
168  Info<< "Constructing fvMesh from IOobject"
169  << endl;
170  }
171 
172  // Check the existance of the cell volumes and read if present
173  // and set the storage of V00
174  if (isFile(time().timePath()/"V0"))
175  {
177  (
178  IOobject
179  (
180  "V0",
181  time().timeName(),
182  *this,
185  ),
186  *this
187  );
188 
189  V00();
190  }
191 
192  // Check the existance of the mesh fluxes, read if present and set the
193  // mesh to be moving
194  if (isFile(time().timePath()/"meshPhi"))
195  {
196  phiPtr_ = new surfaceScalarField
197  (
198  IOobject
199  (
200  "meshPhi",
201  time().timeName(),
202  *this,
205  ),
206  *this
207  );
208 
209  // The mesh is now considered moving so the old-time cell volumes
210  // will be required for the time derivatives so if they haven't been
211  // read initialise to the current cell volumes
212  if (!V0Ptr_)
213  {
215  (
216  IOobject
217  (
218  "V0",
219  time().timeName(),
220  *this,
223  false
224  ),
225  V()
226  );
227  }
228 
229  moving(true);
230  }
231 }
232 
233 
234 Foam::fvMesh::fvMesh
235 (
236  const IOobject& io,
237  const Xfer<pointField>& points,
238  const Xfer<faceList>& faces,
239  const Xfer<labelList>& allOwner,
240  const Xfer<labelList>& allNeighbour,
241  const bool syncPar
242 )
243 :
244  polyMesh(io, points, faces, allOwner, allNeighbour, syncPar),
245  surfaceInterpolation(*this),
246  boundary_(*this),
247  lduPtr_(NULL),
248  curTimeIndex_(time().timeIndex()),
249  VPtr_(NULL),
250  V0Ptr_(NULL),
251  V00Ptr_(NULL),
252  SfPtr_(NULL),
253  magSfPtr_(NULL),
254  CPtr_(NULL),
255  CfPtr_(NULL),
256  phiPtr_(NULL)
257 {
258  if (debug)
259  {
260  Info<< "Constructing fvMesh from components" << endl;
261  }
262 }
263 
264 
265 Foam::fvMesh::fvMesh
266 (
267  const IOobject& io,
268  const Xfer<pointField>& points,
269  const Xfer<faceList>& faces,
270  const Xfer<cellList>& cells,
271  const bool syncPar
272 )
273 :
274  polyMesh(io, points, faces, cells, syncPar),
275  surfaceInterpolation(*this),
276  boundary_(*this),
277  lduPtr_(NULL),
278  curTimeIndex_(time().timeIndex()),
279  VPtr_(NULL),
280  V0Ptr_(NULL),
281  V00Ptr_(NULL),
282  SfPtr_(NULL),
283  magSfPtr_(NULL),
284  CPtr_(NULL),
285  CfPtr_(NULL),
286  phiPtr_(NULL)
287 {
288  if (debug)
289  {
290  Info<< "Constructing fvMesh from components" << endl;
291  }
292 }
293 
294 
295 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
296 
298 {
299  clearOut();
300 }
301 
302 
303 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
304 
306 (
307  const List<polyPatch*> & p,
308  const bool validBoundary
309 )
310 {
311  if (boundary().size())
312  {
314  (
315  "fvMesh::addFvPatches(const List<polyPatch*>&, const bool)"
316  ) << " boundary already exists"
317  << abort(FatalError);
318  }
319 
320  // first add polyPatches
321  addPatches(p, validBoundary);
322  boundary_.addPatches(boundaryMesh());
323 }
324 
325 
327 {
328  if (debug)
329  {
330  Info<< "void fvMesh::removeFvBoundary(): "
331  << "Removing boundary patches."
332  << endl;
333  }
334 
335  // Remove fvBoundaryMesh data first.
336  boundary_.clear();
337  boundary_.setSize(0);
339 
340  clearOut();
341 }
342 
343 
345 {
346  if (debug)
347  {
348  Info<< "polyMesh::readUpdateState fvMesh::readUpdate() : "
349  << "Updating fvMesh. ";
350  }
351 
353 
354  if (state == polyMesh::TOPO_PATCH_CHANGE)
355  {
356  if (debug)
357  {
358  Info << "Boundary and topological update" << endl;
359  }
360 
361  boundary_.readUpdate(boundaryMesh());
362 
363  clearOut();
364 
365  }
366  else if (state == polyMesh::TOPO_CHANGE)
367  {
368  if (debug)
369  {
370  Info << "Topological update" << endl;
371  }
372 
373  clearOut();
374  }
375  else if (state == polyMesh::POINTS_MOVED)
376  {
377  if (debug)
378  {
379  Info << "Point motion update" << endl;
380  }
381 
382  clearGeom();
383  }
384  else
385  {
386  if (debug)
387  {
388  Info << "No update" << endl;
389  }
390  }
391 
392  return state;
393 }
394 
395 
397 {
398  return boundary_;
399 }
400 
401 
403 {
404  if (!lduPtr_)
405  {
406  lduPtr_ = new fvMeshLduAddressing(*this);
407  }
408 
409  return *lduPtr_;
410 }
411 
412 
414 {
415  // Create a mapper
416  const fvMeshMapper mapper(*this, meshMap);
417 
418  // Map all the volFields in the objectRegistry
419  MapGeometricFields<scalar, fvPatchField, fvMeshMapper, volMesh>
420  (mapper);
421  MapGeometricFields<vector, fvPatchField, fvMeshMapper, volMesh>
422  (mapper);
423  MapGeometricFields<sphericalTensor, fvPatchField, fvMeshMapper, volMesh>
424  (mapper);
425  MapGeometricFields<symmTensor, fvPatchField, fvMeshMapper, volMesh>
426  (mapper);
427  MapGeometricFields<tensor, fvPatchField, fvMeshMapper, volMesh>
428  (mapper);
429 
430  // Map all the surfaceFields in the objectRegistry
431  MapGeometricFields<scalar, fvsPatchField, fvMeshMapper, surfaceMesh>
432  (mapper);
433  MapGeometricFields<vector, fvsPatchField, fvMeshMapper, surfaceMesh>
434  (mapper);
435  MapGeometricFields<symmTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
436  (mapper);
437  MapGeometricFields<symmTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
438  (mapper);
439  MapGeometricFields<tensor, fvsPatchField, fvMeshMapper, surfaceMesh>
440  (mapper);
441 
442  // Map all the clouds in the objectRegistry
443  mapClouds(*this, meshMap);
444 
445 
446  const labelList& cellMap = meshMap.cellMap();
447 
448  // Map the old volume. Just map to new cell labels.
449  if (V0Ptr_)
450  {
451  scalarField& V0 = *V0Ptr_;
452 
453  scalarField savedV0(V0);
454  V0.setSize(nCells());
455 
456  forAll(V0, i)
457  {
458  if (cellMap[i] > -1)
459  {
460  V0[i] = savedV0[cellMap[i]];
461  }
462  else
463  {
464  V0[i] = 0.0;
465  }
466  }
467  }
468 
469  // Map the old-old volume. Just map to new cell labels.
470  if (V00Ptr_)
471  {
472  scalarField& V00 = *V00Ptr_;
473 
474  scalarField savedV00(V00);
475  V00.setSize(nCells());
476 
477  forAll(V00, i)
478  {
479  if (cellMap[i] > -1)
480  {
481  V00[i] = savedV00[cellMap[i]];
482  }
483  else
484  {
485  V00[i] = 0.0;
486  }
487  }
488  }
489 }
490 
491 
492 // Temporary helper function to call move points on
493 // MeshObjects
494 template<class Type>
496 {
497  if (mesh.thisDb().foundObject<Type>(Type::typeName))
498  {
499  const_cast<Type&>
500  (
501  mesh.thisDb().lookupObject<Type>
502  (
503  Type::typeName
504  )
505  ).movePoints();
506  }
507 }
508 
509 
511 {
512  // Grab old time volumes if the time has been incremented
513  if (curTimeIndex_ < time().timeIndex())
514  {
515  if (V00Ptr_ && V0Ptr_)
516  {
517  *V00Ptr_ = *V0Ptr_;
518  }
519 
520  if (V0Ptr_)
521  {
522  *V0Ptr_ = V();
523  }
524  else
525  {
527  (
528  IOobject
529  (
530  "V0",
531  time().timeName(),
532  *this,
535  ),
536  V()
537  );
538  }
539 
540  curTimeIndex_ = time().timeIndex();
541  }
542 
543 
544  // delete out of date geometrical information
545  clearGeomNotOldVol();
546 
547 
548  if (!phiPtr_)
549  {
550  // Create mesh motion flux
551  phiPtr_ = new surfaceScalarField
552  (
553  IOobject
554  (
555  "meshPhi",
556  this->time().timeName(),
557  *this,
560  ),
561  *this,
563  );
564  }
565  else
566  {
567  // Grab old time mesh motion fluxes if the time has been incremented
568  if (phiPtr_->timeIndex() != time().timeIndex())
569  {
570  phiPtr_->oldTime();
571  }
572  }
573 
574  surfaceScalarField& phi = *phiPtr_;
575 
576  // Move the polyMesh and set the mesh motion fluxes to the swept-volumes
577 
578  scalar rDeltaT = 1.0/time().deltaT().value();
579 
580  tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
581  scalarField& sweptVols = tsweptVols();
582 
583  phi.internalField() = scalarField::subField(sweptVols, nInternalFaces());
584  phi.internalField() *= rDeltaT;
585 
586  const fvPatchList& patches = boundary();
587 
588  forAll (patches, patchI)
589  {
590  phi.boundaryField()[patchI] = patches[patchI].patchSlice(sweptVols);
591  phi.boundaryField()[patchI] *= rDeltaT;
592  }
593 
594  boundary_.movePoints();
596 
597 
598  // Hack until proper callbacks. Below are all the fvMesh MeshObjects with a
599  // movePoints function.
600  MeshObjectMovePoints<volPointInterpolation>(*this);
601  MeshObjectMovePoints<extendedLeastSquaresVectors>(*this);
602  MeshObjectMovePoints<leastSquaresVectors>(*this);
603  MeshObjectMovePoints<CentredFitData<linearFitPolynomial> >(*this);
604  MeshObjectMovePoints<CentredFitData<quadraticFitPolynomial> >(*this);
605  MeshObjectMovePoints<CentredFitData<quadraticLinearFitPolynomial> >(*this);
606  MeshObjectMovePoints<skewCorrectionVectors>(*this);
607  //MeshObjectMovePoints<quadraticFitSnGradData>(*this);
608 
609  return tsweptVols;
610 }
611 
612 
614 {
615  // Update polyMesh. This needs to keep volume existent!
617 
618  // Clear the sliced fields
619  clearGeomNotOldVol();
620 
621  // Map all fields using current (i.e. not yet mapped) volume
622  mapFields(mpm);
623 
624  // Clear the current volume and other geometry factors
626 
627  clearAddressing();
628 
629  // handleMorph() should also clear out the surfaceInterpolation.
630  // This is a temporary solution
632 }
633 
634 
636 (
640 ) const
641 {
642  return polyMesh::writeObject(fmt, ver, cmp);
643 }
644 
645 
646 //- Write mesh using IO settings from the time
648 {
649  return polyMesh::write();
650 }
651 
652 
653 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
654 
655 bool Foam::fvMesh::operator!=(const fvMesh& bm) const
656 {
657  return &bm != this;
658 }
659 
660 
661 bool Foam::fvMesh::operator==(const fvMesh& bm) const
662 {
663  return &bm == this;
664 }
665 
666 
667 // ************************ vim: set sw=4 sts=4 et: ************************ //