13ebf91b6961b4483d7fd5dcf66d5e0b609cd532
[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 (__sun__) || defined (__sun)) || (!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
373 #       if LBM_INCLUDE_TESTSOLVERS==1
374         // 3d region off... quit
375         if((mUseTestdata)&&(mpTest->mFarfMode>0)) { return; }
376 #       endif // ELBEEM_PLUGIN!=1
377         
378   // main loop region
379         const bool doReduce = true;
380         const int gridLoopBound=1;
381         GRID_REGION_INIT();
382 #if PARALLEL==1
383 #pragma omp parallel default(shared) \
384   reduction(+: \
385           calcCurrentMass,calcCurrentVolume, \
386                 calcCellsFilled,calcCellsEmptied, \
387                 calcNumUsedCells )
388         GRID_REGION_START();
389 #else // PARALLEL==1
390         GRID_REGION_START();
391 #endif // PARALLEL==1
392
393         // local to main
394         CellFlagType nbflag[LBM_DFNUM]; 
395         int oldFlag, newFlag, nbored;
396         LbmFloat m[LBM_DFNUM];
397         LbmFloat rho, ux, uy, uz, tmp, usqr;
398
399         // smago vars
400         LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
401         
402         // ifempty cell conversion flags
403         bool iffilled, ifemptied;
404         LbmFloat nbfracs[LBM_DFNUM]; // ffracs of neighbors
405         int recons[LBM_DFNUM];   // reconstruct this DF?
406         int numRecons;           // how many are reconstructed?
407
408         LbmFloat mass=0., change=0., lcsmqo=0.;
409         rho= ux= uy= uz= usqr= tmp= 0.; 
410         lcsmqadd = lcsmomega = 0.;
411         FORDF0{ lcsmeq[l] = 0.; }
412
413         // ---
414         // now stream etc.
415         // use //template functions for 2D/3D
416
417         GRID_LOOP_START();
418                 // loop start
419                 // stream from current set to other, then collide and store
420                 //errMsg("l2"," at "<<PRINT_IJK<<" id"<<id);
421
422 #               if FSGR_STRICT_DEBUG==1
423                 // safety pointer check
424                 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
425                 if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) || 
426                     (&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) {
427                         errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<<
428                                         RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<<
429                                         (*pFlagSrc)<<","<<(*pFlagDst) <<",  pointers "<<
430           (long)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(long)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
431           (long)(pFlagSrc)<<","<<(long)(pFlagDst)<<" "
432                                         ); 
433                         CAUSE_PANIC;
434                 }       
435                 if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) || 
436                     (&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) {
437                         errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<<
438           (long)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(long)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
439           (long)(ccel)<<","<<(long)(tcel)<<" "
440                                         ); 
441                         CAUSE_PANIC;
442                 }       
443 #               endif
444                 oldFlag = *pFlagSrc;
445                 
446                 // old INTCFCOARSETEST==1
447                 if( (oldFlag & (CFGrFromCoarse)) ) { 
448                         if(( mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
449                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
450                         } else {
451                                 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
452                                 calcNumUsedCells++;
453                         }
454                         continue; // interpolateFineFromCoarse test!
455                 } // interpolateFineFromCoarse test! 
456         
457                 if(oldFlag & (CFMbndInflow)) {
458                         // fluid & if are ok, fill if later on
459                         int isValid = oldFlag & (CFFluid|CFInter);
460                         const LbmFloat iniRho = 1.0;
461                         const int OId = oldFlag>>24;
462                         if(!isValid) {
463                                 // make new if cell
464                                 const LbmVec vel(mObjectSpeeds[OId]);
465                                 // TODO add OPT3D treatment
466                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
467                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
468                                 RAC(tcel, dFlux) = FLUX_INIT;
469                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
470                                 calcCurrentMass += iniRho; 
471                                 calcCurrentVolume += 1.0; 
472                                 calcNumUsedCells++;
473                                 mInitialMass += iniRho;
474                                 // dont treat cell until next step
475                                 continue;
476                         } 
477                 } 
478                 else  // these are exclusive
479                 if(oldFlag & (CFMbndOutflow)) {
480                         int isnotValid = oldFlag & (CFFluid);
481                         if(isnotValid) {
482                                 // remove fluid cells, shouldnt be here anyway
483                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
484                                 mInitialMass -= fluidRho;
485                                 const LbmFloat iniRho = 0.0;
486                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
487                                 RAC(tcel, dFlux) = FLUX_INIT;
488                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
489
490                                 // same as ifemptied for if below
491                                 LbmPoint oemptyp; oemptyp.flag = 0;
492                                 oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
493                                 LIST_EMPTY(oemptyp);
494                                 calcCellsEmptied++;
495                                 continue;
496                         }
497                 }
498
499                 if(oldFlag & (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)) { 
500                         *pFlagDst = oldFlag;
501                         continue;
502                 }
503                 /*if( oldFlag & CFNoBndFluid ) {  // TEST ME FASTER?
504                         OPTIMIZED_STREAMCOLLIDE; PERFORM_USQRMAXCHECK;
505                         RAC(tcel,dFfrac) = 1.0; 
506                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
507                         calcCurrentMass += rho; calcCurrentVolume += 1.0;
508                         calcNumUsedCells++;
509                         continue;
510                 }// TEST ME FASTER? */
511
512                 // only neighbor flags! not own flag
513                 nbored = 0;
514                 
515 #if OPT3D==0
516                 FORDF1 {
517                         nbflag[l] = RFLAG_NB(lev, i,j,k,SRCS(lev),l);
518                         nbored |= nbflag[l];
519                 } 
520 #else
521                 nbflag[dSB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dSB];
522                 nbflag[dWB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWB];
523                 nbflag[ dB] = *(pFlagSrc + (-mLevel[lev].lOffsy)); nbored |= nbflag[dB];
524                 nbflag[dEB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dEB];
525                 nbflag[dNB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNB];
526
527                 nbflag[dSW] = *(pFlagSrc + (-mLevel[lev].lOffsx+-1)); nbored |= nbflag[dSW];
528                 nbflag[ dS] = *(pFlagSrc + (-mLevel[lev].lOffsx)); nbored |= nbflag[dS];
529                 nbflag[dSE] = *(pFlagSrc + (-mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dSE];
530
531                 nbflag[ dW] = *(pFlagSrc + (-1)); nbored |= nbflag[dW];
532                 nbflag[ dE] = *(pFlagSrc + ( 1)); nbored |= nbflag[dE];
533
534                 nbflag[dNW] = *(pFlagSrc + ( mLevel[lev].lOffsx+-1)); nbored |= nbflag[dNW];
535           nbflag[ dN] = *(pFlagSrc + ( mLevel[lev].lOffsx)); nbored |= nbflag[dN];
536                 nbflag[dNE] = *(pFlagSrc + ( mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dNE];
537
538                 nbflag[dST] = *(pFlagSrc + ( mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dST];
539                 nbflag[dWT] = *(pFlagSrc + ( mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWT];
540                 nbflag[ dT] = *(pFlagSrc + ( mLevel[lev].lOffsy)); nbored |= nbflag[dT];
541                 nbflag[dET] = *(pFlagSrc + ( mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dET];
542                 nbflag[dNT] = *(pFlagSrc + ( mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNT];
543                 // */
544 #endif
545
546                 // pointer to destination cell
547                 calcNumUsedCells++;
548
549                 // FLUID cells 
550                 if( oldFlag & CFFluid ) { 
551                         // only standard fluid cells (with nothing except fluid as nbs
552
553                         if(oldFlag&CFMbndInflow) {
554                                 // force velocity for inflow, necessary to have constant direction of flow
555                                 // FIXME , test also set interface cells?
556                                 const int OId = oldFlag>>24;
557                                 //? DEFAULT_STREAM;
558                                 //const LbmFloat fluidRho = 1.0;
559                                 // for submerged inflows, streaming would have to be performed...
560                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
561                                 const LbmVec vel(mObjectSpeeds[OId]);
562                                 ux=vel[0], uy=vel[1], uz=vel[2]; 
563                                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
564                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); }
565                                 rho = fluidRho;
566                                 //errMsg("INFLOW_DEBUG","std at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
567                         } else {
568                                 if(nbored&CFBnd) {
569                                         DEFAULT_STREAM;
570                                         //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
571                                         DEFAULT_COLLIDEG(mLevel[lev].gravity);
572                                         oldFlag &= (~CFNoBndFluid);
573                                 } else {
574                                         // do standard stream/collide
575                                         OPTIMIZED_STREAMCOLLIDE;
576                                         oldFlag |= CFNoBndFluid;
577                                 } 
578                         }
579
580                         PERFORM_USQRMAXCHECK;
581                         // "normal" fluid cells
582                         RAC(tcel,dFfrac) = 1.0; 
583                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
584                         calcCurrentMass += rho; 
585                         calcCurrentVolume += 1.0;
586                         continue;
587                 }
588                 
589                 newFlag  = oldFlag;
590                 // make sure here: always check which flags to really unset...!
591                 newFlag = newFlag & (~( 
592                                         CFNoNbFluid|CFNoNbEmpty| CFNoDelete 
593                                         | CFNoInterpolSrc
594                                         | CFNoBndFluid
595                                         ));
596                 if(!(nbored&CFBndNoslip)) { //NEWSURFT NEWSURFTNOS
597                         newFlag |= CFNoBndFluid;
598                 }
599                 /*if(!(nbored&CFBnd)) { //NEWSURFT NEWSURFTNOS
600                         // explicitly check for noslip neighbors
601                         bool hasnoslipnb = false;
602                         FORDF1 { if((nbflag[l]&CFBnd)&&(nbflag[l]&CFBndNoslip)) hasnoslipnb=true; }
603                         if(!hasnoslipnb) newFlag |= CFNoBndFluid; 
604                 } // */
605
606                 // store own dfs and mass
607                 mass = RAC(ccel,dMass);
608
609                 // WARNING - only interface cells arrive here!
610                 // read distribution funtions of adjacent cells = stream step
611                 DEFAULT_STREAM;
612
613                 if((nbored & CFFluid)==0) { newFlag |= CFNoNbFluid; mNumInvIfCells++; }
614                 if((nbored & CFEmpty)==0) { newFlag |= CFNoNbEmpty; mNumInvIfCells++; }
615
616                 // calculate mass exchange for interface cells 
617                 LbmFloat myfrac = RAC(ccel,dFfrac);
618                 if(myfrac<0.) myfrac=0.; // NEWSURFT
619 #               define nbdf(l) m[ this->dfInv[(l)] ]
620
621                 // update mass 
622                 // only do boundaries for fluid cells, and interface cells without
623                 // any fluid neighbors (assume that interface cells _with_ fluid
624                 // neighbors are affected enough by these) 
625                 // which Df's have to be reconstructed? 
626                 // for fluid cells - just the f_i difference from streaming to empty cells  ----
627                 numRecons = 0;
628                 bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
629                 //onlyBndnb = false; // DEBUG test off
630
631                 FORDF1 { // dfl loop
632                         recons[l] = 0;
633                         nbfracs[l] = 0.0;
634                         // finally, "normal" interface cells ----
635                         if( nbflag[l]&(CFFluid|CFBnd) ) { // NEWTEST! FIXME check!!!!!!!!!!!!!!!!!!
636                                 change = nbdf(l) - MYDF(l);
637                         }
638                         // interface cells - distuingish cells that shouldn't fill/empty 
639                         else if( nbflag[l] & CFInter ) {
640                                 
641                                 LbmFloat mynbfac,nbnbfac;
642                                 // NEW TEST t1
643                                 // t2 -> off
644                                 if((oldFlag&CFNoBndFluid)&&(nbflag[l]&CFNoBndFluid)) {
645                                         mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
646                                         nbnbfac = 1.0/mynbfac;
647                                         onlyBndnb = false;
648                                 } else {
649                                         mynbfac = nbnbfac = 1.0; // switch calc flux off
650                                         goto changeDefault;  // NEWSURFT
651                                         //change = 0.; goto changeDone;  // NEWSURFT
652                                 }
653                                 //mynbfac = nbnbfac = 1.0; // switch calc flux off t3
654
655                                 // perform interface case handling
656                                 if ((oldFlag|nbflag[l])&(CFNoNbFluid|CFNoNbEmpty)) {
657                                 switch (oldFlag&(CFNoNbFluid|CFNoNbEmpty)) {
658                                         case 0: 
659                                                 // we are a normal cell so... 
660                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
661                                                         case CFNoNbFluid: 
662                                                                 // just fill current cell = empty neighbor 
663                                                                 change = nbnbfac*nbdf(l) ; goto changeDone; 
664                                                         case CFNoNbEmpty: 
665                                                                 // just empty current cell = fill neighbor 
666                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
667                                                 }
668                                                 break;
669
670                                         case CFNoNbFluid: 
671                                                 // we dont have fluid nb's so...
672                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
673                                                         case 0: 
674                                                         case CFNoNbEmpty: 
675                                                                 // we have no fluid nb's -> just empty
676                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
677                                                 }
678                                                 break;
679
680                                         case CFNoNbEmpty: 
681                                                 // we dont have empty nb's so...
682                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
683                                                         case 0: 
684                                                         case CFNoNbFluid: 
685                                                                 // we have no empty nb's -> just fill
686                                                                 change = nbnbfac*nbdf(l); goto changeDone; 
687                                                 }
688                                                 break;
689                                 }} // inter-inter exchange
690
691                         changeDefault: ;
692                                 // just do normal mass exchange...
693                                 change = ( nbnbfac*nbdf(l) - mynbfac*MYDF(l) ) ;
694                         changeDone: ;
695                                 nbfracs[l] = QCELL_NB(lev, i,j,k, SRCS(lev),l, dFfrac);
696                                 if(nbfracs[l]<0.) nbfracs[l] = 0.; // NEWSURFT
697                                 change *=  (myfrac + nbfracs[l]) * 0.5;
698                         } // the other cell is interface
699
700                         // last alternative - reconstruction in this direction
701                         else {
702                                 // empty + bnd case
703                                 recons[l] = 1; 
704                                 numRecons++;
705                                 change = 0.0; 
706                         }
707
708                         // modify mass at SRCS
709                         mass += change;
710                 } // l
711                 // normal interface, no if empty/fluid
712
713                 // computenormal
714                 LbmFloat surfaceNormal[3];
715                 computeFluidSurfaceNormal(ccel,pFlagSrc, surfaceNormal);
716
717                 if( (ABS(surfaceNormal[0])+ABS(surfaceNormal[1])+ABS(surfaceNormal[2])) > LBM_EPSILON) {
718                         // normal ok and usable...
719                         FORDF1 {
720                                 if( (this->dfDvecX[l]*surfaceNormal[0] + this->dfDvecY[l]*surfaceNormal[1] + this->dfDvecZ[l]*surfaceNormal[2])  // dot Dvec,norml
721                                                 > LBM_EPSILON) {
722                                         recons[l] = 2; 
723                                         numRecons++;
724                                 } 
725                         }
726                 }
727
728                 // calculate macroscopic cell values
729                 LbmFloat oldUx, oldUy, oldUz;
730                 LbmFloat oldRho; // OLD rho = ccel->rho;
731 #               define REFERENCE_PRESSURE 1.0 // always atmosphere...
732 #               if OPT3D==0
733                 oldRho=RAC(ccel,0);
734                 oldUx = oldUy = oldUz = 0.0;
735                 for(int l=1; l<this->cDfNum; l++) {
736                         oldRho += RAC(ccel,l);
737                         oldUx  += (this->dfDvecX[l]*RAC(ccel,l));
738                         oldUy  += (this->dfDvecY[l]*RAC(ccel,l)); 
739                         oldUz  += (this->dfDvecZ[l]*RAC(ccel,l)); 
740                 } 
741                 // reconstruct dist funcs from empty cells
742                 FORDF1 {
743                         if(recons[ l ]) {
744                                 m[ this->dfInv[l] ] = 
745                                         this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 
746                                         this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 
747                                         - MYDF( l );
748                         }
749                 }
750                 ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
751                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
752 #               else // OPT3D==0
753                 oldRho = + RAC(ccel,dC)  + RAC(ccel,dN )
754                                 + RAC(ccel,dS ) + RAC(ccel,dE )
755                                 + RAC(ccel,dW ) + RAC(ccel,dT )
756                                 + RAC(ccel,dB ) + RAC(ccel,dNE)
757                                 + RAC(ccel,dNW) + RAC(ccel,dSE)
758                                 + RAC(ccel,dSW) + RAC(ccel,dNT)
759                                 + RAC(ccel,dNB) + RAC(ccel,dST)
760                                 + RAC(ccel,dSB) + RAC(ccel,dET)
761                                 + RAC(ccel,dEB) + RAC(ccel,dWT)
762                                 + RAC(ccel,dWB);
763
764                 oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
765                                 + RAC(ccel,dNE) - RAC(ccel,dNW)
766                                 + RAC(ccel,dSE) - RAC(ccel,dSW)
767                                 + RAC(ccel,dET) + RAC(ccel,dEB)
768                                 - RAC(ccel,dWT) - RAC(ccel,dWB);
769
770                 oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
771                                 + RAC(ccel,dNE) + RAC(ccel,dNW)
772                                 - RAC(ccel,dSE) - RAC(ccel,dSW)
773                                 + RAC(ccel,dNT) + RAC(ccel,dNB)
774                                 - RAC(ccel,dST) - RAC(ccel,dSB);
775
776                 oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
777                                 + RAC(ccel,dNT) - RAC(ccel,dNB)
778                                 + RAC(ccel,dST) - RAC(ccel,dSB)
779                                 + RAC(ccel,dET) - RAC(ccel,dEB)
780                                 + RAC(ccel,dWT) - RAC(ccel,dWB);
781
782                 // 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) \
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) \
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         CellFlagType nbored=0;
1176 #       endif // OPT3D==true 
1177
1178         GRID_LOOP_START();
1179                 //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
1180                 const CellFlagType currFlag = *pFlagSrc; //RFLAG(lev, i,j,k,workSet);
1181                 if( (currFlag & (CFEmpty|CFBnd)) ) continue;
1182
1183                 if( (currFlag & (CFInter)) ) {
1184                         // copy all values
1185                         for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
1186                         continue;
1187                 }
1188
1189                 if( (currFlag & CFNoBndFluid)) {
1190                         OPTIMIZED_STREAMCOLLIDE;
1191                 } else {
1192                         FORDF1 {
1193                                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
1194                         } 
1195                         DEFAULT_STREAM;
1196                         DEFAULT_COLLIDEG(mLevel[lev].gravity);
1197                 }
1198                 for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
1199 #if PARALLEL!=1
1200         GRID_LOOPREG_END();
1201 #else // PARALLEL==1
1202 #include "paraloopend.h" // = GRID_LOOPREG_END();
1203 #endif // PARALLEL==1
1204
1205         /* dummy remove warnings */ 
1206         calcCurrentMass = calcCurrentVolume = 0.;
1207         calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
1208         
1209         // change grid
1210   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
1211   mLevel[mMaxRefine].setCurr   ^= 1;
1212 }
1213
1214
1215 /******************************************************************************
1216  * work on lists from updateCellMass to reinit cell flags
1217  *****************************************************************************/
1218
1219 LbmFloat LbmFsgrSolver::getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l) {
1220         //return 0.0; // test
1221         int level = mMaxRefine;
1222         LbmFloat     *ccel  = RACPNT(level, i,j,k, workSet);
1223
1224         // computenormal
1225         CellFlagType *cflag = &RFLAG(level, i,j,k, workSet);
1226         LbmFloat n[3];
1227         computeFluidSurfaceNormal(ccel,cflag, n);
1228         LbmFloat scal = mDvecNrm[l][0]*n[0] + mDvecNrm[l][1]*n[1] + mDvecNrm[l][2]*n[2];
1229
1230         LbmFloat ret = 1.0;
1231         // forward direction, add mass (for filling cells):
1232         if(dirForw) {
1233                 if(scal<LBM_EPSILON) ret = 0.0;
1234                 else ret = scal;
1235         } else {
1236                 // backward for emptying
1237                 if(scal>-LBM_EPSILON) ret = 0.0;
1238                 else ret = scal * -1.0;
1239         }
1240         //errMsg("massd", PRINT_IJK<<" nv"<<nvel<<" : ret="<<ret ); //xit(1); //VECDEB
1241         return ret;
1242 }
1243
1244 // warning - normal compuations are without
1245 //   boundary checks &
1246 //   normalization
1247 void LbmFsgrSolver::computeFluidSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt,LbmFloat *snret) {
1248         const int level = mMaxRefine;
1249         LbmFloat nx,ny,nz, nv1,nv2;
1250         const CellFlagType flagE = *(cflagpnt+1);
1251         const CellFlagType flagW = *(cflagpnt-1);
1252         if(flagE &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } 
1253         else if(flagE &(CFBnd)){ nv1 = 1.; }
1254         else nv1 = 0.0;
1255         if(flagW &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } 
1256         else if(flagW &(CFBnd)){ nv2 = 1.; }
1257         else nv2 = 0.0;
1258         nx = 0.5* (nv2-nv1);
1259
1260         const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
1261         const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
1262         if(flagN &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 
1263         else if(flagN &(CFBnd)){ nv1 = 1.; }
1264         else nv1 = 0.0;
1265         if(flagS &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 
1266         else if(flagS &(CFBnd)){ nv2 = 1.; }
1267         else nv2 = 0.0;
1268         ny = 0.5* (nv2-nv1);
1269
1270 #if LBMDIM==3
1271         const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
1272         const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
1273         if(flagT &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 
1274         else if(flagT &(CFBnd)){ nv1 = 1.; }
1275         else nv1 = 0.0;
1276         if(flagB &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 
1277         else if(flagB &(CFBnd)){ nv2 = 1.; }
1278         else nv2 = 0.0;
1279         nz = 0.5* (nv2-nv1);
1280 #else //LBMDIM==3
1281         nz = 0.0;
1282 #endif //LBMDIM==3
1283
1284         // return vals
1285         snret[0]=nx; snret[1]=ny; snret[2]=nz;
1286 }
1287 void LbmFsgrSolver::computeFluidSurfaceNormalAcc(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
1288         LbmFloat nx=0.,ny=0.,nz=0.;
1289         ccel = NULL; cflagpnt=NULL; // remove warning
1290         snret[0]=nx; snret[1]=ny; snret[2]=nz;
1291 }
1292 void LbmFsgrSolver::computeObstacleSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) {
1293         const int level = mMaxRefine;
1294         LbmFloat nx,ny,nz, nv1,nv2;
1295         ccel = NULL; // remove warning
1296
1297         const CellFlagType flagE = *(cflagpnt+1);
1298         const CellFlagType flagW = *(cflagpnt-1);
1299         if(flagE &(CFBnd)){ nv1 = 1.; }
1300         else nv1 = 0.0;
1301         if(flagW &(CFBnd)){ nv2 = 1.; }
1302         else nv2 = 0.0;
1303         nx = 0.5* (nv2-nv1);
1304
1305         const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx);
1306         const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx);
1307         if(flagN &(CFBnd)){ nv1 = 1.; }
1308         else nv1 = 0.0;
1309         if(flagS &(CFBnd)){ nv2 = 1.; }
1310         else nv2 = 0.0;
1311         ny = 0.5* (nv2-nv1);
1312
1313 #if LBMDIM==3
1314         const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy);
1315         const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy);
1316         if(flagT &(CFBnd)){ nv1 = 1.; }
1317         else nv1 = 0.0;
1318         if(flagB &(CFBnd)){ nv2 = 1.; }
1319         else nv2 = 0.0;
1320         nz = 0.5* (nv2-nv1);
1321 #else //LBMDIM==3
1322         nz = 0.0;
1323 #endif //LBMDIM==3
1324
1325         // return vals
1326         snret[0]=nx; snret[1]=ny; snret[2]=nz;
1327 }
1328 void LbmFsgrSolver::computeObstacleSurfaceNormalAcc(int i,int j,int k, LbmFloat *snret) {
1329         bool nonorm = false;
1330         LbmFloat nx=0.,ny=0.,nz=0.;
1331         if(i<=0)        { nx =  1.; nonorm = true; }
1332         if(i>=mSizex-1) { nx = -1.; nonorm = true; }
1333         if(j<=0)        { ny =  1.; nonorm = true; }
1334         if(j>=mSizey-1) { ny = -1.; nonorm = true; }
1335 #       if LBMDIM==3
1336         if(k<=0)        { nz =  1.; nonorm = true; }
1337         if(k>=mSizez-1) { nz = -1.; nonorm = true; }
1338 #       endif // LBMDIM==3
1339         if(!nonorm) {
1340                 // in domain, revert to helper, use setCurr&mMaxRefine
1341                 LbmVec bnormal;
1342                 CellFlagType *bflag = &RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr);
1343                 LbmFloat     *bcell = RACPNT(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr);
1344                 computeObstacleSurfaceNormal(bcell,bflag, &bnormal[0]);
1345                 // TODO check if there is a normal near here?
1346                 // use wider range otherwise...
1347                 snret[0]=bnormal[0]; snret[1]=bnormal[0]; snret[2]=bnormal[0];
1348                 return;
1349         }
1350         snret[0]=nx; snret[1]=ny; snret[2]=nz;
1351 }
1352
1353 void LbmFsgrSolver::addToNewInterList( int ni, int nj, int nk ) {
1354 #if FSGR_STRICT_DEBUG==10
1355         // dangerous, this can change the simulation...
1356   /*for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1357        iter != mListNewInter.end(); iter++ ) {
1358     if(ni!=iter->x) continue;
1359     if(nj!=iter->y) continue;
1360     if(nk!=iter->z) continue;
1361                 // all 3 values match... skip point
1362                 return;
1363         } */
1364 #endif // FSGR_STRICT_DEBUG==1
1365         // store point
1366         LbmPoint newinter; newinter.flag = 0;
1367         newinter.x = ni; newinter.y = nj; newinter.z = nk;
1368         mListNewInter.push_back(newinter);
1369 }
1370
1371 void LbmFsgrSolver::reinitFlags( int workSet ) { 
1372         // reinitCellFlags OLD mods:
1373         // add all to intel list?
1374         // check ffrac for new cells
1375         // new if cell inits (last loop)
1376         // vweights handling
1377
1378         const int debugFlagreinit = 0;
1379         
1380         // some things need to be read/modified on the other set
1381         int otherSet = (workSet^1);
1382         // fixed level on which to perform 
1383         int workLev = mMaxRefine;
1384
1385   /* modify interface cells from lists */
1386   /* mark filled interface cells as fluid, or emptied as empty */
1387         /* count neighbors and distribute excess mass to interface neighbor cells
1388    * problems arise when there are no interface neighbors anymore
1389          * then just distribute to any fluid neighbors...
1390          */
1391
1392         // for symmetry, first init all neighbor cells */
1393         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1394        iter != mListFull.end(); iter++ ) {
1395     int i=iter->x, j=iter->y, k=iter->z;
1396                 if(debugFlagreinit) errMsg("FULL", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
1397     FORDF1 {
1398                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1399                         //if((LBMDIM>2)&&( (ni<=0) || (nj<=0) || (nk<=0) || (ni>=mLevel[workLev].lSizex-1) || (nj>=mLevel[workLev].lSizey-1) || (nk>=mLevel[workLev].lSizez-1) )) {
1400                         if( (ni<=0) || (nj<=0) || 
1401                                   (ni>=mLevel[workLev].lSizex-1) ||
1402                                   (nj>=mLevel[workLev].lSizey-1) 
1403 #                                       if LBMDIM==3
1404                                   || (nk<=0) || (nk>=mLevel[workLev].lSizez-1) 
1405 #                                       endif // LBMDIM==1
1406                                  ) {
1407                                 continue; } // new bc, dont treat cells on boundary NEWBC
1408       if( RFLAG(workLev, ni,nj,nk, workSet) & CFEmpty ){
1409                                 
1410                                 // preinit speed, get from average surrounding cells
1411                                 // interpolate from non-workset to workset, sets are handled in function
1412
1413                                 // new and empty interface cell, dont change old flag here!
1414                                 addToNewInterList(ni,nj,nk);
1415
1416                                 LbmFloat avgrho = 0.0;
1417                                 LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
1418                                 interpolateCellValues(workLev,ni,nj,nk,workSet, avgrho,avgux,avguy,avguz);
1419
1420                                 // careful with l's...
1421                                 FORDF0M { 
1422                                         QCELL(workLev,ni,nj,nk, workSet, m) = this->getCollideEq( m,avgrho,  avgux, avguy, avguz ); 
1423                                         //QCELL(workLev,ni,nj,nk, workSet, l) = avgnbdf[l]; // CHECK FIXME test?
1424                                 }
1425                                 //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
1426                                 QCELL(workLev,ni,nj,nk, workSet, dMass) = 0.0; //?? new
1427                                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) = 0.0; //?? new
1428                                 //RFLAG(workLev,ni,nj,nk,workSet) = (CellFlagType)(CFInter|CFNoInterpolSrc);
1429                                 changeFlag(workLev,ni,nj,nk,workSet, (CFInter|CFNoInterpolSrc));
1430                                 if(debugFlagreinit) errMsg("NEWE", PRINT_IJK<<" newif "<<PRINT_VEC(ni,nj,nk)<<" rho"<<avgrho<<" vel("<<avgux<<","<<avguy<<","<<avguz<<") " );
1431       }
1432                         /* prevent surrounding interface cells from getting removed as empty cells 
1433                          * (also cells that are not newly inited) */
1434       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1435                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete);
1436                                 changeFlag(workLev,ni,nj,nk, workSet, (RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete));
1437                                 // also add to list...
1438                                 addToNewInterList(ni,nj,nk);
1439                         } // NEW?
1440     }
1441
1442                 // NEW? no extra loop...
1443                 //RFLAG(workLev,i,j,k, workSet) = CFFluid;
1444                 changeFlag(workLev,i,j,k, workSet,CFFluid);
1445         }
1446
1447         /* remove empty interface cells that are not allowed to be removed anyway
1448          * this is important, otherwise the dreaded cell-type-flickering can occur! */
1449   for(int n=0; n<(int)mListEmpty.size(); n++) {
1450     int i=mListEmpty[n].x, j=mListEmpty[n].y, k=mListEmpty[n].z;
1451                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) {
1452                         // treat as "new inter"
1453                         addToNewInterList(i,j,k);
1454                         // remove entry
1455                         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
1456                         if(n<(int)mListEmpty.size()-1) mListEmpty[n] = mListEmpty[mListEmpty.size()-1]; 
1457                         mListEmpty.pop_back();
1458                         n--; 
1459                 }
1460         } 
1461
1462
1463         /* problems arise when adjacent cells empty&fill ->
1464                  let fill cells+surrounding interface cells have the higher importance */
1465   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1466        iter != mListEmpty.end(); iter++ ) {
1467     int i=iter->x, j=iter->y, k=iter->z;
1468                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)){ errMsg("A"," ARGHARGRAG "); } // DEBUG
1469                 if(debugFlagreinit) errMsg("EMPT", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) );
1470
1471                 /* set surrounding fluid cells to interface cells */
1472     FORDF1 {
1473                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1474       if( RFLAG(workLev,ni,nj,nk, workSet) & CFFluid){
1475                                 // init fluid->interface 
1476                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(CFInter); 
1477                                 changeFlag(workLev,ni,nj,nk, workSet, CFInter); 
1478                                 /* new mass = current density */
1479                                 LbmFloat nbrho = QCELL(workLev,ni,nj,nk, workSet, dC);
1480                 for(int rl=1; rl< this->cDfNum ; ++rl) { nbrho += QCELL(workLev,ni,nj,nk, workSet, rl); }
1481                                 QCELL(workLev,ni,nj,nk, workSet, dMass) =  nbrho; 
1482                                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) =  1.0; 
1483
1484                                 // store point
1485                                 addToNewInterList(ni,nj,nk);
1486       }
1487       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter){
1488                                 // test, also add to list...
1489                                 addToNewInterList(ni,nj,nk);
1490                         } // NEW?
1491     }
1492
1493                 /* for symmetry, set our flag right now */
1494                 changeFlag(workLev,i,j,k, workSet, CFEmpty);
1495                 // mark cell not be changed mass... - not necessary, not in list anymore anyway!
1496         } // emptylist
1497
1498
1499         
1500         // precompute weights to get rid of order dependancies
1501         vector<lbmFloatSet> vWeights;
1502         vWeights.resize( mListFull.size() + mListEmpty.size() );
1503         int weightIndex = 0;
1504   int nbCount = 0;
1505         LbmFloat nbWeights[LBM_DFNUM];
1506         LbmFloat nbTotWeights = 0.0;
1507         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1508        iter != mListFull.end(); iter++ ) {
1509     int i=iter->x, j=iter->y, k=iter->z;
1510     nbCount = 0; nbTotWeights = 0.0;
1511     FORDF1 {
1512                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1513       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1514                                 nbCount++;
1515                                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
1516                                 else nbWeights[l] = getMassdWeight(1,i,j,k,workSet,l); // NEWSURFT
1517                                 nbTotWeights += nbWeights[l];
1518       } else {
1519                                 nbWeights[l] = -100.0; // DEBUG;
1520                         }
1521     }
1522                 if(nbCount>0) { 
1523                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
1524         vWeights[weightIndex].val[0] = nbTotWeights;
1525         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
1526         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
1527                 } else { 
1528         vWeights[weightIndex].numNbs = 0.0;
1529                 }
1530                 weightIndex++;
1531         }
1532   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1533        iter != mListEmpty.end(); iter++ ) {
1534     int i=iter->x, j=iter->y, k=iter->z;
1535     nbCount = 0; nbTotWeights = 0.0;
1536     FORDF1 {
1537                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1538       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1539                                 nbCount++;
1540                                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
1541                                 else nbWeights[l] = getMassdWeight(0,i,j,k,workSet,l); // NEWSURFT
1542                                 nbTotWeights += nbWeights[l];
1543       } else {
1544                                 nbWeights[l] = -100.0; // DEBUG;
1545                         }
1546     }
1547                 if(nbCount>0) { 
1548                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
1549         vWeights[weightIndex].val[0] = nbTotWeights;
1550         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
1551         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
1552                 } else { 
1553         vWeights[weightIndex].numNbs = 0.0;
1554                 }
1555                 weightIndex++;
1556         } 
1557         weightIndex = 0;
1558         
1559
1560         /* process full list entries, filled cells are done after this loop */
1561         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1562        iter != mListFull.end(); iter++ ) {
1563     int i=iter->x, j=iter->y, k=iter->z;
1564
1565                 LbmFloat myrho = QCELL(workLev,i,j,k, workSet, dC);
1566     FORDF1 { myrho += QCELL(workLev,i,j,k, workSet, l); } // QCELL.rho
1567
1568     LbmFloat massChange = QCELL(workLev,i,j,k, workSet, dMass) - myrho;
1569                 if(vWeights[weightIndex].numNbs>0.0) {
1570                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
1571                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
1572                         FORDF1 {
1573                                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1574         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1575                                         LbmFloat change = -1.0;
1576                                         if(nbTotWeightsp>0.0) {
1577                                                 //change = massChange * ( nbWeights[l]/nbTotWeightsp );
1578                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
1579                                         } else {
1580                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
1581                                         }
1582                                         QCELL(workLev,ni,nj,nk, workSet, dMass) += change;
1583                                 }
1584                         }
1585                         massChange = 0.0;
1586                 } else {
1587                         // Problem! no interface neighbors
1588                         mFixMass += massChange;
1589                         //TTT mNumProblems++;
1590                         //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<<mFixMass);
1591                 }
1592                 weightIndex++;
1593
1594     // already done? RFLAG(workLev,i,j,k, workSet) = CFFluid;
1595     QCELL(workLev,i,j,k, workSet, dMass) = myrho; // should be rho... but unused?
1596     QCELL(workLev,i,j,k, workSet, dFfrac) = 1.0; // should be rho... but unused?
1597   } // fulllist
1598
1599
1600         /* now, finally handle the empty cells - order is important, has to be after
1601          * full cell handling */
1602   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1603        iter != mListEmpty.end(); iter++ ) {
1604     int i=iter->x, j=iter->y, k=iter->z;
1605
1606     LbmFloat massChange = QCELL(workLev, i,j,k, workSet, dMass);
1607                 if(vWeights[weightIndex].numNbs>0.0) {
1608                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
1609                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
1610                         FORDF1 {
1611                                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1612         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1613                                         LbmFloat change = -1.0;
1614                                         if(nbTotWeightsp>0.0) {
1615                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
1616                                         } else {
1617                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
1618                                         }
1619                                         QCELL(workLev, ni,nj,nk, workSet, dMass) += change;
1620                                 }
1621                         }
1622                         massChange = 0.0;
1623                 } else {
1624                         // Problem! no interface neighbors
1625                         mFixMass += massChange;
1626                         //TTT mNumProblems++;
1627                         //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<mFixMass);
1628                 }
1629                 weightIndex++;
1630                 
1631                 // finally... make it empty 
1632     // already done? RFLAG(workLev,i,j,k, workSet) = CFEmpty;
1633     QCELL(workLev,i,j,k, workSet, dMass) = 0.0;
1634     QCELL(workLev,i,j,k, workSet, dFfrac) = 0.0;
1635         }
1636   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1637        iter != mListEmpty.end(); iter++ ) {
1638     int i=iter->x, j=iter->y, k=iter->z;
1639     changeFlag(workLev,i,j,k, otherSet, CFEmpty);
1640         } 
1641
1642
1643         // check if some of the new interface cells can be removed again 
1644         // never happens !!! not necessary
1645         // calculate ffrac for new IF cells NEW
1646
1647         // how many are really new interface cells?
1648         int numNewIf = 0;
1649   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1650        iter != mListNewInter.end(); iter++ ) {
1651     int i=iter->x, j=iter->y, k=iter->z;
1652                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
1653                         continue; 
1654                         // FIXME remove from list?
1655                 }
1656                 numNewIf++;
1657         }
1658
1659         // redistribute mass, reinit flags
1660         if(debugFlagreinit) errMsg("NEWIF", "total:"<<mListNewInter.size());
1661         float newIfFac = 1.0/(LbmFloat)numNewIf;
1662   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1663        iter != mListNewInter.end(); iter++ ) {
1664     int i=iter->x, j=iter->y, k=iter->z;
1665                 if((i<=0) || (j<=0) || 
1666                          (i>=mLevel[workLev].lSizex-1) ||
1667                          (j>=mLevel[workLev].lSizey-1) ||
1668                          ((LBMDIM==3) && ((k<=0) || (k>=mLevel[workLev].lSizez-1) ) )
1669                          ) {
1670                         continue; } // new bc, dont treat cells on boundary NEWBC
1671                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
1672                         //errMsg("???"," "<<PRINT_IJK);
1673                         continue; 
1674                 } // */
1675
1676     QCELL(workLev,i,j,k, workSet, dMass) += (mFixMass * newIfFac);
1677                 
1678                 int nbored = 0;
1679                 FORDF1 { nbored |= RFLAG_NB(workLev, i,j,k, workSet,l); }
1680                 if(!(nbored & CFBndNoslip)) { RFLAG(workLev,i,j,k, workSet) |= CFNoBndFluid; }
1681                 if(!(nbored & CFFluid))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbFluid; }
1682                 if(!(nbored & CFEmpty))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbEmpty; }
1683
1684                 if(!(RFLAG(workLev,i,j,k, otherSet)&CFInter)) {
1685                         RFLAG(workLev,i,j,k, workSet) = (CellFlagType)(RFLAG(workLev,i,j,k, workSet) | CFNoDelete);
1686                 }
1687                 if(debugFlagreinit) errMsg("NEWIF", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" f"<< convertCellFlagType2String(RFLAG(workLev,i,j,k, workSet))<<" wl"<<workLev );
1688         }
1689
1690         // reinit fill fraction
1691   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1692        iter != mListNewInter.end(); iter++ ) {
1693     int i=iter->x, j=iter->y, k=iter->z;
1694                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { continue; }
1695
1696                 initInterfaceVars(workLev, i,j,k, workSet, false); //int level, int i,int j,int k,int workSet, bool initMass) {
1697                 //LbmFloat nrho = 0.0;
1698                 //FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); }
1699     //QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho;
1700     //QCELL(workLev,i,j,k, workSet, dFlux) = FLUX_INIT;
1701         }
1702
1703         if(mListNewInter.size()>0){ 
1704                 //errMsg("FixMassDisted"," fm:"<<mFixMass<<" nif:"<<mListNewInter.size() );
1705                 mFixMass = 0.0; 
1706         }
1707
1708         // empty lists for next step
1709         mListFull.clear();
1710         mListEmpty.clear();
1711         mListNewInter.clear();
1712 } // reinitFlags
1713
1714
1715