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