applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.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 Application
  25     reconstructParMesh
  26
  27 Description
  28     Reconstructs a mesh using geometric information only.
  29
  30     Writes point/face/cell procAddressing so afterwards reconstructPar can be
  31     used to reconstruct fields.
  32
  33     Note:
  34     - uses geometric matching tolerance (set with -mergeTol (at your option)
  35
  36     If the parallel case does not have correct procBoundaries use the
  37     -fullMatch option which will check all boundary faces (bit slower).
  38
  39 \*---------------------------------------------------------------------------*/
  40
  41 #include "argList.H"
  42 #include "Time.H"
  43 #include "IOobjectList.H"
  44 #include "labelIOList.H"
  45 #include "processorPolyPatch.H"
  46 #include "mapAddedPolyMesh.H"
  47 #include "polyMeshAdder.H"
  48 #include "faceCoupleInfo.H"
  49 #include "fvMeshAdder.H"
  50 #include "polyTopoChange.H"
  51
  52 using namespace Foam;
  53
  54 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  55
  56 // Tolerance (as fraction of the bounding box). Needs to be fairly lax since
  57 // usually meshes get written with limited precision (6 digits)
  58 static const scalar defaultMergeTol = 1E-7;
  59
  60
  61 static void renumber
  62 (
  63     const labelList& map,
  64     labelList& elems
  65 )
  66 {
  67     forAll(elems, i)
  68     {
  69         if (elems[i] >= 0)
  70         {
  71             elems[i] = map[elems[i]];
  72         }
  73     }
  74 }
  75
  76
  77 // Determine which faces are coupled. Uses geometric merge distance.
  78 // Looks either at all boundaryFaces (fullMatch) or only at the
  79 // procBoundaries for procI. Assumes that masterMesh contains already merged
  80 // all the processors < procI.
  81 autoPtr<faceCoupleInfo> determineCoupledFaces
  82 (
  83     const bool fullMatch,
  84     const label procI,
  85     const polyMesh& masterMesh,
  86     const polyMesh& meshToAdd,
  87     const scalar mergeDist
  88 )
  89 {
  90     if (fullMatch || masterMesh.nCells() == 0)
  91     {
  92         return autoPtr<faceCoupleInfo>
  93         (
  94             new faceCoupleInfo
  95             (
  96                 masterMesh,
  97                 meshToAdd,
  98                 mergeDist,      // absolute merging distance
  99                 true            // matching faces identical
 100             )
 101         );
 102     }
 103     else
 104     {
 105         // Pick up all patches on masterMesh ending in "toDDD" where DDD is
 106         // the processor number procI.
 107
 108         const polyBoundaryMesh& masterPatches = masterMesh.boundaryMesh();
 109
 110         const string toProcString("to" + name(procI));
 111
 112         DynamicList<label> masterFaces
 113         (
 114             masterMesh.nFaces()
 115           - masterMesh.nInternalFaces()
 116         );
 117
 118         forAll(masterPatches, patchI)
 119         {
 120             const polyPatch& pp = masterPatches[patchI];
 121
 122             if
 123             (
 124                 isA<processorPolyPatch>(pp)
 125              && (
 126                     pp.name().rfind(toProcString)
 127                  == (pp.name().size()-toProcString.size())
 128                 )
 129             )
 130             {
 131                 label meshFaceI = pp.start();
 132                 forAll(pp, i)
 133                 {
 134                     masterFaces.append(meshFaceI++);
 135                 }
 136             }
 137         }
 138         masterFaces.shrink();
 139
 140
 141         // Pick up all patches on meshToAdd ending in "procBoundaryDDDtoYYY"
 142         // where DDD is the processor number procI and YYY is < procI.
 143
 144         const polyBoundaryMesh& addPatches = meshToAdd.boundaryMesh();
 145
 146         DynamicList<label> addFaces
 147         (
 148             meshToAdd.nFaces()
 149           - meshToAdd.nInternalFaces()
 150         );
 151
 152         forAll(addPatches, patchI)
 153         {
 154             const polyPatch& pp = addPatches[patchI];
 155
 156             if (isA<processorPolyPatch>(pp))
 157             {
 158                 bool isConnected = false;
 159
 160                 for (label mergedProcI = 0; mergedProcI < procI; mergedProcI++)
 161                 {
 162                     const string fromProcString
 163                     (
 164                         "procBoundary"
 165                       + name(procI)
 166                       + "to"
 167                       + name(mergedProcI)
 168                     );
 169
 170                     if (pp.name() == fromProcString)
 171                     {
 172                         isConnected = true;
 173                         break;
 174                     }
 175                 }
 176
 177                 if (isConnected)
 178                 {
 179                     label meshFaceI = pp.start();
 180                     forAll(pp, i)
 181                     {
 182                         addFaces.append(meshFaceI++);
 183                     }
 184                 }
 185             }
 186         }
 187         addFaces.shrink();
 188
 189         return autoPtr<faceCoupleInfo>
 190         (
 191             new faceCoupleInfo
 192             (
 193                 masterMesh,
 194                 masterFaces,
 195                 meshToAdd,
 196                 addFaces,
 197                 mergeDist,      // absolute merging distance
 198                 true,           // matching faces identical?
 199                 false,          // if perfectmatch are faces already ordered
 200                                 // (e.g. processor patches)
 201                 false           // are faces each on separate patch?
 202             )
 203         );
 204     }
 205 }
 206
 207
 208 autoPtr<mapPolyMesh> mergeSharedPoints
 209 (
 210     const scalar mergeDist,
 211     polyMesh& mesh,
 212     labelListList& pointProcAddressing
 213 )
 214 {
 215     // Find out which sets of points get merged and create a map from
 216     // mesh point to unique point.
 217     Map<label> pointToMaster
 218     (
 219         fvMeshAdder::findSharedPoints
 220         (
 221             mesh,
 222             mergeDist
 223         )
 224     );
 225
 226     Info<< "mergeSharedPoints : detected " << pointToMaster.size()
 227         << " points that are to be merged." << endl;
 228
 229     if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
 230     {
 231         return autoPtr<mapPolyMesh>(NULL);
 232     }
 233
 234     polyTopoChange meshMod(mesh);
 235
 236     fvMeshAdder::mergePoints(mesh, pointToMaster, meshMod);
 237
 238     // Change the mesh (no inflation). Note: parallel comms allowed.
 239     autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
 240
 241     // Update fields. No inflation, parallel sync.
 242     mesh.updateMesh(map);
 243
 244     // pointProcAddressing give indices into the master mesh so adapt them
 245     // for changed point numbering.
 246
 247     // Adapt constructMaps for merged points.
 248     forAll(pointProcAddressing, procI)
 249     {
 250         labelList& constructMap = pointProcAddressing[procI];
 251
 252         forAll(constructMap, i)
 253         {
 254             label oldPointI = constructMap[i];
 255
 256             // New label of point after changeMesh.
 257             label newPointI = map().reversePointMap()[oldPointI];
 258
 259             if (newPointI < -1)
 260             {
 261                 constructMap[i] = -newPointI-2;
 262             }
 263             else if (newPointI >= 0)
 264             {
 265                 constructMap[i] = newPointI;
 266             }
 267             else
 268             {
 269                 FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
 270                     << "Problem. oldPointI:" << oldPointI
 271                     << " newPointI:" << newPointI << abort(FatalError);
 272             }
 273         }
 274     }
 275
 276     return map;
 277 }
 278
 279
 280 int main(int argc, char *argv[])
 281 {
 282     argList::addNote
 283     (
 284         "reconstruct a mesh using geometric information only"
 285     );
 286
 287     argList::noParallel();
 288     argList::addOption
 289     (
 290         "mergeTol",
 291         "scalar",
 292         "specify the merge distance relative to the bounding box size "
 293         "(default 1E-7)"
 294     );
 295     argList::addBoolOption
 296     (
 297         "fullMatch",
 298         "do (slower) geometric matching on all boundary faces"
 299     );
 300
 301     #include "addTimeOptions.H"
 302     #include "addRegionOption.H"
 303     #include "setRootCase.H"
 304     #include "createTime.H"
 305
 306     Info<< "This is an experimental tool which tries to merge"
 307         << " individual processor" << nl
 308         << "meshes back into one master mesh. Use it if the original"
 309         << " master mesh has" << nl
 310         << "been deleted or if the processor meshes have been modified"
 311         << " (topology change)." << nl
 312         << "This tool will write the resulting mesh to a new time step"
 313         << " and construct" << nl
 314         << "xxxxProcAddressing files in the processor meshes so"
 315         << " reconstructPar can be" << nl
 316         << "used to regenerate the fields on the master mesh." << nl
 317         << nl
 318         << "Not well tested & use at your own risk!" << nl
 319         << endl;
 320
 321
 322     word regionName = polyMesh::defaultRegion;
 323     word regionDir = word::null;
 324
 325     if (args.optionReadIfPresent("region", regionName))
 326     {
 327         regionDir = regionName;
 328         Info<< "Operating on region " << regionName << nl << endl;
 329     }
 330
 331     scalar mergeTol = defaultMergeTol;
 332     args.optionReadIfPresent("mergeTol", mergeTol);
 333
 334     scalar writeTol = Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
 335
 336     Info<< "Merge tolerance : " << mergeTol << nl
 337         << "Write tolerance : " << writeTol << endl;
 338
 339     if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
 340     {
 341         FatalErrorIn(args.executable())
 342             << "Your current settings specify ASCII writing with "
 343             << IOstream::defaultPrecision() << " digits precision." << endl
 344             << "Your merging tolerance (" << mergeTol << ") is finer than this."
 345             << endl
 346             << "Please change your writeFormat to binary"
 347             << " or increase the writePrecision" << endl
 348             << "or adjust the merge tolerance (-mergeTol)."
 349             << exit(FatalError);
 350     }
 351
 352
 353     const bool fullMatch = args.optionFound("fullMatch");
 354
 355     if (fullMatch)
 356     {
 357         Info<< "Doing geometric matching on all boundary faces." << nl << endl;
 358     }
 359     else
 360     {
 361         Info<< "Doing geometric matching on correct procBoundaries only."
 362             << nl << "This assumes a correct decomposition." << endl;
 363     }
 364
 365
 366
 367     int nProcs = 0;
 368
 369     while
 370     (
 371         isDir
 372         (
 373             args.rootPath()
 374           / args.caseName()
 375           / fileName(word("processor") + name(nProcs))
 376         )
 377     )
 378     {
 379         nProcs++;
 380     }
 381
 382     Info<< "Found " << nProcs << " processor directories" << nl << endl;
 383
 384
 385     // Read all databases.
 386     PtrList<Time> databases(nProcs);
 387
 388     forAll(databases, procI)
 389     {
 390         Info<< "Reading database "
 391             << args.caseName()/fileName(word("processor") + name(procI))
 392             << endl;
 393
 394         databases.set
 395         (
 396             procI,
 397             new Time
 398             (
 399                 Time::controlDictName,
 400                 args.rootPath(),
 401                 args.caseName()/fileName(word("processor") + name(procI))
 402             )
 403         );
 404
 405         Time& procTime = databases[procI];
 406
 407         instantList Times = procTime.times();
 408
 409         // set startTime and endTime depending on -time and -latestTime options
 410 #       include "checkTimeOptions.H"
 411
 412         procTime.setTime(Times[startTime], startTime);
 413
 414         if (procI > 0 && databases[procI-1].value() != procTime.value())
 415         {
 416             FatalErrorIn(args.executable())
 417                 << "Time not equal on processors." << nl
 418                 << "Processor:" << procI-1
 419                 << " time:" << databases[procI-1].value() << nl
 420                 << "Processor:" << procI
 421                 << " time:" << procTime.value()
 422                 << exit(FatalError);
 423         }
 424     }
 425
 426     // Set master time
 427     Info<< "Setting master time to " << databases[0].timeName() << nl << endl;
 428     runTime.setTime(databases[0]);
 429
 430
 431     // Read point on individual processors to determine merge tolerance
 432     // (otherwise single cell domains might give problems)
 433
 434     boundBox bb = boundBox::invertedBox;
 435
 436     for (label procI = 0; procI < nProcs; procI++)
 437     {
 438         fileName pointsInstance
 439         (
 440             databases[procI].findInstance
 441             (
 442                 regionDir/polyMesh::meshSubDir,
 443                 "points"
 444             )
 445         );
 446
 447         if (pointsInstance != databases[procI].timeName())
 448         {
 449             FatalErrorIn(args.executable())
 450                 << "Your time was specified as " << databases[procI].timeName()
 451                 << " but there is no polyMesh/points in that time." << endl
 452                 << "(there is a points file in " << pointsInstance
 453                 << ")" << endl
 454                 << "Please rerun with the correct time specified"
 455                 << " (through the -constant, -time or -latestTime "
 456                 << "(at your option)."
 457                 << endl << exit(FatalError);
 458         }
 459
 460         Info<< "Reading points from "
 461             << databases[procI].caseName()
 462             << " for time = " << databases[procI].timeName()
 463             << nl << endl;
 464
 465         pointIOField points
 466         (
 467             IOobject
 468             (
 469                 "points",
 470                 databases[procI].findInstance
 471                 (
 472                     regionDir/polyMesh::meshSubDir,
 473                     "points"
 474                 ),
 475                 regionDir/polyMesh::meshSubDir,
 476                 databases[procI],
 477                 IOobject::MUST_READ,
 478                 IOobject::NO_WRITE,
 479                 false
 480             )
 481         );
 482
 483         boundBox domainBb(points, false);
 484
 485         bb.min() = min(bb.min(), domainBb.min());
 486         bb.max() = max(bb.max(), domainBb.max());
 487     }
 488     const scalar mergeDist = mergeTol * bb.mag();
 489
 490     Info<< "Overall mesh bounding box  : " << bb << nl
 491         << "Relative tolerance         : " << mergeTol << nl
 492         << "Absolute matching distance : " << mergeDist << nl
 493         << endl;
 494
 495
 496     // Addressing from processor to reconstructed case
 497     labelListList cellProcAddressing(nProcs);
 498     labelListList faceProcAddressing(nProcs);
 499     labelListList pointProcAddressing(nProcs);
 500     labelListList boundaryProcAddressing(nProcs);
 501
 502     // Internal faces on the final reconstructed mesh
 503     label masterInternalFaces;
 504     // Owner addressing on the final reconstructed mesh
 505     labelList masterOwner;
 506
 507     {
 508         // Construct empty mesh.
 509         Info<< "Constructing empty mesh to add to." << nl << endl;
 510         polyMesh masterMesh
 511         (
 512             IOobject
 513             (
 514                 regionName,
 515                 runTime.timeName(),
 516                 runTime,
 517                 IOobject::NO_READ
 518             ),
 519             xferCopy(pointField()),
 520             xferCopy(faceList()),
 521             xferCopy(cellList())
 522         );
 523
 524         for (label procI = 0; procI < nProcs; procI++)
 525         {
 526             Info<< "Reading mesh to add from "
 527                 << databases[procI].caseName()
 528                 << " for time = " << databases[procI].timeName()
 529                 << nl << endl;
 530
 531             polyMesh meshToAdd
 532             (
 533                 IOobject
 534                 (
 535                     regionName,
 536                     databases[procI].timeName(),
 537                     databases[procI]
 538                 )
 539             );
 540
 541             // Initialize its addressing
 542             cellProcAddressing[procI] = identity(meshToAdd.nCells());
 543             faceProcAddressing[procI] = identity(meshToAdd.nFaces());
 544             pointProcAddressing[procI] = identity(meshToAdd.nPoints());
 545             boundaryProcAddressing[procI] =
 546                 identity(meshToAdd.boundaryMesh().size());
 547
 548
 549             // Find geometrically shared points/faces.
 550             autoPtr<faceCoupleInfo> couples = determineCoupledFaces
 551             (
 552                 fullMatch,
 553                 procI,
 554                 masterMesh,
 555                 meshToAdd,
 556                 mergeDist
 557             );
 558
 559
 560             // Add elements to mesh
 561             Info<< "Adding to master mesh" << nl << endl;
 562
 563             autoPtr<mapAddedPolyMesh> map = polyMeshAdder::add
 564             (
 565                 masterMesh,
 566                 meshToAdd,
 567                 couples
 568             );
 569
 570             // Update all addressing so xxProcAddressing points to correct item
 571             // in masterMesh.
 572
 573             // Processors that were already in masterMesh
 574             for (label mergedI = 0; mergedI < procI; mergedI++)
 575             {
 576                 renumber(map().oldCellMap(), cellProcAddressing[mergedI]);
 577                 renumber(map().oldFaceMap(), faceProcAddressing[mergedI]);
 578                 renumber(map().oldPointMap(), pointProcAddressing[mergedI]);
 579                 // Note: boundary is special since can contain -1.
 580                 renumber(map().oldPatchMap(), boundaryProcAddressing[mergedI]);
 581             }
 582
 583             // Added processor
 584             renumber(map().addedCellMap(), cellProcAddressing[procI]);
 585             renumber(map().addedFaceMap(), faceProcAddressing[procI]);
 586             renumber(map().addedPointMap(), pointProcAddressing[procI]);
 587             renumber(map().addedPatchMap(), boundaryProcAddressing[procI]);
 588
 589             Info<< endl;
 590         }
 591
 592         // See if any points on the mastermesh have become connected
 593         // because of connections through processor meshes.
 594         mergeSharedPoints(mergeDist, masterMesh, pointProcAddressing);
 595
 596         // Save some properties on the reconstructed mesh
 597         masterInternalFaces = masterMesh.nInternalFaces();
 598         masterOwner = masterMesh.faceOwner();
 599
 600
 601         Info<< "\nWriting merged mesh to "
 602             << runTime.path()/runTime.timeName()
 603             << nl << endl;
 604
 605         if (!masterMesh.write())
 606         {
 607             FatalErrorIn(args.executable())
 608                 << "Failed writing polyMesh."
 609                 << exit(FatalError);
 610         }
 611     }
 612
 613
 614     // Write the addressing
 615
 616     Info<< "Reconstructing the addressing from the processor meshes"
 617         << " to the newly reconstructed mesh" << nl << endl;
 618
 619     forAll(databases, procI)
 620     {
 621         Info<< "Reading processor " << procI << " mesh from "
 622             << databases[procI].caseName() << endl;
 623
 624         polyMesh procMesh
 625         (
 626             IOobject
 627             (
 628                 regionName,
 629                 databases[procI].timeName(),
 630                 databases[procI]
 631             )
 632         );
 633
 634
 635         // From processor point to reconstructed mesh point
 636
 637         Info<< "Writing pointProcAddressing to "
 638             << databases[procI].caseName()
 639               /procMesh.facesInstance()
 640               /polyMesh::meshSubDir
 641             << endl;
 642
 643         labelIOList
 644         (
 645             IOobject
 646             (
 647                 "pointProcAddressing",
 648                 procMesh.facesInstance(),
 649                 polyMesh::meshSubDir,
 650                 procMesh,
 651                 IOobject::NO_READ,
 652                 IOobject::NO_WRITE,
 653                 false                       // do not register
 654             ),
 655             pointProcAddressing[procI]
 656         ).write();
 657
 658
 659         // From processor face to reconstructed mesh face
 660
 661         Info<< "Writing faceProcAddressing to "
 662             << databases[procI].caseName()
 663               /procMesh.facesInstance()
 664               /polyMesh::meshSubDir
 665             << endl;
 666
 667         labelIOList faceProcAddr
 668         (
 669             IOobject
 670             (
 671                 "faceProcAddressing",
 672                 procMesh.facesInstance(),
 673                 polyMesh::meshSubDir,
 674                 procMesh,
 675                 IOobject::NO_READ,
 676                 IOobject::NO_WRITE,
 677                 false                       // do not register
 678             ),
 679             faceProcAddressing[procI]
 680         );
 681
 682         // Now add turning index to faceProcAddressing.
 683         // See reconstrurPar for meaning of turning index.
 684         forAll(faceProcAddr, procFaceI)
 685         {
 686             label masterFaceI = faceProcAddr[procFaceI];
 687
 688             if
 689             (
 690                !procMesh.isInternalFace(procFaceI)
 691              && masterFaceI < masterInternalFaces
 692             )
 693             {
 694                 // proc face is now external but used to be internal face.
 695                 // Check if we have owner or neighbour.
 696
 697                 label procOwn = procMesh.faceOwner()[procFaceI];
 698                 label masterOwn = masterOwner[masterFaceI];
 699
 700                 if (cellProcAddressing[procI][procOwn] == masterOwn)
 701                 {
 702                     // No turning. Offset by 1.
 703                     faceProcAddr[procFaceI]++;
 704                 }
 705                 else
 706                 {
 707                     // Turned face.
 708                     faceProcAddr[procFaceI] =
 709                         -1 - faceProcAddr[procFaceI];
 710                 }
 711             }
 712             else
 713             {
 714                 // No turning. Offset by 1.
 715                 faceProcAddr[procFaceI]++;
 716             }
 717         }
 718
 719         faceProcAddr.write();
 720
 721
 722         // From processor cell to reconstructed mesh cell
 723
 724         Info<< "Writing cellProcAddressing to "
 725             << databases[procI].caseName()
 726               /procMesh.facesInstance()
 727               /polyMesh::meshSubDir
 728             << endl;
 729
 730         labelIOList
 731         (
 732             IOobject
 733             (
 734                 "cellProcAddressing",
 735                 procMesh.facesInstance(),
 736                 polyMesh::meshSubDir,
 737                 procMesh,
 738                 IOobject::NO_READ,
 739                 IOobject::NO_WRITE,
 740                 false                       // do not register
 741             ),
 742             cellProcAddressing[procI]
 743         ).write();
 744
 745
 746
 747         // From processor patch to reconstructed mesh patch
 748
 749         Info<< "Writing boundaryProcAddressing to "
 750             << databases[procI].caseName()
 751               /procMesh.facesInstance()
 752               /polyMesh::meshSubDir
 753             << endl;
 754
 755         labelIOList
 756         (
 757             IOobject
 758             (
 759                 "boundaryProcAddressing",
 760                 procMesh.facesInstance(),
 761                 polyMesh::meshSubDir,
 762                 procMesh,
 763                 IOobject::NO_READ,
 764                 IOobject::NO_WRITE,
 765                 false                       // do not register
 766             ),
 767             boundaryProcAddressing[procI]
 768         ).write();
 769
 770         Info<< endl;
 771     }
 772
 773     Info<< "End.\n" << endl;
 774
 775     return 0;
 776 }
 777
 778
 779 // ************************************************************************* //