Palabos源码: computeEquilibria

template<typename T, template<typename U> class Descriptor> void Dynamics<T,Descriptor>::computeEquilibria (
        Array<T,Descriptor<T>::q>& fEq,  T rhoBar, Array<T,Descriptor<T>::d> const& j, T jSqr, T thetaBar ) const {
    for (int iPop=0; iPop<Descriptor<T>::q; ++iPop) {
        fEq[iPop] = computeEquilibrium(iPop, rhoBar, j, jSqr, thetaBar);
    } }

dataAnalysisWrapper3D.hh

template<typename T, template<typename U> class Descriptor>
void computeRhoBarJ( MultiBlockLattice3D<T,Descriptor>& lattice,
                     MultiScalarField3D<T>& rhoBar, MultiTensorField3D<T,3>& j, Box3D domain )
{
    std::vector<MultiBlock3D*> fields;
    fields.push_back(&lattice);
    fields.push_back(&rhoBar);
    fields.push_back(&j);
    applyProcessingFunctional (
            new BoxRhoBarJfunctional3D<T,Descriptor>, domain, fields );
}

dataAnalysisFunctional3D.hh

template<typename T, template<typename U> class Descriptor> 
BoxEquilibriumFunctional3D<T,Descriptor>::BoxEquilibriumFunctional3D(plint iComponent_)
    : iComponent(iComponent_)
{ }

template<typename T, template<typename U> class Descriptor> 
void BoxEquilibriumFunctional3D<T,Descriptor>::process (
        Box3D domain, BlockLattice3D<T,Descriptor>& lattice, ScalarField3D<T>& equilibrium)
{
    Dot3D offset = computeRelativeDisplacement(lattice, equilibrium);
    for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
        for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
            for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
                T rhoBar;
                Array<T,Descriptor<T>::d> j;
                Cell<T,Descriptor> const& cell = lattice.get(iX,iY,iZ);
                cell.getDynamics().computeRhoBarJ(cell, rhoBar, j);
                T jSqr = normSqr(j);
                equilibrium.get(iX+offset.x,iY+offset.y,iZ+offset.z) =
                    cell.computeEquilibrium(iComponent, rhoBar, j, jSqr);
            }
        }
    }
}

template<typename T, template<typename U> class Descriptor> 
BoxEquilibriumFunctional3D<T,Descriptor>* BoxEquilibriumFunctional3D<T,Descriptor>::clone() const
{
    return new BoxEquilibriumFunctional3D<T,Descriptor>(*this);
}

template<typename T, template<typename U> class Descriptor> 
void BoxEquilibriumFunctional3D<T,Descriptor>::getTypeOfModification(std::vector<modif::ModifT>& modified) const {
    modified[0] = modif::nothing;
    modified[1] = modif::staticVariables;
}

template<typename T, template<typename U> class Descriptor> 
BlockDomain::DomainT BoxEquilibriumFunctional3D<T,Descriptor>::appliesTo() const {
    return BlockDomain::bulk;
}

Dynamics.h

/// Define density, if possible. Does nothing by default.
virtual void defineDensity(Cell<T,Descriptor>& cell, T density);

/// Define velocity, if possible. Does nothing by default.
virtual void defineVelocity(Cell<T,Descriptor>& cell, Array<T,Descriptor<T>::d> const& u);

Dynamics.hh

template<typename T, template<typename U> class Descriptor>
void Dynamics<T,Descriptor>::defineDensity(Cell<T,Descriptor>& cell, T density) { }

/** This method does nothing by default, unless overloaded in
 *  a deriving class. This is OK, because it's the behavior
 *  users will expect: if they use "defineSomething" on a cell,
 *  they want it to define the boundary condition on boundary
 *  nodes, but nothing to happen in bulk nodes. In this way,
 *  they can be lazy and loop over the whole domain instead
 *  of tracking boundary nodes explicitly.
 **/
template<typename T, template<typename U> class Descriptor>
void Dynamics<T,Descriptor>::defineVelocity(Cell<T,Descriptor>& cell, Array<T,Descriptor<T>::d> const& u) { }

datainitializerGenerics3D.h

template<typename T, template<typename U> class Descriptor, class RhoUFunction>
IniCustomEquilibriumFunctional3D<T,Descriptor,RhoUFunction>::
    IniCustomEquilibriumFunctional3D(RhoUFunction f_)
        : f(f_),
          velocityScale( (T)1 )
{ }

template<typename T, template<typename U> class Descriptor, class RhoUFunction>
void IniCustomEquilibriumFunctional3D<T,Descriptor,RhoUFunction>::execute (
        plint iX, plint iY, plint iZ, Cell<T,Descriptor>& cell ) const
{
    Array<T,Descriptor<T>::d> j;
    T rho;
    f(iX, iY, iZ, rho, j);
    Array<T,Descriptor<T>::d> force;
    force[0] = getExternalForceComponent(cell, 0);
    force[1] = getExternalForceComponent(cell, 1);
    force[2] = getExternalForceComponent(cell, 2);
    j[0] = velocityScale * rho * (j[0] - (T) 0.5 * force[0]);
    j[1] = velocityScale * rho * (j[1] - (T) 0.5 * force[1]);
    j[2] = velocityScale * rho * (j[2] - (T) 0.5 * force[2]);
    T jSqr = VectorTemplate<T,Descriptor>::normSqr(j);
    T rhoBar = Descriptor<T>::rhoBar(rho);
    for (plint iPop=0; iPop<Descriptor<T>::q; ++iPop) {
        cell[iPop] = cell.computeEquilibrium(iPop, rhoBar, j, jSqr);
    }
}

template<typename T, template<typename U> class Descriptor, class RhoUFunction>
IniCustomEquilibriumFunctional3D<T,Descriptor,RhoUFunction>*
    IniCustomEquilibriumFunctional3D<T,Descriptor,RhoUFunction>::clone() const
{
    return new IniCustomEquilibriumFunctional3D<T,Descriptor,RhoUFunction>(*this);
}