amroc/lbm/LBMD3Q19_DR.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_DR_H
00007 #define LBM_D3Q19_DR_H
00008 
00017 #include "LBMD3Q19.h"
00018 #include <cfloat>
00019 
00020 #define LB_FACTOR 1.0e5
00021 
00051 typedef LBMD3Q19<DataType> LBMBaseType;
00052 
00053 template <class DataType>
00054 class LBMD3Q19DR : public LBMBaseType {
00055   typedef LBMBaseType base;
00056 public:
00057   typedef typename base::vec_grid_data_type vec_grid_data_type;
00058   typedef typename base::grid_data_type grid_data_type;
00059   typedef typename base::MicroType MicroType;
00060   typedef typename base::MacroType MacroType;
00061   typedef GridData<MacroType,3> macro_grid_data_type;
00062   typedef typename base::SideName SideName;
00063   typedef typename base::point_type point_type;
00064   typedef Vector<DataType,6> TensorType;
00065 
00066   LBMD3Q19DR() : base() {
00067     DRCoeff=0.005;
00068     DRFlow(0)=1.;
00069     DRFlow(1)=0.;
00070     DRFlow(2)=0.;
00071     DRFlow(3)=0.;
00072     Cs_Smagorinsky = DataType(0.2);
00073   }
00074 
00075   virtual ~LBMD3Q19DR() {}
00076 
00077   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00078     base::register_at(Ctrl,prefix);
00079     RegisterAt(base::LocCtrl,"DR_damping",DRCoeff);
00080     RegisterAt(base::LocCtrl,"DR_rho",DRFlow(0));
00081     RegisterAt(base::LocCtrl,"DR_u",DRFlow(1));
00082     RegisterAt(base::LocCtrl,"DR_v",DRFlow(2));
00083     RegisterAt(base::LocCtrl,"DR_w",DRFlow(3));
00084 
00085     RegisterAt(base::LocCtrl,"DR_xmin",DRmin[0]);
00086     RegisterAt(base::LocCtrl,"DR_xmax",DRmax[0]);
00087     RegisterAt(base::LocCtrl,"DR_ymin",DRmin[1]);
00088     RegisterAt(base::LocCtrl,"DR_ymax",DRmax[1]);
00089     RegisterAt(base::LocCtrl,"DR_zmin",DRmin[2]);
00090     RegisterAt(base::LocCtrl,"DR_zmax",DRmax[2]);
00091 
00092     RegisterAt(base::LocCtrl,"DR_Left",DRbc[0]);
00093     RegisterAt(base::LocCtrl,"DR_Right",DRbc[1]);
00094     RegisterAt(base::LocCtrl,"DR_Bottom",DRbc[2]);
00095     RegisterAt(base::LocCtrl,"DR_Top",DRbc[3]);
00096     RegisterAt(base::LocCtrl,"DR_Front",DRbc[4]);
00097     RegisterAt(base::LocCtrl,"DR_Back",DRbc[5]);
00098   }
00099 
00100   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00101     base::SetupData(gh,ghosts);
00102     base::GH().wlb(WLB);
00103     base::GH().wub(WUB);
00104 
00105     DRminThk[0] = DRmin[0]-WLB[0];
00106     DRminThk[1] = DRmin[1]-WLB[1];
00107     DRminThk[2] = DRmin[2]-WLB[2];
00108     DRmaxThk[0] = WUB[0]-DRmax[0];
00109     DRmaxThk[1] = WUB[1]-DRmax[1];
00110     DRmaxThk[2] = WUB[2]-DRmax[2];
00111 
00112     for (register int i=0; i<3; i++) {
00113       if (DRminThk[i]<DataType(0.)) DRminThk[i] = DataType(0.);
00114       if (DRmaxThk[i]<DataType(0.)) DRmaxThk[i] = DataType(0.);
00115     }
00116 
00117     std::cout << "wholebbox " << base::GH().wholebbox() << WLB[0]<<","<<WLB[1]<<","<<WLB[2] <<") max=(" <<WUB[0]<<","<<WUB[1]<<","<<WUB[2] <<")"<< std::endl;
00118     std::cout << "DR region surrounds box : min=(" << DRmin[0] <<","<< DRmin[1] <<","<< DRmin[2] <<") max=(" <<DRmax[0] <<","<<DRmax[1]<<","<<DRmax[2]<< ") [m]"<<std::endl;
00119     std::cout << "DR region thickness : min=(" << DRminThk[0]<<","<< DRminThk[1]<<","<<DRminThk[2] <<") max=(" <<DRmaxThk[0] <<","<<DRmaxThk[1]<<","<<DRmaxThk[2]<< ") [m]"<<std::endl;
00120     std::cout << "DR region thickness : min=(" << DRminThk[0]/L0() <<","<< DRminThk[1]/L0() <<","<< DRminThk[2]/L0()
00121               <<") max=(" <<DRmaxThk[0]/L0() <<","<<DRmaxThk[1]/L0() <<","<<DRmaxThk[2]/L0() << ") lattice lengths"<<std::endl;
00122     (std::cout << "DR Damping Coefficient " << DRCoeff <<std::endl
00123                << " DR mean flow " << DRFlow << std::endl).flush();
00124 
00125     DRFlow(0) /= R0;
00126     DRFlow(1) /= U0/S0;
00127     DRFlow(2) /= U0/S0;
00128     DRFlow(3) /= U0/S0;
00129     fmeq = Equilibrium(DRFlow);
00130   }
00131 
00132   virtual int DRInside(const DCoords lower, const DCoords upper) const {
00137     int position = 100;
00139     if ( ( lower(0) >= DRmin[0] ) && ( lower(0) <= DRmax[0] ) )
00140       if ( ( lower(1) >= DRmin[1] ) && ( lower(1) <= DRmax[1] ) )
00141         if ( ( lower(2) >= DRmin[2] ) && ( lower(2) <= DRmax[2] ) )
00142           if ( ( upper(0) >= DRmin[0] ) && ( upper(0) <= DRmax[0] ) )
00143             if ( ( upper(1) >= DRmin[1] ) && ( upper(1) <= DRmax[1] ) )
00144               if ( ( upper(2) >= DRmin[2] ) && ( upper(2) <= DRmax[2] ) ) {
00145                 position = 0;
00146               }
00147     return position;
00148   }
00149 
00150   inline virtual void DRCollision(MicroType &f, const MicroType fmeql, const DataType damp, const DataType om) const {
00152     MacroType q = MacroVariables(f);
00153     MicroType feq = Equilibrium(q);
00154     TensorType Sigma = DeviatoricStress(f,feq,om);
00155     DataType S = Strain(q(0),Sigma,om,Cs_Smagorinsky);
00156 #if NUMPLUS >= 3
00157     if ( (turbulence==LES_SmagorinskyMemory || turbulence==LES_dynamicMemory) && f(20)>DataType(1.0) ) {
00158       Sigma = DeviatoricStress(f,feq,f(20));
00159       if (turbulence==LES_SmagorinskyMemory)
00160         S = Strain(q(0),Sigma,f(20),Cs_Smagorinsky);
00161       else
00162         S = Strain(q(0),Sigma,f(20),f(21));
00163     }
00164 #endif
00165     DataType strainDamp = damp*S/S0;
00166     for (register int qi=0; qi<19; qi++) {
00167       f(qi) -= strainDamp*(f(qi)-feq(qi)) +  damp*(feq(qi)-fmeql(qi));
00168     }
00169   }
00170 
00171   virtual void DRDamp(vec_grid_data_type &fvec) const {
00178     int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00179     int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00180     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00181     MicroType *f = (MicroType *)fvec.databuffer();
00182     DCoords dx = base::GH().worldStep(fvec.stepsize());
00183     DataType dt_temp = dx(0)/VelocityScale();
00184     DataType om = Omega(dt_temp);
00185     DCoords wc;
00186     int lc_indx[3] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1), mzs*fvec.stepsize(2) };
00187     DataType WLB[3],  WUB[3];
00188     base::GH().wlb(WLB);
00189     base::GH().wub(WUB);
00190 
00191     for (register int k=mzs; k<=mze; k++) {
00192       lc_indx[2] = k*fvec.stepsize(2)+fvec.lower(2);
00193       for (register int j=mys; j<=mye; j++) {
00194         lc_indx[1] = j*fvec.stepsize(1)+fvec.lower(1);
00195         for (register int i=mxs; i<=mxe; i++) {
00196           lc_indx[0] = i*fvec.stepsize(0)+fvec.lower(0);
00197           wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00198           DataType d1 = (DRmaxThk[0]>DataType(0.) ? (wc(0)-DRmax[0])/DRmaxThk[0] : DataType(0.));
00199           DataType d2 = (DRminThk[0]>DataType(0.) ? (DRmin[0]-wc(0))/DRminThk[0] : DataType(0.));
00200           DataType d3 = (DRmaxThk[1]>DataType(0.) ? (wc(1)-DRmax[1])/DRmaxThk[1] : DataType(0.));
00201           DataType d4 = (DRminThk[1]>DataType(0.) ? (DRmin[1]-wc(1))/DRminThk[1] : DataType(0.));
00202           DataType d5 = (DRmaxThk[2]>DataType(0.) ? (wc(2)-DRmax[2])/DRmaxThk[2] : DataType(0.));
00203           DataType d6 = (DRminThk[2]>DataType(0.) ? (DRmin[2]-wc(2))/DRminThk[2] : DataType(0.));
00204           DataType dist = std::max(d1,std::max(d2,std::max(d3,std::max(d4,std::max(d5,d6)))));
00205           if (dist>DataType(0.)) {
00206             DataType dampVal = DRCoeff*dist*dist;
00207             MicroType feqtmp = fmeq;
00208             if (DRbc[0]==1) { // left no-slip wall
00209               if (wc(0)-WLB[0]<DRminThk[1] && wc(1)-WLB[1]<DRminThk[1]) { // bottom
00210                 dampVal *= (wc(0)-WLB[0])*(wc(0)-WLB[0])/DRminThk[1]/DRminThk[1];
00211               }
00212               else if (wc(0)-WLB[0]<DRmaxThk[1] && WUB[1]-wc(1)<DRmaxThk[1]) { // top
00213                 dampVal *= (wc(0)-WLB[0])*(wc(0)-WLB[0])/DRmaxThk[1]/DRmaxThk[1];
00214               }
00215               else if (wc(0)-WLB[0]<DRminThk[2] && wc(2)-WLB[2]<DRminThk[2]) { // front
00216                 dampVal *= (wc(0)-WLB[0])*(wc(0)-WLB[0])/DRminThk[2]/DRminThk[2];
00217               }
00218               else if (wc(0)-WLB[0]<DRmaxThk[2] && WUB[2]-wc(2)<DRmaxThk[2]) { // back
00219                 dampVal *= (wc(0)-WLB[0])*(wc(0)-WLB[0])/DRmaxThk[2]/DRmaxThk[2];
00220               }
00221             }
00222             if (DRbc[1]==1) { // right no-slip wall
00223               if (WUB[0]-wc(0)<DRminThk[1] && wc(1)-WLB[1]<DRminThk[1]) { // bottom
00224                 dampVal *= (WUB[0]-wc(0))*(WUB[0]-wc(0))/DRminThk[1]/DRminThk[1];
00225               }
00226               else if (WUB[0]-wc(0)<DRmaxThk[1] && WUB[1]-wc(1)<DRmaxThk[1]) { // top
00227                 dampVal *= (WUB[0]-wc(0))*(WUB[0]-wc(0))/DRmaxThk[1]/DRmaxThk[1];
00228               }
00229               else if (WUB[0]-wc(0)<DRminThk[2] && wc(2)-WLB[2]<DRminThk[2]) { // front
00230                 dampVal *= (WUB[0]-wc(0))*(WUB[0]-wc(0))/DRminThk[2]/DRminThk[2];
00231               }
00232               else if (WUB[0]-wc(0)<DRmaxThk[2] && WUB[2]-wc(2)<DRmaxThk[2]) { // back
00233                 dampVal *= (WUB[0]-wc(0))*(WUB[0]-wc(0))/DRmaxThk[2]/DRmaxThk[2];
00234               }
00235             }
00236             if (DRbc[2]==1) { // bottom no-slip wall
00237               if (wc(0)-WLB[0]<DRminThk[0] && wc(1)-WLB[1]<DRminThk[0]) { // left
00238                 dampVal *= (wc(1)-WLB[1])*(wc(1)-WLB[1])/DRminThk[0]/DRminThk[0];
00239               }
00240               else if (WUB[0]-wc(0)<DRmaxThk[0] && wc(1)-WLB[1]<DRmaxThk[0]) { // right
00241                 dampVal *= (wc(1)-WLB[1])*(wc(1)-WLB[1])/DRmaxThk[0]/DRmaxThk[0];
00242               }
00243               else if (wc(2)-WLB[2]<DRminThk[2] && wc(1)-WLB[1]<DRminThk[2]) { // front
00244                 dampVal *= (wc(1)-WLB[1])*(wc(1)-WLB[1])/DRminThk[2]/DRminThk[2];
00245               }
00246               else if (WUB[2]-wc(2)<DRmaxThk[2] && wc(1)-WLB[1]<DRmaxThk[2]) { // back
00247                 dampVal *= (wc(1)-WLB[1])*(wc(1)-WLB[1])/DRmaxThk[2]/DRmaxThk[2];
00248               }
00249             }
00250             if (DRbc[3]==1) { // top no-slip wall
00251               if (wc(0)-WLB[0]<DRminThk[0] && WUB[1]-wc(1)<DRminThk[0]) { // left
00252                 dampVal *= (WUB[1]-wc(1))*(WUB[1]-wc(1))/DRminThk[0]/DRminThk[0];
00253               }
00254               else if (WUB[0]-wc(0)<DRmaxThk[0] && WUB[1]-wc(1)<DRmaxThk[0]) { // right
00255                 dampVal *= (WUB[1]-wc(1))*(WUB[1]-wc(1))/DRmaxThk[0]/DRmaxThk[0];
00256               }
00257               else if (wc(2)-WLB[2]<DRminThk[2] && WUB[1]-wc(1)<DRminThk[2]) { // front
00258                 dampVal *= (WUB[1]-wc(1))*(WUB[1]-wc(1))/DRminThk[2]/DRminThk[2];
00259               }
00260               else if (WUB[2]-wc(2)<DRmaxThk[2] && WUB[1]-wc(1)<DRmaxThk[2]) { // back
00261                 dampVal *= (WUB[1]-wc(1))*(WUB[1]-wc(1))/DRmaxThk[2]/DRmaxThk[2];
00262               }
00263             }
00264             if (DRbc[4]==1) { // front no-slip wall
00265               if (wc(0)-WLB[0]<DRminThk[0] && wc(2)-WLB[2]<DRminThk[0]) { // left
00266                 dampVal *= (wc(2)-WLB[2])*(wc(2)-WLB[2])/DRminThk[0]/DRminThk[0];
00267               }
00268               else if (WUB[0]-wc(0)<DRmaxThk[0] && wc(2)-WLB[2]<DRmaxThk[0]) { // right
00269                 dampVal *= (wc(2)-WLB[2])*(wc(2)-WLB[2])/DRmaxThk[0]/DRmaxThk[0];
00270               }
00271               else if (wc(1)-WLB[1]<DRminThk[1] && wc(2)-WLB[2]<DRminThk[1]) { // bottom
00272                 dampVal *= (wc(2)-WLB[2])*(wc(2)-WLB[2])/DRminThk[1]/DRminThk[1];
00273               }
00274               else if (WUB[1]-wc(1)<DRmaxThk[1] && wc(2)-WLB[2]<DRmaxThk[1]) { // top
00275                 dampVal *= (wc(2)-WLB[2])*(wc(2)-WLB[2])/DRmaxThk[1]/DRmaxThk[1];
00276               }
00277             }
00278             if (DRbc[5]==1) { // back no-slip wall
00279               if (wc(0)-WLB[0]<DRminThk[0] && WUB[2]-wc(2)<DRminThk[0]) { // left
00280                 dampVal *= (WUB[2]-wc(2))*(WUB[2]-wc(2))/DRminThk[0]/DRminThk[0];
00281               }
00282               else if (WUB[0]-wc(0)<DRmaxThk[0] && WUB[2]-wc(2)<DRmaxThk[0]) { // right
00283                 dampVal *= (WUB[2]-wc(2))*(WUB[2]-wc(2))/DRmaxThk[0]/DRmaxThk[0];
00284               }
00285               else if (wc(1)-WLB[1]<DRminThk[1] && WUB[2]-wc(2)<DRminThk[1]) { // bottom
00286                 dampVal *= (WUB[2]-wc(2))*(WUB[2]-wc(2))/DRminThk[1]/DRminThk[1];
00287               }
00288               else if (WUB[1]-wc(1)<DRmaxThk[1] && WUB[2]-wc(2)<DRmaxThk[1]) { // top
00289                 dampVal *= (WUB[2]-wc(2))*(WUB[2]-wc(2))/DRmaxThk[1]/DRmaxThk[1];
00290               }
00291             }
00292             if (DRbc[0]>1 || DRbc[1]>1 || DRbc[2]>1 || DRbc[3]>1 || DRbc[4]>1 || DRbc[5]>1) {
00293               MacroType qtmp = MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
00294               if (wc(0)-WLB[0]<DRminThk[0]) { // left
00295                 switch(DRbc[0]) {
00296                 case 2:
00297                   qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00298                   qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00299                   qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00300                   break;
00301                 case 3:
00302                   qtmp(0) = DRFlow(0);
00303                   break;
00304                 case 4:
00305                 case 5:
00306                   if (DRbc[0]==4) {
00307                     qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00308                     qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00309                     qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00310                   }
00311                   else if (DRbc[0]==5)
00312                     qtmp(0) = DRFlow(0);
00313                   DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,k,Nx,Ny)]), MacroVariables(f[base::idx(i-1,j,k,Nx,Ny)]),
00314                       MacroVariables(f[base::idx(i,j+1,k,Nx,Ny)]), MacroVariables(f[base::idx(i,j-1,k,Nx,Ny)]),
00315                       MacroVariables(f[base::idx(i,j,k+1,Nx,Ny)]), MacroVariables(f[base::idx(i,j,k-1,Nx,Ny)]), dx);
00316                   if (Qcrit>DataType(0.)) {
00317 
00318                   }
00319                   else {
00320                     dampVal *= DataType(0.5);
00321                   }
00322                   break;
00323                 }
00324               }
00325               else if (WUB[0]-wc(0)<DRmaxThk[0]) { // right
00326                 switch(DRbc[1]) {
00327                 case 2:
00328                   qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00329                   qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00330                   qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00331                   break;
00332                 case 3:
00333                   qtmp(0) = DRFlow(0);
00334                   break;
00335                 case 4:
00336                 case 5:
00337                   if (DRbc[0]==4) {
00338                     qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00339                     qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00340                     qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00341                   }
00342                   else if (DRbc[0]==5)
00343                     qtmp(0) = DRFlow(0);
00344                   DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,k,Nx,Ny)]), MacroVariables(f[base::idx(i-1,j,k,Nx,Ny)]),
00345                       MacroVariables(f[base::idx(i,j+1,k,Nx,Ny)]), MacroVariables(f[base::idx(i,j-1,k,Nx,Ny)]),
00346                       MacroVariables(f[base::idx(i,j,k+1,Nx,Ny)]), MacroVariables(f[base::idx(i,j,k-1,Nx,Ny)]), dx);
00347                   if (Qcrit>DataType(0.)) {
00348 
00349                   }
00350                   else {
00351                     dampVal *= DataType(0.5);
00352                   }
00353                   break;
00354                 }
00355               }
00356               else if (wc(1)-WLB[1]<DRminThk[1]) { // bottom
00357                 switch(DRbc[2]) {
00358                 case 2:
00359                   qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00360                   qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00361                   qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00362                   break;
00363                 case 3:
00364                   qtmp(0) = DRFlow(0);
00365                   break;
00366                 case 4:
00367                 case 5:
00368                   if (DRbc[0]==4) {
00369                     qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00370                     qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00371                     qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00372                   }
00373                   else if (DRbc[0]==5)
00374                     qtmp(0) = DRFlow(0);
00375                   DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,k,Nx,Ny)]), MacroVariables(f[base::idx(i-1,j,k,Nx,Ny)]),
00376                       MacroVariables(f[base::idx(i,j+1,k,Nx,Ny)]), MacroVariables(f[base::idx(i,j-1,k,Nx,Ny)]),
00377                       MacroVariables(f[base::idx(i,j,k+1,Nx,Ny)]), MacroVariables(f[base::idx(i,j,k-1,Nx,Ny)]), dx);
00378                   if (Qcrit>DataType(0.)) {
00379 
00380                   }
00381                   else {
00382                     dampVal *= DataType(0.5);
00383                   }
00384                   break;
00385                 }
00386               }
00387               else if (WUB[1]-wc(1)<DRmaxThk[1]) { //top
00388                 switch(DRbc[3]) {
00389                 case 2:
00390                   qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00391                   qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00392                   qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00393                   break;
00394                 case 3:
00395                   qtmp(0) = DRFlow(0);
00396                   break;
00397                 case 4:
00398                 case 5:
00399                   if (DRbc[0]==4) {
00400                     qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00401                     qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00402                     qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00403                   }
00404                   else if (DRbc[0]==5)
00405                     qtmp(0) = DRFlow(0);
00406                   DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,k,Nx,Ny)]), MacroVariables(f[base::idx(i-1,j,k,Nx,Ny)]),
00407                       MacroVariables(f[base::idx(i,j+1,k,Nx,Ny)]), MacroVariables(f[base::idx(i,j-1,k,Nx,Ny)]),
00408                       MacroVariables(f[base::idx(i,j,k+1,Nx,Ny)]), MacroVariables(f[base::idx(i,j,k-1,Nx,Ny)]), dx);
00409                   if (Qcrit>DataType(0.)) {
00410 
00411                   }
00412                   else {
00413                     dampVal *= DataType(0.5);
00414                   }
00415                   break;
00416                 }
00417               }
00418               else if (wc(2)-WLB[2]<DRminThk[2]) { // front
00419                 switch(DRbc[4]) {
00420                 case 2:
00421                   qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00422                   qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00423                   qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00424                   break;
00425                 case 3:
00426                   qtmp(0) = DRFlow(0);
00427                   break;
00428                 case 4:
00429                 case 5:
00430                   if (DRbc[0]==4) {
00431                     qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00432                     qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00433                     qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00434                   }
00435                   else if (DRbc[0]==5)
00436                     qtmp(0) = DRFlow(0);
00437                   DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,k,Nx,Ny)]), MacroVariables(f[base::idx(i-1,j,k,Nx,Ny)]),
00438                       MacroVariables(f[base::idx(i,j+1,k,Nx,Ny)]), MacroVariables(f[base::idx(i,j-1,k,Nx,Ny)]),
00439                       MacroVariables(f[base::idx(i,j,k+1,Nx,Ny)]), MacroVariables(f[base::idx(i,j,k-1,Nx,Ny)]), dx);
00440                   if (Qcrit>DataType(0.)) {
00441 
00442                   }
00443                   else {
00444                     dampVal *= DataType(0.5);
00445                   }
00446                   break;
00447                 }
00448               }
00449               else if (WUB[2]-wc(2)<DRmaxThk[2]) { // back
00450                 switch(DRbc[5]) {
00451                 case 2:
00452                   qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00453                   qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00454                   qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00455                   break;
00456                 case 3:
00457                   qtmp(0) = DRFlow(0);
00458                   break;
00459                 case 4:
00460                 case 5:
00461                   if (DRbc[0]==4) {
00462                     qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00463                     qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00464                     qtmp(3) = DRFlow(3)*qtmp(0)/DRFlow(0);
00465                   }
00466                   else if (DRbc[0]==5)
00467                     qtmp(0) = DRFlow(0);
00468                   DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,k,Nx,Ny)]), MacroVariables(f[base::idx(i-1,j,k,Nx,Ny)]),
00469                       MacroVariables(f[base::idx(i,j+1,k,Nx,Ny)]), MacroVariables(f[base::idx(i,j-1,k,Nx,Ny)]),
00470                       MacroVariables(f[base::idx(i,j,k+1,Nx,Ny)]), MacroVariables(f[base::idx(i,j,k-1,Nx,Ny)]), dx);
00471                   if (Qcrit>DataType(0.)) {
00472 
00473                   }
00474                   else {
00475                     dampVal *= DataType(0.5);
00476                   }
00477                   break;
00478                 }
00479               }
00480               feqtmp = Equilibrium(qtmp);
00481             }
00482             if (dampVal>DataType(0.)) {
00483               DRCollision(f[base::idx(i,j,k,Nx,Ny)],feqtmp,dampVal,om);
00484             }
00485           }
00486         }
00487       }
00488     }
00489 
00490   }
00491 
00492   virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00493                       const double& t, const double& dt, const int& mpass) const {
00494 
00495     base::Step(fvec,ovec,Flux,t,dt,mpass);
00496 
00497     // Damping Region
00498     DCoords lower = base::GH().worldCoords(fvec.lower(), fvec.stepsize());
00499     DCoords upper = base::GH().worldCoords(fvec.upper(), fvec.stepsize());
00500     if (DRInside(lower, upper)!=0) {
00501       DRDamp(fvec);
00502     }
00503 
00504     return DataType(1.);
00505   }
00506 
00507   virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00508                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00509     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00510     int b = fvec.bottom();
00511     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00512             && fvec.stepsize(2)==bb.stepsize(2));
00513     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00514     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00515     int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00516     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00517     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00518     int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00519     MicroType *f = (MicroType *)fvec.databuffer();
00520 
00521     if (type==base::NoSlipWallEntranceExit) {
00522       int lc_indx[3] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1), mzs*fvec.stepsize(2) };
00523       DataType WLB[3],  WUB[3];
00524       base::GH().wlb(WLB);
00525       base::GH().wub(WUB);
00526       DCoords wc;
00527       DataType slip = 0.;
00528       DataType lxs, lxe, lys, lye, lzs, lze;
00529       DataType min[3], max[3];
00530       if (naux==6) {
00531         lxs = aux[0] - WLB[0];
00532         lxe = WUB[0] - aux[1];
00533         lys = aux[2] - WLB[1];
00534         lye = WUB[1] - aux[3];
00535         lzs = aux[4] - WLB[2];
00536         lze = WUB[2] - aux[5];
00537         min[0] = aux[0]; max[0] = aux[1];
00538         min[1] = aux[2]; max[1] = aux[3];
00539         min[2] = aux[4]; max[2] = aux[5];
00540       }
00541       else {
00542         lxs = DRmin[0] - WLB[0];
00543         lxe = WUB[0] - DRmax[0];
00544         lys = DRmin[1] - WLB[1];
00545         lye = WUB[1] - DRmax[1];
00546         lzs = DRmin[2] - WLB[2];
00547         lze = WUB[2] - DRmax[2];
00548         min[0] = DRmin[0]; max[0] = DRmax[0];
00549         min[1] = DRmin[1]; max[1] = DRmax[1];
00550         min[2] = DRmin[2]; max[2] = DRmax[2];
00551       }
00552 
00553       switch (side) {
00554       case base::Left:
00555         lc_indx[0] = mxe*fvec.stepsize(0);
00556         for (register int k=mzs; k<=mze; k++) {
00557           lc_indx[2] = k*fvec.stepsize(2);
00558           for (register int j=mys; j<=mye; j++) {
00559             lc_indx[1] = j*fvec.stepsize(1);
00560             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00561             if (wc(1)>=min[1] && wc(1)<=max[1]) {
00562               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);
00563               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);
00564             }
00565             else {
00566               if (wc(1)<WLB[1] || wc(1) >WUB[1])
00567                 slip = DataType(0.);
00568               else if (wc(1)<min[1])
00569                 slip = (wc(1) - WLB[1])/lys;
00570               else if (wc(1)>max[1])
00571                 slip = (WUB[1] - wc(1))/lye;
00572               if (j>mys && j<mye) {
00573                 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)
00574                     + slip*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
00575                 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)
00576                     + slip*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
00577               }
00578               else if (j==mys) {
00579                 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)
00580                     + slip*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
00581                 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)
00582                     + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);
00583               }
00584               else if (j==mye) {
00585                 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)
00586                     + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
00587                 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)
00588                     + slip*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
00589               }
00590             }
00591             if (wc(2)>=min[2] && wc(2)<=max[2]) {
00592               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);
00593               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);
00594             }
00595             else {
00596               if (wc(2)<WLB[2] || wc(2) >WUB[2])
00597                 slip = DataType(0.);
00598               else if (wc(2)<min[2])
00599                 slip = (wc(2) - WLB[2])/lzs;
00600               else if (wc(2)>max[2])
00601                 slip = (WUB[2] - wc(2))/lze;
00602               if (k>mzs && k<mze) {
00603                 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)
00604                     + slip*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
00605                 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)
00606                     + slip*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
00607               }
00608               else if (k==mzs) {
00609                 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)
00610                     + slip*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
00611                 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)
00612                     + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
00613               }
00614               else if (k==mze) {
00615                 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)
00616                     + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
00617                 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)
00618                     + slip*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
00619               }
00620             }
00621             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00622           }
00623         }
00624         break;
00625       case base::Right:
00626         lc_indx[0] = mxs*fvec.stepsize(0);
00627         for (register int k=mzs; k<=mze; k++) {
00628           lc_indx[2] = k*fvec.stepsize(2);
00629           for (register int j=mys; j<=mye; j++) {
00630             lc_indx[1] = j*fvec.stepsize(1);
00631             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00632             if (wc(1)>=min[1] && wc(1)<=max[1]) {
00633               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);
00634               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);
00635             }
00636 
00637             else {
00638               if (wc(1)<WLB[1] || wc(1) >WUB[1])
00639                 slip = DataType(0.);
00640               else if (wc(1)<min[1])
00641                 slip = (wc(1) - WLB[1])/lys;
00642               else if (wc(1)>max[1])
00643                 slip = (WUB[1] - wc(1))/lye;
00644               if (j>mys && j<mye) {
00645                 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)
00646                     + slip*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
00647                 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)
00648                     + slip*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
00649               }
00650               else if (j==mys) {
00651                 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)
00652                     + slip*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
00653                 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)
00654                     + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
00655               }
00656               else if (j==mye) {
00657                 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)
00658                     + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
00659                 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)
00660                     + slip*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
00661               }
00662             }
00663             if (wc(2)>=min[2] && wc(2)<=max[2]) {
00664               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);
00665               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);
00666             }
00667             else {
00668               if (wc(2)<WLB[2] || wc(2) >WUB[2])
00669                 slip = DataType(0.);
00670               else if (wc(2)<min[2])
00671                 slip = (wc(2) - WLB[2])/lzs;
00672               else if (wc(2)>max[2])
00673                 slip = (WUB[2] - wc(2))/lze;
00674               if (k>mzs && k<mze) {
00675                 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)
00676                     + slip*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
00677                 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)
00678                     + slip*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
00679               }
00680               else if (k==mzs) {
00681                 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)
00682                     + slip*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
00683                 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)
00684                     + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
00685               }
00686               else if (k==mze) {
00687                 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)
00688                     + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
00689                 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)
00690                     + slip*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
00691               }
00692             }
00693             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00694           }
00695         }
00696         break;
00697       case base::Bottom:
00698         lc_indx[1] = mye*fvec.stepsize(1);
00699         for (register int k=mzs; k<=mze; k++) {
00700           lc_indx[2] = k*fvec.stepsize(2);
00701           for (register int i=mxs; i<=mxe; i++) {
00702             lc_indx[0] = i*fvec.stepsize(0);
00703             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00704             if (wc(0)>=min[0] && wc(0)<=max[0]) {
00705               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);
00706               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);
00707             }
00708             else {
00709               if (wc(0)<WLB[0] || wc(0)>WUB[0])
00710                 slip = DataType(0.);
00711               else if (wc(0)<min[0])
00712                 slip = (wc(0) - WLB[0])/lxs;
00713               else if (wc(0)>max[0])
00714                 slip = (WUB[0] - wc(0))/lxe;
00715               if (i>mxs && i<mxe) {
00716                 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)
00717                     + slip*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
00718                 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)
00719                     + slip*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
00720               }
00721               else if (i==mxs) {
00722                 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)
00723                     + slip*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
00724                 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)
00725                     + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
00726               }
00727               else if (i==mxe) {
00728                 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)
00729                     + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
00730                 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)
00731                     + slip*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
00732               }
00733             }
00734             if (wc(2)>=min[2] && wc(2)<=max[2]) {
00735               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);
00736               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);
00737             }
00738             else {
00739               if (wc(2)<WLB[2] || wc(2) >WUB[2])
00740                 slip = DataType(0.);
00741               else if (wc(2)<min[2])
00742                 slip = (wc(2) - WLB[2])/lzs;
00743               else if (wc(2)>max[2])
00744                 slip = (WUB[2] - wc(2))/lze;
00745               if (k>mzs && k<mze) {
00746                 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)
00747                     + slip*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
00748                 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)
00749                     + slip*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
00750               }
00751               else if (k==mzs) {
00752                 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)
00753                     + slip*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
00754                 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)
00755                     + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](14);
00756               }
00757               else if (k==mze) {
00758                 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)
00759                     + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](12);
00760                 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)
00761                     + slip*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
00762               }
00763             }
00764             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00765           }
00766         }
00767         break;
00768       case base::Top:
00769         lc_indx[1] = mys*fvec.stepsize(1);
00770         for (register int k=mzs; k<=mze; k++) {
00771           lc_indx[2] = k*fvec.stepsize(2);
00772           for (register int i=mxs; i<=mxe; i++) {
00773             lc_indx[0] = i*fvec.stepsize(0);
00774             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00775             if (wc(0)>=min[0] && wc(0)<=max[0]) {
00776               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);
00777               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);
00778             }
00779             else {
00780               if (wc(0)<WLB[0] || wc(0) >WUB[0])
00781                 slip = DataType(0.);
00782               else if (wc(0)<min[0])
00783                 slip = (wc(0) - WLB[0])/lxs;
00784               else if (wc(0)>max[0])
00785                 slip = (WUB[0] - wc(0))/lxe;
00786               if (i>mxs && i<mxe) {
00787                 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)
00788                     + slip*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
00789                 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)
00790                     + slip*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
00791               }
00792               else if (i==mxs) {
00793                 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)
00794                     + slip*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
00795                 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)
00796                     + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
00797               }
00798               else if (i==mxe) {
00799                 f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7)
00800                     + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
00801                 f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10)
00802                     + slip*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
00803               }
00804             }
00805             if (wc(2)>=min[2] && wc(2)<=max[2]) {
00806               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);
00807               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);
00808             }
00809             else {
00810               if (wc(2)<WLB[2] || wc(2)>WUB[2])
00811                 slip = DataType(0.);
00812               else if (wc(2)<min[2])
00813                 slip = (wc(2) - WLB[2])/lzs;
00814               else if (wc(2)>max[2])
00815                 slip = (WUB[2] - wc(2))/lze;
00816               if (k>mzs && k<mze) {
00817                 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)
00818                     + slip*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
00819                 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)
00820                     + slip*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
00821               }
00822               else if (k==mzs) {
00823                 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)
00824                     + slip*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
00825                 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)
00826                     + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
00827               }
00828               else if (k==mze) {
00829                 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)
00830                     + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
00831                 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)
00832                     + slip*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
00833               }
00834             }
00835             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00836           }
00837         }
00838         break;
00839       case base::Front:
00840         lc_indx[2] = mze*fvec.stepsize(2);
00841         for (register int j=mys; j<=mye; j++) {
00842           lc_indx[1] = j*fvec.stepsize(1);
00843           for (register int i=mxs; i<=mxe; i++) {
00844             lc_indx[0] = i*fvec.stepsize(0);
00845             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00846             if (wc(0)>=min[0] && wc(0)<=max[0]) {
00847               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);
00848               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);
00849             }
00850             else {
00851               if (wc(0)<WLB[0] || wc(0) >WUB[0])
00852                 slip = DataType(0.);
00853               else if (wc(0)<min[0])
00854                 slip = (wc(0) - WLB[0])/lxs;
00855               else if (wc(0)>max[0])
00856                 slip = (WUB[0] - wc(0))/lxe;
00857               if (i>mxs && i<mxe) {
00858                 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)
00859                     + slip*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
00860                 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)
00861                     + slip*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
00862               }
00863               else if (i==mxs) {
00864                 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)
00865                     + slip*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
00866                 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)
00867                     + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
00868               }
00869               else if (i==mxe) {
00870                 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)
00871                     + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
00872                 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)
00873                     + slip*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
00874               }
00875             }
00876             if (wc(1)>=min[1] && wc(1)<=max[1]) {
00877               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);
00878               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);
00879             }
00880             else {
00881               if (wc(1)<WLB[1] || wc(1) >WUB[1])
00882                 slip = DataType(0.);
00883               else if (wc(1)<min[1])
00884                 slip = (wc(1) - WLB[1])/lys;
00885               else if (wc(1)>max[1])
00886                 slip = (WUB[1] - wc(1))/lye;
00887               if (j>mys && j<mye) {
00888                 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)
00889                     + slip*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
00890                 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)
00891                     + slip*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
00892               }
00893               else if (j==mys) {
00894                 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)
00895                     + slip*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
00896                 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)
00897                     + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
00898               }
00899               else if (j==mye) {
00900                 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)
00901                     + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
00902                 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)
00903                     + slip*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
00904               }
00905             }
00906             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00907           }
00908         }
00909         break;
00910       case base::Back:
00911         lc_indx[2] = mzs*fvec.stepsize(2);
00912         for (register int j=mys; j<=mye; j++) {
00913           lc_indx[1] = j*fvec.stepsize(1);
00914           for (register int i=mxs; i<=mxe; i++) {
00915             lc_indx[0] = i*fvec.stepsize(0);
00916             wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00917             if (wc(0)>=min[0] && wc(0)<=max[0]) {
00918               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);
00919               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);
00920             }
00921             else {
00922               if (wc(0)<WLB[0] || wc(0) >WUB[0])
00923                 slip = DataType(0.);
00924               else if (wc(0)<min[0])
00925                 slip = (wc(0) - WLB[0])/lxs;
00926               else if (wc(0)>max[0])
00927                 slip = (WUB[0] - wc(0))/lxe;
00928               if (i>mxs && i<mxe) {
00929                 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)
00930                     + slip*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
00931                 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)
00932                     + slip*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
00933               }
00934               else if (i==mxs) {
00935                 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)
00936                     + slip*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
00937                 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)
00938                     + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
00939               }
00940               else if (i==mxe) {
00941                 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)
00942                     + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
00943                 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)
00944                     + slip*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
00945               }
00946             }
00947             if (wc(1)>=min[1] && wc(1)<=max[1]) {
00948               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);
00949               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);
00950             }
00951             else {
00952               if (wc(1)<WLB[1] || wc(1) >WUB[1])
00953                 slip = DataType(0.);
00954               else if (wc(1)<min[1])
00955                 slip = (wc(1) - WLB[1])/lys;
00956               else if (wc(1)>max[1])
00957                 slip = (WUB[1] - wc(1))/lye;
00958               if (j>mys && j<mye) {
00959                 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)
00960                     + slip*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
00961                 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)
00962                     + slip*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
00963               }
00964               else if (j==mys) {
00965                 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)
00966                     + slip*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
00967                 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)
00968                     + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
00969               }
00970               else if (j==mye) {
00971                 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)
00972                     + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
00973                 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)
00974                     + slip*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
00975               }
00976             }
00977             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00978           }
00979         }
00980         break;
00981       }
00982     }
00983     else
00984       base::BCStandard(fvec,bb,type,side,aux,naux,scaling);
00985   }
00986 
00987 
00988   inline DataType QcritLocal(const MacroType qxp, const MacroType qxm,
00989                              const MacroType qyp, const MacroType qym,
00990                              const MacroType qzp, const MacroType qzm, const DCoords dx) const {
00991     DataType dxux=(qxp(1)-qxm(1))/(2.*dx(0));
00992     DataType dxuy=(qxp(2)-qxm(2))/(2.*dx(0));
00993     DataType dxuz=(qxp(3)-qxm(3))/(2.*dx(0));
00994     DataType dyux=(qyp(1)-qym(1))/(2.*dx(1));
00995     DataType dyuy=(qyp(2)-qym(2))/(2.*dx(1));
00996     DataType dyuz=(qyp(3)-qym(3))/(2.*dx(1));
00997     DataType dzux=(qyp(1)-qym(1))/(2.*dx(2));
00998     DataType dzuy=(qzp(2)-qzm(2))/(2.*dx(2));
00999     DataType dzuz=(qzp(3)-qzm(3))/(2.*dx(2));
01000     return DataType(-0.5)*(4.*dxux*dxux+4.*dyuy*dyuy+4.*dzuz*dzuz+8.*dxuy*dyux+8.*dxuz*dzux+8.*dyuz*dzuy);
01001   }
01002 
01003   inline virtual const MacroType DRflow() const { return DRFlow; }
01004 
01005 protected:
01006 
01007   int WLBI[3], WUBI[3], DRbc[6];
01008   DataType WLB[3], WUB[3];
01009   DataType DRmin[3], DRmax[3];
01010   DataType DRminThk[3], DRmaxThk[3];
01011   DataType DRCoeff;
01012   MacroType DRFlow;
01013   MicroType fmeq;
01014 };
01015 
01016 
01017 #endif