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