intern/elbeem/intern/elbeem.cpp

   1 /******************************************************************************
   2  *
   3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
   4  * All code distributed as part of El'Beem is covered by the version 2 of the
   5  * GNU General Public License. See the file COPYING for details.
   6  * Copyright 2003-2006 Nils Thuerey
   7  *
   8  * Main program functions
   9  */
  10
  11 #include "elbeem.h"
  12 #include "ntl_blenderdumper.h"
  13 extern "C" void elbeemCheckDebugEnv(void);
  14
  15 #include "ntl_world.h"
  16 #include "ntl_geometrymodel.h"
  17
  18 /*****************************************************************************/
  19 // region of interest global vars
  20 // currently used by e.g. fsgr solver
  21 double guiRoiSX = 0.0;
  22 double guiRoiSY = 0.0;
  23 double guiRoiSZ = 0.0;
  24 double guiRoiEX = 1.0;
  25 double guiRoiEY = 1.0;
  26 double guiRoiEZ = 1.0;
  27 int guiRoiMaxLev=6, guiRoiMinLev=0;
  28
  29 //! global raytracer pointer (=world)
  30 ntlWorld *gpWorld = NULL;
  31
  32
  33
  34 // API
  35
  36 // reset elbeemSimulationSettings struct with defaults
  37 extern "C"
  38 void elbeemResetSettings(elbeemSimulationSettings *set) {
  39         if(!set) return;
  40   set->version = 3;
  41   set->domainId = 0;
  42         for(int i=0 ; i<3; i++) set->geoStart[i] = 0.0;
  43         for(int i=0 ; i<3; i++) set->geoSize[i] = 1.0;
  44   set->resolutionxyz = 64;
  45   set->previewresxyz = 24;
  46   set->realsize = 1.0;
  47   set->viscosity = 0.000001;
  48
  49         for(int i=0 ; i<2; i++) set->gravity[i] = 0.0;
  50         set->gravity[2] = -9.81;
  51
  52   set->animStart = 0;
  53         set->aniFrameTime = 0.01;
  54         set->noOfFrames = 10;
  55   set->gstar = 0.005;
  56   set->maxRefine = -1;
  57   set->generateParticles = 0.0;
  58   set->numTracerParticles = 0;
  59   strcpy(set->outputPath,"./elbeemdata_");
  60
  61         set->channelSizeFrameTime=0;
  62         set->channelFrameTime=NULL;
  63         set->channelSizeViscosity=0;
  64         set->channelViscosity=NULL;
  65         set->channelSizeGravity=0;
  66         set->channelGravity=NULL;
  67
  68         set->domainobsType= FLUIDSIM_OBSTACLE_NOSLIP;
  69         set->domainobsPartslip= 0.;
  70         set->generateVertexVectors = 0;
  71         set->surfaceSmoothing = 1.;
  72         set->surfaceSubdivs = 1;
  73
  74         set->farFieldSize = 0.;
  75         set->runsimCallback = NULL;
  76         set->runsimUserData = NULL;
  77
  78         // init identity
  79         for(int i=0; i<16; i++) set->surfaceTrafo[i] = 0.0;
  80         for(int i=0; i<4; i++) set->surfaceTrafo[i*4+i] = 1.0;
  81 }
  82
  83 // start fluidsim init
  84 extern "C"
  85 int elbeemInit() {
  86         setElbeemState( SIMWORLD_INITIALIZING );
  87         setElbeemErrorString("[none]");
  88         resetGlobalColorSetting();
  89
  90         elbeemCheckDebugEnv();
  91         debMsgStd("performElbeemSimulation",DM_NOTIFY,"El'Beem Simulation Init Start as Plugin, debugLevel:"<<gDebugLevel<<" ...\n", 2);
  92
  93         // create world object with initial settings
  94         ntlBlenderDumper *elbeem = new ntlBlenderDumper();
  95         gpWorld = elbeem;
  96         return 0;
  97 }
  98
  99 // fluidsim end
 100 extern "C"
 101 int elbeemFree() {
 102
 103         return 0;
 104 }
 105
 106 // start fluidsim init
 107 extern "C"
 108 int elbeemAddDomain(elbeemSimulationSettings *settings) {
 109         // has to be inited...
 110         if((getElbeemState() == SIMWORLD_INVALID) && (!gpWorld)) { elbeemInit(); }
 111         if(getElbeemState() != SIMWORLD_INITIALIZING) { errFatal("elbeemAddDomain","Unable to init simulation world",SIMWORLD_INITERROR); }
 112         // create domain with given settings
 113         gpWorld->addDomain(settings);
 114         return 0;
 115 }
 116
 117 // error message access
 118 extern "C"
 119 void elbeemGetErrorString(char *buffer) {
 120         if(!buffer) return;
 121         strncpy(buffer,getElbeemErrorString(),256);
 122 }
 123
 124 // reset elbeemMesh struct with zeroes
 125 extern "C"
 126 void elbeemResetMesh(elbeemMesh *mesh) {
 127         if(!mesh) return;
 128         // init typedef struct elbeemMesh
 129   mesh->type = 0;
 130
 131         mesh->parentDomainId = 0;
 132
 133         /* vertices */
 134   mesh->numVertices = 0;
 135         mesh->vertices = NULL;
 136
 137         mesh->channelSizeVertices = 0;
 138         mesh->channelVertices = NULL;
 139
 140         /* triangles */
 141         mesh->numTriangles = 0;
 142   mesh->triangles = NULL;
 143
 144         /* animation channels */
 145         mesh->channelSizeTranslation = 0;
 146         mesh->channelTranslation = NULL;
 147         mesh->channelSizeRotation = 0;
 148         mesh->channelRotation = NULL;
 149         mesh->channelSizeScale = 0;
 150         mesh->channelScale = NULL;
 151
 152         /* active channel */
 153         mesh->channelSizeActive = 0;
 154         mesh->channelActive = NULL;
 155
 156         mesh->channelSizeInitialVel = 0;
 157         mesh->channelInitialVel = NULL;
 158
 159         mesh->localInivelCoords = 0;
 160
 161         mesh->obstacleType= FLUIDSIM_OBSTACLE_NOSLIP;
 162         mesh->obstaclePartslip= 0.;
 163         mesh->obstacleImpactFactor= 1.;
 164
 165         mesh->volumeInitType= OB_VOLUMEINIT_VOLUME;
 166
 167         /* name of the mesh, mostly for debugging */
 168         mesh->name = "[unnamed]";
 169
 170         /* fluid control settings */
 171         mesh->cpsTimeStart = 0;
 172         mesh->cpsTimeEnd = 0;
 173         mesh->cpsQuality = 0;
 174
 175         mesh->channelSizeAttractforceStrength = 0;
 176         mesh->channelAttractforceStrength = NULL;
 177         mesh->channelSizeAttractforceRadius = 0;
 178         mesh->channelAttractforceRadius = NULL;
 179         mesh->channelSizeVelocityforceStrength = 0;
 180         mesh->channelVelocityforceStrength = NULL;
 181         mesh->channelSizeVelocityforceRadius = 0;
 182         mesh->channelVelocityforceRadius = NULL;
 183 }
 184
 185 int globalMeshCounter = 1;
 186 // add mesh as fluidsim object
 187 extern "C"
 188 int elbeemAddMesh(elbeemMesh *mesh) {
 189         int initType;
 190         if(getElbeemState() != SIMWORLD_INITIALIZING) { errFatal("elbeemAddMesh","World and domain not initialized, call elbeemInit and elbeemAddDomain before...", SIMWORLD_INITERROR); }
 191
 192         switch(mesh->type) {
 193                 case OB_FLUIDSIM_OBSTACLE:
 194                         if     (mesh->obstacleType==FLUIDSIM_OBSTACLE_PARTSLIP) initType = FGI_BNDPART;
 195                         else if(mesh->obstacleType==FLUIDSIM_OBSTACLE_FREESLIP) initType = FGI_BNDFREE;
 196                         else /*if(mesh->obstacleType==FLUIDSIM_OBSTACLE_NOSLIP)*/ initType = FGI_BNDNO;
 197                         break;
 198                 case OB_FLUIDSIM_FLUID: initType = FGI_FLUID; break;
 199                 case OB_FLUIDSIM_INFLOW: initType = FGI_MBNDINFLOW; break;
 200                 case OB_FLUIDSIM_OUTFLOW: initType = FGI_MBNDOUTFLOW; break;
 201                 case OB_FLUIDSIM_CONTROL: initType = FGI_CONTROL; break;
 202                 default: return 1; // invalid type
 203         }
 204
 205         ntlGeometryObjModel *obj = new ntlGeometryObjModel( );
 206         gpWorld->getRenderGlobals()->getSimScene()->addGeoClass( obj );
 207         obj->initModel(
 208                         mesh->numVertices, mesh->vertices, mesh->numTriangles, mesh->triangles,
 209                         mesh->channelSizeVertices, mesh->channelVertices );
 210         if(mesh->name) obj->setName(string(mesh->name));
 211         else {
 212                 char meshname[100];
 213                 snprintf(meshname,100,"mesh%04d",globalMeshCounter);
 214                 obj->setName(string(meshname));
 215         }
 216         globalMeshCounter++;
 217         obj->setGeoInitId( mesh->parentDomainId+1 );
 218         obj->setGeoInitIntersect(true);
 219         obj->setGeoInitType(initType);
 220
 221         // abuse partslip value for control fluid: reverse control keys or not
 222         if(initType == FGI_CONTROL)
 223                 obj->setGeoPartSlipValue(mesh->obstacleType);
 224         else
 225                 obj->setGeoPartSlipValue(mesh->obstaclePartslip);
 226
 227         obj->setGeoImpactFactor(mesh->obstacleImpactFactor);
 228
 229         /* fluid control features */
 230         obj->setCpsTimeStart(mesh->cpsTimeStart);
 231         obj->setCpsTimeEnd(mesh->cpsTimeEnd);
 232         obj->setCpsQuality(mesh->cpsQuality);
 233
 234         if((mesh->volumeInitType<VOLUMEINIT_VOLUME)||(mesh->volumeInitType>VOLUMEINIT_BOTH)) mesh->volumeInitType = VOLUMEINIT_VOLUME;
 235         obj->setVolumeInit(mesh->volumeInitType);
 236         // use channel instead, obj->setInitialVelocity( ntlVec3Gfx(mesh->iniVelocity[0], mesh->iniVelocity[1], mesh->iniVelocity[2]) );
 237
 238         obj->initChannels(
 239                         mesh->channelSizeTranslation, mesh->channelTranslation,
 240                         mesh->channelSizeRotation,    mesh->channelRotation,
 241                         mesh->channelSizeScale,       mesh->channelScale,
 242                         mesh->channelSizeActive,      mesh->channelActive,
 243                         mesh->channelSizeInitialVel,  mesh->channelInitialVel,
 244                         mesh->channelSizeAttractforceStrength,  mesh->channelAttractforceStrength,
 245                         mesh->channelSizeAttractforceRadius,  mesh->channelAttractforceRadius,
 246                         mesh->channelSizeVelocityforceStrength,  mesh->channelVelocityforceStrength,
 247                         mesh->channelSizeVelocityforceRadius,  mesh->channelVelocityforceRadius
 248                 );
 249         obj->setLocalCoordInivel( mesh->localInivelCoords );
 250
 251         debMsgStd("elbeemAddMesh",DM_MSG,"Added elbeem mesh: "<<obj->getName()<<" type="<<initType<<" "<<obj->getIsAnimated(), 9 );
 252         return 0;
 253 }
 254
 255 // do the actual simulation
 256 extern "C"
 257 int elbeemSimulate(void) {
 258         if(!gpWorld) return 1;
 259
 260         gpWorld->finishWorldInit();
 261         if( isSimworldOk() ) {
 262                 myTime_t timestart = getTime();
 263                 gpWorld->renderAnimation();
 264                 myTime_t timeend = getTime();
 265
 266                 if(getElbeemState() != SIMWORLD_STOP) {
 267                         // ok, we're done...
 268                         delete gpWorld;
 269
 270                         gpWorld = NULL;
 271                         debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 );
 272                 } else {
 273                         debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation stopped, time so far: "<<getTimeString(timeend-timestart)<<".", 2 );
 274                 }
 275                 return 0;
 276         }
 277
 278         // failure...
 279         return 1;
 280 }
 281
 282
 283 // continue a previously stopped simulation
 284 extern "C"
 285 int elbeemContinueSimulation(void) {
 286
 287         if(getElbeemState() != SIMWORLD_STOP) {
 288                 errMsg("elbeemContinueSimulation","No running simulation found! Aborting...");
 289                 if(gpWorld) delete gpWorld;
 290                 return 1;
 291         }
 292
 293         myTime_t timestart = getTime();
 294         gpWorld->renderAnimation();
 295         myTime_t timeend = getTime();
 296
 297         if(getElbeemState() != SIMWORLD_STOP) {
 298                 // ok, we're done...
 299                 delete gpWorld;
 300                 gpWorld = NULL;
 301                 debMsgStd("elbeemContinueSimulation",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 );
 302         } else {
 303                 debMsgStd("elbeemContinueSimulation",DM_NOTIFY, "El'Beem simulation stopped, time so far: "<<getTimeString(timeend-timestart)<<".", 2 );
 304         }
 305         return 0;
 306 }
 307
 308
 309 // global vector to flag values to remove
 310 vector<int> gKeepVal;
 311
 312 #define SIMPLIFY_FLOAT_EPSILON (1e-6f)
 313 #define SIMPLIFY_DOUBLE_EPSILON (1e-12f)
 314 #define SFLOATEQ(x,y)  (ABS((x)-(y)) < SIMPLIFY_FLOAT_EPSILON)
 315 #define SDOUBLEEQ(x,y) (ABS((x)-(y)) < SIMPLIFY_DOUBLE_EPSILON)
 316 #define SVECFLOATEQ(x,y)  ( \
 317         (ABS((x)[0]-(y)[0]) < SIMPLIFY_FLOAT_EPSILON) && \
 318         (ABS((x)[1]-(y)[1]) < SIMPLIFY_FLOAT_EPSILON) && \
 319         (ABS((x)[2]-(y)[2]) < SIMPLIFY_FLOAT_EPSILON) )
 320
 321 // helper function - simplify animation channels
 322 extern "C"
 323 int elbeemSimplifyChannelFloat(float *channel, int* size) {
 324         bool changed = false;
 325         int nsize = *size;
 326         int orgsize = *size;
 327         if(orgsize<1) return false;
 328         gKeepVal.resize( orgsize );
 329         for(int i=0; i<orgsize; i++) { gKeepVal[i] = true; }
 330         const bool debugSF = false;
 331
 332         float last = channel[0 + 0];
 333         for(int i=1; i<orgsize; i++) {
 334                 float curr = channel[2*i + 0];
 335                 bool remove = false;
 336                 if(SFLOATEQ(curr,last)) remove = true;
 337                 // dont remove if next value is different
 338                 if((remove)&&(i<orgsize-1)) {
 339                         float next = channel[2*(i+1)+0];
 340                         if(!SFLOATEQ(next,curr)) remove = false;
 341                 }
 342                 if(remove) {
 343                         changed = true;
 344                         gKeepVal[i] = false;
 345                         nsize--;
 346                 }
 347                 if(debugSF) errMsg("elbeemSimplifyChannelFloat","i"<<i<<"/"<<orgsize<<" v"<<channel[ (i*2) + 0 ]<<" t"<<channel[ (i*2) + 1 ]<<"   nsize="<<nsize<<" r"<<remove );
 348                 last = curr;
 349         }
 350
 351         if(changed) {
 352                 nsize = 1;
 353                 for(int i=1; i<orgsize; i++) {
 354                         if(gKeepVal[i]) {
 355                                 channel[ (nsize*2) + 0 ] = channel[ (i*2) + 0 ];
 356                                 channel[ (nsize*2) + 1 ] = channel[ (i*2) + 1 ];
 357                                 nsize++;
 358                         }
 359                 }
 360                 *size = nsize;
 361         }
 362
 363         if(debugSF) for(int i=1; i<nsize; i++) {
 364                 errMsg("elbeemSimplifyChannelFloat","n i"<<i<<"/"<<nsize<<" v"<<channel[ (i*2) + 0 ]<<" t"<<channel[ (i*2) + 1 ] );
 365         }
 366
 367         return changed;
 368 }
 369
 370 extern "C"
 371 int elbeemSimplifyChannelVec3(float *channel, int* size) {
 372         bool changed = false;
 373         int nsize = *size;
 374         int orgsize = *size;
 375         if(orgsize<1) return false;
 376         gKeepVal.resize( orgsize );
 377         for(int i=0; i<orgsize; i++) { gKeepVal[i] = true; }
 378         const bool debugVF = false;
 379
 380         ntlVec3f last( channel[0 + 0], channel[0 + 1], channel[0 + 2] );
 381         for(int i=1; i<orgsize; i++) {
 382                 ntlVec3f curr( channel[4*i + 0], channel[4*i + 1], channel[4*i + 2]);
 383                 bool remove = false;
 384                 if(SVECFLOATEQ(curr,last)) remove = true;
 385                 // dont remove if next value is different
 386                 if((remove)&&(i<orgsize-1)) {
 387                         ntlVec3f next( channel[4*(i+1)+0], channel[4*(i+1)+1], channel[4*(i+1)+2]);
 388                         if(!SVECFLOATEQ(next,curr)) remove = false;
 389                 }
 390                 if(remove) {
 391                         changed = true;
 392                         gKeepVal[i] = false;
 393                         nsize--;
 394                 }
 395                 if(debugVF) errMsg("elbeemSimplifyChannelVec3","i"<<i<<"/"<<orgsize<<" v"<<
 396                                 channel[ (i*4) + 0 ]<<","<< channel[ (i*4) + 1 ]<<","<< channel[ (i*4) + 2 ]<<
 397                                 " t"<<channel[ (i*4) + 3 ]<<"   nsize="<<nsize<<" r"<<remove );
 398                 last = curr;
 399         }
 400
 401         if(changed) {
 402                 nsize = 1;
 403                 for(int i=1; i<orgsize; i++) {
 404                         if(gKeepVal[i]) {
 405                                 for(int j=0; j<4; j++){ channel[ (nsize*4) + j ] = channel[ (i*4) + j ]; }
 406                                 nsize++;
 407                         }
 408                 }
 409                 *size = nsize;
 410         }
 411
 412         if(debugVF) for(int i=1; i<nsize; i++) {
 413                 errMsg("elbeemSimplifyChannelVec3","n i"<<i<<"/"<<nsize<<" v"<<
 414                                 channel[ (i*4) + 0 ]<<","<< channel[ (i*4) + 1 ]<<","<< channel[ (i*4) + 2 ]<<
 415                                 " t"<<channel[ (i*4) + 3 ] );
 416         }
 417
 418         return changed;
 419 }
 420
 421
 422 #undef SIMPLIFY_FLOAT_EPSILON
 423 #undef SIMPLIFY_DOUBLE_EPSILON
 424 #undef SFLOATEQ
 425 #undef SDOUBLEEQ
 426