- forgot to enable mac compile fix
[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,2004,2005 Nils Thuerey
5  *
6  * Standard LBM Factory implementation
7  *
8  *****************************************************************************/
9
10 #include "solver_class.h"
11 #include "solver_relax.h"
12
13 /*****************************************************************************/
14 /*! perform a single LBM step */
15 /*****************************************************************************/
16
17 template<class D>
18 string LbmFsgrSolver<D>::getIdString() { 
19         return string("FsgrSolver[") + D::getIdString(); 
20 }
21 template<class D>
22 void LbmFsgrSolver<D>::step() { 
23         stepMain(); 
24 }
25
26 template<class D>
27 void 
28 LbmFsgrSolver<D>::stepMain()
29 {
30 #if ELBEEM_BLENDER==1
31                 // update gui display
32                 SDL_mutexP(globalBakeLock);
33                 if(globalBakeState<0) {
34                         // this means abort... cause panic
35                         D::mPanic = 1;
36                         errMsg("LbmFsgrSolver::step","Got abort signal from GUI, causing panic, aborting...");
37                 }
38                 SDL_mutexV(globalBakeLock);
39 #endif // ELBEEM_BLENDER==1
40         D::markedClearList(); // DMC clearMarkedCellsList
41
42         // safety check, counter reset
43         D::mNumUsedCells = 0;
44         mNumInterdCells = 0;
45         mNumInvIfCells = 0;
46
47   //debugOutNnl("LbmFsgrSolver::step : "<<D::mStepCnt, 10);
48   if(!D::mSilent){ debMsgNnl("LbmFsgrSolver::step", DM_MSG, D::mName<<" cnt:"<<D::mStepCnt<<"  ", 10); }
49         //debMsgDirect(  "LbmFsgrSolver::step : "<<D::mStepCnt<<" ");
50         myTime_t timestart = getTime();
51         //myTime_t timestart = 0;
52         //if(mStartSymm) { checkSymmetry("step1"); } // DEBUG 
53
54         // important - keep for tadap
55         mCurrentMass = D::mFixMass; // reset here for next step
56         mCurrentVolume = 0.0;
57         
58         //stats
59         mMaxVlen = mMxvz = mMxvy = mMxvx = 0.0;
60
61         //change to single step advance!
62         int levsteps = 0;
63         int dsbits = D::mStepCnt ^ (D::mStepCnt-1);
64         //errMsg("S"," step:"<<D::mStepCnt<<" s-1:"<<(D::mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits));
65         for(int lev=0; lev<=mMaxRefine; lev++) {
66                 //if(! (D::mStepCnt&(1<<lev)) ) {
67                 if( dsbits & (1<<(mMaxRefine-lev)) ) {
68                         //errMsg("S"," l"<<lev);
69
70                         if(lev==mMaxRefine) {
71                                 // always advance fine level...
72                                 fineAdvance();
73                         } else {
74                                 performRefinement(lev);
75                                 performCoarsening(lev);
76                                 coarseRestrictFromFine(lev);
77                                 coarseAdvance(lev);
78                         }
79 #if FSGR_OMEGA_DEBUG==1
80                         errMsg("LbmFsgrSolver::step","LES stats l="<<lev<<" omega="<<mLevel[lev].omega<<" avgOmega="<< (mLevel[lev].avgOmega/mLevel[lev].avgOmegaCnt) );
81                         mLevel[lev].avgOmega = 0.0; mLevel[lev].avgOmegaCnt = 0.0;
82 #endif // FSGR_OMEGA_DEBUG==1
83                         levsteps++;
84                 }
85                 mCurrentMass   += mLevel[lev].lmass;
86                 mCurrentVolume += mLevel[lev].lvolume;
87         }
88
89   // prepare next step
90         D::mStepCnt++;
91
92
93         // some dbugging output follows
94         // calculate MLSUPS
95         myTime_t timeend = getTime();
96
97         D::mNumUsedCells += mNumInterdCells; // count both types for MLSUPS
98         mAvgNumUsedCells += D::mNumUsedCells;
99         D::mMLSUPS = (D::mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000);
100         if(D::mMLSUPS>10000){ D::mMLSUPS = -1; }
101         else { mAvgMLSUPS += D::mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
102         
103         LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1())) / ((timeend-timestart)/(double)1000.0) ) / (1000000);
104         if(totMLSUPS>10000) totMLSUPS = -1;
105         mNumInvIfTotal += mNumInvIfCells; // debug
106
107   // do some formatting 
108   if(!D::mSilent){ 
109                 string sepStr(""); // DEBUG
110 #ifndef USE_MSVC6FIXES
111                 int avgcls = (int)(mAvgNumUsedCells/(long long int)D::mStepCnt);
112 #else
113                 int avgcls = (int)(mAvgNumUsedCells/(_int64)D::mStepCnt);
114 #endif
115                 debMsgDirect( 
116                         "mlsups(curr:"<<D::mMLSUPS<<
117                         " avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< sepStr<<
118                         " totcls:"<<(D::mNumUsedCells)<< sepStr<<
119                         " avgcls:"<< avgcls<< sepStr<<
120                         " intd:"<<mNumInterdCells<< sepStr<<
121                         " invif:"<<mNumInvIfCells<< sepStr<<
122                         " invift:"<<mNumInvIfTotal<< sepStr<<
123                         " fsgrcs:"<<mNumFsgrChanges<< sepStr<<
124                         " filled:"<<D::mNumFilledCells<<", emptied:"<<D::mNumEmptiedCells<< sepStr<<
125                         " mMxv:"<<mMxvx<<","<<mMxvy<<","<<mMxvz<<", tscnts:"<<mTimeSwitchCounts<< sepStr<<
126                         " probs:"<<mNumProblems<< sepStr<<
127                         " simt:"<<mSimulationTime<< sepStr<<
128                         " for '"<<D::mName<<"' " );
129
130                 debMsgDirect(std::endl);
131                 debMsgDirect(D::mStepCnt<<": dccd="<< mCurrentMass<<"/"<<mCurrentVolume<<"(fix="<<D::mFixMass<<",ini="<<mInitialMass<<") ");
132                 debMsgDirect(std::endl);
133
134                 // nicer output
135                 debMsgDirect(std::endl); // 
136                 //debMsgStd(" ",DM_MSG," ",10);
137         } else {
138                 debMsgDirect(".");
139                 //if((mStepCnt%10)==9) debMsgDirect("\n");
140         }
141
142         if(D::mStepCnt==1) {
143                 mMinNoCells = mMaxNoCells = D::mNumUsedCells;
144         } else {
145                 if(D::mNumUsedCells>mMaxNoCells) mMaxNoCells = D::mNumUsedCells;
146                 if(D::mNumUsedCells<mMinNoCells) mMinNoCells = D::mNumUsedCells;
147         }
148         
149         // mass scale test
150         if((mMaxRefine>0)&&(mInitialMass>0.0)) {
151                 LbmFloat mscale = mInitialMass/mCurrentMass;
152
153                 mscale = 1.0;
154                 const LbmFloat dchh = 0.001;
155                 if(mCurrentMass<mInitialMass) mscale = 1.0+dchh;
156                 if(mCurrentMass>mInitialMass) mscale = 1.0-dchh;
157
158                 // use mass rescaling?
159                 // with float precision this seems to be nonsense...
160                 const bool MREnable = false;
161
162                 const int MSInter = 2;
163                 static int mscount = 0;
164                 if( (MREnable) && ((mLevel[0].lsteps%MSInter)== (MSInter-1)) && ( ABS( (mInitialMass/mCurrentMass)-1.0 ) > 0.01) && ( dsbits & (1<<(mMaxRefine-0)) ) ){
165                         // example: FORCE RESCALE MASS! ini:1843.5, cur:1817.6, f=1.01425 step:22153 levstep:5539 msc:37
166                         // mass rescale MASS RESCALE check
167                         errMsg("MDTDD","\n\n");
168                         errMsg("MDTDD","FORCE RESCALE MASS! "
169                                         <<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass)
170                                         <<" step:"<<D::mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
171                                         );
172                         errMsg("MDTDD","\n\n");
173
174                         mscount++;
175                         for(int lev=mMaxRefine; lev>=0 ; lev--) {
176                                 //for(int workSet = 0; workSet<=1; workSet++) {
177                                 int wss = 0;
178                                 int wse = 1;
179 #if COMPRESSGRIDS==1
180                                 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
181 #endif // COMPRESSGRIDS==1
182                                 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
183
184                                         FSGR_FORIJK1(lev) {
185                                                 if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) 
186                                                         ) {
187
188                                                         FORDF0 { QCELL(lev, i,j,k,workSet, l) *= mscale; }
189                                                         QCELL(lev, i,j,k,workSet, dMass) *= mscale;
190                                                         QCELL(lev, i,j,k,workSet, dFfrac) *= mscale;
191
192                                                 } else {
193                                                         continue;
194                                                 }
195                                         }
196                                 }
197                                 mLevel[lev].lmass *= mscale;
198                         }
199                 } 
200
201                 mCurrentMass *= mscale;
202         }// if mass scale test */
203         else {
204                 // use current mass after full step for initial setting
205                 if((mMaxRefine>0)&&(mInitialMass<=0.0) && (levsteps == (mMaxRefine+1))) {
206                         mInitialMass = mCurrentMass;
207                         debMsgStd("MDTDD",DM_NOTIFY,"Second Initial Mass Init: "<<mInitialMass, 2);
208                 }
209         }
210
211         // one of the last things to do - adapt timestep
212         // was in fineAdvance before... 
213         if(mTimeAdap) {
214                 adaptTimestep();
215         } // time adaptivity
216
217         // debug - raw dump of ffrac values
218         /*if((D::mStepCnt%100)==1){
219                 std::ostringstream name;
220                 name <<"fill_" << D::mStepCnt <<".dump";
221                 FILE *file = fopen(name.str().c_str(),"w");
222                 for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
223                  for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
224                         for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
225                                 float val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
226                                 //fwrite( &val, sizeof(val), 1, file); // binary
227                                 fprintf(file, "%f ",val); // text
228                                 //errMsg("W", PRINT_IJK<<" val:"<<val);
229                         }
230                         fprintf(file, "\n"); // text
231                  }
232                         fprintf(file, "\n"); // text
233                 }
234                 fclose(file);
235         } // */
236
237 #if ELBEEM_BLENDER!=1
238         if(mUseTestdata) handleTestdata();
239 #endif // ELBEEM_BLENDER!=1
240 }
241
242 template<class D>
243 void 
244 LbmFsgrSolver<D>::fineAdvance()
245 {
246         // do the real thing...
247         mainLoop( mMaxRefine );
248         if(mUpdateFVHeight) {
249                 // warning assume -Y gravity...
250                 mFVHeight = mCurrentMass*mFVArea/((LbmFloat)(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizez));
251                 if(mFVHeight<1.0) mFVHeight = 1.0;
252                 D::mpParam->setFluidVolumeHeight(mFVHeight);
253         }
254
255         // advance time before timestep change
256         mSimulationTime += D::mpParam->getStepTime();
257         // time adaptivity
258         D::mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) );
259         //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG 
260         if(!D::mSilent){ errMsg("fineAdvance"," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps) ); }
261
262         // update other set
263   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
264   mLevel[mMaxRefine].setCurr   ^= 1;
265   mLevel[mMaxRefine].lsteps++;
266
267         // flag init... (work on current set, to simplify flag checks)
268         reinitFlags( mLevel[mMaxRefine].setCurr );
269         if(!D::mSilent){ errMsg("fineAdvance"," flags reinit on set "<< mLevel[mMaxRefine].setCurr ); }
270 }
271
272
273 /*****************************************************************************/
274 //! coarse/fine step functions
275 /*****************************************************************************/
276
277 // access to own dfs during step (may be changed to local array)
278 #define MYDF(l) RAC(ccel, l)
279
280 template<class D>
281 void 
282 LbmFsgrSolver<D>::mainLoop(int lev)
283 {
284         // loops over _only inner_ cells  -----------------------------------------------------------------------------------
285         LbmFloat calcCurrentMass = 0.0;
286         LbmFloat calcCurrentVolume = 0.0;
287         int      calcCellsFilled = D::mNumFilledCells;
288         int      calcCellsEmptied = D::mNumEmptiedCells;
289         int      calcNumUsedCells = D::mNumUsedCells;
290
291 #if PARALLEL==1
292 #include "paraloop.h"
293 #else // PARALLEL==1
294   { // main loop region
295         int kstart=D::getForZMin1(), kend=D::getForZMax1();
296 #define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz);
297 #endif // PARALLEL==1
298
299
300         // local to loop
301         CellFlagType nbflag[LBM_DFNUM]; 
302         LbmFloat *ccel = NULL;
303         LbmFloat *tcel = NULL;
304         int oldFlag;
305         int newFlag;
306         int nbored;
307         LbmFloat m[LBM_DFNUM];
308         LbmFloat rho, ux, uy, uz, tmp, usqr;
309         LbmFloat mass, change;
310         usqr = tmp = 0.0; 
311 #if OPT3D==1 
312         LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM], lcsmomega;
313 #endif // OPT3D==true 
314
315         
316         // ifempty cell conversion flags
317         bool iffilled, ifemptied;
318         LbmFloat nbfracs[LBM_DFNUM]; // ffracs of neighbors
319         int recons[LBM_DFNUM];   // reconstruct this DF?
320         int numRecons;           // how many are reconstructed?
321
322         // slow surf regions smooth (if below)
323         const LbmFloat smoothStrength = 0.0; //0.01;
324         const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03;
325         const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit;
326         
327         CellFlagType *pFlagSrc;
328         CellFlagType *pFlagDst;
329         pFlagSrc = &RFLAG(lev, 0,1, kstart,SRCS(lev)); // omp
330         pFlagDst = &RFLAG(lev, 0,1, kstart,TSET(lev)); // omp
331         ccel = RACPNT(lev, 0,1, kstart ,SRCS(lev)); // omp
332         tcel = RACPNT(lev, 0,1, kstart ,TSET(lev)); // omp
333         //CellFlagType *pFlagTar = NULL;
334         int pFlagTarOff;
335         if(mLevel[lev].setOther==1) pFlagTarOff = mLevel[lev].lOffsz;
336         else pFlagTarOff = -mLevel[lev].lOffsz;
337 #define ADVANCE_POINTERS(p)     \
338         ccel += (QCELLSTEP*(p));        \
339         tcel += (QCELLSTEP*(p));        \
340         pFlagSrc+= (p); \
341         pFlagDst+= (p); \
342         i+= (p);
343
344
345         // ---
346         // now stream etc.
347
348         // use template functions for 2D/3D
349 #if COMPRESSGRIDS==0
350   for(int k=kstart;k<kend;++k) {
351   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
352   for(int i=0;i<mLevel[lev].lSizex-2;   ) {
353 #else // COMPRESSGRIDS==0
354         int kdir = 1; // COMPRT ON
355         if(mLevel[mMaxRefine].setCurr==1) {
356                 kdir = -1;
357                 int temp = kend;
358                 kend = kstart-1;
359                 kstart = temp-1;
360         } // COMPRT
361
362 #if PARALLEL==0
363   const int id = 0, Nthrds = 1;
364 #endif // PARALLEL==1
365   const int Nj = mLevel[mMaxRefine].lSizey;
366         int jstart = 0+( id * (Nj / Nthrds) );
367         int jend   = 0+( (id+1) * (Nj / Nthrds) );
368   if( ((Nj/Nthrds) *Nthrds) != Nj) {
369     errMsg("LbmFsgrSolver","Invalid domain size Nj="<<Nj<<" Nthrds="<<Nthrds);
370   }
371         // cutoff obstacle boundary
372         if(jstart<1) jstart = 1;
373         if(jend>mLevel[mMaxRefine].lSizey-1) jend = mLevel[mMaxRefine].lSizey-1;
374
375 #if PARALLEL==1
376         errMsg("LbmFsgrSolver::mainLoop","id="<<id<<" js="<<jstart<<" je="<<jend<<" jdir="<<(1) ); // debug
377 #endif // PARALLEL==1
378   for(int k=kstart;k!=kend;k+=kdir) {
379
380         //errMsg("LbmFsgrSolver::mainLoop","k="<<k<<" ks="<<kstart<<" ke="<<kend<<" kdir="<<kdir<<" x*y="<<mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum ); // debug
381   pFlagSrc = &RFLAG(lev, 0, jstart, k, SRCS(lev)); // omp test // COMPRT ON
382   pFlagDst = &RFLAG(lev, 0, jstart, k, TSET(lev)); // omp test
383   ccel = RACPNT(lev,     0, jstart, k, SRCS(lev)); // omp test
384   tcel = RACPNT(lev,     0, jstart, k, TSET(lev)); // omp test // COMPRT ON
385
386   //for(int j=1;j<mLevel[lev].lSizey-1;++j) {
387   for(int j=jstart;j!=jend;++j) {
388   for(int i=0;i<mLevel[lev].lSizex-2;   ) {
389 #endif // COMPRESSGRIDS==0
390
391                 ADVANCE_POINTERS(1); 
392 #if FSGR_STRICT_DEBUG==1
393                 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
394                 if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) || 
395                     (&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) {
396                         errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<<
397                                         RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<<
398                                         (*pFlagSrc)<<","<<(*pFlagDst) <<",  pointers "<<
399           (int)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(int)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
400           (int)(pFlagSrc)<<","<<(int)(pFlagDst)<<" "
401                                         ); 
402                         D::mPanic=1;
403                 }       
404                 if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) || 
405                     (&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) {
406                         errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<<
407           (int)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(int)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
408           (int)(ccel)<<","<<(int)(tcel)<<" "
409                                         ); 
410                         D::mPanic=1;
411                 }       
412 #endif
413                 oldFlag = *pFlagSrc;
414                 // stream from current set to other, then collide and store
415                 
416                 // old INTCFCOARSETEST==1
417                 if( (oldFlag & (CFGrFromCoarse)) ) {  // interpolateFineFromCoarse test!
418                         if(( D::mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
419                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
420                         } else {
421                                 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
422                                 calcNumUsedCells++;
423                         }
424                         continue; // interpolateFineFromCoarse test!
425                 } // interpolateFineFromCoarse test!  old INTCFCOARSETEST==1
426         
427                 if(oldFlag & (CFMbndInflow)) {
428                         // fluid & if are ok, fill if later on
429                         int isValid = oldFlag & (CFFluid|CFInter);
430                         const LbmFloat iniRho = 1.0;
431                         const int OId = oldFlag>>24;
432                         if(!isValid) {
433                                 // make new if cell
434                                 const LbmVec vel(mObjectSpeeds[OId]);
435                                 // TODO add OPT3D treatment
436                                 FORDF0 { RAC(tcel, l) = D::getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
437                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
438                                 RAC(tcel, dFlux) = FLUX_INIT;
439                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
440                                 calcCurrentMass += iniRho; calcCurrentVolume += 1.0; calcNumUsedCells++;
441                                 mInitialMass += iniRho;
442                                 // dont treat cell until next step
443                                 continue;
444                         } 
445                 } 
446                 else  // these are exclusive
447                 if(oldFlag & (CFMbndOutflow)) {
448                         int isnotValid = oldFlag & (CFFluid);
449                         if(isnotValid) {
450                                 // remove fluid cells, shouldnt be here anyway
451                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
452                                 mInitialMass -= fluidRho;
453                                 const LbmFloat iniRho = 0.0;
454                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
455                                 RAC(tcel, dFlux) = FLUX_INIT;
456                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
457
458                                 // same as ifemptied for if below
459                                 LbmPoint emptyp;
460                                 emptyp.x = i; emptyp.y = j; emptyp.z = k;
461 #if PARALLEL==1
462                                 calcListEmpty[id].push_back( emptyp );
463 #else // PARALLEL==1
464                                 mListEmpty.push_back( emptyp );
465 #endif // PARALLEL==1
466                                 calcCellsEmptied++;
467                                 continue;
468                         }
469                 }
470
471                 if(oldFlag & (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)) { 
472                         *pFlagDst = oldFlag;
473                         continue;
474                 }
475                 /*if( oldFlag & CFNoBndFluid ) {  // TEST ME FASTER?
476                         OPTIMIZED_STREAMCOLLIDE; PERFORM_USQRMAXCHECK;
477                         RAC(tcel,dFfrac) = 1.0; 
478                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
479                         calcCurrentMass += rho; calcCurrentVolume += 1.0;
480                         calcNumUsedCells++;
481                         continue;
482                 }// TEST ME FASTER? */
483
484                 // only neighbor flags! not own flag
485                 nbored = 0;
486                 
487 #if OPT3D==0
488                 FORDF1 {
489                         nbflag[l] = RFLAG_NB(lev, i,j,k,SRCS(lev),l);
490                         nbored |= nbflag[l];
491                 } 
492 #else
493                 nbflag[dSB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dSB];
494                 nbflag[dWB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWB];
495                 nbflag[ dB] = *(pFlagSrc + (-mLevel[lev].lOffsy)); nbored |= nbflag[dB];
496                 nbflag[dEB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dEB];
497                 nbflag[dNB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNB];
498
499                 nbflag[dSW] = *(pFlagSrc + (-mLevel[lev].lOffsx+-1)); nbored |= nbflag[dSW];
500                 nbflag[ dS] = *(pFlagSrc + (-mLevel[lev].lOffsx)); nbored |= nbflag[dS];
501                 nbflag[dSE] = *(pFlagSrc + (-mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dSE];
502
503                 nbflag[ dW] = *(pFlagSrc + (-1)); nbored |= nbflag[dW];
504                 nbflag[ dE] = *(pFlagSrc + ( 1)); nbored |= nbflag[dE];
505
506                 nbflag[dNW] = *(pFlagSrc + ( mLevel[lev].lOffsx+-1)); nbored |= nbflag[dNW];
507           nbflag[ dN] = *(pFlagSrc + ( mLevel[lev].lOffsx)); nbored |= nbflag[dN];
508                 nbflag[dNE] = *(pFlagSrc + ( mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dNE];
509
510                 nbflag[dST] = *(pFlagSrc + ( mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dST];
511                 nbflag[dWT] = *(pFlagSrc + ( mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWT];
512                 nbflag[ dT] = *(pFlagSrc + ( mLevel[lev].lOffsy)); nbored |= nbflag[dT];
513                 nbflag[dET] = *(pFlagSrc + ( mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dET];
514                 nbflag[dNT] = *(pFlagSrc + ( mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNT];
515                 // */
516 #endif
517
518                 // pointer to destination cell
519                 calcNumUsedCells++;
520
521                 // FLUID cells 
522                 if( oldFlag & CFFluid ) { 
523                         // only standard fluid cells (with nothing except fluid as nbs
524
525                         if(oldFlag&CFMbndInflow) {
526                                 // force velocity for inflow
527                                 const int OId = oldFlag>>24;
528                                 DEFAULT_STREAM;
529                                 //const LbmFloat fluidRho = 1.0;
530                                 // for submerged inflows, streaming would have to be performed...
531                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
532                                 const LbmVec vel(mObjectSpeeds[OId]);
533                                 ux=vel[0], uy=vel[1], uz=vel[2]; 
534                                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
535                                 FORDF0 { RAC(tcel, l) = D::getCollideEq(l, fluidRho,ux,uy,uz); }
536                                 rho = fluidRho;
537                         } else {
538                                 if(nbored&CFBnd) {
539                                         DEFAULT_STREAM;
540                                         ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
541                                         DEFAULT_COLLIDE;
542                                         oldFlag &= (~CFNoBndFluid);
543                                 } else {
544                                         // do standard stream/collide
545                                         OPTIMIZED_STREAMCOLLIDE;
546                                         // FIXME check for which cells this is executed!
547                                         oldFlag |= CFNoBndFluid;
548                                 } 
549                         }
550
551                         PERFORM_USQRMAXCHECK;
552                         // "normal" fluid cells
553                         RAC(tcel,dFfrac) = 1.0; 
554                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
555                         calcCurrentMass += rho; 
556                         calcCurrentVolume += 1.0;
557                         continue;
558                 }
559                 
560                 newFlag  = oldFlag;
561                 // make sure here: always check which flags to really unset...!
562                 newFlag = newFlag & (~( 
563                                         CFNoNbFluid|CFNoNbEmpty| CFNoDelete 
564                                         | CFNoInterpolSrc
565                                         | CFNoBndFluid
566                                         ));
567                 if(!(nbored&CFBnd)) {
568                         newFlag |= CFNoBndFluid;
569                 }
570
571                 // store own dfs and mass
572                 mass = RAC(ccel,dMass);
573
574                 // WARNING - only interface cells arrive here!
575                 // read distribution funtions of adjacent cells = stream step
576                 DEFAULT_STREAM;
577
578                 if((nbored & CFFluid)==0) { newFlag |= CFNoNbFluid; mNumInvIfCells++; }
579                 if((nbored & CFEmpty)==0) { newFlag |= CFNoNbEmpty; mNumInvIfCells++; }
580
581                 // calculate mass exchange for interface cells 
582                 LbmFloat myfrac = RAC(ccel,dFfrac);
583 #               define nbdf(l) m[ D::dfInv[(l)] ]
584
585                 // update mass 
586                 // only do boundaries for fluid cells, and interface cells without
587                 // any fluid neighbors (assume that interface cells _with_ fluid
588                 // neighbors are affected enough by these) 
589                 // which Df's have to be reconstructed? 
590                 // for fluid cells - just the f_i difference from streaming to empty cells  ----
591                 numRecons = 0;
592
593                 FORDF1 { // dfl loop
594                         recons[l] = 0;
595                         nbfracs[l] = 0.0;
596                         // finally, "normal" interface cells ----
597                         if( nbflag[l]&(CFFluid|CFBnd) ) { // NEWTEST! FIXME check!!!!!!!!!!!!!!!!!!
598                                 change = nbdf(l) - MYDF(l);
599                         }
600                         // interface cells - distuingish cells that shouldn't fill/empty 
601                         else if( nbflag[l] & CFInter ) {
602                                 
603                                 LbmFloat mynbfac,nbnbfac;
604                                 // NEW TEST t1
605                                 // t2 -> off
606                                 if((oldFlag&CFNoBndFluid)&&(nbflag[l]&CFNoBndFluid)) {
607                                         mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
608                                         nbnbfac = 1.0/mynbfac;
609                                 } else {
610                                         mynbfac = nbnbfac = 1.0; // switch calc flux off
611                                 }
612                                 //mynbfac = nbnbfac = 1.0; // switch calc flux off t3
613
614                                 // perform interface case handling
615                                 if ((oldFlag|nbflag[l])&(CFNoNbFluid|CFNoNbEmpty)) {
616                                 switch (oldFlag&(CFNoNbFluid|CFNoNbEmpty)) {
617                                         case 0: 
618                                                 // we are a normal cell so... 
619                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
620                                                         case CFNoNbFluid: 
621                                                                 // just fill current cell = empty neighbor 
622                                                                 change = nbnbfac*nbdf(l) ; goto changeDone; 
623                                                         case CFNoNbEmpty: 
624                                                                 // just empty current cell = fill neighbor 
625                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
626                                                 }
627                                                 break;
628
629                                         case CFNoNbFluid: 
630                                                 // we dont have fluid nb's so...
631                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
632                                                         case 0: 
633                                                         case CFNoNbEmpty: 
634                                                                 // we have no fluid nb's -> just empty
635                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
636                                                 }
637                                                 break;
638
639                                         case CFNoNbEmpty: 
640                                                 // we dont have empty nb's so...
641                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
642                                                         case 0: 
643                                                         case CFNoNbFluid: 
644                                                                 // we have no empty nb's -> just fill
645                                                                 change = nbnbfac*nbdf(l); goto changeDone; 
646                                                 }
647                                                 break;
648                                 }} // inter-inter exchange
649
650                                 // just do normal mass exchange...
651                                 change = ( nbnbfac*nbdf(l) - mynbfac*MYDF(l) ) ;
652                         changeDone: ;
653                                 nbfracs[l] = QCELL_NB(lev, i,j,k, SRCS(lev),l, dFfrac);
654                                 change *=  (myfrac + nbfracs[l]) * 0.5;
655                         } // the other cell is interface
656
657                         // last alternative - reconstruction in this direction
658                         else {
659                                 // empty + bnd case
660                                 recons[l] = 1; 
661                                 numRecons++;
662                                 change = 0.0; 
663                         }
664
665                         // modify mass at SRCS
666                         mass += change;
667                 } // l
668                 // normal interface, no if empty/fluid
669
670                 LbmFloat nv1,nv2;
671                 LbmFloat nx,ny,nz;
672
673                 if(nbflag[dE] &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } else nv1 = 0.0;
674                 if(nbflag[dW] &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } else nv2 = 0.0;
675                 nx = 0.5* (nv2-nv1);
676                 if(nbflag[dN] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
677                 if(nbflag[dS] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
678                 ny = 0.5* (nv2-nv1);
679 #if LBMDIM==3
680                 if(nbflag[dT] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
681                 if(nbflag[dB] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
682                 nz = 0.5* (nv2-nv1);
683 #else // LBMDIM==3
684                 nz = 0.0;
685 #endif // LBMDIM==3
686
687                 if( (ABS(nx)+ABS(ny)+ABS(nz)) > LBM_EPSILON) {
688                         // normal ok and usable...
689                         FORDF1 {
690                                 if( (D::dfDvecX[l]*nx + D::dfDvecY[l]*ny + D::dfDvecZ[l]*nz)  // dot Dvec,norml
691                                                 > LBM_EPSILON) {
692                                         recons[l] = 2; 
693                                         numRecons++;
694                                 } 
695                         }
696                 }
697
698                 // calculate macroscopic cell values
699                 LbmFloat oldUx, oldUy, oldUz;
700                 LbmFloat oldRho; // OLD rho = ccel->rho;
701 #if OPT3D==0
702                         oldRho=RAC(ccel,0);
703                         oldUx = oldUy = oldUz = 0.0;
704                         for(int l=1; l<D::cDfNum; l++) {
705                                 oldRho += RAC(ccel,l);
706                                 oldUx  += (D::dfDvecX[l]*RAC(ccel,l));
707                                 oldUy  += (D::dfDvecY[l]*RAC(ccel,l)); 
708                                 oldUz  += (D::dfDvecZ[l]*RAC(ccel,l)); 
709                         } 
710 #else // OPT3D==0
711                 oldRho = + RAC(ccel,dC)  + RAC(ccel,dN )
712                                 + RAC(ccel,dS ) + RAC(ccel,dE )
713                                 + RAC(ccel,dW ) + RAC(ccel,dT )
714                                 + RAC(ccel,dB ) + RAC(ccel,dNE)
715                                 + RAC(ccel,dNW) + RAC(ccel,dSE)
716                                 + RAC(ccel,dSW) + RAC(ccel,dNT)
717                                 + RAC(ccel,dNB) + RAC(ccel,dST)
718                                 + RAC(ccel,dSB) + RAC(ccel,dET)
719                                 + RAC(ccel,dEB) + RAC(ccel,dWT)
720                                 + RAC(ccel,dWB);
721
722                         oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
723                                 + RAC(ccel,dNE) - RAC(ccel,dNW)
724                                 + RAC(ccel,dSE) - RAC(ccel,dSW)
725                                 + RAC(ccel,dET) + RAC(ccel,dEB)
726                                 - RAC(ccel,dWT) - RAC(ccel,dWB);
727
728                         oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
729                                 + RAC(ccel,dNE) + RAC(ccel,dNW)
730                                 - RAC(ccel,dSE) - RAC(ccel,dSW)
731                                 + RAC(ccel,dNT) + RAC(ccel,dNB)
732                                 - RAC(ccel,dST) - RAC(ccel,dSB);
733
734                         oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
735                                 + RAC(ccel,dNT) - RAC(ccel,dNB)
736                                 + RAC(ccel,dST) - RAC(ccel,dSB)
737                                 + RAC(ccel,dET) - RAC(ccel,dEB)
738                                 + RAC(ccel,dWT) - RAC(ccel,dWB);
739 #endif
740
741                 // now reconstruction
742 #define REFERENCE_PRESSURE 1.0 // always atmosphere...
743 #if OPT3D==0
744                 // NOW - construct dist funcs from empty cells
745                 FORDF1 {
746                         if(recons[ l ]) {
747                                 m[ D::dfInv[l] ] = 
748                                         D::getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 
749                                         D::getCollideEq(D::dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 
750                                         - MYDF( l );
751                         }
752                 }
753 #else
754                 ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
755                 rho = REFERENCE_PRESSURE;
756                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
757                 if(recons[dN ]) { m[dS ] = EQN  + EQS  - MYDF(dN ); }
758                 if(recons[dS ]) { m[dN ] = EQS  + EQN  - MYDF(dS ); }
759                 if(recons[dE ]) { m[dW ] = EQE  + EQW  - MYDF(dE ); }
760                 if(recons[dW ]) { m[dE ] = EQW  + EQE  - MYDF(dW ); }
761                 if(recons[dT ]) { m[dB ] = EQT  + EQB  - MYDF(dT ); }
762                 if(recons[dB ]) { m[dT ] = EQB  + EQT  - MYDF(dB ); }
763                 if(recons[dNE]) { m[dSW] = EQNE + EQSW - MYDF(dNE); }
764                 if(recons[dNW]) { m[dSE] = EQNW + EQSE - MYDF(dNW); }
765                 if(recons[dSE]) { m[dNW] = EQSE + EQNW - MYDF(dSE); }
766                 if(recons[dSW]) { m[dNE] = EQSW + EQNE - MYDF(dSW); }
767                 if(recons[dNT]) { m[dSB] = EQNT + EQSB - MYDF(dNT); }
768                 if(recons[dNB]) { m[dST] = EQNB + EQST - MYDF(dNB); }
769                 if(recons[dST]) { m[dNB] = EQST + EQNB - MYDF(dST); }
770                 if(recons[dSB]) { m[dNT] = EQSB + EQNT - MYDF(dSB); }
771                 if(recons[dET]) { m[dWB] = EQET + EQWB - MYDF(dET); }
772                 if(recons[dEB]) { m[dWT] = EQEB + EQWT - MYDF(dEB); }
773                 if(recons[dWT]) { m[dEB] = EQWT + EQEB - MYDF(dWT); }
774                 if(recons[dWB]) { m[dET] = EQWB + EQET - MYDF(dWB); }
775 #endif          
776
777                 // mass streaming done... do normal collide
778                 ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2];
779                 DEFAULT_COLLIDE;
780                 PERFORM_USQRMAXCHECK;
781                 // rho init from default collide necessary for fill/empty check below
782
783                 // inflow bc handling
784                 if(oldFlag & (CFMbndInflow)) {
785                         // fill if cells in inflow region
786                         if(myfrac<0.5) { 
787                                 mass += 0.25; 
788                                 mInitialMass += 0.25;
789                         }
790                 } 
791
792                 // interface cell filled or emptied?
793                 iffilled = ifemptied = 0;
794                 // interface cells empty/full?, WARNING: to mark these cells, better do it at the end of reinitCellFlags
795                 // interface cell if full?
796                 if( (mass) >= (rho * (1.0+FSGR_MAGICNR)) ) { iffilled = 1; }
797                 // interface cell if empty?
798                 if( (mass) <= (rho * (   -FSGR_MAGICNR)) ) { ifemptied = 1; }
799
800                 if(oldFlag & (CFMbndOutflow)) {
801                         mInitialMass -= mass;
802                         mass = myfrac = 0.0;
803                         iffilled = 0; ifemptied = 1;
804                 }
805
806                 // looks much nicer... LISTTRICK
807 #if FSGR_LISTTRICK==1
808                 if(newFlag&CFNoBndFluid) { // test NEW TEST
809                         if(!iffilled) {
810                                 // remove cells independent from amount of change...
811                                 if( (oldFlag & CFNoNbEmpty)&&(newFlag & CFNoNbEmpty)&&
812                                                 ( (mass>(rho*FSGR_LISTTTHRESHFULL))  || ((nbored&CFInter)==0)  )
813                                         ) { 
814                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","filled "<<PRINT_IJK); };
815                                         iffilled = 1; 
816                                 } 
817                         }
818                         if(!ifemptied) {
819                                 if( (oldFlag & CFNoNbFluid)&&(newFlag & CFNoNbFluid)&&
820                                                 ( (mass<(rho*FSGR_LISTTTHRESHEMPTY)) || ((nbored&CFInter)==0)  )
821                                                 ) 
822                                 { 
823                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","emptied "<<PRINT_IJK); };
824                                         ifemptied = 1; 
825                                 } 
826                         }
827                 } // nobndfluid only */
828 #endif
829
830                 //iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!!
831                 
832
833                 // now that all dfs are known, handle last changes
834                 if(iffilled) {
835                         LbmPoint filledp;
836                         filledp.x = i; filledp.y = j; filledp.z = k;
837 #if PARALLEL==1
838                         calcListFull[id].push_back( filledp );
839 #else // PARALLEL==1
840                         mListFull.push_back( filledp );
841 #endif // PARALLEL==1
842                         //D::mNumFilledCells++; // DEBUG
843                         calcCellsFilled++;
844                 }
845                 else if(ifemptied) {
846                         LbmPoint emptyp;
847                         emptyp.x = i; emptyp.y = j; emptyp.z = k;
848 #if PARALLEL==1
849                         calcListEmpty[id].push_back( emptyp );
850 #else // PARALLEL==1
851                         mListEmpty.push_back( emptyp );
852 #endif // PARALLEL==1
853                         //D::mNumEmptiedCells++; // DEBUG
854                         calcCellsEmptied++;
855                 } else {
856                         // ...
857                 }
858                 
859                 // dont cutoff values -> better cell conversions
860                 RAC(tcel,dFfrac)   = (mass/rho);
861
862                 // init new flux value
863                 float flux = FLUX_INIT; // dxqn on
864                 if(newFlag&CFNoBndFluid) {
865                         //flux = 50.0; // extreme on
866                         for(int nn=1; nn<D::cDfNum; nn++) { 
867                                 if(nbflag[nn] & (CFFluid|CFInter|CFBnd)) { flux += D::dfLength[nn]; }
868                         }
869                         // optical hack - smooth slow moving
870                         // surface regions
871                         if(usqr< sssUsqrLimit) {
872                         for(int nn=1; nn<D::cDfNum; nn++) { 
873                                 if(nbfracs[nn]!=0.0) {
874                                         LbmFloat curSmooth = (sssUsqrLimit-usqr)*sssUsqrLimitInv;
875                                         if(curSmooth>1.0) curSmooth=1.0;
876                                         flux *= (1.0+ smoothStrength*curSmooth * (nbfracs[nn]-myfrac)) ;
877                                 }
878                         } }
879                         // NEW TEST */
880                 }
881                 // flux = FLUX_INIT; // calc flux off
882                 QCELL(lev, i,j,k,TSET(lev), dFlux) = flux; // */
883
884                 // perform mass exchange with streamed values
885                 QCELL(lev, i,j,k,TSET(lev), dMass) = mass; // MASST
886                 // set new flag 
887                 *pFlagDst = (CellFlagType)newFlag;
888                 calcCurrentMass += mass; 
889                 calcCurrentVolume += RAC(tcel,dFfrac);
890
891                 // interface cell handling done...
892         } // i
893         int i=0; //dummy
894         ADVANCE_POINTERS(2);
895         } // j
896
897 #if COMPRESSGRIDS==1
898 #if PARALLEL==1
899         //frintf(stderr," (id=%d k=%d) ",id,k);
900 # pragma omp barrier
901 #endif // PARALLEL==1
902 #else // COMPRESSGRIDS==1
903         int i=0; //dummy
904         ADVANCE_POINTERS(mLevel[lev].lSizex*2);
905 #endif // COMPRESSGRIDS==1
906   } // all cell loop k,j,i
907
908         } // main loop region
909
910         // write vars from parallel computations to class
911         //errMsg("DFINI"," maxr l"<<mMaxRefine<<" cm="<<calcCurrentMass<<" cv="<<calcCurrentVolume );
912         mLevel[lev].lmass    = calcCurrentMass;
913         mLevel[lev].lvolume  = calcCurrentVolume;
914         //mCurrentMass   += calcCurrentMass;
915         //mCurrentVolume += calcCurrentVolume;
916         D::mNumFilledCells  = calcCellsFilled;
917         D::mNumEmptiedCells = calcCellsEmptied;
918         D::mNumUsedCells = calcNumUsedCells;
919 #if PARALLEL==1
920         //errMsg("PARALLELusqrcheck"," curr: "<<mMaxVlen<<"|"<<mMxvx<<","<<mMxvy<<","<<mMxvz);
921         for(int i=0; i<MAX_THREADS; i++) {
922                 for(int j=0; j<calcListFull[i].size() ; j++) mListFull.push_back( calcListFull[i][j] );
923                 for(int j=0; j<calcListEmpty[i].size(); j++) mListEmpty.push_back( calcListEmpty[i][j] );
924                 if(calcMaxVlen[i]>mMaxVlen) { 
925                         mMxvx = calcMxvx[i]; 
926                         mMxvy = calcMxvy[i]; 
927                         mMxvz = calcMxvz[i]; 
928                         mMaxVlen = calcMaxVlen[i]; 
929                 } 
930                 errMsg("PARALLELusqrcheck"," curr: "<<mMaxVlen<<"|"<<mMxvx<<","<<mMxvy<<","<<mMxvz<<
931                                 "      calc["<<i<<": "<<calcMaxVlen[i]<<"|"<<calcMxvx[i]<<","<<calcMxvy[i]<<","<<calcMxvz[i]<<"]  " );
932         }
933 #endif // PARALLEL==1
934
935         // check other vars...?
936 }
937
938 template<class D>
939 void 
940 LbmFsgrSolver<D>::coarseCalculateFluxareas(int lev)
941 {
942         //LbmFloat calcCurrentMass = 0.0;
943         //LbmFloat calcCurrentVolume = 0.0;
944         //LbmFloat *ccel = NULL;
945         //LbmFloat *tcel = NULL;
946         //LbmFloat m[LBM_DFNUM];
947         //LbmFloat rho, ux, uy, uz, tmp, usqr;
948 #if OPT3D==1 
949         //LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM], lcsmomega;
950 #endif // OPT3D==true 
951         //m[0] = tmp = usqr = 0.0;
952
953         //for(int lev=0; lev<mMaxRefine; lev++) { TEST DEBUG
954         FSGR_FORIJK_BOUNDS(lev) {
955                 if( RFLAG(lev, i,j,k,mLevel[lev].setCurr) & CFFluid) {
956                         if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFGrFromCoarse) {
957                                 LbmFloat totArea = mFsgrCellArea[0]; // for l=0
958                                 for(int l=1; l<D::cDirNum; l++) { 
959                                         int ni=(2*i)+D::dfVecX[l], nj=(2*j)+D::dfVecY[l], nk=(2*k)+D::dfVecZ[l];
960                                         if(RFLAG(lev+1, ni,nj,nk, mLevel[lev+1].setCurr)&
961                                                         (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
962                                                         ) { 
963                                                 totArea += mFsgrCellArea[l];
964                                         }
965                                 } // l
966                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = totArea;
967                                 //continue;
968                         } else
969                         if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & (CFEmpty|CFUnused)) {
970                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 1.0;
971                                 //continue;
972                         } else {
973                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 0.0;
974                         }
975                 //errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) ); 
976                 }
977         } // } TEST DEBUG
978         if(!D::mSilent){ debMsgStd("coarseCalculateFluxareas",DM_MSG,"level "<<lev<<" calculated", 7); }
979 }
980         
981 template<class D>
982 void 
983 LbmFsgrSolver<D>::coarseAdvance(int lev)
984 {
985         LbmFloat calcCurrentMass = 0.0;
986         LbmFloat calcCurrentVolume = 0.0;
987
988         LbmFloat *ccel = NULL;
989         LbmFloat *tcel = NULL;
990         LbmFloat m[LBM_DFNUM];
991         LbmFloat rho, ux, uy, uz, tmp, usqr;
992 #if OPT3D==1 
993         LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM], lcsmomega;
994 #endif // OPT3D==true 
995         m[0] = tmp = usqr = 0.0;
996
997         coarseCalculateFluxareas(lev);
998         // copied from fineAdv.
999         CellFlagType *pFlagSrc = &RFLAG(lev, 1,1,getForZMin1(),SRCS(lev));
1000         CellFlagType *pFlagDst = &RFLAG(lev, 1,1,getForZMin1(),TSET(lev));
1001         pFlagSrc -= 1;
1002         pFlagDst -= 1;
1003         ccel = RACPNT(lev, 1,1,getForZMin1() ,SRCS(lev)); // QTEST
1004         ccel -= QCELLSTEP;
1005         tcel = RACPNT(lev, 1,1,getForZMin1() ,TSET(lev)); // QTEST
1006         tcel -= QCELLSTEP;
1007         //if(strstr(D::getName().c_str(),"Debug")){ errMsg("DEBUG","DEBUG!!!!!!!!!!!!!!!!!!!!!!!"); }
1008
1009         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
1010   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
1011   for(int i=1;i<mLevel[lev].lSizex-1;++i) {
1012 #if FSGR_STRICT_DEBUG==1
1013                 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
1014 #endif
1015                 pFlagSrc++;
1016                 pFlagDst++;
1017                 ccel += QCELLSTEP;
1018                 tcel += QCELLSTEP;
1019
1020                 // from coarse cells without unused nbs are not necessary...! -> remove
1021                 if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 
1022                         bool invNb = false;
1023                         FORDF1 { 
1024                                 if(RFLAG_NB(lev, i, j, k, SRCS(lev), l) & CFUnused) { invNb = true; }
1025                         }   
1026                         if(!invNb) {
1027                                 *pFlagSrc = CFFluid|CFGrNorm;
1028 #if ELBEEM_BLENDER!=1
1029                                 errMsg("coarseAdvance","FC2NRM_CHECK Converted CFGrFromCoarse to Norm at "<<lev<<" "<<PRINT_IJK);
1030 #endif // ELBEEM_BLENDER!=1
1031                                 // FIXME add debug check for these types of cells?, move to perform coarsening?
1032                         }
1033                 } // */
1034
1035                 //*(pFlagSrc+pFlagTarOff) = *pFlagSrc; // always set other set...
1036 #if FSGR_STRICT_DEBUG==1
1037                 *pFlagDst = *pFlagSrc; // always set other set...
1038 #else
1039                 *pFlagDst = (*pFlagSrc & (~CFGrCoarseInited)); // always set other set... , remove coarse inited flag
1040 #endif
1041
1042                 // old INTCFCOARSETEST==1
1043                 if((*pFlagSrc) & CFGrFromCoarse) {  // interpolateFineFromCoarse test!
1044                         if(( D::mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
1045                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
1046                         } else {
1047                                 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
1048                                 D::mNumUsedCells++;
1049                         }
1050                         continue; // interpolateFineFromCoarse test!
1051                 } // interpolateFineFromCoarse test! old INTCFCOARSETEST==1
1052
1053                 if( ((*pFlagSrc) & (CFFluid)) ) { 
1054                         ccel = RACPNT(lev, i,j,k ,SRCS(lev)); 
1055                         tcel = RACPNT(lev, i,j,k ,TSET(lev));
1056
1057                         if( ((*pFlagSrc) & (CFGrFromFine)) ) { 
1058                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }    // always copy...?
1059                                 continue; // comes from fine grid
1060                         }
1061                         // also ignore CFGrFromCoarse
1062                         else if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 
1063                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }    // always copy...?
1064                                 continue; 
1065                         }
1066
1067                         OPTIMIZED_STREAMCOLLIDE;
1068                         *pFlagDst |= CFNoBndFluid; // test?
1069                         calcCurrentVolume += RAC(ccel,dFlux); 
1070                         calcCurrentMass   += RAC(ccel,dFlux)*rho;
1071                         //ebugMarkCell(lev+1, 2*i+1,2*j+1,2*k  );
1072 #if FSGR_STRICT_DEBUG==1
1073                         if(rho<-1.0){ debugMarkCell(lev, i,j,k ); 
1074                                 errMsg("INVRHOCELL_CHECK"," l"<<lev<<" "<< PRINT_IJK<<" rho:"<<rho ); 
1075                                 D::mPanic = 1;
1076                         }
1077 #endif // FSGR_STRICT_DEBUG==1
1078                         D::mNumUsedCells++;
1079
1080                 }
1081         } 
1082         pFlagSrc+=2; // after x
1083         pFlagDst+=2; // after x
1084         ccel += (QCELLSTEP*2);
1085         tcel += (QCELLSTEP*2);
1086         } 
1087         pFlagSrc+= mLevel[lev].lSizex*2; // after y
1088         pFlagDst+= mLevel[lev].lSizex*2; // after y
1089         ccel += (QCELLSTEP*mLevel[lev].lSizex*2);
1090         tcel += (QCELLSTEP*mLevel[lev].lSizex*2);
1091         } // all cell loop k,j,i
1092         
1093
1094         //errMsg("coarseAdvance","level "<<lev<<" stepped from "<<mLevel[lev].setCurr<<" to "<<mLevel[lev].setOther);
1095         if(!D::mSilent){ errMsg("coarseAdvance","level "<<lev<<" stepped from "<<SRCS(lev)<<" to "<<TSET(lev)); }
1096         // */
1097
1098         // update other set
1099   mLevel[lev].setOther   = mLevel[lev].setCurr;
1100   mLevel[lev].setCurr   ^= 1;
1101   mLevel[lev].lsteps++;
1102   mLevel[lev].lmass   = calcCurrentMass   * mLevel[lev].lcellfactor;
1103   mLevel[lev].lvolume = calcCurrentVolume * mLevel[lev].lcellfactor;
1104 #ifndef ELBEEM_BLENDER
1105   errMsg("DFINI", " m l"<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<"  lcf="<< mLevel[lev].lcellfactor );
1106   errMsg("DFINI", " v l"<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<"  lcf="<< mLevel[lev].lcellfactor );
1107 #endif // ELBEEM_BLENDER
1108 }
1109
1110 /*****************************************************************************/
1111 //! multi level functions
1112 /*****************************************************************************/
1113
1114
1115 // get dfs from level (lev+1) to (lev) coarse border nodes
1116 template<class D>
1117 void 
1118 LbmFsgrSolver<D>::coarseRestrictFromFine(int lev)
1119 {
1120         if((lev<0) || ((lev+1)>mMaxRefine)) return;
1121 #if FSGR_STRICT_DEBUG==1
1122         // reset all unused cell values to invalid
1123         int unuCnt = 0;
1124         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
1125         for(int j=1;j<mLevel[lev].lSizey-1;++j) {
1126         for(int i=1;i<mLevel[lev].lSizex-1;++i) {
1127                 CellFlagType *pFlagSrc = &RFLAG(lev, i,j,k,mLevel[lev].setCurr);
1128                 if( ((*pFlagSrc) & (CFFluid|CFGrFromFine)) == (CFFluid|CFGrFromFine) ) { 
1129                         FORDF0{ QCELL(lev, i,j,k,mLevel[lev].setCurr, l) = -10000.0;    }
1130                         unuCnt++;
1131                         // set here
1132                 } else if( ((*pFlagSrc) & (CFFluid|CFGrNorm)) == (CFFluid|CFGrNorm) ) { 
1133                         // simulated...
1134                 } else {
1135                         // reset in interpolation
1136                         //errMsg("coarseRestrictFromFine"," reset l"<<lev<<" "<<PRINT_IJK);
1137                 }
1138                 if( ((*pFlagSrc) & (CFEmpty|CFUnused)) ) {  // test, also reset?
1139                         FORDF0{ QCELL(lev, i,j,k,mLevel[lev].setCurr, l) = -10000.0;    }
1140                 } // test
1141         } } }
1142         errMsg("coarseRestrictFromFine"," reset l"<<lev<<" fluid|coarseBorder cells: "<<unuCnt);
1143 #endif // FSGR_STRICT_DEBUG==1
1144         const int srcSet = mLevel[lev+1].setCurr;
1145         const int dstSet = mLevel[lev].setCurr;
1146
1147         LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0;                       
1148         LbmFloat *ccel = NULL;
1149         LbmFloat *tcel = NULL;
1150 #if OPT3D==1 
1151         LbmFloat m[LBM_DFNUM];
1152         // for macro add
1153         LbmFloat usqr;
1154         //LbmFloat *addfcel, *dstcell;
1155         LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM];
1156         LbmFloat lcsmDstOmega, lcsmSrcOmega, lcsmdfscale;
1157 #else // OPT3D==true 
1158         LbmFloat df[LBM_DFNUM];
1159         LbmFloat omegaDst, omegaSrc;
1160         LbmFloat feq[LBM_DFNUM];
1161         LbmFloat dfScale = mDfScaleUp;
1162 #endif // OPT3D==true 
1163
1164         LbmFloat mGaussw[27];
1165         LbmFloat totGaussw = 0.0;
1166         const LbmFloat alpha = 1.0;
1167         const LbmFloat gw = sqrt(2.0*D::cDimension);
1168 #ifndef ELBEEM_BLENDER
1169         errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 test df/dir num!");
1170 #endif
1171         for(int n=0;(n<D::cDirNum); n++) { mGaussw[n] = 0.0; }
1172         //for(int n=0;(n<D::cDirNum); n++) { 
1173         for(int n=0;(n<D::cDfNum); n++) { 
1174                 const LbmFloat d = norm(LbmVec(D::dfVecX[n], D::dfVecY[n], D::dfVecZ[n]));
1175                 LbmFloat w = expf( -alpha*d*d ) - expf( -alpha*gw*gw );
1176                 //errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 cell  n"<<n<<" d"<<d<<" w"<<w);
1177                 mGaussw[n] = w;
1178                 totGaussw += w;
1179         }
1180         for(int n=0;(n<D::cDirNum); n++) { 
1181                 mGaussw[n] = mGaussw[n]/totGaussw;
1182         }
1183
1184         //restrict
1185         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
1186         for(int j=1;j<mLevel[lev].lSizey-1;++j) {
1187         for(int i=1;i<mLevel[lev].lSizex-1;++i) {
1188                 CellFlagType *pFlagSrc = &RFLAG(lev, i,j,k,dstSet);
1189                 if((*pFlagSrc) & (CFFluid)) { 
1190                         if( ((*pFlagSrc) & (CFFluid|CFGrFromFine)) == (CFFluid|CFGrFromFine) ) { 
1191                                 // TODO? optimize?
1192                                 // do resctriction
1193                                 mNumInterdCells++;
1194                                 ccel = RACPNT(lev+1, 2*i,2*j,2*k,srcSet);
1195                                 tcel = RACPNT(lev  , i,j,k      ,dstSet);
1196
1197 #                               if OPT3D==0
1198                                 // add up weighted dfs
1199                                 FORDF0{ df[l] = 0.0;}
1200                                 for(int n=0;(n<D::cDirNum); n++) { 
1201                                         int ni=2*i+1*D::dfVecX[n], nj=2*j+1*D::dfVecY[n], nk=2*k+1*D::dfVecZ[n];
1202                                         ccel = RACPNT(lev+1, ni,nj,nk,srcSet);// CFINTTEST
1203                                         const LbmFloat weight = mGaussw[n];
1204                                         FORDF0{
1205                                                 LbmFloat cdf = weight * RAC(ccel,l);
1206 #                                               if FSGR_STRICT_DEBUG==1
1207                                                 if( cdf<-1.0 ){ errMsg("INVDFCREST_DFCHECK", PRINT_IJK<<" s"<<dstSet<<" from "<<PRINT_VEC(2*i,2*j,2*k)<<" s"<<srcSet<<" df"<<l<<":"<< df[l]); }
1208 #                                               endif
1209                                                 //errMsg("INVDFCREST_DFCHECK", PRINT_IJK<<" s"<<dstSet<<" from "<<PRINT_VEC(2*i,2*j,2*k)<<" s"<<srcSet<<" df"<<l<<":"<< df[l]<<" = "<<cdf<<" , w"<<weight); 
1210                                                 df[l] += cdf;
1211                                         }
1212                                 }
1213
1214                                 // calc rho etc. from weighted dfs
1215                                 rho = ux  = uy  = uz  = 0.0;
1216                                 FORDF0{
1217                                         LbmFloat cdf = df[l];
1218                                         rho += cdf; 
1219                                         ux  += (D::dfDvecX[l]*cdf); 
1220                                         uy  += (D::dfDvecY[l]*cdf);  
1221                                         uz  += (D::dfDvecZ[l]*cdf);  
1222                                 }
1223
1224                                 FORDF0{ feq[l] = D::getCollideEq(l, rho,ux,uy,uz); }
1225                                 if(mLevel[lev  ].lcsmago>0.0) {
1226                                         const LbmFloat Qo = D::getLesNoneqTensorCoeff(df,feq);
1227                                         omegaDst  = D::getLesOmega(mLevel[lev  ].omega,mLevel[lev  ].lcsmago,Qo);
1228                                         omegaSrc = D::getLesOmega(mLevel[lev+1].omega,mLevel[lev+1].lcsmago,Qo);
1229                                 } else {
1230                                         omegaDst = mLevel[lev+0].omega; /* NEWSMAGOT*/ 
1231                                         omegaSrc = mLevel[lev+1].omega;
1232                                 }
1233                                 dfScale   = (mLevel[lev  ].stepsize/mLevel[lev+1].stepsize)* (1.0/omegaDst-1.0)/ (1.0/omegaSrc-1.0); // yu
1234                                 FORDF0{
1235                                         RAC(tcel, l) = feq[l]+ (df[l]-feq[l])*dfScale;
1236                                 } 
1237 #                               else // OPT3D
1238                                 // similar to OPTIMIZED_STREAMCOLLIDE_UNUSED
1239                       
1240                                 //rho = ux = uy = uz = 0.0;
1241                                 MSRC_C  = CCELG_C(0) ;
1242                                 MSRC_N  = CCELG_N(0) ;
1243                                 MSRC_S  = CCELG_S(0) ;
1244                                 MSRC_E  = CCELG_E(0) ;
1245                                 MSRC_W  = CCELG_W(0) ;
1246                                 MSRC_T  = CCELG_T(0) ;
1247                                 MSRC_B  = CCELG_B(0) ;
1248                                 MSRC_NE = CCELG_NE(0);
1249                                 MSRC_NW = CCELG_NW(0);
1250                                 MSRC_SE = CCELG_SE(0);
1251                                 MSRC_SW = CCELG_SW(0);
1252                                 MSRC_NT = CCELG_NT(0);
1253                                 MSRC_NB = CCELG_NB(0);
1254                                 MSRC_ST = CCELG_ST(0);
1255                                 MSRC_SB = CCELG_SB(0);
1256                                 MSRC_ET = CCELG_ET(0);
1257                                 MSRC_EB = CCELG_EB(0);
1258                                 MSRC_WT = CCELG_WT(0);
1259                                 MSRC_WB = CCELG_WB(0);
1260                                 for(int n=1;(n<D::cDirNum); n++) { 
1261                                         ccel = RACPNT(lev+1,  2*i+1*D::dfVecX[n], 2*j+1*D::dfVecY[n], 2*k+1*D::dfVecZ[n]  ,srcSet);
1262                                         MSRC_C  += CCELG_C(n) ;
1263                                         MSRC_N  += CCELG_N(n) ;
1264                                         MSRC_S  += CCELG_S(n) ;
1265                                         MSRC_E  += CCELG_E(n) ;
1266                                         MSRC_W  += CCELG_W(n) ;
1267                                         MSRC_T  += CCELG_T(n) ;
1268                                         MSRC_B  += CCELG_B(n) ;
1269                                         MSRC_NE += CCELG_NE(n);
1270                                         MSRC_NW += CCELG_NW(n);
1271                                         MSRC_SE += CCELG_SE(n);
1272                                         MSRC_SW += CCELG_SW(n);
1273                                         MSRC_NT += CCELG_NT(n);
1274                                         MSRC_NB += CCELG_NB(n);
1275                                         MSRC_ST += CCELG_ST(n);
1276                                         MSRC_SB += CCELG_SB(n);
1277                                         MSRC_ET += CCELG_ET(n);
1278                                         MSRC_EB += CCELG_EB(n);
1279                                         MSRC_WT += CCELG_WT(n);
1280                                         MSRC_WB += CCELG_WB(n);
1281                                 }
1282                                 rho = MSRC_C  + MSRC_N + MSRC_S  + MSRC_E + MSRC_W  + MSRC_T  
1283                                         + MSRC_B  + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT 
1284                                         + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; 
1285                                 ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW 
1286                                         + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB;  
1287                                 uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW 
1288                                         + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB;  
1289                                 uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB 
1290                                         + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB;  
1291                                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
1292                                 \
1293                                 lcsmeq[dC] = EQC ; \
1294                                 COLL_CALCULATE_DFEQ(lcsmeq); \
1295                                 COLL_CALCULATE_NONEQTENSOR(lev+0, MSRC_ )\
1296                                 COLL_CALCULATE_CSMOMEGAVAL(lev+0, lcsmDstOmega); \
1297                                 COLL_CALCULATE_CSMOMEGAVAL(lev+1, lcsmSrcOmega); \
1298                                 \
1299                                 lcsmdfscale   = (mLevel[lev+0].stepsize/mLevel[lev+1].stepsize)* (1.0/lcsmDstOmega-1.0)/ (1.0/lcsmSrcOmega-1.0);  \
1300                                 RAC(tcel, dC ) = (lcsmeq[dC ] + (MSRC_C -lcsmeq[dC ] )*lcsmdfscale);
1301                                 RAC(tcel, dN ) = (lcsmeq[dN ] + (MSRC_N -lcsmeq[dN ] )*lcsmdfscale);
1302                                 RAC(tcel, dS ) = (lcsmeq[dS ] + (MSRC_S -lcsmeq[dS ] )*lcsmdfscale);
1303                                 RAC(tcel, dE ) = (lcsmeq[dE ] + (MSRC_E -lcsmeq[dE ] )*lcsmdfscale);
1304                                 RAC(tcel, dW ) = (lcsmeq[dW ] + (MSRC_W -lcsmeq[dW ] )*lcsmdfscale);
1305                                 RAC(tcel, dT ) = (lcsmeq[dT ] + (MSRC_T -lcsmeq[dT ] )*lcsmdfscale);
1306                                 RAC(tcel, dB ) = (lcsmeq[dB ] + (MSRC_B -lcsmeq[dB ] )*lcsmdfscale);
1307                                 RAC(tcel, dNE) = (lcsmeq[dNE] + (MSRC_NE-lcsmeq[dNE] )*lcsmdfscale);
1308                                 RAC(tcel, dNW) = (lcsmeq[dNW] + (MSRC_NW-lcsmeq[dNW] )*lcsmdfscale);
1309                                 RAC(tcel, dSE) = (lcsmeq[dSE] + (MSRC_SE-lcsmeq[dSE] )*lcsmdfscale);
1310                                 RAC(tcel, dSW) = (lcsmeq[dSW] + (MSRC_SW-lcsmeq[dSW] )*lcsmdfscale);
1311                                 RAC(tcel, dNT) = (lcsmeq[dNT] + (MSRC_NT-lcsmeq[dNT] )*lcsmdfscale);
1312                                 RAC(tcel, dNB) = (lcsmeq[dNB] + (MSRC_NB-lcsmeq[dNB] )*lcsmdfscale);
1313                                 RAC(tcel, dST) = (lcsmeq[dST] + (MSRC_ST-lcsmeq[dST] )*lcsmdfscale);
1314                                 RAC(tcel, dSB) = (lcsmeq[dSB] + (MSRC_SB-lcsmeq[dSB] )*lcsmdfscale);
1315                                 RAC(tcel, dET) = (lcsmeq[dET] + (MSRC_ET-lcsmeq[dET] )*lcsmdfscale);
1316                                 RAC(tcel, dEB) = (lcsmeq[dEB] + (MSRC_EB-lcsmeq[dEB] )*lcsmdfscale);
1317                                 RAC(tcel, dWT) = (lcsmeq[dWT] + (MSRC_WT-lcsmeq[dWT] )*lcsmdfscale);
1318                                 RAC(tcel, dWB) = (lcsmeq[dWB] + (MSRC_WB-lcsmeq[dWB] )*lcsmdfscale);
1319 #                               endif // OPT3D==0
1320
1321                                 //? if((lev<mMaxRefine)&&(D::cDimension==2)) { debugMarkCell(lev,i,j,k); }
1322 #                       if FSGR_STRICT_DEBUG==1
1323                                 //errMsg("coarseRestrictFromFine", "CRFF_DFDEBUG cell  "<<PRINT_IJK<<" rho:"<<rho<<" u:"<<PRINT_VEC(ux,uy,uz)<<" " ); 
1324 #                       endif // FSGR_STRICT_DEBUG==1
1325                                 D::mNumUsedCells++;
1326                         } // from fine & fluid
1327                         else {
1328                                 if(RFLAG(lev+1, 2*i,2*j,2*k,srcSet) & CFGrFromCoarse) {
1329                                         RFLAG(lev, i,j,k,dstSet) |= CFGrToFine;
1330                                 } else {
1331                                         RFLAG(lev, i,j,k,dstSet) &= (~CFGrToFine);
1332                                 }
1333                         }
1334                 } // & fluid
1335         }}}
1336         if(!D::mSilent){ errMsg("coarseRestrictFromFine"," from l"<<(lev+1)<<",s"<<mLevel[lev+1].setCurr<<" to l"<<lev<<",s"<<mLevel[lev].setCurr); }
1337 }
1338
1339 template<class D>
1340 bool 
1341 LbmFsgrSolver<D>::performRefinement(int lev) {
1342         if((lev<0) || ((lev+1)>mMaxRefine)) return false;
1343         bool change = false;
1344         //bool nbsok;
1345         // FIXME remove TIMEINTORDER ?
1346         LbmFloat interTime = 0.0;
1347         // update curr from other, as streaming afterwards works on curr
1348         // thus read only from srcSet, modify other
1349         const int srcSet = mLevel[lev].setOther;
1350         const int dstSet = mLevel[lev].setCurr;
1351         const int srcFineSet = mLevel[lev+1].setCurr;
1352         const bool debugRefinement = false;
1353
1354         // use template functions for 2D/3D
1355         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
1356   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
1357   for(int i=1;i<mLevel[lev].lSizex-1;++i) {
1358
1359                 if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) {
1360                         bool removeFromFine = false;
1361                         const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
1362                         CellFlagType reqType = CFGrNorm;
1363                         if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
1364                         
1365 #if REFINEMENTBORDER==1
1366                         if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet) & reqType) &&
1367                             (!(RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet) & (notAllowed)) )  ){ // ok
1368                         } else {
1369                                 removeFromFine=true;
1370                         }
1371                         /*if(strstr(D::getName().c_str(),"Debug"))
1372                         if(lev+1==mMaxRefine) { // mixborder
1373                                 for(int l=0;((l<D::cDirNum) && (!removeFromFine)); l++) {  // FARBORD
1374                                         int ni=2*i+2*D::dfVecX[l], nj=2*j+2*D::dfVecY[l], nk=2*k+2*D::dfVecZ[l];
1375                                         if(RFLAG(lev+1, ni,nj,nk, srcFineSet)&CFBnd) { // NEWREFT
1376                                                 removeFromFine=true;
1377                                         }
1378                                 }
1379                         } // FARBORD */
1380 #elif REFINEMENTBORDER==2 // REFINEMENTBORDER==1
1381                         FIX
1382                         for(int l=0;((l<D::cDirNum) && (!removeFromFine)); l++) { 
1383                                 int ni=2*i+D::dfVecX[l], nj=2*j+D::dfVecY[l], nk=2*k+D::dfVecZ[l];
1384                                 if(RFLAG(lev+1, ni,nj,nk, srcFineSet)&notSrcAllowed) { // NEWREFT
1385                                         removeFromFine=true;
1386                                 }
1387                         }
1388                         /*for(int l=0;((l<D::cDirNum) && (!removeFromFine)); l++) {  // FARBORD
1389                                 int ni=2*i+2*D::dfVecX[l], nj=2*j+2*D::dfVecY[l], nk=2*k+2*D::dfVecZ[l];
1390                                 if(RFLAG(lev+1, ni,nj,nk, srcFineSet)&notSrcAllowed) { // NEWREFT
1391                                         removeFromFine=true;
1392                                 }
1393                         } // FARBORD */
1394 #elif REFINEMENTBORDER==3 // REFINEMENTBORDER==1
1395                         FIX
1396                         if(lev+1==mMaxRefine) { // mixborder
1397                                 if(RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)&notSrcAllowed) { 
1398                                         removeFromFine=true;
1399                                 }
1400                         } else { // mixborder
1401                                 for(int l=0; l<D::cDirNum; l++) { 
1402                                         int ni=2*i+D::dfVecX[l], nj=2*j+D::dfVecY[l], nk=2*k+D::dfVecZ[l];
1403                                         if(RFLAG(lev+1, ni,nj,nk, srcFineSet)&notSrcAllowed) { // NEWREFT
1404                                                 removeFromFine=true;
1405                                         }
1406                                 }
1407                         } // mixborder
1408                         // also remove from fine cells that are above from fine
1409 #else // REFINEMENTBORDER==1
1410                         ERROR
1411 #endif // REFINEMENTBORDER==1
1412
1413                         if(removeFromFine) {
1414                                 // dont turn CFGrFromFine above interface cells into CFGrNorm
1415                                 //errMsg("performRefinement","Removing CFGrFromFine on lev"<<lev<<" " <<PRINT_IJK<<" srcflag:"<<convertCellFlagType2String(RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet)) <<" set:"<<dstSet );
1416                                 RFLAG(lev, i,j,k, dstSet) = CFEmpty;
1417 #if FSGR_STRICT_DEBUG==1
1418                                 // for interpolation later on during fine grid fixing
1419                                 // these cells are still correctly inited
1420                                 RFLAG(lev, i,j,k, dstSet) |= CFGrCoarseInited;  // remove later on? FIXME?
1421 #endif // FSGR_STRICT_DEBUG==1
1422                                 //RFLAG(lev, i,j,k, mLevel[lev].setOther) = CFEmpty; // FLAGTEST
1423                                 if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev,i,j,k); 
1424                                 change=true;
1425                                 mNumFsgrChanges++;
1426                                 for(int l=1; l<D::cDirNum; l++) { 
1427                                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
1428                                         //errMsg("performRefinement","On lev:"<<lev<<" check: "<<PRINT_VEC(ni,nj,nk)<<" set:"<<dstSet<<" = "<<convertCellFlagType2String(RFLAG(lev, ni,nj,nk, srcSet)) );
1429                                         if( (  RFLAG(lev, ni,nj,nk, srcSet)&CFFluid      ) &&
1430                                                         (!(RFLAG(lev, ni,nj,nk, srcSet)&CFGrFromFine)) ) { // dont change status of nb. from fine cells
1431                                                 // tag as inited for debugging, cell contains fluid DFs anyway
1432                                                 RFLAG(lev, ni,nj,nk, dstSet) = CFFluid|CFGrFromFine|CFGrCoarseInited;
1433                                                 //errMsg("performRefinement","On lev:"<<lev<<" set to from fine: "<<PRINT_VEC(ni,nj,nk)<<" set:"<<dstSet);
1434                                                 //if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 
1435                                         }
1436                                 } // l 
1437
1438                                 // FIXME fix fine level?
1439                         }
1440
1441                         // recheck from fine flag
1442                 }
1443         }}} // TEST
1444
1445
1446         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { // TEST
1447   for(int j=1;j<mLevel[lev].lSizey-1;++j) { // TEST
1448   for(int i=1;i<mLevel[lev].lSizex-1;++i) { // TEST
1449
1450                 // test from coarseAdvance
1451                 // from coarse cells without unused nbs are not necessary...! -> remove
1452                 /*if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 
1453                         bool invNb = false;
1454                         FORDF1 { 
1455                                 if(RFLAG_NB(lev, i, j, k, SRCS(lev), l) & CFUnused) { invNb = true; }
1456                         }   
1457                         if(!invNb) {
1458                                 *pFlagSrc = CFFluid|CFGrNorm;
1459                                 errMsg("coarseAdvance","FC2NRM_CHECK Converted CFGrFromCoarse to Norm at "<<lev<<" "<<PRINT_IJK);
1460                         }
1461                 } // */
1462
1463                 if(RFLAG(lev, i,j,k, srcSet) & CFGrFromCoarse) {
1464
1465                         // from coarse cells without unused nbs are not necessary...! -> remove
1466                         bool invNb = false;
1467                         bool fluidNb = false;
1468                         for(int l=1; l<D::cDirNum; l++) { 
1469                                 if(RFLAG_NB(lev, i, j, k, srcSet, l) & CFUnused) { invNb = true; }
1470                                 if(RFLAG_NB(lev, i, j, k, srcSet, l) & (CFGrNorm)) { fluidNb = true; }
1471                         }   
1472                         if(!invNb) {
1473                                 // no unused cells around -> calculate normally from now on
1474                                 RFLAG(lev, i,j,k, dstSet) = CFFluid|CFGrNorm;
1475                                 if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 
1476                                 change=true;
1477                                 mNumFsgrChanges++;
1478                         } // from advance */
1479                         if(!fluidNb) {
1480                                 // no fluid cells near -> no transfer necessary
1481                                 RFLAG(lev, i,j,k, dstSet) = CFUnused;
1482                                 //RFLAG(lev, i,j,k, mLevel[lev].setOther) = CFUnused; // FLAGTEST
1483                                 if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 
1484                                 change=true;
1485                                 mNumFsgrChanges++;
1486                         } // from advance */
1487
1488
1489                         // dont allow double transfer
1490                         // this might require fixing the neighborhood
1491                         if(RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)&(CFGrFromCoarse)) { 
1492                                 // dont turn CFGrFromFine above interface cells into CFGrNorm
1493                                 //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<lev<<" " <<PRINT_IJK<<" due to finer from coarse cell " );
1494                                 RFLAG(lev, i,j,k, dstSet) = CFFluid|CFGrNorm;
1495                                 if(lev>0) RFLAG(lev-1, i/2,j/2,k/2, mLevel[lev-1].setCurr) &= (~CFGrToFine); // TODO add more of these?
1496                                 if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 
1497                                 change=true;
1498                                 mNumFsgrChanges++;
1499                                 for(int l=1; l<D::cDirNum; l++) { 
1500                                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
1501                                         if(RFLAG(lev, ni,nj,nk, srcSet)&(CFGrNorm)) { //ok
1502                                                 for(int m=1; m<D::cDirNum; m++) { 
1503                                                         int mi=  ni +D::dfVecX[m], mj=  nj +D::dfVecY[m], mk=  nk +D::dfVecZ[m];
1504                                                         if(RFLAG(lev,  mi, mj, mk, srcSet)&CFUnused) {
1505                                                                 // norm cells in neighborhood with unused nbs have to be new border...
1506                                                                 RFLAG(lev, ni,nj,nk, dstSet) = CFFluid|CFGrFromCoarse;
1507                                                                 if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 
1508                                                         }
1509                                                 }
1510                                                 // these alreay have valid values...
1511                                         }
1512                                         else if(RFLAG(lev, ni,nj,nk, srcSet)&(CFUnused)) { //ok
1513                                                 // this should work because we have a valid neighborhood here for now
1514                                                 interpolateCellFromCoarse(lev,  ni, nj, nk, dstSet /*mLevel[lev].setCurr*/, interTime, CFFluid|CFGrFromCoarse, false);
1515                                                 if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 
1516                                                 mNumFsgrChanges++;
1517                                         }
1518                                 } // l 
1519                         } // double transer
1520
1521                 } // from coarse
1522
1523         } } }
1524
1525
1526         // fix dstSet from fine cells here
1527         // warning - checks CFGrFromFine on dstset changed before!
1528         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { // TEST
1529   for(int j=1;j<mLevel[lev].lSizey-1;++j) { // TEST
1530   for(int i=1;i<mLevel[lev].lSizex-1;++i) { // TEST
1531
1532                 //if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) {
1533                 if(RFLAG(lev, i,j,k, dstSet) & CFGrFromFine) {
1534                         // modify finer level border
1535                         if((RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)&(CFGrFromCoarse))) { 
1536                                 //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<(lev+1)<<" from l"<<lev<<" " <<PRINT_IJK );
1537                                 CellFlagType setf = CFFluid;
1538                                 if(lev+1 < mMaxRefine) setf = CFFluid|CFGrNorm;
1539                                 RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)=setf;
1540                                 change=true;
1541                                 mNumFsgrChanges++;
1542                                 for(int l=1; l<D::cDirNum; l++) { 
1543                                         int bi=(2*i)+D::dfVecX[l], bj=(2*j)+D::dfVecY[l], bk=(2*k)+D::dfVecZ[l];
1544                                         if(RFLAG(lev+1,  bi, bj, bk, srcFineSet)&(CFGrFromCoarse)) {
1545                                                 //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<(lev+1)<<" "<<PRINT_VEC(bi,bj,bk) );
1546                                                 RFLAG(lev+1,  bi, bj, bk, srcFineSet) = setf;
1547                                                 if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev+1,bi,bj,bk); 
1548                                         }
1549                                         else if(RFLAG(lev+1,  bi, bj, bk, srcFineSet)&(CFUnused      )) { 
1550                                                 //errMsg("performRefinement","Removing CFUnused on lev"<<(lev+1)<<" "<<PRINT_VEC(bi,bj,bk) );
1551                                                 interpolateCellFromCoarse(lev+1,  bi, bj, bk, srcFineSet, interTime, setf, false);
1552                                                 if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev+1,bi,bj,bk); 
1553                                                 mNumFsgrChanges++;
1554                                         }
1555                                 }
1556                                 for(int l=1; l<D::cDirNum; l++) { 
1557                                         int bi=(2*i)+D::dfVecX[l], bj=(2*j)+D::dfVecY[l], bk=(2*k)+D::dfVecZ[l];
1558                                         if(   (RFLAG(lev+1,  bi, bj, bk, srcFineSet)&CFFluid       ) &&
1559                                                         (!(RFLAG(lev+1,  bi, bj, bk, srcFineSet)&CFGrFromCoarse)) ) {
1560                                                 // all unused nbs now of coarse have to be from coarse
1561                                                 for(int m=1; m<D::cDirNum; m++) { 
1562                                                         int mi=  bi +D::dfVecX[m], mj=  bj +D::dfVecY[m], mk=  bk +D::dfVecZ[m];
1563                                                         if(RFLAG(lev+1,  mi, mj, mk, srcFineSet)&CFUnused) {
1564                                                                 //errMsg("performRefinement","Changing CFUnused on lev"<<(lev+1)<<" "<<PRINT_VEC(mi,mj,mk) );
1565                                                                 interpolateCellFromCoarse(lev+1,  mi, mj, mk, srcFineSet, interTime, CFFluid|CFGrFromCoarse, false);
1566                                                                 if((D::cDimension==2)&&(debugRefinement)) debugMarkCell(lev+1,mi,mj,mk); 
1567                                                                 mNumFsgrChanges++;
1568                                                         }
1569                                                 }
1570                                                 // nbs prepared...
1571                                         }
1572                                 }
1573                         }
1574                         
1575                 } // convert regions of from fine
1576         }}} // TEST
1577
1578         if(!D::mSilent){ errMsg("performRefinement"," for l"<<lev<<" done ("<<change<<") " ); }
1579         return change;
1580 }
1581
1582
1583 // done after refinement
1584 template<class D>
1585 bool 
1586 LbmFsgrSolver<D>::performCoarsening(int lev) {
1587         //if(D::mInitDone){ errMsg("performCoarsening","skip"); return 0;} // DEBUG
1588                                         
1589         if((lev<0) || ((lev+1)>mMaxRefine)) return false;
1590         bool change = false;
1591         bool nbsok;
1592         // hence work on modified curr set
1593         const int srcSet = mLevel[lev].setCurr;
1594         const int dstlev = lev+1;
1595         const int dstFineSet = mLevel[dstlev].setCurr;
1596         const bool debugCoarsening = false;
1597
1598         // use template functions for 2D/3D
1599         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
1600   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
1601   for(int i=1;i<mLevel[lev].lSizex-1;++i) {
1602
1603                         // from coarse cells without unused nbs are not necessary...! -> remove
1604                         // perform check from coarseAdvance here?
1605                         if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) {
1606                                 // remove from fine cells now that are completely in fluid
1607                                 // FIXME? check that new from fine in performRefinement never get deleted here afterwards?
1608                                 // or more general, one cell never changed more than once?
1609                                 const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
1610                                 //const CellFlagType notNbAllowed = (CFInter|CFBnd|CFGrFromFine); unused
1611                                 CellFlagType reqType = CFGrNorm;
1612                                 if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
1613
1614                                 nbsok = true;
1615                                 for(int l=0; l<D::cDirNum && nbsok; l++) { 
1616                                         int ni=(2*i)+D::dfVecX[l], nj=(2*j)+D::dfVecY[l], nk=(2*k)+D::dfVecZ[l];
1617                                         if(   (RFLAG(lev+1, ni,nj,nk, dstFineSet) & reqType) &&
1618                                                         (!(RFLAG(lev+1, ni,nj,nk, dstFineSet) & (notAllowed)) )  ){
1619                                                 // ok
1620                                         } else {
1621                                                 nbsok=false;
1622                                         }
1623                                         /*if(strstr(D::getName().c_str(),"Debug"))
1624                                         if((nbsok)&&(lev+1==mMaxRefine)) { // mixborder
1625                                                 for(int l=0;((l<D::cDirNum) && (nbsok)); l++) {  // FARBORD
1626                                                         int ni=2*i+2*D::dfVecX[l], nj=2*j+2*D::dfVecY[l], nk=2*k+2*D::dfVecZ[l];
1627                                                         if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&CFBnd) { // NEWREFT
1628                                                                 nbsok=false;
1629                                                         }
1630                                                 }
1631                                         } // FARBORD */
1632                                 }
1633                                 // dont turn CFGrFromFine above interface cells into CFGrNorm
1634                                 // now check nbs on same level
1635                                 for(int l=1; l<D::cDirNum && nbsok; l++) { 
1636                                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
1637                                         if(RFLAG(lev, ni,nj,nk, srcSet)&(CFFluid)) { //ok
1638                                         } else {
1639                                                 nbsok = false;
1640                                         }
1641                                 } // l
1642
1643                                 if(nbsok) {
1644                                         // conversion to coarse fluid cell
1645                                         change = true;
1646                                         mNumFsgrChanges++;
1647                                         RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrNorm;
1648                                         // dfs are already ok...
1649                                         //if(D::mInitDone) errMsg("performCoarsening","CFGrFromFine changed to CFGrNorm at lev"<<lev<<" " <<PRINT_IJK );
1650                                         if((D::cDimension==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 
1651
1652                                         // only check complete cubes
1653                                         for(int dx=-1;dx<=1;dx+=2) {
1654                                         for(int dy=-1;dy<=1;dy+=2) {
1655                                         for(int dz=-1*(LBMDIM&1);dz<=1*(LBMDIM&1);dz+=2) { // 2d/3d
1656                                                 // check for norm and from coarse, as the coarse level might just have been refined...
1657                                                 /*if(D::mInitDone) errMsg("performCoarsening","CFGrFromFine subc check "<< "x"<<convertCellFlagType2String( RFLAG(lev, i+dx, j   , k   ,  srcSet))<<" "
1658                                                                         "y"<<convertCellFlagType2String( RFLAG(lev, i   , j+dy, k   ,  srcSet))<<" " "z"<<convertCellFlagType2String( RFLAG(lev, i   , j   , k+dz,  srcSet))<<" "
1659                                                                         "xy"<<convertCellFlagType2String( RFLAG(lev, i+dx, j+dy, k   ,  srcSet))<<" " "xz"<<convertCellFlagType2String( RFLAG(lev, i+dx, j   , k+dz,  srcSet))<<" "
1660                                                                         "yz"<<convertCellFlagType2String( RFLAG(lev, i   , j+dy, k+dz,  srcSet))<<" " "xyz"<<convertCellFlagType2String( RFLAG(lev, i+dx, j+dy, k+dz,  srcSet))<<" " ); // */
1661                                                 if( 
1662                                                                 // we now the flag of the current cell! ( RFLAG(lev, i   , j   , k   ,  srcSet)&(CFGrNorm)) &&
1663                                                                 ( RFLAG(lev, i+dx, j   , k   ,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
1664                                                                 ( RFLAG(lev, i   , j+dy, k   ,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
1665                                                                 ( RFLAG(lev, i   , j   , k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
1666
1667                                                                 ( RFLAG(lev, i+dx, j+dy, k   ,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
1668                                                                 ( RFLAG(lev, i+dx, j   , k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
1669                                                                 ( RFLAG(lev, i   , j+dy, k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
1670                                                                 ( RFLAG(lev, i+dx, j+dy, k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) 
1671                                                         ) {
1672                                                         // middle source node on higher level
1673                                                         int dstx = (2*i)+dx;
1674                                                         int dsty = (2*j)+dy;
1675                                                         int dstz = (2*k)+dz;
1676
1677                                                         mNumFsgrChanges++;
1678                                                         RFLAG(dstlev, dstx,dsty,dstz, dstFineSet) = CFUnused;
1679                                                         RFLAG(dstlev, dstx,dsty,dstz, mLevel[dstlev].setOther) = CFUnused; // FLAGTEST
1680                                                         //if(D::mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init center unused set l"<<dstlev<<" at "<<PRINT_VEC(dstx,dsty,dstz) );
1681
1682                                                         for(int l=1; l<D::cDirNum; l++) { 
1683                                                                 int dstni=dstx+D::dfVecX[l], dstnj=dsty+D::dfVecY[l], dstnk=dstz+D::dfVecZ[l];
1684                                                                 if(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFFluid)) { 
1685                                                                         RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFFluid|CFGrFromCoarse;
1686                                                                 }
1687                                                                 if(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFInter)) { 
1688                                                                         //if(D::mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init CHECK Warning - deleting interface cell...");
1689                                                                         D::mFixMass += QCELL( dstlev, dstni,dstnj,dstnk, dstFineSet, dMass);
1690                                                                         RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFFluid|CFGrFromCoarse;
1691                                                                 }
1692                                                         } // l
1693
1694                                                         // again check nb flags of all surrounding cells to see if any from coarse
1695                                                         // can be convted to unused
1696                                                         for(int l=1; l<D::cDirNum; l++) { 
1697                                                                 int dstni=dstx+D::dfVecX[l], dstnj=dsty+D::dfVecY[l], dstnk=dstz+D::dfVecZ[l];
1698                                                                 // have to be at least from coarse here...
1699                                                                 //errMsg("performCoarsening","CFGrFromFine subcube init unused check l"<<dstlev<<" at "<<PRINT_VEC(dstni,dstnj,dstnk)<<" "<< convertCellFlagType2String(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)) );
1700                                                                 if(!(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFUnused) )) { 
1701                                                                         bool delok = true;
1702                                                                         // careful long range here... check domain bounds?
1703                                                                         for(int m=1; m<D::cDirNum; m++) {                                                                               
1704                                                                                 int chkni=dstni+D::dfVecX[m], chknj=dstnj+D::dfVecY[m], chknk=dstnk+D::dfVecZ[m];
1705                                                                                 if(RFLAG(dstlev, chkni,chknj,chknk, dstFineSet)&(CFUnused|CFGrFromCoarse)) { 
1706                                                                                         // this nb cell is ok for deletion
1707                                                                                 } else { 
1708                                                                                         delok=false; // keep it!
1709                                                                                 }
1710                                                                                 //errMsg("performCoarsening"," CHECK "<<PRINT_VEC(dstni,dstnj,dstnk)<<" to "<<PRINT_VEC( chkni,chknj,chknk )<<" f:"<< convertCellFlagType2String( RFLAG(dstlev, chkni,chknj,chknk, dstFineSet))<<" nbsok"<<delok );
1711                                                                         }
1712                                                                         //errMsg("performCoarsening","CFGrFromFine subcube init unused check l"<<dstlev<<" at "<<PRINT_VEC(dstni,dstnj,dstnk)<<" ok"<<delok );
1713                                                                         if(delok) {
1714                                                                                 mNumFsgrChanges++;
1715                                                                                 RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFUnused;
1716                                                                                 RFLAG(dstlev, dstni,dstnj,dstnk, mLevel[dstlev].setOther) = CFUnused; // FLAGTEST
1717                                                                                 if((D::cDimension==2)&&(debugCoarsening)) debugMarkCell(dstlev,dstni,dstnj,dstnk); 
1718                                                                         }
1719                                                                 }
1720                                                         } // l
1721                                                         // treat subcube
1722                                                         //ebugMarkCell(lev,i+dx,j+dy,k+dz); 
1723                                                         //if(D::mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init, dir:"<<PRINT_VEC(dx,dy,dz) );
1724                                                 }
1725                                         } } }
1726
1727                                 }   // ?
1728                         } // convert regions of from fine
1729         }}} // TEST!
1730
1731                                         // reinit cell area value
1732                                         /*if( RFLAG(lev, i,j,k,srcSet) & CFFluid) {
1733                                                 if( RFLAG(lev+1, i*2,j*2,k*2,dstFineSet) & CFGrFromCoarse) {
1734                                                         LbmFloat totArea = mFsgrCellArea[0]; // for l=0
1735                                                         for(int l=1; l<D::cDirNum; l++) { 
1736                                                                 int ni=(2*i)+D::dfVecX[l], nj=(2*j)+D::dfVecY[l], nk=(2*k)+D::dfVecZ[l];
1737                                                                 if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&
1738                                                                                 (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
1739                                                                                 //(CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
1740                                                                                 ) { 
1741                                                                         //LbmFloat area = 0.25; if(D::dfVecX[l]!=0) area *= 0.5; if(D::dfVecY[l]!=0) area *= 0.5; if(D::dfVecZ[l]!=0) area *= 0.5;
1742                                                                         totArea += mFsgrCellArea[l];
1743                                                                 }
1744                                                         } // l
1745                                                         QCELL(lev, i,j,k,mLevel[lev].setOther, dFlux) = 
1746                                                         QCELL(lev, i,j,k,srcSet, dFlux) = totArea;
1747                                                 } else {
1748                                                         QCELL(lev, i,j,k,mLevel[lev].setOther, dFlux) = 
1749                                                         QCELL(lev, i,j,k,srcSet, dFlux) = 1.0;
1750                                                 }
1751                                                 //errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,srcSet, dFlux) );
1752                                         }
1753                                 // */
1754
1755         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
1756   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
1757   for(int i=1;i<mLevel[lev].lSizex-1;++i) {
1758
1759
1760                         if(RFLAG(lev, i,j,k, srcSet) & CFEmpty) {
1761                                 // check empty -> from fine conversion
1762                                 bool changeToFromFine = false;
1763                                 const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
1764                                 CellFlagType reqType = CFGrNorm;
1765                                 if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
1766
1767 #if REFINEMENTBORDER==1
1768                                 if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & reqType) &&
1769                                     (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) )  ){
1770                                         changeToFromFine=true;
1771                                 }
1772                         /*if(strstr(D::getName().c_str(),"Debug"))
1773                         if((changeToFromFine)&&(lev+1==mMaxRefine)) { // mixborder
1774                                 for(int l=0;((l<D::cDirNum) && (changeToFromFine)); l++) {  // FARBORD
1775                                         int ni=2*i+2*D::dfVecX[l], nj=2*j+2*D::dfVecY[l], nk=2*k+2*D::dfVecZ[l];
1776                                         if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&CFBnd) { // NEWREFT
1777                                                 changeToFromFine=false;
1778                                         }
1779                                 } 
1780                         }// FARBORD */
1781 #elif REFINEMENTBORDER==2 // REFINEMENTBORDER==1
1782                                 if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & reqType) &&
1783                                     (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) )  ){
1784                                         changeToFromFine=true;
1785                                         for(int l=0; ((l<D::cDirNum)&&(changeToFromFine)); l++) { 
1786                                                 int ni=2*i+D::dfVecX[l], nj=2*j+D::dfVecY[l], nk=2*k+D::dfVecZ[l];
1787                                                 if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&(notNbAllowed)) { // NEWREFT
1788                                                         changeToFromFine=false;
1789                                                 }
1790                                         }
1791                                         /*for(int l=0; ((l<D::cDirNum)&&(changeToFromFine)); l++) {  // FARBORD
1792                                                 int ni=2*i+2*D::dfVecX[l], nj=2*j+2*D::dfVecY[l], nk=2*k+2*D::dfVecZ[l];
1793                                                 if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&(notNbAllowed)) { // NEWREFT
1794                                                         changeToFromFine=false;
1795                                                 }
1796                                         } // FARBORD*/
1797                                 }
1798 #elif REFINEMENTBORDER==3 // REFINEMENTBORDER==3
1799                                 FIX!!!
1800                                 if(lev+1==mMaxRefine) { // mixborder
1801                                         if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (CFFluid|CFInter)) &&
1802                                                         (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) )  ){
1803                                                 changeToFromFine=true;
1804                                         }
1805                                 } else {
1806                                 if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (CFFluid)) &&
1807                                     (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) )  ){
1808                                         changeToFromFine=true;
1809                                         for(int l=0; l<D::cDirNum; l++) { 
1810                                                 int ni=2*i+D::dfVecX[l], nj=2*j+D::dfVecY[l], nk=2*k+D::dfVecZ[l];
1811                                                 if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&(notNbAllowed)) { // NEWREFT
1812                                                         changeToFromFine=false;
1813                                                 }
1814                                         }
1815                                 } } // mixborder
1816 #else // REFINEMENTBORDER==3
1817                                 ERROR
1818 #endif // REFINEMENTBORDER==1
1819                                 if(changeToFromFine) {
1820                                         change = true;
1821                                         mNumFsgrChanges++;
1822                                         RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrFromFine;
1823                                         if((D::cDimension==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 
1824                                         // same as restr from fine func! not necessary ?!
1825                                         // coarseRestrictFromFine part */
1826                                 }
1827                         } // only check empty cells
1828
1829         }}} // TEST!
1830
1831         if(!D::mSilent){ errMsg("performCoarsening"," for l"<<lev<<" done " ); }
1832         return change;
1833 }
1834
1835
1836 /*****************************************************************************/
1837 /*! perform a single LBM step */
1838 /*****************************************************************************/
1839 template<class D>
1840 void 
1841 LbmFsgrSolver<D>::adaptTimestep()
1842 {
1843         LbmFloat massTOld=0.0, massTNew=0.0;
1844         LbmFloat volTOld=0.0, volTNew=0.0;
1845
1846         bool rescale = false;  // do any rescale at all?
1847         LbmFloat scaleFac = -1.0; // timestep scaling
1848
1849         LbmFloat levOldOmega[FSGR_MAXNOOFLEVELS];
1850         LbmFloat levOldStepsize[FSGR_MAXNOOFLEVELS];
1851         for(int lev=mMaxRefine; lev>=0 ; lev--) {
1852                 levOldOmega[lev] = mLevel[lev].omega;
1853                 levOldStepsize[lev] = mLevel[lev].stepsize;
1854         }
1855         //if(mTimeSwitchCounts>0){ errMsg("DEB CSKIP",""); return; } // DEBUG
1856
1857         LbmFloat fac = 0.8;          // modify time step by 20%, TODO? do multiple times for large changes?
1858         LbmFloat diffPercent = 0.05; // dont scale if less than 5%
1859         LbmFloat allowMax = D::mpParam->getTadapMaxSpeed();  // maximum allowed velocity
1860         LbmFloat nextmax = D::mpParam->getSimulationMaxSpeed() + norm(mLevel[mMaxRefine].gravity);
1861
1862         //newdt = D::mpParam->getStepTime() * (allowMax/nextmax);
1863         LbmFloat newdt = D::mpParam->getStepTime(); // newtr
1864         if(nextmax>allowMax/fac) {
1865                 newdt = D::mpParam->getStepTime() * fac;
1866         } else {
1867                 if(nextmax<allowMax*fac) {
1868                         newdt = D::mpParam->getStepTime() / fac;
1869                 }
1870         } // newtr
1871         //errMsg("LbmFsgrSolver::adaptTimestep","nextmax="<<nextmax<<" allowMax="<<allowMax<<" fac="<<fac<<" simmaxv="<< D::mpParam->getSimulationMaxSpeed() );
1872
1873         bool minCutoff = false;
1874         LbmFloat desireddt = newdt;
1875         if(newdt>D::mpParam->getMaxStepTime()){ newdt = D::mpParam->getMaxStepTime(); }
1876         if(newdt<D::mpParam->getMinStepTime()){ 
1877                 newdt = D::mpParam->getMinStepTime(); 
1878                 if(nextmax>allowMax/fac){       minCutoff=true; } // only if really large vels...
1879         }
1880
1881         LbmFloat dtdiff = fabs(newdt - D::mpParam->getStepTime());
1882         if(!D::mSilent) {
1883                 debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt<<" max"<<D::mpParam->getMaxStepTime()<<" min"<<D::mpParam->getMinStepTime()<<" diff"<<dtdiff<<
1884                         " simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxStepTime)<<" maxsteps:"<<(mSimulationTime/mMinStepTime) , 10); }
1885
1886         // in range, and more than X% change?
1887         //if( newdt <  D::mpParam->getStepTime() ) // DEBUG
1888         LbmFloat rhoAvg = mCurrentMass/mCurrentVolume;
1889         if( (newdt<=D::mpParam->getMaxStepTime()) && (newdt>=D::mpParam->getMinStepTime()) 
1890                         && (dtdiff>(D::mpParam->getStepTime()*diffPercent)) ) {
1891                 if((newdt>levOldStepsize[mMaxRefine])&&(mTimestepReduceLock)) {
1892                         // wait some more...
1893                         //debMsgNnl("LbmFsgrSolver::TAdp",DM_NOTIFY," Delayed... "<<mTimestepReduceLock<<" ",10);
1894                         debMsgDirect("D");
1895                 } else {
1896                         D::mpParam->setDesiredStepTime( newdt );
1897                         rescale = true;
1898                         if(!D::mSilent) {
1899                                 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"\n\n\n\n",10);
1900                                 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: new="<<newdt<<" old="<<D::mpParam->getStepTime()<<" maxSpeed:"<<D::mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<D::mStepCnt, 10 );
1901                                 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: "<<
1902                                                 "rhoAvg="<<rhoAvg<<" cMass="<<mCurrentMass<<" cVol="<<mCurrentVolume,10);
1903                         }
1904                 } // really change dt
1905         }
1906
1907         if(mTimestepReduceLock>0) mTimestepReduceLock--;
1908
1909         
1910         /*
1911         // forced back and forth switchting (for testing)
1912         const int tadtogInter = 300;
1913         const double tadtogSwitch = 0.66;
1914         errMsg("TIMESWITCHTOGGLETEST","warning enabled "<< tadtogSwitch<<","<<tadtogSwitch<<" !!!!!!!!!!!!!!!!!!!");
1915         if( ((D::mStepCnt% tadtogInter)== (tadtogInter/4*1)-1) ||
1916             ((D::mStepCnt% tadtogInter)== (tadtogInter/4*2)-1) ){
1917                 rescale = true; minCutoff = false;
1918                 newdt = tadtogSwitch * D::mpParam->getStepTime();
1919                 D::mpParam->setDesiredStepTime( newdt );
1920         } else 
1921         if( ((D::mStepCnt% tadtogInter)== (tadtogInter/4*3)-1) ||
1922             ((D::mStepCnt% tadtogInter)== (tadtogInter/4*4)-1) ){
1923                 rescale = true; minCutoff = false;
1924                 newdt = D::mpParam->getStepTime()/tadtogSwitch ;
1925                 D::mpParam->setDesiredStepTime( newdt );
1926         } else {
1927                 rescale = false; minCutoff = false;
1928         }
1929         // */
1930
1931         // test mass rescale
1932
1933         scaleFac = newdt/D::mpParam->getStepTime();
1934         if(rescale) {
1935                 // fixme - warum y, wird jetzt gemittelt...
1936                 mTimestepReduceLock = 4*(mLevel[mMaxRefine].lSizey+mLevel[mMaxRefine].lSizez+mLevel[mMaxRefine].lSizex)/3;
1937
1938                 mTimeSwitchCounts++;
1939                 D::mpParam->calculateAllMissingValues( D::mSilent );
1940                 recalculateObjectSpeeds();
1941                 // calc omega, force for all levels
1942                 initLevelOmegas();
1943                 if(D::mpParam->getStepTime()<mMinStepTime) mMinStepTime = D::mpParam->getStepTime();
1944                 if(D::mpParam->getStepTime()>mMaxStepTime) mMaxStepTime = D::mpParam->getStepTime();
1945
1946                 for(int lev=mMaxRefine; lev>=0 ; lev--) {
1947                         LbmFloat newSteptime = mLevel[lev].stepsize;
1948                         LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]);
1949
1950                         if(!D::mSilent) {
1951                                 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Level: "<<lev<<" Timestep change: "<<
1952                                                 " scaleFac="<<dfScaleFac<<" newDt="<<newSteptime<<" newOmega="<<mLevel[lev].omega,10);
1953                         }
1954                         if(lev!=mMaxRefine) coarseCalculateFluxareas(lev);
1955
1956                         int wss = 0, wse = 1;
1957                         // FIXME always currset!?
1958                         wss = wse = mLevel[lev].setCurr;
1959                         for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
1960                                         // warning - check sets for higher levels...?
1961                                 FSGR_FORIJK1(lev) {
1962                                         if( 
1963                                                         (RFLAG(lev,i,j,k, workSet) & CFFluid) || 
1964                                                         (RFLAG(lev,i,j,k, workSet) & CFInter) ||
1965                                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromCoarse) || 
1966                                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromFine) || 
1967                                                         (RFLAG(lev,i,j,k, workSet) & CFGrNorm) 
1968                                                         ) {
1969                                                 // these cells have to be scaled...
1970                                         } else {
1971                                                 continue;
1972                                         }
1973
1974                                         // collide on current set
1975                                         LbmFloat rhoOld;
1976                                         LbmVec velOld;
1977                                         LbmFloat rho, ux,uy,uz;
1978                                         rho=0.0; ux =  uy = uz = 0.0;
1979                                         for(int l=0; l<D::cDfNum; l++) {
1980                                                 LbmFloat m = QCELL(lev, i, j, k, workSet, l); 
1981                                                 rho += m;
1982                                                 ux  += (D::dfDvecX[l]*m);
1983                                                 uy  += (D::dfDvecY[l]*m); 
1984                                                 uz  += (D::dfDvecZ[l]*m); 
1985                                         } 
1986                                         rhoOld = rho;
1987                                         velOld = LbmVec(ux,uy,uz);
1988
1989                                         LbmFloat rhoNew = (rhoOld-rhoAvg)*scaleFac +rhoAvg;
1990                                         LbmVec velNew = velOld * scaleFac;
1991
1992                                         LbmFloat df[LBM_DFNUM];
1993                                         LbmFloat feqOld[LBM_DFNUM];
1994                                         LbmFloat feqNew[LBM_DFNUM];
1995                                         for(int l=0; l<D::cDfNum; l++) {
1996                                                 feqOld[l] = D::getCollideEq(l,rhoOld, velOld[0],velOld[1],velOld[2] );
1997                                                 feqNew[l] = D::getCollideEq(l,rhoNew, velNew[0],velNew[1],velNew[2] );
1998                                                 df[l] = QCELL(lev, i,j,k,workSet, l);
1999                                         }
2000                                         const LbmFloat Qo = D::getLesNoneqTensorCoeff(df,feqOld);
2001                                         const LbmFloat oldOmega = D::getLesOmega(levOldOmega[lev], mLevel[lev].lcsmago,Qo);
2002                                         const LbmFloat newOmega = D::getLesOmega(mLevel[lev].omega,mLevel[lev].lcsmago,Qo);
2003                                         //newOmega = mLevel[lev].omega; // FIXME debug test
2004
2005                                         //LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]);
2006                                         const LbmFloat dfScale = (newSteptime/newOmega)/(levOldStepsize[lev]/oldOmega);
2007                                         //dfScale = dfScaleFac/newOmega;
2008                                         
2009                                         for(int l=0; l<D::cDfNum; l++) {
2010                                                 // org scaling
2011                                                 //df = eqOld + (df-eqOld)*dfScale; df *= (eqNew/eqOld); // non-eq. scaling, important
2012                                                 // new scaling
2013                                                 LbmFloat dfn = feqNew[l] + (df[l]-feqOld[l])*dfScale*feqNew[l]/feqOld[l]; // non-eq. scaling, important
2014                                                 //df = eqNew + (df-eqOld)*dfScale; // modified ig scaling, no real difference?
2015                                                 QCELL(lev, i,j,k,workSet, l) = dfn;
2016                                         }
2017
2018                                         if(RFLAG(lev,i,j,k, workSet) & CFInter) {
2019                                                 //if(workSet==mLevel[lev].setCurr) 
2020                                                 LbmFloat area = 1.0;
2021                                                 if(lev!=mMaxRefine) area = QCELL(lev, i,j,k,workSet, dFlux);
2022                                                 massTOld += QCELL(lev, i,j,k,workSet, dMass) * area;
2023                                                 volTOld += QCELL(lev, i,j,k,workSet, dFfrac);
2024
2025                                                 // wrong... QCELL(i,j,k,workSet, dMass] = (QCELL(i,j,k,workSet, dFfrac]*rhoNew);
2026                                                 QCELL(lev, i,j,k,workSet, dMass) = (QCELL(lev, i,j,k,workSet, dMass)/rhoOld*rhoNew);
2027                                                 QCELL(lev, i,j,k,workSet, dFfrac) = (QCELL(lev, i,j,k,workSet, dMass)/rhoNew);
2028
2029                                                 //if(workSet==mLevel[lev].setCurr) 
2030                                                 massTNew += QCELL(lev, i,j,k,workSet, dMass);
2031                                                 volTNew += QCELL(lev, i,j,k,workSet, dFfrac);
2032                                         }
2033                                         if(RFLAG(lev,i,j,k, workSet) & CFFluid) { // DEBUG
2034                                                 if(RFLAG(lev,i,j,k, workSet) & (CFGrFromFine|CFGrFromCoarse)) { // DEBUG
2035                                                         // dont include 
2036                                                 } else {
2037                                                         LbmFloat area = 1.0;
2038                                                         if(lev!=mMaxRefine) area = QCELL(lev, i,j,k,workSet, dFlux) * mLevel[lev].lcellfactor;
2039                                                         //if(workSet==mLevel[lev].setCurr) 
2040                                                         massTOld += rhoOld*area;
2041                                                         //if(workSet==mLevel[lev].setCurr) 
2042                                                         massTNew += rhoNew*area;
2043                                                         volTOld += area;
2044                                                         volTNew += area;
2045                                                 }
2046                                         }
2047
2048                                 } // IJK
2049                         } // workSet
2050
2051                 } // lev
2052
2053                 if(!D::mSilent) {
2054                         debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE "<<D::mStepCnt<<
2055                                         " no"<<mTimeSwitchCounts<<" maxdt"<<mMaxStepTime<<
2056                                         " mindt"<<mMinStepTime<<" currdt"<<mLevel[mMaxRefine].stepsize, 10);
2057                         debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE  masst:"<<massTNew<<","<<massTOld<<" org:"<<mCurrentMass<<"; "<<
2058                                         " volt:"<<volTNew<<","<<volTOld<<" org:"<<mCurrentVolume, 10);
2059                 } else {
2060                         debMsgStd("\nLbmOptSolver::step",DM_MSG,"Timestep change by "<< (newdt/levOldStepsize[mMaxRefine]) <<" newDt:"<<newdt
2061                                         <<", oldDt:"<<levOldStepsize[mMaxRefine]<<" newOmega:"<<D::mOmega<<" gStar:"<<D::mpParam->getCurrentGStar() , 10);
2062                 }
2063         } // rescale?
2064         
2065         //errMsg("adaptTimestep","Warning - brute force rescale off!"); minCutoff = false; // DEBUG
2066         if(minCutoff) {
2067                 errMsg("adaptTimestep","Warning - performing Brute-Force rescale... (sim:"<<D::mName<<" step:"<<D::mStepCnt<<" newdt="<<desireddt<<" mindt="<<D::mpParam->getMinStepTime()<<") " );
2068                 //brute force resacle all the time?
2069
2070                 for(int lev=mMaxRefine; lev>=0 ; lev--) {
2071                 int rescs=0;
2072                 int wss = 0, wse = 1;
2073 #if COMPRESSGRIDS==1
2074                 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
2075 #endif // COMPRESSGRIDS==1
2076                 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
2077                 //for(int workSet = 0; workSet<=1; workSet++) {
2078                 FSGR_FORIJK1(lev) {
2079
2080                         //if( (RFLAG(lev, i,j,k, workSet) & CFFluid) || (RFLAG(lev, i,j,k, workSet) & CFInter) ) {
2081                         if( 
2082                                         (RFLAG(lev,i,j,k, workSet) & CFFluid) || 
2083                                         (RFLAG(lev,i,j,k, workSet) & CFInter) ||
2084                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromCoarse) || 
2085                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromFine) || 
2086                                         (RFLAG(lev,i,j,k, workSet) & CFGrNorm) 
2087                                         ) {
2088                                 // these cells have to be scaled...
2089                         } else {
2090                                 continue;
2091                         }
2092
2093                         // collide on current set
2094                         LbmFloat rho, ux,uy,uz;
2095                         rho=0.0; ux =  uy = uz = 0.0;
2096                         for(int l=0; l<D::cDfNum; l++) {
2097                                 LbmFloat m = QCELL(lev, i, j, k, workSet, l); 
2098                                 rho += m;
2099                                 ux  += (D::dfDvecX[l]*m);
2100                                 uy  += (D::dfDvecY[l]*m); 
2101                                 uz  += (D::dfDvecZ[l]*m); 
2102                         } 
2103 #ifndef WIN32
2104                         if (!finite(rho)) {
2105                                 errMsg("adaptTimestep","Brute force non-finite rho at"<<PRINT_IJK);  // DEBUG!
2106                                 rho = 1.0;
2107                                 ux = uy = uz = 0.0;
2108                                 QCELL(lev, i, j, k, workSet, dMass) = 1.0;
2109                                 QCELL(lev, i, j, k, workSet, dFfrac) = 1.0;
2110                         }
2111 #endif // WIN32
2112
2113                         if( (ux*ux+uy*uy+uz*uz)> (allowMax*allowMax) ) {
2114                                 LbmFloat cfac = allowMax/sqrt(ux*ux+uy*uy+uz*uz);
2115                                 ux *= cfac;
2116                                 uy *= cfac;
2117                                 uz *= cfac;
2118                                 for(int l=0; l<D::cDfNum; l++) {
2119                                         QCELL(lev, i, j, k, workSet, l) = D::getCollideEq(l, rho, ux,uy,uz); }
2120                                 rescs++;
2121                                 debMsgDirect("B");
2122                         }
2123
2124                 } } 
2125                         //if(rescs>0) { errMsg("adaptTimestep","!!!!! Brute force rescaling was necessary !!!!!!!"); }
2126                         debMsgStd("adaptTimestep",DM_MSG,"Brute force rescale done. level:"<<lev<<" rescs:"<<rescs, 1);
2127                 //TTT mNumProblems += rescs; // add to problem display...
2128                 } // lev,set,ijk
2129
2130         } // try brute force rescale?
2131
2132         // time adap done...
2133 }
2134
2135
2136
2137
2138 /******************************************************************************
2139  * work on lists from updateCellMass to reinit cell flags
2140  *****************************************************************************/
2141
2142 template<class D>
2143 LbmFloat 
2144 LbmFsgrSolver<D>::getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l) {
2145         //return 0.0; // test
2146         int level = mMaxRefine;
2147         LbmFloat *ccel = RACPNT(level, i,j,k, workSet);
2148
2149         LbmFloat nx,ny,nz, nv1,nv2;
2150         if(RFLAG_NB(level,i,j,k,workSet, dE) &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } else nv1 = 0.0;
2151         if(RFLAG_NB(level,i,j,k,workSet, dW) &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } else nv2 = 0.0;
2152         nx = 0.5* (nv2-nv1);
2153         if(RFLAG_NB(level,i,j,k,workSet, dN) &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
2154         if(RFLAG_NB(level,i,j,k,workSet, dS) &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
2155         ny = 0.5* (nv2-nv1);
2156 #if LBMDIM==3
2157         if(RFLAG_NB(level,i,j,k,workSet, dT) &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
2158         if(RFLAG_NB(level,i,j,k,workSet, dB) &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
2159         nz = 0.5* (nv2-nv1);
2160 #else //LBMDIM==3
2161         nz = 0.0;
2162 #endif //LBMDIM==3
2163         LbmFloat scal = mDvecNrm[l][0]*nx + mDvecNrm[l][1]*ny + mDvecNrm[l][2]*nz;
2164
2165         LbmFloat ret = 1.0;
2166         // forward direction, add mass (for filling cells):
2167         if(dirForw) {
2168                 if(scal<LBM_EPSILON) ret = 0.0;
2169                 else ret = scal;
2170         } else {
2171                 // backward for emptying
2172                 if(scal>-LBM_EPSILON) ret = 0.0;
2173                 else ret = scal * -1.0;
2174         }
2175         //errMsg("massd", PRINT_IJK<<" nv"<<nvel<<" : ret="<<ret ); //xit(1); //VECDEB
2176         return ret;
2177 }
2178
2179 template<class D>
2180 void LbmFsgrSolver<D>::addToNewInterList( int ni, int nj, int nk ) {
2181 #if FSGR_STRICT_DEBUG==10
2182         // dangerous, this can change the simulation...
2183   /*for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
2184        iter != mListNewInter.end(); iter++ ) {
2185     if(ni!=iter->x) continue;
2186     if(nj!=iter->y) continue;
2187     if(nk!=iter->z) continue;
2188                 // all 3 values match... skip point
2189                 return;
2190         } */
2191 #endif // FSGR_STRICT_DEBUG==1
2192         // store point
2193         LbmPoint newinter;
2194         newinter.x = ni; newinter.y = nj; newinter.z = nk;
2195         mListNewInter.push_back(newinter);
2196 }
2197
2198 template<class D>
2199 void LbmFsgrSolver<D>::interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet, bool markNbs) {
2200         LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0;
2201         LbmFloat intDf[19] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
2202
2203 #if OPT3D==1 
2204         // for macro add
2205         LbmFloat addDfFacT, addVal, usqr;
2206         LbmFloat *addfcel, *dstcell;
2207         LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM];
2208         LbmFloat lcsmDstOmega, lcsmSrcOmega, lcsmdfscale;
2209 #endif // OPT3D==true 
2210
2211         // SET required nbs to from coarse (this might overwrite flag several times)
2212         // this is not necessary for interpolateFineFromCoarse
2213         if(markNbs) {
2214         FORDF1{ 
2215                 int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2216                 if(RFLAG(lev,ni,nj,nk,dstSet)&CFUnused) {
2217                         // parents have to be inited!
2218                         interpolateCellFromCoarse(lev, ni, nj, nk, dstSet, t, CFFluid|CFGrFromCoarse, false);
2219                 }
2220         } }
2221
2222         // change flag of cell to be interpolated
2223         RFLAG(lev,i,j,k, dstSet) = flagSet;
2224         mNumInterdCells++;
2225
2226         // interpolation lines...
2227         int betx = i&1;
2228         int bety = j&1;
2229         int betz = k&1;
2230         
2231         if((!betx) && (!bety) && (!betz)) {
2232                 ADD_INT_DFS(lev-1, i/2  ,j/2  ,k/2  , 0.0, 1.0);
2233         }
2234         else if(( betx) && (!bety) && (!betz)) {
2235                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D1);
2236                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D1);
2237         }
2238         else if((!betx) && ( bety) && (!betz)) {
2239                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D1);
2240                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D1);
2241         }
2242         else if((!betx) && (!bety) && ( betz)) {
2243                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D1);
2244                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D1);
2245         }
2246         else if(( betx) && ( bety) && (!betz)) {
2247                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D2);
2248                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D2);
2249                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D2);
2250                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)  , t, WO1D2);
2251         }
2252         else if((!betx) && ( bety) && ( betz)) {
2253                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D2);
2254                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D2);
2255                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D2);
2256                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)+1, t, WO1D2);
2257         }
2258         else if(( betx) && (!bety) && ( betz)) {
2259                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D2);
2260                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D2);
2261                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D2);
2262                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)+1, t, WO1D2);
2263         }
2264         else if(( betx) && ( bety) && ( betz)) {
2265                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D3);
2266                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D3);
2267                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D3);
2268                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)+1, t, WO1D3);
2269                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D3);
2270                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)  , t, WO1D3);
2271                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)+1, t, WO1D3);
2272                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)+1, t, WO1D3);
2273         }
2274         else {
2275                 D::mPanic=1;
2276                 errFatal("interpolateCellFromCoarse","Invalid!?", SIMWORLD_GENERICERROR);       
2277         }
2278
2279         IDF_WRITEBACK;
2280         return;
2281 }
2282
2283 template<class D>
2284 void LbmFsgrSolver<D>::reinitFlags( int workSet )
2285 {
2286         // OLD mods:
2287         // add all to intel list?
2288         // check ffrac for new cells
2289         // new if cell inits (last loop)
2290         // vweights handling
2291
2292 #if ELBEEM_BLENDER==1
2293         const int debugFlagreinit = 0;
2294 #else // ELBEEM_BLENDER==1
2295         const int debugFlagreinit = 0;
2296 #endif // ELBEEM_BLENDER==1
2297         
2298         // some things need to be read/modified on the other set
2299         int otherSet = (workSet^1);
2300         // fixed level on which to perform 
2301         int workLev = mMaxRefine;
2302
2303   /* modify interface cells from lists */
2304   /* mark filled interface cells as fluid, or emptied as empty */
2305         /* count neighbors and distribute excess mass to interface neighbor cells
2306    * problems arise when there are no interface neighbors anymore
2307          * then just distribute to any fluid neighbors...
2308          */
2309
2310         // for symmetry, first init all neighbor cells */
2311         for( vector<LbmPoint>::iterator iter=mListFull.begin();
2312        iter != mListFull.end(); iter++ ) {
2313     int i=iter->x, j=iter->y, k=iter->z;
2314                 if(debugFlagreinit) errMsg("FULL", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
2315     FORDF1 {
2316                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2317       if( RFLAG(workLev, ni,nj,nk, workSet) & CFEmpty ){
2318                                 // new and empty interface cell, dont change old flag here!
2319                                 addToNewInterList(ni,nj,nk);
2320                                 
2321                                 // preinit speed, get from average surrounding cells
2322                                 // interpolate from non-workset to workset, sets are handled in function
2323                                 {
2324                                         // WARNING - other i,j,k than filling cell!
2325                                         int ei=ni; int ej=nj; int ek=nk;
2326                                         LbmFloat avgrho = 0.0;
2327                                         LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
2328                                         LbmFloat cellcnt = 0.0;
2329                                         LbmFloat avgnbdf[LBM_DFNUM];
2330                                         FORDF0M { avgnbdf[m]= 0.0; }
2331
2332                                         for(int nbl=1; nbl< D::cDfNum ; ++nbl) {
2333                                                 if( (RFLAG_NB(workLev,ei,ej,ek,workSet,nbl) & CFFluid) || 
2334                                                         ((!(RFLAG_NB(workLev,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) &&
2335                                                                 (RFLAG_NB(workLev,ei,ej,ek,workSet,nbl) & CFInter) )) { 
2336                                                         cellcnt += 1.0;
2337                                         for(int rl=0; rl< D::cDfNum ; ++rl) { 
2338                                                                 LbmFloat nbdf =  QCELL_NB(workLev,ei,ej,ek, workSet,nbl, rl);
2339                                                                 avgnbdf[rl] += nbdf;
2340                                                                 avgux  += (D::dfDvecX[rl]*nbdf); 
2341                                                                 avguy  += (D::dfDvecY[rl]*nbdf);  
2342                                                                 avguz  += (D::dfDvecZ[rl]*nbdf);  
2343                                                                 avgrho += nbdf;
2344                                                         }
2345                                                 }
2346                                         }
2347
2348                                         if(cellcnt<=0.0) {
2349                                                 // no nbs? just use eq.
2350                                                 //FORDF0 { QCELL(workLev,ei,ej,ek, workSet, l) = D::dfEquil[l]; }
2351                                                 avgrho = 1.0;
2352                                                 avgux = avguy = avguz = 0.0;
2353                                                 //TTT mNumProblems++;
2354 #if ELBEEM_BLENDER!=1
2355                                                 D::mPanic=1; errFatal("NYI2","cellcnt<=0.0",SIMWORLD_GENERICERROR);
2356 #endif // ELBEEM_BLENDER
2357                                         } else {
2358                                                 // init speed
2359                                                 avgux /= cellcnt; avguy /= cellcnt; avguz /= cellcnt;
2360                                                 avgrho /= cellcnt;
2361                                                 FORDF0M { avgnbdf[m] /= cellcnt; } // CHECK FIXME test?
2362                                         }
2363
2364                                         // careful with l's...
2365                                         FORDF0M { 
2366                                                 QCELL(workLev,ei,ej,ek, workSet, m) = D::getCollideEq( m,avgrho,  avgux, avguy, avguz ); 
2367                                                 //QCELL(workLev,ei,ej,ek, workSet, l) = avgnbdf[l]; // CHECK FIXME test?
2368                                         }
2369                                         //errMsg("FNEW", PRINT_VEC(ei,ej,ek)<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<<avgrho<<" vel"<<PRINT_VEC(avgux,avguy,avguz) ); // DEBUG SYMM
2370                                         QCELL(workLev,ei,ej,ek, workSet, dMass) = 0.0; //?? new
2371                                         QCELL(workLev,ei,ej,ek, workSet, dFfrac) = 0.0; //?? new
2372                                         //RFLAG(workLev,ei,ej,ek,workSet) = (CellFlagType)(CFInter|CFNoInterpolSrc);
2373                                         changeFlag(workLev,ei,ej,ek,workSet, (CFInter|CFNoInterpolSrc));
2374                                         if(debugFlagreinit) errMsg("NEWE", PRINT_IJK<<" newif "<<PRINT_VEC(ei,ej,ek)<<" rho"<<avgrho<<" vel("<<avgux<<","<<avguy<<","<<avguz<<") " );
2375                                 } 
2376       }
2377                         /* prevent surrounding interface cells from getting removed as empty cells 
2378                          * (also cells that are not newly inited) */
2379       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
2380                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete);
2381                                 changeFlag(workLev,ni,nj,nk, workSet, (RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete));
2382                                 // also add to list...
2383                                 addToNewInterList(ni,nj,nk);
2384                         } // NEW?
2385     }
2386
2387                 // NEW? no extra loop...
2388                 //RFLAG(workLev,i,j,k, workSet) = CFFluid;
2389                 changeFlag(workLev,i,j,k, workSet,CFFluid);
2390         }
2391
2392         /* remove empty interface cells that are not allowed to be removed anyway
2393          * this is important, otherwise the dreaded cell-type-flickering can occur! */
2394   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
2395        iter != mListEmpty.end(); iter++ ) {
2396     int i=iter->x, j=iter->y, k=iter->z;
2397                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) {
2398                         // remove entry
2399                         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
2400                         iter = mListEmpty.erase(iter); 
2401                         iter--; // and continue with next...
2402
2403                         // treat as "new inter"
2404                         addToNewInterList(i,j,k);
2405                 }
2406         } 
2407
2408
2409         /* problems arise when adjacent cells empty&fill ->
2410                  let fill cells+surrounding interface cells have the higher importance */
2411   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
2412        iter != mListEmpty.end(); iter++ ) {
2413     int i=iter->x, j=iter->y, k=iter->z;
2414                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)){ errMsg("A"," ARGHARGRAG "); } // DEBUG
2415                 if(debugFlagreinit) errMsg("EMPT", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) );
2416
2417                 /* set surrounding fluid cells to interface cells */
2418     FORDF1 {
2419                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2420       if( RFLAG(workLev,ni,nj,nk, workSet) & CFFluid){
2421                                 // init fluid->interface 
2422                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(CFInter); 
2423                                 changeFlag(workLev,ni,nj,nk, workSet, CFInter); 
2424                                 /* new mass = current density */
2425                                 LbmFloat nbrho = QCELL(workLev,ni,nj,nk, workSet, dC);
2426                 for(int rl=1; rl< D::cDfNum ; ++rl) { nbrho += QCELL(workLev,ni,nj,nk, workSet, rl); }
2427                                 QCELL(workLev,ni,nj,nk, workSet, dMass) =  nbrho; 
2428                                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) =  1.0; 
2429
2430                                 // store point
2431                                 addToNewInterList(ni,nj,nk);
2432       }
2433       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter){
2434                                 // test, also add to list...
2435                                 addToNewInterList(ni,nj,nk);
2436                         } // NEW?
2437     }
2438
2439                 /* for symmetry, set our flag right now */
2440                 //RFLAG(workLev,i,j,k, workSet) = CFEmpty;
2441                 changeFlag(workLev,i,j,k, workSet, CFEmpty);
2442                 // mark cell not be changed mass... - not necessary, not in list anymore anyway!
2443         } // emptylist
2444
2445
2446         
2447         // precompute weights to get rid of order dependancies
2448         vector<lbmFloatSet> vWeights;
2449         vWeights.reserve( mListFull.size() + mListEmpty.size() );
2450         int weightIndex = 0;
2451   int nbCount = 0;
2452         LbmFloat nbWeights[LBM_DFNUM];
2453         LbmFloat nbTotWeights = 0.0;
2454         for( vector<LbmPoint>::iterator iter=mListFull.begin();
2455        iter != mListFull.end(); iter++ ) {
2456     int i=iter->x, j=iter->y, k=iter->z;
2457     nbCount = 0; nbTotWeights = 0.0;
2458     FORDF1 {
2459                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2460       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
2461                                 nbCount++;
2462                                 nbWeights[l] = getMassdWeight(1,i,j,k,workSet,l);
2463                                 nbTotWeights += nbWeights[l];
2464       } else {
2465                                 nbWeights[l] = -100.0; // DEBUG;
2466                         }
2467     }
2468                 if(nbCount>0) { 
2469                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
2470         vWeights[weightIndex].val[0] = nbTotWeights;
2471         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
2472         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
2473                 } else { 
2474         vWeights[weightIndex].numNbs = 0.0;
2475                 }
2476                 weightIndex++;
2477         }
2478   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
2479        iter != mListEmpty.end(); iter++ ) {
2480     int i=iter->x, j=iter->y, k=iter->z;
2481     nbCount = 0; nbTotWeights = 0.0;
2482     FORDF1 {
2483                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2484       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
2485                                 nbCount++;
2486                                 nbWeights[l] = getMassdWeight(0,i,j,k,workSet,l);
2487                                 nbTotWeights += nbWeights[l];
2488       } else {
2489                                 nbWeights[l] = -100.0; // DEBUG;
2490                         }
2491     }
2492                 if(nbCount>0) { 
2493                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
2494         vWeights[weightIndex].val[0] = nbTotWeights;
2495         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
2496         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
2497                 } else { 
2498         vWeights[weightIndex].numNbs = 0.0;
2499                 }
2500                 weightIndex++;
2501         } 
2502         weightIndex = 0;
2503         
2504
2505         /* process full list entries, filled cells are done after this loop */
2506         for( vector<LbmPoint>::iterator iter=mListFull.begin();
2507        iter != mListFull.end(); iter++ ) {
2508     int i=iter->x, j=iter->y, k=iter->z;
2509
2510                 LbmFloat myrho = QCELL(workLev,i,j,k, workSet, dC);
2511     FORDF1 { myrho += QCELL(workLev,i,j,k, workSet, l); } // QCELL.rho
2512
2513     LbmFloat massChange = QCELL(workLev,i,j,k, workSet, dMass) - myrho;
2514     /*int nbCount = 0;
2515                 LbmFloat nbWeights[LBM_DFNUM];
2516     FORDF1 {
2517                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2518       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
2519                                 nbCount++;
2520                                 nbWeights[l] = vWeights[weightIndex].val[l];
2521       } else {
2522                         }
2523     }*/
2524
2525                 //errMsg("FDIST", PRINT_IJK<<" mss"<<massChange <<" nb"<< nbCount ); // DEBUG SYMM
2526                 if(vWeights[weightIndex].numNbs>0.0) {
2527                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
2528                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
2529                         FORDF1 {
2530                                 int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2531         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
2532                                         LbmFloat change = -1.0;
2533                                         if(nbTotWeightsp>0.0) {
2534                                                 //change = massChange * ( nbWeights[l]/nbTotWeightsp );
2535                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
2536                                         } else {
2537                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
2538                                         }
2539                                         QCELL(workLev,ni,nj,nk, workSet, dMass) += change;
2540                                 }
2541                         }
2542                         massChange = 0.0;
2543                 } else {
2544                         // Problem! no interface neighbors
2545                         D::mFixMass += massChange;
2546                         //TTT mNumProblems++;
2547                         //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<<D::mFixMass);
2548                 }
2549                 weightIndex++;
2550
2551     // already done? RFLAG(workLev,i,j,k, workSet) = CFFluid;
2552     QCELL(workLev,i,j,k, workSet, dMass) = myrho; // should be rho... but unused?
2553     QCELL(workLev,i,j,k, workSet, dFfrac) = 1.0; // should be rho... but unused?
2554     /*QCELL(workLev,i,j,k, otherSet, dMass) = myrho; // NEW?
2555     QCELL(workLev,i,j,k, otherSet, dFfrac) = 1.0; // NEW? COMPRT */
2556   } // fulllist
2557
2558
2559         /* now, finally handle the empty cells - order is important, has to be after
2560          * full cell handling */
2561   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
2562        iter != mListEmpty.end(); iter++ ) {
2563     int i=iter->x, j=iter->y, k=iter->z;
2564
2565     LbmFloat massChange = QCELL(workLev, i,j,k, workSet, dMass);
2566     /*int nbCount = 0;
2567                 LbmFloat nbWeights[LBM_DFNUM];
2568     FORDF1 {
2569                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2570       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
2571                                 nbCount++;
2572                                 nbWeights[l] = vWeights[weightIndex].val[l];
2573       } else {
2574                                 nbWeights[l] = -100.0; // DEBUG;
2575                         }
2576     }*/
2577
2578                 //errMsg("EDIST", PRINT_IJK<<" mss"<<massChange <<" nb"<< nbCount ); // DEBUG SYMM
2579                 //if(nbCount>0) {
2580                 if(vWeights[weightIndex].numNbs>0.0) {
2581                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
2582                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
2583                         FORDF1 {
2584                                 int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2585         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
2586                                         LbmFloat change = -1.0;
2587                                         if(nbTotWeightsp>0.0) {
2588                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
2589                                         } else {
2590                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
2591                                         }
2592                                         QCELL(workLev, ni,nj,nk, workSet, dMass) += change;
2593                                 }
2594                         }
2595                         massChange = 0.0;
2596                 } else {
2597                         // Problem! no interface neighbors
2598                         D::mFixMass += massChange;
2599                         //TTT mNumProblems++;
2600                         //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<D::mFixMass);
2601                 }
2602                 weightIndex++;
2603                 
2604                 // finally... make it empty 
2605     // already done? RFLAG(workLev,i,j,k, workSet) = CFEmpty;
2606     QCELL(workLev,i,j,k, workSet, dMass) = 0.0;
2607     QCELL(workLev,i,j,k, workSet, dFfrac) = 0.0;
2608         }
2609   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
2610        iter != mListEmpty.end(); iter++ ) {
2611     int i=iter->x, j=iter->y, k=iter->z;
2612     //RFLAG(workLev,i,j,k, otherSet) = CFEmpty;
2613     changeFlag(workLev,i,j,k, otherSet, CFEmpty);
2614     /*QCELL(workLev,i,j,k, otherSet, dMass) = 0.0;
2615     QCELL(workLev,i,j,k, otherSet, dFfrac) = 0.0; // COMPRT OFF */
2616         } 
2617
2618
2619         // check if some of the new interface cells can be removed again 
2620         // never happens !!! not necessary
2621         // calculate ffrac for new IF cells NEW
2622
2623         // how many are really new interface cells?
2624         int numNewIf = 0;
2625   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
2626        iter != mListNewInter.end(); iter++ ) {
2627     int i=iter->x, j=iter->y, k=iter->z;
2628                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
2629                         continue; 
2630                         // FIXME remove from list?
2631                 }
2632                 numNewIf++;
2633         }
2634
2635         // redistribute mass, reinit flags
2636         float newIfFac = 1.0/(LbmFloat)numNewIf;
2637   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
2638        iter != mListNewInter.end(); iter++ ) {
2639     int i=iter->x, j=iter->y, k=iter->z;
2640                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
2641                         //errMsg("???"," "<<PRINT_IJK);
2642                         continue; 
2643                 }
2644
2645     QCELL(workLev,i,j,k, workSet, dMass) += (D::mFixMass * newIfFac);
2646
2647                 int nbored = 0;
2648                 FORDF1 { nbored |= RFLAG_NB(workLev, i,j,k, workSet,l); }
2649                 if((nbored & CFBnd)==0) { RFLAG(workLev,i,j,k, workSet) |= CFNoBndFluid; }
2650                 if((nbored & CFFluid)==0) { RFLAG(workLev,i,j,k, workSet) |= CFNoNbFluid; }
2651                 if((nbored & CFEmpty)==0) { RFLAG(workLev,i,j,k, workSet) |= CFNoNbEmpty; }
2652
2653                 if(!(RFLAG(workLev,i,j,k, otherSet)&CFInter)) {
2654                         RFLAG(workLev,i,j,k, workSet) = (CellFlagType)(RFLAG(workLev,i,j,k, workSet) | CFNoDelete);
2655                 }
2656                 if(debugFlagreinit) errMsg("NEWIF", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" f"<< RFLAG(workLev,i,j,k, workSet)<<" wl"<<workLev );
2657         }
2658
2659         // reinit fill fraction
2660   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
2661        iter != mListNewInter.end(); iter++ ) {
2662     int i=iter->x, j=iter->y, k=iter->z;
2663                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { continue; }
2664
2665                 LbmFloat nrho = 0.0;
2666                 FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); }
2667     QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho;
2668