25673630fc469d497243c26d5e1297f49e70d227
[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
208         int threadval = 1;
209 #if PARALLEL==1
210         threadval = omp_get_max_threads();
211 #endif
212
213         int stepParts = 1;
214         float partSize = _zRes;
215
216 #if PARALLEL==1
217         stepParts = threadval*2;        // Dividing parallelized sections into numOfThreads * 2 sections
218         partSize = (float)_zRes/stepParts;      // Size of one part;
219
220   if (partSize < 4) {stepParts = threadval;                                     // If the slice gets too low (might actually slow things down, change it to larger
221                                         partSize = (float)_zRes/stepParts;}
222   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
223                                         partSize = (float)_zRes/stepParts;}
224 #else
225         int zBegin=0;
226         int zEnd=_zRes;
227 #endif
228
229
230 #if PARALLEL==1
231         #pragma omp parallel
232         {
233         #pragma omp for schedule(static,1)
234         for (int i=0; i<stepParts; i++)
235         {
236                 int zBegin = (int)((float)i*partSize + 0.5f);
237                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
238 #endif
239
240                 wipeBoundariesSL(zBegin, zEnd);
241                 addVorticity(zBegin, zEnd);
242                 addBuoyancy(_heat, _density, zBegin, zEnd);
243                 addForce(zBegin, zEnd);
244
245 #if PARALLEL==1
246         }       // end of parallel
247         #pragma omp barrier
248
249         #pragma omp single
250         {
251 #endif
252         SWAP_POINTERS(_xVelocity, _xVelocityTemp);
253         SWAP_POINTERS(_yVelocity, _yVelocityTemp);
254         SWAP_POINTERS(_zVelocity, _zVelocityTemp);
255 #if PARALLEL==1
256         }       // end of single
257
258         #pragma omp barrier
259
260         #pragma omp for
261         for (int i=0; i<2; i++)
262         {
263                 if (i==0)
264                 {
265 #endif
266                 project();
267 #if PARALLEL==1
268                 }
269                 else
270                 {
271 #endif
272                 diffuseHeat();
273 #if PARALLEL==1
274                 }
275         }
276
277         #pragma omp barrier
278
279         #pragma omp single
280         {
281 #endif
282         SWAP_POINTERS(_xVelocity, _xVelocityOld);
283         SWAP_POINTERS(_yVelocity, _yVelocityOld);
284         SWAP_POINTERS(_zVelocity, _zVelocityOld);
285         SWAP_POINTERS(_density, _densityOld);
286         SWAP_POINTERS(_heat, _heatOld);
287
288         advectMacCormackBegin(0, _zRes);
289
290 #if PARALLEL==1
291         }       // end of single
292
293         #pragma omp barrier
294
295
296         #pragma omp for schedule(static,1)
297         for (int i=0; i<stepParts; i++)
298         {
299
300                 int zBegin = (int)((float)i*partSize + 0.5f);
301                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
302 #endif
303
304         advectMacCormackEnd1(zBegin, zEnd);
305
306 #if PARALLEL==1
307         }       // end of parallel
308
309         #pragma omp barrier
310
311         #pragma omp for schedule(static,1)
312         for (int i=0; i<stepParts; i++)
313         {
314
315                 int zBegin = (int)((float)i*partSize + 0.5f);
316                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
317 #endif
318
319                 advectMacCormackEnd2(zBegin, zEnd);
320
321                 artificialDampingSL(zBegin, zEnd);
322
323                 // Using forces as temp arrays
324
325 #if PARALLEL==1
326         }
327         }
328
329
330
331         for (int i=1; i<stepParts; i++)
332         {
333                 int zPos=(int)((float)i*partSize + 0.5f);
334                 
335                 artificialDampingExactSL(zPos);
336
337         }
338 #endif
339
340         SWAP_POINTERS(_xVelocity, _xForce);
341         SWAP_POINTERS(_yVelocity, _yForce);
342         SWAP_POINTERS(_zVelocity, _zForce);
343
344
345
346
347         _totalTime += _dt;
348         _totalSteps++;          
349
350         for (int i = 0; i < _totalCells; i++)
351         {
352                 _xForce[i] = _yForce[i] = _zForce[i] = 0.0f;
353         }
354
355 }
356
357 //////////////////////////////////////////////////////////////////////
358 // helper function to dampen co-located grid artifacts of given arrays in intervals
359 // (only needed for velocity, strength (w) depends on testcase...
360 //////////////////////////////////////////////////////////////////////
361
362
363 void FLUID_3D::artificialDampingSL(int zBegin, int zEnd) {
364         const float w = 0.9;
365
366         memmove(_xForce+(_slabSize*zBegin), _xVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
367         memmove(_yForce+(_slabSize*zBegin), _yVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
368         memmove(_zForce+(_slabSize*zBegin), _zVelocityTemp+(_slabSize*zBegin), sizeof(float)*_slabSize*(zEnd-zBegin));
369
370
371         if(_totalSteps % 4 == 1) {
372                 for (int z = zBegin+1; z < zEnd-1; z++)
373                         for (int y = 1; y < _res[1]-1; y++)
374                                 for (int x = 1+(y+z)%2; x < _res[0]-1; x+=2) {
375                                         const int index = x + y*_res[0] + z * _slabSize;
376                                         _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
377                                                         _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
378                                                         _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
379                                                         _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
380
381                                         _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
382                                                         _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
383                                                         _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
384                                                         _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
385
386                                         _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
387                                                         _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
388                                                         _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
389                                                         _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
390                                 }
391         }
392
393         if(_totalSteps % 4 == 3) {
394                 for (int z = zBegin+1; z < zEnd-1; z++)
395                         for (int y = 1; y < _res[1]-1; y++)
396                                 for (int x = 1+(y+z+1)%2; x < _res[0]-1; x+=2) {
397                                         const int index = x + y*_res[0] + z * _slabSize;
398                                         _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
399                                                         _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
400                                                         _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
401                                                         _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
402
403                                         _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
404                                                         _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
405                                                         _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
406                                                         _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
407
408                                         _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
409                                                         _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
410                                                         _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
411                                                         _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
412                                 }
413
414         }
415 }
416
417
418
419 void FLUID_3D::artificialDampingExactSL(int pos) {
420         const float w = 0.9;
421         int index, x,y,z;
422         
423
424         size_t posslab;
425
426         for (z=pos-1; z<=pos; z++)
427         {
428         posslab=z * _slabSize;
429
430         if(_totalSteps % 4 == 1) {
431                         for (y = 1; y < _res[1]-1; y++)
432                                 for (x = 1+(y+z)%2; x < _res[0]-1; x+=2) {
433                                         index = x + y*_res[0] + posslab;
434                                         _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
435                                                         _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
436                                                         _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
437                                                         _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
438
439                                         _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
440                                                         _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
441                                                         _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
442                                                         _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
443
444                                         _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
445                                                         _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
446                                                         _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
447                                                         _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
448                                         
449                                 }
450         }
451
452         if(_totalSteps % 4 == 3) {
453                         for (y = 1; y < _res[1]-1; y++)
454                                 for (x = 1+(y+z+1)%2; x < _res[0]-1; x+=2) {
455                                         index = x + y*_res[0] + posslab;
456                                         _xForce[index] = (1-w)*_xVelocityTemp[index] + 1./6. * w*(
457                                                         _xVelocityTemp[index+1] + _xVelocityTemp[index-1] +
458                                                         _xVelocityTemp[index+_res[0]] + _xVelocityTemp[index-_res[0]] +
459                                                         _xVelocityTemp[index+_slabSize] + _xVelocityTemp[index-_slabSize] );
460
461                                         _yForce[index] = (1-w)*_yVelocityTemp[index] + 1./6. * w*(
462                                                         _yVelocityTemp[index+1] + _yVelocityTemp[index-1] +
463                                                         _yVelocityTemp[index+_res[0]] + _yVelocityTemp[index-_res[0]] +
464                                                         _yVelocityTemp[index+_slabSize] + _yVelocityTemp[index-_slabSize] );
465
466                                         _zForce[index] = (1-w)*_zVelocityTemp[index] + 1./6. * w*(
467                                                         _zVelocityTemp[index+1] + _zVelocityTemp[index-1] +
468                                                         _zVelocityTemp[index+_res[0]] + _zVelocityTemp[index-_res[0]] +
469                                                         _zVelocityTemp[index+_slabSize] + _zVelocityTemp[index-_slabSize] );
470                                         
471                                 }
472
473         }
474         }
475 }
476
477 //////////////////////////////////////////////////////////////////////
478 // copy out the boundary in all directions
479 //////////////////////////////////////////////////////////////////////
480 void FLUID_3D::copyBorderAll(float* field, int zBegin, int zEnd)
481 {
482         int index, x, y, z;
483         int zSize = zEnd-zBegin;
484         int _blockTotalCells=_slabSize * zSize;
485
486         if ((zBegin==0))
487         for (int y = 0; y < _yRes; y++)
488                 for (int x = 0; x < _xRes; x++)
489                 {
490                         // front slab
491                         index = x + y * _xRes;
492                         field[index] = field[index + _slabSize];
493     }
494
495         if (zEnd==_zRes)
496         for (y = 0; y < _yRes; y++)
497                 for (x = 0; x < _xRes; x++)
498                 {
499
500                         // back slab
501                         index = x + y * _xRes + _blockTotalCells - _slabSize;
502                         field[index] = field[index - _slabSize];
503     }
504
505         for (z = 0; z < zSize; z++)
506                 for (x = 0; x < _xRes; x++)
507     {
508                         // bottom slab
509                         index = x + z * _slabSize;
510                         field[index] = field[index + _xRes];
511
512                         // top slab
513                         index += _slabSize - _xRes;
514                         field[index] = field[index - _xRes];
515     }
516
517         for (z = 0; z < zSize; z++)
518                 for (y = 0; y < _yRes; y++)
519     {
520                         // left slab
521                         index = y * _xRes + z * _slabSize;
522                         field[index] = field[index + 1];
523
524                         // right slab
525                         index += _xRes - 1;
526                         field[index] = field[index - 1];
527                 }
528 }
529
530 //////////////////////////////////////////////////////////////////////
531 // wipe boundaries of velocity and density
532 //////////////////////////////////////////////////////////////////////
533 void FLUID_3D::wipeBoundaries(int zBegin, int zEnd)
534 {
535         setZeroBorder(_xVelocity, _res, zBegin, zEnd);
536         setZeroBorder(_yVelocity, _res, zBegin, zEnd);
537         setZeroBorder(_zVelocity, _res, zBegin, zEnd);
538         setZeroBorder(_density, _res, zBegin, zEnd);
539 }
540
541 void FLUID_3D::wipeBoundariesSL(int zBegin, int zEnd)
542 {
543         
544         /////////////////////////////////////
545         // setZeroBorder to all:
546         /////////////////////////////////////
547
548         /////////////////////////////////////
549         // setZeroX
550         /////////////////////////////////////
551
552         const int slabSize = _xRes * _yRes;
553         int index, x,y,z;
554
555         for (z = zBegin; z < zEnd; z++)
556                 for (y = 0; y < _yRes; y++)
557                 {
558                         // left slab
559                         index = y * _xRes + z * slabSize;
560                         _xVelocity[index] = 0.0f;
561                         _yVelocity[index] = 0.0f;
562                         _zVelocity[index] = 0.0f;
563                         _density[index] = 0.0f;
564
565                         // right slab
566                         index += _xRes - 1;
567                         _xVelocity[index] = 0.0f;
568                         _yVelocity[index] = 0.0f;
569                         _zVelocity[index] = 0.0f;
570                         _density[index] = 0.0f;
571                 }
572
573         /////////////////////////////////////
574         // setZeroY
575         /////////////////////////////////////
576
577         for (z = zBegin; z < zEnd; z++)
578                 for (x = 0; x < _xRes; x++)
579                 {
580                         // bottom slab
581                         index = x + z * slabSize;
582                         _xVelocity[index] = 0.0f;
583                         _yVelocity[index] = 0.0f;
584                         _zVelocity[index] = 0.0f;
585                         _density[index] = 0.0f;
586
587                         // top slab
588                         index += slabSize - _xRes;
589                         _xVelocity[index] = 0.0f;
590                         _yVelocity[index] = 0.0f;
591                         _zVelocity[index] = 0.0f;
592                         _density[index] = 0.0f;
593
594                 }
595
596         /////////////////////////////////////
597         // setZeroZ
598         /////////////////////////////////////
599
600
601         const int totalCells = _xRes * _yRes * _zRes;
602
603         index = 0;
604         if ((zBegin == 0))
605         for (y = 0; y < _yRes; y++)
606                 for (x = 0; x < _xRes; x++, index++)
607                 {
608                         // front slab
609                         _xVelocity[index] = 0.0f;
610                         _yVelocity[index] = 0.0f;
611                         _zVelocity[index] = 0.0f;
612                         _density[index] = 0.0f;
613     }
614
615         if (zEnd == _zRes)
616         {
617                 index=0;
618                 int indexx=0;
619                 const int cellsslab = totalCells - slabSize;
620
621                 for (y = 0; y < _yRes; y++)
622                         for (x = 0; x < _xRes; x++, index++)
623                         {
624
625                                 // back slab
626                                 indexx = index + cellsslab;
627                                 _xVelocity[indexx] = 0.0f;
628                                 _yVelocity[indexx] = 0.0f;
629                                 _zVelocity[indexx] = 0.0f;
630                                 _density[indexx] = 0.0f;
631                         }
632         }
633
634 }
635 //////////////////////////////////////////////////////////////////////
636 // add forces to velocity field
637 //////////////////////////////////////////////////////////////////////
638 void FLUID_3D::addForce(int zBegin, int zEnd)
639 {
640         int begin=zBegin * _slabSize;
641         int end=begin + (zEnd - zBegin) * _slabSize;
642
643         for (int i = begin; i < end; i++)
644         {
645                 _xVelocityTemp[i] = _xVelocity[i] + _dt * _xForce[i];
646                 _yVelocityTemp[i] = _yVelocity[i] + _dt * _yForce[i];
647                 _zVelocityTemp[i] = _zVelocity[i] + _dt * _zForce[i];
648         }
649 }
650 //////////////////////////////////////////////////////////////////////
651 // project into divergence free field
652 //////////////////////////////////////////////////////////////////////
653 void FLUID_3D::project()
654 {
655         int x, y, z;
656         size_t index;
657
658         float *_pressure = new float[_totalCells];
659         float *_divergence   = new float[_totalCells];
660
661         memset(_pressure, 0, sizeof(float)*_totalCells);
662         memset(_divergence, 0, sizeof(float)*_totalCells);
663         
664         setObstacleBoundaries(_pressure, 0, _zRes);
665         
666         // copy out the boundaries
667         if(DOMAIN_BC_LEFT == 0)  setNeumannX(_xVelocity, _res, 0, _zRes);
668         else setZeroX(_xVelocity, _res, 0, _zRes); 
669
670         if(DOMAIN_BC_TOP == 0)   setNeumannY(_yVelocity, _res, 0, _zRes);
671         else setZeroY(_yVelocity, _res, 0, _zRes); 
672
673         if(DOMAIN_BC_FRONT == 0) setNeumannZ(_zVelocity, _res, 0, _zRes);
674         else setZeroZ(_zVelocity, _res, 0, _zRes);
675
676         // calculate divergence
677         index = _slabSize + _xRes + 1;
678         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
679                 for (y = 1; y < _yRes - 1; y++, index += 2)
680                         for (x = 1; x < _xRes - 1; x++, index++)
681                         {
682                                 float xright = _xVelocity[index + 1];
683                                 float xleft  = _xVelocity[index - 1];
684                                 float yup    = _yVelocity[index + _xRes];
685                                 float ydown  = _yVelocity[index - _xRes];
686                                 float ztop   = _zVelocity[index + _slabSize];
687                                 float zbottom = _zVelocity[index - _slabSize];
688
689                                 if(_obstacles[index+1]) xright = - _xVelocity[index];
690                                 if(_obstacles[index-1]) xleft  = - _xVelocity[index];
691                                 if(_obstacles[index+_xRes]) yup    = - _yVelocity[index];
692                                 if(_obstacles[index-_xRes]) ydown  = - _yVelocity[index];
693                                 if(_obstacles[index+_slabSize]) ztop    = - _zVelocity[index];
694                                 if(_obstacles[index-_slabSize]) zbottom = - _zVelocity[index];
695
696                                 _divergence[index] = -_dx * 0.5f * (
697                                                 xright - xleft +
698                                                 yup - ydown +
699                                                 ztop - zbottom );
700
701                                 // DG: commenting this helps CG to get a better start, 10-20% speed improvement
702                                 // _pressure[index] = 0.0f;
703                         }
704         copyBorderAll(_pressure, 0, _zRes);
705
706         // solve Poisson equation
707         solvePressurePre(_pressure, _divergence, _obstacles);
708
709         setObstaclePressure(_pressure, 0, _zRes);
710
711         // project out solution
712         float invDx = 1.0f / _dx;
713         index = _slabSize + _xRes + 1;
714         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
715                 for (y = 1; y < _yRes - 1; y++, index += 2)
716                         for (x = 1; x < _xRes - 1; x++, index++)
717                         {
718                                 if(!_obstacles[index])
719                                 {
720                                         _xVelocity[index] -= 0.5f * (_pressure[index + 1]     - _pressure[index - 1]) * invDx;
721                                         _yVelocity[index] -= 0.5f * (_pressure[index + _xRes]  - _pressure[index - _xRes]) * invDx;
722                                         _zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
723                                 }
724                         }
725
726         if (_pressure) delete[] _pressure;
727         if (_divergence) delete[] _divergence;
728 }
729
730
731
732
733 //////////////////////////////////////////////////////////////////////
734 // diffuse heat
735 //////////////////////////////////////////////////////////////////////
736 void FLUID_3D::diffuseHeat()
737 {
738         SWAP_POINTERS(_heat, _heatOld);
739
740         copyBorderAll(_heatOld, 0, _zRes);
741         solveHeat(_heat, _heatOld, _obstacles);
742
743         // zero out inside obstacles
744         for (int x = 0; x < _totalCells; x++)
745                 if (_obstacles[x])
746                         _heat[x] = 0.0f;
747
748 }
749
750 //////////////////////////////////////////////////////////////////////
751 // stamp an obstacle in the _obstacles field
752 //////////////////////////////////////////////////////////////////////
753 void FLUID_3D::addObstacle(OBSTACLE* obstacle)
754 {
755         int index = 0;
756         for (int z = 0; z < _zRes; z++)
757                 for (int y = 0; y < _yRes; y++)
758                         for (int x = 0; x < _xRes; x++, index++)
759                                 if (obstacle->inside(x * _dx, y * _dx, z * _dx)) {
760                                         _obstacles[index] = true;
761         }
762 }
763
764 //////////////////////////////////////////////////////////////////////
765 // calculate the obstacle directional types
766 //////////////////////////////////////////////////////////////////////
767 void FLUID_3D::setObstaclePressure(float *_pressure, int zBegin, int zEnd)
768 {
769
770         // compleately TODO
771
772         const size_t index_ = _slabSize + _xRes + 1;
773
774         //int vIndex=_slabSize + _xRes + 1;
775
776         int bb=0;
777         int bt=0;
778
779         if (zBegin == 0) {bb = 1;}
780         if (zEnd == _zRes) {bt = 1;}
781
782         // tag remaining obstacle blocks
783         for (int z = zBegin + bb; z < zEnd - bt; z++)
784         {
785                 size_t index = index_ +(z-1)*_slabSize;
786
787                 for (int y = 1; y < _yRes - 1; y++, index += 2)
788                 {
789                         for (int x = 1; x < _xRes - 1; x++, index++)
790                 {
791                         // could do cascade of ifs, but they are a pain
792                         if (_obstacles[index])
793                         {
794                                 const int top   = _obstacles[index + _slabSize];
795                                 const int bottom= _obstacles[index - _slabSize];
796                                 const int up    = _obstacles[index + _xRes];
797                                 const int down  = _obstacles[index - _xRes];
798                                 const int left  = _obstacles[index - 1];
799                                 const int right = _obstacles[index + 1];
800
801                                 // unused
802                                 // const bool fullz = (top && bottom);
803                                 // const bool fully = (up && down);
804                                 //const bool fullx = (left && right);
805
806                                 _xVelocity[index] =
807                                 _yVelocity[index] =
808                                 _zVelocity[index] = 0.0f;
809                                 _pressure[index] = 0.0f;
810
811                                 // average pressure neighbors
812                                 float pcnt = 0.;
813                                 if (left && !right) {
814                                         _pressure[index] += _pressure[index + 1];
815                                         pcnt += 1.;
816                                 }
817                                 if (!left && right) {
818                                         _pressure[index] += _pressure[index - 1];
819                                         pcnt += 1.;
820                                 }
821                                 if (up && !down) {
822                                         _pressure[index] += _pressure[index - _xRes];
823                                         pcnt += 1.;
824                                 }
825                                 if (!up && down) {
826                                         _pressure[index] += _pressure[index + _xRes];
827                                         pcnt += 1.;
828                                 }
829                                 if (top && !bottom) {
830                                         _pressure[index] += _pressure[index - _slabSize];
831                                         pcnt += 1.;
832                                 }
833                                 if (!top && bottom) {
834                                         _pressure[index] += _pressure[index + _slabSize];
835                                         pcnt += 1.;
836                                 }
837                                 
838                                 if(pcnt > 0.000001f)
839                                         _pressure[index] /= pcnt;
840
841                                 // TODO? set correct velocity bc's
842                                 // velocities are only set to zero right now
843                                 // this means it's not a full no-slip boundary condition
844                                 // but a "half-slip" - still looks ok right now
845                         }
846                         //vIndex++;
847                 }       // x loop
848                 //vIndex += 2;
849                 }       // y loop
850                 //vIndex += 2 * _xRes;
851         }       // z loop
852 }
853
854 void FLUID_3D::setObstacleBoundaries(float *_pressure, int zBegin, int zEnd)
855 {
856         // cull degenerate obstacles , move to addObstacle?
857
858         // r = b - Ax
859         const size_t index_ = _slabSize + _xRes + 1;
860
861         int bb=0;
862         int bt=0;
863
864         if (zBegin == 0) {bb = 1;}
865         if (zEnd == _zRes) {bt = 1;}
866
867         for (int z = zBegin + bb; z < zEnd - bt; z++)
868         {
869                 size_t index = index_ +(z-1)*_slabSize;
870                 
871                 for (int y = 1; y < _yRes - 1; y++, index += 2)
872                 {
873                         for (int x = 1; x < _xRes - 1; x++, index++)
874                         {
875                                 if (_obstacles[index] != EMPTY)
876                                 {
877                                         const int top   = _obstacles[index + _slabSize];
878                                         const int bottom= _obstacles[index - _slabSize];
879                                         const int up    = _obstacles[index + _xRes];
880                                         const int down  = _obstacles[index - _xRes];
881                                         const int left  = _obstacles[index - 1];
882                                         const int right = _obstacles[index + 1];
883
884                                         int counter = 0;
885                                         if (up)    counter++;
886                                         if (down)  counter++;
887                                         if (left)  counter++;
888                                         if (right) counter++;
889                                         if (top)  counter++;
890                                         if (bottom) counter++;
891
892                                         if (counter < 3)
893                                                 _obstacles[index] = EMPTY;
894                                 }
895                                 if (_obstacles[index])
896                                 {
897                                         _xVelocity[index] =
898                                         _yVelocity[index] =
899                                         _zVelocity[index] = 0.0f;
900                                         _pressure[index] = 0.0f;
901                                 }
902                                 //vIndex++;
903                         }       // x-loop
904                         //vIndex += 2;
905                 }       // y-loop
906                 //vIndex += 2* _xRes;
907         }       // z-loop
908 }
909
910 //////////////////////////////////////////////////////////////////////
911 // add buoyancy forces
912 //////////////////////////////////////////////////////////////////////
913 void FLUID_3D::addBuoyancy(float *heat, float *density, int zBegin, int zEnd)
914 {
915         int index = zBegin*_slabSize;
916
917         for (int z = zBegin; z < zEnd; z++)
918                 for (int y = 0; y < _yRes; y++)
919                         for (int x = 0; x < _xRes; x++, index++)
920                         {
921                                 _zForce[index] += *_alpha * density[index] + (*_beta * (heat[index] - _tempAmb)); // DG: was _yForce, changed for Blender
922                         }
923 }
924
925 //////////////////////////////////////////////////////////////////////
926 // add vorticity to the force field
927 //////////////////////////////////////////////////////////////////////
928 void FLUID_3D::addVorticity(int zBegin, int zEnd)
929 {
930         //int x,y,z,index;
931         if(_vorticityEps<=0.) return;
932
933         int _blockSize=zEnd-zBegin;
934         int _blockTotalCells = _slabSize * (_blockSize+2);
935
936         float *_xVorticity, *_yVorticity, *_zVorticity, *_vorticity;
937
938         int _vIndex = _slabSize + _xRes + 1;
939         int bb=0;
940         int bt=0;
941         int bb1=-1;
942         int bt1=-1;
943
944         if (zBegin == 0) {bb1 = 1; bb = 1; _blockTotalCells-=_blockSize;}
945         if (zEnd == _zRes) {bt1 = 1;bt = 1; _blockTotalCells-=_blockSize;}
946
947         _xVorticity = new float[_blockTotalCells];
948         _yVorticity = new float[_blockTotalCells];
949         _zVorticity = new float[_blockTotalCells];
950         _vorticity = new float[_blockTotalCells];
951
952         memset(_xVorticity, 0, sizeof(float)*_blockTotalCells);
953         memset(_yVorticity, 0, sizeof(float)*_blockTotalCells);
954         memset(_zVorticity, 0, sizeof(float)*_blockTotalCells);
955         memset(_vorticity, 0, sizeof(float)*_blockTotalCells);
956
957         //const size_t indexsetupV=_slabSize;
958         const size_t index_ = _slabSize + _xRes + 1;
959
960         // calculate vorticity
961         float gridSize = 0.5f / _dx;
962         //index = _slabSize + _xRes + 1;
963
964
965         size_t vIndex=_xRes + 1;
966
967         for (int z = zBegin + bb1; z < (zEnd - bt1); z++)
968         {
969                 size_t index = index_ +(z-1)*_slabSize;
970                 vIndex = index-(zBegin-1+bb)*_slabSize;
971
972                 for (int y = 1; y < _yRes - 1; y++, index += 2)
973                 {
974                         for (int x = 1; x < _xRes - 1; x++, index++)
975                         {
976
977                                 int up    = _obstacles[index + _xRes] ? index : index + _xRes;
978                                 int down  = _obstacles[index - _xRes] ? index : index - _xRes;
979                                 float dy  = (up == index || down == index) ? 1.0f / _dx : gridSize;
980                                 int out   = _obstacles[index + _slabSize] ? index : index + _slabSize;
981                                 int in    = _obstacles[index - _slabSize] ? index : index - _slabSize;
982                                 float dz  = (out == index || in == index) ? 1.0f / _dx : gridSize;
983                                 int right = _obstacles[index + 1] ? index : index + 1;
984                                 int left  = _obstacles[index - 1] ? index : index - 1;
985                                 float dx  = (right == index || right == index) ? 1.0f / _dx : gridSize;
986
987                                 _xVorticity[vIndex] = (_zVelocity[up] - _zVelocity[down]) * dy + (-_yVelocity[out] + _yVelocity[in]) * dz;
988                                 _yVorticity[vIndex] = (_xVelocity[out] - _xVelocity[in]) * dz + (-_zVelocity[right] + _zVelocity[left]) * dx;
989                                 _zVorticity[vIndex] = (_yVelocity[right] - _yVelocity[left]) * dx + (-_xVelocity[up] + _xVelocity[down])* dy;
990
991                                 _vorticity[vIndex] = sqrtf(_xVorticity[vIndex] * _xVorticity[vIndex] +
992                                                 _yVorticity[vIndex] * _yVorticity[vIndex] +
993                                                 _zVorticity[vIndex] * _zVorticity[vIndex]);
994
995                                 vIndex++;
996                         }
997                         vIndex+=2;
998                 }
999                 //vIndex+=2*_xRes;
1000         }
1001
1002         // calculate normalized vorticity vectors
1003         float eps = _vorticityEps;
1004         
1005         //index = _slabSize + _xRes + 1;
1006         vIndex=_slabSize + _xRes + 1;
1007
1008         for (int z = zBegin + bb; z < (zEnd - bt); z++)
1009         {
1010                 size_t index = index_ +(z-1)*_slabSize;
1011                 vIndex = index-(zBegin-1+bb)*_slabSize;
1012
1013                 for (int y = 1; y < _yRes - 1; y++, index += 2)
1014                 {
1015                         for (int x = 1; x < _xRes - 1; x++, index++)
1016                         {
1017                                 //
1018
1019                                 if (!_obstacles[index])
1020                                 {
1021                                         float N[3];
1022
1023                                         int up    = _obstacles[index + _xRes] ? vIndex : vIndex + _xRes;
1024                                         int down  = _obstacles[index - _xRes] ? vIndex : vIndex - _xRes;
1025                                         float dy  = (up == vIndex || down == vIndex) ? 1.0f / _dx : gridSize;
1026                                         int out   = _obstacles[index + _slabSize] ? vIndex : vIndex + _slabSize;
1027                                         int in    = _obstacles[index - _slabSize] ? vIndex : vIndex - _slabSize;
1028                                         float dz  = (out == vIndex || in == vIndex) ? 1.0f / _dx : gridSize;
1029                                         int right = _obstacles[index + 1] ? vIndex : vIndex + 1;
1030                                         int left  = _obstacles[index - 1] ? vIndex : vIndex - 1;
1031                                         float dx  = (right == vIndex || right == vIndex) ? 1.0f / _dx : gridSize;
1032                                         N[0] = (_vorticity[right] - _vorticity[left]) * dx;
1033                                         N[1] = (_vorticity[up] - _vorticity[down]) * dy;
1034                                         N[2] = (_vorticity[out] - _vorticity[in]) * dz;
1035
1036                                         float magnitude = sqrtf(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
1037                                         if (magnitude > 0.0f)
1038                                         {
1039                                                 magnitude = 1.0f / magnitude;
1040                                                 N[0] *= magnitude;
1041                                                 N[1] *= magnitude;
1042                                                 N[2] *= magnitude;
1043
1044                                                 _xForce[index] += (N[1] * _zVorticity[vIndex] - N[2] * _yVorticity[vIndex]) * _dx * eps;
1045                                                 _yForce[index] -= (N[0] * _zVorticity[vIndex] - N[2] * _xVorticity[vIndex]) * _dx * eps;
1046                                                 _zForce[index] += (N[0] * _yVorticity[vIndex] - N[1] * _xVorticity[vIndex]) * _dx * eps;
1047                                         }
1048                                         }       // if
1049                                         vIndex++;
1050                                         }       // x loop
1051                                 vIndex+=2;
1052                                 }               // y loop
1053                         //vIndex+=2*_xRes;
1054                 }                               // z loop
1055                                 
1056         if (_xVorticity) delete[] _xVorticity;
1057         if (_yVorticity) delete[] _yVorticity;
1058         if (_zVorticity) delete[] _zVorticity;
1059         if (_vorticity) delete[] _vorticity;
1060 }
1061
1062
1063 void FLUID_3D::advectMacCormackBegin(int zBegin, int zEnd)
1064 {
1065         Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
1066
1067         if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocityOld, res, zBegin, zEnd);
1068         else setZeroX(_xVelocityOld, res, zBegin, zEnd);
1069
1070         if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocityOld, res, zBegin, zEnd);
1071         else setZeroY(_yVelocityOld, res, zBegin, zEnd); 
1072
1073         if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocityOld, res, zBegin, zEnd);
1074         else setZeroZ(_zVelocityOld, res, zBegin, zEnd);
1075 }
1076
1077 //////////////////////////////////////////////////////////////////////
1078 // Advect using the MacCormack method from the Selle paper
1079 //////////////////////////////////////////////////////////////////////
1080 void FLUID_3D::advectMacCormackEnd1(int zBegin, int zEnd)
1081 {
1082         Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
1083
1084         const float dt0 = _dt / _dx;
1085
1086         int begin=zBegin * _slabSize;
1087         int end=begin + (zEnd - zBegin) * _slabSize;
1088         for (int x = begin; x < end; x++)
1089                 _xForce[x] = 0.0;
1090
1091         // advectFieldMacCormack1(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res)
1092
1093         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _densityTemp, res, zBegin, zEnd);
1094         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heatTemp, res, zBegin, zEnd);
1095         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocity, res, zBegin, zEnd);
1096         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocity, res, zBegin, zEnd);
1097         advectFieldMacCormack1(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, res, zBegin, zEnd);
1098
1099         // Have to wait untill all the threads are done -> so continuing in step 3
1100 }
1101
1102 //////////////////////////////////////////////////////////////////////
1103 // Advect using the MacCormack method from the Selle paper
1104 //////////////////////////////////////////////////////////////////////
1105 void FLUID_3D::advectMacCormackEnd2(int zBegin, int zEnd)
1106 {
1107         const float dt0 = _dt / _dx;
1108         Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
1109
1110         // use force array as temp arrays
1111         float* t1 = _xForce;
1112
1113         // advectFieldMacCormack2(dt, xVelocity, yVelocity, zVelocity, oldField, newField, tempfield, temp, res, obstacles)
1114
1115         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _density, _densityTemp, t1, res, _obstacles, zBegin, zEnd);
1116         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, _heatTemp, t1, res, _obstacles, zBegin, zEnd);
1117         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocityTemp, _xVelocity, t1, res, _obstacles, zBegin, zEnd);
1118         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocityTemp, _yVelocity, t1, res, _obstacles, zBegin, zEnd);
1119         advectFieldMacCormack2(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocityTemp, _zVelocity, t1, res, _obstacles, zBegin, zEnd);
1120
1121         if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocityTemp, res, zBegin, zEnd);
1122         else setZeroX(_xVelocityTemp, res, zBegin, zEnd);                               
1123
1124         if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocityTemp, res, zBegin, zEnd);
1125         else setZeroY(_yVelocityTemp, res, zBegin, zEnd); 
1126
1127         if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocityTemp, res, zBegin, zEnd);
1128         else setZeroZ(_zVelocityTemp, res, zBegin, zEnd);
1129
1130         setZeroBorder(_density, res, zBegin, zEnd);
1131         setZeroBorder(_heat, res, zBegin, zEnd);
1132
1133
1134         /*int begin=zBegin * _slabSize;
1135         int end=begin + (zEnd - zBegin) * _slabSize;
1136   for (int x = begin; x < end; x++)
1137     _xForce[x] = _yForce[x] = 0.0f;*/
1138 }