src/dynamicMesh/topoChangerFvMesh/mixerFvMesh/mixerFvMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "mixerFvMesh.H"
  27 #include "foamTime.H"
  28 #include "regionSplit.H"
  29 #include "slidingInterface.H"
  30 #include "mapPolyMesh.H"
  31 #include "volMesh.H"
  32 #include "addToRunTimeSelectionTable.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38     defineTypeNameAndDebug(mixerFvMesh, 0);
  39     addToRunTimeSelectionTable(topoChangerFvMesh, mixerFvMesh, IOobject);
  40 }
  41
  42
  43 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  44
  45 void Foam::mixerFvMesh::addZonesAndModifiers()
  46 {
  47     // Add zones and modifiers for motion action
  48     if (cellZones().size() > 0)
  49     {
  50         Info<< "void mixerFvMesh::addZonesAndModifiers() : "
  51             << "Zones and modifiers already present.  Skipping."
  52             << endl;
  53
  54         // Check definition of the modifier
  55         if
  56         (
  57             pointZones().size() > 0
  58          || faceZones().size() > 0
  59         )
  60         {
  61             if (topoChanger_.size() == 0)
  62             {
  63                 FatalErrorIn
  64                 (
  65                     "void mixerFvMesh::addZonesAndModifiers()"
  66                 )   << "Mesh modifiers not read properly.  "
  67                     << "pointZones = " <<  pointZones().size()
  68                     << " faceZones = " << faceZones().size()
  69                     << abort(FatalError);
  70             }
  71         }
  72
  73         return;
  74     }
  75
  76     // Find patches and check sizes
  77
  78     // Moving slider
  79     const word movingSliderName(dict_.subDict("slider").lookup("moving"));
  80     label movingSliderIndex = boundaryMesh().findPatchID(movingSliderName);
  81
  82     if (movingSliderIndex < 0)
  83     {
  84         FatalErrorIn("void mixerFvMesh::addZonesAndModifiers() const")
  85             << "Moving slider patch not found in boundary"
  86             << abort(FatalError);
  87     }
  88
  89     const word staticSliderName(dict_.subDict("slider").lookup("static"));
  90     label staticSliderIndex = boundaryMesh().findPatchID(staticSliderName);
  91
  92     if (staticSliderIndex < 0)
  93     {
  94         FatalErrorIn("void mixerFvMesh::addZonesAndModifiers() const")
  95             << "Static slider patch not found in boundary"
  96             << abort(FatalError);
  97
  98     }
  99
 100     const polyPatch& movingSlider = boundaryMesh()[movingSliderIndex];
 101     const polyPatch& staticSlider = boundaryMesh()[staticSliderIndex];
 102
 103     List<pointZone*> pz;
 104     List<faceZone*> fz;
 105
 106     bool addSlider = false;
 107
 108     if (movingSlider.size() > 0 && staticSlider.size() > 0)
 109     {
 110         addSlider = true;
 111
 112         pz.setSize(1);
 113         fz.setSize(3);
 114
 115         Info<< "Time = " << time().timeName() << endl
 116             << "Adding zones and modifiers to the mesh" << endl;
 117
 118         // Add zones
 119
 120         // Add an empty zone for cut points
 121
 122         pz[0] = new pointZone
 123         (
 124             "cutPointZone",
 125             labelList(0),
 126             0,
 127             pointZones()
 128         );
 129
 130
 131         // Do face zones for slider
 132
 133         // Moving slider
 134         labelList isf(movingSlider.size());
 135
 136         forAll (isf, i)
 137         {
 138             isf[i] = movingSlider.start() + i;
 139         }
 140
 141         fz[0] = new faceZone
 142         (
 143             movingSliderName + "Zone",
 144             isf,
 145             boolList(movingSlider.size(), false),
 146             0,
 147             faceZones()
 148         );
 149
 150         // Static slider
 151         labelList osf(staticSlider.size());
 152
 153         forAll (osf, i)
 154         {
 155             osf[i] = staticSlider.start() + i;
 156         }
 157
 158         fz[1] = new faceZone
 159         (
 160             staticSliderName + "Zone",
 161             osf,
 162             boolList(staticSlider.size(), false),
 163             1,
 164             faceZones()
 165         );
 166
 167         // Add empty zone for cut faces
 168         fz[2] = new faceZone
 169         (
 170             "cutFaceZone",
 171             labelList(0),
 172             boolList(0, false),
 173             2,
 174             faceZones()
 175         );
 176     }
 177
 178     // Add cell zone even if slider does not exist
 179
 180     List<cellZone*> cz(1);
 181
 182     // Mark every cell with its topological region
 183     regionSplit rs(*this);
 184
 185     // Get the region of the cell containing the origin.
 186     label originRegion = rs[findNearestCell(rotatingRegionMarker_)];
 187
 188     labelList movingCells(nCells());
 189     label nMovingCells = 0;
 190
 191     forAll(rs, cellI)
 192     {
 193         if (rs[cellI] == originRegion)
 194         {
 195             movingCells[nMovingCells] = cellI;
 196             nMovingCells++;
 197         }
 198     }
 199
 200     movingCells.setSize(nMovingCells);
 201     Info << "Number of cells in the moving region: " << nMovingCells << endl;
 202
 203     cz[0] = new cellZone
 204     (
 205         "movingCells",
 206         movingCells,
 207         0,
 208         cellZones()
 209     );
 210
 211     Info << "Adding point, face and cell zones" << endl;
 212     addZones(pz, fz, cz);
 213
 214     if (addSlider)
 215     {
 216         // Add a topology modifier
 217         Info << "Adding topology modifiers" << endl;
 218         topoChanger_.setSize(1);
 219
 220         topoChanger_.set
 221         (
 222             0,
 223             new slidingInterface
 224             (
 225                 "mixerSlider",
 226                 0,
 227                 topoChanger_,
 228                 staticSliderName + "Zone",
 229                 movingSliderName + "Zone",
 230                 "cutPointZone",
 231                 "cutFaceZone",
 232                 staticSliderName,
 233                 movingSliderName,
 234                 slidingInterface::INTEGRAL,   // Edge matching algorithm
 235                 attachDetach_,                // Attach-detach action
 236                 intersection::VISIBLE         // Projection algorithm
 237             )
 238         );
 239     }
 240
 241     // Write mesh and modifiers
 242     topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
 243     topoChanger_.write();
 244     write();
 245 }
 246
 247
 248 bool Foam::mixerFvMesh::attached() const
 249 {
 250     bool result = false;
 251
 252     if (topoChanger_.size() > 0)
 253     {
 254         result = refCast<const slidingInterface>(topoChanger_[0]).attached();
 255     }
 256
 257     reduce(result, orOp<bool>());
 258
 259     return result;
 260 }
 261
 262
 263 void Foam::mixerFvMesh::calcMovingMask() const
 264 {
 265     if (debug)
 266     {
 267         Info<< "void mixerFvMesh::calcMovingMask() const : "
 268             << "Calculating point and cell masks"
 269             << endl;
 270     }
 271
 272     if (movingPointsMaskPtr_)
 273     {
 274         FatalErrorIn("void mixerFvMesh::calcMovingMask() const")
 275             << "point mask already calculated"
 276             << abort(FatalError);
 277     }
 278
 279     // Set the point mask
 280     movingPointsMaskPtr_ = new scalarField(allPoints().size(), 0);
 281     scalarField& movingPointsMask = *movingPointsMaskPtr_;
 282
 283     const cellList& c = cells();
 284     const faceList& f = allFaces();
 285
 286     const labelList& cellAddr =
 287         cellZones()[cellZones().findZoneID("movingCells")];
 288
 289     forAll (cellAddr, cellI)
 290     {
 291         const cell& curCell = c[cellAddr[cellI]];
 292
 293         forAll (curCell, faceI)
 294         {
 295             // Mark all the points as moving
 296             const face& curFace = f[curCell[faceI]];
 297
 298             forAll (curFace, pointI)
 299             {
 300                 movingPointsMask[curFace[pointI]] = 1;
 301             }
 302         }
 303     }
 304
 305     if (topoChanger_.size() > 0)
 306     {
 307         // Topo changer present.  Use zones for motion
 308         const word movingSliderZoneName
 309         (
 310             word(dict_.subDict("slider").lookup("moving")) + "Zone"
 311         );
 312
 313         const labelList& movingSliderAddr =
 314             faceZones()[faceZones().findZoneID(movingSliderZoneName)];
 315
 316         forAll (movingSliderAddr, faceI)
 317         {
 318             const face& curFace = f[movingSliderAddr[faceI]];
 319
 320             forAll (curFace, pointI)
 321             {
 322                 movingPointsMask[curFace[pointI]] = 1;
 323             }
 324         }
 325
 326         const word staticSliderZoneName
 327         (
 328             word(dict_.subDict("slider").lookup("static")) + "Zone"
 329         );
 330
 331         const labelList& staticSliderAddr =
 332             faceZones()[faceZones().findZoneID(staticSliderZoneName)];
 333
 334         forAll (staticSliderAddr, faceI)
 335         {
 336             const face& curFace = f[staticSliderAddr[faceI]];
 337
 338             forAll (curFace, pointI)
 339             {
 340                 movingPointsMask[curFace[pointI]] = 0;
 341             }
 342         }
 343     }
 344 }
 345
 346
 347 // Return moving points mask.  Moving points marked with 1
 348 const Foam::scalarField& Foam::mixerFvMesh::movingPointsMask() const
 349 {
 350     if (!movingPointsMaskPtr_)
 351     {
 352         calcMovingMask();
 353     }
 354
 355     return *movingPointsMaskPtr_;
 356 }
 357
 358
 359 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 360
 361 // Construct from components
 362 Foam::mixerFvMesh::mixerFvMesh
 363 (
 364     const IOobject& io
 365 )
 366 :
 367     topoChangerFvMesh(io),
 368     dict_
 369     (
 370         IOdictionary
 371         (
 372             IOobject
 373             (
 374                 "dynamicMeshDict",
 375                 time().constant(),
 376                 *this,
 377                 IOobject::MUST_READ,
 378                 IOobject::NO_WRITE
 379             )
 380         ).subDict(typeName + "Coeffs")
 381     ),
 382     cs_
 383     (
 384         "coordinateSystem",
 385         dict_.subDict("coordinateSystem")
 386     ),
 387     rpm_(readScalar(dict_.lookup("rpm"))),
 388     rotatingRegionMarker_
 389     (
 390         dict_.lookupOrDefault<point>("rotatingRegionMarker", cs_.origin())
 391     ),
 392     attachDetach_(dict_.lookupOrDefault<bool>("attachDetach", false)),
 393     movingPointsMaskPtr_(NULL)
 394 {
 395     // Make sure the coordinate system does not operate in degrees
 396     // Bug fix, HJ, 3/Oct/2011
 397     if (!cs_.inDegrees())
 398     {
 399         WarningIn("mixerFvMesh::mixerFvMesh(const IOobject& io)")
 400             << "Mixer coordinate system is set to operate in radians.  "
 401             << "Changing to rad for correct calculation of angular velocity."
 402             << nl
 403             << "To remove this message please add entry" << nl << nl
 404             << "inDegrees true;" << nl << nl
 405             << "to the specification of the coordinate system"
 406             << endl;
 407
 408         cs_.inDegrees() = true;
 409     }
 410
 411     Info<< "Rotating region marker point: " << rotatingRegionMarker_
 412         << "  Attach-detach action = " << attachDetach_ << endl;
 413
 414     addZonesAndModifiers();
 415
 416     Info<< "Mixer mesh" << nl
 417         << "    origin       : " << cs().origin() << nl
 418         << "    axis         : " << cs().axis() << nl
 419         << "    rpm          : " << rpm_ << nl
 420         << "    marker       : " << rotatingRegionMarker_ << nl
 421         << "    attach-detach: " << attachDetach_
 422         << endl;
 423 }
 424
 425
 426 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 427
 428 Foam::mixerFvMesh::~mixerFvMesh()
 429 {
 430     deleteDemandDrivenData(movingPointsMaskPtr_);
 431 }
 432
 433
 434 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 435
 436 bool Foam::mixerFvMesh::update()
 437 {
 438     if (attached() && attachDetach_)
 439     {
 440         // Detach mesh
 441         Info << "Detaching rotors" << endl;
 442         autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 443
 444         deleteDemandDrivenData(movingPointsMaskPtr_);
 445     }
 446
 447     // Save old points
 448     pointField oldPointsNew = allPoints();
 449
 450     // Move points.  Rotational speed needs to be converted from rpm
 451     movePoints
 452     (
 453         cs_.globalPosition
 454         (
 455             cs_.localPosition(allPoints())
 456           + vector(0, rpm_*360.0*time().deltaT().value()/60.0, 0)
 457             *movingPointsMask()
 458         )
 459     );
 460
 461     autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 462
 463     bool localMorphing = topoChangeMap->morphing();
 464     bool globalMorphing = localMorphing;
 465
 466     reduce(globalMorphing, orOp<bool>());
 467
 468     if (globalMorphing)
 469     {
 470         Info<< "Attaching rotors" << endl;
 471
 472         deleteDemandDrivenData(movingPointsMaskPtr_);
 473
 474         // Move the sliding interface points to correct position
 475         if (localMorphing)
 476         {
 477             pointField mappedOldPointsNew(allPoints().size());
 478             mappedOldPointsNew.map(oldPointsNew, topoChangeMap->pointMap());
 479
 480             // Note: using setOldPoints instead of movePoints.
 481             // HJ, 23/Aug/2015
 482             setOldPoints(mappedOldPointsNew);
 483
 484             resetMotion();
 485             setV0();
 486
 487             // Move the sliding interface points to correct position
 488             movePoints(topoChangeMap->preMotionPoints());
 489         }
 490         else
 491         {
 492             pointField newPoints = allPoints();
 493
 494             // Note: using setOldPoints instead of movePoints.
 495             // HJ, 23/Aug/2015
 496             setOldPoints(oldPointsNew);
 497
 498             resetMotion();
 499             setV0();
 500
 501             // Move the sliding interface points to correct position
 502             movePoints(newPoints);
 503         }
 504     }
 505
 506     return topoChangeMap->morphing();
 507 }
 508
 509
 510 // ************************************************************************* //