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