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