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