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