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