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