amroc/lbm/LBMD3Q19.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2012 Oak Ridge National Laboratory
00004 // Ralf Deiterding, deiterdingr@ornl.gov
00005 
00006 #ifndef LBM_D3Q19_H
00007 #define LBM_D3Q19_H
00008 
00016 #include "LBMBase.h"
00017 #include <cfloat>
00018 
00019 #define LB_FACTOR 1.0e5
00020 
00050 template <class DataType>
00051 class LBMD3Q19 : public LBMBase<Vector<DataType,19>, Vector<DataType,4>, 3> {
00052   typedef LBMBase<Vector<DataType,19>, Vector<DataType,4>, 3> base;
00053 public:
00054   typedef typename base::vec_grid_data_type vec_grid_data_type;
00055   typedef typename base::grid_data_type grid_data_type;
00056   typedef typename base::MicroType MicroType;
00057   typedef typename base::MacroType MacroType;
00058   typedef GridData<MacroType,3> macro_grid_data_type;
00059   typedef typename base::SideName SideName;
00060   typedef typename base::point_type point_type;
00061   typedef Vector<DataType,6> TensorType;
00062 
00063   enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro };
00064   enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall,
00065                       CharacteristicOutlet, CharacteristicInlet, NoSlipWallEntranceExit };
00066   enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall, GFMWallLaw };
00067   enum TurbulenceModel { laminar, LES_Smagorinsky, LES_dynamic };
00068 
00069   LBMD3Q19() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.)  {
00070     cs2  = DataType(1.)/DataType(3.);
00071     cs22 = DataType(2.)*cs2;
00072     cssq = DataType(2.)/DataType(9.);
00073     method[0] = 2; method[1] = 0; method[2] = 0; method[3] = 0;
00074     turbulence = laminar;
00075     mdx[0]=mdx[3]=mdx[4]=mdx[5]=mdx[6]=mdx[11]=mdx[12]=mdx[13]=mdx[14]=0;
00076     mdx[1]=mdx[7]=mdx[9]=mdx[15]=mdx[17]=1;
00077     mdx[2]=mdx[8]=mdx[10]=mdx[16]=mdx[18]=-1;
00078 
00079     mdy[0]=mdy[1]=mdy[2]=mdy[5]=mdy[6]=mdy[15]=mdy[16]=mdy[17]=mdy[18]=0;
00080     mdy[3]=mdy[7]=mdy[10]=mdy[11]=mdy[13]=1;
00081     mdy[4]=mdy[8]=mdy[9]=mdy[12]=mdy[14]=-1;
00082 
00083     mdz[0]=mdz[1]=mdz[2]=mdz[3]=mdz[4]=mdz[7]=mdz[8]=mdz[9]=mdz[10]=0;
00084     mdz[5]=mdz[11]=mdz[14]=mdz[15]=mdz[18]=1;
00085     mdz[6]=mdz[12]=mdz[13]=mdz[16]=mdz[17]=-1;
00086 
00087   }
00088 
00089   virtual ~LBMD3Q19() {}
00090 
00091   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00092     base::LocCtrl = Ctrl.getSubDevice(prefix+"D3Q19");
00093     RegisterAt(base::LocCtrl,"Method(1)",method[0]);
00094     RegisterAt(base::LocCtrl,"Method(2)",method[1]); 
00095     RegisterAt(base::LocCtrl,"Method(3)",method[2]); 
00096     RegisterAt(base::LocCtrl,"Method(4)",method[3]); 
00097     RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00098     RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00099     RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00100     RegisterAt(base::LocCtrl,"Gas_Cs",csp);
00101     RegisterAt(base::LocCtrl,"Gas_nu",nup);
00102     RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00103     RegisterAt(base::LocCtrl,"Gas_W",Wp);
00104     RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00105     RegisterAt(base::LocCtrl,"SpeedUp",S0);
00106   }
00107 
00108   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00109     base::SetupData(gh,ghosts);
00110     if (rhop>0.) R0 = rhop;
00111     if (csp>0.) {
00112       cs2p = csp*csp;
00113       SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00114     }
00115     std::cout << "D3Q19: Gas_rho=" << rhop << "   Gas_Cs=" << csp
00116               << "   Gas_nu=" << nup << std::endl;
00117     std::cout << "D3Q19: L0=" << L0() << "   T0=" << T0() << "   U0=" << U0
00118               << "   R0=" << R0 << std::endl;
00119     if ( TurbulenceType() ==  laminar ) {
00120       std::cout << "LBM Setup() Laminar" << std::endl;
00121     }
00122     if ( TurbulenceType() ==  LES_Smagorinsky ) {
00123       std::cout << "LBM Setup() LES_Smagorinsky : Smagorinsky Constant = " << SmagorinskyConstant()
00124                 << std::endl;
00125 
00126     }
00127     if ( TurbulenceType() ==  LES_dynamic ) {
00128       std::cout << "LBM Setup() LES_dynamic " << std::endl;
00129     }
00130     WriteInit();
00131   }
00132 
00133   virtual void WriteInit() const {
00134     int me = MY_PROC; 
00135     if (me == VizServer) {
00136       std::ofstream ofs("chem.dat", std::ios::out);
00137       ofs << "RU         " << Rp << std::endl;
00138       ofs << "PA         " << BasePressure() << std::endl;
00139       ofs << "Sp(01)     Vicosity" << std::endl;
00140       ofs << "W(01)      " << Wp << std::endl;
00141       ofs << "Sp(02)     |Velocity|" << std::endl;
00142       ofs << "W(02)      " << 1. << std::endl;
00143       ofs << "Sp(03)     Vorticity-x" << std::endl;
00144       ofs << "W(03)      " << 1. << std::endl;
00145       ofs << "Sp(04)     Vorticity-y" << std::endl;
00146       ofs << "W(04)      " << 1. << std::endl;
00147       ofs << "Sp(05)     Vorticity-z" << std::endl;
00148       ofs << "W(05)      " << 1. << std::endl;
00149       ofs << "Sp(06)     |Vorticity|" << std::endl;
00150       ofs << "W(06)      " << 1. << std::endl;
00151       ofs.close();
00152     }    
00153   }
00154 
00155   inline virtual MacroType MacroVariables(const MicroType &f) const {
00156     MacroType q;
00157     q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8)+f(9)
00158         +f(10)+f(11)+f(12)+f(13)+f(14)+f(15)+f(16)+f(17)+f(18);
00159     q(1)=(f(1)-f(2)+f(7)-f(8)+f(9)-f(10)+f(15)-f(16)+f(17)-f(18))/q(0);
00160     q(2)=(f(3)-f(4)+f(7)-f(8)-f(9)+f(10)+f(11)-f(12)+f(13)-f(14))/q(0);
00161     q(3)=(f(5)-f(6)+f(11)-f(12)-f(13)+f(14)+f(15)-f(16)-f(17)+f(18))/q(0);
00162     return q;
00163   }
00164 
00165   inline virtual MicroType Equilibrium(const MacroType &q) const {
00166     MicroType feq;
00167 
00168     DataType rho = q(0);
00169     DataType u   = q(1);
00170     DataType v   = q(2);
00171     DataType w   = q(3);
00172 
00173     DataType ri, rt0, rt1, rt2;
00174     if (method[0]==0) {
00175       ri = DataType(1.);
00176       rt0 = R0 / DataType(3.);
00177       rt1 = R0 / DataType(18.);
00178       rt2 = R0 / DataType(36.);
00179     }
00180     if (method[0]==1) {
00181       ri = rho;
00182       rt0 = DataType(1.) / DataType(3.);
00183       rt1 = DataType(1.) / DataType(18.);
00184       rt2 = DataType(1.) / DataType(36.);
00185     }
00186     else if (method[0]==2) {
00187       ri = DataType(1.);
00188       rt0 = rho / DataType(3.);
00189       rt1 = rho / DataType(18.);
00190       rt2 = rho / DataType(36.);
00191     }
00192 
00193     DataType usq = u * u;
00194     DataType vsq = v * v;
00195     DataType wsq = w * w;
00196     DataType sumsq = (usq + vsq + wsq) / cs22;
00197     DataType sumsqxy = (usq + vsq) / cs22;
00198     DataType sumsqyz = (vsq + wsq) / cs22;
00199     DataType sumsqxz = (usq + wsq) / cs22;
00200     DataType sumsq2xy = sumsqxy * (DataType(1.) - cs2) / cs2;
00201     DataType sumsq2yz = sumsqyz * (DataType(1.) - cs2) / cs2;
00202     DataType sumsq2xz = sumsqxz * (DataType(1.) - cs2) / cs2;
00203     DataType u2 = usq / cssq;
00204     DataType v2 = vsq / cssq;
00205     DataType w2 = wsq / cssq;
00206     DataType u22 = usq / cs22;
00207     DataType v22 = vsq / cs22;
00208     DataType w22 = wsq / cs22;
00209     DataType ui = u / cs2;
00210     DataType vi = v / cs2;
00211     DataType wi = w / cs2;
00212     DataType uv = ui * vi;
00213     DataType uw = ui * wi;
00214     DataType vw = vi * wi;
00215 
00216     feq(0) = rt0 * (ri - sumsq);
00217     //x
00218     feq(1) = rt1 * (ri - sumsq + u2 + ui);
00219     feq(2) = rt1 * (ri - sumsq + u2 - ui);
00220     //y
00221     feq(3) = rt1 * (ri - sumsq + v2 + vi);
00222     feq(4) = rt1 * (ri - sumsq + v2 - vi);
00223     //z
00224     feq(5) = rt1 * (ri - sumsq + w2 + wi);
00225     feq(6) = rt1 * (ri - sumsq + w2 - wi);
00226     //x-y
00227     feq(7) = rt2 * (ri + sumsq2xy -w22 + ui + vi + uv);
00228     feq(8) = rt2 * (ri + sumsq2xy -w22 - ui - vi + uv);
00229     feq(9) = rt2 * (ri + sumsq2xy -w22 + ui - vi - uv);
00230     feq(10) = rt2 * (ri + sumsq2xy -w22 - ui + vi - uv);
00231     //y-z
00232     feq(11) = rt2 * (ri + sumsq2yz -u22 + vi + wi + vw);
00233     feq(12) = rt2 * (ri + sumsq2yz -u22 - vi - wi + vw);
00234     feq(13) = rt2 * (ri + sumsq2yz -u22 + vi - wi - vw);
00235     feq(14) = rt2 * (ri + sumsq2yz -u22 - vi + wi - vw);
00236     //x-z
00237     feq(15) = rt2 * (ri + sumsq2xz -v22 + ui + wi + uw);
00238     feq(16) = rt2 * (ri + sumsq2xz -v22 - ui - wi + uw);
00239     feq(17) = rt2 * (ri + sumsq2xz -v22 + ui - wi - uw);
00240     feq(18) = rt2 * (ri + sumsq2xz -v22 - ui + wi - uw);
00241 
00242     // Third order 
00243     if (method[3]>0) {
00244       DataType cs4 = DataType(1.)/(DataType(4.)*cs2*cs2);
00245       DataType sumsq3 = (usq + vsq + wsq) * cs4;
00246       DataType sumsqxy3 = (usq + vsq) * cs4;
00247       DataType sumsqyz3 = (vsq + wsq) * cs4;
00248       DataType sumsqxz3 = (usq + wsq) * cs4;
00249       DataType upv = (u+v)*(u+v) * cs4;
00250       DataType umv = (u-v)*(u-v) * cs4;
00251       DataType vpw = (v+w)*(v+w) * cs4;
00252       DataType vmw = (v-w)*(v-w) * cs4;
00253       DataType upw = (u+w)*(u+w) * cs4;
00254       DataType umw = (u-w)*(u-w) * cs4;
00255 
00256       //x
00257       feq(1) -= rt1 * u*sumsqyz3;
00258       feq(2) += rt1 * u*sumsqyz3;
00259       //y
00260       feq(3) -= rt1 * v*sumsqxz3;
00261       feq(4) += rt1 * v*sumsqxz3;
00262       //z
00263       feq(5) -= rt1 * w*sumsqxy3;
00264       feq(6) += rt1 * w*sumsqxy3;
00265       //x-y
00266       feq(7) -= rt2 * (u+v)*(sumsq3-upv);
00267       feq(8) += rt2 * (u+v)*(sumsq3-upv);
00268       feq(9) -= rt2 * (u-v)*(sumsq3-umv);
00269       feq(10) += rt2 * (u-v)*(sumsq3-umv);
00270       //y-z
00271       feq(11) -= rt2 * (v+w)*(sumsq3-vpw);
00272       feq(12) += rt2 * (v+w)*(sumsq3-vpw);
00273       feq(13) -= rt2 * (v-w)*(sumsq3-vmw);
00274       feq(14) += rt2 * (v-w)*(sumsq3-vmw);
00275       //x-z
00276       feq(15) -= rt2 * (u+w)*(sumsq3-upw);
00277       feq(16) += rt2 * (u+w)*(sumsq3-upw);
00278       feq(17) -= rt2 * (u-w)*(sumsq3-umw);
00279       feq(18) += rt2 * (u-w)*(sumsq3-umw);
00280     }
00281 
00282     return feq;
00283   }
00284 
00285   inline virtual void Collision(MicroType &f, const DataType dt) const {
00286     MacroType q = MacroVariables(f);
00287     MicroType feq = Equilibrium(q);
00288 
00289     DataType omega;
00290     if (turbulence == laminar) 
00291       omega = Omega(dt);
00292     else if (turbulence == LES_Smagorinsky) 
00293       omega = Omega_LES_Smagorinsky(f,feq,q,dt);
00294     else if (turbulence == LES_dynamic) 
00295       omega = Omega_LES_dynamic(f,feq,q,dt);
00296 
00297     f=(DataType(1.)-omega)*f + omega*feq;
00298   }     
00299 
00300   virtual int IncomingIndices(const int side, int indices[]) const {
00301     switch (side) {
00302     case base::Left:
00303       indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;  
00304       break;
00305     case base::Right:
00306       indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;   
00307       break;
00308     case base::Bottom:
00309       indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;  
00310       break;
00311     case base::Top:
00312       indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;  
00313       break;
00314     case base::Front:
00315       indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;  
00316       break;
00317     case base::Back:
00318       indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;  
00319       break;
00320      }
00321     return 5;
00322   }
00323 
00324   virtual int OutgoingIndices(const int side, int indices[]) const {
00325     switch (side) {
00326     case base::Left:
00327       indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;   
00328       break;
00329     case base::Right:
00330       indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;  
00331       break;
00332     case base::Bottom:
00333       indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;  
00334       break;
00335     case base::Top:
00336       indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;  
00337       break;
00338     case base::Front:
00339       indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;  
00340       break;
00341     case base::Back:
00342       indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;  
00343       break;
00344      }
00345     return 5;
00346   }
00347 
00348   // Revert just one layer of cells at the boundary 
00349   virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb, 
00350                              const int side) const {
00351     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00352     int b = fvec.bottom();
00353     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) 
00354             && fvec.stepsize(2)==bb.stepsize(2));
00355     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00356     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00357     int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00358     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00359     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00360     int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00361     MicroType *f = (MicroType *)fvec.databuffer();
00362 
00363 #ifdef DEBUG_PRINT_FIXUP
00364     ( comm_service::log() << "Reverse streaming at Side: " << side << " of " 
00365       << fvec.bbox() << " using " << bb << "\n").flush();
00366 #endif
00367     
00368     switch (side) {
00369     case base::Left:
00370       for (register int k=mzs; k<=mze; k++) 
00371         for (register int j=mys; j<=mye; j++) {
00372           f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](2);
00373           f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](8);
00374           f[b+base::idx(mxs,j,k,Nx,Ny)](10)= f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](10);
00375           f[b+base::idx(mxs,j,k,Nx,Ny)](16)= f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](16);
00376           f[b+base::idx(mxs,j,k,Nx,Ny)](18)= f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](18);
00377         }
00378       break;
00379     case base::Right:
00380       for (register int k=mzs; k<=mze; k++) 
00381         for (register int j=mys; j<=mye; j++) {
00382           f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](1);
00383           f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](7);
00384           f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](9);
00385           f[b+base::idx(mxe,j,k,Nx,Ny)](15)= f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](15);
00386           f[b+base::idx(mxe,j,k,Nx,Ny)](17)= f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](17);
00387         }
00388       break;
00389     case base::Bottom:
00390       for (register int k=mzs; k<=mze; k++) 
00391         for (register int i=mxs; i<=mxe; i++) {
00392           f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](4);
00393           f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](8);
00394           f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](9);
00395           f[b+base::idx(i,mys,k,Nx,Ny)](12)= f[b+base::idx(i,mys-1,k-1,Nx,Ny)](12);
00396           f[b+base::idx(i,mys,k,Nx,Ny)](14)= f[b+base::idx(i,mys-1,k+1,Nx,Ny)](14);
00397         }
00398       break;
00399     case base::Top:
00400       for (register int k=mzs; k<=mze; k++) 
00401         for (register int i=mxs; i<=mxe; i++) {
00402           f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](3);
00403           f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](7);
00404           f[b+base::idx(i,mye,k,Nx,Ny)](10)= f[b+base::idx(i-1,mye+1,k,Nx,Ny)](10);
00405           f[b+base::idx(i,mye,k,Nx,Ny)](11)= f[b+base::idx(i,mye+1,k+1,Nx,Ny)](11);
00406           f[b+base::idx(i,mye,k,Nx,Ny)](13)= f[b+base::idx(i,mye+1,k-1,Nx,Ny)](13);
00407         }
00408       break;
00409     case base::Front:
00410       for (register int j=mys; j<=mye; j++) 
00411         for (register int i=mxs; i<=mxe; i++) {
00412           f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](6);
00413           f[b+base::idx(i,j,mzs,Nx,Ny)](12)= f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](12);
00414           f[b+base::idx(i,j,mzs,Nx,Ny)](13)= f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](13);
00415           f[b+base::idx(i,j,mzs,Nx,Ny)](16)= f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](16);
00416           f[b+base::idx(i,j,mzs,Nx,Ny)](17)= f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](17);
00417         }
00418       break;
00419     case base::Back:
00420       for (register int j=mys; j<=mye; j++) 
00421         for (register int i=mxs; i<=mxe; i++) {
00422           f[b+base::idx(i,j,mzs,Nx,Ny)](5) = f[b+base::idx(i,j,mzs+1,Nx,Ny)](5);
00423           f[b+base::idx(i,j,mzs,Nx,Ny)](11)= f[b+base::idx(i,j+1,mzs+1,Nx,Ny)](11);
00424           f[b+base::idx(i,j,mzs,Nx,Ny)](14)= f[b+base::idx(i,j-1,mzs+1,Nx,Ny)](14);
00425           f[b+base::idx(i,j,mzs,Nx,Ny)](15)= f[b+base::idx(i+1,j,mzs+1,Nx,Ny)](15);
00426           f[b+base::idx(i,j,mzs,Nx,Ny)](18)= f[b+base::idx(i-1,j,mzs+1,Nx,Ny)](18);
00427         }
00428       break;
00429     }
00430   }   
00431 
00432   virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec, 
00433                          const BBox &bb, const double& dt) const {
00434     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00435     int b = fvec.bottom();
00436     MicroType *f = (MicroType *)fvec.databuffer();
00437     MicroType *of = (MicroType *)ovec.databuffer();
00438 
00439 #ifdef DEBUG_PRINT_FIXUP
00440     ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox() 
00441       << " using " << bb << "\n").flush();
00442 #endif
00443     
00444     assert (fvec.extents(0)==ovec.extents(0));
00445     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) && fvec.stepsize(2)==bb.stepsize(2) &&
00446             fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1) && fvec.stepsize(2)==ovec.stepsize(2));
00447     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00448     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00449     int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00450     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00451     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00452     int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00453 
00454     for (register int k=mzs; k<=mze; k++) 
00455       for (register int j=mys; j<=mye; j++) 
00456         for (register int i=mxs; i<=mxe; i++) {
00457           for (register int n=0; n<base::NMicroVar(); n++) 
00458             f[b+base::idx(i,j,k,Nx,Ny)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],k-mdz[n],Nx,Ny)](n);
00459           Collision(f[b+base::idx(i,j,k,Nx,Ny)],dt);
00460         }
00461   }  
00462 
00463   virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00464                       const double& t, const double& dt, const int& mpass) const {
00465     // Make sure that the physical length and times are scaling consistently
00466     // into lattice quantities
00467     DCoords dx = base::GH().worldStep(fvec.stepsize());
00468     
00469     //std::cout << "STEP   dt=" << dt << "   TO=" << T0() << std::endl;
00470     if ( (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR)
00471          || (std::fabs(dx(1)/L0()-dx(2)/L0()) > DBL_EPSILON*LB_FACTOR) ) {
00472       std::cerr << "LBMD3Q19 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00473                 << " and dx(2)=" << dx(2)/L0()
00474                 << " to be equal." << std::endl << std::flush;
00475       std::exit(-1);
00476     }
00477     if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR)   {
00478       std::cerr << "LBMD3Q19 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
00479                 << " to be equal."
00480                 << "   dt=" << dt << "   TO=" << T0()
00481                 << std::endl << std::flush;
00482       std::exit(-1);
00483     }
00484 
00485     int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00486     MicroType *f = (MicroType *)fvec.databuffer();
00487 
00488     // Streaming - update all cells
00489     for (register int k=Nz-1; k>=1; k--)
00490       for (register int j=Ny-1; j>=1; j--)
00491         for (register int i=0; i<=Nx-2; i++) {
00492           //x-y
00493           f[base::idx(i,j,k,Nx,Ny)](3) = f[base::idx(i,j-1,k,Nx,Ny)](3);
00494           f[base::idx(i,j,k,Nx,Ny)](10) = f[base::idx(i+1,j-1,k,Nx,Ny)](10);
00495           //y-z
00496           //f[base::idx(i,j,k,Nx,Ny)](13) = f[base::idx(i,j-1,k+1,Nx,Ny)](13);//%%%%% -j,+k
00497           //x-z
00498           f[base::idx(i,j,k,Nx,Ny)](5) = f[base::idx(i,j,k-1,Nx,Ny)](5);
00499           f[base::idx(i,j,k,Nx,Ny)](18) = f[base::idx(i+1,j,k-1,Nx,Ny)](18);
00500         }
00501 
00502     for (register int k=0; k<=Nz-2; k++)
00503       for (register int j=Ny-1; j>=1; j--)
00504         for (register int i=0; i<=Nx-2; i++) {
00505           //y-z
00506           f[base::idx(i,j,k,Nx,Ny)](13) = f[base::idx(i,j-1,k+1,Nx,Ny)](13);//%%%%% -j,+k
00507         }
00508     for (register int k=Nz-1; k>=1; k--)
00509       for (register int j=Ny-1; j>=1; j--)
00510         for (register int i=Nx-1; i>=1; i--) {
00511           //x-y
00512           f[base::idx(i,j,k,Nx,Ny)](1) = f[base::idx(i-1,j,k,Nx,Ny)](1);
00513           f[base::idx(i,j,k,Nx,Ny)](7) = f[base::idx(i-1,j-1,k,Nx,Ny)](7);
00514           //y-z
00515           f[base::idx(i,j,k,Nx,Ny)](11) = f[base::idx(i,j-1,k-1,Nx,Ny)](11);
00516           //x-z
00517           f[base::idx(i,j,k,Nx,Ny)](15) = f[base::idx(i-1,j,k-1,Nx,Ny)](15);
00518         }
00519 
00520     for (register int k=0; k<=Nz-2; k++)
00521       for (register int j=0; j<=Ny-2; j++)
00522         for (register int i=Nx-1; i>=1; i--) {
00523           //x-y
00524           f[base::idx(i,j,k,Nx,Ny)](4) = f[base::idx(i,j+1,k,Nx,Ny)](4);
00525           f[base::idx(i,j,k,Nx,Ny)](9) = f[base::idx(i-1,j+1,k,Nx,Ny)](9);
00526           //y-z
00527           //f[base::idx(i,j,k,Nx,Ny)](14) = f[base::idx(i,j+1,k-1,Nx,Ny)](14);//%%%%% +j,-k
00528           //x-z
00529           f[base::idx(i,j,k,Nx,Ny)](6) = f[base::idx(i,j,k+1,Nx,Ny)](6);
00530           f[base::idx(i,j,k,Nx,Ny)](17) = f[base::idx(i-1,j,k+1,Nx,Ny)](17);
00531         }
00532     for (register int k=Nz-1; k>=1; k--)
00533       for (register int j=0; j<=Ny-2; j++)
00534         for (register int i=0; i<=Nx-2; i++) {
00535           //y-z
00536           f[base::idx(i,j,k,Nx,Ny)](14) = f[base::idx(i,j+1,k-1,Nx,Ny)](14);//%%%%% +j,-k
00537         }
00538     for (register int k=0; k<=Nz-2; k++)
00539       for (register int j=0; j<=Ny-2; j++)
00540         for (register int i=0; i<=Nx-2; i++) {
00541           //x-y
00542           f[base::idx(i,j,k,Nx,Ny)](2) = f[base::idx(i+1,j,k,Nx,Ny)](2);
00543           f[base::idx(i,j,k,Nx,Ny)](8) = f[base::idx(i+1,j+1,k,Nx,Ny)](8);
00544           //y-z
00545           f[base::idx(i,j,k,Nx,Ny)](12) = f[base::idx(i,j+1,k+1,Nx,Ny)](12);
00546           //x-z
00547           f[base::idx(i,j,k,Nx,Ny)](16) = f[base::idx(i+1,j,k+1,Nx,Ny)](16);
00548         }
00549 
00550      // Collision - only in the interior region
00551     int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00552     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00553     for (register int k=mzs; k<=mze; k++)
00554       for (register int j=mys; j<=mye; j++)
00555         for (register int i=mxs; i<=mxe; i++)
00556           Collision(f[base::idx(i,j,k,Nx,Ny)],dt);
00557 
00558     return DataType(1.);
00559   }
00560 
00561   virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00562                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00563     int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00564     int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00565     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00566     MicroType *f = (MicroType *)fvec.databuffer();
00567 
00568     //std::cout << "IC Nx=" << Nx << "   Ny=" << Ny << "   Nz=" << Nz << "   type=" << type << std::endl;
00569 
00570     if (type==ConstantMicro) {
00571       //std::cout << "ConstantMicro start\n";
00572       MicroType g;
00573       for (register int i=0; i<base::NMicroVar(); i++)
00574         g(i) = aux[i];
00575 
00576       for (register int k=mzs; k<=mze; k++)
00577         for (register int j=mys; j<=mye; j++)
00578           for (register int i=mxs; i<=mxe; i++)
00579             f[base::idx(i,j,k,Nx,Ny)] = g;
00580       //std::cout << "ConstantMicro end\n";
00581     }
00582     
00583     else if (type==ConstantMacro) {
00584       //std::cout << "ConstantMacro start\n";
00585       MacroType q;
00586       if (scaling==base::Physical) {
00587         q(0) = aux[0]/R0;
00588         q(1) = aux[1]/U0;
00589         q(2) = aux[2]/U0;
00590         q(3) = aux[3]/U0;
00591       }
00592       else
00593         for (register int i=0; i<base::NMacroVar(); i++)
00594           q(i) = aux[i];
00595       q(1) *= S0; q(2) *= S0; q(3) *= S0;
00596       
00597 
00598       /*std::cout << "aux[0]=" << aux[0] << "  aux[1]=" << aux[1]
00599   << "  aux[2]=" << aux[2] << "  aux[3]=" << aux[3]
00600   << "  equlibrium=" << Equilibrium(q)
00601   << "  q=" <<  MacroVariables(Equilibrium(q)) << std::endl;*/
00602 
00603       for (register int k=mzs; k<=mze; k++)
00604         for (register int j=mys; j<=mye; j++)
00605           for (register int i=mxs; i<=mxe; i++)
00606             f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00607       //std::cout << "ConstantMacro end\n";
00608     }
00609     
00610     else if (type==GasAtRest) {
00611       //std::cout << "GasAtRest start\n";
00612       MacroType q;
00613       q(0) = GasDensity()/R0;
00614       q(1) = DataType(0.);
00615       q(2) = DataType(0.);
00616       q(3) = DataType(0.);
00617       for (register int k=mzs; k<=mze; k++)
00618         for (register int j=mys; j<=mye; j++)
00619           for (register int i=mxs; i<=mxe; i++)
00620             f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00621       //std::cout << "GasAtRest end\n";
00622     }
00623   }
00624 
00625   // Set just the one layer of ghost cells at the boundary
00626   virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00627                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00628     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00629     int b = fvec.bottom();
00630     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) 
00631             && fvec.stepsize(2)==bb.stepsize(2));
00632     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00633     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00634     int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00635     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00636     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00637     int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00638     MicroType *f = (MicroType *)fvec.databuffer();
00639     
00640     //std::cout << "BC Nx=" << Nx << "   Ny=" << Ny << "   Nz=" << Nz << std::endl;
00641     if (type==Symmetry) {
00642       switch (side) {
00643       case base::Left:
00644         for (register int k=mzs; k<=mze; k++)
00645           for (register int j=mys; j<=mye; j++) {
00646             if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
00647             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00648             if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);
00649 
00650             if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
00651             if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
00652           }
00653         break;
00654       case base::Right:
00655         for (register int k=mzs; k<=mze; k++)
00656           for (register int j=mys; j<=mye; j++) {
00657             if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
00658             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00659             if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
00660 
00661             if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
00662             if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
00663           }
00664         break;
00665       case base::Bottom:
00666         for (register int k=mzs; k<=mze; k++)
00667           for (register int i=mxs; i<=mxe; i++) {
00668             if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
00669             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00670             if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
00671 
00672             if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
00673             if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
00674           }
00675         break;
00676       case base::Top:
00677         for (register int k=mzs; k<=mze; k++)
00678           for (register int i=mxs; i<=mxe; i++) {
00679             if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
00680             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00681             if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
00682 
00683             if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
00684             if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
00685           }
00686         break;
00687       case base::Front:
00688         //std::cout << "SYMMETRY Front Nz=" << Nz << std::endl;
00689         for (register int j=mys; j<=mye; j++)
00690           for (register int i=mxs; i<=mxe; i++) {
00691             if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
00692             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00693             if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
00694 
00695             if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
00696             if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
00697           }
00698         break;
00699       case base::Back:
00700         //std::cout << "SYMMETRY Back Nz=" << Nz << std::endl;
00701         for (register int j=mys; j<=mye; j++)
00702           for (register int i=mxs; i<=mxe; i++) {
00703             /*std::cout << "NS Back Nz=" << Nz << "  {"<< i << "," << j << "," << mzs << "}"
00704         << " mxs=" << mxs << " mxe=" << mxe
00705         << " mys=" << mys << " mye=" << mye
00706         << " mzs=" << mzs << " mze=" << mze << std::endl;*/
00707             if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
00708             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00709             if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
00710 
00711             if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
00712             if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
00713           }
00714         break;
00715       }
00716     }
00717     
00718     else if (type==SlipWall) {
00719       switch (side) {
00720       case base::Left:
00721         //std::cout << "S Left Nz=" << Nz << std::endl;
00722         for (register int k=mzs; k<=mze; k++)
00723           for (register int j=mys; j<=mye; j++) {
00724             if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
00725             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00726             if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
00727 
00728             if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
00729             if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
00730           }
00731         break;
00732       case base::Right:
00733         //std::cout << "S Right Nz=" << Nz << std::endl;
00734         for (register int k=mzs; k<=mze; k++)
00735           for (register int j=mys; j<=mye; j++) {
00736             if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
00737             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00738             if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
00739 
00740             if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
00741             if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
00742           }
00743         break;
00744       case base::Bottom:
00745         //std::cout << "S Bottom Nz=" << Nz << std::endl;
00746         for (register int k=mzs; k<=mze; k++)
00747           for (register int i=mxs; i<=mxe; i++) {
00748             if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
00749             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00750             if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
00751 
00752             if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
00753             if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
00754           }
00755         break;
00756       case base::Top:
00757         //std::cout << "S Top Nz=" << Nz << std::endl;
00758         for (register int k=mzs; k<=mze; k++)
00759           for (register int i=mxs; i<=mxe; i++) {
00760             if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
00761             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00762             if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
00763 
00764             if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
00765             if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
00766           }
00767         break;
00768       case base::Front:
00769         //std::cout << "S Front Nz=" << Nz << std::endl;
00770         for (register int j=mys; j<=mye; j++)
00771           for (register int i=mxs; i<=mxe; i++) {
00772             if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
00773             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00774             if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
00775 
00776             if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
00777             if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
00778           }
00779         break;
00780       case base::Back:
00781         //std::cout << "S Back Nz=" << Nz << std::endl;
00782         for (register int j=mys; j<=mye; j++)
00783           for (register int i=mxs; i<=mxe; i++) {
00784             /*std::cout << "NS Back Nz=" << Nz << "  {"<< i << "," << j << "," << mzs << "}"
00785         << " mxs=" << mxs << " mxe=" << mxe
00786         << " mys=" << mys << " mye=" << mye
00787         << " mzs=" << mzs << " mze=" << mze << std::endl;*/
00788             if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
00789             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00790             if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
00791 
00792             if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
00793             if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
00794           }
00795         break;
00796       }
00797     }
00798 
00799     else if (type==NoSlipWall) {
00800       switch (side) {
00801       case base::Left:
00802         //std::cout << "NS Left Nz=" << Nz << std::endl;
00803         for (register int k=mzs; k<=mze; k++)
00804           for (register int j=mys; j<=mye; j++) {
00805             if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
00806             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00807             if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
00808 
00809             if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
00810             if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
00811           }
00812         break;
00813       case base::Right:
00814         //std::cout << "NS Right Nz=" << Nz << std::endl;
00815         for (register int k=mzs; k<=mze; k++)
00816           for (register int j=mys; j<=mye; j++) {
00817             if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
00818             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00819             if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
00820 
00821             if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
00822             if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
00823           }
00824         break;
00825       case base::Bottom:
00826         //std::cout << "NS Bottom Nz=" << Nz << std::endl;
00827         for (register int k=mzs; k<=mze; k++)
00828           for (register int i=mxs; i<=mxe; i++) {
00829             if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
00830             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00831             if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
00832 
00833             if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
00834             if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
00835           }
00836         break;
00837       case base::Top:
00838         //std::cout << "NS Top Nz=" << Nz << std::endl;
00839         for (register int k=mzs; k<=mze; k++)
00840           for (register int i=mxs; i<=mxe; i++) {
00841             if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
00842             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00843             if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
00844 
00845             if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
00846             if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
00847           }
00848         break;
00849       case base::Front:
00850         //std::cout << "NS Front Nz=" << Nz << std::endl;
00851         for (register int j=mys; j<=mye; j++)
00852           for (register int i=mxs; i<=mxe; i++) {
00853             if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
00854             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00855             if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
00856 
00857             if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
00858             if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
00859           }
00860         break;
00861       case base::Back:
00862         //std::cout << "NS Back Nz=" << Nz << std::endl;
00863         for (register int j=mys; j<=mye; j++)
00864           for (register int i=mxs; i<=mxe; i++) {
00865             /*std::cout << "NS Back Nz=" << Nz << "  {"<< i << "," << j << "," << mzs << "}"
00866         << " mxs=" << mxs << " mxe=" << mxe
00867         << " mys=" << mys << " mye=" << mye
00868         << " mzs=" << mzs << " mze=" << mze << std::endl;*/
00869             if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
00870             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00871             if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
00872 
00873             if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
00874             if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
00875           }
00876         break;
00877       }
00878     }
00879 
00880     else if (type==Inlet) {
00881       if (aux==0 || naux<=0) return;
00882       DataType a0=S0*aux[0], a1=0., a2=0.;
00883       if (naux>1) {
00884         a1 = S0*aux[1];
00885         a2 = S0*aux[2];
00886       }
00887       if (scaling==base::Physical) {
00888         a0 /= U0;
00889         a1 /= U0;
00890         a2 /= U0;
00891       }
00892       switch (side) {
00893       case base::Left:
00894         for (register int k=mzs; k<=mze; k++)
00895           for (register int j=mys; j<=mye; j++) {
00896             MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00897             q(1) = a0;
00898             if (naux>1) {
00899               q(2) = a1;
00900               q(3) = a2;
00901             }
00902             f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
00903           }
00904         break;
00905       case base::Right:
00906         for (register int k=mzs; k<=mze; k++)
00907           for (register int j=mys; j<=mye; j++) {
00908             MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
00909             q(1) = a0;
00910             if (naux>1) {
00911               q(2) = a1;
00912               q(3) = a2;
00913             }
00914             f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
00915           }
00916         break;
00917       case base::Bottom:
00918         for (register int k=mzs; k<=mze; k++)
00919           for (register int i=mxs; i<=mxe; i++) {
00920             MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
00921             if (naux==1)
00922               q(2) = a0;
00923             else {
00924               q(1) = a0;
00925               q(2) = a1;
00926               q(3) = a2;
00927             }
00928             f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
00929           }
00930         break;
00931       case base::Top:
00932         for (register int k=mzs; k<=mze; k++)
00933           for (register int i=mxs; i<=mxe; i++) {
00934             MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
00935             if (naux==1)
00936               q(2) = a0;
00937             else {
00938               q(1) = a0;
00939               q(2) = a1;
00940               q(3) = a2;
00941             }
00942             f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
00943           }
00944         break;
00945       case base::Front:
00946         for (register int j=mys; j<=mye; j++)
00947           for (register int i=mxs; i<=mxe; i++) {
00948             MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
00949             if (naux==1)
00950               q(3) = a0;
00951             else {
00952               q(1) = a1;
00953               q(2) = a2;
00954               q(3) = a0;
00955             }
00956             f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
00957           }
00958         break;
00959       case base::Back:
00960         for (register int j=mys; j<=mye; j++)
00961           for (register int i=mxs; i<=mxe; i++) {
00962             MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
00963             if (naux==1)
00964               q(3) = a0;
00965             else {
00966               q(1) = a1;
00967               q(2) = a2;
00968               q(3) = a0;
00969             }
00970             f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
00971           }
00972         break;
00973       }
00974     }
00975 
00976     else if (type==Outlet) {
00977       //std::cout << "Outlet" << std::endl;
00978       switch (side) {
00979       case base::Left:
00980         //std::cout << "Outlet Left" << std::endl;
00981         for (register int k=mzs; k<=mze; k++)
00982           for (register int j=mys; j<=mye; j++) {
00983             MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00984             //std::cout << "Outlet Left q(0)=" << q(0) << std::endl;
00985             if (q(0)>0) f[b+base::idx(mxe,j,k,Nx,Ny)] = f[b+base::idx(mxe+1,j,k,Nx,Ny)]/q(0);
00986           }
00987         break;
00988       case base::Right:
00989         //std::cout << "Outlet Right" << std::endl;
00990         for (register int k=mzs; k<=mze; k++)
00991           for (register int j=mys; j<=mye; j++) {
00992             MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
00993             //std::cout << "Outlet Right q(0)=" << q(0) << std::endl;
00994             if (q(0)>0) f[b+base::idx(mxs,j,k,Nx,Ny)] = f[b+base::idx(mxs-1,j,k,Nx,Ny)]/q(0);
00995           }
00996         break;
00997       case base::Bottom:
00998         for (register int k=mzs; k<=mze; k++)
00999           for (register int i=mxs; i<=mxe; i++) {
01000             MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01001             if (q(0)>0) f[b+base::idx(i,mye,k,Nx,Ny)] = f[b+base::idx(i,mye+1,k,Nx,Ny)]/q(0);
01002           }
01003         break;
01004       case base::Top:
01005         for (register int k=mzs; k<=mze; k++)
01006           for (register int i=mxs; i<=mxe; i++) {
01007             MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01008             if (q(0)>0) f[b+base::idx(i,mys,k,Nx,Ny)] = f[b+base::idx(i,mys-1,k,Nx,Ny)]/q(0);
01009           }
01010         break;
01011       case base::Front:
01012         for (register int j=mys; j<=mye; j++)
01013           for (register int i=mxs; i<=mxe; i++) {
01014             MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01015             if (q(0)>0) f[b+base::idx(i,j,mze,Nx,Ny)] = f[b+base::idx(i,j,mze+1,Nx,Ny)]/q(0);
01016           }
01017         break;
01018       case base::Back:
01019         for (register int j=mys; j<=mye; j++)
01020           for (register int i=mxs; i<=mxe; i++) {
01021             MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01022             if (q(0)>0) f[b+base::idx(i,j,mzs,Nx,Ny)] = f[b+base::idx(i,j,mzs-1,Nx,Ny)]/q(0);
01023           }
01024         break;
01025       }
01026     }
01027     
01028     // (Succi pg. 90 and Chen_Martinez_Mei_PhysFluids_8_2527)
01029     else if (type==Pressure) {
01030       if (aux==0 || naux<=0) return;
01031       DataType a0=aux[0];
01032       if (scaling==base::Physical)
01033         a0 /= R0;
01034       switch (side) {
01035       case base::Left:
01036         for (register int k=mzs; k<=mze; k++)
01037           for (register int j=mys; j<=mye; j++) {
01038             MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01039             q(0) = a0;
01040             f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
01041           }
01042         break;
01043       case base::Right:
01044         for (register int k=mzs; k<=mze; k++)
01045           for (register int j=mys; j<=mye; j++) {
01046             MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01047             q(0) = a0;
01048             f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
01049           }
01050         break;
01051       case base::Bottom:
01052         for (register int k=mzs; k<=mze; k++)
01053           for (register int i=mxs; i<=mxe; i++) {
01054             MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01055             q(0) = a0;
01056             f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
01057           }
01058         break;
01059       case base::Top:
01060         for (register int k=mzs; k<=mze; k++)
01061           for (register int i=mxs; i<=mxe; i++) {
01062             MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01063             q(0) = a0;
01064             f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
01065           }
01066         break;
01067       case base::Front:
01068         for (register int j=mys; j<=mye; j++)
01069           for (register int i=mxs; i<=mxe; i++) {
01070             MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01071             q(0) = a0;
01072             f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
01073           }
01074         break;
01075       case base::Back:
01076         for (register int j=mys; j<=mye; j++)
01077           for (register int i=mxs; i<=mxe; i++) {
01078             MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01079             q(0) = a0;
01080             f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
01081           }
01082         break;
01083       }
01084     }
01085     
01086     else if (type==SlidingWall) {
01087       if (aux==0 || naux<=0) return;
01088       if (aux==0 || naux<=0) return;
01089       DataType a0=S0*aux[0];
01090       if (scaling==base::Physical)
01091         a0 /= U0;
01092       switch (side) {
01093       case base::Left:
01094         for (register int k=mzs; k<=mze; k++)
01095           for (register int j=mys; j<=mye; j++) {
01096             MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01097             DataType rd1 = f[b+base::idx(mxe+1,j,k,Nx,Ny)](10), rd2 = f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
01098             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01099             DataType pw = DataType(1.)-qw;
01100             if (j<mye) f[b+base::idx(mxe,j+1,k,Nx,Ny)](9) = pw*rd1+qw*rd2;
01101             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
01102             if (j>mys) f[b+base::idx(mxe,j-1,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01103 
01104             if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01105             if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01106           }
01107         break;
01108       case base::Right:
01109         for (register int k=mzs; k<=mze; k++)
01110           for (register int j=mys; j<=mye; j++) {
01111             MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01112             DataType rd1 = f[b+base::idx(mxs-1,j,k,Nx,Ny)](7), rd2 = f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
01113             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01114             DataType pw = DataType(1.)-qw;
01115             if (j<mye) f[b+base::idx(mxs,j+1,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01116             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
01117             if (j>mys) f[b+base::idx(mxs,j-1,k,Nx,Ny)](10) = qw*rd1+pw*rd2;
01118 
01119             if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
01120             if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
01121           }
01122         break;
01123       case base::Bottom:
01124         for (register int k=mzs; k<=mze; k++)
01125           for (register int i=mxs; i<=mxe; i++) {
01126             MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01127             DataType rd1 = f[b+base::idx(i,mye+1,k,Nx,Ny)](9), rd2 = f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
01128             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01129             DataType pw = DataType(1.)-qw;
01130             if (i<mxe) f[b+base::idx(i+1,mye,k,Nx,Ny)](10) = pw*rd1+qw*rd2;
01131             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
01132             if (i>mxs) f[b+base::idx(i-1,mye,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01133 
01134             if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
01135             if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
01136           }
01137         break;
01138       case base::Top:
01139         for (register int k=mzs; k<=mze; k++)
01140           for (register int i=mxs; i<=mxe; i++) {
01141             MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01142             DataType rd1 = f[b+base::idx(i,mys-1,k,Nx,Ny)](7), rd2 = f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
01143             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01144             DataType pw = DataType(1.)-qw;
01145             if (i<mxe) f[b+base::idx(i+1,mys,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01146             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
01147             if (i>mxs) f[b+base::idx(i-1,mys,k,Nx,Ny)](9) = qw*rd1+pw*rd2;
01148 
01149             if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
01150             if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
01151           }
01152         break;
01153       case base::Front:
01154         for (register int j=mys; j<=mye; j++)
01155           for (register int i=mxs; i<=mxe; i++) {
01156             MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01157             DataType rd1 = f[b+base::idx(i,j,mze+1,Nx,Ny)](17), rd2 = f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
01158             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01159             DataType pw = DataType(1.)-qw;
01160             if (i<mxe) f[b+base::idx(i+1,j,mze,Nx,Ny)](18) = pw*rd1+qw*rd2;
01161             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
01162             if (i>mxs) f[b+base::idx(i-1,j,mze,Nx,Ny)](15) = qw*rd1+pw*rd2;
01163 
01164             if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
01165             if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
01166           }
01167         break;
01168       case base::Back:
01169         for (register int j=mys; j<=mye; j++)
01170           for (register int i=mxs; i<=mxe; i++) {
01171             MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01172             DataType rd1 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](15), rd2 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
01173             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01174             DataType pw = DataType(1.)-qw;
01175             if (i<mxe) f[b+base::idx(i+1,j,mzs,Nx,Ny)](16) = pw*rd1+qw*rd2;
01176             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
01177             if (i>mxs) f[b+base::idx(i-1,j,mzs,Nx,Ny)](17) = qw*rd1+pw*rd2;
01178 
01179             if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
01180             if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
01181           }
01182         break;
01183       }
01184     }
01185 
01189     else if (type==CharacteristicOutlet) {
01190       DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
01191       DCoords dx = base::GH().worldStep(fvec.stepsize());
01192       DataType dt_temp = dx(0)/VelocityScale();
01193       switch (side) {
01194       case base::Left:
01195         for (register int k=mzs; k<=mze; k++)
01196           for (register int j=mys; j<=mye; j++) {
01197             MicroType ftmp = f[b+base::idx(mxe,j,k,Nx,Ny)];
01198             MacroType q = MacroVariables(ftmp);
01199             Collision(ftmp,dt_temp);
01200             MacroType q1 = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01201             MacroType q2 = MacroVariables(f[b+base::idx(mxe+2,j,k,Nx,Ny)]);
01202             MacroType dqdn = NormalDerivative(q,q1,q2);
01203             Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01204 
01205             // Solve LODI equations :
01206             MacroType qn;
01207             qn(0) = q(0) - c1oCsSq2*L(0);
01208             qn(1) = q(1) + c1o2cs*L(0)/q(0);
01209             qn(2) = q(2) - L(1);
01210             qn(3) = q(3) - L(2);
01211 
01212             f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(qn);
01213           }
01214         break;
01215       case base::Right:
01216         for (register int k=mzs; k<=mze; k++)
01217           for (register int j=mys; j<=mye; j++) {
01218             MicroType ftmp = f[b+base::idx(mxs,j,k,Nx,Ny)];
01219             MacroType q = MacroVariables(ftmp);
01220             Collision(ftmp,dt_temp);
01221             MacroType q1 = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01222             MacroType q2 = MacroVariables(f[b+base::idx(mxs-2,j,k,Nx,Ny)]);
01223             MacroType dqdn = NormalDerivative(q,q1,q2);
01224             Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01225 
01226             // Solve LODI equations :
01227             MacroType qn;
01228             qn(0) = q(0) - c1oCsSq2*L(0);
01229             qn(1) = q(1) + c1o2cs*L(0)/q(0);
01230             qn(2) = q(2) - L(1);
01231             qn(3) = q(3) - L(2);
01232 
01233             f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(qn);
01234           }
01235         break;
01236       case base::Bottom:
01237         for (register int k=mzs; k<=mze; k++)
01238           for (register int i=mxs; i<=mxe; i++) {
01239             MicroType ftmp = f[b+base::idx(i,mye,k,Nx,Ny)];
01240             MacroType q = MacroVariables(ftmp);
01241             Collision(ftmp,dt_temp);
01242             MacroType q1 = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01243             MacroType q2 = MacroVariables(f[b+base::idx(i,mye+2,k,Nx,Ny)]);
01244             MacroType dqdn = NormalDerivative(q,q1,q2);
01245             Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01246 
01247             // Solve LODI equations :
01248             MacroType qn;
01249             qn(0) = q(0) - c1oCsSq2*L(0);
01250             qn(1) = q(1) - L(1);
01251             qn(2) = q(2) + c1o2cs*L(0)/q(0);
01252             qn(3) = q(3) - L(2);
01253 
01254             f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(qn);
01255           }
01256         break;
01257       case base::Top:
01258         for (register int k=mzs; k<=mze; k++)
01259           for (register int i=mxs; i<=mxe; i++) {
01260             MicroType ftmp = f[b+base::idx(i,mys,k,Nx,Ny)];
01261             MacroType q = MacroVariables(ftmp);
01262             Collision(ftmp,dt_temp);
01263             MacroType q1 = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01264             MacroType q2 = MacroVariables(f[b+base::idx(i,mys-2,k,Nx,Ny)]);
01265             MacroType dqdn = NormalDerivative(q,q1,q2);
01266             Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01267 
01268             // Solve LODI equations :
01269             MacroType qn;
01270             qn(0) = q(0) - c1oCsSq2*L(0);
01271             qn(1) = q(1) - L(1);
01272             qn(2) = q(2) + c1o2cs*L(0)/q(0);
01273             qn(3) = q(3) - L(2);
01274 
01275             f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(qn);
01276           }
01277         break;
01278       case base::Front:
01279         for (register int j=mys; j<=mye; j++)
01280           for (register int i=mxs; i<=mxe; i++) {
01281             MicroType ftmp = f[b+base::idx(i,j,mze,Nx,Ny)];
01282             MacroType q = MacroVariables(ftmp);
01283             Collision(ftmp,dt_temp);
01284             MacroType q1 = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01285             MacroType q2 = MacroVariables(f[b+base::idx(i,j,mze+2,Nx,Ny)]);
01286             MacroType dqdn = NormalDerivative(q,q1,q2);
01287             Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01288 
01289             // Solve LODI equations :
01290             MacroType qn;
01291             qn(0) = q(0) - c1oCsSq2*L(0);
01292             qn(1) = q(1) - L(1);
01293             qn(2) = q(2) - L(2);
01294             qn(3) = q(3) + c1o2cs*L(0)/q(0);
01295 
01296             f[b+base::idx(i,j,mze,Nx,Ny)] = Equilibrium(qn);
01297           }
01298         break;
01299       case base::Back:
01300         for (register int j=mys; j<=mye; j++)
01301           for (register int i=mxs; i<=mxe; i++) {
01302             MicroType ftmp = f[b+base::idx(i,j,mzs,Nx,Ny)];
01303             MacroType q = MacroVariables(ftmp);
01304             Collision(ftmp,dt_temp);
01305             MacroType q1 = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01306             MacroType q2 = MacroVariables(f[b+base::idx(i,j,mzs-2,Nx,Ny)]);
01307             MacroType dqdn = NormalDerivative(q,q1,q2);
01308             Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01309 
01310             // Solve LODI equations :
01311             MacroType qn;
01312             qn(0) = q(0) - c1oCsSq2*L(0);
01313             qn(1) = q(1) - L(1);
01314             qn(2) = q(2) - L(2);
01315             qn(3) = q(3) + c1o2cs*L(0)/q(0);
01316 
01317             f[b+base::idx(i,j,mzs,Nx,Ny)] = Equilibrium(qn);
01318           }
01319         break;
01320       }
01321     }
01322 
01326     else if (type==CharacteristicInlet) {
01327       if (aux==0 || naux<=0) return;
01328       DataType rho=aux[0];
01329       DataType a0=S0*aux[1];
01330       DataType a1=S0*aux[2];
01331       DataType a2=S0*aux[3];
01332       if (scaling==base::Physical) {
01333         rho /= R0;
01334         a0 /= U0;
01335         a1 /= U0;
01336         a2 /= U0;
01337       }
01338       MacroType q, qn;
01339       q(0) = rho;
01340       q(1) = a0;
01341       q(2) = a1;
01342       q(3) = a2;
01343       DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
01344       switch (side) {
01345       case base::Left:
01346         for (register int k=mzs; k<=mze; k++)
01347           for (register int j=mys; j<=mye; j++) {
01348             MacroType q1 = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01349             MacroType q2 = MacroVariables(f[b+base::idx(mxe+2,j,k,Nx,Ny)]);
01350             Vector<DataType,3> L;
01351             if (naux==1) { // pressure specified
01352               q(1) = q1(1);  q(2) = q1(2);  q(3) = q1(3);
01353               MacroType dqdn = NormalDerivative(q,q1,q2);
01354               L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01355             }
01356             else if (naux==2) { // normal velocity component specified
01357               q(0) = q1(0);  q(2) = q1(2);  q(3) = q1(3);
01358               MacroType dqdn = NormalDerivative(q,q1,q2);
01359               L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01360             }
01361             else if (naux==4) { // normal and transverse velocity components specified
01362               q(0) = q1(0);
01363               MacroType dqdn = NormalDerivative(q,q1,q2);
01364               L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01365             }
01366             MacroType qn;
01367             qn(0) = q(0) - c1oCsSq2*L(0);
01368             qn(1) = q(1) + c1o2cs*L(0)/q(0);
01369             qn(2) = q(2) - L(1);
01370             qn(3) = q(3) - L(2);
01371             f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(qn);
01372           }
01373         break;
01374       case base::Right:
01375         for (register int k=mzs; k<=mze; k++)
01376           for (register int j=mys; j<=mye; j++) {
01377             MacroType q1 = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01378             MacroType q2 = MacroVariables(f[b+base::idx(mxs-2,j,k,Nx,Ny)]);
01379             if (naux==1) { // pressure specified
01380               q(1) = q1(1);  q(2) = q1(2);  q(3) = q1(3);
01381               MacroType dqdn = NormalDerivative(q,q1,q2);
01382               Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01383               MacroType qn;
01384               qn(0) = q(0) - c1oCsSq2*L(0);
01385               qn(1) = q1(1) + c1o2cs*L(0)/q(0);
01386               qn(2) = q1(2) - L(1);
01387               qn(3) = q1(3) - L(2);
01388             }
01389             else if (naux==2) { // normal velocity component specified
01390               q(0) = q1(0);  q(2) = q1(2);  q(3) = q1(3);
01391               MacroType dqdn = NormalDerivative(q,q1,q2);
01392               Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01393               MacroType qn;
01394               qn(0) = q1(0) - c1oCsSq2*L(0);
01395               qn(1) = q(1) + c1o2cs*L(0)/q(0);
01396               qn(2) = q1(2) - L(1);
01397               qn(3) = q1(3) - L(2);
01398             }
01399             else if (naux==4) { // normal and transverse velocity components specified
01400               q(0) = q1(0);
01401               MacroType dqdn = NormalDerivative(q,q1,q2);
01402               Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01403               MacroType qn;
01404               qn(0) = q1(0) - c1oCsSq2*L(0);
01405               qn(1) = q(1) + c1o2cs*L(0)/q(0);
01406               qn(2) = q(2) - L(1);
01407               qn(3) = q(3) - L(2);
01408             }
01409             f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(qn);
01410           }
01411         break;
01412       case base::Bottom:
01413         for (register int k=mzs; k<=mze; k++)
01414           for (register int i=mxs; i<=mxe; i++) {
01415             MacroType q1 = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01416             MacroType q2 = MacroVariables(f[b+base::idx(i,mye+2,k,Nx,Ny)]);
01417             if (naux==1) { // pressure specified
01418               q(1) = q1(1);  q(2) = q1(2);  q(3) = q1(3);
01419               MacroType dqdn = NormalDerivative(q,q1,q2);
01420               Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01421               MacroType qn;
01422               qn(0) = q(0) - c1oCsSq2*L(0);
01423               qn(1) = q1(1) - L(1);
01424               qn(2) = q1(2) + c1o2cs*L(0)/q(0);
01425               qn(3) = q1(3) - L(2);
01426             }
01427             else if (naux==2) { // normal velocity component specified
01428               q(0) = q1(0);  q(2) = q1(2);  q(3) = q1(3);
01429               MacroType dqdn = NormalDerivative(q,q1,q2);
01430               Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01431               MacroType qn;
01432               qn(0) = q1(0) - c1oCsSq2*L(0);
01433               qn(1) = q(1) - L(1);
01434               qn(2) = q1(2) + c1o2cs*L(0)/q(0);
01435               qn(3) = q1(3) - L(2);
01436             }
01437             else if (naux==4) { // normal and transverse velocity components specified
01438               q(0) = q1(0);
01439               MacroType dqdn = NormalDerivative(q,q1,q2);
01440               Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01441               MacroType qn;
01442               qn(0) = q1(0) - c1oCsSq2*L(0);
01443               qn(1) = q(1) - L(1);
01444               qn(2) = q(2) + c1o2cs*L(0)/q(0);
01445               qn(3) = q(3) - L(2);
01446             }
01447             f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(qn);
01448           }
01449         break;
01450       case base::Top:
01451         for (register int k=mzs; k<=mze; k++)
01452           for (register int i=mxs; i<=mxe; i++) {
01453             MacroType q1 = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01454             MacroType q2 = MacroVariables(f[b+base::idx(i,mys-2,k,Nx,Ny)]);
01455             if (naux==1) { // pressure specified
01456               q(1) = q1(1);  q(2) = q1(2);  q(3) = q1(3);
01457               MacroType dqdn = NormalDerivative(q,q1,q2);
01458               Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01459               MacroType qn;
01460               qn(0) = q(0) - c1oCsSq2*L(0);
01461               qn(1) = q1(1) - L(1);
01462               qn(2) = q1(2) + c1o2cs*L(0)/q(0);
01463               qn(3) = q1(3) - L(2);
01464             }
01465             else if (naux==2) { // normal velocity component specified
01466               q(0) = q1(0);  q(2) = q1(2);  q(3) = q1(3);
01467               MacroType dqdn = NormalDerivative(q,q1,q2);
01468               Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01469               MacroType qn;
01470               qn(0) = q1(0) - c1oCsSq2*L(0);
01471               qn(1) = q(1) - L(1);
01472               qn(2) = q1(2) + c1o2cs*L(0)/q(0);
01473               qn(3) = q1(3) - L(2);
01474             }
01475             else if (naux==4) { // normal and transverse velocity components specified
01476               q(0) = q1(0);
01477               MacroType dqdn = NormalDerivative(q,q1,q2);
01478               Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01479               MacroType qn;
01480               qn(0) = q1(0) - c1oCsSq2*L(0);
01481               qn(1) = q(1) - L(1);
01482               qn(2) = q(2) + c1o2cs*L(0)/q(0);
01483               qn(3) = q(3) - L(2);
01484             }
01485             f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(qn);
01486           }
01487         break;
01488       case base::Front:
01489         for (register int j=mys; j<=mye; j++)
01490           for (register int i=mxs; i<=mxe; i++) {
01491             MacroType q1 = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01492             MacroType q2 = MacroVariables(f[b+base::idx(i,j,mze+2,Nx,Ny)]);
01493             if (naux==1) { // pressure specified
01494               q(1) = q1(1);  q(2) = q1(2);  q(3) = q1(3);
01495               MacroType dqdn = NormalDerivative(q,q1,q2);
01496               Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01497               MacroType qn;
01498               qn(0) = q(0) - c1oCsSq2*L(0);
01499               qn(1) = q(1) - L(1);
01500               qn(2) = q(2) - L(2);
01501               qn(3) = q(3) + c1o2cs*L(0)/q(0);
01502             }
01503             else if (naux==2) { // normal velocity component specified
01504               q(0) = q1(0);  q(2) = q1(2);  q(3) = q1(3);
01505               MacroType dqdn = NormalDerivative(q,q1,q2);
01506               Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01507               MacroType qn;
01508               qn(0) = q(0) - c1oCsSq2*L(0);
01509               qn(1) = q(1) - L(1);
01510               qn(2) = q(2) - L(2);
01511               qn(3) = q(3) + c1o2cs*L(0)/q(0);
01512             }
01513             else if (naux==4) { // normal and transverse velocity components specified
01514               q(0) = q1(0);
01515               MacroType dqdn = NormalDerivative(q,q1,q2);
01516               Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01517               MacroType qn;
01518               qn(0) = q(0) - c1oCsSq2*L(0);
01519               qn(1) = q(1) - L(1);
01520               qn(2) = q(2) - L(2);
01521               qn(3) = q(3) + c1o2cs*L(0)/q(0);
01522             }
01523             MacroType dqdn = NormalDerivative(q,q1,q2);
01524             Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01525 
01526             // Solve LODI equations :
01527             MacroType qn;
01528             qn(0) = q(0) - c1oCsSq2*L(0);
01529             qn(1) = q(1) - L(1);
01530             qn(2) = q(2) - L(2);
01531             qn(3) = q(3) + c1o2cs*L(0)/q(0);
01532 
01533             f[b+base::idx(i,j,mze,Nx,Ny)] = Equilibrium(qn);
01534           }
01535         break;
01536       case base::Back:
01537         for (register int j=mys; j<=mye; j++)
01538           for (register int i=mxs; i<=mxe; i++) {
01539             MacroType q1 = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01540             MacroType q2 = MacroVariables(f[b+base::idx(i,j,mzs-2,Nx,Ny)]);
01541             if (naux==1) { // pressure specified
01542               q(1) = q1(1);  q(2) = q1(2);  q(3) = q1(3);
01543               MacroType dqdn = NormalDerivative(q,q1,q2);
01544               Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01545               MacroType qn;
01546               qn(0) = q(0) - c1oCsSq2*L(0);
01547               qn(1) = q(1) - L(1);
01548               qn(2) = q(2) - L(2);
01549               qn(3) = q(3) + c1o2cs*L(0)/q(0);
01550             }
01551             else if (naux==2) { // normal velocity component specified
01552               q(0) = q1(0);  q(2) = q1(2);  q(3) = q1(3);
01553               MacroType dqdn = NormalDerivative(q,q1,q2);
01554               Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01555               MacroType qn;
01556               qn(0) = q(0) - c1oCsSq2*L(0);
01557               qn(1) = q(1) - L(1);
01558               qn(2) = q(2) - L(2);
01559               qn(3) = q(3) + c1o2cs*L(0)/q(0);
01560             }
01561             else if (naux==4) { // normal and transverse velocity components specified
01562               q(0) = q1(0);
01563               MacroType dqdn = NormalDerivative(q,q1,q2);
01564               Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01565               MacroType qn;
01566               qn(0) = q(0) - c1oCsSq2*L(0);
01567               qn(1) = q(1) - L(1);
01568               qn(2) = q(2) - L(2);
01569               qn(3) = q(3) + c1o2cs*L(0)/q(0);
01570             }
01571             f[b+base::idx(i,j,mzs,Nx,Ny)] = Equilibrium(qn);
01572           }
01573         break;
01574       }
01575     }
01576 
01577     else if (type==NoSlipWallEntranceExit) {
01578       int lc_indx[3] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1), mzs*fvec.stepsize(2) };
01579       DataType WLB[3],  WUB[3];
01580       base::GH().wlb(WLB);
01581       base::GH().wub(WUB);
01582       DCoords wc;
01583       DataType slip = 0.;
01584       DataType lxs, lxe, lys, lye, lzs, lze;
01585       DataType min[3], max[3];
01586       if (naux!=6) assert(0);
01587       lxs = aux[0] - WLB[0];
01588       lxe = WUB[0] - aux[1];
01589       lys = aux[2] - WLB[1];
01590       lye = WUB[1] - aux[3];
01591       lzs = aux[4] - WLB[2];
01592       lze = WUB[2] - aux[5];
01593       min[0] = aux[0]; max[0] = aux[1];
01594       min[1] = aux[2]; max[1] = aux[3];
01595       min[2] = aux[4]; max[2] = aux[5];
01596 
01597       switch (side) {
01598       case base::Left:
01599         lc_indx[0] = mxe*fvec.stepsize(0);
01600         for (register int k=mzs; k<=mze; k++) {
01601           lc_indx[2] = k*fvec.stepsize(2);
01602           for (register int j=mys; j<=mye; j++) {
01603             lc_indx[1] = j*fvec.stepsize(1);
01604             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01605             if (wc(1)>=min[1] && wc(1)<=max[1]) {
01606               if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
01607               if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
01608             }
01609             else {
01610               if (wc(1)<WLB[1] || wc(1) >WUB[1])
01611                 slip = DataType(0.);
01612               else if (wc(1)<min[1])
01613                 slip = (wc(1) - WLB[1])/lys;
01614               else if (wc(1)>max[1])
01615                 slip = (WUB[1] - wc(1))/lye;
01616               if (j>mys && j<mye) {
01617                 f[b+base::idx(mxe,j,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10)
01618                     + slip*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
01619                 f[b+base::idx(mxe,j,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8)
01620                     + slip*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
01621               }
01622               else if (j==mys) {
01623                 f[b+base::idx(mxe,j,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k,Nx,Ny)](10)
01624                     + slip*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
01625                 f[b+base::idx(mxe,j,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8)
01626                     + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);
01627               }
01628               else if (j==mye) {
01629                 f[b+base::idx(mxe,j,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10)
01630                     + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
01631                 f[b+base::idx(mxe,j,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k,Nx,Ny)](8)
01632                     + slip*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
01633               }
01634             }
01635             if (wc(2)>=min[2] && wc(2)<=max[2]) {
01636               if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01637               if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01638             }
01639             else {
01640               if (wc(2)<WLB[2] || wc(2) >WUB[2])
01641                 slip = DataType(0.);
01642               else if (wc(2)<min[2])
01643                 slip = (wc(2) - WLB[2])/lzs;
01644               else if (wc(2)>max[2])
01645                 slip = (WUB[2] - wc(2))/lze;
01646               if (k>mzs && k<mze) {
01647                 f[b+base::idx(mxe,j,k,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18)
01648                     + slip*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01649                 f[b+base::idx(mxe,j,k,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16)
01650                     + slip*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01651               }
01652               else if (k==mzs) {
01653                 f[b+base::idx(mxe,j,k,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k,Nx,Ny)](18)
01654                     + slip*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01655                 f[b+base::idx(mxe,j,k,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16)
01656                     + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
01657               }
01658               else if (k==mze) {
01659                 f[b+base::idx(mxe,j,k,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18)
01660                     + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
01661                 f[b+base::idx(mxe,j,k,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k,Nx,Ny)](16)
01662                     + slip*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01663               }
01664             }
01665             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
01666           }
01667         }
01668         break;
01669       case base::Right:
01670         lc_indx[0] = mxs*fvec.stepsize(0);
01671         for (register int k=mzs; k<=mze; k++) {
01672           lc_indx[2] = k*fvec.stepsize(2);
01673           for (register int j=mys; j<=mye; j++) {
01674             lc_indx[1] = j*fvec.stepsize(1);
01675             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01676             if (wc(1)>=min[1] && wc(1)<=max[1]) {
01677               if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
01678               if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
01679             }
01680 
01681             else {
01682               if (wc(1)<WLB[1] || wc(1) >WUB[1])
01683                 slip = DataType(0.);
01684               else if (wc(1)<min[1])
01685                 slip = (wc(1) - WLB[1])/lys;
01686               else if (wc(1)>max[1])
01687                 slip = (WUB[1] - wc(1))/lye;
01688               if (j>mys && j<mye) {
01689                 f[b+base::idx(mxs,j,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7)
01690                     + slip*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
01691                 f[b+base::idx(mxs,j,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9)
01692                     + slip*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
01693               }
01694               else if (j==mys) {
01695                 f[b+base::idx(mxs,j,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k,Nx,Ny)](7)
01696                     + slip*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
01697                 f[b+base::idx(mxs,j,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9)
01698                     + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
01699               }
01700               else if (j==mye) {
01701                 f[b+base::idx(mxs,j,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7)
01702                     + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
01703                 f[b+base::idx(mxs,j,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k,Nx,Ny)](9)
01704                     + slip*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
01705               }
01706             }
01707             if (wc(2)>=min[2] && wc(2)<=max[2]) {
01708               if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
01709               if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
01710             }
01711             else {
01712               if (wc(2)<WLB[2] || wc(2) >WUB[2])
01713                 slip = DataType(0.);
01714               else if (wc(2)<min[2])
01715                 slip = (wc(2) - WLB[2])/lzs;
01716               else if (wc(2)>max[2])
01717                 slip = (WUB[2] - wc(2))/lze;
01718               if (k>mzs && k<mze) {
01719                 f[b+base::idx(mxs,j,k,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15)
01720                     + slip*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
01721                 f[b+base::idx(mxs,j,k,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17)
01722                     + slip*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
01723               }
01724               else if (k==mzs) {
01725                 f[b+base::idx(mxs,j,k,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k,Nx,Ny)](15)
01726                     + slip*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
01727                 f[b+base::idx(mxs,j,k,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17)
01728                     + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
01729               }
01730               else if (k==mze) {
01731                 f[b+base::idx(mxs,j,k,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15)
01732                     + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
01733                 f[b+base::idx(mxs,j,k,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k,Nx,Ny)](17)
01734                     + slip*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
01735               }
01736             }
01737             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
01738           }
01739         }
01740         break;
01741       case base::Bottom:
01742         lc_indx[1] = mye*fvec.stepsize(1);
01743         for (register int k=mzs; k<=mze; k++) {
01744           lc_indx[2] = k*fvec.stepsize(2);
01745           for (register int i=mxs; i<=mxe; i++) {
01746             lc_indx[0] = i*fvec.stepsize(0);
01747             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01748             if (wc(0)>=min[0] && wc(0)<=max[0]) {
01749               if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
01750               if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
01751             }
01752             else {
01753               if (wc(0)<WLB[0] || wc(0)>WUB[0])
01754                 slip = DataType(0.);
01755               else if (wc(0)<min[0])
01756                 slip = (wc(0) - WLB[0])/lxs;
01757               else if (wc(0)>max[0])
01758                 slip = (WUB[0] - wc(0))/lxe;
01759               if (i>mxs && i<mxe) {
01760                 f[b+base::idx(i,mye,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9)
01761                     + slip*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
01762                 f[b+base::idx(i,mye,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8)
01763                     + slip*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
01764               }
01765               else if (i==mxs) {
01766                 f[b+base::idx(i,mye,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k,Nx,Ny)](9)
01767                     + slip*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
01768                 f[b+base::idx(i,mye,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8)
01769                     + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
01770               }
01771               else if (i==mxe) {
01772                 f[b+base::idx(i,mye,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9)
01773                     + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
01774                 f[b+base::idx(i,mye,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k,Nx,Ny)](8)
01775                     + slip*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
01776               }
01777             }
01778             if (wc(2)>=min[2] && wc(2)<=max[2]) {
01779               if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
01780               if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
01781             }
01782             else {
01783               if (wc(2)<WLB[2] || wc(2) >WUB[2])
01784                 slip = DataType(0.);
01785               else if (wc(2)<min[2])
01786                 slip = (wc(2) - WLB[2])/lzs;
01787               else if (wc(2)>max[2])
01788                 slip = (WUB[2] - wc(2))/lze;
01789               if (k>mzs && k<mze) {
01790                 f[b+base::idx(i,mye,k,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14)
01791                     + slip*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
01792                 f[b+base::idx(i,mye,k,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12)
01793                     + slip*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
01794               }
01795               else if (k==mzs) {
01796                 f[b+base::idx(i,mye,k,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k,Nx,Ny)](14)
01797                     + slip*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
01798                 f[b+base::idx(i,mye,k,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12)
01799                     + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](14);
01800               }
01801               else if (k==mze) {
01802                 f[b+base::idx(i,mye,k,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14)
01803                     + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](12);
01804                 f[b+base::idx(i,mye,k,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k,Nx,Ny)](12)
01805                     + slip*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
01806               }
01807             }
01808             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
01809           }
01810         }
01811         break;
01812       case base::Top:
01813         lc_indx[1] = mys*fvec.stepsize(1);
01814         for (register int k=mzs; k<=mze; k++) {
01815           lc_indx[2] = k*fvec.stepsize(2);
01816           for (register int i=mxs; i<=mxe; i++) {
01817             lc_indx[0] = i*fvec.stepsize(0);
01818             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01819             if (wc(0)>=min[0] && wc(0)<=max[0]) {
01820               if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
01821               if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
01822             }
01823             else {
01824               if (wc(0)<WLB[0] || wc(0) >WUB[0])
01825                 slip = DataType(0.);
01826               else if (wc(0)<min[0])
01827                 slip = (wc(0) - WLB[0])/lxs;
01828               else if (wc(0)>max[0])
01829                 slip = (WUB[0] - wc(0))/lxe;
01830               if (i>mxs && i<mxe) {
01831                 f[b+base::idx(i,mys,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7)
01832                     + slip*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
01833                 f[b+base::idx(i,mys,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10)
01834                     + slip*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
01835               }
01836               else if (i==mxs) {
01837                 f[b+base::idx(i,mys,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k,Nx,Ny)](7)
01838                     + slip*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
01839                 f[b+base::idx(i,mys,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10)
01840                     + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
01841               }
01842               else if (i==mxe) {
01843                 f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7)
01844                     + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
01845                 f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10)
01846                     + slip*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
01847               }
01848             }
01849             if (wc(2)>=min[2] && wc(2)<=max[2]) {
01850               if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
01851               if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
01852             }
01853             else {
01854               if (wc(2)<WLB[2] || wc(2)>WUB[2])
01855                 slip = DataType(0.);
01856               else if (wc(2)<min[2])
01857                 slip = (wc(2) - WLB[2])/lzs;
01858               else if (wc(2)>max[2])
01859                 slip = (WUB[2] - wc(2))/lze;
01860               if (k>mzs && k<mze) {
01861                 f[b+base::idx(i,mys,k,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11)
01862                     + slip*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
01863                 f[b+base::idx(i,mys,k,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13)
01864                     + slip*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
01865               }
01866               else if (k==mzs) {
01867                 f[b+base::idx(i,mys,k,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k,Nx,Ny)](11)
01868                     + slip*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
01869                 f[b+base::idx(i,mys,k,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13)
01870                     + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
01871               }
01872               else if (k==mze) {
01873                 f[b+base::idx(i,mys,k,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11)
01874                     + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
01875                 f[b+base::idx(i,mys,k,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k,Nx,Ny)](13)
01876                     + slip*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
01877               }
01878             }
01879             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
01880           }
01881         }
01882         break;
01883       case base::Front:
01884         lc_indx[2] = mze*fvec.stepsize(2);
01885         for (register int j=mys; j<=mye; j++) {
01886           lc_indx[1] = j*fvec.stepsize(1);
01887           for (register int i=mxs; i<=mxe; i++) {
01888             lc_indx[0] = i*fvec.stepsize(0);
01889             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01890             if (wc(0)>=min[0] && wc(0)<=max[0]) {
01891               if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
01892               if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
01893             }
01894             else {
01895               if (wc(0)<WLB[0] || wc(0) >WUB[0])
01896                 slip = DataType(0.);
01897               else if (wc(0)<min[0])
01898                 slip = (wc(0) - WLB[0])/lxs;
01899               else if (wc(0)>max[0])
01900                 slip = (WUB[0] - wc(0))/lxe;
01901               if (i>mxs && i<mxe) {
01902                 f[b+base::idx(i,j,mze,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17)
01903                     + slip*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
01904                 f[b+base::idx(i,j,mze,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16)
01905                     + slip*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
01906               }
01907               else if (i==mxs) {
01908                 f[b+base::idx(i,j,mze,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mze+1,Nx,Ny)](17)
01909                     + slip*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
01910                 f[b+base::idx(i,j,mze,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16)
01911                     + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
01912               }
01913               else if (i==mxe) {
01914                 f[b+base::idx(i,j,mze,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17)
01915                     + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
01916                 f[b+base::idx(i,j,mze,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mze+1,Nx,Ny)](16)
01917                     + slip*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
01918               }
01919             }
01920             if (wc(1)>=min[1] && wc(1)<=max[1]) {
01921               if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
01922               if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
01923             }
01924             else {
01925               if (wc(1)<WLB[1] || wc(1) >WUB[1])
01926                 slip = DataType(0.);
01927               else if (wc(1)<min[1])
01928                 slip = (wc(1) - WLB[1])/lys;
01929               else if (wc(1)>max[1])
01930                 slip = (WUB[1] - wc(1))/lye;
01931               if (j>mys && j<mye) {
01932                 f[b+base::idx(i,j,mze,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13)
01933                     + slip*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
01934                 f[b+base::idx(i,j,mze,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12)
01935                     + slip*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
01936               }
01937               else if (j==mys) {
01938                 f[b+base::idx(i,j,mze,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mze+1,Nx,Ny)](13)
01939                     + slip*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
01940                 f[b+base::idx(i,j,mze,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12)
01941                     + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
01942               }
01943               else if (j==mye) {
01944                 f[b+base::idx(i,j,mze,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13)
01945                     + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
01946                 f[b+base::idx(i,j,mze,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mze+1,Nx,Ny)](12)
01947                     + slip*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
01948               }
01949             }
01950             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
01951           }
01952         }
01953         break;
01954       case base::Back:
01955         lc_indx[2] = mzs*fvec.stepsize(2);
01956         for (register int j=mys; j<=mye; j++) {
01957           lc_indx[1] = j*fvec.stepsize(1);
01958           for (register int i=mxs; i<=mxe; i++) {
01959             lc_indx[0] = i*fvec.stepsize(0);
01960             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01961             if (wc(0)>=min[0] && wc(0)<=max[0]) {
01962               if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
01963               if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
01964             }
01965             else {
01966               if (wc(0)<WLB[0] || wc(0) >WUB[0])
01967                 slip = DataType(0.);
01968               else if (wc(0)<min[0])
01969                 slip = (wc(0) - WLB[0])/lxs;
01970               else if (wc(0)>max[0])
01971                 slip = (WUB[0] - wc(0))/lxe;
01972               if (i>mxs && i<mxe) {
01973                 f[b+base::idx(i,j,mzs,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15)
01974                     + slip*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
01975                 f[b+base::idx(i,j,mzs,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18)
01976                     + slip*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
01977               }
01978               else if (i==mxs) {
01979                 f[b+base::idx(i,j,mzs,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mzs-1,Nx,Ny)](15)
01980                     + slip*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
01981                 f[b+base::idx(i,j,mzs,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18)
01982                     + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
01983               }
01984               else if (i==mxe) {
01985                 f[b+base::idx(i,j,mzs,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15)
01986                     + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
01987                 f[b+base::idx(i,j,mzs,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mzs-1,Nx,Ny)](18)
01988                     + slip*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
01989               }
01990             }
01991             if (wc(1)>=min[1] && wc(1)<=max[1]) {
01992               if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
01993               if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
01994             }
01995             else {
01996               if (wc(1)<WLB[1] || wc(1) >WUB[1])
01997                 slip = DataType(0.);
01998               else if (wc(1)<min[1])
01999                 slip = (wc(1) - WLB[1])/lys;
02000               else if (wc(1)>max[1])
02001                 slip = (WUB[1] - wc(1))/lye;
02002               if (j>mys && j<mye) {
02003                 f[b+base::idx(i,j,mzs,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11)
02004                     + slip*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
02005                 f[b+base::idx(i,j,mzs,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14)
02006                     + slip*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
02007               }
02008               else if (j==mys) {
02009                 f[b+base::idx(i,j,mzs,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mzs-1,Nx,Ny)](11)
02010                     + slip*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
02011                 f[b+base::idx(i,j,mzs,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14)
02012                     + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
02013               }
02014               else if (j==mye) {
02015                 f[b+base::idx(i,j,mzs,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11)
02016                     + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
02017                 f[b+base::idx(i,j,mzs,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mzs-1,Nx,Ny)](14)
02018                     + slip*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
02019               }
02020             }
02021             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
02022           }
02023         }
02024         break;
02025       }
02026     }
02027   }
02028 
02029   inline virtual const MacroType NormalDerivative(const MacroType qa,const MacroType qb, const MacroType qc) const
02030   { return DataType(-1.5)*qa + DataType(2.)*qb - DataType(0.5)*qc; }
02031   inline virtual Vector<DataType,3> WaveAmplitudes(const DataType rho, const DataType vn, const DataType drhodn,
02032                                                    const DataType dvndn, const DataType dvt0dn, const DataType dvt1dn) const {
02033     Vector<DataType,3> L;
02034     L(0) = (vn - LatticeSpeedOfSound())*(cs2*drhodn - rho*LatticeSpeedOfSound()*dvndn);
02035     L(1) = vn*dvt0dn;
02036     L(2) = vn*dvt1dn;
02037     return L;
02038   }
02039 
02040   virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
02041                              const MicroType* f, const point_type* xc, const DataType* distance,
02042                              const point_type* normal, DataType* aux=0, const int naux=0,
02043                              const int scaling=0) const {
02044     if (type==GFMNoSlipWall || type==GFMSlipWall) {
02045       DataType fac = S0;
02046       if (scaling==base::Physical)
02047         fac /= U0;
02048       for (int n=0; n<nc; n++) {
02049         MacroType q=MacroVariables(f[n]);
02050         DataType u=-q(1);
02051         DataType v=-q(2);
02052         DataType w=-q(3);
02053 
02054         // Add boundary velocities if available
02055         if (naux>=3) {
02056           u += fac*aux[naux*n];
02057           v += fac*aux[naux*n+1];
02058           w += fac*aux[naux*n+2];
02059         }
02060         if (type==GFMNoSlipWall) {
02061           // Prescribe entire velocity vector in ghost cells
02062           q(1) += DataType(2.)*u;
02063           q(2) += DataType(2.)*v;
02064           q(3) += DataType(2.)*w;
02065         }
02066         else if (type==GFMSlipWall) {
02067           // Prescribe only normal velocity vector in ghost cells
02068           DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1)+w*normal[n](2));
02069           q(1) += vl*normal[n](0);
02070           q(2) += vl*normal[n](1);
02071           q(3) += vl*normal[n](2);
02072         }
02073         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
02074       }
02075     }
02076 
02077     else if (type==GFMExtrapolation) {
02078       for (int n=0; n<nc; n++) {
02079         MacroType q=MacroVariables(f[n]);
02080         DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1)+q(3)*normal[n](2);
02081         q(1) = vl*normal[n](0);
02082         q(2) = vl*normal[n](1);
02083         q(3) = vl*normal[n](2);
02084         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
02085       }
02086     }
02087 
02088     else if (type==GFMWallLaw ) {
02089       DCoords dx = base::GH().worldStep(fvec.stepsize());
02090       DataType fac = S0;
02091       if (scaling==base::Physical)
02092         fac /= U0;
02093 
02094       // DataType BLH = 10.;
02095       DataType yp_ = 0.,  yp = 0., yMax = 50.;
02096       DataType yPlus[3];
02097       int lc_low[3], lc_up[3], yPlus_lc[3];
02098 
02099       // int mxe = fvec.upper(0)/fvec.stepsize(0);
02100       // int mye = fvec.upper(1)/fvec.stepsize(1);
02101       // int mze = fvec.upper(2)/fvec.stepsize(2);
02102 
02103       DataType shear_vel = 1., karman = 0.41, y0 = 1., law = 1.;
02104 
02105       DataType shearStress;
02106       DataType u, v, w, u_, v_, w_, vel0;
02107       // DataType vel1;
02108       MacroType q, q_;
02109       int kMax = 0;
02110 
02111       for (int n=0; n<nc; n++) {
02112 
02114         q=MacroVariables(f[n]);
02115         q_=MacroVariables(f[n]);
02116 
02117         yPlus[0] = (xc[n](0) - DataType(2.)*dx(0)*normal[n](0));
02118         yPlus[1] = (xc[n](1) - DataType(2.)*dx(1)*normal[n](1));
02119         yPlus[2] = (xc[n](2) - DataType(2.)*dx(2)*normal[n](2));
02120 
02121         yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
02122         yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
02123         yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);
02124 
02125         lc_low[0] = fvec.lower(0);   lc_low[1] = fvec.lower(1);   lc_low[2] = fvec.lower(2);
02126         lc_up[0] = fvec.upper(0);   lc_up[1] = fvec.upper(1);   lc_up[2] = fvec.upper(2);
02127 
02128         //std::cout << "xc " << xc[n] << " yPlus " << yPlus[0] << " " << yPlus[1] << " " << yPlus[2] << std::endl;
02129 
02130         if ( yPlus_lc[0] >  lc_low[0] ) {
02131           if ( yPlus_lc[1] >  lc_low[1] ) {
02132             if ( yPlus_lc[2] >  lc_low[2] ) {
02133 
02134               if ( yPlus_lc[0] < lc_up[0] ) {
02135                 if ( yPlus_lc[1] <  lc_up[1] ) {
02136                   if ( yPlus_lc[2] <  lc_up[2] ) {
02137                     q=MacroVariables( fvec( base::GH().localCoords((const DataType *) yPlus) ) );
02138                   }
02139                 }
02140               }
02141 
02142             }
02143           }
02144         }
02145 
02146         u_=-q(1);
02147         v_=-q(2);
02148         w_=-q(3);
02149         // vel1 = sqrt(u_*u_+v_*v_+w_*w_);
02150 
02151         q=MacroVariables(f[n]);
02152         u=-q(1);
02153         v=-q(2);
02154         w=-q(3);
02155         vel0 = sqrt(u*u+v*v+w*w);
02156 
02157         u_+=q(1);
02158         v_+=q(2);
02159         w_+=q(3);
02160         vel0 = ( u_*normal[n](0) + v_*normal[n](1) + w_*normal[n](2) );
02161 
02162         u_ =  u_ - vel0 * normal[n](0);
02163         v_ =  v_ - vel0 * normal[n](1);
02164         w_ =  w_ - vel0 * normal[n](2);
02165         vel0 = sqrt(u_*u_ + v_*v_ + w_*w_)*U0/S0;
02166 
02167         shearStress = (DataType(0.5)*vel0/dx(0) * q(0)*R0*GasViscosity() );
02168 
02169         // Add boundary velocities if available
02170         if (naux>=3) {
02171           u += fac*aux[naux*n];
02172           v += fac*aux[naux*n+1];
02173           w += fac*aux[naux*n+2];
02174         }
02175         q(1) += DataType(2.)*u;
02176         q(2) += DataType(2.)*v;
02177         q(3) += DataType(2.)*w;
02178         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
02179 
02181         shear_vel = std::sqrt( shearStress/(q(0)*R0) );
02182 
02184         y0 = karman/30.;
02185 
02186         //std::cout << n << " vel0 " << vel0 << " vel1 " << vel1  << " shear_vel = " << shear_vel << std::endl;
02187 
02188         if ( shear_vel > 1.e-5 ) {
02189           kMax = (int) std::ceil( (yMax*GasViscosity())/(dx(0)*shear_vel) );
02190 
02191           for (int k=1; k < kMax; k++) {
02192 
02193             yPlus[0] = (xc[n](0) - k*dx(0)*normal[n](0));
02194             yPlus[1] = (xc[n](1) - k*dx(1)*normal[n](1));
02195             yPlus[2] = (xc[n](2) - k*dx(2)*normal[n](2));
02196 
02197             yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
02198             yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
02199             yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);
02200 
02201             if ( yPlus_lc[0] >  lc_low[0] ) {
02202               if ( yPlus_lc[1] >  lc_low[1] ) {
02203                 if ( yPlus_lc[2] >  lc_low[2] ) {
02204 
02205                   if ( yPlus_lc[0] < lc_up[0] ) {
02206                     if ( yPlus_lc[1] <  lc_up[1] ) {
02207                       if ( yPlus_lc[2] <  lc_up[2] ) {
02208                         yp = k*dx(0);
02209                         yp_ = yp*shear_vel/GasViscosity();
02210 
02213                         if ( yp_ <= 5.0 ) {
02214                           law = yp_;
02215                         }
02217                         else if (yp_ > 5 && yp_ < 30.0 ) {
02218                           //law = (shear_vel/karman*std::log( yp/y0 ) + yp_)*DataType(0.5);
02219                           law = std::min( yp_ , shear_vel/karman*std::log( yp/y0 ) );
02220 
02221                         }
02223                         else if (yp_ >= 30.0 ) {
02224                           law = shear_vel/karman*std::log( yp/y0 );
02225                         }
02226                         //std::cout << "k = " << k << " yp = " << yp << " shear_vel = " << shear_vel << " yp_ = " << yp_ << " law = " << law << std::endl;
02227                         q_=MacroVariables(f[n]);
02228                         u = law*(-q_(1) + fac*aux[naux*n]);
02229                         v = law*(-q_(2) + fac*aux[naux*n+1]);
02230                         w = law*(-q_(3) + fac*aux[naux*n+2]);
02231 
02232                         q_(1) += DataType(2.)*u;
02233                         q_(2) += DataType(2.)*v;
02234                         q_(3) += DataType(2.)*w;
02235                         fvec( base::GH().localCoords((const DataType *) yPlus) ) = Equilibrium(q_);
02236                       }
02237                     }
02238                   }
02239 
02240 
02241                 }
02242               }
02243             }
02244 
02245           }
02246 
02247         }
02248 
02249       }
02250       //std::cout << "====== WallLaw Done\n";
02251     }
02252 
02253   }
02254 
02255   virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
02256                       const int skip_ghosts=1) const {
02257     int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
02258     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
02259     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
02260     MicroType *f = (MicroType *)fvec.databuffer();
02261     DataType *work = (DataType *)workvec.databuffer();
02262     DCoords dx = base::GH().worldStep(fvec.stepsize());
02263     DataType dt = dx(0)*S0/U0;
02264 
02265     if ((cnt>=1 && cnt<=10) || (cnt>=19 && cnt<=24) || (cnt>=28 && cnt<=33)) { 
02266       for (register int k=mzs; k<=mze; k++)
02267         for (register int j=mys; j<=mye; j++)
02268           for (register int i=mxs; i<=mxe; i++) {
02269             MacroType q=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
02270             switch(cnt) {
02271             case 1:
02272               if (method[0]==0) work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0+rhop;
02273               else 
02274                 work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0;
02275               break;
02276             case 2:
02277               work[base::idx(i,j,k,Nx,Ny)]=q(1)*U0/S0;
02278               break;
02279             case 3:
02280               work[base::idx(i,j,k,Nx,Ny)]=q(2)*U0/S0;
02281               break;
02282             case 4:
02283               work[base::idx(i,j,k,Nx,Ny)]=q(3)*U0/S0;
02284               break;
02285             case 6:
02286               work[base::idx(i,j,k,Nx,Ny)]=TempEquation(U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure()); 
02287               break;
02288             case 7:
02289               work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure(); 
02290               break;
02291             case 9: {
02292               MicroType feq = Equilibrium(q);     
02293               work[base::idx(i,j,k,Nx,Ny)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt),
02294                                                         GasSpeedofSound(),dt);
02295               break;
02296             }
02297             case 10:
02298               work[base::idx(i,j,k,Nx,Ny)]=std::sqrt(q(1)*q(1)+q(2)*q(2)+q(3)*q(3))*U0/S0;
02299               break;
02300             case 19:
02301             case 20:
02302             case 21:
02303             case 22:
02304             case 23:
02305             case 24: {
02306               MicroType feq = Equilibrium(q);
02307               DataType omega;
02308               if (turbulence == laminar) 
02309                 omega = Omega(dt);
02310               else if (turbulence == LES_Smagorinsky) 
02311                 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
02312               else if (turbulence == LES_dynamic) 
02313                 omega = Omega_LES_dynamic(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
02314               work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*DeviatoricStress(f[base::idx(i,j,k,Nx,Ny)],feq,omega)(cnt-19);
02315               break;
02316             }
02317             case 28: 
02318             case 29:
02319             case 30: 
02320             case 31: 
02321             case 32: 
02322             case 33: {
02323               DataType omega;
02324               if (turbulence == laminar) 
02325                 omega = Omega(dt);
02326               else if (turbulence == LES_Smagorinsky) {
02327                 MicroType feq = Equilibrium(q);         
02328                 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
02329               }
02330               else if (turbulence == LES_dynamic) {
02331                 MicroType feq = Equilibrium(q);         
02332                 omega = Omega_LES_dynamic(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
02333               }
02334               work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*Stress(f[base::idx(i,j,k,Nx,Ny)],q,omega)(cnt-28);
02335               if (cnt==28 || cnt==30 || cnt==33) work[base::idx(i,j,k,Nx,Ny)]-=BasePressure();
02336               break;
02337             }
02338             }
02339           }
02340     }
02341     else if ((cnt>=11 && cnt<=14) || cnt==18 || (cnt>=40 && cnt<=45)) {
02342       macro_grid_data_type qgrid(fvec.bbox());
02343       MacroType *q = (MacroType *)qgrid.databuffer();
02344       for (register int k=mzs; k<=mze; k++)
02345         for (register int j=mys; j<=mye; j++)
02346           for (register int i=mxs; i<=mxe; i++) {
02347             q[base::idx(i,j,k,Nx,Ny)]=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
02348             q[base::idx(i,j,k,Nx,Ny)](0)*=R0;
02349             q[base::idx(i,j,k,Nx,Ny)](1)*=U0/S0;
02350             q[base::idx(i,j,k,Nx,Ny)](2)*=U0/S0;
02351             q[base::idx(i,j,k,Nx,Ny)](3)*=U0/S0;
02352           }
02353       
02354       if (skip_ghosts==0) {
02355         switch(cnt) {
02356         case 40:
02357           mxs+=1; mxe-=1; 
02358           break;
02359         case 42:
02360           mys+=1; mye-=1; 
02361           break;
02362         case 45:
02363           mzs+=1; mze-=1; 
02364           break;
02365         case 11:
02366         case 44:
02367           mys+=1; mye-=1; mzs+=1; mze-=1;
02368           break;
02369         case 12:
02370         case 43:
02371           mxs+=1; mxe-=1; mzs+=1; mze-=1;
02372           break;
02373         case 13:
02374         case 41:
02375           mxs+=1; mxe-=1; mys+=1; mye-=1; 
02376           break;
02377         case 14:
02378         case 18:
02379           mxs+=1; mxe-=1; mys+=1; mye-=1; mzs+=1; mze-=1;
02380           break;
02381         }
02382       }
02383 
02384       DataType dxux, dxuy, dxuz, dyux, dyuy, dyuz, dzux, dzuy, dzuz; 
02385       for (register int k=mzs; k<=mze; k++)
02386         for (register int j=mys; j<=mye; j++)
02387           for (register int i=mxs; i<=mxe; i++) {
02388             if (cnt==18 || cnt==40)
02389               dxux=(q[base::idx(i+1,j,k,Nx,Ny)](1)-q[base::idx(i-1,j,k,Nx,Ny)](1))/(2.*dx(0));
02390             if (cnt==18 || cnt==42)
02391               dyuy=(q[base::idx(i,j+1,k,Nx,Ny)](2)-q[base::idx(i,j-1,k,Nx,Ny)](2))/(2.*dx(1));
02392             if (cnt==18 || cnt==45)
02393               dzuz=(q[base::idx(i,j,k+1,Nx,Ny)](3)-q[base::idx(i,j,k-1,Nx,Ny)](3))/(2.*dx(2));
02394             if (cnt==11 || cnt==14 || cnt==18 || cnt==44) {
02395               dyuz=(q[base::idx(i,j+1,k,Nx,Ny)](3)-q[base::idx(i,j-1,k,Nx,Ny)](3))/(2.*dx(1));
02396               dzuy=(q[base::idx(i,j,k+1,Nx,Ny)](2)-q[base::idx(i,j,k-1,Nx,Ny)](2))/(2.*dx(2));
02397             }
02398             if (cnt==12 || cnt==14 || cnt==18 || cnt==43) {
02399               dxuz=(q[base::idx(i+1,j,k,Nx,Ny)](3)-q[base::idx(i-1,j,k,Nx,Ny)](3))/(2.*dx(0));
02400               dzux=(q[base::idx(i,j,k+1,Nx,Ny)](1)-q[base::idx(i,j,k-1,Nx,Ny)](1))/(2.*dx(2));
02401             }
02402             if (cnt==13 || cnt==14 || cnt==18 || cnt==41) {         
02403               dxuy=(q[base::idx(i+1,j,k,Nx,Ny)](2)-q[base::idx(i-1,j,k,Nx,Ny)](2))/(2.*dx(0));
02404               dyux=(q[base::idx(i,j+1,k,Nx,Ny)](1)-q[base::idx(i,j-1,k,Nx,Ny)](1))/(2.*dx(1));
02405             }
02406             switch(cnt) {
02407             case 11: 
02408               work[base::idx(i,j,k,Nx,Ny)]=dyuz-dzuy;
02409               break;
02410             case 12: 
02411               work[base::idx(i,j,k,Nx,Ny)]=dzux-dxuz;
02412               break;
02413             case 13: 
02414               work[base::idx(i,j,k,Nx,Ny)]=dxuy-dyux;
02415               break;
02416             case 14: {
02417               DataType ox = dyuz-dzuy, oy = dzux-dxuz, oz = dxuy-dyux;
02418               work[base::idx(i,j,k,Nx,Ny)] = std::sqrt(ox*ox+oy*oy+oz*oz);
02419               break;
02420             }
02421             case 18:
02422               work[base::idx(i,j,k,Nx,Ny)] = DataType(-0.5)*(4.*dxux*dxux+4.*dyuy*dyuy+4.*dzuz*dzuz+4.*dxuy*dyux+4.*dxuz*dzux+4.*dyuz*dzuy);
02423               break;
02424             case 40:
02425               work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dxux; 
02426               break;
02427             case 41:
02428               work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuy+dyux); 
02429               break;
02430             case 42:
02431               work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dyuy; 
02432               break;
02433             case 43:
02434               work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuz+dzux); 
02435               break;
02436             case 44:
02437               work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dyuz+dzuy); 
02438               break;
02439             case 45:
02440               work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dzuz; 
02441               break;
02442             }
02443           }
02444     }
02445   }
02446 
02447   virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
02448                      const int skip_ghosts=1) const {
02449     int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
02450     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
02451     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
02452     MicroType *f = (MicroType *)fvec.databuffer();
02453     DataType *work = (DataType *)workvec.databuffer();
02454     DCoords dx = base::GH().worldStep(fvec.stepsize());
02455 
02456     for (register int k=mzs; k<=mze; k++)
02457       for (register int j=mys; j<=mye; j++)
02458         for (register int i=mxs; i<=mxe; i++) {
02459           switch(cnt) {
02460           case 1:
02461             f[base::idx(i,j,k,Nx,Ny)](0)=work[base::idx(i,j,k,Nx,Ny)]/R0;
02462             break;
02463           case 2:
02464             f[base::idx(i,j,k,Nx,Ny)](1)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
02465             break;
02466           case 3:
02467             f[base::idx(i,j,k,Nx,Ny)](2)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
02468             break;
02469           case 4:
02470             f[base::idx(i,j,k,Nx,Ny)](3)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
02471             f[base::idx(i,j,k,Nx,Ny)] = Equilibrium((const MacroType &)f[base::idx(i,j,k,Nx,Ny)]);
02472             break;
02473           }
02474         }
02475   }
02476 
02477   virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
02478                     const int verbose) const {
02479     int Nx = fvec.extents(0), Ny = fvec.extents(1);
02480     int b = fvec.bottom();
02481     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
02482     BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
02483     int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
02484     int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
02485     int mzs = std::max(bb.lower(2)/bb.stepsize(2),bbmax.lower(2)/bbmax.stepsize(2));
02486     int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
02487     int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
02488     int mze = std::min(bb.upper(2)/bb.stepsize(2),bbmax.upper(2)/bbmax.stepsize(2));
02489     MicroType *f = (MicroType *)fvec.databuffer();
02490 
02491     //std::cout << "CHECK Nx=" << Nx << "   Ny=" << Ny << "   Nz=" << Nz << std::endl;
02492     int result = 1;
02493     for (register int k=mzs; k<=mze; k++)
02494       for (register int j=mys; j<=mye; j++)
02495         for (register int i=mxs; i<=mxe; i++)
02496           for (register int l=0; l<base::NMicroVar(); l++)
02497             if (!(f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l)>-1.e37)) {
02498               result = 0;
02499               if (verbose) {
02500                 DCoords lbcorner = base::GH().worldCoords(fvec.lower(), fvec.stepsize());
02501                 DCoords dx = base::GH().worldStep(fvec.stepsize());
02502                 
02503                 std::cerr << "D3Q19-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << "," << k-mdz[l] << ")(" << l << ")="
02504                           << f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l) << " " << fvec.bbox() 
02505                           << " Coords: (" << (i+0.5)*dx(0)+lbcorner(0) << "," << (j+0.5)*dx(1)+lbcorner(1) << "," 
02506                           << (k+0.5)*dx(2)+lbcorner(2) << ")" << std::endl;
02507               }
02508             }
02509     return result;
02510   }
02511 
02512   inline const TensorType Stress(const MicroType &f, const MacroType &q, const DataType om) const {
02513     TensorType P;
02514     if (method[1]==0) {
02515       P(0) = q(0)*q(1)*q(1)-(f(1)+f(2)+f(7)+f(8)+f(9)+f(10)+f(15)+f(16)+f(17)+f(18))+cs2*q(0);
02516       P(1) = q(0)*q(1)*q(2)-(f(7)+f(8)-f(9)-f(10));
02517       P(2) = q(0)*q(2)*q(2)-(f(3)+f(4)+f(7)+f(8)+f(9)+f(10)+f(11)+f(12)+f(13)+f(14))+cs2*q(0);
02518       P(3) = q(0)*q(1)*q(3)-(f(15)+f(16)-f(17)-f(18));
02519       P(4) = q(0)*q(2)*q(3)-(f(11)+f(12)-f(13)-f(14));
02520       P(5) = q(0)*q(3)*q(3)-(f(5)+f(6)+f(11)+f(12)+f(13)+f(14)+f(15)+f(16)+f(17)+f(18))+cs2*q(0);
02521       P *= -(DataType(1.0)-DataType(0.5)*om);
02522     }
02523     else {
02524       MicroType feq = Equilibrium(q);     
02525       P = DeviatoricStress(f,feq,om);
02526     }
02527     P(0) -= LatticeBasePressure(q(0));  
02528     P(2) -= LatticeBasePressure(q(0));
02529     P(5) -= LatticeBasePressure(q(0));
02530     return P;
02531   }
02532 
02533   inline const TensorType DeviatoricStress(const MicroType &f, const MicroType &feq, const DataType om) const {
02534     TensorType Sigma;
02535     if (method[2]==0) {
02536       MicroType fneq = feq - f;
02537       Sigma(0) = fneq(1)+fneq(2)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
02538       Sigma(1) = fneq(7)+fneq(8)-fneq(9)-fneq(10);
02539       Sigma(2) = fneq(3)+fneq(4)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(11)+fneq(12)+fneq(13)+fneq(14);
02540       Sigma(3) = fneq(15)+fneq(16)-fneq(17)-fneq(18);
02541       Sigma(4) = fneq(11)+fneq(12)-fneq(13)-fneq(14);
02542       Sigma(5) = fneq(5)+fneq(6)+fneq(11)+fneq(12)+fneq(13)+fneq(14)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
02543       Sigma *= -(DataType(1.0)-DataType(0.5)*om);
02544     }
02545     else {
02546       MacroType q = MacroVariables(f);
02547       Sigma = Stress(f,q,om);
02548       Sigma(0) += LatticeBasePressure(q(0));  
02549       Sigma(2) += LatticeBasePressure(q(0));
02550       Sigma(5) += LatticeBasePressure(q(0));
02551     } 
02552     return Sigma;
02553   }
02554 
02555   virtual int NMethodOrder() const { return 2; }
02556 
02557   inline const DataType& L0() const { return base::LengthScale(); }
02558   inline const DataType& T0() const { return base::TimeScale(); }
02559   inline void SetDensityScale(const DataType r0) { R0 = r0; }
02560   inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
02561   inline void SetSpeedUp(const DataType s0) { S0 = s0; }
02562   inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
02563 
02564   inline const DataType& DensityScale() const { return R0; }
02565   inline const DataType VelocityScale() const { return U0/S0; }
02566   // Multiply all velocities and viscosity by S0 to speed up flow artificially, remove this scaling in output
02567   inline const DataType& SpeedUp() const { return S0; }
02568 
02569   // Lattice quantities for the operator
02570   inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
02571   inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
02572   inline virtual DataType LatticeBasePressure(const DataType rho) const {
02573     return (method[0]==0 ? cs2*rho/R0 : cs2*(rho-rhop/R0));
02574   }
02575 
02576   // Physical quantities for the operator
02577   inline void SetGas(DataType rho, DataType nu, DataType cs) {
02578     rhop = rho; nup = nu; csp = cs; cs2p = cs*cs;
02579     SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
02580   }
02581   inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
02582   inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q,
02583                                               const DataType dt) const {
02584     DataType M2 = 0.0;
02585     //x-y
02586     M2 += (f(7) - feq(7))*(f(7) - feq(7));
02587     M2 += (f(8) - feq(8))*(f(8) - feq(8));
02588     M2 += (f(9) - feq(9))*(f(9) - feq(9));
02589     M2 += (f(10) - feq(10))*(f(10) - feq(10));
02590     //y-z
02591     M2 += (f(11) - feq(11))*(f(11) - feq(11));
02592     M2 += (f(12) - feq(12))*(f(12) - feq(12));
02593     M2 += (f(13) - feq(13))*(f(13) - feq(13));
02594     M2 += (f(14) - feq(14))*(f(14) - feq(14));
02595     //x-z
02596     M2 += (f(15) - feq(15))*(f(15) - feq(15));
02597     M2 += (f(16) - feq(16))*(f(16) - feq(16));
02598     M2 += (f(17) - feq(17))*(f(17) - feq(17));
02599     M2 += (f(18) - feq(18))*(f(18) - feq(18));
02600     M2 = std::sqrt(M2);
02601 
02602     DataType tau = DataType(1.0)/Omega(dt);
02603     DataType om_turb =  (DataType(2.0)/( tau
02604                                          + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
02605 
02606     return ( om_turb );
02607   }
02608   inline const DataType Omega_LES_dynamic(const MicroType &f, const MicroType &feq, const MacroType& q,
02609                                           const DataType dt) const {
02610     DataType M2 = 0.0;
02611     //x-y
02612     M2 += (f(7) - feq(7))*(f(7) - feq(7));
02613     M2 += (f(8) - feq(8))*(f(8) - feq(8));
02614     M2 += (f(9) - feq(9))*(f(9) - feq(9));
02615     M2 += (f(10) - feq(10))*(f(10) - feq(10));
02616     //y-z
02617     M2 += (f(11) - feq(11))*(f(11) - feq(11));
02618     M2 += (f(12) - feq(12))*(f(12) - feq(12));
02619     M2 += (f(13) - feq(13))*(f(13) - feq(13));
02620     M2 += (f(14) - feq(14))*(f(14) - feq(14));
02621     //x-z
02622     M2 += (f(15) - feq(15))*(f(15) - feq(15));
02623     M2 += (f(16) - feq(16))*(f(16) - feq(16));
02624     M2 += (f(17) - feq(17))*(f(17) - feq(17));
02625     M2 += (f(18) - feq(18))*(f(18) - feq(18));
02626     M2 = std::sqrt(M2);
02627 
02628     DataType tau = DataType(1.0)/Omega(dt);
02629     DataType vel = sqrt( q(1)*q(1)+q(2)*q(2)+q(3)*q(3) );
02630     DataType Cs_SmagDYN = DataType(0.5)*L0()*L0()*( (vel*vel/csp*(M2-vel)) / ((M2-vel)*(M2-vel))  );
02631     DataType om_turb =  (DataType(2.0)/( tau
02632                                          + sqrt(tau*tau + (DataType(18.)*Cs_SmagDYN*M2)/(R0*q(0)*cs2/S0/S0)) ));
02633 
02634     return ( om_turb );
02635   }
02636   inline const int TurbulenceType() const { return turbulence; }
02637   inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
02638   inline const DataType& GasDensity() const { return rhop; }
02639   inline const DataType& GasViscosity() const { return nup; }
02640   inline const DataType& GasSpeedofSound() const { return csp; }
02641   inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
02642   { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
02643   // Quantities for output only
02644   inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
02645   inline virtual const DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); } 
02646   inline virtual const DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); } 
02647 
02648 protected:
02649   DataType cs2, cs22, cssq;
02650   DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp;
02651   DataType Cs_Smagorinsky, turbulence;
02652   int method[3], mdx[19], mdy[19], mdz[19];
02653 };
02654 
02655 
02656 #endif