applications/solvers/surfaceTracking/freeSurface/makeFreeSurfaceData.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "freeSurface.H"
  27 #include "primitivePatchInterpolation.H"
  28 #include "wedgeFaPatch.H"
  29 #include "wallFvPatch.H"
  30 #include "wedgeFaPatchFields.H"
  31 #include "slipFaPatchFields.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  39
  40 void freeSurface::makeInterpolators()
  41 {
  42     if (debug)
  43     {
  44         Info<< "freeSurface::makeInterpolators() : "
  45             << "making pathc to patch interpolator"
  46             << endl;
  47     }
  48
  49
  50     // It is an error to attempt to recalculate
  51     // if the pointer is already set
  52     if
  53     (
  54         interpolatorBAPtr_ ||
  55         interpolatorABPtr_
  56     )
  57     {
  58         FatalErrorIn("freeSurface::makeInterpolators()")
  59             << "patch to patch interpolators already exists"
  60                 << abort(FatalError);
  61     }
  62
  63
  64     if(aPatchID() == -1)
  65     {
  66         FatalErrorIn("freeSurface::makeInterpolators()")
  67             << "Free surface patch A not defined."
  68             << abort(FatalError);
  69     }
  70
  71
  72     if(bPatchID() == -1)
  73     {
  74         FatalErrorIn("freeSurface::makeInterpolators()")
  75             << "Free surface patch B not defined."
  76             << abort(FatalError);
  77     }
  78
  79 //     patchToPatchInterpolation::setDirectHitTol(1e-2);
  80
  81     interpolatorBAPtr_ = new IOpatchToPatchInterpolation
  82     (
  83         IOobject
  84         (
  85             "baInterpolator",
  86             DB().timeName(),
  87             mesh(),
  88             IOobject::READ_IF_PRESENT,
  89             IOobject::AUTO_WRITE
  90         ),
  91         mesh().boundaryMesh()[bPatchID()],
  92         mesh().boundaryMesh()[aPatchID()],
  93         intersection::VISIBLE
  94         // intersection::HALF_RAY
  95     );
  96
  97
  98     const scalarField& faceDistBA =
  99         interpolatorBAPtr_->faceDistanceToIntersection();
 100
 101     forAll(faceDistBA, faceI)
 102     {
 103         if(mag(faceDistBA[faceI] - GREAT) < SMALL)
 104         {
 105             FatalErrorIn("freeSurface::makeInterpolators()")
 106                 << "Error in B-to-A face patchToPatchInterpolation."
 107                 << abort(FatalError);
 108         }
 109     }
 110
 111     const scalarField& pointDistBA =
 112         interpolatorBAPtr_->pointDistanceToIntersection();
 113
 114     forAll(pointDistBA, pointI)
 115     {
 116         if(mag(pointDistBA[pointI] - GREAT) < SMALL)
 117         {
 118             FatalErrorIn("freeSurface::makeInterpolators()")
 119                 << "Error in B-to-A point patchToPatchInterpolation."
 120                 << abort(FatalError);
 121         }
 122     }
 123
 124
 125     interpolatorABPtr_ = new IOpatchToPatchInterpolation
 126     (
 127         IOobject
 128         (
 129             "abInterpolator",
 130             DB().timeName(),
 131             mesh(),
 132             IOobject::READ_IF_PRESENT,
 133             IOobject::AUTO_WRITE
 134         ),
 135         mesh().boundaryMesh()[aPatchID()],
 136         mesh().boundaryMesh()[bPatchID()],
 137         intersection::VISIBLE
 138         // intersection::HALF_RAY
 139     );
 140
 141
 142     const scalarField& faceDistAB =
 143         interpolatorABPtr_->faceDistanceToIntersection();
 144
 145     forAll(faceDistAB, faceI)
 146     {
 147         if(mag(faceDistAB[faceI] - GREAT) < SMALL)
 148         {
 149             FatalErrorIn("freeSurface::makeInterpolators()")
 150                 << "Error in A-to-B face patchToPatchInterpolation."
 151                 << abort(FatalError);
 152         }
 153     }
 154
 155     const scalarField& pointDistAB =
 156         interpolatorABPtr_->pointDistanceToIntersection();
 157
 158     forAll(pointDistAB, pointI)
 159     {
 160         if(mag(pointDistAB[pointI] - GREAT)<SMALL)
 161         {
 162             FatalErrorIn("freeSurface::makeInterpolators()")
 163                 << "Error in A-to-B point patchToPatchInterpolation."
 164                 << abort(FatalError);
 165         }
 166     }
 167
 168
 169     Info << "\nCheck A-to-B and B-to-A interpolators" << endl;
 170
 171     scalar maxDist = max
 172     (
 173         mag
 174         (
 175             interpolatorABPtr_->faceInterpolate
 176             (
 177                 vectorField(mesh().boundaryMesh()[aPatchID()]
 178                .faceCentres())
 179             )
 180           - mesh().boundaryMesh()[bPatchID()].faceCentres()
 181         )
 182     );
 183
 184     scalar maxDistPt = max
 185     (
 186         mag
 187         (
 188             interpolatorABPtr_->pointInterpolate
 189             (
 190                 vectorField(mesh().boundaryMesh()[aPatchID()]
 191                .localPoints())
 192             )
 193           - mesh().boundaryMesh()[bPatchID()].localPoints()
 194         )
 195     );
 196
 197     Info << "A-to-B interpolation error, face: " << maxDist
 198         << ", point: " << maxDistPt << endl;
 199
 200
 201     maxDist = max
 202     (
 203         mag
 204         (
 205             interpolatorBAPtr_->faceInterpolate
 206             (
 207                 vectorField
 208                 (
 209                     mesh().boundaryMesh()[bPatchID()].faceCentres()
 210                 )
 211             )
 212           - mesh().boundaryMesh()[aPatchID()].faceCentres()
 213         )
 214     );
 215
 216     maxDistPt = max
 217     (
 218         mag
 219         (
 220             interpolatorBAPtr_->pointInterpolate
 221             (
 222                 vectorField
 223                 (
 224                     mesh().boundaryMesh()[bPatchID()].localPoints()
 225                 )
 226             )
 227           - mesh().boundaryMesh()[aPatchID()].localPoints()
 228         )
 229     );
 230
 231     Info << "B-to-A interpolation error, face: " << maxDist
 232         << ", point: " << maxDistPt << endl;
 233 }
 234
 235
 236 void freeSurface::makeControlPoints()
 237 {
 238     if (debug)
 239     {
 240         Info<< "freeSurface::makeControlPoints() : "
 241             << "making control points"
 242             << endl;
 243     }
 244
 245
 246     // It is an error to attempt to recalculate
 247     // if the pointer is already set
 248     if (controlPointsPtr_)
 249     {
 250         FatalErrorIn("freeSurface::makeInterpolators()")
 251             << "patch to patch interpolators already exists"
 252             << abort(FatalError);
 253     }
 254
 255     IOobject controlPointsHeader
 256     (
 257         "controlPoints",
 258         DB().timeName(),
 259         mesh(),
 260         IOobject::MUST_READ
 261     );
 262
 263     if (controlPointsHeader.headerOk())
 264     {
 265         controlPointsPtr_ =
 266             new vectorIOField
 267             (
 268                 IOobject
 269                 (
 270                     "controlPoints",
 271                     DB().timeName(),
 272                     mesh(),
 273                     IOobject::MUST_READ,
 274                     IOobject::AUTO_WRITE
 275                 )
 276             );
 277     }
 278     else
 279     {
 280         controlPointsPtr_ =
 281             new vectorIOField
 282             (
 283                 IOobject
 284                 (
 285                     "controlPoints",
 286                     DB().timeName(),
 287                     mesh(),
 288                     IOobject::NO_READ,
 289                     IOobject::AUTO_WRITE
 290                 ),
 291                 aMesh().areaCentres().internalField()
 292             );
 293
 294         initializeControlPointsPosition();
 295     }
 296 }
 297
 298
 299 void freeSurface::makeMotionPointsMask()
 300 {
 301     if (debug)
 302     {
 303         Info<< "freeSurface::makeMotionPointsMask() : "
 304             << "making motion points mask"
 305             << endl;
 306     }
 307
 308
 309     // It is an error to attempt to recalculate
 310     // if the pointer is already set
 311     if (motionPointsMaskPtr_)
 312     {
 313         FatalErrorIn("freeSurface::motionPointsMask()")
 314             << "motion points mask already exists"
 315             << abort(FatalError);
 316     }
 317
 318
 319     if(aPatchID() == -1)
 320     {
 321         FatalErrorIn("freeSurface::makeMotionPointsMask()")
 322             << "Free surface patch A not defined."
 323             << abort(FatalError);
 324     }
 325
 326
 327     motionPointsMaskPtr_ = new labelList
 328     (
 329         mesh().boundaryMesh()[aPatchID()].nPoints(),
 330         1
 331     );
 332 }
 333
 334
 335 void freeSurface::makeDirections()
 336 {
 337     if (debug)
 338     {
 339         Info<< "freeSurface::makeDirections() : "
 340             << "making displacement directions for points and "
 341             << "control points"
 342             << endl;
 343     }
 344
 345
 346     // It is an error to attempt to recalculate
 347     // if the pointer is already set
 348     if
 349     (
 350         pointsDisplacementDirPtr_ ||
 351         facesDisplacementDirPtr_
 352     )
 353     {
 354         FatalErrorIn("freeSurface::makeDirections()")
 355             << "points and control points displacement directions "
 356             << "already exists"
 357             << abort(FatalError);
 358     }
 359
 360
 361     if(aPatchID() == -1)
 362     {
 363         FatalErrorIn("freeSurface::makeDirections()")
 364             << "Free surface patch A not defined."
 365             << abort(FatalError);
 366     }
 367
 368
 369     pointsDisplacementDirPtr_ =
 370         new vectorField
 371         (
 372             mesh().boundaryMesh()[aPatchID()].nPoints(),
 373             vector::zero
 374         );
 375
 376     facesDisplacementDirPtr_ =
 377         new vectorField
 378         (
 379             mesh().boundaryMesh()[aPatchID()].size(),
 380             vector::zero
 381         );
 382
 383     if(!normalMotionDir())
 384     {
 385         if(mag(motionDir_) < SMALL)
 386         {
 387             FatalErrorIn("freeSurface::makeDirections()")
 388                 << "Zero motion direction"
 389                     << abort(FatalError);
 390         }
 391
 392         facesDisplacementDir() = motionDir_;
 393         pointsDisplacementDir() = motionDir_;
 394     }
 395
 396     updateDisplacementDirections();
 397 }
 398
 399
 400 void freeSurface::makeTotalDisplacement()
 401 {
 402     if (debug)
 403     {
 404         Info<< "freeSurface::makeTotalDisplacement() : "
 405             << "making zero total points displacement"
 406             << endl;
 407     }
 408
 409
 410     // It is an error to attempt to recalculate
 411     // if the pointer is already set
 412     if (totalDisplacementPtr_)
 413     {
 414         FatalErrorIn("freeSurface::makeTotalDisplacement()")
 415             << "total points displacement already exists"
 416             << abort(FatalError);
 417     }
 418
 419     totalDisplacementPtr_ =
 420         new vectorIOField
 421         (
 422             IOobject
 423             (
 424                 "totalDisplacement",
 425                 DB().timeName(),
 426                 mesh(),
 427                 IOobject::NO_READ,
 428                 IOobject::AUTO_WRITE
 429             ),
 430             vectorField
 431             (
 432                 mesh().boundaryMesh()[aPatchID()].nPoints(),
 433                 vector::zero
 434             )
 435         );
 436 }
 437
 438
 439 void freeSurface::readTotalDisplacement()
 440 {
 441     if (debug)
 442     {
 443         Info<< "freeSurface::readTotalDisplacement() : "
 444             << "reading total points displacement if present"
 445             << endl;
 446     }
 447
 448
 449     // It is an error to attempt to recalculate
 450     // if the pointer is already set
 451     if (totalDisplacementPtr_)
 452     {
 453         FatalErrorIn("freeSurface::makeTotalDisplacement()")
 454             << "total points displacement already exists"
 455             << abort(FatalError);
 456     }
 457
 458     if
 459     (
 460         IOobject
 461         (
 462             "totalDisplacement",
 463             DB().timeName(),
 464             mesh(),
 465             IOobject::MUST_READ,
 466             IOobject::AUTO_WRITE
 467         ).headerOk()
 468     )
 469     {
 470         totalDisplacementPtr_ =
 471             new vectorIOField
 472             (
 473                 IOobject
 474                 (
 475                     "totalDisplacement",
 476                     DB().timeName(),
 477                     mesh(),
 478                     IOobject::MUST_READ,
 479                     IOobject::AUTO_WRITE
 480                 )
 481             );
 482     }
 483 }
 484
 485
 486 void freeSurface::makeFaMesh() const
 487 {
 488     if (debug)
 489     {
 490         Info<< "freeSurface::makeFaMesh() : "
 491             << "making finite area mesh"
 492             << endl;
 493     }
 494
 495     // It is an error to attempt to recalculate
 496     // if the pointer is already set
 497     if (aMeshPtr_)
 498     {
 499         FatalErrorIn("freeSurface::makeFaMesh()")
 500             << "finite area mesh already exists"
 501             << abort(FatalError);
 502     }
 503
 504     aMeshPtr_ = new faMesh(mesh());
 505 }
 506
 507 void freeSurface::makeUs() const
 508 {
 509     if (debug)
 510     {
 511         Info<< "freeSurface::makeUs() : "
 512             << "making free-surface velocity field"
 513             << endl;
 514     }
 515
 516
 517     // It is an error to attempt to recalculate
 518     // if the pointer is already set
 519     if (UsPtr_)
 520     {
 521         FatalErrorIn("freeSurface::makeUs()")
 522             << "free-surface velocity field already exists"
 523             << abort(FatalError);
 524     }
 525
 526
 527     wordList patchFieldTypes
 528     (
 529         aMesh().boundary().size(),
 530         zeroGradientFaPatchVectorField::typeName
 531     );
 532
 533     forAll(aMesh().boundary(), patchI)
 534     {
 535         if
 536         (
 537             aMesh().boundary()[patchI].type()
 538          == wedgeFaPatch::typeName
 539         )
 540         {
 541             patchFieldTypes[patchI] =
 542                 wedgeFaPatchVectorField::typeName;
 543         }
 544         else
 545         {
 546             label ngbPolyPatchID =
 547                 aMesh().boundary()[patchI].ngbPolyPatchIndex();
 548
 549             if (ngbPolyPatchID != -1)
 550             {
 551                 if
 552                 (
 553                     mesh().boundary()[ngbPolyPatchID].type()
 554                  == wallFvPatch::typeName
 555                 )
 556                 {
 557                     patchFieldTypes[patchI] =
 558                         slipFaPatchVectorField::typeName;
 559                 }
 560             }
 561         }
 562     }
 563
 564
 565     UsPtr_ = new areaVectorField
 566     (
 567         IOobject
 568         (
 569             "Us",
 570             DB().timeName(),
 571             mesh(),
 572             IOobject::NO_READ,
 573             IOobject::NO_WRITE
 574         ),
 575         aMesh(),
 576         dimensioned<vector>("Us", dimVelocity, vector::zero),
 577         patchFieldTypes
 578     );
 579 }
 580
 581
 582 void freeSurface::makePhis()
 583 {
 584     if (debug)
 585     {
 586         Info<< "freeSurface::makePhis() : "
 587             << "making free-surface fluid flux"
 588             << endl;
 589     }
 590
 591
 592     // It is an error to attempt to recalculate
 593     // if the pointer is already set
 594     if (phisPtr_)
 595     {
 596         FatalErrorIn("freeSurface::makePhis()")
 597             << "free-surface fluid flux already exists"
 598             << abort(FatalError);
 599     }
 600
 601
 602     phisPtr_ = new edgeScalarField
 603     (
 604         IOobject
 605         (
 606             "phis",
 607             DB().timeName(),
 608             mesh(),
 609             IOobject::NO_READ,
 610             IOobject::NO_WRITE
 611         ),
 612         linearEdgeInterpolate(Us()) & aMesh().Le()
 613     );
 614 }
 615
 616
 617 void freeSurface::makeSurfactConc() const
 618 {
 619     if (debug)
 620     {
 621         Info<< "freeSurface::makeSurfactConc() : "
 622             << "making free-surface surfactant concentration field"
 623             << endl;
 624     }
 625
 626
 627     // It is an error to attempt to recalculate
 628     // if the pointer is already set
 629     if (surfactConcPtr_)
 630     {
 631         FatalErrorIn("freeSurface::makeSurfaceConc()")
 632             << "free-surface surfactant concentratio field already exists"
 633             << abort(FatalError);
 634     }
 635
 636     surfactConcPtr_ = new areaScalarField
 637     (
 638         IOobject
 639         (
 640             "Cs",
 641             DB().timeName(),
 642             mesh(),
 643             IOobject::MUST_READ,
 644             IOobject::AUTO_WRITE
 645         ),
 646         aMesh()
 647     );
 648 }
 649
 650
 651 void freeSurface::makeSurfaceTension() const
 652 {
 653     if (debug)
 654     {
 655         Info<< "freeSurface::makeSurfaceTension() : "
 656             << "making surface tension field"
 657             << endl;
 658     }
 659
 660
 661     // It is an error to attempt to recalculate
 662     // if the pointer is already set
 663     if (surfaceTensionPtr_)
 664     {
 665         FatalErrorIn("freeSurface::makeSurfaceTension()")
 666             << "surface tension field already exists"
 667             << abort(FatalError);
 668     }
 669
 670
 671     surfaceTensionPtr_ = new areaScalarField
 672     (
 673         IOobject
 674         (
 675             "surfaceTension",
 676             DB().timeName(),
 677             mesh(),
 678             IOobject::NO_READ,
 679             IOobject::NO_WRITE
 680         ),
 681         cleanInterfaceSurfTension()
 682       + surfactant().surfactR()*
 683         surfactant().surfactT()*
 684         surfactant().surfactSaturatedConc()*
 685         log(1.0 - surfactantConcentration()/
 686         surfactant().surfactSaturatedConc())
 687     );
 688 }
 689
 690
 691 void freeSurface::makeSurfactant() const
 692 {
 693     if (debug)
 694     {
 695         Info<< "freeSurface::makeSurfactant() : "
 696             << "making surfactant properties"
 697             << endl;
 698     }
 699
 700
 701     // It is an error to attempt to recalculate
 702     // if the pointer is already set
 703     if (surfactantPtr_)
 704     {
 705         FatalErrorIn("freeSurface::makeSurfactant()")
 706             << "surfactant properties already exists"
 707             << abort(FatalError);
 708     }
 709
 710
 711     const dictionary& surfactProp =
 712         this->subDict("surfactantProperties");
 713
 714     surfactantPtr_ = new surfactantProperties(surfactProp);
 715 }
 716
 717
 718 void freeSurface::makeFluidIndicator()
 719 {
 720     if (debug)
 721     {
 722         Info<< "freeSurface::makeFluidIndicator() : "
 723             << "making fluid indicator"
 724             << endl;
 725     }
 726
 727
 728     // It is an error to attempt to recalculate
 729     // if the pointer is already set
 730     if (fluidIndicatorPtr_)
 731     {
 732         FatalErrorIn("freeSurface::makeFluidIndicator()")
 733             << "fluid indicator already exists"
 734             << abort(FatalError);
 735     }
 736
 737     fluidIndicatorPtr_ = new volScalarField
 738     (
 739         IOobject
 740         (
 741             "fluidIndicator",
 742             DB().timeName(),
 743             mesh(),
 744             IOobject::NO_READ,
 745             IOobject::NO_WRITE
 746         ),
 747         mesh(),
 748         dimensionedScalar("1", dimless, 1.0),
 749         zeroGradientFvPatchScalarField::typeName
 750     );
 751
 752     volScalarField& fluidIndicator = *fluidIndicatorPtr_;
 753
 754     if (twoFluids())
 755     {
 756         // find start cell
 757         label pointOnShadowPatch =
 758             mesh().boundaryMesh()[bPatchID()][0][0];
 759
 760         label startCell = mesh().pointCells()[pointOnShadowPatch][0];
 761
 762
 763         // get cell-cells addressing
 764         const labelListList& cellCells = mesh().cellCells();
 765
 766         SLList<label> slList(startCell);
 767
 768         while (slList.size())
 769         {
 770             label curCell = slList.removeHead();
 771
 772             if (fluidIndicator[curCell] == 1)
 773             {
 774                 fluidIndicator[curCell] = 0.0;
 775
 776                 for (int i = 0; i < cellCells[curCell].size(); i++)
 777                 {
 778                     slList.append(cellCells[curCell][i]);
 779                 }
 780             }
 781         }
 782     }
 783
 784     fluidIndicator.correctBoundaryConditions();
 785 }
 786
 787
 788 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 789
 790 const IOpatchToPatchInterpolation& freeSurface::interpolatorAB()
 791 {
 792     if (!interpolatorABPtr_)
 793     {
 794         makeInterpolators();
 795     }
 796
 797     return *interpolatorABPtr_;
 798 }
 799
 800
 801 const IOpatchToPatchInterpolation& freeSurface::interpolatorBA()
 802 {
 803     if (!interpolatorBAPtr_)
 804     {
 805         makeInterpolators();
 806     }
 807
 808     return *interpolatorBAPtr_;
 809 }
 810
 811
 812 vectorField& freeSurface::controlPoints()
 813 {
 814     if (!controlPointsPtr_)
 815     {
 816         makeControlPoints();
 817     }
 818
 819     return *controlPointsPtr_;
 820 }
 821
 822
 823 labelList& freeSurface::motionPointsMask()
 824 {
 825     if (!motionPointsMaskPtr_)
 826     {
 827         makeMotionPointsMask();
 828     }
 829
 830     return *motionPointsMaskPtr_;
 831 }
 832
 833
 834 vectorField& freeSurface::pointsDisplacementDir()
 835 {
 836     if (!pointsDisplacementDirPtr_)
 837     {
 838         makeDirections();
 839     }
 840
 841     return *pointsDisplacementDirPtr_;
 842 }
 843
 844
 845 vectorField& freeSurface::facesDisplacementDir()
 846 {
 847     if (!facesDisplacementDirPtr_)
 848     {
 849         makeDirections();
 850     }
 851
 852     return *facesDisplacementDirPtr_;
 853 }
 854
 855
 856 vectorField& freeSurface::totalDisplacement()
 857 {
 858     if (!totalDisplacementPtr_)
 859     {
 860         makeTotalDisplacement();
 861     }
 862
 863     return *totalDisplacementPtr_;
 864 }
 865
 866
 867 faMesh& freeSurface::aMesh()
 868 {
 869     if (!aMeshPtr_)
 870     {
 871         makeFaMesh();
 872     }
 873
 874     return *aMeshPtr_;
 875 }
 876
 877 const faMesh& freeSurface::aMesh() const
 878 {
 879     if (!aMeshPtr_)
 880     {
 881         makeFaMesh();
 882     }
 883
 884     return *aMeshPtr_;
 885 }
 886
 887 areaVectorField& freeSurface::Us()
 888 {
 889     if (!UsPtr_)
 890     {
 891         makeUs();
 892     }
 893
 894     return *UsPtr_;
 895 }
 896
 897 const areaVectorField& freeSurface::Us() const
 898 {
 899     if (!UsPtr_)
 900     {
 901         makeUs();
 902     }
 903
 904     return *UsPtr_;
 905 }
 906
 907 edgeScalarField& freeSurface::Phis()
 908 {
 909     if (!phisPtr_)
 910     {
 911         makePhis();
 912     }
 913
 914     return *phisPtr_;
 915 }
 916
 917 areaScalarField& freeSurface::surfactantConcentration()
 918 {
 919     if (!surfactConcPtr_)
 920     {
 921         makeSurfactConc();
 922     }
 923
 924     return *surfactConcPtr_;
 925 }
 926
 927 const areaScalarField& freeSurface::surfactantConcentration() const
 928 {
 929     if (!surfactConcPtr_)
 930     {
 931         makeSurfactConc();
 932     }
 933
 934     return *surfactConcPtr_;
 935 }
 936
 937 areaScalarField& freeSurface::surfaceTension()
 938 {
 939     if (!surfaceTensionPtr_)
 940     {
 941         makeSurfaceTension();
 942     }
 943
 944     return *surfaceTensionPtr_;
 945 }
 946
 947 const areaScalarField& freeSurface::surfaceTension() const
 948 {
 949     if (!surfaceTensionPtr_)
 950     {
 951         makeSurfaceTension();
 952     }
 953
 954     return *surfaceTensionPtr_;
 955 }
 956
 957 const surfactantProperties& freeSurface::surfactant() const
 958 {
 959     if (!surfactantPtr_)
 960     {
 961         makeSurfactant();
 962     }
 963
 964     return *surfactantPtr_;
 965 }
 966
 967
 968 const volScalarField& freeSurface::fluidIndicator()
 969 {
 970     if (!fluidIndicatorPtr_)
 971     {
 972         makeFluidIndicator();
 973     }
 974
 975     return *fluidIndicatorPtr_;
 976 }
 977
 978
 979 tmp<areaVectorField> freeSurface::surfaceTensionGrad()
 980 {
 981     tmp<areaVectorField> tgrad
 982     (
 983         new areaVectorField
 984         (
 985             IOobject
 986             (
 987                 "surfaceTensionGrad",
 988                 DB().timeName(),
 989                 mesh(),
 990                 IOobject::NO_READ,
 991                 IOobject::NO_WRITE
 992             ),
 993             (-fac::grad(surfactantConcentration())*
 994             surfactant().surfactR()*surfactant().surfactT()/
 995             (1.0 - surfactantConcentration()/
 996             surfactant().surfactSaturatedConc()))()
 997         )
 998     );
 999
1000     return tgrad;
1001 }
1002
1003
1004 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1005
1006 } // End namespace Foam
1007
1008 // ************************************************************************* //