applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.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 anispulation  |
   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 Application
  26     renumberMesh
  27
  28 Description
  29     Renumbers the cell list in order to reduce the bandwidth, reading and
  30     renumbering all fields from all the time directories.
  31
  32 \*---------------------------------------------------------------------------*/
  33
  34 #include "argList.H"
  35 #include "IOobjectList.H"
  36 #include "fvMesh.H"
  37 #include "polyTopoChange.H"
  38 #include "ReadFields.H"
  39 #include "volFields.H"
  40 #include "surfaceFields.H"
  41 #include "bandCompression.H"
  42 #include "faceSet.H"
  43 #include "SortableList.H"
  44 #include "decompositionMethod.H"
  45 #include "fvMeshSubset.H"
  46 #include "zeroGradientFvPatchFields.H"
  47
  48 using namespace Foam;
  49
  50
  51 // Calculate band of matrix
  52 label getBand(const labelList& owner, const labelList& neighbour)
  53 {
  54     label band = 0;
  55
  56     forAll(neighbour, faceI)
  57     {
  58         label diff = neighbour[faceI] - owner[faceI];
  59
  60         if (diff > band)
  61         {
  62             band = diff;
  63         }
  64     }
  65     return band;
  66 }
  67
  68
  69 // Return new to old cell numbering
  70 labelList regionBandCompression
  71 (
  72     const fvMesh& mesh,
  73     const labelList& cellToRegion
  74 )
  75 {
  76     Pout<< "Determining cell order:" << endl;
  77
  78     labelList cellOrder(cellToRegion.size());
  79
  80     label nRegions = max(cellToRegion)+1;
  81
  82     labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
  83
  84     label cellI = 0;
  85
  86     forAll(regionToCells, regionI)
  87     {
  88         Pout<< "    region " << regionI << " starts at " << cellI << endl;
  89
  90         // Per region do a reordering.
  91         fvMeshSubset subsetter(mesh);
  92         subsetter.setLargeCellSubset(cellToRegion, regionI);
  93         const fvMesh& subMesh = subsetter.subMesh();
  94         labelList subCellOrder(bandCompression(subMesh.cellCells()));
  95
  96         const labelList& cellMap = subsetter.cellMap();
  97
  98         forAll(subCellOrder, i)
  99         {
 100             cellOrder[cellI++] = cellMap[subCellOrder[i]];
 101         }
 102     }
 103     Pout<< endl;
 104
 105     return cellOrder;
 106 }
 107
 108
 109 // Determine face order such that inside region faces are sorted
 110 // upper-triangular but inbetween region faces are handled like boundary faces.
 111 labelList regionFaceOrder
 112 (
 113     const primitiveMesh& mesh,
 114     const labelList& cellOrder,     // new to old cell
 115     const labelList& cellToRegion   // old cell to region
 116 )
 117 {
 118     Pout<< "Determining face order:" << endl;
 119
 120     labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
 121
 122     labelList oldToNewFace(mesh.nFaces(), -1);
 123
 124     label newFaceI = 0;
 125
 126     label prevRegion = -1;
 127
 128     forAll (cellOrder, newCellI)
 129     {
 130         label oldCellI = cellOrder[newCellI];
 131
 132         if (cellToRegion[oldCellI] != prevRegion)
 133         {
 134             prevRegion = cellToRegion[oldCellI];
 135             Pout<< "    region " << prevRegion << " internal faces start at "
 136                 << newFaceI << endl;
 137         }
 138
 139         const cell& cFaces = mesh.cells()[oldCellI];
 140
 141         SortableList<label> nbr(cFaces.size());
 142
 143         forAll(cFaces, i)
 144         {
 145             label faceI = cFaces[i];
 146
 147             if (mesh.isInternalFace(faceI))
 148             {
 149                 // Internal face. Get cell on other side.
 150                 label nbrCellI = reverseCellOrder[mesh.faceNeighbour()[faceI]];
 151                 if (nbrCellI == newCellI)
 152                 {
 153                     nbrCellI = reverseCellOrder[mesh.faceOwner()[faceI]];
 154                 }
 155
 156                 if (cellToRegion[oldCellI] != cellToRegion[cellOrder[nbrCellI]])
 157                 {
 158                     // Treat like external face. Do later.
 159                     nbr[i] = -1;
 160                 }
 161                 else if (newCellI < nbrCellI)
 162                 {
 163                     // CellI is master
 164                     nbr[i] = nbrCellI;
 165                 }
 166                 else
 167                 {
 168                     // nbrCell is master. Let it handle this face.
 169                     nbr[i] = -1;
 170                 }
 171             }
 172             else
 173             {
 174                 // External face. Do later.
 175                 nbr[i] = -1;
 176             }
 177         }
 178
 179         nbr.sort();
 180
 181         forAll(nbr, i)
 182         {
 183             if (nbr[i] != -1)
 184             {
 185                 oldToNewFace[cFaces[nbr.indices()[i]]] = newFaceI++;
 186             }
 187         }
 188     }
 189
 190     // Do region interfaces
 191     label nRegions = max(cellToRegion)+1;
 192     {
 193         // Sort in increasing region
 194         SortableList<label> sortKey(mesh.nFaces(), labelMax);
 195
 196         for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
 197         {
 198             label ownRegion = cellToRegion[mesh.faceOwner()[faceI]];
 199             label neiRegion = cellToRegion[mesh.faceNeighbour()[faceI]];
 200
 201             if (ownRegion != neiRegion)
 202             {
 203                 sortKey[faceI] =
 204                     min(ownRegion, neiRegion)*nRegions
 205                    +max(ownRegion, neiRegion);
 206             }
 207         }
 208         sortKey.sort();
 209
 210         // Extract.
 211         label prevKey = -1;
 212         forAll(sortKey, i)
 213         {
 214             label key = sortKey[i];
 215
 216             if (key == labelMax)
 217             {
 218                 break;
 219             }
 220
 221             if (prevKey != key)
 222             {
 223                 Pout<< "    faces inbetween region " << key/nRegions
 224                     << " and " << key%nRegions
 225                     << " start at " << newFaceI << endl;
 226                 prevKey = key;
 227             }
 228
 229             oldToNewFace[sortKey.indices()[i]] = newFaceI++;
 230         }
 231     }
 232
 233     // Leave patch faces intact.
 234     for (label faceI = newFaceI; faceI < mesh.nFaces(); faceI++)
 235     {
 236         oldToNewFace[faceI] = faceI;
 237     }
 238
 239
 240     // Check done all faces.
 241     forAll(oldToNewFace, faceI)
 242     {
 243         if (oldToNewFace[faceI] == -1)
 244         {
 245             FatalErrorIn
 246             (
 247                 "polyDualMesh::getFaceOrder"
 248                 "(const labelList&, const labelList&, const label) const"
 249             )   << "Did not determine new position"
 250                 << " for face " << faceI
 251                 << abort(FatalError);
 252         }
 253     }
 254     Pout<< endl;
 255
 256     return invert(mesh.nFaces(), oldToNewFace);
 257 }
 258
 259
 260 // cellOrder: old cell for every new cell
 261 // faceOrder: old face for every new face. Ordering of boundary faces not
 262 // changed.
 263 autoPtr<mapPolyMesh> reorderMesh
 264 (
 265     polyMesh& mesh,
 266     const labelList& cellOrder,
 267     const labelList& faceOrder
 268 )
 269 {
 270     labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
 271     labelList reverseFaceOrder(invert(faceOrder.size(), faceOrder));
 272
 273     faceList newFaces(reorder(reverseFaceOrder, mesh.faces()));
 274     labelList newOwner
 275     (
 276         renumber
 277         (
 278             reverseCellOrder,
 279             reorder(reverseFaceOrder, mesh.faceOwner())
 280         )
 281     );
 282     labelList newNeighbour
 283     (
 284         renumber
 285         (
 286             reverseCellOrder,
 287             reorder(reverseFaceOrder, mesh.faceNeighbour())
 288         )
 289     );
 290
 291     // Check if any faces need swapping.
 292     forAll(newNeighbour, faceI)
 293     {
 294         label own = newOwner[faceI];
 295         label nei = newNeighbour[faceI];
 296
 297         if (nei < own)
 298         {
 299             newFaces[faceI] = newFaces[faceI].reverseFace();
 300             Swap(newOwner[faceI], newNeighbour[faceI]);
 301         }
 302     }
 303
 304     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 305     labelList patchSizes(patches.size());
 306     labelList patchStarts(patches.size());
 307     labelList oldPatchNMeshPoints(patches.size());
 308     labelListList patchPointMap(patches.size());
 309     forAll(patches, patchI)
 310     {
 311         patchSizes[patchI] = patches[patchI].size();
 312         patchStarts[patchI] = patches[patchI].start();
 313         oldPatchNMeshPoints[patchI] = patches[patchI].nPoints();
 314         patchPointMap[patchI] = identity(patches[patchI].nPoints());
 315     }
 316
 317     mesh.resetPrimitives
 318     (
 319         mesh.nFaces(),
 320         mesh.points(),
 321         newFaces,
 322         newOwner,
 323         newNeighbour,
 324         patchSizes,
 325         patchStarts
 326     );
 327
 328     return autoPtr<mapPolyMesh>
 329     (
 330         new mapPolyMesh
 331         (
 332             mesh,                       //const polyMesh& mesh,
 333             mesh.nPoints(),             // nOldPoints,
 334             mesh.nFaces(),              // nOldFaces,
 335             mesh.nCells(),              // nOldCells,
 336             identity(mesh.nPoints()),   // pointMap,
 337             List<objectMap>(0),         // pointsFromPoints,
 338             faceOrder,                  // faceMap,
 339             List<objectMap>(0),         // facesFromPoints,
 340             List<objectMap>(0),         // facesFromEdges,
 341             List<objectMap>(0),         // facesFromFaces,
 342             cellOrder,                  // cellMap,
 343             List<objectMap>(0),         // cellsFromPoints,
 344             List<objectMap>(0),         // cellsFromEdges,
 345             List<objectMap>(0),         // cellsFromFaces,
 346             List<objectMap>(0),         // cellsFromCells,
 347             identity(mesh.nPoints()),   // reversePointMap,
 348             reverseFaceOrder,           // reverseFaceMap,
 349             reverseCellOrder,           // reverseCellMap,
 350             labelHashSet(0),            // flipFaceFlux,
 351             patchPointMap,              // patchPointMap,
 352             labelListList(0),           // pointZoneMap,
 353             labelListList(0),           // faceZonePointMap,
 354             labelListList(0),           // faceZoneFaceMap,
 355             labelListList(0),           // cellZoneMap,
 356             pointField(0),              // preMotionPoints,
 357             patchStarts,                // oldPatchStarts,
 358             oldPatchNMeshPoints         // oldPatchNMeshPoints
 359         )
 360     );
 361 }
 362
 363
 364 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 365
 366 int main(int argc, char *argv[])
 367 {
 368     argList::validOptions.insert("blockOrder", "");
 369     argList::validOptions.insert("orderPoints", "");
 370     argList::validOptions.insert("writeMaps", "");
 371     argList::validOptions.insert("overwrite", "");
 372
 373 #   include "addTimeOptions.H"
 374
 375 #   include "setRootCase.H"
 376 #   include "createTime.H"
 377
 378     // Get times list
 379     instantList Times = runTime.times();
 380
 381     // set startTime and endTime depending on -time and -latestTime options
 382 #   include "checkTimeOptions.H"
 383
 384     runTime.setTime(Times[startTime], startTime);
 385
 386
 387 #   include "createMesh.H"
 388
 389
 390     const bool blockOrder = args.options().found("blockOrder");
 391
 392     if (blockOrder)
 393     {
 394         Info<< "Ordering cells into regions (using decomposition);"
 395             << " ordering faces into region-internal and region-external." << nl
 396             << endl;
 397     }
 398
 399     const bool orderPoints = args.options().found("orderPoints");
 400
 401     if (orderPoints)
 402     {
 403         Info<< "Ordering points into internal and boundary points." << nl
 404             << endl;
 405     }
 406
 407     const bool writeMaps = args.options().found("writeMaps");
 408
 409     if (writeMaps)
 410     {
 411         Info<< "Writing renumber maps (new to old) to polyMesh." << nl
 412             << endl;
 413     }
 414
 415     bool overwrite = args.options().found("overwrite");
 416
 417     label band = getBand(mesh.faceOwner(), mesh.faceNeighbour());
 418
 419     Info<< "Mesh size: " << returnReduce(mesh.nCells(), sumOp<label>()) << nl
 420         << "Band before renumbering: "
 421         << returnReduce(band, maxOp<label>()) << nl << endl;
 422
 423     // Read objects in time directory
 424     IOobjectList objects(mesh, runTime.timeName());
 425
 426     // Read vol fields.
 427
 428     PtrList<volScalarField> vsFlds;
 429     ReadFields(mesh, objects, vsFlds);
 430
 431     PtrList<volVectorField> vvFlds;
 432     ReadFields(mesh, objects, vvFlds);
 433
 434     PtrList<volSphericalTensorField> vstFlds;
 435     ReadFields(mesh, objects, vstFlds);
 436
 437     PtrList<volSymmTensorField> vsymtFlds;
 438     ReadFields(mesh, objects, vsymtFlds);
 439
 440     PtrList<volTensorField> vtFlds;
 441     ReadFields(mesh, objects, vtFlds);
 442
 443     // Read surface fields.
 444
 445     PtrList<surfaceScalarField> ssFlds;
 446     ReadFields(mesh, objects, ssFlds);
 447
 448     PtrList<surfaceVectorField> svFlds;
 449     ReadFields(mesh, objects, svFlds);
 450
 451     PtrList<surfaceSphericalTensorField> sstFlds;
 452     ReadFields(mesh, objects, sstFlds);
 453
 454     PtrList<surfaceSymmTensorField> ssymtFlds;
 455     ReadFields(mesh, objects, ssymtFlds);
 456
 457     PtrList<surfaceTensorField> stFlds;
 458     ReadFields(mesh, objects, stFlds);
 459
 460
 461     autoPtr<mapPolyMesh> map;
 462
 463     if (blockOrder)
 464     {
 465         // Renumbering in two phases. Should be done in one so mapping of
 466         // fields is done correctly!
 467
 468         // Read decomposePar dictionary
 469         IOdictionary decomposeDict
 470         (
 471             IOobject
 472             (
 473                 "decomposeParDict",
 474                 runTime.system(),
 475                 mesh,
 476                 IOobject::MUST_READ,
 477                 IOobject::NO_WRITE
 478             )
 479         );
 480         autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
 481         (
 482             decomposeDict,
 483             mesh
 484         );
 485
 486         labelList cellToRegion(decomposePtr().decompose(mesh.cellCentres()));
 487
 488         // For debugging: write out region
 489         {
 490             volScalarField cellDist
 491             (
 492                 IOobject
 493                 (
 494                     "cellDist",
 495                     runTime.timeName(),
 496                     mesh,
 497                     IOobject::NO_READ,
 498                     IOobject::NO_WRITE,
 499                     false
 500                 ),
 501                 mesh,
 502                 dimensionedScalar("cellDist", dimless, 0),
 503                 zeroGradientFvPatchScalarField::typeName
 504             );
 505
 506             forAll(cellToRegion, cellI)
 507             {
 508                cellDist[cellI] = cellToRegion[cellI];
 509             }
 510
 511             cellDist.write();
 512
 513             Info<< nl << "Written decomposition as volScalarField to "
 514                 << cellDist.name() << " for use in postprocessing."
 515                 << nl << endl;
 516         }
 517
 518         // Use block based renumbering.
 519         //labelList cellOrder(bandCompression(mesh.cellCells()));
 520         labelList cellOrder(regionBandCompression(mesh, cellToRegion));
 521
 522         // Determine new to old face order with new cell numbering
 523         labelList faceOrder
 524         (
 525             regionFaceOrder
 526             (
 527                 mesh,
 528                 cellOrder,
 529                 cellToRegion
 530             )
 531         );
 532
 533         if (!overwrite)
 534         {
 535             runTime++;
 536         }
 537
 538         // Change the mesh.
 539         map = reorderMesh(mesh, cellOrder, faceOrder);
 540     }
 541     else
 542     {
 543         // Use built-in renumbering.
 544
 545         polyTopoChange meshMod(mesh);
 546
 547         if (!overwrite)
 548         {
 549             runTime++;
 550         }
 551
 552         // Change the mesh.
 553         map = meshMod.changeMesh
 554         (
 555             mesh,
 556             false,      // inflate
 557             true,       // parallel sync
 558             true,       // cell ordering
 559             orderPoints // point ordering
 560         );
 561     }
 562
 563     // Update fields
 564     mesh.updateMesh(map);
 565
 566     // Move mesh (since morphing might not do this)
 567     if (map().hasMotionPoints())
 568     {
 569         mesh.movePoints(map().preMotionPoints());
 570     }
 571
 572
 573     band = getBand(mesh.faceOwner(), mesh.faceNeighbour());
 574
 575     Info<< "Band after renumbering: "
 576         << returnReduce(band, maxOp<label>()) << nl << endl;
 577
 578
 579     if (orderPoints)
 580     {
 581         // Force edge calculation (since only reason that points would need to
 582         // be sorted)
 583         (void)mesh.edges();
 584
 585         label nTotPoints = returnReduce
 586         (
 587             mesh.nPoints(),
 588             sumOp<label>()
 589         );
 590         label nTotIntPoints = returnReduce
 591         (
 592             mesh.nInternalPoints(),
 593             sumOp<label>()
 594         );
 595
 596         label nTotEdges = returnReduce
 597         (
 598             mesh.nEdges(),
 599             sumOp<label>()
 600         );
 601         label nTotIntEdges = returnReduce
 602         (
 603             mesh.nInternalEdges(),
 604             sumOp<label>()
 605         );
 606         label nTotInt0Edges = returnReduce
 607         (
 608             mesh.nInternal0Edges(),
 609             sumOp<label>()
 610         );
 611         label nTotInt1Edges = returnReduce
 612         (
 613             mesh.nInternal1Edges(),
 614             sumOp<label>()
 615         );
 616
 617         Info<< "Points:" << nl
 618             << "    total   : " << nTotPoints << nl
 619             << "    internal: " << nTotIntPoints << nl
 620             << "    boundary: " << nTotPoints-nTotIntPoints << nl
 621             << "Edges:" << nl
 622             << "    total   : " << nTotEdges << nl
 623             << "    internal: " << nTotIntEdges << nl
 624             << "        internal using 0 boundary points: "
 625             << nTotInt0Edges << nl
 626             << "        internal using 1 boundary points: "
 627             << nTotInt1Edges-nTotInt0Edges << nl
 628             << "        internal using 2 boundary points: "
 629             << nTotIntEdges-nTotInt1Edges << nl
 630             << "    boundary: " << nTotEdges-nTotIntEdges << nl
 631             << endl;
 632     }
 633
 634     Info<< "Writing mesh to " << runTime.timeName() << endl;
 635
 636     mesh.write();
 637
 638     if (writeMaps)
 639     {
 640         labelIOList
 641         (
 642             IOobject
 643             (
 644                 "cellMap",
 645                 mesh.facesInstance(),
 646                 polyMesh::meshSubDir,
 647                 mesh,
 648                 IOobject::NO_READ,
 649                 IOobject::NO_WRITE,
 650                 false
 651             ),
 652             map().cellMap()
 653         ).write();
 654         labelIOList
 655         (
 656             IOobject
 657             (
 658                 "faceMap",
 659                 mesh.facesInstance(),
 660                 polyMesh::meshSubDir,
 661                 mesh,
 662                 IOobject::NO_READ,
 663                 IOobject::NO_WRITE,
 664                 false
 665             ),
 666             map().faceMap()
 667         ).write();
 668         labelIOList
 669         (
 670             IOobject
 671             (
 672                 "pointMap",
 673                 mesh.facesInstance(),
 674                 polyMesh::meshSubDir,
 675                 mesh,
 676                 IOobject::NO_READ,
 677                 IOobject::NO_WRITE,
 678                 false
 679             ),
 680             map().pointMap()
 681         ).write();
 682     }
 683
 684     Info<< "\nEnd.\n" << endl;
 685
 686     return 0;
 687 }
 688
 689
 690 // ************************************************************************* //