Re: Slip-free and pressure corner boundary condition

template<typename T, template<typename U> class Lattice>
SlipBoundaryProcessor3D<T,Lattice>::
SlipBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_, int discreteNormalX, int discreteNormalY, int discreteNormalZ)
: x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_)
{
OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1);
reflectionPop[0] = 0;
for (int iPop = 1; iPop < Lattice<T>::q; iPop++) {
reflectionPop[iPop] = 0;
// iPop are the directions which pointing into the fluid, discreteNormal is pointing outwarts
if (Lattice<T>::c[iPop][0]*discreteNormalX + Lattice<T>::c[iPop][1]*discreteNormalY + Lattice<T>::c[iPop][2]*discreteNormalZ < 0) {
// std::cout << “—–” <<s td::endl;
int mirrorDirection0;
int mirrorDirection1;
int mirrorDirection2;
int mult = 2 / (discreteNormalX*discreteNormalX + discreteNormalY*discreteNormalY + discreteNormalZ*discreteNormalZ);

mirrorDirection0 = (Lattice<T>::c[iPop][0] – mult*(Lattice<T>::c[iPop][0]*discreteNormalX + Lattice<T>::c[iPop][1]*discreteNormalY + Lattice<T>::c[iPop][2]*discreteNormalZ)*discreteNormalX);
mirrorDirection1 = (Lattice<T>::c[iPop][1] – mult*(Lattice<T>::c[iPop][0]*discreteNormalX + Lattice<T>::c[iPop][1]*discreteNormalY + Lattice<T>::c[iPop][2]*discreteNormalZ)*discreteNormalY);
mirrorDirection2 = (Lattice<T>::c[iPop][2] – mult*(Lattice<T>::c[iPop][0]*discreteNormalX + Lattice<T>::c[iPop][1]*discreteNormalY + Lattice<T>::c[iPop][2]*discreteNormalZ)*discreteNormalZ);

// bounce back for the case discreteNormalX = discreteNormalY = discreteNormalZ = 1, that is mult=0
if (mult == 0) {
mirrorDirection0 = -Lattice<T>::c[iPop][0];
mirrorDirection1 = -Lattice<T>::c[iPop][1];
mirrorDirection2 = -Lattice<T>::c[iPop][2];
}

// computes mirror jPop
for (reflectionPop[iPop] = 1; reflectionPop[iPop] < Lattice<T>::q ; reflectionPop[iPop]++) {
if (Lattice<T>::c[reflectionPop[iPop]][0]==mirrorDirection0 && Lattice<T>::c[reflectionPop[iPop]][1]==mirrorDirection1 && Lattice<T>::c[reflectionPop[iPop]][2]==mirrorDirection2) {
break;
}
}
//std::cout <<iPop << ” to “<< jPop <<” for discreteNormal= “<< discreteNormalX << “/”<<discreteNormalY <<std::endl;
}
}
}