applications/solvers/surfaceTracking/freeSurface/freeSurface.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 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "freeSurface.H"
  29
  30 #include "volFields.H"
  31 #include "transformField.H"
  32
  33 #include "emptyFaPatch.H"
  34 #include "wedgeFaPatch.H"
  35 #include "wallFvPatch.H"
  36
  37 #include "EulerDdtScheme.H"
  38 #include "CrankNicolsonDdtScheme.H"
  39 #include "backwardDdtScheme.H"
  40
  41 #include "tetFemMatrices.H"
  42 #include "tetPointFields.H"
  43 #include "faceTetPolyPatch.H"
  44 #include "tetPolyPatchInterpolation.H"
  45 #include "fixedValueTetPolyPatchFields.H"
  46 #include "fixedValuePointPatchFields.H"
  47 #include "twoDPointCorrector.H"
  48
  49 #include "slipFvPatchFields.H"
  50 #include "symmetryFvPatchFields.H"
  51 #include "fixedGradientFvPatchFields.H"
  52 #include "zeroGradientCorrectedFvPatchFields.H"
  53 #include "fixedGradientCorrectedFvPatchFields.H"
  54 #include "fixedValueCorrectedFvPatchFields.H"
  55
  56 #include "primitivePatchInterpolation.H"
  57
  58 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  59
  60 namespace Foam
  61 {
  62
  63 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  64
  65 defineTypeNameAndDebug(freeSurface, 0);
  66
  67
  68 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  69
  70 void freeSurface::clearOut()
  71 {
  72     deleteDemandDrivenData(interpolatorABPtr_);
  73     deleteDemandDrivenData(interpolatorBAPtr_);
  74     deleteDemandDrivenData(controlPointsPtr_);
  75     deleteDemandDrivenData(motionPointsMaskPtr_);
  76     deleteDemandDrivenData(pointsDisplacementDirPtr_);
  77     deleteDemandDrivenData(facesDisplacementDirPtr_);
  78     deleteDemandDrivenData(totalDisplacementPtr_);
  79     deleteDemandDrivenData(aMeshPtr_);
  80     deleteDemandDrivenData(UsPtr_);
  81     deleteDemandDrivenData(phisPtr_);
  82     deleteDemandDrivenData(surfactConcPtr_);
  83     deleteDemandDrivenData(surfaceTensionPtr_);
  84     deleteDemandDrivenData(surfactantPtr_);
  85     deleteDemandDrivenData(fluidIndicatorPtr_);
  86 }
  87
  88 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  89
  90 freeSurface::freeSurface
  91 (
  92     dynamicFvMesh& m,
  93     const volScalarField& rho,
  94     volVectorField& Ub,
  95     volScalarField& Pb,
  96     const surfaceScalarField& sfPhi
  97 )
  98 :
  99     IOdictionary
 100     (
 101         IOobject
 102         (
 103             "freeSurfaceProperties",
 104             Ub.mesh().time().constant(),
 105             Ub.mesh(),
 106             IOobject::MUST_READ,
 107             IOobject::NO_WRITE
 108         )
 109     ),
 110     mesh_(m),
 111     rho_(rho),
 112     U_(Ub),
 113     p_(Pb),
 114     phi_(sfPhi),
 115     curTimeIndex_(Ub.mesh().time().timeIndex()),
 116     twoFluids_
 117     (
 118         this->lookup("twoFluids")
 119     ),
 120     normalMotionDir_
 121     (
 122         this->lookup("normalMotionDir")
 123     ),
 124     motionDir_(0, 0, 0),
 125     cleanInterface_
 126     (
 127         this->lookup("cleanInterface")
 128     ),
 129     aPatchID_(-1),
 130     bPatchID_(-1),
 131     muFluidA_
 132     (
 133         this->lookup("muFluidA")
 134     ),
 135     muFluidB_
 136     (
 137         this->lookup("muFluidB")
 138     ),
 139     rhoFluidA_
 140     (
 141         this->lookup("rhoFluidA")
 142     ),
 143     rhoFluidB_
 144     (
 145         this->lookup("rhoFluidB")
 146     ),
 147     g_(this->lookup("g")),
 148     cleanInterfaceSurfTension_
 149     (
 150         this->lookup("surfaceTension")
 151     ),
 152     fixedFreeSurfacePatches_
 153     (
 154         this->lookup("fixedFreeSurfacePatches")
 155     ),
 156     pointNormalsCorrectionPatches_
 157     (
 158         this->lookup("pointNormalsCorrectionPatches")
 159     ),
 160     nFreeSurfCorr_
 161     (
 162         readInt(this->lookup("nFreeSurfaceCorrectors"))
 163     ),
 164     smoothing_(false),
 165     interpolatorABPtr_(NULL),
 166     interpolatorBAPtr_(NULL),
 167     controlPointsPtr_(NULL),
 168     motionPointsMaskPtr_(NULL),
 169     pointsDisplacementDirPtr_(NULL),
 170     facesDisplacementDirPtr_(NULL),
 171     totalDisplacementPtr_(NULL),
 172     aMeshPtr_(NULL),
 173     UsPtr_(NULL),
 174     phisPtr_(NULL),
 175     surfactConcPtr_(NULL),
 176     surfaceTensionPtr_(NULL),
 177     surfactantPtr_(NULL),
 178     fluidIndicatorPtr_(NULL)
 179 {
 180     //Read motion direction
 181     if (!normalMotionDir_)
 182     {
 183         motionDir_ = vector(this->lookup("motionDir"));
 184         motionDir_ /= mag(motionDir_) + SMALL;
 185     }
 186
 187     // Set point normal correction patches
 188     boolList& correction = aMesh().correctPatchPointNormals();
 189
 190     forAll(pointNormalsCorrectionPatches_, patchI)
 191     {
 192         word patchName = pointNormalsCorrectionPatches_[patchI];
 193
 194         label patchID = aMesh().boundary().findPatchID(patchName);
 195
 196         if(patchID == -1)
 197         {
 198             FatalErrorIn
 199             (
 200                 "freeSurface::freeSurface(...)"
 201             )   << "Patch name for point normals correction does not exist"
 202                 << abort(FatalError);
 203         }
 204
 205         correction[patchID] = true;
 206     }
 207
 208     // Clear geometry
 209     aMesh().movePoints();
 210
 211
 212     // Detect the free surface patch
 213     forAll (mesh().boundary(), patchI)
 214     {
 215         if(mesh().boundary()[patchI].name() == "freeSurface")
 216         {
 217             aPatchID_ = patchI;
 218
 219             Info<< "Found free surface patch. ID: " << aPatchID_
 220                 << endl;
 221         }
 222     }
 223
 224     if(aPatchID() == -1)
 225     {
 226         FatalErrorIn("freeSurface::freeSurface(...)")
 227             << "Free surface patch not defined.  Please make sure that "
 228                 << " the free surface patches is named as freeSurface"
 229                 << abort(FatalError);
 230     }
 231
 232
 233     // Detect the free surface shadow patch
 234     if (twoFluids())
 235     {
 236         forAll (mesh().boundary(), patchI)
 237         {
 238             if(mesh().boundary()[patchI].name() == "freeSurfaceShadow")
 239             {
 240                 bPatchID_ = patchI;
 241
 242                 Info<< "Found free surface shadow patch. ID: "
 243                     << bPatchID_ << endl;
 244             }
 245         }
 246
 247         if(bPatchID() == -1)
 248         {
 249             FatalErrorIn("freeSurface::freeSurface(...)")
 250                 << "Free surface shadow patch not defined. "
 251                     << "Please make sure that the free surface shadow patch "
 252                     << "is named as freeSurfaceShadow."
 253                     << abort(FatalError);
 254         }
 255     }
 256
 257
 258     // Mark free surface boundary points
 259     // which belonge to processor patches
 260     forAll(aMesh().boundary(), patchI)
 261     {
 262         if
 263         (
 264             aMesh().boundary()[patchI].type()
 265          == processorFaPatch::typeName
 266         )
 267         {
 268             const labelList& patchPoints =
 269                 aMesh().boundary()[patchI].pointLabels();
 270
 271             forAll(patchPoints, pointI)
 272             {
 273                 motionPointsMask()[patchPoints[pointI]] = -1;
 274             }
 275         }
 276     }
 277
 278
 279     // Mark fixed free surface boundary points
 280     forAll(fixedFreeSurfacePatches_, patchI)
 281     {
 282         label fixedPatchID =
 283             aMesh().boundary().findPatchID
 284             (
 285                 fixedFreeSurfacePatches_[patchI]
 286             );
 287
 288         if(fixedPatchID == -1)
 289         {
 290             FatalErrorIn("freeSurface::freeSurface(...)")
 291                 << "Wrong faPatch name in the fixedFreeSurfacePatches list"
 292                     << " defined in the freeSurfaceProperties dictionary"
 293                     << abort(FatalError);
 294         }
 295
 296         const labelList& patchPoints =
 297             aMesh().boundary()[fixedPatchID].pointLabels();
 298
 299         forAll(patchPoints, pointI)
 300         {
 301             motionPointsMask()[patchPoints[pointI]] = 0;
 302         }
 303     }
 304
 305
 306     // Mark free-surface boundary point
 307     // at the axis of 2-D axisymmetic cases
 308     forAll(aMesh().boundary(), patchI)
 309     {
 310         if
 311         (
 312             aMesh().boundary()[patchI].type()
 313          == wedgeFaPatch::typeName
 314         )
 315         {
 316             const wedgeFaPatch& wedgePatch =
 317                 refCast<const wedgeFaPatch>(aMesh().boundary()[patchI]);
 318
 319             if(wedgePatch.axisPoint() > -1)
 320             {
 321                 motionPointsMask()[wedgePatch.axisPoint()] = 0;
 322
 323                 Info << "Axis point: "
 324                     << wedgePatch.axisPoint()
 325                     << "vector: "
 326                     << aMesh().points()[wedgePatch.axisPoint()] << endl;
 327             }
 328         }
 329     }
 330
 331
 332     // Read free-surface points total displacement if present
 333     readTotalDisplacement();
 334
 335
 336     // Read control points positions if present
 337     controlPoints();
 338
 339
 340     // Check if smoothing switch is set
 341     if (this->found("smoothing"))
 342     {
 343         smoothing_ = Switch(this->lookup("smoothing"));
 344     }
 345 }
 346
 347
 348 // * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * * //
 349
 350 freeSurface::~freeSurface()
 351 {
 352     clearOut();
 353 }
 354
 355
 356 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 357
 358 void freeSurface::updateDisplacementDirections()
 359 {
 360     if(normalMotionDir())
 361     {
 362         // Update point displacement correction
 363         pointsDisplacementDir() = aMesh().pointAreaNormals();
 364
 365         // Correcte point displacement direction
 366         // at the "centerline" symmetryPlane which represents the axis
 367         // of an axisymmetric case
 368         forAll(aMesh().boundary(), patchI)
 369         {
 370             if(aMesh().boundary()[patchI].type() == wedgeFaPatch::typeName)
 371             {
 372                 const wedgeFaPatch& wedgePatch =
 373                     refCast<const wedgeFaPatch>(aMesh().boundary()[patchI]);
 374
 375                 vector axis = wedgePatch.axis();
 376
 377                 label centerLinePatchID =
 378                     aMesh().boundary().findPatchID("centerline");
 379
 380                 if(centerLinePatchID != -1)
 381                 {
 382                     const labelList& pointLabels =
 383                         aMesh().boundary()[centerLinePatchID].pointLabels();
 384
 385                     forAll(pointLabels, pointI)
 386                     {
 387                         vector dir =
 388                             pointsDisplacementDir()[pointLabels[pointI]];
 389
 390                         dir = (dir&axis)*axis;
 391                         dir /= mag(dir);
 392
 393                         pointsDisplacementDir()[pointLabels[pointI]] = dir;
 394                     }
 395                 }
 396                 else
 397                 {
 398                     Info << "Warning: centerline polyPatch does not exist. "
 399                         << "Free surface points displacement directions "
 400                         << "will not be corrected at the axis (centerline)"
 401                         << endl;
 402                 }
 403
 404                 break;
 405             }
 406         }
 407
 408         // Update face displacement direction
 409         facesDisplacementDir() =
 410             aMesh().faceAreaNormals().internalField();
 411
 412         // Correction of control points postion
 413         const vectorField& Cf = aMesh().areaCentres().internalField();
 414
 415         controlPoints() =
 416             facesDisplacementDir()
 417            *(facesDisplacementDir()&(controlPoints() - Cf))
 418           + Cf;
 419     }
 420 }
 421
 422
 423 bool freeSurface::predictPoints()
 424 {
 425     // Smooth interface
 426
 427     if (smoothing_)
 428     {
 429         controlPoints() = aMesh().areaCentres().internalField();
 430         movePoints(scalarField(controlPoints().size(), 0));
 431         movePoints(-fvc::meshPhi(rho(),U())().boundaryField()[aPatchID()]);
 432     }
 433
 434     for
 435     (
 436         int freeSurfCorr=0;
 437         freeSurfCorr<nFreeSurfCorr_;
 438         freeSurfCorr++
 439     )
 440     {
 441         movePoints(phi_.boundaryField()[aPatchID()]);
 442     }
 443
 444     return true;
 445 }
 446
 447
 448 bool freeSurface::correctPoints()
 449 {
 450     for
 451     (
 452         int freeSurfCorr=0;
 453         freeSurfCorr<nFreeSurfCorr_;
 454         freeSurfCorr++
 455     )
 456     {
 457         movePoints(phi_.boundaryField()[aPatchID()]);
 458     }
 459
 460     return true;
 461 }
 462
 463
 464 bool freeSurface::movePoints(const scalarField& interfacePhi)
 465 {
 466     pointField newMeshPoints = mesh().points();
 467
 468     scalarField sweptVolCorr =
 469         interfacePhi
 470       - fvc::meshPhi(rho(),U())().boundaryField()[aPatchID()];
 471
 472     word ddtScheme
 473     (
 474         mesh().schemesDict().ddtScheme
 475         (
 476             "ddt(" + rho().name() + ',' + U().name()+')'
 477         )
 478     );
 479
 480     if
 481     (
 482         ddtScheme
 483      == fv::CrankNicolsonDdtScheme<vector>::typeName
 484     )
 485     {
 486         sweptVolCorr *= (1.0/2.0)*DB().deltaT().value();
 487     }
 488     else if
 489     (
 490         ddtScheme
 491      == fv::EulerDdtScheme<vector>::typeName
 492     )
 493     {
 494         sweptVolCorr *= DB().deltaT().value();
 495     }
 496     else if
 497     (
 498         ddtScheme
 499      == fv::backwardDdtScheme<vector>::typeName
 500     )
 501     {
 502         if (DB().timeIndex() == 1)
 503         {
 504             sweptVolCorr *= DB().deltaT().value();
 505         }
 506         else
 507         {
 508             sweptVolCorr *= (2.0/3.0)*DB().deltaT().value();
 509         }
 510     }
 511     else
 512     {
 513         FatalErrorIn("freeSurface::movePoints()")
 514             << "Unsupported temporal differencing scheme : "
 515                 << ddtScheme
 516                 << abort(FatalError);
 517     }
 518
 519     const scalarField& Sf = aMesh().S();
 520     const vectorField& Nf = aMesh().faceAreaNormals().internalField();
 521
 522     scalarField deltaH =
 523         sweptVolCorr/(Sf*(Nf & facesDisplacementDir()));
 524
 525     pointField displacement = pointDisplacement(deltaH);
 526
 527
 528     // Move only free-surface points
 529
 530     const labelList& meshPointsA =
 531         mesh().boundaryMesh()[aPatchID()].meshPoints();
 532
 533     forAll (displacement, pointI)
 534     {
 535         newMeshPoints[meshPointsA[pointI]] += displacement[pointI];
 536     }
 537
 538     if(twoFluids_)
 539     {
 540         const labelList& meshPointsB =
 541             mesh().boundaryMesh()[bPatchID_].meshPoints();
 542
 543         pointField displacementB =
 544             interpolatorAB().pointInterpolate
 545             (
 546                 displacement
 547             );
 548
 549         forAll (displacementB, pointI)
 550         {
 551             newMeshPoints[meshPointsB[pointI]] += displacementB[pointI];
 552         }
 553     }
 554
 555     // Update total displacement field
 556
 557     if(totalDisplacementPtr_ && (curTimeIndex_ < DB().timeIndex()))
 558     {
 559         FatalErrorIn("freeSurface::movePoints()")
 560             << "Total displacement of free surface points "
 561                 << "from previous time step is not absorbed by the mesh."
 562                 << abort(FatalError);
 563     }
 564     else if (curTimeIndex_ < DB().timeIndex())
 565     {
 566         totalDisplacement() = displacement;
 567
 568         curTimeIndex_ = DB().timeIndex();
 569     }
 570     else
 571     {
 572         totalDisplacement() += displacement;
 573     }
 574
 575     twoDPointCorrector twoDPointCorr(mesh());
 576
 577     twoDPointCorr.correctPoints(newMeshPoints);
 578
 579     mesh().movePoints(newMeshPoints);
 580
 581     // faMesh motion is done automatically, using meshObject
 582     // HJ, 8/Aug/2011
 583 //     aMesh().movePoints(mesh().points());
 584
 585
 586     // Move correctedFvPatchField fvSubMeshes
 587
 588     forAll(U().boundaryField(), patchI)
 589     {
 590         if
 591         (
 592             (
 593                 U().boundaryField()[patchI].type()
 594              == fixedGradientCorrectedFvPatchField<vector>::typeName
 595             )
 596             ||
 597             (
 598                 U().boundaryField()[patchI].type()
 599              == fixedValueCorrectedFvPatchField<vector>::typeName
 600             )
 601             ||
 602             (
 603                 U().boundaryField()[patchI].type()
 604              == zeroGradientCorrectedFvPatchField<vector>::typeName
 605             )
 606         )
 607         {
 608             correctedFvPatchField<vector>& pU =
 609                 refCast<correctedFvPatchField<vector> >
 610                 (
 611                     U().boundaryField()[patchI]
 612                 );
 613
 614             pU.movePatchSubMesh();
 615         }
 616     }
 617
 618     forAll(p().boundaryField(), patchI)
 619     {
 620         if
 621         (
 622             (
 623                 p().boundaryField()[patchI].type()
 624              == fixedGradientCorrectedFvPatchField<scalar>::typeName
 625             )
 626             ||
 627             (
 628                 p().boundaryField()[patchI].type()
 629              == fixedValueCorrectedFvPatchField<scalar>::typeName
 630             )
 631             ||
 632             (
 633                 p().boundaryField()[patchI].type()
 634              == zeroGradientCorrectedFvPatchField<scalar>::typeName
 635             )
 636         )
 637         {
 638             correctedFvPatchField<scalar>& pP =
 639                 refCast<correctedFvPatchField<scalar> >
 640                 (
 641                     p().boundaryField()[patchI]
 642                 );
 643
 644             pP.movePatchSubMesh();
 645         }
 646     }
 647
 648     return true;
 649 }
 650
 651
 652 bool freeSurface::moveMeshPointsForOldFreeSurfDisplacement()
 653 {
 654     if(totalDisplacementPtr_)
 655     {
 656         pointField newPoints = mesh().points();
 657
 658         const labelList& meshPointsA =
 659             mesh().boundaryMesh()[aPatchID()].meshPoints();
 660
 661         forAll (totalDisplacement(), pointI)
 662         {
 663             newPoints[meshPointsA[pointI]] -= totalDisplacement()[pointI];
 664         }
 665
 666
 667         // Check mesh motion solver type
 668         bool feMotionSolver =
 669             mesh().objectRegistry::foundObject<tetPointVectorField>
 670             (
 671                 "motionU"
 672             );
 673         bool fvMotionSolver =
 674             mesh().objectRegistry::foundObject<pointVectorField>
 675             (
 676                 "pointMotionU"
 677             );
 678
 679         if (feMotionSolver)
 680         {
 681             tetPointVectorField& motionU =
 682                 const_cast<tetPointVectorField&>
 683                 (
 684                     mesh().objectRegistry::
 685                     lookupObject<tetPointVectorField>
 686                     (
 687                         "motionU"
 688                     )
 689                 );
 690
 691             fixedValueTetPolyPatchVectorField& motionUaPatch =
 692                 refCast<fixedValueTetPolyPatchVectorField>
 693                 (
 694                     motionU.boundaryField()[aPatchID()]
 695                 );
 696
 697             tetPolyPatchInterpolation tppiAPatch
 698             (
 699                 refCast<const faceTetPolyPatch>
 700                 (
 701                     motionUaPatch.patch()
 702                 )
 703             );
 704
 705             motionUaPatch ==
 706                 tppiAPatch.pointToPointInterpolate
 707                 (
 708                     totalDisplacement()/DB().deltaT().value()
 709                 );
 710
 711             if(twoFluids_)
 712             {
 713                 const labelList& meshPointsB =
 714                     mesh().boundaryMesh()[bPatchID()].meshPoints();
 715
 716                 pointField totDisplacementB =
 717                     interpolatorAB().pointInterpolate
 718                     (
 719                         totalDisplacement()
 720                     );
 721
 722                 forAll (totDisplacementB, pointI)
 723                 {
 724                     newPoints[meshPointsB[pointI]] -=
 725                         totDisplacementB[pointI];
 726                 }
 727
 728                 fixedValueTetPolyPatchVectorField& motionUbPatch =
 729                     refCast<fixedValueTetPolyPatchVectorField>
 730                     (
 731                         motionU.boundaryField()[bPatchID()]
 732                     );
 733
 734                 tetPolyPatchInterpolation tppiBPatch
 735                 (
 736                     refCast<const faceTetPolyPatch>(motionUbPatch.patch())
 737                 );
 738
 739                 motionUbPatch ==
 740                     tppiBPatch.pointToPointInterpolate
 741                     (
 742                         totDisplacementB/DB().deltaT().value()
 743                     );
 744             }
 745         }
 746         else if (fvMotionSolver)
 747         {
 748             pointVectorField& motionU =
 749                 const_cast<pointVectorField&>
 750                 (
 751                     mesh().objectRegistry::
 752                     lookupObject<pointVectorField>
 753                     (
 754                         "pointMotionU"
 755                     )
 756                 );
 757
 758             fixedValuePointPatchVectorField& motionUaPatch =
 759                 refCast<fixedValuePointPatchVectorField>
 760                 (
 761                     motionU.boundaryField()[aPatchID()]
 762                 );
 763
 764             motionUaPatch ==
 765                 totalDisplacement()/DB().deltaT().value();
 766
 767             if(twoFluids_)
 768             {
 769                 const labelList& meshPointsB =
 770                     mesh().boundaryMesh()[bPatchID()].meshPoints();
 771
 772                 pointField totDisplacementB =
 773                     interpolatorAB().pointInterpolate
 774                     (
 775                         totalDisplacement()
 776                     );
 777
 778                 forAll (totDisplacementB, pointI)
 779                 {
 780                     newPoints[meshPointsB[pointI]] -=
 781                         totDisplacementB[pointI];
 782                 }
 783
 784                 fixedValuePointPatchVectorField& motionUbPatch =
 785                     refCast<fixedValuePointPatchVectorField>
 786                     (
 787                         motionU.boundaryField()[bPatchID()]
 788                     );
 789
 790                 motionUbPatch ==
 791                     totDisplacementB/DB().deltaT().value();
 792             }
 793         }
 794
 795         twoDPointCorrector twoDPointCorr(mesh());
 796
 797         twoDPointCorr.correctPoints(newPoints);
 798
 799         mesh().movePoints(newPoints);
 800
 801         deleteDemandDrivenData(totalDisplacementPtr_);
 802
 803         mesh().update();
 804
 805         // faMesh motion is done automatically, using meshObject
 806         // HJ, 8/Aug/2011
 807 //         aMesh().movePoints(mesh().points());
 808
 809         // Move correctedFvPatchField fvSubMeshes
 810
 811         forAll(U().boundaryField(), patchI)
 812         {
 813             if
 814             (
 815                 (
 816                     U().boundaryField()[patchI].type()
 817                  == fixedGradientCorrectedFvPatchField<vector>::typeName
 818                 )
 819                 ||
 820                 (
 821                     U().boundaryField()[patchI].type()
 822                  == fixedValueCorrectedFvPatchField<vector>::typeName
 823                 )
 824                 ||
 825                 (
 826                     U().boundaryField()[patchI].type()
 827                  == zeroGradientCorrectedFvPatchField<vector>::typeName
 828                 )
 829             )
 830             {
 831                 correctedFvPatchField<vector>& aU =
 832                     refCast<correctedFvPatchField<vector> >
 833                     (
 834                         U().boundaryField()[patchI]
 835                     );
 836
 837                 aU.movePatchSubMesh();
 838             }
 839         }
 840
 841         forAll(p().boundaryField(), patchI)
 842         {
 843             if
 844             (
 845                 (
 846                     p().boundaryField()[patchI].type()
 847                  == fixedGradientCorrectedFvPatchField<scalar>::typeName
 848                 )
 849                 ||
 850                 (
 851                     p().boundaryField()[patchI].type()
 852                  == fixedValueCorrectedFvPatchField<scalar>::typeName
 853                 )
 854                 ||
 855                 (
 856                     p().boundaryField()[patchI].type()
 857                  == zeroGradientCorrectedFvPatchField<scalar>::typeName
 858                 )
 859             )
 860             {
 861                 correctedFvPatchField<scalar>& aP =
 862                     refCast<correctedFvPatchField<scalar> >
 863                     (
 864                         p().boundaryField()[patchI]
 865                     );
 866
 867                 aP.movePatchSubMesh();
 868             }
 869         }
 870     }
 871
 872     return true;
 873 }
 874
 875
 876 bool freeSurface::moveMeshPoints()
 877 {
 878         scalarField sweptVolCorr =
 879             phi_.boundaryField()[aPatchID()]
 880           - fvc::meshPhi(rho(),U())().boundaryField()[aPatchID()];
 881
 882         word ddtScheme
 883         (
 884             mesh().schemesDict().ddtScheme
 885             (
 886                 "ddt(" + rho().name() + ',' + U().name()+')'
 887             )
 888         );
 889
 890         if
 891         (
 892             ddtScheme
 893          == fv::CrankNicolsonDdtScheme<vector>::typeName
 894         )
 895         {
 896             sweptVolCorr *= (1.0/2.0)*DB().deltaT().value();
 897         }
 898         else if
 899         (
 900             ddtScheme
 901          == fv::EulerDdtScheme<vector>::typeName
 902         )
 903         {
 904             sweptVolCorr *= DB().deltaT().value();
 905         }
 906         else if
 907         (
 908             ddtScheme
 909          == fv::backwardDdtScheme<vector>::typeName
 910         )
 911         {
 912             sweptVolCorr *= (2.0/3.0)*DB().deltaT().value();
 913         }
 914         else
 915         {
 916             FatalErrorIn("freeSurface::movePoints()")
 917                 << "Unsupported temporal differencing scheme : "
 918                 << ddtScheme
 919                 << abort(FatalError);
 920         }
 921
 922
 923         const scalarField& Sf = aMesh().S();
 924         const vectorField& Nf = aMesh().faceAreaNormals().internalField();
 925
 926         scalarField deltaH =
 927             sweptVolCorr/(Sf*(Nf & facesDisplacementDir()));
 928
 929
 930         pointField displacement = pointDisplacement(deltaH);
 931
 932
 933         //-- Set mesh motion boundary conditions
 934
 935         tetPointVectorField& motionU =
 936             const_cast<tetPointVectorField&>
 937             (
 938                 mesh().objectRegistry::
 939                 lookupObject<tetPointVectorField>
 940                 (
 941                     "motionU"
 942                 )
 943             );
 944
 945         fixedValueTetPolyPatchVectorField& motionUaPatch =
 946             refCast<fixedValueTetPolyPatchVectorField>
 947             (
 948                 motionU.boundaryField()[aPatchID()]
 949             );
 950
 951         tetPolyPatchInterpolation tppiAPatch
 952         (
 953             refCast<const faceTetPolyPatch>
 954             (
 955                 motionUaPatch.patch()
 956             )
 957         );
 958
 959         motionUaPatch ==
 960             tppiAPatch.pointToPointInterpolate
 961             (
 962                 displacement/DB().deltaT().value()
 963             );
 964
 965         if (twoFluids())
 966         {
 967             fixedValueTetPolyPatchVectorField& motionUbPatch =
 968                 refCast<fixedValueTetPolyPatchVectorField>
 969                 (
 970                     motionU.boundaryField()[bPatchID()]
 971                 );
 972
 973             tetPolyPatchInterpolation tppiBPatch
 974             (
 975                 refCast<const faceTetPolyPatch>(motionUbPatch.patch())
 976             );
 977
 978             motionUbPatch ==
 979                 tppiBPatch.pointToPointInterpolate
 980                 (
 981                     interpolatorAB().pointInterpolate
 982                     (
 983                         displacement/DB().deltaT().value()
 984                     )
 985                 );
 986         }
 987
 988         mesh().update();
 989
 990         // faMesh motion is done automatically, using meshObject
 991         // HJ, 8/Aug/2011
 992 //         aMesh().movePoints(mesh().points());
 993
 994
 995         // Move correctedFvPatchField fvSubMeshes
 996
 997         forAll(U().boundaryField(), patchI)
 998         {
 999             if
1000             (
1001                 (
1002                     U().boundaryField()[patchI].type()
1003                  == fixedGradientCorrectedFvPatchField<vector>::typeName
1004                 )
1005                 ||
1006                 (
1007                     U().boundaryField()[patchI].type()
1008                  == fixedValueCorrectedFvPatchField<vector>::typeName
1009                 )
1010                 ||
1011                 (
1012                     U().boundaryField()[patchI].type()
1013                  == zeroGradientCorrectedFvPatchField<vector>::typeName
1014                 )
1015             )
1016             {
1017                 correctedFvPatchField<vector>& aU =
1018                     refCast<correctedFvPatchField<vector> >
1019                     (
1020                         U().boundaryField()[patchI]
1021                     );
1022
1023                 aU.movePatchSubMesh();
1024             }
1025         }
1026
1027         forAll(p().boundaryField(), patchI)
1028         {
1029             if
1030             (
1031                 (
1032                     p().boundaryField()[patchI].type()
1033                  == fixedGradientCorrectedFvPatchField<scalar>::typeName
1034                 )
1035                 ||
1036                 (
1037                     p().boundaryField()[patchI].type()
1038                  == fixedValueCorrectedFvPatchField<scalar>::typeName
1039                 )
1040                 ||
1041                 (
1042                     p().boundaryField()[patchI].type()
1043                  == zeroGradientCorrectedFvPatchField<scalar>::typeName
1044                 )
1045             )
1046             {
1047                 correctedFvPatchField<scalar>& aP =
1048                     refCast<correctedFvPatchField<scalar> >
1049                     (
1050                         p().boundaryField()[patchI]
1051                     );
1052
1053                 aP.movePatchSubMesh();
1054             }
1055         }
1056
1057     return true;
1058 }
1059
1060
1061 void freeSurface::updateBoundaryConditions()
1062 {
1063     updateVelocity();
1064     updateSurfactantConcentration();
1065     updatePressure();
1066 }
1067
1068
1069 void freeSurface::updateVelocity()
1070 {
1071     if(twoFluids())
1072     {
1073         vectorField nA = mesh().boundary()[aPatchID()].nf();
1074
1075         vectorField nB = mesh().boundary()[bPatchID()].nf();
1076
1077         scalarField DnB = interpolatorBA().faceInterpolate
1078         (
1079             mesh().boundary()[bPatchID()].deltaCoeffs()
1080         );
1081
1082         scalarField DnA = mesh().boundary()[aPatchID()].deltaCoeffs();
1083
1084
1085         vectorField UtPA =
1086             U().boundaryField()[aPatchID()].patchInternalField();
1087
1088         if
1089         (
1090             U().boundaryField()[aPatchID()].type()
1091          == fixedGradientCorrectedFvPatchField<vector>::typeName
1092         )
1093         {
1094             fixedGradientCorrectedFvPatchField<vector>& aU =
1095                 refCast<fixedGradientCorrectedFvPatchField<vector> >
1096                 (
1097                     U().boundaryField()[aPatchID()]
1098                 );
1099
1100             UtPA += aU.corrVecGrad();
1101         }
1102
1103         UtPA -= nA*(nA & UtPA);
1104
1105
1106         vectorField UtPB = interpolatorBA().faceInterpolate
1107         (
1108             U().boundaryField()[bPatchID()].patchInternalField()
1109         );
1110
1111         if
1112         (
1113             U().boundaryField()[bPatchID()].type()
1114          == fixedValueCorrectedFvPatchField<vector>::typeName
1115         )
1116         {
1117             fixedValueCorrectedFvPatchField<vector>& bU =
1118                 refCast<fixedValueCorrectedFvPatchField<vector> >
1119                 (
1120                     U().boundaryField()[bPatchID()]
1121                 );
1122
1123             UtPB += interpolatorBA().faceInterpolate(bU.corrVecGrad());
1124         }
1125
1126         UtPB -= nA*(nA & UtPB);
1127
1128         vectorField UtFs = muFluidA().value()*DnA*UtPA
1129           + muFluidB().value()*DnB*UtPB;
1130
1131         vectorField UnFs =
1132             nA*phi_.boundaryField()[aPatchID()]
1133            /mesh().boundary()[aPatchID()].magSf();
1134
1135         Us().internalField() += UnFs - nA*(nA&Us().internalField());
1136         correctUsBoundaryConditions();
1137
1138         UtFs -= (muFluidA().value() - muFluidB().value())*
1139             (fac::grad(Us())&aMesh().faceAreaNormals())().internalField();
1140
1141
1142         vectorField tangentialSurfaceTensionForce(nA.size(), vector::zero);
1143
1144         if(!cleanInterface())
1145         {
1146             tangentialSurfaceTensionForce =
1147                 surfaceTensionGrad()().internalField();
1148         }
1149         else
1150         {
1151             vectorField surfaceTensionForce =
1152                 cleanInterfaceSurfTension().value()
1153                *fac::edgeIntegrate
1154                 (
1155                     aMesh().Le()*aMesh().edgeLengthCorrection()
1156                 )().internalField();
1157
1158             tangentialSurfaceTensionForce =
1159                 surfaceTensionForce
1160               - cleanInterfaceSurfTension().value()
1161                *aMesh().faceCurvatures().internalField()*nA;
1162         }
1163
1164         UtFs += tangentialSurfaceTensionForce;
1165
1166         UtFs /= muFluidA().value()*DnA + muFluidB().value()*DnB + VSMALL;
1167
1168         Us().internalField() = UnFs + UtFs;
1169         correctUsBoundaryConditions();
1170
1171         // Store old-time velocity field U()
1172         U().oldTime();
1173
1174         U().boundaryField()[bPatchID()] ==
1175             interpolatorAB().faceInterpolate(UtFs)
1176           + nB*fvc::meshPhi(rho(),U())().boundaryField()[bPatchID()]/
1177             mesh().boundary()[bPatchID()].magSf();
1178
1179         if
1180         (
1181             p().boundaryField()[bPatchID()].type()
1182          == fixedGradientFvPatchField<scalar>::typeName
1183         )
1184         {
1185             fixedGradientFvPatchField<scalar>& pB =
1186                 refCast<fixedGradientFvPatchField<scalar> >
1187                 (
1188                     p().boundaryField()[bPatchID()]
1189                 );
1190
1191             pB.gradient() =
1192                - rhoFluidB().value()
1193                 *(
1194                      nB&fvc::ddt(U())().boundaryField()[bPatchID()]
1195                  );
1196         }
1197
1198
1199         // Update fixedGradient boundary condition on patch A
1200
1201         vectorField nGradU =
1202             muFluidB().value()*(UtPB - UtFs)*DnA
1203           + tangentialSurfaceTensionForce
1204           - muFluidA().value()*nA*fac::div(Us())().internalField()
1205           + (muFluidB().value() - muFluidA().value())
1206            *(fac::grad(Us())().internalField()&nA);
1207
1208         nGradU /= muFluidA().value() + VSMALL;
1209
1210         if
1211         (
1212             U().boundaryField()[aPatchID()].type()
1213          == fixedGradientCorrectedFvPatchField<vector>::typeName
1214         )
1215         {
1216             fixedGradientCorrectedFvPatchField<vector>& aU =
1217                 refCast<fixedGradientCorrectedFvPatchField<vector> >
1218                 (
1219                     U().boundaryField()[aPatchID()]
1220                 );
1221
1222             aU.gradient() = nGradU;
1223         }
1224         else if
1225         (
1226             U().boundaryField()[aPatchID()].type()
1227          == fixedGradientFvPatchField<vector>::typeName
1228         )
1229         {
1230             fixedGradientFvPatchField<vector>& aU =
1231                 refCast<fixedGradientFvPatchField<vector> >
1232                 (
1233                     U().boundaryField()[aPatchID()]
1234                 );
1235
1236             aU.gradient() = nGradU;
1237         }
1238         else
1239         {
1240             FatalErrorIn("freeSurface::updateVelocity()")
1241                 << "Bounary condition on " << U().name()
1242                     <<  " for freeSurface patch is "
1243                     << U().boundaryField()[aPatchID()].type()
1244                     << ", instead "
1245                     << fixedGradientCorrectedFvPatchField<vector>::typeName
1246                     << " or "
1247                     << fixedGradientFvPatchField<vector>::typeName
1248                     << abort(FatalError);
1249         }
1250     }
1251     else
1252     {
1253         vectorField nA = aMesh().faceAreaNormals().internalField();
1254
1255         vectorField UnFs =
1256             nA*phi_.boundaryField()[aPatchID()]
1257            /mesh().boundary()[aPatchID()].magSf();
1258
1259         // Correct normal component of free-surface velocity
1260         Us().internalField() += UnFs - nA*(nA&Us().internalField());
1261         correctUsBoundaryConditions();
1262
1263         vectorField tangentialSurfaceTensionForce(nA.size(), vector::zero);
1264
1265         if(!cleanInterface())
1266         {
1267             tangentialSurfaceTensionForce =
1268                 surfaceTensionGrad()().internalField();
1269         }
1270         else
1271         {
1272             vectorField surfaceTensionForce =
1273                 cleanInterfaceSurfTension().value()
1274                *fac::edgeIntegrate
1275                 (
1276                     aMesh().Le()*aMesh().edgeLengthCorrection()
1277                 )().internalField();
1278
1279             tangentialSurfaceTensionForce =
1280                 surfaceTensionForce
1281               - cleanInterfaceSurfTension().value()
1282                *aMesh().faceCurvatures().internalField()*nA;
1283
1284             if (muFluidA().value() < SMALL)
1285             {
1286                 tangentialSurfaceTensionForce = vector::zero;
1287             }
1288         }
1289
1290         vectorField tnGradU =
1291             tangentialSurfaceTensionForce/(muFluidA().value() + VSMALL)
1292           - (fac::grad(Us())&aMesh().faceAreaNormals())().internalField();
1293
1294         vectorField UtPA =
1295             U().boundaryField()[aPatchID()].patchInternalField();
1296         UtPA -= nA*(nA & UtPA);
1297
1298         scalarField DnA = mesh().boundary()[aPatchID()].deltaCoeffs();
1299
1300         vectorField UtFs = UtPA + tnGradU/DnA;
1301
1302         Us().internalField() = UtFs + UnFs;
1303         correctUsBoundaryConditions();
1304
1305         vectorField nGradU =
1306             tangentialSurfaceTensionForce/(muFluidA().value() + VSMALL)
1307           - nA*fac::div(Us())().internalField()
1308           - (fac::grad(Us())().internalField()&nA);
1309
1310         if
1311         (
1312             U().boundaryField()[aPatchID()].type()
1313          == fixedGradientCorrectedFvPatchField<vector>::typeName
1314         )
1315         {
1316             fixedGradientCorrectedFvPatchField<vector>& aU =
1317                 refCast<fixedGradientCorrectedFvPatchField<vector> >
1318                 (
1319                     U().boundaryField()[aPatchID()]
1320                 );
1321
1322             aU.gradient() = nGradU;
1323         }
1324         else if
1325         (
1326             U().boundaryField()[aPatchID()].type()
1327          == fixedGradientFvPatchField<vector>::typeName
1328         )
1329         {
1330             fixedGradientFvPatchField<vector>& aU =
1331                 refCast<fixedGradientFvPatchField<vector> >
1332                 (
1333                     U().boundaryField()[aPatchID()]
1334                 );
1335
1336             aU.gradient() = nGradU;
1337         }
1338         else
1339         {
1340             FatalErrorIn("freeSurface::updateVelocity()")
1341                 << "Bounary condition on " << U().name()
1342                     <<  " for freeSurface patch is "
1343                     << U().boundaryField()[aPatchID()].type()
1344                     << ", instead "
1345                     << fixedGradientCorrectedFvPatchField<vector>::typeName
1346                     << " or "
1347                     << fixedGradientFvPatchField<vector>::typeName
1348                     << abort(FatalError);
1349         }
1350     }
1351 }
1352
1353
1354 void freeSurface::updatePressure()
1355 {
1356     // Correct pressure boundary condition at the free-surface
1357
1358     vectorField nA = mesh().boundary()[aPatchID()].nf();
1359
1360     if(twoFluids())
1361     {
1362         scalarField pA =
1363             interpolatorBA().faceInterpolate
1364             (
1365                 p().boundaryField()[bPatchID()]
1366             );
1367
1368         const scalarField& K = aMesh().faceCurvatures().internalField();
1369
1370         Info << "Free surface curvature: min = " << gMin(K)
1371             << ", max = " << gMax(K)
1372             << ", average = " << gAverage(K) << endl << flush;
1373
1374         if(cleanInterface())
1375         {
1376 //             pA -= cleanInterfaceSurfTension().value()*(K - gAverage(K));
1377             pA -= cleanInterfaceSurfTension().value()*K;
1378         }
1379         else
1380         {
1381             scalarField surfTensionK =
1382                 surfaceTension().internalField()*K;
1383
1384             pA -= surfTensionK - gAverage(surfTensionK);
1385         }
1386
1387         pA -= 2.0*(muFluidA().value() - muFluidB().value())
1388             *fac::div(Us())().internalField();
1389
1390 //         vector R0 = gAverage(mesh().C().boundaryField()[aPatchID()]);
1391         vector R0 = vector::zero;
1392
1393         pA -= (rhoFluidA().value() - rhoFluidB().value())*
1394             (
1395                 g_.value()
1396               & (
1397                     mesh().C().boundaryField()[aPatchID()]
1398                   - R0
1399                 )
1400             );
1401
1402         p().boundaryField()[aPatchID()] == pA;
1403     }
1404     else
1405     {
1406 //         vector R0 = gAverage(mesh().C().boundaryField()[aPatchID()]);
1407         vector R0 = vector::zero;
1408
1409         scalarField pA =
1410           - rhoFluidA().value()*
1411             (
1412                 g_.value()
1413               & (
1414                     mesh().C().boundaryField()[aPatchID()]
1415                   - R0
1416                 )
1417             );
1418
1419         const scalarField& K = aMesh().faceCurvatures().internalField();
1420
1421         Info << "Free surface curvature: min = " << gMin(K)
1422             << ", max = " << gMax(K) << ", average = " << gAverage(K)
1423             << endl;
1424
1425         if(cleanInterface())
1426         {
1427 //             pA -= cleanInterfaceSurfTension().value()*(K - gAverage(K));
1428             pA -= cleanInterfaceSurfTension().value()*K;
1429         }
1430         else
1431         {
1432             scalarField surfTensionK =
1433                 surfaceTension().internalField()*K;
1434
1435             pA -= surfTensionK - gAverage(surfTensionK);
1436         }
1437
1438         pA -= 2.0*muFluidA().value()*fac::div(Us())().internalField();
1439
1440         p().boundaryField()[aPatchID()] == pA;
1441     }
1442
1443
1444     // Set modified pressure at patches with fixed apsolute
1445     // pressure
1446
1447 //     vector R0 = gAverage(mesh().C().boundaryField()[aPatchID()]);
1448     vector R0 = vector::zero;
1449
1450     for (int patchI=0; patchI < p().boundaryField().size(); patchI++)
1451     {
1452         if
1453         (
1454             p().boundaryField()[patchI].type()
1455          == fixedValueFvPatchScalarField::typeName
1456         )
1457         {
1458             if (patchI != aPatchID())
1459             {
1460                 p().boundaryField()[patchI] ==
1461                   - rho().boundaryField()[patchI]
1462                    *(g_.value()&(mesh().C().boundaryField()[patchI] - R0));
1463             }
1464         }
1465     }
1466 }
1467
1468
1469 void freeSurface::updateSurfaceFlux()
1470 {
1471     Phis() = fac::interpolate(Us()) & aMesh().Le();
1472 }
1473
1474
1475 void freeSurface::updateSurfactantConcentration()
1476 {
1477     if(!cleanInterface())
1478     {
1479         Info << "Correct surfactant concentration" << endl << flush;
1480
1481         updateSurfaceFlux();
1482
1483         // Crate and solve the surfactanta transport equation
1484         faScalarMatrix CsEqn
1485         (
1486             fam::ddt(surfactantConcentration())
1487           + fam::div(Phis(), surfactantConcentration())
1488           - fam::laplacian
1489             (
1490                 surfactant().surfactDiffusion(),
1491                 surfactantConcentration()
1492             )
1493         );
1494
1495
1496         if(surfactant().soluble())
1497         {
1498             const scalarField& C =
1499                 mesh().boundary()[aPatchID()]
1500                .lookupPatchField<volScalarField, scalar>("C");
1501
1502             areaScalarField Cb
1503             (
1504                 IOobject
1505                 (
1506                     "Cb",
1507                     DB().timeName(),
1508                     mesh(),
1509                     IOobject::NO_READ,
1510                     IOobject::NO_WRITE
1511                 ),
1512                 aMesh(),
1513                 dimensioned<scalar>("Cb", dimMoles/dimVolume, 0),
1514                 zeroGradientFaPatchScalarField::typeName
1515             );
1516
1517             Cb.internalField() = C;
1518             Cb.correctBoundaryConditions();
1519
1520             CsEqn +=
1521                 fam::Sp
1522                 (
1523                     surfactant().surfactAdsorptionCoeff()*Cb
1524                   + surfactant().surfactAdsorptionCoeff()
1525                    *surfactant().surfactDesorptionCoeff(),
1526                     surfactantConcentration()
1527                 )
1528               - surfactant().surfactAdsorptionCoeff()
1529                *Cb*surfactant().surfactSaturatedConc();
1530         }
1531
1532         CsEqn.solve();
1533
1534         Info << "Correct surface tension" << endl;
1535
1536         surfaceTension() =
1537             cleanInterfaceSurfTension()
1538           + surfactant().surfactR()
1539            *surfactant().surfactT()
1540            *surfactant().surfactSaturatedConc()
1541            *log(1.0 - surfactantConcentration()
1542            /surfactant().surfactSaturatedConc());
1543
1544         if(neg(min(surfaceTension().internalField())))
1545         {
1546             FatalErrorIn
1547             (
1548                 "void freeSurface::correctSurfactantConcentration()"
1549             )   << "Surface tension is negative"
1550                     << abort(FatalError);
1551         }
1552     }
1553 }
1554
1555
1556 void freeSurface::correctUsBoundaryConditions()
1557 {
1558     forAll(Us().boundaryField(), patchI)
1559     {
1560         if
1561         (
1562             Us().boundaryField()[patchI].type()
1563          == calculatedFaPatchVectorField::typeName
1564         )
1565         {
1566             vectorField& pUs = Us().boundaryField()[patchI];
1567
1568             pUs = Us().boundaryField()[patchI].patchInternalField();
1569
1570             label ngbPolyPatchID =
1571                 aMesh().boundary()[patchI].ngbPolyPatchIndex();
1572
1573             if(ngbPolyPatchID != -1)
1574             {
1575                 if
1576                 (
1577                     (
1578                         isA<slipFvPatchVectorField>
1579                         (
1580                             U().boundaryField()[ngbPolyPatchID]
1581                         )
1582                     )
1583                  ||
1584                     (
1585                         isA<symmetryFvPatchVectorField>
1586                         (
1587                             U().boundaryField()[ngbPolyPatchID]
1588                         )
1589                     )
1590                 )
1591                 {
1592                     vectorField N =
1593                         aMesh().boundary()[patchI].ngbPolyPatchFaceNormals();
1594
1595                     pUs -= N*(N&pUs);
1596                 }
1597             }
1598         }
1599     }
1600
1601     Us().correctBoundaryConditions();
1602 }
1603
1604
1605 vector freeSurface::totalPressureForce() const
1606 {
1607     const scalarField& S = aMesh().S();
1608
1609     const vectorField& n = aMesh().faceAreaNormals().internalField();
1610
1611     const scalarField& P = p().boundaryField()[aPatchID()];
1612
1613     vectorField pressureForces = S*P*n;
1614
1615     return gSum(pressureForces);
1616 }
1617
1618
1619 vector freeSurface::totalViscousForce() const
1620 {
1621     const scalarField& S = aMesh().S();
1622     const vectorField& n = aMesh().faceAreaNormals().internalField();
1623
1624     vectorField nGradU =
1625         U().boundaryField()[aPatchID()].snGrad();
1626
1627     vectorField viscousForces =
1628       - muFluidA().value()*S
1629        *(
1630             nGradU
1631           + (fac::grad(Us())().internalField()&n)
1632           - (n*fac::div(Us())().internalField())
1633         );
1634
1635     return gSum(viscousForces);
1636 }
1637
1638
1639 vector freeSurface::totalSurfaceTensionForce() const
1640 {
1641     const scalarField& S = aMesh().S();
1642
1643     const vectorField& n = aMesh().faceAreaNormals().internalField();
1644
1645     const scalarField& K = aMesh().faceCurvatures().internalField();
1646
1647     vectorField surfTensionForces(n.size(), vector::zero);
1648
1649     if(cleanInterface())
1650     {
1651         surfTensionForces =
1652             S*cleanInterfaceSurfTension().value()
1653            *fac::edgeIntegrate
1654             (
1655                 aMesh().Le()*aMesh().edgeLengthCorrection()
1656             )().internalField();
1657     }
1658     else
1659     {
1660         surfTensionForces *= surfaceTension().internalField()*K;
1661     }
1662
1663     return gSum(surfTensionForces);
1664 }
1665
1666
1667 void freeSurface::initializeControlPointsPosition()
1668 {
1669     scalarField deltaH = scalarField(aMesh().nFaces(), 0.0);
1670
1671     pointField displacement = pointDisplacement(deltaH);
1672
1673     const faceList& faces = aMesh().faces();
1674     const pointField& points = aMesh().points();
1675
1676
1677     pointField newPoints = points + displacement;
1678
1679     scalarField sweptVol(faces.size(), 0.0);
1680
1681     forAll(faces, faceI)
1682     {
1683         sweptVol[faceI] = -faces[faceI].sweptVol(points, newPoints);
1684     }
1685
1686     vectorField faceArea(faces.size(), vector::zero);
1687
1688     forAll (faceArea, faceI)
1689     {
1690         faceArea[faceI] = faces[faceI].normal(newPoints);
1691     }
1692
1693     forAll(deltaH, faceI)
1694     {
1695         deltaH[faceI] = sweptVol[faceI]/
1696             (faceArea[faceI] & facesDisplacementDir()[faceI]);
1697     }
1698
1699     displacement = pointDisplacement(deltaH);
1700 }
1701
1702
1703 scalar freeSurface::maxCourantNumber()
1704 {
1705     scalar CoNum = 0;
1706
1707     if(cleanInterface())
1708     {
1709         const scalarField& dE =aMesh().lPN();
1710
1711         CoNum = gMax
1712         (
1713             DB().deltaT().value()/
1714             sqrt
1715             (
1716                 rhoFluidA().value()*dE*dE*dE/
1717                 2.0/M_PI/(cleanInterfaceSurfTension().value() + SMALL)
1718             )
1719         );
1720     }
1721     else
1722     {
1723         scalarField sigmaE =
1724             linearEdgeInterpolate(surfaceTension())().internalField()
1725           + SMALL;
1726
1727         const scalarField& dE =aMesh().lPN();
1728
1729         CoNum = gMax
1730         (
1731             DB().deltaT().value()/
1732             sqrt
1733             (
1734                 rhoFluidA().value()*dE*dE*dE/
1735                 2.0/M_PI/sigmaE
1736             )
1737         );
1738     }
1739
1740     return CoNum;
1741 }
1742
1743
1744 void freeSurface::updateProperties()
1745 {
1746     muFluidA_ = dimensionedScalar(this->lookup("muFluidA"));
1747
1748     muFluidB_ = dimensionedScalar(this->lookup("muFluidB"));
1749
1750     rhoFluidA_ = dimensionedScalar(this->lookup("rhoFluidA"));
1751
1752     rhoFluidB_ = dimensionedScalar(this->lookup("rhoFluidB"));
1753
1754     g_ = dimensionedVector(this->lookup("g"));
1755
1756     cleanInterfaceSurfTension_ =
1757         dimensionedScalar(this->lookup("surfaceTension"));
1758 }
1759
1760
1761 void freeSurface::writeVTK() const
1762 {
1763     aMesh().patch().writeVTK
1764     (
1765         DB().timePath()/"freeSurface",
1766         aMesh().patch(),
1767         aMesh().patch().points()
1768     );
1769 }
1770
1771
1772 void freeSurface::writeVTKControlPoints()
1773 {
1774     // Write patch and points into VTK
1775     fileName name(DB().timePath()/"freeSurfaceControlPoints");
1776     OFstream mps(name + ".vtk");
1777
1778     mps << "# vtk DataFile Version 2.0" << nl
1779         << name << ".vtk" << nl
1780         << "ASCII" << nl
1781         << "DATASET POLYDATA" << nl
1782         << "POINTS " << controlPoints().size() << " float" << nl;
1783
1784     forAll(controlPoints(), pointI)
1785     {
1786         mps << controlPoints()[pointI].x() << ' '
1787             << controlPoints()[pointI].y() << ' '
1788             << controlPoints()[pointI].z() << nl;
1789     }
1790
1791     // Write vertices
1792     mps << "VERTICES " << controlPoints().size() << ' '
1793         << controlPoints().size()*2 << nl;
1794
1795     forAll(controlPoints(), pointI)
1796     {
1797         mps << 1 << ' ' << pointI << nl;
1798     }
1799 }
1800
1801
1802 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1803
1804 } // End namespace Foam
1805
1806 // ************************************************************************* //