svn merge -r 12716:12856 https://svn.blender.org/svnroot/bf-blender/trunk/blender
[blender.git] / intern / elbeem / intern / solver_main.cpp
1 /******************************************************************************
2  *
3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4  * Copyright 2003-2006 Nils Thuerey
5  *
6  * Standard LBM Factory implementation
7  *
8  *****************************************************************************/
9
10 #include "solver_class.h"
11 #include "solver_relax.h"
12 #include "particletracer.h"
13 #include "loop_tools.h"
14 #include <stdlib.h>
15
16 #if !defined(linux) && (defined (__sparc) || defined (__sparc__))
17 #include <ieeefp.h>
18 #endif
19
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 extern int glob_mpindex,glob_mpnum;
30 extern bool globOutstrForce;
31
32 // simulation object interface
33 void LbmFsgrSolver::step() { 
34         stepMain();
35 }
36
37 // lbm main step
38 void messageOutputForce(string from);
39 void LbmFsgrSolver::stepMain() { 
40         myTime_t timestart = getTime();
41
42         initLevelOmegas();
43         markedClearList(); // DMC clearMarkedCellsList
44
45         // safety check, counter reset
46         mNumUsedCells = 0;
47         mNumInterdCells = 0;
48         mNumInvIfCells = 0;
49
50   //debugOutNnl("LbmFsgrSolver::step : "<<mStepCnt, 10);
51   if(!mSilent){ debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime, 10); }
52         //debMsgDirect(  "LbmFsgrSolver::step : "<<mStepCnt<<" ");
53         //myTime_t timestart = 0;
54         //if(mStartSymm) { checkSymmetry("step1"); } // DEBUG 
55
56         // time adapt
57         mMaxVlen = mMxvz = mMxvy = mMxvx = 0.0;
58
59         // init moving bc's, can change mMaxVlen
60         initMovingObstacles(false);
61 #if LBM_INCLUDE_TESTSOLVERS==1
62         handleCpdata();
63 #endif
64
65         // important - keep for tadap
66         LbmFloat lastMass = mCurrentMass;
67         mCurrentMass = mFixMass; // reset here for next step
68         mCurrentVolume = 0.0;
69         
70         //change to single step advance!
71         int levsteps = 0;
72         int dsbits = mStepCnt ^ (mStepCnt-1);
73         //errMsg("S"," step:"<<mStepCnt<<" s-1:"<<(mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits));
74         for(int lev=0; lev<=mMaxRefine; lev++) {
75                 //if(! (mStepCnt&(1<<lev)) ) {
76                 if( dsbits & (1<<(mMaxRefine-lev)) ) {
77                         //errMsg("S"," l"<<lev);
78
79                         if(lev==mMaxRefine) {
80                                 // always advance fine level...
81                                 fineAdvance();
82                         } else {
83                                 adaptGrid(lev);
84                                 coarseRestrictFromFine(lev);
85                                 coarseAdvance(lev);
86                         }
87 #if FSGR_OMEGA_DEBUG==1
88                         errMsg("LbmFsgrSolver::step","LES stats l="<<lev<<" omega="<<mLevel[lev].omega<<" avgOmega="<< (mLevel[lev].avgOmega/mLevel[lev].avgOmegaCnt) );
89                         mLevel[lev].avgOmega = 0.0; mLevel[lev].avgOmegaCnt = 0.0;
90 #endif // FSGR_OMEGA_DEBUG==1
91                         levsteps++;
92                 }
93                 mCurrentMass   += mLevel[lev].lmass;
94                 mCurrentVolume += mLevel[lev].lvolume;
95         }
96
97   // prepare next step
98         mStepCnt++;
99
100
101         // some dbugging output follows
102         // calculate MLSUPS
103         myTime_t timeend = getTime();
104
105         mNumUsedCells += mNumInterdCells; // count both types for MLSUPS
106         mAvgNumUsedCells += mNumUsedCells;
107         mMLSUPS = ((double)mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000.0);
108         if(mMLSUPS>10000){ mMLSUPS = -1; }
109         //else { mAvgMLSUPS += mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
110         
111         LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1())) / ((timeend-timestart)/(double)1000.0) ) / (1000000);
112         if(totMLSUPS>10000) totMLSUPS = -1;
113         mNumInvIfTotal += mNumInvIfCells; // debug
114
115   // do some formatting 
116   if(!mSilent){ 
117                 int avgcls = (int)(mAvgNumUsedCells/(LONGINT)mStepCnt);
118         debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime<<
119                         " cur-mlsups:"<<mMLSUPS<< //" avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< 
120                         " totcls:"<<mNumUsedCells<< " avgcls:"<< avgcls<< 
121                         " intd:"<<mNumInterdCells<< " invif:"<<mNumInvIfCells<< 
122                         " invift:"<<mNumInvIfTotal<< " fsgrcs:"<<mNumFsgrChanges<< 
123                         " filled:"<<mNumFilledCells<<", emptied:"<<mNumEmptiedCells<< 
124                         " mMxv:"<<PRINT_VEC(mMxvx,mMxvy,mMxvz)<<", tscnts:"<<mTimeSwitchCounts<< 
125                         //" RWmxv:"<<ntlVec3Gfx(mMxvx,mMxvy,mMxvz)*(mLevel[mMaxRefine].simCellSize / mLevel[mMaxRefine].timestep)<<" "<< /* realworld vel output */
126                         " probs:"<<mNumProblems<< " simt:"<<mSimulationTime<< 
127                         " took:"<< getTimeString(timeend-timestart)<<
128                         " for '"<<mName<<"' " , 10);
129         } else { debMsgDirect("."); }
130
131         if(mStepCnt==1) {
132                 mMinNoCells = mMaxNoCells = mNumUsedCells;
133         } else {
134                 if(mNumUsedCells>mMaxNoCells) mMaxNoCells = mNumUsedCells;
135                 if(mNumUsedCells<mMinNoCells) mMinNoCells = mNumUsedCells;
136         }
137         
138         // mass scale test
139         if((mMaxRefine>0)&&(mInitialMass>0.0)) {
140                 LbmFloat mscale = mInitialMass/mCurrentMass;
141
142                 mscale = 1.0;
143                 const LbmFloat dchh = 0.001;
144                 if(mCurrentMass<mInitialMass) mscale = 1.0+dchh;
145                 if(mCurrentMass>mInitialMass) mscale = 1.0-dchh;
146
147                 // use mass rescaling?
148                 // with float precision this seems to be nonsense...
149                 const bool MREnable = false;
150
151                 const int MSInter = 2;
152                 static int mscount = 0;
153                 if( (MREnable) && ((mLevel[0].lsteps%MSInter)== (MSInter-1)) && ( ABS( (mInitialMass/mCurrentMass)-1.0 ) > 0.01) && ( dsbits & (1<<(mMaxRefine-0)) ) ){
154                         // example: FORCE RESCALE MASS! ini:1843.5, cur:1817.6, f=1.01425 step:22153 levstep:5539 msc:37
155                         // mass rescale MASS RESCALE check
156                         errMsg("MDTDD","\n\n");
157                         errMsg("MDTDD","FORCE RESCALE MASS! "
158                                         <<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass)
159                                         <<" step:"<<mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
160                                         );
161                         errMsg("MDTDD","\n\n");
162
163                         mscount++;
164                         for(int lev=mMaxRefine; lev>=0 ; lev--) {
165                                 //for(int workSet = 0; workSet<=1; workSet++) {
166                                 int wss = 0;
167                                 int wse = 1;
168 #if COMPRESSGRIDS==1
169                                 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
170 #endif // COMPRESSGRIDS==1
171                                 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
172
173                                         FSGR_FORIJK1(lev) {
174                                                 if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) 
175                                                         ) {
176
177                                                         FORDF0 { QCELL(lev, i,j,k,workSet, l) *= mscale; }
178                                                         QCELL(lev, i,j,k,workSet, dMass) *= mscale;
179                                                         QCELL(lev, i,j,k,workSet, dFfrac) *= mscale;
180
181                                                 } else {
182                                                         continue;
183                                                 }
184                                         }
185                                 }
186                                 mLevel[lev].lmass *= mscale;
187                         }
188                 } 
189
190                 mCurrentMass *= mscale;
191         }// if mass scale test */
192         else {
193                 // use current mass after full step for initial setting
194                 if((mMaxRefine>0)&&(mInitialMass<=0.0) && (levsteps == (mMaxRefine+1))) {
195                         mInitialMass = mCurrentMass;
196                         debMsgStd("MDTDD",DM_NOTIFY,"Second Initial Mass Init: "<<mInitialMass, 2);
197                 }
198         }
199
200 #if LBM_INCLUDE_TESTSOLVERS==1
201         if((mUseTestdata)&&(mInitDone)) { handleTestdata(); }
202         mrExchange();
203 #endif
204
205         // advance positions with current grid
206         advanceParticles();
207         if(mpParticles) {
208                 mpParticles->checkDumpTextPositions(mSimulationTime);
209                 mpParticles->checkTrails(mSimulationTime);
210         }
211
212         // one of the last things to do - adapt timestep
213         // was in fineAdvance before... 
214         if(mTimeAdap) {
215                 adaptTimestep();
216         } // time adaptivity
217
218
219 #ifndef WIN32
220         // good indicator for instabilities
221         if( (!finite(mMxvx)) || (!finite(mMxvy)) || (!finite(mMxvz)) ) { CAUSE_PANIC; }
222         if( (!finite(mCurrentMass)) || (!finite(mCurrentVolume)) ) { CAUSE_PANIC; }
223 #endif // WIN32
224
225         // output total step time
226         myTime_t timeend2 = getTime();
227         double mdelta = (lastMass-mCurrentMass);
228         if(ABS(mdelta)<1e-12) mdelta=0.;
229         double effMLSUPS = ((double)mNumUsedCells / ((timeend2-timestart)/(double)1000.0) ) / (1000000.0);
230         if(mInitDone) {
231                 if(effMLSUPS>10000){ effMLSUPS = -1; }
232                 else { mAvgMLSUPS += effMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
233         }
234         
235         debMsgStd("LbmFsgrSolver::stepMain", DM_MSG, "mmpid:"<<glob_mpindex<<" step:"<<mStepCnt
236                         <<" dccd="<< mCurrentMass
237                         //<<",d"<<mdelta
238                         //<<",ds"<<(mCurrentMass-mObjectMassMovnd[1])
239                         //<<"/"<<mCurrentVolume<<"(fix="<<mFixMass<<",ini="<<mInitialMass<<"), "
240                         <<" effMLSUPS=("<< effMLSUPS
241                         <<",avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "
242                         <<" took totst:"<< getTimeString(timeend2-timestart), 3);
243         // nicer output
244         //debMsgDirect(std::endl); 
245         // */
246
247         messageOutputForce("");
248  //#endif // ELBEEM_PLUGIN!=1
249 }
250
251 #define NEWDEBCHECK(str) \
252         if(!this->mPanic){ FSGR_FORIJK_BOUNDS(mMaxRefine) { \
253                 if(RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr)&(CFFluid|CFInter)) { \
254                 for(int l=0;l<dTotalNum;l++) { \
255                         if(!finite(QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr,l))) { errMsg("NNOFIN"," "<<str<<" at "<<PRINT_IJK<<" l"<<l<<" "); }\
256                 }/*for*/ \
257                 }/*if*/ \
258         } }
259
260 void LbmFsgrSolver::fineAdvance()
261 {
262         // do the real thing...
263         //NEWDEBCHECK("t1");
264         mainLoop( mMaxRefine );
265         if(mUpdateFVHeight) {
266                 // warning assume -Y gravity...
267                 mFVHeight = mCurrentMass*mFVArea/((LbmFloat)(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizez));
268                 if(mFVHeight<1.0) mFVHeight = 1.0;
269                 mpParam->setFluidVolumeHeight(mFVHeight);
270         }
271
272         // advance time before timestep change
273         mSimulationTime += mpParam->getTimestep();
274         // time adaptivity
275         mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) );
276         //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG 
277         if(!mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps), 3 ); }
278
279         // update other set
280   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
281   mLevel[mMaxRefine].setCurr   ^= 1;
282   mLevel[mMaxRefine].lsteps++;
283
284         // flag init... (work on current set, to simplify flag checks)
285         reinitFlags( mLevel[mMaxRefine].setCurr );
286         if(!mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," flags reinit on set "<< mLevel[mMaxRefine].setCurr, 3 ); }
287
288         // DEBUG VEL CHECK
289         if(0) {
290                 int lev = mMaxRefine;
291                 int workSet = mLevel[lev].setCurr;
292                 int mi=0,mj=0,mk=0;
293                 LbmFloat compRho=0.;
294                 LbmFloat compUx=0.;
295                 LbmFloat compUy=0.;
296                 LbmFloat compUz=0.;
297                 LbmFloat maxUlen=0.;
298                 LbmVec maxU(0.);
299                 LbmFloat maxRho=-100.;
300                 int ri=0,rj=0,rk=0;
301
302                 FSGR_FORIJK1(lev) {
303                         if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) ) {
304                                 compRho=QCELL(lev, i,j,k,workSet, dC);
305                                 compUx = compUy = compUz = 0.0;
306                                 for(int l=1; l<this->cDfNum; l++) {
307                                         LbmFloat df = QCELL(lev, i,j,k,workSet, l);
308                                         compRho += df;
309                                         compUx  += (this->dfDvecX[l]*df);
310                                         compUy  += (this->dfDvecY[l]*df); 
311                                         compUz  += (this->dfDvecZ[l]*df); 
312                                 } 
313                                 LbmVec u(compUx,compUy,compUz);
314                                 LbmFloat nu = norm(u);
315                                 if(nu>maxUlen) {
316                                         maxU = u;
317                                         maxUlen = nu;
318                                         mi=i; mj=j; mk=k;
319                                 }
320                                 if(compRho>maxRho) {
321                                         maxRho=compRho;
322                                         ri=i; rj=j; rk=k;
323                                 }
324                         } else {
325                                 continue;
326                         }
327                 }
328
329                 errMsg("MAXVELCHECK"," at "<<PRINT_VEC(mi,mj,mk)<<" norm:"<<maxUlen<<" u:"<<maxU);
330                 errMsg("MAXRHOCHECK"," at "<<PRINT_VEC(ri,rj,rk)<<" rho:"<<maxRho);
331                 printLbmCell(lev, 30,36,23, -1);
332         } // DEBUG VEL CHECK
333
334 }
335
336
337
338 // fine step defines
339
340 // access to own dfs during step (may be changed to local array)
341 #define MYDF(l) RAC(ccel, l)
342
343 // drop model definitions
344 #define RWVEL_THRESH 1.5
345 #define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5)
346
347 #if LBMDIM==3
348 // normal
349 #define SLOWDOWNREGION (mSizez/4)
350 #else // LBMDIM==2
351 // off
352 #define SLOWDOWNREGION 10 
353 #endif // LBMDIM==2
354 #define P_LCSMQO 0.01
355
356 /*****************************************************************************/
357 //! fine step function
358 /*****************************************************************************/
359 void 
360 LbmFsgrSolver::mainLoop(int lev)
361 {
362         // loops over _only inner_ cells  -----------------------------------------------------------------------------------
363         
364         // slow surf regions smooth (if below)
365         const LbmFloat smoothStrength = 0.0; //0.01;
366         const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03;
367         const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit;
368
369         const int cutMin  = 1;
370         const int cutConst = mCutoff+2;
371         
372 #       if LBM_INCLUDE_TESTSOLVERS==1
373         // 3d region off... quit
374         if((mUseTestdata)&&(mpTest->mFarfMode>0)) { return; }
375 #       endif // ELBEEM_PLUGIN!=1
376         
377   // main loop region
378         const bool doReduce = true;
379         const int gridLoopBound=1;
380         GRID_REGION_INIT();
381 #if PARALLEL==1
382 #pragma omp parallel default(shared) \
383   reduction(+: \
384           calcCurrentMass,calcCurrentVolume, \
385                 calcCellsFilled,calcCellsEmptied, \
386                 calcNumUsedCells )
387         GRID_REGION_START();
388 #else // PARALLEL==1
389         GRID_REGION_START();
390 #endif // PARALLEL==1
391
392         // local to main
393         CellFlagType nbflag[LBM_DFNUM]; 
394         int oldFlag, newFlag, nbored;
395         LbmFloat m[LBM_DFNUM];
396         LbmFloat rho, ux, uy, uz, tmp, usqr;
397
398         // smago vars
399         LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
400         
401         // ifempty cell conversion flags
402         bool iffilled, ifemptied;
403         LbmFloat nbfracs[LBM_DFNUM]; // ffracs of neighbors
404         int recons[LBM_DFNUM];   // reconstruct this DF?
405         int numRecons;           // how many are reconstructed?
406
407         LbmFloat mass=0., change=0., lcsmqo=0.;
408         rho= ux= uy= uz= usqr= tmp= 0.; 
409         lcsmqadd = lcsmomega = 0.;
410         FORDF0{ lcsmeq[l] = 0.; }
411
412         // ---
413         // now stream etc.
414         // use //template functions for 2D/3D
415
416         GRID_LOOP_START();
417                 // loop start
418                 // stream from current set to other, then collide and store
419                 //errMsg("l2"," at "<<PRINT_IJK<<" id"<<id);
420
421 #               if FSGR_STRICT_DEBUG==1
422                 // safety pointer check
423                 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
424                 if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) || 
425                     (&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) {
426                         errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<<
427                                         RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<<
428                                         (*pFlagSrc)<<","<<(*pFlagDst) <<",  pointers "<<
429           (long)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(long)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
430           (long)(pFlagSrc)<<","<<(long)(pFlagDst)<<" "
431                                         ); 
432                         CAUSE_PANIC;
433                 }       
434                 if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) || 
435                     (&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) {
436                         errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<<
437           (long)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(long)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
438           (long)(ccel)<<","<<(long)(tcel)<<" "
439                                         ); 
440                         CAUSE_PANIC;
441                 }       
442 #               endif
443                 oldFlag = *pFlagSrc;
444                 
445                 // old INTCFCOARSETEST==1
446                 if( (oldFlag & (CFGrFromCoarse)) ) { 
447                         if(( mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
448                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
449                         } else {
450                                 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
451                                 calcNumUsedCells++;
452                         }
453                         continue; // interpolateFineFromCoarse test!
454                 } // interpolateFineFromCoarse test! 
455         
456                 if(oldFlag & (CFMbndInflow)) {
457                         // fluid & if are ok, fill if later on
458                         int isValid = oldFlag & (CFFluid|CFInter);
459                         const LbmFloat iniRho = 1.0;
460                         const int OId = oldFlag>>24;
461                         if(!isValid) {
462                                 // make new if cell
463                                 const LbmVec vel(mObjectSpeeds[OId]);
464                                 // TODO add OPT3D treatment
465                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
466                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
467                                 RAC(tcel, dFlux) = FLUX_INIT;
468                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
469                                 calcCurrentMass += iniRho; 
470                                 calcCurrentVolume += 1.0; 
471                                 calcNumUsedCells++;
472                                 mInitialMass += iniRho;
473                                 // dont treat cell until next step
474                                 continue;
475                         } 
476                 } 
477                 else  // these are exclusive
478                 if(oldFlag & (CFMbndOutflow)) {
479                         int isnotValid = oldFlag & (CFFluid);
480                         if(isnotValid) {
481                                 // remove fluid cells, shouldnt be here anyway
482                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
483                                 mInitialMass -= fluidRho;
484                                 const LbmFloat iniRho = 0.0;
485                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
486                                 RAC(tcel, dFlux) = FLUX_INIT;
487                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
488
489                                 // same as ifemptied for if below
490                                 LbmPoint oemptyp; oemptyp.flag = 0;
491                                 oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
492                                 LIST_EMPTY(oemptyp);
493                                 calcCellsEmptied++;
494                                 continue;
495                         }
496                 }
497
498                 if(oldFlag & (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)) { 
499                         *pFlagDst = oldFlag;
500                         continue;
501                 }
502                 /*if( oldFlag & CFNoBndFluid ) {  // TEST ME FASTER?
503                         OPTIMIZED_STREAMCOLLIDE; PERFORM_USQRMAXCHECK;
504                         RAC(tcel,dFfrac) = 1.0; 
505                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
506                         calcCurrentMass += rho; calcCurrentVolume += 1.0;
507                         calcNumUsedCells++;
508                         continue;
509                 }// TEST ME FASTER? */
510
511                 // only neighbor flags! not own flag
512                 nbored = 0;
513                 
514 #if OPT3D==0
515                 FORDF1 {
516                         nbflag[l] = RFLAG_NB(lev, i,j,k,SRCS(lev),l);
517                         nbored |= nbflag[l];
518                 } 
519 #else
520                 nbflag[dSB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dSB];
521                 nbflag[dWB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWB];
522                 nbflag[ dB] = *(pFlagSrc + (-mLevel[lev].lOffsy)); nbored |= nbflag[dB];
523                 nbflag[dEB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dEB];
524                 nbflag[dNB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNB];
525
526                 nbflag[dSW] = *(pFlagSrc + (-mLevel[lev].lOffsx+-1)); nbored |= nbflag[dSW];
527                 nbflag[ dS] = *(pFlagSrc + (-mLevel[lev].lOffsx)); nbored |= nbflag[dS];
528                 nbflag[dSE] = *(pFlagSrc + (-mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dSE];
529
530                 nbflag[ dW] = *(pFlagSrc + (-1)); nbored |= nbflag[dW];
531                 nbflag[ dE] = *(pFlagSrc + ( 1)); nbored |= nbflag[dE];
532
533                 nbflag[dNW] = *(pFlagSrc + ( mLevel[lev].lOffsx+-1)); nbored |= nbflag[dNW];
534           nbflag[ dN] = *(pFlagSrc + ( mLevel[lev].lOffsx)); nbored |= nbflag[dN];
535                 nbflag[dNE] = *(pFlagSrc + ( mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dNE];
536
537                 nbflag[dST] = *(pFlagSrc + ( mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dST];
538                 nbflag[dWT] = *(pFlagSrc + ( mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWT];
539                 nbflag[ dT] = *(pFlagSrc + ( mLevel[lev].lOffsy)); nbored |= nbflag[dT];
540                 nbflag[dET] = *(pFlagSrc + ( mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dET];
541                 nbflag[dNT] = *(pFlagSrc + ( mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNT];
542                 // */
543 #endif
544
545                 // pointer to destination cell
546                 calcNumUsedCells++;
547
548                 // FLUID cells 
549                 if( oldFlag & CFFluid ) { 
550                         // only standard fluid cells (with nothing except fluid as nbs
551
552                         if(oldFlag&CFMbndInflow) {
553                                 // force velocity for inflow, necessary to have constant direction of flow
554                                 // FIXME , test also set interface cells?
555                                 const int OId = oldFlag>>24;
556                                 //? DEFAULT_STREAM;
557                                 //const LbmFloat fluidRho = 1.0;
558                                 // for submerged inflows, streaming would have to be performed...
559                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
560                                 const LbmVec vel(mObjectSpeeds[OId]);
561                                 ux=vel[0], uy=vel[1], uz=vel[2]; 
562                                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
563                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); }
564                                 rho = fluidRho;
565                                 //errMsg("INFLOW_DEBUG","std at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
566                         } else {
567                                 if(nbored&CFBnd) {
568                                         DEFAULT_STREAM;
569                                         //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
570                                         DEFAULT_COLLIDEG(mLevel[lev].gravity);
571                                         oldFlag &= (~CFNoBndFluid);
572                                 } else {
573                                         // do standard stream/collide
574                                         OPTIMIZED_STREAMCOLLIDE;
575                                         oldFlag |= CFNoBndFluid;
576                                 } 
577                         }
578
579                         PERFORM_USQRMAXCHECK;
580                         // "normal" fluid cells
581                         RAC(tcel,dFfrac) = 1.0; 
582                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
583                         calcCurrentMass += rho; 
584                         calcCurrentVolume += 1.0;
585                         continue;
586                 }
587                 
588                 newFlag  = oldFlag;
589                 // make sure here: always check which flags to really unset...!
590                 newFlag = newFlag & (~( 
591                                         CFNoNbFluid|CFNoNbEmpty| CFNoDelete 
592                                         | CFNoInterpolSrc
593                                         | CFNoBndFluid
594                                         ));
595                 if(!(nbored&CFBndNoslip)) { //NEWSURFT NEWSURFTNOS
596                         newFlag |= CFNoBndFluid;
597                 }
598                 /*if(!(nbored&CFBnd)) { //NEWSURFT NEWSURFTNOS
599                         // explicitly check for noslip neighbors
600                         bool hasnoslipnb = false;
601                         FORDF1 { if((nbflag[l]&CFBnd)&&(nbflag[l]&CFBndNoslip)) hasnoslipnb=true; }
602                         if(!hasnoslipnb) newFlag |= CFNoBndFluid; 
603                 } // */
604
605                 // store own dfs and mass
606                 mass = RAC(ccel,dMass);
607
608                 // WARNING - only interface cells arrive here!
609                 // read distribution funtions of adjacent cells = stream step
610                 DEFAULT_STREAM;
611
612                 if((nbored & CFFluid)==0) { newFlag |= CFNoNbFluid; mNumInvIfCells++; }
613                 if((nbored & CFEmpty)==0) { newFlag |= CFNoNbEmpty; mNumInvIfCells++; }
614
615                 // calculate mass exchange for interface cells 
616                 LbmFloat myfrac = RAC(ccel,dFfrac);
617                 if(myfrac<0.) myfrac=0.; // NEWSURFT
618 #               define nbdf(l) m[ this->dfInv[(l)] ]
619
620                 // update mass 
621                 // only do boundaries for fluid cells, and interface cells without
622                 // any fluid neighbors (assume that interface cells _with_ fluid
623                 // neighbors are affected enough by these) 
624                 // which Df's have to be reconstructed? 
625                 // for fluid cells - just the f_i difference from streaming to empty cells  ----
626                 numRecons = 0;
627                 bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
628                 //onlyBndnb = false; // DEBUG test off
629
630                 FORDF1 { // dfl loop
631                         recons[l] = 0;
632                         nbfracs[l] = 0.0;
633                         // finally, "normal" interface cells ----
634                         if( nbflag[l]&(CFFluid|CFBnd) ) { // NEWTEST! FIXME check!!!!!!!!!!!!!!!!!!
635                                 change = nbdf(l) - MYDF(l);
636                         }
637                         // interface cells - distuingish cells that shouldn't fill/empty 
638                         else if( nbflag[l] & CFInter ) {
639                                 
640                                 LbmFloat mynbfac,nbnbfac;
641                                 // NEW TEST t1
642                                 // t2 -> off
643                                 if((oldFlag&CFNoBndFluid)&&(nbflag[l]&CFNoBndFluid)) {
644                                         mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
645                                         nbnbfac = 1.0/mynbfac;
646                                         onlyBndnb = false;
647                                 } else {
648                                         mynbfac = nbnbfac = 1.0; // switch calc flux off
649                                         goto changeDefault;  // NEWSURFT
650                                         //change = 0.; goto changeDone;  // NEWSURFT
651                                 }
652                                 //mynbfac = nbnbfac = 1.0; // switch calc flux off t3
653
654                                 // perform interface case handling
655                                 if ((oldFlag|nbflag[l])&(CFNoNbFluid|CFNoNbEmpty)) {
656                                 switch (oldFlag&(CFNoNbFluid|CFNoNbEmpty)) {
657                                         case 0: 
658                                                 // we are a normal cell so... 
659                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
660                                                         case CFNoNbFluid: 
661                                                                 // just fill current cell = empty neighbor 
662                                                                 change = nbnbfac*nbdf(l) ; goto changeDone; 
663                                                         case CFNoNbEmpty: 
664                                                                 // just empty current cell = fill neighbor 
665                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
666                                                 }
667                                                 break;
668
669                                         case CFNoNbFluid: 
670                                                 // we dont have fluid nb's so...
671                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
672                                                         case 0: 
673                                                         case CFNoNbEmpty: 
674                                                                 // we have no fluid nb's -> just empty
675                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
676                                                 }
677                                                 break;
678
679                                         case CFNoNbEmpty: 
680                                                 // we dont have empty nb's so...
681                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
682                                                         case 0: 
683                                                         case CFNoNbFluid: 
684                                                                 // we have no empty nb's -> just fill
685                                                                 change = nbnbfac*nbdf(l); goto changeDone; 
686                                                 }
687                                                 break;
688                                 }} // inter-inter exchange
689
690                         changeDefault: ;
691                                 // just do normal mass exchange...
692                                 change = ( nbnbfac*nbdf(l) - mynbfac*MYDF(l) ) ;
693                         changeDone: ;
694                                 nbfracs[l] = QCELL_NB(lev, i,j,k, SRCS(lev),l, dFfrac);
695                                 if(nbfracs[l]<0.) nbfracs[l] = 0.; // NEWSURFT
696                                 change *=  (myfrac + nbfracs[l]) * 0.5;
697                         } // the other cell is interface
698
699                         // last alternative - reconstruction in this direction
700                         else {
701                                 // empty + bnd case
702                                 recons[l] = 1; 
703                                 numRecons++;
704                                 change = 0.0; 
705                         }
706
707                         // modify mass at SRCS
708                         mass += change;
709                 } // l
710                 // normal interface, no if empty/fluid
711
712                 // computenormal
713                 LbmFloat surfaceNormal[3];
714                 computeFluidSurfaceNormal(ccel,pFlagSrc, surfaceNormal);
715
716                 if( (ABS(surfaceNormal[0])+ABS(surfaceNormal[1])+ABS(surfaceNormal[2])) > LBM_EPSILON) {
717                         // normal ok and usable...
718                         FORDF1 {
719                                 if( (this->dfDvecX[l]*surfaceNormal[0] + this->dfDvecY[l]*surfaceNormal[1] + this->dfDvecZ[l]*surfaceNormal[2])  // dot Dvec,norml
720                                                 > LBM_EPSILON) {
721                                         recons[l] = 2; 
722                                         numRecons++;
723                                 } 
724                         }
725                 }
726
727                 // calculate macroscopic cell values
728                 LbmFloat oldUx, oldUy, oldUz;
729                 LbmFloat oldRho; // OLD rho = ccel->rho;
730 #               define REFERENCE_PRESSURE 1.0 // always atmosphere...
731 #               if OPT3D==0
732                 oldRho=RAC(ccel,0);
733                 oldUx = oldUy = oldUz = 0.0;
734                 for(int l=1; l<this->cDfNum; l++) {
735                         oldRho += RAC(ccel,l);
736                         oldUx  += (this->dfDvecX[l]*RAC(ccel,l));
737                         oldUy  += (this->dfDvecY[l]*RAC(ccel,l)); 
738                         oldUz  += (this->dfDvecZ[l]*RAC(ccel,l)); 
739                 } 
740                 // reconstruct dist funcs from empty cells
741                 FORDF1 {
742                         if(recons[ l ]) {
743                                 m[ this->dfInv[l] ] = 
744                                         this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 
745                                         this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 
746                                         - MYDF( l );
747                         }
748                 }
749                 ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
750                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
751 #               else // OPT3D==0
752                 oldRho = + RAC(ccel,dC)  + RAC(ccel,dN )
753                                 + RAC(ccel,dS ) + RAC(ccel,dE )
754                                 + RAC(ccel,dW ) + RAC(ccel,dT )
755                                 + RAC(ccel,dB ) + RAC(ccel,dNE)
756                                 + RAC(ccel,dNW) + RAC(ccel,dSE)
757                                 + RAC(ccel,dSW) + RAC(ccel,dNT)
758                                 + RAC(ccel,dNB) + RAC(ccel,dST)
759                                 + RAC(ccel,dSB) + RAC(ccel,dET)
760                                 + RAC(ccel,dEB) + RAC(ccel,dWT)
761                                 + RAC(ccel,dWB);
762
763                 oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
764                                 + RAC(ccel,dNE) - RAC(ccel,dNW)
765                                 + RAC(ccel,dSE) - RAC(ccel,dSW)
766                                 + RAC(ccel,dET) + RAC(ccel,dEB)
767                                 - RAC(ccel,dWT) - RAC(ccel,dWB);
768
769                 oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
770                                 + RAC(ccel,dNE) + RAC(ccel,dNW)
771                                 - RAC(ccel,dSE) - RAC(ccel,dSW)
772                                 + RAC(ccel,dNT) + RAC(ccel,dNB)
773                                 - RAC(ccel,dST) - RAC(ccel,dSB);
774
775                 oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
776                                 + RAC(ccel,dNT) - RAC(ccel,dNB)
777                                 + RAC(ccel,dST) - RAC(ccel,dSB)
778                                 + RAC(ccel,dET) - RAC(ccel,dEB)
779                                 + RAC(ccel,dWT) - RAC(ccel,dWB);
780
781                 // now reconstruction
782                 ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
783                 rho = REFERENCE_PRESSURE;
784                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
785                 if(recons[dN ]) { m[dS ] = EQN  + EQS  - MYDF(dN ); }
786                 if(recons[dS ]) { m[dN ] = EQS  + EQN  - MYDF(dS ); }
787                 if(recons[dE ]) { m[dW ] = EQE  + EQW  - MYDF(dE ); }
788                 if(recons[dW ]) { m[dE ] = EQW  + EQE  - MYDF(dW ); }
789                 if(recons[dT ]) { m[dB ] = EQT  + EQB  - MYDF(dT ); }
790                 if(recons[dB ]) { m[dT ] = EQB  + EQT  - MYDF(dB ); }
791                 if(recons[dNE]) { m[dSW] = EQNE + EQSW - MYDF(dNE); }
792                 if(recons[dNW]) { m[dSE] = EQNW + EQSE - MYDF(dNW); }
793                 if(recons[dSE]) { m[dNW] = EQSE + EQNW - MYDF(dSE); }
794                 if(recons[dSW]) { m[dNE] = EQSW + EQNE - MYDF(dSW); }
795                 if(recons[dNT]) { m[dSB] = EQNT + EQSB - MYDF(dNT); }
796                 if(recons[dNB]) { m[dST] = EQNB + EQST - MYDF(dNB); }
797                 if(recons[dST]) { m[dNB] = EQST + EQNB - MYDF(dST); }
798                 if(recons[dSB]) { m[dNT] = EQSB + EQNT - MYDF(dSB); }
799                 if(recons[dET]) { m[dWB] = EQET + EQWB - MYDF(dET); }
800                 if(recons[dEB]) { m[dWT] = EQEB + EQWT - MYDF(dEB); }
801                 if(recons[dWT]) { m[dEB] = EQWT + EQEB - MYDF(dWT); }
802                 if(recons[dWB]) { m[dET] = EQWB + EQET - MYDF(dWB); }
803 #               endif           
804
805
806                 // inflow bc handling
807                 if(oldFlag & (CFMbndInflow)) {
808                         // fill if cells in inflow region
809                         if(myfrac<0.5) { 
810                                 mass += 0.25; 
811                                 mInitialMass += 0.25;
812                         }
813                         const int OId = oldFlag>>24;
814                         const LbmVec vel(mObjectSpeeds[OId]);
815                         ux=vel[0], uy=vel[1], uz=vel[2]; 
816                         //? usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
817                         //FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); } rho = fluidRho;
818                         rho = REFERENCE_PRESSURE;
819                         FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
820                         //errMsg("INFLOW_DEBUG","if at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
821                 }  else { 
822                         // NEWSURFT, todo optimize!
823                         if(onlyBndnb) { //if(usqr<0.001*0.001) {
824                                 rho = ux = uy = uz = 0.;
825                                 FORDF0{ 
826                                         rho += m[l]; 
827                                         ux  += (this->dfDvecX[l]*m[l]); 
828                                         uy  += (this->dfDvecY[l]*m[l]);  
829                                         uz  += (this->dfDvecZ[l]*m[l]);  
830                                 }
831                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
832                         } else {// NEWSURFT */
833                                 if(usqr>0.3*0.3) { 
834                                         // force reset! , warning causes distortions...
835                                         FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,0.,0.,0.); }
836                                 } else {
837                                 // normal collide
838                                 // mass streaming done... do normal collide
839                                 LbmVec grav = mLevel[lev].gravity*mass;
840                                 DEFAULT_COLLIDEG(grav);
841                                 PERFORM_USQRMAXCHECK; }
842                                 // rho init from default collide necessary for fill/empty check below
843                         } // test
844                 }
845
846                 // testing..., particle generation
847                 // also check oldFlag for CFNoNbFluid, test
848                 // for inflow no pargen test // NOBUBBB!
849                 if((mInitDone) 
850                                 // dont allow new if cells, or submerged ones
851                                 && (!((oldFlag|newFlag)& (CFNoDelete|CFNoNbEmpty) )) 
852                                 // dont try to subtract from empty cells
853                                 && (mass>0.) && (mPartGenProb>0.0)) {
854                         bool doAdd = true;
855                         bool bndOk=true;
856                         if( (i<cutMin)||(i>mSizex-cutMin)||
857                                         (j<cutMin)||(j>mSizey-cutMin)||
858                                         (k<cutMin)||(k>mSizez-cutMin) ) { bndOk=false; }
859                         if(!bndOk) doAdd=false;
860                         
861                         LbmFloat prob = (rand()/(RAND_MAX+1.0));
862                         LbmFloat basethresh = mPartGenProb*lcsmqo*(LbmFloat)(mSizez+mSizey+mSizex)*0.5*0.333;
863
864                         // physical drop model
865                         if(mPartUsePhysModel) {
866                                 LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep);
867                                 LbmFloat rux = (ux * realWorldFac);
868                                 LbmFloat ruy = (uy * realWorldFac);
869                                 LbmFloat ruz = (uz * realWorldFac);
870                                 LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz));
871                                 basethresh *= rl;
872
873                                 // reduce probability in outer region?
874                                 const int pibord = mLevel[mMaxRefine].lSizex/2-cutConst;
875                                 const int pjbord = mLevel[mMaxRefine].lSizey/2-cutConst;
876                                 LbmFloat pifac = 1.-(LbmFloat)(ABS(i-pibord)) / (LbmFloat)(pibord);
877                                 LbmFloat pjfac = 1.-(LbmFloat)(ABS(j-pjbord)) / (LbmFloat)(pjbord);
878                                 if(pifac<0.) pifac=0.;
879                                 if(pjfac<0.) pjfac=0.;
880
881                                 //if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
882                                 if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
883                                         // add
884                                 } else {
885                                         doAdd = false; // dont...
886                                 }
887
888                                 // "wind" disturbance
889                                 // use realworld relative velocity here instead?
890                                 if( (doAdd && 
891                                                 ((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) )// normal checks
892                                                 ||(k>mSizez-SLOWDOWNREGION)   ) {
893                                         LbmFloat nuz = uz;
894                                         if(k>mSizez-SLOWDOWNREGION) {
895                                                 // special case
896                                                 LbmFloat zfac = (LbmFloat)( k-(mSizez-SLOWDOWNREGION) );
897                                                 zfac /= (LbmFloat)(SLOWDOWNREGION);
898                                                 nuz += (1.0) * zfac; // check max speed? OFF?
899                                                 //errMsg("TOPT"," at "<<PRINT_IJK<<" zfac"<<zfac);
900                                         } else {
901                                                 // normal probability
902                                                 //? LbmFloat fac = P_LCSMQO-lcsmqo;
903                                                 //? jdf *= fac;
904                                         }
905                                         FORDF1 {
906                                                 const LbmFloat jdf = 0.05 * (rand()/(RAND_MAX+1.0));
907                                                 // TODO  use wind velocity?
908                                                 if(jdf>0.025) {
909                                                 const LbmFloat add =  this->dfLength[l]*(-ux*this->dfDvecX[l]-uy*this->dfDvecY[l]-nuz*this->dfDvecZ[l])*jdf;
910                                                 RAC(tcel,l) += add; }
911                                         }
912                                         //errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) );
913                                 } // wind disturbance
914                         } // mPartUsePhysModel
915                         else {
916                                 // empirical model
917                                 //if((prob<basethresh) && (lcsmqo>0.0095)) { // add
918                                 if((prob<basethresh) && (lcsmqo>0.012)) { // add
919                                 } else { doAdd = false; }// dont...
920                         } 
921
922
923                         // remove noise
924                         if(usqr<0.0001) doAdd=false;   // TODO check!?
925
926                         // dont try to subtract from empty cells
927                         // ensure cell has enough mass for new drop
928                         LbmFloat newPartsize = 1.0;
929                         if(mPartUsePhysModel) {
930                                 // 1-10
931                                 newPartsize += 9.0* (rand()/(RAND_MAX+1.0));
932                         } else {
933                                 // 1-5, overall size has to be less than
934                                 // .62 (ca. 0.5) to make drops significantly smaller 
935                                 // than a full cell!
936                                 newPartsize += 4.0* (rand()/(RAND_MAX+1.0));
937                         }
938                         LbmFloat dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize); //PARTMASS(mPartDropMassSub*newPartsize); // mass: 4/3 pi r^3 rho
939                         while(dropmass>mass) {
940                                 newPartsize -= 0.2;
941                                 dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize);
942                         }
943                         if(newPartsize<=1.) doAdd=false;
944
945                         if( (doAdd)  ) { // init new particle
946                                 for(int s=0; s<1; s++) { // one part!
947                                 const LbmFloat posjitter = 0.05;
948                                 const LbmFloat posjitteroffs = posjitter*-0.5;
949                                 LbmFloat jpx = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
950                                 LbmFloat jpy = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
951                                 LbmFloat jpz = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
952
953                                 const LbmFloat jitterstr = 1.0;
954                                 const LbmFloat jitteroffs = jitterstr*-0.5;
955                                 LbmFloat jx = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
956                                 LbmFloat jy = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
957                                 LbmFloat jz = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
958
959                                 // average normal & velocity 
960                                 // -> mostly along velocity dir, many into surface
961                                 // fluid velocity (not normalized!)
962                                 LbmVec flvelVel = LbmVec(ux,uy,uz);
963                                 LbmFloat flvelLen = norm(flvelVel);
964                                 // surface normal
965                                 LbmVec normVel = LbmVec(surfaceNormal[0],surfaceNormal[1],surfaceNormal[2]);
966                                 normalize(normVel);
967                                 LbmFloat normScale = (0.01+flvelLen);
968                                 // jitter vector, 0.2 * flvel
969                                 LbmVec jittVel = LbmVec(jx,jy,jz)*(0.05+flvelLen)*0.1;
970                                 // weighten velocities
971                                 const LbmFloat flvelWeight = 0.9;
972                                 LbmVec newpartVel = normVel*normScale*(1.-flvelWeight) + flvelVel*(flvelWeight) + jittVel; 
973
974                                 // offset towards surface (hide popping)
975                                 jpx += -normVel[0]*0.4;
976                                 jpy += -normVel[1]*0.4;
977                                 jpz += -normVel[2]*0.4;
978
979                                 LbmFloat srci=i+0.5+jpx, srcj=j+0.5+jpy, srck=k+0.5+jpz;
980                                 int type=0;
981                                 type = PART_DROP;
982
983 #                               if LBMDIM==2
984                                 newpartVel[2]=0.; srck=0.5;
985 #                               endif // LBMDIM==2
986                                 // subtract drop mass
987                                 mass -= dropmass;
988                                 // init new particle
989                                 {
990                                         ParticleObject np( ntlVec3Gfx(srci,srcj,srck) );
991                                         np.setVel(newpartVel[0],newpartVel[1],newpartVel[2]);
992                                         np.setStatus(PART_IN);
993                                         np.setType(type);
994                                         //if((s%3)==2) np.setType(PART_FLOAT);
995                                         np.setSize(newPartsize);
996                                         //errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<norm(LbmVec(ux,uy,uz)) <<" add"<<doAdd<<" pvel="<<norm(newpartVel)<<" size="<<newPartsize );
997                                         //errMsg("NEWPT","u="<<newpartVel<<" norm="<<normVel<<" flvel="<<flvelVel<<" jitt="<<jittVel );
998                                         FSGR_ADDPART(np);
999                                 } // new part
1000                 //errMsg("PIT","NEW pit"<<np.getId()<<" pos:"<< np.getPos()<<" status:"<<convertFlags2String(np.getFlags())<<" vel:"<< np.getVel()  );
1001                                 } // multiple parts
1002                         } // doAdd
1003                 } // */
1004
1005
1006                 // interface cell filled or emptied?
1007                 iffilled = ifemptied = 0;
1008                 // interface cells empty/full?, WARNING: to mark these cells, better do it at the end of reinitCellFlags
1009                 // interface cell if full?
1010                 if( (mass) >= (rho * (1.0+FSGR_MAGICNR)) ) { iffilled = 1; }
1011                 // interface cell if empty?
1012                 if( (mass) <= (rho * (   -FSGR_MAGICNR)) ) { ifemptied = 1; }
1013
1014                 if(oldFlag & (CFMbndOutflow)) {
1015                         mInitialMass -= mass;
1016                         mass = myfrac = 0.0;
1017                         iffilled = 0; ifemptied = 1;
1018                 }
1019
1020                 // looks much nicer... LISTTRICK
1021 #               if FSGR_LISTTRICK==1
1022                 //if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
1023                 if(newFlag&CFNoBndFluid) { // test NEW TEST
1024                         if(!iffilled) {
1025                                 // remove cells independent from amount of change...
1026                                 if( (oldFlag & CFNoNbEmpty)&&(newFlag & CFNoNbEmpty)&&
1027                                                 ( (mass>(rho*FSGR_LISTTTHRESHFULL))  || ((nbored&CFInter)==0)  )) { 
1028                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","filled "<<PRINT_IJK); };
1029                                         iffilled = 1; 
1030                                 } 
1031                         }
1032                         if(!ifemptied) {
1033                                 if( (oldFlag & CFNoNbFluid)&&(newFlag & CFNoNbFluid)&&
1034                                                 ( (mass<(rho*FSGR_LISTTTHRESHEMPTY)) || ((nbored&CFInter)==0)  )) { 
1035                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","emptied "<<PRINT_IJK); };
1036                                         ifemptied = 1; 
1037                                 } 
1038                         }
1039                 } // nobndfluid only */
1040 #               endif
1041                 //iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!!
1042                 
1043
1044                 // now that all dfs are known, handle last changes
1045                 if(iffilled) {
1046                         LbmPoint filledp; filledp.flag=0;
1047                         if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1;  // NEWSURFT
1048                         filledp.x = i; filledp.y = j; filledp.z = k;
1049                         LIST_FULL(filledp);
1050                         //mNumFilledCells++; // DEBUG
1051                         calcCellsFilled++;
1052                 }
1053                 else if(ifemptied) {
1054                         LbmPoint emptyp; emptyp.flag=0;
1055                         if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; //  NEWSURFT
1056                         emptyp.x = i; emptyp.y = j; emptyp.z = k;
1057                         LIST_EMPTY(emptyp);
1058                         //mNumEmptiedCells++; // DEBUG
1059                         calcCellsEmptied++;
1060                 } 
1061                 // dont cutoff values -> better cell conversions
1062                 RAC(tcel,dFfrac)   = (mass/rho);
1063
1064                 // init new flux value
1065                 float flux = FLUX_INIT; // dxqn on
1066                 if(newFlag&CFNoBndFluid) {
1067                         //flux = 50.0; // extreme on
1068                         for(int nn=1; nn<this->cDfNum; nn++) { 
1069                                 if(nbflag[nn] & (CFFluid|CFInter|CFBnd)) { flux += this->dfLength[nn]; }
1070                         }
1071                         // optical hack - smooth slow moving
1072                         // surface regions
1073                         if(usqr< sssUsqrLimit) {
1074                         for(int nn=1; nn<this->cDfNum; nn++) { 
1075                                 if(nbfracs[nn]!=0.0) {
1076                                         LbmFloat curSmooth = (sssUsqrLimit-usqr)*sssUsqrLimitInv;
1077                                         if(curSmooth>1.0) curSmooth=1.0;
1078                                         flux *= (1.0+ smoothStrength*curSmooth * (nbfracs[nn]-myfrac)) ;
1079                                 }
1080                         } }
1081                         // NEW TEST */
1082                 }
1083                 // flux = FLUX_INIT; // calc flux off
1084                 QCELL(lev, i,j,k,TSET(lev), dFlux) = flux; // */
1085
1086                 // perform mass exchange with streamed values
1087                 QCELL(lev, i,j,k,TSET(lev), dMass) = mass; // MASST
1088                 // set new flag 
1089                 *pFlagDst = (CellFlagType)newFlag;
1090                 calcCurrentMass += mass; 
1091                 calcCurrentVolume += RAC(tcel,dFfrac);
1092
1093                 // interface cell handling done...
1094
1095 #if PARALLEL!=1
1096         GRID_LOOPREG_END();
1097 #else // PARALLEL==1
1098 #include "paraloopend.h" // = GRID_LOOPREG_END();
1099 #endif // PARALLEL==1
1100
1101         // write vars from computations to class
1102         mLevel[lev].lmass    = calcCurrentMass;
1103         mLevel[lev].lvolume  = calcCurrentVolume;
1104         mNumFilledCells  = calcCellsFilled;
1105         mNumEmptiedCells = calcCellsEmptied;
1106         mNumUsedCells = calcNumUsedCells;
1107 }
1108
1109
1110
1111 void 
1112 LbmFsgrSolver::preinitGrids()
1113 {
1114         const int lev = mMaxRefine;
1115         const bool doReduce = false;
1116         const int gridLoopBound=0;
1117
1118         // preinit both grids
1119         for(int s=0; s<2; s++) {
1120         
1121                 GRID_REGION_INIT();
1122 #if PARALLEL==1
1123 #pragma omp parallel default(shared) \
1124   reduction(+: \
1125           calcCurrentMass,calcCurrentVolume, \
1126                 calcCellsFilled,calcCellsEmptied, \
1127                 calcNumUsedCells )
1128 #endif // PARALLEL==1
1129                 GRID_REGION_START();
1130                 GRID_LOOP_START();
1131                         for(int l=0; l<dTotalNum; l++) { RAC(ccel,l) = 0.; }
1132                         *pFlagSrc =0;
1133                         *pFlagDst =0;
1134                         //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
1135 #if PARALLEL!=1
1136                 GRID_LOOPREG_END();
1137 #else // PARALLEL==1
1138 #include "paraloopend.h" // = GRID_LOOPREG_END();
1139 #endif // PARALLEL==1
1140
1141                 /* dummy, remove warnings */ 
1142                 calcCurrentMass = calcCurrentVolume = 0.;
1143                 calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
1144                 
1145                 // change grid
1146                 mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
1147                 mLevel[mMaxRefine].setCurr   ^= 1;
1148         }
1149 }
1150
1151 void 
1152 LbmFsgrSolver::standingFluidPreinit()
1153 {
1154         const int lev = mMaxRefine;
1155         const bool doReduce = false;
1156         const int gridLoopBound=1;
1157
1158         GRID_REGION_INIT();
1159 #if PARALLEL==1
1160 #pragma omp parallel default(shared) \
1161   reduction(+: \
1162           calcCurrentMass,calcCurrentVolume, \
1163                 calcCellsFilled,calcCellsEmptied, \
1164                 calcNumUsedCells )
1165 #endif // PARALLEL==1
1166         GRID_REGION_START();
1167
1168         LbmFloat rho, ux,uy,uz, usqr; 
1169         CellFlagType nbflag[LBM_DFNUM];
1170         LbmFloat m[LBM_DFNUM];
1171         LbmFloat lcsmqo;
1172 #       if OPT3D==1 
1173         LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
1174         CellFlagType nbored=0;
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