src/engine/engineTopoChangerMesh/layerARGambit/addLayerARGambitMeshModifiers.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 "layerARGambit.H"
  28 #include "slidingInterface.H"
  29 #include "layerAdditionRemoval.H"
  30 #include "surfaceFields.H"
  31
  32 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  33
  34 void Foam::layerARGambit::addZonesAndModifiers()
  35 {
  36     // Add the zones and mesh modifiers to operate piston motion
  37
  38     if
  39     (
  40         pointZones().size() > 0
  41      || faceZones().size() > 0
  42      || cellZones().size() > 0
  43     )
  44     {
  45         Info<< "void layerARGambit::addZonesAndModifiers() : "
  46             << "Zones and modifiers already present.  Skipping."
  47             << endl;
  48
  49         if (topoChanger_.size() == 0)
  50         {
  51             FatalErrorIn
  52             (
  53                 "void layerARGambit::addZonesAndModifiers()"
  54             )   << "Mesh modifiers not read properly"
  55                 << abort(FatalError);
  56         }
  57
  58         setVirtualPistonPosition();
  59         checkAndCalculate();
  60
  61         return;
  62     }
  63
  64     checkAndCalculate();
  65
  66
  67     Info<< "Time = " << engTime().theta() << endl
  68         << "Adding zones to the engine mesh" << endl;
  69
  70     //fz = 1: faces where layer are added/removed
  71     //pz = 2: points below the virtual piston faces and head points
  72
  73     List<pointZone*> pz(2);
  74     List<faceZone*> fz(1);
  75     List<cellZone*> cz(0);
  76
  77     label nPointZones = 0;
  78     label nFaceZones = 0;
  79
  80     // Add the piston zone
  81     if (piston().patchID().active() && offSet() > SMALL)
  82     {
  83
  84         // Piston position
  85
  86         label pistonPatchID = piston().patchID().index();
  87
  88         scalar zPist = max(boundary()[pistonPatchID].patch().localPoints()).z();
  89
  90         scalar zPistV = zPist + offSet();
  91
  92         labelList zone1(faceCentres().size());
  93         boolList flipZone1(faceCentres().size(), false);
  94         label nZoneFaces1 = 0;
  95
  96         bool foundAtLeastOne = false;
  97         scalar zHigher = GREAT;
  98         scalar zLower = GREAT;
  99         scalar dh = GREAT;
 100         scalar dl = GREAT;
 101
 102         forAll (faceCentres(), faceI)
 103         {
 104             scalar zc = faceCentres()[faceI].z();
 105             vector n = faceAreas()[faceI]/mag(faceAreas()[faceI]);
 106             scalar dd = n & vector(0,0,1);
 107
 108             if (mag(dd) > 0.1)
 109             {
 110                 if (zPistV - zc > 0 && zPistV - zc < dl)
 111                 {
 112                     zLower = zc;
 113                     dl = zPistV - zc;
 114                 }
 115
 116                 if (zc - zPistV > 0 && zc - zPistV < dh)
 117                 {
 118                     zHigher = zc;
 119                     dh = zc - zHigher;
 120                 }
 121
 122                 if
 123                 (
 124                     zc > zPistV - delta()
 125                     && zc < zPistV + delta()
 126                 )
 127                 {
 128                     foundAtLeastOne = true;
 129                     if ((faceAreas()[faceI] & vector(0,0,1)) < 0)
 130                     {
 131                         flipZone1[nZoneFaces1] = true;
 132                     }
 133
 134                     zone1[nZoneFaces1] = faceI;
 135                     nZoneFaces1++;
 136                 }
 137             }
 138         }
 139
 140         // if no cut was found use the layer above
 141         if (!foundAtLeastOne)
 142         {
 143
 144             zPistV = zHigher;
 145
 146             forAll (faceCentres(), faceI)
 147             {
 148                 scalar zc = faceCentres()[faceI].z();
 149                 vector n = faceAreas()[faceI]/mag(faceAreas()[faceI]);
 150                 scalar dd = n & vector(0,0,1);
 151                 if (dd > 0.1)
 152                 {
 153
 154                     if
 155                     (
 156                         zc > zPistV - delta()
 157                         && zc < zPistV + delta()
 158                     )
 159                     {
 160                         if ((faceAreas()[faceI] & vector(0,0,1)) < 0)
 161                         {
 162                             flipZone1[nZoneFaces1] = true;
 163                         }
 164
 165                         zone1[nZoneFaces1] = faceI;
 166                         nZoneFaces1++;
 167                     }
 168                 }
 169             }
 170
 171         }
 172
 173         zone1.setSize(nZoneFaces1);
 174         flipZone1.setSize(nZoneFaces1);
 175
 176         fz[nFaceZones]=
 177             new faceZone
 178             (
 179                 "pistonLayerFaces",
 180                 zone1,
 181                 flipZone1,
 182                 nFaceZones,
 183                 faceZones()
 184             );
 185
 186         nFaceZones++;
 187
 188
 189         // Construct point zones
 190
 191
 192         // Points which don't move (= cylinder head)
 193         DynamicList<label> headPoints(nPoints() / 10);
 194
 195         // Points below the piston which moves with the piston displacement
 196         DynamicList<label> pistonPoints(nPoints() / 10);
 197
 198         label nHeadPoints = 0;
 199
 200         forAll (points(), pointI)
 201         {
 202             scalar zCoord = points()[pointI].z();
 203
 204             if (zCoord > deckHeight() - delta())
 205             {
 206                 headPoints.append(pointI);
 207                 nHeadPoints++;
 208             }
 209             else if (zCoord < zPistV + delta())
 210             {
 211                 pistonPoints.append(pointI);
 212             }
 213         }
 214
 215         Info << "Number of head points = " << nHeadPoints << endl;
 216
 217
 218         pz[nPointZones] =
 219             new pointZone
 220             (
 221                 "headPoints",
 222                 headPoints.shrink(),
 223                 nPointZones,
 224                 pointZones()
 225             );
 226
 227         nPointZones++;
 228
 229         pz[nPointZones] =
 230             new pointZone
 231             (
 232                 "pistonPoints",
 233                 pistonPoints.shrink(),
 234                 nPointZones,
 235                 pointZones()
 236             );
 237
 238         nPointZones++;
 239
 240     }
 241     else if(piston().patchID().active() && offSet() <= SMALL)
 242     {
 243         label pistonPatchID = piston().patchID().index();
 244
 245         const polyPatch& pistonPatch =
 246             boundaryMesh()[piston().patchID().index()];
 247
 248         labelList pistonPatchLabels(pistonPatch.size(), pistonPatch.start());
 249
 250         forAll (pistonPatchLabels, i)
 251         {
 252             pistonPatchLabels[i] += i;
 253         }
 254
 255         fz[nFaceZones] =
 256             new faceZone
 257             (
 258                 "pistonLayerFaces",
 259                 pistonPatchLabels,
 260                 boolList(pistonPatchLabels.size(), true),
 261                 nFaceZones,
 262                 faceZones()
 263             );
 264         nFaceZones++;
 265         // Construct point zones
 266
 267         scalar zPistV = max(boundary()[pistonPatchID].patch().localPoints()).z();
 268
 269         // Points which don't move (= cylinder head)
 270         DynamicList<label> headPoints(nPoints() / 10);
 271
 272         // Points below the piston which moves with the piston displacement
 273         DynamicList<label> pistonPoints(nPoints() / 10);
 274
 275         label nHeadPoints = 0;
 276
 277         forAll (points(), pointI)
 278         {
 279             scalar zCoord = points()[pointI].z();
 280
 281             if (zCoord > deckHeight() - delta())
 282             {
 283                 headPoints.append(pointI);
 284                 nHeadPoints++;
 285                 //Info<< "HeadPoint:" << pointI << " coord:" << points[pointI]
 286                 //   << endl;
 287             }
 288             else if (zCoord < zPistV + delta())
 289             {
 290                 pistonPoints.append(pointI);
 291                 //Info<< "PistonPoint:" << pointI << " coord:" << points[pointI]
 292                 //    << endl;
 293             }
 294         }
 295
 296         Info << "Number of head points = " << nHeadPoints << endl;
 297
 298
 299         pz[nPointZones] =
 300             new pointZone
 301             (
 302                 "headPoints",
 303                 headPoints.shrink(),
 304                 nPointZones,
 305                 pointZones()
 306             );
 307
 308         nPointZones++;
 309
 310         pz[nPointZones] =
 311             new pointZone
 312             (
 313                 "pistonPoints",
 314                 pistonPoints.shrink(),
 315                 nPointZones,
 316                 pointZones()
 317             );
 318
 319         nPointZones++;
 320     }
 321
 322     Info<< "Adding " << nPointZones << " point and "
 323         << nFaceZones << " face zones" << endl;
 324
 325     pz.setSize(nPointZones);
 326     fz.setSize(nFaceZones);
 327     addZones(pz, fz, cz);
 328
 329     List<polyMeshModifier*> tm(1);
 330     label nMods = 0;
 331
 332     // Add piston layer addition
 333     Info << "Adding Layer Addition/Removal Mesh Modifier" << endl;
 334
 335     if (piston().patchID().active())
 336     {
 337         topoChanger_.setSize(1);
 338         topoChanger_.set
 339         (
 340             0,
 341             new layerAdditionRemoval
 342             (
 343                 "pistonLayer",
 344                 nMods,
 345                 topoChanger_,
 346                 "pistonLayerFaces",
 347                 piston().minLayer(),
 348                 piston().maxLayer()
 349             )
 350         );
 351     }
 352
 353     // Write mesh and modifiers
 354     topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
 355     topoChanger_.write();
 356     write();
 357
 358     // Calculating the virtual piston position
 359     setVirtualPistonPosition();
 360
 361     Info << "virtualPistonPosition = " << virtualPistonPosition() << endl;
 362     Info << "piston position = " << pistonPosition() << endl;
 363 }
 364
 365
 366 // ************************************************************************* //