src/dynamicMesh/layerAdditionRemoval/layerAdditionRemoval.C

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