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