ENH: autoLayerDriver: better layering information message
[OpenFOAM-2.0.x.git] / src / finiteVolume / fvMesh / fvMesh.C
blob51ded1992f336227d735535fdc977e25f79c955e
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
4    \\    /   O peration     |
5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
6      \\/     M anipulation  |
7 -------------------------------------------------------------------------------
8 License
9     This file is part of OpenFOAM.
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.
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.
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/>.
24 \*---------------------------------------------------------------------------*/
26 #include "fvMesh.H"
27 #include "volFields.H"
28 #include "surfaceFields.H"
29 #include "slicedVolFields.H"
30 #include "slicedSurfaceFields.H"
31 #include "SubField.H"
32 #include "demandDrivenData.H"
33 #include "fvMeshLduAddressing.H"
34 #include "emptyPolyPatch.H"
35 #include "mapPolyMesh.H"
36 #include "MapFvFields.H"
37 #include "fvMeshMapper.H"
38 #include "mapClouds.H"
40 #include "volPointInterpolation.H"
41 #include "extendedLeastSquaresVectors.H"
42 #include "extendedLeastSquaresVectors.H"
43 #include "leastSquaresVectors.H"
44 #include "CentredFitData.H"
45 #include "linearFitPolynomial.H"
46 #include "quadraticFitPolynomial.H"
47 #include "quadraticLinearFitPolynomial.H"
48 //#include "quadraticFitSnGradData.H"
49 #include "skewCorrectionVectors.H"
52 #include "centredCECCellToFaceStencilObject.H"
53 #include "centredCFCCellToFaceStencilObject.H"
54 #include "centredCPCCellToFaceStencilObject.H"
55 #include "centredFECCellToFaceStencilObject.H"
56 #include "upwindCECCellToFaceStencilObject.H"
57 #include "upwindCFCCellToFaceStencilObject.H"
58 #include "upwindCPCCellToFaceStencilObject.H"
59 #include "upwindFECCellToFaceStencilObject.H"
61 #include "centredCFCFaceToCellStencilObject.H"
63 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
65 defineTypeNameAndDebug(Foam::fvMesh, 0);
67 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
69 void Foam::fvMesh::clearGeomNotOldVol()
71     slicedVolScalarField::DimensionedInternalField* VPtr =
72         static_cast<slicedVolScalarField::DimensionedInternalField*>(VPtr_);
73     deleteDemandDrivenData(VPtr);
74     VPtr_ = NULL;
76     deleteDemandDrivenData(SfPtr_);
77     deleteDemandDrivenData(magSfPtr_);
78     deleteDemandDrivenData(CPtr_);
79     deleteDemandDrivenData(CfPtr_);
83 void Foam::fvMesh::clearGeom()
85     clearGeomNotOldVol();
87     deleteDemandDrivenData(V0Ptr_);
88     deleteDemandDrivenData(V00Ptr_);
90     // Mesh motion flux cannot be deleted here because the old-time flux
91     // needs to be saved.
93     // Things geometry dependent that are not updated.
94     volPointInterpolation::Delete(*this);
95     extendedLeastSquaresVectors::Delete(*this);
96     leastSquaresVectors::Delete(*this);
97     CentredFitData<linearFitPolynomial>::Delete(*this);
98     CentredFitData<quadraticFitPolynomial>::Delete(*this);
99     CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
100     skewCorrectionVectors::Delete(*this);
101     //quadraticFitSnGradData::Delete(*this);
105 void Foam::fvMesh::clearAddressing()
107     deleteDemandDrivenData(lduPtr_);
109     // Hack until proper callbacks. Below are all the fvMesh-MeshObjects.
111     volPointInterpolation::Delete(*this);
112     extendedLeastSquaresVectors::Delete(*this);
113     leastSquaresVectors::Delete(*this);
114     CentredFitData<linearFitPolynomial>::Delete(*this);
115     CentredFitData<quadraticFitPolynomial>::Delete(*this);
116     CentredFitData<quadraticLinearFitPolynomial>::Delete(*this);
117     skewCorrectionVectors::Delete(*this);
118     //quadraticFitSnGradData::Delete(*this);
120     centredCECCellToFaceStencilObject::Delete(*this);
121     centredCFCCellToFaceStencilObject::Delete(*this);
122     centredCPCCellToFaceStencilObject::Delete(*this);
123     centredFECCellToFaceStencilObject::Delete(*this);
124     // Is this geometry related - cells distorting to upwind direction?
125     upwindCECCellToFaceStencilObject::Delete(*this);
126     upwindCFCCellToFaceStencilObject::Delete(*this);
127     upwindCPCCellToFaceStencilObject::Delete(*this);
128     upwindFECCellToFaceStencilObject::Delete(*this);
130     centredCFCFaceToCellStencilObject::Delete(*this);
134 void Foam::fvMesh::clearOut()
136     clearGeom();
137     surfaceInterpolation::clearOut();
139     clearAddressing();
141     // Clear mesh motion flux
142     deleteDemandDrivenData(phiPtr_);
144     polyMesh::clearOut();
148 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
150 Foam::fvMesh::fvMesh(const IOobject& io)
152     polyMesh(io),
153     surfaceInterpolation(*this),
154     fvSchemes(static_cast<const objectRegistry&>(*this)),
155     fvSolution(static_cast<const objectRegistry&>(*this)),
156     data(static_cast<const objectRegistry&>(*this)),
157     boundary_(*this, boundaryMesh()),
158     lduPtr_(NULL),
159     curTimeIndex_(time().timeIndex()),
160     VPtr_(NULL),
161     V0Ptr_(NULL),
162     V00Ptr_(NULL),
163     SfPtr_(NULL),
164     magSfPtr_(NULL),
165     CPtr_(NULL),
166     CfPtr_(NULL),
167     phiPtr_(NULL)
169     if (debug)
170     {
171         Info<< "Constructing fvMesh from IOobject"
172             << endl;
173     }
175     // Check the existance of the cell volumes and read if present
176     // and set the storage of V00
177     if (isFile(time().timePath()/"V0"))
178     {
179         V0Ptr_ = new DimensionedField<scalar, volMesh>
180         (
181             IOobject
182             (
183                 "V0",
184                 time().timeName(),
185                 *this,
186                 IOobject::MUST_READ,
187                 IOobject::NO_WRITE
188             ),
189             *this
190         );
192         V00();
193     }
195     // Check the existance of the mesh fluxes, read if present and set the
196     // mesh to be moving
197     if (isFile(time().timePath()/"meshPhi"))
198     {
199         phiPtr_ = new surfaceScalarField
200         (
201             IOobject
202             (
203                 "meshPhi",
204                 time().timeName(),
205                 *this,
206                 IOobject::MUST_READ,
207                 IOobject::AUTO_WRITE
208             ),
209             *this
210         );
212         // The mesh is now considered moving so the old-time cell volumes
213         // will be required for the time derivatives so if they haven't been
214         // read initialise to the current cell volumes
215         if (!V0Ptr_)
216         {
217             V0Ptr_ = new DimensionedField<scalar, volMesh>
218             (
219                 IOobject
220                 (
221                     "V0",
222                     time().timeName(),
223                     *this,
224                     IOobject::NO_READ,
225                     IOobject::NO_WRITE,
226                     false
227                 ),
228                 V()
229             );
230         }
232         moving(true);
233     }
237 Foam::fvMesh::fvMesh
239     const IOobject& io,
240     const Xfer<pointField>& points,
241     const Xfer<faceList>& faces,
242     const Xfer<labelList>& allOwner,
243     const Xfer<labelList>& allNeighbour,
244     const bool syncPar
247     polyMesh(io, points, faces, allOwner, allNeighbour, syncPar),
248     surfaceInterpolation(*this),
249     fvSchemes(static_cast<const objectRegistry&>(*this)),
250     fvSolution(static_cast<const objectRegistry&>(*this)),
251     data(static_cast<const objectRegistry&>(*this)),
252     boundary_(*this),
253     lduPtr_(NULL),
254     curTimeIndex_(time().timeIndex()),
255     VPtr_(NULL),
256     V0Ptr_(NULL),
257     V00Ptr_(NULL),
258     SfPtr_(NULL),
259     magSfPtr_(NULL),
260     CPtr_(NULL),
261     CfPtr_(NULL),
262     phiPtr_(NULL)
264     if (debug)
265     {
266         Info<< "Constructing fvMesh from components" << endl;
267     }
271 Foam::fvMesh::fvMesh
273     const IOobject& io,
274     const Xfer<pointField>& points,
275     const Xfer<faceList>& faces,
276     const Xfer<cellList>& cells,
277     const bool syncPar
280     polyMesh(io, points, faces, cells, syncPar),
281     surfaceInterpolation(*this),
282     fvSchemes(static_cast<const objectRegistry&>(*this)),
283     fvSolution(static_cast<const objectRegistry&>(*this)),
284     data(static_cast<const objectRegistry&>(*this)),
285     boundary_(*this),
286     lduPtr_(NULL),
287     curTimeIndex_(time().timeIndex()),
288     VPtr_(NULL),
289     V0Ptr_(NULL),
290     V00Ptr_(NULL),
291     SfPtr_(NULL),
292     magSfPtr_(NULL),
293     CPtr_(NULL),
294     CfPtr_(NULL),
295     phiPtr_(NULL)
297     if (debug)
298     {
299         Info<< "Constructing fvMesh from components" << endl;
300     }
304 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
306 Foam::fvMesh::~fvMesh()
308     clearOut();
312 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
314 void Foam::fvMesh::addFvPatches
316     const List<polyPatch*> & p,
317     const bool validBoundary
320     if (boundary().size())
321     {
322         FatalErrorIn
323         (
324             "fvMesh::addFvPatches(const List<polyPatch*>&, const bool)"
325         )   << " boundary already exists"
326             << abort(FatalError);
327     }
329     // first add polyPatches
330     addPatches(p, validBoundary);
331     boundary_.addPatches(boundaryMesh());
335 void Foam::fvMesh::removeFvBoundary()
337     if (debug)
338     {
339         Info<< "void fvMesh::removeFvBoundary(): "
340             << "Removing boundary patches."
341             << endl;
342     }
344     // Remove fvBoundaryMesh data first.
345     boundary_.clear();
346     boundary_.setSize(0);
347     polyMesh::removeBoundary();
349     clearOut();
353 Foam::polyMesh::readUpdateState Foam::fvMesh::readUpdate()
355     if (debug)
356     {
357         Info<< "polyMesh::readUpdateState fvMesh::readUpdate() : "
358             << "Updating fvMesh.  ";
359     }
361     polyMesh::readUpdateState state = polyMesh::readUpdate();
363     if (state == polyMesh::TOPO_PATCH_CHANGE)
364     {
365         if (debug)
366         {
367             Info<< "Boundary and topological update" << endl;
368         }
370         boundary_.readUpdate(boundaryMesh());
372         clearOut();
374     }
375     else if (state == polyMesh::TOPO_CHANGE)
376     {
377         if (debug)
378         {
379             Info<< "Topological update" << endl;
380         }
382         clearOut();
383     }
384     else if (state == polyMesh::POINTS_MOVED)
385     {
386         if (debug)
387         {
388             Info<< "Point motion update" << endl;
389         }
391         clearGeom();
392     }
393     else
394     {
395         if (debug)
396         {
397             Info<< "No update" << endl;
398         }
399     }
401     return state;
405 const Foam::fvBoundaryMesh& Foam::fvMesh::boundary() const
407     return boundary_;
411 const Foam::lduAddressing& Foam::fvMesh::lduAddr() const
413     if (!lduPtr_)
414     {
415         lduPtr_ = new fvMeshLduAddressing(*this);
416     }
418     return *lduPtr_;
422 void Foam::fvMesh::mapFields(const mapPolyMesh& meshMap)
424     // Create a mapper
425     const fvMeshMapper mapper(*this, meshMap);
427     // Map all the volFields in the objectRegistry
428     MapGeometricFields<scalar, fvPatchField, fvMeshMapper, volMesh>
429     (mapper);
430     MapGeometricFields<vector, fvPatchField, fvMeshMapper, volMesh>
431     (mapper);
432     MapGeometricFields<sphericalTensor, fvPatchField, fvMeshMapper, volMesh>
433     (mapper);
434     MapGeometricFields<symmTensor, fvPatchField, fvMeshMapper, volMesh>
435     (mapper);
436     MapGeometricFields<tensor, fvPatchField, fvMeshMapper, volMesh>
437     (mapper);
439     // Map all the surfaceFields in the objectRegistry
440     MapGeometricFields<scalar, fvsPatchField, fvMeshMapper, surfaceMesh>
441     (mapper);
442     MapGeometricFields<vector, fvsPatchField, fvMeshMapper, surfaceMesh>
443     (mapper);
444     MapGeometricFields<symmTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
445     (mapper);
446     MapGeometricFields<symmTensor, fvsPatchField, fvMeshMapper, surfaceMesh>
447     (mapper);
448     MapGeometricFields<tensor, fvsPatchField, fvMeshMapper, surfaceMesh>
449     (mapper);
451     // Map all the clouds in the objectRegistry
452     mapClouds(*this, meshMap);
455     const labelList& cellMap = meshMap.cellMap();
457     // Map the old volume. Just map to new cell labels.
458     if (V0Ptr_)
459     {
460         scalarField& V0 = *V0Ptr_;
462         scalarField savedV0(V0);
463         V0.setSize(nCells());
465         forAll(V0, i)
466         {
467             if (cellMap[i] > -1)
468             {
469                 V0[i] = savedV0[cellMap[i]];
470             }
471             else
472             {
473                 V0[i] = 0.0;
474             }
475         }
476     }
478     // Map the old-old volume. Just map to new cell labels.
479     if (V00Ptr_)
480     {
481         scalarField& V00 = *V00Ptr_;
483         scalarField savedV00(V00);
484         V00.setSize(nCells());
486         forAll(V00, i)
487         {
488             if (cellMap[i] > -1)
489             {
490                 V00[i] = savedV00[cellMap[i]];
491             }
492             else
493             {
494                 V00[i] = 0.0;
495             }
496         }
497     }
501 // Temporary helper function to call move points on
502 // MeshObjects
503 template<class Type>
504 void MeshObjectMovePoints(const Foam::fvMesh& mesh)
506     if (mesh.thisDb().foundObject<Type>(Type::typeName))
507     {
508         const_cast<Type&>
509         (
510             mesh.thisDb().lookupObject<Type>
511             (
512                 Type::typeName
513             )
514         ).movePoints();
515     }
519 Foam::tmp<Foam::scalarField> Foam::fvMesh::movePoints(const pointField& p)
521     // Grab old time volumes if the time has been incremented
522     if (curTimeIndex_ < time().timeIndex())
523     {
524         if (V00Ptr_ && V0Ptr_)
525         {
526             *V00Ptr_ = *V0Ptr_;
527         }
529         if (V0Ptr_)
530         {
531             *V0Ptr_ = V();
532         }
533         else
534         {
535             V0Ptr_ = new DimensionedField<scalar, volMesh>
536             (
537                 IOobject
538                 (
539                     "V0",
540                     time().timeName(),
541                     *this,
542                     IOobject::NO_READ,
543                     IOobject::NO_WRITE
544                 ),
545                 V()
546             );
547         }
549         curTimeIndex_ = time().timeIndex();
550     }
553     // delete out of date geometrical information
554     clearGeomNotOldVol();
557     if (!phiPtr_)
558     {
559         // Create mesh motion flux
560         phiPtr_ = new surfaceScalarField
561         (
562             IOobject
563             (
564                 "meshPhi",
565                 this->time().timeName(),
566                 *this,
567                 IOobject::NO_READ,
568                 IOobject::AUTO_WRITE
569             ),
570             *this,
571             dimVolume/dimTime
572         );
573     }
574     else
575     {
576         // Grab old time mesh motion fluxes if the time has been incremented
577         if (phiPtr_->timeIndex() != time().timeIndex())
578         {
579             phiPtr_->oldTime();
580         }
581     }
583     surfaceScalarField& phi = *phiPtr_;
585     // Move the polyMesh and set the mesh motion fluxes to the swept-volumes
587     scalar rDeltaT = 1.0/time().deltaTValue();
589     tmp<scalarField> tsweptVols = polyMesh::movePoints(p);
590     scalarField& sweptVols = tsweptVols();
592     phi.internalField() = scalarField::subField(sweptVols, nInternalFaces());
593     phi.internalField() *= rDeltaT;
595     const fvPatchList& patches = boundary();
597     forAll(patches, patchI)
598     {
599         phi.boundaryField()[patchI] = patches[patchI].patchSlice(sweptVols);
600         phi.boundaryField()[patchI] *= rDeltaT;
601     }
603     boundary_.movePoints();
604     surfaceInterpolation::movePoints();
607     // Hack until proper callbacks. Below are all the fvMesh MeshObjects with a
608     // movePoints function.
609     MeshObjectMovePoints<volPointInterpolation>(*this);
610     MeshObjectMovePoints<extendedLeastSquaresVectors>(*this);
611     MeshObjectMovePoints<leastSquaresVectors>(*this);
612     MeshObjectMovePoints<CentredFitData<linearFitPolynomial> >(*this);
613     MeshObjectMovePoints<CentredFitData<quadraticFitPolynomial> >(*this);
614     MeshObjectMovePoints<CentredFitData<quadraticLinearFitPolynomial> >(*this);
615     MeshObjectMovePoints<skewCorrectionVectors>(*this);
616     //MeshObjectMovePoints<quadraticFitSnGradData>(*this);
618     return tsweptVols;
622 void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm)
624     // Update polyMesh. This needs to keep volume existent!
625     polyMesh::updateMesh(mpm);
627     // Clear the sliced fields
628     clearGeomNotOldVol();
630     // Map all fields using current (i.e. not yet mapped) volume
631     mapFields(mpm);
633     // Clear the current volume and other geometry factors
634     surfaceInterpolation::clearOut();
636     clearAddressing();
638     // handleMorph() should also clear out the surfaceInterpolation.
639     // This is a temporary solution
640     surfaceInterpolation::movePoints();
644 bool Foam::fvMesh::writeObjects
646     IOstream::streamFormat fmt,
647     IOstream::versionNumber ver,
648     IOstream::compressionType cmp
649 ) const
651     return polyMesh::writeObject(fmt, ver, cmp);
655 //- Write mesh using IO settings from the time
656 bool Foam::fvMesh::write() const
658     return polyMesh::write();
662 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
664 bool Foam::fvMesh::operator!=(const fvMesh& bm) const
666     return &bm != this;
670 bool Foam::fvMesh::operator==(const fvMesh& bm) const
672     return &bm == this;
676 // ************************************************************************* //