fix some complier warnings and add -Wundef to CMake's default GCC warnings.
[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)
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.5); // 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./(float)_resBig[0], 1./(float)_resBig[1], 1./(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./(float)_resSm[0], 1./(float)_resSm[1], 1./(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         // allocate & init texture coordinates
98         _tcU = new float[_totalCellsSm];
99         _tcV = new float[_totalCellsSm];
100         _tcW = new float[_totalCellsSm];
101         _tcTemp = new float[_totalCellsSm];
102         
103         // map all 
104         const float dx = 1./(float)(_resSm[0]);
105         const float dy = 1./(float)(_resSm[1]);
106         const float dz = 1./(float)(_resSm[2]);
107         int index = 0;
108         for (int z = 0; z < _zResSm; z++) 
109         for (int y = 0; y < _yResSm; y++) 
110                 for (int x = 0; x < _xResSm; x++, index++)
111                 {
112                 _tcU[index] = x*dx;
113                 _tcV[index] = y*dy;
114                 _tcW[index] = z*dz;
115                 _tcTemp[index] = 0.;
116                 }
117         
118         // noise tiles
119         _noiseTile = new float[noiseTileSize * noiseTileSize * noiseTileSize];
120         /*
121         std::string noiseTileFilename = std::string("noise.wavelets");
122         generateTile_WAVELET(_noiseTile, noiseTileFilename);
123         */
124         setNoise(noisetype);
125         /*
126         std::string noiseTileFilename = std::string("noise.fft");
127         generatTile_FFT(_noiseTile, noiseTileFilename);
128         */
129 }
130
131 //////////////////////////////////////////////////////////////////////
132 // destructor
133 //////////////////////////////////////////////////////////////////////
134 WTURBULENCE::~WTURBULENCE() {
135   delete[] _densityBig;
136   delete[] _densityBigOld;
137
138   delete[] _tcU;
139   delete[] _tcV;
140   delete[] _tcW;
141   delete[] _tcTemp;
142
143   delete[] _noiseTile;
144 }
145
146 //////////////////////////////////////////////////////////////////////
147 // Change noise type
148 //
149 // type (1<<0) = wavelet / 2
150 // type (1<<1) = FFT / 4
151 // type (1<<2) = curl / 8
152 //////////////////////////////////////////////////////////////////////
153 void WTURBULENCE::setNoise(int type)
154 {
155         if(type == (1<<1)) // FFT
156         {
157                 // needs fft
158                 #if FFTW3==1
159                 std::string noiseTileFilename = std::string("noise.fft");
160                 generatTile_FFT(_noiseTile, noiseTileFilename);
161                 #endif
162         }
163         else if(type == (1<<2)) // curl
164         {
165                 // TODO: not supported yet
166         }
167         else // standard - wavelet
168         {
169                 std::string noiseTileFilename = std::string("noise.wavelets");
170                 generateTile_WAVELET(_noiseTile, noiseTileFilename);
171         }
172 }
173
174 // init direct access functions from blender
175 void WTURBULENCE::initBlenderRNA(float *strength)
176 {
177         _strength = strength;
178 }
179
180 //////////////////////////////////////////////////////////////////////
181 // Get the smallest valid x derivative
182 //
183 // Takes the one-sided finite difference in both directions and
184 // selects the smaller of the two
185 //////////////////////////////////////////////////////////////////////
186 static float minDx(int x, int y, int z, float* input, Vec3Int res)
187 {
188   const int index = x + y * res[0] + z * res[0] * res[1];
189   const int maxx = res[0]-2;
190
191   // get grid values
192   float center = input[index];
193   float left  = (x <= 1)    ? FLT_MAX : input[index - 1];
194   float right = (x >= maxx) ? FLT_MAX : input[index + 1];
195
196   const float dx = res[0];
197
198   // get all the derivative estimates
199   float dLeft   = (x <= 1)     ? FLT_MAX : (center - left) * dx;
200   float dRight  = (x >= maxx)  ? FLT_MAX : (right - center) * dx;
201   float dCenter = (x <= 1 || x >= maxx) ? FLT_MAX : (right - left) * dx * 0.5f;
202
203   // if it's on a boundary, only one estimate is valid
204   if (x <= 1) return dRight;
205   if (x >= maxx) return dLeft;
206
207   // if it's not on a boundary, get the smallest one
208   float finalD;
209   finalD = (fabs(dCenter) < fabs(dRight)) ? dCenter : dRight;
210   finalD = (fabs(finalD)  < fabs(dLeft))  ? finalD  : dLeft;
211
212   return finalD;
213 }
214
215 //////////////////////////////////////////////////////////////////////
216 // get the smallest valid y derivative
217 //
218 // Takes the one-sided finite difference in both directions and
219 // selects the smaller of the two
220 //////////////////////////////////////////////////////////////////////
221 static float minDy(int x, int y, int z, float* input, Vec3Int res)
222 {
223   const int index = x + y * res[0] + z * res[0] * res[1];
224   const int maxy = res[1]-2;
225
226   // get grid values
227   float center = input[index];
228   float down  = (y <= 1) ? FLT_MAX : input[index - res[0]];
229   float up = (y >= maxy) ? FLT_MAX : input[index + res[0]];
230
231   const float dx = res[1]; // only for square domains
232
233   // get all the derivative estimates
234   float dDown   = (y <= 1)  ? FLT_MAX : (center - down) * dx;
235   float dUp  = (y >= maxy)  ? FLT_MAX : (up - center) * dx;
236   float dCenter = (y <= 1 || y >= maxy) ? FLT_MAX : (up - down) * dx * 0.5f;
237
238   // if it's on a boundary, only one estimate is valid
239   if (y <= 1) return dUp;
240   if (y >= maxy) return dDown;
241
242   // if it's not on a boundary, get the smallest one
243   float finalD = (fabs(dCenter) < fabs(dUp)) ? dCenter : dUp;
244   finalD = (fabs(finalD) < fabs(dDown)) ? finalD : dDown;
245
246   return finalD;
247 }
248
249 //////////////////////////////////////////////////////////////////////
250 // get the smallest valid z derivative
251 //
252 // Takes the one-sided finite difference in both directions and
253 // selects the smaller of the two
254 //////////////////////////////////////////////////////////////////////
255 static float minDz(int x, int y, int z, float* input, Vec3Int res)
256 {
257   const int slab = res[0]*res[1];
258   const int index = x + y * res[0] + z * slab;
259   const int maxz = res[2]-2;
260
261   // get grid values
262   float center = input[index];
263   float front  = (z <= 1) ? FLT_MAX : input[index - slab];
264   float back = (z >= maxz) ? FLT_MAX : input[index + slab];
265
266   const float dx = res[2]; // only for square domains
267
268   // get all the derivative estimates
269   float dfront   = (z <= 1)  ? FLT_MAX : (center - front) * dx;
270   float dback  = (z >= maxz)  ? FLT_MAX : (back - center) * dx;
271   float dCenter = (z <= 1 || z >= maxz) ? FLT_MAX : (back - front) * dx * 0.5f;
272
273   // if it's on a boundary, only one estimate is valid
274   if (z <= 1) return dback;
275   if (z >= maxz) return dfront;
276
277   // if it's not on a boundary, get the smallest one
278   float finalD = (fabs(dCenter) < fabs(dback)) ? dCenter : dback;
279   finalD = (fabs(finalD) < fabs(dfront)) ? finalD : dfront;
280
281   return finalD;
282 }
283
284 //////////////////////////////////////////////////////////////////////
285 // handle texture coordinates (advection, reset, eigenvalues), 
286 // Beware -- uses big density maccormack as temporary arrays
287 ////////////////////////////////////////////////////////////////////// 
288 void WTURBULENCE::advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *tempBig1, float *tempBig2) {
289
290   // advection
291   SWAP_POINTERS(_tcTemp, _tcU);
292   FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
293   FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
294   FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
295   FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel, 
296       _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
297   FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel, 
298       _tcTemp, _tcU, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
299
300   SWAP_POINTERS(_tcTemp, _tcV);
301   FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
302   FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
303   FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
304   FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel, 
305       _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
306   FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel, 
307       _tcTemp, _tcV, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
308
309   SWAP_POINTERS(_tcTemp, _tcW);
310   FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
311   FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
312   FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
313   FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel, 
314       _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
315   FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel, 
316       _tcTemp, _tcW, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
317 }
318
319 //////////////////////////////////////////////////////////////////////
320 // Compute the eigenvalues of the advected texture
321 ////////////////////////////////////////////////////////////////////// 
322 void WTURBULENCE::computeEigenvalues(float *_eigMin, float *_eigMax) {
323   // stats
324   float maxeig = -1.;
325   float mineig = 10.;
326
327   // texture coordinate eigenvalues
328   for (int z = 1; z < _zResSm-1; z++) {
329     for (int y = 1; y < _yResSm-1; y++) 
330       for (int x = 1; x < _xResSm-1; x++)
331       {
332         const int index = x+ y *_resSm[0] + z*_slabSizeSm;
333
334         // compute jacobian
335         float jacobian[3][3] = {
336           { minDx(x, y, z, _tcU, _resSm), minDx(x, y, z, _tcV, _resSm), minDx(x, y, z, _tcW, _resSm) } ,
337           { minDy(x, y, z, _tcU, _resSm), minDy(x, y, z, _tcV, _resSm), minDy(x, y, z, _tcW, _resSm) } ,
338           { minDz(x, y, z, _tcU, _resSm), minDz(x, y, z, _tcV, _resSm), minDz(x, y, z, _tcW, _resSm) }
339         };
340
341         // ONLY compute the eigenvalues after checking that the matrix
342         // is nonsingular
343         sLU LU = computeLU(jacobian);
344
345         if (isNonsingular(LU))
346         {
347           // get the analytic eigenvalues, quite slow right now...
348           Vec3 eigenvalues = Vec3(1.);
349           computeEigenvalues3x3( &eigenvalues[0], jacobian);
350           _eigMax[index] = MAX3V(eigenvalues);
351           _eigMin[index] = MIN3V(eigenvalues);
352           maxeig = MAX(_eigMax[index],maxeig);
353           mineig = MIN(_eigMin[index],mineig);
354         }
355         else
356         {
357           _eigMax[index] = 10.0f;
358           _eigMin[index] = 0.1;
359         }
360       }
361   }
362 }
363
364 //////////////////////////////////////////////////////////////////////
365 // advect & reset texture coordinates based on eigenvalues
366 ////////////////////////////////////////////////////////////////////// 
367 void WTURBULENCE::resetTextureCoordinates(float *_eigMin, float *_eigMax) 
368 {
369   // allowed deformation of the textures
370   const float limit = 2.f;
371   const float limitInv = 1./limit;
372
373   // standard reset
374   int resets = 0;
375   const float dx = 1./(float)(_resSm[0]);
376   const float dy = 1./(float)(_resSm[1]);
377   const float dz = 1./(float)(_resSm[2]);
378
379   for (int z = 1; z < _zResSm-1; z++)
380     for (int y = 1; y < _yResSm-1; y++)
381       for (int x = 1; x < _xResSm-1; x++)
382       {
383         const int index = x+ y *_resSm[0] + z*_slabSizeSm;
384         if (_eigMax[index] > limit || _eigMin[index] < limitInv)
385         {
386           _tcU[index] = (float)x * dx;
387           _tcV[index] = (float)y * dy;
388           _tcW[index] = (float)z * dz;
389           resets++;
390         }
391       }
392 }
393
394 //////////////////////////////////////////////////////////////////////
395 // Compute the highest frequency component of the wavelet
396 // decomposition
397 //////////////////////////////////////////////////////////////////////
398 void WTURBULENCE::decomposeEnergy(float *_energy, float *_highFreqEnergy)
399 {
400   // do the decomposition -- the goal here is to have
401   // the energy with the high frequency component stomped out
402   // stored in _tcTemp when it is done. _highFreqEnergy is only used
403   // as an additional temp array
404   
405   // downsample input
406   downsampleXNeumann(_highFreqEnergy, _energy, _xResSm, _yResSm, _zResSm);
407   downsampleYNeumann(_tcTemp, _highFreqEnergy, _xResSm, _yResSm, _zResSm);
408   downsampleZNeumann(_highFreqEnergy, _tcTemp, _xResSm, _yResSm, _zResSm);
409
410   // upsample input
411   upsampleZNeumann(_tcTemp, _highFreqEnergy, _xResSm, _yResSm, _zResSm);
412   upsampleYNeumann(_highFreqEnergy, _tcTemp, _xResSm, _yResSm, _zResSm);
413   upsampleXNeumann(_tcTemp, _highFreqEnergy, _xResSm, _yResSm, _zResSm);
414
415   // subtract the down and upsampled field from the original field -- 
416   // what should be left over is solely the high frequency component
417         int index = 0;
418   for (int z = 0; z < _zResSm; z++) 
419     for (int y = 0; y < _yResSm; y++) {
420       for (int x = 0; x < _xResSm; x++, index++) {
421         // brute force reset of boundaries
422         if(z >= _zResSm - 1 || x >= _xResSm - 1 || y >= _yResSm - 1 || z <= 0 || y <= 0 || x <= 0) 
423           _highFreqEnergy[index] = 0.; 
424         else 
425           _highFreqEnergy[index] = _energy[index] - _tcTemp[index];
426     }
427   }
428 }
429
430 //////////////////////////////////////////////////////////////////////
431 // compute velocity from energies and march into obstacles
432 // for wavelet decomposition
433 ////////////////////////////////////////////////////////////////////// 
434 void WTURBULENCE::computeEnergy(float *_energy, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) 
435 {
436   // compute everywhere
437   for (int x = 0; x < _totalCellsSm; x++) 
438     _energy[x] = 0.5f * (xvel[x] * xvel[x] + yvel[x] * yvel[x] + zvel[x] * zvel[x]);
439
440   FLUID_3D::copyBorderX(_energy, _resSm, 0 , _resSm[2]);
441   FLUID_3D::copyBorderY(_energy, _resSm, 0 , _resSm[2]);
442   FLUID_3D::copyBorderZ(_energy, _resSm, 0 , _resSm[2]);
443
444   // pseudo-march the values into the obstacles
445   // the wavelet upsampler only uses a 3x3 support neighborhood, so
446   // propagating the values in by 4 should be sufficient
447   int index;
448
449   // iterate
450   for (int iter = 0; iter < 4; iter++)
451   {
452     index = _slabSizeSm + _xResSm + 1;
453     for (int z = 1; z < _zResSm - 1; z++, index += 2 * _xResSm)
454       for (int y = 1; y < _yResSm - 1; y++, index += 2)
455         for (int x = 1; x < _xResSm - 1; x++, index++)
456           if (obstacles[index] && obstacles[index] != RETIRED)
457           {
458             float sum = 0.0f;
459             int valid = 0;
460
461             if (!obstacles[index + 1] || obstacles[index + 1] == RETIRED)
462             {
463               sum += _energy[index + 1];
464               valid++;
465             }
466             if (!obstacles[index - 1] || obstacles[index - 1] == RETIRED)
467             {
468               sum += _energy[index - 1];
469               valid++;
470             }
471             if (!obstacles[index + _xResSm] || obstacles[index + _xResSm] == RETIRED)
472             {
473               sum += _energy[index + _xResSm];
474               valid++;
475             }
476             if (!obstacles[index - _xResSm] || obstacles[index - _xResSm] == RETIRED)
477             {
478               sum += _energy[index - _xResSm];
479               valid++;
480             }
481             if (!obstacles[index + _slabSizeSm] || obstacles[index + _slabSizeSm] == RETIRED)
482             {
483               sum += _energy[index + _slabSizeSm];
484               valid++;
485             }
486             if (!obstacles[index - _slabSizeSm] || obstacles[index - _slabSizeSm] == RETIRED)
487             {
488               sum += _energy[index - _slabSizeSm];
489               valid++;
490             }
491             if (valid > 0)
492             {
493               _energy[index] = sum / (float)valid;
494               obstacles[index] = MARCHED;
495             }
496           }
497     index = _slabSizeSm + _xResSm + 1;
498     for (int z = 1; z < _zResSm - 1; z++, index += 2 * _xResSm)
499       for (int y = 1; y < _yResSm - 1; y++, index += 2)
500         for (int x = 1; x < _xResSm - 1; x++, index++)
501           if (obstacles[index] == MARCHED)
502             obstacles[index] = RETIRED;
503   }
504   index = _slabSizeSm + _xResSm + 1;
505   for (int z = 1; z < _zResSm - 1; z++, index += 2 * _xResSm)
506     for (int y = 1; y < _yResSm - 1; y++, index += 2)
507       for (int x = 1; x < _xResSm - 1; x++, index++)
508         if (obstacles[index])
509           obstacles[index] = 1;
510 }
511
512 //////////////////////////////////////////////////////////////////////////////////////////
513 // Evaluate derivatives
514 //////////////////////////////////////////////////////////////////////////////////////////
515 Vec3 WTURBULENCE::WVelocity(Vec3 orgPos)
516 {
517   // arbitrarily offset evaluation points
518   const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0);
519   const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0);
520   const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0);
521
522   const float f1y = WNoiseDy(p1, _noiseTile);
523   const float f1z = WNoiseDz(p1, _noiseTile);
524
525   const float f2x = WNoiseDx(p2, _noiseTile);
526   const float f2z = WNoiseDz(p2, _noiseTile);
527
528   const float f3x = WNoiseDx(p3, _noiseTile);
529   const float f3y = WNoiseDy(p3, _noiseTile);
530
531   Vec3 ret = Vec3( 
532       f3y - f2z,
533       f1z - f3x,
534       f2x - f1y ); 
535   return ret;
536 }
537
538 //////////////////////////////////////////////////////////////////////////////////////////
539 // Evaluate derivatives with Jacobian
540 //////////////////////////////////////////////////////////////////////////////////////////
541 Vec3 WTURBULENCE::WVelocityWithJacobian(Vec3 orgPos, float* xUnwarped, float* yUnwarped, float* zUnwarped)
542 {
543   // arbitrarily offset evaluation points
544   const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0);
545   const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0);
546   const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0);
547
548   Vec3 final;
549   final[0] = WNoiseDx(p1, _noiseTile);
550   final[1] = WNoiseDy(p1, _noiseTile);
551   final[2] = WNoiseDz(p1, _noiseTile);
552   // UNUSED const float f1x = xUnwarped[0] * final[0] + xUnwarped[1] * final[1] + xUnwarped[2] * final[2];
553   const float f1y = yUnwarped[0] * final[0] + yUnwarped[1] * final[1] + yUnwarped[2] * final[2];
554   const float f1z = zUnwarped[0] * final[0] + zUnwarped[1] * final[1] + zUnwarped[2] * final[2];
555
556   final[0] = WNoiseDx(p2, _noiseTile);
557   final[1] = WNoiseDy(p2, _noiseTile);
558   final[2] = WNoiseDz(p2, _noiseTile);
559   const float f2x = xUnwarped[0] * final[0] + xUnwarped[1] * final[1] + xUnwarped[2] * final[2];
560   // UNUSED const float f2y = yUnwarped[0] * final[0] + yUnwarped[1] * final[1] + yUnwarped[2] * final[2];
561   const float f2z = zUnwarped[0] * final[0] + zUnwarped[1] * final[1] + zUnwarped[2] * final[2];
562
563   final[0] = WNoiseDx(p3, _noiseTile);
564   final[1] = WNoiseDy(p3, _noiseTile);
565   final[2] = WNoiseDz(p3, _noiseTile);
566   const float f3x = xUnwarped[0] * final[0] + xUnwarped[1] * final[1] + xUnwarped[2] * final[2];
567   const float f3y = yUnwarped[0] * final[0] + yUnwarped[1] * final[1] + yUnwarped[2] * final[2];
568   // UNUSED const float f3z = zUnwarped[0] * final[0] + zUnwarped[1] * final[1] + zUnwarped[2] * final[2];
569
570   Vec3 ret = Vec3( 
571       f3y - f2z,
572       f1z - f3x,
573       f2x - f1y ); 
574   return ret;
575 }
576
577
578 //////////////////////////////////////////////////////////////////////
579 // perform an actual noise advection step
580 //////////////////////////////////////////////////////////////////////
581 /*void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) 
582 {
583         // enlarge timestep to match grid
584         const float dt = dtOrg * _amplify;
585         const float invAmp = 1.0f / _amplify;
586         float *tempBig1 = new float[_totalCellsBig];
587         float *tempBig2 = new float[_totalCellsBig];
588         float *bigUx = new float[_totalCellsBig];
589         float *bigUy = new float[_totalCellsBig];
590         float *bigUz = new float[_totalCellsBig]; 
591         float *_energy = new float[_totalCellsSm];
592         float *highFreqEnergy = new float[_totalCellsSm];
593         float *eigMin  = new float[_totalCellsSm];
594         float *eigMax  = new float[_totalCellsSm];
595
596         memset(tempBig1, 0, sizeof(float)*_totalCellsBig);
597         memset(tempBig2, 0, sizeof(float)*_totalCellsBig);
598         memset(highFreqEnergy, 0, sizeof(float)*_totalCellsSm);
599         memset(eigMin, 0, sizeof(float)*_totalCellsSm);
600         memset(eigMax, 0, sizeof(float)*_totalCellsSm);
601
602         // prepare textures
603         advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2);
604
605   // compute eigenvalues of the texture coordinates
606   computeEigenvalues(eigMin, eigMax);
607
608   // do wavelet decomposition of energy
609   computeEnergy(_energy, xvel, yvel, zvel, obstacles);
610   decomposeEnergy(_energy, highFreqEnergy);
611
612   // zero out coefficients inside of the obstacle
613   for (int x = 0; x < _totalCellsSm; x++)
614     if (obstacles[x]) _energy[x] = 0.f;
615
616   float maxVelocity = 0.;
617   for (int z = 1; z < _zResBig - 1; z++) 
618     for (int y = 1; y < _yResBig - 1; y++) 
619       for (int x = 1; x < _xResBig - 1; x++)
620       {
621         // get unit position for both fine and coarse grid
622         const Vec3 pos = Vec3(x,y,z);
623         const Vec3 posSm = pos * invAmp;
624         
625         // get grid index for both fine and coarse grid
626         const int index = x + y *_xResBig + z *_slabSizeBig;
627         const int indexSmall = (int)posSm[0] + (int)posSm[1] * _xResSm + (int)posSm[2] * _slabSizeSm;
628         
629         // get a linearly interpolated velocity and texcoords
630         // from the coarse grid
631         Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel, 
632             posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
633         Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW, 
634             posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
635
636         // multiply the texture coordinate by _resSm so that turbulence
637         // synthesis begins at the first octave that the coarse grid 
638         // cannot capture
639         Vec3 texCoord = Vec3(uvw[0] * _resSm[0], 
640                              uvw[1] * _resSm[1],
641                              uvw[2] * _resSm[2]); 
642
643         // retrieve wavelet energy at highest frequency
644         float energy = INTERPOLATE::lerp3d(
645             highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
646
647         // base amplitude for octave 0
648         float coefficient = sqrtf(2.0f * fabs(energy));
649         const float amplitude = *_strength * fabs(0.5 * coefficient) * persistence;
650
651         // add noise to velocity, but only if the turbulence is
652         // sufficiently undeformed, and the energy is large enough
653         // to make a difference
654         const bool addNoise = eigMax[indexSmall] < 2. && 
655                               eigMin[indexSmall] > 0.5;
656         if (addNoise && amplitude > _cullingThreshold) {
657           // base amplitude for octave 0
658           float amplitudeScaled = amplitude;
659
660           for (int octave = 0; octave < _octaves; octave++)
661           {
662             // multiply the vector noise times the maximum allowed
663             // noise amplitude at this octave, and add it to the total
664             vel += WVelocity(texCoord) * amplitudeScaled;
665
666             // scale coefficient for next octave
667             amplitudeScaled *= persistence;
668             texCoord *= 2.0f;
669           }
670         }
671
672         // Store velocity + turbulence in big grid for maccormack step
673         //
674         // If you wanted to save memory, you would instead perform a 
675         // semi-Lagrangian backtrace for the current grid cell here. Then
676         // you could just throw the velocity away.
677         bigUx[index] = vel[0];
678         bigUy[index] = vel[1];
679         bigUz[index] = vel[2];
680
681         // compute the velocity magnitude for substepping later
682         const float velMag = bigUx[index] * bigUx[index] + 
683                              bigUy[index] * bigUy[index] + 
684                              bigUz[index] * bigUz[index];
685         if (velMag > maxVelocity) maxVelocity = velMag;
686
687         // zero out velocity inside obstacles
688         float obsCheck = INTERPOLATE::lerp3dToFloat(
689             obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); 
690         if (obsCheck > 0.95) 
691           bigUx[index] = bigUy[index] = bigUz[index] = 0.;
692       }
693
694   // prepare density for an advection
695   SWAP_POINTERS(_densityBig, _densityBigOld);
696
697   // based on the maximum velocity present, see if we need to substep,
698   // but cap the maximum number of substeps to 5
699   const int maxSubSteps = 5; 
700   maxVelocity = sqrt(maxVelocity) * dt;
701   int totalSubsteps = (int)(maxVelocity / (float)maxSubSteps);
702   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
703   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
704   const float dtSubdiv = dt / (float)totalSubsteps;
705
706   // set boundaries of big velocity grid
707   FLUID_3D::setZeroX(bigUx, _resBig, 0, _resBig[2]); 
708   FLUID_3D::setZeroY(bigUy, _resBig, 0, _resBig[2]); 
709   FLUID_3D::setZeroZ(bigUz, _resBig, 0, _resBig[2]);
710
711   // do the MacCormack advection, with substepping if necessary
712   for(int substep = 0; substep < totalSubsteps; substep++)
713   {
714     FLUID_3D::advectFieldMacCormack(dtSubdiv, bigUx, bigUy, bigUz, 
715         _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL);
716
717     if (substep < totalSubsteps - 1) 
718       SWAP_POINTERS(_densityBig, _densityBigOld);
719   } // substep
720   
721   // wipe the density borders
722   FLUID_3D::setZeroBorder(_densityBig, _resBig, 0, _resBig[2]);
723     
724   // reset texture coordinates now in preparation for next timestep
725   // Shouldn't do this before generating the noise because then the 
726   // eigenvalues stored do not reflect the underlying texture coordinates
727   resetTextureCoordinates(eigMin, eigMax);
728   
729   delete[] tempBig1;
730   delete[] tempBig2;
731   delete[] bigUx;
732   delete[] bigUy;
733   delete[] bigUz;
734   delete[] _energy;
735   delete[] highFreqEnergy;
736
737   delete[] eigMin;
738   delete[] eigMax;
739   
740
741   _totalStepsBig++;
742 }*/
743
744 //struct
745
746 //////////////////////////////////////////////////////////////////////
747 // perform the full turbulence algorithm, including OpenMP 
748 // if available
749 //////////////////////////////////////////////////////////////////////
750 void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles)
751 {
752         // enlarge timestep to match grid
753         const float dt = dtOrg * _amplify;
754         const float invAmp = 1.0f / _amplify;
755         float *tempBig1 = (float *)calloc(_totalCellsBig, sizeof(float));
756         float *tempBig2 = (float *)calloc(_totalCellsBig, sizeof(float));
757         float *bigUx = (float *)calloc(_totalCellsBig, sizeof(float));
758         float *bigUy = (float *)calloc(_totalCellsBig, sizeof(float));
759         float *bigUz = (float *)calloc(_totalCellsBig, sizeof(float)); 
760         float *_energy = (float *)calloc(_totalCellsSm, sizeof(float));
761         float *highFreqEnergy = (float *)calloc(_totalCellsSm, sizeof(float));
762         float *eigMin  = (float *)calloc(_totalCellsSm, sizeof(float));
763         float *eigMax  = (float *)calloc(_totalCellsSm, sizeof(float));
764
765         memset(_tcTemp, 0, sizeof(float)*_totalCellsSm);
766
767
768         // prepare textures
769         advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2);
770
771         // do wavelet decomposition of energy
772         computeEnergy(_energy, xvel, yvel, zvel, obstacles);
773
774         for (int x = 0; x < _totalCellsSm; x++)
775                 if (obstacles[x]) _energy[x] = 0.f;
776
777         decomposeEnergy(_energy, highFreqEnergy);
778
779         // zero out coefficients inside of the obstacle
780         for (int x = 0; x < _totalCellsSm; x++)
781                 if (obstacles[x]) highFreqEnergy[x] = 0.f;
782
783         Vec3Int ressm(_xResSm, _yResSm, _zResSm);
784         FLUID_3D::setNeumannX(highFreqEnergy, ressm, 0 , ressm[2]);
785         FLUID_3D::setNeumannY(highFreqEnergy, ressm, 0 , ressm[2]);
786         FLUID_3D::setNeumannZ(highFreqEnergy, ressm, 0 , ressm[2]);
787
788
789    int threadval = 1;
790 #if PARALLEL==1
791   threadval = omp_get_max_threads();
792 #endif
793
794
795   // parallel region setup
796   // Uses omp_get_max_trheads to get number of required cells.
797   float* maxVelMagThreads = new float[threadval];
798
799   for (int i=0; i<threadval; i++) maxVelMagThreads[i] = -1.0f;
800
801 #if PARALLEL==1
802
803 #pragma omp parallel
804 #endif
805   { float maxVelMag1 = 0.;
806 #if PARALLEL==1
807     const int id  = omp_get_thread_num(); /*, num = omp_get_num_threads(); */
808 #endif
809
810   // vector noise main loop
811 #if PARALLEL==1
812 #pragma omp for schedule(static,1)
813 #endif
814   for (int zSmall = 0; zSmall < _zResSm; zSmall++)
815   {
816   for (int ySmall = 0; ySmall < _yResSm; ySmall++) 
817   for (int xSmall = 0; xSmall < _xResSm; xSmall++)
818   {
819     const int indexSmall = xSmall + ySmall * _xResSm + zSmall * _slabSizeSm;
820
821     // compute jacobian
822     float jacobian[3][3] = {
823       { minDx(xSmall, ySmall, zSmall, _tcU, _resSm), minDx(xSmall, ySmall, zSmall, _tcV, _resSm), minDx(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
824       { minDy(xSmall, ySmall, zSmall, _tcU, _resSm), minDy(xSmall, ySmall, zSmall, _tcV, _resSm), minDy(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
825       { minDz(xSmall, ySmall, zSmall, _tcU, _resSm), minDz(xSmall, ySmall, zSmall, _tcV, _resSm), minDz(xSmall, ySmall, zSmall, _tcW, _resSm) }
826     };
827
828     // get LU factorization of texture jacobian and apply 
829     // it to unit vectors
830     sLU LU = computeLU(jacobian);
831     float xUnwarped[3], yUnwarped[3], zUnwarped[3];
832     float xWarped[3], yWarped[3], zWarped[3];
833     bool nonSingular = isNonsingular(LU);
834
835         xUnwarped[0] = 1.0f; xUnwarped[1] = 0.0f; xUnwarped[2] = 0.0f;
836         yUnwarped[0] = 0.0f; yUnwarped[1] = 1.0f; yUnwarped[2] = 0.0f;
837         zUnwarped[0] = 0.0f; zUnwarped[1] = 0.0f; zUnwarped[2] = 1.0f;
838
839         xWarped[0] = 1.0f; xWarped[1] = 0.0f; xWarped[2] = 0.0f;
840         yWarped[0] = 0.0f; yWarped[1] = 1.0f; yWarped[2] = 0.0f;
841         zWarped[0] = 0.0f; zWarped[1] = 0.0f; zWarped[2] = 1.0f;
842
843 #if 0
844         // UNUSED
845     float eigMax = 10.0f;
846     float eigMin = 0.1f;
847 #endif
848     if (nonSingular)
849     {
850       solveLU3x3(LU, xUnwarped, xWarped);
851       solveLU3x3(LU, yUnwarped, yWarped);
852       solveLU3x3(LU, zUnwarped, zWarped);
853
854       // compute the eigenvalues while we have the Jacobian available
855       Vec3 eigenvalues = Vec3(1.);
856       computeEigenvalues3x3( &eigenvalues[0], jacobian);
857       eigMax[indexSmall] = MAX3V(eigenvalues);
858       eigMin[indexSmall] = MIN3V(eigenvalues);
859     }
860     
861     // make sure to skip one on the beginning and end
862     int xStart = (xSmall == 0) ? 1 : 0;
863     int xEnd   = (xSmall == _xResSm - 1) ? _amplify - 1 : _amplify;
864     int yStart = (ySmall == 0) ? 1 : 0;
865     int yEnd   = (ySmall == _yResSm - 1) ? _amplify - 1 : _amplify;
866     int zStart = (zSmall == 0) ? 1 : 0;
867     int zEnd   = (zSmall == _zResSm - 1) ? _amplify - 1 : _amplify;
868       
869     for (int zBig = zStart; zBig < zEnd; zBig++) 
870     for (int yBig = yStart; yBig < yEnd; yBig++) 
871     for (int xBig = xStart; xBig < xEnd; xBig++)
872     {
873       const int x = xSmall * _amplify + xBig;
874       const int y = ySmall * _amplify + yBig;
875       const int z = zSmall * _amplify + zBig;
876       
877       // get unit position for both fine and coarse grid
878       const Vec3 pos = Vec3(x,y,z);
879       const Vec3 posSm = pos * invAmp;
880       
881       // get grid index for both fine and coarse grid
882       const int index = x + y *_xResBig + z *_slabSizeBig;
883       
884       // get a linearly interpolated velocity and texcoords
885       // from the coarse grid
886       Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel, 
887           posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
888       Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW, 
889           posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
890
891       // multiply the texture coordinate by _resSm so that turbulence
892       // synthesis begins at the first octave that the coarse grid 
893       // cannot capture
894       Vec3 texCoord = Vec3(uvw[0] * _resSm[0], 
895                            uvw[1] * _resSm[1],
896                            uvw[2] * _resSm[2]); 
897
898       // retrieve wavelet energy at highest frequency
899       float energy = INTERPOLATE::lerp3d(
900           highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
901
902       // base amplitude for octave 0
903       float coefficient = sqrtf(2.0f * fabs(energy));
904       const float amplitude = *_strength * fabs(0.5 * coefficient) * persistence;
905
906       // add noise to velocity, but only if the turbulence is
907       // sufficiently undeformed, and the energy is large enough
908       // to make a difference
909       const bool addNoise = eigMax[indexSmall] < 2. && 
910                             eigMin[indexSmall] > 0.5;
911       if (addNoise && amplitude > _cullingThreshold) {
912         // base amplitude for octave 0
913         float amplitudeScaled = amplitude;
914
915         for (int octave = 0; octave < _octaves; octave++)
916         {
917           // multiply the vector noise times the maximum allowed
918           // noise amplitude at this octave, and add it to the total
919           vel += WVelocityWithJacobian(texCoord, &xUnwarped[0], &yUnwarped[0], &zUnwarped[0]) * amplitudeScaled;
920
921           // scale coefficient for next octave
922           amplitudeScaled *= persistence;
923           texCoord *= 2.0f;
924         }
925       }
926
927       // Store velocity + turbulence in big grid for maccormack step
928       //
929       // If you wanted to save memory, you would instead perform a 
930       // semi-Lagrangian backtrace for the current grid cell here. Then
931       // you could just throw the velocity away.
932       bigUx[index] = vel[0];
933       bigUy[index] = vel[1];
934       bigUz[index] = vel[2];
935
936       // compute the velocity magnitude for substepping later
937       const float velMag = bigUx[index] * bigUx[index] + 
938                            bigUy[index] * bigUy[index] + 
939                            bigUz[index] * bigUz[index];
940       if (velMag > maxVelMag1) maxVelMag1 = velMag;
941
942       // zero out velocity inside obstacles
943       float obsCheck = INTERPOLATE::lerp3dToFloat(
944           obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); 
945       if (obsCheck > 0.95) 
946         bigUx[index] = bigUy[index] = bigUz[index] = 0.;
947     } // xyz*/
948
949 #if PARALLEL==1
950     maxVelMagThreads[id] = maxVelMag1;
951 #else
952     maxVelMagThreads[0] = maxVelMag1;
953 #endif
954   }
955   }
956   } // omp
957   
958   // compute maximum over threads
959   float maxVelMag = maxVelMagThreads[0];
960 #if PARALLEL==1
961   for (int i = 1; i < threadval; i++) 
962     if (maxVelMag < maxVelMagThreads[i]) 
963       maxVelMag = maxVelMagThreads[i];
964 #endif
965   delete [] maxVelMagThreads;
966
967
968   // prepare density for an advection
969   SWAP_POINTERS(_densityBig, _densityBigOld);
970
971   // based on the maximum velocity present, see if we need to substep,
972   // but cap the maximum number of substeps to 5
973   const int maxSubSteps = 25;
974   const int maxVel = 5;
975   maxVelMag = sqrt(maxVelMag) * dt;
976   int totalSubsteps = (int)(maxVelMag / (float)maxVel);
977   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
978   // printf("totalSubsteps: %d\n", totalSubsteps);
979   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
980   const float dtSubdiv = dt / (float)totalSubsteps;
981
982   // set boundaries of big velocity grid
983   FLUID_3D::setZeroX(bigUx, _resBig, 0 , _resBig[2]); 
984   FLUID_3D::setZeroY(bigUy, _resBig, 0 , _resBig[2]); 
985   FLUID_3D::setZeroZ(bigUz, _resBig, 0 , _resBig[2]);
986
987 #if PARALLEL==1
988   int stepParts = threadval*2;  // Dividing parallelized sections into numOfThreads * 2 sections
989   float partSize = (float)_zResBig/stepParts;   // Size of one part;
990
991   if (partSize < 4) {stepParts = threadval;                                     // If the slice gets too low (might actually slow things down, change it to larger
992                                         partSize = (float)_zResBig/stepParts;}
993   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
994                                         partSize = (float)_zResBig/stepParts;}
995 #else
996   int zBegin=0;
997   int zEnd=_resBig[2];
998 #endif
999
1000   // do the MacCormack advection, with substepping if necessary
1001   for(int substep = 0; substep < totalSubsteps; substep++)
1002   {
1003
1004 #if PARALLEL==1
1005         #pragma omp parallel
1006         {
1007
1008         #pragma omp for schedule(static,1)
1009         for (int i=0; i<stepParts; i++)
1010         {
1011                 int zBegin = (int)((float)i*partSize + 0.5f);
1012                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
1013 #endif
1014                 FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz, 
1015                     _densityBigOld, tempBig1, _resBig, zBegin, zEnd);
1016 #if PARALLEL==1
1017         }
1018
1019         #pragma omp barrier
1020
1021         #pragma omp for schedule(static,1)
1022         for (int i=0; i<stepParts; i++)
1023         {
1024                 int zBegin = (int)((float)i*partSize + 0.5f);
1025                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
1026 #endif
1027                 FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz, 
1028                     _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL, zBegin, zEnd);
1029 #if PARALLEL==1
1030         }
1031         }
1032 #endif
1033
1034     if (substep < totalSubsteps - 1) 
1035       SWAP_POINTERS(_densityBig, _densityBigOld);
1036   } // substep
1037
1038   free(tempBig1);
1039   free(tempBig2);
1040   free(bigUx);
1041   free(bigUy);
1042   free(bigUz);
1043   free(_energy);
1044   free(highFreqEnergy);
1045   
1046   // wipe the density borders
1047   FLUID_3D::setZeroBorder(_densityBig, _resBig, 0 , _resBig[2]);
1048     
1049   // reset texture coordinates now in preparation for next timestep
1050   // Shouldn't do this before generating the noise because then the 
1051   // eigenvalues stored do not reflect the underlying texture coordinates
1052   resetTextureCoordinates(eigMin, eigMax);
1053
1054   free(eigMin);
1055   free(eigMax);
1056   
1057   // output files
1058   // string prefix = string("./amplified.preview/density_bigxy_");
1059   // FLUID_3D::writeImageSliceXY(_densityBig, _resBig, _resBig[2]/2, prefix, _totalStepsBig, 1.0f);
1060   //string df3prefix = string("./df3/density_big_");
1061   //IMAGE::dumpDF3(_totalStepsBig, df3prefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
1062   // string pbrtPrefix = string("./pbrt/density_big_");
1063   // IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
1064   
1065   _totalStepsBig++;
1066 }