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