src/dynamicMesh/polyMeshModifiers/layerAdditionRemoval/layerAdditionRemoval.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 Description
  26     Cell layer addition/removal mesh modifier
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "layerAdditionRemoval.H"
  31 #include "polyTopoChanger.H"
  32 #include "polyMesh.H"
  33 #include "Time.H"
  34 #include "primitiveMesh.H"
  35 #include "polyTopoChange.H"
  36 #include "addToRunTimeSelectionTable.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 namespace Foam
  41 {
  42     defineTypeNameAndDebug(layerAdditionRemoval, 0);
  43     addToRunTimeSelectionTable
  44     (
  45         polyMeshModifier,
  46         layerAdditionRemoval,
  47         dictionary
  48     );
  49 }
  50
  51
  52 const Foam::scalar Foam::layerAdditionRemoval::addDelta_ = 0.3;
  53 const Foam::scalar Foam::layerAdditionRemoval::removeDelta_ = 0.1;
  54
  55
  56 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  57
  58 void Foam::layerAdditionRemoval::checkDefinition()
  59 {
  60     if (!faceZoneID_.active())
  61     {
  62         FatalErrorIn
  63         (
  64             "void Foam::layerAdditionRemoval::checkDefinition()"
  65         )   << "Master face zone named " << faceZoneID_.name()
  66             << " cannot be found."
  67             << abort(FatalError);
  68     }
  69
  70     if
  71     (
  72         minLayerThickness_ < VSMALL
  73      || maxLayerThickness_ < minLayerThickness_
  74     )
  75     {
  76         FatalErrorIn
  77         (
  78             "void Foam::layerAdditionRemoval::checkDefinition()"
  79         )   << "Incorrect layer thickness definition."
  80             << abort(FatalError);
  81     }
  82
  83     if (topoChanger().mesh().faceZones()[faceZoneID_.index()].size() == 0)
  84     {
  85         FatalErrorIn
  86         (
  87             "void Foam::layerAdditionRemoval::checkDefinition()"
  88         )   << "Face extrusion zone contains no faces.  Please check your "
  89             << "mesh definition."
  90             << abort(FatalError);
  91     }
  92
  93     if (debug)
  94     {
  95         Pout<< "Cell layer addition/removal object " << name() << " :" << nl
  96             << "    faceZoneID: " << faceZoneID_ << endl;
  97     }
  98 }
  99
 100 Foam::scalar Foam::layerAdditionRemoval::readOldThickness
 101 (
 102     const dictionary& dict
 103 )
 104 {
 105     if (dict.found("oldLayerThickness"))
 106     {
 107         return readScalar(dict.lookup("oldLayerThickness"));
 108     }
 109     else
 110     {
 111         return -1.0;
 112     }
 113 }
 114
 115
 116 void Foam::layerAdditionRemoval::clearAddressing() const
 117 {
 118     // Layer removal data
 119     deleteDemandDrivenData(pointsPairingPtr_);
 120     deleteDemandDrivenData(facesPairingPtr_);
 121 }
 122
 123
 124 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 125
 126 // Construct from components
 127 Foam::layerAdditionRemoval::layerAdditionRemoval
 128 (
 129     const word& name,
 130     const label index,
 131     const polyTopoChanger& mme,
 132     const word& zoneName,
 133     const scalar minThickness,
 134     const scalar maxThickness
 135 )
 136 :
 137     polyMeshModifier(name, index, mme, true),
 138     faceZoneID_(zoneName, mme.mesh().faceZones()),
 139     minLayerThickness_(minThickness),
 140     maxLayerThickness_(maxThickness),
 141     oldLayerThickness_(-1.0),
 142     pointsPairingPtr_(NULL),
 143     facesPairingPtr_(NULL),
 144     triggerRemoval_(-1),
 145     triggerAddition_(-1)
 146 {
 147     checkDefinition();
 148 }
 149
 150
 151 // Construct from dictionary
 152 Foam::layerAdditionRemoval::layerAdditionRemoval
 153 (
 154     const word& name,
 155     const dictionary& dict,
 156     const label index,
 157     const polyTopoChanger& mme
 158 )
 159 :
 160     polyMeshModifier(name, index, mme, Switch(dict.lookup("active"))),
 161     faceZoneID_(dict.lookup("faceZoneName"), mme.mesh().faceZones()),
 162     minLayerThickness_(readScalar(dict.lookup("minLayerThickness"))),
 163     maxLayerThickness_(readScalar(dict.lookup("maxLayerThickness"))),
 164     oldLayerThickness_(readOldThickness(dict)),
 165     pointsPairingPtr_(NULL),
 166     facesPairingPtr_(NULL),
 167     triggerRemoval_(-1),
 168     triggerAddition_(-1)
 169 {
 170     checkDefinition();
 171 }
 172
 173
 174 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 175
 176 Foam::layerAdditionRemoval::~layerAdditionRemoval()
 177 {
 178     clearAddressing();
 179 }
 180
 181
 182 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 183
 184 bool Foam::layerAdditionRemoval::changeTopology() const
 185 {
 186     // Protect from multiple calculation in the same time-step
 187     if (triggerRemoval_ > -1 || triggerAddition_ > -1)
 188     {
 189         return true;
 190     }
 191
 192     // Go through all the cells in the master layer and calculate
 193     // approximate layer thickness as the ratio of the cell volume and
 194     // face area in the face zone.
 195     // Layer addition:
 196     //     When the max thickness exceeds the threshold, trigger refinement.
 197     // Layer removal:
 198     //     When the min thickness falls below the threshold, trigger removal.
 199
 200     const faceZone& fz = topoChanger().mesh().faceZones()[faceZoneID_.index()];
 201     const labelList& mc = fz.masterCells();
 202
 203     const scalarField& V = topoChanger().mesh().cellVolumes();
 204     const vectorField& S = topoChanger().mesh().faceAreas();
 205
 206     if (min(V) < -VSMALL)
 207     {
 208         FatalErrorIn("bool layerAdditionRemoval::changeTopology() const")
 209             << "negative cell volume. Error in mesh motion before "
 210             << "topological change."
 211             << abort(FatalError);
 212     }
 213
 214     scalar avgDelta = 0;
 215     scalar minDelta = GREAT;
 216     scalar maxDelta = 0;
 217
 218     forAll (fz, faceI)
 219     {
 220         scalar curDelta = V[mc[faceI]]/mag(S[fz[faceI]]);
 221         avgDelta += curDelta;
 222         minDelta = min(minDelta, curDelta);
 223         maxDelta = max(maxDelta, curDelta);
 224     }
 225
 226     avgDelta /= fz.size();
 227
 228     ////MJ Alternative thickness determination
 229     //{
 230     //    // Edges on layer.
 231     //    const Map<label>& zoneMeshPointMap = fz().meshPointMap();
 232     //
 233     //    label nDelta = 0;
 234     //
 235     //    // Edges with only one point on zone
 236     //    const polyMesh& mesh = topoChanger().mesh();
 237     //
 238     //    forAll(mc, faceI)
 239     //    {
 240     //        const cell& cFaces = mesh.cells()[mc[faceI]];
 241     //        const edgeList cellEdges(cFaces.edges(mesh.faces()));
 242     //
 243     //        forAll(cellEdges, i)
 244     //        {
 245     //            const edge& e = cellEdges[i];
 246     //
 247     //            if (zoneMeshPointMap.found(e[0]))
 248     //            {
 249     //                if (!zoneMeshPointMap.found(e[1]))
 250     //                {
 251     //                    scalar curDelta = e.mag(mesh.points());
 252     //                    avgDelta += curDelta;
 253     //                    nDelta++;
 254     //                    minDelta = min(minDelta, curDelta);
 255     //                    maxDelta = max(maxDelta, curDelta);
 256     //                }
 257     //            }
 258     //            else
 259     //            {
 260     //                if (zoneMeshPointMap.found(e[1]))
 261     //                {
 262     //                    scalar curDelta = e.mag(mesh.points());
 263     //                    avgDelta += curDelta;
 264     //                    nDelta++;
 265     //                    minDelta = min(minDelta, curDelta);
 266     //                    maxDelta = max(maxDelta, curDelta);
 267     //                }
 268     //            }
 269     //        }
 270     //    }
 271     //
 272     //    avgDelta /= nDelta;
 273     //}
 274
 275     if (debug)
 276     {
 277         Pout<< "bool layerAdditionRemoval::changeTopology() const "
 278             << " for object " << name() << " : " << nl
 279             << "Layer thickness: min: " << minDelta
 280             << " max: " << maxDelta << " avg: " << avgDelta
 281             << " old thickness: " << oldLayerThickness_ << nl
 282             << "Removal threshold: " << minLayerThickness_
 283             << " addition threshold: " << maxLayerThickness_ << endl;
 284     }
 285
 286     bool topologicalChange = false;
 287
 288     // If the thickness is decreasing and crosses the min thickness,
 289     // trigger removal
 290     if (oldLayerThickness_ < 0)
 291     {
 292         if (debug)
 293         {
 294             Pout << "First step. No addition/removal" << endl;
 295         }
 296
 297         // No topological changes allowed before first mesh motion
 298         //
 299         oldLayerThickness_ = avgDelta;
 300
 301         topologicalChange = false;
 302     }
 303     else if (avgDelta < oldLayerThickness_)
 304     {
 305         // Layers moving towards removal
 306         if (minDelta < minLayerThickness_)
 307         {
 308             // Check layer pairing
 309             if (setLayerPairing())
 310             {
 311                 // A mesh layer detected.  Check that collapse is valid
 312                 if (validCollapse())
 313                 {
 314                     // At this point, info about moving the old mesh
 315                     // in a way to collapse the cells in the removed
 316                     // layer is available.  Not sure what to do with
 317                     // it.
 318
 319                     if (debug)
 320                     {
 321                         Pout<< "bool layerAdditionRemoval::changeTopology() "
 322                             << " const for object " << name() << " : "
 323                             << "Triggering layer removal" << endl;
 324                     }
 325
 326                     triggerRemoval_ = topoChanger().morphIndex();
 327
 328                     // Old thickness looses meaning.
 329                     // Set it up to indicate layer removal
 330                     oldLayerThickness_ = GREAT;
 331
 332                     topologicalChange = true;
 333                 }
 334                 else
 335                 {
 336                     // No removal, clear addressing
 337                     clearAddressing();
 338                 }
 339             }
 340         }
 341         else
 342         {
 343             oldLayerThickness_ = avgDelta;
 344         }
 345     }
 346     else
 347     {
 348         // Layers moving towards addition
 349         if (maxDelta > maxLayerThickness_)
 350         {
 351             if (debug)
 352             {
 353                 Pout<< "bool layerAdditionRemoval::changeTopology() const "
 354                     << " for object " << name() << " : "
 355                     << "Triggering layer addition" << endl;
 356             }
 357
 358             triggerAddition_ = topoChanger().morphIndex();
 359
 360             // Old thickness looses meaning.
 361             // Set it up to indicate layer removal
 362             oldLayerThickness_ = 0;
 363
 364             topologicalChange = true;
 365         }
 366         else
 367         {
 368             oldLayerThickness_ = avgDelta;
 369         }
 370     }
 371
 372     return topologicalChange;
 373 }
 374
 375
 376 void Foam::layerAdditionRemoval::setRefinement(polyTopoChange& ref) const
 377 {
 378     // Insert the layer addition/removal instructions
 379     // into the topological change
 380
 381     if (triggerRemoval_ == topoChanger().morphIndex())
 382     {
 383         removeCellLayer(ref);
 384
 385         // Clear addressing.  This also resets the addition/removal data
 386         if (debug)
 387         {
 388             Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
 389                 << " for object " << name() << " : "
 390                 << "Clearing addressing after layer removal. " << endl;
 391         }
 392
 393         triggerRemoval_ = -1;
 394         clearAddressing();
 395     }
 396
 397     if (triggerAddition_ == topoChanger().morphIndex())
 398     {
 399         addCellLayer(ref);
 400
 401         // Clear addressing.  This also resets the addition/removal data
 402         if (debug)
 403         {
 404             Pout<< "layerAdditionRemoval::setRefinement(polyTopoChange& ref) "
 405                 << " for object " << name() << " : "
 406                 << "Clearing addressing after layer addition. " << endl;
 407         }
 408
 409         triggerAddition_ = -1;
 410         clearAddressing();
 411     }
 412 }
 413
 414
 415 void Foam::layerAdditionRemoval::updateMesh(const mapPolyMesh&)
 416 {
 417     if (debug)
 418     {
 419         Pout<< "layerAdditionRemoval::updateMesh(const mapPolyMesh&) "
 420             << " for object " << name() << " : "
 421             << "Clearing addressing on external request. ";
 422
 423         if (pointsPairingPtr_ || facesPairingPtr_)
 424         {
 425             Pout << "Pointers set." << endl;
 426         }
 427         else
 428         {
 429             Pout << "Pointers not set." << endl;
 430         }
 431     }
 432
 433     // Mesh has changed topologically.  Update local topological data
 434     faceZoneID_.update(topoChanger().mesh().faceZones());
 435
 436     clearAddressing();
 437 }
 438
 439
 440 void Foam::layerAdditionRemoval::setMinLayerThickness(const scalar t) const
 441 {
 442     if
 443     (
 444         t < VSMALL
 445      || maxLayerThickness_ < t
 446     )
 447     {
 448         FatalErrorIn
 449         (
 450             "void layerAdditionRemoval::setMinLayerThickness("
 451             "const scalar t) const"
 452         )   << "Incorrect layer thickness definition."
 453             << abort(FatalError);
 454     }
 455
 456     minLayerThickness_ = t;
 457 }
 458
 459
 460 void Foam::layerAdditionRemoval::setMaxLayerThickness(const scalar t) const
 461 {
 462     if (t < minLayerThickness_)
 463     {
 464         FatalErrorIn
 465         (
 466             "void layerAdditionRemoval::setMaxLayerThickness("
 467             "const scalar t) const"
 468         )   << "Incorrect layer thickness definition."
 469             << abort(FatalError);
 470     }
 471
 472     maxLayerThickness_ = t;
 473 }
 474
 475
 476 void Foam::layerAdditionRemoval::write(Ostream& os) const
 477 {
 478     os  << nl << type() << nl
 479         << name()<< nl
 480         << faceZoneID_ << nl
 481         << minLayerThickness_ << nl
 482         << oldLayerThickness_ << nl
 483         << maxLayerThickness_ << endl;
 484 }
 485
 486
 487 void Foam::layerAdditionRemoval::writeDict(Ostream& os) const
 488 {
 489     os  << nl << name() << nl << token::BEGIN_BLOCK << nl
 490         << "    type " << type()
 491         << token::END_STATEMENT << nl
 492         << "    faceZoneName " << faceZoneID_.name()
 493         << token::END_STATEMENT << nl
 494         << "    minLayerThickness " << minLayerThickness_
 495         << token::END_STATEMENT << nl
 496         << "    maxLayerThickness " << maxLayerThickness_
 497         << token::END_STATEMENT << nl
 498         << "    oldLayerThickness " << oldLayerThickness_
 499         << token::END_STATEMENT << nl
 500         << "    active " << active()
 501         << token::END_STATEMENT << nl
 502         << token::END_BLOCK << endl;
 503 }
 504
 505
 506 // ************************************************************************* //