src/engine/engineTopoChangerMesh/engineValveSliding/engineValveSlidingMove.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 Class
  25     verticalValvesGambit
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "engineValveSliding.H"
  30 #include "slidingInterface.H"
  31 #include "layerAdditionRemoval.H"
  32 #include "attachDetach.H"
  33 #include "surfaceFields.H"
  34 #include "regionSplit.H"
  35 #include "componentMixedTetPolyPatchVectorField.H"
  36 #include "mapPolyMesh.H"
  37
  38 #include "tetPolyMesh.H"
  39 #include "tetPointFields.H"
  40 #include "elementFields.H"
  41 #include "fixedValueTetPolyPatchFields.H"
  42 #include "slipTetPolyPatchFields.H"
  43
  44 #include "tetFem.H"
  45
  46 #include "symmetryFvPatch.H"
  47 #include "wedgeFvPatch.H"
  48 #include "emptyFvPatch.H"
  49 #include "zeroGradientTetPolyPatchFields.H"
  50 #include "tetMotionSolver.H"
  51
  52 #include "fixedValueTetPolyPatchFields.H"
  53 #include "mixedTetPolyPatchFields.H"
  54 #include "slipTetPolyPatchFields.H"
  55 #include "zeroGradientTetPolyPatchFields.H"
  56
  57 #include "zeroGradientFvPatchFields.H"
  58 #include "fvCFD.H"
  59
  60 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  61
  62 void Foam::engineValveSliding::makeLayersLive()
  63 {
  64     // Enable layering
  65     forAll (topoChanger_, modI)
  66     {
  67         if (isA<layerAdditionRemoval>(topoChanger_[modI]))
  68         {
  69             topoChanger_[modI].enable();
  70         }
  71         else if (isA<slidingInterface>(topoChanger_[modI]))
  72         {
  73             topoChanger_[modI].disable();
  74         }
  75         else if (isA<attachDetach>(topoChanger_[modI]))
  76         {
  77             topoChanger_[modI].enable();
  78         }
  79         else
  80         {
  81             FatalErrorIn("void Foam::engineTopoFvMesh::makeLayersLive()")
  82                 << "Don't know what to do with mesh modifier "
  83                 << modI << " of type " << topoChanger_[modI].type()
  84                 << abort(FatalError);
  85         }
  86     }
  87 }
  88
  89 void Foam::engineValveSliding::makeSlidersLive()
  90 {
  91     // Enable sliding interface
  92     forAll (topoChanger_, modI)
  93     {
  94         if (isA<layerAdditionRemoval>(topoChanger_[modI]))
  95         {
  96             topoChanger_[modI].disable();
  97         }
  98         else if (isA<slidingInterface>(topoChanger_[modI]))
  99         {
 100             topoChanger_[modI].enable();
 101         }
 102         else if (isA<attachDetach>(topoChanger_[modI]))
 103         {
 104             topoChanger_[modI].enable();
 105         }
 106         else
 107         {
 108             FatalErrorIn("void Foam::engineTopoFvMesh::makeLayersLive()")
 109                 << "Don't know what to do with mesh modifier "
 110                 << modI << " of type " << topoChanger_[modI].type()
 111                 << abort(FatalError);
 112         }
 113     }
 114
 115 }
 116
 117 bool Foam::engineValveSliding::attached() const
 118 {
 119     const polyTopoChanger& morphs = topoChanger_;
 120
 121     bool result = false;
 122
 123     forAll (morphs, modI)
 124     {
 125         if (typeid(morphs[modI]) == typeid(slidingInterface))
 126         {
 127             result =
 128                 result
 129              || refCast<const slidingInterface>(morphs[modI]).attached();
 130         }
 131     }
 132
 133     // Check thal all sliders are in sync (debug only)
 134     forAll (morphs, modI)
 135     {
 136         if (typeid(morphs[modI]) == typeid(slidingInterface))
 137         {
 138             if
 139             (
 140                 result
 141              != refCast<const slidingInterface>(morphs[modI]).attached()
 142             )
 143             {
 144                 FatalErrorIn("bool movingSquaresTM::attached() const")
 145                     << "Slider " << modI << " named " << morphs[modI].name()
 146                     << " out of sync: Should be" << result
 147                     << abort(FatalError);
 148             }
 149         }
 150     }
 151
 152     return result;
 153 }
 154
 155
 156 void Foam::engineValveSliding::valveDetach()
 157 {
 158     // Enable sliding interface
 159     forAll (topoChanger_, modI)
 160     {
 161         if (isA<attachDetach>(topoChanger_[modI]))
 162         {
 163             const attachDetach& ad =
 164                 refCast<const attachDetach>(topoChanger_[modI]);
 165
 166             const word masterName = ad.masterPatchID().name();
 167
 168             // Find the valve with that name
 169             label valveIndex = -1;
 170
 171             forAll (valves_, valveI)
 172             {
 173                 if
 174                 (
 175                     valves_[valveI].detachInCylinderPatchID().name()
 176                  == masterName
 177                 )
 178                 {
 179                     valveIndex = valveI;
 180                     break;
 181                 }
 182             }
 183
 184             if (valveIndex < 0)
 185             {
 186                 FatalErrorIn
 187                 (
 188                     "void Foam::engineTopoFvMesh::prepareValveDetach()"
 189                 )   << "Cannot match patch for attach/detach " << modI
 190                     << abort(FatalError);
 191             }
 192
 193             if (debug)
 194             {
 195                 Info<< " valveI: " << valveIndex << " attached: "
 196                     << ad.attached()
 197                     << " valve open: " << valves_[valveIndex].isOpen()
 198                     << endl;
 199             }
 200
 201             ad.setDetach();
 202         }
 203     }
 204 }
 205
 206 void Foam::engineValveSliding::valveAttach()
 207 {
 208     // Enable sliding interface
 209     forAll (topoChanger_, modI)
 210     {
 211         if (isA<attachDetach>(topoChanger_[modI]))
 212         {
 213             const attachDetach& ad =
 214                 refCast<const attachDetach>(topoChanger_[modI]);
 215
 216             const word masterName = ad.masterPatchID().name();
 217
 218             // Find the valve with that name
 219             label valveIndex = -1;
 220
 221             forAll (valves_, valveI)
 222             {
 223                 if
 224                 (
 225                     valves_[valveI].detachInCylinderPatchID().name()
 226                  == masterName
 227                 )
 228                 {
 229                     valveIndex = valveI;
 230                     break;
 231                 }
 232             }
 233
 234             if (valveIndex < 0)
 235             {
 236                 FatalErrorIn
 237                 (
 238                     "void Foam::engineTopoFvMesh::prepareValveDetach()"
 239                 )   << "Cannot match patch for attach/detach " << modI
 240                     << abort(FatalError);
 241             }
 242
 243             if (debug)
 244             {
 245                 Info<< " valveI: " << valveIndex << " attached: "
 246                     << ad.attached()
 247                     << " valve open: " << valves_[valveIndex].isOpen()
 248                     << endl;
 249             }
 250
 251             ad.setAttach();
 252         }
 253     }
 254 }
 255
 256
 257 void Foam::engineValveSliding::prepareValveDetach()
 258 {
 259     // Enable sliding interface
 260     forAll (topoChanger_, modI)
 261     {
 262         if (isA<attachDetach>(topoChanger_[modI]))
 263         {
 264             const attachDetach& ad =
 265                 refCast<const attachDetach>(topoChanger_[modI]);
 266
 267             const word masterName = ad.masterPatchID().name();
 268
 269             // Find the valve with that name
 270             label valveIndex = -1;
 271
 272             forAll (valves_, valveI)
 273             {
 274                 if
 275                 (
 276                     valves_[valveI].detachInCylinderPatchID().name()
 277                  == masterName
 278                 )
 279                 {
 280                     valveIndex = valveI;
 281                     break;
 282                 }
 283             }
 284
 285             if (valveIndex < 0)
 286             {
 287                 FatalErrorIn
 288                 (
 289                     "void Foam::engineTopoFvMesh::prepareValveDetach()"
 290                 )   << "Cannot match patch for attach/detach " << modI
 291                     << abort(FatalError);
 292             }
 293
 294             if (debug)
 295             {
 296                 Info<< " valveI: " << valveIndex << " attached: "
 297                     << ad.attached()
 298                     << " valve open: " << valves_[valveIndex].isOpen()
 299                     << endl;
 300             }
 301
 302             if (valves_[valveIndex].isOpen())
 303             {
 304                 ad.setAttach();
 305             }
 306             else
 307             {
 308                 ad.setDetach();
 309             }
 310         }
 311     }
 312 }
 313
 314
 315 bool Foam::engineValveSliding::update()
 316 {
 317
 318     tetMotionSolver& mSolver =
 319         refCast<tetMotionSolver>(msPtr_());
 320
 321     // Detaching the interfacethobois2DSlidingDeform
 322     if (attached())
 323     {
 324         Info << "Decoupling sliding interfaces" << endl;
 325         makeSlidersLive();
 326         valveDetach();
 327         autoPtr<mapPolyMesh> topoChangeMap1 = topoChanger_.changeMesh();
 328
 329         Info << "sliding interfaces successfully decoupled!!!" << endl;
 330         if (topoChangeMap1->morphing())
 331         {
 332             mSolver.updateMesh(topoChangeMap1());
 333         }
 334     }
 335     else
 336     {
 337         Info << "Sliding interfaces decoupled" << endl;
 338         valveDetach();
 339     }
 340
 341     Info << "Executing layer action" << endl;
 342
 343     scalar deltaZ = engTime().pistonDisplacement().value();
 344
 345     // Piston Layering
 346     makeLayersLive();
 347
 348     {
 349         forAll(valves_, valveI)
 350         {
 351             scalar valveDisplacement = valves_[valveI].curVelocity()*
 352                 valves_[valveI].cs().axis().z()*engTime().deltaT().value();
 353
 354             Info << "valveDisplacement = " << valveDisplacement << nl
 355                 << "valve velocity =" << valves_[valveI].curVelocity() << nl
 356                 << "valve lift =" << valves_[valveI].curLift() << nl
 357                 << "valve deformation lift ="
 358                 << valves_[valveI].deformationLift() << endl;
 359         }
 360
 361     }
 362
 363     forAll(valves_,valveI)
 364     {
 365         if(valves_[valveI].curLift() >= valves_[valveI].deformationLift())
 366         {
 367             if(valves_[valveI].poppetPatchID().active())
 368             {
 369                 // Find valve layer mesh modifier
 370                 const label valveLayerID =
 371                     topoChanger_.findModifierID
 372                     (
 373                         "valvePoppetLayer"
 374                       + Foam::name(valveI + 1)
 375                     );
 376
 377                 if (valveLayerID < 0)
 378                 {
 379                         FatalErrorIn("void engineFvMesh::moveAndMorph()")
 380                            << "valve modifier not found."
 381                            << abort(FatalError);
 382                 }
 383
 384                 topoChanger_[valveLayerID].enable();
 385             }
 386         }
 387         else
 388         {
 389             if(valves_[valveI].poppetPatchID().active())
 390             {
 391                 // Find valve layer mesh modifier
 392                 const label valveLayerID =
 393                     topoChanger_.findModifierID
 394                     (
 395                         "valvePoppetLayer"
 396                       + Foam::name(valveI + 1)
 397                     );
 398
 399                 if (valveLayerID < 0)
 400                 {
 401                         FatalErrorIn("void engineFvMesh::moveAndMorph()")
 402                            << "valve modifier not found."
 403                            << abort(FatalError);
 404                 }
 405
 406                 topoChanger_[valveLayerID].disable();
 407             }
 408         }
 409     }
 410
 411     pointField newPoints = points();
 412
 413     autoPtr<mapPolyMesh> topoChangeMap2 = topoChanger_.changeMesh();
 414
 415     // Changing topology by hand
 416     if (topoChangeMap2->morphing())
 417     {
 418         mSolver.updateMesh(topoChangeMap2());
 419
 420         if (topoChangeMap2->hasMotionPoints())
 421         {
 422             movePoints(topoChangeMap2->preMotionPoints());
 423             newPoints = topoChangeMap2->preMotionPoints();
 424         }
 425         setV0();
 426         resetMotion();
 427     }
 428
 429     // Work array for new points position.
 430
 431 // NEW !
 432
 433     // Reset the position of layered interfaces
 434
 435     {
 436
 437
 438 #        include "setMotionBoundaryConditionEngineValveSliding.H"
 439
 440         DynamicList<label> constrainedPoints(mSolver.curPoints()().size()/100);
 441         DynamicList<vector> constrainedVelocity
 442         (
 443             mSolver.curPoints()().size()/100
 444         );
 445
 446 #       include "setEngineValveSlidingConstraint.H"
 447
 448         forAll (constrainedPoints, i)
 449         {
 450             mSolver.setConstraint
 451             (
 452                 constrainedPoints[i],
 453                 constrainedVelocity[i]
 454             );
 455         }
 456
 457 //        Info << mSolver.motionU() << endl;
 458
 459         Info << "pre solve" << endl;
 460
 461         mSolver.solve();
 462
 463         newPoints = mSolver.curPoints();
 464         movePoints(newPoints);
 465
 466         Info << "max mesh phi, 1 = " << max(phi()) << endl;
 467         Info << "min mesh phi, 1 = " << min(phi()) << endl;
 468
 469         setVirtualPositions();
 470         mSolver.clearConstraints();
 471
 472
 473
 474 //      layering
 475
 476 #       include "moveValvePointsEngineValveSliding.H"
 477         movePoints(newPoints);
 478 //        newPoints = points();
 479         newPoints = allPoints();
 480         setVirtualPositions();
 481
 482         Info << "max mesh phi, 2 = " << max(phi()) << endl;
 483         Info << "min mesh phi, 2 = " << min(phi()) << endl;
 484
 485     }
 486
 487     Info << ", deckHeight = " << deckHeight()
 488         << ", pistonPosition = " << pistonPosition()    << endl;
 489
 490     pistonPosition() += deltaZ;
 491
 492     Info << "max V()-V0() before = " << max(mag(V() - V0())) << endl;
 493     Info << "max V() before = " << max(mag(V())) << endl;
 494
 495     Info<< "BEFORE mag(points - oldPoints)/deltaT = "
 496         << max(mag(points() - oldPoints()))/engTime().deltaT().value() << endl;
 497
 498 /*
 499     forAll(points(), pp)
 500     {
 501         Info<< "OLD " <<  "point " << pp << " " << points()[pp]
 502             << " old point = " << oldPoints()[pp] << mag(points()[pp] -
 503         oldPoints()[pp])/engTime().deltaT().value() << endl;
 504     }
 505
 506 */
 507
 508 //*/ //Tommaso, 23/5/2008
 509
 510     {
 511         // Grab old points to correct the motion
 512         pointField oldPointsNew = oldAllPoints();
 513 //        pointField oldPointsNew = oldPoints();
 514
 515         // Attach the interface
 516         Info << "Coupling sliding interfaces" << endl;
 517         makeSlidersLive();
 518
 519         // Changing topology by hand
 520         autoPtr<mapPolyMesh> topoChangeMap4 = topoChanger_.changeMesh();
 521
 522         // Changing topology by hand
 523         if(topoChangeMap4->morphing())
 524         {
 525             mSolver.updateMesh(topoChangeMap4());
 526
 527             if (topoChangeMap4->hasMotionPoints())
 528             {
 529 //                 Info<< "Topology change; executing pre-motion "
 530 //                     << "after sliding attach" << endl;
 531
 532 //                 Info<< "topoChangeMap4->preMotionPoints().size() = "
 533 //                     << topoChangeMap4->preMotionPoints().size() << endl;
 534 //                 Info << "allPoints.size() = " << allPoints().size() << endl;
 535 //                 Info << "points.size() = " << points().size() << endl;
 536
 537                 movePoints(topoChangeMap4->preMotionPoints());
 538 //            newPoints = topoChangeMap4->preMotionPoints();
 539 //            newPoints = allPoints();
 540                 newPoints = allPoints();
 541             }
 542         }
 543
 544         surfaceScalarField correctedMeshPhi = phi();
 545
 546
 547         Info << "max V()-V0() after = " << max(mag(V() - V0())) << endl;
 548         Info << "max V() before = " << max(mag(V())) << endl;
 549
 550         {
 551             Info << "AFTER mag(points - oldPoints)/deltaT = "
 552                 << max(mag(points() - oldPoints()))/engTime().deltaT().value()
 553                 << endl;
 554             Info << "max V()-V0() after movePoints?!? = "
 555                 << max(mag(V() - V0())) << endl;
 556
 557             pointField mappedOldPointsNew(allPoints().size());
 558
 559             mappedOldPointsNew.map(oldPointsNew, topoChangeMap4->pointMap());
 560
 561             forAll(valves(), valveI)
 562             {
 563                 const labelList& cutPointsAddressing =
 564                     pointZones()
 565                     [
 566                         pointZones().findZoneID
 567                         (
 568                             "cutPointsV"
 569                             + Foam::name(valveI + 1)
 570                         )
 571                     ];
 572
 573                 forAll(cutPointsAddressing, i)
 574                 {
 575                     mappedOldPointsNew[cutPointsAddressing[i]] =
 576                         newPoints[cutPointsAddressing[i]];
 577                 }
 578             }
 579
 580             pointField newPoints = allPoints();
 581
 582             movePoints(mappedOldPointsNew);
 583
 584             resetMotion();
 585             setV0();
 586             movePoints(newPoints);
 587
 588             Info<< "max mesh phi, sliding mesh attached = "
 589                 << max(mag(phi().internalField())) << endl;
 590             Info<< "max mesh phi, sliding mesh attached = "
 591                 << min(mag(phi().internalField())) << endl;
 592
 593             Info<< "AFTER 2 mag(points - oldPoints)/deltaT = "
 594                 << max(mag(points() - oldPoints()))/engTime().deltaT().value()
 595                 << endl;
 596             Info << "AFTER 2 max V()-V0() after = "
 597                 << max(mag(V() - V0())) << endl;
 598
 599         }
 600     }
 601
 602     return true;
 603 }
 604
 605