src/engine/engineTopoChangerMesh/thoboisMesh/thoboisMeshMove.C

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