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