5f40a823238cc9233b85fd93b8bf66e8ae3216d3
[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
23 #include "FLUID_3D.h"
24 #include "IMAGE.h"
25 #include <INTERPOLATE.h>
26 #include "SPHERE.h"
27 #include <zlib.h>
28
29 // boundary conditions of the fluid domain
30 #define DOMAIN_BC_FRONT  1
31 #define DOMAIN_BC_TOP    0
32 #define DOMAIN_BC_LEFT   1
33 #define DOMAIN_BC_BACK   DOMAIN_BC_FRONT
34 #define DOMAIN_BC_BOTTOM DOMAIN_BC_TOP
35 #define DOMAIN_BC_RIGHT  DOMAIN_BC_LEFT
36
37 //////////////////////////////////////////////////////////////////////
38 // Construction/Destruction
39 //////////////////////////////////////////////////////////////////////
40
41 FLUID_3D::FLUID_3D(int *res, int amplify, float *p0, float dt) :
42         _xRes(res[0]), _yRes(res[1]), _zRes(res[2]), _res(0.0f), _dt(dt)
43 {
44         // set simulation consts
45         // _dt = dt; // 0.10
46         
47         // start point of array
48         _p0[0] = p0[0];
49         _p0[1] = p0[1];
50         _p0[2] = p0[2];
51
52         _iterations = 100;
53         _tempAmb = 0; 
54         _heatDiffusion = 1e-3;
55         _vorticityEps = 2.0;
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         if(amplify)
63                 _wTurbulence = new WTURBULENCE(_res[0],_res[1],_res[2], amplify);
64         else
65                 _wTurbulence = NULL;
66         
67         // scale the constants according to the refinement of the grid
68         _dx = 1.0f / (float)_maxRes;
69         float scaling = 64.0f / _maxRes;
70         scaling = (scaling < 1.0f) ? 1.0f : scaling;
71         _vorticityEps /= scaling;
72
73         // allocate arrays
74         _totalCells   = _xRes * _yRes * _zRes;
75         _slabSize = _xRes * _yRes;
76         _divergence   = new float[_totalCells];
77         _pressure     = new float[_totalCells];
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         _vorticity    = new float[_totalCells];
88         _density      = new float[_totalCells];
89         _densityOld   = new float[_totalCells];
90         _heat         = new float[_totalCells];
91         _heatOld      = new float[_totalCells];
92         _residual     = new float[_totalCells];
93         _direction    = new float[_totalCells];
94         _q            = new float[_totalCells];
95         _obstacles    = new unsigned char[_totalCells];
96         _xVorticity   = new float[_totalCells];
97         _yVorticity   = new float[_totalCells];
98         _zVorticity   = new float[_totalCells];
99         _h                        = new float[_totalCells];
100         _Precond          = new float[_totalCells];
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                 _divergence[x]   = 0.0f;
109                 _pressure[x]     = 0.0f;
110                 _xVelocity[x]    = 0.0f;
111                 _yVelocity[x]    = 0.0f;
112                 _zVelocity[x]    = 0.0f;
113                 _xVelocityOld[x] = 0.0f;
114                 _yVelocityOld[x] = 0.0f;
115                 _zVelocityOld[x] = 0.0f;
116                 _xForce[x]       = 0.0f;
117                 _yForce[x]       = 0.0f;
118                 _zForce[x]       = 0.0f;
119                 _xVorticity[x]   = 0.0f;
120                 _yVorticity[x]   = 0.0f;
121                 _zVorticity[x]   = 0.0f;
122                 _residual[x]     = 0.0f;
123                 _h[x]                    = 0.0f;
124                 _Precond[x]              = 0.0f;
125                 _obstacles[x]    = false;
126         }
127
128         // set side obstacles
129   int index;
130   for (int y = 0; y < _yRes; y++) // z
131     for (int x = 0; x < _xRes; x++)
132     {
133       // front slab
134       index = x + y * _xRes;
135       if(DOMAIN_BC_BOTTOM==1) _obstacles[index] = 1;
136
137       // back slab
138       index += _totalCells - _slabSize;
139       if(DOMAIN_BC_TOP==1) _obstacles[index] = 1;
140     }
141   for (int z = 0; z < _zRes; z++) // y
142     for (int x = 0; x < _xRes; x++)
143     {
144       // bottom slab
145       index = x + z * _slabSize;
146       if(DOMAIN_BC_FRONT==1) _obstacles[index] = 1;
147
148       // top slab
149       index += _slabSize - _xRes;
150       if(DOMAIN_BC_BACK==1) _obstacles[index] = 1;
151     }
152   for (int z = 0; z < _zRes; z++) // x
153     for (int y = 0; y < _yRes; y++)
154     {
155       // left slab
156       index = y * _xRes + z * _slabSize;
157       if(DOMAIN_BC_LEFT==1) _obstacles[index] = 1;
158
159       // right slab
160       index += _xRes - 1;
161       if(DOMAIN_BC_RIGHT==1) _obstacles[index] = 1;
162     }
163
164         /*
165         SPHERE *obsSphere = NULL;
166         obsSphere = new SPHERE(0.375,0.5,0.375, 0.1); // for 4 to 3 domain
167         addObstacle(obsSphere);
168         delete obsSphere;
169         */
170 }
171
172 FLUID_3D::~FLUID_3D()
173 {
174         if (_divergence) delete[] _divergence;
175         if (_pressure) delete[] _pressure;
176         if (_xVelocity) delete[] _xVelocity;
177         if (_yVelocity) delete[] _yVelocity;
178         if (_zVelocity) delete[] _zVelocity;
179         if (_xVelocityOld) delete[] _xVelocityOld;
180         if (_yVelocityOld) delete[] _yVelocityOld;
181         if (_zVelocityOld) delete[] _zVelocityOld;
182         if (_xForce) delete[] _xForce;
183         if (_yForce) delete[] _yForce;
184         if (_zForce) delete[] _zForce;
185         if (_residual) delete[] _residual;
186         if (_direction) delete[] _direction;
187         if (_q)       delete[] _q;
188         if (_density) delete[] _density;
189         if (_densityOld) delete[] _densityOld;
190         if (_heat) delete[] _heat;
191         if (_heatOld) delete[] _heatOld;
192         if (_xVorticity) delete[] _xVorticity;
193         if (_yVorticity) delete[] _yVorticity;
194         if (_zVorticity) delete[] _zVorticity;
195         if (_vorticity) delete[] _vorticity;
196         if (_h) delete[] _h;
197         if (_Precond) delete[] _Precond;
198         if (_obstacles) delete[] _obstacles;
199     if (_wTurbulence) delete _wTurbulence;
200
201     printf("deleted fluid\n");
202 }
203
204 // init direct access functions from blender
205 void FLUID_3D::initBlenderRNA(float *alpha, float *beta)
206 {
207         _alpha = alpha;
208         _beta = beta;
209         
210         // XXX TODO DEBUG
211         // *_alpha = 0;
212         // *_beta = 0;
213 }
214
215 //////////////////////////////////////////////////////////////////////
216 // step simulation once
217 //////////////////////////////////////////////////////////////////////
218 void FLUID_3D::step()
219 {
220         // wipe forces
221         for (int i = 0; i < _totalCells; i++)
222                 _xForce[i] = _yForce[i] = _zForce[i] = 0.0f;
223
224         wipeBoundaries();
225
226         // run the solvers
227         addVorticity();
228         addBuoyancy(_heat, _density);
229         addForce();
230         project();
231         diffuseHeat();
232
233         // advect everything
234         advectMacCormack();
235
236         if(_wTurbulence) {
237                 _wTurbulence->stepTurbulenceFull(_dt/_dx,
238                                 _xVelocity, _yVelocity, _zVelocity, _obstacles);
239                 // _wTurbulence->stepTurbulenceReadable(_dt/_dx,
240                 //  _xVelocity, _yVelocity, _zVelocity, _obstacles);
241         }
242 /*
243  // no file output
244   float *src = _density;
245         string prefix = string("./original.preview/density_fullxy_");
246         writeImageSliceXY(src,_res, _res[2]/2, prefix, _totalSteps);
247 */
248         // artificial damping -- this is necessary because we use a
249   // collated grid, and at very coarse grid resolutions, banding
250   // artifacts can occur
251         artificialDamping(_xVelocity);
252         artificialDamping(_yVelocity);
253         artificialDamping(_zVelocity);
254 /*
255 // no file output
256   string pbrtPrefix = string("./pbrt/density_small_");
257   IMAGE::dumpPBRT(_totalSteps, pbrtPrefix, _density, _res[0],_res[1],_res[2]);
258   */
259         _totalTime += _dt;
260         _totalSteps++;
261
262         // clear obstacles
263         for (int i = 0; i < _totalCells; i++)
264         {
265                 if(_obstacles[i])
266                 {
267                         _xVelocity[i] =
268                         _yVelocity[i] =
269                         _zVelocity[i] = 0.0f;
270                         _pressure[i] = 0.0f;
271                 }
272                 _obstacles[i] = 0;
273         }
274 }
275
276 //////////////////////////////////////////////////////////////////////
277 // helper function to dampen co-located grid artifacts of given arrays in intervals
278 // (only needed for velocity, strength (w) depends on testcase...
279 //////////////////////////////////////////////////////////////////////
280 void FLUID_3D::artificialDamping(float* field) {
281         const float w = 0.9;
282         if(_totalSteps % 4 == 1) {
283                 for (int z = 1; z < _res[2]-1; z++)
284                         for (int y = 1; y < _res[1]-1; y++)
285                                 for (int x = 1+(y+z)%2; x < _res[0]-1; x+=2) {
286                                         const int index = x + y*_res[0] + z * _slabSize;
287                                         field[index] = (1-w)*field[index] + 1./6. * w*(
288                                                         field[index+1] + field[index-1] +
289                                                         field[index+_res[0]] + field[index-_res[0]] +
290                                                         field[index+_slabSize] + field[index-_slabSize] );
291                                 }
292         }
293         if(_totalSteps % 4 == 3) {
294                 for (int z = 1; z < _res[2]-1; z++)
295                         for (int y = 1; y < _res[1]-1; y++)
296                                 for (int x = 1+(y+z+1)%2; x < _res[0]-1; x+=2) {
297                                         const int index = x + y*_res[0] + z * _slabSize;
298                                         field[index] = (1-w)*field[index] + 1./6. * w*(
299                                                         field[index+1] + field[index-1] +
300                                                         field[index+_res[0]] + field[index-_res[0]] +
301                                                         field[index+_slabSize] + field[index-_slabSize] );
302                                 }
303         }
304 }
305
306 //////////////////////////////////////////////////////////////////////
307 // copy out the boundary in all directions
308 //////////////////////////////////////////////////////////////////////
309 void FLUID_3D::copyBorderAll(float* field)
310 {
311         int index;
312         for (int y = 0; y < _yRes; y++)
313                 for (int x = 0; x < _xRes; x++)
314                 {
315                         // front slab
316                         index = x + y * _xRes;
317                         field[index] = field[index + _slabSize];
318
319                         // back slab
320                         index += _totalCells - _slabSize;
321                         field[index] = field[index - _slabSize];
322     }
323
324         for (int z = 0; z < _zRes; z++)
325                 for (int x = 0; x < _xRes; x++)
326     {
327                         // bottom slab
328                         index = x + z * _slabSize;
329                         field[index] = field[index + _xRes];
330
331                         // top slab
332                         index += _slabSize - _xRes;
333                         field[index] = field[index - _xRes];
334     }
335
336         for (int z = 0; z < _zRes; z++)
337                 for (int y = 0; y < _yRes; y++)
338     {
339                         // left slab
340                         index = y * _xRes + z * _slabSize;
341                         field[index] = field[index + 1];
342
343                         // right slab
344                         index += _xRes - 1;
345                         field[index] = field[index - 1];
346                 }
347 }
348
349 //////////////////////////////////////////////////////////////////////
350 // wipe boundaries of velocity and density
351 //////////////////////////////////////////////////////////////////////
352 void FLUID_3D::wipeBoundaries()
353 {
354         setZeroBorder(_xVelocity, _res);
355         setZeroBorder(_yVelocity, _res);
356         setZeroBorder(_zVelocity, _res);
357         setZeroBorder(_density, _res);
358 }
359
360 //////////////////////////////////////////////////////////////////////
361 // add forces to velocity field
362 //////////////////////////////////////////////////////////////////////
363 void FLUID_3D::addForce()
364 {
365         for (int i = 0; i < _totalCells; i++)
366         {
367                 _xVelocity[i] += _dt * _xForce[i];
368                 _yVelocity[i] += _dt * _yForce[i];
369                 _zVelocity[i] += _dt * _zForce[i];
370         }
371 }
372
373 //////////////////////////////////////////////////////////////////////
374 // project into divergence free field
375 //////////////////////////////////////////////////////////////////////
376 void FLUID_3D::project()
377 {
378         int index, x, y, z;
379
380         setObstacleBoundaries();
381
382         // copy out the boundaries
383         if(DOMAIN_BC_LEFT == 0)  setNeumannX(_xVelocity, _res);
384         else setZeroX(_xVelocity, _res);
385
386         if(DOMAIN_BC_TOP == 0)   setNeumannZ(_zVelocity, _res);
387         else setZeroZ(_zVelocity, _res);
388
389         if(DOMAIN_BC_FRONT == 0) setNeumannY(_yVelocity, _res);
390         else setZeroY(_yVelocity, _res);
391
392         // calculate divergence
393         index = _slabSize + _xRes + 1;
394         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
395                 for (y = 1; y < _yRes - 1; y++, index += 2)
396                         for (x = 1; x < _xRes - 1; x++, index++)
397                         {
398                                 float xright = _xVelocity[index + 1];
399                                 float xleft  = _xVelocity[index - 1];
400                                 float yup    = _yVelocity[index + _xRes];
401                                 float ydown  = _yVelocity[index - _xRes];
402                                 float ztop   = _zVelocity[index + _slabSize];
403                                 float zbottom = _zVelocity[index - _slabSize];
404
405                                 if(_obstacles[index+1]) xright = - _xVelocity[index];
406                                 if(_obstacles[index-1]) xleft  = - _xVelocity[index];
407                                 if(_obstacles[index+_xRes]) yup    = - _yVelocity[index];
408                                 if(_obstacles[index-_xRes]) ydown  = - _yVelocity[index];
409                                 if(_obstacles[index+_slabSize]) ztop    = - _zVelocity[index];
410                                 if(_obstacles[index-_slabSize]) zbottom = - _zVelocity[index];
411
412                                 _divergence[index] = -_dx * 0.5f * (
413                                                 xright - xleft +
414                                                 yup - ydown +
415                                                 ztop - zbottom );
416
417                                 // DG: commenting this helps CG to get a better start, 10-20% speed improvement
418                                 // _pressure[index] = 0.0f;
419                         }
420         copyBorderAll(_pressure);
421
422         // solve Poisson equation
423         solvePressurePre(_pressure, _divergence, _obstacles);
424
425         setObstaclePressure();
426
427         // project out solution
428         float invDx = 1.0f / _dx;
429         index = _slabSize + _xRes + 1;
430         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
431                 for (y = 1; y < _yRes - 1; y++, index += 2)
432                         for (x = 1; x < _xRes - 1; x++, index++)
433                         {
434                                 if(!_obstacles[index])
435                                 {
436                                         _xVelocity[index] -= 0.5f * (_pressure[index + 1]     - _pressure[index - 1]) * invDx;
437                                         _yVelocity[index] -= 0.5f * (_pressure[index + _xRes]  - _pressure[index - _xRes]) * invDx;
438                                         _zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
439                                 }
440                         }
441 }
442
443 //////////////////////////////////////////////////////////////////////
444 // diffuse heat
445 //////////////////////////////////////////////////////////////////////
446 void FLUID_3D::diffuseHeat()
447 {
448         SWAP_POINTERS(_heat, _heatOld);
449
450         copyBorderAll(_heatOld);
451         solveHeat(_heat, _heatOld, _obstacles);
452
453         // zero out inside obstacles
454         for (int x = 0; x < _totalCells; x++)
455                 if (_obstacles[x])
456                         _heat[x] = 0.0f;
457 }
458
459 //////////////////////////////////////////////////////////////////////
460 // stamp an obstacle in the _obstacles field
461 //////////////////////////////////////////////////////////////////////
462 void FLUID_3D::addObstacle(OBSTACLE* obstacle)
463 {
464         int index = 0;
465         for (int z = 0; z < _zRes; z++)
466                 for (int y = 0; y < _yRes; y++)
467                         for (int x = 0; x < _xRes; x++, index++)
468                                 if (obstacle->inside(x * _dx, y * _dx, z * _dx)) {
469                                         _obstacles[index] = true;
470         }
471 }
472
473 //////////////////////////////////////////////////////////////////////
474 // calculate the obstacle directional types
475 //////////////////////////////////////////////////////////////////////
476 void FLUID_3D::setObstaclePressure()
477 {
478         // tag remaining obstacle blocks
479         for (int z = 1, index = _slabSize + _xRes + 1;
480                         z < _zRes - 1; z++, index += 2 * _xRes)
481                 for (int y = 1; y < _yRes - 1; y++, index += 2)
482                         for (int x = 1; x < _xRes - 1; x++, index++)
483                 {
484                         // could do cascade of ifs, but they are a pain
485                         if (_obstacles[index])
486                         {
487                                 const int top   = _obstacles[index + _slabSize];
488                                 const int bottom= _obstacles[index - _slabSize];
489                                 const int up    = _obstacles[index + _xRes];
490                                 const int down  = _obstacles[index - _xRes];
491                                 const int left  = _obstacles[index - 1];
492                                 const int right = _obstacles[index + 1];
493
494                                 // unused
495                                 // const bool fullz = (top && bottom);
496                                 // const bool fully = (up && down);
497                                 //const bool fullx = (left && right);
498
499                                 _xVelocity[index] =
500                                 _yVelocity[index] =
501                                 _zVelocity[index] = 0.0f;
502                                 _pressure[index] = 0.0f;
503
504                                 // average pressure neighbors
505                                 float pcnt = 0., vp = 0.;
506                                 if (left && !right) {
507                                         _pressure[index] += _pressure[index + 1];
508                                         pcnt += 1.;
509                                 }
510                                 if (!left && right) {
511                                         _pressure[index] += _pressure[index - 1];
512                                         pcnt += 1.;
513                                 }
514                                 if (up && !down) {
515                                         _pressure[index] += _pressure[index - _xRes];
516                                         pcnt += 1.;
517                                 }
518                                 if (!up && down) {
519                                         _pressure[index] += _pressure[index + _xRes];
520                                         pcnt += 1.;
521                                 }
522                                 if (top && !bottom) {
523                                         _pressure[index] += _pressure[index - _slabSize];
524                                         pcnt += 1.;
525                                         // _zVelocity[index] +=  - _zVelocity[index - _slabSize];
526                                         // vp += 1.0;
527                                 }
528                                 if (!top && bottom) {
529                                         _pressure[index] += _pressure[index + _slabSize];
530                                         pcnt += 1.;
531                                         // _zVelocity[index] +=  - _zVelocity[index + _slabSize];
532                                         // vp += 1.0;
533                                 }
534                                 
535                                 if(pcnt > 0.000001f)
536                                         _pressure[index] /= pcnt;
537
538                                 // test - dg
539                                 // if(vp > 0.000001f)
540                                 //      _zVelocity[index] /= vp;
541
542                                 // TODO? set correct velocity bc's
543                                 // velocities are only set to zero right now
544                                 // this means it's not a full no-slip boundary condition
545                                 // but a "half-slip" - still looks ok right now
546                         }
547                 }
548 }
549
550 void FLUID_3D::setObstacleBoundaries()
551 {
552         // cull degenerate obstacles , move to addObstacle?
553         for (int z = 1, index = _slabSize + _xRes + 1;
554                         z < _zRes - 1; z++, index += 2 * _xRes)
555                 for (int y = 1; y < _yRes - 1; y++, index += 2)
556                         for (int x = 1; x < _xRes - 1; x++, index++)
557                         {
558                                 if (_obstacles[index] != EMPTY)
559                                 {
560                                         const int top   = _obstacles[index + _slabSize];
561                                         const int bottom= _obstacles[index - _slabSize];
562                                         const int up    = _obstacles[index + _xRes];
563                                         const int down  = _obstacles[index - _xRes];
564                                         const int left  = _obstacles[index - 1];
565                                         const int right = _obstacles[index + 1];
566
567                                         int counter = 0;
568                                         if (up)    counter++;
569                                         if (down)  counter++;
570                                         if (left)  counter++;
571                                         if (right) counter++;
572                                         if (top)  counter++;
573                                         if (bottom) counter++;
574
575                                         if (counter < 3)
576                                                 _obstacles[index] = EMPTY;
577                                 }
578                                 if (_obstacles[index])
579                                 {
580                                         _xVelocity[index] =
581                                         _yVelocity[index] =
582                                         _zVelocity[index] = 0.0f;
583                                         _pressure[index] = 0.0f;
584                                 }
585                         }
586 }
587
588 //////////////////////////////////////////////////////////////////////
589 // add buoyancy forces
590 //////////////////////////////////////////////////////////////////////
591 void FLUID_3D::addBuoyancy(float *heat, float *density)
592 {
593         int index = 0;
594
595         for (int z = 0; z < _zRes; z++)
596                 for (int y = 0; y < _yRes; y++)
597                         for (int x = 0; x < _xRes; x++, index++)
598                         {
599                                 _zForce[index] += *_alpha * density[index] + (*_beta * (heat[index] - _tempAmb)); // DG: was _yForce, changed for Blender
600                         }
601 }
602
603 //////////////////////////////////////////////////////////////////////
604 // add vorticity to the force field
605 //////////////////////////////////////////////////////////////////////
606 void FLUID_3D::addVorticity()
607 {
608         int x,y,z,index;
609         if(_vorticityEps<=0.) return;
610
611         // calculate vorticity
612         float gridSize = 0.5f / _dx;
613         index = _slabSize + _xRes + 1;
614         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
615                 for (y = 1; y < _yRes - 1; y++, index += 2)
616                         for (x = 1; x < _xRes - 1; x++, index++)
617                         {
618                                 int up    = _obstacles[index + _xRes] ? index : index + _xRes;
619                                 int down  = _obstacles[index - _xRes] ? index : index - _xRes;
620                                 float dy  = (up == index || down == index) ? 1.0f / _dx : gridSize;
621                                 int out   = _obstacles[index + _slabSize] ? index : index + _slabSize;
622                                 int in    = _obstacles[index - _slabSize] ? index : index - _slabSize;
623                                 float dz  = (out == index || in == index) ? 1.0f / _dx : gridSize;
624                                 int right = _obstacles[index + 1] ? index : index + 1;
625                                 int left  = _obstacles[index - 1] ? index : index - 1;
626                                 float dx  = (right == index || right == index) ? 1.0f / _dx : gridSize;
627
628                                 _xVorticity[index] = (_zVelocity[up] - _zVelocity[down]) * dy + (-_yVelocity[out] + _yVelocity[in]) * dz;
629                                 _yVorticity[index] = (_xVelocity[out] - _xVelocity[in]) * dz + (-_zVelocity[right] + _zVelocity[left]) * dx;
630                                 _zVorticity[index] = (_yVelocity[right] - _yVelocity[left]) * dx + (-_xVelocity[up] + _xVelocity[down])* dy;
631
632                                 _vorticity[index] = sqrtf(_xVorticity[index] * _xVorticity[index] +
633                                                 _yVorticity[index] * _yVorticity[index] +
634                                                 _zVorticity[index] * _zVorticity[index]);
635                         }
636
637         // calculate normalized vorticity vectors
638         float eps = _vorticityEps;
639         index = _slabSize + _xRes + 1;
640         for (z = 1; z < _zRes - 1; z++, index += 2 * _xRes)
641                 for (y = 1; y < _yRes - 1; y++, index += 2)
642                         for (x = 1; x < _xRes - 1; x++, index++)
643                                 if (!_obstacles[index])
644                                 {
645                                         float N[3];
646
647                                         int up    = _obstacles[index + _xRes] ? index : index + _xRes;
648                                         int down  = _obstacles[index - _xRes] ? index : index - _xRes;
649                                         float dy  = (up == index || down == index) ? 1.0f / _dx : gridSize;
650                                         int out   = _obstacles[index + _slabSize] ? index : index + _slabSize;
651                                         int in    = _obstacles[index - _slabSize] ? index : index - _slabSize;
652                                         float dz  = (out == index || in == index) ? 1.0f / _dx : gridSize;
653                                         int right = _obstacles[index + 1] ? index : index + 1;
654                                         int left  = _obstacles[index - 1] ? index : index - 1;
655                                         float dx  = (right == index || right == index) ? 1.0f / _dx : gridSize;
656                                         N[0] = (_vorticity[right] - _vorticity[left]) * dx;
657                                         N[1] = (_vorticity[up] - _vorticity[down]) * dy;
658                                         N[2] = (_vorticity[out] - _vorticity[in]) * dz;
659
660                                         float magnitude = sqrtf(N[0] * N[0] + N[1] * N[1] + N[2] * N[2]);
661                                         if (magnitude > 0.0f)
662                                         {
663                                                 magnitude = 1.0f / magnitude;
664                                                 N[0] *= magnitude;
665                                                 N[1] *= magnitude;
666                                                 N[2] *= magnitude;
667
668                                                 _xForce[index] += (N[1] * _zVorticity[index] - N[2] * _yVorticity[index]) * _dx * eps;
669                                                 _yForce[index] -= (N[0] * _zVorticity[index] - N[2] * _xVorticity[index]) * _dx * eps;
670                                                 _zForce[index] += (N[0] * _yVorticity[index] - N[1] * _xVorticity[index]) * _dx * eps;
671                                         }
672                                 }
673 }
674
675 //////////////////////////////////////////////////////////////////////
676 // Advect using the MacCormack method from the Selle paper
677 //////////////////////////////////////////////////////////////////////
678 void FLUID_3D::advectMacCormack()
679 {
680         Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
681
682         if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocity, res);
683         else setZeroX(_xVelocity, res);
684
685         if(DOMAIN_BC_TOP == 0) copyBorderZ(_zVelocity, res);
686         else setZeroZ(_zVelocity, res);
687
688         if(DOMAIN_BC_FRONT == 0) copyBorderY(_yVelocity, res);
689         else setZeroY(_yVelocity, res);
690
691         SWAP_POINTERS(_xVelocity, _xVelocityOld);
692         SWAP_POINTERS(_yVelocity, _yVelocityOld);
693         SWAP_POINTERS(_zVelocity, _zVelocityOld);
694         SWAP_POINTERS(_density, _densityOld);
695         SWAP_POINTERS(_heat, _heatOld);
696
697         const float dt0 = _dt / _dx;
698         // use force arrays as temp arrays
699   for (int x = 0; x < _totalCells; x++)
700     _xForce[x] = _yForce[x] = 0.0;
701         float* t1 = _xForce;
702         float* t2 = _yForce;
703
704         advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _densityOld, _density, t1,t2, res, NULL);
705         advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _heatOld, _heat, t1,t2, res, NULL);
706         advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _xVelocityOld, _xVelocity, t1,t2, res, NULL);
707         advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _yVelocityOld, _yVelocity, t1,t2, res, NULL);
708         advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, t1,t2, res, NULL);
709
710         if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocity, res);
711         else setZeroX(_xVelocity, res);
712
713         if(DOMAIN_BC_TOP == 0) copyBorderZ(_zVelocity, res);
714         else setZeroZ(_zVelocity, res);
715
716         if(DOMAIN_BC_FRONT == 0) copyBorderY(_yVelocity, res);
717         else setZeroY(_yVelocity, res);
718
719         setZeroBorder(_density, res);
720         setZeroBorder(_heat, res);
721
722   for (int x = 0; x < _totalCells; x++)
723     t1[x] = t2[x] = 0.0;
724 }