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