applications/utilities/mesh/advanced/refinementLevel/refinementLevel.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     Tries to figure out what the refinement level is on refined cartesian
  26     meshes. Run BEFORE snapping.
  27
  28     Writes
  29     - volScalarField 'refinementLevel' with current refinement level.
  30     - cellSet 'refCells' which are the cells that need to be refined to satisfy
  31       2:1 refinement.
  32
  33     Works by dividing cells into volume bins.
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "argList.H"
  38 #include "Time.H"
  39 #include "polyMesh.H"
  40 #include "cellSet.H"
  41 #include "SortableList.H"
  42 #include "labelIOList.H"
  43 #include "fvMesh.H"
  44 #include "volFields.H"
  45
  46 using namespace Foam;
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49
  50 // Return true if any cells had to be split to keep a difference between
  51 // neighbouring refinement levels < limitDiff. Puts cells into refCells and
  52 // update refLevel to account for refinement.
  53 bool limitRefinementLevel
  54 (
  55     const primitiveMesh& mesh,
  56     labelList& refLevel,
  57     cellSet& refCells
  58 )
  59 {
  60     const labelListList& cellCells = mesh.cellCells();
  61
  62     label oldNCells = refCells.size();
  63
  64     forAll(cellCells, cellI)
  65     {
  66         const labelList& cCells = cellCells[cellI];
  67
  68         forAll(cCells, i)
  69         {
  70             if (refLevel[cCells[i]] > (refLevel[cellI]+1))
  71             {
  72                 // Found neighbour with >=2 difference in refLevel.
  73                 refCells.insert(cellI);
  74                 refLevel[cellI]++;
  75                 break;
  76             }
  77         }
  78     }
  79
  80     if (refCells.size() > oldNCells)
  81     {
  82         Info<< "Added an additional " << refCells.size() - oldNCells
  83             << " cells to satisfy 1:2 refinement level"
  84             << endl;
  85
  86         return true;
  87     }
  88     else
  89     {
  90         return false;
  91     }
  92 }
  93
  94
  95 // Main program:
  96
  97 int main(int argc, char *argv[])
  98 {
  99     argList::addBoolOption
 100     (
 101         "readLevel",
 102         "read level from refinementLevel file"
 103     );
 104
 105     #include "setRootCase.H"
 106     #include "createTime.H"
 107     #include "createPolyMesh.H"
 108
 109     Info<< "Dividing cells into bins depending on cell volume.\nThis will"
 110         << " correspond to refinement levels for a mesh with only 2x2x2"
 111         << " refinement\n"
 112         << "The upper range for every bin is always 1.1 times the lower range"
 113         << " to allow for some truncation error."
 114         << nl << endl;
 115
 116     const bool readLevel = args.optionFound("readLevel");
 117
 118     const scalarField& vols = mesh.cellVolumes();
 119
 120     SortableList<scalar> sortedVols(vols);
 121
 122     // All cell labels, sorted per bin.
 123     DynamicList<DynamicList<label> > bins;
 124
 125     // Lower/upper limits
 126     DynamicList<scalar> lowerLimits;
 127     DynamicList<scalar> upperLimits;
 128
 129     // Create bin0. Have upperlimit as factor times lowerlimit.
 130     bins.append(DynamicList<label>());
 131     lowerLimits.append(sortedVols[0]);
 132     upperLimits.append(1.1 * lowerLimits.last());
 133
 134     forAll(sortedVols, i)
 135     {
 136         if (sortedVols[i] > upperLimits.last())
 137         {
 138             // New value outside of current bin
 139
 140             // Shrink old bin.
 141             DynamicList<label>& bin = bins.last();
 142
 143             bin.shrink();
 144
 145             Info<< "Collected " << bin.size() << " elements in bin "
 146                 << lowerLimits.last() << " .. "
 147                 << upperLimits.last() << endl;
 148
 149             // Create new bin.
 150             bins.append(DynamicList<label>());
 151             lowerLimits.append(sortedVols[i]);
 152             upperLimits.append(1.1 * lowerLimits.last());
 153
 154             Info<< "Creating new bin " << lowerLimits.last()
 155                 << " .. " << upperLimits.last()
 156                 << endl;
 157         }
 158
 159         // Append to current bin.
 160         DynamicList<label>& bin = bins.last();
 161
 162         bin.append(sortedVols.indices()[i]);
 163     }
 164     Info<< endl;
 165
 166     bins.last().shrink();
 167     bins.shrink();
 168     lowerLimits.shrink();
 169     upperLimits.shrink();
 170
 171
 172     //
 173     // Write to cellSets.
 174     //
 175
 176     Info<< "Volume bins:" << nl;
 177     forAll(bins, binI)
 178     {
 179         const DynamicList<label>& bin = bins[binI];
 180
 181         cellSet cells(mesh, "vol" + name(binI), bin.size());
 182
 183         forAll(bin, i)
 184         {
 185             cells.insert(bin[i]);
 186         }
 187
 188         Info<< "    " << lowerLimits[binI] << " .. " << upperLimits[binI]
 189             << "  : writing " << bin.size() << " cells to cellSet "
 190             << cells.name() << endl;
 191
 192         cells.write();
 193     }
 194
 195
 196
 197     //
 198     // Convert bins into refinement level.
 199     //
 200
 201
 202     // Construct fvMesh to be able to construct volScalarField
 203
 204     fvMesh fMesh
 205     (
 206         IOobject
 207         (
 208             fvMesh::defaultRegion,
 209             runTime.timeName(),
 210             runTime
 211         ),
 212         xferCopy(mesh.points()),   // could we safely re-use the data?
 213         xferCopy(mesh.faces()),
 214         xferCopy(mesh.cells())
 215     );
 216
 217     // Add the boundary patches
 218     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 219
 220     List<polyPatch*> p(patches.size());
 221
 222     forAll(p, patchI)
 223     {
 224         p[patchI] = patches[patchI].clone(fMesh.boundaryMesh()).ptr();
 225     }
 226
 227     fMesh.addFvPatches(p);
 228
 229
 230     // Refinement level
 231     IOobject refHeader
 232     (
 233         "refinementLevel",
 234         runTime.timeName(),
 235         polyMesh::defaultRegion,
 236         runTime
 237     );
 238
 239     if (!readLevel && refHeader.headerOk())
 240     {
 241         WarningIn(args.executable())
 242             << "Detected " << refHeader.name() << " file in "
 243             << polyMesh::defaultRegion <<  " directory. Please remove to"
 244             << " recreate it or use the -readLevel option to use it"
 245             << endl;
 246         return 1;
 247     }
 248
 249
 250     labelIOList refLevel
 251     (
 252         IOobject
 253         (
 254             "refinementLevel",
 255             runTime.timeName(),
 256             mesh,
 257             IOobject::NO_READ,
 258             IOobject::AUTO_WRITE
 259         ),
 260         labelList(mesh.nCells(), 0)
 261     );
 262
 263     if (readLevel)
 264     {
 265         refLevel = labelIOList(refHeader);
 266     }
 267
 268     // Construct volScalarField with same info for post processing
 269     volScalarField postRefLevel
 270     (
 271         IOobject
 272         (
 273             "refinementLevel",
 274             runTime.timeName(),
 275             mesh,
 276             IOobject::NO_READ,
 277             IOobject::NO_WRITE
 278         ),
 279         fMesh,
 280         dimensionedScalar("zero", dimless/dimTime, 0)
 281     );
 282
 283     // Set cell values
 284     forAll(bins, binI)
 285     {
 286         const DynamicList<label>& bin = bins[binI];
 287
 288         forAll(bin, i)
 289         {
 290             refLevel[bin[i]] = bins.size() - binI - 1;
 291             postRefLevel[bin[i]] = refLevel[bin[i]];
 292         }
 293     }
 294
 295
 296     // For volScalarField: set boundary values to same as cell.
 297     // Note: could also put
 298     // zeroGradient b.c. on postRefLevel and do evaluate.
 299     forAll(postRefLevel.boundaryField(), patchI)
 300     {
 301         const polyPatch& pp = patches[patchI];
 302
 303         fvPatchScalarField& bField = postRefLevel.boundaryField()[patchI];
 304
 305         Info<< "Setting field for patch "<< endl;
 306
 307         forAll(bField, faceI)
 308         {
 309             label own = mesh.faceOwner()[pp.start() + faceI];
 310
 311             bField[faceI] = postRefLevel[own];
 312         }
 313     }
 314
 315     Info<< "Determined current refinement level and writing to "
 316         << postRefLevel.name() << " (as volScalarField; for post processing)"
 317         << nl
 318         << polyMesh::defaultRegion/refLevel.name()
 319         << " (as labelIOList; for meshing)" << nl
 320         << endl;
 321
 322     refLevel.write();
 323     postRefLevel.write();
 324
 325
 326     // Find out cells to refine to keep to 2:1 refinement level restriction
 327
 328     // Cells to refine
 329     cellSet refCells(mesh, "refCells", 100);
 330
 331     while
 332     (
 333         limitRefinementLevel
 334         (
 335             mesh,
 336             refLevel,       // current refinement level
 337             refCells        // cells to refine
 338         )
 339     )
 340     {}
 341
 342     if (refCells.size())
 343     {
 344         Info<< "Collected " << refCells.size() << " cells that need to be"
 345             << " refined to get closer to overall 2:1 refinement level limit"
 346             << nl
 347             << "Written cells to be refined to cellSet " << refCells.name()
 348             << nl << endl;
 349
 350         refCells.write();
 351
 352         Info<< "After refinement this tool can be run again to see if the 2:1"
 353             << " limit is observed all over the mesh" << nl << endl;
 354     }
 355     else
 356     {
 357         Info<< "All cells in the mesh observe the 2:1 refinement level limit"
 358             << nl << endl;
 359     }
 360
 361     Info<< "\nEnd\n" << endl;
 362     return 0;
 363 }
 364
 365
 366 // ************************************************************************* //