00001
00002
00003
00004
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) {
00209 if (wc(0)-WLB[0]<DRminThk[1] && wc(1)-WLB[1]<DRminThk[1]) {
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]) {
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]) {
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]) {
00219 dampVal *= (wc(0)-WLB[0])*(wc(0)-WLB[0])/DRmaxThk[2]/DRmaxThk[2];
00220 }
00221 }
00222 if (DRbc[1]==1) {
00223 if (WUB[0]-wc(0)<DRminThk[1] && wc(1)-WLB[1]<DRminThk[1]) {
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]) {
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]) {
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]) {
00233 dampVal *= (WUB[0]-wc(0))*(WUB[0]-wc(0))/DRmaxThk[2]/DRmaxThk[2];
00234 }
00235 }
00236 if (DRbc[2]==1) {
00237 if (wc(0)-WLB[0]<DRminThk[0] && wc(1)-WLB[1]<DRminThk[0]) {
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]) {
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]) {
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]) {
00247 dampVal *= (wc(1)-WLB[1])*(wc(1)-WLB[1])/DRmaxThk[2]/DRmaxThk[2];
00248 }
00249 }
00250 if (DRbc[3]==1) {
00251 if (wc(0)-WLB[0]<DRminThk[0] && WUB[1]-wc(1)<DRminThk[0]) {
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]) {
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]) {
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]) {
00261 dampVal *= (WUB[1]-wc(1))*(WUB[1]-wc(1))/DRmaxThk[2]/DRmaxThk[2];
00262 }
00263 }
00264 if (DRbc[4]==1) {
00265 if (wc(0)-WLB[0]<DRminThk[0] && wc(2)-WLB[2]<DRminThk[0]) {
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]) {
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]) {
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]) {
00275 dampVal *= (wc(2)-WLB[2])*(wc(2)-WLB[2])/DRmaxThk[1]/DRmaxThk[1];
00276 }
00277 }
00278 if (DRbc[5]==1) {
00279 if (wc(0)-WLB[0]<DRminThk[0] && WUB[2]-wc(2)<DRminThk[0]) {
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]) {
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]) {
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]) {
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]) {
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]) {
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]) {
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]) {
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]) {
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]) {
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
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