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