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