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