applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

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