|
Blender
V2.59
|
00001 00004 /****************************************************************************** 00005 * 00006 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method 00007 * Copyright 2003-2006 Nils Thuerey 00008 * 00009 * Particle Viewer/Tracer 00010 * 00011 *****************************************************************************/ 00012 00013 #include <stdio.h> 00014 //#include "../libs/my_gl.h" 00015 //#include "../libs/my_glu.h" 00016 00017 /* own lib's */ 00018 #include "particletracer.h" 00019 #include "ntl_matrices.h" 00020 #include "ntl_ray.h" 00021 #include "ntl_matrices.h" 00022 00023 #include <zlib.h> 00024 00025 00026 // particle object id counter 00027 int ParticleObjectIdCnt = 1; 00028 00029 /****************************************************************************** 00030 * Standard constructor 00031 *****************************************************************************/ 00032 ParticleTracer::ParticleTracer() : 00033 ntlGeometryObject(), 00034 mParts(), 00035 //mTrailLength(1), mTrailInterval(1),mTrailIntervalCounter(0), 00036 mPartSize(0.01), 00037 mStart(-1.0), mEnd(1.0), 00038 mSimStart(-1.0), mSimEnd(1.0), 00039 mPartScale(0.1) , mPartHeadDist( 0.1 ), mPartTailDist( -0.1 ), mPartSegments( 4 ), 00040 mValueScale(0), 00041 mValueCutoffTop(0.0), mValueCutoffBottom(0.0), 00042 mDumpParts(0), //mDumpText(0), 00043 mDumpTextFile(""), 00044 mDumpTextInterval(0.), mDumpTextLastTime(0.), mDumpTextCount(0), 00045 mShowOnly(0), 00046 mNumInitialParts(0), mpTrafo(NULL), 00047 mInitStart(-1.), mInitEnd(-1.), 00048 mPrevs(), mTrailTimeLast(0.), mTrailInterval(-1.), mTrailLength(0) 00049 { 00050 debMsgStd("ParticleTracer::ParticleTracer",DM_MSG,"inited",10); 00051 }; 00052 00053 ParticleTracer::~ParticleTracer() { 00054 debMsgStd("ParticleTracer::~ParticleTracer",DM_MSG,"destroyed",10); 00055 } 00056 00057 /*****************************************************************************/ 00059 /*****************************************************************************/ 00060 void ParticleTracer::parseAttrList(AttributeList *att) 00061 { 00062 AttributeList *tempAtt = mpAttrs; 00063 mpAttrs = att; 00064 00065 mNumInitialParts = mpAttrs->readInt("particles",mNumInitialParts, "ParticleTracer","mNumInitialParts", false); 00066 //errMsg(" NUMP"," "<<mNumInitialParts); 00067 mPartScale = mpAttrs->readFloat("part_scale",mPartScale, "ParticleTracer","mPartScale", false); 00068 mPartHeadDist = mpAttrs->readFloat("part_headdist",mPartHeadDist, "ParticleTracer","mPartHeadDist", false); 00069 mPartTailDist = mpAttrs->readFloat("part_taildist",mPartTailDist, "ParticleTracer","mPartTailDist", false); 00070 mPartSegments = mpAttrs->readInt ("part_segments",mPartSegments, "ParticleTracer","mPartSegments", false); 00071 mValueScale = mpAttrs->readInt ("part_valscale",mValueScale, "ParticleTracer","mValueScale", false); 00072 mValueCutoffTop = mpAttrs->readFloat("part_valcutofftop",mValueCutoffTop, "ParticleTracer","mValueCutoffTop", false); 00073 mValueCutoffBottom = mpAttrs->readFloat("part_valcutoffbottom",mValueCutoffBottom, "ParticleTracer","mValueCutoffBottom", false); 00074 00075 mDumpParts = mpAttrs->readInt ("part_dump",mDumpParts, "ParticleTracer","mDumpParts", false); 00076 // mDumpText deprecatd, use mDumpTextInterval>0. instead 00077 mShowOnly = mpAttrs->readInt ("part_showonly",mShowOnly, "ParticleTracer","mShowOnly", false); 00078 mDumpTextFile= mpAttrs->readString("part_textdumpfile",mDumpTextFile, "ParticleTracer","mDumpTextFile", false); 00079 mDumpTextInterval= mpAttrs->readFloat("part_textdumpinterval",mDumpTextInterval, "ParticleTracer","mDumpTextInterval", false); 00080 00081 string matPart; 00082 matPart = mpAttrs->readString("material_part", "default", "ParticleTracer","material", false); 00083 setMaterialName( matPart ); 00084 00085 mInitStart = mpAttrs->readFloat("part_initstart",mInitStart, "ParticleTracer","mInitStart", false); 00086 mInitEnd = mpAttrs->readFloat("part_initend", mInitEnd, "ParticleTracer","mInitEnd", false); 00087 00088 // unused... 00089 //int mTrailLength = 0; // UNUSED 00090 //int mTrailInterval= 0; // UNUSED 00091 mTrailLength = mpAttrs->readInt("traillength",mTrailLength, "ParticleTracer","mTrailLength", false); 00092 mTrailInterval= mpAttrs->readFloat("trailinterval",mTrailInterval, "ParticleTracer","mTrailInterval", false); 00093 00094 // restore old list 00095 mpAttrs = tempAtt; 00096 } 00097 00098 /****************************************************************************** 00099 * draw the particle array 00100 *****************************************************************************/ 00101 void ParticleTracer::draw() 00102 { 00103 } 00104 00105 /****************************************************************************** 00106 * init trafo matrix 00107 *****************************************************************************/ 00108 void ParticleTracer::initTrafoMatrix() { 00109 ntlVec3Gfx scale = ntlVec3Gfx( 00110 (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]), 00111 (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]), 00112 (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2]) 00113 ); 00114 ntlVec3Gfx trans = mStart; 00115 if(!mpTrafo) mpTrafo = new ntlMat4Gfx(0.0); 00116 mpTrafo->initId(); 00117 for(int i=0; i<3; i++) { mpTrafo->value[i][i] = scale[i]; } 00118 for(int i=0; i<3; i++) { mpTrafo->value[i][3] = trans[i]; } 00119 } 00120 00121 /****************************************************************************** 00122 * adapt time step by rescaling velocities 00123 *****************************************************************************/ 00124 void ParticleTracer::adaptPartTimestep(float factor) { 00125 for(size_t i=0; i<mParts.size(); i++) { 00126 mParts[i].setVel( mParts[i].getVel() * factor ); 00127 } 00128 } 00129 00130 00131 /****************************************************************************** 00132 * add a particle at this position 00133 *****************************************************************************/ 00134 void ParticleTracer::addParticle(float x, float y, float z) { 00135 ntlVec3Gfx p(x,y,z); 00136 ParticleObject part( p ); 00137 mParts.push_back( part ); 00138 } 00139 void ParticleTracer::addFullParticle(ParticleObject &np) { 00140 mParts.push_back( np ); 00141 } 00142 00143 00144 void ParticleTracer::cleanup() { 00145 // cleanup 00146 int last = (int)mParts.size()-1; 00147 if(mDumpTextInterval>0.) { errMsg("ParticleTracer::cleanup","Skipping cleanup due to text dump..."); return; } 00148 00149 for(int i=0; i<=last; i++) { 00150 if( mParts[i].getActive()==false ) { 00151 ParticleObject *p = &mParts[i]; 00152 ParticleObject *p2 = &mParts[last]; 00153 *p = *p2; last--; mParts.pop_back(); 00154 } 00155 } 00156 } 00157 00158 extern bool glob_mpactive; 00159 extern int glob_mpindex,glob_mpnum; 00160 00161 /****************************************************************************** 00162 *! dump particles if desired 00163 *****************************************************************************/ 00164 void ParticleTracer::notifyOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename, double simtime) { 00165 debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"obj:"<<this->getName()<<" frame:"<<frameNrStr<<" dumpp"<<mDumpParts<<" t"<<simtime, 10); // DEBUG 00166 00167 if( 00168 (dumptype==DUMP_FULLGEOMETRY)&& 00169 (mDumpParts>0)) { 00170 // dump to binary file 00171 std::ostringstream boutfilename(""); 00172 boutfilename << outfilename <<"_particles_" << frameNrStr; 00173 if(glob_mpactive) { 00174 if(glob_mpindex>0) { boutfilename << "mp"<<glob_mpindex; } 00175 } 00176 boutfilename << ".gz"; 00177 debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" #"<<frameNr , 7); 00178 //debMsgStd("ParticleTracer::notifyOfDump",DM_MSG,"B-Dumping: partgeodeb sim:"<<mSimStart<<","<<mSimEnd<<" geosize:"<<mStart<<","<<mEnd,2 ); 00179 00180 // output to zipped file 00181 gzFile gzf; 00182 gzf = gzopen(boutfilename.str().c_str(), "wb1"); 00183 if(gzf) { 00184 int numParts; 00185 if(sizeof(numParts)!=4) { errMsg("ParticleTracer::notifyOfDump","Invalid int size"); return; } 00186 // only dump active particles 00187 numParts = 0; 00188 for(size_t i=0; i<mParts.size(); i++) { 00189 if(!mParts[i].getActive()) continue; 00190 numParts++; 00191 } 00192 gzwrite(gzf, &numParts, sizeof(numParts)); 00193 for(size_t i=0; i<mParts.size(); i++) { 00194 if(!mParts[i].getActive()) { continue; } 00195 ParticleObject *p = &mParts[i]; 00196 //int type = p->getType(); // export whole type info 00197 int type = p->getFlags(); // debug export whole type & status info 00198 ntlVec3Gfx pos = p->getPos(); 00199 float size = p->getSize(); 00200 00201 if(type&PART_FLOAT) { // WARNING same handling for dump! 00202 // add one gridcell offset 00203 //pos[2] += 1.0; 00204 } 00205 // display as drop for now externally 00206 //else if(type&PART_TRACER) { type |= PART_DROP; } 00207 00208 pos = (*mpTrafo) * pos; 00209 00210 ntlVec3Gfx v = p->getVel(); 00211 v[0] *= mpTrafo->value[0][0]; 00212 v[1] *= mpTrafo->value[1][1]; 00213 v[2] *= mpTrafo->value[2][2]; 00214 // FIXME check: pos = (*mpTrafo) * pos; 00215 gzwrite(gzf, &type, sizeof(type)); 00216 gzwrite(gzf, &size, sizeof(size)); 00217 for(int j=0; j<3; j++) { gzwrite(gzf, &pos[j], sizeof(float)); } 00218 for(int j=0; j<3; j++) { gzwrite(gzf, &v[j], sizeof(float)); } 00219 } 00220 gzclose( gzf ); 00221 } 00222 } // dump? 00223 } 00224 00225 void ParticleTracer::checkDumpTextPositions(double simtime) { 00226 // dfor partial & full dump 00227 if(mDumpTextInterval>0.) { 00228 debMsgStd("ParticleTracer::checkDumpTextPositions",DM_MSG,"t="<<simtime<<" last:"<<mDumpTextLastTime<<" inter:"<<mDumpTextInterval,7); 00229 } 00230 00231 if((mDumpTextInterval>0.) && (simtime>mDumpTextLastTime+mDumpTextInterval)) { 00232 // dump to binary file 00233 std::ostringstream boutfilename(""); 00234 if(mDumpTextFile.length()>1) { 00235 boutfilename << mDumpTextFile << ".cpart2"; 00236 } else { 00237 boutfilename << boutfilename <<"_particles" << ".cpart2"; 00238 } 00239 debMsgStd("ParticleTracer::checkDumpTextPositions",DM_MSG,"T-Dumping: "<< this->getName() <<", particles:"<<mParts.size()<<" "<< " to "<<boutfilename.str()<<" " , 7); 00240 00241 int numParts = 0; 00242 // only dump bubble particles 00243 for(size_t i=0; i<mParts.size(); i++) { 00244 //if(!mParts[i].getActive()) continue; 00245 //if(!(mParts[i].getType()&PART_BUBBLE)) continue; 00246 numParts++; 00247 } 00248 00249 // output to text file 00250 //gzFile gzf; 00251 FILE *stf; 00252 if(mDumpTextCount==0) { 00253 //gzf = gzopen(boutfilename.str().c_str(), "w0"); 00254 stf = fopen(boutfilename.str().c_str(), "w"); 00255 00256 fprintf( stf, "\n\n# cparts generated by elbeem \n# no. of parts \nN %d \n\n",numParts); 00257 // fixed time scale for now 00258 fprintf( stf, "T %f \n\n", 1.0); 00259 } else { 00260 //gzf = gzopen(boutfilename.str().c_str(), "a+0"); 00261 stf = fopen(boutfilename.str().c_str(), "a+"); 00262 } 00263 00264 // add current set 00265 if(stf) { 00266 fprintf( stf, "\n\n# new set at frame %d,t%f,p%d --------------------------------- \n\n", mDumpTextCount, simtime, numParts ); 00267 fprintf( stf, "S %f \n\n", simtime ); 00268 00269 for(size_t i=0; i<mParts.size(); i++) { 00270 ParticleObject *p = &mParts[i]; 00271 ntlVec3Gfx pos = p->getPos(); 00272 float size = p->getSize(); 00273 float infl = 1.; 00274 //if(!mParts[i].getActive()) { size=0.; } // switch "off" 00275 if(!mParts[i].getActive()) { infl=0.; } // switch "off" 00276 if(!mParts[i].getInFluid()) { infl=0.; } // switch "off" 00277 if(mParts[i].getLifeTime()<0.) { infl=0.; } // not yet active... 00278 00279 pos = (*mpTrafo) * pos; 00280 ntlVec3Gfx v = p->getVel(); 00281 v[0] *= mpTrafo->value[0][0]; 00282 v[1] *= mpTrafo->value[1][1]; 00283 v[2] *= mpTrafo->value[2][2]; 00284 00285 fprintf( stf, "P %f %f %f \n", pos[0],pos[1],pos[2] ); 00286 if(size!=1.0) fprintf( stf, "s %f \n", size ); 00287 if(infl!=1.0) fprintf( stf, "i %f \n", infl ); 00288 fprintf( stf, "\n" ); 00289 } 00290 00291 fprintf( stf, "# %d end ", mDumpTextCount ); 00292 //gzclose( gzf ); 00293 fclose( stf ); 00294 00295 mDumpTextCount++; 00296 } 00297 00298 mDumpTextLastTime += mDumpTextInterval; 00299 } 00300 00301 } 00302 00303 00304 void ParticleTracer::checkTrails(double time) { 00305 if(mTrailLength<1) return; 00306 if(time-mTrailTimeLast > mTrailInterval) { 00307 00308 if( (int)mPrevs.size() < mTrailLength) mPrevs.resize( mTrailLength ); 00309 for(int i=mPrevs.size()-1; i>0; i--) { 00310 mPrevs[i] = mPrevs[i-1]; 00311 //errMsg("TRAIL"," from "<<i<<" to "<<(i-1) ); 00312 } 00313 mPrevs[0] = mParts; 00314 00315 mTrailTimeLast += mTrailInterval; 00316 } 00317 } 00318 00319 00320 /****************************************************************************** 00321 * Get triangles for rendering 00322 *****************************************************************************/ 00323 void ParticleTracer::getTriangles(double time, vector<ntlTriangle> *triangles, 00324 vector<ntlVec3Gfx> *vertices, 00325 vector<ntlVec3Gfx> *normals, int objectId ) 00326 { 00327 #ifdef ELBEEM_PLUGIN 00328 // suppress warnings... 00329 vertices = NULL; triangles = NULL; 00330 normals = NULL; objectId = 0; 00331 time = 0.; 00332 #else // ELBEEM_PLUGIN 00333 int pcnt = 0; 00334 // currently not used in blender 00335 objectId = 0; // remove, deprecated 00336 if(mDumpParts>1) { 00337 return; // only dump, no tri-gen 00338 } 00339 00340 const bool debugParts = false; 00341 int tris = 0; 00342 int segments = mPartSegments; 00343 ntlVec3Gfx scale = ntlVec3Gfx( (mEnd[0]-mStart[0])/(mSimEnd[0]-mSimStart[0]), (mEnd[1]-mStart[1])/(mSimEnd[1]-mSimStart[1]), (mEnd[2]-mStart[2])/(mSimEnd[2]-mSimStart[2])); 00344 ntlVec3Gfx trans = mStart; 00345 time = 0.; // doesnt matter 00346 00347 for(size_t t=0; t<mPrevs.size()+1; t++) { 00348 vector<ParticleObject> *dparts; 00349 if(t==0) { 00350 dparts = &mParts; 00351 } else { 00352 dparts = &mPrevs[t-1]; 00353 } 00354 //errMsg("TRAILT","prevs"<<t<<"/"<<mPrevs.size()<<" parts:"<<dparts->size() ); 00355 00356 gfxReal partscale = mPartScale; 00357 if(t>1) { 00358 partscale *= (gfxReal)(mPrevs.size()+1-t) / (gfxReal)(mPrevs.size()+1); 00359 } 00360 gfxReal partNormSize = 0.01 * partscale; 00361 //for(size_t i=0; i<mParts.size(); i++) { 00362 for(size_t i=0; i<dparts->size(); i++) { 00363 ParticleObject *p = &( (*dparts)[i] ); // mParts[i]; 00364 00365 if(mShowOnly!=10) { 00366 // 10=show only deleted 00367 if( p->getActive()==false ) continue; 00368 } else { 00369 if( p->getActive()==true ) continue; 00370 } 00371 int type = p->getType(); 00372 if(mShowOnly>0) { 00373 switch(mShowOnly) { 00374 case 1: if(!(type&PART_BUBBLE)) continue; break; 00375 case 2: if(!(type&PART_DROP)) continue; break; 00376 case 3: if(!(type&PART_INTER)) continue; break; 00377 case 4: if(!(type&PART_FLOAT)) continue; break; 00378 case 5: if(!(type&PART_TRACER)) continue; break; 00379 } 00380 } else { 00381 // by default dont display inter 00382 if(type&PART_INTER) continue; 00383 } 00384 00385 pcnt++; 00386 ntlVec3Gfx pnew = p->getPos(); 00387 if(type&PART_FLOAT) { // WARNING same handling for dump! 00388 if(p->getStatus()&PART_IN) { pnew[2] += 0.8; } // offset for display 00389 // add one gridcell offset 00390 //pnew[2] += 1.0; 00391 } 00392 #if LBMDIM==2 00393 pnew[2] += 0.001; // DEBUG 00394 pnew[2] += 0.009; // DEBUG 00395 #endif 00396 00397 ntlVec3Gfx pdir = p->getVel(); 00398 gfxReal plen = normalize( pdir ); 00399 if( plen < 1e-05) pdir = ntlVec3Gfx(-1.0 ,0.0 ,0.0); 00400 ntlVec3Gfx pos = (*mpTrafo) * pnew; 00401 gfxReal partsize = 0.0; 00402 if(debugParts) errMsg("DebugParts"," i"<<i<<" new"<<pnew<<" vel"<<pdir<<" pos="<<pos ); 00403 //if(i==0 &&(debugParts)) errMsg("DebugParts"," i"<<i<<" new"<<pnew[0]<<" pos="<<pos[0]<<" scale="<<scale[0]<<" t="<<trans[0] ); 00404 00405 // value length scaling? 00406 if(mValueScale==1) { 00407 partsize = partscale * plen; 00408 } else if(mValueScale==2) { 00409 // cut off scaling 00410 if(plen > mValueCutoffTop) continue; 00411 if(plen < mValueCutoffBottom) continue; 00412 partsize = partscale * plen; 00413 } else { 00414 partsize = partscale; // no length scaling 00415 } 00416 //if(type&(PART_DROP|PART_BUBBLE)) 00417 partsize *= p->getSize()/5.0; 00418 00419 ntlVec3Gfx pstart( mPartHeadDist *partsize, 0.0, 0.0 ); 00420 ntlVec3Gfx pend ( mPartTailDist *partsize, 0.0, 0.0 ); 00421 gfxReal phi = 0.0; 00422 gfxReal phiD = 2.0*M_PI / (gfxReal)segments; 00423 00424 ntlMat4Gfx cvmat; 00425 cvmat.initId(); 00426 pdir *= -1.0; 00427 ntlVec3Gfx cv1 = pdir; 00428 ntlVec3Gfx cv2 = ntlVec3Gfx(pdir[1], -pdir[0], 0.0); 00429 ntlVec3Gfx cv3 = cross( cv1, cv2); 00430 //? for(int l=0; l<3; l++) { cvmat.value[l][0] = cv1[l]; cvmat.value[l][1] = cv2[l]; cvmat.value[l][2] = cv3[l]; } 00431 pstart = (cvmat * pstart); 00432 pend = (cvmat * pend); 00433 00434 for(int s=0; s<segments; s++) { 00435 ntlVec3Gfx p1( 0.0 ); 00436 ntlVec3Gfx p2( 0.0 ); 00437 00438 gfxReal radscale = partNormSize; 00439 radscale = (partsize+partNormSize)*0.5; 00440 p1[1] += cos(phi) * radscale; 00441 p1[2] += sin(phi) * radscale; 00442 p2[1] += cos(phi + phiD) * radscale; 00443 p2[2] += sin(phi + phiD) * radscale; 00444 ntlVec3Gfx n1 = ntlVec3Gfx( 0.0, cos(phi), sin(phi) ); 00445 ntlVec3Gfx n2 = ntlVec3Gfx( 0.0, cos(phi + phiD), sin(phi + phiD) ); 00446 ntlVec3Gfx ns = n1*0.5 + n2*0.5; 00447 00448 p1 = (cvmat * p1); 00449 p2 = (cvmat * p2); 00450 00451 sceneAddTriangle( pos+pstart, pos+p1, pos+p2, 00452 ns,n1,n2, ntlVec3Gfx(0.0), 1, triangles,vertices,normals ); 00453 sceneAddTriangle( pos+pend , pos+p2, pos+p1, 00454 ns,n2,n1, ntlVec3Gfx(0.0), 1, triangles,vertices,normals ); 00455 00456 phi += phiD; 00457 tris += 2; 00458 } 00459 } 00460 00461 } // t 00462 00463 debMsgStd("ParticleTracer::getTriangles",DM_MSG,"Dumped "<<pcnt<<"/"<<mParts.size()<<" parts, tris:"<<tris<<", showonly:"<<mShowOnly,10); 00464 return; // DEBUG 00465 00466 #endif // ELBEEM_PLUGIN 00467 } 00468 00469 00470 00471