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 (!isA<wedgePolyPatch>(patch))
 272         {
 273             const vectorField& n = patch.faceAreas();
 274
 275             scalarField cosAngle = mag(n/mag(n) & cellPlane.normal());
 276
 277             if (mag(min(cosAngle) - max(cosAngle)) > 1E-6)
 278             {
 279                 // cosAngle should be either ~1 over all faces (2D front and
 280                 // back) or ~0 (all other patches perp to 2D)
 281                 return -1;
 282             }
 283         }
 284     }
 285
 286     return axisIndex;
 287 }
 288
 289
 290 // Main program:
 291
 292 int main(int argc, char *argv[])
 293 {
 294     Foam::argList::validOptions.insert("dict", "");
 295     Foam::argList::validOptions.insert("overwrite", "");
 296
 297 #   include "setRootCase.H"
 298 #   include "createTime.H"
 299     runTime.functionObjects().off();
 300 #   include "createPolyMesh.H"
 301     const word oldInstance = mesh.pointsInstance();
 302
 303     printEdgeStats(mesh);
 304
 305
 306     //
 307     // Read/construct control dictionary
 308     //
 309
 310     bool readDict = args.optionFound("dict");
 311     bool overwrite = args.optionFound("overwrite");
 312
 313     // List of cells to refine
 314     labelList refCells;
 315
 316     // Dictionary to control refinement
 317     dictionary refineDict;
 318
 319     if (readDict)
 320     {
 321         Info<< "Refining according to refineMeshDict" << nl << endl;
 322
 323         refineDict =
 324             IOdictionary
 325             (
 326                 IOobject
 327                 (
 328                     "refineMeshDict",
 329                     runTime.system(),
 330                     mesh,
 331                     IOobject::MUST_READ,
 332                     IOobject::NO_WRITE
 333                 )
 334             );
 335
 336         word setName(refineDict.lookup("set"));
 337
 338         cellSet cells(mesh, setName);
 339
 340         Pout<< "Read " << cells.size() << " cells from cellSet "
 341             << cells.instance()/cells.local()/cells.name()
 342             << endl << endl;
 343
 344         refCells = cells.toc();
 345     }
 346     else
 347     {
 348         Info<< "Refining all cells" << nl << endl;
 349
 350         // Select all cells
 351         refCells.setSize(mesh.nCells());
 352
 353         forAll(mesh.cells(), cellI)
 354         {
 355             refCells[cellI] = cellI;
 356         }
 357
 358
 359         // Set refinement directions based on 2D/3D
 360         label axisIndex = twoDNess(mesh);
 361
 362         if (axisIndex == -1)
 363         {
 364             Info<< "3D case; refining all directions" << nl << endl;
 365
 366             wordList directions(3);
 367             directions[0] = "tan1";
 368             directions[1] = "tan2";
 369             directions[2] = "normal";
 370             refineDict.add("directions", directions);
 371
 372             // Use hex cutter
 373             refineDict.add("useHexTopology", "true");
 374         }
 375         else
 376         {
 377             wordList directions(2);
 378
 379             if (axisIndex == 0)
 380             {
 381                 Info<< "2D case; refining in directions y,z\n" << endl;
 382                 directions[0] = "tan2";
 383                 directions[1] = "normal";
 384             }
 385             else if (axisIndex == 1)
 386             {
 387                 Info<< "2D case; refining in directions x,z\n" << endl;
 388                 directions[0] = "tan1";
 389                 directions[1] = "normal";
 390             }
 391             else
 392             {
 393                 Info<< "2D case; refining in directions x,y\n" << endl;
 394                 directions[0] = "tan1";
 395                 directions[1] = "tan2";
 396             }
 397
 398             refineDict.add("directions", directions);
 399
 400             // Use standard cutter
 401             refineDict.add("useHexTopology", "false");
 402         }
 403
 404         refineDict.add("coordinateSystem", "global");
 405
 406         dictionary coeffsDict;
 407         coeffsDict.add("tan1", vector(1, 0, 0));
 408         coeffsDict.add("tan2", vector(0, 1, 0));
 409         refineDict.add("globalCoeffs", coeffsDict);
 410
 411         refineDict.add("geometricCut", "false");
 412         refineDict.add("writeMesh", "false");
 413     }
 414
 415
 416     string oldTimeName(runTime.timeName());
 417
 418     if (!overwrite)
 419     {
 420         runTime++;
 421     }
 422
 423
 424     // Multi-directional refinement (does multiple iterations)
 425     multiDirRefinement multiRef(mesh, refCells, refineDict);
 426
 427
 428     // Write resulting mesh
 429     if (overwrite)
 430     {
 431         mesh.setInstance(oldInstance);
 432     }
 433     mesh.write();
 434
 435
 436     // Get list of cell splits.
 437     // (is for every cell in old mesh the cells they have been split into)
 438     const labelListList& oldToNew = multiRef.addedCells();
 439
 440
 441     // Create cellSet with added cells for easy inspection
 442     cellSet newCells(mesh, "refinedCells", refCells.size());
 443
 444     forAll(oldToNew, oldCellI)
 445     {
 446         const labelList& added = oldToNew[oldCellI];
 447
 448         forAll(added, i)
 449         {
 450             newCells.insert(added[i]);
 451         }
 452     }
 453
 454     Pout<< "Writing refined cells (" << newCells.size() << ") to cellSet "
 455         << newCells.instance()/newCells.local()/newCells.name()
 456         << endl << endl;
 457
 458     newCells.write();
 459
 460
 461
 462
 463     //
 464     // Invert cell split to construct map from new to old
 465     //
 466
 467     labelIOList newToOld
 468     (
 469         IOobject
 470         (
 471             "cellMap",
 472             runTime.timeName(),
 473             polyMesh::meshSubDir,
 474             mesh,
 475             IOobject::NO_READ,
 476             IOobject::AUTO_WRITE
 477         ),
 478         mesh.nCells()
 479     );
 480     newToOld.note() =
 481         "From cells in mesh at "
 482       + runTime.timeName()
 483       + " to cells in mesh at "
 484       + oldTimeName;
 485
 486
 487     forAll(oldToNew, oldCellI)
 488     {
 489         const labelList& added = oldToNew[oldCellI];
 490
 491         if (added.size())
 492         {
 493             forAll(added, i)
 494             {
 495                 newToOld[added[i]] = oldCellI;
 496             }
 497         }
 498         else
 499         {
 500             // Unrefined cell
 501             newToOld[oldCellI] = oldCellI;
 502         }
 503     }
 504
 505     Info<< "Writing map from new to old cell to "
 506         << newToOld.objectPath() << nl << endl;
 507
 508     newToOld.write();
 509
 510
 511     // Some statistics.
 512
 513     printEdgeStats(mesh);
 514
 515     Info<< "End\n" << endl;
 516
 517     return 0;
 518 }
 519
 520
 521 // ************************************************************************* //