src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "movingConeTopoFvMesh.H"
  28 #include "Time.H"
  29 #include "mapPolyMesh.H"
  30 #include "layerAdditionRemoval.H"
  31 #include "addToRunTimeSelectionTable.H"
  32 #include "volMesh.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38     defineTypeNameAndDebug(movingConeTopoFvMesh, 0);
  39
  40     addToRunTimeSelectionTable
  41     (
  42         topoChangerFvMesh,
  43         movingConeTopoFvMesh,
  44         IOobject
  45     );
  46 }
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51 Foam::tmp<Foam::scalarField> Foam::movingConeTopoFvMesh::vertexMarkup
  52 (
  53     const pointField& p,
  54     const scalar& curLeft,
  55     const scalar& curRight
  56 ) const
  57 {
  58     Info<< "Updating vertex markup.  curLeft: "
  59         << curLeft << " curRight: " << curRight << endl;
  60
  61     tmp<scalarField> tvertexMarkup(new scalarField(p.size()));
  62     scalarField& vertexMarkup = tvertexMarkup();
  63
  64     forAll (p, pI)
  65     {
  66         if (p[pI].x() < curLeft - 1e-10)
  67         {
  68             vertexMarkup[pI] = -1;
  69         }
  70         else if (p[pI].x() > curRight + 1e-10)
  71         {
  72             vertexMarkup[pI] = 1;
  73         }
  74         else
  75         {
  76             vertexMarkup[pI] = 0;
  77         }
  78     }
  79
  80     return tvertexMarkup;
  81 }
  82
  83
  84 void Foam::movingConeTopoFvMesh::addZonesAndModifiers()
  85 {
  86     // Add zones and modifiers for motion action
  87
  88     if
  89     (
  90         pointZones().size() > 0
  91      || faceZones().size() > 0
  92      || cellZones().size() > 0
  93     )
  94     {
  95         Info<< "void movingConeTopoFvMesh::addZonesAndModifiers() : "
  96             << "Zones and modifiers already present.  Skipping."
  97             << endl;
  98
  99         if (topoChanger_.size() == 0)
 100         {
 101             FatalErrorIn
 102             (
 103                 "void movingConeTopoFvMesh::addZonesAndModifiers()"
 104             )   << "Mesh modifiers not read properly"
 105                 << abort(FatalError);
 106         }
 107
 108         return;
 109     }
 110
 111     Info<< "Time = " << time().timeName() << endl
 112         << "Adding zones and modifiers to the mesh" << endl;
 113
 114     const vectorField& fc = faceCentres();
 115     const vectorField& fa = faceAreas();
 116
 117     labelList zone1(fc.size());
 118     boolList flipZone1(fc.size(), false);
 119     label nZoneFaces1 = 0;
 120
 121     labelList zone2(fc.size());
 122     boolList flipZone2(fc.size(), false);
 123     label nZoneFaces2 = 0;
 124
 125     forAll (fc, faceI)
 126     {
 127         if
 128         (
 129             fc[faceI].x() > -0.003501
 130          && fc[faceI].x() < -0.003499
 131         )
 132         {
 133             if ((fa[faceI] & vector(1, 0, 0)) < 0)
 134             {
 135                 flipZone1[nZoneFaces1] = true;
 136             }
 137
 138             zone1[nZoneFaces1] = faceI;
 139             Info<< "face " << faceI << " for zone 1.  Flip: "
 140                 << flipZone1[nZoneFaces1] << endl;
 141             nZoneFaces1++;
 142         }
 143         else if
 144         (
 145             fc[faceI].x() > -0.00701
 146          && fc[faceI].x() < -0.00699
 147         )
 148         {
 149             zone2[nZoneFaces2] = faceI;
 150
 151             if ((fa[faceI] & vector(1, 0, 0)) > 0)
 152             {
 153                 flipZone2[nZoneFaces2] = true;
 154             }
 155
 156             Info << "face " << faceI << " for zone 2.  Flip: "
 157                 << flipZone2[nZoneFaces2] << endl;
 158             nZoneFaces2++;
 159         }
 160     }
 161
 162     zone1.setSize(nZoneFaces1);
 163     flipZone1.setSize(nZoneFaces1);
 164
 165     zone2.setSize(nZoneFaces2);
 166     flipZone2.setSize(nZoneFaces2);
 167
 168     Info << "zone: " << zone1 << endl;
 169     Info << "zone: " << zone2 << endl;
 170
 171     List<pointZone*> pz(0);
 172     List<faceZone*> fz(2);
 173     List<cellZone*> cz(0);
 174
 175     label nFz = 0;
 176
 177     fz[nFz] =
 178         new faceZone
 179         (
 180             "rightExtrusionFaces",
 181             zone1,
 182             flipZone1,
 183             nFz,
 184             faceZones()
 185         );
 186     nFz++;
 187
 188     fz[nFz] =
 189         new faceZone
 190         (
 191             "leftExtrusionFaces",
 192             zone2,
 193             flipZone2,
 194             nFz,
 195             faceZones()
 196         );
 197     nFz++;
 198
 199     fz.setSize(nFz);
 200
 201     Info << "Adding mesh zones." << endl;
 202     addZones(pz, fz, cz);
 203
 204
 205     // Add layer addition/removal interfaces
 206
 207     topoChanger_.setSize(2);
 208     label nMods = 0;
 209
 210     topoChanger_.set
 211     (
 212         0,
 213         new layerAdditionRemoval
 214         (
 215             "right",
 216             nMods,
 217             topoChanger_,
 218             "rightExtrusionFaces",
 219             readScalar
 220             (
 221                 motionDict_.subDict("right").lookup("minThickness")
 222             ),
 223             readScalar
 224             (
 225                 motionDict_.subDict("right").lookup("maxThickness")
 226             )
 227         )
 228     );
 229     nMods++;
 230
 231     topoChanger_.set
 232     (
 233         1,
 234         new layerAdditionRemoval
 235         (
 236             "left",
 237             nMods,
 238             topoChanger_,
 239             "leftExtrusionFaces",
 240             readScalar
 241             (
 242                 motionDict_.subDict("left").lookup("minThickness")
 243             ),
 244             readScalar
 245             (
 246                 motionDict_.subDict("left").lookup("maxThickness")
 247             )
 248         )
 249     );
 250     nMods++;
 251
 252     Info << "Adding " << nMods << " mesh modifiers" << endl;
 253
 254     // Write mesh and modifiers
 255     topoChanger_.write();
 256     write();
 257 }
 258
 259
 260 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 261
 262 // Construct from components
 263 Foam::movingConeTopoFvMesh::movingConeTopoFvMesh(const IOobject& io)
 264 :
 265     topoChangerFvMesh(io),
 266     motionDict_
 267     (
 268         IOdictionary
 269         (
 270             IOobject
 271             (
 272                 "dynamicMeshDict",
 273                 time().constant(),
 274                 *this,
 275                 IOobject::MUST_READ,
 276                 IOobject::NO_WRITE
 277             )
 278         ).subDict(typeName + "Coeffs")
 279     ),
 280     motionVelAmplitude_(motionDict_.lookup("motionVelAmplitude")),
 281     motionVelPeriod_(readScalar(motionDict_.lookup("motionVelPeriod"))),
 282     curMotionVel_
 283     (
 284         motionVelAmplitude_*
 285         Foam::sin(time().value()*M_PI/motionVelPeriod_)
 286     ),
 287     leftEdge_(readScalar(motionDict_.lookup("leftEdge"))),
 288     curLeft_(readScalar(motionDict_.lookup("leftObstacleEdge"))),
 289     curRight_(readScalar(motionDict_.lookup("rightObstacleEdge"))),
 290     motionMask_
 291     (
 292         vertexMarkup
 293         (
 294             points(),
 295             curLeft_,
 296             curRight_
 297         )
 298     )
 299 {
 300     addZonesAndModifiers();
 301
 302     curLeft_ = average
 303     (
 304         faceZones()
 305         [
 306             faceZones().findZoneID("leftExtrusionFaces")
 307         ]().localPoints()
 308     ).x();
 309
 310     curRight_ = average
 311     (
 312         faceZones()
 313         [
 314             faceZones().findZoneID("rightExtrusionFaces")
 315         ]().localPoints()
 316     ).x();
 317 }
 318
 319
 320 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 321
 322 Foam::movingConeTopoFvMesh::~movingConeTopoFvMesh()
 323 {}
 324
 325
 326 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 327
 328 bool Foam::movingConeTopoFvMesh::update()
 329 {
 330     autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 331
 332     // Calculate the new point positions depending on whether the
 333     // topological change has happened or not
 334     pointField newPoints;
 335
 336     vector curMotionVel_ =
 337         motionVelAmplitude_*
 338         Foam::sin(time().value()*M_PI/motionVelPeriod_);
 339
 340     bool meshChanged = topoChangeMap->morphing();
 341
 342     if (meshChanged)
 343     {
 344         Info << "Topology change. Calculating motion points" << endl;
 345
 346         if (topoChangeMap().hasMotionPoints())
 347         {
 348             motionMask_ =
 349                 vertexMarkup
 350                 (
 351                     topoChangeMap().preMotionPoints(),
 352                     curLeft_,
 353                     curRight_
 354                 );
 355         }
 356         else
 357         {
 358             motionMask_ =
 359                 vertexMarkup
 360                 (
 361                     points(),
 362                     curLeft_,
 363                     curRight_
 364                 );
 365         }
 366
 367         // Create new points by moving old points but using the
 368         // pre-motion points at the motion selector for the moving
 369         // region
 370         newPoints =
 371             points()
 372           + (
 373                 pos(0.5 - mag(motionMask_)) // cells above the body
 374 //               + pos(motionMask_ - 0.5)*      // cells in front of the body
 375 //                 (
 376 //                     points().component(vector::X)/curRight
 377 //                 )
 378 //               + pos(-motionMask_ - 0.5)*      // cells behind the body
 379 //                 (
 380 //                     (
 381 //                         points().component(vector::X)
 382 //                       - leftEdge
 383 //                     )/
 384 //                     (curLeft_ - leftEdge_)
 385 //                 )
 386             )*curMotionVel_*time().deltaT().value();
 387
 388         // Correct mesh motion for correct volume continuity
 389         movePoints(topoChangeMap().preMotionPoints());
 390         resetMotion();
 391         setV0();
 392     }
 393     else
 394     {
 395         Info << "No topology change" << endl;
 396         // Set the mesh motion
 397         newPoints =
 398             points()
 399           + (
 400                 pos(0.5 - mag(motionMask_)) // cells above the body
 401 //               + pos(motionMask_ - 0.5)*  // cells in front of the body
 402 //                 (
 403 //                     points().component(vector::X)/curRight
 404 //                 )
 405 //               + pos(-motionMask_ - 0.5)*  // cells behind the body
 406 //                 (
 407 //                     (
 408 //                         points().component(vector::X)
 409 //                       - leftEdge
 410 //                     )/
 411 //                     (curLeft_ - leftEdge_)
 412 //                 )
 413             )*curMotionVel_*time().deltaT().value();
 414     }
 415
 416     curLeft_ += curMotionVel_.x()*time().deltaT().value();
 417     curRight_ += curMotionVel_.x()*time().deltaT().value();
 418
 419     // The mesh now contains the cells with zero volume
 420
 421     Info << "Executing mesh motion" << endl;
 422     movePoints(newPoints);
 423
 424     //  The mesh now has got non-zero volume cells
 425
 426     return meshChanged;
 427 }
 428
 429
 430 // ************************************************************************* //