src/engine/engineTopoChangerMesh/layerARGambit/addLayerARGambitMeshModifiers.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   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 "layerARGambit.H"
  27 #include "slidingInterface.H"
  28 #include "layerAdditionRemoval.H"
  29 #include "surfaceFields.H"
  30
  31 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  32
  33 void Foam::layerARGambit::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 layerARGambit::addZonesAndModifiers() : "
  45             << "Zones and modifiers already present.  Skipping."
  46             << endl;
  47
  48         if (topoChanger_.size() == 0)
  49         {
  50             FatalErrorIn
  51             (
  52                 "void layerARGambit::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
  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 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 (mag(dd) > 0.1)
 107             {
 108                 if (zPistV - zc > 0 && zPistV - zc < dl)
 109                 {
 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
 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
 215         pz[nPointZones] =
 216             new pointZone
 217             (
 218                 "headPoints",
 219                 headPoints.shrink(),
 220                 nPointZones,
 221                 pointZones()
 222             );
 223
 224         nPointZones++;
 225
 226         pz[nPointZones] =
 227             new pointZone
 228             (
 229                 "pistonPoints",
 230                 pistonPoints.shrink(),
 231                 nPointZones,
 232                 pointZones()
 233             );
 234
 235         nPointZones++;
 236
 237     }
 238     else if(piston().patchID().active() && offSet() <= SMALL)
 239     {
 240         label pistonPatchID = piston().patchID().index();
 241
 242         const polyPatch& pistonPatch =
 243             boundaryMesh()[piston().patchID().index()];
 244
 245         labelList pistonPatchLabels(pistonPatch.size(), pistonPatch.start());
 246
 247         forAll (pistonPatchLabels, i)
 248         {
 249             pistonPatchLabels[i] += i;
 250         }
 251
 252         fz[nFaceZones] =
 253             new faceZone
 254             (
 255                 "pistonLayerFaces",
 256                 pistonPatchLabels,
 257                 boolList(pistonPatchLabels.size(), true),
 258                 nFaceZones,
 259                 faceZones()
 260             );
 261         nFaceZones++;
 262         // Construct point zones
 263
 264         scalar zPistV = 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                 //Info<< "HeadPoint:" << pointI << " coord:" << points[pointI]
 283                 //   << endl;
 284             }
 285             else if (zCoord < zPistV + delta())
 286             {
 287                 pistonPoints.append(pointI);
 288                 //Info<< "PistonPoint:" << pointI << " coord:" << points[pointI]
 289                 //    << endl;
 290             }
 291         }
 292
 293         Info << "Number of head points = " << nHeadPoints << endl;
 294
 295
 296         pz[nPointZones] =
 297             new pointZone
 298             (
 299                 "headPoints",
 300                 headPoints.shrink(),
 301                 nPointZones,
 302                 pointZones()
 303             );
 304
 305         nPointZones++;
 306
 307         pz[nPointZones] =
 308             new pointZone
 309             (
 310                 "pistonPoints",
 311                 pistonPoints.shrink(),
 312                 nPointZones,
 313                 pointZones()
 314             );
 315
 316         nPointZones++;
 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     // Write mesh and modifiers
 351     topoChanger_.writeOpt() = IOobject::AUTO_WRITE;
 352     topoChanger_.write();
 353     write();
 354
 355     // Calculating the virtual piston position
 356     setVirtualPistonPosition();
 357
 358     Info << "virtualPistonPosition = " << virtualPistonPosition() << endl;
 359     Info << "piston position = " << pistonPosition() << endl;
 360 }
 361
 362
 363 // ************************************************************************* //