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