applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Application
  26     reconstructParMesh
  27
  28 Author
  29     Hrvoje Jasak, Wikki Ltd.  All rights reserved
  30
  31 Description
  32     Reconstructs a mesh using geometrical matching and catenation.
  33     Use following topological changes in parallel to create global mesh
  34     and xxxxProcAddressing files in the processor meshes.
  35
  36 \*---------------------------------------------------------------------------*/
  37
  38 #include "fvCFD.H"
  39 #include "processorMeshesReconstructor.H"
  40 #include "fvFieldReconstructor.H"
  41 #include "pointFieldReconstructor.H"
  42 #include "tetPointFieldReconstructor.H"
  43 #include "reconstructLagrangian.H"
  44
  45 #include "faCFD.H"
  46 #include "faMesh.H"
  47 #include "processorFaMeshes.H"
  48 #include "faFieldReconstructor.H"
  49
  50 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  51
  52 int main(int argc, char *argv[])
  53 {
  54     // enable -constant ... if someone really wants it
  55     // enable -zeroTime to prevent accidentally trashing the initial fields
  56     timeSelector::addOptions(false, true);
  57     argList::noParallel();
  58 #   include "addRegionOption.H"
  59     argList::validOptions.insert("cellDist", "");
  60     argList::validOptions.insert("fields", "\"(list of fields)\"");
  61     argList::validOptions.insert("noLagrangian", "");
  62
  63 #   include "setRootCase.H"
  64
  65     bool writeCellDist = args.optionFound("cellDist");
  66
  67 #   include "createTime.H"
  68
  69     HashSet<word> selectedFields;
  70     if (args.optionFound("fields"))
  71     {
  72         args.optionLookup("fields")() >> selectedFields;
  73     }
  74
  75     bool noLagrangian = args.optionFound("noLagrangian");
  76
  77     // Determine the processor count directly
  78     label nProcs = 0;
  79     while (isDir(args.path()/(word("processor") + name(nProcs))))
  80     {
  81         ++nProcs;
  82     }
  83
  84     if (!nProcs)
  85     {
  86         FatalErrorIn(args.executable())
  87             << "No processor* directories found"
  88             << exit(FatalError);
  89     }
  90
  91     // Create the processor databases
  92     PtrList<Time> databases(nProcs);
  93
  94     forAll (databases, procI)
  95     {
  96         databases.set
  97         (
  98             procI,
  99             new Time
 100             (
 101                 Time::controlDictName,
 102                 args.rootPath(),
 103                 args.caseName()/fileName(word("processor") + name(procI))
 104             )
 105         );
 106     }
 107
 108     // use the times list from the master processor
 109     // and select a subset based on the command-line options
 110     instantList timeDirs = timeSelector::select
 111     (
 112         databases[0].times(),
 113         args
 114     );
 115
 116     if (timeDirs.empty())
 117     {
 118         FatalErrorIn(args.executable())
 119             << "No times selected"
 120             << exit(FatalError);
 121     }
 122
 123     Foam::word regionName = polyMesh::defaultRegion;
 124
 125     if (args.optionReadIfPresent("region", regionName))
 126     {
 127         Info<< "Selecting region " << regionName << " for time = "
 128             << runTime.timeName() << Foam::nl << Foam::endl;
 129     }
 130
 131     // Set all times on processor meshes equal to reconstructed mesh
 132     forAll (databases, procI)
 133     {
 134         Info<< "Reading database for processor " << procI << endl;
 135
 136         databases[procI].setTime(runTime.timeName(), runTime.timeIndex());
 137     }
 138
 139     // Read all meshes and addressing to reconstructed mesh
 140     processorMeshesReconstructor procMeshes(databases, regionName);
 141
 142     autoPtr<fvMesh> meshPtr = procMeshes.reconstructMesh(runTime);
 143
 144
 145     // Mesh write will be controlled by hand
 146     meshPtr->write();
 147     procMeshes.writeAddressing();
 148     meshPtr->setMotionWriteOpt(IOobject::NO_WRITE);
 149     meshPtr->setTopoWriteOpt(IOobject::NO_WRITE);
 150
 151     // Write cell decomposition
 152     if (writeCellDist)
 153     {
 154         // Write as volScalarField for post-processing
 155         Info<< "Writing cellDist to time " << runTime.timeName()
 156             << endl;
 157         volScalarField cellDist
 158         (
 159             IOobject
 160             (
 161                 "cellDist",
 162                 runTime.timeName(),
 163                 meshPtr(),
 164                 IOobject::NO_READ,
 165                 IOobject::NO_WRITE
 166             ),
 167             meshPtr(),
 168             dimensionedScalar("cellDist", dimless, 0),
 169             zeroGradientFvPatchScalarField::typeName
 170         );
 171         scalarField& cellDistIn = cellDist.internalField();
 172
 173         label cellI = 0;
 174
 175         forAll (procMeshes.meshes(), procI)
 176         {
 177             for
 178             (
 179                 label i = 0;
 180                 i < procMeshes.meshes()[procI].nCells();
 181                 i++
 182             )
 183             {
 184                 cellDistIn[cellI] = procI;
 185                 cellI++;
 186             }
 187         }
 188
 189         cellDist.write();
 190     }
 191
 192     // Get region prefix for lagrangian
 193     fileName regionPrefix = "";
 194     if (regionName != fvMesh::defaultRegion)
 195     {
 196         regionPrefix = regionName;
 197     }
 198
 199
 200     // Loop over all times
 201     forAll (timeDirs, timeI)
 202     {
 203         // Set time for global database
 204         runTime.setTime(timeDirs[timeI], timeI);
 205
 206         Info << "Time = " << runTime.timeName() << endl << endl;
 207
 208         // Set time for all databases
 209         forAll (databases, procI)
 210         {
 211             databases[procI].setTime(timeDirs[timeI], timeI);
 212         }
 213
 214         polyMesh::readUpdateState procStat = procMeshes.readUpdate();
 215
 216         if (procStat == polyMesh::UNCHANGED)
 217         {
 218             Info<< "Mesh unchanged" << endl;
 219
 220             meshPtr->setMotionWriteOpt(IOobject::NO_WRITE);
 221             meshPtr->setTopoWriteOpt(IOobject::NO_WRITE);
 222         }
 223         else if (procStat == polyMesh::POINTS_MOVED)
 224         {
 225             Info<< "Mesh motion detected.  Reconstruct motion points"
 226                 << endl;
 227
 228             // Reconstruct the points for moving mesh cases and write them out
 229             procMeshes.reconstructPoints(meshPtr());
 230
 231             // Set write options
 232             meshPtr->setMotionWriteOpt(IOobject::AUTO_WRITE);
 233             meshPtr->setTopoWriteOpt(IOobject::NO_WRITE);
 234
 235             // Global mesh write
 236             meshPtr->write();
 237         }
 238         else if
 239         (
 240             procStat == polyMesh::TOPO_CHANGE
 241          || procStat == polyMesh::TOPO_PATCH_CHANGE
 242         )
 243         {
 244             Info<< "Topological change detected.  Reconstructing mesh"
 245                 << endl;
 246
 247             // Reconstruct mesh
 248             meshPtr = procMeshes.reconstructMesh(runTime);
 249
 250             // Set write options
 251             meshPtr->setMotionWriteOpt(IOobject::AUTO_WRITE);
 252             meshPtr->setTopoWriteOpt(IOobject::AUTO_WRITE);
 253             procMeshes.writeAddressing();
 254
 255             // Global mesh write
 256             meshPtr->write();
 257
 258             // Write out mapping in processor directories
 259             forAll (databases, procI)
 260             {
 261                 databases[procI].write();
 262             }
 263
 264             if (writeCellDist)
 265             {
 266                 // Write as volScalarField for post-processing
 267                 Info<< "Writing cellDist to time " << runTime.timeName()
 268                     << endl;
 269                 volScalarField cellDist
 270                 (
 271                     IOobject
 272                     (
 273                         "cellDist",
 274                         runTime.timeName(),
 275                         meshPtr(),
 276                         IOobject::NO_READ,
 277                         IOobject::NO_WRITE
 278                     ),
 279                     meshPtr(),
 280                     dimensionedScalar("cellDist", dimless, 0),
 281                     zeroGradientFvPatchScalarField::typeName
 282                 );
 283                 scalarField& cellDistIn = cellDist.internalField();
 284
 285                 label cellI = 0;
 286
 287                 forAll (procMeshes.meshes(), procI)
 288                 {
 289                     for
 290                     (
 291                         label i = 0;
 292                         i < procMeshes.meshes()[procI].nCells();
 293                         i++
 294                     )
 295                     {
 296                         cellDistIn[cellI] = procI;
 297                         cellI++;
 298                     }
 299                 }
 300
 301                 cellDist.write();
 302             }
 303         }
 304         else
 305         {
 306             FatalErrorIn(args.executable())
 307                 << "Unknown readUpdate state"
 308                 << abort(FatalError);
 309         }
 310
 311         fvMesh& mesh = meshPtr();
 312
 313         // Get list of objects from processor0 database
 314         IOobjectList objects(procMeshes.meshes()[0], databases[0].timeName());
 315
 316
 317         // If there are any FV fields, reconstruct them
 318
 319         if
 320         (
 321             objects.lookupClass(volScalarField::typeName).size()
 322          || objects.lookupClass(volVectorField::typeName).size()
 323          || objects.lookupClass(volSphericalTensorField::typeName).size()
 324          || objects.lookupClass(volSymmTensorField::typeName).size()
 325          || objects.lookupClass(volTensorField::typeName).size()
 326          || objects.lookupClass(surfaceScalarField::typeName).size()
 327          || objects.lookupClass(surfaceVectorField::typeName).size()
 328          || objects.lookupClass(surfaceSphericalTensorField::typeName).size()
 329          || objects.lookupClass(surfaceSymmTensorField::typeName).size()
 330          || objects.lookupClass(surfaceTensorField::typeName).size()
 331         )
 332         {
 333             Info << "Reconstructing FV fields" << nl << endl;
 334
 335             fvFieldReconstructor fvReconstructor
 336             (
 337                 mesh,
 338                 procMeshes.meshes(),
 339                 procMeshes.faceProcAddressing(),
 340                 procMeshes.cellProcAddressing(),
 341                 procMeshes.boundaryProcAddressing()
 342             );
 343
 344             fvReconstructor.reconstructFvVolumeFields<scalar>
 345             (
 346                 objects,
 347                 selectedFields
 348             );
 349             fvReconstructor.reconstructFvVolumeFields<vector>
 350             (
 351                 objects,
 352                 selectedFields
 353             );
 354             fvReconstructor.reconstructFvVolumeFields<sphericalTensor>
 355             (
 356                 objects,
 357                 selectedFields
 358             );
 359             fvReconstructor.reconstructFvVolumeFields<symmTensor>
 360             (
 361                 objects,
 362                 selectedFields
 363             );
 364             fvReconstructor.reconstructFvVolumeFields<tensor>
 365             (
 366                 objects,
 367                 selectedFields
 368             );
 369
 370             fvReconstructor.reconstructFvSurfaceFields<scalar>
 371             (
 372                 objects,
 373                 selectedFields
 374             );
 375             fvReconstructor.reconstructFvSurfaceFields<vector>
 376             (
 377                 objects,
 378                 selectedFields
 379             );
 380             fvReconstructor.reconstructFvSurfaceFields<sphericalTensor>
 381             (
 382                 objects,
 383                 selectedFields
 384             );
 385             fvReconstructor.reconstructFvSurfaceFields<symmTensor>
 386             (
 387                 objects,
 388                 selectedFields
 389             );
 390             fvReconstructor.reconstructFvSurfaceFields<tensor>
 391             (
 392                 objects,
 393                 selectedFields
 394             );
 395         }
 396         else
 397         {
 398             Info << "No FV fields" << nl << endl;
 399         }
 400
 401
 402         // If there are any point fields, reconstruct them
 403         if
 404         (
 405             objects.lookupClass(pointScalarField::typeName).size()
 406          || objects.lookupClass(pointVectorField::typeName).size()
 407          || objects.lookupClass(pointSphericalTensorField::typeName).size()
 408          || objects.lookupClass(pointSymmTensorField::typeName).size()
 409          || objects.lookupClass(pointTensorField::typeName).size()
 410         )
 411         {
 412             Info << "Reconstructing point fields" << nl << endl;
 413
 414             pointMesh pMesh(mesh);
 415             PtrList<pointMesh> pMeshes(procMeshes.meshes().size());
 416
 417             forAll (pMeshes, procI)
 418             {
 419                 pMeshes.set(procI, new pointMesh(procMeshes.meshes()[procI]));
 420             }
 421
 422             pointFieldReconstructor pointReconstructor
 423             (
 424                 pMesh,
 425                 pMeshes,
 426                 procMeshes.pointProcAddressing(),
 427                 procMeshes.boundaryProcAddressing()
 428             );
 429
 430             pointReconstructor.reconstructFields<scalar>(objects);
 431             pointReconstructor.reconstructFields<vector>(objects);
 432             pointReconstructor.reconstructFields<sphericalTensor>(objects);
 433             pointReconstructor.reconstructFields<symmTensor>(objects);
 434             pointReconstructor.reconstructFields<tensor>(objects);
 435         }
 436         else
 437         {
 438             Info << "No point fields" << nl << endl;
 439         }
 440
 441         // If there are any tetFem fields, reconstruct them
 442         if
 443         (
 444             objects.lookupClass(tetPointScalarField::typeName).size()
 445          || objects.lookupClass(tetPointVectorField::typeName).size()
 446          || objects.lookupClass(tetPointSphericalTensorField::typeName).size()
 447          || objects.lookupClass(tetPointSymmTensorField::typeName).size()
 448          || objects.lookupClass(tetPointTensorField::typeName).size()
 449
 450          || objects.lookupClass(elementScalarField::typeName).size()
 451          || objects.lookupClass(elementVectorField::typeName).size()
 452         )
 453         {
 454             Info << "Reconstructing tet point fields" << nl << endl;
 455
 456             tetPolyMesh tetMesh(mesh);
 457             PtrList<tetPolyMesh> tetMeshes(procMeshes.meshes().size());
 458
 459             forAll (tetMeshes, procI)
 460             {
 461                 tetMeshes.set
 462                 (
 463                     procI,
 464                     new tetPolyMesh(procMeshes.meshes()[procI])
 465                 );
 466             }
 467
 468             tetPointFieldReconstructor tetPointReconstructor
 469             (
 470                 tetMesh,
 471                 tetMeshes,
 472                 procMeshes.pointProcAddressing(),
 473                 procMeshes.faceProcAddressing(),
 474                 procMeshes.cellProcAddressing(),
 475                 procMeshes.boundaryProcAddressing()
 476             );
 477
 478             // Reconstruct tet point fields
 479             tetPointReconstructor.reconstructTetPointFields<scalar>(objects);
 480             tetPointReconstructor.reconstructTetPointFields<vector>(objects);
 481             tetPointReconstructor.
 482                 reconstructTetPointFields<sphericalTensor>(objects);
 483             tetPointReconstructor.
 484                 reconstructTetPointFields<symmTensor>(objects);
 485             tetPointReconstructor.reconstructTetPointFields<tensor>(objects);
 486
 487             tetPointReconstructor.reconstructElementFields<scalar>(objects);
 488             tetPointReconstructor.reconstructElementFields<vector>(objects);
 489         }
 490         else
 491         {
 492             Info << "No tetFem fields" << nl << endl;
 493         }
 494
 495
 496         // If there are any clouds, reconstruct them.
 497         // The problem is that a cloud of size zero will not get written so
 498         // in pass 1 we determine the cloud names and per cloud name the
 499         // fields. Note that the fields are stored as IOobjectList from
 500         // the first processor that has them. They are in pass2 only used
 501         // for name and type (scalar, vector etc).
 502
 503         if (!noLagrangian)
 504         {
 505             HashTable<IOobjectList> cloudObjects;
 506
 507             forAll (databases, procI)
 508             {
 509                 fileNameList cloudDirs
 510                 (
 511                     readDir
 512                     (
 513                         databases[procI].timePath()/regionPrefix/cloud::prefix,
 514                         fileName::DIRECTORY
 515                     )
 516                 );
 517
 518                 forAll (cloudDirs, i)
 519                 {
 520                     // Check if we already have cloud objects for
 521                     // this cloudname
 522                     HashTable<IOobjectList>::const_iterator iter =
 523                         cloudObjects.find(cloudDirs[i]);
 524
 525                     if (iter == cloudObjects.end())
 526                     {
 527                         // Do local scan for valid cloud objects
 528                         IOobjectList sprayObjs
 529                         (
 530                             procMeshes.meshes()[procI],
 531                             databases[procI].timeName(),
 532                             cloud::prefix/cloudDirs[i]
 533                         );
 534
 535                         IOobject* positionsPtr = sprayObjs.lookup("positions");
 536
 537                         if (positionsPtr)
 538                         {
 539                             cloudObjects.insert(cloudDirs[i], sprayObjs);
 540                         }
 541                     }
 542                 }
 543             }
 544
 545
 546             if (cloudObjects.size())
 547             {
 548                 // Pass2: reconstruct the cloud
 549                 forAllConstIter(HashTable<IOobjectList>, cloudObjects, iter)
 550                 {
 551                     const word cloudName = string::validate<word>(iter.key());
 552
 553                     // Objects (on arbitrary processor)
 554                     const IOobjectList& sprayObjs = iter();
 555
 556                     Info<< "Reconstructing lagrangian fields for cloud "
 557                         << cloudName << nl << endl;
 558
 559                     reconstructLagrangianPositions
 560                     (
 561                         mesh,
 562                         cloudName,
 563                         procMeshes.meshes(),
 564                         procMeshes.faceProcAddressing(),
 565                         procMeshes.cellProcAddressing()
 566                     );
 567                     reconstructLagrangianFields<label>
 568                     (
 569                         cloudName,
 570                         mesh,
 571                         procMeshes.meshes(),
 572                         sprayObjs
 573                     );
 574                     reconstructLagrangianFields<scalar>
 575                     (
 576                         cloudName,
 577                         mesh,
 578                         procMeshes.meshes(),
 579                         sprayObjs
 580                     );
 581                     reconstructLagrangianFields<vector>
 582                     (
 583                         cloudName,
 584                         mesh,
 585                         procMeshes.meshes(),
 586                         sprayObjs
 587                     );
 588                     reconstructLagrangianFields<sphericalTensor>
 589                     (
 590                         cloudName,
 591                         mesh,
 592                         procMeshes.meshes(),
 593                         sprayObjs
 594                     );
 595                     reconstructLagrangianFields<symmTensor>
 596                     (
 597                         cloudName,
 598                         mesh,
 599                         procMeshes.meshes(),
 600                         sprayObjs
 601                     );
 602                     reconstructLagrangianFields<tensor>
 603                     (
 604                         cloudName,
 605                         mesh,
 606                         procMeshes.meshes(),
 607                         sprayObjs
 608                     );
 609                 }
 610             }
 611             else
 612             {
 613                 Info << "No lagrangian fields" << nl << endl;
 614             }
 615         }
 616
 617         // If there are any FA fields, reconstruct them
 618
 619         if
 620         (
 621             objects.lookupClass(areaScalarField::typeName).size()
 622          || objects.lookupClass(areaVectorField::typeName).size()
 623          || objects.lookupClass(areaSphericalTensorField::typeName).size()
 624          || objects.lookupClass(areaSymmTensorField::typeName).size()
 625          || objects.lookupClass(areaTensorField::typeName).size()
 626          || objects.lookupClass(edgeScalarField::typeName).size()
 627         )
 628         {
 629             Info << "Reconstructing FA fields" << nl << endl;
 630
 631             faMesh aMesh(mesh);
 632
 633             processorFaMeshes procFaMeshes(procMeshes.meshes());
 634
 635             faFieldReconstructor faReconstructor
 636             (
 637                 aMesh,
 638                 procFaMeshes.meshes(),
 639                 procFaMeshes.edgeProcAddressing(),
 640                 procFaMeshes.faceProcAddressing(),
 641                 procFaMeshes.boundaryProcAddressing()
 642             );
 643
 644             faReconstructor.reconstructFaAreaFields<scalar>(objects);
 645             faReconstructor.reconstructFaAreaFields<vector>(objects);
 646             faReconstructor
 647                .reconstructFaAreaFields<sphericalTensor>(objects);
 648             faReconstructor.reconstructFaAreaFields<symmTensor>(objects);
 649             faReconstructor.reconstructFaAreaFields<tensor>(objects);
 650
 651             faReconstructor.reconstructFaEdgeFields<scalar>(objects);
 652         }
 653         else
 654         {
 655             Info << "No FA fields" << nl << endl;
 656         }
 657
 658         // If there are any "uniform" directories copy them from
 659         // the master processor
 660
 661         fileName uniformDir0 = databases[0].timePath()/"uniform";
 662         if (isDir(uniformDir0))
 663         {
 664             cp(uniformDir0, runTime.timePath());
 665         }
 666     }
 667
 668     Info<< "End.\n" << endl;
 669
 670     return 0;
 671 }
 672
 673
 674 // ************************************************************************* //