src/topoChangerFvMesh/movingConeTopoFvMesh/movingConeTopoFvMesh.C

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