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