Blender V2.61 - r43446

solver_main.cpp

Go to the documentation of this file.
00001 
00004 /******************************************************************************
00005  *
00006  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
00007  * Copyright 2003-2006 Nils Thuerey
00008  *
00009  * Standard LBM Factory implementation
00010  *
00011  *****************************************************************************/
00012 
00013 #include "solver_class.h"
00014 #include "solver_relax.h"
00015 #include "particletracer.h"
00016 #include "loop_tools.h"
00017 #include <stdlib.h>
00018 
00019 /*****************************************************************************/
00021 /*****************************************************************************/
00022 
00023 double globdfcnt;
00024 double globdfavg[19];
00025 double globdfmax[19];
00026 double globdfmin[19];
00027 extern int glob_mpindex,glob_mpnum;
00028 extern bool globOutstrForce;
00029 
00030 // simulation object interface
00031 void LbmFsgrSolver::step() { 
00032     stepMain();
00033 }
00034 
00035 // lbm main step
00036 void messageOutputForce(string from);
00037 void LbmFsgrSolver::stepMain() { 
00038     myTime_t timestart = getTime();
00039 
00040     initLevelOmegas();
00041     markedClearList(); // DMC clearMarkedCellsList
00042 
00043     // safety check, counter reset
00044     mNumUsedCells = 0;
00045     mNumInterdCells = 0;
00046     mNumInvIfCells = 0;
00047 
00048   //debugOutNnl("LbmFsgrSolver::step : "<<mStepCnt, 10);
00049   if(!mSilent){ debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime, 10); }
00050     //debMsgDirect(  "LbmFsgrSolver::step : "<<mStepCnt<<" ");
00051     //myTime_t timestart = 0;
00052     //if(mStartSymm) { checkSymmetry("step1"); } // DEBUG 
00053 
00054     // time adapt
00055     mMaxVlen = mMxvz = mMxvy = mMxvx = 0.0;
00056 
00057     // init moving bc's, can change mMaxVlen
00058     initMovingObstacles(false);
00059     
00060     // handle fluid control 
00061     handleCpdata();
00062 
00063     // important - keep for tadap
00064     LbmFloat lastMass = mCurrentMass;
00065     mCurrentMass = mFixMass; // reset here for next step
00066     mCurrentVolume = 0.0;
00067     
00068     //change to single step advance!
00069     int levsteps = 0;
00070     int dsbits = mStepCnt ^ (mStepCnt-1);
00071     //errMsg("S"," step:"<<mStepCnt<<" s-1:"<<(mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits));
00072     for(int lev=0; lev<=mMaxRefine; lev++) {
00073         //if(! (mStepCnt&(1<<lev)) ) {
00074         if( dsbits & (1<<(mMaxRefine-lev)) ) {
00075             //errMsg("S"," l"<<lev);
00076 
00077             if(lev==mMaxRefine) {
00078                 // always advance fine level...
00079                 fineAdvance();
00080             } else {
00081                 adaptGrid(lev);
00082                 coarseRestrictFromFine(lev);
00083                 coarseAdvance(lev);
00084             }
00085 #if FSGR_OMEGA_DEBUG==1
00086             errMsg("LbmFsgrSolver::step","LES stats l="<<lev<<" omega="<<mLevel[lev].omega<<" avgOmega="<< (mLevel[lev].avgOmega/mLevel[lev].avgOmegaCnt) );
00087             mLevel[lev].avgOmega = 0.0; mLevel[lev].avgOmegaCnt = 0.0;
00088 #endif // FSGR_OMEGA_DEBUG==1
00089             levsteps++;
00090         }
00091         mCurrentMass   += mLevel[lev].lmass;
00092         mCurrentVolume += mLevel[lev].lvolume;
00093     }
00094 
00095   // prepare next step
00096     mStepCnt++;
00097 
00098 
00099     // some dbugging output follows
00100     // calculate MLSUPS
00101     myTime_t timeend = getTime();
00102 
00103     mNumUsedCells += mNumInterdCells; // count both types for MLSUPS
00104     mAvgNumUsedCells += mNumUsedCells;
00105     mMLSUPS = ((double)mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000.0);
00106     if(mMLSUPS>10000){ mMLSUPS = -1; }
00107     //else { mAvgMLSUPS += mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
00108     
00109     LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1())) / ((timeend-timestart)/(double)1000.0) ) / (1000000);
00110     if(totMLSUPS>10000) totMLSUPS = -1;
00111     mNumInvIfTotal += mNumInvIfCells; // debug
00112 
00113   // do some formatting 
00114   if(!mSilent){ 
00115         int avgcls = (int)(mAvgNumUsedCells/(LONGINT)mStepCnt);
00116     debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime<<
00117             " cur-mlsups:"<<mMLSUPS<< //" avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< 
00118             " totcls:"<<mNumUsedCells<< " avgcls:"<< avgcls<< 
00119             " intd:"<<mNumInterdCells<< " invif:"<<mNumInvIfCells<< 
00120             " invift:"<<mNumInvIfTotal<< " fsgrcs:"<<mNumFsgrChanges<< 
00121             " filled:"<<mNumFilledCells<<", emptied:"<<mNumEmptiedCells<< 
00122             " mMxv:"<<PRINT_VEC(mMxvx,mMxvy,mMxvz)<<", tscnts:"<<mTimeSwitchCounts<< 
00123             //" RWmxv:"<<ntlVec3Gfx(mMxvx,mMxvy,mMxvz)*(mLevel[mMaxRefine].simCellSize / mLevel[mMaxRefine].timestep)<<" "<< /* realworld vel output */
00124             " probs:"<<mNumProblems<< " simt:"<<mSimulationTime<< 
00125             " took:"<< getTimeString(timeend-timestart)<<
00126             " for '"<<mName<<"' " , 10);
00127     } else { debMsgDirect("."); }
00128 
00129     if(mStepCnt==1) {
00130         mMinNoCells = mMaxNoCells = mNumUsedCells;
00131     } else {
00132         if(mNumUsedCells>mMaxNoCells) mMaxNoCells = mNumUsedCells;
00133         if(mNumUsedCells<mMinNoCells) mMinNoCells = mNumUsedCells;
00134     }
00135     
00136     // mass scale test
00137     if((mMaxRefine>0)&&(mInitialMass>0.0)) {
00138         LbmFloat mscale = mInitialMass/mCurrentMass;
00139 
00140         mscale = 1.0;
00141         const LbmFloat dchh = 0.001;
00142         if(mCurrentMass<mInitialMass) mscale = 1.0+dchh;
00143         if(mCurrentMass>mInitialMass) mscale = 1.0-dchh;
00144 
00145         // use mass rescaling?
00146         // with float precision this seems to be nonsense...
00147         const bool MREnable = false;
00148 
00149         const int MSInter = 2;
00150         static int mscount = 0;
00151         if( (MREnable) && ((mLevel[0].lsteps%MSInter)== (MSInter-1)) && ( ABS( (mInitialMass/mCurrentMass)-1.0 ) > 0.01) && ( dsbits & (1<<(mMaxRefine-0)) ) ){
00152             // example: FORCE RESCALE MASS! ini:1843.5, cur:1817.6, f=1.01425 step:22153 levstep:5539 msc:37
00153             // mass rescale MASS RESCALE check
00154             errMsg("MDTDD","\n\n");
00155             errMsg("MDTDD","FORCE RESCALE MASS! "
00156                     <<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass)
00157                     <<" step:"<<mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
00158                     );
00159             errMsg("MDTDD","\n\n");
00160 
00161             mscount++;
00162             for(int lev=mMaxRefine; lev>=0 ; lev--) {
00163                 //for(int workSet = 0; workSet<=1; workSet++) {
00164                 int wss = 0;
00165                 int wse = 1;
00166 #if COMPRESSGRIDS==1
00167                 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
00168 #endif // COMPRESSGRIDS==1
00169                 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
00170 
00171                     FSGR_FORIJK1(lev) {
00172                         if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) 
00173                             ) {
00174 
00175                             FORDF0 { QCELL(lev, i,j,k,workSet, l) *= mscale; }
00176                             QCELL(lev, i,j,k,workSet, dMass) *= mscale;
00177                             QCELL(lev, i,j,k,workSet, dFfrac) *= mscale;
00178 
00179                         } else {
00180                             continue;
00181                         }
00182                     }
00183                 }
00184                 mLevel[lev].lmass *= mscale;
00185             }
00186         } 
00187 
00188         mCurrentMass *= mscale;
00189     }// if mass scale test */
00190     else {
00191         // use current mass after full step for initial setting
00192         if((mMaxRefine>0)&&(mInitialMass<=0.0) && (levsteps == (mMaxRefine+1))) {
00193             mInitialMass = mCurrentMass;
00194             debMsgStd("MDTDD",DM_NOTIFY,"Second Initial Mass Init: "<<mInitialMass, 2);
00195         }
00196     }
00197 
00198 #if LBM_INCLUDE_TESTSOLVERS==1
00199     if((mUseTestdata)&&(mInitDone)) { handleTestdata(); }
00200     mrExchange();
00201 #endif
00202 
00203     // advance positions with current grid
00204     advanceParticles();
00205     if(mpParticles) {
00206         mpParticles->checkDumpTextPositions(mSimulationTime);
00207         mpParticles->checkTrails(mSimulationTime);
00208     }
00209 
00210     // one of the last things to do - adapt timestep
00211     // was in fineAdvance before... 
00212     if(mTimeAdap) {
00213         adaptTimestep();
00214     } // time adaptivity
00215 
00216 
00217 #ifndef WIN32
00218     // good indicator for instabilities
00219     if( (!finite(mMxvx)) || (!finite(mMxvy)) || (!finite(mMxvz)) ) { CAUSE_PANIC; }
00220     if( (!finite(mCurrentMass)) || (!finite(mCurrentVolume)) ) { CAUSE_PANIC; }
00221 #endif // WIN32
00222 
00223     // output total step time
00224     myTime_t timeend2 = getTime();
00225     double mdelta = (lastMass-mCurrentMass);
00226     if(ABS(mdelta)<1e-12) mdelta=0.;
00227     double effMLSUPS = ((double)mNumUsedCells / ((timeend2-timestart)/(double)1000.0) ) / (1000000.0);
00228     if(mInitDone) {
00229         if(effMLSUPS>10000){ effMLSUPS = -1; }
00230         else { mAvgMLSUPS += effMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
00231     }
00232     
00233     debMsgStd("LbmFsgrSolver::stepMain", DM_MSG, "mmpid:"<<glob_mpindex<<" step:"<<mStepCnt
00234             <<" dccd="<< mCurrentMass
00235             //<<",d"<<mdelta
00236             //<<",ds"<<(mCurrentMass-mObjectMassMovnd[1])
00237             //<<"/"<<mCurrentVolume<<"(fix="<<mFixMass<<",ini="<<mInitialMass<<"), "
00238             <<" effMLSUPS=("<< effMLSUPS
00239             <<",avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "
00240             <<" took totst:"<< getTimeString(timeend2-timestart), 3);
00241     // nicer output
00242     //debMsgDirect(std::endl); 
00243     // */
00244 
00245     messageOutputForce("");
00246  //#endif // ELBEEM_PLUGIN!=1
00247 }
00248 
00249 #define NEWDEBCHECK(str) \
00250     if(!this->mPanic){ FSGR_FORIJK_BOUNDS(mMaxRefine) { \
00251         if(RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr)&(CFFluid|CFInter)) { \
00252         for(int l=0;l<dTotalNum;l++) { \
00253             if(!finite(QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr,l))) { errMsg("NNOFIN"," "<<str<<" at "<<PRINT_IJK<<" l"<<l<<" "); }\
00254         }/*for*/ \
00255         }/*if*/ \
00256     } }
00257 
00258 void LbmFsgrSolver::fineAdvance()
00259 {
00260     // do the real thing...
00261     //NEWDEBCHECK("t1");
00262     mainLoop( mMaxRefine );
00263     if(mUpdateFVHeight) {
00264         // warning assume -Y gravity...
00265         mFVHeight = mCurrentMass*mFVArea/((LbmFloat)(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizez));
00266         if(mFVHeight<1.0) mFVHeight = 1.0;
00267         mpParam->setFluidVolumeHeight(mFVHeight);
00268     }
00269 
00270     // advance time before timestep change
00271     mSimulationTime += mpParam->getTimestep();
00272     // time adaptivity
00273     mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) );
00274     //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG 
00275     if(!mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps), 3 ); }
00276 
00277     // update other set
00278   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
00279   mLevel[mMaxRefine].setCurr   ^= 1;
00280   mLevel[mMaxRefine].lsteps++;
00281 
00282     // flag init... (work on current set, to simplify flag checks)
00283     reinitFlags( mLevel[mMaxRefine].setCurr );
00284     if(!mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," flags reinit on set "<< mLevel[mMaxRefine].setCurr, 3 ); }
00285 
00286     // DEBUG VEL CHECK
00287     if(0) {
00288         int lev = mMaxRefine;
00289         int workSet = mLevel[lev].setCurr;
00290         int mi=0,mj=0,mk=0;
00291         LbmFloat compRho=0.;
00292         LbmFloat compUx=0.;
00293         LbmFloat compUy=0.;
00294         LbmFloat compUz=0.;
00295         LbmFloat maxUlen=0.;
00296         LbmVec maxU(0.);
00297         LbmFloat maxRho=-100.;
00298         int ri=0,rj=0,rk=0;
00299 
00300         FSGR_FORIJK1(lev) {
00301             if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) ) {
00302                 compRho=QCELL(lev, i,j,k,workSet, dC);
00303                 compUx = compUy = compUz = 0.0;
00304                 for(int l=1; l<this->cDfNum; l++) {
00305                     LbmFloat df = QCELL(lev, i,j,k,workSet, l);
00306                     compRho += df;
00307                     compUx  += (this->dfDvecX[l]*df);
00308                     compUy  += (this->dfDvecY[l]*df); 
00309                     compUz  += (this->dfDvecZ[l]*df); 
00310                 } 
00311                 LbmVec u(compUx,compUy,compUz);
00312                 LbmFloat nu = norm(u);
00313                 if(nu>maxUlen) {
00314                     maxU = u;
00315                     maxUlen = nu;
00316                     mi=i; mj=j; mk=k;
00317                 }
00318                 if(compRho>maxRho) {
00319                     maxRho=compRho;
00320                     ri=i; rj=j; rk=k;
00321                 }
00322             } else {
00323                 continue;
00324             }
00325         }
00326 
00327         errMsg("MAXVELCHECK"," at "<<PRINT_VEC(mi,mj,mk)<<" norm:"<<maxUlen<<" u:"<<maxU);
00328         errMsg("MAXRHOCHECK"," at "<<PRINT_VEC(ri,rj,rk)<<" rho:"<<maxRho);
00329         printLbmCell(lev, 30,36,23, -1);
00330     } // DEBUG VEL CHECK
00331 
00332 }
00333 
00334 
00335 
00336 // fine step defines
00337 
00338 // access to own dfs during step (may be changed to local array)
00339 #define MYDF(l) RAC(ccel, l)
00340 
00341 // drop model definitions
00342 #define RWVEL_THRESH 1.5
00343 #define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5)
00344 
00345 #if LBMDIM==3
00346 // normal
00347 #define SLOWDOWNREGION (mSizez/4)
00348 #else // LBMDIM==2
00349 // off
00350 #define SLOWDOWNREGION 10 
00351 #endif // LBMDIM==2
00352 #define P_LCSMQO 0.01
00353 
00354 /*****************************************************************************/
00356 /*****************************************************************************/
00357 void 
00358 LbmFsgrSolver::mainLoop(int lev)
00359 {
00360     // loops over _only inner_ cells  -----------------------------------------------------------------------------------
00361     
00362     // slow surf regions smooth (if below)
00363     const LbmFloat smoothStrength = 0.0; //0.01;
00364     const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03;
00365     const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit;
00366 
00367     const int cutMin  = 1;
00368     const int cutConst = mCutoff+2;
00369 
00370 
00371 #   if LBM_INCLUDE_TESTSOLVERS==1
00372     // 3d region off... quit
00373     if((mUseTestdata)&&(mpTest->mFarfMode>0)) { return; }
00374 #   endif // ELBEEM_PLUGIN!=1
00375     
00376   // main loop region
00377     const bool doReduce = true;
00378     const int gridLoopBound=1;
00379     GRID_REGION_INIT();
00380 #if PARALLEL==1
00381 #pragma omp parallel default(shared) \
00382   reduction(+: \
00383       calcCurrentMass,calcCurrentVolume, \
00384         calcCellsFilled,calcCellsEmptied, \
00385         calcNumUsedCells )
00386     GRID_REGION_START();
00387 #else // PARALLEL==1
00388     GRID_REGION_START();
00389 #endif // PARALLEL==1
00390 
00391     // local to main
00392     CellFlagType nbflag[LBM_DFNUM]; 
00393     int oldFlag, newFlag, nbored;
00394     LbmFloat m[LBM_DFNUM];
00395     LbmFloat rho, ux, uy, uz, tmp, usqr;
00396 
00397     // smago vars
00398     LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
00399     
00400     // ifempty cell conversion flags
00401     bool iffilled, ifemptied;
00402     LbmFloat nbfracs[LBM_DFNUM]; // ffracs of neighbors
00403     int recons[LBM_DFNUM];   // reconstruct this DF?
00404     int numRecons;           // how many are reconstructed?
00405 
00406     LbmFloat mass=0., change=0., lcsmqo=0.;
00407     rho= ux= uy= uz= usqr= tmp= 0.; 
00408     lcsmqadd = lcsmomega = 0.;
00409     FORDF0{ lcsmeq[l] = 0.; }
00410 
00411     // ---
00412     // now stream etc.
00413     // use //template functions for 2D/3D
00414 
00415     GRID_LOOP_START();
00416         // loop start
00417         // stream from current set to other, then collide and store
00418         //errMsg("l2"," at "<<PRINT_IJK<<" id"<<id);
00419 
00420 #       if FSGR_STRICT_DEBUG==1
00421         // safety pointer check
00422         rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
00423         if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) || 
00424             (&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) {
00425             errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<<
00426                     RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<<
00427                     (*pFlagSrc)<<","<<(*pFlagDst) <<",  pointers "<<
00428           (long)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(long)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
00429           (long)(pFlagSrc)<<","<<(long)(pFlagDst)<<" "
00430                     ); 
00431             CAUSE_PANIC;
00432         }   
00433         if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) || 
00434             (&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) {
00435             errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<<
00436           (long)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(long)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
00437           (long)(ccel)<<","<<(long)(tcel)<<" "
00438                     ); 
00439             CAUSE_PANIC;
00440         }   
00441 #       endif
00442         oldFlag = *pFlagSrc;
00443         
00444         // old INTCFCOARSETEST==1
00445         if( (oldFlag & (CFGrFromCoarse)) ) { 
00446             if(( mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
00447                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
00448             } else {
00449                 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
00450                 calcNumUsedCells++;
00451             }
00452             continue; // interpolateFineFromCoarse test!
00453         } // interpolateFineFromCoarse test! 
00454     
00455         if(oldFlag & (CFMbndInflow)) {
00456             // fluid & if are ok, fill if later on
00457             int isValid = oldFlag & (CFFluid|CFInter);
00458             const LbmFloat iniRho = 1.0;
00459             const int OId = oldFlag>>24;
00460             if(!isValid) {
00461                 // make new if cell
00462                 const LbmVec vel(mObjectSpeeds[OId]);
00463                 // TODO add OPT3D treatment
00464                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
00465                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
00466                 RAC(tcel, dFlux) = FLUX_INIT;
00467                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
00468                 calcCurrentMass += iniRho; 
00469                 calcCurrentVolume += 1.0; 
00470                 calcNumUsedCells++;
00471                 mInitialMass += iniRho;
00472                 // dont treat cell until next step
00473                 continue;
00474             } 
00475         } 
00476         else  // these are exclusive
00477         if(oldFlag & (CFMbndOutflow)) {
00478             int isnotValid = oldFlag & (CFFluid);
00479             if(isnotValid) {
00480                 // remove fluid cells, shouldnt be here anyway
00481                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
00482                 mInitialMass -= fluidRho;
00483                 const LbmFloat iniRho = 0.0;
00484                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
00485                 RAC(tcel, dFlux) = FLUX_INIT;
00486                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
00487 
00488                 // same as ifemptied for if below
00489                 LbmPoint oemptyp; oemptyp.flag = 0;
00490                 oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
00491                 LIST_EMPTY(oemptyp);
00492                 calcCellsEmptied++;
00493                 continue;
00494             }
00495         }
00496 
00497         if(oldFlag & (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)) { 
00498             *pFlagDst = oldFlag;
00499             continue;
00500         }
00501         /*if( oldFlag & CFNoBndFluid ) {  // TEST ME FASTER?
00502             OPTIMIZED_STREAMCOLLIDE; PERFORM_USQRMAXCHECK;
00503             RAC(tcel,dFfrac) = 1.0; 
00504             *pFlagDst = (CellFlagType)oldFlag; // newFlag;
00505             calcCurrentMass += rho; calcCurrentVolume += 1.0;
00506             calcNumUsedCells++;
00507             continue;
00508         }// TEST ME FASTER? */
00509 
00510         // only neighbor flags! not own flag
00511         nbored = 0;
00512         
00513 #if OPT3D==0
00514         FORDF1 {
00515             nbflag[l] = RFLAG_NB(lev, i,j,k,SRCS(lev),l);
00516             nbored |= nbflag[l];
00517         } 
00518 #else
00519         nbflag[dSB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dSB];
00520         nbflag[dWB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWB];
00521         nbflag[ dB] = *(pFlagSrc + (-mLevel[lev].lOffsy)); nbored |= nbflag[dB];
00522         nbflag[dEB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dEB];
00523         nbflag[dNB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNB];
00524 
00525         nbflag[dSW] = *(pFlagSrc + (-mLevel[lev].lOffsx+-1)); nbored |= nbflag[dSW];
00526         nbflag[ dS] = *(pFlagSrc + (-mLevel[lev].lOffsx)); nbored |= nbflag[dS];
00527         nbflag[dSE] = *(pFlagSrc + (-mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dSE];
00528 
00529         nbflag[ dW] = *(pFlagSrc + (-1)); nbored |= nbflag[dW];
00530         nbflag[ dE] = *(pFlagSrc + ( 1)); nbored |= nbflag[dE];
00531 
00532         nbflag[dNW] = *(pFlagSrc + ( mLevel[lev].lOffsx+-1)); nbored |= nbflag[dNW];
00533       nbflag[ dN] = *(pFlagSrc + ( mLevel[lev].lOffsx)); nbored |= nbflag[dN];
00534         nbflag[dNE] = *(pFlagSrc + ( mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dNE];
00535 
00536         nbflag[dST] = *(pFlagSrc + ( mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dST];
00537         nbflag[dWT] = *(pFlagSrc + ( mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWT];
00538         nbflag[ dT] = *(pFlagSrc + ( mLevel[lev].lOffsy)); nbored |= nbflag[dT];
00539         nbflag[dET] = *(pFlagSrc + ( mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dET];
00540         nbflag[dNT] = *(pFlagSrc + ( mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNT];
00541         // */
00542 #endif
00543 
00544         // pointer to destination cell
00545         calcNumUsedCells++;
00546 
00547         // FLUID cells 
00548         if( oldFlag & CFFluid ) { 
00549             // only standard fluid cells (with nothing except fluid as nbs
00550 
00551             if(oldFlag&CFMbndInflow) {
00552                 // force velocity for inflow, necessary to have constant direction of flow
00553                 // FIXME , test also set interface cells?
00554                 const int OId = oldFlag>>24;
00555                 //? DEFAULT_STREAM;
00556                 //const LbmFloat fluidRho = 1.0;
00557                 // for submerged inflows, streaming would have to be performed...
00558                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
00559                 const LbmVec vel(mObjectSpeeds[OId]);
00560                 ux=vel[0], uy=vel[1], uz=vel[2]; 
00561                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
00562                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); }
00563                 rho = fluidRho;
00564                 //errMsg("INFLOW_DEBUG","std at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
00565             } else {
00566                 if(nbored&CFBnd) {
00567                     DEFAULT_STREAM;
00568                     //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
00569                     DEFAULT_COLLIDEG(mLevel[lev].gravity);
00570                     oldFlag &= (~CFNoBndFluid);
00571                 } else {
00572                     // do standard stream/collide
00573                     OPTIMIZED_STREAMCOLLIDE;
00574                     oldFlag |= CFNoBndFluid;
00575                 } 
00576             }
00577 
00578             PERFORM_USQRMAXCHECK;
00579             // "normal" fluid cells
00580             RAC(tcel,dFfrac) = 1.0; 
00581             *pFlagDst = (CellFlagType)oldFlag; // newFlag;
00582             calcCurrentMass += rho; 
00583             calcCurrentVolume += 1.0;
00584             continue;
00585         }
00586         
00587         newFlag  = oldFlag;
00588         // make sure here: always check which flags to really unset...!
00589         newFlag = newFlag & (~( 
00590                     CFNoNbFluid|CFNoNbEmpty| CFNoDelete 
00591                     | CFNoInterpolSrc
00592                     | CFNoBndFluid
00593                     ));
00594         if(!(nbored&CFBndNoslip)) { //NEWSURFT NEWSURFTNOS
00595             newFlag |= CFNoBndFluid;
00596         }
00597         /*if(!(nbored&CFBnd)) { //NEWSURFT NEWSURFTNOS
00598             // explicitly check for noslip neighbors
00599             bool hasnoslipnb = false;
00600             FORDF1 { if((nbflag[l]&CFBnd)&&(nbflag[l]&CFBndNoslip)) hasnoslipnb=true; }
00601             if(!hasnoslipnb) newFlag |= CFNoBndFluid; 
00602         } // */
00603 
00604         // store own dfs and mass
00605         mass = RAC(ccel,dMass);
00606 
00607         // WARNING - only interface cells arrive here!
00608         // read distribution funtions of adjacent cells = stream step
00609         DEFAULT_STREAM;
00610 
00611         if((nbored & CFFluid)==0) { newFlag |= CFNoNbFluid; mNumInvIfCells++; }
00612         if((nbored & CFEmpty)==0) { newFlag |= CFNoNbEmpty; mNumInvIfCells++; }
00613 
00614         // calculate mass exchange for interface cells 
00615         LbmFloat myfrac = RAC(ccel,dFfrac);
00616         if(myfrac<0.) myfrac=0.; // NEWSURFT
00617 #       define nbdf(l) m[ this->dfInv[(l)] ]
00618 
00619         // update mass 
00620         // only do boundaries for fluid cells, and interface cells without
00621         // any fluid neighbors (assume that interface cells _with_ fluid
00622         // neighbors are affected enough by these) 
00623         // which Df's have to be reconstructed? 
00624         // for fluid cells - just the f_i difference from streaming to empty cells  ----
00625         numRecons = 0;
00626         bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
00627         //onlyBndnb = false; // DEBUG test off
00628 
00629         FORDF1 { // dfl loop
00630             recons[l] = 0;
00631             nbfracs[l] = 0.0;
00632             // finally, "normal" interface cells ----
00633             if( nbflag[l]&(CFFluid|CFBnd) ) { // NEWTEST! FIXME check!!!!!!!!!!!!!!!!!!
00634                 change = nbdf(l) - MYDF(l);
00635             }
00636             // interface cells - distuingish cells that shouldn't fill/empty 
00637             else if( nbflag[l] & CFInter ) {
00638                 
00639                 LbmFloat mynbfac,nbnbfac;
00640                 // NEW TEST t1
00641                 // t2 -> off
00642                 if((oldFlag&CFNoBndFluid)&&(nbflag[l]&CFNoBndFluid)) {
00643                     mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
00644                     nbnbfac = 1.0/mynbfac;
00645                     onlyBndnb = false;
00646                 } else {
00647                     mynbfac = nbnbfac = 1.0; // switch calc flux off
00648                     goto changeDefault;  // NEWSURFT
00649                     //change = 0.; goto changeDone;  // NEWSURFT
00650                 }
00651                 //mynbfac = nbnbfac = 1.0; // switch calc flux off t3
00652 
00653                 // perform interface case handling
00654                 if ((oldFlag|nbflag[l])&(CFNoNbFluid|CFNoNbEmpty)) {
00655                 switch (oldFlag&(CFNoNbFluid|CFNoNbEmpty)) {
00656                     case 0: 
00657                         // we are a normal cell so... 
00658                         switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
00659                             case CFNoNbFluid: 
00660                                 // just fill current cell = empty neighbor 
00661                                 change = nbnbfac*nbdf(l) ; goto changeDone; 
00662                             case CFNoNbEmpty: 
00663                                 // just empty current cell = fill neighbor 
00664                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
00665                         }
00666                         break;
00667 
00668                     case CFNoNbFluid: 
00669                         // we dont have fluid nb's so...
00670                         switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
00671                             case 0: 
00672                             case CFNoNbEmpty: 
00673                                 // we have no fluid nb's -> just empty
00674                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
00675                         }
00676                         break;
00677 
00678                     case CFNoNbEmpty: 
00679                         // we dont have empty nb's so...
00680                         switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
00681                             case 0: 
00682                             case CFNoNbFluid: 
00683                                 // we have no empty nb's -> just fill
00684                                 change = nbnbfac*nbdf(l); goto changeDone; 
00685                         }
00686                         break;
00687                 }} // inter-inter exchange
00688 
00689             changeDefault: ;
00690                 // just do normal mass exchange...
00691                 change = ( nbnbfac*nbdf(l) - mynbfac*MYDF(l) ) ;
00692             changeDone: ;
00693                 nbfracs[l] = QCELL_NB(lev, i,j,k, SRCS(lev),l, dFfrac);
00694                 if(nbfracs[l]<0.) nbfracs[l] = 0.; // NEWSURFT
00695                 change *=  (myfrac + nbfracs[l]) * 0.5;
00696             } // the other cell is interface
00697 
00698             // last alternative - reconstruction in this direction
00699             else {
00700                 // empty + bnd case
00701                 recons[l] = 1; 
00702                 numRecons++;
00703                 change = 0.0; 
00704             }
00705 
00706             // modify mass at SRCS
00707             mass += change;
00708         } // l
00709         // normal interface, no if empty/fluid
00710 
00711         // computenormal
00712         LbmFloat surfaceNormal[3];
00713         computeFluidSurfaceNormal(ccel,pFlagSrc, surfaceNormal);
00714 
00715         if( (ABS(surfaceNormal[0])+ABS(surfaceNormal[1])+ABS(surfaceNormal[2])) > LBM_EPSILON) {
00716             // normal ok and usable...
00717             FORDF1 {
00718                 if( (this->dfDvecX[l]*surfaceNormal[0] + this->dfDvecY[l]*surfaceNormal[1] + this->dfDvecZ[l]*surfaceNormal[2])  // dot Dvec,norml
00719                         > LBM_EPSILON) {
00720                     recons[l] = 2; 
00721                     numRecons++;
00722                 } 
00723             }
00724         }
00725 
00726         // calculate macroscopic cell values
00727         LbmFloat oldUx, oldUy, oldUz;
00728         LbmFloat oldRho; // OLD rho = ccel->rho;
00729 #       define REFERENCE_PRESSURE 1.0 // always atmosphere...
00730 #       if OPT3D==0
00731         oldRho=RAC(ccel,0);
00732         oldUx = oldUy = oldUz = 0.0;
00733         for(int l=1; l<this->cDfNum; l++) {
00734             oldRho += RAC(ccel,l);
00735             oldUx  += (this->dfDvecX[l]*RAC(ccel,l));
00736             oldUy  += (this->dfDvecY[l]*RAC(ccel,l)); 
00737             oldUz  += (this->dfDvecZ[l]*RAC(ccel,l)); 
00738         } 
00739         // reconstruct dist funcs from empty cells
00740         FORDF1 {
00741             if(recons[ l ]) {
00742                 m[ this->dfInv[l] ] = 
00743                     this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 
00744                     this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 
00745                     - MYDF( l );
00746             }
00747         }
00748         ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
00749         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
00750 #       else // OPT3D==0
00751         oldRho = + RAC(ccel,dC)  + RAC(ccel,dN )
00752                 + RAC(ccel,dS ) + RAC(ccel,dE )
00753                 + RAC(ccel,dW ) + RAC(ccel,dT )
00754                 + RAC(ccel,dB ) + RAC(ccel,dNE)
00755                 + RAC(ccel,dNW) + RAC(ccel,dSE)
00756                 + RAC(ccel,dSW) + RAC(ccel,dNT)
00757                 + RAC(ccel,dNB) + RAC(ccel,dST)
00758                 + RAC(ccel,dSB) + RAC(ccel,dET)
00759                 + RAC(ccel,dEB) + RAC(ccel,dWT)
00760                 + RAC(ccel,dWB);
00761 
00762         oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
00763                 + RAC(ccel,dNE) - RAC(ccel,dNW)
00764                 + RAC(ccel,dSE) - RAC(ccel,dSW)
00765                 + RAC(ccel,dET) + RAC(ccel,dEB)
00766                 - RAC(ccel,dWT) - RAC(ccel,dWB);
00767 
00768         oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
00769                 + RAC(ccel,dNE) + RAC(ccel,dNW)
00770                 - RAC(ccel,dSE) - RAC(ccel,dSW)
00771                 + RAC(ccel,dNT) + RAC(ccel,dNB)
00772                 - RAC(ccel,dST) - RAC(ccel,dSB);
00773 
00774         oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
00775                 + RAC(ccel,dNT) - RAC(ccel,dNB)
00776                 + RAC(ccel,dST) - RAC(ccel,dSB)
00777                 + RAC(ccel,dET) - RAC(ccel,dEB)
00778                 + RAC(ccel,dWT) - RAC(ccel,dWB);
00779 
00780         // now reconstruction
00781         ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
00782         rho = REFERENCE_PRESSURE;
00783         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
00784         if(recons[dN ]) { m[dS ] = EQN  + EQS  - MYDF(dN ); }
00785         if(recons[dS ]) { m[dN ] = EQS  + EQN  - MYDF(dS ); }
00786         if(recons[dE ]) { m[dW ] = EQE  + EQW  - MYDF(dE ); }
00787         if(recons[dW ]) { m[dE ] = EQW  + EQE  - MYDF(dW ); }
00788         if(recons[dT ]) { m[dB ] = EQT  + EQB  - MYDF(dT ); }
00789         if(recons[dB ]) { m[dT ] = EQB  + EQT  - MYDF(dB ); }
00790         if(recons[dNE]) { m[dSW] = EQNE + EQSW - MYDF(dNE); }
00791         if(recons[dNW]) { m[dSE] = EQNW + EQSE - MYDF(dNW); }
00792         if(recons[dSE]) { m[dNW] = EQSE + EQNW - MYDF(dSE); }
00793         if(recons[dSW]) { m[dNE] = EQSW + EQNE - MYDF(dSW); }
00794         if(recons[dNT]) { m[dSB] = EQNT + EQSB - MYDF(dNT); }
00795         if(recons[dNB]) { m[dST] = EQNB + EQST - MYDF(dNB); }
00796         if(recons[dST]) { m[dNB] = EQST + EQNB - MYDF(dST); }
00797         if(recons[dSB]) { m[dNT] = EQSB + EQNT - MYDF(dSB); }
00798         if(recons[dET]) { m[dWB] = EQET + EQWB - MYDF(dET); }
00799         if(recons[dEB]) { m[dWT] = EQEB + EQWT - MYDF(dEB); }
00800         if(recons[dWT]) { m[dEB] = EQWT + EQEB - MYDF(dWT); }
00801         if(recons[dWB]) { m[dET] = EQWB + EQET - MYDF(dWB); }
00802 #       endif       
00803 
00804 
00805         // inflow bc handling
00806         if(oldFlag & (CFMbndInflow)) {
00807             // fill if cells in inflow region
00808             if(myfrac<0.5) { 
00809                 mass += 0.25; 
00810                 mInitialMass += 0.25;
00811             }
00812             const int OId = oldFlag>>24;
00813             const LbmVec vel(mObjectSpeeds[OId]);
00814             ux=vel[0], uy=vel[1], uz=vel[2]; 
00815             //? usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
00816             //FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); } rho = fluidRho;
00817             rho = REFERENCE_PRESSURE;
00818             FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
00819             //errMsg("INFLOW_DEBUG","if at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
00820         }  else { 
00821             // NEWSURFT, todo optimize!
00822             if(onlyBndnb) { //if(usqr<0.001*0.001) {
00823                 rho = ux = uy = uz = 0.;
00824                 FORDF0{ 
00825                     rho += m[l]; 
00826                     ux  += (this->dfDvecX[l]*m[l]); 
00827                     uy  += (this->dfDvecY[l]*m[l]);  
00828                     uz  += (this->dfDvecZ[l]*m[l]);  
00829                 }
00830                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
00831             } else {// NEWSURFT */
00832                 if(usqr>0.3*0.3) { 
00833                     // force reset! , warning causes distortions...
00834                     FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,0.,0.,0.); }
00835                 } else {
00836                 // normal collide
00837                 // mass streaming done... do normal collide
00838                 LbmVec grav = mLevel[lev].gravity*mass;
00839                 DEFAULT_COLLIDEG(grav);
00840                 PERFORM_USQRMAXCHECK; }
00841                 // rho init from default collide necessary for fill/empty check below
00842             } // test
00843         }
00844 
00845         // testing..., particle generation
00846         // also check oldFlag for CFNoNbFluid, test
00847         // for inflow no pargen test // NOBUBBB!
00848         if((mInitDone) 
00849                 // dont allow new if cells, or submerged ones
00850                 && (!((oldFlag|newFlag)& (CFNoDelete|CFNoNbEmpty) )) 
00851                 // dont try to subtract from empty cells
00852                 && (mass>0.) && (mPartGenProb>0.0)) {
00853             bool doAdd = true;
00854             bool bndOk=true;
00855             if( (i<cutMin)||(i>mSizex-cutMin)||
00856                     (j<cutMin)||(j>mSizey-cutMin)||
00857                     (k<cutMin)||(k>mSizez-cutMin) ) { bndOk=false; }
00858             if(!bndOk) doAdd=false;
00859             
00860             LbmFloat prob = (rand()/(RAND_MAX+1.0));
00861             LbmFloat basethresh = mPartGenProb*lcsmqo*(LbmFloat)(mSizez+mSizey+mSizex)*0.5*0.333;
00862 
00863             // physical drop model
00864             if(mPartUsePhysModel) {
00865                 LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep);
00866                 LbmFloat rux = (ux * realWorldFac);
00867                 LbmFloat ruy = (uy * realWorldFac);
00868                 LbmFloat ruz = (uz * realWorldFac);
00869                 LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz));
00870                 basethresh *= rl;
00871 
00872                 // reduce probability in outer region?
00873                 const int pibord = mLevel[mMaxRefine].lSizex/2-cutConst;
00874                 const int pjbord = mLevel[mMaxRefine].lSizey/2-cutConst;
00875                 LbmFloat pifac = 1.-(LbmFloat)(ABS(i-pibord)) / (LbmFloat)(pibord);
00876                 LbmFloat pjfac = 1.-(LbmFloat)(ABS(j-pjbord)) / (LbmFloat)(pjbord);
00877                 if(pifac<0.) pifac=0.;
00878                 if(pjfac<0.) pjfac=0.;
00879 
00880                 //if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
00881                 if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
00882                     // add
00883                 } else {
00884                     doAdd = false; // dont...
00885                 }
00886 
00887                 // "wind" disturbance
00888                 // use realworld relative velocity here instead?
00889                 if( (doAdd && 
00890                         ((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) )// normal checks
00891                         ||(k>mSizez-SLOWDOWNREGION)   ) {
00892                     LbmFloat nuz = uz;
00893                     if(k>mSizez-SLOWDOWNREGION) {
00894                         // special case
00895                         LbmFloat zfac = (LbmFloat)( k-(mSizez-SLOWDOWNREGION) );
00896                         zfac /= (LbmFloat)(SLOWDOWNREGION);
00897                         nuz += (1.0) * zfac; // check max speed? OFF?
00898                         //errMsg("TOPT"," at "<<PRINT_IJK<<" zfac"<<zfac);
00899                     } else {
00900                         // normal probability
00901                         //? LbmFloat fac = P_LCSMQO-lcsmqo;
00902                         //? jdf *= fac;
00903                     }
00904                     FORDF1 {
00905                         const LbmFloat jdf = 0.05 * (rand()/(RAND_MAX+1.0));
00906                         // TODO  use wind velocity?
00907                         if(jdf>0.025) {
00908                         const LbmFloat add =  this->dfLength[l]*(-ux*this->dfDvecX[l]-uy*this->dfDvecY[l]-nuz*this->dfDvecZ[l])*jdf;
00909                         RAC(tcel,l) += add; }
00910                     }
00911                     //errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) );
00912                 } // wind disturbance
00913             } // mPartUsePhysModel
00914             else {
00915                 // empirical model
00916                 //if((prob<basethresh) && (lcsmqo>0.0095)) { // add
00917                 if((prob<basethresh) && (lcsmqo>0.012)) { // add
00918                 } else { doAdd = false; }// dont...
00919             } 
00920 
00921 
00922             // remove noise
00923             if(usqr<0.0001) doAdd=false;   // TODO check!?
00924 
00925             // dont try to subtract from empty cells
00926             // ensure cell has enough mass for new drop
00927             LbmFloat newPartsize = 1.0;
00928             if(mPartUsePhysModel) {
00929                 // 1-10
00930                 newPartsize += 9.0* (rand()/(RAND_MAX+1.0));
00931             } else {
00932                 // 1-5, overall size has to be less than
00933                 // .62 (ca. 0.5) to make drops significantly smaller 
00934                 // than a full cell!
00935                 newPartsize += 4.0* (rand()/(RAND_MAX+1.0));
00936             }
00937             LbmFloat dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize); //PARTMASS(mPartDropMassSub*newPartsize); // mass: 4/3 pi r^3 rho
00938             while(dropmass>mass) {
00939                 newPartsize -= 0.2;
00940                 dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize);
00941             }
00942             if(newPartsize<=1.) doAdd=false;
00943 
00944             if( (doAdd)  ) { // init new particle
00945                 for(int s=0; s<1; s++) { // one part!
00946                 const LbmFloat posjitter = 0.05;
00947                 const LbmFloat posjitteroffs = posjitter*-0.5;
00948                 LbmFloat jpx = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
00949                 LbmFloat jpy = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
00950                 LbmFloat jpz = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
00951 
00952                 const LbmFloat jitterstr = 1.0;
00953                 const LbmFloat jitteroffs = jitterstr*-0.5;
00954                 LbmFloat jx = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
00955                 LbmFloat jy = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
00956                 LbmFloat jz = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
00957 
00958                 // average normal & velocity 
00959                 // -> mostly along velocity dir, many into surface
00960                 // fluid velocity (not normalized!)
00961                 LbmVec flvelVel = LbmVec(ux,uy,uz);
00962                 LbmFloat flvelLen = norm(flvelVel);
00963                 // surface normal
00964                 LbmVec normVel = LbmVec(surfaceNormal[0],surfaceNormal[1],surfaceNormal[2]);
00965                 normalize(normVel);
00966                 LbmFloat normScale = (0.01+flvelLen);
00967                 // jitter vector, 0.2 * flvel
00968                 LbmVec jittVel = LbmVec(jx,jy,jz)*(0.05+flvelLen)*0.1;
00969                 // weighten velocities
00970                 const LbmFloat flvelWeight = 0.9;
00971                 LbmVec newpartVel = normVel*normScale*(1.-flvelWeight) + flvelVel*(flvelWeight) + jittVel; 
00972 
00973                 // offset towards surface (hide popping)
00974                 jpx += -normVel[0]*0.4;
00975                 jpy += -normVel[1]*0.4;
00976                 jpz += -normVel[2]*0.4;
00977 
00978                 LbmFloat srci=i+0.5+jpx, srcj=j+0.5+jpy, srck=k+0.5+jpz;
00979                 int type=0;
00980                 type = PART_DROP;
00981 
00982 #               if LBMDIM==2
00983                 newpartVel[2]=0.; srck=0.5;
00984 #               endif // LBMDIM==2
00985                 // subtract drop mass
00986                 mass -= dropmass;
00987                 // init new particle
00988                 {
00989                     ParticleObject np( ntlVec3Gfx(srci,srcj,srck) );
00990                     np.setVel(newpartVel[0],newpartVel[1],newpartVel[2]);
00991                     np.setStatus(PART_IN);
00992                     np.setType(type);
00993                     //if((s%3)==2) np.setType(PART_FLOAT);
00994                     np.setSize(newPartsize);
00995                     //errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<norm(LbmVec(ux,uy,uz)) <<" add"<<doAdd<<" pvel="<<norm(newpartVel)<<" size="<<newPartsize );
00996                     //errMsg("NEWPT","u="<<newpartVel<<" norm="<<normVel<<" flvel="<<flvelVel<<" jitt="<<jittVel );
00997                     FSGR_ADDPART(np);
00998                 } // new part
00999             //errMsg("PIT","NEW pit"<<np.getId()<<" pos:"<< np.getPos()<<" status:"<<convertFlags2String(np.getFlags())<<" vel:"<< np.getVel()  );
01000                 } // multiple parts
01001             } // doAdd
01002         } // */
01003 
01004 
01005         // interface cell filled or emptied?
01006         iffilled = ifemptied = 0;
01007         // interface cells empty/full?, WARNING: to mark these cells, better do it at the end of reinitCellFlags
01008         // interface cell if full?
01009         if( (mass) >= (rho * (1.0+FSGR_MAGICNR)) ) { iffilled = 1; }
01010         // interface cell if empty?
01011         if( (mass) <= (rho * (   -FSGR_MAGICNR)) ) { ifemptied = 1; }
01012 
01013         if(oldFlag & (CFMbndOutflow)) {
01014             mInitialMass -= mass;
01015             mass = myfrac = 0.0;
01016             iffilled = 0; ifemptied = 1;
01017         }
01018 
01019         // looks much nicer... LISTTRICK
01020 #       if FSGR_LISTTRICK==1
01021         //if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
01022         if(newFlag&CFNoBndFluid) { // test NEW TEST
01023             if(!iffilled) {
01024                 // remove cells independent from amount of change...
01025                 if( (oldFlag & CFNoNbEmpty)&&(newFlag & CFNoNbEmpty)&&
01026                         ( (mass>(rho*FSGR_LISTTTHRESHFULL))  || ((nbored&CFInter)==0)  )) { 
01027                     //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","filled "<<PRINT_IJK); };
01028                     iffilled = 1; 
01029                 } 
01030             }
01031             if(!ifemptied) {
01032                 if( (oldFlag & CFNoNbFluid)&&(newFlag & CFNoNbFluid)&&
01033                         ( (mass<(rho*FSGR_LISTTTHRESHEMPTY)) || ((nbored&CFInter)==0)  )) { 
01034                     //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","emptied "<<PRINT_IJK); };
01035                     ifemptied = 1; 
01036                 } 
01037             }
01038         } // nobndfluid only */
01039 #       endif
01040         //iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!!
01041         
01042 
01043         // now that all dfs are known, handle last changes
01044         if(iffilled) {
01045             LbmPoint filledp; filledp.flag=0;
01046             if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1;  // NEWSURFT
01047             filledp.x = i; filledp.y = j; filledp.z = k;
01048             LIST_FULL(filledp);
01049             //mNumFilledCells++; // DEBUG
01050             calcCellsFilled++;
01051         }
01052         else if(ifemptied) {
01053             LbmPoint emptyp; emptyp.flag=0;
01054             if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; //  NEWSURFT
01055             emptyp.x = i; emptyp.y = j; emptyp.z = k;
01056             LIST_EMPTY(emptyp);
01057             //mNumEmptiedCells++; // DEBUG
01058             calcCellsEmptied++;
01059         } 
01060         // dont cutoff values -> better cell conversions
01061         RAC(tcel,dFfrac)   = (mass/rho);
01062 
01063         // init new flux value
01064         float flux = FLUX_INIT; // dxqn on
01065         if(newFlag&CFNoBndFluid) {
01066             //flux = 50.0; // extreme on
01067             for(int nn=1; nn<this->cDfNum; nn++) { 
01068                 if(nbflag[nn] & (CFFluid|CFInter|CFBnd)) { flux += this->dfLength[nn]; }
01069             }
01070             // optical hack - smooth slow moving
01071             // surface regions
01072             if(usqr< sssUsqrLimit) {
01073             for(int nn=1; nn<this->cDfNum; nn++) { 
01074                 if(nbfracs[nn]!=0.0) {
01075                     LbmFloat curSmooth = (sssUsqrLimit-usqr)*sssUsqrLimitInv;
01076                     if(curSmooth>1.0) curSmooth=1.0;
01077                     flux *= (1.0+ smoothStrength*curSmooth * (nbfracs[nn]-myfrac)) ;
01078                 }
01079             } }
01080             // NEW TEST */
01081         }
01082         // flux = FLUX_INIT; // calc flux off
01083         QCELL(lev, i,j,k,TSET(lev), dFlux) = flux; // */
01084 
01085         // perform mass exchange with streamed values
01086         QCELL(lev, i,j,k,TSET(lev), dMass) = mass; // MASST
01087         // set new flag 
01088         *pFlagDst = (CellFlagType)newFlag;
01089         calcCurrentMass += mass; 
01090         calcCurrentVolume += RAC(tcel,dFfrac);
01091 
01092         // interface cell handling done...
01093 
01094 #if PARALLEL!=1
01095     GRID_LOOPREG_END();
01096 #else // PARALLEL==1
01097 #include "paraloopend.h" // = GRID_LOOPREG_END();
01098 #endif // PARALLEL==1
01099 
01100     // write vars from computations to class
01101     mLevel[lev].lmass    = calcCurrentMass;
01102     mLevel[lev].lvolume  = calcCurrentVolume;
01103     mNumFilledCells  = calcCellsFilled;
01104     mNumEmptiedCells = calcCellsEmptied;
01105     mNumUsedCells = calcNumUsedCells;
01106 }
01107 
01108 
01109 
01110 void 
01111 LbmFsgrSolver::preinitGrids()
01112 {
01113     const int lev = mMaxRefine;
01114     const bool doReduce = false;
01115     const int gridLoopBound=0;
01116 
01117     // preinit both grids
01118     for(int s=0; s<2; s++) {
01119     
01120         GRID_REGION_INIT();
01121 #if PARALLEL==1
01122 #pragma omp parallel default(shared) \
01123   reduction(+: \
01124       calcCurrentMass,calcCurrentVolume, \
01125         calcCellsFilled,calcCellsEmptied, \
01126         calcNumUsedCells )
01127 #endif // PARALLEL==1
01128         GRID_REGION_START();
01129         GRID_LOOP_START();
01130             for(int l=0; l<dTotalNum; l++) { RAC(ccel,l) = 0.; }
01131             *pFlagSrc =0;
01132             *pFlagDst =0;
01133             //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
01134 #if PARALLEL!=1
01135         GRID_LOOPREG_END();
01136 #else // PARALLEL==1
01137 #include "paraloopend.h" // = GRID_LOOPREG_END();
01138 #endif // PARALLEL==1
01139 
01140         /* dummy, remove warnings */ 
01141         calcCurrentMass = calcCurrentVolume = 0.;
01142         calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
01143         
01144         // change grid
01145         mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
01146         mLevel[mMaxRefine].setCurr   ^= 1;
01147     }
01148 }
01149 
01150 void 
01151 LbmFsgrSolver::standingFluidPreinit()
01152 {
01153     const int lev = mMaxRefine;
01154     const bool doReduce = false;
01155     const int gridLoopBound=1;
01156 
01157     GRID_REGION_INIT();
01158 #if PARALLEL==1
01159 #pragma omp parallel default(shared) \
01160   reduction(+: \
01161       calcCurrentMass,calcCurrentVolume, \
01162         calcCellsFilled,calcCellsEmptied, \
01163         calcNumUsedCells )
01164 #endif // PARALLEL==1
01165     GRID_REGION_START();
01166 
01167     LbmFloat rho, ux,uy,uz, usqr; 
01168     CellFlagType nbflag[LBM_DFNUM];
01169     LbmFloat m[LBM_DFNUM];
01170     LbmFloat lcsmqo;
01171 #   if OPT3D==1 
01172     LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
01173     CellFlagType nbored=0;
01174 #   endif // OPT3D==true 
01175 
01176     GRID_LOOP_START();
01177         //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
01178         const CellFlagType currFlag = *pFlagSrc; //RFLAG(lev, i,j,k,workSet);
01179         if( (currFlag & (CFEmpty|CFBnd)) ) continue;
01180 
01181         if( (currFlag & (CFInter)) ) {
01182             // copy all values
01183             for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
01184             continue;
01185         }
01186 
01187         if( (currFlag & CFNoBndFluid)) {
01188             OPTIMIZED_STREAMCOLLIDE;
01189         } else {
01190             FORDF1 {
01191                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
01192             } 
01193             DEFAULT_STREAM;
01194             DEFAULT_COLLIDEG(mLevel[lev].gravity);
01195         }
01196         for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
01197 #if PARALLEL!=1
01198     GRID_LOOPREG_END();
01199 #else // PARALLEL==1
01200 #include "paraloopend.h" // = GRID_LOOPREG_END();
01201 #endif // PARALLEL==1
01202 
01203     /* dummy remove warnings */ 
01204     calcCurrentMass = calcCurrentVolume = 0.;
01205     calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
01206     
01207     // change grid
01208   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
01209   mLevel[mMaxRefine].setCurr   ^= 1;
01210 }
01211 
01212 
01213 /******************************************************************************
01214  * work on lists from updateCellMass to reinit cell flags
01215  *****************************************************************************/
01216 
01217 LbmFloat LbmFsgrSolver::getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l) {
01218     //return 0.0; // test
01219     int level = mMaxRefine;
01220     LbmFloat     *ccel  = RACPNT(level, i,j,k, workSet);
01221 
01222     // computenormal
01223     CellFlagType *cflag = &RFLAG(level, i,j,k, workSet);
01224     LbmFloat n[3];
01225     computeFluidSurfaceNormal(ccel,cflag, n);
01226     LbmFloat scal = mDvecNrm[l][0]*n[0] + mDvecNrm[l][1]*n[1] + mDvecNrm[l][2]*n[2];
01227 
01228     LbmFloat ret = 1.0;
01229     // forward direction, add mass (for filling cells):
01230     if(dirForw) {
01231         if(scal<LBM_EPSILON) ret = 0.0;
01232         else ret = scal;
01233     } else {
01234         // backward for emptying
01235         if(scal>-LBM_EPSILON) ret = 0.0;
01236         else ret = scal * -1.0;
01237     }
01238     //errMsg("massd", PRINT_IJK<<" nv"<<nvel<<" : ret="<<ret ); //xit(1); //VECDEB
01239     return ret;
01240 }
01241 
01242 // warning - normal compuations are without
01243 //   boundary checks &
01244 //   normalization
01245 void LbmFsgrSolver::computeFluidSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt,LbmFloat *snret) {
01246     const int level = mMaxRefine;
01247     LbmFloat nx,ny,nz, nv1,nv2;
01248     const CellFlagType flagE = *(cflagpnt+1);
01249     const CellFlagType flagW = *(cflagpnt-1);
01250     if(flagE &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } 
01251     else if(flagE &(CFBnd)){ nv1 = 1.; }
01252     else nv1 = 0.0;
01253     if(flagW &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } 
01254     else if(flagW &(CFBnd)){ nv2 = 1.; }
01255     else nv2 = 0.0;
01256     nx = 0.5* (nv2-nv1);
01257 
01258     const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
01259     const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
01260     if(flagN &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 
01261     else if(flagN &(CFBnd)){ nv1 = 1.; }
01262     else nv1 = 0.0;
01263     if(flagS &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 
01264     else if(flagS &(CFBnd)){ nv2 = 1.; }
01265     else nv2 = 0.0;
01266     ny = 0.5* (nv2-nv1);
01267 
01268 #if LBMDIM==3
01269     const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
01270     const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
01271     if(flagT &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 
01272     else if(flagT &(CFBnd)){ nv1 = 1.; }
01273     else nv1 = 0.0;
01274     if(flagB &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 
01275     else if(flagB &(CFBnd)){ nv2 = 1.; }
01276     else nv2 = 0.0;
01277     nz = 0.5* (nv2-nv1);
01278 #else //LBMDIM==3
01279     nz = 0.0;
01280 #endif //LBMDIM==3
01281 
01282     // return vals
01283     snret[0]=nx; snret[1]=ny; snret[2]=nz;
01284 }
01285 void LbmFsgrSolver::computeFluidSurfaceNormalAcc(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
01286     LbmFloat nx=0.,ny=0.,nz=0.;
01287     ccel = NULL; cflagpnt=NULL; // remove warning
01288     snret[0]=nx; snret[1]=ny; snret[2]=nz;
01289 }
01290 void LbmFsgrSolver::computeObstacleSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
01291     const int level = mMaxRefine;
01292     LbmFloat nx,ny,nz, nv1,nv2;
01293     ccel = NULL; // remove warning
01294 
01295     const CellFlagType flagE = *(cflagpnt+1);
01296     const CellFlagType flagW = *(cflagpnt-1);
01297     if(flagE &(CFBnd)){ nv1 = 1.; }
01298     else nv1 = 0.0;
01299     if(flagW &(CFBnd)){ nv2 = 1.; }
01300     else nv2 = 0.0;
01301     nx = 0.5* (nv2-nv1);
01302 
01303     const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
01304     const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
01305     if(flagN &(CFBnd)){ nv1 = 1.; }
01306     else nv1 = 0.0;
01307     if(flagS &(CFBnd)){ nv2 = 1.; }
01308     else nv2 = 0.0;
01309     ny = 0.5* (nv2-nv1);
01310 
01311 #if LBMDIM==3
01312     const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
01313     const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
01314     if(flagT &(CFBnd)){ nv1 = 1.; }
01315     else nv1 = 0.0;
01316     if(flagB &(CFBnd)){ nv2 = 1.; }
01317     else nv2 = 0.0;
01318     nz = 0.5* (nv2-nv1);
01319 #else //LBMDIM==3
01320     nz = 0.0;
01321 #endif //LBMDIM==3
01322 
01323     // return vals
01324     snret[0]=nx; snret[1]=ny; snret[2]=nz;
01325 }
01326 void LbmFsgrSolver::computeObstacleSurfaceNormalAcc(int i,int j,int k, LbmFloat *snret) {
01327     bool nonorm = false;
01328     LbmFloat nx=0.,ny=0.,nz=0.;
01329     if(i<=0)        { nx =  1.; nonorm = true; }
01330     if(i>=mSizex-1) { nx = -1.; nonorm = true; }
01331     if(j<=0)        { ny =  1.; nonorm = true; }
01332     if(j>=mSizey-1) { ny = -1.; nonorm = true; }
01333 #   if LBMDIM==3
01334     if(k<=0)        { nz =  1.; nonorm = true; }
01335     if(k>=mSizez-1) { nz = -1.; nonorm = true; }
01336 #   endif // LBMDIM==3
01337     if(!nonorm) {
01338         // in domain, revert to helper, use setCurr&mMaxRefine
01339         LbmVec bnormal;
01340         CellFlagType *bflag = &RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr);
01341         LbmFloat     *bcell = RACPNT(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr);
01342         computeObstacleSurfaceNormal(bcell,bflag, &bnormal[0]);
01343         // TODO check if there is a normal near here?
01344         // use wider range otherwise...
01345         snret[0]=bnormal[0]; snret[1]=bnormal[0]; snret[2]=bnormal[0];
01346         return;
01347     }
01348     snret[0]=nx; snret[1]=ny; snret[2]=nz;
01349 }
01350 
01351 void LbmFsgrSolver::addToNewInterList( int ni, int nj, int nk ) {
01352 #if FSGR_STRICT_DEBUG==10
01353     // dangerous, this can change the simulation...
01354   /*for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
01355        iter != mListNewInter.end(); iter++ ) {
01356     if(ni!=iter->x) continue;
01357     if(nj!=iter->y) continue;
01358     if(nk!=iter->z) continue;
01359         // all 3 values match... skip point
01360         return;
01361     } */
01362 #endif // FSGR_STRICT_DEBUG==1
01363     // store point
01364     LbmPoint newinter; newinter.flag = 0;
01365     newinter.x = ni; newinter.y = nj; newinter.z = nk;
01366     mListNewInter.push_back(newinter);
01367 }
01368 
01369 void LbmFsgrSolver::reinitFlags( int workSet ) { 
01370     // reinitCellFlags OLD mods:
01371     // add all to intel list?
01372     // check ffrac for new cells
01373     // new if cell inits (last loop)
01374     // vweights handling
01375 
01376     const int debugFlagreinit = 0;
01377     
01378     // some things need to be read/modified on the other set
01379     int otherSet = (workSet^1);
01380     // fixed level on which to perform 
01381     int workLev = mMaxRefine;
01382 
01383   /* modify interface cells from lists */
01384   /* mark filled interface cells as fluid, or emptied as empty */
01385     /* count neighbors and distribute excess mass to interface neighbor cells
01386    * problems arise when there are no interface neighbors anymore
01387      * then just distribute to any fluid neighbors...
01388      */
01389 
01390     // for symmetry, first init all neighbor cells */
01391     for( vector<LbmPoint>::iterator iter=mListFull.begin();
01392        iter != mListFull.end(); iter++ ) {
01393     int i=iter->x, j=iter->y, k=iter->z;
01394         if(debugFlagreinit) errMsg("FULL", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
01395     FORDF1 {
01396             int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
01397             //if((LBMDIM>2)&&( (ni<=0) || (nj<=0) || (nk<=0) || (ni>=mLevel[workLev].lSizex-1) || (nj>=mLevel[workLev].lSizey-1) || (nk>=mLevel[workLev].lSizez-1) )) {
01398             if( (ni<=0) || (nj<=0) || 
01399                   (ni>=mLevel[workLev].lSizex-1) ||
01400                   (nj>=mLevel[workLev].lSizey-1) 
01401 #                   if LBMDIM==3
01402                   || (nk<=0) || (nk>=mLevel[workLev].lSizez-1) 
01403 #                   endif // LBMDIM==1
01404                  ) {
01405                 continue; } // new bc, dont treat cells on boundary NEWBC
01406       if( RFLAG(workLev, ni,nj,nk, workSet) & CFEmpty ){
01407                 
01408                 // preinit speed, get from average surrounding cells
01409                 // interpolate from non-workset to workset, sets are handled in function
01410 
01411                 // new and empty interface cell, dont change old flag here!
01412                 addToNewInterList(ni,nj,nk);
01413 
01414                 LbmFloat avgrho = 0.0;
01415                 LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
01416                 interpolateCellValues(workLev,ni,nj,nk,workSet, avgrho,avgux,avguy,avguz);
01417 
01418                 // careful with l's...
01419                 FORDF0M { 
01420                     QCELL(workLev,ni,nj,nk, workSet, m) = this->getCollideEq( m,avgrho,  avgux, avguy, avguz ); 
01421                     //QCELL(workLev,ni,nj,nk, workSet, l) = avgnbdf[l]; // CHECK FIXME test?
01422                 }
01423                 //errMsg("FNEW", PRINT_VEC(ni,nj,nk)<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<<avgrho<<" vel"<<PRINT_VEC(avgux,avguy,avguz) ); // DEBUG SYMM
01424                 QCELL(workLev,ni,nj,nk, workSet, dMass) = 0.0; //?? new
01425                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) = 0.0; //?? new
01426                 //RFLAG(workLev,ni,nj,nk,workSet) = (CellFlagType)(CFInter|CFNoInterpolSrc);
01427                 changeFlag(workLev,ni,nj,nk,workSet, (CFInter|CFNoInterpolSrc));
01428                 if(debugFlagreinit) errMsg("NEWE", PRINT_IJK<<" newif "<<PRINT_VEC(ni,nj,nk)<<" rho"<<avgrho<<" vel("<<avgux<<","<<avguy<<","<<avguz<<") " );
01429       }
01430             /* prevent surrounding interface cells from getting removed as empty cells 
01431              * (also cells that are not newly inited) */
01432       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
01433                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete);
01434                 changeFlag(workLev,ni,nj,nk, workSet, (RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete));
01435                 // also add to list...
01436                 addToNewInterList(ni,nj,nk);
01437             } // NEW?
01438     }
01439 
01440         // NEW? no extra loop...
01441         //RFLAG(workLev,i,j,k, workSet) = CFFluid;
01442         changeFlag(workLev,i,j,k, workSet,CFFluid);
01443     }
01444 
01445     /* remove empty interface cells that are not allowed to be removed anyway
01446      * this is important, otherwise the dreaded cell-type-flickering can occur! */
01447   for(int n=0; n<(int)mListEmpty.size(); n++) {
01448     int i=mListEmpty[n].x, j=mListEmpty[n].y, k=mListEmpty[n].z;
01449         if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) {
01450             // treat as "new inter"
01451             addToNewInterList(i,j,k);
01452             // remove entry
01453             if(debugFlagreinit) errMsg("EMPT REMOVED!!!", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
01454             if(n<(int)mListEmpty.size()-1) mListEmpty[n] = mListEmpty[mListEmpty.size()-1]; 
01455             mListEmpty.pop_back();
01456             n--; 
01457         }
01458     } 
01459 
01460 
01461     /* problems arise when adjacent cells empty&fill ->
01462          let fill cells+surrounding interface cells have the higher importance */
01463   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
01464        iter != mListEmpty.end(); iter++ ) {
01465     int i=iter->x, j=iter->y, k=iter->z;
01466         if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)){ errMsg("A"," ARGHARGRAG "); } // DEBUG
01467         if(debugFlagreinit) errMsg("EMPT", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) );
01468 
01469         /* set surrounding fluid cells to interface cells */
01470     FORDF1 {
01471             int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
01472       if( RFLAG(workLev,ni,nj,nk, workSet) & CFFluid){
01473                 // init fluid->interface 
01474                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(CFInter); 
01475                 changeFlag(workLev,ni,nj,nk, workSet, CFInter); 
01476                 /* new mass = current density */
01477                 LbmFloat nbrho = QCELL(workLev,ni,nj,nk, workSet, dC);
01478             for(int rl=1; rl< this->cDfNum ; ++rl) { nbrho += QCELL(workLev,ni,nj,nk, workSet, rl); }
01479                 QCELL(workLev,ni,nj,nk, workSet, dMass) =  nbrho; 
01480                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) =  1.0; 
01481 
01482                 // store point
01483                 addToNewInterList(ni,nj,nk);
01484       }
01485       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter){
01486                 // test, also add to list...
01487                 addToNewInterList(ni,nj,nk);
01488             } // NEW?
01489     }
01490 
01491         /* for symmetry, set our flag right now */
01492         changeFlag(workLev,i,j,k, workSet, CFEmpty);
01493         // mark cell not be changed mass... - not necessary, not in list anymore anyway!
01494     } // emptylist
01495 
01496 
01497     
01498     // precompute weights to get rid of order dependancies
01499     vector<lbmFloatSet> vWeights;
01500     vWeights.resize( mListFull.size() + mListEmpty.size() );
01501     int weightIndex = 0;
01502   int nbCount = 0;
01503     LbmFloat nbWeights[LBM_DFNUM];
01504     LbmFloat nbTotWeights = 0.0;
01505     for( vector<LbmPoint>::iterator iter=mListFull.begin();
01506        iter != mListFull.end(); iter++ ) {
01507     int i=iter->x, j=iter->y, k=iter->z;
01508     nbCount = 0; nbTotWeights = 0.0;
01509     FORDF1 {
01510             int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
01511       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
01512                 nbCount++;
01513                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
01514                 else nbWeights[l] = getMassdWeight(1,i,j,k,workSet,l); // NEWSURFT
01515                 nbTotWeights += nbWeights[l];
01516       } else {
01517                 nbWeights[l] = -100.0; // DEBUG;
01518             }
01519     }
01520         if(nbCount>0) { 
01521             //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
01522         vWeights[weightIndex].val[0] = nbTotWeights;
01523         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
01524         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
01525         } else { 
01526         vWeights[weightIndex].numNbs = 0.0;
01527         }
01528         weightIndex++;
01529     }
01530   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
01531        iter != mListEmpty.end(); iter++ ) {
01532     int i=iter->x, j=iter->y, k=iter->z;
01533     nbCount = 0; nbTotWeights = 0.0;
01534     FORDF1 {
01535             int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
01536       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
01537                 nbCount++;
01538                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
01539                 else nbWeights[l] = getMassdWeight(0,i,j,k,workSet,l); // NEWSURFT
01540                 nbTotWeights += nbWeights[l];
01541       } else {
01542                 nbWeights[l] = -100.0; // DEBUG;
01543             }
01544     }
01545         if(nbCount>0) { 
01546             //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
01547         vWeights[weightIndex].val[0] = nbTotWeights;
01548         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
01549         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
01550         } else { 
01551         vWeights[weightIndex].numNbs = 0.0;
01552         }
01553         weightIndex++;
01554     } 
01555     weightIndex = 0;
01556     
01557 
01558     /* process full list entries, filled cells are done after this loop */
01559     for( vector<LbmPoint>::iterator iter=mListFull.begin();
01560        iter != mListFull.end(); iter++ ) {
01561     int i=iter->x, j=iter->y, k=iter->z;
01562 
01563         LbmFloat myrho = QCELL(workLev,i,j,k, workSet, dC);
01564     FORDF1 { myrho += QCELL(workLev,i,j,k, workSet, l); } // QCELL.rho
01565 
01566     LbmFloat massChange = QCELL(workLev,i,j,k, workSet, dMass) - myrho;
01567         if(vWeights[weightIndex].numNbs>0.0) {
01568             const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
01569             //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
01570             FORDF1 {
01571                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
01572         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
01573                     LbmFloat change = -1.0;
01574                     if(nbTotWeightsp>0.0) {
01575                         //change = massChange * ( nbWeights[l]/nbTotWeightsp );
01576                         change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
01577                     } else {
01578                         change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
01579                     }
01580                     QCELL(workLev,ni,nj,nk, workSet, dMass) += change;
01581                 }
01582             }
01583             massChange = 0.0;
01584         } else {
01585             // Problem! no interface neighbors
01586             mFixMass += massChange;
01587             //TTT mNumProblems++;
01588             //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<<mFixMass);
01589         }
01590         weightIndex++;
01591 
01592     // already done? RFLAG(workLev,i,j,k, workSet) = CFFluid;
01593     QCELL(workLev,i,j,k, workSet, dMass) = myrho; // should be rho... but unused?
01594     QCELL(workLev,i,j,k, workSet, dFfrac) = 1.0; // should be rho... but unused?
01595   } // fulllist
01596 
01597 
01598     /* now, finally handle the empty cells - order is important, has to be after
01599      * full cell handling */
01600   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
01601        iter != mListEmpty.end(); iter++ ) {
01602     int i=iter->x, j=iter->y, k=iter->z;
01603 
01604     LbmFloat massChange = QCELL(workLev, i,j,k, workSet, dMass);
01605         if(vWeights[weightIndex].numNbs>0.0) {
01606             const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
01607             //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
01608             FORDF1 {
01609                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
01610         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
01611                     LbmFloat change = -1.0;
01612                     if(nbTotWeightsp>0.0) {
01613                         change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
01614                     } else {
01615                         change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
01616                     }
01617                     QCELL(workLev, ni,nj,nk, workSet, dMass) += change;
01618                 }
01619             }
01620             massChange = 0.0;
01621         } else {
01622             // Problem! no interface neighbors
01623             mFixMass += massChange;
01624             //TTT mNumProblems++;
01625             //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<mFixMass);
01626         }
01627         weightIndex++;
01628         
01629         // finally... make it empty 
01630     // already done? RFLAG(workLev,i,j,k, workSet) = CFEmpty;
01631     QCELL(workLev,i,j,k, workSet, dMass) = 0.0;
01632     QCELL(workLev,i,j,k, workSet, dFfrac) = 0.0;
01633     }
01634   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
01635        iter != mListEmpty.end(); iter++ ) {
01636     int i=iter->x, j=iter->y, k=iter->z;
01637     changeFlag(workLev,i,j,k, otherSet, CFEmpty);
01638     } 
01639 
01640 
01641     // check if some of the new interface cells can be removed again 
01642     // never happens !!! not necessary
01643     // calculate ffrac for new IF cells NEW
01644 
01645     // how many are really new interface cells?
01646     int numNewIf = 0;
01647   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
01648        iter != mListNewInter.end(); iter++ ) {
01649     int i=iter->x, j=iter->y, k=iter->z;
01650         if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
01651             continue; 
01652             // FIXME remove from list?
01653         }
01654         numNewIf++;
01655     }
01656 
01657     // redistribute mass, reinit flags
01658     if(debugFlagreinit) errMsg("NEWIF", "total:"<<mListNewInter.size());
01659     float newIfFac = 1.0/(LbmFloat)numNewIf;
01660   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
01661        iter != mListNewInter.end(); iter++ ) {
01662     int i=iter->x, j=iter->y, k=iter->z;
01663         if((i<=0) || (j<=0) || 
01664              (i>=mLevel[workLev].lSizex-1) ||
01665              (j>=mLevel[workLev].lSizey-1) ||
01666              ((LBMDIM==3) && ((k<=0) || (k>=mLevel[workLev].lSizez-1) ) )
01667              ) {
01668             continue; } // new bc, dont treat cells on boundary NEWBC
01669         if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
01670             //errMsg("???"," "<<PRINT_IJK);
01671             continue; 
01672         } // */
01673 
01674     QCELL(workLev,i,j,k, workSet, dMass) += (mFixMass * newIfFac);
01675         
01676         int nbored = 0;
01677         FORDF1 { nbored |= RFLAG_NB(workLev, i,j,k, workSet,l); }
01678         if(!(nbored & CFBndNoslip)) { RFLAG(workLev,i,j,k, workSet) |= CFNoBndFluid; }
01679         if(!(nbored & CFFluid))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbFluid; }
01680         if(!(nbored & CFEmpty))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbEmpty; }
01681 
01682         if(!(RFLAG(workLev,i,j,k, otherSet)&CFInter)) {
01683             RFLAG(workLev,i,j,k, workSet) = (CellFlagType)(RFLAG(workLev,i,j,k, workSet) | CFNoDelete);
01684         }
01685         if(debugFlagreinit) errMsg("NEWIF", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" f"<< convertCellFlagType2String(RFLAG(workLev,i,j,k, workSet))<<" wl"<<workLev );
01686     }
01687 
01688     // reinit fill fraction
01689   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
01690        iter != mListNewInter.end(); iter++ ) {
01691     int i=iter->x, j=iter->y, k=iter->z;
01692         if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { continue; }
01693 
01694         initInterfaceVars(workLev, i,j,k, workSet, false); //int level, int i,int j,int k,int workSet, bool initMass) {
01695         //LbmFloat nrho = 0.0;
01696         //FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); }
01697     //QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho;
01698     //QCELL(workLev,i,j,k, workSet, dFlux) = FLUX_INIT;
01699     }
01700 
01701     if(mListNewInter.size()>0){ 
01702         //errMsg("FixMassDisted"," fm:"<<mFixMass<<" nif:"<<mListNewInter.size() );
01703         mFixMass = 0.0; 
01704     }
01705 
01706     // empty lists for next step
01707     mListFull.clear();
01708     mListEmpty.clear();
01709     mListNewInter.clear();
01710 } // reinitFlags
01711 
01712 
01713