 |
Blender V2.61 - r43446
|
solver_util.cpp
Go to the documentation of this file.
00001 00004 /****************************************************************************** 00005 * 00006 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method 00007 * Copyright 2003-2006 Nils Thuerey 00008 * 00009 * Standard LBM Factory implementation 00010 * 00011 *****************************************************************************/ 00012 00013 #include "solver_class.h" 00014 #include "solver_relax.h" 00015 #include "particletracer.h" 00016 00017 // MPT 00018 #include "ntl_world.h" 00019 #include "simulation_object.h" 00020 00021 #include <stdlib.h> 00022 #include <zlib.h> 00023 #ifndef sqrtf 00024 #define sqrtf sqrt 00025 #endif 00026 00027 /****************************************************************************** 00028 * helper functions 00029 *****************************************************************************/ 00030 00031 // try to enhance surface? 00032 #define SURFACE_ENH 2 00033 00034 extern bool glob_mpactive; 00035 extern bool glob_mpnum; 00036 extern bool glob_mpindex; 00037 00039 void LbmFsgrSolver::prepareVisualization( void ) { 00040 int lev = mMaxRefine; 00041 int workSet = mLevel[lev].setCurr; 00042 00043 int mainGravDir=6; // if normalizing fails, we asume z-direction gravity 00044 LbmFloat mainGravLen = 0.; 00045 FORDF1{ 00046 LbmFloat thisGravLen = dot(LbmVec(dfVecX[l],dfVecY[l],dfVecZ[l]), mLevel[mMaxRefine].gravity ); 00047 if(thisGravLen>mainGravLen) { 00048 mainGravLen = thisGravLen; 00049 mainGravDir = l; 00050 } 00051 } 00052 00053 #if LBMDIM==2 00054 // 2d, place in the middle of isofield slice (k=2) 00055 # define ZKD1 0 00056 // 2d z offset = 2, lbmGetData adds 1, so use one here 00057 # define ZKOFF 1 00058 // reset all values... 00059 for(int k= 0; k< 5; ++k) 00060 for(int j=0;j<mLevel[lev].lSizey-0;j++) 00061 for(int i=0;i<mLevel[lev].lSizex-0;i++) { 00062 *mpIso->lbmGetData(i,j,ZKOFF)=0.0; 00063 } 00064 #else // LBMDIM==2 00065 // 3d, use normal bounds 00066 # define ZKD1 1 00067 # define ZKOFF k 00068 // reset all values... 00069 for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 00070 for(int j=0;j<mLevel[lev].lSizey-0;j++) 00071 for(int i=0;i<mLevel[lev].lSizex-0;i++) { 00072 *mpIso->lbmGetData(i,j,ZKOFF)=0.0; 00073 } 00074 #endif // LBMDIM==2 00075 00076 // MPT, ignore 00077 if((glob_mpactive) && (glob_mpnum>1) && (glob_mpindex==0)) { 00078 mpIso->resetAll(0.); 00079 } 00080 00081 00082 LbmFloat minval = mIsoValue*1.05; // / mIsoWeight[13]; 00083 // add up... 00084 float val = 0.0; 00085 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) 00086 for(int j=1;j<mLevel[lev].lSizey-1;j++) 00087 for(int i=1;i<mLevel[lev].lSizex-1;i++) { 00088 const CellFlagType cflag = RFLAG(lev, i,j,k,workSet); 00089 //if(cflag&(CFBnd|CFEmpty)) { 00090 00091 #if SURFACE_ENH==0 00092 00093 // no enhancements... 00094 if( (cflag&(CFFluid|CFUnused)) ) { 00095 val = 1.; 00096 } else if( (cflag&CFInter) ) { 00097 val = (QCELL(lev, i,j,k,workSet, dFfrac)); 00098 } else { 00099 continue; 00100 } 00101 00102 #else // SURFACE_ENH!=1 00103 if(cflag&CFBnd) { 00104 // treated in second loop 00105 continue; 00106 } else if(cflag&CFUnused) { 00107 val = 1.; 00108 } else if( (cflag&CFFluid) && (cflag&CFNoBndFluid)) { 00109 // optimized fluid 00110 val = 1.; 00111 } else if( (cflag&(CFEmpty|CFInter|CFFluid)) ) { 00112 int noslipbnd = 0; 00113 int intercnt = 0; 00114 FORDF1 { 00115 const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l); 00116 if(nbflag&CFInter){ intercnt++; } 00117 00118 // check all directions otherwise we get bugs with splashes on obstacles 00119 if(l!=mainGravDir) continue; // only check bnd along main grav. dir 00120 //if((nbflag&CFBnd)&&(nbflag&CFBndNoslip)){ noslipbnd=1; } 00121 if((nbflag&CFBnd)){ noslipbnd=1; } 00122 } 00123 00124 if(cflag&CFEmpty) { 00125 // special empty treatment 00126 if((noslipbnd)&&(intercnt>6)) { 00127 *mpIso->lbmGetData(i,j,ZKOFF) += minval; 00128 } else if((noslipbnd)&&(intercnt>0)) { 00129 // necessary? 00130 *mpIso->lbmGetData(i,j,ZKOFF) += mIsoValue*0.9; 00131 } else { 00132 // nothing to do... 00133 } 00134 00135 continue; 00136 } else if(cflag&(CFNoNbEmpty|CFFluid)) { 00137 // no empty nb interface cells are treated as full 00138 val=1.0; 00139 } else { 00140 val = (QCELL(lev, i,j,k,workSet, dFfrac)); 00141 } 00142 00143 if(noslipbnd) { 00144 if(val<minval) val = minval; 00145 *mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] ); 00146 } 00147 } else { // all others, unused? 00148 continue; 00149 } 00150 #endif // SURFACE_ENH>0 00151 00152 *mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] ); 00153 *mpIso->lbmGetData( i , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[1] ); 00154 *mpIso->lbmGetData( i+1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[2] ); 00155 00156 *mpIso->lbmGetData( i-1 , j ,ZKOFF-ZKD1) += ( val * mIsoWeight[3] ); 00157 *mpIso->lbmGetData( i , j ,ZKOFF-ZKD1) += ( val * mIsoWeight[4] ); 00158 *mpIso->lbmGetData( i+1 , j ,ZKOFF-ZKD1) += ( val * mIsoWeight[5] ); 00159 00160 *mpIso->lbmGetData( i-1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[6] ); 00161 *mpIso->lbmGetData( i , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[7] ); 00162 *mpIso->lbmGetData( i+1 , j+1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[8] ); 00163 00164 00165 *mpIso->lbmGetData( i-1 , j-1 ,ZKOFF ) += ( val * mIsoWeight[9] ); 00166 *mpIso->lbmGetData( i , j-1 ,ZKOFF ) += ( val * mIsoWeight[10] ); 00167 *mpIso->lbmGetData( i+1 , j-1 ,ZKOFF ) += ( val * mIsoWeight[11] ); 00168 00169 *mpIso->lbmGetData( i-1 , j ,ZKOFF ) += ( val * mIsoWeight[12] ); 00170 *mpIso->lbmGetData( i , j ,ZKOFF ) += ( val * mIsoWeight[13] ); 00171 *mpIso->lbmGetData( i+1 , j ,ZKOFF ) += ( val * mIsoWeight[14] ); 00172 00173 *mpIso->lbmGetData( i-1 , j+1 ,ZKOFF ) += ( val * mIsoWeight[15] ); 00174 *mpIso->lbmGetData( i , j+1 ,ZKOFF ) += ( val * mIsoWeight[16] ); 00175 *mpIso->lbmGetData( i+1 , j+1 ,ZKOFF ) += ( val * mIsoWeight[17] ); 00176 00177 00178 *mpIso->lbmGetData( i-1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[18] ); 00179 *mpIso->lbmGetData( i , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[19] ); 00180 *mpIso->lbmGetData( i+1 , j-1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[20] ); 00181 00182 *mpIso->lbmGetData( i-1 , j ,ZKOFF+ZKD1) += ( val * mIsoWeight[21] ); 00183 *mpIso->lbmGetData( i , j ,ZKOFF+ZKD1)+= ( val * mIsoWeight[22] ); 00184 *mpIso->lbmGetData( i+1 , j ,ZKOFF+ZKD1) += ( val * mIsoWeight[23] ); 00185 00186 *mpIso->lbmGetData( i-1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[24] ); 00187 *mpIso->lbmGetData( i , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[25] ); 00188 *mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] ); 00189 } 00190 00191 // TEST!? 00192 #if SURFACE_ENH>=2 00193 00194 if(mFsSurfGenSetting&fssgNoObs) { 00195 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) 00196 for(int j=1;j<mLevel[lev].lSizey-1;j++) 00197 for(int i=1;i<mLevel[lev].lSizex-1;i++) { 00198 const CellFlagType cflag = RFLAG(lev, i,j,k,workSet); 00199 if(cflag&(CFBnd)) { 00200 CellFlagType nbored=0; 00201 LbmFloat avgfill=0.,avgfcnt=0.; 00202 FORDF1 { 00203 const int ni = i+dfVecX[l]; 00204 const int nj = j+dfVecY[l]; 00205 const int nk = ZKOFF+dfVecZ[l]; 00206 00207 const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet); 00208 nbored |= nbflag; 00209 if(nbflag&CFInter) { 00210 avgfill += QCELL(lev, ni,nj,nk, workSet,dFfrac); avgfcnt += 1.; 00211 } else if(nbflag&CFFluid) { 00212 avgfill += 1.; avgfcnt += 1.; 00213 } else if(nbflag&CFEmpty) { 00214 avgfill += 0.; avgfcnt += 1.; 00215 } 00216 00217 //if( (ni<0) || (nj<0) || (nk<0) 00218 //|| (ni>=mLevel[mMaxRefine].lSizex) 00219 //|| (nj>=mLevel[mMaxRefine].lSizey) 00220 //|| (nk>=mLevel[mMaxRefine].lSizez) ) continue; 00221 } 00222 00223 if(nbored&CFInter) { 00224 if(avgfcnt>0.) avgfill/=avgfcnt; 00225 *mpIso->lbmGetData(i,j,ZKOFF) = avgfill; continue; 00226 } 00227 else if(nbored&CFFluid) { 00228 *mpIso->lbmGetData(i,j,ZKOFF) = 1.; continue; 00229 } 00230 00231 } 00232 } 00233 00234 // move surface towards inner "row" of obstacle 00235 // cells if necessary (all obs cells without fluid/inter 00236 // nbs (=iso==0) next to obstacles...) 00237 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) 00238 for(int j=1;j<mLevel[lev].lSizey-1;j++) 00239 for(int i=1;i<mLevel[lev].lSizex-1;i++) { 00240 const CellFlagType cflag = RFLAG(lev, i,j,k,workSet); 00241 if( (cflag&(CFBnd)) && (*mpIso->lbmGetData(i,j,ZKOFF)==0.)) { 00242 int bndnbcnt=0; 00243 FORDF1 { 00244 const int ni = i+dfVecX[l]; 00245 const int nj = j+dfVecY[l]; 00246 const int nk = ZKOFF+dfVecZ[l]; 00247 const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet); 00248 if(nbflag&CFBnd) bndnbcnt++; 00249 } 00250 if(bndnbcnt>0) *mpIso->lbmGetData(i,j,ZKOFF)=mIsoValue*0.95; 00251 } 00252 } 00253 } 00254 // */ 00255 00256 if(mFsSurfGenSetting&fssgNoNorth) 00257 for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 00258 for(int j=0;j<mLevel[lev].lSizey-0;j++) { 00259 *mpIso->lbmGetData(0, j,ZKOFF) = *mpIso->lbmGetData(1, j,ZKOFF); 00260 } 00261 if(mFsSurfGenSetting&fssgNoEast) 00262 for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 00263 for(int i=0;i<mLevel[lev].lSizex-0;i++) { 00264 *mpIso->lbmGetData(i,0, ZKOFF) = *mpIso->lbmGetData(i,1, ZKOFF); 00265 } 00266 if(mFsSurfGenSetting&fssgNoSouth) 00267 for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 00268 for(int j=0;j<mLevel[lev].lSizey-0;j++) { 00269 *mpIso->lbmGetData(mLevel[lev].lSizex-1,j,ZKOFF) = *mpIso->lbmGetData(mLevel[lev].lSizex-2,j,ZKOFF); 00270 } 00271 if(mFsSurfGenSetting&fssgNoWest) 00272 for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k) 00273 for(int i=0;i<mLevel[lev].lSizex-0;i++) { 00274 *mpIso->lbmGetData(i,mLevel[lev].lSizey-1,ZKOFF) = *mpIso->lbmGetData(i,mLevel[lev].lSizey-2,ZKOFF); 00275 } 00276 if(LBMDIM>2) { 00277 if(mFsSurfGenSetting&fssgNoBottom) 00278 for(int j=0;j<mLevel[lev].lSizey-0;j++) 00279 for(int i=0;i<mLevel[lev].lSizex-0;i++) { 00280 *mpIso->lbmGetData(i,j,0 ) = *mpIso->lbmGetData(i,j,1 ); 00281 } 00282 if(mFsSurfGenSetting&fssgNoTop) 00283 for(int j=0;j<mLevel[lev].lSizey-0;j++) 00284 for(int i=0;i<mLevel[lev].lSizex-0;i++) { 00285 *mpIso->lbmGetData(i,j,mLevel[lev].lSizez-1) = *mpIso->lbmGetData(i,j,mLevel[lev].lSizez-2); 00286 } 00287 } 00288 #endif // SURFACE_ENH>=2 00289 00290 00291 // update preview, remove 2d? 00292 if((mOutputSurfacePreview)&&(LBMDIM==3)) { 00293 int pvsx = (int)(mPreviewFactor*mSizex); 00294 int pvsy = (int)(mPreviewFactor*mSizey); 00295 int pvsz = (int)(mPreviewFactor*mSizez); 00296 //float scale = (float)mSizex / previewSize; 00297 LbmFloat scalex = (LbmFloat)mSizex/(LbmFloat)pvsx; 00298 LbmFloat scaley = (LbmFloat)mSizey/(LbmFloat)pvsy; 00299 LbmFloat scalez = (LbmFloat)mSizez/(LbmFloat)pvsz; 00300 for(int k= 0; k< (pvsz-1); ++k) 00301 for(int j=0;j< pvsy;j++) 00302 for(int i=0;i< pvsx;i++) { 00303 *mpPreviewSurface->lbmGetData(i,j,k) = *mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley), (int)(k*scalez) ); 00304 } 00305 // set borders again... 00306 for(int k= 0; k< (pvsz-1); ++k) { 00307 for(int j=0;j< pvsy;j++) { 00308 *mpPreviewSurface->lbmGetData(0,j,k) = *mpIso->lbmGetData( 0, (int)(j*scaley), (int)(k*scalez) ); 00309 *mpPreviewSurface->lbmGetData(pvsx-1,j,k) = *mpIso->lbmGetData( mSizex-1, (int)(j*scaley), (int)(k*scalez) ); 00310 } 00311 for(int i=0;i< pvsx;i++) { 00312 *mpPreviewSurface->lbmGetData(i,0,k) = *mpIso->lbmGetData( (int)(i*scalex), 0, (int)(k*scalez) ); 00313 *mpPreviewSurface->lbmGetData(i,pvsy-1,k) = *mpIso->lbmGetData( (int)(i*scalex), mSizey-1, (int)(k*scalez) ); 00314 } 00315 } 00316 for(int j=0;j<pvsy;j++) 00317 for(int i=0;i<pvsx;i++) { 00318 *mpPreviewSurface->lbmGetData(i,j,0) = *mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , 0); 00319 *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = *mpIso->lbmGetData( (int)(i*scalex), (int)(j*scaley) , mSizez-1); 00320 } 00321 00322 if(mFarFieldSize>=1.2) { 00323 // also remove preview border 00324 for(int k= 0; k< (pvsz-1); ++k) { 00325 for(int j=0;j< pvsy;j++) { 00326 *mpPreviewSurface->lbmGetData(0,j,k) = 00327 *mpPreviewSurface->lbmGetData(1,j,k) = 00328 *mpPreviewSurface->lbmGetData(2,j,k); 00329 *mpPreviewSurface->lbmGetData(pvsx-1,j,k) = 00330 *mpPreviewSurface->lbmGetData(pvsx-2,j,k) = 00331 *mpPreviewSurface->lbmGetData(pvsx-3,j,k); 00332 //0.0; 00333 } 00334 for(int i=0;i< pvsx;i++) { 00335 *mpPreviewSurface->lbmGetData(i,0,k) = 00336 *mpPreviewSurface->lbmGetData(i,1,k) = 00337 *mpPreviewSurface->lbmGetData(i,2,k); 00338 *mpPreviewSurface->lbmGetData(i,pvsy-1,k) = 00339 *mpPreviewSurface->lbmGetData(i,pvsy-2,k) = 00340 *mpPreviewSurface->lbmGetData(i,pvsy-3,k); 00341 //0.0; 00342 } 00343 } 00344 for(int j=0;j<pvsy;j++) 00345 for(int i=0;i<pvsx;i++) { 00346 *mpPreviewSurface->lbmGetData(i,j,0) = 00347 *mpPreviewSurface->lbmGetData(i,j,1) = 00348 *mpPreviewSurface->lbmGetData(i,j,2); 00349 *mpPreviewSurface->lbmGetData(i,j,pvsz-1) = 00350 *mpPreviewSurface->lbmGetData(i,j,pvsz-2) = 00351 *mpPreviewSurface->lbmGetData(i,j,pvsz-3); 00352 //0.0; 00353 } 00354 } 00355 } 00356 00357 // MPT 00358 #if LBM_INCLUDE_TESTSOLVERS==1 00359 mrIsoExchange(); 00360 #endif // LBM_INCLUDE_TESTSOLVERS==1 00361 00362 return; 00363 } 00364 00366 void LbmFsgrSolver::recalculateObjectSpeeds() { 00367 const bool debugRecalc = false; 00368 int numobjs = (int)(this->mpGiObjects->size()); 00369 // note - (numobjs + 1) is entry for domain settings 00370 00371 if(debugRecalc) errMsg("recalculateObjectSpeeds","start, #obj:"<<numobjs); 00372 if(numobjs>255-1) { 00373 errFatal("LbmFsgrSolver::recalculateObjectSpeeds","More than 256 objects currently not supported...",SIMWORLD_INITERROR); 00374 return; 00375 } 00376 mObjectSpeeds.resize(numobjs+1); 00377 for(int i=0; i<(int)(numobjs+0); i++) { 00378 mObjectSpeeds[i] = vec2L(this->mpParam->calculateLattVelocityFromRw( vec2P( (*this->mpGiObjects)[i]->getInitialVelocity(mSimulationTime) ))); 00379 if(debugRecalc) errMsg("recalculateObjectSpeeds","id"<<i<<" set to "<< mObjectSpeeds[i]<<", unscaled:"<< (*this->mpGiObjects)[i]->getInitialVelocity(mSimulationTime) ); 00380 } 00381 00382 // also reinit part slip values here 00383 mObjectPartslips.resize(numobjs+1); 00384 for(int i=0; i<=(int)(numobjs+0); i++) { 00385 if(i<numobjs) { 00386 mObjectPartslips[i] = (LbmFloat)(*this->mpGiObjects)[i]->getGeoPartSlipValue(); 00387 } else { 00388 // domain setting 00389 mObjectPartslips[i] = this->mDomainPartSlipValue; 00390 } 00391 LbmFloat set = mObjectPartslips[i]; 00392 00393 // as in setInfluenceVelocity 00394 const LbmFloat dt = mLevel[mMaxRefine].timestep; 00395 const LbmFloat dtInter = 0.01; 00396 //LbmFloat facFv = 1.-set; 00397 // mLevel[mMaxRefine].timestep 00398 LbmFloat facNv = (LbmFloat)( 1.-pow( (double)(set), (double)(dt/dtInter)) ); 00399 errMsg("mObjectPartslips","id:"<<i<<"/"<<numobjs<<" ts:"<<dt<< " its:"<<(dt/dtInter) <<" set"<<set<<" nv"<<facNv<<" test:"<< 00400 pow( (double)(1.-facNv),(double)(dtInter/dt)) ); 00401 mObjectPartslips[i] = facNv; 00402 00403 if(debugRecalc) errMsg("recalculateObjectSpeeds","id"<<i<<" parts "<< mObjectPartslips[i] ); 00404 } 00405 00406 if(debugRecalc) errMsg("recalculateObjectSpeeds","done, domain:"<<mObjectPartslips[numobjs]<<" n"<<numobjs); 00407 } 00408 00409 00410 00411 /*****************************************************************************/ 00413 /*****************************************************************************/ 00414 vector<ntlGeometryObject*> LbmFsgrSolver::getDebugObjects() { 00415 vector<ntlGeometryObject*> debo; 00416 if(this->mOutputSurfacePreview) { 00417 debo.push_back( mpPreviewSurface ); 00418 } 00419 #if LBM_INCLUDE_TESTSOLVERS==1 00420 if(mUseTestdata) { 00421 vector<ntlGeometryObject*> tdebo; 00422 tdebo = mpTest->getDebugObjects(); 00423 for(size_t i=0; i<tdebo.size(); i++) debo.push_back( tdebo[i] ); 00424 } 00425 #endif // ELBEEM_PLUGIN 00426 return debo; 00427 } 00428 00429 /****************************************************************************** 00430 * particle handling 00431 *****************************************************************************/ 00432 00434 int LbmFsgrSolver::initParticles() { 00435 int workSet = mLevel[mMaxRefine].setCurr; 00436 int tries = 0; 00437 int num = 0; 00438 ParticleTracer *partt = mpParticles; 00439 00440 partt->setStart( this->mvGeoStart + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) ); 00441 partt->setEnd ( this->mvGeoEnd + ntlVec3Gfx(mLevel[mMaxRefine].nodeSize*0.5) ); 00442 00443 partt->setSimStart( ntlVec3Gfx(0.0) ); 00444 partt->setSimEnd ( ntlVec3Gfx(mSizex, mSizey, getForZMaxBnd(mMaxRefine)) ); 00445 00446 while( (num<partt->getNumInitialParticles()) && (tries<100*partt->getNumInitialParticles()) ) { 00447 LbmFloat x,y,z,t; 00448 x = 1.0+(( (LbmFloat)(mSizex-3.) ) * (rand()/(RAND_MAX+1.0)) ); 00449 y = 1.0+(( (LbmFloat)(mSizey-3.) ) * (rand()/(RAND_MAX+1.0)) ); 00450 z = 1.0+(( (LbmFloat) getForZMax1(mMaxRefine)-2. )* (rand()/(RAND_MAX+1.0)) ); 00451 int i = (int)(x+0.5); 00452 int j = (int)(y+0.5); 00453 int k = (int)(z+0.5); 00454 if(LBMDIM==2) { 00455 k = 0; z = 0.5; // place in the middle of domain 00456 } 00457 00458 //if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFFluid) ) 00459 //&& ( RFLAG(mMaxRefine, i,j,k, workSet)& CFNoNbFluid ) 00460 //if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter|CFMbndInflow) ) { 00461 if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFNoBndFluid|CFUnused) ) { 00462 bool cellOk = true; 00463 //? FORDF1 { if(!(RFLAG_NB(mMaxRefine,i,j,k,workSet, l) & CFFluid)) cellOk = false; } 00464 if(!cellOk) continue; 00465 // in fluid... 00466 partt->addParticle(x,y,z); 00467 partt->getLast()->setStatus(PART_IN); 00468 partt->getLast()->setType(PART_TRACER); 00469 00470 partt->getLast()->setSize(1.); 00471 // randomize size 00472 partt->getLast()->setSize(0.5 + (rand()/(RAND_MAX+1.0))); 00473 00474 if( ( partt->getInitStart()>0.) 00475 && ( partt->getInitEnd()>0.) 00476 && ( partt->getInitEnd()>partt->getInitStart() )) { 00477 t = partt->getInitStart()+ (partt->getInitEnd()-partt->getInitStart())*(rand()/(RAND_MAX+1.0)); 00478 partt->getLast()->setLifeTime( -t ); 00479 } 00480 num++; 00481 } 00482 tries++; 00483 } // */ 00484 00485 /*FSGR_FORIJK1(mMaxRefine) { 00486 if( (RFLAG(mMaxRefine,i,j,k, workSet) & (CFNoBndFluid)) ) { 00487 LbmFloat rndn = (rand()/(RAND_MAX+1.0)); 00488 if(rndn>0.0) { 00489 ntlVec3Gfx pos( (LbmFloat)(i)-0.5, (LbmFloat)(j)-0.5, (LbmFloat)(k)-0.5 ); 00490 if(LBMDIM==2) { pos[2]=0.5; } 00491 partt->addParticle( pos[0],pos[1],pos[2] ); 00492 partt->getLast()->setStatus(PART_IN); 00493 partt->getLast()->setType(PART_TRACER); 00494 partt->getLast()->setSize(1.0); 00495 } 00496 } 00497 } // */ 00498 00499 00500 // DEBUG TEST 00501 #if LBM_INCLUDE_TESTSOLVERS==1 00502 if(mUseTestdata) { 00503 const bool partDebug=false; 00504 if(mpTest->mPartTestcase==0){ errMsg("LbmTestdata"," part init "<<mpTest->mPartTestcase); } 00505 if(mpTest->mPartTestcase==-12){ 00506 const int lev = mMaxRefine; 00507 for(int i=5;i<15;i++) { 00508 LbmFloat x,y,z; 00509 y = 0.5+(LbmFloat)(i); 00510 x = mLevel[lev].lSizex/20.0*10.0; 00511 z = mLevel[lev].lSizez/20.0*2.0; 00512 partt->addParticle(x,y,z); 00513 partt->getLast()->setStatus(PART_IN); 00514 partt->getLast()->setType(PART_BUBBLE); 00515 partt->getLast()->setSize( (-4.0+(LbmFloat)i)/1.0 ); 00516 if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() ); 00517 } 00518 } 00519 if(mpTest->mPartTestcase==-11){ 00520 const int lev = mMaxRefine; 00521 for(int i=5;i<15;i++) { 00522 LbmFloat x,y,z; 00523 y = 10.5+(LbmFloat)(i); 00524 x = mLevel[lev].lSizex/20.0*10.0; 00525 z = mLevel[lev].lSizez/20.0*40.0; 00526 partt->addParticle(x,y,z); 00527 partt->getLast()->setStatus(PART_IN); 00528 partt->getLast()->setType(PART_DROP); 00529 partt->getLast()->setSize( (-4.0+(LbmFloat)i)/1.0 ); 00530 if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() ); 00531 } 00532 } 00533 // place floats on rectangular region FLOAT_JITTER_BND 00534 if(mpTest->mPartTestcase==-10){ 00535 const int lev = mMaxRefine; 00536 const int sx = mLevel[lev].lSizex; 00537 const int sy = mLevel[lev].lSizey; 00538 //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) { 00539 //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) { 00540 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) { 00541 //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) { 00542 LbmFloat x,y,z; 00543 x = 0.0+(LbmFloat)(i); 00544 y = 0.0+(LbmFloat)(j); 00545 //z = mLevel[lev].lSizez/10.0*2.5 - 1.0; 00546 z = mLevel[lev].lSizez/20.0*9.5 - 1.0; 00547 //z = mLevel[lev].lSizez/20.0*4.5 - 1.0; 00548 partt->addParticle(x,y,z); 00549 //if( (i>0)&&(i<sx) && (j>0)&&(j<sy) ) { partt->getLast()->setStatus(PART_IN); } else { partt->getLast()->setStatus(PART_OUT); } 00550 partt->getLast()->setStatus(PART_IN); 00551 partt->getLast()->setType(PART_FLOAT); 00552 partt->getLast()->setSize( 15.0 ); 00553 if(partDebug) errMsg("PARTTT","SET "<<PRINT_VEC(x,y,z)<<" p"<<partt->getLast()->getPos() <<" s"<<partt->getLast()->getSize() ); 00554 } 00555 } } 00556 } 00557 // DEBUG TEST 00558 #endif // LBM_INCLUDE_TESTSOLVERS 00559 00560 00561 debMsgStd("LbmFsgrSolver::initParticles",DM_MSG,"Added "<<num<<" particles, genProb:"<<this->mPartGenProb<<", tries:"<<tries, 10); 00562 if(num != partt->getNumParticles()) return 1; 00563 00564 return 0; 00565 } 00566 00567 // helper function for particle debugging 00568 /*static string getParticleStatusString(int state) { 00569 std::ostringstream out; 00570 if(state&PART_DROP) out << "dropp "; 00571 if(state&PART_TRACER) out << "tracr "; 00572 if(state&PART_BUBBLE) out << "bubbl "; 00573 if(state&PART_FLOAT) out << "float "; 00574 if(state&PART_INTER) out << "inter "; 00575 00576 if(state&PART_IN) out << "inn "; 00577 if(state&PART_OUT) out << "out "; 00578 if(state&PART_INACTIVE) out << "INACT "; 00579 if(state&PART_OUTFLUID) out << "outfluid "; 00580 return out.str(); 00581 } // */ 00582 00583 #define P_CHANGETYPE(p, newtype) \ 00584 p->setLifeTime(0.); \ 00585 /* errMsg("PIT","U pit"<<(p)->getId()<<" pos:"<< (p)->getPos()<<" status:"<<convertFlags2String((p)->getFlags())<<" to "<< (newtype) ); */ \ 00586 p->setType(newtype); 00587 00588 // tracer defines 00589 #define TRACE_JITTER 0.025 00590 #define TRACE_RAND (rand()/(RAND_MAX+1.0))*TRACE_JITTER-(TRACE_JITTER*0.5) 00591 #define FFGET_NORM(var,dl) \ 00592 if(RFLAG_NB(lev,i,j,k,workSet, dl) &(CFInter)){ (var) = QCELL_NB(lev,i,j,k,workSet,dl,dFfrac); } \ 00593 else if(RFLAG_NB(lev,i,j,k,workSet, dl) &(CFFluid|CFUnused)){ (var) = 1.; } else (var) = 0.0; 00594 00595 // float jitter 00596 #define FLOAT_JITTER_BND (FLOAT_JITTER*2.0) 00597 #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)) 00598 00599 #define DEL_PART { \ 00600 /*errMsg("PIT","DEL AT "<< __LINE__<<" type:"<<p->getType()<<" "); */ \ 00601 p->setActive( false ); \ 00602 continue; } 00603 00604 void LbmFsgrSolver::advanceParticles() { 00605 const int level = mMaxRefine; 00606 const int workSet = mLevel[level].setCurr; 00607 LbmFloat vx=0.0,vy=0.0,vz=0.0; 00608 //int debugOutCounter=0; // debug output counter 00609 00610 myTime_t parttstart = getTime(); 00611 const LbmFloat cellsize = this->mpParam->getCellSize(); 00612 const LbmFloat timestep = this->mpParam->getTimestep(); 00613 //const LbmFloat viscAir = 1.79 * 1e-5; // RW2L kin. viscosity, mu 00614 //const LbmFloat viscWater = 1.0 * 1e-6; // RW2L kin. viscosity, mu 00615 const LbmFloat rhoAir = 1.2; // [kg m^-3] RW2L 00616 const LbmFloat rhoWater = 1000.0; // RW2L 00617 const LbmFloat minDropSize = 0.0005; // [m], = 2mm RW2L 00618 const LbmVec velAir(0.); // [m / s] 00619 00620 const LbmFloat r1 = 0.005; // r max 00621 const LbmFloat r2 = 0.0005; // r min 00622 const LbmFloat v1 = 9.0; // v max 00623 const LbmFloat v2 = 2.0; // v min 00624 const LbmVec rwgrav = vec2L( this->mpParam->getGravity(mSimulationTime) ); 00625 const bool useff = (mFarFieldSize>1.2); // if(mpTest->mFarfMode>0){ 00626 00627 // TODO scale bubble/part damping dep. on timestep, also drop bnd rev damping 00628 const int cutval = mCutoff; // use full border!? 00629 if(this->mStepCnt%50==49) { mpParticles->cleanup(); } 00630 for(vector<ParticleObject>::iterator pit= mpParticles->getParticlesBegin(); 00631 pit!= mpParticles->getParticlesEnd(); pit++) { 00632 //if((*pit).getPos()[2]>10.) errMsg("PIT"," pit"<<(*pit).getId()<<" pos:"<< (*pit).getPos()<<" status:["<<getParticleStatusString((*pit).getFlags())<<"] vel:"<< (*pit).getVel() ); 00633 if( (*pit).getActive()==false ) continue; 00634 // skip until reached 00635 ParticleObject *p = &(*pit); 00636 if(p->getLifeTime()<0.){ 00637 if(p->getLifeTime() < -mSimulationTime) continue; 00638 else p->setLifeTime(-mLevel[level].timestep); // zero for following update 00639 } 00640 int i,j,k; 00641 p->setLifeTime(p->getLifeTime()+mLevel[level].timestep); 00642 00643 // nearest neighbor, particle positions don't include empty bounds 00644 ntlVec3Gfx pos = p->getPos(); 00645 i= (int)pos[0]; j= (int)pos[1]; k= (int)pos[2];// no offset necessary 00646 if(LBMDIM==2) { k = 0; } 00647 00648 // only testdata handling, all for sws 00649 #if LBM_INCLUDE_TESTSOLVERS==1 00650 if(useff && (mpTest->mFarfMode>0)) { 00651 p->setStatus(PART_OUT); 00652 mpTest->handleParticle(p, i,j,k); continue; 00653 } 00654 #endif // LBM_INCLUDE_TESTSOLVERS==1 00655 00656 // in out tests 00657 if(p->getStatus()&PART_IN) { // IN 00658 if( (i<cutval)||(i>mSizex-1-cutval)|| 00659 (j<cutval)||(j>mSizey-1-cutval) 00660 //||(k<cutval)||(k>mSizez-1-cutval) 00661 ) { 00662 if(!useff) { DEL_PART; 00663 } else { 00664 p->setStatus(PART_OUT); 00665 } 00666 } 00667 } else { // OUT rough check 00668 // check in again? 00669 if( (i>=cutval)&&(i<=mSizex-1-cutval)&& 00670 (j>=cutval)&&(j<=mSizey-1-cutval) 00671 ) { 00672 p->setStatus(PART_IN); 00673 } 00674 } 00675 00676 if( (p->getType()==PART_BUBBLE) || 00677 (p->getType()==PART_TRACER) ) { 00678 00679 // no interpol 00680 vx = vy = vz = 0.0; 00681 if(p->getStatus()&PART_IN) { // IN 00682 if(k>=cutval) { 00683 if(k>mSizez-1-cutval) DEL_PART; 00684 00685 if( RFLAG(level, i,j,k, workSet)&(CFFluid|CFUnused) ) { 00686 // still ok 00687 int partLev = level; 00688 int si=i, sj=j, sk=k; 00689 while(partLev>0 && RFLAG(partLev, si,sj,sk, workSet)&(CFUnused)) { 00690 partLev--; 00691 si/=2; 00692 sj/=2; 00693 sk/=2; 00694 } 00695 // get velocity from fluid cell 00696 if( RFLAG(partLev, si,sj,sk, workSet)&(CFFluid) ) { 00697 LbmFloat *ccel = RACPNT(partLev, si,sj,sk, workSet); 00698 FORDF1{ 00699 LbmFloat cdf = RAC(ccel, l); 00700 // TODO update below 00701 vx += (this->dfDvecX[l]*cdf); 00702 vy += (this->dfDvecY[l]*cdf); 00703 vz += (this->dfDvecZ[l]*cdf); 00704 } 00705 // remove gravity influence 00706 const LbmFloat lesomega = mLevel[level].omega; // no les 00707 vx -= mLevel[level].gravity[0] * lesomega*0.5; 00708 vy -= mLevel[level].gravity[1] * lesomega*0.5; 00709 vz -= mLevel[level].gravity[2] * lesomega*0.5; 00710 } // fluid vel 00711 00712 } else { // OUT 00713 // out of bounds, deactivate... 00714 // FIXME make fsgr treatment 00715 if(p->getType()==PART_BUBBLE) { P_CHANGETYPE(p, PART_FLOAT ); continue; } 00716 } 00717 } else { 00718 // below 3d region, just rise 00719 } 00720 } else { // OUT 00721 # if LBM_INCLUDE_TESTSOLVERS==1 00722 if(useff) { mpTest->handleParticle(p, i,j,k); } 00723 else DEL_PART; 00724 # else // LBM_INCLUDE_TESTSOLVERS==1 00725 DEL_PART; 00726 # endif // LBM_INCLUDE_TESTSOLVERS==1 00727 // TODO use x,y vel...? 00728 } 00729 00730 ntlVec3Gfx v = p->getVel(); // dampen... 00731 if( (useff)&& (p->getType()==PART_BUBBLE) ) { 00732 // test rise 00733 00734 if(mPartUsePhysModel) { 00735 LbmFloat radius = p->getSize() * minDropSize; 00736 LbmVec velPart = vec2L(p->getVel()) *cellsize/timestep; // L2RW, lattice velocity 00737 LbmVec velWater = LbmVec(vx,vy,vz) *cellsize/timestep;// L2RW, fluid velocity 00738 LbmVec velRel = velWater - velPart; 00739 //LbmFloat velRelNorm = norm(velRel); 00740 LbmFloat pvolume = rhoAir * 4.0/3.0 * M_PI* radius*radius*radius; // volume: 4/3 pi r^3 00741 00742 LbmVec fb = -rwgrav* pvolume *rhoWater; 00743 LbmVec fd = velRel*6.0*M_PI*radius* (1e-3); //viscWater; 00744 LbmVec change = (fb+fd) *10.0*timestep *(timestep/cellsize); 00745 /*if(debugOutCounter<0) { 00746 errMsg("PIT","BTEST1 vol="<<pvolume<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel); 00747 errMsg("PIT","BTEST2 cellsize="<<cellsize<<" timestep="<<timestep<<" viscW="<<viscWater<<" ss/mb="<<(timestep/(pvolume*rhoAir))); 00748 errMsg("PIT","BTEST2 grav="<<rwgrav<<" " ); 00749 errMsg("PIT","BTEST2 change="<<(change)<<" fb="<<(fb)<<" fd="<<(fd)<<" "); 00750 errMsg("PIT","BTEST2 change="<<norm(change)<<" fb="<<norm(fb)<<" fd="<<norm(fd)<<" "); 00751 } // DEBUG */ 00752 00753 LbmVec fd2 = (LbmVec(vx,vy,vz)-vec2L(p->getVel())) * 6.0*M_PI*radius* (1e-3); //viscWater; 00754 LbmFloat w = 0.99; 00755 vz = (1.0-w)*vz + w*(p->getVel()[2]-0.5*(p->getSize()/5.0)*mLevel[level].gravity[2]); 00756 v = ntlVec3Gfx(vx,vy,vz)+vec2G(fd2); 00757 p->setVel( v ); 00758 } else { 00759 // non phys, half old, half fluid, use slightly slower acc 00760 v = v*0.5 + ntlVec3Gfx(vx,vy,vz)* 0.5-vec2G(mLevel[level].gravity)*0.5; 00761 p->setVel( v * 0.99 ); 00762 } 00763 p->advanceVel(); 00764 00765 } else if(p->getType()==PART_TRACER) { 00766 v = ntlVec3Gfx(vx,vy,vz); 00767 CellFlagType fflag = RFLAG(level, i,j,k, workSet); 00768 00769 if(fflag&(CFFluid|CFInter) ) { p->setInFluid(true); 00770 } else { p->setInFluid(false); } 00771 00772 if( (( fflag&CFFluid ) && ( fflag&CFNoBndFluid )) || 00773 (( fflag&CFInter ) && (!(fflag&CFNoNbFluid)) ) ) { 00774 // only real fluid 00775 # if LBMDIM==3 00776 p->advance( TRACE_RAND,TRACE_RAND,TRACE_RAND); 00777 # else 00778 p->advance( TRACE_RAND,TRACE_RAND, 0.); 00779 # endif 00780 00781 } else { 00782 // move inwards along normal, make sure normal is valid first 00783 // todo use class funcs! 00784 const int lev = level; 00785 LbmFloat nx=0.,ny=0.,nz=0., nv1,nv2; 00786 bool nonorm = false; 00787 if(i<=0) { nx = -1.; nonorm = true; } 00788 if(i>=mSizex-1) { nx = 1.; nonorm = true; } 00789 if(j<=0) { ny = -1.; nonorm = true; } 00790 if(j>=mSizey-1) { ny = 1.; nonorm = true; } 00791 # if LBMDIM==3 00792 if(k<=0) { nz = -1.; nonorm = true; } 00793 if(k>=mSizez-1) { nz = 1.; nonorm = true; } 00794 # endif // LBMDIM==3 00795 if(!nonorm) { 00796 FFGET_NORM(nv1,dE); FFGET_NORM(nv2,dW); 00797 nx = 0.5* (nv2-nv1); 00798 FFGET_NORM(nv1,dN); FFGET_NORM(nv2,dS); 00799 ny = 0.5* (nv2-nv1); 00800 # if LBMDIM==3 00801 FFGET_NORM(nv1,dT); FFGET_NORM(nv2,dB); 00802 nz = 0.5* (nv2-nv1); 00803 # else // LBMDIM==3 00804 nz = 0.; 00805 # endif // LBMDIM==3 00806 } else { 00807 v = p->getVel() + vec2G(mLevel[level].gravity); 00808 } 00809 p->advanceVec( (ntlVec3Gfx(nx,ny,nz)) * -0.1 ); // + vec2G(mLevel[level].gravity); 00810 } 00811 } 00812 00813 p->setVel( v ); 00814 p->advanceVel(); 00815 } 00816 00817 // drop handling 00818 else if(p->getType()==PART_DROP) { 00819 ntlVec3Gfx v = p->getVel(); // dampen... 00820 00821 if(mPartUsePhysModel) { 00822 LbmFloat radius = p->getSize() * minDropSize; 00823 LbmVec velPart = vec2L(p->getVel()) *cellsize /timestep; // * cellsize / timestep; // L2RW, lattice velocity 00824 LbmVec velRel = velAir - velPart; 00825 //LbmVec velRelLat = velRel /cellsize*timestep; // L2RW 00826 LbmFloat velRelNorm = norm(velRel); 00827 // TODO calculate values in lattice units, compute CD?!??! 00828 LbmFloat mb = rhoWater * 4.0/3.0 * M_PI* radius*radius*radius; // mass: 4/3 pi r^3 rho 00829 const LbmFloat rw = (r1-radius)/(r1-r2); 00830 const LbmFloat rmax = (0.5 + 0.5*rw); 00831 const LbmFloat vmax = (v2 + (v1-v2)* (1.0-rw) ); 00832 const LbmFloat cd = (rmax) * (velRelNorm)/(vmax); 00833 00834 LbmVec fg = rwgrav * mb;// * (1.0-rhoAir/rhoWater); 00835 LbmVec fd = velRel* velRelNorm* cd*M_PI *rhoAir *0.5 *radius*radius; 00836 LbmVec change = (fg+ fd ) *timestep / mb *(timestep/cellsize); 00837 //if(k>0) { errMsg("\nPIT","NTEST1 mb="<<mb<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel<<" pgetVel="<<p->getVel() ); } 00838 00839 v += vec2G(change); 00840 p->setVel(v); 00841 // NEW 00842 } else { 00843 p->setVel( v ); 00844 int gravk = (int)(p->getPos()[2]+mLevel[level].gravity[2]); 00845 if(gravk>=0 && gravk<mSizez && RFLAG(level, i,j,gravk, workSet)&CFBnd) { 00846 // dont add for "resting" parts 00847 v[2] = 0.; 00848 p->setVel( v*0.9 ); // restdamping 00849 } else { 00850 p->addToVel( vec2G(mLevel[level].gravity) ); 00851 } 00852 } // OLD 00853 p->advanceVel(); 00854 00855 if(p->getStatus()&PART_IN) { // IN 00856 if(k<cutval) { DEL_PART; continue; } 00857 if(k<=mSizez-1-cutval){ 00858 CellFlagType pflag = RFLAG(level, i,j,k, workSet); 00859 //errMsg("PIT move"," at "<<PRINT_IJK<<" flag"<<convertCellFlagType2String(pflag) ); 00860 if(pflag & (CFBnd)) { 00861 handleObstacleParticle(p); 00862 continue; 00863 } else if(pflag & (CFEmpty)) { 00864 // still ok 00865 } else if((pflag & CFInter) 00866 //&&(!(RFLAG(level, i,j,k, workSet)& CFNoNbFluid)) 00867 ) { 00868 // add to no nb fluid i.f.'s, so skip if interface with fluid nb 00869 } else if(pflag & (CFFluid|CFUnused|CFInter) ){ // interface cells ignored here due to previous check! 00870 // add dropmass again, (these are only interf. with nonbfl.) 00871 int oi= (int)(pos[0]-1.25*v[0]+0.5); 00872 int oj= (int)(pos[1]-1.25*v[1]+0.5); 00873 int ok= (int)(pos[2]-1.25*v[2]+0.5); 00874 const LbmFloat size = p->getSize(); 00875 const LbmFloat dropmass = ParticleObject::getMass(mPartDropMassSub*size); 00876 bool orgcellok = false; 00877 if( (oi<0)||(oi>mSizex-1)|| 00878 (oj<0)||(oj>mSizey-1)|| 00879 (ok<0)||(ok>mSizez-1) ) { 00880 // org cell not ok! 00881 } else if( RFLAG(level, oi,oj,ok, workSet) & (CFInter) ){ 00882 orgcellok = true; 00883 } else { 00884 // search upward for interface 00885 oi=i; oj=j; ok=k; 00886 for(int kk=0; kk<5 && ok<=mSizez-2; kk++) { 00887 ok++; // check sizez-2 due to this increment! 00888 if( RFLAG(level, oi,oj,ok, workSet) & (CFInter) ){ 00889 kk = 5; orgcellok = true; 00890 } 00891 } 00892 } 00893 00894 //errMsg("PTIMPULSE"," new v"<<v<<" at "<<PRINT_VEC(oi,oj,ok)<<" , was "<<PRINT_VEC(i,j,k)<<" ok "<<orgcellok ); 00895 if(orgcellok) { 00896 QCELL(level, oi,oj,ok, workSet, dMass) += dropmass; 00897 QCELL(level, oi,oj,ok, workSet, dFfrac) += dropmass; // assume rho=1? 00898 00899 if(RFLAG(level, oi,oj,ok, workSet) & CFNoBndFluid){ 00900 // check speed, perhaps normalize 00901 gfxReal vlensqr = normNoSqrt(v); 00902 if(vlensqr > 0.166*0.166) { 00903 v *= 1./sqrtf((float)vlensqr)*0.166; 00904 } 00905 // compute cell velocity 00906 LbmFloat *tcel = RACPNT(level, oi,oj,ok, workSet); 00907 LbmFloat velUx=0., velUy=0., velUz=0.; 00908 FORDF0 { 00909 velUx += (this->dfDvecX[l]*RAC(tcel,l)); 00910 velUy += (this->dfDvecY[l]*RAC(tcel,l)); 00911 velUz += (this->dfDvecZ[l]*RAC(tcel,l)); 00912 } 00913 // add impulse 00914 /* 00915 LbmFloat cellVelSqr = velUx*velUx+ velUy*velUy+ velUz*velUz; 00916 //errMsg("PTIMPULSE"," new v"<<v<<" cvs"<<cellVelSqr<<"="<<sqrt(cellVelSqr)); 00917 if(cellVelSqr< 0.166*0.166) { 00918 FORDF1 { 00919 const LbmFloat add = 3. * dropmass * this->dfLength[l]*(v[0]*this->dfDvecX[l]+v[1]*this->dfDvecY[l]+v[2]*this->dfDvecZ[l]); 00920 RAC(tcel,l) += add; 00921 } } // */ 00922 } // only add impulse away from obstacles! 00923 } // orgcellok 00924 00925 // FIXME make fsgr treatment 00926 P_CHANGETYPE(p, PART_FLOAT ); continue; 00927 // jitter in cell to prevent stacking when hitting a steep surface 00928 ntlVec3Gfx cpos = p->getPos(); 00929 cpos[0] += (rand()/(RAND_MAX+1.0))-0.5; 00930 cpos[1] += (rand()/(RAND_MAX+1.0))-0.5; 00931 cpos[2] += (rand()/(RAND_MAX+1.0))-0.5; 00932 p->setPos(cpos); 00933 } else { 00934 DEL_PART; 00935 this->mNumParticlesLost++; 00936 } 00937 } 00938 } else { // OUT 00939 # if LBM_INCLUDE_TESTSOLVERS==1 00940 if(useff) { mpTest->handleParticle(p, i,j,k); } 00941 else{ DEL_PART; } 00942 # else // LBM_INCLUDE_TESTSOLVERS==1 00943 DEL_PART; 00944 # endif // LBM_INCLUDE_TESTSOLVERS==1 00945 } 00946 00947 } // air particle 00948 00949 // inter particle 00950 else if(p->getType()==PART_INTER) { 00951 // unused!? 00952 if(p->getStatus()&PART_IN) { // IN 00953 if((k<cutval)||(k>mSizez-1-cutval)) { 00954 // undecided particle above or below... remove? 00955 DEL_PART; 00956 } 00957 00958 CellFlagType pflag = RFLAG(level, i,j,k, workSet); 00959 if(pflag& CFInter ) { 00960 // still ok 00961 } else if(pflag& (CFFluid|CFUnused) ) { 00962 P_CHANGETYPE(p, PART_FLOAT ); continue; 00963 } else if(pflag& CFEmpty ) { 00964 P_CHANGETYPE(p, PART_DROP ); continue; 00965 } else if(pflag& CFBnd ) { 00966 P_CHANGETYPE(p, PART_FLOAT ); continue; 00967 } 00968 } else { // OUT 00969 // undecided particle outside... remove? 00970 DEL_PART; 00971 } 00972 } 00973 00974 // float particle 00975 else if(p->getType()==PART_FLOAT) { 00976 00977 if(p->getStatus()&PART_IN) { // IN 00978 if(k<cutval) DEL_PART; 00979 // not valid for mass... 00980 vx = vy = vz = 0.0; 00981 00982 // define from particletracer.h 00983 #if MOVE_FLOATS==1 00984 const int DEPTH_AVG=3; // only average interface vels 00985 int ccnt=0; 00986 for(int kk=0;kk<DEPTH_AVG;kk+=1) { 00987 if((k-kk)<1) continue; 00988 if(RFLAG(level, i,j,k, workSet)&(CFInter)) {} else continue; 00989 ccnt++; 00990 FORDF1{ 00991 LbmFloat cdf = QCELL(level, i,j,k-kk, workSet, l); 00992 vx += (this->dfDvecX[l]*cdf); 00993 vy += (this->dfDvecY[l]*cdf); 00994 vz += (this->dfDvecZ[l]*cdf); 00995 } 00996 } 00997 if(ccnt) { 00998 // use halved surface velocity (todo, use omega instead) 00999 vx /=(LbmFloat)(ccnt * 2.0); // half xy speed! value2 01000 vy /=(LbmFloat)(ccnt * 2.0); 01001 vz /=(LbmFloat)(ccnt); } 01002 #else // MOVE_FLOATS==1 01003 vx=vy=0.; //p->setVel(ntlVec3Gfx(0.) ); // static_float 01004 #endif // MOVE_FLOATS==1 01005 vx += (rand()/(RAND_MAX+1.0))*(FLOAT_JITTER*0.2)-(FLOAT_JITTER*0.2*0.5); 01006 vy += (rand()/(RAND_MAX+1.0))*(FLOAT_JITTER*0.2)-(FLOAT_JITTER*0.2*0.5); 01007 01008 //bool delfloat = false; 01009 if( ( RFLAG(level, i,j,k, workSet)& (CFFluid|CFUnused) ) ) { 01010 // in fluid cell 01011 vz = p->getVel()[2]-1.0*mLevel[level].gravity[2]; // simply rise... 01012 if(vz<0.) vz=0.; 01013 } else if( ( RFLAG(level, i,j,k, workSet)& CFBnd ) ) { 01014 // force downwards movement, move below obstacle... 01015 //vz = p->getVel()[2]+1.0*mLevel[level].gravity[2]; // fall... 01016 //if(vz>0.) vz=0.; 01017 DEL_PART; 01018 } else if( ( RFLAG(level, i,j,k, workSet)& CFInter ) ) { 01019 // keep in interface , one grid cell offset is added in part. gen 01020 } else { // all else... 01021 if( ( RFLAG(level, i,j,k-1, workSet)& (CFFluid|CFInter) ) ) { 01022 vz = p->getVel()[2]+2.0*mLevel[level].gravity[2]; // fall... 01023 if(vz>0.) vz=0.; } 01024 else { DEL_PART; } 01025 } 01026 01027 p->setVel( vec2G( ntlVec3Gfx(vx,vy,vz) ) ); //? 01028 p->advanceVel(); 01029 } else { 01030 #if LBM_INCLUDE_TESTSOLVERS==1 01031 if(useff) { mpTest->handleParticle(p, i,j,k); } 01032 else DEL_PART; 01033 #else // LBM_INCLUDE_TESTSOLVERS==1 01034 DEL_PART; 01035 #endif // LBM_INCLUDE_TESTSOLVERS==1 01036 } 01037 01038 // additional bnd jitter 01039 if((0) && (useff) && (p->getLifeTime()<3.*mLevel[level].timestep)) { 01040 // use half butoff border 1/8 01041 int maxdw = (int)(mLevel[level].lSizex*0.125*0.5); 01042 if(maxdw<3) maxdw=3; 01043 if((j>=0)&&(j<=mSizey-1)) { 01044 if(ABS(i-( cutval))<maxdw) { p->advance( FLOAT_JITTBNDRAND( ABS(i-( cutval))), 0.,0.); } 01045 if(ABS(i-(mSizex-1-cutval))<maxdw) { p->advance( FLOAT_JITTBNDRAND( ABS(i-(mSizex-1-cutval))), 0.,0.); } 01046 } 01047 } 01048 } // PART_FLOAT 01049 01050 // unknown particle type 01051 else { 01052 errMsg("LbmFsgrSolver::advanceParticles","PIT pit invalid type!? "<<p->getStatus() ); 01053 } 01054 } 01055 myTime_t parttend = getTime(); 01056 debMsgStd("LbmFsgrSolver::advanceParticles",DM_MSG,"Time for particle update:"<< getTimeString(parttend-parttstart)<<", #particles:"<<mpParticles->getNumParticles() , 10 ); 01057 } 01058 01059 void LbmFsgrSolver::notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename) { 01060 int workSet = mLevel[mMaxRefine].setCurr; 01061 std::ostringstream name; 01062 01063 // debug - raw dump of ffrac values, as text! 01064 if(mDumpRawText) { 01065 name << outfilename<< frameNrStr <<".dump"; 01066 FILE *file = fopen(name.str().c_str(),"w"); 01067 if(file) { 01068 01069 for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k) { 01070 for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++) { 01071 for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) { 01072 float val = 0.; 01073 if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) { 01074 val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac); 01075 if(val<0.) val=0.; 01076 if(val>1.) val=1.; 01077 } 01078 if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.; 01079 fprintf(file, "%f ",val); // text 01080 //errMsg("W", PRINT_IJK<<" val:"<<val); 01081 } 01082 fprintf(file, "\n"); // text 01083 } 01084 fprintf(file, "\n"); // text 01085 } 01086 fclose(file); 01087 01088 } // file 01089 } // */ 01090 01091 if(mDumpRawBinary) { 01092 if(!mDumpRawBinaryZip) { 01093 // unzipped, only fill 01094 name << outfilename<< frameNrStr <<".bdump"; 01095 FILE *file = fopen(name.str().c_str(),"w"); 01096 if(file) { 01097 for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k) { 01098 for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++) { 01099 for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) { 01100 float val = 0.; 01101 if(RFLAG(mMaxRefine, i,j,k, workSet) & CFInter) { 01102 val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac); 01103 if(val<0.) val=0.; 01104 if(val>1.) val=1.; 01105 } 01106 if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) val = 1.; 01107 fwrite( &val, sizeof(val), 1, file); // binary 01108 } 01109 } 01110 } 01111 fclose(file); 01112 } // file 01113 } // unzipped 01114 else { 01115 // zipped, use iso values 01116 prepareVisualization(); 01117 name << outfilename<< frameNrStr <<".bdump.gz"; 01118 gzFile gzf = gzopen(name.str().c_str(),"wb9"); 01119 if(gzf) { 01120 // write size 01121 int s; 01122 s=mSizex; gzwrite(gzf, &s, sizeof(s)); 01123 s=mSizey; gzwrite(gzf, &s, sizeof(s)); 01124 s=mSizez; gzwrite(gzf, &s, sizeof(s)); 01125 01126 // write isovalues 01127 for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k) { 01128 for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) { 01129 for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) { 01130 float val = 0.; 01131 val = *mpIso->lbmGetData( i,j,k ); 01132 gzwrite(gzf, &val, sizeof(val)); 01133 } 01134 } 01135 } 01136 gzclose(gzf); 01137 } // gzf 01138 } // zip 01139 } // bin dump 01140 01141 dumptype = 0; frameNr = 0; // get rid of warning 01142 } 01143 01145 void LbmFsgrSolver::handleObstacleParticle(ParticleObject *p) { 01146 //if(normNoSqrt(v)<=0.) continue; // skip stuck 01147 /* 01148 p->setVel( v * -1. ); // revert 01149 p->advanceVel(); // move back twice... 01150 if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFBndNoslip)) { 01151 p->setVel( v * -0.5 ); // revert & dampen 01152 } 01153 p->advanceVel(); 01154 // */ 01155 // TODO mark/remove stuck parts!? 01156 01157 const int level = mMaxRefine; 01158 const int workSet = mLevel[level].setCurr; 01159 LbmVec v = vec2L( p->getVel() ); 01160 if(normNoSqrt(v)<=0.) { 01161 p->setVel(vec2G(mLevel[level].gravity)); 01162 } 01163 01164 CellFlagType pflag = CFBnd; 01165 ntlVec3Gfx posOrg(p->getPos()); 01166 ntlVec3Gfx npos(0.); 01167 int ni=1,nj=1,nk=1; 01168 int tries = 0; 01169 01170 // try to undo movement 01171 p->advanceVec( (p->getVel()-vec2G(mLevel[level].gravity)) * -2.); 01172 01173 npos = p->getPos(); ni= (int)npos[0]; 01174 nj= (int)npos[1]; nk= (int)npos[2]; 01175 if(LBMDIM==2) { nk = 0; } 01176 //errMsg("BOUNDCPAR"," t"<<PRINT_VEC(ni,nj,nk)<<" v"<<v<<" p"<<npos); 01177 01178 // delete out of domain 01179 if(!checkDomainBounds(level,ni,nj,nk)) { 01180 //errMsg("BOUNDCPAR"," DEL! "); 01181 p->setActive( false ); 01182 return; 01183 } 01184 pflag = RFLAG(level, ni,nj,nk, workSet); 01185 01186 // try to force particle out of boundary 01187 bool haveNorm = false; 01188 LbmVec bnormal; 01189 if(pflag&CFBnd) { 01190 npos = posOrg; ni= (int)npos[0]; 01191 nj= (int)npos[1]; nk= (int)npos[2]; 01192 if(LBMDIM==2) { nk = 0; } 01193 01194 computeObstacleSurfaceNormalAcc(ni,nj,nk, &bnormal[0]); 01195 haveNorm = true; 01196 normalize(bnormal); 01197 bnormal *= 0.25; 01198 01199 tries = 1; 01200 while(pflag&CFBnd && tries<=5) { 01201 // use increasing step sizes 01202 p->advanceVec( vec2G( bnormal *0.5 *(gfxReal)tries ) ); 01203 npos = p->getPos(); 01204 ni= (int)npos[0]; 01205 nj= (int)npos[1]; 01206 nk= (int)npos[2]; 01207 01208 // delete out of domain 01209 if(!checkDomainBounds(level,ni,nj,nk)) { 01210 //errMsg("BOUNDCPAR"," DEL! "); 01211 p->setActive( false ); 01212 return; 01213 } 01214 pflag = RFLAG(level, ni,nj,nk, workSet); 01215 tries++; 01216 } 01217 01218 // really stuck, delete... 01219 if(pflag&CFBnd) { 01220 p->setActive( false ); 01221 return; 01222 } 01223 } 01224 01225 // not in bound anymore! 01226 if(!haveNorm) { 01227 CellFlagType *bflag = &RFLAG(level, ni,nj,nk, workSet); 01228 LbmFloat *bcell = RACPNT(level, ni,nj,nk, workSet); 01229 computeObstacleSurfaceNormal(bcell,bflag, &bnormal[0]); 01230 } 01231 normalize(bnormal); 01232 LbmVec normComp = bnormal * dot(vec2L(v),bnormal); 01233 //errMsg("BOUNDCPAR","bnormal"<<bnormal<<" normComp"<<normComp<<" newv"<<(v-normComp) ); 01234 v = (v-normComp)*0.9; // only move tangential 01235 v *= 0.9; // restdamping , todo use timestep 01236 p->setVel(vec2G(v)); 01237 p->advanceVel(); 01238 } 01239 01240 /*****************************************************************************/ 01242 /*****************************************************************************/ 01243 void 01244 LbmFsgrSolver::printLbmCell(int level, int i, int j, int k, int set) { 01245 stdCellId *newcid = new stdCellId; 01246 newcid->level = level; 01247 newcid->x = i; 01248 newcid->y = j; 01249 newcid->z = k; 01250 01251 // this function is not called upon clicking, then its from setMouseClick 01252 debugPrintNodeInfo( newcid, set ); 01253 delete newcid; 01254 } 01255 void 01256 LbmFsgrSolver::debugMarkCellCall(int level, int vi,int vj,int vk) { 01257 stdCellId *newcid = new stdCellId; 01258 newcid->level = level; 01259 newcid->x = vi; 01260 newcid->y = vj; 01261 newcid->z = vk; 01262 this->addCellToMarkedList( newcid ); 01263 } 01264 01265 01266 /*****************************************************************************/ 01267 // implement CellIterator<UniformFsgrCellIdentifier> interface 01268 /*****************************************************************************/ 01269 01270 01271 01272 // values from guiflkt.cpp 01273 extern double guiRoiSX, guiRoiSY, guiRoiSZ, guiRoiEX, guiRoiEY, guiRoiEZ; 01274 extern int guiRoiMaxLev, guiRoiMinLev; 01275 #define CID_SX (int)( (mLevel[cid->level].lSizex-1) * guiRoiSX ) 01276 #define CID_SY (int)( (mLevel[cid->level].lSizey-1) * guiRoiSY ) 01277 #define CID_SZ (int)( (mLevel[cid->level].lSizez-1) * guiRoiSZ ) 01278 01279 #define CID_EX (int)( (mLevel[cid->level].lSizex-1) * guiRoiEX ) 01280 #define CID_EY (int)( (mLevel[cid->level].lSizey-1) * guiRoiEY ) 01281 #define CID_EZ (int)( (mLevel[cid->level].lSizez-1) * guiRoiEZ ) 01282 01283 CellIdentifierInterface* 01284 LbmFsgrSolver::getFirstCell( ) { 01285 int level = mMaxRefine; 01286 01287 #if LBMDIM==3 01288 if(mMaxRefine>0) { level = mMaxRefine-1; } // NO1HIGHESTLEV DEBUG 01289 #endif 01290 level = guiRoiMaxLev; 01291 if(level>mMaxRefine) level = mMaxRefine; 01292 01293 //errMsg("LbmFsgrSolver::getFirstCell","Celliteration started..."); 01294 stdCellId *cid = new stdCellId; 01295 cid->level = level; 01296 cid->x = CID_SX; 01297 cid->y = CID_SY; 01298 cid->z = CID_SZ; 01299 return cid; 01300 } 01301 01302 LbmFsgrSolver::stdCellId* 01303 LbmFsgrSolver::convertBaseCidToStdCid( CellIdentifierInterface* basecid) { 01304 //stdCellId *cid = dynamic_cast<stdCellId*>( basecid ); 01305 stdCellId *cid = (stdCellId*)( basecid ); 01306 return cid; 01307 } 01308 01309 void LbmFsgrSolver::advanceCell( CellIdentifierInterface* basecid) { 01310 stdCellId *cid = convertBaseCidToStdCid(basecid); 01311 if(cid->getEnd()) return; 01312 01313 //debugOut(" ADb "<<cid->x<<","<<cid->y<<","<<cid->z<<" e"<<cid->getEnd(), 10); 01314 cid->x++; 01315 if(cid->x > CID_EX){ cid->x = CID_SX; cid->y++; 01316 if(cid->y > CID_EY){ cid->y = CID_SY; cid->z++; 01317 if(cid->z > CID_EZ){ 01318 cid->level--; 01319 cid->x = CID_SX; 01320 cid->y = CID_SY; 01321 cid->z = CID_SZ; 01322 if(cid->level < guiRoiMinLev) { 01323 cid->level = guiRoiMaxLev; 01324 cid->setEnd( true ); 01325 } 01326 } 01327 } 01328 } 01329 //debugOut(" ADa "<<cid->x<<","<<cid->y<<","<<cid->z<<" e"<<cid->getEnd(), 10); 01330 } 01331 01332 bool LbmFsgrSolver::noEndCell( CellIdentifierInterface* basecid) { 01333 stdCellId *cid = convertBaseCidToStdCid(basecid); 01334 return (!cid->getEnd()); 01335 } 01336 01337 void LbmFsgrSolver::deleteCellIterator( CellIdentifierInterface** cid ) { 01338 delete *cid; 01339 *cid = NULL; 01340 } 01341 01342 CellIdentifierInterface* LbmFsgrSolver::getCellAt( ntlVec3Gfx pos ) { 01343 //int cellok = false; 01344 pos -= (this->mvGeoStart); 01345 01346 LbmFloat mmaxsize = mLevel[mMaxRefine].nodeSize; 01347 for(int level=mMaxRefine; level>=0; level--) { // finest first 01348 //for(int level=0; level<=mMaxRefine; level++) { // coarsest first 01349 LbmFloat nsize = mLevel[level].nodeSize; 01350 int x,y,z; 01351 // CHECK +- maxsize? 01352 x = (int)((pos[0]+0.5*mmaxsize) / nsize ); 01353 y = (int)((pos[1]+0.5*mmaxsize) / nsize ); 01354 z = (int)((pos[2]+0.5*mmaxsize) / nsize ); 01355 if(LBMDIM==2) z = 0; 01356 01357 // double check... 01358 if(x<0) continue; 01359 if(y<0) continue; 01360 if(z<0) continue; 01361 if(x>=mLevel[level].lSizex) continue; 01362 if(y>=mLevel[level].lSizey) continue; 01363 if(z>=mLevel[level].lSizez) continue; 01364 01365 // return fluid/if/border cells 01366 if( ( (RFLAG(level, x,y,z, mLevel[level].setCurr)&(CFUnused)) ) || 01367 ( (level<mMaxRefine) && (RFLAG(level, x,y,z, mLevel[level].setCurr)&(CFUnused|CFEmpty)) ) ) { 01368 continue; 01369 } // */ 01370 01371 stdCellId *newcid = new stdCellId; 01372 newcid->level = level; 01373 newcid->x = x; 01374 newcid->y = y; 01375 newcid->z = z; 01376 //errMsg("cellAt",this->mName<<" "<<pos<<" l"<<level<<":"<<x<<","<<y<<","<<z<<" "<<convertCellFlagType2String(RFLAG(level, x,y,z, mLevel[level].setCurr)) ); 01377 return newcid; 01378 } 01379 01380 return NULL; 01381 } 01382 01383 01384 // INFO functions 01385 01386 int LbmFsgrSolver::getCellSet ( CellIdentifierInterface* basecid) { 01387 stdCellId *cid = convertBaseCidToStdCid(basecid); 01388 return mLevel[cid->level].setCurr; 01389 //return mLevel[cid->level].setOther; 01390 } 01391 01392 int LbmFsgrSolver::getCellLevel ( CellIdentifierInterface* basecid) { 01393 stdCellId *cid = convertBaseCidToStdCid(basecid); 01394 return cid->level; 01395 } 01396 01397 ntlVec3Gfx LbmFsgrSolver::getCellOrigin ( CellIdentifierInterface* basecid) { 01398 ntlVec3Gfx ret; 01399 01400 stdCellId *cid = convertBaseCidToStdCid(basecid); 01401 ntlVec3Gfx cs( mLevel[cid->level].nodeSize ); 01402 if(LBMDIM==2) { cs[2] = 0.0; } 01403 01404 if(LBMDIM==2) { 01405 ret =(this->mvGeoStart + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], (this->mvGeoEnd[2]-this->mvGeoStart[2])*0.5 ) 01406 + ntlVec3Gfx(0.0,0.0,cs[1]*-0.25)*cid->level ) 01407 +getCellSize(basecid); 01408 } else { 01409 ret =(this->mvGeoStart + ntlVec3Gfx( cid->x *cs[0], cid->y *cs[1], cid->z *cs[2] )) 01410 +getCellSize(basecid); 01411 } 01412 return (ret); 01413 } 01414 01415 ntlVec3Gfx LbmFsgrSolver::getCellSize ( CellIdentifierInterface* basecid) { 01416 // return half size 01417 stdCellId *cid = convertBaseCidToStdCid(basecid); 01418 ntlVec3Gfx retvec( mLevel[cid->level].nodeSize * 0.5 ); 01419 // 2d display as rectangles 01420 if(LBMDIM==2) { retvec[2] = 0.0; } 01421 return (retvec); 01422 } 01423 01424 LbmFloat LbmFsgrSolver::getCellDensity ( CellIdentifierInterface* basecid,int set) { 01425 stdCellId *cid = convertBaseCidToStdCid(basecid); 01426 01427 // skip non-fluid cells 01428 if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&(CFFluid|CFInter)) { 01429 // ok go on... 01430 } else { 01431 return 0.; 01432 } 01433 01434 LbmFloat rho = 0.0; 01435 FORDF0 { rho += QCELL(cid->level, cid->x,cid->y,cid->z, set, l); } // ORG 01436 return ((rho-1.0) * mLevel[cid->level].simCellSize / mLevel[cid->level].timestep) +1.0; // ORG 01437 /*if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFInter) { // test 01438 LbmFloat ux,uy,uz; 01439 ux=uy=uz= 0.0; 01440 int lev = cid->level; 01441 LbmFloat df[27], feqOld[27]; 01442 FORDF0 { 01443 rho += QCELL(lev, cid->x,cid->y,cid->z, set, l); 01444 ux += this->dfDvecX[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l); 01445 uy += this->dfDvecY[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l); 01446 uz += this->dfDvecZ[l]* QCELL(lev, cid->x,cid->y,cid->z, set, l); 01447 df[l] = QCELL(lev, cid->x,cid->y,cid->z, set, l); 01448 } 01449 FORDF0 { 01450 feqOld[l] = getCollideEq(l, rho,ux,uy,uz); 01451 } 01452 // debugging mods 01453 //const LbmFloat Qo = this->getLesNoneqTensorCoeff(df,feqOld); 01454 //const LbmFloat modOmega = this->getLesOmega(mLevel[lev].omega, mLevel[lev].lcsmago,Qo); 01455 //rho = (2.0-modOmega) *25.0; 01456 //rho = Qo*100.0; 01457 //if(cid->x==24){ errMsg("MODOMT"," at "<<PRINT_VEC(cid->x,cid->y,cid->z)<<" = "<<rho<<" "<<Qo); } 01458 //else{ rho=0.0; } 01459 } // test 01460 return rho; // test */ 01461 } 01462 01463 LbmVec LbmFsgrSolver::getCellVelocity ( CellIdentifierInterface* basecid,int set) { 01464 stdCellId *cid = convertBaseCidToStdCid(basecid); 01465 01466 // skip non-fluid cells 01467 if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&(CFFluid|CFInter)) { 01468 // ok go on... 01469 } else { 01470 return LbmVec(0.0); 01471 } 01472 01473 LbmFloat ux,uy,uz; 01474 ux=uy=uz= 0.0; 01475 FORDF0 { 01476 ux += this->dfDvecX[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l); 01477 uy += this->dfDvecY[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l); 01478 uz += this->dfDvecZ[l]* QCELL(cid->level, cid->x,cid->y,cid->z, set, l); 01479 } 01480 LbmVec vel(ux,uy,uz); 01481 // TODO fix... 01482 return (vel * mLevel[cid->level].simCellSize / mLevel[cid->level].timestep * this->mDebugVelScale); // normal 01483 } 01484 01485 LbmFloat LbmFsgrSolver::getCellDf( CellIdentifierInterface* basecid,int set, int dir) { 01486 stdCellId *cid = convertBaseCidToStdCid(basecid); 01487 return QCELL(cid->level, cid->x,cid->y,cid->z, set, dir); 01488 } 01489 LbmFloat LbmFsgrSolver::getCellMass( CellIdentifierInterface* basecid,int set) { 01490 stdCellId *cid = convertBaseCidToStdCid(basecid); 01491 return QCELL(cid->level, cid->x,cid->y,cid->z, set, dMass); 01492 } 01493 LbmFloat LbmFsgrSolver::getCellFill( CellIdentifierInterface* basecid,int set) { 01494 stdCellId *cid = convertBaseCidToStdCid(basecid); 01495 if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFInter) return QCELL(cid->level, cid->x,cid->y,cid->z, set, dFfrac); 01496 if(RFLAG(cid->level, cid->x,cid->y,cid->z, set)&CFFluid) return 1.0; 01497 return 0.0; 01498 //return QCELL(cid->level, cid->x,cid->y,cid->z, set, dFfrac); 01499 } 01500 CellFlagType LbmFsgrSolver::getCellFlag( CellIdentifierInterface* basecid,int set) { 01501 stdCellId *cid = convertBaseCidToStdCid(basecid); 01502 return RFLAG(cid->level, cid->x,cid->y,cid->z, set); 01503 } 01504 01505 LbmFloat LbmFsgrSolver::getEquilDf( int l ) { 01506 return this->dfEquil[l]; 01507 } 01508 01509 01510 ntlVec3Gfx LbmFsgrSolver::getVelocityAt (float xp, float yp, float zp) { 01511 ntlVec3Gfx avgvel(0.0); 01512 LbmFloat avgnum = 0.; 01513 01514 // taken from getCellAt! 01515 const int level = mMaxRefine; 01516 const int workSet = mLevel[level].setCurr; 01517 const LbmFloat nsize = mLevel[level].nodeSize; 01518 const int x = (int)((-this->mvGeoStart[0]+xp-0.5*nsize) / nsize ); 01519 const int y = (int)((-this->mvGeoStart[1]+yp-0.5*nsize) / nsize ); 01520 int z = (int)((-this->mvGeoStart[2]+zp-0.5*nsize) / nsize ); 01521 if(LBMDIM==2) z=0; 01522 //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) ); 01523 01524 // return fluid/if/border cells 01525 // search neighborhood, do smoothing 01526 FORDF0{ 01527 const int i = x+this->dfVecX[l]; 01528 const int j = y+this->dfVecY[l]; 01529 const int k = z+this->dfVecZ[l]; 01530 01531 if( (i<0) || (j<0) || (k<0) 01532 || (i>=mLevel[level].lSizex) 01533 || (j>=mLevel[level].lSizey) 01534 || (k>=mLevel[level].lSizez) ) continue; 01535 01536 if( (RFLAG(level, i,j,k, mLevel[level].setCurr)&(CFFluid|CFInter)) ) { 01537 ntlVec3Gfx vel(0.0); 01538 LbmFloat *ccel = RACPNT(level, i,j,k ,workSet); // omp 01539 for(int n=1; n<this->cDfNum; n++) { 01540 vel[0] += (this->dfDvecX[n]*RAC(ccel,n)); 01541 vel[1] += (this->dfDvecY[n]*RAC(ccel,n)); 01542 vel[2] += (this->dfDvecZ[n]*RAC(ccel,n)); 01543 } 01544 01545 avgvel += vel; 01546 avgnum += 1.0; 01547 if(l==0) { // center slightly more weight 01548 avgvel += vel; avgnum += 1.0; 01549 } 01550 } // */ 01551 } 01552 01553 if(avgnum>0.) { 01554 ntlVec3Gfx retv = avgvel / avgnum; 01555 retv *= nsize/mLevel[level].timestep; 01556 // scale for current animation settings (frame time) 01557 retv *= mpParam->getCurrentAniFrameTime(); 01558 //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() ); 01559 return retv; 01560 } 01561 // no cells here...? 01562 //errMsg("DUMPVEL"," at "<<PRINT_VEC(xp,yp,zp)<<" v"<<avgvel<<" n"<<avgnum<<" no vel !?"); 01563 return ntlVec3Gfx(0.); 01564 } 01565 01566 #if LBM_USE_GUI==1 01567 01568 void 01569 LbmFsgrSolver::debugDisplay(int set){ 01570 //lbmDebugDisplay< LbmFsgrSolver >( set, this ); 01571 lbmDebugDisplay( set ); 01572 } 01573 #endif 01574 01575 /*****************************************************************************/ 01576 // strict debugging functions 01577 /*****************************************************************************/ 01578 #if FSGR_STRICT_DEBUG==1 01579 #define STRICT_EXIT *((int *)0)=0; 01580 01581 int LbmFsgrSolver::debLBMGI(int level, int ii,int ij,int ik, int is) { 01582 if(level < 0){ errMsg("LbmStrict::debLBMGI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01583 if(level > mMaxRefine){ errMsg("LbmStrict::debLBMGI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01584 01585 if((ii==-1)&&(ij==0)) { 01586 // special case for main loop, ok 01587 } else { 01588 if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01589 if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01590 if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01591 if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01592 } 01593 if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01594 if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01595 if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01596 if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01597 return _LBMGI(level, ii,ij,ik, is); 01598 }; 01599 01600 CellFlagType& LbmFsgrSolver::debRFLAG(int level, int xx,int yy,int zz,int set){ 01601 return _RFLAG(level, xx,yy,zz,set); 01602 }; 01603 01604 CellFlagType& LbmFsgrSolver::debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir) { 01605 if(dir<0) { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; } 01606 // warning might access all spatial nbs 01607 if(dir>this->cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; } 01608 return _RFLAG_NB(level, xx,yy,zz,set, dir); 01609 }; 01610 01611 CellFlagType& LbmFsgrSolver::debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir) { 01612 if(dir<0) { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; } 01613 if(dir>this->cDirNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; } 01614 return _RFLAG_NBINV(level, xx,yy,zz,set, dir); 01615 }; 01616 01617 int LbmFsgrSolver::debLBMQI(int level, int ii,int ij,int ik, int is, int l) { 01618 if(level < 0){ errMsg("LbmStrict::debLBMQI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01619 if(level > mMaxRefine){ errMsg("LbmStrict::debLBMQI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01620 01621 if((ii==-1)&&(ij==0)) { 01622 // special case for main loop, ok 01623 } else { 01624 if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01625 if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01626 if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01627 if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01628 } 01629 if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01630 if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01631 if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01632 if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; } 01633 if(l<0) { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; } 01634 if(l>this->cDfNum){ // dFfrac is an exception 01635 if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } } 01636 #if COMPRESSGRIDS==1 01637 //if((!this->mInitDone) && (is!=mLevel[level].setCurr)){ STRICT_EXIT; } // COMPRT debug 01638 #endif // COMPRESSGRIDS==1 01639 return _LBMQI(level, ii,ij,ik, is, l); 01640 }; 01641 01642 LbmFloat& LbmFsgrSolver::debQCELL(int level, int xx,int yy,int zz,int set,int l) { 01643 //errMsg("LbmStrict","debQCELL debug: l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" l"<<l<<" index"<<LBMGI(level, xx,yy,zz,set)); 01644 return _QCELL(level, xx,yy,zz,set,l); 01645 }; 01646 01647 LbmFloat& LbmFsgrSolver::debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l) { 01648 if(dir<0) { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; } 01649 if(dir>this->cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; } 01650 return _QCELL_NB(level, xx,yy,zz,set, dir,l); 01651 }; 01652 01653 LbmFloat& LbmFsgrSolver::debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l) { 01654 if(dir<0) { errMsg("LbmStrict"," invD- l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; } 01655 if(dir>this->cDfNum){ errMsg("LbmStrict"," invD+ l"<<level<<"|"<<xx<<","<<yy<<","<<zz<<" s"<<set<<" d"<<dir); STRICT_EXIT; } 01656 return _QCELL_NBINV(level, xx,yy,zz,set, dir,l); 01657 }; 01658 01659 LbmFloat* LbmFsgrSolver::debRACPNT(int level, int ii,int ij,int ik, int is ) { 01660 return _RACPNT(level, ii,ij,ik, is ); 01661 }; 01662 01663 LbmFloat& LbmFsgrSolver::debRAC(LbmFloat* s,int l) { 01664 if(l<0) { errMsg("LbmStrict"," invD- "<<" l"<<l); STRICT_EXIT; } 01665 if(l>dTotalNum){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } 01666 //if(l>this->cDfNum){ // dFfrac is an exception 01667 //if((l != dMass) && (l != dFfrac) && (l != dFlux)){ errMsg("LbmStrict"," invD+ "<<" l"<<l); STRICT_EXIT; } } 01668 return _RAC(s,l); 01669 }; 01670 01671 #endif // FSGR_STRICT_DEBUG==1 01672 01673 01674 /****************************************************************************** 01675 * GUI&debugging functions 01676 *****************************************************************************/ 01677 01678 01679 #if LBM_USE_GUI==1 01680 #define USE_GLUTILITIES 01681 #include "../gui/gui_utilities.h" 01682 01684 void LbmFsgrSolver::debugDisplayNode(int dispset, CellIdentifierInterface* cell ) { 01685 //debugOut(" DD: "<<cell->getAsString() , 10); 01686 ntlVec3Gfx org = this->getCellOrigin( cell ); 01687 ntlVec3Gfx halfsize = this->getCellSize( cell ); 01688 int set = this->getCellSet( cell ); 01689 //debugOut(" DD: "<<cell->getAsString()<<" "<< (dispset->type) , 10); 01690 01691 bool showcell = true; 01692 int linewidth = 1; 01693 ntlColor col(0.5); 01694 LbmFloat cscale = 1.0; //dispset->scale; 01695 01696 #define DRAWDISPCUBE(col,scale) \ 01697 { glLineWidth( linewidth ); \ 01698 glColor3f( (col)[0], (col)[1], (col)[2]); \ 01699 ntlVec3Gfx s = org-(halfsize * (scale)); \ 01700 ntlVec3Gfx e = org+(halfsize * (scale)); \ 01701 drawCubeWire( s,e ); } 01702 01703 CellFlagType flag = this->getCellFlag(cell, set ); 01704 // always check types 01705 if(flag& CFInvalid ) { if(!guiShowInvalid ) return; } 01706 if(flag& CFUnused ) { if(!guiShowInvalid ) return; } 01707 if(flag& CFEmpty ) { if(!guiShowEmpty ) return; } 01708 if(flag& CFInter ) { if(!guiShowInterface) return; } 01709 if(flag& CFNoDelete ) { if(!guiShowNoDelete ) return; } 01710 if(flag& CFBnd ) { if(!guiShowBnd ) return; } 01711 01712 // only dismiss one of these types 01713 if(flag& CFGrFromCoarse) { if(!guiShowCoarseInner ) return; } // inner not really interesting 01714 else 01715 if(flag& CFGrFromFine) { if(!guiShowCoarseBorder ) return; } 01716 else 01717 if(flag& CFFluid ) { if(!guiShowFluid ) return; } 01718 01719 switch(dispset) { 01720 case FLUIDDISPNothing: { 01721 showcell = false; 01722 } break; 01723 case FLUIDDISPCelltypes: { 01724 cscale = 0.5; 01725 01726 if(flag& CFNoDelete) { // debug, mark nodel cells 01727 ntlColor ccol(0.7,0.0,0.0); 01728 DRAWDISPCUBE(ccol, 0.1); 01729 } 01730 if(flag& CFPersistMask) { // mark persistent flags 01731 ntlColor ccol(0.5); 01732 DRAWDISPCUBE(ccol, 0.125); 01733 } 01734 if(flag& CFNoBndFluid) { // mark persistent flags 01735 ntlColor ccol(0,0,1); 01736 DRAWDISPCUBE(ccol, 0.075); 01737 } 01738 01739 if(flag& CFInvalid) { 01740 cscale = 0.50; 01741 col = ntlColor(0.0,0,0.0); 01742 } 01743 else if(flag& CFBnd) { 01744 cscale = 0.59; 01745 col = ntlColor(0.4); 01746 } 01747 01748 else if(flag& CFInter) { 01749 cscale = 0.55; 01750 col = ntlColor(0,1,1); 01751 01752 } else if(flag& CFGrFromCoarse) { 01753 // draw as - with marker 01754 ntlColor col2(0.0,1.0,0.3); 01755 DRAWDISPCUBE(col2, 0.1); 01756 cscale = 0.5; 01757 showcell=false; // DEBUG 01758 } 01759 else if(flag& CFFluid) { 01760 cscale = 0.5; 01761 if(flag& CFGrToFine) { 01762 ntlColor col2(0.5,0.0,0.5); 01763 DRAWDISPCUBE(col2, 0.1); 01764 col = ntlColor(0,0,1); 01765 } 01766 if(flag& CFGrFromFine) { 01767 ntlColor col2(1.0,1.0,0.0); 01768 DRAWDISPCUBE(col2, 0.1); 01769 col = ntlColor(0,0,1); 01770 } else if(flag& CFGrFromCoarse) { 01771 // draw as fluid with marker 01772 ntlColor col2(0.0,1.0,0.3); 01773 DRAWDISPCUBE(col2, 0.1); 01774 col = ntlColor(0,0,1); 01775 } else { 01776 col = ntlColor(0,0,1); 01777 } 01778 } 01779 else if(flag& CFEmpty) { 01780 showcell=false; 01781 } 01782 01783 } break; 01784 case FLUIDDISPVelocities: { 01785 // dont use cube display 01786 LbmVec vel = this->getCellVelocity( cell, set ); 01787 glBegin(GL_LINES); 01788 glColor3f( 0.0,0.0,0.0 ); 01789 glVertex3f( org[0], org[1], org[2] ); 01790 org += vec2G(vel * 10.0 * cscale); 01791 glColor3f( 1.0,1.0,1.0 ); 01792 glVertex3f( org[0], org[1], org[2] ); 01793 glEnd(); 01794 showcell = false; 01795 } break; 01796 case FLUIDDISPCellfills: { 01797 cscale = 0.5; 01798 if(flag& CFFluid) { 01799 cscale = 0.75; 01800 col = ntlColor(0,0,0.5); 01801 } 01802 else if(flag& CFInter) { 01803 cscale = 0.75 * this->getCellMass(cell,set); 01804 col = ntlColor(0,1,1); 01805 } 01806 else { 01807 showcell=false; 01808 } 01809 01810 if( ABS(this->getCellMass(cell,set)) < 10.0 ) { 01811 cscale = 0.75 * this->getCellMass(cell,set); 01812 } else { 01813 showcell = false; 01814 } 01815 if(cscale>0.0) { 01816 col = ntlColor(0,1,1); 01817 } else { 01818 col = ntlColor(1,1,0); 01819 } 01820 // TODO 01821 } break; 01822 case FLUIDDISPDensity: { 01823 LbmFloat rho = this->getCellDensity(cell,set); 01824 cscale = rho*rho * 0.25; 01825 col = ntlColor( MIN(0.5+cscale,1.0) , MIN(0.0+cscale,1.0), MIN(0.0+cscale,1.0) ); 01826 cscale *= 2.0; 01827 } break; 01828 case FLUIDDISPGrid: { 01829 cscale = 0.59; 01830 col = ntlColor(1.0); 01831 } break; 01832 default: { 01833 cscale = 0.5; 01834 col = ntlColor(1.0,0.0,0.0); 01835 } break; 01836 } 01837 01838 if(!showcell) return; 01839 if(cscale==0.0) return; // dont draw zero values 01840 DRAWDISPCUBE(col, cscale); 01841 } 01842 01844 // D has to implement the CellIterator interface 01845 void LbmFsgrSolver::lbmDebugDisplay(int dispset) { 01846 // DEBUG always display testdata 01847 #if LBM_INCLUDE_TESTSOLVERS==1 01848 if(mUseTestdata){ 01849 cpDebugDisplay(dispset); 01850 mpTest->testDebugDisplay(dispset); 01851 } 01852 #endif // LBM_INCLUDE_TESTSOLVERS==1 01853 if(dispset<=FLUIDDISPNothing) return; 01854 //if(!dispset->on) return; 01855 glDisable( GL_LIGHTING ); // dont light lines 01856 01857 #if LBM_INCLUDE_TESTSOLVERS==1 01858 if((!mUseTestdata)|| (mUseTestdata)&&(mpTest->mFarfMode<=0)) { 01859 #endif // LBM_INCLUDE_TESTSOLVERS==1 01860 01861 LbmFsgrSolver::CellIdentifier cid = this->getFirstCell(); 01862 for(; this->noEndCell( cid ); 01863 this->advanceCell( cid ) ) { 01864 this->debugDisplayNode(dispset, cid ); 01865 } 01866 delete cid; 01867 01868 #if LBM_INCLUDE_TESTSOLVERS==1 01869 } // 3d check 01870 #endif // LBM_INCLUDE_TESTSOLVERS==1 01871 01872 glEnable( GL_LIGHTING ); // dont light lines 01873 } 01874 01876 // D has to implement the CellIterator interface 01877 void LbmFsgrSolver::lbmMarkedCellDisplay() { 01878 //fluidDispSettings dispset; 01879 // trick - display marked cells as grid displa -> white, big 01880 int dispset = FLUIDDISPGrid; 01881 glDisable( GL_LIGHTING ); // dont light lines 01882 01883 LbmFsgrSolver::CellIdentifier cid = this->markedGetFirstCell(); 01884 while(cid) { 01885 this->debugDisplayNode(dispset, cid ); 01886 cid = this->markedAdvanceCell(); 01887 } 01888 delete cid; 01889 01890 glEnable( GL_LIGHTING ); // dont light lines 01891 } 01892 01893 #endif // LBM_USE_GUI==1 01894 01896 void LbmFsgrSolver::debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet) { 01897 //string printInfo, 01898 // force printing of one set? default = -1 = off 01899 bool printDF = false; 01900 bool printRho = false; 01901 bool printVel = false; 01902 bool printFlag = false; 01903 bool printGeom = false; 01904 bool printMass=false; 01905 bool printBothSets = false; 01906 string printInfo = this->getNodeInfoString(); 01907 01908 for(size_t i=0; i<printInfo.length()-0; i++) { 01909 char what = printInfo[i]; 01910 switch(what) { 01911 case '+': // all on 01912 printDF = true; printRho = true; printVel = true; printFlag = true; printGeom = true; printMass = true ; 01913 printBothSets = true; break; 01914 case '-': // all off 01915 printDF = false; printRho = false; printVel = false; printFlag = false; printGeom = false; printMass = false; 01916 printBothSets = false; break; 01917 case 'd': printDF = true; break; 01918 case 'r': printRho = true; break; 01919 case 'v': printVel = true; break; 01920 case 'f': printFlag = true; break; 01921 case 'g': printGeom = true; break; 01922 case 'm': printMass = true; break; 01923 case 's': printBothSets = true; break; 01924 default: 01925 errFatal("debugPrintNodeInfo","Invalid node info id "<<what,SIMWORLD_GENERICERROR); return; 01926 } 01927 } 01928 01929 ntlVec3Gfx org = this->getCellOrigin( cell ); 01930 ntlVec3Gfx halfsize = this->getCellSize( cell ); 01931 int set = this->getCellSet( cell ); 01932 debMsgStd("debugPrintNodeInfo",DM_NOTIFY, "Printing cell info '"<<printInfo<<"' for node: "<<cell->getAsString()<<" from "<<this->getName()<<" currSet:"<<set , 1); 01933 if(printGeom) debMsgStd(" ",DM_MSG, "Org:"<<org<<" Halfsize:"<<halfsize<<" ", 1); 01934 01935 int setmax = 2; 01936 if(!printBothSets) setmax = 1; 01937 if(forceSet>=0) setmax = 1; 01938 01939 for(int s=0; s<setmax; s++) { 01940 int workset = set; 01941 if(s==1){ workset = (set^1); } 01942 if(forceSet>=0) workset = forceSet; 01943 debMsgStd(" ",DM_MSG, "Printing set:"<<workset<<" orgSet:"<<set, 1); 01944 01945 if(printDF) { 01946 for(int l=0; l<LBM_DFNUM; l++) { // FIXME ?? 01947 debMsgStd(" ",DM_MSG, " Df"<<l<<": "<<this->getCellDf(cell,workset,l), 1); 01948 } 01949 } 01950 if(printRho) { 01951 debMsgStd(" ",DM_MSG, " Rho: "<<this->getCellDensity(cell,workset), 1); 01952 } 01953 if(printVel) { 01954 debMsgStd(" ",DM_MSG, " Vel: "<<this->getCellVelocity(cell,workset), 1); 01955 } 01956 if(printFlag) { 01957 CellFlagType flag = this->getCellFlag(cell,workset); 01958 debMsgStd(" ",DM_MSG, " Flg: "<< flag<<" "<<convertFlags2String( flag ) <<" "<<convertCellFlagType2String( flag ), 1); 01959 } 01960 if(printMass) { 01961 debMsgStd(" ",DM_MSG, " Mss: "<<this->getCellMass(cell,workset), 1); 01962 } 01963 // dirty... TODO fixme 01964 debMsgStd(" ",DM_MSG, " Flx: "<<this->getCellDf(cell,workset,dFlux), 1); 01965 } 01966 } 01967 01968