src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C

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