Blender  V2.59
elbeem.cpp
Go to the documentation of this file.
00001 
00004 /******************************************************************************
00005  *
00006  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
00007  * All code distributed as part of El'Beem is covered by the version 2 of the 
00008  * GNU General Public License. See the file COPYING for details.
00009  * Copyright 2003-2006 Nils Thuerey
00010  *
00011  * Main program functions
00012  */
00013 
00014 #include "elbeem.h"
00015 #include "ntl_blenderdumper.h"
00016 extern "C" void elbeemCheckDebugEnv(void);
00017 
00018 #include "ntl_world.h"
00019 #include "ntl_geometrymodel.h"
00020 
00021 /*****************************************************************************/
00022 // region of interest global vars
00023 // currently used by e.g. fsgr solver
00024 double guiRoiSX = 0.0;
00025 double guiRoiSY = 0.0;
00026 double guiRoiSZ = 0.0;
00027 double guiRoiEX = 1.0;
00028 double guiRoiEY = 1.0;
00029 double guiRoiEZ = 1.0;
00030 int guiRoiMaxLev=6, guiRoiMinLev=0;
00031 
00033 ntlWorld *gpWorld = NULL;
00034 
00035 
00036 
00037 // API
00038 
00039 // reset elbeemSimulationSettings struct with defaults
00040 extern "C" 
00041 void elbeemResetSettings(elbeemSimulationSettings *set) {
00042         if(!set) return;
00043   set->version = 3;
00044   set->domainId = 0;
00045         for(int i=0 ; i<3; i++) set->geoStart[i] = 0.0;
00046         for(int i=0 ; i<3; i++) set->geoSize[i] = 1.0;
00047   set->resolutionxyz = 64;
00048   set->previewresxyz = 24;
00049   set->realsize = 1.0;
00050   set->viscosity = 0.000001;
00051 
00052         for(int i=0 ; i<2; i++) set->gravity[i] = 0.0;
00053         set->gravity[2] = -9.81;
00054 
00055   set->animStart = 0; 
00056         set->aniFrameTime = 0.01;
00057         set->noOfFrames = 10;
00058   set->gstar = 0.005;
00059   set->maxRefine = -1;
00060   set->generateParticles = 0.0;
00061   set->numTracerParticles = 0;
00062   strcpy(set->outputPath,"./elbeemdata_");
00063 
00064         set->channelSizeFrameTime=0;
00065         set->channelFrameTime=NULL;
00066         set->channelSizeViscosity=0;
00067         set->channelViscosity=NULL;
00068         set->channelSizeGravity=0;
00069         set->channelGravity=NULL; 
00070 
00071         set->domainobsType= FLUIDSIM_OBSTACLE_NOSLIP;
00072         set->domainobsPartslip= 0.;
00073         set->generateVertexVectors = 0;
00074         set->surfaceSmoothing = 1.;
00075         set->surfaceSubdivs = 1;
00076 
00077         set->farFieldSize = 0.;
00078         set->runsimCallback = NULL;
00079         set->runsimUserData = NULL;
00080 
00081         // init identity
00082         for(int i=0; i<16; i++) set->surfaceTrafo[i] = 0.0;
00083         for(int i=0; i<4; i++) set->surfaceTrafo[i*4+i] = 1.0;
00084 }
00085 
00086 // start fluidsim init
00087 extern "C" 
00088 int elbeemInit() {
00089         setElbeemState( SIMWORLD_INITIALIZING );
00090         setElbeemErrorString("[none]");
00091         resetGlobalColorSetting();
00092 
00093         elbeemCheckDebugEnv();
00094         debMsgStd("performElbeemSimulation",DM_NOTIFY,"El'Beem Simulation Init Start as Plugin, debugLevel:"<<gDebugLevel<<" ...\n", 2);
00095         
00096         // create world object with initial settings
00097         ntlBlenderDumper *elbeem = new ntlBlenderDumper(); 
00098         gpWorld = elbeem;
00099         return 0;
00100 }
00101 
00102 // fluidsim end
00103 extern "C" 
00104 int elbeemFree() {
00105         
00106         return 0;
00107 }
00108 
00109 // start fluidsim init
00110 extern "C" 
00111 int elbeemAddDomain(elbeemSimulationSettings *settings) {
00112         // has to be inited...
00113         if((getElbeemState() == SIMWORLD_INVALID) && (!gpWorld)) { elbeemInit(); }
00114         if(getElbeemState() != SIMWORLD_INITIALIZING) { errFatal("elbeemAddDomain","Unable to init simulation world",SIMWORLD_INITERROR); }
00115         // create domain with given settings
00116         gpWorld->addDomain(settings);
00117         return 0;
00118 }
00119 
00120 // error message access
00121 extern "C" 
00122 void elbeemGetErrorString(char *buffer) {
00123         if(!buffer) return;
00124         strncpy(buffer,getElbeemErrorString(),256);
00125 }
00126 
00127 // reset elbeemMesh struct with zeroes
00128 extern "C" 
00129 void elbeemResetMesh(elbeemMesh *mesh) {
00130         if(!mesh) return;
00131         // init typedef struct elbeemMesh
00132   mesh->type = 0;
00133 
00134         mesh->parentDomainId = 0;
00135 
00136         /* vertices */
00137   mesh->numVertices = 0;
00138         mesh->vertices = NULL;
00139 
00140         mesh->channelSizeVertices = 0;
00141         mesh->channelVertices = NULL;
00142 
00143         /* triangles */
00144         mesh->numTriangles = 0;
00145   mesh->triangles = NULL;
00146 
00147         /* animation channels */
00148         mesh->channelSizeTranslation = 0;
00149         mesh->channelTranslation = NULL;
00150         mesh->channelSizeRotation = 0;
00151         mesh->channelRotation = NULL;
00152         mesh->channelSizeScale = 0;
00153         mesh->channelScale = NULL;
00154 
00155         /* active channel */
00156         mesh->channelSizeActive = 0;
00157         mesh->channelActive = NULL;
00158 
00159         mesh->channelSizeInitialVel = 0;
00160         mesh->channelInitialVel = NULL;
00161 
00162         mesh->localInivelCoords = 0;
00163 
00164         mesh->obstacleType= FLUIDSIM_OBSTACLE_NOSLIP;
00165         mesh->obstaclePartslip= 0.;
00166         mesh->obstacleImpactFactor= 1.;
00167 
00168         mesh->volumeInitType= OB_VOLUMEINIT_VOLUME;
00169 
00170         /* name of the mesh, mostly for debugging */
00171         mesh->name = "[unnamed]";
00172         
00173         /* fluid control settings */
00174         mesh->cpsTimeStart = 0;
00175         mesh->cpsTimeEnd = 0;
00176         mesh->cpsQuality = 0;
00177         
00178         mesh->channelSizeAttractforceStrength = 0;
00179         mesh->channelAttractforceStrength = NULL;
00180         mesh->channelSizeAttractforceRadius = 0;
00181         mesh->channelAttractforceRadius = NULL;
00182         mesh->channelSizeVelocityforceStrength = 0;
00183         mesh->channelVelocityforceStrength = NULL;
00184         mesh->channelSizeVelocityforceRadius = 0;
00185         mesh->channelVelocityforceRadius = NULL;
00186 }
00187 
00188 int globalMeshCounter = 1;
00189 // add mesh as fluidsim object
00190 extern "C" 
00191 int elbeemAddMesh(elbeemMesh *mesh) {
00192         int initType;
00193         if(getElbeemState() != SIMWORLD_INITIALIZING) { errFatal("elbeemAddMesh","World and domain not initialized, call elbeemInit and elbeemAddDomain before...", SIMWORLD_INITERROR); }
00194 
00195         switch(mesh->type) {
00196                 case OB_FLUIDSIM_OBSTACLE: 
00197                         if     (mesh->obstacleType==FLUIDSIM_OBSTACLE_PARTSLIP) initType = FGI_BNDPART; 
00198                         else if(mesh->obstacleType==FLUIDSIM_OBSTACLE_FREESLIP) initType = FGI_BNDFREE; 
00199                         else /*if(mesh->obstacleType==FLUIDSIM_OBSTACLE_NOSLIP)*/ initType = FGI_BNDNO; 
00200                         break;
00201                 case OB_FLUIDSIM_FLUID: initType = FGI_FLUID; break;
00202                 case OB_FLUIDSIM_INFLOW: initType = FGI_MBNDINFLOW; break;
00203                 case OB_FLUIDSIM_OUTFLOW: initType = FGI_MBNDOUTFLOW; break;
00204                 case OB_FLUIDSIM_CONTROL: initType = FGI_CONTROL; break;
00205                 default: return 1; // invalid type
00206         }
00207         
00208         ntlGeometryObjModel *obj = new ntlGeometryObjModel( );
00209         gpWorld->getRenderGlobals()->getSimScene()->addGeoClass( obj );
00210         gpWorld->getRenderGlobals()->getRenderScene()->addGeoClass(obj);
00211         obj->initModel(
00212                         mesh->numVertices, mesh->vertices, mesh->numTriangles, mesh->triangles,
00213                         mesh->channelSizeVertices, mesh->channelVertices );
00214         if(mesh->name) obj->setName(string(mesh->name));
00215         else {
00216                 char meshname[100];
00217                 snprintf(meshname,100,"mesh%04d",globalMeshCounter);
00218                 obj->setName(string(meshname));
00219         }
00220         globalMeshCounter++;
00221         obj->setGeoInitId( mesh->parentDomainId+1 );
00222         obj->setGeoInitIntersect(true);
00223         obj->setGeoInitType(initType);
00224         
00225         // abuse partslip value for control fluid: reverse control keys or not
00226         if(initType == FGI_CONTROL)
00227                 obj->setGeoPartSlipValue(mesh->obstacleType);
00228         else
00229                 obj->setGeoPartSlipValue(mesh->obstaclePartslip);
00230         
00231         obj->setGeoImpactFactor(mesh->obstacleImpactFactor);
00232         
00233         /* fluid control features */
00234         obj->setCpsTimeStart(mesh->cpsTimeStart);
00235         obj->setCpsTimeEnd(mesh->cpsTimeEnd);
00236         obj->setCpsQuality(mesh->cpsQuality);
00237         
00238         if((mesh->volumeInitType<VOLUMEINIT_VOLUME)||(mesh->volumeInitType>VOLUMEINIT_BOTH)) mesh->volumeInitType = VOLUMEINIT_VOLUME;
00239         obj->setVolumeInit(mesh->volumeInitType);
00240         // use channel instead, obj->setInitialVelocity( ntlVec3Gfx(mesh->iniVelocity[0], mesh->iniVelocity[1], mesh->iniVelocity[2]) );
00241         
00242         obj->initChannels(
00243                         mesh->channelSizeTranslation, mesh->channelTranslation, 
00244                         mesh->channelSizeRotation,    mesh->channelRotation, 
00245                         mesh->channelSizeScale,       mesh->channelScale,
00246                         mesh->channelSizeActive,      mesh->channelActive,
00247                         mesh->channelSizeInitialVel,  mesh->channelInitialVel,
00248                         mesh->channelSizeAttractforceStrength,  mesh->channelAttractforceStrength,
00249                         mesh->channelSizeAttractforceRadius,  mesh->channelAttractforceRadius,
00250                         mesh->channelSizeVelocityforceStrength,  mesh->channelVelocityforceStrength,
00251                         mesh->channelSizeVelocityforceRadius,  mesh->channelVelocityforceRadius
00252                 );
00253         obj->setLocalCoordInivel( mesh->localInivelCoords );
00254 
00255         debMsgStd("elbeemAddMesh",DM_MSG,"Added elbeem mesh: "<<obj->getName()<<" type="<<initType<<" "<<obj->getIsAnimated(), 9 );
00256         return 0;
00257 }
00258 
00259 // do the actual simulation
00260 extern "C" 
00261 int elbeemSimulate(void) {
00262         if(!gpWorld) return 1;
00263 
00264         gpWorld->finishWorldInit();
00265         if( isSimworldOk() ) {
00266                 myTime_t timestart = getTime();
00267                 gpWorld->renderAnimation();
00268                 myTime_t timeend = getTime();
00269 
00270                 if(getElbeemState() != SIMWORLD_STOP) {
00271                         // ok, we're done...
00272                         delete gpWorld;
00273                         
00274                         gpWorld = NULL;
00275                         debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 ); 
00276                 } else {
00277                         debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation stopped, time so far: "<<getTimeString(timeend-timestart)<<".", 2 ); 
00278                 }
00279                 return 0;
00280         } 
00281 
00282         // failure...
00283         return 1;
00284 }
00285 
00286 
00287 // continue a previously stopped simulation
00288 extern "C" 
00289 int elbeemContinueSimulation(void) {
00290 
00291         if(getElbeemState() != SIMWORLD_STOP) {
00292                 errMsg("elbeemContinueSimulation","No running simulation found! Aborting...");
00293                 if(gpWorld) delete gpWorld;
00294                 return 1;
00295         }
00296 
00297         myTime_t timestart = getTime();
00298         gpWorld->renderAnimation();
00299         myTime_t timeend = getTime();
00300 
00301         if(getElbeemState() != SIMWORLD_STOP) {
00302                 // ok, we're done...
00303                 delete gpWorld;
00304                 gpWorld = NULL;
00305                 debMsgStd("elbeemContinueSimulation",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 ); 
00306         } else {
00307                 debMsgStd("elbeemContinueSimulation",DM_NOTIFY, "El'Beem simulation stopped, time so far: "<<getTimeString(timeend-timestart)<<".", 2 ); 
00308         }
00309         return 0;
00310 } 
00311 
00312 
00313 // global vector to flag values to remove
00314 vector<int> gKeepVal;
00315 
00316 #define SIMPLIFY_FLOAT_EPSILON (1e-6f)
00317 #define SIMPLIFY_DOUBLE_EPSILON (1e-12f)
00318 #define SFLOATEQ(x,y)  (ABS((x)-(y)) < SIMPLIFY_FLOAT_EPSILON)
00319 #define SDOUBLEEQ(x,y) (ABS((x)-(y)) < SIMPLIFY_DOUBLE_EPSILON)
00320 #define SVECFLOATEQ(x,y)  ( \
00321         (ABS((x)[0]-(y)[0]) < SIMPLIFY_FLOAT_EPSILON) && \
00322         (ABS((x)[1]-(y)[1]) < SIMPLIFY_FLOAT_EPSILON) && \
00323         (ABS((x)[2]-(y)[2]) < SIMPLIFY_FLOAT_EPSILON) )
00324 
00325 // helper function - simplify animation channels
00326 extern "C" 
00327 int elbeemSimplifyChannelFloat(float *channel, int* size) {
00328         bool changed = false;
00329         int nsize = *size;
00330         int orgsize = *size;
00331         if(orgsize<1) return false;
00332         gKeepVal.resize( orgsize );
00333         for(int i=0; i<orgsize; i++) { gKeepVal[i] = true; }
00334         const bool debugSF = false;
00335 
00336         float last = channel[0 + 0];
00337         for(int i=1; i<orgsize; i++) {
00338                 float curr = channel[2*i + 0];
00339                 bool remove = false;
00340                 if(SFLOATEQ(curr,last)) remove = true;
00341                 // dont remove if next value is different
00342                 if((remove)&&(i<orgsize-1)) {
00343                         float next = channel[2*(i+1)+0];
00344                         if(!SFLOATEQ(next,curr)) remove = false;
00345                 }
00346                 if(remove) {
00347                         changed = true;
00348                         gKeepVal[i] = false;
00349                         nsize--;
00350                 }
00351                 if(debugSF) errMsg("elbeemSimplifyChannelFloat","i"<<i<<"/"<<orgsize<<" v"<<channel[ (i*2) + 0 ]<<" t"<<channel[ (i*2) + 1 ]<<"   nsize="<<nsize<<" r"<<remove );
00352                 last = curr;
00353         }
00354 
00355         if(changed) {
00356                 nsize = 1;
00357                 for(int i=1; i<orgsize; i++) {
00358                         if(gKeepVal[i]) {
00359                                 channel[ (nsize*2) + 0 ] = channel[ (i*2) + 0 ];
00360                                 channel[ (nsize*2) + 1 ] = channel[ (i*2) + 1 ];
00361                                 nsize++;
00362                         }
00363                 }
00364                 *size = nsize;
00365         }
00366 
00367         if(debugSF) for(int i=1; i<nsize; i++) {
00368                 errMsg("elbeemSimplifyChannelFloat","n i"<<i<<"/"<<nsize<<" v"<<channel[ (i*2) + 0 ]<<" t"<<channel[ (i*2) + 1 ] );
00369         }
00370 
00371         return changed;
00372 }
00373 
00374 extern "C" 
00375 int elbeemSimplifyChannelVec3(float *channel, int* size) {
00376         bool changed = false;
00377         int nsize = *size;
00378         int orgsize = *size;
00379         if(orgsize<1) return false;
00380         gKeepVal.resize( orgsize );
00381         for(int i=0; i<orgsize; i++) { gKeepVal[i] = true; }
00382         const bool debugVF = false;
00383 
00384         ntlVec3f last( channel[0 + 0], channel[0 + 1], channel[0 + 2] );
00385         for(int i=1; i<orgsize; i++) {
00386                 ntlVec3f curr( channel[4*i + 0], channel[4*i + 1], channel[4*i + 2]);
00387                 bool remove = false;
00388                 if(SVECFLOATEQ(curr,last)) remove = true;
00389                 // dont remove if next value is different
00390                 if((remove)&&(i<orgsize-1)) {
00391                         ntlVec3f next( channel[4*(i+1)+0], channel[4*(i+1)+1], channel[4*(i+1)+2]);
00392                         if(!SVECFLOATEQ(next,curr)) remove = false;
00393                 }
00394                 if(remove) {
00395                         changed = true;
00396                         gKeepVal[i] = false;
00397                         nsize--;
00398                 }
00399                 if(debugVF) errMsg("elbeemSimplifyChannelVec3","i"<<i<<"/"<<orgsize<<" v"<<
00400                                 channel[ (i*4) + 0 ]<<","<< channel[ (i*4) + 1 ]<<","<< channel[ (i*4) + 2 ]<<
00401                                 " t"<<channel[ (i*4) + 3 ]<<"   nsize="<<nsize<<" r"<<remove );
00402                 last = curr;
00403         }
00404 
00405         if(changed) {
00406                 nsize = 1;
00407                 for(int i=1; i<orgsize; i++) {
00408                         if(gKeepVal[i]) {
00409                                 for(int j=0; j<4; j++){ channel[ (nsize*4) + j ] = channel[ (i*4) + j ]; }
00410                                 nsize++;
00411                         }
00412                 }
00413                 *size = nsize;
00414         }
00415 
00416         if(debugVF) for(int i=1; i<nsize; i++) {
00417                 errMsg("elbeemSimplifyChannelVec3","n i"<<i<<"/"<<nsize<<" v"<<
00418                                 channel[ (i*4) + 0 ]<<","<< channel[ (i*4) + 1 ]<<","<< channel[ (i*4) + 2 ]<<
00419                                 " t"<<channel[ (i*4) + 3 ] );
00420         }
00421 
00422         return changed;
00423 }
00424 
00425 
00426 #undef SIMPLIFY_FLOAT_EPSILON
00427 #undef SIMPLIFY_DOUBLE_EPSILON
00428 #undef SFLOATEQ
00429 #undef SDOUBLEEQ
00430