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

tmp<errorEstimate<Type> >
div
(
    const surfaceScalarField& flux,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = vf.mesh();

    const scalarField& vols = mesh.V();
    const surfaceVectorField& faceCentres = mesh.Cf();
    const volVectorField& cellCentres = mesh.C();
    const fvPatchList& patches = mesh.boundary();
    const unallocLabelList& owner = mesh.owner();
    const unallocLabelList& neighbour = mesh.neighbour();

    Field<Type> res(vols.size(), pTraits<Type>::zero);
    scalarField aNorm(vols.size(), 0.0);

    // Get sign of flux
    const surfaceScalarField signF = pos(flux);

    // Calculate gradient of the solution
    GeometricField
    <
        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
    > gradVf = fvc::grad(vf);

    // Internal faces
    forAll (owner, faceI)
    {
        // Calculate the centre of the face
        const vector& curFaceCentre = faceCentres[faceI];

        // Owner
        vector ownD = curFaceCentre - cellCentres[owner[faceI]];

        // Subtract convection
        res[owner[faceI]] -=
        (
            vf[owner[faceI]]
          + (ownD & gradVf[owner[faceI]])
        )*flux[faceI];

        aNorm[owner[faceI]] += signF[faceI]*flux[faceI];

        // Neighbour
        vector neiD = curFaceCentre - cellCentres[neighbour[faceI]];

        // Subtract convection
        res[neighbour[faceI]] +=
        (
            vf[neighbour[faceI]]
          + (neiD & gradVf[neighbour[faceI]])
        )*flux[faceI];

        aNorm[neighbour[faceI]] -= (1.0 - signF[faceI])*flux[faceI];
    }

    forAll (patches, patchI)
    {
        const vectorField& patchFaceCentres =
            faceCentres.boundaryField()[patchI];

        const scalarField& patchFlux = flux.boundaryField()[patchI];
        const scalarField& patchSignFlux = signF.boundaryField()[patchI];

        const labelList& fCells = patches[patchI].faceCells();

        forAll (fCells, faceI)
        {
            vector d =
                patchFaceCentres[faceI] - cellCentres[fCells[faceI]];

            // Subtract convection
            res[fCells[faceI]] -=
            (
                vf[fCells[faceI]]
              + (d & gradVf[fCells[faceI]])
            )*patchFlux[faceI];

            aNorm[fCells[faceI]] += patchSignFlux[faceI]*patchFlux[faceI];
        }
    }

scalar fieldOperations::
getScalingFlowRate(const surfaceScalarField& phi)
{
  return mag( gSum(phi.boundaryField()[scalingPatchID_]) );
}

void Foam::MULES::limiter
(
    scalarField& allLambda,
    const RdeltaTType& rDeltaT,
    const RhoType& rho,
    const volScalarField& psi,
    const surfaceScalarField& phiBD,
    const surfaceScalarField& phiCorr,
    const SpType& Sp,
    const SuType& Su,
    const scalar psiMax,
    const scalar psiMin
)
{
    const scalarField& psiIf = psi;
    const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();

    const fvMesh& mesh = psi.mesh();

    const dictionary& MULEScontrols = mesh.solverDict(psi.name());

    label nLimiterIter
    (
        MULEScontrols.lookupOrDefault<label>("nLimiterIter", 3)
    );

    scalar smoothLimiter
    (
        MULEScontrols.lookupOrDefault<scalar>("smoothLimiter", 0)
    );

    const scalarField& psi0 = psi.oldTime();

    const labelUList& owner = mesh.owner();
    const labelUList& neighb = mesh.neighbour();
    tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
    const scalarField& V = tVsc();

    const scalarField& phiBDIf = phiBD;
    const surfaceScalarField::GeometricBoundaryField& phiBDBf =
        phiBD.boundaryField();

    const scalarField& phiCorrIf = phiCorr;
    const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
        phiCorr.boundaryField();

    slicedSurfaceScalarField lambda
    (
        IOobject
        (
            "lambda",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimless,
        allLambda,
        false   // Use slices for the couples
    );

    scalarField& lambdaIf = lambda;
    surfaceScalarField::GeometricBoundaryField& lambdaBf =
        lambda.boundaryField();

    scalarField psiMaxn(psiIf.size(), psiMin);
    scalarField psiMinn(psiIf.size(), psiMax);

    scalarField sumPhiBD(psiIf.size(), 0.0);

    scalarField sumPhip(psiIf.size(), VSMALL);
    scalarField mSumPhim(psiIf.size(), VSMALL);

    forAll(phiCorrIf, facei)
    {
        label own = owner[facei];
        label nei = neighb[facei];

        psiMaxn[own] = max(psiMaxn[own], psiIf[nei]);
        psiMinn[own] = min(psiMinn[own], psiIf[nei]);

        psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]);
        psiMinn[nei] = min(psiMinn[nei], psiIf[own]);

        sumPhiBD[own] += phiBDIf[facei];
        sumPhiBD[nei] -= phiBDIf[facei];

        scalar phiCorrf = phiCorrIf[facei];

        if (phiCorrf > 0.0)
        {
            sumPhip[own] += phiCorrf;
            mSumPhim[nei] += phiCorrf;
        }
        else
        {
            mSumPhim[own] -= phiCorrf;
            sumPhip[nei] -= phiCorrf;
//.........这里部分代码省略.........

void Foam::MULES::limiter
(
    scalargpuField& allLambda,
    const RdeltaTType& rDeltaT,
    const RhoType& rho,
    const volScalarField& psi,
    const surfaceScalarField& phiBD,
    const surfaceScalarField& phiCorr,
    const SpType& Sp,
    const SuType& Su,
    const scalar psiMax,
    const scalar psiMin,
    const label nLimiterIter
)
{
    const scalargpuField& psiIf = psi.getField();
    const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField();

    const scalargpuField& psi0 = psi.oldTime();

    const fvMesh& mesh = psi.mesh();

    const labelgpuList& owner = mesh.owner();
    const labelgpuList& neighb = mesh.neighbour();
    const labelgpuList& losort = mesh.lduAddr().losortAddr();

    const labelgpuList& ownStart = mesh.lduAddr().ownerStartAddr();
    const labelgpuList& losortStart = mesh.lduAddr().losortStartAddr();

    tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc();
    const scalargpuField& V = tVsc().getField();

    const scalargpuField& phiBDIf = phiBD;
    const surfaceScalarField::GeometricBoundaryField& phiBDBf =
        phiBD.boundaryField();

    const scalargpuField& phiCorrIf = phiCorr;
    const surfaceScalarField::GeometricBoundaryField& phiCorrBf =
        phiCorr.boundaryField();

    slicedSurfaceScalarField lambda
    (
        IOobject
        (
            "lambda",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimless,
        allLambda,
        false   // Use slices for the couples
    );

    scalargpuField& lambdaIf = lambda;
    surfaceScalarField::GeometricBoundaryField& lambdaBf =
        lambda.boundaryField();

    scalargpuField psiMaxn(psiIf.size(), psiMin);
    scalargpuField psiMinn(psiIf.size(), psiMax);

    scalargpuField sumPhiBD(psiIf.size(), 0.0);

    scalargpuField sumPhip(psiIf.size(), VSMALL);
    scalargpuField mSumPhim(psiIf.size(), VSMALL);

    thrust::for_each
    (
        thrust::make_counting_iterator(0),
        thrust::make_counting_iterator(0)+psiIf.size(),
        limiterMULESFunctor
        (
            owner.data(),
            neighb.data(),
            ownStart.data(),
            losortStart.data(),
            losort.data(),
            psiIf.data(),
            phiBDIf.data(),
            phiCorrIf.data(),
            psiMaxn.data(),
            psiMinn.data(),
            sumPhiBD.data(),
            sumPhip.data(),
            mSumPhim.data()
        )
    );


    forAll(phiCorrBf, patchi)
    {
        const fvPatchScalarField& psiPf = psiBf[patchi];
        const scalargpuField& phiBDPf = phiBDBf[patchi];
        const scalargpuField& phiCorrPf = phiCorrBf[patchi];

        const labelgpuList& pcells = mesh.lduAddr().patchSortCells(patchi);

//.........这里部分代码省略.........