NLD_boundaryDynamics3D.hh

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef NLD_BOUNDARY_DYNAMICS_3D_HH
#define NLD_BOUNDARY_DYNAMICS_3D_HH

#include "boundaryCondition/NLD_boundaryDynamics3D.h"
#include "boundaryCondition/boundaryTemplates.h"
#include "core/cell.h"
#include "finiteDifference/fdStencils1D.h"

namespace plb {

namespace flatWall {

template<typename T, template<typename U> class Descriptor>
T computeRhoBar(Cell<T,Descriptor>& cell, int direction, int orientation, T u, T f)
{
    T rhoOnWall = T();
    T rhoNormal = T();
    for (int iPop=0; iPop<Descriptor<T>::q; ++iPop) {
        if (Descriptor<T>::c[iPop][direction]==0) {
            rhoOnWall += cell[iPop];
        }
        else if (Descriptor<T>::c[iPop][direction]==orientation) {
            rhoNormal += cell[iPop];
        }
    }

    T velNormal = (T)orientation * (u-(T)0.5*f);
    T rhoBar = T();
    if (cell.getDynamics().velIsJ()) {
        rhoBar = (T)2*rhoNormal+rhoOnWall-velNormal;
    }
    else {
        rhoBar =((T)2*rhoNormal+rhoOnWall-Descriptor<T>::SkordosFactor()*velNormal)
                    / ((T)1+velNormal);
    }
    return rhoBar;
}

template<typename T, template<typename U> class Descriptor>
void compute_PiNeq (
     Dynamics<T,Descriptor> const& dynamics,
     Cell<T,Descriptor> const& cell, int direction, int orientation,
     T rhoBar, Array<T,Descriptor<T>::d> const& j, T jSqr,
     Array<T,SymmetricTensor<T,Descriptor>::n>& PiNeq )
{
    typedef Descriptor<T> L;

    // Compute off-equilibrium for known particle populations.
    Array<T,Descriptor<T>::q> fNeq;

    for (int iPop=0; iPop<L::q; ++iPop) {
        if (iPop==0) {
            fNeq[iPop] = T();  // fNeq[0] will not be used anyway
        }
        else {
            int cNormal = L::c[iPop][direction];
            if (cNormal==0 || cNormal==orientation) {
                fNeq[iPop] = cell[iPop] - dynamics.computeEquilibrium(iPop, rhoBar, j, jSqr);
            }
        }
    }

    // Compute PiNeq from fNeq, by using "bounce-back of off-equilibrium part" rule.
    int iPi = 0;
    for (int iAlpha=0; iAlpha<L::d; ++iAlpha) {
        for (int iBeta=iAlpha; iBeta<L::d; ++iBeta) {
            PiNeq[iPi] = T();
            for (int iPop=0; iPop<L::q; ++iPop) {
                int cNormal = L::c[iPop][direction];
                if (cNormal==0) {
                    PiNeq[iPi] += L::c[iPop][iAlpha]*L::c[iPop][iBeta] *fNeq[iPop];
                }
                else if (cNormal==orientation) {
                    PiNeq[iPi] += (T)2 * L::c[iPop][iAlpha]*L::c[iPop][iBeta] *fNeq[iPop];
                }
            }
            ++iPi;
        }
    }
}

template<typename T, template<typename U> class Descriptor>
Array<T,3> computeJ (
        Cell<T,Descriptor> const& cell, int direction,
        int orientation, Array<T,3> j, T rhoBar )
{
    T rhoOnWall=T(), rhoNormal=T();
    for (int iPop=0; iPop<Descriptor<T>::q; ++iPop) {
        if (Descriptor<T>::c[iPop][direction]==0) {
            rhoOnWall += cell[iPop];
        }
        else if (Descriptor<T>::c[iPop][direction]==orientation) {
            rhoNormal += cell[iPop];
        }
    }

    j[direction] = (T)orientation*((T)2*rhoNormal+rhoOnWall-rhoBar);
    return j;
}

}  // namespace flatWall


/* *************** Class NonLocalBoundaryDynamics3D ****** */

template<typename T, template<typename U> class Descriptor>
NonLocalBoundaryDynamics3D<T,Descriptor>::NonLocalBoundaryDynamics3D(Dynamics<T,Descriptor>* baseDynamics_)
    : CompositeDynamics<T,Descriptor>(baseDynamics_, false)
{ }

template<typename T, template<typename U> class Descriptor>
bool NonLocalBoundaryDynamics3D<T,Descriptor>::isBoundary() const {
    return true;
}

template<typename T, template<typename U> class Descriptor>
bool NonLocalBoundaryDynamics3D<T,Descriptor>::isNonLocal() const {
    return true;
}




/* *************** Class NLD_VelocityCornerUtil3D ****** */

template<typename T, template<typename U> class Descriptor>
class NLD_VelocityCornerUtil3D {
public:
    NLD_VelocityCornerUtil3D (
            int xNormal_, int yNormal_, int zNormal_,
            plint iX_, plint iY_, plint iZ_,
            BlockLattice3D<T,Descriptor>& lattice_,
            NonLocalBoundaryDynamics3D<T,Descriptor>& dynamics_)
        : xNormal(xNormal_),
          yNormal(yNormal_),
          zNormal(zNormal_),
          iX(iX_), iY(iY_), iZ(iZ_),
          lattice(lattice_),
          dynamics(dynamics_)
    { }
    void computeVelocityGradient(int deriveDirection, int orientation, Array<T,3>& gradient)
    {
        Array<T,3> u0, u1;
        lattice.get(iX,iY,iZ).computeVelocity(u0);
        lattice.get (
            iX+(deriveDirection==0 ? (-orientation):0),
            iY+(deriveDirection==1 ? (-orientation):0),
            iZ+(deriveDirection==2 ? (-orientation):0)  ).computeVelocity(u1);
        for (int iD=0; iD<Descriptor<T>::d; ++iD) {
            gradient[iD] = -orientation * fd::o1_fwd_diff(u0[iD], u1[iD]);
        }
    }
    void computeVelocityGradientX(Array<T,3> const& u0, Array<T,3>& gradient) {
        Array<T,3> u1;
        T rhoBar, dummyTheta;
        NonLocalBoundaryDynamics3D<T,Descriptor>::staticEdgeComputeMacroscopic (
                0, yNormal, zNormal, iX -1*xNormal, iY -0*yNormal, iZ -0*zNormal, lattice, rhoBar, u1, dummyTheta );
        if (!dynamics.velIsJ()) {
            u1 /= Descriptor<T>::fullRho(rhoBar);
        }
        for (int iD=0; iD<Descriptor<T>::d; ++iD) {
            gradient[iD] = -xNormal * fd::o1_fwd_diff(u0[iD], u1[iD]);
        }
    }
    void computeVelocityGradientY(Array<T,3> const& u0, Array<T,3>& gradient) {
        Array<T,3> u1;
        T rhoBar, dummyTheta;
        NonLocalBoundaryDynamics3D<T,Descriptor>::staticEdgeComputeMacroscopic (
                1, zNormal, xNormal, iX -0*xNormal, iY -1*yNormal, iZ -0*zNormal, lattice, rhoBar, u1, dummyTheta );
        if (!dynamics.velIsJ()) {
            u1 /= Descriptor<T>::fullRho(rhoBar);
        }
        for (int iD=0; iD<Descriptor<T>::d; ++iD) {
            gradient[iD] = -yNormal * fd::o1_fwd_diff(u0[iD], u1[iD]);
        }
    }
    void computeVelocityGradientZ(Array<T,3> const& u0, Array<T,3>& gradient) {
        Array<T,3> u1;
        T rhoBar, dummyTheta;
        NonLocalBoundaryDynamics3D<T,Descriptor>::staticEdgeComputeMacroscopic (
                2, xNormal, yNormal, iX -0*xNormal, iY -0*yNormal, iZ -1*zNormal, lattice, rhoBar, u1, dummyTheta );
        if (!dynamics.velIsJ()) {
            u1 /= Descriptor<T>::fullRho(rhoBar);
        }
        for (int iD=0; iD<Descriptor<T>::d; ++iD) {
            gradient[iD] = -zNormal * fd::o1_fwd_diff(u0[iD], u1[iD]);
        }
    }
    void computeMacroscopic(Array<T,3> const& u, T& rhoBar, Array<T,3>& j)
    {
        T rho100, rho010, rho001;
        T dummyTheta;
        Array<T,3> dummyJ;
        NonLocalBoundaryDynamics3D<T,Descriptor>::staticEdgeComputeMacroscopic (
                0, yNormal, zNormal, iX -1*xNormal, iY -0*yNormal, iZ -0*zNormal, lattice, rho100, dummyJ, dummyTheta );
        NonLocalBoundaryDynamics3D<T,Descriptor>::staticEdgeComputeMacroscopic (
                1, zNormal, xNormal, iX -0*xNormal, iY -1*yNormal, iZ -0*zNormal, lattice, rho010, dummyJ, dummyTheta );
        NonLocalBoundaryDynamics3D<T,Descriptor>::staticEdgeComputeMacroscopic (
                2, xNormal, yNormal, iX -0*xNormal, iY -0*yNormal, iZ -1*zNormal, lattice, rho001, dummyJ, dummyTheta );
        rhoBar = (T)1/(T)3 * (rho001 + rho010 + rho100);
        j = u;
        if (dynamics.velIsJ()) {
            j *= Descriptor<T>::fullRho(rhoBar);
        }
    }
    void boundaryCompletion(Array<T,3> const& u)
    {
        typedef SymmetricTensorImpl<T,Descriptor<T>::d> S;
        T rhoBar;
        Array<T,3> j;
        computeMacroscopic(u, rhoBar, j);
        T rho = Descriptor<T>::fullRho(rhoBar);
        Array<T,3> dx_u, dy_u, dz_u;
        computeVelocityGradientX(u, dx_u);
        computeVelocityGradientY(u, dy_u);
        computeVelocityGradientZ(u, dz_u);

        T dx_ux = dx_u[0];
        T dy_ux = dy_u[0];
        T dz_ux = dz_u[0];
        T dx_uy = dx_u[1];
        T dy_uy = dy_u[1];
        T dz_uy = dz_u[1];
        T dx_uz = dx_u[2];
        T dy_uz = dy_u[2];
        T dz_uz = dz_u[2];
        T omega = dynamics.getOmega();
        T sToPi = - rho / Descriptor<T>::invCs2 / omega;
        Array<T,SymmetricTensor<T,Descriptor>::n> pi;
        pi[S::xx] = (T)2 * dx_ux * sToPi;
        pi[S::yy] = (T)2 * dy_uy * sToPi;
        pi[S::zz] = (T)2 * dz_uz * sToPi;
        pi[S::xy] = (dx_uy + dy_ux) * sToPi;
        pi[S::xz] = (dx_uz + dz_ux) * sToPi;
        pi[S::yz] = (dy_uz + dz_uy) * sToPi;

        T jSqr = VectorTemplate<T,Descriptor>::normSqr(j);

        Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
        for (plint iPop = 0; iPop < Descriptor<T>::q; ++iPop) {
            cell[iPop] = dynamics.computeEquilibrium(iPop,rhoBar,j,jSqr) +
                             offEquilibriumTemplates<T,Descriptor>::fromPiToFneq(iPop, pi);
        }
    }
private:
    int xNormal, yNormal, zNormal;
    plint iX,iY,iZ;
    BlockLattice3D<T,Descriptor>& lattice;
    NonLocalBoundaryDynamics3D<T,Descriptor>& dynamics;
};


/* *************** Class NLD_VelocityPlaneUtil3D ****** */

template<typename T, template<typename U> class Descriptor>
class NLD_VelocityPlaneUtil3D {
public:
    NLD_VelocityPlaneUtil3D (
            int direction_, int orientation_,
            plint iX_, plint iY_, plint iZ_,
            BlockLattice3D<T,Descriptor>& lattice_,
            NonLocalBoundaryDynamics3D<T,Descriptor>& dynamics_ )
        : direction(direction_),
          orientation(orientation_),
          iX(iX_), iY(iY_), iZ(iZ_),
          lattice(lattice_),
          dynamics(dynamics_)
    {
        direction1 = (direction+1)%3;
        direction2 = (direction+2)%3;
    }
    void computeMacroscopic(Array<T,3> const& u, T& rhoBar, Array<T,3>& j) const
    {
        Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
        T uNormal = u[direction];
        T fNormal = getExternalForceComponent(cell,direction);
        rhoBar = flatWall::computeRhoBar<T,Descriptor> (
                cell, direction, orientation, uNormal, fNormal );
        j = u;
        if (!dynamics.velIsJ()) {
            j *= Descriptor<T>::fullRho(rhoBar);
        }
    }
    Array<T,3> extrapolateVelocity() const {
        Array<T,3> uTot; uTot.resetToZero();
        static T weight0 = 1./9.;
        static T weight1 = 1./36.;
        static T weight2 = 4./9.;
        static T weights[] = { weight1, weight0, weight1,
                               weight0, weight2, weight0,
                               weight1, weight0, weight1 };
        plint i=0;
        plint numCells=0;
        T sumWeights = T();
        for (plint i1=-1; i1<=+1; ++i1) {
            for (plint i2=-1; i2<=+1; ++i2) {
                Array<plint,3> neighb(iX,iY,iZ);
                neighb[direction] -= orientation;
                neighb[direction1] += i1;
                neighb[direction2] += i2;
                Cell<T,Descriptor>& neighborCell = lattice.get(neighb[0],neighb[1],neighb[2]);
                if (!neighborCell.getDynamics().isBoundary())  {
                    Array<T,3> u; neighborCell.computeVelocity(u);
                    uTot += weights[i]*u;
                    ++numCells;
                    sumWeights += weights[i];
                }
                ++i;
            }
        }
        if (numCells != 9) {
            uTot /= sumWeights;
        }
        return uTot;
    }
private:
    int direction, orientation;
    int direction1, direction2;
    int iX,iY,iZ;
    BlockLattice3D<T,Descriptor>& lattice;
    NonLocalBoundaryDynamics3D<T,Descriptor>& dynamics;
};


/* *************** Class NLD_VelocityEdgeUtil3D ****** */

template<typename T, template<typename U> class Descriptor>
class NLD_VelocityEdgeUtil3D {
public:
    NLD_VelocityEdgeUtil3D (
            int plane_, int normal1_, int normal2_,
            plint iX_, plint iY_, plint iZ_,
            BlockLattice3D<T,Descriptor>& lattice_,
            NonLocalBoundaryDynamics3D<T,Descriptor>& dynamics_ )
        : plane(plane_),
          normal1(normal1_),
          normal2(normal2_),
          iX(iX_), iY(iY_), iZ(iZ_),
          lattice(lattice_),
          dynamics(dynamics_)
    {
        direction1 = (plane+1)%3;
        direction2 = (plane+2)%3;
    }
    T getNeighborRhoBar(plint step1, plint step2)
    {
        PLB_ASSERT(step1==0 || step2==0);
        Array<T,3> dummyJ;
        T dummyTheta;
        T rhoBar;
        Array<int,3> coords (iX, iY, iZ);
        coords[direction1] += -normal1*step1;
        coords[direction2] += -normal2*step2;
        if (step1==0) {
            NonLocalBoundaryDynamics3D<T,Descriptor>::staticPlaneComputeMacroscopic (
                    direction1, normal1, coords[0], coords[1], coords[2], lattice, rhoBar, dummyJ, dummyTheta );
            PLB_ASSERT(step2 != 0);
        }
        else {  // step2==0.
            NonLocalBoundaryDynamics3D<T,Descriptor>::staticPlaneComputeMacroscopic (
                    direction2, normal2, coords[0], coords[1], coords[2], lattice, rhoBar, dummyJ, dummyTheta );
        }
        return rhoBar;
    }
    void computeMacroscopic(Array<T,3> const& u, T& rhoBar, Array<T,3>& j)
    {
        T rho10 = getNeighborRhoBar(1,0);
        T rho01 = getNeighborRhoBar(0,1);
        T rho20 = getNeighborRhoBar(2,0);
        T rho02 = getNeighborRhoBar(0,2);
        rhoBar = (T)2/(T)3*(rho01+rho10)-(T)1/(T)6*(rho02+rho20);
        j = u;
        if (dynamics.velIsJ()) {
            j *= Descriptor<T>::fullRho(rhoBar);
        }
    }
    void computeEdgeGradient(Array<T,3>& gradient) {
        Array<T,3> u_p1, u_m1;
        Array<plint,3> posPlus(iX,iY,iZ), posMinus(iX,iY,iZ);
        ++posPlus[plane];
        --posMinus[plane];

        NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dynPlus = 
            dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                    &lattice.get(posPlus[0], posPlus[1], posPlus[2]).getDynamics() );
        NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dynMinus = 
            dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                    &lattice.get(posMinus[0], posMinus[1], posMinus[2]).getDynamics() );
        if (dynPlus && dynMinus) {
            u_p1 = dynPlus->get_u();
            u_m1 = dynMinus->get_u();
            for (int iD=0; iD<Descriptor<T>::d; ++iD) {
                gradient[iD] = fd::ctl_diff(u_p1[iD],u_m1[iD]);
            }
        }
        else {
            gradient.resetToZero();
        }
    }
    void computeGradientIntoPlane(Array<T,3> const& u0, int deriveDirection, int deriveOrientation,
                                  int normalDirection, int normalOrientation, Array<T,3>& gradient)
    {
        Array<T,3> u1;
        Array<plint,3> pos1(iX,iY,iZ);
        pos1[deriveDirection] -= deriveOrientation;
        
        T rhoBar, dummyTheta;
        NonLocalBoundaryDynamics3D<T,Descriptor>::staticPlaneComputeMacroscopic (
                normalDirection, normalOrientation, pos1[0], pos1[1], pos1[2], lattice, rhoBar, u1, dummyTheta );
        if (!dynamics.velIsJ()) {
            u1 /= Descriptor<T>::fullRho(rhoBar);
        }
        for (int iD=0; iD<Descriptor<T>::d; ++iD) {
            gradient[iD] = -deriveOrientation * fd::o1_fwd_diff(u0[iD], u1[iD]);
        }
    }
    void boundaryCompletion(Array<T,3> const& u)
    {
        typedef SymmetricTensorImpl<T,Descriptor<T>::d> S;
        T rhoBar;
        Array<T,3> j;
        computeMacroscopic(u, rhoBar, j);
        T rho = Descriptor<T>::fullRho(rhoBar);
        std::vector<Array<T,3> > dA_uB_(3);
        computeEdgeGradient(dA_uB_[0]);
        computeGradientIntoPlane(u, direction1, normal1, direction2, normal2, dA_uB_[1]);
        computeGradientIntoPlane(u, direction2, normal2, direction1, normal1, dA_uB_[2]);

        std::vector<Array<T,3> > dA_uB(3);
        for (int iBeta=0; iBeta<3; ++iBeta) {
            dA_uB[plane][iBeta]      = dA_uB_[0][iBeta];
            dA_uB[direction1][iBeta] = dA_uB_[1][iBeta];
            dA_uB[direction2][iBeta] = dA_uB_[2][iBeta];
        }

        T sToPi = - rho / Descriptor<T>::invCs2 / dynamics.getOmega();
        Array<T,SymmetricTensor<T,Descriptor>::n> pi;
        pi[S::xx] = (T)2 * dA_uB[0][0] * sToPi;
        pi[S::yy] = (T)2 * dA_uB[1][1] * sToPi;
        pi[S::zz] = (T)2 * dA_uB[2][2] * sToPi;
        pi[S::xy] = (dA_uB[0][1]+dA_uB[1][0]) * sToPi;
        pi[S::xz] = (dA_uB[0][2]+dA_uB[2][0]) * sToPi;
        pi[S::yz] = (dA_uB[1][2]+dA_uB[2][1]) * sToPi;

        T jSqr = VectorTemplate<T,Descriptor>::normSqr(j);

        Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
        for (plint iPop = 0; iPop < Descriptor<T>::q; ++iPop) {
            cell[iPop] = dynamics.computeEquilibrium(iPop,rhoBar,j,jSqr) +
                             offEquilibriumTemplates<T,Descriptor>::fromPiToFneq(iPop, pi);
        }
    }
private:
    int plane, normal1, normal2, direction1, direction2;
    plint iX,iY,iZ;
    BlockLattice3D<T,Descriptor>& lattice;
    NonLocalBoundaryDynamics3D<T,Descriptor>& dynamics;
};


/* *************** Class NLD_VelocityBoundaryDynamics3D ****** */

template<typename T, template<typename U> class Descriptor>
int NLD_VelocityBoundaryDynamics3D<T,Descriptor>::id =
    meta::registerGeneralDynamics<T,Descriptor, NLD_VelocityBoundaryDynamics3D<T,Descriptor> >
            ( std::string("NLD_Velocity") );

template<typename T, template<typename U> class Descriptor>
NLD_VelocityBoundaryDynamics3D<T,Descriptor>::NLD_VelocityBoundaryDynamics3D(Dynamics<T,Descriptor>* baseDynamics_)
    : NonLocalBoundaryDynamics3D<T,Descriptor>(baseDynamics_)
{ }

template<typename T, template<typename U> class Descriptor>
NLD_VelocityBoundaryDynamics3D<T,Descriptor>::NLD_VelocityBoundaryDynamics3D(HierarchicUnserializer& unserializer)
    : NonLocalBoundaryDynamics3D<T,Descriptor>(0)
{
    this->unserialize(unserializer);
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::defineVelocity (
        Cell<T,Descriptor>& cell, Array<T,Descriptor<T>::d> const& velocity_ )
{
    u = velocity_;
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::computeVelocity (
        Cell<T,Descriptor> const& cell, Array<T,Descriptor<T>::d>& velocity_ ) const
{
    velocity_ = u;
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::serialize(HierarchicSerializer& serializer) const {
    serializer.addValues<T,3>(u);
    NonLocalBoundaryDynamics3D<T,Descriptor>::serialize(serializer);
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::unserialize(HierarchicUnserializer& unserializer) {
    unserializer.readValues<T,3>(u);
    NonLocalBoundaryDynamics3D<T,Descriptor>::unserialize(unserializer);
}

template<typename T, template<typename U> class Descriptor>
int NLD_VelocityBoundaryDynamics3D<T,Descriptor>::getId() const
{
    return id;
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::prepareCollision(Cell<T,Descriptor>& cell)
{ }

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


template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::planeComputeMacroscopic (
            int direction, int orientation, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice,
            T& rhoBar, Array<T,3>& j, T& thetaBar )
{
    NLD_VelocityPlaneUtil3D<T,Descriptor> (
            direction, orientation, iX,iY,iZ, lattice, *this ).computeMacroscopic(this->u, rhoBar, j);
    thetaBar = T();
}


template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::edgeComputeMacroscopic (
            int plane, int normal1, int normal2, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice,
            T& rhoBar, Array<T,3>& j, T& thetaBar )
{
    NLD_VelocityEdgeUtil3D<T,Descriptor> edgeUtil (
            plane, normal1, normal2, iX, iY, iZ, lattice, *this );
    edgeUtil.computeMacroscopic(this->u, rhoBar, j);
    thetaBar = T();
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::cornerComputeMacroscopic (
            int xNormal, int yNormal, int zNormal, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice,
            T& rhoBar, Array<T,3>& j, T& thetaBar )
{
    NLD_VelocityCornerUtil3D<T,Descriptor> cornerUtil(xNormal, yNormal, zNormal, iX,iY,iZ, lattice, *this);
    cornerUtil.computeMacroscopic(this->u, rhoBar, j);
    thetaBar = T();
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::planeBoundaryCompletion (
            int direction, int orientation, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice )
{

    Array<T,3> j;
    T rhoBar;
    NLD_VelocityPlaneUtil3D<T,Descriptor> (
            direction, orientation, iX,iY,iZ, lattice, *this ).computeMacroscopic(this->u, rhoBar, j);
    T jSqr = VectorTemplate<T,Descriptor>::normSqr(j);

    Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
    Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
    flatWall::compute_PiNeq<T,Descriptor> (
            this->getBaseDynamics(), cell, direction, orientation, rhoBar, j, jSqr, PiNeq );

    this->getBaseDynamics().regularize(cell, rhoBar, j, jSqr, PiNeq);
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::edgeBoundaryCompletion (
            int plane, int normal1, int normal2, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice )
{
    int direction1 = (plane+1)%3;
    int direction2 = (plane+2)%3;
    Array<plint,3> nbCell1(iX,iY,iZ), nbCell2(iX,iY,iZ);
    nbCell1[direction1] -= normal1;
    nbCell2[direction2] -= normal2;

    // 1. Detect case where a direction is periodic.
    bool isWall1 =
        lattice.get(nbCell1[0],nbCell1[1],nbCell1[2]).getDynamics().isBoundary();
    bool isWall2 =
        lattice.get(nbCell2[0],nbCell2[1],nbCell2[2]).getDynamics().isBoundary();

    // If both directions are periodic, there is no need for a boundary condition.
    if (!isWall1 && !isWall2) return;
    // If exactly one direction is periodic, then the edge behaves like a plane.
    if (!isWall1) {
        planeBoundaryCompletion(direction1, normal1, iX,iY,iZ, lattice);
        return;
    }
    else if (!isWall2) {
        planeBoundaryCompletion(direction2, normal2, iX,iY,iZ, lattice);
        return;
    }

    // 2. Implement general edge algorithm.
    NLD_VelocityEdgeUtil3D<T,Descriptor> edgeUtil (
            plane, normal1, normal2, iX, iY, iZ, lattice, *this );
    edgeUtil.boundaryCompletion(this->u);
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityBoundaryDynamics3D<T,Descriptor>::cornerBoundaryCompletion (
            int xNormal, int yNormal, int zNormal, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice )
{
    Array<plint,3> middleNbCell(iX,iY,iZ), nbCellX(iX,iY,iZ), nbCellY(iX,iY,iZ), nbCellZ(iX,iY,iZ);
    middleNbCell -= Array<plint,3>(xNormal,yNormal,zNormal);
    nbCellX[0] -= xNormal;
    nbCellY[1] -= yNormal;
    nbCellZ[2] -= zNormal;

    // 1. To start width, detect periodic directions where the corner actually
    //    becomes edge, plane, or non-boundary.

    bool xPeriodic = ! lattice.get(iX,iY-yNormal,iZ-zNormal).getDynamics().isBoundary();
    bool yPeriodic = ! lattice.get(iX-xNormal,iY,iZ-zNormal).getDynamics().isBoundary();
    bool zPeriodic = ! lattice.get(iX-xNormal,iY-yNormal,iZ).getDynamics().isBoundary();

    // Note that we only need to consider cases where corners become edges,
    //   because further degenerate cases (two or three periodic directions)
    //   are treated in the edge algorithm.
    if (xPeriodic) {
        edgeBoundaryCompletion(0, yNormal, zNormal, iX,iY,iZ, lattice);
        return;
    }
    if (yPeriodic) {
        edgeBoundaryCompletion(1, zNormal, xNormal, iX,iY,iZ, lattice);
        return;
    }
    if (zPeriodic) {
        edgeBoundaryCompletion(2, xNormal, yNormal, iX,iY,iZ, lattice);
        return;
    }

    // 2. Implement generic case.
    NLD_VelocityCornerUtil3D<T,Descriptor> cornerUtil(xNormal, yNormal, zNormal, iX,iY,iZ, lattice, *this);
    cornerUtil.boundaryCompletion(this->u);
}


/* *************** Class NLD_VelocityNeumannBoundaryDynamics3D ****** */

template<typename T, template<typename U> class Descriptor>
int NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::id =
    meta::registerGeneralDynamics<T,Descriptor, NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor> >
            ( std::string("NLD_VelocityNeumann") );

template<typename T, template<typename U> class Descriptor>
NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::NLD_VelocityNeumannBoundaryDynamics3D (
        Dynamics<T,Descriptor>* baseDynamics_ )
    : NonLocalBoundaryDynamics3D<T,Descriptor>(baseDynamics_)
{ }

template<typename T, template<typename U> class Descriptor>
NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::NLD_VelocityNeumannBoundaryDynamics3D (
        HierarchicUnserializer& unserializer )
    : NonLocalBoundaryDynamics3D<T,Descriptor>(0)
{
    this->unserialize(unserializer);
}

template<typename T, template<typename U> class Descriptor>
int NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::getId() const
{
    return id;
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::prepareCollision(Cell<T,Descriptor>& cell)
{ }

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


template<typename T, template<typename U> class Descriptor>
void NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::planeComputeMacroscopic (
            int direction, int orientation, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice,
            T& rhoBar, Array<T,3>& j, T& thetaBar )
{
    NLD_VelocityPlaneUtil3D<T,Descriptor> planeUtil (
            direction, orientation, iX,iY,iZ, lattice, *this );
    Array<T,3> u(planeUtil.extrapolateVelocity());
    planeUtil.computeMacroscopic(u, rhoBar, j);
    thetaBar = T();
}


template<typename T, template<typename U> class Descriptor>
void NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::edgeComputeMacroscopic (
            int plane, int normal1, int normal2, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice,
            T& rhoBar, Array<T,3>& j, T& thetaBar )
{
    int direction1 = (plane+1)%3;
    int direction2 = (plane+2)%3;

    Array<plint,3> nbCell1(iX,iY,iZ), nbCell2(iX,iY,iZ);
    nbCell1[direction1] -= normal1;
    nbCell2[direction2] -= normal2;
    Cell<T,Descriptor>& neighbor1 = lattice.get(nbCell1[0], nbCell1[1], nbCell1[2]);
    Cell<T,Descriptor>& neighbor2 = lattice.get(nbCell2[0], nbCell2[1], nbCell2[2]);
    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn1 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor1.getDynamics() );
    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn2 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor2.getDynamics() );
    if (!(dyn1 || dyn2)) {
        Array<plint,3> middleNbCell(iX,iY,iZ);
        middleNbCell[direction1] -= normal1;
        middleNbCell[direction2] -= normal2;
        Cell<T,Descriptor>& middleNeighbor = lattice.get(middleNbCell[0], middleNbCell[1], middleNbCell[2]);
        Array<T,3> u;
        middleNeighbor.computeVelocity(u);
        NLD_VelocityEdgeUtil3D<T,Descriptor> edgeUtil (
                plane, normal1, normal2, iX, iY, iZ, lattice, *this );
        edgeUtil.computeMacroscopic(u, rhoBar, j);
        thetaBar = T();
    }
    else {
        T rhoBar1, rhoBar2;
        Array<T,3> j1, j2;
        if (dyn1) {
            dyn1->edgeComputeMacroscopic (
                    plane, normal1, normal2, iX,iY,iZ, lattice, rhoBar1, j1, thetaBar );
            rhoBar = rhoBar1;
            j = j1;
        }
        if (dyn2) {
            dyn2->edgeComputeMacroscopic (
                    plane, normal1, normal2, iX,iY,iZ, lattice, rhoBar2, j2, thetaBar );
            rhoBar = rhoBar2;
            j = j2;
        }
        if (dyn1 && dyn2) {
            j = (j1+j2)/(T)2;
            rhoBar = (rhoBar1+rhoBar2)/(T)2;
        }
    }
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::cornerComputeMacroscopic (
            int xNormal, int yNormal, int zNormal, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice,
            T& rhoBar, Array<T,3>& j, T& thetaBar )
{
    Array<plint,3> nbCell1(iX-xNormal,iY,iZ), nbCell2(iX,iY-yNormal,iZ), nbCell3(iX,iY,iZ-zNormal);
    Cell<T,Descriptor>& neighbor1 = lattice.get(nbCell1[0], nbCell1[1], nbCell1[2]);
    Cell<T,Descriptor>& neighbor2 = lattice.get(nbCell2[0], nbCell2[1], nbCell2[2]);
    Cell<T,Descriptor>& neighbor3 = lattice.get(nbCell3[0], nbCell3[1], nbCell3[2]);

    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn1 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor1.getDynamics() );
    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn2 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor2.getDynamics() );
    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn3 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor3.getDynamics() );

    if (!(dyn1 || dyn2 || dyn3)) {
        Array<plint,3> nbCell(iX,iY,iZ);
        nbCell -= Array<plint,3>(xNormal,yNormal,zNormal);
        Cell<T,Descriptor>& middleNeighbor = lattice.get(nbCell[0], nbCell[1], nbCell[2]);
        Array<T,3> u; middleNeighbor.computeVelocity(u);

        NLD_VelocityCornerUtil3D<T,Descriptor> cornerUtil(xNormal, yNormal, zNormal, iX,iY,iZ, lattice, *this);
        cornerUtil.computeMacroscopic(u, rhoBar, j);
        thetaBar = T();
    }
    else {
        rhoBar = T();
        j.resetToZero();
        plint numNeighbors=0;
        if (dyn1) {
            T rhoBar_;
            Array<T,3> j_;
            dyn1->cornerComputeMacroscopic (
                    xNormal,yNormal,zNormal, iX,iY,iZ, lattice, rhoBar_, j_, thetaBar );
            rhoBar += rhoBar_;
            j += j_;
            ++numNeighbors;
        }
        if (dyn2) {
            T rhoBar_;
            Array<T,3> j_;
            dyn2->cornerComputeMacroscopic (
                    xNormal,yNormal,zNormal, iX,iY,iZ, lattice, rhoBar_, j_, thetaBar );
            rhoBar += rhoBar_;
            j += j_;
            ++numNeighbors;
        }
        if (dyn3) {
            T rhoBar_;
            Array<T,3> j_;
            dyn3->cornerComputeMacroscopic (
                    xNormal,yNormal,zNormal, iX,iY,iZ, lattice, rhoBar_, j_, thetaBar );
            rhoBar += rhoBar_;
            j += j_;
            ++numNeighbors;
        }
        if (numNeighbors>1) {
            rhoBar /= (T)numNeighbors;
            j /= (T)numNeighbors;
        }
    }
}


template<typename T, template<typename U> class Descriptor>
void NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::planeBoundaryCompletion (
            int direction, int orientation, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice )
{
    NLD_VelocityPlaneUtil3D<T,Descriptor> planeUtil (
            direction, orientation, iX,iY,iZ, lattice, *this );
    Array<T,3> u(planeUtil.extrapolateVelocity());
    T rhoBar;
    Array<T,3> j;
    planeUtil.computeMacroscopic(u, rhoBar, j);
    T jSqr = VectorTemplate<T,Descriptor>::normSqr(j);

    Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
    Array<T,SymmetricTensor<T,Descriptor>::n> PiNeq;
    flatWall::compute_PiNeq<T,Descriptor> (
            this->getBaseDynamics(), cell, direction, orientation, rhoBar, j, jSqr, PiNeq );

    this->getBaseDynamics().regularize(cell, rhoBar, j, jSqr, PiNeq);
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::edgeBoundaryCompletion (
            int plane, int normal1, int normal2, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice )
{
    int direction1 = (plane+1)%3;
    int direction2 = (plane+2)%3;

    Array<plint,3> nbCell1(iX,iY,iZ), nbCell2(iX,iY,iZ), middleNbCell(iX,iY,iZ);
    nbCell1[direction1] -= normal1;
    nbCell2[direction2] -= normal2;
    middleNbCell[direction1] -= normal1;
    middleNbCell[direction2] -= normal2;

    // 1. Detect case where a direction is periodic.
    bool isWall1 =
        lattice.get(nbCell1[0],nbCell1[1],nbCell1[2]).getDynamics().isBoundary();
    bool isWall2 =
        lattice.get(nbCell2[0],nbCell2[1],nbCell2[2]).getDynamics().isBoundary();

    // If both directions are periodic, there is no need for a boundary condition.
    if (!isWall1 && !isWall2) return;
    // If exactly one direction is periodic, then the edge behaves like a plane.
    if (!isWall1) {
        planeBoundaryCompletion(direction1, normal1, iX,iY,iZ, lattice);
        return;
    }
    else if (!isWall2) {
        planeBoundaryCompletion(direction2, normal2, iX,iY,iZ, lattice);
        return;
    }

    Cell<T,Descriptor>& neighbor1 = lattice.get(nbCell1[0], nbCell1[1], nbCell1[2]);
    Cell<T,Descriptor>& neighbor2 = lattice.get(nbCell2[0], nbCell2[1], nbCell2[2]);
    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn1 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor1.getDynamics() );
    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn2 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor2.getDynamics() );

    if (!(dyn1 || dyn2)) {
        Cell<T,Descriptor>& middleNeighbor = lattice.get(middleNbCell[0], middleNbCell[1], middleNbCell[2]);
        Array<T,3> u; middleNeighbor.computeVelocity(u);
        NLD_VelocityEdgeUtil3D<T,Descriptor> edgeUtil (
                plane, normal1, normal2, iX, iY, iZ, lattice, *this );
        edgeUtil.boundaryCompletion(u);
    }
    else {
        if (dyn1) {
            dyn1->edgeBoundaryCompletion (
                        plane, normal1, normal2, iX, iY, iZ, lattice );
        }
        if (dyn2) {
            Cell<T,Descriptor> saveCell(lattice.get(iX,iY,iZ));
            dyn2->edgeBoundaryCompletion (
                        plane, normal1, normal2, iX, iY, iZ, lattice );
            if (dyn1) {
                Cell<T,Descriptor>& cell(lattice.get(iX,iY,iZ));
                for (plint iPop=0; iPop<Descriptor<T>::q; ++iPop) {
                    cell[iPop] = (cell[iPop]+saveCell[iPop])/2.;
                }
            }
        }
    }
}

template<typename T, template<typename U> class Descriptor>
void NLD_VelocityNeumannBoundaryDynamics3D<T,Descriptor>::cornerBoundaryCompletion (
            int xNormal, int yNormal, int zNormal, plint iX, plint iY, plint iZ,
            BlockLattice3D<T,Descriptor>& lattice )
{
    Cell<T,Descriptor>& cell = lattice.get(iX,iY,iZ);
    Array<plint,3> nbCell(iX,iY,iZ);
    nbCell -= Array<plint,3>(xNormal,yNormal,zNormal);

    // 1. To start width, detect periodic directions where the corner actually
    //    becomes edge, plane, or non-boundary.

    bool xPeriodic = ! lattice.get(iX,iY-yNormal,iZ-zNormal).getDynamics().isBoundary();
    bool yPeriodic = ! lattice.get(iX-xNormal,iY,iZ-zNormal).getDynamics().isBoundary();
    bool zPeriodic = ! lattice.get(iX-xNormal,iY-yNormal,iZ).getDynamics().isBoundary();

    // Note that we only need to consider cases where corners become edges,
    //   because further degenerate cases (two or three periodic directions)
    //   are treated in the edge algorithm.
    if (xPeriodic) {
        edgeBoundaryCompletion(0, yNormal, zNormal, iX,iY,iZ, lattice);
        return;
    }
    if (yPeriodic) {
        edgeBoundaryCompletion(1, zNormal, xNormal, iX,iY,iZ, lattice);
        return;
    }
    if (zPeriodic) {
        edgeBoundaryCompletion(2, xNormal, yNormal, iX,iY,iZ, lattice);
        return;
    }


    Array<plint,3> nbCell1(iX-xNormal,iY,iZ), nbCell2(iX,iY-yNormal,iZ), nbCell3(iX,iY,iZ-zNormal);
    Cell<T,Descriptor>& neighbor1 = lattice.get(nbCell1[0], nbCell1[1], nbCell1[2]);
    Cell<T,Descriptor>& neighbor2 = lattice.get(nbCell2[0], nbCell2[1], nbCell2[2]);
    Cell<T,Descriptor>& neighbor3 = lattice.get(nbCell3[0], nbCell3[1], nbCell3[2]);

    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn1 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor1.getDynamics() );
    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn2 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor2.getDynamics() );
    NLD_VelocityBoundaryDynamics3D<T,Descriptor>* dyn3 = 
        dynamic_cast<NLD_VelocityBoundaryDynamics3D<T,Descriptor>*> (
                &neighbor3.getDynamics() );

    if (!(dyn1 || dyn2 || dyn3)) {
        Cell<T,Descriptor>& middleNeighbor = lattice.get(nbCell[0], nbCell[1], nbCell[2]);
        Array<T,3> u; middleNeighbor.computeVelocity(u);
        NLD_VelocityCornerUtil3D<T,Descriptor> cornerUtil(xNormal, yNormal, zNormal, iX,iY,iZ, lattice, *this);
        cornerUtil.boundaryCompletion(u);
    }
    else {
        Cell<T,Descriptor> tmpCell(lattice.get(iX,iY,iZ));
        for (plint iPop=0; iPop<Descriptor<T>::q; ++iPop) {
            tmpCell[iPop] = T();
        }
        plint numNeighbors=0;
        if (dyn1) {
            dyn1->cornerBoundaryCompletion (
                    xNormal,yNormal,zNormal, iX,iY,iZ, lattice);
            for (plint iPop=0; iPop<Descriptor<T>::q; ++iPop) {
                tmpCell[iPop] += cell[iPop];
            }
            ++numNeighbors;
        }
        if (dyn2) {
            dyn2->cornerBoundaryCompletion (
                    xNormal,yNormal,zNormal, iX,iY,iZ, lattice);
            for (plint iPop=0; iPop<Descriptor<T>::q; ++iPop) {
                tmpCell[iPop] += cell[iPop];
            }
            ++numNeighbors;
        }
        if (dyn3) {
            dyn3->cornerBoundaryCompletion (
                    xNormal,yNormal,zNormal, iX,iY,iZ, lattice);
            for (plint iPop=0; iPop<Descriptor<T>::q; ++iPop) {
                tmpCell[iPop] += cell[iPop];
            }
            ++numNeighbors;
        }
        if (numNeighbors>1) {
            for (plint iPop=0; iPop<Descriptor<T>::q; ++iPop) {
                cell[iPop] = tmpCell[iPop] / (T)numNeighbors;
            }
        }
    }
}


}  // namespace plb

// namespace plb

#endif  // NLD_BOUNDARY_DYNAMICS_3D_HH