C++ volScalarField::boundaryField方法代码示例

void mixtureKEpsilon<BasicTurbulenceModel>::correctInletOutlet
(
    volScalarField& vsf,
    const volScalarField& refVsf
) const
{
    volScalarField::GeometricBoundaryField& bf = vsf.boundaryField();
    const volScalarField::GeometricBoundaryField& refBf =
        refVsf.boundaryField();

    forAll(bf, patchi)
    {
        if
        (
            isA<inletOutletFvPatchScalarField>(bf[patchi])
         && isA<inletOutletFvPatchScalarField>(refBf[patchi])
        )
        {
            refCast<inletOutletFvPatchScalarField>
            (bf[patchi]).refValue() =
            refCast<const inletOutletFvPatchScalarField>
            (refBf[patchi]).refValue();
        }
    }
}

void Foam::bound(volScalarField& vsf, const dimensionedScalar& vsf0)
{
    scalar minVsf = min(vsf).value();

    if (minVsf < vsf0.value())
    {
        Info<< "bounding " << vsf.name()
            << ", min: " << gMin(vsf.internalField())
            << " max: " << gMax(vsf.internalField())
            << " average: " << gAverage(vsf.internalField())
            << endl;

        vsf.internalField() = max
        (
            max
            (
                vsf.internalField(),
                fvc::average(max(vsf, vsf0))().internalField()
                // Bug fix: was assuming bound on zero.  HJ, 25/Nov/2008
                *pos(vsf0.value() - vsf.internalField())
            ),
            vsf0.value()
        );

        vsf.correctBoundaryConditions();
        vsf.boundaryField() = max(vsf.boundaryField(), vsf0.value());
    }
}

void calcYPlus
(
    const TurbulenceModel& turbulenceModel,
    const fvMesh& mesh,
    const volVectorField& U,
    volScalarField& yPlus
)
{
    volScalarField::GeometricBoundaryField d = nearWallDist(mesh).y();

    const volScalarField::GeometricBoundaryField nutBf =
        turbulenceModel->nut()().boundaryField();

    const volScalarField::GeometricBoundaryField nuEffBf =
        turbulenceModel->nuEff()().boundaryField();

    const volScalarField::GeometricBoundaryField nuBf =
        turbulenceModel->nu()().boundaryField();

    const fvPatchList& patches = mesh.boundary();

    forAll(patches, patchi)
    {
        const fvPatch& patch = patches[patchi];

        if (isA<nutWallFunctionFvPatchScalarField>(nutBf[patchi]))
        {
            const nutWallFunctionFvPatchScalarField& nutPf =
                dynamic_cast<const nutWallFunctionFvPatchScalarField&>
                (
                    nutBf[patchi]
                );

            yPlus.boundaryField()[patchi] = nutPf.yPlus();
            const scalarField& Yp = yPlus.boundaryField()[patchi];

            Info<< "Patch " << patchi
                << " named " << nutPf.patch().name()
                << ", wall-function " << nutPf.type()
                << ", y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
                << " average: " << gAverage(Yp) << nl << endl;
        }
        else if (isA<wallFvPatch>(patch))
        {
            yPlus.boundaryField()[patchi] =
                d[patchi]
               *sqrt
                (
                    nuEffBf[patchi]
                   *mag(U.boundaryField()[patchi].snGrad())
                )/nuBf[patchi];
            const scalarField& Yp = yPlus.boundaryField()[patchi];

            Info<< "Patch " << patchi
                << " named " << patch.name()
                << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
                << " average: " << gAverage(Yp) << nl << endl;
        }
    }
}

void Foam::boundMinMax
(
    volScalarField& vsf,
    const dimensionedScalar& vsf0,
    const dimensionedScalar& vsf1
)
{
    scalar minVsf = min(vsf).value();
    scalar maxVsf = max(vsf).value();

    if (minVsf < vsf0.value() || maxVsf > vsf1.value())
    {
        Info<< "bounding " << vsf.name()
            << ", min: " << gMin(vsf.internalField())
            << " max: " << gMax(vsf.internalField())
            << " average: " << gAverage(vsf.internalField())
            << endl;
    }

    if (minVsf < vsf0.value())
    {
        vsf.internalField() = max
        (
            max
            (
                vsf.internalField(),
                fvc::average(max(vsf, vsf0))().internalField()
                *pos(vsf0.value() - vsf.internalField())
            ),
            vsf0.value()
        );

        vsf.correctBoundaryConditions();
        vsf.boundaryField() = max(vsf.boundaryField(), vsf0.value());
    }

    if (maxVsf > vsf1.value())
    {
        vsf.internalField() = min
        (
            min
            (
                vsf.internalField(),
                fvc::average(min(vsf, vsf1))().internalField()
                *neg(vsf1.value() - vsf.internalField())
                // This is needed when all values are above max
                // HJ, 18/Apr/2009
              + pos(vsf1.value() - vsf.internalField())*vsf1.value()
            ),
            vsf1.value()
        );

        vsf.correctBoundaryConditions();
        vsf.boundaryField() = min(vsf.boundaryField(), vsf1.value());
    }
}

void Foam::yPlusRAS::calcIncompressibleYPlus
(
    const fvMesh& mesh,
    volScalarField& yPlus
)
{
    typedef incompressible::nutWallFunctionFvPatchScalarField
        wallFunctionPatchField;

    const incompressible::RASModel& model =
        mesh.lookupObject<incompressible::RASModel>("RASProperties");

    const volScalarField nut(model.nut());
    const volScalarField::GeometricBoundaryField& nutPatches =
        nut.boundaryField();

    bool foundPatch = false;
    forAll(nutPatches, patchI)
    {
        if (isA<wallFunctionPatchField>(nutPatches[patchI]))
        {
            foundPatch = true;

            const wallFunctionPatchField& nutPw =
                dynamic_cast<const wallFunctionPatchField&>(nutPatches[patchI]);

            yPlus.boundaryField()[patchI] = nutPw.yPlus();
            const scalarField& Yp = yPlus.boundaryField()[patchI];

            scalar minYp = gMin(Yp);
            scalar maxYp = gMax(Yp);
            scalar avgYp = gAverage(Yp);

            if (log_)
            {
                Info<< "    patch " << nutPw.patch().name()
                    << " y+ : min = " << minYp << ", max = " << maxYp
                    << ", average = " << avgYp << nl;
            }

            if (Pstream::master())
            {
                file() << obr_.time().value() << token::TAB
                    << nutPw.patch().name() << token::TAB
                    << minYp << token::TAB << maxYp << token::TAB
                    << avgYp << endl;
            }
        }
    }

    if (log_ && !foundPatch)
    {
        Info<< "    no " << wallFunctionPatchField::typeName << " patches"
            << endl;
    }
}

void calcIncompressibleYPlus
(
    const fvMesh& mesh,
    const Time& runTime,
    const volVectorField& U,
    volScalarField& yPlus
)
{
    typedef incompressible::RASModels::nutWallFunctionFvPatchScalarField
        wallFunctionPatchField;

    #include "createPhi.H"

    singlePhaseTransportModel laminarTransport(U, phi);

    autoPtr<incompressible::RASModel> RASModel
    (
        incompressible::RASModel::New(U, phi, laminarTransport)
    );

    const volScalarField::GeometricBoundaryField nutPatches =
        RASModel->nut()().boundaryField();

    bool foundNutPatch = false;
    forAll(nutPatches, patchi)
    {
        if (isA<wallFunctionPatchField>(nutPatches[patchi]))
        {
            foundNutPatch = true;

            const wallFunctionPatchField& nutPw =
                dynamic_cast<const wallFunctionPatchField&>
                    (nutPatches[patchi]);

            yPlus.boundaryField()[patchi] = nutPw.yPlus();
            const scalarField& Yp = yPlus.boundaryField()[patchi];

            Info<< "Patch " << patchi
                << " named " << nutPw.patch().name()
                << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
                << " average: " << gAverage(Yp) << nl << endl;
        }
    }

    if (!foundNutPatch)
    {
        Info<< "    no " << wallFunctionPatchField::typeName << " patches"
            << endl;
    }
}

vectorField fieldOperations::
getWallPointMotion
(
  const fvMesh& mesh,
  const volScalarField& C, 
  const label movingPatchID
)
{
  // interpolate the concentration from cells to wall faces
  coupledPatchInterpolation patchInterpolator
  (
    mesh.boundaryMesh()[movingPatchID], mesh 
  );

  // concentration and normals on the faces
  scalarField pointCface = -C.boundaryField()[movingPatchID].snGrad();
  vectorField pointNface = mesh.boundaryMesh()[movingPatchID].faceNormals();
  
  scalarField motionC = patchInterpolator.faceToPointInterpolate(pointCface);
  vectorField motionN = patchInterpolator.faceToPointInterpolate(pointNface);
  
  // normalize point normals to 1
  forAll(motionN, ii) motionN[ii]/=mag(motionN[ii]);
  
  return motionC*motionN;
}

Foam::tmp<Foam::volScalarField>
Foam::laminarModels::generalizedNewtonianViscosityModels::strainRateFunction::
nu
(
    const volScalarField& nu0,
    const volScalarField& strainRate
) const
{
    tmp<volScalarField> tnu
    (
        volScalarField::New
        (
            IOobject::groupName(type() + ":nu", nu0.group()),
            nu0.mesh(),
            dimensionedScalar(dimViscosity, 0)
        )
    );

    tnu.ref().primitiveFieldRef() = strainRateFunction_->value(strainRate);

    volScalarField::Boundary& nuBf = tnu.ref().boundaryFieldRef();
    const volScalarField::Boundary& sigmaBf = strainRate.boundaryField();

    forAll(nuBf, patchi)
    {
        nuBf[patchi] = strainRateFunction_->value(sigmaBf[patchi]);
    }

tmp<GeometricField<Type, fvPatchField, volMesh> >
EulerLocalDdtScheme<Type>::fvcDdt
(
    const volScalarField& rho,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const objectRegistry& registry = this->mesh();

    // get access to the scalar beta[i]
    const scalarField& beta =
        registry.lookupObject<scalarField>(deltaTName_);

    volScalarField rDeltaT =
        1.0/(beta[0]*registry.lookupObject<volScalarField>(deltaTauName_));

    IOobject ddtIOobject
    (
        "ddt("+rho.name()+','+vf.name()+')',
        mesh().time().timeName(),
        mesh()
    );

    if (mesh().moving())
    {
        return tmp<GeometricField<Type, fvPatchField, volMesh> >
        (
            new GeometricField<Type, fvPatchField, volMesh>
            (
                ddtIOobject,
                mesh(),
                rDeltaT.dimensions()*rho.dimensions()*vf.dimensions(),
                rDeltaT.internalField()*
                (
                    rho.internalField()*vf.internalField()
                  - rho.oldTime().internalField()
                   *vf.oldTime().internalField()*mesh().V0()/mesh().V()
                ),
                rDeltaT.boundaryField()*
                (
                    rho.boundaryField()*vf.boundaryField()
                  - rho.oldTime().boundaryField()
                   *vf.oldTime().boundaryField()
                )
            )
        );
    }
    else
    {
        return tmp<GeometricField<Type, fvPatchField, volMesh> >
        (
            new GeometricField<Type, fvPatchField, volMesh>
            (
                ddtIOobject,
                rDeltaT*(rho*vf - rho.oldTime()*vf.oldTime())
            )
        );
    }
}

void Foam::compressibleTwoPhaseMixtureThermo::heBoundaryCorrection(volScalarField& h)
{
    volScalarField::GeometricBoundaryField& hbf = h.boundaryField();

    forAll(hbf, patchi)
    {
        if (isA<gradientEnergyFvPatchScalarField>(hbf[patchi]))
        {
            refCast<gradientEnergyFvPatchScalarField>(hbf[patchi]).gradient()
                = hbf[patchi].fvPatchField::snGrad();
        }
        else if (isA<mixedEnergyFvPatchScalarField>(hbf[patchi]))
        {
            refCast<mixedEnergyFvPatchScalarField>(hbf[patchi]).refGrad()
                = hbf[patchi].fvPatchField::snGrad();
        }
    }
}

void Foam::basicThermo::eBoundaryCorrection(volScalarField& e)
{
    volScalarField::GeometricBoundaryField& ebf = e.boundaryField();

    forAll(ebf, patchi)
    {
        if (isA<gradientInternalEnergyFvPatchScalarField>(ebf[patchi]))
        {
            refCast<gradientInternalEnergyFvPatchScalarField>(ebf[patchi])
                .gradient() = ebf[patchi].fvPatchField::snGrad();
        }
        else if (isA<mixedInternalEnergyFvPatchScalarField>(ebf[patchi]))
        {
            refCast<mixedInternalEnergyFvPatchScalarField>(ebf[patchi])
                .refGrad() = ebf[patchi].fvPatchField::snGrad();
        }
    }
}

void calcCompressibleYPlus
(
    const fvMesh& mesh,
    const Time& runTime,
    const volVectorField& U,
    volScalarField& yPlus
)
{
    typedef compressible::RASModels::mutkWallFunctionFvPatchScalarField
        wallFunctionPatchField;

    IOobject rhoHeader
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    );

    if (!rhoHeader.headerOk())
    {
        Info<< "    no rho field" << endl;
        return;
    }

    Info<< "Reading field rho\n" << endl;
    volScalarField rho(rhoHeader, mesh);

    #include "compressibleCreatePhi.H"

    autoPtr<basicThermo> pThermo
    (
        basicThermo::New(mesh)
    );
    basicThermo& thermo = pThermo();

    autoPtr<compressible::RASModel> RASModel
    (
        compressible::RASModel::New
        (
            rho,
            U,
            phi,
            thermo
        )
    );

    const volScalarField::GeometricBoundaryField mutPatches =
        RASModel->mut()().boundaryField();

    bool foundMutPatch = false;
    forAll(mutPatches, patchi)
    {
        if (isA<wallFunctionPatchField>(mutPatches[patchi]))
        {
            foundMutPatch = true;

            const wallFunctionPatchField& mutPw =
                dynamic_cast<const wallFunctionPatchField&>
                    (mutPatches[patchi]);

            yPlus.boundaryField()[patchi] = mutPw.yPlus();
            const scalarField& Yp = yPlus.boundaryField()[patchi];

            Info<< "Patch " << patchi
                << " named " << mutPw.patch().name()
                << " y+ : min: " << gMin(Yp) << " max: " << gMax(Yp)
                << " average: " << gAverage(Yp) << nl << endl;
        }
    }

    if (!foundMutPatch)
    {
        Info<< "    no " << wallFunctionPatchField::typeName << " patches"
            << endl;
    }
}

void Foam::yPlus::calcYPlus
(
    const TurbulenceModel& turbulenceModel,
    const fvMesh& mesh,
    volScalarField& yPlus
)
{
    volScalarField::GeometricBoundaryField d = nearWallDist(mesh).y();

    const volScalarField::GeometricBoundaryField nutBf =
        turbulenceModel.nut()().boundaryField();

    const volScalarField::GeometricBoundaryField nuEffBf =
        turbulenceModel.nuEff()().boundaryField();

    const volScalarField::GeometricBoundaryField nuBf =
        turbulenceModel.nu()().boundaryField();

    const fvPatchList& patches = mesh.boundary();

    forAll(patches, patchi)
    {
        const fvPatch& patch = patches[patchi];

        if (isA<nutWallFunctionFvPatchScalarField>(nutBf[patchi]))
        {
            const nutWallFunctionFvPatchScalarField& nutPf =
                dynamic_cast<const nutWallFunctionFvPatchScalarField&>
                (
                    nutBf[patchi]
                );

            yPlus.boundaryField()[patchi] = nutPf.yPlus();
            const scalarField& yPlusp = yPlus.boundaryField()[patchi];

            const scalar minYplus = gMin(yPlusp);
            const scalar maxYplus = gMax(yPlusp);
            const scalar avgYplus = gAverage(yPlusp);

            if (Pstream::master())
            {
                if (log_) Info
                    << "    patch " << patch.name()
                    << " y+ : min = " << minYplus << ", max = " << maxYplus
                    << ", average = " << avgYplus << nl;

                file() << obr_.time().value()
                    << token::TAB << patch.name()
                    << token::TAB << minYplus
                    << token::TAB << maxYplus
                    << token::TAB << avgYplus
                    << endl;
            }
        }
        else if (isA<wallFvPatch>(patch))
        {
            yPlus.boundaryField()[patchi] =
                d[patchi]
               *sqrt
                (
                    nuEffBf[patchi]
                   *mag(turbulenceModel.U().boundaryField()[patchi].snGrad())
                )/nuBf[patchi];
            const scalarField& yPlusp = yPlus.boundaryField()[patchi];

            const scalar minYplus = gMin(yPlusp);
            const scalar maxYplus = gMax(yPlusp);
            const scalar avgYplus = gAverage(yPlusp);

            if (Pstream::master())
            {
                if (log_) Info
                    << "    patch " << patch.name()
                    << " y+ : min = " << minYplus << ", max = " << maxYplus
                    << ", average = " << avgYplus << nl;

                file() << obr_.time().value()
                    << token::TAB << patch.name()
                    << token::TAB << minYplus
                    << token::TAB << maxYplus
                    << token::TAB << avgYplus
                    << endl;
            }
        }
    }
}

void Foam::yPlusLES::calcIncompressibleYPlus
(
    const fvMesh& mesh,
    const volVectorField& U,
    volScalarField& yPlus
)
{
    const incompressible::LESModel& model =
        mesh.lookupObject<incompressible::LESModel>("LESProperties");

    volScalarField::GeometricBoundaryField d = nearWallDist(mesh).y();
    volScalarField nuEff(model.nuEff());

    const fvPatchList& patches = mesh.boundary();

    const volScalarField nuLam(model.nu());

    bool foundPatch = false;
    forAll(patches, patchI)
    {
        const fvPatch& currPatch = patches[patchI];

        if (isA<wallFvPatch>(currPatch))
        {
            foundPatch = true;
            yPlus.boundaryField()[patchI] =
                d[patchI]
               *sqrt
                (
                    nuEff.boundaryField()[patchI]
                   *mag(U.boundaryField()[patchI].snGrad())
                )
               /nuLam.boundaryField()[patchI];

            const scalarField& Yp = yPlus.boundaryField()[patchI];

            scalar minYp = gMin(Yp);
            scalar maxYp = gMax(Yp);
            scalar avgYp = gAverage(Yp);

            if (log_)
            {
                Info<< "    patch " << currPatch.name()
                    << " y+ : min = " << minYp << ", max = " << maxYp
                    << ", average = " << avgYp << nl;
            }

            if (Pstream::master())
            {
                file() << obr_.time().value() << token::TAB
                    << currPatch.name() << token::TAB
                    << minYp << token::TAB << maxYp << token::TAB
                    << avgYp << endl;
            }
        }
    }

    if (log_ && !foundPatch)
    {
        Info<< "    no " << wallFvPatch::typeName << " patches" << endl;
    }
}

tmp<GeometricField<Type, fvPatchField, volMesh>>
EulerD2dt2Scheme<Type>::fvcD2dt2
(
    const volScalarField& rho,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    dimensionedScalar rDeltaT2 =
        4.0/sqr(mesh().time().deltaT() + mesh().time().deltaT0());

    IOobject d2dt2IOobject
    (
        "d2dt2("+rho.name()+','+vf.name()+')',
        mesh().time().timeName(),
        mesh(),
        IOobject::NO_READ,
        IOobject::NO_WRITE
    );

    scalar deltaT = mesh().time().deltaTValue();
    scalar deltaT0 = mesh().time().deltaT0Value();

    scalar coefft   = (deltaT + deltaT0)/(2*deltaT);
    scalar coefft00 = (deltaT + deltaT0)/(2*deltaT0);

    if (mesh().moving())
    {
        scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
        scalar quarterRdeltaT2 = 0.25*rDeltaT2.value();

        const scalarField VV0rhoRho0
        (
            (mesh().V() + mesh().V0())
          * (rho.primitiveField() + rho.oldTime().primitiveField())
        );

        const scalarField V0V00rho0Rho00
        (
            (mesh().V0() + mesh().V00())
          * (
                rho.oldTime().primitiveField()
              + rho.oldTime().oldTime().primitiveField()
            )
        );

        return tmp<GeometricField<Type, fvPatchField, volMesh>>
        (
            new GeometricField<Type, fvPatchField, volMesh>
            (
                d2dt2IOobject,
                mesh(),
                rDeltaT2.dimensions()*rho.dimensions()*vf.dimensions(),
                quarterRdeltaT2*
                (
                    coefft*VV0rhoRho0*vf.primitiveField()

                  - (coefft*VV0rhoRho0 + coefft00*V0V00rho0Rho00)
                   *vf.oldTime().primitiveField()

                  + (coefft00*V0V00rho0Rho00)
                   *vf.oldTime().oldTime().primitiveField()
                )/mesh().V(),
                halfRdeltaT2*
                (
                    coefft
                   *(rho.boundaryField() + rho.oldTime().boundaryField())
                   *vf.boundaryField()

                  - (
                        coefft
                       *(
                           rho.boundaryField()
                         + rho.oldTime().boundaryField()
                        )
                      + coefft00
                       *(
                           rho.oldTime().boundaryField()
                         + rho.oldTime().oldTime().boundaryField()
                        )
                    )*vf.oldTime().boundaryField()

                  + coefft00
                   *(
                       rho.oldTime().boundaryField()
                     + rho.oldTime().oldTime().boundaryField()
                    )*vf.oldTime().oldTime().boundaryField()
                )
            )
        );
    }
    else
    {
        dimensionedScalar halfRdeltaT2 = 0.5*rDeltaT2;

        const volScalarField rhoRho0(rho + rho.oldTime());
        const volScalarField rho0Rho00(rho.oldTime() +rho.oldTime().oldTime());

        return tmp<GeometricField<Type, fvPatchField, volMesh>>
        (
            new GeometricField<Type, fvPatchField, volMesh>
//.........这里部分代码省略.........