src/postProcessing/functionObjects/field/streamLine/streamLine.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "Pstream.H"
  27 #include "functionObjectList.H"
  28 #include "streamLine.H"
  29 #include "fvMesh.H"
  30 #include "streamLineParticleCloud.H"
  31 #include "ReadFields.H"
  32 #include "meshSearch.H"
  33 #include "sampledSet.H"
  34 #include "globalIndex.H"
  35 #include "mapDistribute.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 defineTypeNameAndDebug(Foam::streamLine, 0);
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 void Foam::streamLine::track()
  44 {
  45     const Time& runTime = obr_.time();
  46     const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
  47
  48     IDLList<streamLineParticle> initialParticles;
  49     streamLineParticleCloud particles
  50     (
  51         mesh,
  52         cloudName_,
  53         initialParticles
  54     );
  55
  56     const sampledSet& seedPoints = sampledSetPtr_();
  57
  58     forAll(seedPoints, i)
  59     {
  60         particles.addParticle
  61         (
  62             new streamLineParticle
  63             (
  64                 mesh,
  65                 seedPoints[i],
  66                 seedPoints.cells()[i],
  67                 lifeTime_               // lifetime
  68             )
  69         );
  70     }
  71
  72     label nSeeds = returnReduce(particles.size(), sumOp<label>());
  73
  74     Info<< "Seeded " << nSeeds << " particles." << endl;
  75
  76     // Read or lookup fields
  77     PtrList<volScalarField> vsFlds;
  78     PtrList<interpolationCellPoint<scalar> > vsInterp;
  79     PtrList<volVectorField> vvFlds;
  80     PtrList<interpolationCellPoint<vector> > vvInterp;
  81
  82     label UIndex = -1;
  83
  84     if (loadFromFiles_)
  85     {
  86         IOobjectList allObjects(mesh, runTime.timeName());
  87
  88         IOobjectList objects(2*fields_.size());
  89         forAll(fields_, i)
  90         {
  91             objects.add(*allObjects[fields_[i]]);
  92         }
  93
  94         ReadFields(mesh, objects, vsFlds);
  95         vsInterp.setSize(vsFlds.size());
  96         forAll(vsFlds, i)
  97         {
  98             vsInterp.set(i, new interpolationCellPoint<scalar>(vsFlds[i]));
  99         }
 100         ReadFields(mesh, objects, vvFlds);
 101         vvInterp.setSize(vvFlds.size());
 102         forAll(vvFlds, i)
 103         {
 104             vvInterp.set(i, new interpolationCellPoint<vector>(vvFlds[i]));
 105         }
 106     }
 107     else
 108     {
 109         label nScalar = 0;
 110         label nVector = 0;
 111
 112         forAll(fields_, i)
 113         {
 114             if (mesh.foundObject<volScalarField>(fields_[i]))
 115             {
 116                 nScalar++;
 117             }
 118             else if (mesh.foundObject<volVectorField>(fields_[i]))
 119             {
 120                 nVector++;
 121             }
 122             else
 123             {
 124                 FatalErrorIn("streamLine::execute()")
 125                     << "Cannot find field " << fields_[i] << endl
 126                     << "Valid scalar fields are:"
 127                     << mesh.names(volScalarField::typeName) << endl
 128                     << "Valid vector fields are:"
 129                     << mesh.names(volVectorField::typeName)
 130                     << exit(FatalError);
 131             }
 132         }
 133         vsInterp.setSize(nScalar);
 134         nScalar = 0;
 135         vvInterp.setSize(nVector);
 136         nVector = 0;
 137
 138         forAll(fields_, i)
 139         {
 140             if (mesh.foundObject<volScalarField>(fields_[i]))
 141             {
 142                 const volScalarField& f = mesh.lookupObject<volScalarField>
 143                 (
 144                     fields_[i]
 145                 );
 146                 vsInterp.set
 147                 (
 148                     nScalar++,
 149                     new interpolationCellPoint<scalar>(f)
 150                 );
 151             }
 152             else if (mesh.foundObject<volVectorField>(fields_[i]))
 153             {
 154                 const volVectorField& f = mesh.lookupObject<volVectorField>
 155                 (
 156                     fields_[i]
 157                 );
 158
 159                 if (f.name() == UName_)
 160                 {
 161                     UIndex = nVector;
 162                 }
 163
 164                 vvInterp.set
 165                 (
 166                     nVector++,
 167                     new interpolationCellPoint<vector>(f)
 168                 );
 169             }
 170         }
 171     }
 172
 173     // Store the names
 174     scalarNames_.setSize(vsInterp.size());
 175     forAll(vsInterp, i)
 176     {
 177         scalarNames_[i] = vsInterp[i].psi().name();
 178     }
 179     vectorNames_.setSize(vvInterp.size());
 180     forAll(vvInterp, i)
 181     {
 182         vectorNames_[i] = vvInterp[i].psi().name();
 183     }
 184
 185     // Check that we know the index of U in the interpolators.
 186
 187     if (UIndex == -1)
 188     {
 189         FatalErrorIn("streamLine::execute()")
 190             << "Cannot find field to move particles with : " << UName_
 191             << endl
 192             << "This field has to be present in the sampled fields "
 193             << fields_
 194             << " and in the objectRegistry." << endl
 195             << exit(FatalError);
 196     }
 197
 198     // Sampled data
 199     // ~~~~~~~~~~~~
 200
 201     // Size to maximum expected sizes.
 202     allTracks_.clear();
 203     allTracks_.setCapacity(nSeeds);
 204     allScalars_.setSize(vsInterp.size());
 205     forAll(allScalars_, i)
 206     {
 207         allScalars_[i].clear();
 208         allScalars_[i].setCapacity(nSeeds);
 209     }
 210     allVectors_.setSize(vvInterp.size());
 211     forAll(allVectors_, i)
 212     {
 213         allVectors_[i].clear();
 214         allVectors_[i].setCapacity(nSeeds);
 215     }
 216
 217     // additional particle info
 218     streamLineParticle::trackingData td
 219     (
 220         particles,
 221         vsInterp,
 222         vvInterp,
 223         UIndex,         // index of U in vvInterp
 224         trackForward_,  // track in +u direction
 225         nSubCycle_,     // step through cells in steps?
 226         allTracks_,
 227         allScalars_,
 228         allVectors_
 229     );
 230
 231
 232     // Set very large dt. Note: cannot use GREAT since 1/GREAT is SMALL
 233     // which is a trigger value for the tracking...
 234     const scalar trackTime = Foam::sqrt(GREAT);
 235
 236     // Track
 237     particles.move(td, trackTime);
 238 }
 239
 240
 241 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 242
 243 Foam::streamLine::streamLine
 244 (
 245     const word& name,
 246     const objectRegistry& obr,
 247     const dictionary& dict,
 248     const bool loadFromFiles
 249 )
 250 :
 251     dict_(dict),
 252     name_(name),
 253     obr_(obr),
 254     loadFromFiles_(loadFromFiles),
 255     active_(true)
 256 {
 257     // Only active if a fvMesh is available
 258     if (isA<fvMesh>(obr_))
 259     {
 260         read(dict_);
 261     }
 262     else
 263     {
 264         active_ = false;
 265         WarningIn
 266         (
 267             "streamLine::streamLine\n"
 268             "(\n"
 269                 "const word&,\n"
 270                 "const objectRegistry&,\n"
 271                 "const dictionary&,\n"
 272                 "const bool\n"
 273             ")"
 274         )   << "No fvMesh available, deactivating."
 275             << nl << endl;
 276     }
 277 }
 278
 279
 280 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 281
 282 Foam::streamLine::~streamLine()
 283 {}
 284
 285
 286 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 287
 288 void Foam::streamLine::read(const dictionary& dict)
 289 {
 290     if (active_)
 291     {
 292         //dict_ = dict;
 293         dict.lookup("fields") >> fields_;
 294         if (dict.found("UName"))
 295         {
 296             dict.lookup("UName") >> UName_;
 297         }
 298         else
 299         {
 300             UName_ = dict.lookupOrDefault<word>("U", "U");
 301         }
 302         dict.lookup("trackForward") >> trackForward_;
 303         dict.lookup("lifeTime") >> lifeTime_;
 304         if (lifeTime_ < 1)
 305         {
 306             FatalErrorIn(":streamLine::read(const dictionary&)")
 307                 << "Illegal value " << lifeTime_ << " for lifeTime"
 308                 << exit(FatalError);
 309         }
 310
 311         dict.lookup("nSubCycle") >> nSubCycle_;
 312         if (nSubCycle_ < 1)
 313         {
 314             nSubCycle_ = 1;
 315         }
 316
 317         cloudName_ = dict.lookupOrDefault<word>("cloudName", "streamLine");
 318         dict.lookup("seedSampleSet") >> seedSet_;
 319
 320         const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
 321
 322         meshSearchPtr_.reset(new meshSearch(mesh, false));
 323
 324         const dictionary& coeffsDict = dict.subDict(seedSet_ + "Coeffs");
 325         sampledSetPtr_ = sampledSet::New
 326         (
 327             seedSet_,
 328             mesh,
 329             meshSearchPtr_(),
 330             coeffsDict
 331         );
 332         coeffsDict.lookup("axis") >> sampledSetAxis_;
 333
 334         scalarFormatterPtr_ = writer<scalar>::New(dict.lookup("setFormat"));
 335         vectorFormatterPtr_ = writer<vector>::New(dict.lookup("setFormat"));
 336     }
 337 }
 338
 339
 340 void Foam::streamLine::execute()
 341 {
 342 //    const Time& runTime = obr_.time();
 343 //    Pout<< "**streamLine::execute : time:" << runTime.timeName() << endl;
 344 //
 345 //    bool isOutputTime = false;
 346 //
 347 //    const functionObjectList& fobs = runTime.functionObjects();
 348 //
 349 //    forAll(fobs, i)
 350 //    {
 351 //        if (isA<streamLineFunctionObject>(fobs[i]))
 352 //        {
 353 //            const streamLineFunctionObject& fo =
 354 //                dynamic_cast<const streamLineFunctionObject&>(fobs[i]);
 355 //
 356 //            if (fo.name() == name_)
 357 //            {
 358 //                Pout<< "found me:" << i << endl;
 359 //                if (fo.outputControl().output())
 360 //                {
 361 //                    isOutputTime = true;
 362 //                    break;
 363 //                }
 364 //            }
 365 //        }
 366 //    }
 367 //
 368 //
 369 //    if (active_ && isOutputTime)
 370 //    {
 371 //        track();
 372 //    }
 373 }
 374
 375
 376 void Foam::streamLine::end()
 377 {}
 378
 379
 380 void Foam::streamLine::write()
 381 {
 382     if (active_)
 383     {
 384         const Time& runTime = obr_.time();
 385         const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
 386
 387
 388         // Do all injection and tracking
 389         track();
 390
 391
 392         if (Pstream::parRun())
 393         {
 394             // Append slave tracks to master ones
 395             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 396
 397             globalIndex globalTrackIDs(allTracks_.size());
 398
 399             // Construct a distribution map to pull all to the master.
 400             labelListList sendMap(Pstream::nProcs());
 401             labelListList recvMap(Pstream::nProcs());
 402
 403             if (Pstream::master())
 404             {
 405                 // Master: receive all. My own first, then consecutive
 406                 // processors.
 407                 label trackI = 0;
 408
 409                 forAll(recvMap, procI)
 410                 {
 411                     labelList& fromProc = recvMap[procI];
 412                     fromProc.setSize(globalTrackIDs.localSize(procI));
 413                     forAll(fromProc, i)
 414                     {
 415                         fromProc[i] = trackI++;
 416                     }
 417                 }
 418             }
 419
 420             labelList& toMaster = sendMap[0];
 421             toMaster.setSize(globalTrackIDs.localSize());
 422             forAll(toMaster, i)
 423             {
 424                 toMaster[i] = i;
 425             }
 426
 427             const mapDistribute distMap
 428             (
 429                 globalTrackIDs.size(),
 430                 sendMap.xfer(),
 431                 recvMap.xfer()
 432             );
 433
 434
 435             // Distribute the track positions. Note: use scheduled comms
 436             // to prevent buffering.
 437             mapDistribute::distribute
 438             (
 439                 Pstream::scheduled,
 440                 distMap.schedule(),
 441                 distMap.constructSize(),
 442                 distMap.subMap(),
 443                 distMap.constructMap(),
 444                 allTracks_
 445             );
 446
 447             // Distribute the scalars
 448             forAll(allScalars_, scalarI)
 449             {
 450                 mapDistribute::distribute
 451                 (
 452                     Pstream::scheduled,
 453                     distMap.schedule(),
 454                     distMap.constructSize(),
 455                     distMap.subMap(),
 456                     distMap.constructMap(),
 457                     allScalars_[scalarI]
 458                 );
 459             }
 460             // Distribute the vectors
 461             forAll(allVectors_, vectorI)
 462             {
 463                 mapDistribute::distribute
 464                 (
 465                     Pstream::scheduled,
 466                     distMap.schedule(),
 467                     distMap.constructSize(),
 468                     distMap.subMap(),
 469                     distMap.constructMap(),
 470                     allVectors_[vectorI]
 471                 );
 472             }
 473         }
 474
 475
 476         label n = 0;
 477         forAll(allTracks_, trackI)
 478         {
 479             n += allTracks_[trackI].size();
 480         }
 481
 482         Info<< "Tracks:" << allTracks_.size()
 483             << "  total samples:" << n << endl;
 484
 485
 486         // Massage into form suitable for writers
 487         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 488
 489         if (Pstream::master() && allTracks_.size())
 490         {
 491             // Make output directory
 492
 493             fileName vtkPath
 494             (
 495                 Pstream::parRun()
 496               ? runTime.path()/".."/"sets"/name()
 497               : runTime.path()/"sets"/name()
 498             );
 499             if (mesh.name() != fvMesh::defaultRegion)
 500             {
 501                 vtkPath = vtkPath/mesh.name();
 502             }
 503             vtkPath = vtkPath/mesh.time().timeName();
 504
 505             mkDir(vtkPath);
 506
 507             // Convert track positions
 508
 509             PtrList<coordSet> tracks(allTracks_.size());
 510             forAll(allTracks_, trackI)
 511             {
 512                 tracks.set
 513                 (
 514                     trackI,
 515                     new coordSet
 516                     (
 517                         "track" + Foam::name(trackI),
 518                         sampledSetAxis_                 //"xyz"
 519                     )
 520                 );
 521                 tracks[trackI].transfer(allTracks_[trackI]);
 522             }
 523
 524             // Convert scalar values
 525
 526             if (allScalars_.size() > 0)
 527             {
 528                 List<List<scalarField> > scalarValues(allScalars_.size());
 529
 530                 forAll(allScalars_, scalarI)
 531                 {
 532                     DynamicList<scalarList>& allTrackVals =
 533                         allScalars_[scalarI];
 534                     scalarValues[scalarI].setSize(allTrackVals.size());
 535
 536                     forAll(allTrackVals, trackI)
 537                     {
 538                         scalarList& trackVals = allTrackVals[trackI];
 539                         scalarValues[scalarI][trackI].transfer(trackVals);
 540                     }
 541                 }
 542
 543                 fileName vtkFile
 544                 (
 545                     vtkPath
 546                   / scalarFormatterPtr_().getFileName
 547                     (
 548                         tracks[0],
 549                         scalarNames_
 550                     )
 551                 );
 552
 553                 Info<< "Writing data to " << vtkFile.path() << endl;
 554
 555                 scalarFormatterPtr_().write
 556                 (
 557                     true,           // writeTracks
 558                     tracks,
 559                     scalarNames_,
 560                     scalarValues,
 561                     OFstream(vtkFile)()
 562                 );
 563             }
 564
 565             // Convert vector values
 566
 567             if (allVectors_.size() > 0)
 568             {
 569                 List<List<vectorField> > vectorValues(allVectors_.size());
 570
 571                 forAll(allVectors_, vectorI)
 572                 {
 573                     DynamicList<vectorList>& allTrackVals =
 574                         allVectors_[vectorI];
 575                     vectorValues[vectorI].setSize(allTrackVals.size());
 576
 577                     forAll(allTrackVals, trackI)
 578                     {
 579                         vectorList& trackVals = allTrackVals[trackI];
 580                         vectorValues[vectorI][trackI].transfer(trackVals);
 581                     }
 582                 }
 583
 584                 fileName vtkFile
 585                 (
 586                     vtkPath
 587                   / vectorFormatterPtr_().getFileName
 588                     (
 589                         tracks[0],
 590                         vectorNames_
 591                     )
 592                 );
 593
 594                 //Info<< "Writing vector data to " << vtkFile << endl;
 595
 596                 vectorFormatterPtr_().write
 597                 (
 598                     true,           // writeTracks
 599                     tracks,
 600                     vectorNames_,
 601                     vectorValues,
 602                     OFstream(vtkFile)()
 603                 );
 604             }
 605         }
 606     }
 607 }
 608
 609
 610 void Foam::streamLine::updateMesh(const mapPolyMesh&)
 611 {
 612     read(dict_);
 613 }
 614
 615
 616 void Foam::streamLine::movePoints(const pointField&)
 617 {
 618     // Moving mesh affects the search tree
 619     read(dict_);
 620 }
 621
 622
 623 //void Foam::streamLine::readUpdate(const polyMesh::readUpdateState state)
 624 //{
 625 //    if (state != UNCHANGED)
 626 //    {
 627 //        read(dict_);
 628 //    }
 629 //}
 630
 631
 632 // ************************************************************************* //