src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
   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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "movingConeTopoFvMesh.H"
  27 #include "Time.H"
  28 #include "mapPolyMesh.H"
  29 #include "layerAdditionRemoval.H"
  30 #include "addToRunTimeSelectionTable.H"
  31 #include "meshTools.H"
  32 #include "OFstream.H"
  33
  34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38     defineTypeNameAndDebug(movingConeTopoFvMesh, 0);
  39
  40     addToRunTimeSelectionTable
  41     (
  42         topoChangerFvMesh,
  43         movingConeTopoFvMesh,
  44         IOobject
  45     );
  46 }
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51 Foam::tmp<Foam::scalarField> Foam::movingConeTopoFvMesh::vertexMarkup
  52 (
  53     const pointField& p,
  54     const scalar curLeft,
  55     const scalar curRight
  56 ) const
  57 {
  58     Info<< "Updating vertex markup.  curLeft: "
  59         << curLeft << " curRight: " << curRight << endl;
  60
  61     tmp<scalarField> tvertexMarkup(new scalarField(p.size()));
  62     scalarField& vertexMarkup = tvertexMarkup();
  63
  64     forAll(p, pI)
  65     {
  66         if (p[pI].x() < curLeft - SMALL)
  67         {
  68             vertexMarkup[pI] = -1;
  69         }
  70         else if (p[pI].x() > curRight + SMALL)
  71         {
  72             vertexMarkup[pI] = 1;
  73         }
  74         else
  75         {
  76             vertexMarkup[pI] = 0;
  77         }
  78     }
  79
  80     return tvertexMarkup;
  81 }
  82
  83
  84 void Foam::movingConeTopoFvMesh::addZonesAndModifiers()
  85 {
  86     // Add zones and modifiers for motion action
  87
  88     if
  89     (
  90         pointZones().size()
  91      || faceZones().size()
  92      || cellZones().size()
  93      || topoChanger_.size()
  94     )
  95     {
  96         Info<< "void movingConeTopoFvMesh::addZonesAndModifiers() : "
  97             << "Zones and modifiers already present.  Skipping."
  98             << endl;
  99
 100         return;
 101     }
 102
 103     Info<< "Time = " << time().timeName() << endl
 104         << "Adding zones and modifiers to the mesh" << endl;
 105
 106     const vectorField& fc = faceCentres();
 107     const vectorField& fa = faceAreas();
 108
 109     labelList zone1(fc.size());
 110     boolList flipZone1(fc.size(), false);
 111     label nZoneFaces1 = 0;
 112
 113     labelList zone2(fc.size());
 114     boolList flipZone2(fc.size(), false);
 115     label nZoneFaces2 = 0;
 116
 117     forAll(fc, faceI)
 118     {
 119         if
 120         (
 121             fc[faceI].x() > -0.003501
 122          && fc[faceI].x() < -0.003499
 123         )
 124         {
 125             if ((fa[faceI] & vector(1, 0, 0)) < 0)
 126             {
 127                 flipZone1[nZoneFaces1] = true;
 128             }
 129
 130             zone1[nZoneFaces1] = faceI;
 131             Info<< "face " << faceI << " for zone 1.  Flip: "
 132                 << flipZone1[nZoneFaces1] << endl;
 133             nZoneFaces1++;
 134         }
 135         else if
 136         (
 137             fc[faceI].x() > -0.00701
 138          && fc[faceI].x() < -0.00699
 139         )
 140         {
 141             zone2[nZoneFaces2] = faceI;
 142
 143             if ((fa[faceI] & vector(1, 0, 0)) > 0)
 144             {
 145                 flipZone2[nZoneFaces2] = true;
 146             }
 147
 148             Info<< "face " << faceI << " for zone 2.  Flip: "
 149                 << flipZone2[nZoneFaces2] << endl;
 150             nZoneFaces2++;
 151         }
 152     }
 153
 154     zone1.setSize(nZoneFaces1);
 155     flipZone1.setSize(nZoneFaces1);
 156
 157     zone2.setSize(nZoneFaces2);
 158     flipZone2.setSize(nZoneFaces2);
 159
 160     Info<< "zone: " << zone1 << endl;
 161     Info<< "zone: " << zone2 << endl;
 162
 163     List<pointZone*> pz(0);
 164     List<faceZone*> fz(2);
 165     List<cellZone*> cz(0);
 166
 167     label nFz = 0;
 168
 169     fz[nFz] =
 170         new faceZone
 171         (
 172             "rightExtrusionFaces",
 173             zone1,
 174             flipZone1,
 175             nFz,
 176             faceZones()
 177         );
 178     nFz++;
 179
 180     fz[nFz] =
 181         new faceZone
 182         (
 183             "leftExtrusionFaces",
 184             zone2,
 185             flipZone2,
 186             nFz,
 187             faceZones()
 188         );
 189     nFz++;
 190
 191     fz.setSize(nFz);
 192
 193     Info<< "Adding mesh zones." << endl;
 194     addZones(pz, fz, cz);
 195
 196
 197     // Add layer addition/removal interfaces
 198
 199     List<polyMeshModifier*> tm(2);
 200     label nMods = 0;
 201
 202     tm[nMods] =
 203         new layerAdditionRemoval
 204         (
 205             "right",
 206             nMods,
 207             topoChanger_,
 208             "rightExtrusionFaces",
 209             readScalar
 210             (
 211                 motionDict_.subDict("right").lookup("minThickness")
 212             ),
 213             readScalar
 214             (
 215                 motionDict_.subDict("right").lookup("maxThickness")
 216             )
 217         );
 218     nMods++;
 219
 220     tm[nMods] = new layerAdditionRemoval
 221     (
 222         "left",
 223         nMods,
 224         topoChanger_,
 225         "leftExtrusionFaces",
 226         readScalar
 227         (
 228             motionDict_.subDict("left").lookup("minThickness")
 229         ),
 230         readScalar
 231         (
 232             motionDict_.subDict("left").lookup("maxThickness")
 233         )
 234     );
 235     nMods++;
 236     tm.setSize(nMods);
 237
 238     Info<< "Adding " << nMods << " mesh modifiers" << endl;
 239     topoChanger_.addTopologyModifiers(tm);
 240
 241     write();
 242 }
 243
 244
 245 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 246
 247 // Construct from components
 248 Foam::movingConeTopoFvMesh::movingConeTopoFvMesh(const IOobject& io)
 249 :
 250     topoChangerFvMesh(io),
 251     motionDict_
 252     (
 253         IOdictionary
 254         (
 255             IOobject
 256             (
 257                 "dynamicMeshDict",
 258                 time().constant(),
 259                 *this,
 260                 IOobject::MUST_READ_IF_MODIFIED,
 261                 IOobject::NO_WRITE,
 262                 false
 263             )
 264         ).subDict(typeName + "Coeffs")
 265     ),
 266     motionVelAmplitude_(motionDict_.lookup("motionVelAmplitude")),
 267     motionVelPeriod_(readScalar(motionDict_.lookup("motionVelPeriod"))),
 268     curMotionVel_
 269     (
 270         motionVelAmplitude_*
 271         Foam::sin(time().value()*M_PI/motionVelPeriod_)
 272     ),
 273     leftEdge_(readScalar(motionDict_.lookup("leftEdge"))),
 274     curLeft_(readScalar(motionDict_.lookup("leftObstacleEdge"))),
 275     curRight_(readScalar(motionDict_.lookup("rightObstacleEdge")))
 276 {
 277     Pout<< "Initial time:" << time().value()
 278         << " Initial curMotionVel_:" << curMotionVel_
 279         << endl;
 280
 281     addZonesAndModifiers();
 282
 283     curLeft_ = average
 284     (
 285         faceZones()
 286         [
 287             faceZones().findZoneID("leftExtrusionFaces")
 288         ]().localPoints()
 289     ).x() - SMALL;
 290
 291     curRight_ = average
 292     (
 293         faceZones()
 294         [
 295             faceZones().findZoneID("rightExtrusionFaces")
 296         ]().localPoints()
 297     ).x() + SMALL;
 298
 299     motionMask_ = vertexMarkup
 300     (
 301         points(),
 302         curLeft_,
 303         curRight_
 304     );
 305 }
 306
 307
 308 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 309
 310 Foam::movingConeTopoFvMesh::~movingConeTopoFvMesh()
 311 {}
 312
 313
 314 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 315
 316 bool Foam::movingConeTopoFvMesh::update()
 317 {
 318     // Do mesh changes (use inflation - put new points in topoChangeMap)
 319     autoPtr<mapPolyMesh> topoChangeMap = topoChanger_.changeMesh(true);
 320
 321     // Calculate the new point positions depending on whether the
 322     // topological change has happened or not
 323     pointField newPoints;
 324
 325     vector curMotionVel_ =
 326         motionVelAmplitude_*
 327         Foam::sin(time().value()*M_PI/motionVelPeriod_);
 328
 329     Pout<< "time:" << time().value() << " curMotionVel_:" << curMotionVel_
 330         << " curLeft:" << curLeft_ << " curRight:" << curRight_
 331         << endl;
 332
 333     if (topoChangeMap.valid())
 334     {
 335         Info<< "Topology change. Calculating motion points" << endl;
 336
 337         if (topoChangeMap().hasMotionPoints())
 338         {
 339             Info<< "Topology change. Has premotion points" << endl;
 340             //Info<< "preMotionPoints:" << topoChangeMap().preMotionPoints()
 341             //    << endl;
 342
 343             //mkDir(time().timePath());
 344             //{
 345             //    OFstream str(time().timePath()/"meshPoints.obj");
 346             //    Pout<< "Writing mesh with meshPoints to " << str.name()
 347             //        << endl;
 348             //
 349             //    const pointField& currentPoints = points();
 350             //    label vertI = 0;
 351             //    forAll(currentPoints, pointI)
 352             //    {
 353             //        meshTools::writeOBJ(str, currentPoints[pointI]);
 354             //        vertI++;
 355             //    }
 356             //    forAll(edges(), edgeI)
 357             //    {
 358             //        const edge& e = edges()[edgeI];
 359             //        str << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
 360             //    }
 361             //}
 362             //{
 363             //    OFstream str(time().timePath()/"preMotionPoints.obj");
 364             //    Pout<< "Writing mesh with preMotionPoints to " << str.name()
 365             //        << endl;
 366             //
 367             //    const pointField& newPoints =
 368             //        topoChangeMap().preMotionPoints();
 369             //    label vertI = 0;
 370             //    forAll(newPoints, pointI)
 371             //    {
 372             //        meshTools::writeOBJ(str, newPoints[pointI]);
 373             //        vertI++;
 374             //    }
 375             //    forAll(edges(), edgeI)
 376             //    {
 377             //        const edge& e = edges()[edgeI];
 378             //        str << "l " << e[0]+1 << ' ' << e[1]+1 << nl;
 379             //    }
 380             //}
 381
 382
 383             motionMask_ =
 384                 vertexMarkup
 385                 (
 386                     topoChangeMap().preMotionPoints(),
 387                     curLeft_,
 388                     curRight_
 389                 );
 390
 391             // Move points inside the motionMask
 392             newPoints =
 393                 topoChangeMap().preMotionPoints()
 394               + (
 395                     pos(0.5 - mag(motionMask_)) // cells above the body
 396                 )*curMotionVel_*time().deltaT().value();
 397         }
 398         else
 399         {
 400             Info<< "Topology change. Already set mesh points" << endl;
 401
 402             motionMask_ =
 403                 vertexMarkup
 404                 (
 405                     points(),
 406                     curLeft_,
 407                     curRight_
 408                 );
 409
 410             // Move points inside the motionMask
 411             newPoints =
 412                 points()
 413               + (
 414                     pos(0.5 - mag(motionMask_)) // cells above the body
 415                 )*curMotionVel_*time().deltaT().value();
 416         }
 417     }
 418     else
 419     {
 420         Info<< "No topology change" << endl;
 421         // Set the mesh motion
 422         newPoints =
 423             points()
 424           + (
 425                 pos(0.5 - mag(motionMask_)) // cells above the body
 426            )*curMotionVel_*time().deltaT().value();
 427     }
 428
 429     // The mesh now contains the cells with zero volume
 430     Info << "Executing mesh motion" << endl;
 431     movePoints(newPoints);
 432     //  The mesh now has got non-zero volume cells
 433
 434     curLeft_ = average
 435     (
 436         faceZones()
 437         [
 438             faceZones().findZoneID("leftExtrusionFaces")
 439         ]().localPoints()
 440     ).x() - SMALL;
 441
 442     curRight_ = average
 443     (
 444         faceZones()
 445         [
 446             faceZones().findZoneID("rightExtrusionFaces")
 447         ]().localPoints()
 448     ).x() + SMALL;
 449
 450
 451     return true;
 452 }
 453
 454
 455 // ************************************************************************* //