decomposePar.C 55.1 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6 7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2017 OpenFOAM Foundation
9
    Copyright (C) 2016-2020 OpenCFD Ltd.
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

14 15 16 17
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
18 19 20 21 22 23 24

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
26 27 28 29

Application
    decomposePar

30 31 32
Group
    grpParallelUtilities

33
Description
34 35
    Automatically decomposes a mesh and fields of a case for parallel
    execution of OpenFOAM.
36

37
Usage
38
    \b decomposePar [OPTIONS]
39

40
    Options:
41 42 43 44
      - \par -allRegions
        Decompose all regions in regionProperties. Does not check for
        existence of processor*.

45
      - \par -case \<dir\>
46
        Specify case directory to use (instead of the cwd).
47

48 49
      - \par -cellDist
        Write the cell distribution as a labelList, for use with 'manual'
50
        decomposition method and as a volScalarField for visualization.
51

52 53
      - \par -constant
        Include the 'constant/' dir in the times list.
54

55 56
      - \par -copyUniform
        Copy any \a uniform directories too.
mattijs's avatar
mattijs committed
57

58
      - \par -copyZero
59 60
        Copy \a 0 directory to processor* rather than decompose the fields.

61
      - \par -debug-switch \<name=val\>
62 63
        Specify the value of a registered debug switch. Default is 1
        if the value is omitted. (Can be used multiple times)
64

65
      - \par -decomposeParDict \<file\>
66
        Use specified file for decomposePar dictionary.
67

68 69 70
      - \par -dry-run
        Test without writing the decomposition. Changes -cellDist to
        only write volScalarField.
71

72
      - \par -fields
73
        Use existing geometry decomposition and convert fields only.
74

75
      - \par fileHandler \<handler\>
76
        Override the file handler type.
77

78
      - \par -force
79 80
        Remove any existing \a processor subdirectories before decomposing the
        geometry.
Mark Olesen's avatar
Mark Olesen committed
81

82
      - \par -ifRequired
83 84 85 86 87 88
        Only decompose the geometry if the number of domains has changed from a
        previous decomposition. No \a processor subdirectories will be removed
        unless the \a -force option is also specified. This option can be used
        to avoid redundant geometry decomposition (eg, in scripts), but should
        be used with caution when the underlying (serial) geometry or the
        decomposition method etc. have been changed between decompositions.
89

90
      - \par -info-switch \<name=val\>
91 92 93 94 95 96
        Specify the value of a registered info switch. Default is 1
        if the value is omitted. (Can be used multiple times)

      - \par -latestTime
        Select the latest time.

97
      - \par -lib \<name\>
98 99 100 101 102 103 104 105 106 107 108
        Additional library or library list to load (can be used multiple times).

      - \par -noFunctionObjects
        Do not execute function objects.

      - \par -noSets
        Skip decomposing cellSets, faceSets, pointSets.

      - \par -noZero
        Exclude the \a 0 dir from the times list.

109
      - \par -opt-switch \<name=val\>
110 111 112 113 114 115
        Specify the value of a registered optimisation switch (int/bool).
        Default is 1 if the value is omitted. (Can be used multiple times)

      - \par -region \<regionName\>
        Decompose named region. Does not check for existence of processor*.

116
      - \par -time \<ranges\>
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
        Override controlDict settings and decompose selected times. Does not
        re-decompose the mesh i.e. does not handle moving mesh or changing
        mesh cases. Eg, ':10,20 40:70 1000:', 'none', etc.

      - \par -verbose
        Additional verbosity.

      - \par -doc
        Display documentation in browser.

      - \par -doc-source
        Display source code in browser.

      - \par -help
        Display short help and exit.

      - \par -help-man
        Display full help (manpage format) and exit.

      - \par -help-notes
        Display help notes (description) and exit.

      - \par -help-full
        Display full help and exit.

142 143 144 145 146 147
\*---------------------------------------------------------------------------*/

#include "OSspecific.H"
#include "fvCFD.H"
#include "IOobjectList.H"
#include "domainDecomposition.H"
148
#include "domainDecompositionDryRun.H"
149
#include "labelIOField.H"
150
#include "labelFieldIOField.H"
151
#include "scalarIOField.H"
152
#include "scalarFieldIOField.H"
153
#include "vectorIOField.H"
154
#include "vectorFieldIOField.H"
155
#include "sphericalTensorIOField.H"
156
#include "sphericalTensorFieldIOField.H"
157
#include "symmTensorIOField.H"
158
#include "symmTensorFieldIOField.H"
159
#include "tensorIOField.H"
160
#include "tensorFieldIOField.H"
161
#include "pointFields.H"
162
#include "regionProperties.H"
163 164

#include "readFields.H"
165
#include "dimFieldDecomposer.H"
166 167 168
#include "fvFieldDecomposer.H"
#include "pointFieldDecomposer.H"
#include "lagrangianFieldDecomposer.H"
169
#include "decompositionModel.H"
170

Hrvoje Jasak's avatar
Hrvoje Jasak committed
171 172 173 174 175
#include "faCFD.H"
#include "emptyFaPatch.H"
#include "faMeshDecomposition.H"
#include "faFieldDecomposer.H"

176 177
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

178 179 180
namespace Foam
{

181 182 183
const labelIOList& procAddressing
(
    const PtrList<fvMesh>& procMeshList,
184
    const label proci,
185 186 187 188
    const word& name,
    PtrList<labelIOList>& procAddressingList
)
{
189
    const fvMesh& procMesh = procMeshList[proci];
190

191
    if (!procAddressingList.set(proci))
192 193 194
    {
        procAddressingList.set
        (
195
            proci,
196 197 198 199 200 201 202 203 204
            new labelIOList
            (
                IOobject
                (
                    name,
                    procMesh.facesInstance(),
                    procMesh.meshSubDir,
                    procMesh,
                    IOobject::MUST_READ,
205 206
                    IOobject::NO_WRITE,
                    false
207 208 209 210
                )
            )
        );
    }
211
    return procAddressingList[proci];
212 213 214
}


215 216 217 218 219 220 221 222 223 224 225 226 227
void decomposeUniform
(
    const bool copyUniform,
    const domainDecomposition& mesh,
    const Time& processorDb,
    const word& regionDir = word::null
)
{
    const Time& runTime = mesh.time();

    // Any uniform data to copy/link?
    const fileName uniformDir(regionDir/"uniform");

228
    if (fileHandler().isDir(runTime.timePath()/uniformDir))
229 230 231 232 233
    {
        Info<< "Detected additional non-decomposed files in "
            << runTime.timePath()/uniformDir
            << endl;

234 235
        const fileName timePath =
            fileHandler().filePath(processorDb.timePath());
236

237 238 239 240
        // If no fields have been decomposed the destination
        // directory will not have been created so make sure.
        mkDir(timePath);

241 242
        if (copyUniform || mesh.distributed())
        {
243 244 245 246 247 248 249 250
            if (!fileHandler().exists(timePath/uniformDir))
            {
                fileHandler().cp
                (
                    runTime.timePath()/uniformDir,
                    timePath/uniformDir
                );
            }
251 252 253
        }
        else
        {
254
            // Link with relative paths
255 256 257 258 259 260 261 262 263
            string parentPath = string("..")/"..";

            if (regionDir != word::null)
            {
                parentPath = parentPath/"..";
            }

            fileName currentDir(cwd());
            chDir(timePath);
264 265 266 267 268 269 270 271 272

            if (!fileHandler().exists(uniformDir))
            {
                fileHandler().ln
                (
                    parentPath/runTime.timeName()/uniformDir,
                    uniformDir
                );
            }
273 274 275 276 277
            chDir(currentDir);
        }
    }
}

278 279 280
}


281
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
282

283 284
int main(int argc, char *argv[])
{
285 286
    argList::addNote
    (
287
        "Decompose a mesh and fields of a case for parallel execution"
288 289
    );

290
    argList::noParallel();
291
    argList::addOption
292 293 294
    (
        "decomposeParDict",
        "file",
295
        "Use specified file for decomposePar dictionary"
296
    );
297
    #include "addRegionOption.H"
298
    argList::addBoolOption
299 300
    (
        "allRegions",
301
        "Operate on all regions in regionProperties"
302 303
    );
    argList::addBoolOption
304 305 306 307 308 309 310 311 312 313 314
    (
        "dry-run",
        "Test without writing the decomposition. "
        "Changes -cellDist to only write volScalarField."
    );
    argList::addBoolOption
    (
        "verbose",
        "Additional verbosity"
    );
    argList::addBoolOption
315 316
    (
        "cellDist",
317 318
        "Write cell distribution as a labelList - for use with 'manual' "
        "decomposition method and as a volScalarField for visualization."
319
    );
320
    argList::addBoolOption
321 322
    (
        "copyZero",
323
        "Copy 0/ directory to processor*/ rather than decompose the fields"
324 325
    );
    argList::addBoolOption
326 327
    (
        "copyUniform",
328
        "Copy any uniform/ directories too"
329 330 331 332
    );
    argList::addBoolOption
    (
        "fields",
333
        "Use existing geometry decomposition and convert fields only"
334 335
    );
    argList::addBoolOption
336
    (
mattijs's avatar
mattijs committed
337
        "noSets",
338
        "Skip decomposing cellSets, faceSets, pointSets"
339 340
    );
    argList::addBoolOption
341 342
    (
        "force",
343
        "Remove existing processor*/ subdirs before decomposing the geometry"
344 345 346 347
    );
    argList::addBoolOption
    (
        "ifRequired",
348
        "Only decompose geometry if the number of domains has changed"
349
    );
mattijs's avatar
mattijs committed
350

351 352
    // Allow explicit -constant, have zero from time range
    timeSelector::addOptions(true, false);  // constant(true), zero(false)
353

354
    #include "setRootCase.H"
355

356
    const bool dryrun           = args.found("dry-run");
357 358 359
    const bool optRegion        = args.found("region");
    const bool allRegions       = args.found("allRegions");
    const bool writeCellDist    = args.found("cellDist");
360 361 362
    const bool verbose          = args.found("verbose");

    // Most of these are ignored for dry-run (not triggered anywhere)
363 364 365
    const bool copyZero         = args.found("copyZero");
    const bool copyUniform      = args.found("copyUniform");
    const bool decomposeSets    = !args.found("noSets");
366
    const bool decomposeIfRequired = args.found("ifRequired");
367 368 369 370

    bool decomposeFieldsOnly = args.found("fields");
    bool forceOverwrite      = args.found("force");

371

mattijs's avatar
mattijs committed
372
    // Set time from database
373
    #include "createTime.H"
374 375 376 377 378 379 380 381 382 383 384 385

    // Allow override of time (unless dry-run)
    instantList times;
    if (dryrun)
    {
        Info<< "\ndry-run: ignoring -copy*, -fields, -force, time selection"
            << nl;
    }
    else
    {
        times = timeSelector::selectIfPresent(runTime, args);
    }
386

387
    // Allow override of decomposeParDict location
388 389 390 391 392
    fileName decompDictFile(args.get<fileName>("decomposeParDict", ""));
    if (!decompDictFile.empty() && !decompDictFile.isAbsolute())
    {
        decompDictFile = runTime.globalPath()/decompDictFile;
    }
393

394
    // Get all region names
395 396
    wordList regionNames;
    if (allRegions)
397
    {
398
        regionNames = regionProperties(runTime).names();
399

400 401
        Info<< "Decomposing all regions in regionProperties" << nl
            << "    " << flatOutput(regionNames) << nl << endl;
402
    }
403
    else
404
    {
405
        regionNames.resize(1);
406 407
        regionNames.first() =
            args.getOrDefault<word>("region", fvMesh::defaultRegion);
408 409
    }

410
    forAll(regionNames, regioni)
411
    {
412
        const word& regionName = regionNames[regioni];
413
        const word& regionDir =
414
            (regionName == fvMesh::defaultRegion ? word::null : regionName);
415

416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
        if (dryrun)
        {
            Info<< "dry-run: decomposing mesh " << regionName << nl << nl
                << "Create mesh..." << flush;

            domainDecompositionDryRun decompTest
            (
                IOobject
                (
                    regionName,
                    runTime.timeName(),
                    runTime,
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE,
                    false
                ),
                decompDictFile
            );

            decompTest.execute(writeCellDist, verbose);
            continue;
        }

439
        Info<< "\n\nDecomposing mesh " << regionName << nl << endl;
440

441 442
        // Determine the existing processor count directly
        label nProcs = fileHandler().nProcs(runTime.path(), regionDir);
443

444 445 446
        // Get requested numberOfSubdomains directly from the dictionary.
        // Note: have no mesh yet so cannot use decompositionModel::New
        const label nDomains = decompositionMethod::nDomains
447 448
        (
            IOdictionary
449
            (
450
                IOobject::selectIO
451
                (
452 453
                    IOobject
                    (
454
                        decompositionModel::canonicalName,
455
                        runTime.time().system(),
456
                        regionDir,  // region (if non-default)
457
                        runTime,
458
                        IOobject::MUST_READ,
459 460 461 462
                        IOobject::NO_WRITE,
                        false
                    ),
                    decompDictFile
463
                )
464
            )
465
        );
466

467 468 469
        // Give file handler a chance to determine the output directory
        const_cast<fileOperation&>(fileHandler()).setNProcs(nDomains);

470 471 472 473
        if (decomposeFieldsOnly)
        {
            // Sanity check on previously decomposed case
            if (nProcs != nDomains)
474
            {
475
                FatalErrorInFunction
476 477 478 479
                    << "Specified -fields, but the case was decomposed with "
                    << nProcs << " domains"
                    << nl
                    << "instead of " << nDomains
480
                    << " domains as specified in decomposeParDict" << nl
481
                    << exit(FatalError);
482 483
            }
        }
484 485 486
        else if (nProcs)
        {
            bool procDirsProblem = true;
487

488
            if (decomposeIfRequired && nProcs == nDomains)
489
            {
490
                // We can reuse the decomposition
491 492 493
                decomposeFieldsOnly = true;
                procDirsProblem = false;
                forceOverwrite = false;
494

495 496
                Info<< "Using existing processor directories" << nl;
            }
497

498
            if (allRegions || optRegion)
499 500 501 502 503
            {
                procDirsProblem = false;
                forceOverwrite = false;
            }

504 505 506 507
            if (forceOverwrite)
            {
                Info<< "Removing " << nProcs
                    << " existing processor directories" << endl;
508

509 510
                // Remove existing processors directory
                fileNameList dirs
511
                (
512 513 514 515 516
                    fileHandler().readDir
                    (
                        runTime.path(),
                        fileName::Type::DIRECTORY
                    )
517
                );
518
                forAllReverse(dirs, diri)
519
                {
520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544
                    const fileName& d = dirs[diri];

                    // Starts with 'processors'
                    if (d.find("processors") == 0)
                    {
                        if (fileHandler().exists(d))
                        {
                            fileHandler().rmDir(d);
                        }
                    }

                    // Starts with 'processor'
                    if (d.find("processor") == 0)
                    {
                        // Check that integer after processor
                        fileName num(d.substr(9));
                        label proci = -1;
                        if (Foam::read(num.c_str(), proci))
                        {
                            if (fileHandler().exists(d))
                            {
                                fileHandler().rmDir(d);
                            }
                        }
                    }
545
                }
546

547 548
                procDirsProblem = false;
            }
549

550 551
            if (procDirsProblem)
            {
552
                FatalErrorInFunction
553 554 555 556 557 558 559 560 561
                    << "Case is already decomposed with " << nProcs
                    << " domains, use the -force option or manually" << nl
                    << "remove processor directories before decomposing. e.g.,"
                    << nl
                    << "    rm -rf " << runTime.path().c_str() << "/processor*"
                    << nl
                    << exit(FatalError);
            }
        }
562

563 564
        Info<< "Create mesh" << endl;
        domainDecomposition mesh
mattijs's avatar
mattijs committed
565
        (
566 567 568 569
            IOobject
            (
                regionName,
                runTime.timeName(),
570 571 572 573
                runTime,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
574 575
            ),
            decompDictFile
mattijs's avatar
mattijs committed
576
        );
577

578 579
        // Decompose the mesh
        if (!decomposeFieldsOnly)
mattijs's avatar
mattijs committed
580
        {
581 582
            mesh.decomposeMesh();

583
            mesh.writeDecomposition(decomposeSets);
584

585
            if (writeCellDist)
mattijs's avatar
mattijs committed
586
            {
587
                const labelList& procIds = mesh.cellToProc();
588

589
                // Write decomposition as volScalarField for visualization
590
                volScalarField cellDist
mattijs's avatar
mattijs committed
591
                (
592
                    IOobject
mattijs's avatar
mattijs committed
593
                    (
594 595 596 597 598 599 600
                        "cellDist",
                        runTime.timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                    ),
                    mesh,
601 602
                    dimensionedScalar("cellDist", dimless, -1),
                    zeroGradientFvPatchScalarField::typeName
mattijs's avatar
mattijs committed
603
                );
604

605
                forAll(procIds, celli)
606
                {
607 608
                   cellDist[celli] = procIds[celli];
                }
mattijs's avatar
mattijs committed
609

610
                cellDist.correctBoundaryConditions();
611
                cellDist.write();
mattijs's avatar
mattijs committed
612

613
                Info<< nl << "Wrote decomposition as volScalarField to "
614
                    << cellDist.name() << " for visualization."
615
                    << endl;
616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636

                // Write decomposition as labelList for use with 'manual'
                // decomposition method.
                labelIOList cellDecomposition
                (
                    IOobject
                    (
                        "cellDecomposition",
                        mesh.facesInstance(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::NO_WRITE,
                        false
                    ),
                    procIds
                );
                cellDecomposition.write();

                Info<< nl << "Wrote decomposition to "
                    << cellDecomposition.objectPath()
                    << " for use in manual decomposition." << endl;
637
            }
638

639
            fileHandler().flush();
640
        }
mattijs's avatar
mattijs committed
641 642


643
        if (copyZero)
644
        {
645 646
            // Copy the 0 directory into each of the processor directories
            fileName prevTimePath;
647
            for (label proci = 0; proci < mesh.nProcs(); ++proci)
648
            {
649
                Time processorDb
mattijs's avatar
mattijs committed
650
                (
651 652
                    Time::controlDictName,
                    args.rootPath(),
653
                    args.caseName()/("processor" + Foam::name(proci))
mattijs's avatar
mattijs committed
654
                );
655
                processorDb.setTime(runTime);
656

657
                if (fileHandler().isDir(runTime.timePath()))
658
                {
659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679
                    // Get corresponding directory name (to handle processors/)
                    const fileName timePath
                    (
                        fileHandler().objectPath
                        (
                            IOobject
                            (
                                "",
                                processorDb.timeName(),
                                processorDb
                            ),
                            word::null
                        )
                    );

                    if (timePath != prevTimePath)
                    {
                        Info<< "Processor " << proci
                            << ": copying " << runTime.timePath() << nl
                            << " to " << timePath << endl;
                        fileHandler().cp(runTime.timePath(), timePath);
680

681 682
                        prevTimePath = timePath;
                    }
683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703
                }
            }
        }
        else
        {
            // Decompose the field files

            // Cached processor meshes and maps. These are only preserved if
            // running with multiple times.
            PtrList<Time> processorDbList(mesh.nProcs());
            PtrList<fvMesh> procMeshList(mesh.nProcs());
            PtrList<labelIOList> faceProcAddressingList(mesh.nProcs());
            PtrList<labelIOList> cellProcAddressingList(mesh.nProcs());
            PtrList<labelIOList> boundaryProcAddressingList(mesh.nProcs());
            PtrList<fvFieldDecomposer> fieldDecomposerList(mesh.nProcs());
            PtrList<dimFieldDecomposer> dimFieldDecomposerList(mesh.nProcs());
            PtrList<labelIOList> pointProcAddressingList(mesh.nProcs());
            PtrList<pointFieldDecomposer> pointFieldDecomposerList
            (
                mesh.nProcs()
            );
704 705


706 707 708 709 710 711 712 713 714 715 716 717 718 719
            // Loop over all times
            forAll(times, timeI)
            {
                runTime.setTime(times[timeI], timeI);

                Info<< "Time = " << runTime.timeName() << endl;

                // Search for list of objects for this time
                IOobjectList objects(mesh, runTime.timeName());


                // Construct the vol fields
                // ~~~~~~~~~~~~~~~~~~~~~~~~
                PtrList<volScalarField> volScalarFields;
720
                readFields(mesh, objects, volScalarFields, false);
721
                PtrList<volVectorField> volVectorFields;
722
                readFields(mesh, objects, volVectorFields, false);
723
                PtrList<volSphericalTensorField> volSphericalTensorFields;
724
                readFields(mesh, objects, volSphericalTensorFields, false);
725
                PtrList<volSymmTensorField> volSymmTensorFields;
726
                readFields(mesh, objects, volSymmTensorFields, false);
727
                PtrList<volTensorField> volTensorFields;
728
                readFields(mesh, objects, volTensorFields, false);
729 730 731 732 733


                // Construct the dimensioned fields
                // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                PtrList<DimensionedField<scalar, volMesh>> dimScalarFields;
734
                readFields(mesh, objects, dimScalarFields);
735
                PtrList<DimensionedField<vector, volMesh>> dimVectorFields;
736
                readFields(mesh, objects, dimVectorFields);
737 738
                PtrList<DimensionedField<sphericalTensor, volMesh>>
                    dimSphericalTensorFields;
739
                readFields(mesh, objects, dimSphericalTensorFields);
740 741
                PtrList<DimensionedField<symmTensor, volMesh>>
                    dimSymmTensorFields;
742
                readFields(mesh, objects, dimSymmTensorFields);
743
                PtrList<DimensionedField<tensor, volMesh>> dimTensorFields;
744
                readFields(mesh, objects, dimTensorFields);
745 746 747 748 749


                // Construct the surface fields
                // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                PtrList<surfaceScalarField> surfaceScalarFields;
750
                readFields(mesh, objects, surfaceScalarFields, false);
751
                PtrList<surfaceVectorField> surfaceVectorFields;
752
                readFields(mesh, objects, surfaceVectorFields, false);
753 754
                PtrList<surfaceSphericalTensorField>
                    surfaceSphericalTensorFields;
755
                readFields(mesh, objects, surfaceSphericalTensorFields, false);
756
                PtrList<surfaceSymmTensorField> surfaceSymmTensorFields;
757
                readFields(mesh, objects, surfaceSymmTensorFields, false);
758
                PtrList<surfaceTensorField> surfaceTensorFields;
759
                readFields(mesh, objects, surfaceTensorFields, false);
760 761 762 763 764 765 766


                // Construct the point fields
                // ~~~~~~~~~~~~~~~~~~~~~~~~~~
                const pointMesh& pMesh = pointMesh::New(mesh);

                PtrList<pointScalarField> pointScalarFields;
767
                readFields(pMesh, objects, pointScalarFields, false);
768
                PtrList<pointVectorField> pointVectorFields;
769
                readFields(pMesh, objects, pointVectorFields, false);
770
                PtrList<pointSphericalTensorField> pointSphericalTensorFields;
771
                readFields(pMesh, objects, pointSphericalTensorFields, false);
772
                PtrList<pointSymmTensorField> pointSymmTensorFields;
773
                readFields(pMesh, objects, pointSymmTensorFields, false);
774
                PtrList<pointTensorField> pointTensorFields;
775
                readFields(pMesh, objects, pointTensorFields, false);
776 777 778 779 780 781 782


                // Construct the Lagrangian fields
                // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

                fileNameList cloudDirs
                (
783
                    fileHandler().readDir
784
                    (
785 786
                        runTime.timePath()/cloud::prefix,
                        fileName::DIRECTORY
787
                    )
788
                );
789

790 791 792 793 794 795 796 797 798 799
                // Particles
                PtrList<Cloud<indexedParticle>> lagrangianPositions
                (
                    cloudDirs.size()
                );
                // Particles per cell
                PtrList<List<SLList<indexedParticle*>*>> cellParticles
                (
                    cloudDirs.size()
                );
800

801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852
                PtrList<PtrList<labelIOField>> lagrangianLabelFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<labelFieldCompactIOField>>
                lagrangianLabelFieldFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<scalarIOField>> lagrangianScalarFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<scalarFieldCompactIOField>>
                lagrangianScalarFieldFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<vectorIOField>> lagrangianVectorFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<vectorFieldCompactIOField>>
                lagrangianVectorFieldFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<sphericalTensorIOField>>
                lagrangianSphericalTensorFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<sphericalTensorFieldCompactIOField>>
                    lagrangianSphericalTensorFieldFields(cloudDirs.size());
                PtrList<PtrList<symmTensorIOField>> lagrangianSymmTensorFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<symmTensorFieldCompactIOField>>
                lagrangianSymmTensorFieldFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<tensorIOField>> lagrangianTensorFields
                (
                    cloudDirs.size()
                );
                PtrList<PtrList<tensorFieldCompactIOField>>
                lagrangianTensorFieldFields
                (
                    cloudDirs.size()
                );
853

854
                label cloudI = 0;
Mark Olesen's avatar
Mark Olesen committed
855

856
                for (const fileName& cloudDir : cloudDirs)
857
                {
858
                    IOobjectList cloudObjects
859 860 861
                    (
                        mesh,
                        runTime.timeName(),
862
                        cloud::prefix/cloudDir,
863 864 865
                        IOobject::MUST_READ,
                        IOobject::NO_WRITE,
                        false
866
                    );
867

868 869 870 871 872 873
                    // Note: look up "positions" for backwards compatibility
                    if
                    (
                        cloudObjects.found("coordinates")
                     || cloudObjects.found("positions")
                    )
874 875 876
                    {
                        // Read lagrangian particles
                        // ~~~~~~~~~~~~~~~~~~~~~~~~~
877

878
                        Info<< "Identified lagrangian data set: "
879
                            << cloudDir << endl;
880

881 882 883 884 885 886
                        lagrangianPositions.set
                        (
                            cloudI,
                            new Cloud<indexedParticle>
                            (
                                mesh,
887
                                cloudDir,
888 889 890
                                false
                            )
                        );
891 892


893 894
                        // Sort particles per cell
                        // ~~~~~~~~~~~~~~~~~~~~~~~
895

896 897 898 899 900 901 902 903 904
                        cellParticles.set
                        (
                            cloudI,
                            new List<SLList<indexedParticle*>*>
                            (
                                mesh.nCells(),
                                static_cast<SLList<indexedParticle*>*>(nullptr)
                            )
                        );
905

906
                        label i = 0;
907

908
                        for (indexedParticle& p : lagrangianPositions[cloudI])
909
                        {
910
                            p.index() = i++;
911

912
                            label celli = p.cell();
913 914 915 916 917 918 919

                            // Check
                            if (celli < 0 || celli >= mesh.nCells())
                            {
                                FatalErrorInFunction
                                    << "Illegal cell number " << celli
                                    << " for particle with index "
920
                                    << p.index()
921
                                    << " at position "
922
                                    << p.position() << nl
923 924 925 926
                                    << "Cell number should be between 0 and "
                                    << mesh.nCells()-1 << nl
                                    << "On this mesh the particle should"
                                    << " be in cell "
927
                                    << mesh.findCell(p.position())
928 929 930 931 932 933 934