elbeem update:
[blender.git] / intern / elbeem / intern / solver_util.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 #include "particletracer.h"
13
14 // MPT
15 #include "ntl_world.h"
16 #include "simulation_object.h"
17
18
19 /******************************************************************************
20  * helper functions
21  *****************************************************************************/
22
23 // try to enhance surface?
24 #define SURFACE_ENH 2
25
26 //! for raytracing
27 void LbmFsgrSolver::prepareVisualization( void ) {
28         int lev = mMaxRefine;
29         int workSet = mLevel[lev].setCurr;
30
31         int mainGravDir=0;
32         LbmFloat mainGravLen = 0.;
33         FORDF1{
34                 LbmFloat thisGravLen = dot(LbmVec(dfVecX[l],dfVecY[l],dfVecZ[l]), getNormalized(mLevel[mMaxRefine].gravity) );
35                 if(thisGravLen>mainGravLen) {
36                         mainGravLen = thisGravLen;
37                         mainGravDir = l;
38                 }
39                 //errMsg("GRAVT","l"<<l<<" dfv"<<LbmVec(dfVecX[l],dfVecY[l],dfVecZ[l])<<" g"<< mLevel[mMaxRefine].gravity<< " mgl"<<mainGravLen<<" mgd"<<mainGravDir);
40         }
41
42         //make same prepareVisualization and getIsoSurface...
43 #if LBMDIM==2
44         // 2d, place in the middle of isofield slice (k=2)
45 #  define ZKD1 0
46         // 2d z offset = 2, lbmGetData adds 1, so use one here
47 #  define ZKOFF 1
48         // reset all values...
49         for(int k= 0; k< 5; ++k) 
50    for(int j=0;j<mLevel[lev].lSizey-0;j++) 
51     for(int i=0;i<mLevel[lev].lSizex-0;i++) {
52                 *this->mpIso->lbmGetData(i,j,ZKOFF)=0.0;
53         }
54 #else // LBMDIM==2
55         // 3d, use normal bounds
56 #  define ZKD1 1
57 #  define ZKOFF k
58         // reset all values...
59         for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
60    for(int j=0;j<mLevel[lev].lSizey-0;j++) 
61     for(int i=0;i<mLevel[lev].lSizex-0;i++) {
62                 *this->mpIso->lbmGetData(i,j,ZKOFF)=0.0;
63         }
64 #endif // LBMDIM==2
65
66         // MPT, ignore
67         //if( strstr( this->getName().c_str(), "mpfluid2" ) != NULL) {
68                 //errMsg("MPT TESTX","skipping mpfluid2");
69                 //return;
70         //}
71         if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
72                 mpIso->resetAll(0.);
73         }
74
75
76         LbmFloat minval = this->mIsoValue*1.05; // / mIsoWeight[13]; 
77         // add up...
78         float val = 0.0;
79         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) 
80    for(int j=1;j<mLevel[lev].lSizey-1;j++) 
81     for(int i=1;i<mLevel[lev].lSizex-1;i++) {
82                         const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
83                         //if(cflag&(CFBnd|CFEmpty)) {
84
85 #if SURFACE_ENH==0
86
87                         // no enhancements...
88                         if( (cflag&(CFFluid|CFUnused)) ) {
89                                 val = 1.;
90                         } else if( (cflag&CFInter) ) {
91                                 val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
92                         } else {
93                                 continue;
94                         }
95
96 #else // SURFACE_ENH!=1
97                         if(cflag&CFBnd) {
98                                 // treated in second loop
99                                 continue;
100                         } else if(cflag&CFUnused) {
101                                 val = 1.;
102                         } else if( (cflag&CFFluid) && (cflag&CFNoBndFluid)) {
103                                 // optimized fluid
104                                 val = 1.;
105                         } else if( (cflag&(CFEmpty|CFInter|CFFluid)) ) {
106                                 int noslipbnd = 0;
107                                 int intercnt = 0;
108                                 FORDF1 { 
109                                         const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
110                                         if(nbflag&CFInter){ intercnt++; }
111
112                                         if(l!=mainGravDir) continue; // only check bnd along main grav. dir
113                                         //if((nbflag&CFBnd)&&(nbflag&CFBndNoslip)){ noslipbnd=1; }
114                                         if((nbflag&CFBnd)){ noslipbnd=1; }
115                                 }
116
117                                 if(cflag&CFEmpty) {
118                                         // special empty treatment
119                                         if((noslipbnd)&&(intercnt>6)) {
120                                                 *this->mpIso->lbmGetData(i,j,ZKOFF) += minval;
121                                         } else if((noslipbnd)&&(intercnt>0)) {
122                                                 // necessary?
123                                                 *this->mpIso->lbmGetData(i,j,ZKOFF) += this->mIsoValue*0.9;
124                                         } else {
125                                                 // nothing to do...
126                                         }
127
128                                         continue;
129                                 } else if(cflag&(CFNoNbEmpty|CFFluid)) {
130                                         // no empty nb interface cells are treated as full
131                                         val=1.0;
132                                 } else {
133                                         val = (QCELL(lev, i,j,k,workSet, dFfrac)); 
134                                 }
135                                 
136                                 if(noslipbnd) {
137                                         if(val<minval) val = minval; 
138                                         *this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); 
139                                 }
140 #endif // SURFACE_ENH>0
141
142                         } else { // all others, unused?
143                                 continue;
144                         } 
145
146                         *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); 
147                         *this->mpIso->lbmGetData( i   , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[1] ); 
148                         *this->mpIso->lbmGetData( i+1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[2] ); 
149
150                         *this->mpIso->lbmGetData( i-1 , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[3] ); 
151                         *this->mpIso->lbmGetData( i   , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[4] ); 
152                         *this->mpIso->lbmGetData( i+1 , j   ,ZKOFF-ZKD1) += ( val * mIsoWeight[5] ); 
153
154                         *this->mpIso->lbmGetData( i-1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[6] ); 
155                         *this->mpIso->lbmGetData( i   , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[7] ); 
156                         *this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[8] ); 
157
158
159                         *this->mpIso->lbmGetData( i-1 , j-1  ,ZKOFF  ) += ( val * mIsoWeight[9] ); 
160                         *this->mpIso->lbmGetData( i   , j-1  ,ZKOFF  ) += ( val * mIsoWeight[10] ); 
161                         *this->mpIso->lbmGetData( i+1 , j-1  ,ZKOFF  ) += ( val * mIsoWeight[11] ); 
162
163                         *this->mpIso->lbmGetData( i-1 , j    ,ZKOFF  ) += ( val * mIsoWeight[12] ); 
164                         *this->mpIso->lbmGetData( i   , j    ,ZKOFF  ) += ( val * mIsoWeight[13] ); 
165                         *this->mpIso->lbmGetData( i+1 , j    ,ZKOFF  ) += ( val * mIsoWeight[14] ); 
166
167                         *this->mpIso->lbmGetData( i-1 , j+1  ,ZKOFF  ) += ( val * mIsoWeight[15] ); 
168                         *this->mpIso->lbmGetData( i   , j+1  ,ZKOFF  ) += ( val * mIsoWeight[16] ); 
169                         *this->mpIso->lbmGetData( i+1 , j+1  ,ZKOFF  ) += ( val * mIsoWeight[17] ); 
170
171
172                         *this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[18] ); 
173                         *this->mpIso->lbmGetData( i   , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[19] ); 
174                         *this->mpIso->lbmGetData( i+1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[20] ); 
175
176                         *this->mpIso->lbmGetData( i-1 , j   ,ZKOFF+ZKD1) += ( val * mIsoWeight[21] ); 
177                         *this->mpIso->lbmGetData( i   , j   ,ZKOFF+ZKD1)+= ( val * mIsoWeight[22] ); 
178                         *this->mpIso->lbmGetData( i+1 , j   ,ZKOFF+ZKD1) += ( val * mIsoWeight[23] ); 
179
180                         *this->mpIso->lbmGetData( i-1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[24] ); 
181                         *this->mpIso->lbmGetData( i   , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[25] ); 
182                         *this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); 
183         }
184
185         // TEST!?
186 #if SURFACE_ENH>=2
187
188         if(mFsSurfGenSetting&fssgNoObs) {
189                 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) 
190                  for(int j=1;j<mLevel[lev].lSizey-1;j++) 
191                         for(int i=1;i<mLevel[lev].lSizex-1;i++) {
192                                 const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
193                                 if(cflag&(CFBnd)) {
194                                         CellFlagType nbored=0;
195                                         LbmFloat avgfill=0.,avgfcnt=0.;
196                                         FORDF1 { 
197                                                 const int ni = i+this->dfVecX[l];
198                                                 const int nj = j+this->dfVecY[l];
199                                                 const int nk = ZKOFF+this->dfVecZ[l];
200
201                                                 const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet);
202                                                 nbored |= nbflag;
203                                                 if(nbflag&CFInter) {
204                                                         avgfill += QCELL(lev, ni,nj,nk, workSet,dFfrac); avgfcnt += 1.;
205                                                 } else if(nbflag&CFFluid) {
206                                                         avgfill += 1.; avgfcnt += 1.;
207                                                 } else if(nbflag&CFEmpty) {
208                                                         avgfill += 0.; avgfcnt += 1.;
209                                                 }
210
211                                                 //if( (ni<0) || (nj<0) || (nk<0) 
212                                                  //|| (ni>=mLevel[mMaxRefine].lSizex) 
213                                                  //|| (nj>=mLevel[mMaxRefine].lSizey) 
214                                                  //|| (nk>=mLevel[mMaxRefine].lSizez) ) continue;
215                                         }
216
217                                         if(nbored&CFInter) {
218                                                 if(avgfcnt>0.) avgfill/=avgfcnt;
219                                                 *this->mpIso->lbmGetData(i,j,ZKOFF) = avgfill; continue;
220                                         } 
221                                         else if(nbored&CFFluid) {
222                                                 *this->mpIso->lbmGetData(i,j,ZKOFF) = 1.; continue;
223                                         }
224
225                                 }
226                         }
227
228                 // move surface towards inner "row" of obstacle
229                 // cells if necessary (all obs cells without fluid/inter
230                 // nbs (=iso==0) next to obstacles...)
231                 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) 
232                         for(int j=1;j<mLevel[lev].lSizey-1;j++) 
233                                 for(int i=1;i<mLevel[lev].lSizex-1;i++) {
234                                         const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
235                                         if( (cflag&(CFBnd)) && (*this->mpIso->lbmGetData(i,j,ZKOFF)==0.)) {
236                                                 int bndnbcnt=0;
237                                                 FORDF1 { 
238                                                         const int ni = i+this->dfVecX[l];
239                                                         const int nj = j+this->dfVecY[l];
240                                                         const int nk = ZKOFF+this->dfVecZ[l];
241                                                         const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet);
242                                                         if(nbflag&CFBnd) bndnbcnt++;
243                                                 }
244                                                 if(bndnbcnt>0) *this->mpIso->lbmGetData(i,j,ZKOFF)=this->mIsoValue*0.95; 
245                                         }
246                                 }
247         }
248         // */
249
250         if(mFsSurfGenSetting&fssgNoNorth) 
251                 for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
252                         for(int j=0;j<mLevel[lev].lSizey-0;j++) {
253                                 *this->mpIso->lbmGetData(0,                   j,ZKOFF) = *this->mpIso->lbmGetData(1,                   j,ZKOFF);
254                         }
255         if(mFsSurfGenSetting&fssgNoEast) 
256                 for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
257                         for(int i=0;i<mLevel[lev].lSizex-0;i++) {
258                                 *this->mpIso->lbmGetData(i,0,                   ZKOFF) = *this->mpIso->lbmGetData(i,1,                   ZKOFF);
259                         }
260         if(mFsSurfGenSetting&fssgNoSouth) 
261                 for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
262                         for(int j=0;j<mLevel[lev].lSizey-0;j++) {
263                                 *this->mpIso->lbmGetData(mLevel[lev].lSizex-1,j,ZKOFF) = *this->mpIso->lbmGetData(mLevel[lev].lSizex-2,j,ZKOFF);
264                         }
265         if(mFsSurfGenSetting&fssgNoWest) 
266                 for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 
267                         for(int i=0;i<mLevel[lev].lSizex-0;i++) {
268                                 *this->mpIso->lbmGetData(i,mLevel[lev].lSizey-1,ZKOFF) = *this->mpIso->lbmGetData(i,mLevel[lev].lSizey-2,ZKOFF);
269                         }
270         if(LBMDIM>2) {
271                 if(mFsSurfGenSetting&fssgNoBottom) 
272                         for(int j=0;j<mLevel[lev].lSizey-0;j++) 
273                                 for(int i=0;i<mLevel[lev].lSizex-0;i++) {
274                                         *this->mpIso->lbmGetData(i,j,0                   ) = *this->mpIso->lbmGetData(i,j,1                   );
275                                 } 
276                 if(mFsSurfGenSetting&fssgNoTop) 
277                         for(int j=0;j<mLevel[lev].lSizey-0;j++) 
278                                 for(int i=0;i<mLevel[lev].lSizex-0;i++) {
279                                         *this->mpIso->lbmGetData(i,j,mLevel[lev].lSizez-1) = *this->mpIso->lbmGetData(i,j,mLevel[lev].lSizez-2);
280                                 } 
281         }
282 #endif // SURFACE_ENH>=2
283
284
285         // update preview, remove 2d?
286         if((this->mOutputSurfacePreview)&&(LBMDIM==3)) {
287                 int pvsx = (int)(this->mPreviewFactor*this->mSizex);
288                 int pvsy = (int)(this->mPreviewFactor*this->mSizey);
289                 int pvsz = (int)(this->mPreviewFactor*this->mSizez);
290                 //float scale = (float)this->mSizex / previewSize;
291                 LbmFloat scalex = (LbmFloat)this->mSizex/(LbmFloat)pvsx;
292                 LbmFloat scaley = (LbmFloat)this->mSizey/(LbmFloat)pvsy;
293                 LbmFloat scalez = (LbmFloat)this->mSizez/(LbmFloat)pvsz;
294                 for(int k= 0; k< (pvsz-1); ++k) 
295                 for(int j=0;j< pvsy;j++) 
296                 for(int i=0;i< pvsx;i++) {
297                                         *mpPreviewSurface->lbmGetData(i,j,k) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley), (int)(k*scalez) );
298                                 }
299                 // set borders again...
300                 for(int k= 0; k< (pvsz-1); ++k) {
301                         for(int j=0;j< pvsy;j++) {
302                                 *mpPreviewSurface->lbmGetData(0,j,k) = *this->mpIso->lbmGetData( 0, (int)(j*scaley), (int)(k*scalez) );
303                                 *mpPreviewSurface->lbmGetData(pvsx-1,j,k) = *this->mpIso->lbmGetData( this->mSizex-1, (int)(j*scaley), (int)(k*scalez) );
304                         }
305                         for(int i=0;i< pvsx;i++) {
306                                 *mpPreviewSurface->lbmGetData(i,0,k) = *this->mpIso->lbmGetData( (int)(i*scalex), 0, (int)(k*scalez) );
307                                 *mpPreviewSurface->lbmGetData(i,pvsy-1,k) = *this->mpIso->lbmGetData( (int)(i*scalex), this->mSizey-1, (int)(k*scalez) );
308                         }
309                 }
310                 for(int j=0;j<pvsy;j++)
311                         for(int i=0;i<pvsx;i++) {      
312                                 *mpPreviewSurface->lbmGetData(i,j,0) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , 0);
313                                 *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = *this->mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , this->mSizez-1);
314                         }
315
316                 if(mFarFieldSize>=2.) {
317                         // also remove preview border
318                         for(int k= 0; k< (pvsz-1); ++k) {
319                                 for(int j=0;j< pvsy;j++) {
320                                         *mpPreviewSurface->lbmGetData(0,j,k) = 
321                                         *mpPreviewSurface->lbmGetData(1,j,k) =  
322                                         *mpPreviewSurface->lbmGetData(2,j,k);
323                                         *mpPreviewSurface->lbmGetData(pvsx-1,j,k) = 
324                                         *mpPreviewSurface->lbmGetData(pvsx-2,j,k) = 
325                                         *mpPreviewSurface->lbmGetData(pvsx-3,j,k);
326                                         //0.0;
327                                 }
328                                 for(int i=0;i< pvsx;i++) {
329                                         *mpPreviewSurface->lbmGetData(i,0,k) = 
330                                         *mpPreviewSurface->lbmGetData(i,1,k) = 
331                                         *mpPreviewSurface->lbmGetData(i,2,k);
332                                         *mpPreviewSurface->lbmGetData(i,pvsy-1,k) = 
333                                         *mpPreviewSurface->lbmGetData(i,pvsy-2,k) = 
334                                         *mpPreviewSurface->lbmGetData(i,pvsy-3,k);
335                                         //0.0;
336                                 }
337                         }
338                         for(int j=0;j<pvsy;j++)
339                                 for(int i=0;i<pvsx;i++) {      
340                                         *mpPreviewSurface->lbmGetData(i,j,0) = 
341                                         *mpPreviewSurface->lbmGetData(i,j,1) = 
342                                         *mpPreviewSurface->lbmGetData(i,j,2);
343                                         *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = 
344                                         *mpPreviewSurface->lbmGetData(i,j,pvsz-2) = 
345                                         *mpPreviewSurface->lbmGetData(i,j,pvsz-3);
346                                         //0.0;
347                         }
348                 }
349         }
350
351         // MPT
352 #if LBM_INCLUDE_TESTSOLVERS==1
353         if( ( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) && (mInitDone) ) {
354                 errMsg("MPT TESTX","special mpfluid1");
355                 vector<SimulationObject*> *sims = mpGlob->getSims();
356                 for(int simi=0; simi<(int)(sims->size()); simi++) {
357                         SimulationObject *sim = (*sims)[simi];
358                         if( strstr( sim->getName().c_str(), "mpfluid2" ) != NULL) {
359                                 LbmFsgrSolver *fsgr = (LbmFsgrSolver *)(sim->getSolver());
360                                 int msyLev = mMaxRefine;
361                                 int simLev = fsgr->mMaxRefine;
362
363                                 IsoSurface* simIso = fsgr->mpIso;
364                                 int idst = mLevel[msyLev].lSizex-2 -1; // start at boundary
365                                 int isrc = 0                    +2;
366                                 if((simIso)&&(fsgr->mInitDone)) {
367                                 fsgr->prepareVisualization();
368
369                                 for(int k=0;k<=mLevel[simLev].lSizez-1;++k) {
370                                         for(int j=0;j<=mLevel[simLev].lSizey-1;++j) {
371                                                 for(int i=isrc; i<mLevel[simLev].lSizex-1; i++) {
372                                                         //LbmFloat *cceldst = &QCELL(          msyLev ,idst+i,j,k, msySet,0);
373                                                         //LbmFloat *ccelsrc = &SIM_QCELL(fsgr, simLev ,isrc+i,j,k, simSet,0);
374                                                         //for(int l=0; l<dTotalNum; l++) {
375                                                                 //RAC(cceldst,l) = RAC(ccelsrc,l);
376                                                         //}
377                                                         // both flag sets, make sure curr/other are same!?
378                                                         //RFLAG(msyLev ,idst+i,j,k, 0) = SIM_RFLAG(fsgr, simLev ,isrc+i,j,k, 0);
379                                                         //RFLAG(msyLev ,idst+i,j,k, 1) = SIM_RFLAG(fsgr, simLev ,isrc+i,j,k, 1);
380                                                         errMsg("TESTX"," "<<PRINT_VEC(idst+i,j,k)<<" < "<<PRINT_VEC(i,j,k) );
381                                                         *mpIso->lbmGetData(idst+i,j,k) = *simIso->lbmGetData(i,j,k);
382                                                 }
383                                         }
384                                 }
385
386                                 } // simIso
387                         }
388                 }
389         }
390 #endif // LBM_INCLUDE_TESTSOLVERS==1
391
392         return;
393 }
394
395 /*! calculate speeds of fluid objects (or inflow) */
396 void LbmFsgrSolver::recalculateObjectSpeeds() {
397         const bool debugRecalc = false;
398         int numobjs = (int)(this->mpGiObjects->size());
399         // note - (numobjs + 1) is entry for domain settings
400
401         if(debugRecalc) errMsg("recalculateObjectSpeeds","start, #obj:"<<numobjs);
402         if(numobjs>255-1) {
403                 errFatal("LbmFsgrSolver::recalculateObjectSpeeds","More than 256 objects currently not supported...",SIMWORLD_INITERROR);
404                 return;
405         }
406         mObjectSpeeds.resize(numobjs+1);
407         for(int i=0; i<(int)(numobjs+0); i++) {
408                 mObjectSpeeds[i] = vec2L(this->mpParam->calculateLattVelocityFromRw( vec2P( (*this->mpGiObjects)[i]->getInitialVelocity(mSimulationTime) )));
409                 if(debugRecalc) errMsg("recalculateObjectSpeeds","id"<<i<<" set to "<< mObjectSpeeds[i]<<", unscaled:"<< (*this->mpGiObjects)[i]->getInitialVelocity(mSimulationTime) );
410         }
411
412         // also reinit part slip values here
413         mObjectPartslips.resize(numobjs+1);
414         for(int i=0; i<=(int)(numobjs+0); i++) {
415                 if(i<numobjs) {
416                         mObjectPartslips[i] = (LbmFloat)(*this->mpGiObjects)[i]->getGeoPartSlipValue();
417                 } else {
418                         // domain setting
419                         mObjectPartslips[i] = this->mDomainPartSlipValue;
420                 }
421                 LbmFloat set = mObjectPartslips[i];
422
423                 // as in setInfluenceVelocity
424                 const LbmFloat dt = mLevel[mMaxRefine].timestep;
425                 const LbmFloat dtInter = 0.01;
426                 //LbmFloat facFv = 1.-set;
427                 // mLevel[mMaxRefine].timestep
428                 LbmFloat facNv = (LbmFloat)( 1.-pow( (double)(set), (double)(dt/dtInter)) );
429                 errMsg("mObjectPartslips","id:"<<i<<"/"<<numobjs<<"  ts:"<<dt<< " its:"<<(dt/dtInter) <<" set"<<set<<" nv"<<facNv<<" test:"<<
430                          pow( (double)(1.-facNv),(double)(dtInter/dt))  );
431                 mObjectPartslips[i] = facNv;
432
433                 if(debugRecalc) errMsg("recalculateObjectSpeeds","id"<<i<<" parts "<< mObjectPartslips[i] );
434         }
435
436         if(debugRecalc) errMsg("recalculateObjectSpeeds","done, domain:"<<mObjectPartslips[numobjs]<<" n"<<numobjs);
437 }
438
439
440
441 /*****************************************************************************/
442 /*! debug object display */
443 /*****************************************************************************/
444 vector<ntlGeometryObject*> LbmFsgrSolver::getDebugObjects() { 
445         vector<ntlGeometryObject*> debo; 
446         if(this->mOutputSurfacePreview) {
447                 debo.push_back( mpPreviewSurface );
448         }
449 #if LBM_INCLUDE_TESTSOLVERS==1
450         if(mUseTestdata) {
451                 vector<ntlGeometryObject*> tdebo; 
452                 tdebo = mpTest->getDebugObjects();
453                 for(size_t i=0; i<tdebo.size(); i++) debo.push_back( tdebo[i] );
454         }
455 #endif // ELBEEM_PLUGIN
456         return debo; 
457 }
458
459 /******************************************************************************
460  * particle handling
461  *****************************************************************************/
462
463 /*! init particle positions */
464 int LbmFsgrSolver::initParticles() { 
465   int workSet = mLevel[mMaxRefine].setCurr;
466   int tries = 0;
467   int num = 0;
468         ParticleTracer *partt = mpParticles;
469
470   partt->setStart( this->mvGeoStart + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) );
471   partt->setEnd  ( this->mvGeoEnd   + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) );
472
473   partt->setSimStart( ntlVec3Gfx(0.0) );
474   partt->setSimEnd  ( ntlVec3Gfx(this->mSizex,   this->mSizey,   getForZMaxBnd(mMaxRefine)) );
475   
476   while( (num<partt->getNumInitialParticles()) && (tries<100*partt->getNumInitialParticles()) ) {
477     LbmFloat x,y,z,t;
478     x = 1.0+(( (LbmFloat)(this->mSizex-3.) )          * (rand()/(RAND_MAX+1.0)) );
479     y = 1.0+(( (LbmFloat)(this->mSizey-3.) )          * (rand()/(RAND_MAX+1.0)) );
480     z = 1.0+(( (LbmFloat) getForZMax1(mMaxRefine)-2. )* (rand()/(RAND_MAX+1.0)) );
481     int i = (int)(x+0.5);
482     int j = (int)(y+0.5);
483     int k = (int)(z+0.5);
484     if(LBMDIM==2) {
485       k = 0; z = 0.5; // place in the middle of domain
486     }
487
488     //if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid) ) 
489     //&& TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFNoNbFluid ) 
490     //if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter|CFMbndInflow) ) { 
491     if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFNoBndFluid|CFUnused) ) { 
492                         bool cellOk = true;
493                         //? FORDF1 { if(!(RFLAG_NB(mMaxRefine,i,j,k,workSet, l) & CFFluid)) cellOk = false; }
494                         if(!cellOk) continue;
495       // in fluid...
496       partt->addParticle(x,y,z);
497                         partt->getLast()->setStatus(PART_IN);
498                         partt->getLast()->setType(PART_TRACER);
499
500                         partt->getLast()->setSize(1.);
501                         // randomize size
502                         partt->getLast()->setSize(0.5 + (rand()/(RAND_MAX+1.0)));
503
504                         if( ( partt->getInitStart()>0.)
505                                         && ( partt->getInitEnd()>0.)
506                                         && ( partt->getInitEnd()>partt->getInitStart() )) {
507                 t = partt->getInitStart()+ (partt->getInitEnd()-partt->getInitStart())*(rand()/(RAND_MAX+1.0));
508                                 partt->getLast()->setLifeTime( -t );
509                         }
510       num++;
511     }
512     tries++;
513   } // */
514
515         /*FSGR_FORIJK1(mMaxRefine) {
516                 if( (RFLAG(mMaxRefine,i,j,k, workSet) & (CFNoBndFluid)) ) {
517         LbmFloat rndn = (rand()/(RAND_MAX+1.0));
518                         if(rndn>0.0) {
519                                 ntlVec3Gfx pos( (LbmFloat)(i)-0.5, (LbmFloat)(j)-0.5, (LbmFloat)(k)-0.5 );
520                                 if(LBMDIM==2) { pos[2]=0.5; }
521                                 partt->addParticle( pos[0],pos[1],pos[2] );
522                                 partt->getLast()->setStatus(PART_IN);
523                                 partt->getLast()->setType(PART_TRACER);
524                                 partt->getLast()->setSize(1.0);
525                         }
526                 }
527         } // */
528
529
530         // DEBUG TEST
531 #if LBM_INCLUDE_TESTSOLVERS==1
532         if(mUseTestdata) { 
533                 const bool partDebug=false;
534                 if(mpTest->mDebugvalue2!=0.0){ errMsg("LbmTestdata"," part init "<<mpTest->mDebugvalue2); }
535                 if(mpTest->mDebugvalue2==-12.0){ 
536                         const int lev = mMaxRefine;
537                         for(int i=5;i<15;i++) {
538                                 LbmFloat x,y,z;
539                                 y = 0.5+(LbmFloat)(i);
540                                 x = mLevel[lev].lSizex/20.0*10.0;
541                                 z = mLevel[lev].lSizez/20.0*2.0;
542                                 partt->addParticle(x,y,z);
543                                 partt->getLast()->setStatus(PART_IN);
544                                 partt->getLast()->setType(PART_BUBBLE);
545                                 partt->getLast()->setSize(  (-4.0+(LbmFloat)i)/1.0  );
546                                 if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() );
547                         }
548                 }
549                 if(mpTest->mDebugvalue2==-11.0){ 
550                         const int lev = mMaxRefine;
551                         for(int i=5;i<15;i++) {
552                                 LbmFloat x,y,z;
553                                 y = 10.5+(LbmFloat)(i);
554                                 x = mLevel[lev].lSizex/20.0*10.0;
555                                 z = mLevel[lev].lSizez/20.0*40.0;
556                                 partt->addParticle(x,y,z);
557                                 partt->getLast()->setStatus(PART_IN);
558                                 partt->getLast()->setType(PART_DROP);
559                                 partt->getLast()->setSize(  (-4.0+(LbmFloat)i)/1.0  );
560                                 if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() );
561                         }
562                 }
563                 // place floats on rectangular region FLOAT_JITTER_BND
564                 if(mpTest->mDebugvalue2==-10.0){ 
565                         const int lev = mMaxRefine;
566                         const int sx = mLevel[lev].lSizex;
567                         const int sy = mLevel[lev].lSizey;
568                         //for(int j=-(int)(sy*0.25);j<-(int)(sy*0.25)+2;++j) { for(int i=-(int)(sx*0.25);i<-(int)(sy*0.25)+2;++i) {
569                         //for(int j=-(int)(sy*1.25);j<(int)(2.25*sy);++j) { for(int i=-(int)(sx*1.25);i<(int)(2.25*sx);++i) {
570                         for(int j=-(int)(sy*0.3);j<(int)(1.3*sy);++j) { for(int i=-(int)(sx*0.3);i<(int)(1.3*sx);++i) {
571                         //for(int j=-(int)(sy*0.2);j<(int)(0.2*sy);++j) { for(int i= (int)(sx*0.5);i<= (int)(0.51*sx);++i) {
572                                         LbmFloat x,y,z;
573                                         x = 0.0+(LbmFloat)(i);
574                                         y = 0.0+(LbmFloat)(j);
575                                         //z = mLevel[lev].lSizez/10.0*2.5 - 1.0;
576                                         z = mLevel[lev].lSizez/20.0*9.5 - 1.0;
577                                         //z = mLevel[lev].lSizez/20.0*4.5 - 1.0;
578                                         partt->addParticle(x,y,z);
579                                         //if( (i>0)&&(i<sx) && (j>0)&&(j<sy) ) { partt->getLast()->setStatus(PART_IN); } else { partt->getLast()->setStatus(PART_OUT); }
580                                         partt->getLast()->setStatus(PART_IN);
581                                         partt->getLast()->setType(PART_FLOAT);
582                                         partt->getLast()->setSize( 15.0 );
583                                         if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() );
584                          }
585                 }       }
586         } 
587         // DEBUG TEST
588 #endif // LBM_INCLUDE_TESTSOLVERS
589
590         
591   debMsgStd("LbmFsgrSolver::initParticles",DM_MSG,"Added "<<num<<" particles, genProb:"<<this->mPartGenProb<<", tries:"<<tries, 10);
592   if(num != partt->getNumParticles()) return 1;
593
594         return 0;
595 }
596
597 #define P_CHANGETYPE(p, newtype) \
598                 p->setLifeTime(0.); \
599     /* errMsg("PIT","U pit"<<(p)->getId()<<" pos:"<< (p)->getPos()<<" status:"<<convertFlags2String((p)->getFlags())<<" to "<< (newtype) ); */ \
600                 p->setType(newtype); 
601
602 // tracer defines
603 #define TRACE_JITTER 0.025
604 #define TRACE_RAND (rand()/(RAND_MAX+1.0))*TRACE_JITTER-(TRACE_JITTER*0.5)
605 #define FFGET_NORM(var,dl) \
606                                                         if(RFLAG_NB(lev,i,j,k,workSet, dl) &(CFInter)){ (var) = QCELL_NB(lev,i,j,k,workSet,dl,dFfrac); } \
607                                                         else if(RFLAG_NB(lev,i,j,k,workSet, dl) &(CFFluid|CFUnused)){ (var) = 1.; } else (var) = 0.0;
608
609 // float jitter
610 #define FLOAT_JITTER_BND (FLOAT_JITTER*2.0)
611 #define FLOAT_JITTBNDRAND(x) ((rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND*(1.-(x/(LbmFloat)maxdw))-(FLOAT_JITTER_BND*(1.-(x)/(LbmFloat)maxdw)*0.5)) 
612
613
614 void LbmFsgrSolver::advanceParticles() { 
615   int workSet = mLevel[mMaxRefine].setCurr;
616         LbmFloat vx=0.0,vy=0.0,vz=0.0;
617         //LbmFloat df[27]; //feq[27];
618 #define DEL_PART { \
619         /*errMsg("PIT","DEL AT "<< __LINE__<<" type:"<<p->getType()<<"  ");  */ \
620         p->setActive( false ); \
621         continue; }
622
623         myTime_t parttstart = getTime(); 
624         const LbmFloat cellsize = this->mpParam->getCellSize();
625         const LbmFloat timestep = this->mpParam->getTimestep();
626         //const LbmFloat viscAir = 1.79 * 1e-5; // RW2L kin. viscosity, mu
627         const LbmFloat viscWater = 1.0 * 1e-6; // RW2L kin. viscosity, mu
628         const LbmFloat rhoAir = 1.2;  // [kg m^-3] RW2L
629         const LbmFloat rhoWater = 1000.0; // RW2L
630         const LbmFloat minDropSize = 0.0005; // [m], = 2mm  RW2L
631         const LbmVec   velAir(0.); // [m / s]
632
633         const LbmFloat r1 = 0.005;  // r max
634         const LbmFloat r2 = 0.0005; // r min
635         const LbmFloat v1 = 9.0; // v max
636         const LbmFloat v2 = 2.0; // v min
637         const LbmVec rwgrav = vec2L( this->mpParam->getGravity(mSimulationTime) );
638         const bool useff = (mFarFieldSize>2.); // if(mpTest->mDebugvalue1>0.0){ 
639
640         // TODO use timestep size
641         const int cutval = mCutoff; // use full border!?
642         int actCnt=0;
643         if(this->mStepCnt%50==49) { mpParticles->cleanup(); }
644   for(vector<ParticleObject>::iterator pit= mpParticles->getParticlesBegin();
645       pit!= mpParticles->getParticlesEnd(); pit++) {
646     //errMsg("PIT"," pit"<<(*pit).getId()<<" pos:"<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())<<" vel:"<< (*pit).getVel()  );
647     //errMsg("PIT"," pit pos:"<< (*pit)->getPos()<<" vel:"<< (*pit)->getVel()<<" status:"<<convertFlags2String((*pit)->getFlags()) <<" " <<mpParticles->getStart()<<" "<<mpParticles->getEnd() );
648                 //int flag = (*pit).getFlags();
649     if( (*pit).getActive()==false ) continue;
650                 // skip until reached
651                 ParticleObject *p = &(*pit);
652                 if(p->getLifeTime()<0.){ 
653                         if(p->getLifeTime()<-this->mSimulationTime) continue; 
654                         else p->setLifeTime(-mLevel[mMaxRefine].timestep); // zero for following update
655                 }
656     int i,j,k;
657                 p->setLifeTime(p->getLifeTime()+mLevel[mMaxRefine].timestep);
658
659                 // nearest neighbor, particle positions don't include empty bounds
660                 ntlVec3Gfx pos = p->getPos();
661                 i= (int)(pos[0]+0.5);
662                 j= (int)(pos[1]+0.5);
663                 k= (int)(pos[2]+0.5);
664                 if(LBMDIM==2) { k = 0; }
665
666                 // only testdata handling, all for sws
667 #if LBM_INCLUDE_TESTSOLVERS==1
668                 if(useff && (mpTest->mDebugvalue1>0.)) {
669                         p->setStatus(PART_OUT);
670                         mpTest->handleParticle(p, i,j,k); continue;
671                 } 
672 #endif // LBM_INCLUDE_TESTSOLVERS==1
673
674                 // FIXME , add k tests again, remove per type ones...
675                 if(p->getStatus()&PART_IN) { // IN
676                         if( (i<cutval)||(i>this->mSizex-1-cutval)||
677                                         (j<cutval)||(j>this->mSizey-1-cutval)
678                                         //||(k<cutval)||(k>this->mSizez-1-cutval) 
679                                         ) {
680                                 if(!useff) { DEL_PART;
681                                 } else { 
682                                         p->setStatus(PART_OUT); 
683                                         /* del? */ //if((rand()/(RAND_MAX+1.0))<0.5) DEL_PART;
684                                 }
685                         } 
686                 } else { // OUT rough check
687                         // check in again?
688                         if( (i>=cutval)&&(i<=this->mSizex-1-cutval)&&
689                                         (j>=cutval)&&(j<=this->mSizey-1-cutval)
690                                         //&&(k>=cutval)&&(k<=this->mSizez-1-cutval) 
691                                         ) {
692                                 p->setStatus(PART_IN);
693                                 /* del? */ //if((rand()/(RAND_MAX+1.0))<0.5) DEL_PART;
694                         }
695                 }
696                 //p->setStatus(PART_OUT);// DEBUG always out!
697
698                 if( (p->getType()==PART_BUBBLE) ||
699                     (p->getType()==PART_TRACER) ) {
700
701                         // no interpol
702                         vx = vy = vz = 0.0;
703                         if(p->getStatus()&PART_IN) { // IN
704                                 if(k>=cutval) {
705                                         if(k>this->mSizez-1-cutval) DEL_PART; 
706
707                                         if( RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFUnused) ) {
708                                                 // still ok
709                                                 int partLev = mMaxRefine;
710                                                 int si=i, sj=j, sk=k;
711                                                 while(partLev>0 && RFLAG(partLev, si,sj,sk, workSet)&(CFUnused)) {
712                                                         partLev--;
713                                                         si/=2;
714                                                         sj/=2;
715                                                         sk/=2;
716                                                 }
717                                                 //LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
718                                                 LbmFloat *ccel = RACPNT(partLev, si,sj,sk, workSet);
719                                                 FORDF1{
720                                                         //LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
721                                                         LbmFloat cdf = RAC(ccel, l);
722                                                         // TODO update below
723                                                         //df[l] = cdf;
724                                                         vx  += (this->dfDvecX[l]*cdf); 
725                                                         vy  += (this->dfDvecY[l]*cdf);  
726                                                         vz  += (this->dfDvecZ[l]*cdf);  
727                                                 }
728
729                                                 // remove gravity influence
730                                                 const LbmFloat lesomega = mLevel[mMaxRefine].omega; // no les
731                                                 vx -= mLevel[mMaxRefine].gravity[0] * lesomega*0.5;
732                                                 vy -= mLevel[mMaxRefine].gravity[1] * lesomega*0.5;
733                                                 vz -= mLevel[mMaxRefine].gravity[2] * lesomega*0.5;
734
735                                         } else { // OUT
736                                                 // out of bounds, deactivate...
737                                                 // FIXME make fsgr treatment
738                                                 if(p->getType()==PART_BUBBLE) { P_CHANGETYPE(p, PART_FLOAT ); continue; }
739                                         }
740                                 } else {
741                                         // below 3d region, just rise
742                                 }
743                         } else { // OUT
744 #                               if LBM_INCLUDE_TESTSOLVERS==1
745                                 if(useff) { mpTest->handleParticle(p, i,j,k); }
746                                 else DEL_PART;
747 #                               else // LBM_INCLUDE_TESTSOLVERS==1
748                                 DEL_PART;
749 #                               endif // LBM_INCLUDE_TESTSOLVERS==1
750                                 // TODO use x,y vel...?
751                         }
752
753                         ntlVec3Gfx v = p->getVel(); // dampen...
754                         if( (useff)&& (p->getType()==PART_BUBBLE) ) {
755                                 // test rise
756                                 //O vz = p->getVel()[2]-0.5*mLevel[mMaxRefine].gravity[2];
757
758                                 LbmFloat radius = p->getSize() * minDropSize;
759                                 LbmVec   velPart = vec2L(p->getVel()) *cellsize/timestep; // L2RW, lattice velocity
760                                 LbmVec   velWater = LbmVec(vx,vy,vz) *cellsize/timestep;// L2RW, fluid velocity
761                                 LbmVec   velRel = velWater - velPart;
762                                 LbmFloat velRelNorm = norm(velRel);
763                                 // TODO calculate values in lattice units, compute CD?!??!
764                                 LbmFloat pvolume = rhoAir * 4.0/3.0 * M_PI* radius*radius*radius; // volume: 4/3 pi r^3
765                                 //const LbmFloat cd = 
766
767                                 LbmVec fb = -rwgrav* pvolume *rhoWater;
768                                 LbmVec fd = velRel*6.0*M_PI*radius* (1e-3); //viscWater;
769                                 LbmVec change = (fb+fd) *10.0*timestep  *(timestep/cellsize);
770                                 //LbmVec change = (fb+fd) *timestep / (pvolume*rhoAir)  *(timestep/cellsize);
771                                 //actCnt++; // should be after active test
772                                 if(actCnt<0) {
773                                         errMsg("PIT","BTEST1   vol="<<pvolume<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel);
774                                         errMsg("PIT","BTEST2        cellsize="<<cellsize<<" timestep="<<timestep<<" viscW="<<viscWater<<" ss/mb="<<(timestep/(pvolume*rhoAir)));
775                                         errMsg("PIT","BTEST2        grav="<<rwgrav<<"  " );
776                                         errMsg("PIT","BTEST2        change="<<(change)<<" fb="<<(fb)<<" fd="<<(fd)<<" ");
777                                         errMsg("PIT","BTEST2        change="<<norm(change)<<" fb="<<norm(fb)<<" fd="<<norm(fd)<<" ");
778                                 }
779                                         
780                                 //v += change;
781                                 //v += ntlVec3Gfx(vx,vy,vz);
782                                 LbmVec fd2 = (LbmVec(vx,vy,vz)-vec2L(p->getVel())) * 6.0*M_PI*radius* (1e-3); //viscWater;
783                                 LbmFloat w = 0.99;
784                                 vz = (1.0-w)*vz + w*(p->getVel()[2]-0.5*(p->getSize()/5.0)*mLevel[mMaxRefine].gravity[2]);
785                                 v = ntlVec3Gfx(vx,vy,vz)+vec2G(fd2);
786                         } else if(p->getType()==PART_TRACER) {
787                                 v = ntlVec3Gfx(vx,vy,vz);
788                                 CellFlagType fflag = RFLAG(mMaxRefine, i,j,k, workSet);
789
790                                 if(fflag&(CFFluid|CFInter) ) { p->setInFluid(true);
791                                 } else { p->setInFluid(false); }
792
793                                 if( (( fflag&CFFluid ) && ( fflag&CFNoBndFluid )) ||
794                                                 (( fflag&CFInter ) && (!(fflag&CFNoNbFluid)) ) ) {
795                                         // only real fluid
796 #                                       if LBMDIM==3
797                                         p->advance( TRACE_RAND,TRACE_RAND,TRACE_RAND);
798 #                                       else
799                                         p->advance( TRACE_RAND,TRACE_RAND, 0.);
800 #                                       endif
801
802                                 } else {
803                                         // move inwards along normal, make sure normal is valid first
804                                         const int lev = mMaxRefine;
805                                         LbmFloat nx=0.,ny=0.,nz=0., nv1,nv2;
806                                         bool nonorm = false;
807                                         if(i<=0)              { nx = -1.; nonorm = true; }
808                                         if(i>=this->mSizex-1) { nx =  1.; nonorm = true; }
809                                         if(j<=0)              { ny = -1.; nonorm = true; }
810                                         if(j>=this->mSizey-1) { ny =  1.; nonorm = true; }
811 #                                       if LBMDIM==3
812                                         if(k<=0)              { nz = -1.; nonorm = true; }
813                                         if(k>=this->mSizez-1) { nz =  1.; nonorm = true; }
814 #                                       endif // LBMDIM==3
815                                         if(!nonorm) {
816                                                 FFGET_NORM(nv1,dE); FFGET_NORM(nv2,dW);
817                                                 nx = 0.5* (nv2-nv1);
818                                                 FFGET_NORM(nv1,dN); FFGET_NORM(nv2,dS);
819                                                 ny = 0.5* (nv2-nv1);
820 #                                               if LBMDIM==3
821                                                 FFGET_NORM(nv1,dT); FFGET_NORM(nv2,dB);
822                                                 nz = 0.5* (nv2-nv1);
823 #                                               else // LBMDIM==3
824                                                 nz = 0.;
825 #                                               endif // LBMDIM==3
826                                         } else {
827                                                 v = p->getVel() + vec2G(mLevel[mMaxRefine].gravity);
828                                         }
829
830                                         p->advanceVec( (ntlVec3Gfx(nx,ny,nz)) * -0.1 ); // + vec2G(mLevel[mMaxRefine].gravity);
831                                         //if(normNoSqrt(v)<LBM_EPSILON) v = mLevel[mMaxRefine].gravity;
832                                 }
833                         }
834
835                         p->setVel( v );
836                         p->advanceVel();
837                         //errMsg("PPPP"," pos"<<p->getPos()<<" "<<p->getVel() );
838                 } 
839
840                 // drop handling
841                 else if(p->getType()==PART_DROP) {
842                         ntlVec3Gfx v = p->getVel(); // dampen...
843
844                         if(1) {
845                                 LbmFloat radius = p->getSize() * minDropSize;
846                                 LbmVec   velPart = vec2L(p->getVel()) *cellsize /timestep; // * cellsize / timestep; // L2RW, lattice velocity
847                                 LbmVec   velRel = velAir - velPart;
848                                 //LbmVec   velRelLat = velRel /cellsize*timestep; // L2RW
849                                 LbmFloat velRelNorm = norm(velRel);
850                                 // TODO calculate values in lattice units, compute CD?!??!
851                                 LbmFloat mb = rhoWater * 4.0/3.0 * M_PI* radius*radius*radius; // mass: 4/3 pi r^3 rho
852                                 const LbmFloat rw = (r1-radius)/(r1-r2);
853                                 const LbmFloat rmax = (0.5 + 0.5*rw);
854                                 const LbmFloat vmax = (v2 + (v1-v2)* (1.0-rw) );
855                                 const LbmFloat cd = (rmax) * (velRelNorm)/(vmax);
856
857                                 LbmVec fg = rwgrav * mb;//  * (1.0-rhoAir/rhoWater);
858                                 LbmVec fd = velRel* velRelNorm* cd*M_PI *rhoAir *0.5 *radius*radius;
859                                 LbmVec change = (fg+   fd ) *timestep / mb  *(timestep/cellsize);
860
861                                 //actCnt++; // should be after active test
862                                 if(actCnt<0) {
863                                         errMsg("\nPIT","NTEST1   mb="<<mb<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel<<" pgetVel="<<p->getVel() );
864                                         //errMsg("PIT","NTEST2        cellsize="<<cellsize<<" timestep="<<timestep<<" viscAir="<<viscAir<<" ss/mb="<<(timestep/mb));
865                                         //errMsg("PIT","NTEST2        grav="<<rwgrav<<" mb="<<mb<<" "<<" cd="<<cd );
866                                         //errMsg("PIT","NTEST2        change="<<norm(change)<<" fg="<<norm(fg)<<" fd="<<norm(fd)<<" ");
867                                 }
868
869                                 v += vec2G(change);
870                                 p->setVel(v); 
871                                 // NEW
872                         } else {
873                                 p->setVel( v * 0.999 ); // dampen...
874                                 p->setVel( v ); // DEBUG!
875                                 p->addToVel( vec2G(mLevel[mMaxRefine].gravity) );\
876                         } // OLD
877                         p->advanceVel();
878
879                         if(p->getStatus()&PART_IN) { // IN
880                                 if(k<cutval) { DEL_PART; continue; }
881                                 if(k<=this->mSizez-1-cutval){ 
882                                         //if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter)) {
883                                         if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter|CFBnd)) {
884                                                 // still ok
885                                         // shipt3 } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter) ){
886                                         } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFUnused) ){
887                                                 // FIXME make fsgr treatment
888                                                 P_CHANGETYPE(p, PART_FLOAT ); continue; 
889                                                 // jitter in cell to prevent stacking when hitting a steep surface
890                                                 ntlVec3Gfx cpos = p->getPos(); cpos[0] += (rand()/(RAND_MAX+1.0))-0.5;
891                                                 cpos[1] += (rand()/(RAND_MAX+1.0))-0.5; p->setPos(cpos);
892                                                 //} else DEL_PART;
893                                         } else {
894                                                 DEL_PART;
895                                                 this->mNumParticlesLost++;
896                                         }
897                                 }
898                         } else { // OUT
899 #                               if LBM_INCLUDE_TESTSOLVERS==1
900                                 if(useff) { mpTest->handleParticle(p, i,j,k); }
901                                 else{ DEL_PART; }
902 #                               else // LBM_INCLUDE_TESTSOLVERS==1
903                                 { DEL_PART; }
904 #                               endif // LBM_INCLUDE_TESTSOLVERS==1
905                         }
906
907                 } // air particle
908
909                 // inter particle
910                 else if(p->getType()==PART_INTER) {
911                         if(p->getStatus()&PART_IN) { // IN
912                                 if((k<cutval)||(k>this->mSizez-1-cutval)) {
913                                         // undecided particle above or below... remove?
914                                         DEL_PART; 
915                                 }
916
917                                 if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) {
918                                         // still ok
919                                 } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid|CFUnused) ) ) {
920                         //errMsg("PIT","NEWBUB pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
921
922                                         //P_CHANGETYPE(p, PART_BUBBLE ); continue;
923                                         // currently bubbles off! DEBUG!
924                                         //DEL_PART; // DEBUG bubbles off for siggraph
925                                         P_CHANGETYPE(p, PART_FLOAT ); continue;
926                                 } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) {
927                         //errMsg("PIT","NEWDROP pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
928                                         P_CHANGETYPE(p, PART_DROP ); continue;
929                                 } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFBnd ) ) {
930                                         // 060423 new test
931                                         P_CHANGETYPE(p, PART_FLOAT ); continue;
932                                 }
933                         } else { // OUT
934                                 // undecided particle outside... remove?
935                                 DEL_PART; 
936                                 //P_CHANGETYPE(p, PART_FLOAT ); continue;
937                         }
938                 }
939
940                 // float particle
941                 else if(p->getType()==PART_FLOAT) {
942                         //  test - delte on boundary!?
943                         //if( (i<=cutval)||(i>=this->mSizex-1-cutval)|| (j<=cutval)||(j>=this->mSizey-1-cutval)) { DEL_PART; } // DEBUG TEST
944
945 #                               if LBM_INCLUDE_TESTSOLVERS==1
946                         /*LbmFloat prob = 1.0;
947                         // vanishing
948                         prob = (rand()/(RAND_MAX+1.0));
949                         // increse del prob. up to max height by given factor
950                         const int fhCheckh = (int)(mpTest->mFluidHeight*1.25);
951                         const LbmFloat fhProbh = 25.;
952                         if((useff)&&(k>fhCheckh)) {
953                                 LbmFloat fac = (LbmFloat)(k-fhCheckh)/(LbmFloat)(fhProbh*(mLevel[mMaxRefine].lSizez-fhCheckh));
954                                 prob /= fac; 
955                         }
956                         //if(prob<mLevel[mMaxRefine].timestep*0.1) DEL_PART;
957                         // forced vanishing
958                         //? if(k>this->mSizez*3/4) {    if(prob<3.0*mLevel[mMaxRefine].timestep*0.1) DEL_PART;}
959                         // */
960
961 #                               endif // LBM_INCLUDE_TESTSOLVERS==1
962
963                         if(p->getStatus()&PART_IN) { // IN
964                                 //if((k<cutval)||(k>this->mSizez-1-cutval)) DEL_PART; 
965                                 if(k<cutval) DEL_PART; 
966                                 if(k>this->mSizez-1-cutval) { P_CHANGETYPE(p, PART_DROP ); continue; } // create new drop
967                                 //ntlVec3Gfx v = getVelocityAt(i,j,k);
968                                 vx = vy = vz = 0.0;
969
970                                 // define from particletracer.h
971 #if MOVE_FLOATS==1
972                                 const int DEPTH_AVG=3; // only average interface vels
973                                 int ccnt=0;
974                                 for(int kk=0;kk<DEPTH_AVG;kk+=1) {
975                                         if((k-kk)<1) continue;
976                                         if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFInter)) {} else continue;
977                                         ccnt++;
978                                         FORDF1{
979                                                 LbmFloat cdf = QCELL(mMaxRefine, i,j,k-kk, workSet, l);
980                                                 //df[l] = cdf;
981                                                 vx  += (this->dfDvecX[l]*cdf); 
982                                                 vy  += (this->dfDvecY[l]*cdf);  
983                                                 vz  += (this->dfDvecZ[l]*cdf);  
984                                         }
985                                 }
986                                 if(ccnt) {
987                                 // use halved surface velocity (todo, use omega instead)
988                                 vx /=(LbmFloat)(ccnt * 2.0); // half xy speed! value2
989                                 vy /=(LbmFloat)(ccnt * 2.0);
990                                 vz /=(LbmFloat)(ccnt); }
991 #else // MOVE_FLOATS==1
992                                 vx=vy=0.; //p->setVel(ntlVec3Gfx(0.) ); // static_float
993 #endif // MOVE_FLOATS==1
994                                 vx += (rand()/(RAND_MAX+1.0))*(FLOAT_JITTER*0.2)-(FLOAT_JITTER*0.2*0.5);
995                                 vy += (rand()/(RAND_MAX+1.0))*(FLOAT_JITTER*0.2)-(FLOAT_JITTER*0.2*0.5);
996
997                                 //bool delfloat = false;
998                                 if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid|CFUnused) ) ) {
999                                         // in fluid cell
1000                                         /*if((1) && (k<this->mSizez-3) && 
1001                                                         (
1002                                                           TESTFLAG( RFLAG(mMaxRefine, i,j,k+1, workSet), CFInter ) ||
1003                                                           TESTFLAG( RFLAG(mMaxRefine, i,j,k+2, workSet), CFInter ) 
1004                                                         ) ) {
1005                                                 vz = p->getVel()[2]-0.5*mLevel[mMaxRefine].gravity[2];
1006                                                 if(vz<0.0) vz=0.0;
1007                                         } else delfloat=true;
1008                                         // keep below obstacles
1009                                         if((delfloat) && (TESTFLAG( RFLAG(mMaxRefine, i,j,k+1, workSet), CFBnd )) ) {
1010                                                 delfloat=false; vz=0.0;
1011                                         } // */
1012                                         vz = p->getVel()[2]-1.0*mLevel[mMaxRefine].gravity[2]; // simply rise...
1013                                         if(vz<0.) vz=0.;
1014                                 } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFBnd ) ) {
1015                                         //vz = p->getVel()[2]-0.2; // force downwards movement, move below obstacle...
1016                                         vz = p->getVel()[2]+1.0*mLevel[mMaxRefine].gravity[2]; // fall...
1017                                         if(vz>0.) vz=0.;
1018                                 } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) {
1019                                         // keep in interface , one grid cell offset is added in part. gen
1020                                 } else { //if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) { // shipt?, was all else before
1021                                         vz = p->getVel()[2]+1.0*mLevel[mMaxRefine].gravity[2]; // fall...
1022                                         if(vz>0.) vz=0.;
1023                                         // check if above inter, remove otherwise
1024                                         /*if((1) && (k>2) && (
1025                                                                 TESTFLAG( RFLAG(mMaxRefine, i,j,k-1, workSet), CFInter ) ||
1026                                                                 TESTFLAG( RFLAG(mMaxRefine, i,j,k-2, workSet), CFInter ) 
1027                                                                 ) ) {
1028                                                 vz = p->getVel()[2]+0.5*mLevel[mMaxRefine].gravity[2];
1029                                                 if(vz>0.0) vz=0.0;
1030                                         } else delfloat=true; // */
1031                                         //P_CHANGETYPE(p, PART_DROP ); continue; // create new drop
1032                                 }
1033                                 //if(delfloat) DEL_PART;
1034                                 /*
1035                                 // move down from empty
1036                                 else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFEmpty ) ) {
1037                                         vz = p->getVel()[2]+0.5*mLevel[mMaxRefine].gravity[2];
1038                                         if(vz>0.0) vz=0.0;
1039                                         //DEL_PART; // ????
1040                                 } else  {        DEL_PART; } // */
1041                                 //vz = 0.0; // DEBUG
1042
1043                                 p->setVel( vec2G( ntlVec3Gfx(vx,vy,vz) ) ); //?
1044                                 //p->setVel( vec2G(v)*0.75 + p->getVel()*0.25 ); //?
1045                                 p->advanceVel();
1046                                 //errMsg("PIT","IN pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
1047                         } else {
1048 #if LBM_INCLUDE_TESTSOLVERS==1
1049                                 if(useff) { mpTest->handleParticle(p, i,j,k); }
1050                                 else DEL_PART; 
1051 #else // LBM_INCLUDE_TESTSOLVERS==1
1052                                 DEL_PART; 
1053 #endif // LBM_INCLUDE_TESTSOLVERS==1
1054                                 //errMsg("PIT","OUT pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags())  );
1055                         }
1056                                 
1057                         // additional bnd jitter
1058                         if((0) && (useff) && (p->getLifeTime()<3.*mLevel[mMaxRefine].timestep)) {
1059                                 // use half butoff border 1/8
1060                                 int maxdw = (int)(mLevel[mMaxRefine].lSizex*0.125*0.5);
1061                                 if(maxdw<3) maxdw=3;
1062                                 if((j>=0)&&(j<=this->mSizey-1)) {
1063                                         if(ABS(i-(               cutval))<maxdw) { p->advance(  FLOAT_JITTBNDRAND( ABS(i-(               cutval))), 0.,0.); }
1064                                         if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance(  FLOAT_JITTBNDRAND( ABS(i-(this->mSizex-1-cutval))), 0.,0.); }
1065                                 }
1066                                 //if( (i<cutval)||(i>this->mSizex-1-cutval)|| //(j<cutval)||(j>this->mSizey-1-cutval)
1067                         }
1068                 }  // PART_FLOAT
1069                 
1070                 // unknown particle type        
1071                 else {
1072                         errMsg("LbmFsgrSolver::advanceParticles","PIT pit invalid type!? "<<p->getStatus() );
1073                 }
1074   }
1075         myTime_t parttend = getTime(); 
1076         debMsgStd("LbmFsgrSolver::advanceParticles",DM_MSG,"Time for particle update:"<< getTimeString(parttend-parttstart)<<", #particles:"<<mpParticles->getNumParticles() , 10 );
1077 }
1078
1079 void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) {
1080         int workSet = mLevel[mMaxRefine].setCurr;
1081         std::ostringstream name;
1082
1083         // debug - raw dump of ffrac values
1084         if(getenv("ELBEEM_RAWDEBUGDUMP")) {
1085                 //name <<"fill_" << this->mStepCnt <<".dump";
1086                 name << outfilename<< frameNrStr <<".dump";
1087                 FILE *file = fopen(name.str().c_str(),"w");
1088                 if(file) {
1089
1090                         for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
1091                                 for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
1092                                         for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
1093                                                 float val = 0.;
1094                                                 if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) {
1095                                                         val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
1096                                                         if(val<0.) val=0.;
1097                                                         if(val>1.) val=1.;
1098                                                 }
1099                                                 if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.;
1100                                                 fprintf(file, "%f ",val); // text
1101                                                 //errMsg("W", PRINT_IJK<<" val:"<<val);
1102                                         }
1103                                         fprintf(file, "\n"); // text
1104                                 }
1105                                 fprintf(file, "\n"); // text
1106                         }
1107                         fclose(file);
1108
1109                 } // file
1110         } // */
1111         if(getenv("ELBEEM_BINDEBUGDUMP")) {
1112                 name << outfilename<< frameNrStr <<".bdump";
1113                 FILE *file = fopen(name.str().c_str(),"w");
1114                 if(file) {
1115                         for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
1116                                 for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
1117                                         for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
1118                                                 float val = 0.;
1119                                                 if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) {
1120                                                         val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
1121                                                         if(val<0.) val=0.;
1122                                                         if(val>1.) val=1.;
1123                                                 }
1124                                                 if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.;
1125                                                 fwrite( &val, sizeof(val), 1, file); // binary
1126                                         }
1127                                 }
1128                         }
1129                         fclose(file);
1130                 } // file
1131         }
1132
1133         dumptype = 0; frameNr = 0; // get rid of warning
1134 }
1135
1136 /*****************************************************************************/
1137 /*! internal quick print function (for debugging) */
1138 /*****************************************************************************/
1139 void 
1140 LbmFsgrSolver::printLbmCell(int level, int i, int j, int k, int set) {
1141         stdCellId *newcid = new stdCellId;
1142         newcid->level = level;
1143         newcid->x = i;
1144         newcid->y = j;
1145         newcid->z = k;
1146
1147         // this function is not called upon clicking, then its from setMouseClick
1148         debugPrintNodeInfo( newcid, set );
1149         delete newcid;
1150 }
1151 void 
1152 LbmFsgrSolver::debugMarkCellCall(int level, int vi,int vj,int vk) {
1153         stdCellId *newcid = new stdCellId;
1154         newcid->level = level;
1155         newcid->x = vi;
1156         newcid->y = vj;
1157         newcid->z = vk;
1158         this->addCellToMarkedList( newcid );
1159 }
1160
1161                 
1162 /*****************************************************************************/
1163 // implement CellIterator<UniformFsgrCellIdentifier> interface
1164 /*****************************************************************************/
1165
1166
1167
1168 // values from guiflkt.cpp
1169 extern double guiRoiSX, guiRoiSY, guiRoiSZ, guiRoiEX, guiRoiEY, guiRoiEZ;
1170 extern int guiRoiMaxLev, guiRoiMinLev;
1171 #define CID_SX (int)( (mLevel[cid->level].lSizex-1) * guiRoiSX )
1172 #define CID_SY (int)( (mLevel[cid->level].lSizey-1) * guiRoiSY )
1173 #define CID_SZ (int)( (mLevel[cid->level].lSizez-1) * guiRoiSZ )
1174
1175 #define CID_EX (int)( (mLevel[cid->level].lSizex-1) * guiRoiEX )
1176 #define CID_EY (int)( (mLevel[cid->level].lSizey-1) * guiRoiEY )
1177 #define CID_EZ (int)( (mLevel[cid->level].lSizez-1) * guiRoiEZ )
1178
1179 CellIdentifierInterface* 
1180 LbmFsgrSolver::getFirstCell( ) {
1181         int level = mMaxRefine;
1182
1183 #if LBMDIM==3
1184         if(mMaxRefine>0) { level = mMaxRefine-1; } // NO1HIGHESTLEV DEBUG
1185 #endif
1186         level = guiRoiMaxLev;
1187         if(level>mMaxRefine) level = mMaxRefine;
1188         
1189         //errMsg("LbmFsgrSolver::getFirstCell","Celliteration started...");
1190         stdCellId *cid = new stdCellId;
1191         cid->level = level;
1192         cid->x = CID_SX;
1193         cid->y = CID_SY;
1194         cid->z = CID_SZ;
1195         return cid;
1196 }
1197
1198 LbmFsgrSolver::stdCellId* 
1199 LbmFsgrSolver::convertBaseCidToStdCid( CellIdentifierInterface* basecid) {
1200         //stdCellId *cid = dynamic_cast<stdCellId*>( basecid );
1201         stdCellId *cid = (stdCellId*)( basecid );
1202         return cid;
1203 }
1204
1205 void LbmFsgrSolver::advanceCell( CellIdentifierInterface* basecid) {
1206         stdCellId *cid = convertBaseCidToStdCid(basecid);
1207         if(cid->getEnd()) return;
1208
1209         //debugOut(" ADb "<<cid->x<<","<<cid->y<<","<<cid->z<<" e"<<cid->getEnd(), 10);
1210         cid->x++;
1211         if(cid->x > CID_EX){ cid->x = CID_SX; cid->y++; 
1212                 if(cid->y > CID_EY){ cid->y = CID_SY; cid->z++; 
1213                         if(cid->z > CID_EZ){ 
1214                                 cid->level--;
1215                                 cid->x = CID_SX; 
1216                                 cid->y = CID_SY; 
1217                                 cid->z = CID_SZ; 
1218                                 if(cid->level < guiRoiMinLev) {
1219                                         cid->level = guiRoiMaxLev;
1220                                         cid->setEnd( true );
1221                                 }
1222                         }
1223                 }
1224         }
1225         //debugOut(" ADa "<<cid->x<<","<<cid->y<<","<<cid->z<<" e"<<cid->getEnd(), 10);
1226 }
1227
1228 bool LbmFsgrSolver::noEndCell( CellIdentifierInterface* basecid) {
1229         stdCellId *cid = convertBaseCidToStdCid(basecid);
1230         return (!cid->getEnd());
1231 }
1232
1233 void LbmFsgrSolver::deleteCellIterator( CellIdentifierInterface** cid ) {
1234         delete *cid;
1235         *cid = NULL;
1236 }
1237
1238 CellIdentifierInterface* LbmFsgrSolver::getCellAt( ntlVec3Gfx pos ) {
1239         //int cellok = false;
1240         pos -= (this->mvGeoStart);
1241
1242         LbmFloat mmaxsize = mLevel[mMaxRefine].nodeSize;
1243         for(int level=mMaxRefine; level>=0; level--) { // finest first
1244         //for(int level=0; level<=mMaxRefine; level++) { // coarsest first
1245                 LbmFloat nsize = mLevel[level].nodeSize;
1246                 int x,y,z;
1247                 // CHECK +- maxsize?
1248                 x = (int)((pos[0]+0.5*mmaxsize) / nsize );
1249                 y = (int)((pos[1]+0.5*mmaxsize) / nsize );
1250                 z = (int)((pos[2]+0.5*mmaxsize) / nsize );
1251                 if(LBMDIM==2) z = 0;
1252
1253                 // double check...
1254                 if(x<0) continue;
1255                 if(y<0) continue;
1256                 if(z<0) continue;
1257                 if(x>=mLevel[level].lSizex) continue;
1258                 if(y>=mLevel[level].lSizey) continue;
1259                 if(z>=mLevel[level].lSizez) continue;
1260
1261                 // return fluid/if/border cells
1262                 if( ( (RFLAG(level, x,y,z, mLevel[level].setCurr)&(CFUnused)) ) ||
1263                           ( (level<mMaxRefine) && (RFLAG(level, x,y,z, mLevel[level].setCurr)&(CFUnused|CFEmpty)) ) ) {
1264                         continue;
1265                 } // */
1266
1267                 stdCellId *newcid = new stdCellId;
1268                 newcid->level = level;
1269                 newcid->x = x;
1270                 newcid->y = y;
1271                 newcid->z = z;
1272                 //errMsg("cellAt",this->mName<<" "<<pos<<" l"<<level<<":"<<x<<","<<y<<","<<z<<" "<<convertCellFlagType2String(RFLAG(level, x,y,z, mLevel[level].setCurr)) );
1273                 return newcid;
1274         }
1275
1276         return NULL;
1277 }
1278
1279
1280 // INFO functions
1281
1282 int      LbmFsgrSolver::getCellSet      ( CellIdentifierInterface* basecid) {
1283         stdCellId *cid = convertBaseCidToStdCid(basecid);
1284         return mLevel[cid->level].setCurr;
1285         //return mLevel[cid->level].setOther;
1286 }
1287
1288 int      LbmFsgrSolver::getCellLevel    ( CellIdentifierInterface* basecid) {
1289         stdCellId *cid = convertBaseCidToStdCid(basecid);
1290         return cid->level;
1291 }
1292
1293 ntlVec3Gfx   LbmFsgrSolver::getCellOrigin   ( CellIdentifierInterface* basecid) {
1294         ntlVec3Gfx ret;
1295
1296         stdCellId *cid = convertBaseCidToStdCid(basecid);
1297         ntlVec3Gfx cs( mLevel[cid->level].nodeSize );
1298         if(LBMDIM==2) { cs[2] = 0.0; }
1299
1300         if(LBMDIM==2) {
1301                 ret =(this->mvGeoStart + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], (this->mvGeoEnd[2]-this->mvGeoStart[2])*0.5 )
1302                                 + ntlVec3Gfx(0.0,0.0,cs[1]*-0.25)*cid->level )
1303                         +getCellSize(basecid);
1304         } else {
1305                 ret =(this->mvGeoStart + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], cid->z *cs[2] ))
1306                         +getCellSize(basecid);
1307         }
1308         return (ret);
1309 }
1310
1311 ntlVec3Gfx   LbmFsgrSolver::getCellSize     ( CellIdentifierInterface* basecid) {
1312         // return half size
1313         stdCellId *cid = convertBaseCidToStdCid(basecid);
1314         ntlVec3Gfx retvec( mLevel[cid->level].nodeSize * 0.5 );
1315         // 2d display as rectangles
1316         if(LBMDIM==2) { retvec[2] = 0.0; }
1317         return (retvec);
1318 }
1319
1320 LbmFloat LbmFsgrSolver::getCellDensity  ( CellIdentifierInterface* basecid,int set) {
1321         stdCellId *cid = convertBaseCidToStdCid(basecid);
1322
1323         LbmFloat rho = 0.0;
1324         FORDF0 { rho += QCELL(cid->level, cid->x,cid->y,cid->z, set, l); } // ORG
1325         return ((rho-1.0) * mLevel[cid->level].simCellSize / mLevel[cid->level].timestep) +1.0; // ORG
1326         /*if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFInter) { // test
1327                 LbmFloat ux,uy,uz;
1328                 ux=uy=uz= 0.0;
1329                 int lev = cid->level;
1330                 LbmFloat df[27], feqOld[27];
1331                 FORDF0 {
1332                         rho += QCELL(lev, cid->x,cid->y,cid->z, set, l);
1333                         ux += this->dfDvecX[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l);
1334                         uy += this->dfDvecY[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l);
1335                         uz += this->dfDvecZ[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l);
1336                         df[l] = QCELL(lev, cid->x,cid->y,cid->z, set, l);
1337                 }
1338                 FORDF0 {
1339                         feqOld[l] = getCollideEq(l, rho,ux,uy,uz); 
1340                 }
1341                 // debugging mods
1342                 //const LbmFloat Qo = this->getLesNoneqTensorCoeff(df,feqOld);
1343                 //const LbmFloat modOmega = this->getLesOmega(mLevel[lev].omega, mLevel[lev].lcsmago,Qo);
1344                 //rho = (2.0-modOmega) *25.0;
1345                 //rho = Qo*100.0;
1346                 //if(cid->x==24){ errMsg("MODOMT"," at "<<PRINT_VEC(cid->x,cid->y,cid->z)<<" = "<<rho<<" "<<Qo); }
1347                 //else{ rho=0.0; }
1348         } // test 
1349         return rho; // test */
1350 }
1351
1352 LbmVec   LbmFsgrSolver::getCellVelocity ( CellIdentifierInterface* basecid,int set) {
1353         stdCellId *cid = convertBaseCidToStdCid(basecid);
1354
1355         // skip non-fluid cells
1356         if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&(CFFluid|CFInter)) {
1357                 // ok go on...
1358         } else {
1359                 return LbmVec(0.0);
1360         }
1361
1362         LbmFloat ux,uy,uz;
1363         ux=uy=uz= 0.0;
1364         FORDF0 {
1365                 ux += this->dfDvecX[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
1366                 uy += this->dfDvecY[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
1367                 uz += this->dfDvecZ[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l);
1368         }
1369         LbmVec vel(ux,uy,uz);
1370         // TODO fix...
1371         return (vel * mLevel[cid->level].simCellSize / mLevel[cid->level].timestep * this->mDebugVelScale); // normal
1372 }
1373
1374 LbmFloat   LbmFsgrSolver::getCellDf( CellIdentifierInterface* basecid,int set, int dir) {
1375         stdCellId *cid = convertBaseCidToStdCid(basecid);
1376         return QCELL(cid->level, cid->x,cid->y,cid->z, set, dir);
1377 }
1378 LbmFloat   LbmFsgrSolver::getCellMass( CellIdentifierInterface* basecid,int set) {
1379         stdCellId *cid = convertBaseCidToStdCid(basecid);
1380         return QCELL(cid->level, cid->x,cid->y,cid->z, set, dMass);
1381 }
1382 LbmFloat   LbmFsgrSolver::getCellFill( CellIdentifierInterface* basecid,int set) {
1383         stdCellId *cid = convertBaseCidToStdCid(basecid);
1384         if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFInter) return QCELL(cid->level, cid->x,cid->y,cid->z, set, dFfrac);
1385         if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFFluid) return 1.0;
1386         return 0.0;
1387         //return QCELL(cid->level, cid->x,cid->y,cid->z, set, dFfrac);
1388 }
1389 CellFlagType LbmFsgrSolver::getCellFlag( CellIdentifierInterface* basecid,int set) {
1390         stdCellId *cid = convertBaseCidToStdCid(basecid);
1391         return RFLAG(cid->level, cid->x,cid->y,cid->z, set);
1392 }
1393
1394 LbmFloat LbmFsgrSolver::getEquilDf( int l ) {
1395         return this->dfEquil[l];
1396 }
1397
1398
1399 ntlVec3Gfx LbmFsgrSolver::getVelocityAt   (float xp, float yp, float zp) {
1400         ntlVec3Gfx avgvel(0.0);
1401         LbmFloat   avgnum = 0.;
1402
1403         // taken from getCellAt!
1404         const int level = mMaxRefine;
1405         const int workSet = mLevel[level].setCurr;
1406         const LbmFloat nsize = mLevel[level].nodeSize;
1407         const int x = (int)((-this->mvGeoStart[0]+xp-0.5*nsize) / nsize );
1408         const int y = (int)((-this->mvGeoStart[1]+yp-0.5*nsize) / nsize );
1409         int       z = (int)((-this->mvGeoStart[2]+zp-0.5*nsize) / nsize );
1410         if(LBMDIM==2) z=0;
1411         //errMsg("DUMPVEL","p"<<PRINT_VEC(xp,yp,zp)<<" at "<<PRINT_VEC(x,y,z)<<" max"<<PRINT_VEC(mLevel[level].lSizex,mLevel[level].lSizey,mLevel[level].lSizez) );
1412
1413         // return fluid/if/border cells
1414         // search neighborhood, do smoothing
1415         FORDF0{ 
1416                 const int i = x+this->dfVecX[l];
1417                 const int j = y+this->dfVecY[l];
1418                 const int k = z+this->dfVecZ[l];
1419
1420                 if( (i<0) || (j<0) || (k<0) 
1421                  || (i>=mLevel[level].lSizex) 
1422                  || (j>=mLevel[level].lSizey) 
1423                  || (k>=mLevel[level].lSizez) ) continue;
1424
1425                 if( (RFLAG(level, i,j,k, mLevel[level].setCurr)&(CFFluid|CFInter)) ) {
1426                         ntlVec3Gfx vel(0.0);
1427                         LbmFloat *ccel = RACPNT(level, i,j,k ,workSet); // omp
1428                         for(int n=1; n<this->cDfNum; n++) {
1429                                 vel[0]  += (this->dfDvecX[n]*RAC(ccel,n));
1430                                 vel[1]  += (this->dfDvecY[n]*RAC(ccel,n)); 
1431                                 vel[2]  += (this->dfDvecZ[n]*RAC(ccel,n)); 
1432                         } 
1433
1434                         avgvel += vel;
1435                         avgnum += 1.0;
1436                         if(l==0) { // center slightly more weight
1437                                 avgvel += vel; avgnum += 1.0;
1438                         }
1439                 } // */
1440         }
1441
1442         if(avgnum>0.) {
1443                 ntlVec3Gfx retv = avgvel / avgnum;
1444                 retv *= nsize/mLevel[level].timestep;
1445                 // scale for current animation settings (frame time)
1446                 retv *= mpParam->getCurrentAniFrameTime();
1447                 //errMsg("DUMPVEL","t"<<mSimulationTime<<" at "<<PRINT_VEC(xp,yp,zp)<<" ret:"<<retv<<", avgv:"<<avgvel<<" n"<<avgnum<<" nsize"<<nsize<<" ts"<<mLevel[level].timestep<<" fr"<<mpParam->getCurrentAniFrameTime() );
1448                 return retv;
1449         }
1450         // no cells here...?
1451         //errMsg("DUMPVEL"," at "<<PRINT_VEC(xp,yp,zp)<<" v"<<avgvel<<" n"<<avgnum<<" no vel !?");
1452         return ntlVec3Gfx(0.);
1453 }
1454
1455 #if LBM_USE_GUI==1
1456 //! show simulation info (implement SimulationObject pure virtual func)
1457 void 
1458 LbmFsgrSolver::debugDisplay(int set){ 
1459         //lbmDebugDisplay< LbmFsgrSolver >( set, this ); 
1460         lbmDebugDisplay( set ); 
1461 }
1462 #endif
1463
1464 /*****************************************************************************/
1465 // strict debugging functions
1466 /*****************************************************************************/
1467 #if FSGR_STRICT_DEBUG==1
1468 #define STRICT_EXIT *((int *)0)=0;
1469
1470 int LbmFsgrSolver::debLBMGI(int level, int ii,int ij,int ik, int is) {
1471         if(level <  0){ errMsg("LbmStrict::debLBMGI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
1472         if(level >  mMaxRefine){ errMsg("LbmStrict::debLBMGI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
1473
1474         if((ii==-1)&&(ij==0)) {
1475                 // special case for main loop, ok
1476         } else {
1477                 if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1478                 if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1479                 if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1480                 if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1481         }
1482         if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1483         if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1484         if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1485         if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1486         return _LBMGI(level, ii,ij,ik, is);
1487 };
1488
1489 CellFlagType& LbmFsgrSolver::debRFLAG(int level, int xx,int yy,int zz,int set){
1490         return _RFLAG(level, xx,yy,zz,set);   
1491 };
1492
1493 CellFlagType& LbmFsgrSolver::debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir) {
1494         if(dir<0)         { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1495         // warning might access all spatial nbs
1496         if(dir>this->cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1497         return _RFLAG_NB(level, xx,yy,zz,set, dir);
1498 };
1499
1500 CellFlagType& LbmFsgrSolver::debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir) {
1501         if(dir<0)         { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1502         if(dir>this->cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1503         return _RFLAG_NBINV(level, xx,yy,zz,set, dir);
1504 };
1505
1506 int LbmFsgrSolver::debLBMQI(int level, int ii,int ij,int ik, int is, int l) {
1507         if(level <  0){ errMsg("LbmStrict::debLBMQI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
1508         if(level >  mMaxRefine){ errMsg("LbmStrict::debLBMQI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 
1509
1510         if((ii==-1)&&(ij==0)) {
1511                 // special case for main loop, ok
1512         } else {
1513                 if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1514                 if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1515                 if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1516                 if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1517         }
1518         if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1519         if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1520         if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1521         if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
1522         if(l<0)        { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; }
1523         if(l>this->cDfNum){  // dFfrac is an exception
1524                 if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } }
1525 #if COMPRESSGRIDS==1
1526         //if((!this->mInitDone) && (is!=mLevel[level].setCurr)){ STRICT_EXIT; } // COMPRT debug
1527 #endif // COMPRESSGRIDS==1
1528         return _LBMQI(level, ii,ij,ik, is, l);
1529 };
1530
1531 LbmFloat& LbmFsgrSolver::debQCELL(int level, int xx,int yy,int zz,int set,int l) {
1532         //errMsg("LbmStrict","debQCELL debug: l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" l"<<l<<" index"<<LBMGI(level, xx,yy,zz,set)); 
1533         return _QCELL(level, xx,yy,zz,set,l);
1534 };
1535
1536 LbmFloat& LbmFsgrSolver::debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l) {
1537         if(dir<0)        { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1538         if(dir>this->cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1539         return _QCELL_NB(level, xx,yy,zz,set, dir,l);
1540 };
1541
1542 LbmFloat& LbmFsgrSolver::debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l) {
1543         if(dir<0)        { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1544         if(dir>this->cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; }
1545         return _QCELL_NBINV(level, xx,yy,zz,set, dir,l);
1546 };
1547
1548 LbmFloat* LbmFsgrSolver::debRACPNT(int level,  int ii,int ij,int ik, int is ) {
1549         return _RACPNT(level, ii,ij,ik, is );
1550 };
1551
1552 LbmFloat& LbmFsgrSolver::debRAC(LbmFloat* s,int l) {
1553         if(l<0)        { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; }
1554         if(l>dTotalNum){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } 
1555         //if(l>this->cDfNum){ // dFfrac is an exception 
1556         //if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } }
1557         return _RAC(s,l);
1558 };
1559
1560 #endif // FSGR_STRICT_DEBUG==1
1561
1562
1563 /******************************************************************************
1564  * GUI&debugging functions
1565  *****************************************************************************/
1566
1567
1568 #if LBM_USE_GUI==1
1569 #define USE_GLUTILITIES
1570 #include "../gui/gui_utilities.h"
1571
1572 //! display a single node
1573 void LbmFsgrSolver::debugDisplayNode(int dispset, CellIdentifierInterface* cell ) {
1574         //debugOut(" DD: "<<cell->getAsString() , 10);
1575         ntlVec3Gfx org      = this->getCellOrigin( cell );
1576         ntlVec3Gfx halfsize = this->getCellSize( cell );
1577         int    set      = this->getCellSet( cell );
1578         //debugOut(" DD: "<<cell->getAsString()<<" "<< (dispset->type) , 10);
1579
1580         bool     showcell = true;
1581         int      linewidth = 1;
1582         ntlColor col(0.5);
1583         LbmFloat cscale = 1.0; //dispset->scale;
1584
1585 #define DRAWDISPCUBE(col,scale) \
1586         {       glLineWidth( linewidth ); \
1587           glColor3f( (col)[0], (col)[1], (col)[2]); \
1588                 ntlVec3Gfx s = org-(halfsize * (scale)); \
1589                 ntlVec3Gfx e = org+(halfsize * (scale)); \
1590                 drawCubeWire( s,e ); }
1591
1592         CellFlagType flag = this->getCellFlag(cell, set );
1593         // always check types
1594         if(flag& CFInvalid  ) { if(!guiShowInvalid  ) return; }
1595         if(flag& CFUnused   ) { if(!guiShowInvalid  ) return; }
1596         if(flag& CFEmpty    ) { if(!guiShowEmpty    ) return; }
1597         if(flag& CFInter    ) { if(!guiShowInterface) return; }
1598         if(flag& CFNoDelete ) { if(!guiShowNoDelete ) return; }
1599         if(flag& CFBnd      ) { if(!guiShowBnd      ) return; }
1600
1601         // only dismiss one of these types 
1602         if(flag& CFGrFromCoarse)  { if(!guiShowCoarseInner  ) return; } // inner not really interesting
1603         else
1604         if(flag& CFGrFromFine) { if(!guiShowCoarseBorder ) return; }
1605         else
1606         if(flag& CFFluid    )    { if(!guiShowFluid    ) return; }
1607
1608         switch(dispset) {
1609                 case FLUIDDISPNothing: {
1610                                 showcell = false;
1611                         } break;
1612                 case FLUIDDISPCelltypes: {
1613                                 cscale = 0.5;
1614
1615                                 if(flag& CFNoDelete) { // debug, mark nodel cells
1616                                         ntlColor ccol(0.7,0.0,0.0);
1617                                         DRAWDISPCUBE(ccol, 0.1);
1618                                 }
1619                                 if(flag& CFPersistMask) { // mark persistent flags
1620                                         ntlColor ccol(0.5);
1621                                         DRAWDISPCUBE(ccol, 0.125);
1622                                 }
1623                                 if(flag& CFNoBndFluid) { // mark persistent flags
1624                                         ntlColor ccol(0,0,1);
1625                                         DRAWDISPCUBE(ccol, 0.075);
1626                                 }
1627
1628                                 if(flag& CFInvalid) {
1629                                         cscale = 0.50;
1630                                         col = ntlColor(0.0,0,0.0);
1631                                 }
1632                                 else if(flag& CFBnd) {
1633                                         cscale = 0.59;
1634                                         col = ntlColor(0.4);
1635                                 }
1636
1637                                 else if(flag& CFInter) {
1638                                         cscale = 0.55;
1639                                         col = ntlColor(0,1,1);
1640
1641                                 } else if(flag& CFGrFromCoarse) {
1642                                         // draw as - with marker
1643                                         ntlColor col2(0.0,1.0,0.3);
1644                                         DRAWDISPCUBE(col2, 0.1);
1645                                         cscale = 0.5;
1646                                         showcell=false; // DEBUG
1647                                 }
1648                                 else if(flag& CFFluid) {
1649                                         cscale = 0.5;
1650                                         if(flag& CFGrToFine) {
1651                                                 ntlColor col2(0.5,0.0,0.5);
1652                                                 DRAWDISPCUBE(col2, 0.1);
1653                                                 col = ntlColor(0,0,1);
1654                                         }
1655                                         if(flag& CFGrFromFine) {
1656                                                 ntlColor col2(1.0,1.0,0.0);
1657                                                 DRAWDISPCUBE(col2, 0.1);
1658                                                 col = ntlColor(0,0,1);
1659                                         } else if(flag& CFGrFromCoarse) {
1660                                                 // draw as fluid with marker
1661                                                 ntlColor col2(0.0,1.0,0.3);
1662                                                 DRAWDISPCUBE(col2, 0.1);
1663                                                 col = ntlColor(0,0,1);
1664                                         } else {
1665                                                 col = ntlColor(0,0,1);
1666                                         }
1667                                 }
1668                                 else if(flag& CFEmpty) {
1669                                         showcell=false;
1670                                 }
1671
1672                         } break;
1673                 case FLUIDDISPVelocities: {
1674                                 // dont use cube display
1675                                 LbmVec vel = this->getCellVelocity( cell, set );
1676                                 glBegin(GL_LINES);
1677                                 glColor3f( 0.0,0.0,0.0 );
1678                                 glVertex3f( org[0], org[1], org[2] );
1679                                 org += vec2G(vel * 10.0 * cscale);
1680                                 glColor3f( 1.0,1.0,1.0 );
1681                                 glVertex3f( org[0], org[1], org[2] );
1682                                 glEnd();
1683                                 showcell = false;
1684                         } break;
1685                 case FLUIDDISPCellfills: {
1686                                 cscale = 0.5;
1687                                 if(flag& CFFluid) {
1688                                         cscale = 0.75;
1689                                         col = ntlColor(0,0,0.5);
1690                                 }
1691                                 else if(flag& CFInter) {
1692                                         cscale = 0.75 * this->getCellMass(cell,set);
1693                                         col = ntlColor(0,1,1);
1694                                 }
1695                                 else {
1696                                         showcell=false;
1697                                 }
1698
1699                                         if( ABS(this->getCellMass(cell,set)) < 10.0 ) {
1700                                                 cscale = 0.75 * this->getCellMass(cell,set);
1701                                         } else {
1702                                                 showcell = false;
1703                                         }
1704                                         if(cscale>0.0) {
1705                                                 col = ntlColor(0,1,1);
1706                                         } else {
1707                                                 col = ntlColor(1,1,0);
1708                                         }
1709                         // TODO
1710                         } break;
1711                 case FLUIDDISPDensity: {
1712                                 LbmFloat rho = this->getCellDensity(cell,set);
1713                                 cscale = rho*rho * 0.25;
1714                                 col = ntlColor( MIN(0.5+cscale,1.0) , MIN(0.0+cscale,1.0), MIN(0.0+cscale,1.0) );
1715                                 cscale *= 2.0;
1716                         } break;
1717                 case FLUIDDISPGrid: {
1718                                 cscale = 0.59;
1719                                 col = ntlColor(1.0);
1720                         } break;
1721                 default: {
1722                                 cscale = 0.5;
1723                                 col = ntlColor(1.0,0.0,0.0);
1724                         } break;
1725         }
1726
1727         if(!showcell) return;
1728         if(cscale==0.0) return; // dont draw zero values
1729         DRAWDISPCUBE(col, cscale);
1730 }
1731
1732 //! debug display function
1733 //  D has to implement the CellIterator interface
1734 void LbmFsgrSolver::lbmDebugDisplay(int dispset) {
1735         // DEBUG always display testdata
1736 #if LBM_INCLUDE_TESTSOLVERS==1
1737         if(mUseTestdata){ 
1738                 cpDebugDisplay(dispset); 
1739                 mpTest->testDebugDisplay(dispset); 
1740         }
1741 #endif // LBM_INCLUDE_TESTSOLVERS==1
1742         if(dispset<=FLUIDDISPNothing) return;
1743         //if(!dispset->on) return;
1744         glDisable( GL_LIGHTING ); // dont light lines
1745
1746 #if LBM_INCLUDE_TESTSOLVERS==1
1747         if((!mUseTestdata)|| (mUseTestdata)&&(mpTest->mDebugvalue1<=0.0)) {
1748 #endif // LBM_INCLUDE_TESTSOLVERS==1
1749
1750         LbmFsgrSolver::CellIdentifier cid = this->getFirstCell();
1751         for(; this->noEndCell( cid );
1752               this->advanceCell( cid ) ) {
1753                 this->debugDisplayNode(dispset, cid );
1754         }
1755         delete cid;
1756
1757 #if LBM_INCLUDE_TESTSOLVERS==1
1758         } // 3d check
1759 #endif // LBM_INCLUDE_TESTSOLVERS==1
1760
1761         glEnable( GL_LIGHTING ); // dont light lines
1762 }
1763
1764 //! debug display function
1765 //  D has to implement the CellIterator interface
1766 void LbmFsgrSolver::lbmMarkedCellDisplay() {
1767         //fluidDispSettings dispset;
1768         // trick - display marked cells as grid displa -> white, big
1769         int dispset = FLUIDDISPGrid;
1770         glDisable( GL_LIGHTING ); // dont light lines
1771         
1772         LbmFsgrSolver::CellIdentifier cid = this->markedGetFirstCell();
1773         while(cid) {
1774                 this->debugDisplayNode(dispset, cid );
1775                 cid = this->markedAdvanceCell();
1776         }
1777         delete cid;
1778
1779         glEnable( GL_LIGHTING ); // dont light lines
1780 }
1781
1782 #endif // LBM_USE_GUI==1
1783
1784 //! display a single node
1785 void LbmFsgrSolver::debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet) {
1786                 //string printInfo,
1787                 // force printing of one set? default = -1 = off
1788   bool printDF     = false;
1789   bool printRho    = false;
1790   bool printVel    = false;
1791   bool printFlag   = false;
1792   bool printGeom   = false;
1793   bool printMass=false;
1794         bool printBothSets = false;
1795         string printInfo = this->getNodeInfoString();
1796
1797         for(size_t i=0; i<printInfo.length()-0; i++) {
1798                 char what = printInfo[i];
1799                 switch(what) {
1800                         case '+': // all on
1801                                                                 printDF = true; printRho = true; printVel = true; printFlag = true; printGeom = true; printMass = true ;
1802                                                                 printBothSets = true; break;
1803                         case '-': // all off
1804                                                                 printDF = false; printRho = false; printVel = false; printFlag = false; printGeom = false; printMass = false; 
1805                                                                 printBothSets = false; break;
1806                         case 'd': printDF = true; break;
1807                         case 'r': printRho = true; break;
1808                         case 'v': printVel = true; break;
1809                         case 'f': printFlag = true; break;
1810                         case 'g': printGeom = true; break;
1811                         case 'm': printMass = true; break;
1812                         case 's': printBothSets = true; break;
1813                         default: 
1814                                 errFatal("debugPrintNodeInfo","Invalid node info id "<<what,SIMWORLD_GENERICERROR); return;
1815                 }
1816         }
1817
1818         ntlVec3Gfx org      = this->getCellOrigin( cell );
1819         ntlVec3Gfx halfsize = this->getCellSize( cell );
1820         int    set      = this->getCellSet( cell );
1821         debMsgStd("debugPrintNodeInfo",DM_NOTIFY, "Printing cell info '"<<printInfo<<"' for node: "<<cell->getAsString()<<" from "<<this->getName()<<" currSet:"<<set , 1);
1822         if(printGeom) debMsgStd("                  ",DM_MSG, "Org:"<<org<<" Halfsize:"<<halfsize<<" ", 1);
1823
1824         int setmax = 2;
1825         if(!printBothSets) setmax = 1;
1826         if(forceSet>=0) setmax = 1;
1827
1828         for(int s=0; s<setmax; s++) {
1829                 int workset = set;
1830                 if(s==1){ workset = (set^1); }          
1831                 if(forceSet>=0) workset = forceSet;
1832                 debMsgStd("                  ",DM_MSG, "Printing set:"<<workset<<" orgSet:"<<set, 1);
1833                 
1834                 if(printDF) {
1835                         for(int l=0; l<LBM_DFNUM; l++) { // FIXME ??
1836                                 debMsgStd("                  ",DM_MSG, "  Df"<<l<<": "<<this->getCellDf(cell,workset,l), 1);
1837                         }
1838                 }
1839                 if(printRho) {
1840                         debMsgStd("                  ",DM_MSG, "  Rho: "<<this->getCellDensity(cell,workset), 1);
1841                 }
1842                 if(printVel) {
1843                         debMsgStd("                  ",DM_MSG, "  Vel: "<<this->getCellVelocity(cell,workset), 1);
1844                 }
1845                 if(printFlag) {
1846                         CellFlagType flag = this->getCellFlag(cell,workset);
1847                         debMsgStd("                  ",DM_MSG, "  Flg: "<< flag<<" "<<convertFlags2String( flag ) <<" "<<convertCellFlagType2String( flag ), 1);
1848                 }
1849                 if(printMass) {
1850                         debMsgStd("                  ",DM_MSG, "  Mss: "<<this->getCellMass(cell,workset), 1);
1851                 }
1852                 // dirty... TODO fixme
1853                 debMsgStd("                  ",DM_MSG, "  Flx: "<<this->getCellDf(cell,workset,dFlux), 1);
1854         }
1855 }
1856
1857