src/engine/engineTopoChangerMesh/twoStrokeEngine/twoStrokeEngineMove.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 "twoStrokeEngine.H"
  28 #include "slidingInterface.H"
  29 #include "layerAdditionRemoval.H"
  30 #include "surfaceFields.H"
  31 #include "regionSplit.H"
  32 #include "componentMixedTetPolyPatchVectorField.H"
  33 #include "mapPolyMesh.H"
  34
  35 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  36
  37 void Foam::twoStrokeEngine::makeLayersLive()
  38 {
  39     const polyTopoChanger& morphs = topoChanger_;
  40
  41     // Enable layering
  42     forAll (morphs, modI)
  43     {
  44         if (typeid(morphs[modI]) == typeid(layerAdditionRemoval))
  45         {
  46             morphs[modI].enable();
  47         }
  48         else if (typeid(morphs[modI]) == typeid(slidingInterface))
  49         {
  50             morphs[modI].disable();
  51         }
  52         else
  53         {
  54             FatalErrorIn("void Foam::twoStrokeEngine::makeLayersLive()")
  55                 << "Don't know what to do with mesh modifier "
  56                 << modI << " of type " << morphs[modI].type()
  57                 << abort(FatalError);
  58         }
  59     }
  60 }
  61
  62 void Foam::twoStrokeEngine::makeSlidersLive()
  63 {
  64     const polyTopoChanger& morphs = topoChanger_;
  65
  66     // Enable sliding interface
  67     forAll (morphs, modI)
  68     {
  69         if (typeid(morphs[modI]) == typeid(layerAdditionRemoval))
  70         {
  71             morphs[modI].disable();
  72         }
  73         else if (typeid(morphs[modI]) == typeid(slidingInterface))
  74         {
  75             morphs[modI].enable();
  76         }
  77         else
  78         {
  79             FatalErrorIn("void movingSquaresTM::makeLayersLive()")
  80                 << "Don't know what to do with mesh modifier "
  81                 << modI << " of type " << morphs[modI].type()
  82                 << abort(FatalError);
  83         }
  84     }
  85
  86 }
  87
  88 bool Foam::twoStrokeEngine::attached() const
  89 {
  90     const polyTopoChanger& morphs = topoChanger_;
  91
  92     bool result = false;
  93
  94     forAll (morphs, modI)
  95     {
  96         if (typeid(morphs[modI]) == typeid(slidingInterface))
  97         {
  98             result =
  99                 result
 100              || refCast<const slidingInterface>(morphs[modI]).attached();
 101         }
 102     }
 103
 104     // Check thal all sliders are in sync (debug only)
 105     forAll (morphs, modI)
 106     {
 107         if (typeid(morphs[modI]) == typeid(slidingInterface))
 108         {
 109             if
 110             (
 111                 result
 112              != refCast<const slidingInterface>(morphs[modI]).attached()
 113             )
 114             {
 115                 FatalErrorIn("bool movingSquaresTM::attached() const")
 116                     << "Slider " << modI << " named " << morphs[modI].name()
 117                     << " out of sync: Should be" << result
 118                     << abort(FatalError);
 119             }
 120         }
 121     }
 122
 123     return result;
 124 }
 125
 126
 127 bool Foam::twoStrokeEngine::update()
 128 {
 129     // Detaching the interface
 130     if (attached())
 131     {
 132         Info << "Decoupling sliding interfaces" << endl;
 133         makeSlidersLive();
 134
 135         // Changing topology by hand
 136         autoPtr<mapPolyMesh> topoChangeMap5 = topoChanger_.changeMesh();
 137
 138         if (topoChangeMap5->hasMotionPoints() && topoChangeMap5->morphing())
 139         {
 140             Info << "Topology change; executing pre-motion after "
 141                 << "sliding detach" << endl;
 142             movePoints(topoChangeMap5->preMotionPoints());
 143         }
 144
 145         Info << "sliding interfaces successfully decoupled!!!" << endl;
 146     }
 147     else
 148     {
 149         Info << "Sliding interfaces decoupled" << endl;
 150     }
 151
 152     Info << "Executing layer action" << endl;
 153
 154     // Piston Layering
 155
 156     makeLayersLive();
 157     // Changing topology by hand
 158
 159
 160 // /* Tommaso, 23/5/2008
 161
 162     // Find piston mesh modifier
 163     const label pistonLayerID =
 164         topoChanger_.findModifierID("pistonLayer");
 165
 166     if (pistonLayerID < 0)
 167     {
 168         FatalErrorIn("void engineFvMesh::moveAndMorph()")
 169             << "Piston modifier not found."
 170             << abort(FatalError);
 171     }
 172
 173     scalar minLayerThickness = piston().minLayer();
 174     scalar deltaZ = engTime().pistonDisplacement().value();
 175     virtualPistonPosition() += deltaZ;
 176
 177     Info << "virtualPistonPosition = " << virtualPistonPosition()
 178     << ", deckHeight = " << deckHeight()
 179     << ", pistonPosition = " << pistonPosition() << endl;
 180
 181     if (realDeformation())
 182     {
 183         // Dectivate piston layer
 184         Info << "Mesh deformation mode" << endl;
 185         topoChanger_[pistonLayerID].disable();
 186     }
 187     else
 188     {
 189         // Activate piston layer
 190         Info << "Piston layering mode" << endl;
 191         topoChanger_[pistonLayerID].enable();
 192     }
 193
 194
 195     // Changing topology by hand
 196     autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh();
 197
 198     // Work array for new points position.
 199     pointField newPoints = allPoints();
 200     const pointField& refPoints = allPoints();
 201
 202     if (topoChangeMap->morphing())
 203     {
 204         if (topoChangeMap->hasMotionPoints())
 205         {
 206             Info<< "Topology change; executing pre-motion after "
 207                 << "dynamic layering" << endl;
 208             movePoints(topoChangeMap->preMotionPoints());
 209             newPoints = topoChangeMap->preMotionPoints();
 210         }
 211
 212         setV0();
 213         resetMotion();
 214     }
 215
 216     // Reset the position of layered interfaces
 217
 218     boolList scaleDisp(nPoints(), true);
 219
 220     boolList pistonPoint(newPoints.size(), false);
 221     boolList headPoint(newPoints.size(), false);
 222
 223     // Make a list of piston and head points. HJ, 2/Dec/2009
 224
 225     labelList pistonPoints;
 226     labelList headPoints;
 227     {
 228         label pistonCellIndex = cellZones().findZoneID("pistonCells");
 229
 230         if (pistonCellIndex < 0)
 231         {
 232             FatalErrorIn("bool twoStrokeEngine::update()")
 233                 << "Cannot find cell zone pistonCells"
 234                 << abort(FatalError);
 235         }
 236
 237
 238         const labelList& pistonCells = cellZones()[pistonCellIndex];
 239
 240         label headCellIndex = cellZones().findZoneID("headCells");
 241
 242         if (headCellIndex < 0)
 243         {
 244             FatalErrorIn("bool twoStrokeEngine::update()")
 245                 << "Cannot find cell zone headCells"
 246                 << abort(FatalError);
 247         }
 248
 249         const labelList& headCells = cellZones()[headCellIndex];
 250
 251         const labelListList& cp = cellPoints();
 252
 253         boolList count(newPoints.size(), false);
 254
 255         forAll (pistonCells, cellI)
 256         {
 257             const labelList& curCellPoints = cp[pistonCells[cellI]];
 258
 259             forAll (curCellPoints, i)
 260             {
 261                 count[curCellPoints[i]] = true;
 262             }
 263         }
 264
 265         // Count the points
 266         label nCounted = 0;
 267         forAll (count, pointI)
 268         {
 269             if (count[pointI] == true)
 270             {
 271                 nCounted++;
 272             }
 273         }
 274
 275         pistonPoints.setSize(nCounted);
 276
 277         // Collect the points
 278         nCounted = 0;
 279         forAll (count, pointI)
 280         {
 281             if (count[pointI] == true)
 282             {
 283                 pistonPoints[nCounted] = pointI;
 284                 nCounted++;
 285             }
 286         }
 287
 288         // Repeat for head points
 289         count = false;
 290
 291         forAll (headCells, cellI)
 292         {
 293             const labelList& curCellPoints = cp[pistonCells[cellI]];
 294
 295             forAll (curCellPoints, i)
 296             {
 297                 count[curCellPoints[i]] = true;
 298             }
 299         }
 300
 301         // Count the points
 302         nCounted = 0;
 303         forAll (count, pointI)
 304         {
 305             if (count[pointI] == true)
 306             {
 307                 nCounted++;
 308             }
 309         }
 310
 311         headPoints.setSize(nCounted);
 312
 313         // Collect the points
 314         nCounted = 0;
 315         forAll (count, pointI)
 316         {
 317             if (count[pointI] == true)
 318             {
 319                 headPoints[nCounted] = pointI;
 320                 nCounted++;
 321             }
 322         }
 323     }
 324
 325
 326     label nScaled = nPoints();
 327
 328 //     label pistonPtsIndex = pointZones().findZoneID("pistonPoints");
 329 //     const labelList& pistonPoints = pointZones()[pistonPtsIndex];
 330
 331 //     label headPtsIndex = pointZones().findZoneID("headPoints");
 332 //     const labelList& headPoints = pointZones()[headPtsIndex];
 333
 334     const scalarField& movingPointsM = movingPointsMask();
 335
 336     forAll(pistonPoints, i)
 337     {
 338         label pointI = pistonPoints[i];
 339         pistonPoint[pointI] = true;
 340         point& p = newPoints[pointI];
 341
 342         if (p.z() < pistonPosition() - 1.0e-6)
 343         {
 344             scaleDisp[pointI] = false;
 345             nScaled--;
 346         }
 347     }
 348
 349     forAll(headPoints, i)
 350     {
 351         headPoint[headPoints[i]] = true;
 352         scaleDisp[headPoints[i]] = false;
 353         nScaled--;
 354     }
 355
 356     if (realDeformation())
 357     {
 358         forAll(scaleDisp, pointI)
 359         {
 360             point& p = newPoints[pointI];
 361
 362             if (scaleDisp[pointI])
 363             {
 364                 p.z() += movingPointsM[pointI]*deltaZ*
 365                     (deckHeight() - p.z())/(deckHeight() - pistonPosition());
 366             }
 367             else
 368             {
 369                 if (!headPoint[pointI])
 370                 {
 371                     p.z() += movingPointsM[pointI]*deltaZ;
 372                 }
 373             }
 374         }
 375     }
 376     else
 377     {
 378         // Always move piston
 379         scalar pistonTopZ = -GREAT;
 380         forAll(pistonPoints, i)
 381         {
 382             point& p = newPoints[pistonPoints[i]];
 383             p.z() += deltaZ*movingPointsM[pistonPoints[i]];
 384             pistonTopZ = max(pistonTopZ, p.z());
 385         }
 386
 387
 388         // NN! fix. only needed for compression
 389         if (deltaZ > 0.0)
 390         {
 391             // check if piston-points have moved beyond the layer above
 392             forAll(newPoints, i)
 393             {
 394                 if (!pistonPoint[i])
 395                 {
 396                     if (virtualPistonPosition() > newPoints[i].z())
 397                     {
 398                         newPoints[i].z() =
 399                         (1.0 - movingPointsM[i])*newPoints[i].z()
 400                         +
 401                         movingPointsM[i] *
 402                         (pistonTopZ + 0.9*minLayerThickness);
 403                     }
 404                 }
 405             }
 406         }
 407     }
 408
 409
 410
 411     movePoints(newPoints);
 412     deleteDemandDrivenData(movingPointsMaskPtr_);
 413
 414     pistonPosition() += deltaZ;
 415
 416 //*/ //Tommaso, 23/5/2008
 417
 418     {
 419         // Grab old points to correct the motion
 420         pointField oldPointsNew = oldAllPoints();
 421
 422         // Attach the interface
 423         Info << "Coupling sliding interfaces" << endl;
 424         makeSlidersLive();
 425
 426         // Changing topology by hand
 427         autoPtr<mapPolyMesh> topoChangeMap4 = topoChanger_.changeMesh();
 428
 429         if (topoChangeMap4->morphing())
 430         {
 431             if (topoChangeMap4->hasMotionPoints())
 432             {
 433                 Info<< "Topology change; executing pre-motion after "
 434                     << "sliding attach" << endl;
 435
 436 //                 Info<< "topoChangeMap4->preMotionPoints().size() = "
 437 //                     << topoChangeMap4->preMotionPoints().size() << nl
 438 //                     << "allPoints.size() = " << allPoints().size() << nl
 439 //                     << "points.size() = " << points().size() << endl;
 440
 441                 movePoints(topoChangeMap4->preMotionPoints());
 442                 newPoints = points();
 443             }
 444
 445             {
 446                 pointField mappedOldPointsNew(allPoints().size());
 447
 448                 mappedOldPointsNew.map
 449                 (
 450                     oldPointsNew, topoChangeMap4->pointMap()
 451                 );
 452
 453                 forAll(scavInPortPatches_, patchi)
 454                 {
 455                     const labelList& cutPointsAddressing =
 456                         pointZones()
 457                         [
 458                             pointZones().findZoneID
 459                             (
 460                                 "cutPointZone" + Foam::name(patchi + 1)
 461                             )
 462                         ];
 463
 464                     forAll(cutPointsAddressing, i)
 465                     {
 466                         mappedOldPointsNew[cutPointsAddressing[i]] =
 467                             newPoints[cutPointsAddressing[i]];
 468                     }
 469
 470                     forAll(cutPointsAddressing, i)
 471                     {
 472                         if
 473                         (
 474                             newPoints[cutPointsAddressing[i]].z()
 475                           > virtualPistonPosition()
 476                         )
 477                         {
 478                             mappedOldPointsNew[cutPointsAddressing[i]].z() =
 479                                 newPoints[cutPointsAddressing[i]].z();
 480                         }
 481                         else
 482                         {
 483                             mappedOldPointsNew[cutPointsAddressing[i]].z() =
 484                                 newPoints[cutPointsAddressing[i]].z() - deltaZ;
 485                         }
 486                     }
 487                 }
 488                 pointField newPoints = allPoints();
 489
 490                 movePoints(mappedOldPointsNew);
 491
 492                 resetMotion();
 493                 setV0();
 494                 movePoints(newPoints);
 495             }
 496         }
 497     }
 498
 499     return true;
 500 }
 501
 502
 503 // ************************************************************************* //