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