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) 1991-2010 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
  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 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::noParallel();
 283     argList::validOptions.insert("mergeTol", "relative merge distance");
 284     argList::validOptions.insert("fullMatch", "");
 285
 286 #   include "addTimeOptions.H"
 287 #   include "addRegionOption.H"
 288 #   include "setRootCase.H"
 289 #   include "createTime.H"
 290
 291     Info<< "This is an experimental tool which tries to merge"
 292         << " individual processor" << nl
 293         << "meshes back into one master mesh. Use it if the original"
 294         << " master mesh has" << nl
 295         << "been deleted or if the processor meshes have been modified"
 296         << " (topology change)." << nl
 297         << "This tool will write the resulting mesh to a new time step"
 298         << " and construct" << nl
 299         << "xxxxProcAddressing files in the processor meshes so"
 300         << " reconstructPar can be" << nl
 301         << "used to regenerate the fields on the master mesh." << nl
 302         << nl
 303         << "Not well tested & use at your own risk!" << nl
 304         << endl;
 305
 306
 307     word regionName = polyMesh::defaultRegion;
 308     fileName regionPrefix = "";
 309     if (args.optionFound("region"))
 310     {
 311         regionName = args.option("region");
 312         regionPrefix = regionName;
 313         Info<< "Operating on region " << regionName << nl << endl;
 314     }
 315
 316     scalar mergeTol = defaultMergeTol;
 317     args.optionReadIfPresent("mergeTol", mergeTol);
 318
 319     scalar writeTol = Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
 320
 321     Info<< "Merge tolerance : " << mergeTol << nl
 322         << "Write tolerance : " << writeTol << endl;
 323
 324     if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
 325     {
 326         FatalErrorIn(args.executable())
 327             << "Your current settings specify ASCII writing with "
 328             << IOstream::defaultPrecision() << " digits precision." << endl
 329             << "Your merging tolerance (" << mergeTol << ") is finer than this."
 330             << endl
 331             << "Please change your writeFormat to binary"
 332             << " or increase the writePrecision" << endl
 333             << "or adjust the merge tolerance (-mergeTol)."
 334             << exit(FatalError);
 335     }
 336
 337
 338     const bool fullMatch = args.optionFound("fullMatch");
 339
 340     if (fullMatch)
 341     {
 342         Info<< "Doing geometric matching on all boundary faces." << nl << endl;
 343     }
 344     else
 345     {
 346         Info<< "Doing geometric matching on correct procBoundaries only."
 347             << nl << "This assumes a correct decomposition." << endl;
 348     }
 349
 350
 351
 352     int nProcs = 0;
 353
 354     while
 355     (
 356         isDir
 357         (
 358             args.rootPath()
 359           / args.caseName()
 360           / fileName(word("processor") + name(nProcs))
 361         )
 362     )
 363     {
 364         nProcs++;
 365     }
 366
 367     Info<< "Found " << nProcs << " processor directories" << nl << endl;
 368
 369
 370     // Read all databases.
 371     PtrList<Time> databases(nProcs);
 372
 373     forAll (databases, procI)
 374     {
 375         Info<< "Reading database "
 376             << args.caseName()/fileName(word("processor") + name(procI))
 377             << endl;
 378
 379         databases.set
 380         (
 381             procI,
 382             new Time
 383             (
 384                 Time::controlDictName,
 385                 args.rootPath(),
 386                 args.caseName()/fileName(word("processor") + name(procI))
 387             )
 388         );
 389
 390         Time& procTime = databases[procI];
 391
 392         instantList Times = procTime.times();
 393
 394         // set startTime and endTime depending on -time and -latestTime options
 395 #       include "checkTimeOptions.H"
 396
 397         procTime.setTime(Times[startTime], startTime);
 398
 399         if (procI > 0 && databases[procI-1].value() != procTime.value())
 400         {
 401             FatalErrorIn(args.executable())
 402                 << "Time not equal on processors." << nl
 403                 << "Processor:" << procI-1
 404                 << " time:" << databases[procI-1].value() << nl
 405                 << "Processor:" << procI
 406                 << " time:" << procTime.value()
 407                 << exit(FatalError);
 408         }
 409     }
 410
 411     // Set master time
 412     Info<< "Setting master time to " << databases[0].timeName() << nl << endl;
 413     runTime.setTime(databases[0]);
 414
 415
 416     // Read point on individual processors to determine merge tolerance
 417     // (otherwise single cell domains might give problems)
 418
 419     boundBox bb = boundBox::invertedBox;
 420
 421     for (label procI = 0; procI < nProcs; procI++)
 422     {
 423         fileName pointsInstance
 424         (
 425             databases[procI].findInstance
 426             (
 427                 regionPrefix/polyMesh::meshSubDir,
 428                 "points"
 429             )
 430         );
 431
 432         if (pointsInstance != databases[procI].timeName())
 433         {
 434             FatalErrorIn(args.executable())
 435                 << "Your time was specified as " << databases[procI].timeName()
 436                 << " but there is no polyMesh/points in that time." << endl
 437                 << "(there is a points file in " << pointsInstance
 438                 << ")" << endl
 439                 << "Please rerun with the correct time specified"
 440                 << " (through the -constant, -time or -latestTime option)."
 441                 << endl << exit(FatalError);
 442         }
 443
 444         Info<< "Reading points from "
 445             << databases[procI].caseName()
 446             << " for time = " << databases[procI].timeName()
 447             << nl << endl;
 448
 449         pointIOField points
 450         (
 451             IOobject
 452             (
 453                 "points",
 454                 databases[procI].findInstance
 455                 (
 456                     regionPrefix/polyMesh::meshSubDir,
 457                     "points"
 458                 ),
 459                 regionPrefix/polyMesh::meshSubDir,
 460                 databases[procI],
 461                 IOobject::MUST_READ,
 462                 IOobject::NO_WRITE,
 463                 false
 464             )
 465         );
 466
 467         boundBox domainBb(points, false);
 468
 469         bb.min() = min(bb.min(), domainBb.min());
 470         bb.max() = max(bb.max(), domainBb.max());
 471     }
 472     const scalar mergeDist = mergeTol * bb.mag();
 473
 474     Info<< "Overall mesh bounding box  : " << bb << nl
 475         << "Relative tolerance         : " << mergeTol << nl
 476         << "Absolute matching distance : " << mergeDist << nl
 477         << endl;
 478
 479
 480     // Addressing from processor to reconstructed case
 481     labelListList cellProcAddressing(nProcs);
 482     labelListList faceProcAddressing(nProcs);
 483     labelListList pointProcAddressing(nProcs);
 484     labelListList boundaryProcAddressing(nProcs);
 485
 486     // Internal faces on the final reconstructed mesh
 487     label masterInternalFaces;
 488     // Owner addressing on the final reconstructed mesh
 489     labelList masterOwner;
 490
 491     {
 492         // Construct empty mesh.
 493         Info<< "Constructing empty mesh to add to." << nl << endl;
 494         polyMesh masterMesh
 495         (
 496             IOobject
 497             (
 498                 regionName,
 499                 runTime.timeName(),
 500                 runTime,
 501                 IOobject::NO_READ
 502             ),
 503             xferCopy(pointField()),
 504             xferCopy(faceList()),
 505             xferCopy(cellList())
 506         );
 507
 508         for (label procI = 0; procI < nProcs; procI++)
 509         {
 510             Info<< "Reading mesh to add from "
 511                 << databases[procI].caseName()
 512                 << " for time = " << databases[procI].timeName()
 513                 << nl << endl;
 514
 515             polyMesh meshToAdd
 516             (
 517                 IOobject
 518                 (
 519                     regionName,
 520                     databases[procI].timeName(),
 521                     databases[procI]
 522                 )
 523             );
 524
 525             // Initialize its addressing
 526             cellProcAddressing[procI] = identity(meshToAdd.nCells());
 527             faceProcAddressing[procI] = identity(meshToAdd.nFaces());
 528             pointProcAddressing[procI] = identity(meshToAdd.nPoints());
 529             boundaryProcAddressing[procI] =
 530                 identity(meshToAdd.boundaryMesh().size());
 531
 532
 533             // Find geometrically shared points/faces.
 534             autoPtr<faceCoupleInfo> couples = determineCoupledFaces
 535             (
 536                 fullMatch,
 537                 procI,
 538                 masterMesh,
 539                 meshToAdd,
 540                 mergeDist
 541             );
 542
 543
 544             // Add elements to mesh
 545             Info<< "Adding to master mesh" << nl << endl;
 546
 547             autoPtr<mapAddedPolyMesh> map = polyMeshAdder::add
 548             (
 549                 masterMesh,
 550                 meshToAdd,
 551                 couples
 552             );
 553
 554             // Update all addressing so xxProcAddressing points to correct item
 555             // in masterMesh.
 556
 557             // Processors that were already in masterMesh
 558             for (label mergedI = 0; mergedI < procI; mergedI++)
 559             {
 560                 renumber(map().oldCellMap(), cellProcAddressing[mergedI]);
 561                 renumber(map().oldFaceMap(), faceProcAddressing[mergedI]);
 562                 renumber(map().oldPointMap(), pointProcAddressing[mergedI]);
 563                 // Note: boundary is special since can contain -1.
 564                 renumber(map().oldPatchMap(), boundaryProcAddressing[mergedI]);
 565             }
 566
 567             // Added processor
 568             renumber(map().addedCellMap(), cellProcAddressing[procI]);
 569             renumber(map().addedFaceMap(), faceProcAddressing[procI]);
 570             renumber(map().addedPointMap(), pointProcAddressing[procI]);
 571             renumber(map().addedPatchMap(), boundaryProcAddressing[procI]);
 572
 573             Info<< endl;
 574         }
 575
 576         // See if any points on the mastermesh have become connected
 577         // because of connections through processor meshes.
 578         mergeSharedPoints(mergeDist, masterMesh, pointProcAddressing);
 579
 580         // Save some properties on the reconstructed mesh
 581         masterInternalFaces = masterMesh.nInternalFaces();
 582         masterOwner = masterMesh.faceOwner();
 583
 584
 585         Info<< "\nWriting merged mesh to "
 586             << runTime.path()/runTime.timeName()
 587             << nl << endl;
 588
 589         if (!masterMesh.write())
 590         {
 591             FatalErrorIn(args.executable())
 592                 << "Failed writing polyMesh."
 593                 << exit(FatalError);
 594         }
 595     }
 596
 597
 598     // Write the addressing
 599
 600     Info<< "Reconstructing the addressing from the processor meshes"
 601         << " to the newly reconstructed mesh" << nl << endl;
 602
 603     forAll(databases, procI)
 604     {
 605         Info<< "Reading processor " << procI << " mesh from "
 606             << databases[procI].caseName() << endl;
 607
 608         polyMesh procMesh
 609         (
 610             IOobject
 611             (
 612                 regionName,
 613                 databases[procI].timeName(),
 614                 databases[procI]
 615             )
 616         );
 617
 618
 619         // From processor point to reconstructed mesh point
 620
 621         Info<< "Writing pointProcAddressing to "
 622             << databases[procI].caseName()
 623               /procMesh.facesInstance()
 624               /polyMesh::meshSubDir
 625             << endl;
 626
 627         labelIOList
 628         (
 629             IOobject
 630             (
 631                 "pointProcAddressing",
 632                 procMesh.facesInstance(),
 633                 polyMesh::meshSubDir,
 634                 procMesh,
 635                 IOobject::NO_READ,
 636                 IOobject::NO_WRITE,
 637                 false                       // do not register
 638             ),
 639             pointProcAddressing[procI]
 640         ).write();
 641
 642
 643         // From processor face to reconstructed mesh face
 644
 645         Info<< "Writing faceProcAddressing to "
 646             << databases[procI].caseName()
 647               /procMesh.facesInstance()
 648               /polyMesh::meshSubDir
 649             << endl;
 650
 651         labelIOList faceProcAddr
 652         (
 653             IOobject
 654             (
 655                 "faceProcAddressing",
 656                 procMesh.facesInstance(),
 657                 polyMesh::meshSubDir,
 658                 procMesh,
 659                 IOobject::NO_READ,
 660                 IOobject::NO_WRITE,
 661                 false                       // do not register
 662             ),
 663             faceProcAddressing[procI]
 664         );
 665
 666         // Now add turning index to faceProcAddressing.
 667         // See reconstrurPar for meaning of turning index.
 668         forAll(faceProcAddr, procFaceI)
 669         {
 670             label masterFaceI = faceProcAddr[procFaceI];
 671
 672             if
 673             (
 674                !procMesh.isInternalFace(procFaceI)
 675              && masterFaceI < masterInternalFaces
 676             )
 677             {
 678                 // proc face is now external but used to be internal face.
 679                 // Check if we have owner or neighbour.
 680
 681                 label procOwn = procMesh.faceOwner()[procFaceI];
 682                 label masterOwn = masterOwner[masterFaceI];
 683
 684                 if (cellProcAddressing[procI][procOwn] == masterOwn)
 685                 {
 686                     // No turning. Offset by 1.
 687                     faceProcAddr[procFaceI]++;
 688                 }
 689                 else
 690                 {
 691                     // Turned face.
 692                     faceProcAddr[procFaceI] =
 693                         -1 - faceProcAddr[procFaceI];
 694                 }
 695             }
 696             else
 697             {
 698                 // No turning. Offset by 1.
 699                 faceProcAddr[procFaceI]++;
 700             }
 701         }
 702
 703         faceProcAddr.write();
 704
 705
 706         // From processor cell to reconstructed mesh cell
 707
 708         Info<< "Writing cellProcAddressing to "
 709             << databases[procI].caseName()
 710               /procMesh.facesInstance()
 711               /polyMesh::meshSubDir
 712             << endl;
 713
 714         labelIOList
 715         (
 716             IOobject
 717             (
 718                 "cellProcAddressing",
 719                 procMesh.facesInstance(),
 720                 polyMesh::meshSubDir,
 721                 procMesh,
 722                 IOobject::NO_READ,
 723                 IOobject::NO_WRITE,
 724                 false                       // do not register
 725             ),
 726             cellProcAddressing[procI]
 727         ).write();
 728
 729
 730
 731         // From processor patch to reconstructed mesh patch
 732
 733         Info<< "Writing boundaryProcAddressing to "
 734             << databases[procI].caseName()
 735               /procMesh.facesInstance()
 736               /polyMesh::meshSubDir
 737             << endl;
 738
 739         labelIOList
 740         (
 741             IOobject
 742             (
 743                 "boundaryProcAddressing",
 744                 procMesh.facesInstance(),
 745                 polyMesh::meshSubDir,
 746                 procMesh,
 747                 IOobject::NO_READ,
 748                 IOobject::NO_WRITE,
 749                 false                       // do not register
 750             ),
 751             boundaryProcAddressing[procI]
 752         ).write();
 753
 754         Info<< endl;
 755     }
 756
 757     Info<< "End.\n" << endl;
 758
 759     return 0;
 760 }
 761
 762
 763 // ************************************************************************* //