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) 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     Utility to refine cells in multiple directions.
  26
  27     Either supply -all option to refine all cells (3D refinement for 3D
  28     cases; 2D for 2D cases) or reads a refineMeshDict with
  29     - cellSet to refine
  30     - directions to refine
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "argList.H"
  35 #include "polyMesh.H"
  36 #include "Time.H"
  37 #include "undoableMeshCutter.H"
  38 #include "hexCellLooper.H"
  39 #include "cellSet.H"
  40 #include "twoDPointCorrector.H"
  41 #include "directions.H"
  42 #include "OFstream.H"
  43 #include "multiDirRefinement.H"
  44 #include "labelIOList.H"
  45 #include "wedgePolyPatch.H"
  46 #include "plane.H"
  47 #include "SubField.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             const 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     argList::addNote
 295     (
 296         "refine cells in multiple directions"
 297     );
 298
 299     #include "addOverwriteOption.H"
 300     #include "addRegionOption.H"
 301     argList::addBoolOption
 302     (
 303         "dict",
 304         "refine according to system/refineMeshDict"
 305     );
 306
 307     #include "setRootCase.H"
 308     #include "createTime.H"
 309     runTime.functionObjects().off();
 310     #include "createNamedPolyMesh.H"
 311     const word oldInstance = mesh.pointsInstance();
 312
 313     printEdgeStats(mesh);
 314
 315     //
 316     // Read/construct control dictionary
 317     //
 318
 319     const bool readDict = args.optionFound("dict");
 320     const bool overwrite = args.optionFound("overwrite");
 321
 322     // List of cells to refine
 323     labelList refCells;
 324
 325     // Dictionary to control refinement
 326     dictionary refineDict;
 327
 328     if (readDict)
 329     {
 330         Info<< "Refining according to refineMeshDict" << nl << endl;
 331
 332         refineDict = IOdictionary
 333         (
 334             IOobject
 335             (
 336                 "refineMeshDict",
 337                 runTime.system(),
 338                 mesh,
 339                 IOobject::MUST_READ_IF_MODIFIED,
 340                 IOobject::NO_WRITE
 341             )
 342         );
 343
 344         const word setName(refineDict.lookup("set"));
 345
 346         cellSet cells(mesh, setName);
 347
 348         Pout<< "Read " << cells.size() << " cells from cellSet "
 349             << cells.instance()/cells.local()/cells.name()
 350             << endl << endl;
 351
 352         refCells = cells.toc();
 353     }
 354     else
 355     {
 356         Info<< "Refining all cells" << nl << endl;
 357
 358         // Select all cells
 359         refCells.setSize(mesh.nCells());
 360
 361         forAll(mesh.cells(), cellI)
 362         {
 363             refCells[cellI] = cellI;
 364         }
 365
 366
 367         // Set refinement directions based on 2D/3D
 368         label axisIndex = twoDNess(mesh);
 369
 370         if (axisIndex == -1)
 371         {
 372             Info<< "3D case; refining all directions" << nl << endl;
 373
 374             wordList directions(3);
 375             directions[0] = "tan1";
 376             directions[1] = "tan2";
 377             directions[2] = "normal";
 378             refineDict.add("directions", directions);
 379
 380             // Use hex cutter
 381             refineDict.add("useHexTopology", "true");
 382         }
 383         else
 384         {
 385             wordList directions(2);
 386
 387             if (axisIndex == 0)
 388             {
 389                 Info<< "2D case; refining in directions y,z\n" << endl;
 390                 directions[0] = "tan2";
 391                 directions[1] = "normal";
 392             }
 393             else if (axisIndex == 1)
 394             {
 395                 Info<< "2D case; refining in directions x,z\n" << endl;
 396                 directions[0] = "tan1";
 397                 directions[1] = "normal";
 398             }
 399             else
 400             {
 401                 Info<< "2D case; refining in directions x,y\n" << endl;
 402                 directions[0] = "tan1";
 403                 directions[1] = "tan2";
 404             }
 405
 406             refineDict.add("directions", directions);
 407
 408             // Use standard cutter
 409             refineDict.add("useHexTopology", "false");
 410         }
 411
 412         refineDict.add("coordinateSystem", "global");
 413
 414         dictionary coeffsDict;
 415         coeffsDict.add("tan1", vector(1, 0, 0));
 416         coeffsDict.add("tan2", vector(0, 1, 0));
 417         refineDict.add("globalCoeffs", coeffsDict);
 418
 419         refineDict.add("geometricCut", "false");
 420         refineDict.add("writeMesh", "false");
 421     }
 422
 423
 424     string oldTimeName(runTime.timeName());
 425
 426     if (!overwrite)
 427     {
 428         runTime++;
 429     }
 430
 431
 432     // Multi-directional refinement (does multiple iterations)
 433     multiDirRefinement multiRef(mesh, refCells, refineDict);
 434
 435
 436     // Write resulting mesh
 437     if (overwrite)
 438     {
 439         mesh.setInstance(oldInstance);
 440     }
 441     mesh.write();
 442
 443
 444     // Get list of cell splits.
 445     // (is for every cell in old mesh the cells they have been split into)
 446     const labelListList& oldToNew = multiRef.addedCells();
 447
 448
 449     // Create cellSet with added cells for easy inspection
 450     cellSet newCells(mesh, "refinedCells", refCells.size());
 451
 452     forAll(oldToNew, oldCellI)
 453     {
 454         const labelList& added = oldToNew[oldCellI];
 455
 456         forAll(added, i)
 457         {
 458             newCells.insert(added[i]);
 459         }
 460     }
 461
 462     Pout<< "Writing refined cells (" << newCells.size() << ") to cellSet "
 463         << newCells.instance()/newCells.local()/newCells.name()
 464         << endl << endl;
 465
 466     newCells.write();
 467
 468
 469
 470
 471     //
 472     // Invert cell split to construct map from new to old
 473     //
 474
 475     labelIOList newToOld
 476     (
 477         IOobject
 478         (
 479             "cellMap",
 480             runTime.timeName(),
 481             polyMesh::meshSubDir,
 482             mesh,
 483             IOobject::NO_READ,
 484             IOobject::AUTO_WRITE
 485         ),
 486         mesh.nCells()
 487     );
 488     newToOld.note() =
 489         "From cells in mesh at "
 490       + runTime.timeName()
 491       + " to cells in mesh at "
 492       + oldTimeName;
 493
 494
 495     forAll(oldToNew, oldCellI)
 496     {
 497         const labelList& added = oldToNew[oldCellI];
 498
 499         if (added.size())
 500         {
 501             forAll(added, i)
 502             {
 503                 newToOld[added[i]] = oldCellI;
 504             }
 505         }
 506         else
 507         {
 508             // Unrefined cell
 509             newToOld[oldCellI] = oldCellI;
 510         }
 511     }
 512
 513     Info<< "Writing map from new to old cell to "
 514         << newToOld.objectPath() << nl << endl;
 515
 516     newToOld.write();
 517
 518
 519     // Some statistics.
 520
 521     printEdgeStats(mesh);
 522
 523     Info<< "End\n" << endl;
 524
 525     return 0;
 526 }
 527
 528
 529 // ************************************************************************* //