Filling in branch from trunk
[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 // 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