applications/utilities/mesh/manipulation/refineMesh/refineMesh.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 Description
  26     Utility to refine cells in multiple directions.
  27
  28     Either supply -all option to refine all cells (3D refinement for 3D
  29     cases; 2D for 2D cases) or reads a refineMeshDict with
  30     - cellSet to refine
  31     - directions to refine
  32
  33 \*---------------------------------------------------------------------------*/
  34
  35 #include "argList.H"
  36 #include "polyMesh.H"
  37 #include "Time.H"
  38 #include "undoableMeshCutter.H"
  39 #include "hexCellLooper.H"
  40 #include "cellSet.H"
  41 #include "twoDPointCorrector.H"
  42 #include "directions.H"
  43 #include "OFstream.H"
  44 #include "multiDirRefinement.H"
  45 #include "labelIOList.H"
  46 #include "wedgePolyPatch.H"
  47 #include "plane.H"
  48
  49 using namespace Foam;
  50
  51 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  52
  53
  54 // Max cos angle for edges to be considered aligned with axis.
  55 static const scalar edgeTol = 1E-3;
  56
  57
  58 // Calculate some edge statistics on mesh.
  59 void printEdgeStats(const primitiveMesh& mesh)
  60 {
  61     label nX = 0;
  62     label nY = 0;
  63     label nZ = 0;
  64
  65     scalar minX = GREAT;
  66     scalar maxX = -GREAT;
  67     vector x(1, 0, 0);
  68
  69     scalar minY = GREAT;
  70     scalar maxY = -GREAT;
  71     vector y(0, 1, 0);
  72
  73     scalar minZ = GREAT;
  74     scalar maxZ = -GREAT;
  75     vector z(0, 0, 1);
  76
  77     scalar minOther = GREAT;
  78     scalar maxOther = -GREAT;
  79
  80     const edgeList& edges = mesh.edges();
  81
  82     forAll(edges, edgeI)
  83     {
  84         const edge& e = edges[edgeI];
  85
  86         vector eVec(e.vec(mesh.points()));
  87
  88         scalar eMag = mag(eVec);
  89
  90         eVec /= eMag;
  91
  92         if (mag(eVec & x) > 1-edgeTol)
  93         {
  94             minX = min(minX, eMag);
  95             maxX = max(maxX, eMag);
  96             nX++;
  97         }
  98         else if (mag(eVec & y) > 1-edgeTol)
  99         {
 100             minY = min(minY, eMag);
 101             maxY = max(maxY, eMag);
 102             nY++;
 103         }
 104         else if (mag(eVec & z) > 1-edgeTol)
 105         {
 106             minZ = min(minZ, eMag);
 107             maxZ = max(maxZ, eMag);
 108             nZ++;
 109         }
 110         else
 111         {
 112             minOther = min(minOther, eMag);
 113             maxOther = max(maxOther, eMag);
 114         }
 115     }
 116
 117     Pout<< "Mesh edge statistics:" << endl
 118         << "    x aligned :  number:" << nX << "\tminLen:" << minX
 119         << "\tmaxLen:" << maxX << endl
 120         << "    y aligned :  number:" << nY << "\tminLen:" << minY
 121         << "\tmaxLen:" << maxY << endl
 122         << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
 123         << "\tmaxLen:" << maxZ << endl
 124         << "    other     :  number:" << mesh.nEdges() - nX - nY - nZ
 125         << "\tminLen:" << minOther
 126         << "\tmaxLen:" << maxOther << endl << endl;
 127 }
 128
 129
 130 // Return index of coordinate axis.
 131 label axis(const vector& normal)
 132 {
 133     label axisIndex = -1;
 134
 135     if (mag(normal & point(1, 0, 0)) > (1-edgeTol))
 136     {
 137         axisIndex = 0;
 138     }
 139     else if (mag(normal & point(0, 1, 0)) > (1-edgeTol))
 140     {
 141         axisIndex = 1;
 142     }
 143     else if (mag(normal & point(0, 0, 1)) > (1-edgeTol))
 144     {
 145         axisIndex = 2;
 146     }
 147
 148     return axisIndex;
 149 }
 150
 151
 152 //- Returns -1 or cartesian coordinate component (0=x, 1=y, 2=z) of normal
 153 //  in case of 2D mesh
 154 label twoDNess(const polyMesh& mesh)
 155 {
 156     const pointField& ctrs = mesh.cellCentres();
 157
 158     if (ctrs.size() < 2)
 159     {
 160         return -1;
 161     }
 162
 163
 164     //
 165     // 1. All cell centres on single plane aligned with x, y or z
 166     //
 167
 168     // Determine 3 points to base plane on.
 169     vector vec10 = ctrs[1] - ctrs[0];
 170     vec10 /= mag(vec10);
 171
 172     label otherCellI = -1;
 173
 174     for (label cellI = 2; cellI < ctrs.size(); cellI++)
 175     {
 176         vector vec(ctrs[cellI] - ctrs[0]);
 177         vec /= mag(vec);
 178
 179         if (mag(vec & vec10) < 0.9)
 180         {
 181             // ctrs[cellI] not in line with n
 182             otherCellI = cellI;
 183
 184             break;
 185         }
 186     }
 187
 188     if (otherCellI == -1)
 189     {
 190         // Cannot find cell to make decent angle with cell0-cell1 vector.
 191         // Note: what to do here? All cells (almost) in one line. Maybe 1D case?
 192         return -1;
 193     }
 194
 195     plane cellPlane(ctrs[0], ctrs[1], ctrs[otherCellI]);
 196
 197
 198     forAll(ctrs, cellI)
 199     {
 200         const labelList& cEdges = mesh.cellEdges()[cellI];
 201
 202         scalar minLen = GREAT;
 203
 204         forAll(cEdges, i)
 205         {
 206             minLen = min(minLen, mesh.edges()[cEdges[i]].mag(mesh.points()));
 207         }
 208
 209         if (cellPlane.distance(ctrs[cellI]) > 1E-6*minLen)
 210         {
 211             // Centres not in plane
 212             return  -1;
 213         }
 214     }
 215
 216     label axisIndex = axis(cellPlane.normal());
 217
 218     if (axisIndex == -1)
 219     {
 220         return axisIndex;
 221     }
 222
 223
 224     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 225
 226
 227     //
 228     // 2. No edges without points on boundary
 229     //
 230
 231     // Mark boundary points
 232     boolList boundaryPoint(mesh.points().size(), false);
 233
 234     forAll(patches, patchI)
 235     {
 236         const polyPatch& patch = patches[patchI];
 237
 238         forAll(patch, patchFaceI)
 239         {
 240             const face& f = patch[patchFaceI];
 241
 242             forAll(f, fp)
 243             {
 244                 boundaryPoint[f[fp]] = true;
 245             }
 246         }
 247     }
 248
 249
 250     const edgeList& edges = mesh.edges();
 251
 252     forAll(edges, edgeI)
 253     {
 254         const edge& e = edges[edgeI];
 255
 256         if (!boundaryPoint[e.start()] && !boundaryPoint[e.end()])
 257         {
 258             // Edge has no point on boundary.
 259             return -1;
 260         }
 261     }
 262
 263
 264     // 3. For all non-wedge patches: all faces either perp or aligned with
 265     //    cell-plane normal. (wedge patches already checked upon construction)
 266
 267     forAll(patches, patchI)
 268     {
 269         const polyPatch& patch = patches[patchI];
 270
 271         if
 272         (
 273             typeid(patch) != typeid(wedgePolyPatch)
 274         )
 275         {
 276             const vectorField& n = patch.faceAreas();
 277
 278             scalarField cosAngle = mag(n/mag(n) & cellPlane.normal());
 279
 280             if (mag(min(cosAngle) - max(cosAngle)) > 1E-6)
 281             {
 282                 // cosAngle should be either ~1 over all faces (2D front and
 283                 // back) or ~0 (all other patches perp to 2D)
 284                 return -1;
 285             }
 286         }
 287     }
 288
 289     return axisIndex;
 290 }
 291
 292
 293 // Main program:
 294
 295 int main(int argc, char *argv[])
 296 {
 297     Foam::argList::validOptions.insert("dict", "");
 298     Foam::argList::validOptions.insert("overwrite", "");
 299
 300 #   include "setRootCase.H"
 301 #   include "createTime.H"
 302     runTime.functionObjects().off();
 303 #   include "createPolyMesh.H"
 304     const word oldInstance = mesh.pointsInstance();
 305
 306     printEdgeStats(mesh);
 307
 308
 309     //
 310     // Read/construct control dictionary
 311     //
 312
 313     bool readDict = args.optionFound("dict");
 314     bool overwrite = args.optionFound("overwrite");
 315
 316     // List of cells to refine
 317     labelList refCells;
 318
 319     // Dictionary to control refinement
 320     dictionary refineDict;
 321
 322     if (readDict)
 323     {
 324         Info<< "Refining according to refineMeshDict" << nl << endl;
 325
 326         refineDict =
 327             IOdictionary
 328             (
 329                 IOobject
 330                 (
 331                     "refineMeshDict",
 332                     runTime.system(),
 333                     mesh,
 334                     IOobject::MUST_READ,
 335                     IOobject::NO_WRITE
 336                 )
 337             );
 338
 339         word setName(refineDict.lookup("set"));
 340
 341         cellSet cells(mesh, setName);
 342
 343         Pout<< "Read " << cells.size() << " cells from cellSet "
 344             << cells.instance()/cells.local()/cells.name()
 345             << endl << endl;
 346
 347         refCells = cells.toc();
 348     }
 349     else
 350     {
 351         Info<< "Refining all cells" << nl << endl;
 352
 353         // Select all cells
 354         refCells.setSize(mesh.nCells());
 355
 356         forAll(mesh.cells(), cellI)
 357         {
 358             refCells[cellI] = cellI;
 359         }
 360
 361
 362         // Set refinement directions based on 2D/3D
 363         label axisIndex = twoDNess(mesh);
 364
 365         if (axisIndex == -1)
 366         {
 367             Info<< "3D case; refining all directions" << nl << endl;
 368
 369             wordList directions(3);
 370             directions[0] = "tan1";
 371             directions[1] = "tan2";
 372             directions[2] = "normal";
 373             refineDict.add("directions", directions);
 374
 375             // Use hex cutter
 376             refineDict.add("useHexTopology", "true");
 377         }
 378         else
 379         {
 380             wordList directions(2);
 381
 382             if (axisIndex == 0)
 383             {
 384                 Info<< "2D case; refining in directions y,z\n" << endl;
 385                 directions[0] = "tan2";
 386                 directions[1] = "normal";
 387             }
 388             else if (axisIndex == 1)
 389             {
 390                 Info<< "2D case; refining in directions x,z\n" << endl;
 391                 directions[0] = "tan1";
 392                 directions[1] = "normal";
 393             }
 394             else
 395             {
 396                 Info<< "2D case; refining in directions x,y\n" << endl;
 397                 directions[0] = "tan1";
 398                 directions[1] = "tan2";
 399             }
 400
 401             refineDict.add("directions", directions);
 402
 403             // Use standard cutter
 404             refineDict.add("useHexTopology", "false");
 405         }
 406
 407         refineDict.add("coordinateSystem", "global");
 408
 409         dictionary coeffsDict;
 410         coeffsDict.add("tan1", vector(1, 0, 0));
 411         coeffsDict.add("tan2", vector(0, 1, 0));
 412         refineDict.add("globalCoeffs", coeffsDict);
 413
 414         refineDict.add("geometricCut", "false");
 415         refineDict.add("writeMesh", "false");
 416     }
 417
 418
 419     string oldTimeName(runTime.timeName());
 420
 421     if (!overwrite)
 422     {
 423         runTime++;
 424     }
 425
 426
 427     // Multi-directional refinement (does multiple iterations)
 428     multiDirRefinement multiRef(mesh, refCells, refineDict);
 429
 430
 431     // Write resulting mesh
 432     if (overwrite)
 433     {
 434         mesh.setInstance(oldInstance);
 435     }
 436     mesh.write();
 437
 438
 439     // Get list of cell splits.
 440     // (is for every cell in old mesh the cells they have been split into)
 441     const labelListList& oldToNew = multiRef.addedCells();
 442
 443
 444     // Create cellSet with added cells for easy inspection
 445     cellSet newCells(mesh, "refinedCells", refCells.size());
 446
 447     forAll(oldToNew, oldCellI)
 448     {
 449         const labelList& added = oldToNew[oldCellI];
 450
 451         forAll(added, i)
 452         {
 453             newCells.insert(added[i]);
 454         }
 455     }
 456
 457     Pout<< "Writing refined cells (" << newCells.size() << ") to cellSet "
 458         << newCells.instance()/newCells.local()/newCells.name()
 459         << endl << endl;
 460
 461     newCells.write();
 462
 463
 464
 465
 466     //
 467     // Invert cell split to construct map from new to old
 468     //
 469
 470     labelIOList newToOld
 471     (
 472         IOobject
 473         (
 474             "cellMap",
 475             runTime.timeName(),
 476             polyMesh::meshSubDir,
 477             mesh,
 478             IOobject::NO_READ,
 479             IOobject::AUTO_WRITE
 480         ),
 481         mesh.nCells()
 482     );
 483     newToOld.note() =
 484         "From cells in mesh at "
 485       + runTime.timeName()
 486       + " to cells in mesh at "
 487       + oldTimeName;
 488
 489
 490     forAll(oldToNew, oldCellI)
 491     {
 492         const labelList& added = oldToNew[oldCellI];
 493
 494         if (added.size())
 495         {
 496             forAll(added, i)
 497             {
 498                 newToOld[added[i]] = oldCellI;
 499             }
 500         }
 501         else
 502         {
 503             // Unrefined cell
 504             newToOld[oldCellI] = oldCellI;
 505         }
 506     }
 507
 508     Info<< "Writing map from new to old cell to "
 509         << newToOld.objectPath() << nl << endl;
 510
 511     newToOld.write();
 512
 513
 514     // Some statistics.
 515
 516     printEdgeStats(mesh);
 517
 518     Info<< "End\n" << endl;
 519
 520     return 0;
 521 }
 522
 523
 524 // ************************************************************************* //