Foam::S_MDCPP::S_MDCPP
(
    const word& name,
    const volVectorField& U,
    const surfaceScalarField& phi,
    const dictionary& dict
)
    :
    viscoelasticLaw(name, U, phi),
    tau_
    (
       IOobject
       (
           "tau" + name,
           U.time().timeName(),
           U.mesh(),
           IOobject::MUST_READ,
           IOobject::AUTO_WRITE
       ),
       U.mesh()
    ),
    I_
    (
       dimensionedSymmTensor
       (
           "I",
           dimensionSet(0, 0, 0, 0, 0, 0, 0),
           symmTensor
           (
               1, 0, 0,
               1, 0,
               1
           )
       )
    ),
    rho_(dict.lookup("rho")),
    etaS_(dict.lookup("etaS")),
    etaP_(dict.lookup("etaP")),
    zeta_(dict.lookup("zeta")),
    lambdaOb_(dict.lookup("lambdaOb")),
    lambdaOs_(dict.lookup("lambdaOs")),
    q_(dict.lookup("q"))
{}

continuousGasKEpsilon<BasicTurbulenceModel>::continuousGasKEpsilon
(
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName,
    const word& type
)
:
    kEpsilon<BasicTurbulenceModel>
    (
        alpha,
        rho,
        U,
        alphaPhi,
        phi,
        transport,
        propertiesName,
        type
    ),

    liquidTurbulencePtr_(NULL),

    nutEff_
    (
        IOobject
        (
            IOobject::groupName("nutEff", U.group()),
            this->runTime_.timeName(),
            this->mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        this->nut_
    ),

    alphaInversion_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "alphaInversion",
            this->coeffDict_,
            0.7
        )
    )
{
    if (type == typeName)
    {
        kEpsilon<BasicTurbulenceModel>::correctNut();
        this->printCoeffs(type);
    }
}

Foam::autoPtr<Foam::SRF::SRFModel> Foam::SRF::SRFModel::New
(
    const volVectorField& Urel
)
{
    // get model name, but do not register the dictionary
    // otherwise it is registered in the database twice
    const word modelType
    (
        IOdictionary
        (
            IOobject
            (
                "SRFProperties",
                Urel.time().constant(),
                Urel.db(),
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE,
                false
            )
        ).lookup("SRFModel")
    );

    Info<< "Selecting SRFModel " << modelType << endl;

    dictionaryConstructorTable::iterator cstrIter =
        dictionaryConstructorTablePtr_->find(modelType);

    if (cstrIter == dictionaryConstructorTablePtr_->end())
    {
        FatalErrorIn
        (
            "SRFModel::New(const fvMesh&)"
        )   << "Unknown SRFModel type "
            << modelType << nl << nl
            << "Valid SRFModel types are :" << nl
            << dictionaryConstructorTablePtr_->sortedToc()
            << exit(FatalError);
    }

    return autoPtr<SRFModel>(cstrIter()(Urel));
}

Foam::scalar Foam::fv::patchMeanVelocityForce::magUbarAve
(
    const volVectorField& U
) const
{
    vector2D sumAmagUsumA
    (
        sum
        (
            (flowDir_ & U.boundaryField()[patchi_])
           *mesh_.boundary()[patchi_].magSf()
        ),
        sum(mesh_.boundary()[patchi_].magSf())
    );


    // If the mean velocity force is applied to a cyclic patch
    // for parallel runs include contributions from processorCyclic patches
    // generated from the decomposition of the cyclic patch
    const polyBoundaryMesh& patches = mesh_.boundaryMesh();

    if (Pstream::parRun() && isA<cyclicPolyPatch>(patches[patchi_]))
    {
        labelList processorCyclicPatches
        (
            processorCyclicPolyPatch::patchIDs(patch_, patches)
        );

        forAll(processorCyclicPatches, pcpi)
        {
            const label patchi = processorCyclicPatches[pcpi];

            sumAmagUsumA.x() +=
                sum
                (
                    (flowDir_ & U.boundaryField()[patchi])
                   *mesh_.boundary()[patchi].magSf()
                );

            sumAmagUsumA.y() += sum(mesh_.boundary()[patchi].magSf());
        }
    }

viscoelasticModel::viscoelasticModel
(
    const volScalarField& alpha,
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    IOdictionary
    (
        IOobject
        (
            "viscoelasticProperties",
            U.time().constant(),
            U.db(),
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    ),
    lawPtr_(viscoelasticLaw::New(word::null, alpha, U, phi, subDict("rheology")))
{}

void dilute::setVectorAverage
(
    volVectorField& field,
    double**& value,
    double**& weight,
    volScalarField& weightField,
    double**const& mask,
    double**const& weight2,         //allows the specification of a 2nd weight field
    bool weightWithWeight2    //switch to activate 2nd weight field
) const
{
    label cellI;
    vector valueVec;
    scalar weightP;

    if(weightWithWeight2) //use weight2, e.g., mass-averaged - has no effect, just weight is DIFFERENT!
    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        for(int subCell=0;subCell<particleCloud_.cellsPerParticle()[index][0];subCell++)
        {
            cellI = particleCloud_.cellIDs()[index][subCell];
            if (cellI >= 0)
            {
                for(int i=0;i<3;i++) valueVec[i] = value[index][i];
                weightP = weight[index][subCell]*weight2[index][subCell];
                weightField[cellI] += weightP;
                if(weightP > 0) field[cellI] = valueVec; //field[cellI] = valueVec/weightP;
            }
        }
    }
    else //standard, i.e., volume-averaged - has no effect, just weight is DIFFERENT!
    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        for(int subCell=0;subCell<particleCloud_.cellsPerParticle()[index][0];subCell++)
        {
            //Info << "subCell=" << subCell << endl;
            cellI = particleCloud_.cellIDs()[index][subCell];

            if (cellI >= 0)
            {
                for(int i=0;i<3;i++) valueVec[i] = value[index][i];
                weightP = weight[index][subCell];
                weightField[cellI] += weightP;
                if(weightP > 0) field[cellI] = valueVec; //field[cellI] = valueVec/weightP;
                else Warning << "!!! W A R N I N G --- weightP <= 0" << endl;
            }
        }
    }

    // correct cell values to patches
    field.correctBoundaryConditions();
}

scaleSimilarity::scaleSimilarity
(
    const volVectorField& U,
    const surfaceScalarField& phi,
    transportModel& transport
)
:
    LESModel(typeName, U, phi, transport),
    filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
    filter_(filterPtr_())
{
    printCoeffs();
}

LESModel::LESModel
(
    const word& type,
    const volScalarField& rho,
    const volVectorField& U,
    const surfaceScalarField& phi,
    const basicThermo& thermoPhysicalModel,
    const word& turbulenceModelName
)
:
    turbulenceModel(rho, U, phi, thermoPhysicalModel, turbulenceModelName),

    IOdictionary
    (
        IOobject
        (
            "LESProperties",
            U.time().constant(),
            U.db(),
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    ),

    printCoeffs_(lookupOrDefault<Switch>("printCoeffs", false)),
    coeffDict_(subOrEmptyDict(type + "Coeffs")),

    kMin_("kMin", sqr(dimVelocity), SMALL),

    delta_(LESdelta::New("delta", U.mesh(), *this))
{
    kMin_.readIfPresent(*this);

    // Force the construction of the mesh deltaCoeffs which may be needed
    // for the construction of the derived models and BCs
    mesh_.deltaCoeffs();
}

void Foam::twoStrokeEngine::setBoundaryVelocity(volVectorField& U)
{
    vector pistonVel = piston().cs().axis()*engTime().pistonSpeed().value();

    //  On the piston movingWallVelocity is used.
    // There is no need to update the piston velocity

    forAll (scavInPortPatches_, patchi)
    {
        const label curPatchID =
            boundaryMesh().findPatchID(scavInPortPatches_[patchi]);

        U.boundaryField()[curPatchID] == pistonVel;
    }
}

void Foam::porosityModels::powerLaw::calcForce
(
    const volVectorField& U,
    const volScalarField& rho,
    const volScalarField& mu,
    vectorField& force
) const
{
    scalarField Udiag(U.size(), 0.0);
    const scalarField& V = mesh_.V();

    apply(Udiag, V, rho, U);

    force = Udiag*U;
}

dynOneEqEddy::dynOneEqEddy
(
    const volVectorField& U,
    const surfaceScalarField& phi,
    transportModel& transport,
    const word& turbulenceModelName,
    const word& modelName
)
:
    LESModel(modelName, U, phi, transport, turbulenceModelName),
    GenEddyVisc(U, phi, transport),

    k_
    (
        IOobject
        (
            "k",
            runTime_.timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),

    simpleFilter_(U.mesh()),
    filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
    filter_(filterPtr_())
{
    bound(k_, kMin_);

    const volScalarField KK(0.5*(filter_(magSqr(U)) - magSqr(filter_(U))));
    updateSubGridScaleFields(symm(fvc::grad(U)), KK);

    printCoeffs();
}

void dense::setVectorAverage
(
    volVectorField& field,
    double**& value,
    double**& weight,
    volScalarField& weightField,
    double**const& mask
) const
{
    label cellI;
    vector valueVec;
    scalar weightP;

    for(int index=0; index< particleCloud_.numberOfParticles(); index++)
    {
        if(mask[index][0])
        {
            for(int subCell=0;subCell<particleCloud_.voidFractionM().cellsPerParticle()[index][0];subCell++)
            {
                cellI = particleCloud_.cellIDs()[index][subCell];

                if (cellI >= 0)
                {
                    for(int i=0;i<3;i++) valueVec[i] = value[index][i];
                    weightP = weight[index][subCell];

                    // first entry in this cell
                    if(weightField[cellI] == 0)
                    {
                        field[cellI] = valueVec;
                        weightField[cellI] = weightP;
                    }
                    else
                    {
                        field[cellI] = (field[cellI]*weightField[cellI]+valueVec*weightP)/(weightField[cellI]+weightP);
                        weightField[cellI] += weightP;
                    }
                }
            }//forAllSubPoints
        }
    }

    // correct cell values to patches
    field.correctBoundaryConditions();
}

dynOneEqEddy::dynOneEqEddy
(
    const volScalarField& rho,
    const volVectorField& U,
    const surfaceScalarField& phi,
    const basicThermo& thermoPhysicalModel
)
:
    LESModel(typeName, rho, U, phi, thermoPhysicalModel),
    GenEddyVisc(rho, U, phi, thermoPhysicalModel),

    filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
    filter_(filterPtr_())
{
    updateSubGridScaleFields(dev(symm(fvc::grad(U))));

    printCoeffs();
}

Foam::nonlinearEddyViscosity<BasicTurbulenceModel>::nonlinearEddyViscosity
(
    const word& modelName,
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName
)
:
    eddyViscosity<BasicTurbulenceModel>
    (
        modelName,
        alpha,
        rho,
        U,
        alphaRhoPhi,
        phi,
        transport,
        propertiesName
    ),

    nonlinearStress_
    (
        IOobject
        (
            IOobject::groupName("nonlinearStress", U.group()),
            this->runTime_.timeName(),
            this->mesh_
        ),
        this->mesh_,
        dimensionedSymmTensor
        (
            "nonlinearStress",
            sqr(dimVelocity),
            symmTensor::zero
        )
    )
{}

void Foam::porosityModels::solidification::apply
(
    tensorField& AU,
    const RhoFieldType& rho,
    const volVectorField& U
) const
{
    if (alphaName_ == "none")
    {
        return apply(AU, geometricOneField(), rho, U);
    }
    else
    {
        const volScalarField& alpha = mesh_.lookupObject<volScalarField>
        (
            IOobject::groupName(alphaName_, U.group())
        );

        return apply(AU, alpha, rho, U);
    }
}

Foam::mixtureViscosityModels::slurry::slurry
(
    const word& name,
    const dictionary& viscosityProperties,
    const volVectorField& U,
    const surfaceScalarField& phi,
    const word modelName
)
:
    mixtureViscosityModel(name, viscosityProperties, U, phi),
    alpha_
    (
        U.mesh().lookupObject<volScalarField>
        (
            IOobject::groupName
            (
                viscosityProperties.lookupOrDefault<word>("alpha", "alpha"),
                viscosityProperties.dictName()
            )
        )
    )
{}

Foam::eddyViscosity<BasicTurbulenceModel>::eddyViscosity
(
    const word& type,
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName
)
:
    linearViscousStress<BasicTurbulenceModel>
    (
        type,
        alpha,
        rho,
        U,
        alphaRhoPhi,
        phi,
        transport,
        propertiesName
    ),

    nut_
    (
        IOobject
        (
            IOobject::groupName("nut", U.group()),
            this->runTime_.timeName(),
            this->mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        this->mesh_
    )
{}

dynLagrangianCsBound::dynLagrangianCsBound
(
    const volVectorField& U,
    const surfaceScalarField& phi,
    transportModel& transport,
    const word& turbulenceModelName,
    const word& modelName
)
:
    LESModel(modelName, U, phi, transport, turbulenceModelName),
    GenEddyVisc(U, phi, transport),
    CsMin(0.07),
    CsMax(0.14),
    flm_
    (
        IOobject
        (
            "flm",
            runTime_.timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),
    fmm_
    (
        IOobject
        (
            "fmm",
            runTime_.timeName(),
            mesh_,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_
    ),
    Cs_
    (
        IOobject
        (
            "Cs",
            runTime_.timeName(),
            mesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        Foam::sqrt(flm_/fmm_)
    ),
    theta_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "theta",
            coeffDict_,
            1.5
        )
    ),
    simpleFilter_(U.mesh()),
    filterPtr_(LESfilter::New(U.mesh(), coeffDict())),
    filter_(filterPtr_()),
    flm0_("flm0", flm_.dimensions(), 0.0),
    fmm0_("fmm0", fmm_.dimensions(), VSMALL)
{
    updateSubGridScaleFields(fvc::grad(U));

    printCoeffs();
}

kOmegaSSTSASnew<BasicTurbulenceModel>::kOmegaSSTSASnew
(
    const alphaField& alpha,
    const rhoField& rho,
    const volVectorField& U,
    const surfaceScalarField& alphaRhoPhi,
    const surfaceScalarField& phi,
    const transportModel& transport,
    const word& propertiesName,
    const word& type
)
:
    kOmegaSST<BasicTurbulenceModel>
    (
        alpha,
        rho,
        U,
        alphaRhoPhi,
        phi,
        transport,
        propertiesName
    ),

    Cs_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "Cs",
            this->coeffDict_,
            0.262
        )
    ),
    kappa_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "kappa",
            this->coeffDict_,
            0.41
        )
    ),
    zeta2_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "zeta2",
            this->coeffDict_,
            1.47
        )
    ),
    sigmaPhi_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "sigmaPhi",
            this->coeffDict_,
            2.0/3.0
        )
    ),
    C_
    (
        dimensioned<scalar>::lookupOrAddToDict
        (
            "C",
            this->coeffDict_,
            2
        )
    ),

    delta_
    (
        LESdelta::New
        (
            IOobject::groupName("delta", U.group()),
            *this,
            this->coeffDict_
        )
    )
{}

// Construct from components
Foam::spray::spray
(
    const volVectorField& U,
    const volScalarField& rho,
    const volScalarField& p,
    const volScalarField& T,
    const basicMultiComponentMixture& composition,
    const PtrList<gasThermoPhysics>& gasProperties,
    const dictionary&,
    const dimensionedVector& g,
    bool readFields
)
:
    Cloud<parcel>(U.mesh(), false), // suppress className checking on positions
    runTime_(U.time()),
    time0_(runTime_.value()),
    mesh_(U.mesh()),
    rndGen_(label(0)),
    g_(g.value()),

    U_(U),
    rho_(rho),
    p_(p),
    T_(T),

    sprayProperties_
    (
        IOobject
        (
            "sprayProperties",
            U.time().constant(),
            U.db(),
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    ),

    ambientPressure_(p_.average().value()),
    ambientTemperature_(T_.average().value()),

    injectors_
    (
        IOobject
        (
            "injectorProperties",
            U.time().constant(),
            U.db(),
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        ),
        injector::iNew(U.time())
    ),
    atomization_
    (
        atomizationModel::New
        (
            sprayProperties_,
            *this
        )
    ),
    drag_
    (
        dragModel::New
        (
            sprayProperties_
        )
    ),
    evaporation_
    (
        evaporationModel::New
        (
            sprayProperties_
        )
    ),
    heatTransfer_
    (
        heatTransferModel::New
        (
            sprayProperties_
        )
    ),
    wall_
    (
        wallModel::New
        (
            sprayProperties_,
            U,
            *this
        )
    ),
    breakupModel_
    (
        breakupModel::New
        (
            sprayProperties_,
            *this
        )
    ),
    collisionModel_
//.........这里部分代码省略.........