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