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