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