2888c7e172ea3f7b2a825a6183dfbde8b1ed1594
[blender.git] / intern / elbeem / intern / solver_main.cpp
1 /******************************************************************************
2  *
3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4  * Copyright 2003,2004,2005 Nils Thuerey
5  *
6  * Standard LBM Factory implementation
7  *
8  *****************************************************************************/
9
10 #include "solver_class.h"
11 #include "solver_relax.h"
12 #include "particletracer.h"
13
14 /*****************************************************************************/
15 /*! perform a single LBM step */
16 /*****************************************************************************/
17
18
19 void LbmFsgrSolver::step() { 
20         initLevelOmegas();
21         stepMain(); 
22 }
23
24 void LbmFsgrSolver::stepMain()
25 {
26         this->markedClearList(); // DMC clearMarkedCellsList
27
28         // safety check, counter reset
29         this->mNumUsedCells = 0;
30         mNumInterdCells = 0;
31         mNumInvIfCells = 0;
32
33   //debugOutNnl("LbmFsgrSolver::step : "<<this->mStepCnt, 10);
34   if(!this->mSilent){ debMsgStd("LbmFsgrSolver::step", DM_MSG, this->mName<<" cnt:"<<this->mStepCnt<<" t:"<<mSimulationTime, 10); }
35         //debMsgDirect(  "LbmFsgrSolver::step : "<<this->mStepCnt<<" ");
36         myTime_t timestart = getTime();
37         //myTime_t timestart = 0;
38         //if(mStartSymm) { checkSymmetry("step1"); } // DEBUG 
39
40         // time adapt
41         mMaxVlen = mMxvz = mMxvy = mMxvx = 0.0;
42
43         // init moving bc's, can change mMaxVlen
44         initMovingObstacles(false);
45 #if LBM_INCLUDE_TESTSOLVERS==1
46         handleCpdata();
47 #endif
48
49         // important - keep for tadap
50         LbmFloat lastMass = mCurrentMass;
51         mCurrentMass = this->mFixMass; // reset here for next step
52         mCurrentVolume = 0.0;
53         
54         //change to single step advance!
55         int levsteps = 0;
56         int dsbits = this->mStepCnt ^ (this->mStepCnt-1);
57         //errMsg("S"," step:"<<this->mStepCnt<<" s-1:"<<(this->mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits));
58         for(int lev=0; lev<=mMaxRefine; lev++) {
59                 //if(! (this->mStepCnt&(1<<lev)) ) {
60                 if( dsbits & (1<<(mMaxRefine-lev)) ) {
61                         //errMsg("S"," l"<<lev);
62
63                         if(lev==mMaxRefine) {
64                                 // always advance fine level...
65                                 fineAdvance();
66                         } else {
67                                 adaptGrid(lev);
68                                 coarseRestrictFromFine(lev);
69                                 coarseAdvance(lev);
70                         }
71 #if FSGR_OMEGA_DEBUG==1
72                         errMsg("LbmFsgrSolver::step","LES stats l="<<lev<<" omega="<<mLevel[lev].omega<<" avgOmega="<< (mLevel[lev].avgOmega/mLevel[lev].avgOmegaCnt) );
73                         mLevel[lev].avgOmega = 0.0; mLevel[lev].avgOmegaCnt = 0.0;
74 #endif // FSGR_OMEGA_DEBUG==1
75                         levsteps++;
76                 }
77                 mCurrentMass   += mLevel[lev].lmass;
78                 mCurrentVolume += mLevel[lev].lvolume;
79         }
80
81   // prepare next step
82         this->mStepCnt++;
83
84
85         // some dbugging output follows
86         // calculate MLSUPS
87         myTime_t timeend = getTime();
88
89         this->mNumUsedCells += mNumInterdCells; // count both types for MLSUPS
90         mAvgNumUsedCells += this->mNumUsedCells;
91         this->mMLSUPS = (this->mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000);
92         if(this->mMLSUPS>10000){ this->mMLSUPS = -1; }
93         else { mAvgMLSUPS += this->mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
94         
95         LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1())) / ((timeend-timestart)/(double)1000.0) ) / (1000000);
96         if(totMLSUPS>10000) totMLSUPS = -1;
97         mNumInvIfTotal += mNumInvIfCells; // debug
98
99   // do some formatting 
100   if(!this->mSilent){ 
101                 string sepStr(""); // DEBUG
102                 int avgcls = (int)(mAvgNumUsedCells/(LONGINT)this->mStepCnt);
103         debMsgStd("LbmFsgrSolver::step", DM_MSG, this->mName<<" cnt:"<<this->mStepCnt<<" t:"<<mSimulationTime<<
104                 //debMsgDirect( 
105                         " mlsups(curr:"<<this->mMLSUPS<<
106                         " avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< sepStr<<
107                         " totcls:"<<(this->mNumUsedCells)<< sepStr<<
108                         " avgcls:"<< avgcls<< sepStr<<
109                         " intd:"<<mNumInterdCells<< sepStr<<
110                         " invif:"<<mNumInvIfCells<< sepStr<<
111                         " invift:"<<mNumInvIfTotal<< sepStr<<
112                         " fsgrcs:"<<mNumFsgrChanges<< sepStr<<
113                         " filled:"<<this->mNumFilledCells<<", emptied:"<<this->mNumEmptiedCells<< sepStr<<
114                         " mMxv:"<<mMxvx<<","<<mMxvy<<","<<mMxvz<<", tscnts:"<<mTimeSwitchCounts<< sepStr<<
115                         " RWmxv:"<<ntlVec3Gfx(mMxvx,mMxvy,mMxvz)*(mLevel[mMaxRefine].simCellSize / mLevel[mMaxRefine].timestep)<<" "<< /* realworld vel output */
116                         " probs:"<<mNumProblems<< sepStr<<
117                         " simt:"<<mSimulationTime<< sepStr<<
118                         " for '"<<this->mName<<"' " , 10);
119         } else { debMsgDirect("."); }
120
121         if(this->mStepCnt==1) {
122                 mMinNoCells = mMaxNoCells = this->mNumUsedCells;
123         } else {
124                 if(this->mNumUsedCells>mMaxNoCells) mMaxNoCells = this->mNumUsedCells;
125                 if(this->mNumUsedCells<mMinNoCells) mMinNoCells = this->mNumUsedCells;
126         }
127         
128         // mass scale test
129         if((mMaxRefine>0)&&(mInitialMass>0.0)) {
130                 LbmFloat mscale = mInitialMass/mCurrentMass;
131
132                 mscale = 1.0;
133                 const LbmFloat dchh = 0.001;
134                 if(mCurrentMass<mInitialMass) mscale = 1.0+dchh;
135                 if(mCurrentMass>mInitialMass) mscale = 1.0-dchh;
136
137                 // use mass rescaling?
138                 // with float precision this seems to be nonsense...
139                 const bool MREnable = false;
140
141                 const int MSInter = 2;
142                 static int mscount = 0;
143                 if( (MREnable) && ((mLevel[0].lsteps%MSInter)== (MSInter-1)) && ( ABS( (mInitialMass/mCurrentMass)-1.0 ) > 0.01) && ( dsbits & (1<<(mMaxRefine-0)) ) ){
144                         // example: FORCE RESCALE MASS! ini:1843.5, cur:1817.6, f=1.01425 step:22153 levstep:5539 msc:37
145                         // mass rescale MASS RESCALE check
146                         errMsg("MDTDD","\n\n");
147                         errMsg("MDTDD","FORCE RESCALE MASS! "
148                                         <<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass)
149                                         <<" step:"<<this->mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
150                                         );
151                         errMsg("MDTDD","\n\n");
152
153                         mscount++;
154                         for(int lev=mMaxRefine; lev>=0 ; lev--) {
155                                 //for(int workSet = 0; workSet<=1; workSet++) {
156                                 int wss = 0;
157                                 int wse = 1;
158 #if COMPRESSGRIDS==1
159                                 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
160 #endif // COMPRESSGRIDS==1
161                                 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
162
163                                         FSGR_FORIJK1(lev) {
164                                                 if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) 
165                                                         ) {
166
167                                                         FORDF0 { QCELL(lev, i,j,k,workSet, l) *= mscale; }
168                                                         QCELL(lev, i,j,k,workSet, dMass) *= mscale;
169                                                         QCELL(lev, i,j,k,workSet, dFfrac) *= mscale;
170
171                                                 } else {
172                                                         continue;
173                                                 }
174                                         }
175                                 }
176                                 mLevel[lev].lmass *= mscale;
177                         }
178                 } 
179
180                 mCurrentMass *= mscale;
181         }// if mass scale test */
182         else {
183                 // use current mass after full step for initial setting
184                 if((mMaxRefine>0)&&(mInitialMass<=0.0) && (levsteps == (mMaxRefine+1))) {
185                         mInitialMass = mCurrentMass;
186                         debMsgStd("MDTDD",DM_NOTIFY,"Second Initial Mass Init: "<<mInitialMass, 2);
187                 }
188         }
189
190 #if LBM_INCLUDE_TESTSOLVERS==1
191         if((mUseTestdata)&&(this->mInitDone)) { handleTestdata(); }
192 #endif
193
194         // advance positions with current grid
195         advanceParticles();
196         if(mpParticles) {
197                 mpParticles->checkDumpTextPositions(mSimulationTime);
198                 mpParticles->checkTrails(mSimulationTime);
199         }
200
201         // one of the last things to do - adapt timestep
202         // was in fineAdvance before... 
203         if(mTimeAdap) {
204                 adaptTimestep();
205         } // time adaptivity
206
207
208 #ifndef WIN32
209         // good indicator for instabilities
210         if( (!finite(mMxvx)) || (!finite(mMxvy)) || (!finite(mMxvz)) ) { CAUSE_PANIC; }
211         if( (!finite(mCurrentMass)) || (!finite(mCurrentVolume)) ) { CAUSE_PANIC; }
212 #endif // WIN32
213
214         // output total step time
215         timeend = getTime();
216         debMsgStd("LbmFsgrSolver::stepMain",DM_MSG,"step:"<<this->mStepCnt
217                         <<": dccd="<< mCurrentMass<<",d"<<(lastMass-mCurrentMass)<<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<"), "
218                         <<" totst:"<<getTimeString(timeend-timestart), 3);
219         // nicer output
220         debMsgDirect(std::endl); 
221
222  //#endif // ELBEEM_PLUGIN!=1
223 }
224
225 #define NEWDEBCHECK(str) \
226         if(!this->mPanic){ FSGR_FORIJK_BOUNDS(mMaxRefine) { \
227                 if(RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr)&(CFFluid|CFInter)) { \
228                 for(int l=0;l<dTotalNum;l++) { \
229                         if(!finite(QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr,l))) { errMsg("NNOFIN"," "<<str<<" at "<<PRINT_IJK<<" l"<<l<<" "); }\
230                 }/*for*/ \
231                 }/*if*/ \
232         } }
233
234 void LbmFsgrSolver::fineAdvance()
235 {
236         // do the real thing...
237         //NEWDEBCHECK("t1");
238         mainLoop( mMaxRefine );
239         if(mUpdateFVHeight) {
240                 // warning assume -Y gravity...
241                 mFVHeight = mCurrentMass*mFVArea/((LbmFloat)(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizez));
242                 if(mFVHeight<1.0) mFVHeight = 1.0;
243                 this->mpParam->setFluidVolumeHeight(mFVHeight);
244         }
245
246         // advance time before timestep change
247         mSimulationTime += this->mpParam->getTimestep();
248         // time adaptivity
249         this->mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) );
250         //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG 
251         if(!this->mSilent){ errMsg("fineAdvance"," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps) ); }
252
253         // update other set
254   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
255   mLevel[mMaxRefine].setCurr   ^= 1;
256   mLevel[mMaxRefine].lsteps++;
257
258         // flag init... (work on current set, to simplify flag checks)
259         reinitFlags( mLevel[mMaxRefine].setCurr );
260         if(!this->mSilent){ errMsg("fineAdvance"," flags reinit on set "<< mLevel[mMaxRefine].setCurr ); }
261 }
262
263
264 /*****************************************************************************/
265 //! fine step function
266 /*****************************************************************************/
267
268
269 // access to own dfs during step (may be changed to local array)
270 #define MYDF(l) RAC(ccel, l)
271
272 void 
273 LbmFsgrSolver::mainLoop(int lev)
274 {
275         // loops over _only inner_ cells  -----------------------------------------------------------------------------------
276         LbmFloat calcCurrentMass = 0.0;
277         LbmFloat calcCurrentVolume = 0.0;
278         int      calcCellsFilled = this->mNumFilledCells;
279         int      calcCellsEmptied = this->mNumEmptiedCells;
280         int      calcNumUsedCells = this->mNumUsedCells;
281         const int cutMin  = 1;
282         const int cutConst = mCutoff+2;
283
284 #       if LBM_INCLUDE_TESTSOLVERS==1
285         // 3d region off... quit
286         if((mUseTestdata)&&(mpTest->mDebugvalue1>0.0)) { return; }
287 #endif // ELBEEM_PLUGIN!=1
288         //printLbmCell(lev, 6,6,16, mLevel[lev].setCurr ); // DEBUG
289         
290 #if PARALLEL==1
291 #include "paraloop.h"
292 #else // PARALLEL==1
293   { // main loop region
294         int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine);
295         //{ errMsg("LbmFsgrSolver::mainLoop","Err MAINADVANCE0 ks:"<< kstart<<" ke:"<<kend<<" dim:"<<LBMDIMcDimension<<" mlsz:"<< mLevel[mMaxRefine].lSizez<<" zmax1:"<<getForZMax1(mMaxRefine) ); } // DEBUG
296 #define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz);
297 #endif // PARALLEL==1
298
299         // local to loop
300         CellFlagType nbflag[LBM_DFNUM]; 
301         LbmFloat *ccel = NULL;
302         LbmFloat *tcel = NULL;
303         int oldFlag;
304         int newFlag;
305         int nbored;
306         LbmFloat m[LBM_DFNUM];
307         LbmFloat rho, ux, uy, uz, tmp, usqr;
308         LbmFloat mass, change, lcsmqo=0.0;
309         usqr = tmp = 0.0; 
310 #if OPT3D==1 
311         LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
312 #endif // OPT3D==true 
313
314         
315         // ifempty cell conversion flags
316         bool iffilled, ifemptied;
317         LbmFloat nbfracs[LBM_DFNUM]; // ffracs of neighbors
318         int recons[LBM_DFNUM];   // reconstruct this DF?
319         int numRecons;           // how many are reconstructed?
320
321         // slow surf regions smooth (if below)
322         const LbmFloat smoothStrength = 0.0; //0.01;
323         const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03;
324         const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit;
325         
326         CellFlagType *pFlagSrc;
327         CellFlagType *pFlagDst;
328         pFlagSrc = &RFLAG(lev, 0,1, kstart,SRCS(lev)); // omp
329         pFlagDst = &RFLAG(lev, 0,1, kstart,TSET(lev)); // omp
330         ccel = RACPNT(lev, 0,1, kstart ,SRCS(lev)); // omp
331         tcel = RACPNT(lev, 0,1, kstart ,TSET(lev)); // omp
332         //CellFlagType *pFlagTar = NULL;
333         int pFlagTarOff;
334         if(mLevel[lev].setOther==1) pFlagTarOff = mLevel[lev].lOffsz;
335         else pFlagTarOff = -mLevel[lev].lOffsz;
336 #define ADVANCE_POINTERS(p)     \
337         ccel += (QCELLSTEP*(p));        \
338         tcel += (QCELLSTEP*(p));        \
339         pFlagSrc+= (p); \
340         pFlagDst+= (p); \
341         i+= (p);
342
343         // ---
344         // now stream etc.
345
346         //{ errMsg("LbmFsgrSolver::mainLoop","Err MAINADVANCE0 ks:"<< kstart<<" ke:"<<kend<<" dim:"<<LBMDIM<<" mlsz:"<<mLevel[mMaxRefine].lSizez ); } // DEBUG
347
348         // use //template functions for 2D/3D
349 #if COMPRESSGRIDS==0
350   for(int k=kstart;k<kend;++k) {
351   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
352   for(int i=0;i<mLevel[lev].lSizex-2;   ) {
353 #else // COMPRESSGRIDS==0
354         int kdir = 1; // COMPRT ON
355         if(mLevel[mMaxRefine].setCurr==1) {
356                 kdir = -1;
357                 int temp = kend;
358                 kend = kstart-1;
359                 kstart = temp-1;
360         } // COMPRT
361
362 #if PARALLEL==0
363   const int id = 0, Nthrds = 1;
364 #endif // PARALLEL==1
365   const int Nj = mLevel[mMaxRefine].lSizey;
366         int jstart = 0+( id * (Nj / Nthrds) );
367         int jend   = 0+( (id+1) * (Nj / Nthrds) );
368   if( ((Nj/Nthrds) *Nthrds) != Nj) {
369     errMsg("LbmFsgrSolver","Invalid domain size Nj="<<Nj<<" Nthrds="<<Nthrds);
370   }
371         // cutoff obstacle boundary
372         if(jstart<1) jstart = 1;
373         if(jend>mLevel[mMaxRefine].lSizey-1) jend = mLevel[mMaxRefine].lSizey-1;
374
375 #if PARALLEL==1
376         PARA_INITIALIZE();
377         //errMsg("LbmFsgrSolver::mainLoop","id="<<id<<" js="<<jstart<<" je="<<jend<<" jdir="<<(1) ); // debug
378 #endif // PARALLEL==1
379   for(int k=kstart;k!=kend;k+=kdir) {
380
381         //errMsg("LbmFsgrSolver::mainLoop","k="<<k<<" ks="<<kstart<<" ke="<<kend<<" kdir="<<kdir<<" x*y="<<mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum ); // debug
382   pFlagSrc = &RFLAG(lev, 0, jstart, k, SRCS(lev)); // omp test // COMPRT ON
383   pFlagDst = &RFLAG(lev, 0, jstart, k, TSET(lev)); // omp test
384   ccel = RACPNT(lev,     0, jstart, k, SRCS(lev)); // omp test
385   tcel = RACPNT(lev,     0, jstart, k, TSET(lev)); // omp test // COMPRT ON
386
387   //for(int j=1;j<mLevel[lev].lSizey-1;++j) {
388   for(int j=jstart;j!=jend;++j) {
389   for(int i=0;i<mLevel[lev].lSizex-2;   ) {
390 #endif // COMPRESSGRIDS==0
391
392                 ADVANCE_POINTERS(1); 
393 #if FSGR_STRICT_DEBUG==1
394                 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
395                 if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) || 
396                     (&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) {
397                         errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<<
398                                         RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<<
399                                         (*pFlagSrc)<<","<<(*pFlagDst) <<",  pointers "<<
400           (int)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(int)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
401           (int)(pFlagSrc)<<","<<(int)(pFlagDst)<<" "
402                                         ); 
403                         CAUSE_PANIC;
404                 }       
405                 if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) || 
406                     (&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) {
407                         errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<<
408           (int)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(int)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
409           (int)(ccel)<<","<<(int)(tcel)<<" "
410                                         ); 
411                         CAUSE_PANIC;
412                 }       
413 #endif
414                 oldFlag = *pFlagSrc;
415                 //DEBUG if(LBMDIM==2) { errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<< RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<" "); }
416                 // stream from current set to other, then collide and store
417                 
418                 // old INTCFCOARSETEST==1
419                 if( (oldFlag & (CFGrFromCoarse)) ) {  // interpolateFineFromCoarse test!
420                         if(( this->mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
421                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
422                         } else {
423                                 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
424                                 calcNumUsedCells++;
425                         }
426                         continue; // interpolateFineFromCoarse test!
427                 } // interpolateFineFromCoarse test!  old INTCFCOARSETEST==1
428         
429                 if(oldFlag & (CFMbndInflow)) {
430                         // fluid & if are ok, fill if later on
431                         int isValid = oldFlag & (CFFluid|CFInter);
432                         const LbmFloat iniRho = 1.0;
433                         const int OId = oldFlag>>24;
434                         if(!isValid) {
435                                 // make new if cell
436                                 const LbmVec vel(mObjectSpeeds[OId]);
437                                 // TODO add OPT3D treatment
438                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
439                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
440                                 RAC(tcel, dFlux) = FLUX_INIT;
441                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
442                                 calcCurrentMass += iniRho; calcCurrentVolume += 1.0; calcNumUsedCells++;
443                                 mInitialMass += iniRho;
444                                 //errMsg("INFLOW_DEBUG","ini at "<<PRINT_IJK<<" v="<<vel<<" inirho="<<iniRho);
445                                 // dont treat cell until next step
446                                 continue;
447                         } 
448                 } 
449                 else  // these are exclusive
450                 if(oldFlag & (CFMbndOutflow)) {
451                         int isnotValid = oldFlag & (CFFluid);
452                         if(isnotValid) {
453                                 // remove fluid cells, shouldnt be here anyway
454                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
455                                 mInitialMass -= fluidRho;
456                                 const LbmFloat iniRho = 0.0;
457                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
458                                 RAC(tcel, dFlux) = FLUX_INIT;
459                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
460
461                                 // same as ifemptied for if below
462                                 LbmPoint emptyp; emptyp.flag = 0;
463                                 emptyp.x = i; emptyp.y = j; emptyp.z = k;
464 #if PARALLEL==1
465                                 calcListEmpty[id].push_back( emptyp );
466 #else // PARALLEL==1
467                                 mListEmpty.push_back( emptyp );
468 #endif // PARALLEL==1
469                                 calcCellsEmptied++;
470                                 continue;
471                         }
472                 }
473
474                 if(oldFlag & (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)) { 
475                         *pFlagDst = oldFlag;
476                         continue;
477                 }
478                 /*if( oldFlag & CFNoBndFluid ) {  // TEST ME FASTER?
479                         OPTIMIZED_STREAMCOLLIDE; PERFORM_USQRMAXCHECK;
480                         RAC(tcel,dFfrac) = 1.0; 
481                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
482                         calcCurrentMass += rho; calcCurrentVolume += 1.0;
483                         calcNumUsedCells++;
484                         continue;
485                 }// TEST ME FASTER? */
486
487                 // only neighbor flags! not own flag
488                 nbored = 0;
489                 
490 #if OPT3D==0
491                 FORDF1 {
492                         nbflag[l] = RFLAG_NB(lev, i,j,k,SRCS(lev),l);
493                         nbored |= nbflag[l];
494                 } 
495 #else
496                 nbflag[dSB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dSB];
497                 nbflag[dWB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWB];
498                 nbflag[ dB] = *(pFlagSrc + (-mLevel[lev].lOffsy)); nbored |= nbflag[dB];
499                 nbflag[dEB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dEB];
500                 nbflag[dNB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNB];
501
502                 nbflag[dSW] = *(pFlagSrc + (-mLevel[lev].lOffsx+-1)); nbored |= nbflag[dSW];
503                 nbflag[ dS] = *(pFlagSrc + (-mLevel[lev].lOffsx)); nbored |= nbflag[dS];
504                 nbflag[dSE] = *(pFlagSrc + (-mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dSE];
505
506                 nbflag[ dW] = *(pFlagSrc + (-1)); nbored |= nbflag[dW];
507                 nbflag[ dE] = *(pFlagSrc + ( 1)); nbored |= nbflag[dE];
508
509                 nbflag[dNW] = *(pFlagSrc + ( mLevel[lev].lOffsx+-1)); nbored |= nbflag[dNW];
510           nbflag[ dN] = *(pFlagSrc + ( mLevel[lev].lOffsx)); nbored |= nbflag[dN];
511                 nbflag[dNE] = *(pFlagSrc + ( mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dNE];
512
513                 nbflag[dST] = *(pFlagSrc + ( mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dST];
514                 nbflag[dWT] = *(pFlagSrc + ( mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWT];
515                 nbflag[ dT] = *(pFlagSrc + ( mLevel[lev].lOffsy)); nbored |= nbflag[dT];
516                 nbflag[dET] = *(pFlagSrc + ( mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dET];
517                 nbflag[dNT] = *(pFlagSrc + ( mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNT];
518                 // */
519 #endif
520
521                 // pointer to destination cell
522                 calcNumUsedCells++;
523
524                 // FLUID cells 
525                 if( oldFlag & CFFluid ) { 
526                         // only standard fluid cells (with nothing except fluid as nbs
527
528                         if(oldFlag&CFMbndInflow) {
529                                 // force velocity for inflow, necessary to have constant direction of flow
530                                 // FIXME , test also set interface cells?
531                                 const int OId = oldFlag>>24;
532                                 //? DEFAULT_STREAM;
533                                 //const LbmFloat fluidRho = 1.0;
534                                 // for submerged inflows, streaming would have to be performed...
535                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
536                                 const LbmVec vel(mObjectSpeeds[OId]);
537                                 ux=vel[0], uy=vel[1], uz=vel[2]; 
538                                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
539                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); }
540                                 rho = fluidRho;
541                                 //errMsg("INFLOW_DEBUG","std at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
542                         } else {
543                                 if(nbored&CFBnd) {
544                                         DEFAULT_STREAM;
545                                         //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
546                                         DEFAULT_COLLIDEG(mLevel[lev].gravity);
547                                         oldFlag &= (~CFNoBndFluid);
548                                 } else {
549                                         // do standard stream/collide
550                                         OPTIMIZED_STREAMCOLLIDE;
551                                         oldFlag |= CFNoBndFluid;
552                                 } 
553                         }
554
555                         PERFORM_USQRMAXCHECK;
556                         // "normal" fluid cells
557                         RAC(tcel,dFfrac) = 1.0; 
558                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
559                         calcCurrentMass += rho; 
560                         calcCurrentVolume += 1.0;
561                         continue;
562                 }
563                 
564                 newFlag  = oldFlag;
565                 // make sure here: always check which flags to really unset...!
566                 newFlag = newFlag & (~( 
567                                         CFNoNbFluid|CFNoNbEmpty| CFNoDelete 
568                                         | CFNoInterpolSrc
569                                         | CFNoBndFluid
570                                         ));
571                 if(!(nbored&CFBndNoslip)) { //NEWSURFT NEWSURFTNOS
572                         newFlag |= CFNoBndFluid;
573                 }
574                 /*if(!(nbored&CFBnd)) { //NEWSURFT NEWSURFTNOS
575                         // explicitly check for noslip neighbors
576                         bool hasnoslipnb = false;
577                         FORDF1 { if((nbflag[l]&CFBnd)&&(nbflag[l]&CFBndNoslip)) hasnoslipnb=true; }
578                         if(!hasnoslipnb) newFlag |= CFNoBndFluid; 
579                 } // */
580
581                 // store own dfs and mass
582                 mass = RAC(ccel,dMass);
583
584                 // WARNING - only interface cells arrive here!
585                 // read distribution funtions of adjacent cells = stream step
586                 DEFAULT_STREAM;
587
588                 if((nbored & CFFluid)==0) { newFlag |= CFNoNbFluid; mNumInvIfCells++; }
589                 if((nbored & CFEmpty)==0) { newFlag |= CFNoNbEmpty; mNumInvIfCells++; }
590
591                 // calculate mass exchange for interface cells 
592                 LbmFloat myfrac = RAC(ccel,dFfrac);
593                 if(myfrac<0.) myfrac=0.; // NEWSURFT
594 #               define nbdf(l) m[ this->dfInv[(l)] ]
595
596                 // update mass 
597                 // only do boundaries for fluid cells, and interface cells without
598                 // any fluid neighbors (assume that interface cells _with_ fluid
599                 // neighbors are affected enough by these) 
600                 // which Df's have to be reconstructed? 
601                 // for fluid cells - just the f_i difference from streaming to empty cells  ----
602                 numRecons = 0;
603                 bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
604
605                 FORDF1 { // dfl loop
606                         recons[l] = 0;
607                         nbfracs[l] = 0.0;
608                         // finally, "normal" interface cells ----
609                         if( nbflag[l]&(CFFluid|CFBnd) ) { // NEWTEST! FIXME check!!!!!!!!!!!!!!!!!!
610                                 change = nbdf(l) - MYDF(l);
611                         }
612                         // interface cells - distuingish cells that shouldn't fill/empty 
613                         else if( nbflag[l] & CFInter ) {
614                                 
615                                 LbmFloat mynbfac,nbnbfac;
616                                 // NEW TEST t1
617                                 // t2 -> off
618                                 if((oldFlag&CFNoBndFluid)&&(nbflag[l]&CFNoBndFluid)) {
619                                         mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
620                                         nbnbfac = 1.0/mynbfac;
621                                         onlyBndnb = false;
622                                 } else {
623                                         mynbfac = nbnbfac = 1.0; // switch calc flux off
624                                         goto changeDefault;  // NEWSURFT
625                                         //change = 0.; goto changeDone;  // NEWSURFT
626                                 }
627                                 //mynbfac = nbnbfac = 1.0; // switch calc flux off t3
628
629                                 // perform interface case handling
630                                 if ((oldFlag|nbflag[l])&(CFNoNbFluid|CFNoNbEmpty)) {
631                                 switch (oldFlag&(CFNoNbFluid|CFNoNbEmpty)) {
632                                         case 0: 
633                                                 // we are a normal cell so... 
634                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
635                                                         case CFNoNbFluid: 
636                                                                 // just fill current cell = empty neighbor 
637                                                                 change = nbnbfac*nbdf(l) ; goto changeDone; 
638                                                         case CFNoNbEmpty: 
639                                                                 // just empty current cell = fill neighbor 
640                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
641                                                 }
642                                                 break;
643
644                                         case CFNoNbFluid: 
645                                                 // we dont have fluid nb's so...
646                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
647                                                         case 0: 
648                                                         case CFNoNbEmpty: 
649                                                                 // we have no fluid nb's -> just empty
650                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
651                                                 }
652                                                 break;
653
654                                         case CFNoNbEmpty: 
655                                                 // we dont have empty nb's so...
656                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
657                                                         case 0: 
658                                                         case CFNoNbFluid: 
659                                                                 // we have no empty nb's -> just fill
660                                                                 change = nbnbfac*nbdf(l); goto changeDone; 
661                                                 }
662                                                 break;
663                                 }} // inter-inter exchange
664
665                         changeDefault: ;
666                                 // just do normal mass exchange...
667                                 change = ( nbnbfac*nbdf(l) - mynbfac*MYDF(l) ) ;
668                         changeDone: ;
669                                 nbfracs[l] = QCELL_NB(lev, i,j,k, SRCS(lev),l, dFfrac);
670                                 if(nbfracs[l]<0.) nbfracs[l] = 0.; // NEWSURFT
671                                 change *=  (myfrac + nbfracs[l]) * 0.5;
672                         } // the other cell is interface
673
674                         // last alternative - reconstruction in this direction
675                         else {
676                                 // empty + bnd case
677                                 recons[l] = 1; 
678                                 numRecons++;
679                                 change = 0.0; 
680                         }
681
682                         // modify mass at SRCS
683                         mass += change;
684                 } // l
685                 // normal interface, no if empty/fluid
686
687                 LbmFloat nv1,nv2;
688                 LbmFloat nx,ny,nz;
689
690                 if(nbflag[dE] &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } else nv1 = 0.0;
691                 if(nbflag[dW] &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } else nv2 = 0.0;
692                 nx = 0.5* (nv2-nv1);
693                 if(nbflag[dN] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
694                 if(nbflag[dS] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
695                 ny = 0.5* (nv2-nv1);
696 #if LBMDIM==3
697                 if(nbflag[dT] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
698                 if(nbflag[dB] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
699                 nz = 0.5* (nv2-nv1);
700 #else // LBMDIM==3
701                 nz = 0.0;
702 #endif // LBMDIM==3
703
704                 if( (ABS(nx)+ABS(ny)+ABS(nz)) > LBM_EPSILON) {
705                         // normal ok and usable...
706                         FORDF1 {
707                                 if( (this->dfDvecX[l]*nx + this->dfDvecY[l]*ny + this->dfDvecZ[l]*nz)  // dot Dvec,norml
708                                                 > LBM_EPSILON) {
709                                         recons[l] = 2; 
710                                         numRecons++;
711                                 } 
712                         }
713                 }
714
715                 // calculate macroscopic cell values
716                 LbmFloat oldUx, oldUy, oldUz;
717                 LbmFloat oldRho; // OLD rho = ccel->rho;
718 #if OPT3D==0
719                         oldRho=RAC(ccel,0);
720                         oldUx = oldUy = oldUz = 0.0;
721                         for(int l=1; l<this->cDfNum; l++) {
722                                 oldRho += RAC(ccel,l);
723                                 oldUx  += (this->dfDvecX[l]*RAC(ccel,l));
724                                 oldUy  += (this->dfDvecY[l]*RAC(ccel,l)); 
725                                 oldUz  += (this->dfDvecZ[l]*RAC(ccel,l)); 
726                         } 
727 #else // OPT3D==0
728                 oldRho = + RAC(ccel,dC)  + RAC(ccel,dN )
729                                 + RAC(ccel,dS ) + RAC(ccel,dE )
730                                 + RAC(ccel,dW ) + RAC(ccel,dT )
731                                 + RAC(ccel,dB ) + RAC(ccel,dNE)
732                                 + RAC(ccel,dNW) + RAC(ccel,dSE)
733                                 + RAC(ccel,dSW) + RAC(ccel,dNT)
734                                 + RAC(ccel,dNB) + RAC(ccel,dST)
735                                 + RAC(ccel,dSB) + RAC(ccel,dET)
736                                 + RAC(ccel,dEB) + RAC(ccel,dWT)
737                                 + RAC(ccel,dWB);
738
739                         oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
740                                 + RAC(ccel,dNE) - RAC(ccel,dNW)
741                                 + RAC(ccel,dSE) - RAC(ccel,dSW)
742                                 + RAC(ccel,dET) + RAC(ccel,dEB)
743                                 - RAC(ccel,dWT) - RAC(ccel,dWB);
744
745                         oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
746                                 + RAC(ccel,dNE) + RAC(ccel,dNW)
747                                 - RAC(ccel,dSE) - RAC(ccel,dSW)
748                                 + RAC(ccel,dNT) + RAC(ccel,dNB)
749                                 - RAC(ccel,dST) - RAC(ccel,dSB);
750
751                         oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
752                                 + RAC(ccel,dNT) - RAC(ccel,dNB)
753                                 + RAC(ccel,dST) - RAC(ccel,dSB)
754                                 + RAC(ccel,dET) - RAC(ccel,dEB)
755                                 + RAC(ccel,dWT) - RAC(ccel,dWB);
756 #endif
757
758                 // now reconstruction
759 #define REFERENCE_PRESSURE 1.0 // always atmosphere...
760 #if OPT3D==0
761                 // NOW - construct dist funcs from empty cells
762                 FORDF1 {
763                         if(recons[ l ]) {
764                                 m[ this->dfInv[l] ] = 
765                                         this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 
766                                         this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 
767                                         - MYDF( l );
768                         }
769                 }
770                 usqr = 1.5 * (oldUx*oldUx + oldUy*oldUy + oldUz*oldUz); // needed later on
771 #else
772                 ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
773                 rho = REFERENCE_PRESSURE;
774                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
775                 if(recons[dN ]) { m[dS ] = EQN  + EQS  - MYDF(dN ); }
776                 if(recons[dS ]) { m[dN ] = EQS  + EQN  - MYDF(dS ); }
777                 if(recons[dE ]) { m[dW ] = EQE  + EQW  - MYDF(dE ); }
778                 if(recons[dW ]) { m[dE ] = EQW  + EQE  - MYDF(dW ); }
779                 if(recons[dT ]) { m[dB ] = EQT  + EQB  - MYDF(dT ); }
780                 if(recons[dB ]) { m[dT ] = EQB  + EQT  - MYDF(dB ); }
781                 if(recons[dNE]) { m[dSW] = EQNE + EQSW - MYDF(dNE); }
782                 if(recons[dNW]) { m[dSE] = EQNW + EQSE - MYDF(dNW); }
783                 if(recons[dSE]) { m[dNW] = EQSE + EQNW - MYDF(dSE); }
784                 if(recons[dSW]) { m[dNE] = EQSW + EQNE - MYDF(dSW); }
785                 if(recons[dNT]) { m[dSB] = EQNT + EQSB - MYDF(dNT); }
786                 if(recons[dNB]) { m[dST] = EQNB + EQST - MYDF(dNB); }
787                 if(recons[dST]) { m[dNB] = EQST + EQNB - MYDF(dST); }
788                 if(recons[dSB]) { m[dNT] = EQSB + EQNT - MYDF(dSB); }
789                 if(recons[dET]) { m[dWB] = EQET + EQWB - MYDF(dET); }
790                 if(recons[dEB]) { m[dWT] = EQEB + EQWT - MYDF(dEB); }
791                 if(recons[dWT]) { m[dEB] = EQWT + EQEB - MYDF(dWT); }
792                 if(recons[dWB]) { m[dET] = EQWB + EQET - MYDF(dWB); }
793 #endif          
794
795
796                 // inflow bc handling
797                 if(oldFlag & (CFMbndInflow)) {
798                         // fill if cells in inflow region
799                         if(myfrac<0.5) { 
800                                 mass += 0.25; 
801                                 mInitialMass += 0.25;
802                         }
803                         const int OId = oldFlag>>24;
804                         const LbmVec vel(mObjectSpeeds[OId]);
805                         ux=vel[0], uy=vel[1], uz=vel[2]; 
806                         //? usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
807                         //FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); } rho = fluidRho;
808                         rho = REFERENCE_PRESSURE;
809                         FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
810                         //errMsg("INFLOW_DEBUG","if at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
811                 }  else { 
812                         // NEWSURFT, todo optimize!
813                         if(onlyBndnb) { //if(usqr<0.001*0.001) {
814                                 rho = ux = uy = uz = 0.;
815                                 FORDF0{ 
816                                         rho += m[l]; 
817                                         ux  += (this->dfDvecX[l]*m[l]); 
818                                         uy  += (this->dfDvecY[l]*m[l]);  
819                                         uz  += (this->dfDvecZ[l]*m[l]);  
820                                 }
821                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
822                         } else {// NEWSURFT */
823                                 // normal collide
824                                 // mass streaming done... do normal collide
825                                 LbmVec grav = mLevel[lev].gravity*mass;
826                                 DEFAULT_COLLIDEG(grav);
827                                 PERFORM_USQRMAXCHECK;
828                                 // rho init from default collide necessary for fill/empty check below
829                         } // test
830                 }
831
832                 // testing..., particle generation
833                 // also check oldFlag for CFNoNbFluid, test
834                 // for inflow no pargen test // NOBUBBB!
835                 if((this->mInitDone) //&&(mUseTestdata) 
836                                 //&& (!((oldFlag|newFlag)&CFNoNbEmpty)) 
837                                 && (!((oldFlag|newFlag)&CFNoDelete)) 
838                                 && (this->mPartGenProb>0.0)) {
839                         bool doAdd = true;
840                         bool bndOk=true;
841                         //if( (i<1+mCutoff)||(i>this->mSizex-1-mCutoff)||
842                                         //(j<1+mCutoff)||(j>this->mSizey-1-mCutoff)||
843                                         //(k<1+mCutoff)||(k>this->mSizez-1-mCutoff) ) { bndOk=false; }
844                         if( (i<cutMin)||(i>this->mSizex-cutMin)||
845                                         (j<cutMin)||(j>this->mSizey-cutMin)||
846                                         (k<cutMin)||(k>this->mSizez-cutMin) ) { bndOk=false; }
847                         if(!bndOk) doAdd=false;
848                         
849                         LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep);
850                         LbmFloat rux = (ux * realWorldFac);
851                         LbmFloat ruy = (uy * realWorldFac);
852                         LbmFloat ruz = (uz * realWorldFac);
853                         LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz));
854                         // WHMOD
855
856                         LbmFloat prob = (rand()/(RAND_MAX+1.0));
857                         LbmFloat basethresh = this->mPartGenProb*lcsmqo*rl;
858
859                         // reduce probability in outer region
860                         //if( (i<cutConst)||(i>this->mSizex-cutConst) ){ prob *= 0.5; }
861                         //if( (j<cutConst)||(j>this->mSizey-cutConst) ){ prob *= 0.5; }
862                         //if( (k<cutConst)||(k>this->mSizez-cutConst) ){ prob *= 0.5; }
863                         // NEW TEST 040627
864                         //if( (i<cutConst)||(i>this->mSizex-cutConst) ){ doAdd=false; }
865                         //if( (j<cutConst)||(j>this->mSizey-cutConst) ){ doAdd=false; }
866                         //if( (k<cutConst)||(k>this->mSizez-cutConst) ){ doAdd=false; }
867                         const int pibord = mLevel[mMaxRefine].lSizex/2-cutConst;
868                         const int pjbord = mLevel[mMaxRefine].lSizey/2-cutConst;
869                         LbmFloat pifac = 1.-(LbmFloat)(ABS(i-pibord)) / (LbmFloat)(pibord);
870                         LbmFloat pjfac = 1.-(LbmFloat)(ABS(j-pjbord)) / (LbmFloat)(pjbord);
871                         if(pifac<0.) pifac=0.;
872                         if(pjfac<0.) pjfac=0.;
873                         //const LbmFloat pkfac = 1.-(LbmFloat)(ABS(k-mLevel[mMaxRefine].lSizez/2))/(LbmFloat)(mLevel[mMaxRefine].lSizez/2);
874                         //errMsg("PROBTTT"," at "<<PRINT_IJK<<" prob"<<prob<<" pifac"<<pifac<<" pjfac"<<pjfac<<" "<<(basethresh*rl*pifac*pjfac));
875                         //prob *= pifac*pjfac; //*pkfac;
876
877 //#define RWVEL_THRESH 1.0
878 #define RWVEL_THRESH 1.5
879 #define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5)
880                         //if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
881                         if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
882                                 // add
883                         } else {
884                                 doAdd = false; // dont...
885                         }
886
887 //#define SLOWDOWNREGION (2*mCutoff)
888 #if LBMDIM==3
889                         // normal
890 #define SLOWDOWNREGION (this->mSizez/4)
891 #else // LBMDIM==2
892                         // off
893 #define SLOWDOWNREGION 10 
894 #endif // LBMDIM==2
895 #define P_LCSMQO 0.01
896
897                         // "wind" disturbance
898                         // use realworld relative velocity here instead?
899                         if( (doAdd && 
900                                         ((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) )// normal checks
901                                         ||(k>this->mSizez-SLOWDOWNREGION)   ) {
902                                 LbmFloat nuz = uz;
903                                 //? LbmFloat jdf; // = 0.05 * (rand()/(RAND_MAX+1.0));
904                                 //? if(rl>RWVEL_WINDTHRESH) jdf *= (rl-RWVEL_WINDTHRESH);
905                                 if(k>this->mSizez-SLOWDOWNREGION) {
906                                         // special case
907                                         LbmFloat zfac = (LbmFloat)( k-(this->mSizez-SLOWDOWNREGION) );
908                                         zfac /= (LbmFloat)(SLOWDOWNREGION);
909                                         nuz += (1.0) * zfac; // check max speed? OFF?
910                                         //errMsg("TOPT"," at "<<PRINT_IJK<<" zfac"<<zfac);
911                                 } else {
912                                         // normal probability
913                                         //? LbmFloat fac = P_LCSMQO-lcsmqo;
914                                         //? jdf *= fac;
915                                 }
916                                 FORDF1 {
917                                         const LbmFloat jdf = 0.05 * (rand()/(RAND_MAX+1.0));
918                                         // TODO  use wind velocity?
919                                         if(jdf>0.025) {
920                                         const LbmFloat add =  this->dfLength[l]*(-ux*this->dfDvecX[l]-uy*this->dfDvecY[l]-nuz*this->dfDvecZ[l])*jdf;
921                                         RAC(tcel,l) += add; }
922                                 }
923                                 //errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<PRINT_VEC(ux,uy,nuz)<<" rwv"<<PRINT_VEC(rux,ruy,ruz) );
924                                 //errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) );
925                         } // wind disturbance
926
927                         if(usqr<0.0001) doAdd=false;   // TODO check!?
928                         // if outside, and 20% above sea level, delete, TODO really check level?
929                         //if((!bndOk)&&((LbmFloat)k>pTest->mFluidHeight*1.5)) { doAdd=true; } // FORCEDISSOLVE
930                         //if(this->mStepCnt>700) errMsg("DFJITT"," at "<<PRINT_IJK<<"rwl:"<<rl<<"  usqr:"<<usqr <<" qo:"<<lcsmqo<<" add="<<doAdd );
931
932                         if( (doAdd)  ) { // ADD DROP
933                                 LbmFloat len = norm(LbmVec(ux,uy,uz));
934                                 // WHMOD
935                                 //for(int s=0; s<10; s++) { // multiple parts
936                                 for(int s=0; s<1; s++) { // one part!
937                                 //LbmFloat prob = this->mPartGenProb * 0.02* (rand()/(RAND_MAX+1.0));
938                                 const LbmFloat posjitter = 1.0;
939                                 const LbmFloat posjitteroffs = posjitter*-0.5;
940                                 LbmFloat jpx = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
941                                 LbmFloat jpy = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
942                                 LbmFloat jpz = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
943
944                                 const LbmFloat jitterstr = 0.1;
945                                 const LbmFloat jitteroffs = jitterstr*-0.5;
946                                 LbmFloat jx = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
947                                 LbmFloat jy = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
948                                 LbmFloat jz = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
949
950                                 // average normal & velocity 
951                                 // -> mostly along velocity dir, many into surface
952                                 LbmVec pv = (LbmVec(nx+jx,ny+jy,nz+jz)*0.75 + getNormalized(LbmVec(ux,uy,uz)) )*0.35; 
953                                 normalize(pv);
954
955                                 LbmFloat srci = i+0.5+jpx; // TEST? + (pv[0]*1.41);
956                                 LbmFloat srcj = j+0.5+jpy; // TEST? + (pv[1]*1.41);
957                                 LbmFloat srck = k+0.5+jpz; // TEST? + (pv[2]*1.41);
958                                 int type=0;
959                                 //if((s%3)!=2) {} else { type=PART_FLOAT; }
960                                 //type = PART_DROP;
961                                 type = PART_INTER;
962                                 // drop
963                                 /*srci += (pv[0]*1.41);
964                                 srcj += (pv[1]*1.41);
965                                 srck += (pv[2]*1.41);
966                                 if(!(RFLAG(lev, (int)(srci),(int)(srcj),(int)(srck),SRCS(lev)) &CFEmpty)) continue; // only add in good direction */
967
968                                 pv *= len;
969                                 LbmFloat size = 1.0+ 9.0* (rand()/(RAND_MAX+1.0));
970
971                                 mpParticles->addParticle(srci, srcj, srck); //i+0.5+jpx,j+0.5+jpy,k+0.5+jpz);
972                                 mpParticles->getLast()->setVel(pv[0],pv[1],pv[2]);
973                                 //? mpParticles->getLast()->advanceVel(); // advance a bit outwards
974                                 mpParticles->getLast()->setStatus(PART_IN);
975                                 mpParticles->getLast()->setType(type);
976                                 //if((s%3)==2) mpParticles->getLast()->setType(PART_FLOAT);
977                                 mpParticles->getLast()->setSize(size);
978                                 //errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<PRINT_VEC(ux,uy,uz) <<" RWu="<<PRINT_VEC(rux,ruy,ruz)<<" add"<<doAdd<<" pvel="<<pv );
979                                 //errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<norm(LbmVec(ux,uy,uz)) <<" RWu="<<norm(LbmVec(rux,ruy,ruz))<<" add"<<doAdd<<" pvel="<<norm(pv) );
980                                 //if(!bndOk){ mass -= val2fac*size*0.02; } // FORCEDISSOLVE
981                                 //const LbmFloat val2fac = 1.0; //? TODO N test? /(this->mPartGenProb); // FIXME remove factor!
982                                 //mass -= val2fac*size*0.0015; // NTEST!
983                                 //mass -= val2fac*size*0.001; // NTEST!
984                                 mass -= size*0.0020; // NTEST!
985 #if LBMDIM==2
986                                 mpParticles->getLast()->setVel(pv[0],pv[1],0.0);
987                                 mpParticles->getLast()->setPos(ntlVec3Gfx(srci,srcj,0.5));
988 #endif // LBMDIM==2
989                 //errMsg("PIT","NEW pit"<<mpParticles->getLast()->getId()<<" pos:"<< mpParticles->getLast()->getPos()<<" status:"<<convertFlags2String(mpParticles->getLast()->getFlags())<<" vel:"<< mpParticles->getLast()->getVel()  );
990                                 } // multiple parts
991                         } // doAdd
992                 } // */
993
994
995                 // interface cell filled or emptied?
996                 iffilled = ifemptied = 0;
997                 // interface cells empty/full?, WARNING: to mark these cells, better do it at the end of reinitCellFlags
998                 // interface cell if full?
999                 if( (mass) >= (rho * (1.0+FSGR_MAGICNR)) ) { iffilled = 1; }
1000                 // interface cell if empty?
1001                 if( (mass) <= (rho * (   -FSGR_MAGICNR)) ) { ifemptied = 1; }
1002
1003                 if(oldFlag & (CFMbndOutflow)) {
1004                         mInitialMass -= mass;
1005                         mass = myfrac = 0.0;
1006                         iffilled = 0; ifemptied = 1;
1007                 }
1008
1009                 // looks much nicer... LISTTRICK
1010 #if FSGR_LISTTRICK==1
1011                 if(newFlag&CFNoBndFluid) { // test NEW TEST
1012                         if(!iffilled) {
1013                                 // remove cells independent from amount of change...
1014                                 if( (oldFlag & CFNoNbEmpty)&&(newFlag & CFNoNbEmpty)&&
1015                                                 ( (mass>(rho*FSGR_LISTTTHRESHFULL))  || ((nbored&CFInter)==0)  )) { 
1016                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","filled "<<PRINT_IJK); };
1017                                         iffilled = 1; 
1018                                 } 
1019                         }
1020                         if(!ifemptied) {
1021                                 if( (oldFlag & CFNoNbFluid)&&(newFlag & CFNoNbFluid)&&
1022                                                 ( (mass<(rho*FSGR_LISTTTHRESHEMPTY)) || ((nbored&CFInter)==0)  )) { 
1023                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","emptied "<<PRINT_IJK); };
1024                                         ifemptied = 1; 
1025                                 } 
1026                         }
1027                 } // nobndfluid only */
1028 #endif
1029                 //iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!!
1030                 
1031
1032                 // now that all dfs are known, handle last changes
1033                 if(iffilled) {
1034                         LbmPoint filledp; filledp.flag=0;
1035                         if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1;  // NEWSURFT
1036                         filledp.x = i; filledp.y = j; filledp.z = k;
1037 #if PARALLEL==1
1038                         calcListFull[id].push_back( filledp );
1039 #else // PARALLEL==1
1040                         mListFull.push_back( filledp );
1041 #endif // PARALLEL==1
1042                         //this->mNumFilledCells++; // DEBUG
1043                         calcCellsFilled++;
1044                 }
1045                 else if(ifemptied) {
1046                         LbmPoint emptyp; emptyp.flag=0;
1047                         if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; //  NEWSURFT
1048                         emptyp.x = i; emptyp.y = j; emptyp.z = k;
1049 #if PARALLEL==1
1050                         calcListEmpty[id].push_back( emptyp );
1051 #else // PARALLEL==1
1052                         mListEmpty.push_back( emptyp );
1053 #endif // PARALLEL==1
1054                         //this->mNumEmptiedCells++; // DEBUG
1055                         calcCellsEmptied++;
1056                 } 
1057                 // dont cutoff values -> better cell conversions
1058                 RAC(tcel,dFfrac)   = (mass/rho);
1059
1060                 // init new flux value
1061                 float flux = FLUX_INIT; // dxqn on
1062                 if(newFlag&CFNoBndFluid) {
1063                         //flux = 50.0; // extreme on
1064                         for(int nn=1; nn<this->cDfNum; nn++) { 
1065                                 if(nbflag[nn] & (CFFluid|CFInter|CFBnd)) { flux += this->dfLength[nn]; }
1066                         }
1067                         // optical hack - smooth slow moving
1068                         // surface regions
1069                         if(usqr< sssUsqrLimit) {
1070                         for(int nn=1; nn<this->cDfNum; nn++) { 
1071                                 if(nbfracs[nn]!=0.0) {
1072                                         LbmFloat curSmooth = (sssUsqrLimit-usqr)*sssUsqrLimitInv;
1073                                         if(curSmooth>1.0) curSmooth=1.0;
1074                                         flux *= (1.0+ smoothStrength*curSmooth * (nbfracs[nn]-myfrac)) ;
1075                                 }
1076                         } }
1077                         // NEW TEST */
1078                 }
1079                 // flux = FLUX_INIT; // calc flux off
1080                 QCELL(lev, i,j,k,TSET(lev), dFlux) = flux; // */
1081
1082                 // perform mass exchange with streamed values
1083                 QCELL(lev, i,j,k,TSET(lev), dMass) = mass; // MASST
1084                 // set new flag 
1085                 *pFlagDst = (CellFlagType)newFlag;
1086                 calcCurrentMass += mass; 
1087                 calcCurrentVolume += RAC(tcel,dFfrac);
1088
1089                 // interface cell handling done...
1090         } // i
1091         int i=0; //dummy
1092         ADVANCE_POINTERS(2);
1093         } // j
1094
1095 #if COMPRESSGRIDS==1
1096 #if PARALLEL==1
1097         //frintf(stderr," (id=%d k=%d) ",id,k);
1098 # pragma omp barrier
1099 #endif // PARALLEL==1
1100 #else // COMPRESSGRIDS==1
1101         int i=0; //dummy
1102         ADVANCE_POINTERS(mLevel[lev].lSizex*2);
1103 #endif // COMPRESSGRIDS==1
1104   } // all cell loop k,j,i
1105
1106         } // main loop region
1107
1108         // write vars from parallel computations to class
1109         //errMsg("DFINI"," maxr l"<<mMaxRefine<<" cm="<<calcCurrentMass<<" cv="<<calcCurrentVolume );
1110         mLevel[lev].lmass    = calcCurrentMass;
1111         mLevel[lev].lvolume  = calcCurrentVolume;
1112         this->mNumFilledCells  = calcCellsFilled;
1113         this->mNumEmptiedCells = calcCellsEmptied;
1114         this->mNumUsedCells = calcNumUsedCells;
1115 #if PARALLEL==1
1116         PARA_FINISH();
1117 #endif // PARALLEL==1
1118
1119         // check other vars...?
1120 }
1121
1122
1123
1124 /******************************************************************************
1125  * work on lists from updateCellMass to reinit cell flags
1126  *****************************************************************************/
1127
1128 LbmFloat LbmFsgrSolver::getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l) {
1129         //return 0.0; // test
1130         int level = mMaxRefine;
1131         LbmFloat *ccel = RACPNT(level, i,j,k, workSet);
1132
1133         LbmFloat nx,ny,nz, nv1,nv2;
1134         if(RFLAG_NB(level,i,j,k,workSet, dE) &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } else nv1 = 0.0;
1135         if(RFLAG_NB(level,i,j,k,workSet, dW) &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } else nv2 = 0.0;
1136         nx = 0.5* (nv2-nv1);
1137         if(RFLAG_NB(level,i,j,k,workSet, dN) &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
1138         if(RFLAG_NB(level,i,j,k,workSet, dS) &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
1139         ny = 0.5* (nv2-nv1);
1140 #if LBMDIM==3
1141         if(RFLAG_NB(level,i,j,k,workSet, dT) &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
1142         if(RFLAG_NB(level,i,j,k,workSet, dB) &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
1143         nz = 0.5* (nv2-nv1);
1144 #else //LBMDIM==3
1145         nz = 0.0;
1146 #endif //LBMDIM==3
1147         LbmFloat scal = mDvecNrm[l][0]*nx + mDvecNrm[l][1]*ny + mDvecNrm[l][2]*nz;
1148
1149         LbmFloat ret = 1.0;
1150         // forward direction, add mass (for filling cells):
1151         if(dirForw) {
1152                 if(scal<LBM_EPSILON) ret = 0.0;
1153                 else ret = scal;
1154         } else {
1155                 // backward for emptying
1156                 if(scal>-LBM_EPSILON) ret = 0.0;
1157                 else ret = scal * -1.0;
1158         }
1159         //errMsg("massd", PRINT_IJK<<" nv"<<nvel<<" : ret="<<ret ); //xit(1); //VECDEB
1160         return ret;
1161 }
1162
1163 void LbmFsgrSolver::addToNewInterList( int ni, int nj, int nk ) {
1164 #if FSGR_STRICT_DEBUG==10
1165         // dangerous, this can change the simulation...
1166   /*for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1167        iter != mListNewInter.end(); iter++ ) {
1168     if(ni!=iter->x) continue;
1169     if(nj!=iter->y) continue;
1170     if(nk!=iter->z) continue;
1171                 // all 3 values match... skip point
1172                 return;
1173         } */
1174 #endif // FSGR_STRICT_DEBUG==1
1175         // store point
1176         LbmPoint newinter; newinter.flag = 0;
1177         newinter.x = ni; newinter.y = nj; newinter.z = nk;
1178         mListNewInter.push_back(newinter);
1179 }
1180
1181 void LbmFsgrSolver::reinitFlags( int workSet ) { 
1182         // reinitCellFlags OLD mods:
1183         // add all to intel list?
1184         // check ffrac for new cells
1185         // new if cell inits (last loop)
1186         // vweights handling
1187
1188         const int debugFlagreinit = 0;
1189         
1190         // some things need to be read/modified on the other set
1191         int otherSet = (workSet^1);
1192         // fixed level on which to perform 
1193         int workLev = mMaxRefine;
1194
1195   /* modify interface cells from lists */
1196   /* mark filled interface cells as fluid, or emptied as empty */
1197         /* count neighbors and distribute excess mass to interface neighbor cells
1198    * problems arise when there are no interface neighbors anymore
1199          * then just distribute to any fluid neighbors...
1200          */
1201
1202         // for symmetry, first init all neighbor cells */
1203         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1204        iter != mListFull.end(); iter++ ) {
1205     int i=iter->x, j=iter->y, k=iter->z;
1206                 if(debugFlagreinit) errMsg("FULL", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
1207     FORDF1 {
1208                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1209                         //if((LBMDIM>2)&&( (ni<=0) || (nj<=0) || (nk<=0) || (ni>=mLevel[workLev].lSizex-1) || (nj>=mLevel[workLev].lSizey-1) || (nk>=mLevel[workLev].lSizez-1) )) {
1210                         if( (ni<=0) || (nj<=0) || 
1211                                   (ni>=mLevel[workLev].lSizex-1) ||
1212                                   (nj>=mLevel[workLev].lSizey-1) 
1213 #                                       if LBMDIM==3
1214                                   || (nk<=0) || (nk>=mLevel[workLev].lSizez-1) 
1215 #                                       endif // LBMDIM==1
1216                                  ) {
1217                                 continue; } // new bc, dont treat cells on boundary NEWBC
1218       if( RFLAG(workLev, ni,nj,nk, workSet) & CFEmpty ){
1219                                 
1220                                 // preinit speed, get from average surrounding cells
1221                                 // interpolate from non-workset to workset, sets are handled in function
1222
1223                                 // new and empty interface cell, dont change old flag here!
1224                                 addToNewInterList(ni,nj,nk);
1225
1226                                 LbmFloat avgrho = 0.0;
1227                                 LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
1228                                 interpolateCellValues(workLev,ni,nj,nk,workSet, avgrho,avgux,avguy,avguz);
1229
1230                                 // careful with l's...
1231                                 FORDF0M { 
1232                                         QCELL(workLev,ni,nj,nk, workSet, m) = this->getCollideEq( m,avgrho,  avgux, avguy, avguz ); 
1233                                         //QCELL(workLev,ni,nj,nk, workSet, l) = avgnbdf[l]; // CHECK FIXME test?
1234                                 }
1235                                 //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
1236                                 QCELL(workLev,ni,nj,nk, workSet, dMass) = 0.0; //?? new
1237                                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) = 0.0; //?? new
1238                                 //RFLAG(workLev,ni,nj,nk,workSet) = (CellFlagType)(CFInter|CFNoInterpolSrc);
1239                                 changeFlag(workLev,ni,nj,nk,workSet, (CFInter|CFNoInterpolSrc));
1240                                 if(debugFlagreinit) errMsg("NEWE", PRINT_IJK<<" newif "<<PRINT_VEC(ni,nj,nk)<<" rho"<<avgrho<<" vel("<<avgux<<","<<avguy<<","<<avguz<<") " );
1241       }
1242                         /* prevent surrounding interface cells from getting removed as empty cells 
1243                          * (also cells that are not newly inited) */
1244       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1245                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete);
1246                                 changeFlag(workLev,ni,nj,nk, workSet, (RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete));
1247                                 // also add to list...
1248                                 addToNewInterList(ni,nj,nk);
1249                         } // NEW?
1250     }
1251
1252                 // NEW? no extra loop...
1253                 //RFLAG(workLev,i,j,k, workSet) = CFFluid;
1254                 changeFlag(workLev,i,j,k, workSet,CFFluid);
1255         }
1256
1257         /* remove empty interface cells that are not allowed to be removed anyway
1258          * this is important, otherwise the dreaded cell-type-flickering can occur! */
1259   for(int n=0; n<(int)mListEmpty.size(); n++) {
1260     int i=mListEmpty[n].x, j=mListEmpty[n].y, k=mListEmpty[n].z;
1261                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) {
1262                         // treat as "new inter"
1263                         addToNewInterList(i,j,k);
1264                         // remove entry
1265                         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
1266                         if(n<(int)mListEmpty.size()-1) mListEmpty[n] = mListEmpty[mListEmpty.size()-1]; 
1267                         mListEmpty.pop_back();
1268                         n--; 
1269                 }
1270         } 
1271
1272
1273         /* problems arise when adjacent cells empty&fill ->
1274                  let fill cells+surrounding interface cells have the higher importance */
1275   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1276        iter != mListEmpty.end(); iter++ ) {
1277     int i=iter->x, j=iter->y, k=iter->z;
1278                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)){ errMsg("A"," ARGHARGRAG "); } // DEBUG
1279                 if(debugFlagreinit) errMsg("EMPT", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) );
1280
1281                 /* set surrounding fluid cells to interface cells */
1282     FORDF1 {
1283                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1284       if( RFLAG(workLev,ni,nj,nk, workSet) & CFFluid){
1285                                 // init fluid->interface 
1286                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(CFInter); 
1287                                 changeFlag(workLev,ni,nj,nk, workSet, CFInter); 
1288                                 /* new mass = current density */
1289                                 LbmFloat nbrho = QCELL(workLev,ni,nj,nk, workSet, dC);
1290                 for(int rl=1; rl< this->cDfNum ; ++rl) { nbrho += QCELL(workLev,ni,nj,nk, workSet, rl); }
1291                                 QCELL(workLev,ni,nj,nk, workSet, dMass) =  nbrho; 
1292                                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) =  1.0; 
1293
1294                                 // store point
1295                                 addToNewInterList(ni,nj,nk);
1296       }
1297       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter){
1298                                 // test, also add to list...
1299                                 addToNewInterList(ni,nj,nk);
1300                         } // NEW?
1301     }
1302
1303                 /* for symmetry, set our flag right now */
1304                 changeFlag(workLev,i,j,k, workSet, CFEmpty);
1305                 // mark cell not be changed mass... - not necessary, not in list anymore anyway!
1306         } // emptylist
1307
1308
1309         
1310         // precompute weights to get rid of order dependancies
1311         vector<lbmFloatSet> vWeights;
1312         vWeights.resize( mListFull.size() + mListEmpty.size() );
1313         int weightIndex = 0;
1314   int nbCount = 0;
1315         LbmFloat nbWeights[LBM_DFNUM];
1316         LbmFloat nbTotWeights = 0.0;
1317         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1318        iter != mListFull.end(); iter++ ) {
1319     int i=iter->x, j=iter->y, k=iter->z;
1320     nbCount = 0; nbTotWeights = 0.0;
1321     FORDF1 {
1322                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1323       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1324                                 nbCount++;
1325                                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
1326                                 else nbWeights[l] = getMassdWeight(1,i,j,k,workSet,l); // NEWSURFT
1327                                 nbTotWeights += nbWeights[l];
1328       } else {
1329                                 nbWeights[l] = -100.0; // DEBUG;
1330                         }
1331     }
1332                 if(nbCount>0) { 
1333                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
1334         vWeights[weightIndex].val[0] = nbTotWeights;
1335         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
1336         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
1337                 } else { 
1338         vWeights[weightIndex].numNbs = 0.0;
1339                 }
1340                 weightIndex++;
1341         }
1342   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1343        iter != mListEmpty.end(); iter++ ) {
1344     int i=iter->x, j=iter->y, k=iter->z;
1345     nbCount = 0; nbTotWeights = 0.0;
1346     FORDF1 {
1347                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1348       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1349                                 nbCount++;
1350                                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
1351                                 else nbWeights[l] = getMassdWeight(0,i,j,k,workSet,l); // NEWSURFT
1352                                 nbTotWeights += nbWeights[l];
1353       } else {
1354                                 nbWeights[l] = -100.0; // DEBUG;
1355                         }
1356     }
1357                 if(nbCount>0) { 
1358                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
1359         vWeights[weightIndex].val[0] = nbTotWeights;
1360         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
1361         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
1362                 } else { 
1363         vWeights[weightIndex].numNbs = 0.0;
1364                 }
1365                 weightIndex++;
1366         } 
1367         weightIndex = 0;
1368         
1369
1370         /* process full list entries, filled cells are done after this loop */
1371         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1372        iter != mListFull.end(); iter++ ) {
1373     int i=iter->x, j=iter->y, k=iter->z;
1374
1375                 LbmFloat myrho = QCELL(workLev,i,j,k, workSet, dC);
1376     FORDF1 { myrho += QCELL(workLev,i,j,k, workSet, l); } // QCELL.rho
1377
1378     LbmFloat massChange = QCELL(workLev,i,j,k, workSet, dMass) - myrho;
1379                 if(vWeights[weightIndex].numNbs>0.0) {
1380                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
1381                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
1382                         FORDF1 {
1383                                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1384         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1385                                         LbmFloat change = -1.0;
1386                                         if(nbTotWeightsp>0.0) {
1387                                                 //change = massChange * ( nbWeights[l]/nbTotWeightsp );
1388                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
1389                                         } else {
1390                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
1391                                         }
1392                                         QCELL(workLev,ni,nj,nk, workSet, dMass) += change;
1393                                 }
1394                         }
1395                         massChange = 0.0;
1396                 } else {
1397                         // Problem! no interface neighbors
1398                         this->mFixMass += massChange;
1399                         //TTT mNumProblems++;
1400                         //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<<this->mFixMass);
1401                 }
1402                 weightIndex++;
1403
1404     // already done? RFLAG(workLev,i,j,k, workSet) = CFFluid;
1405     QCELL(workLev,i,j,k, workSet, dMass) = myrho; // should be rho... but unused?
1406     QCELL(workLev,i,j,k, workSet, dFfrac) = 1.0; // should be rho... but unused?
1407   } // fulllist
1408
1409
1410         /* now, finally handle the empty cells - order is important, has to be after
1411          * full cell handling */
1412   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1413        iter != mListEmpty.end(); iter++ ) {
1414     int i=iter->x, j=iter->y, k=iter->z;
1415
1416     LbmFloat massChange = QCELL(workLev, i,j,k, workSet, dMass);
1417                 if(vWeights[weightIndex].numNbs>0.0) {
1418                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
1419                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
1420                         FORDF1 {
1421                                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1422         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1423                                         LbmFloat change = -1.0;
1424                                         if(nbTotWeightsp>0.0) {
1425                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
1426                                         } else {
1427                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
1428                                         }
1429                                         QCELL(workLev, ni,nj,nk, workSet, dMass) += change;
1430                                 }
1431                         }
1432                         massChange = 0.0;
1433                 } else {
1434                         // Problem! no interface neighbors
1435                         this->mFixMass += massChange;
1436                         //TTT mNumProblems++;
1437                         //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<this->mFixMass);
1438                 }
1439                 weightIndex++;
1440                 
1441                 // finally... make it empty 
1442     // already done? RFLAG(workLev,i,j,k, workSet) = CFEmpty;
1443     QCELL(workLev,i,j,k, workSet, dMass) = 0.0;
1444     QCELL(workLev,i,j,k, workSet, dFfrac) = 0.0;
1445         }
1446   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1447        iter != mListEmpty.end(); iter++ ) {
1448     int i=iter->x, j=iter->y, k=iter->z;
1449     changeFlag(workLev,i,j,k, otherSet, CFEmpty);
1450         } 
1451
1452
1453         // check if some of the new interface cells can be removed again 
1454         // never happens !!! not necessary
1455         // calculate ffrac for new IF cells NEW
1456
1457         // how many are really new interface cells?
1458         int numNewIf = 0;
1459   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1460        iter != mListNewInter.end(); iter++ ) {
1461     int i=iter->x, j=iter->y, k=iter->z;
1462                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
1463                         continue; 
1464                         // FIXME remove from list?
1465                 }
1466                 numNewIf++;
1467         }
1468
1469         // redistribute mass, reinit flags
1470         if(debugFlagreinit) errMsg("NEWIF", "total:"<<mListNewInter.size());
1471         float newIfFac = 1.0/(LbmFloat)numNewIf;
1472   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1473        iter != mListNewInter.end(); iter++ ) {
1474     int i=iter->x, j=iter->y, k=iter->z;
1475                 if((i<=0) || (j<=0) || 
1476                          (i>=mLevel[workLev].lSizex-1) ||
1477                          (j>=mLevel[workLev].lSizey-1) ||
1478                          ((LBMDIM==3) && ((k<=0) || (k>=mLevel[workLev].lSizez-1) ) )
1479                          ) {
1480                         continue; } // new bc, dont treat cells on boundary NEWBC
1481                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
1482                         //errMsg("???"," "<<PRINT_IJK);
1483                         continue; 
1484                 } // */
1485
1486     QCELL(workLev,i,j,k, workSet, dMass) += (this->mFixMass * newIfFac);
1487                 
1488                 int nbored = 0;
1489                 FORDF1 { nbored |= RFLAG_NB(workLev, i,j,k, workSet,l); }
1490                 if(!(nbored & CFBndNoslip)) { RFLAG(workLev,i,j,k, workSet) |= CFNoBndFluid; }
1491                 if(!(nbored & CFFluid))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbFluid; }
1492                 if(!(nbored & CFEmpty))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbEmpty; }
1493
1494                 if(!(RFLAG(workLev,i,j,k, otherSet)&CFInter)) {
1495                         RFLAG(workLev,i,j,k, workSet) = (CellFlagType)(RFLAG(workLev,i,j,k, workSet) | CFNoDelete);
1496                 }
1497                 if(debugFlagreinit) errMsg("NEWIF", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" f"<< convertCellFlagType2String(RFLAG(workLev,i,j,k, workSet))<<" wl"<<workLev );
1498         }
1499
1500         // reinit fill fraction
1501   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1502        iter != mListNewInter.end(); iter++ ) {
1503     int i=iter->x, j=iter->y, k=iter->z;
1504                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { continue; }
1505
1506                 LbmFloat nrho = 0.0;
1507                 FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); }
1508     QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho;
1509     QCELL(workLev,i,j,k, workSet, dFlux) = FLUX_INIT;
1510         }
1511
1512         if(mListNewInter.size()>0){ 
1513                 //errMsg("FixMassDisted"," fm:"<<this->mFixMass<<" nif:"<<mListNewInter.size() );
1514                 this->mFixMass = 0.0; 
1515         }
1516
1517         // empty lists for next step
1518         mListFull.clear();
1519         mListEmpty.clear();
1520         mListNewInter.clear();
1521 } // reinitFlags
1522
1523
1524