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) 2011 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 "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 bool Foam::regionModels::regionModel1D::read(const dictionary& dict)
 163 {
 164     if (regionModel::read(dict))
 165     {
 166         moveMesh_.readIfPresent("moveMesh", coeffs_);
 167
 168         return true;
 169     }
 170     else
 171     {
 172         return false;
 173     }
 174 }
 175
 176
 177 Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
 178 (
 179     const scalarList& deltaV,
 180     const scalar minDelta
 181 )
 182 {
 183     tmp<labelField> tcellMoveMap(new labelField(regionMesh().nCells(), 0));
 184     labelField& cellMoveMap = tcellMoveMap();
 185
 186     if (!moveMesh_)
 187     {
 188         return cellMoveMap;
 189     }
 190
 191     pointField oldPoints = regionMesh().points();
 192     pointField newPoints = oldPoints;
 193
 194     const polyBoundaryMesh& bm = regionMesh().boundaryMesh();
 195
 196     forAll(intCoupledPatchIDs_, localPatchI)
 197     {
 198         label patchI = intCoupledPatchIDs_[localPatchI];
 199         const polyPatch pp = bm[patchI];
 200         const vectorField& cf = regionMesh().Cf().boundaryField()[patchI];
 201
 202         forAll(pp, patchFaceI)
 203         {
 204             const labelList& faces = boundaryFaceFaces_[patchFaceI];
 205             const labelList& cells = boundaryFaceCells_[patchFaceI];
 206             const vector n = pp.faceNormals()[patchFaceI];
 207             const vector sf = pp.faceAreas()[patchFaceI];
 208
 209             List<point> oldCf(faces.size() + 1);
 210             oldCf[0] = cf[patchFaceI];
 211             forAll(faces, i)
 212             {
 213                 oldCf[i + 1] = regionMesh().faceCentres()[faces[i]];
 214             }
 215
 216             vector newDelta = vector::zero;
 217             point nbrCf = oldCf[0];
 218
 219             forAll(faces, i)
 220             {
 221                 const label faceI = faces[i];
 222                 const label cellI = cells[i];
 223
 224                 const face f = regionMesh().faces()[faceI];
 225
 226                 newDelta += (deltaV[cellI]/mag(sf))*n;
 227
 228                 vector localDelta = vector::zero;
 229                 forAll(f, pti)
 230                 {
 231                     const label pointI = f[pti];
 232
 233                     if
 234                     (
 235                         ((nbrCf - (oldPoints[pointI] + newDelta)) & n)
 236                       > minDelta
 237                     )
 238                     {
 239                         newPoints[pointI] = oldPoints[pointI] + newDelta;
 240                         localDelta = newDelta;
 241                         cellMoveMap[cellI] = 1;
 242                     }
 243                 }
 244                 nbrCf = oldCf[i + 1] + localDelta;
 245             }
 246
 247             // Modify boundary
 248             const label bFaceI = boundaryFaceOppositeFace_[patchFaceI];
 249             const face f = regionMesh().faces()[bFaceI];
 250             const label cellI = cells[cells.size() - 1];
 251             newDelta += (deltaV[cellI]/mag(sf))*n;
 252             forAll(f, pti)
 253             {
 254                 const label pointI = f[pti];
 255
 256                 if
 257                 (
 258                     ((nbrCf - (oldPoints[pointI] + newDelta)) & n)
 259                   > minDelta
 260                 )
 261                 {
 262                     newPoints[pointI] = oldPoints[pointI] + newDelta;
 263                     cellMoveMap[cellI] = 1;
 264                 }
 265             }
 266         }
 267     }
 268
 269     // Move points
 270     regionMesh().movePoints(newPoints);
 271
 272     return tcellMoveMap;
 273 }
 274
 275
 276 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 277
 278 Foam::regionModels::regionModel1D::regionModel1D(const fvMesh& mesh)
 279 :
 280     regionModel(mesh),
 281     boundaryFaceFaces_(),
 282     boundaryFaceCells_(),
 283     boundaryFaceOppositeFace_(),
 284     nLayers_(0),
 285     nMagSfPtr_(NULL),
 286     moveMesh_(false)
 287 {}
 288
 289
 290 Foam::regionModels::regionModel1D::regionModel1D
 291 (
 292     const fvMesh& mesh,
 293     const word& regionType,
 294     const word& modelName,
 295     bool readFields
 296 )
 297 :
 298     regionModel(mesh, regionType, modelName, false),
 299     boundaryFaceFaces_(regionMesh().nCells()),
 300     boundaryFaceCells_(regionMesh().nCells()),
 301     boundaryFaceOppositeFace_(regionMesh().nCells()),
 302     nLayers_(0),
 303     nMagSfPtr_(NULL),
 304     moveMesh_(true)
 305 {
 306     if (active_)
 307     {
 308         constructMeshObjects();
 309         initialise();
 310
 311         if (readFields)
 312         {
 313             read();
 314         }
 315     }
 316 }
 317
 318
 319 Foam::regionModels::regionModel1D::regionModel1D
 320 (
 321     const fvMesh& mesh,
 322     const word& regionType,
 323     const word& modelName,
 324     const dictionary& dict,
 325     bool readFields
 326 )
 327 :
 328     regionModel(mesh, regionType, modelName, dict, readFields),
 329     boundaryFaceFaces_(regionMesh().nCells()),
 330     boundaryFaceCells_(regionMesh().nCells()),
 331     boundaryFaceOppositeFace_(regionMesh().nCells()),
 332     nLayers_(0),
 333     nMagSfPtr_(NULL),
 334     moveMesh_(false)
 335 {
 336     if (active_)
 337     {
 338         constructMeshObjects();
 339         initialise();
 340
 341         if (readFields)
 342         {
 343             read(dict);
 344         }
 345     }
 346 }
 347
 348
 349 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 350
 351 Foam::regionModels::regionModel1D::~regionModel1D()
 352 {}
 353
 354
 355 // ************************************************************************* //