svn merge -r 16866:17042 https://svn.blender.org/svnroot/bf-blender/trunk/blender
[blender.git] / 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         gpWorld->getRenderGlobals()->getRenderScene()->addGeoClass(obj);
208         obj->initModel(
209                         mesh->numVertices, mesh->vertices, mesh->numTriangles, mesh->triangles,
210                         mesh->channelSizeVertices, mesh->channelVertices );
211         if(mesh->name) obj->setName(string(mesh->name));
212         else {
213                 char meshname[100];
214                 snprintf(meshname,100,"mesh%04d",globalMeshCounter);
215                 obj->setName(string(meshname));
216         }
217         globalMeshCounter++;
218         obj->setGeoInitId( mesh->parentDomainId+1 );
219         obj->setGeoInitIntersect(true);
220         obj->setGeoInitType(initType);
221         
222         // abuse partslip value for control fluid: reverse control keys or not
223         if(initType == FGI_CONTROL)
224                 obj->setGeoPartSlipValue(mesh->obstacleType);
225         else
226                 obj->setGeoPartSlipValue(mesh->obstaclePartslip);
227         
228         obj->setGeoImpactFactor(mesh->obstacleImpactFactor);
229         
230         /* fluid control features */
231         obj->setCpsTimeStart(mesh->cpsTimeStart);
232         obj->setCpsTimeEnd(mesh->cpsTimeEnd);
233         obj->setCpsQuality(mesh->cpsQuality);
234         
235         if((mesh->volumeInitType<VOLUMEINIT_VOLUME)||(mesh->volumeInitType>VOLUMEINIT_BOTH)) mesh->volumeInitType = VOLUMEINIT_VOLUME;
236         obj->setVolumeInit(mesh->volumeInitType);
237         // use channel instead, obj->setInitialVelocity( ntlVec3Gfx(mesh->iniVelocity[0], mesh->iniVelocity[1], mesh->iniVelocity[2]) );
238         
239         obj->initChannels(
240                         mesh->channelSizeTranslation, mesh->channelTranslation, 
241                         mesh->channelSizeRotation,    mesh->channelRotation, 
242                         mesh->channelSizeScale,       mesh->channelScale,
243                         mesh->channelSizeActive,      mesh->channelActive,
244                         mesh->channelSizeInitialVel,  mesh->channelInitialVel,
245                         mesh->channelSizeAttractforceStrength,  mesh->channelAttractforceStrength,
246                         mesh->channelSizeAttractforceRadius,  mesh->channelAttractforceRadius,
247                         mesh->channelSizeVelocityforceStrength,  mesh->channelVelocityforceStrength,
248                         mesh->channelSizeVelocityforceRadius,  mesh->channelVelocityforceRadius
249                 );
250         obj->setLocalCoordInivel( mesh->localInivelCoords );
251
252         debMsgStd("elbeemAddMesh",DM_MSG,"Added elbeem mesh: "<<obj->getName()<<" type="<<initType<<" "<<obj->getIsAnimated(), 9 );
253         return 0;
254 }
255
256 // do the actual simulation
257 extern "C" 
258 int elbeemSimulate(void) {
259         if(!gpWorld) return 1;
260
261         gpWorld->finishWorldInit();
262         if( isSimworldOk() ) {
263                 myTime_t timestart = getTime();
264                 gpWorld->renderAnimation();
265                 myTime_t timeend = getTime();
266
267                 if(getElbeemState() != SIMWORLD_STOP) {
268                         // ok, we're done...
269                         delete gpWorld;
270                         
271                         gpWorld = NULL;
272                         debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 ); 
273                 } else {
274                         debMsgStd("elbeemSimulate",DM_NOTIFY, "El'Beem simulation stopped, time so far: "<<getTimeString(timeend-timestart)<<".", 2 ); 
275                 }
276                 return 0;
277         } 
278
279         // failure...
280         return 1;
281 }
282
283
284 // continue a previously stopped simulation
285 extern "C" 
286 int elbeemContinueSimulation(void) {
287
288         if(getElbeemState() != SIMWORLD_STOP) {
289                 errMsg("elbeemContinueSimulation","No running simulation found! Aborting...");
290                 if(gpWorld) delete gpWorld;
291                 return 1;
292         }
293
294         myTime_t timestart = getTime();
295         gpWorld->renderAnimation();
296         myTime_t timeend = getTime();
297
298         if(getElbeemState() != SIMWORLD_STOP) {
299                 // ok, we're done...
300                 delete gpWorld;
301                 gpWorld = NULL;
302                 debMsgStd("elbeemContinueSimulation",DM_NOTIFY, "El'Beem simulation done, time: "<<getTimeString(timeend-timestart)<<".\n", 2 ); 
303         } else {
304                 debMsgStd("elbeemContinueSimulation",DM_NOTIFY, "El'Beem simulation stopped, time so far: "<<getTimeString(timeend-timestart)<<".", 2 ); 
305         }
306         return 0;
307
308
309
310 // global vector to flag values to remove
311 vector<int> gKeepVal;
312
313 #define SIMPLIFY_FLOAT_EPSILON (1e-6f)
314 #define SIMPLIFY_DOUBLE_EPSILON (1e-12f)
315 #define SFLOATEQ(x,y)  (ABS((x)-(y)) < SIMPLIFY_FLOAT_EPSILON)
316 #define SDOUBLEEQ(x,y) (ABS((x)-(y)) < SIMPLIFY_DOUBLE_EPSILON)
317 #define SVECFLOATEQ(x,y)  ( \
318         (ABS((x)[0]-(y)[0]) < SIMPLIFY_FLOAT_EPSILON) && \
319         (ABS((x)[1]-(y)[1]) < SIMPLIFY_FLOAT_EPSILON) && \
320         (ABS((x)[2]-(y)[2]) < SIMPLIFY_FLOAT_EPSILON) )
321
322 // helper function - simplify animation channels
323 extern "C" 
324 int elbeemSimplifyChannelFloat(float *channel, int* size) {
325         bool changed = false;
326         int nsize = *size;
327         int orgsize = *size;
328         if(orgsize<1) return false;
329         gKeepVal.resize( orgsize );
330         for(int i=0; i<orgsize; i++) { gKeepVal[i] = true; }
331         const bool debugSF = false;
332
333         float last = channel[0 + 0];
334         for(int i=1; i<orgsize; i++) {
335                 float curr = channel[2*i + 0];
336                 bool remove = false;
337                 if(SFLOATEQ(curr,last)) remove = true;
338                 // dont remove if next value is different
339                 if((remove)&&(i<orgsize-1)) {
340                         float next = channel[2*(i+1)+0];
341                         if(!SFLOATEQ(next,curr)) remove = false;
342                 }
343                 if(remove) {
344                         changed = true;
345                         gKeepVal[i] = false;
346                         nsize--;
347                 }
348                 if(debugSF) errMsg("elbeemSimplifyChannelFloat","i"<<i<<"/"<<orgsize<<" v"<<channel[ (i*2) + 0 ]<<" t"<<channel[ (i*2) + 1 ]<<"   nsize="<<nsize<<" r"<<remove );
349                 last = curr;
350         }
351
352         if(changed) {
353                 nsize = 1;
354                 for(int i=1; i<orgsize; i++) {
355                         if(gKeepVal[i]) {
356                                 channel[ (nsize*2) + 0 ] = channel[ (i*2) + 0 ];
357                                 channel[ (nsize*2) + 1 ] = channel[ (i*2) + 1 ];
358                                 nsize++;
359                         }
360                 }
361                 *size = nsize;
362         }
363
364         if(debugSF) for(int i=1; i<nsize; i++) {
365                 errMsg("elbeemSimplifyChannelFloat","n i"<<i<<"/"<<nsize<<" v"<<channel[ (i*2) + 0 ]<<" t"<<channel[ (i*2) + 1 ] );
366         }
367
368         return changed;
369 }
370
371 extern "C" 
372 int elbeemSimplifyChannelVec3(float *channel, int* size) {
373         bool changed = false;
374         int nsize = *size;
375         int orgsize = *size;
376         if(orgsize<1) return false;
377         gKeepVal.resize( orgsize );
378         for(int i=0; i<orgsize; i++) { gKeepVal[i] = true; }
379         const bool debugVF = false;
380
381         ntlVec3f last( channel[0 + 0], channel[0 + 1], channel[0 + 2] );
382         for(int i=1; i<orgsize; i++) {
383                 ntlVec3f curr( channel[4*i + 0], channel[4*i + 1], channel[4*i + 2]);
384                 bool remove = false;
385                 if(SVECFLOATEQ(curr,last)) remove = true;
386                 // dont remove if next value is different
387                 if((remove)&&(i<orgsize-1)) {
388                         ntlVec3f next( channel[4*(i+1)+0], channel[4*(i+1)+1], channel[4*(i+1)+2]);
389                         if(!SVECFLOATEQ(next,curr)) remove = false;
390                 }
391                 if(remove) {
392                         changed = true;
393                         gKeepVal[i] = false;
394                         nsize--;
395                 }
396                 if(debugVF) errMsg("elbeemSimplifyChannelVec3","i"<<i<<"/"<<orgsize<<" v"<<
397                                 channel[ (i*4) + 0 ]<<","<< channel[ (i*4) + 1 ]<<","<< channel[ (i*4) + 2 ]<<
398                                 " t"<<channel[ (i*4) + 3 ]<<"   nsize="<<nsize<<" r"<<remove );
399                 last = curr;
400         }
401
402         if(changed) {
403                 nsize = 1;
404                 for(int i=1; i<orgsize; i++) {
405                         if(gKeepVal[i]) {
406                                 for(int j=0; j<4; j++){ channel[ (nsize*4) + j ] = channel[ (i*4) + j ]; }
407                                 nsize++;
408                         }
409                 }
410                 *size = nsize;
411         }
412
413         if(debugVF) for(int i=1; i<nsize; i++) {
414                 errMsg("elbeemSimplifyChannelVec3","n i"<<i<<"/"<<nsize<<" v"<<
415                                 channel[ (i*4) + 0 ]<<","<< channel[ (i*4) + 1 ]<<","<< channel[ (i*4) + 2 ]<<
416                                 " t"<<channel[ (i*4) + 3 ] );
417         }
418
419         return changed;
420 }
421
422
423 #undef SIMPLIFY_FLOAT_EPSILON
424 #undef SIMPLIFY_DOUBLE_EPSILON
425 #undef SFLOATEQ
426 #undef SDOUBLEEQ
427