src/fvMotionSolver/fvMotionSolvers/displacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "displacementComponentLaplacianFvMotionSolver.H"
  28 #include "motionDiffusivity.H"
  29 #include "fvmLaplacian.H"
  30 #include "addToRunTimeSelectionTable.H"
  31 #include "mapPolyMesh.H"
  32 #include "volPointInterpolation.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38     defineTypeNameAndDebug(displacementComponentLaplacianFvMotionSolver, 0);
  39
  40     addToRunTimeSelectionTable
  41     (
  42         fvMotionSolver,
  43         displacementComponentLaplacianFvMotionSolver,
  44         dictionary
  45     );
  46 }
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51 Foam::direction Foam::displacementComponentLaplacianFvMotionSolver::cmpt
  52 (
  53     const word& cmptName
  54 ) const
  55 {
  56     if (cmptName == "x")
  57     {
  58         return vector::X;
  59     }
  60     else if (cmptName == "y")
  61     {
  62         return vector::Y;
  63     }
  64     else if (cmptName == "z")
  65     {
  66         return vector::Z;
  67     }
  68     else
  69     {
  70         FatalErrorIn
  71         (
  72             "displacementComponentLaplacianFvMotionSolver::"
  73             "displacementComponentLaplacianFvMotionSolver"
  74             "(const polyMesh& mesh, Istream& msData)"
  75         )   << "Given component name " << cmptName << " should be x, y or z"
  76             << exit(FatalError);
  77
  78         return 0;
  79     }
  80 }
  81
  82
  83 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  84
  85 Foam::displacementComponentLaplacianFvMotionSolver::
  86 displacementComponentLaplacianFvMotionSolver
  87 (
  88     const polyMesh& mesh,
  89     Istream& msData
  90 )
  91 :
  92     fvMotionSolver(mesh),
  93     cmptName_(msData),
  94     cmpt_(cmpt(cmptName_)),
  95     points0_
  96     (
  97         pointIOField
  98         (
  99             IOobject
 100             (
 101                 "points",
 102                 time().constant(),
 103                 polyMesh::meshSubDir,
 104                 mesh,
 105                 IOobject::MUST_READ,
 106                 IOobject::NO_WRITE,
 107                 false
 108             )
 109         ).component(cmpt_)
 110     ),
 111     pointDisplacement_
 112     (
 113         IOobject
 114         (
 115             "pointDisplacement" + cmptName_,
 116             fvMesh_.time().timeName(),
 117             fvMesh_,
 118             IOobject::MUST_READ,
 119             IOobject::AUTO_WRITE
 120         ),
 121         pointMesh::New(fvMesh_)
 122     ),
 123     cellDisplacement_
 124     (
 125         IOobject
 126         (
 127             "cellDisplacement" + cmptName_,
 128             mesh.time().timeName(),
 129             mesh,
 130             IOobject::READ_IF_PRESENT,
 131             IOobject::AUTO_WRITE
 132         ),
 133         fvMesh_,
 134         dimensionedScalar
 135         (
 136             "cellDisplacement",
 137             pointDisplacement_.dimensions(),
 138             0
 139         ),
 140         cellMotionBoundaryTypes<scalar>(pointDisplacement_.boundaryField())
 141     ),
 142     pointLocation_(NULL),
 143     diffusivityPtr_
 144     (
 145         motionDiffusivity::New(*this, lookup("diffusivity"))
 146     ),
 147     frozenPointsZone_
 148     (
 149         found("frozenPointsZone")
 150       ? fvMesh_.pointZones().findZoneID(lookup("frozenPointsZone"))
 151       : -1
 152     )
 153 {
 154     IOobject io
 155     (
 156         "pointLocation",
 157         fvMesh_.time().timeName(),
 158         fvMesh_,
 159         IOobject::MUST_READ,
 160         IOobject::AUTO_WRITE
 161     );
 162
 163     if (debug)
 164     {
 165         Info<< "displacementComponentLaplacianFvMotionSolver:" << nl
 166             << "    diffusivity       : " << diffusivityPtr_().type() << nl
 167             << "    frozenPoints zone : " << frozenPointsZone_ << endl;
 168     }
 169
 170     if (io.headerOk())
 171     {
 172         pointLocation_.reset
 173         (
 174             new pointVectorField
 175             (
 176                 io,
 177                 pointMesh::New(fvMesh_)
 178             )
 179         );
 180
 181         if (debug)
 182         {
 183             Info<< "displacementComponentLaplacianFvMotionSolver :"
 184                 << " Read pointVectorField "
 185                 << io.name() << " to be used for boundary conditions on points."
 186                 << nl
 187                 << "Boundary conditions:"
 188                 << pointLocation_().boundaryField().types() << endl;
 189         }
 190     }
 191 }
 192
 193
 194 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 195
 196 Foam::displacementComponentLaplacianFvMotionSolver::
 197 ~displacementComponentLaplacianFvMotionSolver()
 198 {}
 199
 200
 201 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 202
 203 Foam::tmp<Foam::pointField>
 204 Foam::displacementComponentLaplacianFvMotionSolver::curPoints() const
 205 {
 206     volPointInterpolation::New(fvMesh_).interpolate
 207     (
 208         cellDisplacement_,
 209         pointDisplacement_
 210     );
 211
 212     if (pointLocation_.valid())
 213     {
 214         if (debug)
 215         {
 216             Info<< "displacementComponentLaplacianFvMotionSolver : applying "
 217                 << " boundary conditions on " << pointLocation_().name()
 218                 << " to new point location."
 219                 << endl;
 220         }
 221
 222         // Apply pointLocation_ b.c. to mesh points.
 223
 224         pointLocation_().internalField() = fvMesh_.points();
 225
 226         pointLocation_().internalField().replace
 227         (
 228             cmpt_,
 229             points0_ + pointDisplacement_.internalField()
 230         );
 231
 232         pointLocation_().correctBoundaryConditions();
 233
 234         // Implement frozen points
 235         if (frozenPointsZone_ != -1)
 236         {
 237             const pointZone& pz = fvMesh_.pointZones()[frozenPointsZone_];
 238
 239             forAll(pz, i)
 240             {
 241                 label pointI = pz[i];
 242
 243                 pointLocation_()[pointI][cmpt_] = points0_[pointI];
 244             }
 245         }
 246
 247         twoDCorrectPoints(pointLocation_().internalField());
 248
 249         return tmp<pointField>(pointLocation_().internalField());
 250     }
 251     else
 252     {
 253         tmp<pointField> tcurPoints(new pointField(fvMesh_.points()));
 254
 255         tcurPoints().replace
 256         (
 257             cmpt_,
 258             points0_ + pointDisplacement_.internalField()
 259         );
 260
 261         // Implement frozen points
 262         if (frozenPointsZone_ != -1)
 263         {
 264             const pointZone& pz = fvMesh_.pointZones()[frozenPointsZone_];
 265
 266             forAll(pz, i)
 267             {
 268                 label pointI = pz[i];
 269
 270                 tcurPoints()[pointI][cmpt_] = points0_[pointI];
 271             }
 272         }
 273
 274         twoDCorrectPoints(tcurPoints());
 275
 276         return tcurPoints;
 277     }
 278 }
 279
 280
 281 void Foam::displacementComponentLaplacianFvMotionSolver::solve()
 282 {
 283     // The points have moved so before interpolation update
 284     // the fvMotionSolver accordingly
 285     movePoints(fvMesh_.points());
 286
 287     diffusivityPtr_->correct();
 288     pointDisplacement_.boundaryField().updateCoeffs();
 289
 290     Foam::solve
 291     (
 292         fvm::laplacian
 293         (
 294             diffusivityPtr_->operator()(),
 295             cellDisplacement_,
 296             "laplacian(diffusivity,cellDisplacement)"
 297         )
 298     );
 299 }
 300
 301
 302 void Foam::displacementComponentLaplacianFvMotionSolver::updateMesh
 303 (
 304     const mapPolyMesh& mpm
 305 )
 306 {
 307     fvMotionSolver::updateMesh(mpm);
 308
 309     // Map points0_. Bit special since we somehow have to come up with
 310     // a sensible points0 position for introduced points.
 311     // Find out scaling between points0 and current points
 312
 313     // Get the new points either from the map or the mesh
 314     const scalarField points
 315     (
 316         mpm.hasMotionPoints()
 317       ? mpm.preMotionPoints().component(cmpt_)
 318       : fvMesh_.points().component(cmpt_)
 319     );
 320
 321     // Get extents of points0 and points and determine scale
 322     const scalar scale =
 323         (gMax(points0_)-gMin(points0_))
 324        /(gMax(points)-gMin(points));
 325
 326     scalarField newPoints0(mpm.pointMap().size());
 327
 328     forAll(newPoints0, pointI)
 329     {
 330         label oldPointI = mpm.pointMap()[pointI];
 331
 332         if (oldPointI >= 0)
 333         {
 334             label masterPointI = mpm.reversePointMap()[oldPointI];
 335
 336             if (masterPointI == pointI)
 337             {
 338                 newPoints0[pointI] = points0_[oldPointI];
 339             }
 340             else
 341             {
 342                 // New point. Assume motion is scaling.
 343                 newPoints0[pointI] =
 344                     points0_[oldPointI]
 345                   + scale*(points[pointI]-points[masterPointI]);
 346             }
 347         }
 348         else
 349         {
 350             FatalErrorIn
 351             (
 352                 "displacementLaplacianFvMotionSolver::updateMesh"
 353                 "(const mapPolyMesh& mpm)"
 354             )   << "Cannot work out coordinates of introduced vertices."
 355                 << " New vertex " << pointI << " at coordinate "
 356                 << points[pointI] << exit(FatalError);
 357         }
 358     }
 359     points0_.transfer(newPoints0);
 360
 361     // Update diffusivity. Note two stage to make sure old one is de-registered
 362     // before creating/registering new one.
 363     diffusivityPtr_.reset(NULL);
 364     diffusivityPtr_ = motionDiffusivity::New(*this, lookup("diffusivity"));
 365 }
 366
 367
 368 // ************************************************************************* //