Warning fixes, one actual bug found in sequencer sound wave drawing. Also
[blender-staging.git] / intern / smoke / intern / FLUID_3D.cpp
1 //////////////////////////////////////////////////////////////////////
2 // This file is part of Wavelet Turbulence.
3 //
4 // Wavelet Turbulence is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 //
9 // Wavelet Turbulence is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
16 //
17 // Copyright 2008 Theodore Kim and Nils Thuerey
18 //
19 // FLUID_3D.cpp: implementation of the FLUID_3D class.
20 //
21 //////////////////////////////////////////////////////////////////////
22 // Heavy parallel optimization done. Many of the old functions now
23 // take begin and end parameters and process only specified part of the data.
24 // Some functions were divided into multiple ones.
25 //              - MiikaH
26 //////////////////////////////////////////////////////////////////////
27
28 #include "FLUID_3D.h"
29 #include "IMAGE.h"
30 #include <INTERPOLATE.h>
31 #include "SPHERE.h"
32 #include <zlib.h>
33
34 #if PARALLEL==1
35 #include <omp.h>
36 #endif // PARALLEL 
37
38 // boundary conditions of the fluid domain
39 #define DOMAIN_BC_FRONT  0 // z
40 #define DOMAIN_BC_TOP    1 // y
41 #define DOMAIN_BC_LEFT   1 // x
42 #define DOMAIN_BC_BACK   DOMAIN_BC_FRONT
43 #define DOMAIN_BC_BOTTOM DOMAIN_BC_TOP
44 #define DOMAIN_BC_RIGHT  DOMAIN_BC_LEFT
45
46 //////////////////////////////////////////////////////////////////////
47 // Construction/Destruction
48 //////////////////////////////////////////////////////////////////////
49
50 FLUID_3D::FLUID_3D(int *res, float *p0, float dt) :
51         _xRes(res[0]), _yRes(res[1]), _zRes(res[2]), _res(0.0f), _dt(dt)
52 {
53         // set simulation consts
54         // _dt = dt; // 0.10
55         
56         // start point of array
57         _p0[0] = p0[0];
58         _p0[1] = p0[1];
59         _p0[2] = p0[2];
60
61         _iterations = 100;
62         _tempAmb = 0; 
63         _heatDiffusion = 1e-3;
64         _vorticityEps = 2.0;
65         _totalTime = 0.0f;
66         _totalSteps = 0;
67         _res = Vec3Int(_xRes,_yRes,_zRes);
68         _maxRes = MAX3(_xRes, _yRes, _zRes);
69         
70         // initialize wavelet turbulence
71         /*
72         if(amplify)
73                 _wTurbulence = new WTURBULENCE(_res[0],_res[1],_res[2], amplify, noisetype);
74         else
75                 _wTurbulence = NULL;
76         */
77         
78         // scale the constants according to the refinement of the grid
79         _dx = 1.0f / (float)_maxRes;
80         float scaling = 64.0f / _maxRes;
81         scaling = (scaling < 1.0f) ? 1.0f : scaling;
82         _vorticityEps /= scaling;
83
84         // allocate arrays
85         _totalCells   = _xRes * _yRes * _zRes;
86         _slabSize = _xRes * _yRes;
87         _xVelocity    = new float[_totalCells];
88         _yVelocity    = new float[_totalCells];
89         _zVelocity    = new float[_totalCells];
90         _xVelocityOld = new float[_totalCells];
91         _yVelocityOld = new float[_totalCells];
92         _zVelocityOld = new float[_totalCells];
93         _xForce       = new float[_totalCells];
94         _yForce       = new float[_totalCells];
95         _zForce       = new float[_totalCells];
96         _density      = new float[_totalCells];
97         _densityOld   = new float[_totalCells];
98         _heat         = new float[_totalCells];
99         _heatOld      = new float[_totalCells];
100         _obstacles    = new unsigned char[_totalCells]; // set 0 at end of step
101
102         // For threaded version:
103         _xVelocityTemp = new float[_totalCells];
104         _yVelocityTemp = new float[_totalCells];
105         _zVelocityTemp = new float[_totalCells];
106         _densityTemp   = new float[_totalCells];
107         _heatTemp      = new float[_totalCells];
108
109         // DG TODO: check if alloc went fine
110
111         for (int x = 0; x < _totalCells; x++)
112         {
113                 _density[x]      = 0.0f;
114                 _densityOld[x]   = 0.0f;
115                 _heat[x]         = 0.0f;
116                 _heatOld[x]      = 0.0f;
117                 _xVelocity[x]    = 0.0f;
118                 _yVelocity[x]    = 0.0f;
119                 _zVelocity[x]    = 0.0f;
120                 _xVelocityOld[x] = 0.0f;
121                 _yVelocityOld[x] = 0.0f;
122                 _zVelocityOld[x] = 0.0f;
123                 _xForce[x]       = 0.0f;
124                 _yForce[x]       = 0.0f;
125                 _zForce[x]       = 0.0f;
126                 _obstacles[x]    = false;
127         }
128
129         // set side obstacles
130         int index;
131         for (int y = 0; y < _yRes; y++)
132         for (int x = 0; x < _xRes; x++)
133         {
134                 // front slab
135                 index = x + y * _xRes;
136                 if(DOMAIN_BC_FRONT==1) _obstacles[index] = 1;
137
138                 // back slab
139                 index += _totalCells - _slabSize;
140                 if(DOMAIN_BC_BACK==1) _obstacles[index] = 1;
141         }
142
143         for (int z = 0; z < _zRes; z++)
144         for (int x = 0; x < _xRes; x++)
145         {
146                 // bottom slab
147                 index = x + z * _slabSize;
148                 if(DOMAIN_BC_BOTTOM==1) _obstacles[index] = 1;
149
150                 // top slab
151                 index += _slabSize - _xRes;
152                 if(DOMAIN_BC_TOP==1) _obstacles[index] = 1;
153         }
154
155         for (int z = 0; z < _zRes; z++)
156         for (int y = 0; y < _yRes; y++)
157         {
158                 // left slab
159                 index = y * _xRes + z * _slabSize;
160                 if(DOMAIN_BC_LEFT==1) _obstacles[index] = 1;
161
162                 // right slab
163                 index += _xRes - 1;
164                 if(DOMAIN_BC_RIGHT==1) _obstacles[index] = 1;
165         }
166 }
167
168 FLUID_3D::~FLUID_3D()
169 {
170         if (_xVelocity) delete[] _xVelocity;
171         if (_yVelocity) delete[] _yVelocity;
172         if (_zVelocity) delete[] _zVelocity;
173         if (_xVelocityOld) delete[] _xVelocityOld;
174         if (_yVelocityOld) delete[] _yVelocityOld;
175         if (_zVelocityOld) delete[] _zVelocityOld;
176         if (_xForce) delete[] _xForce;
177         if (_yForce) delete[] _yForce;
178         if (_zForce) delete[] _zForce;
179         if (_density) delete[] _density;
180         if (_densityOld) delete[] _densityOld;
181         if (_heat) delete[] _heat;
182         if (_heatOld) delete[] _heatOld;
183         if (_obstacles) delete[] _obstacles;
184     // if (_wTurbulence) delete _wTurbulence;
185
186         if (_xVelocityTemp) delete[] _xVelocityTemp;
187         if (_yVelocityTemp) delete[] _yVelocityTemp;
188         if (_zVelocityTemp) delete[] _zVelocityTemp;
189         if (_densityTemp) delete[] _densityTemp;
190         if (_heatTemp) delete[] _heatTemp;
191
192     // printf("deleted fluid\n");
193 }
194
195 // init direct access functions from blender
196 void FLUID_3D::initBlenderRNA(float *alpha, float *beta)
197 {
198         _alpha = alpha;
199         _beta = beta;
200 }
201
202 //////////////////////////////////////////////////////////////////////
203 // step simulation once
204 //////////////////////////////////////////////////////////////////////
205 void FLUID_3D::step()
206 {
207 #if PARALLEL==1
208         int threadval = 1;
209         threadval = omp_get_max_threads();
210
211         int stepParts = 1;
212         float partSize = _zRes;
213
214         stepParts = threadval*2;        // Dividing parallelized sections into numOfThreads * 2 sections
215         partSize = (float)_zRes/stepParts;      // Size of one part;
216
217   if (partSize < 4) {stepParts = threadval;                                     // If the slice gets too low (might actually slow things down, change it to larger
218                                         partSize = (float)_zRes/stepParts;}
219   if (partSize < 4) {stepParts = (int)(ceil((float)_zRes/4.0f));        // If it's still too low (only possible on future systems with +24 cores), change it to 4
220                                         partSize = (float)_zRes/stepParts;}
221 #else
222         int zBegin=0;
223         int zEnd=_zRes;
224 #endif
225
226
227 #if PARALLEL==1
228         #pragma omp parallel
229         {
230         #pragma omp for schedule(static,1)
231         for (int i=0; i<stepParts; i++)
232         {
233                 int zBegin = (int)((float)i*partSize + 0.5f);
234                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
235 #endif
236
237                 wipeBoundariesSL(zBegin, zEnd);
238                 addVorticity(zBegin, zEnd);
239                 addBuoyancy(_heat, _density, zBegin, zEnd);
240                 addForce(zBegin, zEnd);
241
242 #if PARALLEL==1
243         }       // end of parallel
244         #pragma omp barrier
245
246         #pragma omp single
247         {
248 #endif
249         SWAP_POINTERS(_xVelocity, _xVelocityTemp);
250         SWAP_POINTERS(_yVelocity, _yVelocityTemp);
251         SWAP_POINTERS(_zVelocity, _zVelocityTemp);
252 #if PARALLEL==1
253         }       // end of single
254
255         #pragma omp barrier
256
257         #pragma omp for
258         for (int i=0; i<2; i++)
259         {
260                 if (i==0)
261                 {
262 #endif
263                 project();
264 #if PARALLEL==1
265                 }
266                 else
267                 {
268 #endif
269                 diffuseHeat();
270 #if PARALLEL==1
271                 }
272         }
273
274         #pragma omp barrier
275
276         #pragma omp single
277         {
278 #endif
279         SWAP_POINTERS(_xVelocity, _xVelocityOld);
280         SWAP_POINTERS(_yVelocity, _yVelocityOld);
281         SWAP_POINTERS(_zVelocity, _zVelocityOld);
282         SWAP_POINTERS(_density, _densityOld);
283         SWAP_POINTERS(_heat, _heatOld);
284
285         advectMacCormackBegin(0, _zRes);
286
287 #if PARALLEL==1
288         }       // end of single
289
290         #pragma omp barrier
291
292
293         #pragma omp for schedule(static,1)
294         for (int i=0; i<stepParts; i++)
295         {
296
297                 int zBegin = (int)((float)i*partSize + 0.5f);
298                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
299 #endif
300
301         advectMacCormackEnd1(zBegin, zEnd);
302
303 #if PARALLEL==1
304         }       // end of parallel
305
306         #pragma omp barrier
307
308         #pragma omp for schedule(static,1)
309         for (int i=0; i<stepParts; i++)
310         {
311
312                 int zBegin = (int)((float)i*partSize + 0.5f);
313                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
314 #endif
315
316                 advectMacCormackEnd2(zBegin, zEnd);
317
318                 artificialDampingSL(zBegin, zEnd);
319
320                 // Using forces as temp arrays
321
322 #if PARALLEL==1
323         }
324         }
325
326
327
328         for (int i=1; i<stepParts; i++)
329         {
330                 int zPos=(int)((float)i*partSize + 0.5f);
331                 
332                 artificialDampingExactSL(zPos);
333
334         }
335 #endif
336
337         SWAP_POINTERS(_xVelocity, _xForce);
338         SWAP_POINTERS(_yVelocity, _yForce);
339         SWAP_POINTERS(_zVelocity, _zForce);
340
341
342
343
344         _totalTime += _dt;
345         _totalSteps++;          
346
347         for (int i = 0; i < _totalCells; i++)
348         {
349                 _xForce[i] = _yForce[i] = _zForce[i] = 0.0f;
350         }
351
352 }
353
354 //////////////////////////////////////////////////////////////////////
355 // helper function to dampen co-located grid artifacts of given arrays in intervals
356 // (only needed for velocity, strength (w) depends on testcase...
357 //////////////////////////////////////////////////////////////////////
358
359
360 void FLUID_3D::artificialDampingSL(int zBegin, int zEnd) {
361         const float w = 0.9;
362
363         memmove(_xForce+(_slabSize*zBegin), _xVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
364         memmove(_yForce+(_slabSize*zBegin), _yVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
365         memmove(_zForce+(_slabSize*zBegin), _zVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
366
367
368         if(_totalSteps % 4 == 1) {
369                 for (int z = zBegin+1; z < zEnd-1; z++)
370                         for (int y = 1; y < _res[1]-1; y++)
371                                 for (int x = 1+(y+z)%2; x < _res[0]-1; x+=2) {
372                                         const int index = x + y*_res[0] + z * _slabSize;
373                                         _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
374                                                         _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
375                                                         _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
376                                                         _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
377
378                                         _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
379                                                         _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
380                                                         _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
381                                                         _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
382
383                                         _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
384                                                         _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
385                                                         _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
386                                                         _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
387                                 }
388         }
389
390         if(_totalSteps % 4 == 3) {
391                 for (int z = zBegin+1; z < zEnd-1; z++)
392                         for (int y = 1; y < _res[1]-1; y++)
393                                 for (int x = 1+(y+z+1)%2; x < _res[0]-1; x+=2) {
394                                         const int index = x + y*_res[0] + z * _slabSize;
395                                         _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
396                                                         _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
397                                                         _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
398                                                         _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
399
400                                         _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
401                                                         _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
402                                                         _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
403                                                         _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
404
405                                         _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
406                                                         _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
407                                                         _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
408                                                         _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
409                                 }
410
411         }
412 }
413
414
415
416 void FLUID_3D::artificialDampingExactSL(int pos) {
417         const float w = 0.9;
418         int index, x,y,z;
419         
420
421         size_t posslab;
422
423         for (z=pos-1; z<=pos; z++)
424         {
425         posslab=z * _slabSize;
426
427         if(_totalSteps % 4 == 1) {
428                         for (y = 1; y < _res[1]-1; y++)
429                                 for (x = 1+(y+z)%2; x < _res[0]-1; x+=2) {
430                                         index = x + y*_res[0] + posslab;
431                                         _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
432                                                         _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
433                                                         _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
434                                                         _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
435
436                                         _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
437                                                         _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
438                                                         _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
439                                                         _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
440
441                                         _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
442                                                         _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
443                                                         _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
444                                                         _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
445                                         
446                                 }
447         }
448
449         if(_totalSteps % 4 == 3) {
450                         for (y = 1; y < _res[1]-1; y++)
451                                 for (x = 1+(y+z+1)%2; x < _res[0]-1; x+=2) {
452                                         index = x + y*_res[0] + posslab;
453                                         _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
454                                                         _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
455                                                         _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
456                                                         _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
457
458                                         _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
459                                                         _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
460                                                         _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
461                                                         _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
462
463                                         _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
464                                                         _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
465                                                         _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
466                                                         _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
467                                         
468                                 }
469
470         }
471         }
472 }
473
474 //////////////////////////////////////////////////////////////////////
475 // copy out the boundary in all directions
476 //////////////////////////////////////////////////////////////////////
477 void FLUID_3D::copyBorderAll(float* field, int zBegin, int zEnd)
478 {
479         int index, x, y, z;
480         int zSize = zEnd-zBegin;
481         int _blockTotalCells=_slabSize * zSize;
482
483         if ((zBegin==0))
484         for (int y = 0; y < _yRes; y++)
485                 for (int x = 0; x < _xRes; x++)
486                 {
487                         // front slab
488                         index = x + y * _xRes;
489                         field[index] = field[index + _slabSize];
490     }
491
492         if (zEnd==_zRes)
493         for (y = 0; y < _yRes; y++)
494                 for (x = 0; x < _xRes; x++)
495                 {
496
497                         // back slab
498                         index = x + y * _xRes + _blockTotalCells - _slabSize;
499                         field[index] = field[index - _slabSize];
500     }
501
502         for (z = 0; z < zSize; z++)
503                 for (x = 0; x < _xRes; x++)
504     {
505                         // bottom slab
506                         index = x + z * _slabSize;
507                         field[index] = field[index + _xRes];
508
509                         // top slab
510                         index += _slabSize - _xRes;
511                         field[index] = field[index - _xRes];
512     }
513
514         for (z = 0; z < zSize; z++)
515                 for (y = 0; y < _yRes; y++)
516     {
517                         // left slab
518                         index = y * _xRes + z * _slabSize;
519                         field[index] = field[index + 1];
520
521                         // right slab
522                         index += _xRes - 1;
523                         field[index] = field[index - 1];
524                 }
525 }
526
527 //////////////////////////////////////////////////////////////////////
528 // wipe boundaries of velocity and density
529 //////////////////////////////////////////////////////////////////////
530 void FLUID_3D::wipeBoundaries(int zBegin, int zEnd)
531 {
532         setZeroBorder(_xVelocity, _res, zBegin, zEnd);
533         setZeroBorder(_yVelocity, _res, zBegin, zEnd);
534         setZeroBorder(_zVelocity, _res, zBegin, zEnd);
535         setZeroBorder(_density, _res, zBegin, zEnd);
536 }
537
538 void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
539 {
540         
541         /////////////////////////////////////
542         // setZeroBorder to all:
543         /////////////////////////////////////
544
545         /////////////////////////////////////
546         // setZeroX
547         /////////////////////////////////////
548
549         const int slabSize = _xRes * _yRes;
550         int index, x,y,z;
551
552         for (z = zBegin; z < zEnd; z++)
553                 for (y = 0; y < _yRes; y++)
554                 {
555                         // left slab
556                         index = y * _xRes + z * slabSize;
557                         _xVelocity[index] = 0.0f;
558                         _yVelocity[index] = 0.0f;
559                         _zVelocity[index] = 0.0f;
560                         _density[index] = 0.0f;
561
562                         // right slab
563                         index += _xRes - 1;
564                         _xVelocity[index] = 0.0f;
565                         _yVelocity[index] = 0.0f;
566                         _zVelocity[index] = 0.0f;
567                         _density[index] = 0.0f;
568                 }
569
570         /////////////////////////////////////
571         // setZeroY
572         /////////////////////////////////////
573
574         for (z = zBegin; z < zEnd; z++)
575                 for (x = 0; x < _xRes; x++)
576                 {
577                         // bottom slab
578                         index = x + z * slabSize;
579                         _xVelocity[index] = 0.0f;
580                         _yVelocity[index] = 0.0f;
581                         _zVelocity[index] = 0.0f;
582                         _density[index] = 0.0f;
583
584                         // top slab
585                         index += slabSize - _xRes;
586                         _xVelocity[index] = 0.0f;
587                         _yVelocity[index] = 0.0f;
588                         _zVelocity[index] = 0.0f;
589                         _density[index] = 0.0f;
590
591                 }
592
593         /////////////////////////////////////
594         // setZeroZ
595         /////////////////////////////////////
596
597
598         const int totalCells = _xRes * _yRes * _zRes;
599
600         index = 0;
601         if ((zBegin == 0))
602         for (y = 0; y < _yRes; y++)
603                 for (x = 0; x < _xRes; x++, index++)
604                 {
605                         // front slab
606                         _xVelocity[index] = 0.0f;
607                         _yVelocity[index] = 0.0f;
608                         _zVelocity[index] = 0.0f;
609                         _density[index] = 0.0f;
610     }
611
612         if (zEnd == _zRes)
613         {
614                 index=0;
615                 int indexx=0;
616                 const int cellsslab = totalCells - slabSize;
617
618                 for (y = 0; y < _yRes; y++)
619                         for (x = 0; x < _xRes; x++, index++)
620                         {
621
622                                 // back slab
623                                 indexx = index + cellsslab;
624                                 _xVelocity[indexx] = 0.0f;
625                                 _yVelocity[indexx] = 0.0f;
626                                 _zVelocity[indexx] = 0.0f;
627                                 _density[indexx] = 0.0f;
628                         }
629         }
630
631 }
632 //////////////////////////////////////////////////////////////////////
633 // add forces to velocity field
634 //////////////////////////////////////////////////////////////////////
635 void FLUID_3D::addForce(int zBegin, int zEnd)
636 {
637         int begin=zBegin * _slabSize;
638         int end=begin + (zEnd - zBegin) * _slabSize;
639
640         for (int i = begin; i < end; i++)
641         {
642                 _xVelocityTemp[i] = _xVelocity[i] + _dt * _xForce[i];
643                 _yVelocityTemp[i] = _yVelocity[i] + _dt * _yForce[i];
644                 _zVelocityTemp[i] = _zVelocity[i] + _dt * _zForce[i];
645         }
646 }
647 //////////////////////////////////////////////////////////////////////
648 // project into divergence free field
649 //////////////////////////////////////////////////////////////////////
650 void FLUID_3D::project()
651 {
652         int x, y, z;
653         size_t index;
654
655         float *_pressure = new float[_totalCells];
656         float *_divergence   = new float[_totalCells];
657
658         memset(_pressure, 0, sizeof(float)*_totalCells);
659         memset(_divergence, 0, sizeof(float)*_totalCells);
660         
661         setObstacleBoundaries(_pressure, 0, _zRes);
662         
663         // copy out the boundaries
664         if(DOMAIN_BC_LEFT == 0)  setNeumannX(_xVelocity, _res, 0, _zRes);
665         else setZeroX(_xVelocity, _res, 0, _zRes); 
666
667         if(DOMAIN_BC_TOP == 0)   setNeumannY(_yVelocity, _res, 0, _zRes);
668         else setZeroY(_yVelocity, _res, 0, _zRes); 
669
670         if(DOMAIN_BC_FRONT == 0) setNeumannZ(_zVelocity, _res, 0, _zRes);
671         else setZeroZ(_zVelocity, _res, 0, _zRes);
672
673         // calculate divergence
674         index = _slabSize + _xRes + 1;
675         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
676                 for (y = 1; y < _yRes - 1; y++, index += 2)
677                         for (x = 1; x < _xRes - 1; x++, index++)
678                         {
679                                 float xright = _xVelocity[index + 1];
680                                 float xleft  = _xVelocity[index - 1];
681                                 float yup    = _yVelocity[index + _xRes];
682                                 float ydown  = _yVelocity[index - _xRes];
683                                 float ztop   = _zVelocity[index + _slabSize];
684                                 float zbottom = _zVelocity[index - _slabSize];
685
686                                 if(_obstacles[index+1]) xright = - _xVelocity[index];
687                                 if(_obstacles[index-1]) xleft  = - _xVelocity[index];
688                                 if(_obstacles[index+_xRes]) yup    = - _yVelocity[index];
689                                 if(_obstacles[index-_xRes]) ydown  = - _yVelocity[index];
690                                 if(_obstacles[index+_slabSize]) ztop    = - _zVelocity[index];
691                                 if(_obstacles[index-_slabSize]) zbottom = - _zVelocity[index];
692
693                                 _divergence[index] = -_dx * 0.5f * (
694                                                 xright - xleft +
695                                                 yup - ydown +
696                                                 ztop - zbottom );
697
698                                 // DG: commenting this helps CG to get a better start, 10-20% speed improvement
699                                 // _pressure[index] = 0.0f;
700                         }
701         copyBorderAll(_pressure, 0, _zRes);
702
703         // solve Poisson equation
704         solvePressurePre(_pressure, _divergence, _obstacles);
705
706         setObstaclePressure(_pressure, 0, _zRes);
707
708         // project out solution
709         float invDx = 1.0f / _dx;
710         index = _slabSize + _xRes + 1;
711         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
712                 for (y = 1; y < _yRes - 1; y++, index += 2)
713                         for (x = 1; x < _xRes - 1; x++, index++)
714                         {
715                                 if(!_obstacles[index])
716                                 {
717                                         _xVelocity[index] -= 0.5f * (_pressure[index + 1]     - _pressure[index - 1]) * invDx;
718                                         _yVelocity[index] -= 0.5f * (_pressure[index + _xRes]  - _pressure[index - _xRes]) * invDx;
719                                         _zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
720                                 }
721                         }
722
723         if (_pressure) delete[] _pressure;
724         if (_divergence) delete[] _divergence;
725 }
726
727
728
729
730 //////////////////////////////////////////////////////////////////////
731 // diffuse heat
732 //////////////////////////////////////////////////////////////////////
733 void FLUID_3D::diffuseHeat()
734 {
735         SWAP_POINTERS(_heat, _heatOld);
736
737         copyBorderAll(_heatOld, 0, _zRes);
738         solveHeat(_heat, _heatOld, _obstacles);
739
740         // zero out inside obstacles
741         for (int x = 0; x < _totalCells; x++)
742                 if (_obstacles[x])
743                         _heat[x] = 0.0f;
744
745 }
746
747 //////////////////////////////////////////////////////////////////////
748 // stamp an obstacle in the _obstacles field
749 //////////////////////////////////////////////////////////////////////
750 void FLUID_3D::addObstacle(OBSTACLE* obstacle)
751 {
752         int index = 0;
753         for (int z = 0; z < _zRes; z++)
754                 for (int y = 0; y < _yRes; y++)
755                         for (int x = 0; x < _xRes; x++, index++)
756                                 if (obstacle->inside(x * _dx, y * _dx, z * _dx)) {
757                                         _obstacles[index] = true;
758         }
759 }
760
761 //////////////////////////////////////////////////////////////////////
762 // calculate the obstacle directional types
763 //////////////////////////////////////////////////////////////////////
764 void FLUID_3D::setObstaclePressure(float *_pressure, int zBegin, int zEnd)
765 {
766
767         // compleately TODO
768
769         const size_t index_ = _slabSize + _xRes + 1;
770
771         //int vIndex=_slabSize + _xRes + 1;
772
773         int bb=0;
774         int bt=0;
775
776         if (zBegin == 0) {bb = 1;}
777         if (zEnd == _zRes) {bt = 1;}
778
779         // tag remaining obstacle blocks
780         for (int z = zBegin + bb; z < zEnd - bt; z++)
781         {
782                 size_t index = index_ +(z-1)*_slabSize;
783
784                 for (int y = 1; y < _yRes - 1; y++, index += 2)
785                 {
786                         for (int x = 1; x < _xRes - 1; x++, index++)
787                 {
788                         // could do cascade of ifs, but they are a pain
789                         if (_obstacles[index])
790                         {
791                                 const int top   = _obstacles[index + _slabSize];
792                                 const int bottom= _obstacles[index - _slabSize];
793                                 const int up    = _obstacles[index + _xRes];
794                                 const int down  = _obstacles[index - _xRes];
795                                 const int left  = _obstacles[index - 1];
796                                 const int right = _obstacles[index + 1];
797
798                                 // unused
799                                 // const bool fullz = (top && bottom);
800                                 // const bool fully = (up && down);
801                                 //const bool fullx = (left && right);
802
803                                 _xVelocity[index] =
804                                 _yVelocity[index] =
805                                 _zVelocity[index] = 0.0f;
806                                 _pressure[index] = 0.0f;
807
808                                 // average pressure neighbors
809                                 float pcnt = 0.;
810                                 if (left && !right) {
811                                         _pressure[index] += _pressure[index + 1];
812                                         pcnt += 1.;
813                                 }
814                                 if (!left && right) {
815                                         _pressure[index] += _pressure[index - 1];
816                                         pcnt += 1.;
817                                 }
818                                 if (up && !down) {
819                                         _pressure[index] += _pressure[index - _xRes];
820                                         pcnt += 1.;
821                                 }
822                                 if (!up && down) {
823                                         _pressure[index] += _pressure[index + _xRes];
824                                         pcnt += 1.;
825                                 }
826                                 if (top && !bottom) {
827                                         _pressure[index] += _pressure[index - _slabSize];
828                                         pcnt += 1.;
829                                 }
830                                 if (!top && bottom) {
831                                         _pressure[index] += _pressure[index + _slabSize];
832                                         pcnt += 1.;
833                                 }
834                                 
835                                 if(pcnt > 0.000001f)
836                                         _pressure[index] /= pcnt;
837
838                                 // TODO? set correct velocity bc's
839                                 // velocities are only set to zero right now
840                                 // this means it's not a full no-slip boundary condition
841                                 // but a "half-slip" - still looks ok right now
842                         }
843                         //vIndex++;
844                 }       // x loop
845                 //vIndex += 2;
846                 }       // y loop
847                 //vIndex += 2 * _xRes;
848         }       // z loop
849 }
850
851 void FLUID_3D::setObstacleBoundaries(float *_pressure, int zBegin, int zEnd)
852 {
853         // cull degenerate obstacles , move to addObstacle?
854
855         // r = b - Ax
856         const size_t index_ = _slabSize + _xRes + 1;
857
858         int bb=0;
859         int bt=0;
860
861         if (zBegin == 0) {bb = 1;}
862         if (zEnd == _zRes) {bt = 1;}
863
864         for (int z = zBegin + bb; z < zEnd - bt; z++)
865         {
866                 size_t index = index_ +(z-1)*_slabSize;
867                 
868                 for (int y = 1; y < _yRes - 1; y++, index += 2)
869                 {
870                         for (int x = 1; x < _xRes - 1; x++, index++)
871                         {
872                                 if (_obstacles[index] != EMPTY)
873                                 {
874                                         const int top   = _obstacles[index + _slabSize];
875                                         const int bottom= _obstacles[index - _slabSize];
876                                         const int up    = _obstacles[index + _xRes];
877                                         const int down  = _obstacles[index - _xRes];
878                                         const int left  = _obstacles[index - 1];
879                                         const int right = _obstacles[index + 1];
880
881                                         int counter = 0;
882                                         if (up)    counter++;
883                                         if (down)  counter++;
884                                         if (left)  counter++;
885                                         if (right) counter++;
886                                         if (top)  counter++;
887                                         if (bottom) counter++;
888
889                                         if (counter < 3)
890                                                 _obstacles[index] = EMPTY;
891                                 }
892                                 if (_obstacles[index])
893                                 {
894                                         _xVelocity[index] =
895                                         _yVelocity[index] =
896                                         _zVelocity[index] = 0.0f;
897                                         _pressure[index] = 0.0f;
898                                 }
899                                 //vIndex++;
900                         }       // x-loop
901                         //vIndex += 2;
902                 }       // y-loop
903                 //vIndex += 2* _xRes;
904         }       // z-loop
905 }
906
907 //////////////////////////////////////////////////////////////////////
908 // add buoyancy forces
909 //////////////////////////////////////////////////////////////////////
910 void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd)
911 {
912         int index = zBegin*_slabSize;
913
914         for (int z = zBegin; z < zEnd; z++)
915                 for (int y = 0; y < _yRes; y++)
916                         for (int x = 0; x < _xRes; x++, index++)
917                         {
918                                 _zForce[index] += *_alpha * density[index] + (*_beta * (heat[index] - _tempAmb)); // DG: was _yForce, changed for Blender
919                         }
920 }
921
922 //////////////////////////////////////////////////////////////////////
923 // add vorticity to the force field
924 //////////////////////////////////////////////////////////////////////
925 void FLUID_3D::addVorticity(int zBegin, int zEnd)
926 {
927         //int x,y,z,index;
928         if(_vorticityEps<=0.) return;
929
930         int _blockSize=zEnd-zBegin;
931         int _blockTotalCells = _slabSize * (_blockSize+2);
932
933         float *_xVorticity, *_yVorticity, *_zVorticity, *_vorticity;
934
935         int bb=0;
936         int bt=0;
937         int bb1=-1;
938         int bt1=-1;
939
940         if (zBegin == 0) {bb1 = 1; bb = 1; _blockTotalCells-=_blockSize;}
941         if (zEnd == _zRes) {bt1 = 1;bt = 1; _blockTotalCells-=_blockSize;}
942
943         _xVorticity = new float[_blockTotalCells];
944         _yVorticity = new float[_blockTotalCells];
945         _zVorticity = new float[_blockTotalCells];
946         _vorticity = new float[_blockTotalCells];
947
948         memset(_xVorticity, 0, sizeof(float)*_blockTotalCells);
949         memset(_yVorticity, 0, sizeof(float)*_blockTotalCells);
950         memset(_zVorticity, 0, sizeof(float)*_blockTotalCells);
951         memset(_vorticity, 0, sizeof(float)*_blockTotalCells);
952
953         //const size_t indexsetupV=_slabSize;
954         const size_t index_ = _slabSize + _xRes + 1;
955
956         // calculate vorticity
957         float gridSize = 0.5f / _dx;
958         //index = _slabSize + _xRes + 1;
959
960
961         size_t vIndex=_xRes + 1;
962
963         for (int z = zBegin + bb1; z < (zEnd - bt1); z++)
964         {
965                 size_t index = index_ +(z-1)*_slabSize;
966                 vIndex = index-(zBegin-1+bb)*_slabSize;
967
968                 for (int y = 1; y < _yRes - 1; y++, index += 2)
969                 {
970                         for (int x = 1; x < _xRes - 1; x++, index++)
971                         {
972
973                                 int up    = _obstacles[index + _xRes] ? index : index + _xRes;
974                                 int down  = _obstacles[index - _xRes] ? index : index - _xRes;
975                                 float dy  = (up == index || down == index) ? 1.0f / _dx : gridSize;
976                                 int out   = _obstacles[index + _slabSize] ? index : index + _slabSize;
977                                 int in    = _obstacles[index - _slabSize] ? index : index - _slabSize;
978                                 float dz  = (out == index || in == index) ? 1.0f / _dx : gridSize;
979                                 int right = _obstacles[index + 1] ? index : index + 1;
980                                 int left  = _obstacles[index - 1] ? index : index - 1;
981                                 float dx  = (right == index || right == index) ? 1.0f / _dx : gridSize;
982
983                                 _xVorticity[vIndex] = (_zVelocity[up] - _zVelocity[down]) * dy + (-_yVelocity[out] + _yVelocity[in]) * dz;
984                                 _yVorticity[vIndex] = (_xVelocity[out] - _xVelocity[in]) * dz + (-_zVelocity[right] + _zVelocity[left]) * dx;
985                                 _zVorticity[vIndex] = (_yVelocity[right] - _yVelocity[left]) * dx + (-_xVelocity[up] + _xVelocity[down])* dy;
986
987                                 _vorticity[vIndex] = sqrtf(_xVorticity[vIndex] * _xVorticity[vIndex] +
988                                                 _yVorticity[vIndex] * _yVorticity[vIndex] +
989                                                 _zVorticity[vIndex] * _zVorticity[vIndex]);
990
991                                 vIndex++;
992                         }
993                         vIndex+=2;
994                 }
995                 //vIndex+=2*_xRes;
996         }
997
998         // calculate normalized vorticity vectors
999         float eps = _vorticityEps;
1000         
1001         //index = _slabSize + _xRes + 1;
1002         vIndex=_slabSize + _xRes + 1;
1003
1004         for (int z = zBegin + bb; z < (zEnd - bt); z++)
1005         {
1006                 size_t index = index_ +(z-1)*_slabSize;
1007                 vIndex = index-(zBegin-1+bb)*_slabSize;
1008
1009                 for (int y = 1; y < _yRes - 1; y++, index += 2)
1010                 {
1011                         for (int x = 1; x < _xRes - 1; x++, index++)
1012                         {
1013                                 //
1014
1015                                 if (!_obstacles[index])
1016                                 {
1017                                         float N[3];
1018
1019                                         int up    = _obstacles[index + _xRes] ? vIndex : vIndex + _xRes;
1020                                         int down  = _obstacles[index - _xRes] ? vIndex : vIndex - _xRes;
1021                                         float dy  = (up == vIndex || down == vIndex) ? 1.0f / _dx : gridSize;
1022                                         int out   = _obstacles[index + _slabSize] ? vIndex : vIndex + _slabSize;
1023                                         int in    = _obstacles[index - _slabSize] ? vIndex : vIndex - _slabSize;
1024                                         float dz  = (out == vIndex || in == vIndex) ? 1.0f / _dx : gridSize;
1025                                         int right = _obstacles[index + 1] ? vIndex : vIndex + 1;
1026                                         int left  = _obstacles[index - 1] ? vIndex : vIndex - 1;
1027                                         float dx  = (right == vIndex || right == vIndex) ? 1.0f / _dx : gridSize;
1028                                         N[0] = (_vorticity[right] - _vorticity[left]) * dx;
1029                                         N[1] = (_vorticity[up] - _vorticity[down]) * dy;
1030                                         N[2] = (_vorticity[out] - _vorticity[in]) * dz;
1031
1032                                         float magnitude = sqrtf(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
1033                                         if (magnitude > 0.0f)
1034                                         {
1035                                                 magnitude = 1.0f / magnitude;
1036                                                 N[0] *= magnitude;
1037                                                 N[1] *= magnitude;
1038                                                 N[2] *= magnitude;
1039
1040                                                 _xForce[index] += (N[1] * _zVorticity[vIndex] - N[2] * _yVorticity[vIndex]) * _dx * eps;
1041                                                 _yForce[index] -= (N[0] * _zVorticity[vIndex] - N[2] * _xVorticity[vIndex]) * _dx * eps;
1042                                                 _zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * eps;
1043                                         }
1044                                         }       // if
1045                                         vIndex++;
1046                                         }       // x loop
1047                                 vIndex+=2;
1048                                 }               // y loop
1049                         //vIndex+=2*_xRes;
1050                 }                               // z loop
1051                                 
1052         if (_xVorticity) delete[] _xVorticity;
1053         if (_yVorticity) delete[] _yVorticity;
1054         if (_zVorticity) delete[] _zVorticity;
1055         if (_vorticity) delete[] _vorticity;
1056 }
1057
1058
1059 void FLUID_3D::advectMacCormackBegin(int zBegin, int zEnd)
1060 {
1061         Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
1062
1063         if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocityOld, res, zBegin, zEnd);
1064         else setZeroX(_xVelocityOld, res, zBegin, zEnd);
1065
1066         if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocityOld, res, zBegin, zEnd);
1067         else setZeroY(_yVelocityOld, res, zBegin, zEnd); 
1068
1069         if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocityOld, res, zBegin, zEnd);
1070         else setZeroZ(_zVelocityOld, res, zBegin, zEnd);
1071 }
1072
1073 //////////////////////////////////////////////////////////////////////
1074 // Advect using the MacCormack method from the Selle paper
1075 //////////////////////////////////////////////////////////////////////
1076 void FLUID_3D::advectMacCormackEnd1(int zBegin, int zEnd)
1077 {
1078         Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
1079
1080         const float dt0 = _dt / _dx;
1081
1082         int begin=zBegin * _slabSize;
1083         int end=begin + (zEnd - zBegin) * _slabSize;
1084         for (int x = begin; x < end; x++)
1085                 _xForce[x] = 0.0;
1086
1087         // advectFieldMacCormack1(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res)
1088
1089         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _densityTemp, res, zBegin, zEnd);
1090         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, res, zBegin, zEnd);
1091         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocity, res, zBegin, zEnd);
1092         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocity, res, zBegin, zEnd);
1093         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, res, zBegin, zEnd);
1094
1095         // Have to wait untill all the threads are done -> so continuing in step 3
1096 }
1097
1098 //////////////////////////////////////////////////////////////////////
1099 // Advect using the MacCormack method from the Selle paper
1100 //////////////////////////////////////////////////////////////////////
1101 void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd)
1102 {
1103         const float dt0 = _dt / _dx;
1104         Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
1105
1106         // use force array as temp arrays
1107         float* t1 = _xForce;
1108
1109         // advectFieldMacCormack2(dt, xVelocity, yVelocity, zVelocity, oldField, newField, tempfield, temp, res, obstacles)
1110
1111         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _density, _densityTemp, t1, res, _obstacles, zBegin, zEnd);
1112         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, t1, res, _obstacles, zBegin, zEnd);
1113         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocityTemp, _xVelocity, t1, res, _obstacles, zBegin, zEnd);
1114         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocityTemp, _yVelocity, t1, res, _obstacles, zBegin, zEnd);
1115         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocityTemp, _zVelocity, t1, res, _obstacles, zBegin, zEnd);
1116
1117         if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocityTemp, res, zBegin, zEnd);
1118         else setZeroX(_xVelocityTemp, res, zBegin, zEnd);                               
1119
1120         if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocityTemp, res, zBegin, zEnd);
1121         else setZeroY(_yVelocityTemp, res, zBegin, zEnd); 
1122
1123         if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd);
1124         else setZeroZ(_zVelocityTemp, res, zBegin, zEnd);
1125
1126         setZeroBorder(_density, res, zBegin, zEnd);
1127         setZeroBorder(_heat, res, zBegin, zEnd);
1128
1129
1130         /*int begin=zBegin * _slabSize;
1131         int end=begin + (zEnd - zBegin) * _slabSize;
1132   for (int x = begin; x < end; x++)
1133     _xForce[x] = _yForce[x] = 0.0f;*/
1134 }