Silence some annoying warnings when doing full build with strict flags
[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(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         GRID_REGION_INIT();
380 #if PARALLEL==1
381 #pragma omp parallel default(shared) num_threads(mNumOMPThreads) \
382   reduction(+: \
383           calcCurrentMass,calcCurrentVolume, \
384                 calcCellsFilled,calcCellsEmptied, \
385                 calcNumUsedCells )
386         GRID_REGION_START();
387 #else // PARALLEL==1
388         GRID_REGION_START();
389 #endif // PARALLEL==1
390
391         // local to main
392         CellFlagType nbflag[LBM_DFNUM]; 
393         int oldFlag, newFlag, nbored;
394         LbmFloat m[LBM_DFNUM];
395         LbmFloat rho, ux, uy, uz, tmp, usqr;
396
397         // smago vars
398         LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
399         
400         // ifempty cell conversion flags
401         bool iffilled, ifemptied;
402         LbmFloat nbfracs[LBM_DFNUM]; // ffracs of neighbors
403         int recons[LBM_DFNUM];   // reconstruct this DF?
404         int numRecons;           // how many are reconstructed?
405
406         LbmFloat mass=0., change=0., lcsmqo=0.;
407         rho= ux= uy= uz= usqr= tmp= 0.; 
408         lcsmqadd = lcsmomega = 0.;
409         FORDF0{ lcsmeq[l] = 0.; }
410
411         // ---
412         // now stream etc.
413         // use //template functions for 2D/3D
414
415         GRID_LOOP_START();
416                 // loop start
417                 // stream from current set to other, then collide and store
418                 //errMsg("l2"," at "<<PRINT_IJK<<" id"<<id);
419
420 #               if FSGR_STRICT_DEBUG==1
421                 // safety pointer check
422                 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
423                 if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) || 
424                     (&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) {
425                         errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<<
426                                         RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<<
427                                         (*pFlagSrc)<<","<<(*pFlagDst) <<",  pointers "<<
428           (long)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(long)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
429           (long)(pFlagSrc)<<","<<(long)(pFlagDst)<<" "
430                                         ); 
431                         CAUSE_PANIC;
432                 }       
433                 if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) || 
434                     (&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) {
435                         errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<<
436           (long)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(long)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
437           (long)(ccel)<<","<<(long)(tcel)<<" "
438                                         ); 
439                         CAUSE_PANIC;
440                 }       
441 #               endif
442                 oldFlag = *pFlagSrc;
443                 
444                 // old INTCFCOARSETEST==1
445                 if( (oldFlag & (CFGrFromCoarse)) ) { 
446                         if(( mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
447                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
448                         } else {
449                                 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
450                                 calcNumUsedCells++;
451                         }
452                         continue; // interpolateFineFromCoarse test!
453                 } // interpolateFineFromCoarse test! 
454         
455                 if(oldFlag & (CFMbndInflow)) {
456                         // fluid & if are ok, fill if later on
457                         int isValid = oldFlag & (CFFluid|CFInter);
458                         const LbmFloat iniRho = 1.0;
459                         const int OId = oldFlag>>24;
460                         if(!isValid) {
461                                 // make new if cell
462                                 const LbmVec vel(mObjectSpeeds[OId]);
463                                 // TODO add OPT3D treatment
464                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
465                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
466                                 RAC(tcel, dFlux) = FLUX_INIT;
467                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
468                                 calcCurrentMass += iniRho; 
469                                 calcCurrentVolume += 1.0; 
470                                 calcNumUsedCells++;
471                                 mInitialMass += iniRho;
472                                 // dont treat cell until next step
473                                 continue;
474                         } 
475                 } 
476                 else  // these are exclusive
477                 if(oldFlag & (CFMbndOutflow)) {
478                         int isnotValid = oldFlag & (CFFluid);
479                         if(isnotValid) {
480                                 // remove fluid cells, shouldnt be here anyway
481                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
482                                 mInitialMass -= fluidRho;
483                                 const LbmFloat iniRho = 0.0;
484                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
485                                 RAC(tcel, dFlux) = FLUX_INIT;
486                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
487
488                                 // same as ifemptied for if below
489                                 LbmPoint oemptyp; oemptyp.flag = 0;
490                                 oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
491                                 LIST_EMPTY(oemptyp);
492                                 calcCellsEmptied++;
493                                 continue;
494                         }
495                 }
496
497                 if(oldFlag & (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)) { 
498                         *pFlagDst = oldFlag;
499                         continue;
500                 }
501                 /*if( oldFlag & CFNoBndFluid ) {  // TEST ME FASTER?
502                         OPTIMIZED_STREAMCOLLIDE; PERFORM_USQRMAXCHECK;
503                         RAC(tcel,dFfrac) = 1.0; 
504                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
505                         calcCurrentMass += rho; calcCurrentVolume += 1.0;
506                         calcNumUsedCells++;
507                         continue;
508                 }// TEST ME FASTER? */
509
510                 // only neighbor flags! not own flag
511                 nbored = 0;
512                 
513 #if OPT3D==0
514                 FORDF1 {
515                         nbflag[l] = RFLAG_NB(lev, i,j,k,SRCS(lev),l);
516                         nbored |= nbflag[l];
517                 } 
518 #else
519                 nbflag[dSB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dSB];
520                 nbflag[dWB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWB];
521                 nbflag[ dB] = *(pFlagSrc + (-mLevel[lev].lOffsy)); nbored |= nbflag[dB];
522                 nbflag[dEB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dEB];
523                 nbflag[dNB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNB];
524
525                 nbflag[dSW] = *(pFlagSrc + (-mLevel[lev].lOffsx+-1)); nbored |= nbflag[dSW];
526                 nbflag[ dS] = *(pFlagSrc + (-mLevel[lev].lOffsx)); nbored |= nbflag[dS];
527                 nbflag[dSE] = *(pFlagSrc + (-mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dSE];
528
529                 nbflag[ dW] = *(pFlagSrc + (-1)); nbored |= nbflag[dW];
530                 nbflag[ dE] = *(pFlagSrc + ( 1)); nbored |= nbflag[dE];
531
532                 nbflag[dNW] = *(pFlagSrc + ( mLevel[lev].lOffsx+-1)); nbored |= nbflag[dNW];
533           nbflag[ dN] = *(pFlagSrc + ( mLevel[lev].lOffsx)); nbored |= nbflag[dN];
534                 nbflag[dNE] = *(pFlagSrc + ( mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dNE];
535
536                 nbflag[dST] = *(pFlagSrc + ( mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dST];
537                 nbflag[dWT] = *(pFlagSrc + ( mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWT];
538                 nbflag[ dT] = *(pFlagSrc + ( mLevel[lev].lOffsy)); nbored |= nbflag[dT];
539                 nbflag[dET] = *(pFlagSrc + ( mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dET];
540                 nbflag[dNT] = *(pFlagSrc + ( mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNT];
541                 // */
542 #endif
543
544                 // pointer to destination cell
545                 calcNumUsedCells++;
546
547                 // FLUID cells 
548                 if( oldFlag & CFFluid ) { 
549                         // only standard fluid cells (with nothing except fluid as nbs
550
551                         if(oldFlag&CFMbndInflow) {
552                                 // force velocity for inflow, necessary to have constant direction of flow
553                                 // FIXME , test also set interface cells?
554                                 const int OId = oldFlag>>24;
555                                 //? DEFAULT_STREAM;
556                                 //const LbmFloat fluidRho = 1.0;
557                                 // for submerged inflows, streaming would have to be performed...
558                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
559                                 const LbmVec vel(mObjectSpeeds[OId]);
560                                 ux=vel[0], uy=vel[1], uz=vel[2]; 
561                                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
562                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); }
563                                 rho = fluidRho;
564                                 //errMsg("INFLOW_DEBUG","std at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
565                         } else {
566                                 if(nbored&CFBnd) {
567                                         DEFAULT_STREAM;
568                                         //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
569                                         DEFAULT_COLLIDEG(mLevel[lev].gravity);
570                                         oldFlag &= (~CFNoBndFluid);
571                                 } else {
572                                         // do standard stream/collide
573                                         OPTIMIZED_STREAMCOLLIDE;
574                                         oldFlag |= CFNoBndFluid;
575                                 } 
576                         }
577
578                         PERFORM_USQRMAXCHECK;
579                         // "normal" fluid cells
580                         RAC(tcel,dFfrac) = 1.0; 
581                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
582                         calcCurrentMass += rho; 
583                         calcCurrentVolume += 1.0;
584                         continue;
585                 }
586                 
587                 newFlag  = oldFlag;
588                 // make sure here: always check which flags to really unset...!
589                 newFlag = newFlag & (~( 
590                                         CFNoNbFluid|CFNoNbEmpty| CFNoDelete 
591                                         | CFNoInterpolSrc
592                                         | CFNoBndFluid
593                                         ));
594                 if(!(nbored&CFBndNoslip)) { //NEWSURFT NEWSURFTNOS
595                         newFlag |= CFNoBndFluid;
596                 }
597                 /*if(!(nbored&CFBnd)) { //NEWSURFT NEWSURFTNOS
598                         // explicitly check for noslip neighbors
599                         bool hasnoslipnb = false;
600                         FORDF1 { if((nbflag[l]&CFBnd)&&(nbflag[l]&CFBndNoslip)) hasnoslipnb=true; }
601                         if(!hasnoslipnb) newFlag |= CFNoBndFluid; 
602                 } // */
603
604                 // store own dfs and mass
605                 mass = RAC(ccel,dMass);
606
607                 // WARNING - only interface cells arrive here!
608                 // read distribution funtions of adjacent cells = stream step
609                 DEFAULT_STREAM;
610
611                 if((nbored & CFFluid)==0) { newFlag |= CFNoNbFluid; mNumInvIfCells++; }
612                 if((nbored & CFEmpty)==0) { newFlag |= CFNoNbEmpty; mNumInvIfCells++; }
613
614                 // calculate mass exchange for interface cells 
615                 LbmFloat myfrac = RAC(ccel,dFfrac);
616                 if(myfrac<0.) myfrac=0.; // NEWSURFT
617 #               define nbdf(l) m[ this->dfInv[(l)] ]
618
619                 // update mass 
620                 // only do boundaries for fluid cells, and interface cells without
621                 // any fluid neighbors (assume that interface cells _with_ fluid
622                 // neighbors are affected enough by these) 
623                 // which Df's have to be reconstructed? 
624                 // for fluid cells - just the f_i difference from streaming to empty cells  ----
625                 numRecons = 0;
626                 bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
627                 //onlyBndnb = false; // DEBUG test off
628
629                 FORDF1 { // dfl loop
630                         recons[l] = 0;
631                         nbfracs[l] = 0.0;
632                         // finally, "normal" interface cells ----
633                         if( nbflag[l]&(CFFluid|CFBnd) ) { // NEWTEST! FIXME check!!!!!!!!!!!!!!!!!!
634                                 change = nbdf(l) - MYDF(l);
635                         }
636                         // interface cells - distuingish cells that shouldn't fill/empty 
637                         else if( nbflag[l] & CFInter ) {
638                                 
639                                 LbmFloat mynbfac,nbnbfac;
640                                 // NEW TEST t1
641                                 // t2 -> off
642                                 if((oldFlag&CFNoBndFluid)&&(nbflag[l]&CFNoBndFluid)) {
643                                         mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
644                                         nbnbfac = 1.0/mynbfac;
645                                         onlyBndnb = false;
646                                 } else {
647                                         mynbfac = nbnbfac = 1.0; // switch calc flux off
648                                         goto changeDefault;  // NEWSURFT
649                                         //change = 0.; goto changeDone;  // NEWSURFT
650                                 }
651                                 //mynbfac = nbnbfac = 1.0; // switch calc flux off t3
652
653                                 // perform interface case handling
654                                 if ((oldFlag|nbflag[l])&(CFNoNbFluid|CFNoNbEmpty)) {
655                                 switch (oldFlag&(CFNoNbFluid|CFNoNbEmpty)) {
656                                         case 0: 
657                                                 // we are a normal cell so... 
658                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
659                                                         case CFNoNbFluid: 
660                                                                 // just fill current cell = empty neighbor 
661                                                                 change = nbnbfac*nbdf(l) ; goto changeDone; 
662                                                         case CFNoNbEmpty: 
663                                                                 // just empty current cell = fill neighbor 
664                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
665                                                 }
666                                                 break;
667
668                                         case CFNoNbFluid: 
669                                                 // we dont have fluid nb's so...
670                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
671                                                         case 0: 
672                                                         case CFNoNbEmpty: 
673                                                                 // we have no fluid nb's -> just empty
674                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
675                                                 }
676                                                 break;
677
678                                         case CFNoNbEmpty: 
679                                                 // we dont have empty nb's so...
680                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
681                                                         case 0: 
682                                                         case CFNoNbFluid: 
683                                                                 // we have no empty nb's -> just fill
684                                                                 change = nbnbfac*nbdf(l); goto changeDone; 
685                                                 }
686                                                 break;
687                                 }} // inter-inter exchange
688
689                         changeDefault: ;
690                                 // just do normal mass exchange...
691                                 change = ( nbnbfac*nbdf(l) - mynbfac*MYDF(l) ) ;
692                         changeDone: ;
693                                 nbfracs[l] = QCELL_NB(lev, i,j,k, SRCS(lev),l, dFfrac);
694                                 if(nbfracs[l]<0.) nbfracs[l] = 0.; // NEWSURFT
695                                 change *=  (myfrac + nbfracs[l]) * 0.5;
696                         } // the other cell is interface
697
698                         // last alternative - reconstruction in this direction
699                         else {
700                                 // empty + bnd case
701                                 recons[l] = 1; 
702                                 numRecons++;
703                                 change = 0.0; 
704                         }
705
706                         // modify mass at SRCS
707                         mass += change;
708                 } // l
709                 // normal interface, no if empty/fluid
710
711                 // computenormal
712                 LbmFloat surfaceNormal[3];
713                 computeFluidSurfaceNormal(ccel,pFlagSrc, surfaceNormal);
714
715                 if( (ABS(surfaceNormal[0])+ABS(surfaceNormal[1])+ABS(surfaceNormal[2])) > LBM_EPSILON) {
716                         // normal ok and usable...
717                         FORDF1 {
718                                 if( (this->dfDvecX[l]*surfaceNormal[0] + this->dfDvecY[l]*surfaceNormal[1] + this->dfDvecZ[l]*surfaceNormal[2])  // dot Dvec,norml
719                                                 > LBM_EPSILON) {
720                                         recons[l] = 2; 
721                                         numRecons++;
722                                 } 
723                         }
724                 }
725
726                 // calculate macroscopic cell values
727                 LbmFloat oldUx, oldUy, oldUz;
728                 LbmFloat oldRho; // OLD rho = ccel->rho;
729 #               define REFERENCE_PRESSURE 1.0 // always atmosphere...
730 #               if OPT3D==0
731                 oldRho=RAC(ccel,0);
732                 oldUx = oldUy = oldUz = 0.0;
733                 for(int l=1; l<this->cDfNum; l++) {
734                         oldRho += RAC(ccel,l);
735                         oldUx  += (this->dfDvecX[l]*RAC(ccel,l));
736                         oldUy  += (this->dfDvecY[l]*RAC(ccel,l)); 
737                         oldUz  += (this->dfDvecZ[l]*RAC(ccel,l)); 
738                 } 
739                 // reconstruct dist funcs from empty cells
740                 FORDF1 {
741                         if(recons[ l ]) {
742                                 m[ this->dfInv[l] ] = 
743                                         this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 
744                                         this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 
745                                         - MYDF( l );
746                         }
747                 }
748                 ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
749                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
750 #               else // OPT3D==0
751                 oldRho = + RAC(ccel,dC)  + RAC(ccel,dN )
752                                 + RAC(ccel,dS ) + RAC(ccel,dE )
753                                 + RAC(ccel,dW ) + RAC(ccel,dT )
754                                 + RAC(ccel,dB ) + RAC(ccel,dNE)
755                                 + RAC(ccel,dNW) + RAC(ccel,dSE)
756                                 + RAC(ccel,dSW) + RAC(ccel,dNT)
757                                 + RAC(ccel,dNB) + RAC(ccel,dST)
758                                 + RAC(ccel,dSB) + RAC(ccel,dET)
759                                 + RAC(ccel,dEB) + RAC(ccel,dWT)
760                                 + RAC(ccel,dWB);
761
762                 oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
763                                 + RAC(ccel,dNE) - RAC(ccel,dNW)
764                                 + RAC(ccel,dSE) - RAC(ccel,dSW)
765                                 + RAC(ccel,dET) + RAC(ccel,dEB)
766                                 - RAC(ccel,dWT) - RAC(ccel,dWB);
767
768                 oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
769                                 + RAC(ccel,dNE) + RAC(ccel,dNW)
770                                 - RAC(ccel,dSE) - RAC(ccel,dSW)
771                                 + RAC(ccel,dNT) + RAC(ccel,dNB)
772                                 - RAC(ccel,dST) - RAC(ccel,dSB);
773
774                 oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
775                                 + RAC(ccel,dNT) - RAC(ccel,dNB)
776                                 + RAC(ccel,dST) - RAC(ccel,dSB)
777                                 + RAC(ccel,dET) - RAC(ccel,dEB)
778                                 + RAC(ccel,dWT) - RAC(ccel,dWB);
779
780                 (void)oldRho;
781
782                 // now reconstruction
783                 ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
784                 rho = REFERENCE_PRESSURE;
785                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
786                 if(recons[dN ]) { m[dS ] = EQN  + EQS  - MYDF(dN ); }
787                 if(recons[dS ]) { m[dN ] = EQS  + EQN  - MYDF(dS ); }
788                 if(recons[dE ]) { m[dW ] = EQE  + EQW  - MYDF(dE ); }
789                 if(recons[dW ]) { m[dE ] = EQW  + EQE  - MYDF(dW ); }
790                 if(recons[dT ]) { m[dB ] = EQT  + EQB  - MYDF(dT ); }
791                 if(recons[dB ]) { m[dT ] = EQB  + EQT  - MYDF(dB ); }
792                 if(recons[dNE]) { m[dSW] = EQNE + EQSW - MYDF(dNE); }
793                 if(recons[dNW]) { m[dSE] = EQNW + EQSE - MYDF(dNW); }
794                 if(recons[dSE]) { m[dNW] = EQSE + EQNW - MYDF(dSE); }
795                 if(recons[dSW]) { m[dNE] = EQSW + EQNE - MYDF(dSW); }
796                 if(recons[dNT]) { m[dSB] = EQNT + EQSB - MYDF(dNT); }
797                 if(recons[dNB]) { m[dST] = EQNB + EQST - MYDF(dNB); }
798                 if(recons[dST]) { m[dNB] = EQST + EQNB - MYDF(dST); }
799                 if(recons[dSB]) { m[dNT] = EQSB + EQNT - MYDF(dSB); }
800                 if(recons[dET]) { m[dWB] = EQET + EQWB - MYDF(dET); }
801                 if(recons[dEB]) { m[dWT] = EQEB + EQWT - MYDF(dEB); }
802                 if(recons[dWT]) { m[dEB] = EQWT + EQEB - MYDF(dWT); }
803                 if(recons[dWB]) { m[dET] = EQWB + EQET - MYDF(dWB); }
804 #               endif           
805
806
807                 // inflow bc handling
808                 if(oldFlag & (CFMbndInflow)) {
809                         // fill if cells in inflow region
810                         if(myfrac<0.5) { 
811                                 mass += 0.25; 
812                                 mInitialMass += 0.25;
813                         }
814                         const int OId = oldFlag>>24;
815                         const LbmVec vel(mObjectSpeeds[OId]);
816                         ux=vel[0], uy=vel[1], uz=vel[2]; 
817                         //? usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
818                         //FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); } rho = fluidRho;
819                         rho = REFERENCE_PRESSURE;
820                         FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
821                         //errMsg("INFLOW_DEBUG","if at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
822                 }  else { 
823                         // NEWSURFT, todo optimize!
824                         if(onlyBndnb) { //if(usqr<0.001*0.001) {
825                                 rho = ux = uy = uz = 0.;
826                                 FORDF0{ 
827                                         rho += m[l]; 
828                                         ux  += (this->dfDvecX[l]*m[l]); 
829                                         uy  += (this->dfDvecY[l]*m[l]);  
830                                         uz  += (this->dfDvecZ[l]*m[l]);  
831                                 }
832                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
833                         } else {// NEWSURFT */
834                                 if(usqr>0.3*0.3) { 
835                                         // force reset! , warning causes distortions...
836                                         FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,0.,0.,0.); }
837                                 } else {
838                                 // normal collide
839                                 // mass streaming done... do normal collide
840                                 LbmVec grav = mLevel[lev].gravity*mass;
841                                 DEFAULT_COLLIDEG(grav);
842                                 PERFORM_USQRMAXCHECK; }
843                                 // rho init from default collide necessary for fill/empty check below
844                         } // test
845                 }
846
847                 // testing..., particle generation
848                 // also check oldFlag for CFNoNbFluid, test
849                 // for inflow no pargen test // NOBUBBB!
850                 if((mInitDone) 
851                                 // dont allow new if cells, or submerged ones
852                                 && (!((oldFlag|newFlag)& (CFNoDelete|CFNoNbEmpty) )) 
853                                 // dont try to subtract from empty cells
854                                 && (mass>0.) && (mPartGenProb>0.0)) {
855                         bool doAdd = true;
856                         bool bndOk=true;
857                         if( (i<cutMin)||(i>mSizex-cutMin)||
858                                         (j<cutMin)||(j>mSizey-cutMin)||
859                                         (k<cutMin)||(k>mSizez-cutMin) ) { bndOk=false; }
860                         if(!bndOk) doAdd=false;
861                         
862                         LbmFloat prob = (rand()/(RAND_MAX+1.0));
863                         LbmFloat basethresh = mPartGenProb*lcsmqo*(LbmFloat)(mSizez+mSizey+mSizex)*0.5*0.333;
864
865                         // physical drop model
866                         if(mPartUsePhysModel) {
867                                 LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep);
868                                 LbmFloat rux = (ux * realWorldFac);
869                                 LbmFloat ruy = (uy * realWorldFac);
870                                 LbmFloat ruz = (uz * realWorldFac);
871                                 LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz));
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 = LbmVec(ux,uy,uz);
964                                 LbmFloat flvelLen = norm(flvelVel);
965                                 // surface normal
966                                 LbmVec normVel = LbmVec(surfaceNormal[0],surfaceNormal[1],surfaceNormal[2]);
967                                 normalize(normVel);
968                                 LbmFloat normScale = (0.01+flvelLen);
969                                 // jitter vector, 0.2 * flvel
970                                 LbmVec jittVel = LbmVec(jx,jy,jz)*(0.05+flvelLen)*0.1;
971                                 // weighten velocities
972                                 const LbmFloat flvelWeight = 0.9;
973                                 LbmVec newpartVel = normVel*normScale*(1.-flvelWeight) + flvelVel*(flvelWeight) + jittVel; 
974
975                                 // offset towards surface (hide popping)
976                                 jpx += -normVel[0]*0.4;
977                                 jpy += -normVel[1]*0.4;
978                                 jpz += -normVel[2]*0.4;
979
980                                 LbmFloat srci=i+0.5+jpx, srcj=j+0.5+jpy, srck=k+0.5+jpz;
981                                 int type=0;
982                                 type = PART_DROP;
983
984 #                               if LBMDIM==2
985                                 newpartVel[2]=0.; srck=0.5;
986 #                               endif // LBMDIM==2
987                                 // subtract drop mass
988                                 mass -= dropmass;
989                                 // init new particle
990                                 {
991                                         ParticleObject np( ntlVec3Gfx(srci,srcj,srck) );
992                                         np.setVel(newpartVel[0],newpartVel[1],newpartVel[2]);
993                                         np.setStatus(PART_IN);
994                                         np.setType(type);
995                                         //if((s%3)==2) np.setType(PART_FLOAT);
996                                         np.setSize(newPartsize);
997                                         //errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<norm(LbmVec(ux,uy,uz)) <<" add"<<doAdd<<" pvel="<<norm(newpartVel)<<" size="<<newPartsize );
998                                         //errMsg("NEWPT","u="<<newpartVel<<" norm="<<normVel<<" flvel="<<flvelVel<<" jitt="<<jittVel );
999                                         FSGR_ADDPART(np);
1000                                 } // new part
1001                 //errMsg("PIT","NEW pit"<<np.getId()<<" pos:"<< np.getPos()<<" status:"<<convertFlags2String(np.getFlags())<<" vel:"<< np.getVel()  );
1002                                 } // multiple parts
1003                         } // doAdd
1004                 } // */
1005
1006
1007                 // interface cell filled or emptied?
1008                 iffilled = ifemptied = 0;
1009                 // interface cells empty/full?, WARNING: to mark these cells, better do it at the end of reinitCellFlags
1010                 // interface cell if full?
1011                 if( (mass) >= (rho * (1.0+FSGR_MAGICNR)) ) { iffilled = 1; }
1012                 // interface cell if empty?
1013                 if( (mass) <= (rho * (   -FSGR_MAGICNR)) ) { ifemptied = 1; }
1014
1015                 if(oldFlag & (CFMbndOutflow)) {
1016                         mInitialMass -= mass;
1017                         mass = myfrac = 0.0;
1018                         iffilled = 0; ifemptied = 1;
1019                 }
1020
1021                 // looks much nicer... LISTTRICK
1022 #               if FSGR_LISTTRICK==1
1023                 //if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
1024                 if(newFlag&CFNoBndFluid) { // test NEW TEST
1025                         if(!iffilled) {
1026                                 // remove cells independent from amount of change...
1027                                 if( (oldFlag & CFNoNbEmpty)&&(newFlag & CFNoNbEmpty)&&
1028                                                 ( (mass>(rho*FSGR_LISTTTHRESHFULL))  || ((nbored&CFInter)==0)  )) { 
1029                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","filled "<<PRINT_IJK); };
1030                                         iffilled = 1; 
1031                                 } 
1032                         }
1033                         if(!ifemptied) {
1034                                 if( (oldFlag & CFNoNbFluid)&&(newFlag & CFNoNbFluid)&&
1035                                                 ( (mass<(rho*FSGR_LISTTTHRESHEMPTY)) || ((nbored&CFInter)==0)  )) { 
1036                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","emptied "<<PRINT_IJK); };
1037                                         ifemptied = 1; 
1038                                 } 
1039                         }
1040                 } // nobndfluid only */
1041 #               endif
1042                 //iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!!
1043                 
1044
1045                 // now that all dfs are known, handle last changes
1046                 if(iffilled) {
1047                         LbmPoint filledp; filledp.flag=0;
1048                         if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1;  // NEWSURFT
1049                         filledp.x = i; filledp.y = j; filledp.z = k;
1050                         LIST_FULL(filledp);
1051                         //mNumFilledCells++; // DEBUG
1052                         calcCellsFilled++;
1053                 }
1054                 else if(ifemptied) {
1055                         LbmPoint emptyp; emptyp.flag=0;
1056                         if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; //  NEWSURFT
1057                         emptyp.x = i; emptyp.y = j; emptyp.z = k;
1058                         LIST_EMPTY(emptyp);
1059                         //mNumEmptiedCells++; // DEBUG
1060                         calcCellsEmptied++;
1061                 } 
1062                 // dont cutoff values -> better cell conversions
1063                 RAC(tcel,dFfrac)   = (mass/rho);
1064
1065                 // init new flux value
1066                 float flux = FLUX_INIT; // dxqn on
1067                 if(newFlag&CFNoBndFluid) {
1068                         //flux = 50.0; // extreme on
1069                         for(int nn=1; nn<this->cDfNum; nn++) { 
1070                                 if(nbflag[nn] & (CFFluid|CFInter|CFBnd)) { flux += this->dfLength[nn]; }
1071                         }
1072                         // optical hack - smooth slow moving
1073                         // surface regions
1074                         if(usqr< sssUsqrLimit) {
1075                         for(int nn=1; nn<this->cDfNum; nn++) { 
1076                                 if(nbfracs[nn]!=0.0) {
1077                                         LbmFloat curSmooth = (sssUsqrLimit-usqr)*sssUsqrLimitInv;
1078                                         if(curSmooth>1.0) curSmooth=1.0;
1079                                         flux *= (1.0+ smoothStrength*curSmooth * (nbfracs[nn]-myfrac)) ;
1080                                 }
1081                         } }
1082                         // NEW TEST */
1083                 }
1084                 // flux = FLUX_INIT; // calc flux off
1085                 QCELL(lev, i,j,k,TSET(lev), dFlux) = flux; // */
1086
1087                 // perform mass exchange with streamed values
1088                 QCELL(lev, i,j,k,TSET(lev), dMass) = mass; // MASST
1089                 // set new flag 
1090                 *pFlagDst = (CellFlagType)newFlag;
1091                 calcCurrentMass += mass; 
1092                 calcCurrentVolume += RAC(tcel,dFfrac);
1093
1094                 // interface cell handling done...
1095
1096 #if PARALLEL!=1
1097         GRID_LOOPREG_END();
1098 #else // PARALLEL==1
1099 #include "paraloopend.h" // = GRID_LOOPREG_END();
1100 #endif // PARALLEL==1
1101
1102         // write vars from computations to class
1103         mLevel[lev].lmass    = calcCurrentMass;
1104         mLevel[lev].lvolume  = calcCurrentVolume;
1105         mNumFilledCells  = calcCellsFilled;
1106         mNumEmptiedCells = calcCellsEmptied;
1107         mNumUsedCells = calcNumUsedCells;
1108 }
1109
1110
1111
1112 void 
1113 LbmFsgrSolver::preinitGrids()
1114 {
1115         const int lev = mMaxRefine;
1116         const bool doReduce = false;
1117         const int gridLoopBound=0;
1118
1119         // preinit both grids
1120         for(int s=0; s<2; s++) {
1121         
1122                 GRID_REGION_INIT();
1123 #if PARALLEL==1
1124 #pragma omp parallel default(shared) num_threads(mNumOMPThreads) \
1125   reduction(+: \
1126           calcCurrentMass,calcCurrentVolume, \
1127                 calcCellsFilled,calcCellsEmptied, \
1128                 calcNumUsedCells )
1129 #endif // PARALLEL==1
1130                 GRID_REGION_START();
1131                 GRID_LOOP_START();
1132                         for(int l=0; l<dTotalNum; l++) { RAC(ccel,l) = 0.; }
1133                         *pFlagSrc =0;
1134                         *pFlagDst =0;
1135                         //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
1136 #if PARALLEL!=1
1137                 GRID_LOOPREG_END();
1138 #else // PARALLEL==1
1139 #include "paraloopend.h" // = GRID_LOOPREG_END();
1140 #endif // PARALLEL==1
1141
1142                 /* dummy, remove warnings */ 
1143                 calcCurrentMass = calcCurrentVolume = 0.;
1144                 calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
1145                 
1146                 // change grid
1147                 mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
1148                 mLevel[mMaxRefine].setCurr   ^= 1;
1149         }
1150 }
1151
1152 void 
1153 LbmFsgrSolver::standingFluidPreinit()
1154 {
1155         const int lev = mMaxRefine;
1156         const bool doReduce = false;
1157         const int gridLoopBound=1;
1158
1159         GRID_REGION_INIT();
1160 #if PARALLEL==1
1161 #pragma omp parallel default(shared) num_threads(mNumOMPThreads) \
1162   reduction(+: \
1163           calcCurrentMass,calcCurrentVolume, \
1164                 calcCellsFilled,calcCellsEmptied, \
1165                 calcNumUsedCells )
1166 #endif // PARALLEL==1
1167         GRID_REGION_START();
1168
1169         LbmFloat rho, ux,uy,uz, usqr; 
1170         CellFlagType nbflag[LBM_DFNUM];
1171         LbmFloat m[LBM_DFNUM];
1172         LbmFloat lcsmqo;
1173 #       if OPT3D==1 
1174         LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
1175 #       endif // OPT3D==true 
1176
1177         GRID_LOOP_START();
1178                 //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
1179                 const CellFlagType currFlag = *pFlagSrc; //RFLAG(lev, i,j,k,workSet);
1180                 if( (currFlag & (CFEmpty|CFBnd)) ) continue;
1181
1182                 if( (currFlag & (CFInter)) ) {
1183                         // copy all values
1184                         for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
1185                         continue;
1186                 }
1187
1188                 if( (currFlag & CFNoBndFluid)) {
1189                         OPTIMIZED_STREAMCOLLIDE;
1190                 } else {
1191                         FORDF1 {
1192                                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
1193                         } 
1194                         DEFAULT_STREAM;
1195                         DEFAULT_COLLIDEG(mLevel[lev].gravity);
1196                 }
1197                 for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
1198 #if PARALLEL!=1
1199         GRID_LOOPREG_END();
1200 #else // PARALLEL==1
1201 #include "paraloopend.h" // = GRID_LOOPREG_END();
1202 #endif // PARALLEL==1
1203
1204         /* dummy remove warnings */ 
1205         calcCurrentMass = calcCurrentVolume = 0.;
1206         calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
1207         
1208         // change grid
1209   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
1210   mLevel[mMaxRefine].setCurr   ^= 1;
1211 }
1212
1213
1214 /******************************************************************************
1215  * work on lists from updateCellMass to reinit cell flags
1216  *****************************************************************************/
1217
1218 LbmFloat LbmFsgrSolver::getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l) {
1219         //return 0.0; // test
1220         int level = mMaxRefine;
1221         LbmFloat     *ccel  = RACPNT(level, i,j,k, workSet);
1222
1223         // computenormal
1224         CellFlagType *cflag = &RFLAG(level, i,j,k, workSet);
1225         LbmFloat n[3];
1226         computeFluidSurfaceNormal(ccel,cflag, n);
1227         LbmFloat scal = mDvecNrm[l][0]*n[0] + mDvecNrm[l][1]*n[1] + mDvecNrm[l][2]*n[2];
1228
1229         LbmFloat ret = 1.0;
1230         // forward direction, add mass (for filling cells):
1231         if(dirForw) {
1232                 if(scal<LBM_EPSILON) ret = 0.0;
1233                 else ret = scal;
1234         } else {
1235                 // backward for emptying
1236                 if(scal>-LBM_EPSILON) ret = 0.0;
1237                 else ret = scal * -1.0;
1238         }
1239         //errMsg("massd", PRINT_IJK<<" nv"<<nvel<<" : ret="<<ret ); //xit(1); //VECDEB
1240         return ret;
1241 }
1242
1243 // warning - normal compuations are without
1244 //   boundary checks &
1245 //   normalization
1246 void LbmFsgrSolver::computeFluidSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt,LbmFloat *snret) {
1247         const int level = mMaxRefine;
1248         LbmFloat nx,ny,nz, nv1,nv2;
1249         const CellFlagType flagE = *(cflagpnt+1);
1250         const CellFlagType flagW = *(cflagpnt-1);
1251         if(flagE &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } 
1252         else if(flagE &(CFBnd)){ nv1 = 1.; }
1253         else nv1 = 0.0;
1254         if(flagW &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } 
1255         else if(flagW &(CFBnd)){ nv2 = 1.; }
1256         else nv2 = 0.0;
1257         nx = 0.5* (nv2-nv1);
1258
1259         const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
1260         const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
1261         if(flagN &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 
1262         else if(flagN &(CFBnd)){ nv1 = 1.; }
1263         else nv1 = 0.0;
1264         if(flagS &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 
1265         else if(flagS &(CFBnd)){ nv2 = 1.; }
1266         else nv2 = 0.0;
1267         ny = 0.5* (nv2-nv1);
1268
1269 #if LBMDIM==3
1270         const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
1271         const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
1272         if(flagT &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 
1273         else if(flagT &(CFBnd)){ nv1 = 1.; }
1274         else nv1 = 0.0;
1275         if(flagB &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 
1276         else if(flagB &(CFBnd)){ nv2 = 1.; }
1277         else nv2 = 0.0;
1278         nz = 0.5* (nv2-nv1);
1279 #else //LBMDIM==3
1280         nz = 0.0;
1281 #endif //LBMDIM==3
1282
1283         // return vals
1284         snret[0]=nx; snret[1]=ny; snret[2]=nz;
1285 }
1286 void LbmFsgrSolver::computeFluidSurfaceNormalAcc(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
1287         LbmFloat nx=0.,ny=0.,nz=0.;
1288         ccel = NULL; cflagpnt=NULL; // remove warning
1289         snret[0]=nx; snret[1]=ny; snret[2]=nz;
1290 }
1291 void LbmFsgrSolver::computeObstacleSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
1292         const int level = mMaxRefine;
1293         LbmFloat nx,ny,nz, nv1,nv2;
1294         ccel = NULL; // remove warning
1295
1296         const CellFlagType flagE = *(cflagpnt+1);
1297         const CellFlagType flagW = *(cflagpnt-1);
1298         if(flagE &(CFBnd)){ nv1 = 1.; }
1299         else nv1 = 0.0;
1300         if(flagW &(CFBnd)){ nv2 = 1.; }
1301         else nv2 = 0.0;
1302         nx = 0.5* (nv2-nv1);
1303
1304         const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
1305         const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
1306         if(flagN &(CFBnd)){ nv1 = 1.; }
1307         else nv1 = 0.0;
1308         if(flagS &(CFBnd)){ nv2 = 1.; }
1309         else nv2 = 0.0;
1310         ny = 0.5* (nv2-nv1);
1311
1312 #if LBMDIM==3
1313         const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
1314         const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
1315         if(flagT &(CFBnd)){ nv1 = 1.; }
1316         else nv1 = 0.0;
1317         if(flagB &(CFBnd)){ nv2 = 1.; }
1318         else nv2 = 0.0;
1319         nz = 0.5* (nv2-nv1);
1320 #else //LBMDIM==3
1321         nz = 0.0;
1322 #endif //LBMDIM==3
1323
1324         // return vals
1325         snret[0]=nx; snret[1]=ny; snret[2]=nz;
1326 }
1327 void LbmFsgrSolver::computeObstacleSurfaceNormalAcc(int i,int j,int k, LbmFloat *snret) {
1328         bool nonorm = false;
1329         LbmFloat nx=0.,ny=0.,nz=0.;
1330         if(i<=0)        { nx =  1.; nonorm = true; }
1331         if(i>=mSizex-1) { nx = -1.; nonorm = true; }
1332         if(j<=0)        { ny =  1.; nonorm = true; }
1333         if(j>=mSizey-1) { ny = -1.; nonorm = true; }
1334 #       if LBMDIM==3
1335         if(k<=0)        { nz =  1.; nonorm = true; }
1336         if(k>=mSizez-1) { nz = -1.; nonorm = true; }
1337 #       endif // LBMDIM==3
1338         if(!nonorm) {
1339                 // in domain, revert to helper, use setCurr&mMaxRefine
1340                 LbmVec bnormal;
1341                 CellFlagType *bflag = &RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr);
1342                 LbmFloat     *bcell = RACPNT(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr);
1343                 computeObstacleSurfaceNormal(bcell,bflag, &bnormal[0]);
1344                 // TODO check if there is a normal near here?
1345                 // use wider range otherwise...
1346                 snret[0]=bnormal[0]; snret[1]=bnormal[0]; snret[2]=bnormal[0];
1347                 return;
1348         }
1349         snret[0]=nx; snret[1]=ny; snret[2]=nz;
1350 }
1351
1352 void LbmFsgrSolver::addToNewInterList( int ni, int nj, int nk ) {
1353 #if FSGR_STRICT_DEBUG==10
1354         // dangerous, this can change the simulation...
1355   /*for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1356        iter != mListNewInter.end(); iter++ ) {
1357     if(ni!=iter->x) continue;
1358     if(nj!=iter->y) continue;
1359     if(nk!=iter->z) continue;
1360                 // all 3 values match... skip point
1361                 return;
1362         } */
1363 #endif // FSGR_STRICT_DEBUG==1
1364         // store point
1365         LbmPoint newinter; newinter.flag = 0;
1366         newinter.x = ni; newinter.y = nj; newinter.z = nk;
1367         mListNewInter.push_back(newinter);
1368 }
1369
1370 void LbmFsgrSolver::reinitFlags( int workSet ) { 
1371         // reinitCellFlags OLD mods:
1372         // add all to intel list?
1373         // check ffrac for new cells
1374         // new if cell inits (last loop)
1375         // vweights handling
1376
1377         const int debugFlagreinit = 0;
1378         
1379         // some things need to be read/modified on the other set
1380         int otherSet = (workSet^1);
1381         // fixed level on which to perform 
1382         int workLev = mMaxRefine;
1383
1384   /* modify interface cells from lists */
1385   /* mark filled interface cells as fluid, or emptied as empty */
1386         /* count neighbors and distribute excess mass to interface neighbor cells
1387    * problems arise when there are no interface neighbors anymore
1388          * then just distribute to any fluid neighbors...
1389          */
1390
1391         // for symmetry, first init all neighbor cells */
1392         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1393        iter != mListFull.end(); iter++ ) {
1394     int i=iter->x, j=iter->y, k=iter->z;
1395                 if(debugFlagreinit) errMsg("FULL", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
1396     FORDF1 {
1397                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1398                         //if((LBMDIM>2)&&( (ni<=0) || (nj<=0) || (nk<=0) || (ni>=mLevel[workLev].lSizex-1) || (nj>=mLevel[workLev].lSizey-1) || (nk>=mLevel[workLev].lSizez-1) )) {
1399                         if( (ni<=0) || (nj<=0) || 
1400                                   (ni>=mLevel[workLev].lSizex-1) ||
1401                                   (nj>=mLevel[workLev].lSizey-1) 
1402 #                                       if LBMDIM==3
1403                                   || (nk<=0) || (nk>=mLevel[workLev].lSizez-1) 
1404 #                                       endif // LBMDIM==1
1405                                  ) {
1406                                 continue; } // new bc, dont treat cells on boundary NEWBC
1407       if( RFLAG(workLev, ni,nj,nk, workSet) & CFEmpty ){
1408                                 
1409                                 // preinit speed, get from average surrounding cells
1410                                 // interpolate from non-workset to workset, sets are handled in function
1411
1412                                 // new and empty interface cell, dont change old flag here!
1413                                 addToNewInterList(ni,nj,nk);
1414
1415                                 LbmFloat avgrho = 0.0;
1416                                 LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
1417                                 interpolateCellValues(workLev,ni,nj,nk,workSet, avgrho,avgux,avguy,avguz);
1418
1419                                 // careful with l's...
1420                                 FORDF0M { 
1421                                         QCELL(workLev,ni,nj,nk, workSet, m) = this->getCollideEq( m,avgrho,  avgux, avguy, avguz ); 
1422                                         //QCELL(workLev,ni,nj,nk, workSet, l) = avgnbdf[l]; // CHECK FIXME test?
1423                                 }
1424                                 //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
1425                                 QCELL(workLev,ni,nj,nk, workSet, dMass) = 0.0; //?? new
1426                                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) = 0.0; //?? new
1427                                 //RFLAG(workLev,ni,nj,nk,workSet) = (CellFlagType)(CFInter|CFNoInterpolSrc);
1428                                 changeFlag(workLev,ni,nj,nk,workSet, (CFInter|CFNoInterpolSrc));
1429                                 if(debugFlagreinit) errMsg("NEWE", PRINT_IJK<<" newif "<<PRINT_VEC(ni,nj,nk)<<" rho"<<avgrho<<" vel("<<avgux<<","<<avguy<<","<<avguz<<") " );
1430       }
1431                         /* prevent surrounding interface cells from getting removed as empty cells 
1432                          * (also cells that are not newly inited) */
1433       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1434                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete);
1435                                 changeFlag(workLev,ni,nj,nk, workSet, (RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete));
1436                                 // also add to list...
1437                                 addToNewInterList(ni,nj,nk);
1438                         } // NEW?
1439     }
1440
1441                 // NEW? no extra loop...
1442                 //RFLAG(workLev,i,j,k, workSet) = CFFluid;
1443                 changeFlag(workLev,i,j,k, workSet,CFFluid);
1444         }
1445
1446         /* remove empty interface cells that are not allowed to be removed anyway
1447          * this is important, otherwise the dreaded cell-type-flickering can occur! */
1448   for(int n=0; n<(int)mListEmpty.size(); n++) {
1449     int i=mListEmpty[n].x, j=mListEmpty[n].y, k=mListEmpty[n].z;
1450                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) {
1451                         // treat as "new inter"
1452                         addToNewInterList(i,j,k);
1453                         // remove entry
1454                         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
1455                         if(n<(int)mListEmpty.size()-1) mListEmpty[n] = mListEmpty[mListEmpty.size()-1]; 
1456                         mListEmpty.pop_back();
1457                         n--; 
1458                 }
1459         } 
1460
1461
1462         /* problems arise when adjacent cells empty&fill ->
1463                  let fill cells+surrounding interface cells have the higher importance */
1464   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1465        iter != mListEmpty.end(); iter++ ) {
1466     int i=iter->x, j=iter->y, k=iter->z;
1467                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)){ errMsg("A"," ARGHARGRAG "); } // DEBUG
1468                 if(debugFlagreinit) errMsg("EMPT", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) );
1469
1470                 /* set surrounding fluid cells to interface cells */
1471     FORDF1 {
1472                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1473       if( RFLAG(workLev,ni,nj,nk, workSet) & CFFluid){
1474                                 // init fluid->interface 
1475                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(CFInter); 
1476                                 changeFlag(workLev,ni,nj,nk, workSet, CFInter); 
1477                                 /* new mass = current density */
1478                                 LbmFloat nbrho = QCELL(workLev,ni,nj,nk, workSet, dC);
1479                 for(int rl=1; rl< this->cDfNum ; ++rl) { nbrho += QCELL(workLev,ni,nj,nk, workSet, rl); }
1480                                 QCELL(workLev,ni,nj,nk, workSet, dMass) =  nbrho; 
1481                                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) =  1.0; 
1482
1483                                 // store point
1484                                 addToNewInterList(ni,nj,nk);
1485       }
1486       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter){
1487                                 // test, also add to list...
1488                                 addToNewInterList(ni,nj,nk);
1489                         } // NEW?
1490     }
1491
1492                 /* for symmetry, set our flag right now */
1493                 changeFlag(workLev,i,j,k, workSet, CFEmpty);
1494                 // mark cell not be changed mass... - not necessary, not in list anymore anyway!
1495         } // emptylist
1496
1497
1498         
1499         // precompute weights to get rid of order dependancies
1500         vector<lbmFloatSet> vWeights;
1501         vWeights.resize( mListFull.size() + mListEmpty.size() );
1502         int weightIndex = 0;
1503   int nbCount = 0;
1504         LbmFloat nbWeights[LBM_DFNUM];
1505         LbmFloat nbTotWeights = 0.0;
1506         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1507        iter != mListFull.end(); iter++ ) {
1508     int i=iter->x, j=iter->y, k=iter->z;
1509     nbCount = 0; nbTotWeights = 0.0;
1510     FORDF1 {
1511                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1512       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1513                                 nbCount++;
1514                                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
1515                                 else nbWeights[l] = getMassdWeight(1,i,j,k,workSet,l); // NEWSURFT
1516                                 nbTotWeights += nbWeights[l];
1517       } else {
1518                                 nbWeights[l] = -100.0; // DEBUG;
1519                         }
1520     }
1521                 if(nbCount>0) { 
1522                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
1523         vWeights[weightIndex].val[0] = nbTotWeights;
1524         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
1525         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
1526                 } else { 
1527         vWeights[weightIndex].numNbs = 0.0;
1528                 }
1529                 weightIndex++;
1530         }
1531   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1532        iter != mListEmpty.end(); iter++ ) {
1533     int i=iter->x, j=iter->y, k=iter->z;
1534     nbCount = 0; nbTotWeights = 0.0;
1535     FORDF1 {
1536                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1537       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1538                                 nbCount++;
1539                                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
1540                                 else nbWeights[l] = getMassdWeight(0,i,j,k,workSet,l); // NEWSURFT
1541                                 nbTotWeights += nbWeights[l];
1542       } else {
1543                                 nbWeights[l] = -100.0; // DEBUG;
1544                         }
1545     }
1546                 if(nbCount>0) { 
1547                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
1548         vWeights[weightIndex].val[0] = nbTotWeights;
1549         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
1550         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
1551                 } else { 
1552         vWeights[weightIndex].numNbs = 0.0;
1553                 }
1554                 weightIndex++;
1555         } 
1556         weightIndex = 0;
1557         
1558
1559         /* process full list entries, filled cells are done after this loop */
1560         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1561        iter != mListFull.end(); iter++ ) {
1562     int i=iter->x, j=iter->y, k=iter->z;
1563
1564                 LbmFloat myrho = QCELL(workLev,i,j,k, workSet, dC);
1565     FORDF1 { myrho += QCELL(workLev,i,j,k, workSet, l); } // QCELL.rho
1566
1567     LbmFloat massChange = QCELL(workLev,i,j,k, workSet, dMass) - myrho;
1568                 if(vWeights[weightIndex].numNbs>0.0) {
1569                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
1570                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
1571                         FORDF1 {
1572                                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1573         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1574                                         LbmFloat change = -1.0;
1575                                         if(nbTotWeightsp>0.0) {
1576                                                 //change = massChange * ( nbWeights[l]/nbTotWeightsp );
1577                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
1578                                         } else {
1579                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
1580                                         }
1581                                         QCELL(workLev,ni,nj,nk, workSet, dMass) += change;
1582                                 }
1583                         }
1584                         massChange = 0.0;
1585                 } else {
1586                         // Problem! no interface neighbors
1587                         mFixMass += massChange;
1588                         //TTT mNumProblems++;
1589                         //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<<mFixMass);
1590                 }
1591                 weightIndex++;
1592
1593     // already done? RFLAG(workLev,i,j,k, workSet) = CFFluid;
1594     QCELL(workLev,i,j,k, workSet, dMass) = myrho; // should be rho... but unused?
1595     QCELL(workLev,i,j,k, workSet, dFfrac) = 1.0; // should be rho... but unused?
1596   } // fulllist
1597
1598
1599         /* now, finally handle the empty cells - order is important, has to be after
1600          * full cell handling */
1601   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1602        iter != mListEmpty.end(); iter++ ) {
1603     int i=iter->x, j=iter->y, k=iter->z;
1604
1605     LbmFloat massChange = QCELL(workLev, i,j,k, workSet, dMass);
1606                 if(vWeights[weightIndex].numNbs>0.0) {
1607                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
1608                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
1609                         FORDF1 {
1610                                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1611         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1612                                         LbmFloat change = -1.0;
1613                                         if(nbTotWeightsp>0.0) {
1614                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
1615                                         } else {
1616                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
1617                                         }
1618                                         QCELL(workLev, ni,nj,nk, workSet, dMass) += change;
1619                                 }
1620                         }
1621                         massChange = 0.0;
1622                 } else {
1623                         // Problem! no interface neighbors
1624                         mFixMass += massChange;
1625                         //TTT mNumProblems++;
1626                         //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<mFixMass);
1627                 }
1628                 weightIndex++;
1629                 
1630                 // finally... make it empty 
1631     // already done? RFLAG(workLev,i,j,k, workSet) = CFEmpty;
1632     QCELL(workLev,i,j,k, workSet, dMass) = 0.0;
1633     QCELL(workLev,i,j,k, workSet, dFfrac) = 0.0;
1634         }
1635   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1636        iter != mListEmpty.end(); iter++ ) {
1637     int i=iter->x, j=iter->y, k=iter->z;
1638     changeFlag(workLev,i,j,k, otherSet, CFEmpty);
1639         } 
1640
1641
1642         // check if some of the new interface cells can be removed again 
1643         // never happens !!! not necessary
1644         // calculate ffrac for new IF cells NEW
1645
1646         // how many are really new interface cells?
1647         int numNewIf = 0;
1648   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1649        iter != mListNewInter.end(); iter++ ) {
1650     int i=iter->x, j=iter->y, k=iter->z;
1651                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
1652                         continue; 
1653                         // FIXME remove from list?
1654                 }
1655                 numNewIf++;
1656         }
1657
1658         // redistribute mass, reinit flags
1659         if(debugFlagreinit) errMsg("NEWIF", "total:"<<mListNewInter.size());
1660         float newIfFac = 1.0/(LbmFloat)numNewIf;
1661   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1662        iter != mListNewInter.end(); iter++ ) {
1663     int i=iter->x, j=iter->y, k=iter->z;
1664                 if((i<=0) || (j<=0) || 
1665                          (i>=mLevel[workLev].lSizex-1) ||
1666                          (j>=mLevel[workLev].lSizey-1) ||
1667                          ((LBMDIM==3) && ((k<=0) || (k>=mLevel[workLev].lSizez-1) ) )
1668                          ) {
1669                         continue; } // new bc, dont treat cells on boundary NEWBC
1670                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
1671                         //errMsg("???"," "<<PRINT_IJK);
1672                         continue; 
1673                 } // */
1674
1675     QCELL(workLev,i,j,k, workSet, dMass) += (mFixMass * newIfFac);
1676                 
1677                 int nbored = 0;
1678                 FORDF1 { nbored |= RFLAG_NB(workLev, i,j,k, workSet,l); }
1679                 if(!(nbored & CFBndNoslip)) { RFLAG(workLev,i,j,k, workSet) |= CFNoBndFluid; }
1680                 if(!(nbored & CFFluid))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbFluid; }
1681                 if(!(nbored & CFEmpty))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbEmpty; }
1682
1683                 if(!(RFLAG(workLev,i,j,k, otherSet)&CFInter)) {
1684                         RFLAG(workLev,i,j,k, workSet) = (CellFlagType)(RFLAG(workLev,i,j,k, workSet) | CFNoDelete);
1685                 }
1686                 if(debugFlagreinit) errMsg("NEWIF", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" f"<< convertCellFlagType2String(RFLAG(workLev,i,j,k, workSet))<<" wl"<<workLev );
1687         }
1688
1689         // reinit fill fraction
1690   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1691        iter != mListNewInter.end(); iter++ ) {
1692     int i=iter->x, j=iter->y, k=iter->z;
1693                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { continue; }
1694
1695                 initInterfaceVars(workLev, i,j,k, workSet, false); //int level, int i,int j,int k,int workSet, bool initMass) {
1696                 //LbmFloat nrho = 0.0;
1697                 //FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); }
1698     //QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho;
1699     //QCELL(workLev,i,j,k, workSet, dFlux) = FLUX_INIT;
1700         }
1701
1702         if(mListNewInter.size()>0){ 
1703                 //errMsg("FixMassDisted"," fm:"<<mFixMass<<" nif:"<<mListNewInter.size() );
1704                 mFixMass = 0.0; 
1705         }
1706
1707         // empty lists for next step
1708         mListFull.clear();
1709         mListEmpty.clear();
1710         mListNewInter.clear();
1711 } // reinitFlags
1712
1713
1714