src/engine/engineTopoChangerMesh/thoboisMesh/thoboisMeshMove.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 "thoboisMesh.H"
  28 #include "slidingInterface.H"
  29 #include "layerAdditionRemoval.H"
  30 #include "attachDetach.H"
  31 #include "surfaceFields.H"
  32 #include "regionSplit.H"
  33 #include "componentMixedTetPolyPatchVectorField.H"
  34 #include "mapPolyMesh.H"
  35
  36 #include "tetPolyMesh.H"
  37 #include "tetPointFields.H"
  38 #include "elementFields.H"
  39 #include "fixedValueTetPolyPatchFields.H"
  40 #include "slipTetPolyPatchFields.H"
  41 #include "symmetryTetPolyPatch.H"
  42
  43 #include "tetFem.H"
  44
  45 #include "symmetryFvPatch.H"
  46 #include "wedgeFvPatch.H"
  47 #include "emptyFvPatch.H"
  48 #include "zeroGradientTetPolyPatchFields.H"
  49 #include "tetDecompositionMotionSolver.H"
  50
  51 #include "fixedValueTetPolyPatchFields.H"
  52 #include "mixedTetPolyPatchFields.H"
  53 #include "slipTetPolyPatchFields.H"
  54 #include "zeroGradientTetPolyPatchFields.H"
  55
  56
  57
  58 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  59 void Foam::thoboisMesh::makeLayersLive()
  60 {
  61     // Enable layering
  62     forAll (topoChanger_, modI)
  63     {
  64         if (isA<layerAdditionRemoval>(topoChanger_[modI]))
  65         {
  66             topoChanger_[modI].enable();
  67         }
  68         else if (isA<slidingInterface>(topoChanger_[modI]))
  69         {
  70             topoChanger_[modI].disable();
  71         }
  72         else if (isA<attachDetach>(topoChanger_[modI]))
  73         {
  74             topoChanger_[modI].enable();
  75         }
  76         else
  77         {
  78             FatalErrorIn("void Foam::engineTopoFvMesh::makeLayersLive()")
  79                 << "Don't know what to do with mesh modifier "
  80                 << modI << " of type " << topoChanger_[modI].type()
  81                 << abort(FatalError);
  82         }
  83     }
  84 }
  85
  86
  87
  88 void Foam::thoboisMesh::valveDetach()
  89 {
  90     // Enable sliding interface
  91     forAll (topoChanger_, modI)
  92     {
  93         if (isA<attachDetach>(topoChanger_[modI]))
  94         {
  95             const attachDetach& ad =
  96                 refCast<const attachDetach>(topoChanger_[modI]);
  97
  98             const word masterName = ad.masterPatchID().name();
  99
 100             // Find the valve with that name
 101             label valveIndex = -1;
 102
 103             forAll (valves_, valveI)
 104             {
 105                 if
 106                 (
 107                     valves_[valveI].detachInCylinderPatchID().name()
 108                  == masterName
 109                 )
 110                 {
 111                     valveIndex = valveI;
 112                     break;
 113                 }
 114             }
 115
 116             if (valveIndex < 0)
 117             {
 118                 FatalErrorIn
 119                 (
 120                     "void Foam::engineTopoFvMesh::prepareValveDetach()"
 121                 )   << "Cannot match patch for attach/detach " << modI
 122                     << abort(FatalError);
 123             }
 124
 125             if (debug)
 126             {
 127                 Info<< " valveI: " << valveIndex << " attached: "
 128                     << ad.attached()
 129                     << " valve open: " << valves_[valveIndex].isOpen()
 130                     << endl;
 131             }
 132
 133             ad.setDetach();
 134
 135         }
 136     }
 137 }
 138
 139 void Foam::thoboisMesh::valveAttach()
 140 {
 141     // Enable sliding interface
 142     forAll (topoChanger_, modI)
 143     {
 144         if (isA<attachDetach>(topoChanger_[modI]))
 145         {
 146             const attachDetach& ad =
 147                 refCast<const attachDetach>(topoChanger_[modI]);
 148
 149             const word masterName = ad.masterPatchID().name();
 150
 151             // Find the valve with that name
 152             label valveIndex = -1;
 153
 154             forAll (valves_, valveI)
 155             {
 156                 if
 157                 (
 158                     valves_[valveI].detachInCylinderPatchID().name()
 159                  == masterName
 160                 )
 161                 {
 162                     valveIndex = valveI;
 163                     break;
 164                 }
 165             }
 166
 167             if (valveIndex < 0)
 168             {
 169                 FatalErrorIn
 170                 (
 171                     "void Foam::engineTopoFvMesh::prepareValveDetach()"
 172                 )   << "Cannot match patch for attach/detach " << modI
 173                     << abort(FatalError);
 174             }
 175
 176             if (debug)
 177             {
 178                 Info<< " valveI: " << valveIndex << " attached: "
 179                     << ad.attached()
 180                     << " valve open: " << valves_[valveIndex].isOpen()
 181                     << endl;
 182             }
 183
 184             ad.setAttach();
 185
 186         }
 187     }
 188
 189 }
 190
 191 void Foam::thoboisMesh::prepareValveDetach()
 192 {
 193     // Enable sliding interface
 194     forAll (topoChanger_, modI)
 195     {
 196         if (isA<attachDetach>(topoChanger_[modI]))
 197         {
 198             const attachDetach& ad =
 199                 refCast<const attachDetach>(topoChanger_[modI]);
 200
 201             const word masterName = ad.masterPatchID().name();
 202
 203             // Find the valve with that name
 204             label valveIndex = -1;
 205
 206             forAll (valves_, valveI)
 207             {
 208                 if
 209                 (
 210                     valves_[valveI].detachInCylinderPatchID().name()
 211                  == masterName
 212                 )
 213                 {
 214                     valveIndex = valveI;
 215                     break;
 216                 }
 217             }
 218
 219             if (valveIndex < 0)
 220             {
 221                 FatalErrorIn
 222                 (
 223                     "void Foam::engineTopoFvMesh::prepareValveDetach()"
 224                 )   << "Cannot match patch for attach/detach " << modI
 225                     << abort(FatalError);
 226             }
 227
 228             if (debug)
 229             {
 230                 Info<< " valveI: " << valveIndex << " attached: "
 231                     << ad.attached()
 232                     << " valve open: " << valves_[valveIndex].isOpen()
 233                     << endl;
 234             }
 235
 236             if (valves_[valveIndex].isOpen())
 237             {
 238                 ad.setAttach();
 239             }
 240             else
 241             {
 242                 ad.setDetach();
 243             }
 244         }
 245     }
 246 }
 247
 248
 249 bool Foam::thoboisMesh::update()
 250 {
 251     Info << "bool Foam::layerSmooth::update()" << endl;
 252     tetDecompositionMotionSolver& mSolver =
 253         refCast<tetDecompositionMotionSolver>(msPtr_());
 254
 255     tetPointVectorField& motionU = mSolver.motionU();
 256
 257     // Find piston mesh modifier
 258
 259     const label pistonLayerID =
 260         topoChanger_.findModifierID("pistonLayer");
 261
 262     {
 263         valveDetach();
 264
 265         topoChanger_[pistonLayerID].disable();
 266
 267         autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 268
 269         if (topoChangeMap->morphing())
 270         {
 271             mSolver.updateMesh(topoChangeMap());
 272             Info << "mSolver.updateMesh(topoChangeMap())" << endl;
 273         }
 274
 275     }
 276
 277     scalar deltaZ = engTime().pistonDisplacement().value();
 278
 279     // deltaZ set to zero, FIXED PISTON POSITION
 280     deltaZ = 0.0;
 281
 282     Info<< "deltaZ = " << deltaZ << " Piston at:" << pistonPosition()
 283         << endl;
 284     virtualPistonPosition() += deltaZ;
 285
 286     Info << "pistonLayerID: " << pistonLayerID << endl;
 287
 288     bool realDeformation = deformation();
 289
 290     if
 291     (
 292         virtualPistonPosition() + deltaZ
 293       > deckHeight() - engTime().clearance().value() - SMALL
 294     )
 295     {
 296         realDeformation = true;
 297     }
 298
 299     if (realDeformation)
 300     {
 301         // Disable layer addition
 302         Info << "**Mesh deformation mode" << endl;
 303         topoChanger_[pistonLayerID].disable();
 304     }
 305     else
 306     {
 307         // Enable layer addition
 308         Info << "**Piston layering mode" << endl;
 309         topoChanger_[pistonLayerID].enable();
 310     }
 311
 312     autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 313
 314     bool meshChanged = topoChangeMap->morphing();
 315
 316     pointField newPoints = allPoints();
 317
 318     if (meshChanged)
 319     {
 320         mSolver.updateMesh(topoChangeMap());
 321
 322         if (topoChangeMap().hasMotionPoints())
 323         {
 324             movePoints(topoChangeMap().preMotionPoints());
 325             newPoints = topoChangeMap().preMotionPoints();
 326         }
 327         setV0();
 328         resetMotion();
 329     }
 330
 331     if (!deformation())
 332     {
 333 //#       include "movePistonPointsLayering.H"
 334         movePoints(newPoints);
 335
 336 #       include "setValveMotionBoundaryConditionThobois.H"
 337
 338         // Set piston velocity
 339         if (piston().patchID().active())
 340         {
 341
 342             componentMixedTetPolyPatchVectorField& pp =
 343                 refCast<componentMixedTetPolyPatchVectorField>
 344                 (
 345                     motionU.boundaryField()[piston().patchID().index()]
 346                 );
 347
 348             pp.refValue() = vector::zero;
 349         }
 350
 351         motionU.correctBoundaryConditions();
 352
 353
 354         DynamicList<label> constrainedPoints(mSolver.curPoints()().size()/100);
 355         DynamicList<vector> constrainedVelocity
 356         (
 357             mSolver.curPoints()().size()/100
 358         );
 359
 360 #       include "setThoboisMeshConstraintsNoDeformation.H"
 361
 362
 363         // Set piston velocity
 364         if (piston().patchID().active())
 365         {
 366             componentMixedTetPolyPatchVectorField& pp =
 367                 refCast<componentMixedTetPolyPatchVectorField>
 368                 (
 369                     motionU.boundaryField()[piston().patchID().index()]
 370                 );
 371
 372             pp.refValue() = vector::zero;
 373         }
 374
 375         motionU.correctBoundaryConditions();
 376
 377         labelList constrainedPointsList(constrainedPoints.shrink());
 378         vectorField constrainedVelocityField(constrainedVelocity.shrink());
 379
 380         forAll (constrainedPointsList, i)
 381         {
 382             mSolver.setConstraint
 383             (
 384                 constrainedPointsList[i],
 385                 constrainedVelocityField[i]
 386             );
 387         }
 388
 389 //         mSolver.solve
 390 //         (
 391 //             labelList(constrainedPoints.shrink()),
 392 //             constrainedVelocity.shrink()
 393 //         );
 394
 395         mSolver.solve();
 396
 397         //set to zero the motion U along the x and y directions
 398
 399         newPoints = mSolver.curPoints();
 400         movePoints(newPoints);
 401         setVirtualPositions();
 402         mSolver.clearConstraints();
 403     }
 404     else
 405     {
 406 #       include "setValveMotionBoundaryConditionThobois.H"
 407 #       include "setPistonMotionBoundaryConditionThobois.H"
 408         DynamicList<label> constrainedPoints(mSolver.curPoints()().size()/100);
 409         DynamicList<vector> constrainedVelocity
 410         (
 411             mSolver.curPoints()().size()/100
 412         );
 413
 414 #       include "setThoboisMeshConstraints.H"
 415
 416         // Set piston velocity
 417         if (piston().patchID().active())
 418         {
 419
 420             componentMixedTetPolyPatchVectorField& pp =
 421                 refCast<componentMixedTetPolyPatchVectorField>
 422                 (
 423                     motionU.boundaryField()[piston().patchID().index()]
 424                 );
 425
 426             pp.refValue() = vector::zero;
 427         }
 428
 429 //        motionU.correctBoundaryConditions();
 430
 431         labelList constrainedPointsList(constrainedPoints.shrink());
 432         vectorField constrainedVelocityField(constrainedVelocity.shrink());
 433
 434         forAll (constrainedPointsList, i)
 435         {
 436             mSolver.setConstraint
 437             (
 438                 constrainedPointsList[i],
 439                 constrainedVelocityField[i]
 440             );
 441         }
 442
 443 //        mSolver.solve
 444 //        (
 445 //            labelList(constrainedPoints.shrink()),
 446 //            constrainedVelocity.shrink()
 447 //        );
 448
 449         mSolver.solve();
 450
 451         //set to zero the motion U along the x and y directions
 452
 453         newPoints = mSolver.curPoints();
 454         movePoints(newPoints);
 455         setVirtualPositions();
 456         mSolver.clearConstraints();
 457     }
 458
 459     {
 460         pointField oldPointsNew = oldPoints();
 461         pointField newPointsNew = points();
 462
 463         prepareValveDetach();
 464         topoChanger_[pistonLayerID].disable();
 465         autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 466
 467 //        Info << motionU << endl;
 468
 469         if (topoChangeMap->morphing())
 470         {
 471             mSolver.updateMesh(topoChangeMap());
 472
 473            {
 474                // correct the motion after attaching the sliding interface
 475
 476                pointField mappedOldPointsNew(allPoints().size());
 477
 478                mappedOldPointsNew.map(oldPointsNew, topoChangeMap->pointMap());
 479                pointField newPoints = allPoints();
 480
 481                movePoints(mappedOldPointsNew);
 482                resetMotion();
 483                setV0();
 484                movePoints(newPoints);
 485            }
 486         }
 487     }
 488
 489     Info << "CHANGED LAST" << endl;
 490
 491 /*
 492     if(moving() && meshChanged)
 493     {
 494         Info << "MOOOOOOOOOOVING!!!!!!!" << endl;
 495         Info << "min V0() post motion = " << min(V0()) << endl;
 496     }
 497
 498     Info << "min V() post-motion = " << min(V()) << endl;
 499     Info << "max phi() post-motion = " << max(phi()) << endl;
 500     Info << "min phi() post-motion = " << min(phi()) << endl;
 501 */
 502
 503     Info<< "Total cylinder volume at CA " << engTime().timeName() << " = "
 504         << sum(V()) << endl;
 505
 506     return meshChanged;
 507 }