src/dynamicMesh/topoChangerFvMesh/linearValveLayersFvMesh/linearValveLayersFvMesh.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 "linearValveLayersFvMesh.H"
  27 #include "foamTime.H"
  28 #include "slidingInterface.H"
  29 #include "layerAdditionRemoval.H"
  30 #include "pointField.H"
  31 #include "mapPolyMesh.H"
  32 #include "polyTopoChange.H"
  33 #include "pointZone.H"
  34 #include "volMesh.H"
  35 #include "addToRunTimeSelectionTable.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41     defineTypeNameAndDebug(linearValveLayersFvMesh, 0);
  42
  43     addToRunTimeSelectionTable
  44     (
  45         topoChangerFvMesh,
  46         linearValveLayersFvMesh,
  47         IOobject
  48     );
  49 }
  50
  51
  52 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  53
  54 void Foam::linearValveLayersFvMesh::addZonesAndModifiers()
  55 {
  56     // Inner slider
  57     const word innerSliderName(motionDict_.subDict("slider").lookup("inside"));
  58
  59     // Outer slider
  60     const word outerSliderName
  61     (
  62         motionDict_.subDict("slider").lookup("outside")
  63     );
  64
  65     bool initialised = false;
  66
  67     // Check if zones and modifiers for motion action are present
  68     label insideZoneID = faceZones().findZoneID(innerSliderName + "Zone");
  69     label outsideZoneID = faceZones().findZoneID(outerSliderName + "Zone");
  70
  71     if
  72     (
  73         insideZoneID > -1
  74      || outsideZoneID > -1
  75     )
  76     {
  77         // Zones found.  Check topo changer
  78
  79         if (topoChanger_.empty())
  80         {
  81             FatalErrorIn
  82             (
  83                 "void linearValveLayersFvMesh::addZonesAndModifiers()"
  84             )   << "Mesh modifiers not read properly"
  85                 << abort(FatalError);
  86         }
  87
  88         initialised = true;
  89     }
  90
  91     // Check if slider has been initialised on any of the processors
  92     reduce(initialised, orOp<bool>());
  93
  94     if (initialised)
  95     {
  96         InfoIn("void linearValveLayersFvMesh::addZonesAndModifiers()")
  97             << "Zones and modifiers already present.  Skipping."
  98             << endl;
  99
 100         return;
 101     }
 102
 103     // Add zones and modifiers for motion action
 104     Info<< "Time = " << time().timeName() << endl
 105         << "Adding zones and modifiers to the mesh" << endl;
 106
 107     // Add zones
 108     label nPz = 0;
 109     label nFz = 0;
 110     label nCz = 0;
 111     List<pointZone*> pz(pointZones().size() + 1);
 112     List<faceZone*> fz(faceZones().size() + 4);
 113     List<cellZone*> cz(cellZones().size());
 114
 115     // Add a topology modifier
 116     topoChanger_.setSize(2);
 117     label nTc = 0;
 118
 119     // Copy existing point zones
 120     forAll (pointZones(), zoneI)
 121     {
 122         pz[nPz] = pointZones()[zoneI].clone(pointZones()).ptr();
 123         nPz++;
 124     }
 125
 126     // Copy existing face zones
 127     forAll (faceZones(), zoneI)
 128     {
 129         fz[nFz] = faceZones()[zoneI].clone(faceZones()).ptr();
 130         nFz++;
 131     }
 132
 133     // Copy existing cell zones
 134     forAll (cellZones(), zoneI)
 135     {
 136         cz[nCz] = cellZones()[zoneI].clone(cellZones()).ptr();
 137         nCz++;
 138     }
 139
 140
 141
 142     // Do face zones for slider
 143
 144     // Inner slider
 145     const polyPatch& innerSlider =
 146         boundaryMesh()[boundaryMesh().findPatchID(innerSliderName)];
 147
 148     // Outer slider
 149     const polyPatch& outerSlider =
 150         boundaryMesh()[boundaryMesh().findPatchID(outerSliderName)];
 151
 152     if (!innerSlider.empty() && !outerSlider.empty())
 153     {
 154         Pout<< "Adding sliding interface between patches "
 155             << innerSliderName << " and " << outerSliderName << endl;
 156
 157
 158         // Add an empty zone for cut points
 159         pz[nPz] = new pointZone
 160         (
 161             "cutPointZone",
 162             labelList(0),
 163             nPz,
 164             pointZones()
 165         );
 166         nPz++;
 167
 168         labelList isf(innerSlider.size());
 169
 170         forAll (isf, i)
 171         {
 172             isf[i] = innerSlider.start() + i;
 173         }
 174
 175         fz[nFz] = new faceZone
 176         (
 177             innerSliderName + "Zone",
 178             isf,
 179             boolList(innerSlider.size(), false),
 180             nFz,
 181             faceZones()
 182         );
 183         nFz++;
 184
 185         labelList osf(outerSlider.size());
 186
 187         forAll (osf, i)
 188         {
 189             osf[i] = outerSlider.start() + i;
 190         }
 191
 192         fz[nFz] = new faceZone
 193         (
 194             outerSliderName + "Zone",
 195             osf,
 196             boolList(outerSlider.size(), false),
 197             nFz,
 198             faceZones()
 199         );
 200         nFz++;
 201
 202         // Add empty zone for cut faces
 203         fz[nFz] = new faceZone
 204         (
 205             "cutFaceZone",
 206             labelList(0),
 207             boolList(0, false),
 208             nFz,
 209             faceZones()
 210         );
 211         nFz++;
 212     }
 213
 214     // Add face zone for layer addition.  This is present on all processors
 215     const word layerPatchName
 216     (
 217         motionDict_.subDict("layer").lookup("patch")
 218     );
 219
 220     const polyPatch& layerPatch =
 221         boundaryMesh()[boundaryMesh().findPatchID(layerPatchName)];
 222
 223     labelList lpf(layerPatch.size());
 224
 225     forAll (lpf, i)
 226     {
 227         lpf[i] = layerPatch.start() + i;
 228     }
 229
 230     fz[nFz] = new faceZone
 231     (
 232         "valveLayerZone",
 233         lpf,
 234         boolList(layerPatch.size(), true),
 235         nFz,
 236         faceZones()
 237     );
 238     nFz++;
 239
 240     // Resize the number of live zones
 241     pz.setSize(nPz);
 242     fz.setSize(nFz);
 243     // Cell zones remain unchanged
 244
 245     Info << "Adding point and face zones" << endl;
 246     removeZones();
 247     addZones(pz, fz, cz);
 248
 249     if (!innerSlider.empty() && !outerSlider.empty())
 250     {
 251         // Set the topo changer for sliding interface
 252         topoChanger_.set
 253         (
 254             0,
 255             new slidingInterface
 256             (
 257                 "valveSlider",
 258                 0,
 259                 topoChanger_,
 260                 outerSliderName + "Zone",
 261                 innerSliderName + "Zone",
 262                 "cutPointZone",
 263                 "cutFaceZone",
 264                 outerSliderName,
 265                 innerSliderName,
 266                 slidingInterface::INTEGRAL,   // Edge matching algorithm
 267                 true,                         // Attach-detach action
 268                 intersection::VISIBLE         // Projection algorithm
 269             )
 270         );
 271
 272         // Record one added topo modifyer
 273         nTc = 1;
 274     }
 275
 276     // Add the layer addition-removal topo modifyer
 277     topoChanger_.set
 278     (
 279         nTc,
 280         new layerAdditionRemoval
 281         (
 282             "valveLayer",
 283             nTc,
 284             topoChanger_,
 285             "valveLayerZone",
 286             readScalar
 287             (
 288                 motionDict_.subDict("layer").lookup("minThickness")
 289             ),
 290             readScalar
 291             (
 292                 motionDict_.subDict("layer").lookup("maxThickness")
 293             )
 294         )
 295     );
 296     nTc++;
 297
 298     // Reset the size of mesh modifiers
 299     topoChanger_.setSize(nTc);
 300     Pout<< nTc << " topology modifiers on processor" << endl;
 301     // Write mesh and modifiers
 302     topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
 303     topoChanger_.write();
 304     write();
 305
 306     // Update the mesh for changes in zones.  This needs to
 307     // happen after write, because mesh instance will be changed
 308     // HJ and OP, 20/Nov/2013
 309     syncUpdateMesh();
 310 }
 311
 312
 313 void Foam::linearValveLayersFvMesh::makeLayersLive()
 314 {
 315     const polyTopoChanger& topoChanges = topoChanger_;
 316
 317     // Enable layering
 318     forAll (topoChanges, modI)
 319     {
 320         if (isA<layerAdditionRemoval>(topoChanges[modI]))
 321         {
 322             topoChanges[modI].enable();
 323         }
 324         else if (isA<slidingInterface>(topoChanges[modI]))
 325         {
 326             topoChanges[modI].disable();
 327         }
 328         else
 329         {
 330             FatalErrorIn("void linearValveLayersFvMesh::makeLayersLive()")
 331                 << "Don't know what to do with mesh modifier "
 332                 << modI << " of type " << topoChanges[modI].type()
 333                 << abort(FatalError);
 334         }
 335     }
 336 }
 337
 338
 339 void Foam::linearValveLayersFvMesh::makeSlidersLive()
 340 {
 341     const polyTopoChanger& topoChanges = topoChanger_;
 342
 343     // Enable sliding interface
 344     forAll (topoChanges, modI)
 345     {
 346         if (isA<layerAdditionRemoval>(topoChanges[modI]))
 347         {
 348             topoChanges[modI].disable();
 349         }
 350         else if (isA<slidingInterface>(topoChanges[modI]))
 351         {
 352             topoChanges[modI].enable();
 353         }
 354         else
 355         {
 356             FatalErrorIn("void linearValveLayersFvMesh::makeLayersLive()")
 357                 << "Don't know what to do with mesh modifier "
 358                 << modI << " of type " << topoChanges[modI].type()
 359                 << abort(FatalError);
 360         }
 361     }
 362 }
 363
 364
 365 bool Foam::linearValveLayersFvMesh::attached() const
 366 {
 367     const polyTopoChanger& topoChanges = topoChanger_;
 368
 369     bool result = false;
 370
 371     forAll (topoChanges, modI)
 372     {
 373         if (isA<slidingInterface>(topoChanges[modI]))
 374         {
 375             result =
 376                 result
 377              || refCast<const slidingInterface>(topoChanges[modI]).attached();
 378         }
 379     }
 380
 381     // Check that all sliders are in sync (debug only)
 382     forAll (topoChanges, modI)
 383     {
 384         if (isA<slidingInterface>(topoChanges[modI]))
 385         {
 386             if
 387             (
 388                 result
 389              != refCast<const slidingInterface>(topoChanges[modI]).attached()
 390             )
 391             {
 392                 FatalErrorIn("bool linearValveLayersFvMesh::attached() const")
 393                     << "Slider " << modI << " named "
 394                     << topoChanges[modI].name()
 395                     << " out of sync: Should be" << result
 396                     << abort(FatalError);
 397             }
 398         }
 399     }
 400
 401     // Sync across processors
 402     reduce(result, orOp<bool>());
 403
 404     return result;
 405 }
 406
 407
 408 Foam::tmp<Foam::pointField>
 409 Foam::linearValveLayersFvMesh::newLayerPoints() const
 410 {
 411     tmp<pointField> tnewLayerPoints
 412     (
 413         new pointField(allPoints())
 414     );
 415
 416     pointField& np = tnewLayerPoints();
 417
 418     const word layerPatchName
 419     (
 420         motionDict_.subDict("layer").lookup("patch")
 421     );
 422
 423     const polyPatch& layerPatch =
 424         boundaryMesh()[boundaryMesh().findPatchID(layerPatchName)];
 425
 426     const labelList& patchPoints = layerPatch.meshPoints();
 427
 428     const vector vel
 429     (
 430         motionDict_.lookup("pistonVelocity")
 431     );
 432
 433     forAll (patchPoints, ppI)
 434     {
 435         np[patchPoints[ppI]] += vel*time().deltaT().value();
 436     }
 437
 438     return tnewLayerPoints;
 439 }
 440
 441
 442
 443 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 444
 445 // Construct from components
 446 Foam::linearValveLayersFvMesh::linearValveLayersFvMesh(const IOobject& io)
 447 :
 448     topoChangerFvMesh(io),
 449     motionDict_
 450     (
 451         IOdictionary
 452         (
 453             IOobject
 454             (
 455                 "dynamicMeshDict",
 456                 time().constant(),
 457                 *this,
 458                 IOobject::MUST_READ,
 459                 IOobject::NO_WRITE
 460             )
 461         ).subDict(typeName + "Coeffs")
 462     )
 463 {
 464     addZonesAndModifiers();
 465 }
 466
 467
 468 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 469
 470 Foam::linearValveLayersFvMesh::~linearValveLayersFvMesh()
 471 {}
 472
 473
 474 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 475
 476 bool Foam::linearValveLayersFvMesh::update()
 477 {
 478     // Detaching the interface
 479     if (attached())
 480     {
 481         Info<< "Decoupling sliding interfaces" << endl;
 482         makeSlidersLive();
 483
 484         // Changing topology by hand
 485         autoPtr<mapPolyMesh> topoChangeMap1 = topoChanger_.changeMesh();
 486
 487         bool localMorphing1 = topoChangeMap1->morphing();
 488
 489         // Note: Since we are detaching, global morphing is always true
 490         // HJ, 7/Mar/2011
 491
 492         if (localMorphing1)
 493         {
 494             Info << "Topology change; executing pre-motion after "
 495                 << "sliding detach" << endl;
 496             movePoints(topoChangeMap1->preMotionPoints());
 497         }
 498         else
 499         {
 500             pointField newPoints = allPoints();
 501
 502             // Dummy motion
 503             movePoints(newPoints);
 504         }
 505
 506         Info<< "sliding interfaces successfully decoupled!!!" << endl;
 507     }
 508     else
 509     {
 510         Info<< "Sliding interfaces decoupled" << endl;
 511     }
 512
 513     // Perform layer action and mesh motion
 514     makeLayersLive();
 515
 516     // Changing topology by hand
 517     autoPtr<mapPolyMesh> topoChangeMap2 = topoChanger_.changeMesh();
 518
 519     bool localMorphing2 = topoChangeMap2->morphing();
 520     bool globalMorphing2 = localMorphing2;
 521
 522     reduce(globalMorphing2, orOp<bool>());
 523
 524     // Work array for new points position.
 525     pointField newPoints = allPoints();
 526
 527     if (globalMorphing2)
 528     {
 529         Info<< "Topology change; executing pre-motion after "
 530             << "dynamic layering" << endl;
 531
 532         if (localMorphing2)
 533         {
 534             Info << "Topology change; executing pre-motion" << endl;
 535             // Note: using setOldPoints instead of movePoints.
 536             // HJ, 23/Aug/2015
 537             setOldPoints(topoChangeMap2->preMotionPoints());
 538             newPoints = topoChangeMap2->preMotionPoints();
 539         }
 540         else
 541         {
 542             // Note: using setOldPoints instead of movePoints.
 543             // HJ, 23/Aug/2015
 544             setOldPoints(newPoints);
 545         }
 546
 547         setV0();
 548         resetMotion();
 549     }
 550
 551     // Move points to change layer thickness
 552     movePoints(newLayerPoints());
 553
 554     // Changing topology by hand
 555     {
 556         // Grab old points to correct the motion
 557         pointField oldPointsNew = oldAllPoints();
 558
 559         // Attach the interface
 560         Info << "Coupling sliding interfaces" << endl;
 561         makeSlidersLive();
 562
 563         autoPtr<mapPolyMesh> topoChangeMap3 = topoChanger_.changeMesh();
 564
 565         bool localMorphing3 = topoChangeMap3->morphing();
 566         bool globalMorphing3 = localMorphing3;
 567
 568         reduce(globalMorphing3, orOp<bool>());
 569
 570         if (globalMorphing3)
 571         {
 572             Info<< "Topology change; executing pre-motion after "
 573                 << "sliding attach" << endl;
 574
 575             // Grab points
 576             newPoints = allPoints();
 577
 578             if (localMorphing3)
 579             {
 580                 // If there is layering, pick up correct points
 581                 if (topoChangeMap3->hasMotionPoints())
 582                 {
 583                     newPoints = topoChangeMap3->preMotionPoints();
 584                 }
 585
 586                 pointField mappedOldPointsNew(newPoints.size());
 587
 588                 mappedOldPointsNew.map
 589                 (
 590                     oldPointsNew,
 591                     topoChangeMap3->pointMap()
 592                 );
 593
 594                 // Solve the correct mesh motion to make sure motion fluxes
 595                 // are solved for and not mapped
 596                 // Note: using setOldPoints instead of movePoints.
 597                 // HJ, 23/Aug/2015
 598                 setOldPoints(mappedOldPointsNew);
 599
 600                 resetMotion();
 601                 setV0();
 602
 603                 // Set new point motion
 604                 movePoints(newPoints);
 605             }
 606             else
 607             {
 608                 // No local topological change.  Execute double motion for
 609                 // sync with topological changes
 610                 // Note: using setOldPoints instead of movePoints.
 611                 // HJ, 23/Aug/2015
 612                 setOldPoints(oldPointsNew);
 613
 614                 resetMotion();
 615                 setV0();
 616
 617                 // Set new point motion
 618                 movePoints(newPoints);
 619             }
 620         }
 621     }
 622
 623     return true;
 624 }
 625
 626
 627 // ************************************************************************* //