Fix: smoke noise tile was saved in Blender executable directory, which is often write...
[blender.git] / intern / smoke / intern / WTURBULENCE.cpp
1 /** \file smoke/intern/WTURBULENCE.cpp
2  *  \ingroup smoke
3  */
4 //////////////////////////////////////////////////////////////////////
5 // This file is part of Wavelet Turbulence.
6 // 
7 // Wavelet Turbulence is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 // 
12 // Wavelet Turbulence is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 // GNU General Public License for more details.
16 // 
17 // You should have received a copy of the GNU General Public License
18 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
19 // 
20 // Copyright 2008 Theodore Kim and Nils Thuerey
21 // 
22 // WTURBULENCE handling
23 ///////////////////////////////////////////////////////////////////////////////////
24 // Parallelized turbulence even further. TNT matrix library functions
25 // rewritten to improve performance.
26 //              - MiikaH
27 //////////////////////////////////////////////////////////////////////
28
29 #include "WTURBULENCE.h"
30 #include "INTERPOLATE.h"
31 #include "IMAGE.h"
32 #include <MERSENNETWISTER.h>
33 #include "WAVELET_NOISE.h"
34 #include "FFT_NOISE.h"
35 #include "EIGENVALUE_HELPER.h"
36 #include "LU_HELPER.h"
37 #include "SPHERE.h"
38 #include <zlib.h>
39 #include <math.h>
40
41 // needed to access static advection functions
42 #include "FLUID_3D.h"
43
44 #if PARALLEL==1
45 #include <omp.h>
46 #endif // PARALLEL 
47
48 // 2^ {-5/6}
49 static const float persistence = 0.56123f;
50
51 //////////////////////////////////////////////////////////////////////
52 // constructor
53 //////////////////////////////////////////////////////////////////////
54 WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype, const char *noisefile_path, int init_fire, int init_colors)
55 {
56         // if noise magnitude is below this threshold, its contribution
57         // is negilgible, so stop evaluating new octaves
58         _cullingThreshold = 1e-3;
59         
60         // factor by which to increase the simulation resolution
61         _amplify = amplify;
62         
63         // manually adjust the overall amount of turbulence
64         // DG - RNA-fied _strength = 2.;
65         
66         // add the corresponding octaves of noise
67         _octaves = (int)(log((float)_amplify) / log(2.0f) + 0.5f); // XXX DEBUG/ TODO: int casting correct? - dg
68         
69         // noise resolution
70         _xResBig = _amplify * xResSm;
71         _yResBig = _amplify * yResSm;
72         _zResBig = _amplify * zResSm;
73         _resBig = Vec3Int(_xResBig, _yResBig, _zResBig);
74         _invResBig = Vec3(1.0f/(float)_resBig[0], 1.0f/(float)_resBig[1], 1.0f/(float)_resBig[2]);
75         _slabSizeBig = _xResBig*_yResBig;
76         _totalCellsBig = _slabSizeBig * _zResBig;
77         
78         // original / small resolution
79         _xResSm = xResSm;
80         _yResSm = yResSm;
81         _zResSm = zResSm;
82         _resSm = Vec3Int(xResSm, yResSm, zResSm);
83         _invResSm = Vec3(1.0f/(float)_resSm[0], 1.0f/(float)_resSm[1], 1.0f/(float)_resSm[2] );
84         _slabSizeSm = _xResSm*_yResSm;
85         _totalCellsSm = _slabSizeSm * _zResSm;
86         
87         // allocate high resolution density field
88         _totalStepsBig = 0;
89         _densityBig = new float[_totalCellsBig];
90         _densityBigOld = new float[_totalCellsBig];
91         
92         for(int i = 0; i < _totalCellsBig; i++) {
93                 _densityBig[i] = 
94                 _densityBigOld[i] = 0.;
95         }
96
97         /* fire */
98         _flameBig = _fuelBig = _fuelBigOld = NULL;
99         _reactBig = _reactBigOld = NULL;
100         if (init_fire) {
101                 initFire();
102         }
103         /* colors */
104         _color_rBig = _color_rBigOld = NULL;
105         _color_gBig = _color_gBigOld = NULL;
106         _color_bBig = _color_bBigOld = NULL;
107         if (init_colors) {
108                 initColors(0.0f, 0.0f, 0.0f);
109         }
110         
111         // allocate & init texture coordinates
112         _tcU = new float[_totalCellsSm];
113         _tcV = new float[_totalCellsSm];
114         _tcW = new float[_totalCellsSm];
115         _tcTemp = new float[_totalCellsSm];
116         
117         // map all 
118         const float dx = 1.0f/(float)(_resSm[0]);
119         const float dy = 1.0f/(float)(_resSm[1]);
120         const float dz = 1.0f/(float)(_resSm[2]);
121         int index = 0;
122         for (int z = 0; z < _zResSm; z++) 
123         for (int y = 0; y < _yResSm; y++) 
124                 for (int x = 0; x < _xResSm; x++, index++)
125                 {
126                 _tcU[index] = x*dx;
127                 _tcV[index] = y*dy;
128                 _tcW[index] = z*dz;
129                 _tcTemp[index] = 0.;
130                 }
131         
132         // noise tiles
133         _noiseTile = new float[noiseTileSize * noiseTileSize * noiseTileSize];
134         setNoise(noisetype, noisefile_path);
135 }
136
137 void WTURBULENCE::initFire()
138 {
139         if (!_fuelBig) {
140                 _flameBig = new float[_totalCellsBig];
141                 _fuelBig = new float[_totalCellsBig];
142                 _fuelBigOld = new float[_totalCellsBig];
143                 _reactBig = new float[_totalCellsBig];
144                 _reactBigOld = new float[_totalCellsBig];
145
146                 for(int i = 0; i < _totalCellsBig; i++) {
147                         _flameBig[i] = 
148                         _fuelBig[i] = 
149                         _fuelBigOld[i] = 0.;
150                         _reactBig[i] = 
151                         _reactBigOld[i] = 0.;
152                 }
153         }
154 }
155
156 void WTURBULENCE::initColors(float init_r, float init_g, float init_b)
157 {
158         if (!_color_rBig) {
159                 _color_rBig = new float[_totalCellsBig];
160                 _color_rBigOld = new float[_totalCellsBig];
161                 _color_gBig = new float[_totalCellsBig];
162                 _color_gBigOld = new float[_totalCellsBig];
163                 _color_bBig = new float[_totalCellsBig];
164                 _color_bBigOld = new float[_totalCellsBig];
165
166                 for(int i = 0; i < _totalCellsBig; i++) {
167                         _color_rBig[i] = _densityBig[i] * init_r;
168                         _color_rBigOld[i] = 0.0f;
169                         _color_gBig[i] = _densityBig[i] * init_g;
170                         _color_gBigOld[i] = 0.0f;
171                         _color_bBig[i] = _densityBig[i] * init_b;
172                         _color_bBigOld[i] = 0.0f;
173                 }
174         }
175 }
176
177 //////////////////////////////////////////////////////////////////////
178 // destructor
179 //////////////////////////////////////////////////////////////////////
180 WTURBULENCE::~WTURBULENCE() {
181   delete[] _densityBig;
182   delete[] _densityBigOld;
183   if (_flameBig) delete[] _flameBig;
184   if (_fuelBig) delete[] _fuelBig;
185   if (_fuelBigOld) delete[] _fuelBigOld;
186   if (_reactBig) delete[] _reactBig;
187   if (_reactBigOld) delete[] _reactBigOld;
188
189   if (_color_rBig) delete[] _color_rBig;
190   if (_color_rBigOld) delete[] _color_rBigOld;
191   if (_color_gBig) delete[] _color_gBig;
192   if (_color_gBigOld) delete[] _color_gBigOld;
193   if (_color_bBig) delete[] _color_bBig;
194   if (_color_bBigOld) delete[] _color_bBigOld;
195
196   delete[] _tcU;
197   delete[] _tcV;
198   delete[] _tcW;
199   delete[] _tcTemp;
200
201   delete[] _noiseTile;
202 }
203
204 //////////////////////////////////////////////////////////////////////
205 // Change noise type
206 //
207 // type (1<<0) = wavelet / 2
208 // type (1<<1) = FFT / 4
209 // type (1<<2) = curl / 8
210 //////////////////////////////////////////////////////////////////////
211 void WTURBULENCE::setNoise(int type, const char *noisefile_path)
212 {
213         if(type == (1<<1)) // FFT
214         {
215 #ifdef WITH_FFTW3
216                 // needs fft
217                 std::string noiseTileFilename = std::string(noisefile_path) + std::string("noise.fft");
218                 generatTile_FFT(_noiseTile, noiseTileFilename);
219                 return;
220 #else
221                 fprintf(stderr, "FFTW not enabled, falling back to wavelet noise.\n");
222 #endif
223         }
224 #if 0
225         if(type == (1<<2)) // curl
226         {
227                 // TODO: not supported yet
228                 return;
229         }
230 #endif
231
232         std::string noiseTileFilename = std::string(noisefile_path) + std::string("noise.wavelets");
233         generateTile_WAVELET(_noiseTile, noiseTileFilename);
234 }
235
236 // init direct access functions from blender
237 void WTURBULENCE::initBlenderRNA(float *strength)
238 {
239         _strength = strength;
240 }
241
242 //////////////////////////////////////////////////////////////////////
243 // Get the smallest valid x derivative
244 //
245 // Takes the one-sided finite difference in both directions and
246 // selects the smaller of the two
247 //////////////////////////////////////////////////////////////////////
248 static float minDx(int x, int y, int z, float* input, Vec3Int res)
249 {
250   const int index = x + y * res[0] + z * res[0] * res[1];
251   const int maxx = res[0]-2;
252
253   // get grid values
254   float center = input[index];
255   float left  = (x <= 1)    ? FLT_MAX : input[index - 1];
256   float right = (x >= maxx) ? FLT_MAX : input[index + 1];
257
258   const float dx = res[0];
259
260   // get all the derivative estimates
261   float dLeft   = (x <= 1)     ? FLT_MAX : (center - left) * dx;
262   float dRight  = (x >= maxx)  ? FLT_MAX : (right - center) * dx;
263   float dCenter = (x <= 1 || x >= maxx) ? FLT_MAX : (right - left) * dx * 0.5f;
264
265   // if it's on a boundary, only one estimate is valid
266   if (x <= 1) return dRight;
267   if (x >= maxx) return dLeft;
268
269   // if it's not on a boundary, get the smallest one
270   float finalD;
271   finalD = (fabs(dCenter) < fabs(dRight)) ? dCenter : dRight;
272   finalD = (fabs(finalD)  < fabs(dLeft))  ? finalD  : dLeft;
273
274   return finalD;
275 }
276
277 //////////////////////////////////////////////////////////////////////
278 // get the smallest valid y derivative
279 //
280 // Takes the one-sided finite difference in both directions and
281 // selects the smaller of the two
282 //////////////////////////////////////////////////////////////////////
283 static float minDy(int x, int y, int z, float* input, Vec3Int res)
284 {
285   const int index = x + y * res[0] + z * res[0] * res[1];
286   const int maxy = res[1]-2;
287
288   // get grid values
289   float center = input[index];
290   float down  = (y <= 1) ? FLT_MAX : input[index - res[0]];
291   float up = (y >= maxy) ? FLT_MAX : input[index + res[0]];
292
293   const float dx = res[1]; // only for square domains
294
295   // get all the derivative estimates
296   float dDown   = (y <= 1)  ? FLT_MAX : (center - down) * dx;
297   float dUp  = (y >= maxy)  ? FLT_MAX : (up - center) * dx;
298   float dCenter = (y <= 1 || y >= maxy) ? FLT_MAX : (up - down) * dx * 0.5f;
299
300   // if it's on a boundary, only one estimate is valid
301   if (y <= 1) return dUp;
302   if (y >= maxy) return dDown;
303
304   // if it's not on a boundary, get the smallest one
305   float finalD = (fabs(dCenter) < fabs(dUp)) ? dCenter : dUp;
306   finalD = (fabs(finalD) < fabs(dDown)) ? finalD : dDown;
307
308   return finalD;
309 }
310
311 //////////////////////////////////////////////////////////////////////
312 // get the smallest valid z derivative
313 //
314 // Takes the one-sided finite difference in both directions and
315 // selects the smaller of the two
316 //////////////////////////////////////////////////////////////////////
317 static float minDz(int x, int y, int z, float* input, Vec3Int res)
318 {
319   const int slab = res[0]*res[1];
320   const int index = x + y * res[0] + z * slab;
321   const int maxz = res[2]-2;
322
323   // get grid values
324   float center = input[index];
325   float front  = (z <= 1) ? FLT_MAX : input[index - slab];
326   float back = (z >= maxz) ? FLT_MAX : input[index + slab];
327
328   const float dx = res[2]; // only for square domains
329
330   // get all the derivative estimates
331   float dfront   = (z <= 1)  ? FLT_MAX : (center - front) * dx;
332   float dback  = (z >= maxz)  ? FLT_MAX : (back - center) * dx;
333   float dCenter = (z <= 1 || z >= maxz) ? FLT_MAX : (back - front) * dx * 0.5f;
334
335   // if it's on a boundary, only one estimate is valid
336   if (z <= 1) return dback;
337   if (z >= maxz) return dfront;
338
339   // if it's not on a boundary, get the smallest one
340   float finalD = (fabs(dCenter) < fabs(dback)) ? dCenter : dback;
341   finalD = (fabs(finalD) < fabs(dfront)) ? finalD : dfront;
342
343   return finalD;
344 }
345
346 //////////////////////////////////////////////////////////////////////
347 // handle texture coordinates (advection, reset, eigenvalues), 
348 // Beware -- uses big density maccormack as temporary arrays
349 ////////////////////////////////////////////////////////////////////// 
350 void WTURBULENCE::advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *tempBig1, float *tempBig2) {
351
352   // advection
353   SWAP_POINTERS(_tcTemp, _tcU);
354   FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
355   FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
356   FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
357   FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel, 
358       _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
359   FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel, 
360       _tcTemp, _tcU, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
361
362   SWAP_POINTERS(_tcTemp, _tcV);
363   FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
364   FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
365   FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
366   FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel, 
367       _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
368   FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel, 
369       _tcTemp, _tcV, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
370
371   SWAP_POINTERS(_tcTemp, _tcW);
372   FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
373   FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
374   FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
375   FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel, 
376       _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
377   FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel, 
378       _tcTemp, _tcW, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
379 }
380
381 //////////////////////////////////////////////////////////////////////
382 // Compute the eigenvalues of the advected texture
383 ////////////////////////////////////////////////////////////////////// 
384 void WTURBULENCE::computeEigenvalues(float *_eigMin, float *_eigMax) {
385   // stats
386   float maxeig = -1.;
387   float mineig = 10.;
388
389   // texture coordinate eigenvalues
390   for (int z = 1; z < _zResSm-1; z++) {
391     for (int y = 1; y < _yResSm-1; y++) 
392       for (int x = 1; x < _xResSm-1; x++)
393       {
394         const int index = x+ y *_resSm[0] + z*_slabSizeSm;
395
396         // compute jacobian
397         float jacobian[3][3] = {
398           { minDx(x, y, z, _tcU, _resSm), minDx(x, y, z, _tcV, _resSm), minDx(x, y, z, _tcW, _resSm) } ,
399           { minDy(x, y, z, _tcU, _resSm), minDy(x, y, z, _tcV, _resSm), minDy(x, y, z, _tcW, _resSm) } ,
400           { minDz(x, y, z, _tcU, _resSm), minDz(x, y, z, _tcV, _resSm), minDz(x, y, z, _tcW, _resSm) }
401         };
402
403         // ONLY compute the eigenvalues after checking that the matrix
404         // is nonsingular
405         sLU LU = computeLU(jacobian);
406
407         if (isNonsingular(LU))
408         {
409           // get the analytic eigenvalues, quite slow right now...
410           Vec3 eigenvalues = Vec3(1.);
411           computeEigenvalues3x3( &eigenvalues[0], jacobian);
412           _eigMax[index] = MAX3V(eigenvalues);
413           _eigMin[index] = MIN3V(eigenvalues);
414           maxeig = MAX(_eigMax[index],maxeig);
415           mineig = MIN(_eigMin[index],mineig);
416         }
417         else
418         {
419           _eigMax[index] = 10.0f;
420           _eigMin[index] = 0.1;
421         }
422       }
423   }
424 }
425
426 //////////////////////////////////////////////////////////////////////
427 // advect & reset texture coordinates based on eigenvalues
428 ////////////////////////////////////////////////////////////////////// 
429 void WTURBULENCE::resetTextureCoordinates(float *_eigMin, float *_eigMax) 
430 {
431   // allowed deformation of the textures
432   const float limit = 2.f;
433   const float limitInv = 1.0f/limit;
434
435   // standard reset
436   int resets = 0;
437   const float dx = 1.0f/(float)(_resSm[0]);
438   const float dy = 1.0f/(float)(_resSm[1]);
439   const float dz = 1.0f/(float)(_resSm[2]);
440
441   for (int z = 1; z < _zResSm-1; z++)
442     for (int y = 1; y < _yResSm-1; y++)
443       for (int x = 1; x < _xResSm-1; x++)
444       {
445         const int index = x+ y *_resSm[0] + z*_slabSizeSm;
446         if (_eigMax[index] > limit || _eigMin[index] < limitInv)
447         {
448           _tcU[index] = (float)x * dx;
449           _tcV[index] = (float)y * dy;
450           _tcW[index] = (float)z * dz;
451           resets++;
452         }
453       }
454 }
455
456 //////////////////////////////////////////////////////////////////////
457 // Compute the highest frequency component of the wavelet
458 // decomposition
459 //////////////////////////////////////////////////////////////////////
460 void WTURBULENCE::decomposeEnergy(float *_energy, float *_highFreqEnergy)
461 {
462   // do the decomposition -- the goal here is to have
463   // the energy with the high frequency component stomped out
464   // stored in _tcTemp when it is done. _highFreqEnergy is only used
465   // as an additional temp array
466   
467   // downsample input
468   downsampleXNeumann(_highFreqEnergy, _energy, _xResSm, _yResSm, _zResSm);
469   downsampleYNeumann(_tcTemp, _highFreqEnergy, _xResSm, _yResSm, _zResSm);
470   downsampleZNeumann(_highFreqEnergy, _tcTemp, _xResSm, _yResSm, _zResSm);
471
472   // upsample input
473   upsampleZNeumann(_tcTemp, _highFreqEnergy, _xResSm, _yResSm, _zResSm);
474   upsampleYNeumann(_highFreqEnergy, _tcTemp, _xResSm, _yResSm, _zResSm);
475   upsampleXNeumann(_tcTemp, _highFreqEnergy, _xResSm, _yResSm, _zResSm);
476
477   // subtract the down and upsampled field from the original field -- 
478   // what should be left over is solely the high frequency component
479         int index = 0;
480   for (int z = 0; z < _zResSm; z++) 
481     for (int y = 0; y < _yResSm; y++) {
482       for (int x = 0; x < _xResSm; x++, index++) {
483         // brute force reset of boundaries
484         if(z >= _zResSm - 1 || x >= _xResSm - 1 || y >= _yResSm - 1 || z <= 0 || y <= 0 || x <= 0) 
485           _highFreqEnergy[index] = 0.; 
486         else 
487           _highFreqEnergy[index] = _energy[index] - _tcTemp[index];
488     }
489   }
490 }
491
492 //////////////////////////////////////////////////////////////////////
493 // compute velocity from energies and march into obstacles
494 // for wavelet decomposition
495 ////////////////////////////////////////////////////////////////////// 
496 void WTURBULENCE::computeEnergy(float *_energy, float* xvel, float* yvel, float* zvel, unsigned char *origObstacles) 
497 {
498   unsigned char *obstacles = new unsigned char[_totalCellsSm];
499   memcpy(obstacles, origObstacles, sizeof(unsigned char) * _totalCellsSm);
500
501   // compute everywhere
502   for (int x = 0; x < _totalCellsSm; x++) 
503     _energy[x] = 0.5f * (xvel[x] * xvel[x] + yvel[x] * yvel[x] + zvel[x] * zvel[x]);
504
505   FLUID_3D::copyBorderX(_energy, _resSm, 0 , _resSm[2]);
506   FLUID_3D::copyBorderY(_energy, _resSm, 0 , _resSm[2]);
507   FLUID_3D::copyBorderZ(_energy, _resSm, 0 , _resSm[2]);
508
509   // pseudo-march the values into the obstacles
510   // the wavelet upsampler only uses a 3x3 support neighborhood, so
511   // propagating the values in by 4 should be sufficient
512   int index;
513
514   // iterate
515   for (int iter = 0; iter < 4; iter++)
516   {
517     index = _slabSizeSm + _xResSm + 1;
518     for (int z = 1; z < _zResSm - 1; z++, index += 2 * _xResSm)
519       for (int y = 1; y < _yResSm - 1; y++, index += 2)
520         for (int x = 1; x < _xResSm - 1; x++, index++)
521           if (obstacles[index] && obstacles[index] != RETIRED)
522           {
523             float sum = 0.0f;
524             int valid = 0;
525
526             if (!obstacles[index + 1] || obstacles[index + 1] == RETIRED)
527             {
528               sum += _energy[index + 1];
529               valid++;
530             }
531             if (!obstacles[index - 1] || obstacles[index - 1] == RETIRED)
532             {
533               sum += _energy[index - 1];
534               valid++;
535             }
536             if (!obstacles[index + _xResSm] || obstacles[index + _xResSm] == RETIRED)
537             {
538               sum += _energy[index + _xResSm];
539               valid++;
540             }
541             if (!obstacles[index - _xResSm] || obstacles[index - _xResSm] == RETIRED)
542             {
543               sum += _energy[index - _xResSm];
544               valid++;
545             }
546             if (!obstacles[index + _slabSizeSm] || obstacles[index + _slabSizeSm] == RETIRED)
547             {
548               sum += _energy[index + _slabSizeSm];
549               valid++;
550             }
551             if (!obstacles[index - _slabSizeSm] || obstacles[index - _slabSizeSm] == RETIRED)
552             {
553               sum += _energy[index - _slabSizeSm];
554               valid++;
555             }
556             if (valid > 0)
557             {
558               _energy[index] = sum / (float)valid;
559               obstacles[index] = MARCHED;
560             }
561           }
562     index = _slabSizeSm + _xResSm + 1;
563     for (int z = 1; z < _zResSm - 1; z++, index += 2 * _xResSm)
564       for (int y = 1; y < _yResSm - 1; y++, index += 2)
565         for (int x = 1; x < _xResSm - 1; x++, index++)
566           if (obstacles[index] == MARCHED)
567             obstacles[index] = RETIRED;
568   }
569   index = _slabSizeSm + _xResSm + 1;
570   for (int z = 1; z < _zResSm - 1; z++, index += 2 * _xResSm)
571     for (int y = 1; y < _yResSm - 1; y++, index += 2)
572       for (int x = 1; x < _xResSm - 1; x++, index++)
573         if (obstacles[index])
574           obstacles[index] = 1; // DG TODO ? animated obstacle flag?
575
576   delete [] obstacles;
577 }
578
579 //////////////////////////////////////////////////////////////////////////////////////////
580 // Evaluate derivatives
581 //////////////////////////////////////////////////////////////////////////////////////////
582 Vec3 WTURBULENCE::WVelocity(Vec3 orgPos)
583 {
584   // arbitrarily offset evaluation points
585   const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0);
586   const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0);
587   const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0);
588
589   const float f1y = WNoiseDy(p1, _noiseTile);
590   const float f1z = WNoiseDz(p1, _noiseTile);
591
592   const float f2x = WNoiseDx(p2, _noiseTile);
593   const float f2z = WNoiseDz(p2, _noiseTile);
594
595   const float f3x = WNoiseDx(p3, _noiseTile);
596   const float f3y = WNoiseDy(p3, _noiseTile);
597
598   Vec3 ret = Vec3( 
599       f3y - f2z,
600       f1z - f3x,
601       f2x - f1y ); 
602   return ret;
603 }
604
605 //////////////////////////////////////////////////////////////////////////////////////////
606 // Evaluate derivatives with Jacobian
607 //////////////////////////////////////////////////////////////////////////////////////////
608 Vec3 WTURBULENCE::WVelocityWithJacobian(Vec3 orgPos, float* xUnwarped, float* yUnwarped, float* zUnwarped)
609 {
610   // arbitrarily offset evaluation points
611   const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0);
612   const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0);
613   const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0);
614
615   Vec3 final;
616   final[0] = WNoiseDx(p1, _noiseTile);
617   final[1] = WNoiseDy(p1, _noiseTile);
618   final[2] = WNoiseDz(p1, _noiseTile);
619   // UNUSED const float f1x = xUnwarped[0] * final[0] + xUnwarped[1] * final[1] + xUnwarped[2] * final[2];
620   const float f1y = yUnwarped[0] * final[0] + yUnwarped[1] * final[1] + yUnwarped[2] * final[2];
621   const float f1z = zUnwarped[0] * final[0] + zUnwarped[1] * final[1] + zUnwarped[2] * final[2];
622
623   final[0] = WNoiseDx(p2, _noiseTile);
624   final[1] = WNoiseDy(p2, _noiseTile);
625   final[2] = WNoiseDz(p2, _noiseTile);
626   const float f2x = xUnwarped[0] * final[0] + xUnwarped[1] * final[1] + xUnwarped[2] * final[2];
627   // UNUSED const float f2y = yUnwarped[0] * final[0] + yUnwarped[1] * final[1] + yUnwarped[2] * final[2];
628   const float f2z = zUnwarped[0] * final[0] + zUnwarped[1] * final[1] + zUnwarped[2] * final[2];
629
630   final[0] = WNoiseDx(p3, _noiseTile);
631   final[1] = WNoiseDy(p3, _noiseTile);
632   final[2] = WNoiseDz(p3, _noiseTile);
633   const float f3x = xUnwarped[0] * final[0] + xUnwarped[1] * final[1] + xUnwarped[2] * final[2];
634   const float f3y = yUnwarped[0] * final[0] + yUnwarped[1] * final[1] + yUnwarped[2] * final[2];
635   // UNUSED const float f3z = zUnwarped[0] * final[0] + zUnwarped[1] * final[1] + zUnwarped[2] * final[2];
636
637   Vec3 ret = Vec3( 
638       f3y - f2z,
639       f1z - f3x,
640       f2x - f1y ); 
641   return ret;
642 }
643
644
645 //////////////////////////////////////////////////////////////////////
646 // perform an actual noise advection step
647 //////////////////////////////////////////////////////////////////////
648 /*void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) 
649 {
650         // enlarge timestep to match grid
651         const float dt = dtOrg * _amplify;
652         const float invAmp = 1.0f / _amplify;
653         float *tempBig1 = new float[_totalCellsBig];
654         float *tempBig2 = new float[_totalCellsBig];
655         float *bigUx = new float[_totalCellsBig];
656         float *bigUy = new float[_totalCellsBig];
657         float *bigUz = new float[_totalCellsBig]; 
658         float *_energy = new float[_totalCellsSm];
659         float *highFreqEnergy = new float[_totalCellsSm];
660         float *eigMin  = new float[_totalCellsSm];
661         float *eigMax  = new float[_totalCellsSm];
662
663         memset(tempBig1, 0, sizeof(float)*_totalCellsBig);
664         memset(tempBig2, 0, sizeof(float)*_totalCellsBig);
665         memset(highFreqEnergy, 0, sizeof(float)*_totalCellsSm);
666         memset(eigMin, 0, sizeof(float)*_totalCellsSm);
667         memset(eigMax, 0, sizeof(float)*_totalCellsSm);
668
669         // prepare textures
670         advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2);
671
672   // compute eigenvalues of the texture coordinates
673   computeEigenvalues(eigMin, eigMax);
674
675   // do wavelet decomposition of energy
676   computeEnergy(_energy, xvel, yvel, zvel, obstacles);
677   decomposeEnergy(_energy, highFreqEnergy);
678
679   // zero out coefficients inside of the obstacle
680   for (int x = 0; x < _totalCellsSm; x++)
681     if (obstacles[x]) _energy[x] = 0.f;
682
683   float maxVelocity = 0.;
684   for (int z = 1; z < _zResBig - 1; z++) 
685     for (int y = 1; y < _yResBig - 1; y++) 
686       for (int x = 1; x < _xResBig - 1; x++)
687       {
688         // get unit position for both fine and coarse grid
689         const Vec3 pos = Vec3(x,y,z);
690         const Vec3 posSm = pos * invAmp;
691         
692         // get grid index for both fine and coarse grid
693         const int index = x + y *_xResBig + z *_slabSizeBig;
694         const int indexSmall = (int)posSm[0] + (int)posSm[1] * _xResSm + (int)posSm[2] * _slabSizeSm;
695         
696         // get a linearly interpolated velocity and texcoords
697         // from the coarse grid
698         Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel, 
699             posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
700         Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW, 
701             posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
702
703         // multiply the texture coordinate by _resSm so that turbulence
704         // synthesis begins at the first octave that the coarse grid 
705         // cannot capture
706         Vec3 texCoord = Vec3(uvw[0] * _resSm[0], 
707                              uvw[1] * _resSm[1],
708                              uvw[2] * _resSm[2]); 
709
710         // retrieve wavelet energy at highest frequency
711         float energy = INTERPOLATE::lerp3d(
712             highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
713
714         // base amplitude for octave 0
715         float coefficient = sqrtf(2.0f * fabs(energy));
716         const float amplitude = *_strength * fabs(0.5 * coefficient) * persistence;
717
718         // add noise to velocity, but only if the turbulence is
719         // sufficiently undeformed, and the energy is large enough
720         // to make a difference
721         const bool addNoise = eigMax[indexSmall] < 2. && 
722                               eigMin[indexSmall] > 0.5;
723         if (addNoise && amplitude > _cullingThreshold) {
724           // base amplitude for octave 0
725           float amplitudeScaled = amplitude;
726
727           for (int octave = 0; octave < _octaves; octave++)
728           {
729             // multiply the vector noise times the maximum allowed
730             // noise amplitude at this octave, and add it to the total
731             vel += WVelocity(texCoord) * amplitudeScaled;
732
733             // scale coefficient for next octave
734             amplitudeScaled *= persistence;
735             texCoord *= 2.0f;
736           }
737         }
738
739         // Store velocity + turbulence in big grid for maccormack step
740         //
741         // If you wanted to save memory, you would instead perform a 
742         // semi-Lagrangian backtrace for the current grid cell here. Then
743         // you could just throw the velocity away.
744         bigUx[index] = vel[0];
745         bigUy[index] = vel[1];
746         bigUz[index] = vel[2];
747
748         // compute the velocity magnitude for substepping later
749         const float velMag = bigUx[index] * bigUx[index] + 
750                              bigUy[index] * bigUy[index] + 
751                              bigUz[index] * bigUz[index];
752         if (velMag > maxVelocity) maxVelocity = velMag;
753
754         // zero out velocity inside obstacles
755         float obsCheck = INTERPOLATE::lerp3dToFloat(
756             obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); 
757         if (obsCheck > 0.95) 
758           bigUx[index] = bigUy[index] = bigUz[index] = 0.;
759       }
760
761   // prepare density for an advection
762   SWAP_POINTERS(_densityBig, _densityBigOld);
763
764   // based on the maximum velocity present, see if we need to substep,
765   // but cap the maximum number of substeps to 5
766   const int maxSubSteps = 5; 
767   maxVelocity = sqrt(maxVelocity) * dt;
768   int totalSubsteps = (int)(maxVelocity / (float)maxSubSteps);
769   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
770   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
771   const float dtSubdiv = dt / (float)totalSubsteps;
772
773   // set boundaries of big velocity grid
774   FLUID_3D::setZeroX(bigUx, _resBig, 0, _resBig[2]); 
775   FLUID_3D::setZeroY(bigUy, _resBig, 0, _resBig[2]); 
776   FLUID_3D::setZeroZ(bigUz, _resBig, 0, _resBig[2]);
777
778   // do the MacCormack advection, with substepping if necessary
779   for(int substep = 0; substep < totalSubsteps; substep++)
780   {
781     FLUID_3D::advectFieldMacCormack(dtSubdiv, bigUx, bigUy, bigUz, 
782         _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL);
783
784     if (substep < totalSubsteps - 1) 
785       SWAP_POINTERS(_densityBig, _densityBigOld);
786   } // substep
787   
788   // wipe the density borders
789   FLUID_3D::setZeroBorder(_densityBig, _resBig, 0, _resBig[2]);
790     
791   // reset texture coordinates now in preparation for next timestep
792   // Shouldn't do this before generating the noise because then the 
793   // eigenvalues stored do not reflect the underlying texture coordinates
794   resetTextureCoordinates(eigMin, eigMax);
795   
796   delete[] tempBig1;
797   delete[] tempBig2;
798   delete[] bigUx;
799   delete[] bigUy;
800   delete[] bigUz;
801   delete[] _energy;
802   delete[] highFreqEnergy;
803
804   delete[] eigMin;
805   delete[] eigMax;
806   
807
808   _totalStepsBig++;
809 }*/
810
811 //struct
812
813 //////////////////////////////////////////////////////////////////////
814 // perform the full turbulence algorithm, including OpenMP 
815 // if available
816 //////////////////////////////////////////////////////////////////////
817 void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles)
818 {
819         // enlarge timestep to match grid
820         const float dt = dtOrg * _amplify;
821         const float invAmp = 1.0f / _amplify;
822         float *tempFuelBig = NULL, *tempReactBig = NULL;
823         float *tempColor_rBig = NULL, *tempColor_gBig = NULL, *tempColor_bBig = NULL;
824         float *tempDensityBig = (float *)calloc(_totalCellsBig, sizeof(float));
825         float *tempBig = (float *)calloc(_totalCellsBig, sizeof(float));
826         float *bigUx = (float *)calloc(_totalCellsBig, sizeof(float));
827         float *bigUy = (float *)calloc(_totalCellsBig, sizeof(float));
828         float *bigUz = (float *)calloc(_totalCellsBig, sizeof(float)); 
829         float *_energy = (float *)calloc(_totalCellsSm, sizeof(float));
830         float *highFreqEnergy = (float *)calloc(_totalCellsSm, sizeof(float));
831         float *eigMin  = (float *)calloc(_totalCellsSm, sizeof(float));
832         float *eigMax  = (float *)calloc(_totalCellsSm, sizeof(float));
833
834         if (_fuelBig) {
835                 tempFuelBig = (float *)calloc(_totalCellsBig, sizeof(float));
836                 tempReactBig = (float *)calloc(_totalCellsBig, sizeof(float));
837         }
838         if (_color_rBig) {
839                 tempColor_rBig = (float *)calloc(_totalCellsBig, sizeof(float));
840                 tempColor_gBig = (float *)calloc(_totalCellsBig, sizeof(float));
841                 tempColor_bBig = (float *)calloc(_totalCellsBig, sizeof(float));
842         }
843
844         memset(_tcTemp, 0, sizeof(float)*_totalCellsSm);
845
846
847         // prepare textures
848         advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempDensityBig, tempBig);
849
850         // do wavelet decomposition of energy
851         computeEnergy(_energy, xvel, yvel, zvel, obstacles);
852
853         for (int x = 0; x < _totalCellsSm; x++)
854                 if (obstacles[x]) _energy[x] = 0.f;
855
856         decomposeEnergy(_energy, highFreqEnergy);
857
858         // zero out coefficients inside of the obstacle
859         for (int x = 0; x < _totalCellsSm; x++)
860                 if (obstacles[x]) highFreqEnergy[x] = 0.f;
861
862         Vec3Int ressm(_xResSm, _yResSm, _zResSm);
863         FLUID_3D::setNeumannX(highFreqEnergy, ressm, 0 , ressm[2]);
864         FLUID_3D::setNeumannY(highFreqEnergy, ressm, 0 , ressm[2]);
865         FLUID_3D::setNeumannZ(highFreqEnergy, ressm, 0 , ressm[2]);
866
867
868    int threadval = 1;
869 #if PARALLEL==1
870   threadval = omp_get_max_threads();
871 #endif
872
873
874   // parallel region setup
875   // Uses omp_get_max_trheads to get number of required cells.
876   float* maxVelMagThreads = new float[threadval];
877
878   for (int i=0; i<threadval; i++) maxVelMagThreads[i] = -1.0f;
879
880 #if PARALLEL==1
881
882 #pragma omp parallel
883 #endif
884   { float maxVelMag1 = 0.;
885 #if PARALLEL==1
886     const int id  = omp_get_thread_num(); /*, num = omp_get_num_threads(); */
887 #endif
888
889   // vector noise main loop
890 #if PARALLEL==1
891 #pragma omp for schedule(static,1)
892 #endif
893   for (int zSmall = 0; zSmall < _zResSm; zSmall++)
894   {
895   for (int ySmall = 0; ySmall < _yResSm; ySmall++) 
896   for (int xSmall = 0; xSmall < _xResSm; xSmall++)
897   {
898     const int indexSmall = xSmall + ySmall * _xResSm + zSmall * _slabSizeSm;
899
900     // compute jacobian
901     float jacobian[3][3] = {
902       { minDx(xSmall, ySmall, zSmall, _tcU, _resSm), minDx(xSmall, ySmall, zSmall, _tcV, _resSm), minDx(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
903       { minDy(xSmall, ySmall, zSmall, _tcU, _resSm), minDy(xSmall, ySmall, zSmall, _tcV, _resSm), minDy(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
904       { minDz(xSmall, ySmall, zSmall, _tcU, _resSm), minDz(xSmall, ySmall, zSmall, _tcV, _resSm), minDz(xSmall, ySmall, zSmall, _tcW, _resSm) }
905     };
906
907     // get LU factorization of texture jacobian and apply 
908     // it to unit vectors
909     sLU LU = computeLU(jacobian);
910     float xUnwarped[3], yUnwarped[3], zUnwarped[3];
911     float xWarped[3], yWarped[3], zWarped[3];
912     bool nonSingular = isNonsingular(LU);
913
914         xUnwarped[0] = 1.0f; xUnwarped[1] = 0.0f; xUnwarped[2] = 0.0f;
915         yUnwarped[0] = 0.0f; yUnwarped[1] = 1.0f; yUnwarped[2] = 0.0f;
916         zUnwarped[0] = 0.0f; zUnwarped[1] = 0.0f; zUnwarped[2] = 1.0f;
917
918         xWarped[0] = 1.0f; xWarped[1] = 0.0f; xWarped[2] = 0.0f;
919         yWarped[0] = 0.0f; yWarped[1] = 1.0f; yWarped[2] = 0.0f;
920         zWarped[0] = 0.0f; zWarped[1] = 0.0f; zWarped[2] = 1.0f;
921
922 #if 0
923         // UNUSED
924     float eigMax = 10.0f;
925     float eigMin = 0.1f;
926 #endif
927     if (nonSingular)
928     {
929       solveLU3x3(LU, xUnwarped, xWarped);
930       solveLU3x3(LU, yUnwarped, yWarped);
931       solveLU3x3(LU, zUnwarped, zWarped);
932
933       // compute the eigenvalues while we have the Jacobian available
934       Vec3 eigenvalues = Vec3(1.);
935       computeEigenvalues3x3( &eigenvalues[0], jacobian);
936       eigMax[indexSmall] = MAX3V(eigenvalues);
937       eigMin[indexSmall] = MIN3V(eigenvalues);
938     }
939     
940     // make sure to skip one on the beginning and end
941     int xStart = (xSmall == 0) ? 1 : 0;
942     int xEnd   = (xSmall == _xResSm - 1) ? _amplify - 1 : _amplify;
943     int yStart = (ySmall == 0) ? 1 : 0;
944     int yEnd   = (ySmall == _yResSm - 1) ? _amplify - 1 : _amplify;
945     int zStart = (zSmall == 0) ? 1 : 0;
946     int zEnd   = (zSmall == _zResSm - 1) ? _amplify - 1 : _amplify;
947       
948     for (int zBig = zStart; zBig < zEnd; zBig++) 
949     for (int yBig = yStart; yBig < yEnd; yBig++) 
950     for (int xBig = xStart; xBig < xEnd; xBig++)
951     {
952       const int x = xSmall * _amplify + xBig;
953       const int y = ySmall * _amplify + yBig;
954       const int z = zSmall * _amplify + zBig;
955       
956       // get unit position for both fine and coarse grid
957       const Vec3 pos = Vec3(x,y,z);
958       const Vec3 posSm = pos * invAmp;
959       
960       // get grid index for both fine and coarse grid
961       const int index = x + y *_xResBig + z *_slabSizeBig;
962       
963       // get a linearly interpolated velocity and texcoords
964       // from the coarse grid
965       Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel, 
966           posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
967       Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW, 
968           posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
969
970       // multiply the texture coordinate by _resSm so that turbulence
971       // synthesis begins at the first octave that the coarse grid 
972       // cannot capture
973       Vec3 texCoord = Vec3(uvw[0] * _resSm[0], 
974                            uvw[1] * _resSm[1],
975                            uvw[2] * _resSm[2]); 
976
977       // retrieve wavelet energy at highest frequency
978       float energy = INTERPOLATE::lerp3d(
979           highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
980
981       // base amplitude for octave 0
982       float coefficient = sqrtf(2.0f * fabs(energy));
983       const float amplitude = *_strength * fabs(0.5f * coefficient) * persistence;
984
985       // add noise to velocity, but only if the turbulence is
986       // sufficiently undeformed, and the energy is large enough
987       // to make a difference
988       const bool addNoise = eigMax[indexSmall] < 2.0f &&
989                             eigMin[indexSmall] > 0.5f;
990       if (addNoise && amplitude > _cullingThreshold) {
991         // base amplitude for octave 0
992         float amplitudeScaled = amplitude;
993
994         for (int octave = 0; octave < _octaves; octave++)
995         {
996           // multiply the vector noise times the maximum allowed
997           // noise amplitude at this octave, and add it to the total
998           vel += WVelocityWithJacobian(texCoord, &xUnwarped[0], &yUnwarped[0], &zUnwarped[0]) * amplitudeScaled;
999
1000           // scale coefficient for next octave
1001           amplitudeScaled *= persistence;
1002           texCoord *= 2.0f;
1003         }
1004       }
1005
1006       // Store velocity + turbulence in big grid for maccormack step
1007       //
1008       // If you wanted to save memory, you would instead perform a 
1009       // semi-Lagrangian backtrace for the current grid cell here. Then
1010       // you could just throw the velocity away.
1011       bigUx[index] = vel[0];
1012       bigUy[index] = vel[1];
1013       bigUz[index] = vel[2];
1014
1015       // compute the velocity magnitude for substepping later
1016       const float velMag = bigUx[index] * bigUx[index] + 
1017                            bigUy[index] * bigUy[index] + 
1018                            bigUz[index] * bigUz[index];
1019       if (velMag > maxVelMag1) maxVelMag1 = velMag;
1020
1021       // zero out velocity inside obstacles
1022       float obsCheck = INTERPOLATE::lerp3dToFloat(
1023           obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); 
1024       if (obsCheck > 0.95f)
1025         bigUx[index] = bigUy[index] = bigUz[index] = 0.;
1026     } // xyz*/
1027
1028 #if PARALLEL==1
1029     maxVelMagThreads[id] = maxVelMag1;
1030 #else
1031     maxVelMagThreads[0] = maxVelMag1;
1032 #endif
1033   }
1034   }
1035   } // omp
1036   
1037   // compute maximum over threads
1038   float maxVelMag = maxVelMagThreads[0];
1039 #if PARALLEL==1
1040   for (int i = 1; i < threadval; i++) 
1041     if (maxVelMag < maxVelMagThreads[i]) 
1042       maxVelMag = maxVelMagThreads[i];
1043 #endif
1044   delete [] maxVelMagThreads;
1045
1046
1047   // prepare density for an advection
1048   SWAP_POINTERS(_densityBig, _densityBigOld);
1049   SWAP_POINTERS(_fuelBig, _fuelBigOld);
1050   SWAP_POINTERS(_reactBig, _reactBigOld);
1051   SWAP_POINTERS(_color_rBig, _color_rBigOld);
1052   SWAP_POINTERS(_color_gBig, _color_gBigOld);
1053   SWAP_POINTERS(_color_bBig, _color_bBigOld);
1054
1055   // based on the maximum velocity present, see if we need to substep,
1056   // but cap the maximum number of substeps to 5
1057   const int maxSubSteps = 25;
1058   const int maxVel = 5;
1059   maxVelMag = sqrt(maxVelMag) * dt;
1060   int totalSubsteps = (int)(maxVelMag / (float)maxVel);
1061   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
1062   // printf("totalSubsteps: %d\n", totalSubsteps);
1063   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
1064   const float dtSubdiv = dt / (float)totalSubsteps;
1065
1066   // set boundaries of big velocity grid
1067   FLUID_3D::setZeroX(bigUx, _resBig, 0 , _resBig[2]); 
1068   FLUID_3D::setZeroY(bigUy, _resBig, 0 , _resBig[2]); 
1069   FLUID_3D::setZeroZ(bigUz, _resBig, 0 , _resBig[2]);
1070
1071 #if PARALLEL==1
1072   int stepParts = threadval*2;  // Dividing parallelized sections into numOfThreads * 2 sections
1073   float partSize = (float)_zResBig/stepParts;   // Size of one part;
1074
1075   if (partSize < 4) {stepParts = threadval;                                     // If the slice gets too low (might actually slow things down, change it to larger
1076                                         partSize = (float)_zResBig/stepParts;}
1077   if (partSize < 4) {stepParts = (int)(ceil((float)_zResBig/4.0f));     // If it's still too low (only possible on future systems with +24 cores), change it to 4
1078                                         partSize = (float)_zResBig/stepParts;}
1079 #else
1080   int zBegin=0;
1081   int zEnd=_resBig[2];
1082 #endif
1083
1084   // do the MacCormack advection, with substepping if necessary
1085   for(int substep = 0; substep < totalSubsteps; substep++)
1086   {
1087
1088 #if PARALLEL==1
1089         #pragma omp parallel
1090         {
1091
1092         #pragma omp for schedule(static,1)
1093         for (int i=0; i<stepParts; i++)
1094         {
1095                 int zBegin = (int)((float)i*partSize + 0.5f);
1096                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
1097 #endif
1098                 FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
1099                     _densityBigOld, tempDensityBig, _resBig, zBegin, zEnd);
1100                 if (_fuelBig) {
1101                         FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
1102                                 _fuelBigOld, tempFuelBig, _resBig, zBegin, zEnd);
1103                         FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
1104                                 _reactBigOld, tempReactBig, _resBig, zBegin, zEnd);
1105                 }
1106                 if (_color_rBig) {
1107                         FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
1108                                 _color_rBigOld, tempColor_rBig, _resBig, zBegin, zEnd);
1109                         FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
1110                                 _color_gBigOld, tempColor_gBig, _resBig, zBegin, zEnd);
1111                         FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
1112                                 _color_bBigOld, tempColor_bBig, _resBig, zBegin, zEnd);
1113                 }
1114 #if PARALLEL==1
1115         }
1116
1117         #pragma omp barrier
1118
1119         #pragma omp for schedule(static,1)
1120         for (int i=0; i<stepParts; i++)
1121         {
1122                 int zBegin = (int)((float)i*partSize + 0.5f);
1123                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
1124 #endif
1125                 FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
1126                     _densityBigOld, _densityBig, tempDensityBig, tempBig, _resBig, NULL, zBegin, zEnd);
1127                 if (_fuelBig) {
1128                         FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
1129                                 _fuelBigOld, _fuelBig, tempFuelBig, tempBig, _resBig, NULL, zBegin, zEnd);
1130                         FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
1131                                 _reactBigOld, _reactBig, tempReactBig, tempBig, _resBig, NULL, zBegin, zEnd);
1132                 }
1133                 if (_color_rBig) {
1134                         FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
1135                                 _color_rBigOld, _color_rBig, tempColor_rBig, tempBig, _resBig, NULL, zBegin, zEnd);
1136                         FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
1137                                 _color_gBigOld, _color_gBig, tempColor_gBig, tempBig, _resBig, NULL, zBegin, zEnd);
1138                         FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
1139                                 _color_bBigOld, _color_bBig, tempColor_bBig, tempBig, _resBig, NULL, zBegin, zEnd);
1140                 }
1141 #if PARALLEL==1
1142         }
1143         }
1144 #endif
1145
1146         if (substep < totalSubsteps - 1) {
1147       SWAP_POINTERS(_densityBig, _densityBigOld);
1148           SWAP_POINTERS(_fuelBig, _fuelBigOld);
1149           SWAP_POINTERS(_reactBig, _reactBigOld);
1150           SWAP_POINTERS(_color_rBig, _color_rBigOld);
1151           SWAP_POINTERS(_color_gBig, _color_gBigOld);
1152           SWAP_POINTERS(_color_bBig, _color_bBigOld);
1153         }
1154   } // substep
1155
1156   free(tempDensityBig);
1157   if (tempFuelBig) free(tempFuelBig);
1158   if (tempReactBig) free(tempReactBig);
1159   if (tempColor_rBig) free(tempColor_rBig);
1160   if (tempColor_gBig) free(tempColor_gBig);
1161   if (tempColor_bBig) free(tempColor_bBig);
1162   free(tempBig);
1163   free(bigUx);
1164   free(bigUy);
1165   free(bigUz);
1166   free(_energy);
1167   free(highFreqEnergy);
1168   
1169   // wipe the density borders
1170   FLUID_3D::setZeroBorder(_densityBig, _resBig, 0 , _resBig[2]);
1171   if (_fuelBig) {
1172         FLUID_3D::setZeroBorder(_fuelBig, _resBig, 0 , _resBig[2]);
1173         FLUID_3D::setZeroBorder(_reactBig, _resBig, 0 , _resBig[2]);
1174   }
1175   if (_color_rBig) {
1176           FLUID_3D::setZeroBorder(_color_rBig, _resBig, 0 , _resBig[2]);
1177           FLUID_3D::setZeroBorder(_color_gBig, _resBig, 0 , _resBig[2]);
1178           FLUID_3D::setZeroBorder(_color_bBig, _resBig, 0 , _resBig[2]);
1179   }
1180     
1181   // reset texture coordinates now in preparation for next timestep
1182   // Shouldn't do this before generating the noise because then the 
1183   // eigenvalues stored do not reflect the underlying texture coordinates
1184   resetTextureCoordinates(eigMin, eigMax);
1185
1186   free(eigMin);
1187   free(eigMax);
1188   
1189   // output files
1190   // string prefix = string("./amplified.preview/density_bigxy_");
1191   // FLUID_3D::writeImageSliceXY(_densityBig, _resBig, _resBig[2]/2, prefix, _totalStepsBig, 1.0f);
1192   //string df3prefix = string("./df3/density_big_");
1193   //IMAGE::dumpDF3(_totalStepsBig, df3prefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
1194   // string pbrtPrefix = string("./pbrt/density_big_");
1195   // IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
1196   
1197   _totalStepsBig++;
1198 }