src/regionModels/regionModel/regionModel1D/regionModel1D.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2010-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 "regionModel1D.H"
  27
  28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  29
  30 namespace Foam
  31 {
  32 namespace regionModels
  33 {
  34     defineTypeNameAndDebug(regionModel1D, 0);
  35 }
  36 }
  37
  38 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  39
  40 void Foam::regionModels::regionModel1D::constructMeshObjects()
  41 {
  42     const fvMesh& regionMesh = regionMeshPtr_();
  43
  44     nMagSfPtr_.reset
  45     (
  46         new surfaceScalarField
  47         (
  48             IOobject
  49             (
  50                 "nMagSf",
  51                 time().timeName(),
  52                 regionMesh,
  53                 IOobject::NO_READ,
  54                 IOobject::NO_WRITE
  55             ),
  56             regionMesh,
  57             dimensionedScalar("zero", dimArea, 0.0)
  58         )
  59     );
  60 }
  61
  62
  63 void Foam::regionModels::regionModel1D::initialise()
  64 {
  65     if (debug)
  66     {
  67         Pout<< "regionModel1D::initialise()" << endl;
  68     }
  69
  70     // Calculate boundaryFaceFaces and boundaryFaceCells
  71
  72     DynamicList<label> faceIDs;
  73     DynamicList<label> cellIDs;
  74
  75     label localPyrolysisFaceI = 0;
  76
  77     const polyBoundaryMesh& rbm = regionMesh().boundaryMesh();
  78
  79     forAll(intCoupledPatchIDs_, i)
  80     {
  81         const label patchI = intCoupledPatchIDs_[i];
  82         const polyPatch& ppCoupled = rbm[patchI];
  83         forAll(ppCoupled, localFaceI)
  84         {
  85             label faceI = ppCoupled.start() + localFaceI;
  86             label cellI = -1;
  87             label nFaces = 0;
  88             label nCells = 0;
  89             do
  90             {
  91                 label ownCellI = regionMesh().faceOwner()[faceI];
  92                 if (ownCellI != cellI)
  93                 {
  94                     cellI = ownCellI;
  95                 }
  96                 else
  97                 {
  98                     cellI = regionMesh().faceNeighbour()[faceI];
  99                 }
 100                 nCells++;
 101                 cellIDs.append(cellI);
 102                 const cell& cFaces = regionMesh().cells()[cellI];
 103                 faceI = cFaces.opposingFaceLabel(faceI, regionMesh().faces());
 104                 faceIDs.append(faceI);
 105                 nFaces++;
 106             } while (regionMesh().isInternalFace(faceI));
 107
 108             boundaryFaceOppositeFace_[localPyrolysisFaceI] = faceI;
 109             faceIDs.remove(); //remove boundary face.
 110             nFaces--;
 111
 112             boundaryFaceFaces_[localPyrolysisFaceI].transfer(faceIDs);
 113             boundaryFaceCells_[localPyrolysisFaceI].transfer(cellIDs);
 114
 115             localPyrolysisFaceI++;
 116             nLayers_ = nCells;
 117         }
 118     }
 119
 120     boundaryFaceOppositeFace_.setSize(localPyrolysisFaceI);
 121
 122     surfaceScalarField& nMagSf = nMagSfPtr_();
 123
 124     forAll(intCoupledPatchIDs_, i)
 125     {
 126         const label patchI = intCoupledPatchIDs_[i];
 127         const polyPatch& ppCoupled = rbm[patchI];
 128         const vectorField& pNormals = ppCoupled.faceNormals();
 129         nMagSf.boundaryField()[patchI] =
 130             regionMesh().Sf().boundaryField()[patchI] & pNormals;
 131         forAll(pNormals, localFaceI)
 132         {
 133             const vector& n = pNormals[localFaceI];
 134             const labelList& faces = boundaryFaceFaces_[localFaceI];
 135             forAll (faces, faceI)
 136             {
 137                 const label faceID = faces[faceI];
 138                 nMagSf[faceID] = regionMesh().Sf()[faceID] & n;
 139             }
 140         }
 141     }
 142 }
 143
 144
 145 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 146
 147 bool Foam::regionModels::regionModel1D::read()
 148 {
 149     if (regionModel::read())
 150     {
 151         moveMesh_.readIfPresent("moveMesh", coeffs_);
 152
 153         return true;
 154     }
 155     else
 156     {
 157         return false;
 158     }
 159 }
 160
 161
 162 Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
 163 (
 164     const scalarList& deltaV,
 165     const scalar minDelta
 166 )
 167 {
 168     tmp<labelField> tcellMoveMap(new labelField(regionMesh().nCells(), 0));
 169     labelField& cellMoveMap = tcellMoveMap();
 170
 171     if (!moveMesh_)
 172     {
 173         return cellMoveMap;
 174     }
 175
 176     pointField oldPoints = regionMesh().points();
 177     pointField newPoints = oldPoints;
 178
 179     const polyBoundaryMesh& bm = regionMesh().boundaryMesh();
 180
 181     forAll(intCoupledPatchIDs_, localPatchI)
 182     {
 183         label patchI = intCoupledPatchIDs_[localPatchI];
 184         const polyPatch pp = bm[patchI];
 185         const vectorField& cf = regionMesh().Cf().boundaryField()[patchI];
 186
 187         forAll(pp, patchFaceI)
 188         {
 189             const labelList& faces = boundaryFaceFaces_[patchFaceI];
 190             const labelList& cells = boundaryFaceCells_[patchFaceI];
 191             const vector n = pp.faceNormals()[patchFaceI];
 192             const vector sf = pp.faceAreas()[patchFaceI];
 193
 194             List<point> oldCf(faces.size() + 1);
 195             oldCf[0] = cf[patchFaceI];
 196             forAll(faces, i)
 197             {
 198                 oldCf[i + 1] = regionMesh().faceCentres()[faces[i]];
 199             }
 200
 201             vector newDelta = vector::zero;
 202             point nbrCf = oldCf[0];
 203
 204             forAll(faces, i)
 205             {
 206                 const label faceI = faces[i];
 207                 const label cellI = cells[i];
 208
 209                 const face f = regionMesh().faces()[faceI];
 210
 211                 newDelta += (deltaV[cellI]/mag(sf))*n;
 212
 213                 vector localDelta = vector::zero;
 214                 forAll(f, pti)
 215                 {
 216                     const label pointI = f[pti];
 217
 218                     if
 219                     (
 220                         ((nbrCf - (oldPoints[pointI] + newDelta)) & n)
 221                       > minDelta
 222                     )
 223                     {
 224                         newPoints[pointI] = oldPoints[pointI] + newDelta;
 225                         localDelta = newDelta;
 226                         cellMoveMap[cellI] = 1;
 227                     }
 228                 }
 229                 nbrCf = oldCf[i + 1] + localDelta;
 230             }
 231
 232             // Modify boundary
 233             const label bFaceI = boundaryFaceOppositeFace_[patchFaceI];
 234             const face f = regionMesh().faces()[bFaceI];
 235             const label cellI = cells[cells.size() - 1];
 236             newDelta += (deltaV[cellI]/mag(sf))*n;
 237             forAll(f, pti)
 238             {
 239                 const label pointI = f[pti];
 240
 241                 if
 242                 (
 243                     ((nbrCf - (oldPoints[pointI] + newDelta)) & n)
 244                   > minDelta
 245                 )
 246                 {
 247                     newPoints[pointI] = oldPoints[pointI] + newDelta;
 248                     cellMoveMap[cellI] = 1;
 249                 }
 250             }
 251         }
 252     }
 253
 254     // Move points
 255     regionMesh().movePoints(newPoints);
 256
 257     return tcellMoveMap;
 258 }
 259
 260
 261 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 262
 263 Foam::regionModels::regionModel1D::regionModel1D(const fvMesh& mesh)
 264 :
 265     regionModel(mesh),
 266     boundaryFaceFaces_(),
 267     boundaryFaceCells_(),
 268     boundaryFaceOppositeFace_(),
 269     nLayers_(0),
 270     nMagSfPtr_(NULL),
 271     moveMesh_(false)
 272 {}
 273
 274
 275 Foam::regionModels::regionModel1D::regionModel1D
 276 (
 277     const fvMesh& mesh,
 278     const word& regionType,
 279     const word& modelName,
 280     bool readFields
 281 )
 282 :
 283     regionModel(mesh, regionType, modelName, false),
 284     boundaryFaceFaces_(regionMesh().nCells()),
 285     boundaryFaceCells_(regionMesh().nCells()),
 286     boundaryFaceOppositeFace_(regionMesh().nCells()),
 287     nLayers_(0),
 288     nMagSfPtr_(NULL),
 289     moveMesh_(true)
 290 {
 291     if (active_)
 292     {
 293         constructMeshObjects();
 294         initialise();
 295
 296         if (readFields)
 297         {
 298             read();
 299         }
 300     }
 301 }
 302
 303
 304 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 305
 306 Foam::regionModels::regionModel1D::~regionModel1D()
 307 {}
 308
 309
 310 // ************************************************************************* //