Smoke:
authorDaniel Genrich <daniel.genrich@gmx.net>
Wed, 9 Sep 2009 18:39:40 +0000 (18:39 +0000)
committerDaniel Genrich <daniel.genrich@gmx.net>
Wed, 9 Sep 2009 18:39:40 +0000 (18:39 +0000)
* Enable cache for high res + new preview
* Bugfix for smoke banding (in cooperation with N_T)

Hint: Work-in-progress regarding collision objects so can be broken, didn't test

Hint2: jahka enabled a general particle panel but
* bake button doesn't work
* step is not supported for cloth
* several other things there ;)

21 files changed:
intern/smoke/extern/smoke_API.h
intern/smoke/intern/FLUID_3D.cpp
intern/smoke/intern/FLUID_3D_STATIC.cpp
intern/smoke/intern/WTURBULENCE.cpp
intern/smoke/intern/WTURBULENCE.h
intern/smoke/intern/smoke_API.cpp
release/ui/buttons_physics_smoke.py
source/blender/blenkernel/BKE_pointcache.h
source/blender/blenkernel/BKE_smoke.h
source/blender/blenkernel/intern/pointcache.c
source/blender/blenkernel/intern/smoke.c
source/blender/blenloader/intern/readfile.c
source/blender/editors/space_view3d/drawobject.c
source/blender/editors/space_view3d/drawvolume.c
source/blender/editors/space_view3d/view3d_intern.h
source/blender/gpu/CMakeLists.txt
source/blender/gpu/GPU_draw.h
source/blender/gpu/intern/gpu_draw.c
source/blender/gpu/intern/gpu_extensions.c
source/blender/makesdna/DNA_smoke_types.h
source/blender/makesrna/intern/rna_smoke.c

index 1a3edce2344ef554aeaee1e848e37da625cf6aea..5607df70cf3b8273254fa006e627b3a8afd9110f 100644 (file)
@@ -63,11 +63,11 @@ void smoke_turbulence_free(struct WTURBULENCE *wt);
 void smoke_turbulence_step(struct WTURBULENCE *wt, struct FLUID_3D *fluid);
 
 float *smoke_turbulence_get_density(struct WTURBULENCE *wt);
-void smoke_turbulence_get_res(struct WTURBULENCE *wt, unsigned int *res);
+void smoke_turbulence_get_res(struct WTURBULENCE *wt, int *res);
 void smoke_turbulence_set_noise(struct WTURBULENCE *wt, int type);
-void smoke_turbulence_initBlenderRNA(struct WTURBULENCE *wt, float *strength);
+void smoke_initWaveletBlenderRNA(struct WTURBULENCE *wt, float *strength);
 
-void smoke_turbulence_dissolve(struct WTURBULENCE *wt, int speed, int log);
+void smoke_dissolve_wavelet(struct WTURBULENCE *wt, int speed, int log);
 
 // export
 void smoke_turbulence_export(struct WTURBULENCE *wt, float **dens, float **densold, float **tcu, float **tcv, float **tcw);
index 89dd893198b3546c31dbe7491afbff4d04148b93..7574a0ec16e45083b475d4a3a7b63399451bfae1 100644 (file)
@@ -27,9 +27,9 @@
 #include <zlib.h>
 
 // boundary conditions of the fluid domain
-#define DOMAIN_BC_FRONT  1
-#define DOMAIN_BC_TOP    0
-#define DOMAIN_BC_LEFT   1
+#define DOMAIN_BC_FRONT  0 // z
+#define DOMAIN_BC_TOP    1 // y
+#define DOMAIN_BC_LEFT   1 // x
 #define DOMAIN_BC_BACK   DOMAIN_BC_FRONT
 #define DOMAIN_BC_BOTTOM DOMAIN_BC_TOP
 #define DOMAIN_BC_RIGHT  DOMAIN_BC_LEFT
@@ -111,47 +111,42 @@ FLUID_3D::FLUID_3D(int *res, float *p0, float dt) :
        }
 
        // set side obstacles
-  size_t index;
-  for (int y = 0; y < _yRes; y++) // z
-    for (int x = 0; x < _xRes; x++)
-    {
-      // front slab
-      index = x + y * _xRes;
-      if(DOMAIN_BC_BOTTOM==1) _obstacles[index] = 1;
+       int index;
+       for (int y = 0; y < _yRes; y++)
+       for (int x = 0; x < _xRes; x++)
+       {
+               // front slab
+               index = x + y * _xRes;
+               if(DOMAIN_BC_FRONT==1) _obstacles[index] = 1;
 
-      // back slab
-      index += _totalCells - _slabSize;
-      if(DOMAIN_BC_TOP==1) _obstacles[index] = 1;
-    }
-  for (int z = 0; z < _zRes; z++) // y
-    for (int x = 0; x < _xRes; x++)
-    {
-      // bottom slab
-      index = x + z * _slabSize;
-      if(DOMAIN_BC_FRONT==1) _obstacles[index] = 1;
+               // back slab
+               index += _totalCells - _slabSize;
+               if(DOMAIN_BC_BACK==1) _obstacles[index] = 1;
+       }
 
-      // top slab
-      index += _slabSize - _xRes;
-      if(DOMAIN_BC_BACK==1) _obstacles[index] = 1;
-    }
-  for (int z = 0; z < _zRes; z++) // x
-    for (int y = 0; y < _yRes; y++)
-    {
-      // left slab
-      index = y * _xRes + z * _slabSize;
-      if(DOMAIN_BC_LEFT==1) _obstacles[index] = 1;
+       for (int z = 0; z < _zRes; z++)
+       for (int x = 0; x < _xRes; x++)
+       {
+               // bottom slab
+               index = x + z * _slabSize;
+               if(DOMAIN_BC_BOTTOM==1) _obstacles[index] = 1;
 
-      // right slab
-      index += _xRes - 1;
-      if(DOMAIN_BC_RIGHT==1) _obstacles[index] = 1;
-    }
+               // top slab
+               index += _slabSize - _xRes;
+               if(DOMAIN_BC_TOP==1) _obstacles[index] = 1;
+       }
 
-       /*
-       SPHERE *obsSphere = NULL;
-       obsSphere = new SPHERE(0.375,0.5,0.375, 0.1); // for 4 to 3 domain
-       addObstacle(obsSphere);
-       delete obsSphere;
-       */
+       for (int z = 0; z < _zRes; z++)
+       for (int y = 0; y < _yRes; y++)
+       {
+               // left slab
+               index = y * _xRes + z * _slabSize;
+               if(DOMAIN_BC_LEFT==1) _obstacles[index] = 1;
+
+               // right slab
+               index += _xRes - 1;
+               if(DOMAIN_BC_RIGHT==1) _obstacles[index] = 1;
+       }
 }
 
 FLUID_3D::~FLUID_3D()
@@ -191,7 +186,7 @@ void FLUID_3D::step()
        for (int i = 0; i < _totalCells; i++)
        {
                _xForce[i] = _yForce[i] = _zForce[i] = 0.0f;
-               _obstacles[i] &= ~2;
+               // _obstacles[i] &= ~2;
        }
 
        wipeBoundaries();
@@ -232,7 +227,8 @@ void FLUID_3D::step()
        _totalTime += _dt;
        _totalSteps++;  
 
-       memset(_obstacles, 0, sizeof(unsigned char)*_xRes*_yRes*_zRes);
+       // todo xxx dg: only clear obstacles, not boundaries
+       // memset(_obstacles, 0, sizeof(unsigned char)*_xRes*_yRes*_zRes);
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -270,7 +266,7 @@ void FLUID_3D::artificialDamping(float* field) {
 //////////////////////////////////////////////////////////////////////
 void FLUID_3D::copyBorderAll(float* field)
 {
-       size_t index;
+       int index;
        for (int y = 0; y < _yRes; y++)
                for (int x = 0; x < _xRes; x++)
                {
@@ -350,13 +346,13 @@ void FLUID_3D::project()
 
        // copy out the boundaries
        if(DOMAIN_BC_LEFT == 0)  setNeumannX(_xVelocity, _res);
-       else setZeroX(_xVelocity, _res);
+       else setZeroX(_xVelocity, _res); 
 
-       if(DOMAIN_BC_TOP == 0)   setNeumannZ(_zVelocity, _res);
-       else setZeroZ(_zVelocity, _res);
+       if(DOMAIN_BC_TOP == 0)   setNeumannY(_yVelocity, _res);
+       else setZeroY(_yVelocity, _res); 
 
-       if(DOMAIN_BC_FRONT == 0) setNeumannY(_yVelocity, _res);
-       else setZeroY(_yVelocity, _res);
+       if(DOMAIN_BC_FRONT == 0) setNeumannZ(_zVelocity, _res);
+       else setZeroZ(_zVelocity, _res);
 
        // calculate divergence
        index = _slabSize + _xRes + 1;
@@ -400,12 +396,16 @@ void FLUID_3D::project()
                for (y = 1; y < _yRes - 1; y++, index += 2)
                        for (x = 1; x < _xRes - 1; x++, index++)
                        {
-                               if(!_obstacles[index])
+                               // if(!_obstacles[index])
                                {
                                        _xVelocity[index] -= 0.5f * (_pressure[index + 1]     - _pressure[index - 1]) * invDx;
                                        _yVelocity[index] -= 0.5f * (_pressure[index + _xRes]  - _pressure[index - _xRes]) * invDx;
                                        _zVelocity[index] -= 0.5f * (_pressure[index + _slabSize] - _pressure[index - _slabSize]) * invDx;
-                               }
+                               }/*
+                               else
+                               {
+                                       _xVelocity[index] = _yVelocity[index] = _zVelocity[index] = 0.0f;
+                               }*/
                        }
 
        if (_pressure) delete[] _pressure;
@@ -669,13 +669,13 @@ void FLUID_3D::advectMacCormack()
        Vec3Int res = Vec3Int(_xRes,_yRes,_zRes);
 
        if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocity, res);
-       else setZeroX(_xVelocity, res);
+       else setZeroX(_xVelocity, res); 
 
-       if(DOMAIN_BC_TOP == 0) copyBorderZ(_zVelocity, res);
-       else setZeroZ(_zVelocity, res);
+       if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocity, res);
+       else setZeroY(_yVelocity, res); 
 
-       if(DOMAIN_BC_FRONT == 0) copyBorderY(_yVelocity, res);
-       else setZeroY(_yVelocity, res);
+       if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocity, res);
+       else setZeroZ(_zVelocity, res);
 
        SWAP_POINTERS(_xVelocity, _xVelocityOld);
        SWAP_POINTERS(_yVelocity, _yVelocityOld);
@@ -698,13 +698,13 @@ void FLUID_3D::advectMacCormack()
        advectFieldMacCormack(dt0, _xVelocityOld, _yVelocityOld, _zVelocityOld, _zVelocityOld, _zVelocity, t1,t2, res, _obstacles);
 
        if(DOMAIN_BC_LEFT == 0) copyBorderX(_xVelocity, res);
-       else setZeroX(_xVelocity, res);
+       else setZeroX(_xVelocity, res); 
 
-       if(DOMAIN_BC_TOP == 0) copyBorderZ(_zVelocity, res);
-       else setZeroZ(_zVelocity, res);
+       if(DOMAIN_BC_TOP == 0) copyBorderY(_yVelocity, res);
+       else setZeroY(_yVelocity, res); 
 
-       if(DOMAIN_BC_FRONT == 0) copyBorderY(_yVelocity, res);
-       else setZeroY(_yVelocity, res);
+       if(DOMAIN_BC_FRONT == 0) copyBorderZ(_zVelocity, res);
+       else setZeroZ(_zVelocity, res);
 
        setZeroBorder(_density, res);
        setZeroBorder(_heat, res);
index 4474129beea7a89f6e7f51572855ef9960fc42b6..4909c071c3dd3c66dcebfd09b7956501c2b77aaa 100644 (file)
@@ -80,7 +80,7 @@ void FLUID_3D::addSmokeTestCase(float* field, Vec3Int res, float value)
 void FLUID_3D::setNeumannX(float* field, Vec3Int res)
 {
        const int slabSize = res[0] * res[1];
-       size_t index;
+       int index;
        for (int z = 0; z < res[2]; z++)
                for (int y = 0; y < res[1]; y++)
                {
@@ -92,6 +92,18 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res)
                        index += res[0] - 1;
                        field[index] = field[index - 2];
                }
+
+       // fix, force top slab to only allow outwards flux
+       for (int y = 0; y < res[1]; y++)
+               for (int z = 0; z < res[2]; z++)
+               {
+                       // top slab
+                       int index = y * res[0] + z * slabSize;
+                       index += res[0] - 1;
+                       if(field[index]<0.) field[index] = 0.;
+                       index -= 1;
+                       if(field[index]<0.) field[index] = 0.;
+               }
  }
 
 //////////////////////////////////////////////////////////////////////
@@ -100,18 +112,31 @@ void FLUID_3D::setNeumannX(float* field, Vec3Int res)
 void FLUID_3D::setNeumannY(float* field, Vec3Int res)
 {
        const int slabSize = res[0] * res[1];
-       size_t index;
+       int index;
        for (int z = 0; z < res[2]; z++)
                for (int x = 0; x < res[0]; x++)
                {
-                       // front slab
+                       // bottom slab
                        index = x + z * slabSize;
                        field[index] = field[index + 2 * res[0]];
 
-                       // back slab
+                       // top slab
                        index += slabSize - res[0];
                        field[index] = field[index - 2 * res[0]];
                }
+
+       // fix, force top slab to only allow outwards flux
+       for (int z = 0; z < res[2]; z++)
+               for (int x = 0; x < res[0]; x++)
+               {
+                       // top slab
+                       int index = x + z * slabSize;
+                       index += slabSize - res[0];
+                       if(field[index]<0.) field[index] = 0.;
+                       index -= res[0];
+                       if(field[index]<0.) field[index] = 0.;
+               }
+               
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -121,15 +146,15 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res)
 {
        const int slabSize = res[0] * res[1];
        const int totalCells = res[0] * res[1] * res[2];
-       size_t index;
+       int index;
        for (int y = 0; y < res[1]; y++)
                for (int x = 0; x < res[0]; x++)
                {
-                       // bottom slab
+                       // front slab
                        index = x + y * res[0];
                        field[index] = field[index + 2 * slabSize];
 
-                       // top slab
+                       // back slab
                        index += totalCells - slabSize;
                        field[index] = field[index - 2 * slabSize];
                }
@@ -139,11 +164,11 @@ void FLUID_3D::setNeumannZ(float* field, Vec3Int res)
                for (int x = 0; x < res[0]; x++)
                {
                        // top slab
-                       index = x + y * res[0];
+                       int index = x + y * res[0];
                        index += totalCells - slabSize;
-                       if(field[index]<0.) field[index] = 0.0f;
+                       if(field[index]<0.) field[index] = 0.;
                        index -= slabSize;
-                       if(field[index]<0.) field[index] = 0.0f;
+                       if(field[index]<0.) field[index] = 0.;
                }
                
 }
@@ -231,13 +256,14 @@ void FLUID_3D::copyBorderX(float* field, Vec3Int res)
 void FLUID_3D::copyBorderY(float* field, Vec3Int res)
 {
        const int slabSize = res[0] * res[1];
+       const int totalCells = res[0] * res[1] * res[2];
        int index;
        for (int z = 0; z < res[2]; z++)
                for (int x = 0; x < res[0]; x++)
                {
                        // bottom slab
                        index = x + z * slabSize;
-                       field[index] = field[index + res[0]];
+                       field[index] = field[index + res[0]]; 
                        // top slab
                        index += slabSize - res[0];
                        field[index] = field[index - res[0]];
@@ -253,7 +279,7 @@ void FLUID_3D::copyBorderZ(float* field, Vec3Int res)
                {
                        // front slab
                        index = x + y * res[0];
-                       field[index] = field[index + slabSize];
+                       field[index] = field[index + slabSize]; 
                        // back slab
                        index += totalCells - slabSize;
                        field[index] = field[index - slabSize];
@@ -269,18 +295,21 @@ void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const
        const int xres = res[0];
        const int yres = res[1];
        const int zres = res[2];
+       static int hits = 0;
+       static int total = 0;
        const int slabSize = res[0] * res[1];
 
        // scale dt up to grid resolution
-#if PARALLEL==1
-#pragma omp parallel for schedule(static)
+#if PARALLEL==1 && !_WIN32
+#pragma omp parallel
+#pragma omp for  schedule(static)
 #endif
        for (int z = 0; z < zres; z++)
                for (int y = 0; y < yres; y++)
                        for (int x = 0; x < xres; x++)
                        {
                                const int index = x + y * xres + z * xres*yres;
-
+                               
         // backtrace
                                float xTrace = x - dt * velx[index];
                                float yTrace = y - dt * vely[index];
@@ -337,8 +366,8 @@ void FLUID_3D::advectFieldSemiLagrange(const float dt, const float* velx, const
 //
 // comments are the pseudocode from selle's paper
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity,
-               float* oldField, float* newField, float* temp1, float* temp2, Vec3Int res, const unsigned char* obstacles)
+void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, const float* yVelocity, const float* zVelocity, 
+                               float* oldField, float* newField, float* temp1, float* temp2, Vec3Int res, const unsigned char* obstacles)
 {
        float* phiHatN  = temp1;
        float* phiHatN1 = temp2;
@@ -359,7 +388,7 @@ void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, con
        advectFieldSemiLagrange( -1.0*dt, xVelocity, yVelocity, zVelocity, phiHatN1, phiHatN, res);
 
        // phiN1 = phiHatN1 + (phiN - phiHatN) / 2
-       const int border = 0;
+       const int border = 0; 
        for (int z = border; z < sz-border; z++)
                for (int y = border; y < sy-border; y++)
                        for (int x = border; x < sx-border; x++) {
@@ -376,7 +405,7 @@ void FLUID_3D::advectFieldMacCormack(const float dt, const float* xVelocity, con
 
        // if the error estimate was bad, revert to first order
        clampOutsideRays(dt, xVelocity, yVelocity, zVelocity, oldField, newField, res, obstacles, phiHatN1);
-}
+} 
 
 
 //////////////////////////////////////////////////////////////////////
@@ -454,17 +483,18 @@ void FLUID_3D::clampExtrema(const float dt, const float* velx, const float* vely
 }
 
 //////////////////////////////////////////////////////////////////////
-// Reverts any backtraces that go into boundaries back to first
+// Reverts any backtraces that go into boundaries back to first 
 // order -- in this case the error correction term was totally
 // incorrect
 //////////////////////////////////////////////////////////////////////
-void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely, const float* velz,
-               float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection)
+void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float* vely,  const float* velz,
+                               float* oldField, float* newField, Vec3Int res, const unsigned char* obstacles, const float *oldAdvection)
 {
        const int sx= res[0];
        const int sy= res[1];
        const int sz= res[2];
        const int slabSize = res[0] * res[1];
+
        for (int z = 1; z < sz-1; z++)
                for (int y = 1; y < sy-1; y++)
                        for (int x = 1; x < sx-1; x++)
@@ -479,7 +509,7 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float*
                                float zTrace    = z - dt * velz[index];
 
                                // see if it goes outside the boundaries
-                               bool hasObstacle =
+                               bool hasObstacle = 
                                        (zTrace < 1.0f)    || (zTrace > sz - 2.0f) ||
                                        (yTrace < 1.0f)    || (yTrace > sy - 2.0f) ||
                                        (xTrace < 1.0f)    || (xTrace > sx - 2.0f) ||
@@ -515,7 +545,7 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float*
                                int z0 = (int)zBackward;
                                int z1 = z0 + 1;
                                if(obstacles && !hasObstacle) {
-                                       hasObstacle = hasObstacle ||
+                                       hasObstacle = hasObstacle || 
                                                obstacles[x0 + y0 * sx + z0*slabSize] ||
                                                obstacles[x0 + y1 * sx + z0*slabSize] ||
                                                obstacles[x1 + y0 * sx + z0*slabSize] ||
@@ -535,7 +565,7 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float*
                                z0 = (int)zTrace;
                                z1 = z0 + 1;
                                if(obstacles && !hasObstacle) {
-                                       hasObstacle = hasObstacle ||
+                                       hasObstacle = hasObstacle || 
                                                obstacles[x0 + y0 * sx + z0*slabSize] ||
                                                obstacles[x0 + y1 * sx + z0*slabSize] ||
                                                obstacles[x1 + y0 * sx + z0*slabSize] ||
@@ -577,84 +607,8 @@ void FLUID_3D::clampOutsideRays(const float dt, const float* velx, const float*
                                                u1 * (s0 * (t0 * oldField[i001] +
                                                                        t1 * oldField[i011]) +
                                                                s1 * (t0 * oldField[i101] +
-                                                                       t1 * oldField[i111]));
+                                                                       t1 * oldField[i111])); 
                                }
                        } // xyz
 }
 
-//////////////////////////////////////////////////////////////////////
-// image output
-//////////////////////////////////////////////////////////////////////
-/*
-void FLUID_3D::writeImageSliceXY(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale) {
-  writeProjectedIntern(field, res, 0,1, prefix, picCnt, scale);
-}
-void FLUID_3D::writeImageSliceYZ(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale) {
-  writeProjectedIntern(field, res, 1,2, prefix, picCnt, scale);
-}
-void FLUID_3D::writeImageSliceXZ(const float *field, Vec3Int res, int slice, string prefix, int picCnt, float scale) {
-  writeProjectedIntern(field, res, 0,2, prefix, picCnt, scale);
-}
-*/
-
-//////////////////////////////////////////////////////////////////////
-// Helper function for projecting densities along a dimension
-//////////////////////////////////////////////////////////////////////
-/*
-static int getOtherDir(int dir1, int dir2) {
-       switch(dir1) {
-               case 0:
-                       switch(dir2) {
-                               case 1: return 2;
-                               case 2: return 1; }
-                       break;
-               case 1:
-                       switch(dir2) {
-                               case 0: return 2;
-                               case 2: return 0; }
-                       break;
-               case 2:
-                       switch(dir2) {
-                               case 0: return 1;
-                               case 1: return 0; }
-                       break;
-               default:
-                       return 0;
-       }
-       return 0;
-}
-*/
-
-//////////////////////////////////////////////////////////////////////
-// average densities along third spatial direction
-//////////////////////////////////////////////////////////////////////
-/*
-void FLUID_3D::writeProjectedIntern(const float *field, Vec3Int res,
-               int dir1, int dir2, string prefix, int picCnt, float scale) {
-       const int nitems = res[dir1]*res[dir2];
-       const int otherDir = getOtherDir(dir1,dir2);
-       float *buf = new float[nitems];
-       Vec3Int min = Vec3Int(0);
-       Vec3Int max = res;
-
-       min[otherDir] = 0;
-       max[otherDir] = res[otherDir];
-       float div = 1./(float)MIN3V(res); // normalize for shorter sides, old: res[otherDir];
-       div *= 4.; //slightly increase contrast
-       for(int i=0; i<nitems; i++) buf[i] = 0.;
-
-       Vec3Int cnt = 0;
-       for (cnt[2] = min[2]; cnt[2] < max[2]; cnt[2]++) {
-               for (cnt[1] = min[1]; cnt[1] < max[1]; cnt[1]++)
-                       for (cnt[0] = min[0]; cnt[0] < max[0]; cnt[0]++)
-                       {
-                               const int index = cnt[0] + cnt[1] * res[0] + cnt[2] * res[0]*res[1];
-                               const int bufindex = cnt[dir1] + cnt[dir2] * res[dir1];
-                               buf[bufindex] += field[index] * scale *div;
-                       }
-       }
-       // IMAGE::dumpNumberedPNG(picCnt, prefix, buf, res[dir1], res[dir2]);
-       delete[] buf;
-}
-*/
-
index dd092d4f0cc3e72386192f962538313728bb3695..a1b2aaf30f2292f17df164c7e9c4ad6555d7eb58 100644 (file)
@@ -92,10 +92,6 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int no
        _tcW = new float[_totalCellsSm];
        _tcTemp = new float[_totalCellsSm];
        
-       // allocate & init energy terms
-       _energy = new float[_totalCellsSm];
-       _highFreqEnergy = new float[_totalCellsSm];
-       
        // map all 
        const float dx = 1./(float)(_resSm[0]);
        const float dy = 1./(float)(_resSm[1]);
@@ -109,15 +105,8 @@ WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int no
                _tcV[index] = y*dy;
                _tcW[index] = z*dz;
                _tcTemp[index] = 0.;
-               _energy[index] = 0.;
                }
        
-       // allocate eigenvalue arrays
-       _eigMin = new float[_totalCellsSm];
-       _eigMax = new float[_totalCellsSm];
-       for(int i=0; i < _totalCellsSm; i++) 
-               _eigMin[i] = _eigMax[i] =  0.;
-       
        // noise tiles
        _noiseTile = new float[noiseTileSize * noiseTileSize * noiseTileSize];
        /*
@@ -143,12 +132,7 @@ WTURBULENCE::~WTURBULENCE() {
   delete[] _tcW;
   delete[] _tcTemp;
 
-  delete[] _eigMin;
-  delete[] _eigMax;
   delete[] _noiseTile;
-
-  delete[] _energy;
-  delete[] _highFreqEnergy;
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -293,34 +277,34 @@ static float minDz(int x, int y, int z, float* input, Vec3Int res)
 // handle texture coordinates (advection, reset, eigenvalues), 
 // Beware -- uses big density maccormack as temporary arrays
 ////////////////////////////////////////////////////////////////////// 
-void WTURBULENCE::advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *_tempBig1, float *_tempBig2) {
+void WTURBULENCE::advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *tempBig1, float *tempBig2) {
   // advection
   SWAP_POINTERS(_tcTemp, _tcU);
   FLUID_3D::copyBorderX(_tcTemp, _resSm);
   FLUID_3D::copyBorderY(_tcTemp, _resSm);
   FLUID_3D::copyBorderZ(_tcTemp, _resSm);
   FLUID_3D::advectFieldMacCormack(dtOrg, xvel, yvel, zvel, 
-      _tcTemp, _tcU, _tempBig1, _tempBig2, _resSm, NULL);
+      _tcTemp, _tcU, tempBig1, tempBig2, _resSm, NULL);
 
   SWAP_POINTERS(_tcTemp, _tcV);
   FLUID_3D::copyBorderX(_tcTemp, _resSm);
   FLUID_3D::copyBorderY(_tcTemp, _resSm);
   FLUID_3D::copyBorderZ(_tcTemp, _resSm);
   FLUID_3D::advectFieldMacCormack(dtOrg, xvel, yvel, zvel, 
-      _tcTemp, _tcV, _tempBig1, _tempBig2, _resSm, NULL);
+      _tcTemp, _tcV, tempBig1, tempBig2, _resSm, NULL);
 
   SWAP_POINTERS(_tcTemp, _tcW);
   FLUID_3D::copyBorderX(_tcTemp, _resSm);
   FLUID_3D::copyBorderY(_tcTemp, _resSm);
   FLUID_3D::copyBorderZ(_tcTemp, _resSm);
   FLUID_3D::advectFieldMacCormack(dtOrg, xvel, yvel, zvel, 
-      _tcTemp, _tcW, _tempBig1, _tempBig2, _resSm, NULL);
+      _tcTemp, _tcW, tempBig1, tempBig2, _resSm, NULL);
 }
 
 //////////////////////////////////////////////////////////////////////
 // Compute the eigenvalues of the advected texture
 ////////////////////////////////////////////////////////////////////// 
-void WTURBULENCE::computeEigenvalues() {
+void WTURBULENCE::computeEigenvalues(float *_eigMin, float *_eigMax) {
   // stats
   float maxeig = -1.;
   float mineig = 10.;
@@ -365,7 +349,7 @@ void WTURBULENCE::computeEigenvalues() {
 //////////////////////////////////////////////////////////////////////
 // advect & reset texture coordinates based on eigenvalues
 ////////////////////////////////////////////////////////////////////// 
-void WTURBULENCE::resetTextureCoordinates() 
+void WTURBULENCE::resetTextureCoordinates(float *_eigMin, float *_eigMax
 {
   // allowed deformation of the textures
   const float limit = 2.f;
@@ -396,7 +380,7 @@ void WTURBULENCE::resetTextureCoordinates()
 // Compute the highest frequency component of the wavelet
 // decomposition
 //////////////////////////////////////////////////////////////////////
-void WTURBULENCE::decomposeEnergy()
+void WTURBULENCE::decomposeEnergy(float *_energy, float *_highFreqEnergy)
 {
   // do the decomposition -- the goal here is to have
   // the energy with the high frequency component stomped out
@@ -432,7 +416,7 @@ void WTURBULENCE::decomposeEnergy()
 // compute velocity from energies and march into obstacles
 // for wavelet decomposition
 ////////////////////////////////////////////////////////////////////// 
-void WTURBULENCE::computeEnergy(float* xvel, float* yvel, float* zvel, unsigned char *obstacles) 
+void WTURBULENCE::computeEnergy(float *_energy, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) 
 {
   // compute everywhere
   for (int x = 0; x < _totalCellsSm; x++) 
@@ -491,7 +475,7 @@ void WTURBULENCE::computeEnergy(float* xvel, float* yvel, float* zvel, unsigned
             }
             if (valid > 0)
             {
-              _energy[index] = sum / valid;
+              _energy[index] = sum / (float)valid;
               obstacles[index] = MARCHED;
             }
           }
@@ -516,9 +500,9 @@ void WTURBULENCE::computeEnergy(float* xvel, float* yvel, float* zvel, unsigned
 Vec3 WTURBULENCE::WVelocity(Vec3 orgPos)
 {
   // arbitrarily offset evaluation points
-  const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2,0,0);
-  const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2,0);
-  const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2);
+  const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0);
+  const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0);
+  const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0);
 
   const float f1y = WNoiseDy(p1, _noiseTile);
   const float f1z = WNoiseDz(p1, _noiseTile);
@@ -542,9 +526,9 @@ Vec3 WTURBULENCE::WVelocity(Vec3 orgPos)
 Vec3 WTURBULENCE::WVelocityWithJacobian(Vec3 orgPos, float* xUnwarped, float* yUnwarped, float* zUnwarped)
 {
   // arbitrarily offset evaluation points
-  const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2,0,0);
-  const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2,0);
-  const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2);
+  const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0);
+  const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0);
+  const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0);
 
   Vec3 final;
   final[0] = WNoiseDx(p1, _noiseTile);
@@ -575,44 +559,40 @@ Vec3 WTURBULENCE::WVelocityWithJacobian(Vec3 orgPos, float* xUnwarped, float* yU
   return ret;
 }
 
+
 //////////////////////////////////////////////////////////////////////
 // perform an actual noise advection step
 //////////////////////////////////////////////////////////////////////
 void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles) 
 {
-  // big velocity macCormack fields
-       float* _bigUx;
-       float* _bigUy;
-       float* _bigUz;
-
-       // temp arrays for BFECC and MacCormack - they have more convenient
-       // names in the actual implementations
-       float* _tempBig1;
-       float* _tempBig2;
-
-       // allocate high resolution velocity field. Note that this is only
-       // necessary because we use MacCormack advection. For semi-Lagrangian
-       // advection, these arrays are not necessary.
-       _tempBig1 = new float[_totalCellsBig];
-       _tempBig2 = new float[_totalCellsBig];
-
-  // enlarge timestep to match grid
-  const float dt = dtOrg * _amplify;
-  const float invAmp = 1.0f / _amplify;
-
-  _bigUx = new float[_totalCellsBig];
-  _bigUy = new float[_totalCellsBig];
-  _bigUz = new float[_totalCellsBig]; 
-
-  // prepare textures
-  advectTextureCoordinates(dtOrg, xvel,yvel,zvel, _tempBig1, _tempBig2);
+       // enlarge timestep to match grid
+       const float dt = dtOrg * _amplify;
+       const float invAmp = 1.0f / _amplify;
+       float *tempBig1 = new float[_totalCellsBig];
+       float *tempBig2 = new float[_totalCellsBig];
+       float *bigUx = new float[_totalCellsBig];
+       float *bigUy = new float[_totalCellsBig];
+       float *bigUz = new float[_totalCellsBig]; 
+       float *_energy = new float[_totalCellsSm];
+       float *highFreqEnergy = new float[_totalCellsSm];
+       float *eigMin  = new float[_totalCellsSm];
+       float *eigMax  = new float[_totalCellsSm];
+
+       memset(tempBig1, 0, sizeof(float)*_totalCellsBig);
+       memset(tempBig2, 0, sizeof(float)*_totalCellsBig);
+       memset(highFreqEnergy, 0, sizeof(float)*_totalCellsSm);
+       memset(eigMin, 0, sizeof(float)*_totalCellsSm);
+       memset(eigMax, 0, sizeof(float)*_totalCellsSm);
+
+       // prepare textures
+       advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2);
 
   // compute eigenvalues of the texture coordinates
-  computeEigenvalues();
+  computeEigenvalues(eigMin, eigMax);
 
   // do wavelet decomposition of energy
-  computeEnergy(xvel, yvel, zvel, obstacles);
-  decomposeEnergy();
+  computeEnergy(_energy, xvel, yvel, zvel, obstacles);
+  decomposeEnergy(_energy, highFreqEnergy);
 
   // zero out coefficients inside of the obstacle
   for (int x = 0; x < _totalCellsSm; x++)
@@ -647,7 +627,7 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
 
         // retrieve wavelet energy at highest frequency
         float energy = INTERPOLATE::lerp3d(
-            _highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
+            highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
 
         // base amplitude for octave 0
         float coefficient = sqrtf(2.0f * fabs(energy));
@@ -656,8 +636,8 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
         // add noise to velocity, but only if the turbulence is
         // sufficiently undeformed, and the energy is large enough
         // to make a difference
-        const bool addNoise = _eigMax[indexSmall] < 2. && 
-                              _eigMin[indexSmall] > 0.5;
+        const bool addNoise = eigMax[indexSmall] < 2. && 
+                              eigMin[indexSmall] > 0.5;
         if (addNoise && amplitude > _cullingThreshold) {
           // base amplitude for octave 0
           float amplitudeScaled = amplitude;
@@ -679,21 +659,21 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
         // If you wanted to save memory, you would instead perform a 
         // semi-Lagrangian backtrace for the current grid cell here. Then
         // you could just throw the velocity away.
-        _bigUx[index] = vel[0];
-        _bigUy[index] = vel[1];
-        _bigUz[index] = vel[2];
+        bigUx[index] = vel[0];
+        bigUy[index] = vel[1];
+        bigUz[index] = vel[2];
 
         // compute the velocity magnitude for substepping later
-        const float velMag = _bigUx[index] * _bigUx[index] + 
-                             _bigUy[index] * _bigUy[index] + 
-                             _bigUz[index] * _bigUz[index];
+        const float velMag = bigUx[index] * bigUx[index] + 
+                             bigUy[index] * bigUy[index] + 
+                             bigUz[index] * bigUz[index];
         if (velMag > maxVelocity) maxVelocity = velMag;
 
         // zero out velocity inside obstacles
         float obsCheck = INTERPOLATE::lerp3dToFloat(
             obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); 
         if (obsCheck > 0.95) 
-          _bigUx[index] = _bigUy[index] = _bigUz[index] = 0.;
+          bigUx[index] = bigUy[index] = bigUz[index] = 0.;
       }
 
   // prepare density for an advection
@@ -701,24 +681,23 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
 
   // based on the maximum velocity present, see if we need to substep,
   // but cap the maximum number of substeps to 5
-  const int maxSubSteps = 25;
-  const int maxVel = 5;
+  const int maxSubSteps = 5; 
   maxVelocity = sqrt(maxVelocity) * dt;
-  int totalSubsteps = (int)(maxVelocity / (float)maxVel);
+  int totalSubsteps = (int)(maxVelocity / (float)maxSubSteps);
   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
   const float dtSubdiv = dt / (float)totalSubsteps;
 
   // set boundaries of big velocity grid
-  FLUID_3D::setZeroX(_bigUx, _resBig); 
-  FLUID_3D::setZeroY(_bigUy, _resBig); 
-  FLUID_3D::setZeroZ(_bigUz, _resBig);
+  FLUID_3D::setZeroX(bigUx, _resBig); 
+  FLUID_3D::setZeroY(bigUy, _resBig); 
+  FLUID_3D::setZeroZ(bigUz, _resBig);
 
   // do the MacCormack advection, with substepping if necessary
   for(int substep = 0; substep < totalSubsteps; substep++)
   {
-    FLUID_3D::advectFieldMacCormack(dtSubdiv, _bigUx, _bigUy, _bigUz, 
-        _densityBigOld, _densityBig, _tempBig1, _tempBig2, _resBig, NULL);
+    FLUID_3D::advectFieldMacCormack(dtSubdiv, bigUx, bigUy, bigUz, 
+        _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL);
 
     if (substep < totalSubsteps - 1) 
       SWAP_POINTERS(_densityBig, _densityBigOld);
@@ -730,25 +709,21 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
   // reset texture coordinates now in preparation for next timestep
   // Shouldn't do this before generating the noise because then the 
   // eigenvalues stored do not reflect the underlying texture coordinates
-  resetTextureCoordinates();
+  resetTextureCoordinates(eigMin, eigMax);
   
-  // output files
-  /*
-  string prefix = string("./amplified.preview/density_bigxy_");
-  FLUID_3D::writeImageSliceXY(_densityBig, _resBig, _resBig[2]/2, prefix, _totalStepsBig, 1.0f);
-  //string df3Prefix = string("./df3/density_big_");
-  //IMAGE::dumpDF3(_totalStepsBig, df3Prefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
-  string pbrtPrefix = string("./pbrt/density_big_");
-  IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
-       */
-  _totalStepsBig++;
+  delete[] tempBig1;
+  delete[] tempBig2;
+  delete[] bigUx;
+  delete[] bigUy;
+  delete[] bigUz;
+  delete[] _energy;
+  delete[] highFreqEnergy;
 
-       delete[] _bigUx;
-       delete[] _bigUy;
-       delete[] _bigUz;
+  delete[] eigMin;
+  delete[] eigMax;
+  
 
-       delete[] _tempBig1;
-    delete[] _tempBig2;
+  _totalStepsBig++;
 }
 
 //////////////////////////////////////////////////////////////////////
@@ -757,257 +732,258 @@ void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel,
 //////////////////////////////////////////////////////////////////////
 void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles)
 {
-       // big velocity macCormack fields
-       float* _bigUx;
-       float* _bigUy;
-       float* _bigUz;
-
-       // temp arrays for BFECC and MacCormack - they have more convenient
-       // names in the actual implementations
-       float* _tempBig1;
-       float* _tempBig2;
-
-       // allocate high resolution velocity field. Note that this is only
-       // necessary because we use MacCormack advection. For semi-Lagrangian
-       // advection, these arrays are not necessary.
-       _tempBig1 = new float[_totalCellsBig];
-       _tempBig2 = new float[_totalCellsBig];
-
        // enlarge timestep to match grid
        const float dt = dtOrg * _amplify;
        const float invAmp = 1.0f / _amplify;
-
-       _bigUx = new float[_totalCellsBig];
-       _bigUy = new float[_totalCellsBig];
-       _bigUz = new float[_totalCellsBig]; 
+       float *tempBig1 = new float[_totalCellsBig];
+       float *tempBig2 = new float[_totalCellsBig];
+       float *bigUx = new float[_totalCellsBig];
+       float *bigUy = new float[_totalCellsBig];
+       float *bigUz = new float[_totalCellsBig]; 
+       float *_energy = new float[_totalCellsSm];
+       float *highFreqEnergy = new float[_totalCellsSm];
+       float *eigMin  = new float[_totalCellsSm];
+       float *eigMax  = new float[_totalCellsSm];
+
+       memset(highFreqEnergy, 0, sizeof(float)*_totalCellsSm);
+       memset(eigMin, 0, sizeof(float)*_totalCellsSm);
+       memset(eigMax, 0, sizeof(float)*_totalCellsSm);
 
        // prepare textures
-       advectTextureCoordinates(dtOrg, xvel,yvel,zvel, _tempBig1, _tempBig2);
+       advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2);
 
        // do wavelet decomposition of energy
-       computeEnergy(xvel, yvel, zvel, obstacles);
-       decomposeEnergy();
+       computeEnergy(_energy, xvel, yvel, zvel, obstacles);
 
-       // zero out coefficients inside of the obstacle
        for (int x = 0; x < _totalCellsSm; x++)
                if (obstacles[x]) _energy[x] = 0.f;
 
-       // parallel region setup
-       float maxVelMagThreads[8] = { -1., -1., -1., -1., -1., -1., -1., -1. };
-
-       #if PARALLEL==1
-               #pragma omp parallel
-       #endif
-         { float maxVelMag1 = 0.;
-       #if PARALLEL==1
-               const int id  = omp_get_thread_num(); /*, num = omp_get_num_threads(); */
-       #endif
-
-       // vector noise main loop
-       #if PARALLEL==1
-       #pragma omp for  schedule(static)
-       #endif
-       for (int zSmall = 0; zSmall < _zResSm; zSmall++) 
-       for (int ySmall = 0; ySmall < _yResSm; ySmall++) 
-       for (int xSmall = 0; xSmall < _xResSm; xSmall++)
-       {
-               const int indexSmall = xSmall + ySmall * _xResSm + zSmall * _slabSizeSm;
-
-               // compute jacobian
-               float jacobian[3][3] = {
-                 { minDx(xSmall, ySmall, zSmall, _tcU, _resSm), minDx(xSmall, ySmall, zSmall, _tcV, _resSm), minDx(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
-                 { minDy(xSmall, ySmall, zSmall, _tcU, _resSm), minDy(xSmall, ySmall, zSmall, _tcV, _resSm), minDy(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
-                 { minDz(xSmall, ySmall, zSmall, _tcU, _resSm), minDz(xSmall, ySmall, zSmall, _tcV, _resSm), minDz(xSmall, ySmall, zSmall, _tcW, _resSm) }
-               };
-
-               // get LU factorization of texture jacobian and apply 
-               // it to unit vectors
-               JAMA::LU<float> LU = computeLU3x3(jacobian);
-               float xUnwarped[] = {1.0f, 0.0f, 0.0f};
-               float yUnwarped[] = {0.0f, 1.0f, 0.0f};
-               float zUnwarped[] = {0.0f, 0.0f, 1.0f};
-               float xWarped[] = {1.0f, 0.0f, 0.0f};
-               float yWarped[] = {0.0f, 1.0f, 0.0f};
-               float zWarped[] = {0.0f, 0.0f, 1.0f};
-               bool nonSingular = LU.isNonsingular();
-
-               #if 0
-                       // UNUSED
-                       float eigMax = 10.0f;
-                       float eigMin = 0.1f;
-               #endif
+       decomposeEnergy(_energy, highFreqEnergy);
 
-               if (nonSingular)
-               {
-                       solveLU3x3(LU, xUnwarped, xWarped);
-                       solveLU3x3(LU, yUnwarped, yWarped);
-                       solveLU3x3(LU, zUnwarped, zWarped);
-
-                       // compute the eigenvalues while we have the Jacobian available
-                       Vec3 eigenvalues = Vec3(1.);
-                       computeEigenvalues3x3( &eigenvalues[0], jacobian);
-                       _eigMax[indexSmall] = MAX3V(eigenvalues);
-                       _eigMin[indexSmall] = MIN3V(eigenvalues);
-               }
+       // zero out coefficients inside of the obstacle
+       for (int x = 0; x < _totalCellsSm; x++)
+               if (obstacles[x]) highFreqEnergy[x] = 0.f;
 
-               // make sure to skip one on the beginning and end
-               int xStart = (xSmall == 0) ? 1 : 0;
-               int xEnd   = (xSmall == _xResSm - 1) ? _amplify - 1 : _amplify;
-               int yStart = (ySmall == 0) ? 1 : 0;
-               int yEnd   = (ySmall == _yResSm - 1) ? _amplify - 1 : _amplify;
-               int zStart = (zSmall == 0) ? 1 : 0;
-               int zEnd   = (zSmall == _zResSm - 1) ? _amplify - 1 : _amplify;
-                 
-               for (int zBig = zStart; zBig < zEnd; zBig++) 
-               for (int yBig = yStart; yBig < yEnd; yBig++) 
-               for (int xBig = xStart; xBig < xEnd; xBig++)
-               {
-                       const int x = xSmall * _amplify + xBig;
-                       const int y = ySmall * _amplify + yBig;
-                       const int z = zSmall * _amplify + zBig;
-
-                       // get unit position for both fine and coarse grid
-                       const Vec3 pos = Vec3(x,y,z);
-                       const Vec3 posSm = pos * invAmp;
-
-                       // get grid index for both fine and coarse grid
-                       const int index = x + y *_xResBig + z *_slabSizeBig;
-
-                       // get a linearly interpolated velocity and texcoords
-                       // from the coarse grid
-                       Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel, 
-                         posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
-                       Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW, 
-                         posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
-
-                       // multiply the texture coordinate by _resSm so that turbulence
-                       // synthesis begins at the first octave that the coarse grid 
-                       // cannot capture
-                       Vec3 texCoord = Vec3(uvw[0] * _resSm[0], 
-                                                          uvw[1] * _resSm[1],
-                                                          uvw[2] * _resSm[2]); 
-
-                       // retrieve wavelet energy at highest frequency
-                       float energy = INTERPOLATE::lerp3d(
-                         _highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
-
-                       // base amplitude for octave 0
-                       float coefficient = sqrtf(2.0f * fabs(energy));
-                       const float amplitude = *_strength * fabs(0.5 * coefficient) * persistence;
-
-                       // add noise to velocity, but only if the turbulence is
-                       // sufficiently undeformed, and the energy is large enough
-                       // to make a difference
-                       const bool addNoise = _eigMax[indexSmall] < 2. && 
-                                                               _eigMin[indexSmall] > 0.5;
-
-                       if (addNoise && amplitude > _cullingThreshold) {
-                               // base amplitude for octave 0
-                               float amplitudeScaled = amplitude;
-
-                               for (int octave = 0; octave < _octaves; octave++)
-                               {
-                                       // multiply the vector noise times the maximum allowed
-                                       // noise amplitude at this octave, and add it to the total
-                                       vel += WVelocityWithJacobian(texCoord, &xUnwarped[0], &yUnwarped[0], &zUnwarped[0]) * amplitudeScaled;
-
-                                       // scale coefficient for next octave
-                                       amplitudeScaled *= persistence;
-                                       texCoord *= 2.0f;
-                               }
-                       }
-
-                       // Store velocity + turbulence in big grid for maccormack step
-                       //
-                       // If you wanted to save memory, you would instead perform a 
-                       // semi-Lagrangian backtrace for the current grid cell here. Then
-                       // you could just throw the velocity away.
-                       _bigUx[index] = vel[0];
-                       _bigUy[index] = vel[1];
-                       _bigUz[index] = vel[2];
-
-                       // compute the velocity magnitude for substepping later
-                       const float velMag = _bigUx[index] * _bigUx[index] + 
-                                                          _bigUy[index] * _bigUy[index] + 
-                                                          _bigUz[index] * _bigUz[index];
-                       if (velMag > maxVelMag1) maxVelMag1 = velMag;
-
-                       // zero out velocity inside obstacles
-                       float obsCheck = INTERPOLATE::lerp3dToFloat(
-                         obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); 
-
-                       if (obsCheck > 0.95) 
-                               _bigUx[index] = _bigUy[index] = _bigUz[index] = 0.;
-               } // xyz
-
-               #if PARALLEL==1
-               maxVelMagThreads[id] = maxVelMag1;
-               #else
-               maxVelMagThreads[0] = maxVelMag1;
-               #endif
-               }
-       } // omp
-
-       // compute maximum over threads
-       float maxVelMag = maxVelMagThreads[0];
-       #if PARALLEL==1
-       for (int i = 1; i < 8; i++) 
-       if (maxVelMag < maxVelMagThreads[i]) 
-         maxVelMag = maxVelMagThreads[i];
-       #endif
-
-       // prepare density for an advection
-       SWAP_POINTERS(_densityBig, _densityBigOld);
-
-       // based on the maximum velocity present, see if we need to substep,
-       // but cap the maximum number of substeps to 5
-       const int maxSubSteps = 25;
-       const int maxVel = 5;
-       maxVelMag = sqrt(maxVelMag) * dt;
-       int totalSubsteps = (int)(maxVelMag / (float)maxVel);
-       totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
-       // printf("totalSubsteps: %d\n", totalSubsteps);
-       totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
-       const float dtSubdiv = dt / (float)totalSubsteps;
-
-       // set boundaries of big velocity grid
-       FLUID_3D::setZeroX(_bigUx, _resBig); 
-       FLUID_3D::setZeroY(_bigUy, _resBig); 
-       FLUID_3D::setZeroZ(_bigUz, _resBig);
-
-       // do the MacCormack advection, with substepping if necessary
-       for(int substep = 0; substep < totalSubsteps; substep++)
-       {
-       FLUID_3D::advectFieldMacCormack(dtSubdiv, _bigUx, _bigUy, _bigUz, 
-               _densityBigOld, _densityBig, _tempBig1, _tempBig2, _resBig, NULL);
+       Vec3Int ressm(_xResSm, _yResSm, _zResSm);
+       FLUID_3D::setNeumannX(highFreqEnergy, ressm);
+       FLUID_3D::setNeumannY(highFreqEnergy, ressm);
+       FLUID_3D::setNeumannZ(highFreqEnergy, ressm);
+
+  // parallel region setup
+  float maxVelMagThreads[8] = { -1., -1., -1., -1., -1., -1., -1., -1. };
+#if PARALLEL==1
+#pragma omp parallel
+#endif
+  { float maxVelMag1 = 0.;
+#if PARALLEL==1
+    const int id  = omp_get_thread_num(); /*, num = omp_get_num_threads(); */
+#endif
+
+  // vector noise main loop
+#if PARALLEL==1
+#pragma omp for  schedule(static)
+#endif
+  for (int zSmall = 0; zSmall < _zResSm; zSmall++) 
+  for (int ySmall = 0; ySmall < _yResSm; ySmall++) 
+  for (int xSmall = 0; xSmall < _xResSm; xSmall++)
+  {
+    const int indexSmall = xSmall + ySmall * _xResSm + zSmall * _slabSizeSm;
+
+    // compute jacobian
+    float jacobian[3][3] = {
+      { minDx(xSmall, ySmall, zSmall, _tcU, _resSm), minDx(xSmall, ySmall, zSmall, _tcV, _resSm), minDx(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
+      { minDy(xSmall, ySmall, zSmall, _tcU, _resSm), minDy(xSmall, ySmall, zSmall, _tcV, _resSm), minDy(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
+      { minDz(xSmall, ySmall, zSmall, _tcU, _resSm), minDz(xSmall, ySmall, zSmall, _tcV, _resSm), minDz(xSmall, ySmall, zSmall, _tcW, _resSm) }
+    };
+
+    // get LU factorization of texture jacobian and apply 
+    // it to unit vectors
+    JAMA::LU<float> LU = computeLU3x3(jacobian);
+    float xUnwarped[] = {1.0f, 0.0f, 0.0f};
+    float yUnwarped[] = {0.0f, 1.0f, 0.0f};
+    float zUnwarped[] = {0.0f, 0.0f, 1.0f};
+    float xWarped[] = {1.0f, 0.0f, 0.0f};
+    float yWarped[] = {0.0f, 1.0f, 0.0f};
+    float zWarped[] = {0.0f, 0.0f, 1.0f};
+    bool nonSingular = LU.isNonsingular();
+#if 0
+       // UNUSED
+    float eigMax = 10.0f;
+    float eigMin = 0.1f;
+#endif
+    if (nonSingular)
+    {
+      solveLU3x3(LU, xUnwarped, xWarped);
+      solveLU3x3(LU, yUnwarped, yWarped);
+      solveLU3x3(LU, zUnwarped, zWarped);
+
+      // compute the eigenvalues while we have the Jacobian available
+      Vec3 eigenvalues = Vec3(1.);
+      computeEigenvalues3x3( &eigenvalues[0], jacobian);
+      eigMax[indexSmall] = MAX3V(eigenvalues);
+      eigMin[indexSmall] = MIN3V(eigenvalues);
+    }
+    
+    // make sure to skip one on the beginning and end
+    int xStart = (xSmall == 0) ? 1 : 0;
+    int xEnd   = (xSmall == _xResSm - 1) ? _amplify - 1 : _amplify;
+    int yStart = (ySmall == 0) ? 1 : 0;
+    int yEnd   = (ySmall == _yResSm - 1) ? _amplify - 1 : _amplify;
+    int zStart = (zSmall == 0) ? 1 : 0;
+    int zEnd   = (zSmall == _zResSm - 1) ? _amplify - 1 : _amplify;
+      
+    for (int zBig = zStart; zBig < zEnd; zBig++) 
+    for (int yBig = yStart; yBig < yEnd; yBig++) 
+    for (int xBig = xStart; xBig < xEnd; xBig++)
+    {
+      const int x = xSmall * _amplify + xBig;
+      const int y = ySmall * _amplify + yBig;
+      const int z = zSmall * _amplify + zBig;
+      
+      // get unit position for both fine and coarse grid
+      const Vec3 pos = Vec3(x,y,z);
+      const Vec3 posSm = pos * invAmp;
+      
+      // get grid index for both fine and coarse grid
+      const int index = x + y *_xResBig + z *_slabSizeBig;
+      
+      // get a linearly interpolated velocity and texcoords
+      // from the coarse grid
+      Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel, 
+          posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
+      Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW, 
+          posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
+
+      // multiply the texture coordinate by _resSm so that turbulence
+      // synthesis begins at the first octave that the coarse grid 
+      // cannot capture
+      Vec3 texCoord = Vec3(uvw[0] * _resSm[0], 
+                           uvw[1] * _resSm[1],
+                           uvw[2] * _resSm[2]); 
+
+      // retrieve wavelet energy at highest frequency
+      float energy = INTERPOLATE::lerp3d(
+          highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
+
+      // base amplitude for octave 0
+      float coefficient = sqrtf(2.0f * fabs(energy));
+      const float amplitude = *_strength * fabs(0.5 * coefficient) * persistence;
+
+      // add noise to velocity, but only if the turbulence is
+      // sufficiently undeformed, and the energy is large enough
+      // to make a difference
+      const bool addNoise = eigMax[indexSmall] < 2. && 
+                            eigMin[indexSmall] > 0.5;
+      if (addNoise && amplitude > _cullingThreshold) {
+        // base amplitude for octave 0
+        float amplitudeScaled = amplitude;
 
-       if (substep < totalSubsteps - 1) 
-         SWAP_POINTERS(_densityBig, _densityBigOld);
-       } // substep
+        for (int octave = 0; octave < _octaves; octave++)
+        {
+          // multiply the vector noise times the maximum allowed
+          // noise amplitude at this octave, and add it to the total
+          vel += WVelocityWithJacobian(texCoord, &xUnwarped[0], &yUnwarped[0], &zUnwarped[0]) * amplitudeScaled;
 
-       // wipe the density borders
-       FLUID_3D::setZeroBorder(_densityBig, _resBig);
+          // scale coefficient for next octave
+          amplitudeScaled *= persistence;
+          texCoord *= 2.0f;
+        }
+      }
 
-       // reset texture coordinates now in preparation for next timestep
-       // Shouldn't do this before generating the noise because then the 
-       // eigenvalues stored do not reflect the underlying texture coordinates
-       resetTextureCoordinates();
+      // Store velocity + turbulence in big grid for maccormack step
+      //
+      // If you wanted to save memory, you would instead perform a 
+      // semi-Lagrangian backtrace for the current grid cell here. Then
+      // you could just throw the velocity away.
+      bigUx[index] = vel[0];
+      bigUy[index] = vel[1];
+      bigUz[index] = vel[2];
+
+      // compute the velocity magnitude for substepping later
+      const float velMag = bigUx[index] * bigUx[index] + 
+                           bigUy[index] * bigUy[index] + 
+                           bigUz[index] * bigUz[index];
+      if (velMag > maxVelMag1) maxVelMag1 = velMag;
+
+      // zero out velocity inside obstacles
+      float obsCheck = INTERPOLATE::lerp3dToFloat(
+          obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm); 
+      if (obsCheck > 0.95) 
+        bigUx[index] = bigUy[index] = bigUz[index] = 0.;
+    } // xyz
 
-       // output files
-       // string prefix = string("./amplified.preview/density_bigxy_");
-       // FLUID_3D::writeImageSliceXY(_densityBig, _resBig, _resBig[2]/2, prefix, _totalStepsBig, 1.0f);
-       //string df3prefix = string("./df3/density_big_");
-       //IMAGE::dumpDF3(_totalStepsBig, df3prefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
-       // string pbrtPrefix = string("./pbrt/density_big_");
-       // IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
+#if PARALLEL==1
+    maxVelMagThreads[id] = maxVelMag1;
+#else
+    maxVelMagThreads[0] = maxVelMag1;
+#endif
+  }
+  } // omp
+  
+  // compute maximum over threads
+  float maxVelMag = maxVelMagThreads[0];
+#if PARALLEL==1
+  for (int i = 1; i < 8; i++) 
+    if (maxVelMag < maxVelMagThreads[i]) 
+      maxVelMag = maxVelMagThreads[i];
+#endif
 
-       _totalStepsBig++;
+  // prepare density for an advection
+  SWAP_POINTERS(_densityBig, _densityBigOld);
 
-       delete[] _bigUx;
-       delete[] _bigUy;
-       delete[] _bigUz;
+  // based on the maximum velocity present, see if we need to substep,
+  // but cap the maximum number of substeps to 5
+  const int maxSubSteps = 25;
+  const int maxVel = 5;
+  maxVelMag = sqrt(maxVelMag) * dt;
+  int totalSubsteps = (int)(maxVelMag / (float)maxVel);
+  totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
+  // printf("totalSubsteps: %d\n", totalSubsteps);
+  totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
+  const float dtSubdiv = dt / (float)totalSubsteps;
 
-       delete[] _tempBig1;
-       delete[] _tempBig2;
+  // set boundaries of big velocity grid
+  FLUID_3D::setZeroX(bigUx, _resBig); 
+  FLUID_3D::setZeroY(bigUy, _resBig); 
+  FLUID_3D::setZeroZ(bigUz, _resBig);
+
+  // do the MacCormack advection, with substepping if necessary
+  for(int substep = 0; substep < totalSubsteps; substep++)
+  {
+    FLUID_3D::advectFieldMacCormack(dtSubdiv, bigUx, bigUy, bigUz, 
+        _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL);
+
+    if (substep < totalSubsteps - 1) 
+      SWAP_POINTERS(_densityBig, _densityBigOld);
+  } // substep
+
+  delete[] tempBig1;
+  delete[] tempBig2;
+  delete[] bigUx;
+  delete[] bigUy;
+  delete[] bigUz;
+  delete[] _energy;
+  delete[] highFreqEnergy;
+  
+  // wipe the density borders
+  FLUID_3D::setZeroBorder(_densityBig, _resBig);
+    
+  // reset texture coordinates now in preparation for next timestep
+  // Shouldn't do this before generating the noise because then the 
+  // eigenvalues stored do not reflect the underlying texture coordinates
+  resetTextureCoordinates(eigMin, eigMax);
+
+  delete[] eigMin;
+  delete[] eigMax;
+  
+  // output files
+  // string prefix = string("./amplified.preview/density_bigxy_");
+  // FLUID_3D::writeImageSliceXY(_densityBig, _resBig, _resBig[2]/2, prefix, _totalStepsBig, 1.0f);
+  //string df3prefix = string("./df3/density_big_");
+  //IMAGE::dumpDF3(_totalStepsBig, df3prefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
+  // string pbrtPrefix = string("./pbrt/density_big_");
+  // IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
+  
+  _totalStepsBig++;
 }
 
index c21e002ad4873f128bfbd668cd5031eb28f07c06..0aa978e9e52fcdce586b6fbd85d436ed39389b8d 100644 (file)
@@ -49,10 +49,10 @@ class WTURBULENCE
                void stepTurbulenceFull(float dt, float* xvel, float* yvel, float* zvel, unsigned char *obstacles);
        
                // texcoord functions
-               void advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *_tempBig1, float *_tempBig2);
-               void resetTextureCoordinates();
+               void advectTextureCoordinates(float dtOrg, float* xvel, float* yvel, float* zvel, float *tempBig1, float *tempBig2);
+               void resetTextureCoordinates(float *_eigMin, float *_eigMax);
 
-               void computeEnergy(float* xvel, float* yvel, float* zvel, unsigned char *obstacles);
+               void computeEnergy(float *energy, float* xvel, float* yvel, float* zvel, unsigned char *obstacles);
 
                // evaluate wavelet noise function
                Vec3 WVelocity(Vec3 p);
@@ -63,8 +63,6 @@ class WTURBULENCE
                inline float* getArrayTcU() { return _tcU; }
                inline float* getArrayTcV() { return _tcV; }
                inline float* getArrayTcW() { return _tcW; }
-               inline float* getArrayEigMin() { return _eigMin; }
-               inline float* getArrayEigMax() { return _eigMax; }
 
                inline Vec3Int getResSm() { return _resSm; } // small resolution
                inline Vec3Int getResBig() { return _resBig; }
@@ -81,15 +79,15 @@ class WTURBULENCE
                
                // noise settings
                float _cullingThreshold;
-               float _noiseStrength;
-               float _noiseSizeScale;
-               bool _uvwAdvection;
-               bool _uvwReset;
-               float _noiseTimeanimSpeed;
-               int _dumpInterval;
-               int _noiseControlType;
+               // float _noiseStrength;
+               // float _noiseSizeScale;
+               // bool _uvwAdvection;
+               // bool _uvwReset;
+               // float _noiseTimeanimSpeed;
+               // int _dumpInterval;
+               // nt _noiseControlType;
                // debug, scale density for projections output images
-               float _outputScale;
+               // float _outputScale;
 
                // noise resolution
                int _xResBig;
@@ -117,24 +115,15 @@ class WTURBULENCE
                float* _tcW;
                float* _tcTemp;
 
-               float* _eigMin; // no save -dg
-               float* _eigMax; // no save -dg
-
-               // wavelet decomposition of velocity energies
-               float* _energy; // no save -dg
-
                // noise data
                float* _noiseTile;
                //float* _noiseTileExt;
 
                // step counter
                int _totalStepsBig;
-
-               // highest frequency component of wavelet decomposition
-               float* _highFreqEnergy; // no save -dg
                
-               void computeEigenvalues();
-               void decomposeEnergy();
+               void computeEigenvalues(float *_eigMin, float *_eigMax);
+               void decomposeEnergy(float *energy, float *_highFreqEnergy);
 };
 
 #endif // WTURBULENCE_H
index 5a016b51bbeb846ffa78da9aa49977fadc0e1335..058831dddbbe868975e25c82fa81eb6e9ba336e9 100644 (file)
@@ -137,7 +137,7 @@ extern "C" void smoke_dissolve(FLUID_3D *fluid, int speed, int log)
        }
 }
 
-extern "C" void smoke_turbulence_dissolve(WTURBULENCE *wt, int speed, int log)
+extern "C" void smoke_dissolve_wavelet(WTURBULENCE *wt, int speed, int log)
 {
        float *density = wt->getDensityBig();
        Vec3Int r = wt->getResBig();
@@ -172,7 +172,7 @@ extern "C" void smoke_turbulence_dissolve(WTURBULENCE *wt, int speed, int log)
        }
 }
 
-extern "C" void smoke_turbulence_initBlenderRNA(WTURBULENCE *wt, float *strength)
+extern "C" void smoke_initWaveletBlenderRNA(WTURBULENCE *wt, float *strength)
 {
        wt->initBlenderRNA(strength);
 }
@@ -241,13 +241,14 @@ extern "C" float *smoke_turbulence_get_density(WTURBULENCE *wt)
        return wt ? wt->getDensityBig() : NULL;
 }
 
-extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, unsigned int *res)
+extern "C" void smoke_turbulence_get_res(WTURBULENCE *wt, int *res)
 {
        if(wt)
        {
-               res[0] = wt->_resBig[0];
-               res[1] = wt->_resBig[1];
-               res[2] = wt->_resBig[2];
+               Vec3Int r = wt->getResBig();
+               res[0] = r[0];
+               res[1] = r[1];
+               res[2] = r[2];
        }
 }
 
index c87f71bff4249069634094ab528fa445442bbdd0..83c1ffc2e9afafe654b47e76519d2ab59d184830 100644 (file)
@@ -3,11 +3,6 @@ import bpy
 
 from buttons_particle import point_cache_ui
 
-def smoke_panel_enabled_low(smd):
-       if smd.smoke_type == 'TYPE_DOMAIN':
-               return smd.domain.point_cache.baked==False
-       return True
-
 class PhysicButtonsPanel(bpy.types.Panel):
        __space_type__ = 'PROPERTIES'
        __region_type__ = 'WINDOW'
@@ -45,8 +40,6 @@ class PHYSICS_PT_smoke(PhysicButtonsPanel):
                        split.itemL()
 
                if md:
-               
-                       # layout.enabled = smoke_panel_enabled(md)
                        layout.itemR(md, "smoke_type", expand=True)
                
                        if md.smoke_type == 'TYPE_DOMAIN':
@@ -66,7 +59,7 @@ class PHYSICS_PT_smoke(PhysicButtonsPanel):
                                col.itemR(domain, "dissolve_smoke", text="Dissolve")
                                sub = col.column()
                                sub.active = domain.dissolve_smoke
-                               sub.itemR(domain, "dissolve_speed", text="Speed")
+                               sub.itemR(domain, "dissolve_speed", text="Time")
                                sub.itemR(domain, "dissolve_smoke_log", text="Slow")
                                
                        elif md.smoke_type == 'TYPE_FLOW':
@@ -90,14 +83,17 @@ class PHYSICS_PT_smoke(PhysicButtonsPanel):
                                        
                        #elif md.smoke_type == 'TYPE_COLL':
                        #       layout.itemS()
-                                       
+
 class PHYSICS_PT_smoke_groups(PhysicButtonsPanel):
        __label__ = "Smoke Groups"
        __default_closed__ = True
        
        def poll(self, context):
                md = context.smoke
-               return md and (md.smoke_type == 'TYPE_DOMAIN')
+               if md:
+                               return (md.smoke_type == 'TYPE_DOMAIN')
+               
+               return False
 
        def draw(self, context):
                layout = self.layout
@@ -116,7 +112,7 @@ class PHYSICS_PT_smoke_groups(PhysicButtonsPanel):
                col = split.column()
                col.itemL(text="Collision Group:")
                col.itemR(group, "coll_group", text="")
-               
+
 class PHYSICS_PT_smoke_cache(PhysicButtonsPanel):
        __label__ = "Smoke Cache"
        __default_closed__ = True
@@ -129,7 +125,7 @@ class PHYSICS_PT_smoke_cache(PhysicButtonsPanel):
                layout = self.layout
 
                md = context.smoke.domain_settings
-               cache = md.point_cache
+               cache = md.point_cache_low
                        
                point_cache_ui(self, cache, cache.baked==False, 0, 1)
                                        
@@ -152,7 +148,7 @@ class PHYSICS_PT_smoke_highres(PhysicButtonsPanel):
                md = context.smoke.domain_settings
 
                split = layout.split()
-                       
+
                col = split.column()
                col.itemL(text="Resolution:")
                col.itemR(md, "amplify", text="Divisions")
@@ -161,66 +157,26 @@ class PHYSICS_PT_smoke_highres(PhysicButtonsPanel):
                col.itemL(text="Noise Method:")
                col.row().itemR(md, "noise_type", text="")
                col.itemR(md, "strength")
-               col.itemR(md, "show_highres")
+               sub.itemR(md, "viewhighres")
                
 class PHYSICS_PT_smoke_cache_highres(PhysicButtonsPanel):
-       __label__ = "Smoke Cache"
+       __label__ = "Smoke High Resolution Cache"
        __default_closed__ = True
 
        def poll(self, context):
-               return (context.smoke)
+               md = context.smoke
+               return md and (md.smoke_type == 'TYPE_DOMAIN') and md.domain_settings.highres
 
        def draw(self, context):
                layout = self.layout
 
-               md = context.smoke
-
-               cache = md.point_cache
-                       
-               layout.set_context_pointer("PointCache", cache)
-                       
-               row = layout.row()
-               row.template_list(cache, "point_cache_list", cache, "active_point_cache_index")
-               col = row.column(align=True)
-               col.itemO("ptcache.add_new", icon='ICON_ZOOMIN', text="")
-               col.itemO("ptcache.remove", icon='ICON_ZOOMOUT', text="")
-                       
-               row = layout.row()
-               row.itemR(cache, "name")
-                       
-               row = layout.row()
-               row.itemR(cache, "start_frame")
-               row.itemR(cache, "end_frame")
-                       
-               row = layout.row()
-                       
-               if cache.baked == True:
-                       row.itemO("ptcache.free_bake", text="Free Bake")
-               else:
-                       row.item_booleanO("ptcache.bake", "bake", True, text="Bake")
-                       
-               subrow = row.row()
-               subrow.enabled = cache.frames_skipped or cache.outdated
-               subrow.itemO("ptcache.bake", "bake", False, text="Calculate to Current Frame")
-                               
-               row = layout.row()
-               #row.enabled = smoke_panel_enabled(psys)
-               row.itemO("ptcache.bake_from_cache", text="Current Cache to Bake")
-               
-               row = layout.row()
-               #row.enabled = smoke_panel_enabled(psys)
-                       
-               layout.itemL(text=cache.info)
-                       
-               layout.itemS()
+               md = context.smoke.domain_settings
+               cache = md.point_cache_high
                        
-               row = layout.row()
-               row.itemO("ptcache.bake_all", "bake", True, text="Bake All Dynamics")
-               row.itemO("ptcache.free_bake_all", text="Free All Bakes")
-               layout.itemO("ptcache.bake_all", "bake", False, text="Update All Dynamics to current frame")
-
+               point_cache_ui(self, cache, cache.baked==False, 0, 1)
+                                       
 bpy.types.register(PHYSICS_PT_smoke)
 bpy.types.register(PHYSICS_PT_smoke_cache)
+bpy.types.register(PHYSICS_PT_smoke_highres)
 bpy.types.register(PHYSICS_PT_smoke_groups)
-#bpy.types.register(PHYSICS_PT_smoke_highres)
-#bpy.types.register(PHYSICS_PT_smoke_cache_highres)
+bpy.types.register(PHYSICS_PT_smoke_cache_highres)
index b7ab07b0f917ebfc17d6c49d8e1c97f67d4f9bed..5ae10d736fde178027bdac58c73c68a4827baa4a 100644 (file)
@@ -232,6 +232,7 @@ void BKE_ptcache_id_from_softbody(PTCacheID *pid, struct Object *ob, struct Soft
 void BKE_ptcache_id_from_particles(PTCacheID *pid, struct Object *ob, struct ParticleSystem *psys);
 void BKE_ptcache_id_from_cloth(PTCacheID *pid, struct Object *ob, struct ClothModifierData *clmd);
 void BKE_ptcache_id_from_smoke(PTCacheID *pid, struct Object *ob, struct SmokeModifierData *smd);
+void BKE_ptcache_id_from_smoke_turbulence(PTCacheID *pid, struct Object *ob, struct SmokeModifierData *smd);
 
 void BKE_ptcache_ids_from_object(struct ListBase *lb, struct Object *ob);
 
index fddcf0fea839af286f676fd17ff062eb72d8a0e7..0f8e9c5edf5d61f4ec23b3e3ecbc5e1ccb3660f8 100644 (file)
 #ifndef BKE_SMOKE_H_
 #define BKE_SMOKE_H_
 
-typedef int (*bresenham_callback) (float *input, int res[3], int *pixel, float *tRay);
+typedef float (*bresenham_callback) (float *result, float *input, int res[3], int *pixel, float *tRay, float correct);
 
 void smokeModifier_do(struct SmokeModifierData *smd, struct Scene *scene, struct Object *ob, struct DerivedMesh *dm, int useRenderParams, int isFinalCalc);
 
 void smokeModifier_free (struct SmokeModifierData *smd);
 void smokeModifier_reset(struct SmokeModifierData *smd);
+void smokeModifier_reset_turbulence(struct SmokeModifierData *smd);
 void smokeModifier_createType(struct SmokeModifierData *smd);
 
+long long smoke_get_mem_req(int xres, int yres, int zres, int amplify);
 
 #endif /* BKE_SMOKE_H_ */
index 912acc4786f3b6d58b68dc35ea8dcf906963bc29..1ab1daf1782618502a5cf9be6f8d606e6b9135bd 100644 (file)
@@ -471,6 +471,19 @@ static int ptcache_totpoint_smoke(void *smoke_v)
                return 0;
 }
 
+/* Smoke functions */
+static int ptcache_totpoint_smoke_turbulence(void *smoke_v)
+{
+       SmokeModifierData *smd= (SmokeModifierData *)smoke_v;
+       SmokeDomainSettings *sds = smd->domain;
+       
+       if(sds->wt) {
+               return sds->res_wt[0]*sds->res_wt[1]*sds->res_wt[2];
+       }
+       else
+               return 0;
+}
+
 // forward decleration
 static int ptcache_file_write(PTCacheFile *pf, void *f, size_t tot, int size);
 
@@ -521,7 +534,7 @@ static int ptcache_compress_write(PTCacheFile *pf, unsigned char *in, unsigned i
 }
 
 static int ptcache_write_smoke(PTCacheFile *pf, void *smoke_v)
-{
+{      
        SmokeModifierData *smd= (SmokeModifierData *)smoke_v;
        SmokeDomainSettings *sds = smd->domain;
        
@@ -535,7 +548,7 @@ static int ptcache_write_smoke(PTCacheFile *pf, void *smoke_v)
 
                smoke_export(sds->fluid, &dt, &dx, &dens, &densold, &heat, &heatold, &vx, &vy, &vz, &vxold, &vyold, &vzold, &obstacles);
 
-               ptcache_compress_write(pf, (unsigned char *)sds->view3d, in_len*4, out, mode);
+               ptcache_compress_write(pf, (unsigned char *)sds->shadow, in_len, out, mode);
                ptcache_compress_write(pf, (unsigned char *)dens, in_len, out, mode);
                ptcache_compress_write(pf, (unsigned char *)densold, in_len, out, mode);        
                ptcache_compress_write(pf, (unsigned char *)heat, in_len, out, mode);
@@ -554,36 +567,36 @@ static int ptcache_write_smoke(PTCacheFile *pf, void *smoke_v)
                
                return 1;
        }
-
        return 0;
 }
 
-/*
 static int ptcache_write_smoke_turbulence(PTCacheFile *pf, void *smoke_v)
 {
        SmokeModifierData *smd= (SmokeModifierData *)smoke_v;
        SmokeDomainSettings *sds = smd->domain;
        
        if(sds->wt) {
-               unsigned int res_big[3];
-               size_t res = sds->res[0]*sds->res[1]*sds->res[2];
+               unsigned int res_big_array[3];
+               unsigned int res_big;
+               unsigned int res = sds->res[0]*sds->res[1]*sds->res[2];
                float *dens, *densold, *tcu, *tcv, *tcw;
                unsigned int in_len = sizeof(float)*(unsigned int)res;
-               unsigned int in_len_big = sizeof(float) * (unsigned int)res_big;
+               unsigned int in_len_big;
                unsigned char *out;
                int mode;
 
-               smoke_turbulence_get_res(sds->wt, res_big);
-               mode = res_big[0]*res_big[1]*res_big[2] >= 1000000 ? 2 : 1;
+               smoke_turbulence_get_res(sds->wt, res_big_array);
+               res_big = res_big_array[0]*res_big_array[1]*res_big_array[2];
+               mode =  res_big >= 1000000 ? 2 : 1;
+               in_len_big = sizeof(float) * (unsigned int)res_big;
 
                smoke_turbulence_export(sds->wt, &dens, &densold, &tcu, &tcv, &tcw);
 
                out = (unsigned char *)MEM_callocN(LZO_OUT_LEN(in_len_big), "pointcache_lzo_buffer");
-
                ptcache_compress_write(pf, (unsigned char *)dens, in_len_big, out, mode);
                ptcache_compress_write(pf, (unsigned char *)densold, in_len_big, out, mode);    
-
                MEM_freeN(out);
+
                out = (unsigned char *)MEM_callocN(LZO_OUT_LEN(in_len), "pointcache_lzo_buffer");
                ptcache_compress_write(pf, (unsigned char *)tcu, in_len, out, mode);
                ptcache_compress_write(pf, (unsigned char *)tcv, in_len, out, mode);
@@ -594,7 +607,6 @@ static int ptcache_write_smoke_turbulence(PTCacheFile *pf, void *smoke_v)
        }
        return 0;
 }
-*/
 
 // forward decleration
 static int ptcache_file_read(PTCacheFile *pf, void *f, size_t tot, int size);
@@ -649,7 +661,7 @@ static void ptcache_read_smoke(PTCacheFile *pf, void *smoke_v)
                
                smoke_export(sds->fluid, &dt, &dx, &dens, &densold, &heat, &heatold, &vx, &vy, &vz, &vxold, &vyold, &vzold, &obstacles);
 
-               ptcache_compress_read(pf, (unsigned char *)sds->view3d, out_len*4);
+               ptcache_compress_read(pf, (unsigned char *)sds->shadow, out_len);
                ptcache_compress_read(pf, (unsigned char*)dens, out_len);
                ptcache_compress_read(pf, (unsigned char*)densold, out_len);
                ptcache_compress_read(pf, (unsigned char*)heat, out_len);
@@ -666,26 +678,32 @@ static void ptcache_read_smoke(PTCacheFile *pf, void *smoke_v)
        }
 }
 
-/*
 static void ptcache_read_smoke_turbulence(PTCacheFile *pf, void *smoke_v)
 {
        SmokeModifierData *smd= (SmokeModifierData *)smoke_v;
        SmokeDomainSettings *sds = smd->domain;
        
        if(sds->fluid) {
-               unsigned int res[3];
+               unsigned int res = sds->res[0]*sds->res[1]*sds->res[2];
+               unsigned int res_big, res_big_array[3];
                float *dens, *densold, *tcu, *tcv, *tcw;
                unsigned int out_len = sizeof(float)*(unsigned int)res;
+               unsigned int out_len_big;
 
-               smoke_turbulence_get_res(sds->wt, res);
+               smoke_turbulence_get_res(sds->wt, res_big_array);
+               res_big = res_big_array[0]*res_big_array[1]*res_big_array[2];
+               out_len_big = sizeof(float) * (unsigned int)res_big;
 
                smoke_turbulence_export(sds->wt, &dens, &densold, &tcu, &tcv, &tcw);
 
-               ptcache_compress_read(pf, (unsigned char*)dens, out_len);
-               
+               ptcache_compress_read(pf, (unsigned char*)dens, out_len_big);
+               ptcache_compress_read(pf, (unsigned char*)densold, out_len_big);
+
+               ptcache_compress_read(pf, (unsigned char*)tcu, out_len);
+               ptcache_compress_read(pf, (unsigned char*)tcv, out_len);
+               ptcache_compress_read(pf, (unsigned char*)tcw, out_len);                
        }
 }
-*/
 
 void BKE_ptcache_id_from_smoke(PTCacheID *pid, struct Object *ob, struct SmokeModifierData *smd)
 {
@@ -716,7 +734,7 @@ void BKE_ptcache_id_from_smoke(PTCacheID *pid, struct Object *ob, struct SmokeMo
        pid->write_header= ptcache_write_basic_header;
        pid->read_header= ptcache_read_basic_header;
 
-       pid->data_types= (1<<BPHYS_DATA_LOCATION); // bogus values tot make pointcache happy
+       pid->data_types= (1<<BPHYS_DATA_LOCATION); // bogus values to make pointcache happy
        pid->info_types= 0;
 }
 
@@ -736,13 +754,13 @@ void BKE_ptcache_id_from_smoke_turbulence(PTCacheID *pid, struct Object *ob, str
        pid->cache_ptr= &sds->point_cache[1];
        pid->ptcaches= &sds->ptcaches[1];
 
-       pid->totpoint= pid->totwrite= ptcache_totpoint_smoke;
+       pid->totpoint= pid->totwrite= ptcache_totpoint_smoke_turbulence;
 
        pid->write_elem= NULL;
        pid->read_elem= NULL;
 
-       pid->read_stream = ptcache_read_smoke;
-       pid->write_stream = ptcache_write_smoke;
+       pid->read_stream = ptcache_read_smoke_turbulence;
+       pid->write_stream = ptcache_write_smoke_turbulence;
        
        pid->interpolate_elem= NULL;
 
@@ -820,6 +838,10 @@ void BKE_ptcache_ids_from_object(ListBase *lb, Object *ob)
                                pid= MEM_callocN(sizeof(PTCacheID), "PTCacheID");
                                BKE_ptcache_id_from_smoke(pid, ob, (SmokeModifierData*)md);
                                BLI_addtail(lb, pid);
+
+                               pid= MEM_callocN(sizeof(PTCacheID), "PTCacheID");
+                               BKE_ptcache_id_from_smoke_turbulence(pid, ob, (SmokeModifierData*)md);
+                               BLI_addtail(lb, pid);
                        }
                }
        }
@@ -1824,6 +1846,8 @@ int BKE_ptcache_id_reset(Scene *scene, PTCacheID *pid, int mode)
                        psys_reset(pid->calldata, PSYS_RESET_DEPSGRAPH);
                else if(pid->type == PTCACHE_TYPE_SMOKE_DOMAIN)
                        smokeModifier_reset(pid->calldata);
+               else if(pid->type == PTCACHE_TYPE_SMOKE_HIGHRES)
+                       smokeModifier_reset(pid->calldata);
        }
        if(clear)
                BKE_ptcache_id_clear(pid, PTCACHE_CLEAR_ALL, 0);
@@ -1878,6 +1902,9 @@ int BKE_ptcache_object_reset(Scene *scene, Object *ob, int mode)
                        {
                                BKE_ptcache_id_from_smoke(&pid, ob, (SmokeModifierData*)md);
                                reset |= BKE_ptcache_id_reset(scene, &pid, mode);
+
+                               BKE_ptcache_id_from_smoke_turbulence(&pid, ob, (SmokeModifierData*)md);
+                               reset |= BKE_ptcache_id_reset(scene, &pid, mode);
                        }
                }
        }
index 223d48012df1f55d43da4bbe56313bcec9fa236e..026ccf4fc7d458c3958bb3a9bf1652e99c162231 100644 (file)
@@ -92,10 +92,10 @@ static void tend ( void )
 {
        QueryPerformanceCounter ( &liCurrentTime );
 }
-//static double tval()
-//{
-//     return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart)* (double)1000.0/(double)liFrequency.QuadPart ));
-//}
+static double tval()
+{
+       return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart)* (double)1000.0/(double)liFrequency.QuadPart ));
+}
 #else
 #include <sys/time.h>
 static struct timeval _tstart, _tend;
@@ -125,7 +125,6 @@ struct SmokeModifierData;
 // forward declerations
 static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int correct);
 void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *tris, int numfaces, int numtris, int **tridivs, float cell_len);
-void smoke_prepare_View(SmokeModifierData *smd, float framenr, float *light, int have_light);
 
 #define TRI_UVOFFSET (1./4.)
 
@@ -167,6 +166,7 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive
                // calc other res with max_res provided
                VECSUB(size, max, min);
 
+               // prevent crash when initializing a plane as domain
                if((size[0] < FLT_EPSILON) || (size[1] < FLT_EPSILON) || (size[2] < FLT_EPSILON))
                        return 0;
 
@@ -215,22 +215,28 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive
                // dt max is 0.1
                smd->domain->fluid = smoke_init(smd->domain->res, smd->domain->p0, 0.1);
                smd->time = scene->r.cfra;
-               smd->domain->firstframe = smd->time;
 
-               if(!smd->domain->wt && (smd->domain->flags & MOD_SMOKE_HIGHRES))
+               if(smd->domain->flags & MOD_SMOKE_HIGHRES)
                {
-                       smd->domain->wt = smoke_turbulence_init(smd->domain->res,  smd->domain->amplify + 1, smd->domain->noise);
-                       smoke_turbulence_initBlenderRNA(smd->domain->wt, &smd->domain->strength);
+                       smd->domain->wt = smoke_turbulence_init(smd->domain->res, smd->domain->amplify + 1, smd->domain->noise);
+                       smd->domain->res_wt[0] = smd->domain->res[0] * (smd->domain->amplify + 1);
+                       smd->domain->res_wt[1] = smd->domain->res[1] * (smd->domain->amplify + 1);                      
+                       smd->domain->res_wt[2] = smd->domain->res[2] * (smd->domain->amplify + 1);                      
+                       smd->domain->dx_wt = smd->domain->dx / (smd->domain->amplify + 1);              
+                       printf("smd->domain->amplify: %d\n",  smd->domain->amplify);
+                       printf("(smd->domain->flags & MOD_SMOKE_HIGHRES)\n");
                }
 
-               if(!smd->domain->view3d)
-               {
-                       // TVox is for transparency
-                       smd->domain->view3d = MEM_callocN(sizeof(float)*smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2]*4, "Smoke_tVox");
-               }
+               if(!smd->domain->shadow)
+                       smd->domain->shadow = MEM_callocN(sizeof(float) * smd->domain->res[0] * smd->domain->res[1] * smd->domain->res[2], "SmokeDomainShadow");
 
                smoke_initBlenderRNA(smd->domain->fluid, &(smd->domain->alpha), &(smd->domain->beta));
 
+               if(smd->domain->wt)     
+               {
+                       smoke_initWaveletBlenderRNA(smd->domain->wt, &(smd->domain->strength));
+                       printf("smoke_initWaveletBlenderRNA\n");
+               }
                return 1;
        }
        else if((smd->type & MOD_SMOKE_TYPE_FLOW) && smd->flow)
@@ -270,13 +276,11 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive
                        SmokeCollSettings *scs = smd->coll;
                        MVert *mvert = dm->getVertArray(dm);
                        MFace *mface = dm->getFaceArray(dm);
-                       size_t i = 0;
-                       int divs = 0;
+                       int i = 0, divs = 0;
                        int *tridivs = NULL;
                        float cell_len = 1.0 / 50.0; // for res = 50
-                       size_t newdivs = 0;
-                       //size_t max_points = 0;
-                       size_t quads = 0, facecounter = 0;
+                       int newdivs = 0;
+                       int quads = 0, facecounter = 0;
 
                        // copy obmat
                        Mat4CpyMat4(scs->mat, ob->obmat);
@@ -314,7 +318,7 @@ int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, Derive
                                int again = 0;
                                do
                                {
-                                       size_t j, k;
+                                       int j, k;
                                        int divs1 = tridivs[3 * facecounter + 0];
                                        int divs2 = tridivs[3 * facecounter + 1];
                                        //int divs3 = tridivs[3 * facecounter + 2];
@@ -521,10 +525,10 @@ void smokeModifier_freeDomain(SmokeModifierData *smd)
 {
        if(smd->domain)
        {
-               // free visualisation buffers
-               if(smd->domain->view3d)
-                       MEM_freeN(smd->domain->view3d);
-               
+               if(smd->domain->shadow)
+                               MEM_freeN(smd->domain->shadow);
+                       smd->domain->shadow = NULL;
+
                if(smd->domain->fluid)
                        smoke_free(smd->domain->fluid);
 
@@ -583,41 +587,46 @@ void smokeModifier_freeCollision(SmokeModifierData *smd)
        }
 }
 
+void smokeModifier_reset_turbulence(struct SmokeModifierData *smd)
+{
+       if(smd && smd->domain && smd->domain->wt)
+       {
+               smoke_turbulence_free(smd->domain->wt);
+               smd->domain->wt = NULL;
+       }
+
+       smd->domain->point_cache[1]->flag &= ~PTCACHE_SIMULATION_VALID;
+       smd->domain->point_cache[1]->flag |= PTCACHE_OUTDATED;
+       smd->domain->point_cache[1]->simframe= 0;
+       smd->domain->point_cache[1]->last_exact= 0;
+}
+
 void smokeModifier_reset(struct SmokeModifierData *smd)
 {
        if(smd)
        {
                if(smd->domain)
                {
-                       if(smd->domain->view3d)
-                               MEM_freeN(smd->domain->view3d);
-                       smd->domain->view3d = NULL;
-
-                       smd->domain->tex = NULL;
+                       if(smd->domain->shadow)
+                               MEM_freeN(smd->domain->shadow);
+                       smd->domain->shadow = NULL;
 
                        if(smd->domain->fluid)
                        {
                                smoke_free(smd->domain->fluid);
                                smd->domain->fluid = NULL;
                        }
-
-                       if(smd->domain->wt)
-                       {
-                               smoke_turbulence_free(smd->domain->wt);
-                               smd->domain->wt = NULL;
-                       }
-               
+                                       
                        smd->domain->point_cache[0]->flag &= ~PTCACHE_SIMULATION_VALID;
                        smd->domain->point_cache[0]->flag |= PTCACHE_OUTDATED;
                        smd->domain->point_cache[0]->simframe= 0;
                        smd->domain->point_cache[0]->last_exact= 0;
 
-                       smd->domain->point_cache[1]->flag &= ~PTCACHE_SIMULATION_VALID;
-                       smd->domain->point_cache[1]->flag |= PTCACHE_OUTDATED;
-                       smd->domain->point_cache[1]->simframe= 0;
-                       smd->domain->point_cache[1]->last_exact= 0;
+                       smokeModifier_reset_turbulence(smd);
 
-                       // printf("reset_domain\n");
+                       smd->time = -1;
+
+                       // printf("reset domain end\n");
                }
                else if(smd->flow)
                {
@@ -685,22 +694,21 @@ void smokeModifier_createType(struct SmokeModifierData *smd)
 
                        /* set some standard values */
                        smd->domain->fluid = NULL;
+                       smd->domain->wt = NULL;                 
                        smd->domain->eff_group = NULL;
                        smd->domain->fluid_group = NULL;
                        smd->domain->coll_group = NULL;
                        smd->domain->maxres = 32;
+                       smd->domain->amplify = 1;                       
+                       smd->domain->omega = 1.0;                       
                        smd->domain->alpha = -0.001;
                        smd->domain->beta = 0.1;
                        smd->domain->flags = MOD_SMOKE_DISSOLVE_LOG;
-                       smd->domain->diss_speed = 5;
-                       smd->domain->strength = 2.0f;
-                       smd->domain->amplify = 1;
+                       smd->domain->strength = 2.0;
                        smd->domain->noise = MOD_SMOKE_NOISEWAVE;
-                       smd->domain->wt = NULL;
-
+                       smd->domain->diss_speed = 5;
                        // init 3dview buffer
-                       smd->domain->view3d = NULL;
-                       smd->domain->tex = NULL;
+                       smd->domain->viewsettings = 0;
                }
                else if(smd->type & MOD_SMOKE_TYPE_FLOW)
                {
@@ -734,9 +742,314 @@ void smokeModifier_createType(struct SmokeModifierData *smd)
        }
 }
 
-// forward declaration
-void smoke_simulate_domain(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm);
+// forward decleration
+void smoke_calc_transparency(float *result, float *input, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb, float correct);
+static float calc_voxel_transp(float *result, float *input, int res[3], int *pixel, float *tRay, float correct);
+static int get_lamp(Scene *scene, float *light)
+{      
+       Base *base_tmp = NULL;  
+       for(base_tmp = scene->base.first; base_tmp; base_tmp= base_tmp->next)   
+       {               
+               if(base_tmp->object->type == OB_LAMP)           
+               {                       
+                       Lamp *la = (Lamp *)base_tmp->object->data;      
+
+                       if(la->type == LA_LOCAL)                        
+                       {                               
+                               VECCOPY(light, base_tmp->object->obmat[3]);                             
+                               return 1;                       
+                       }               
+               }       
+       }       
+       return 0;
+}
+
+static void smoke_calc_domain(Scene *scene, Object *ob, SmokeModifierData *smd)
+{
+       SmokeDomainSettings *sds = smd->domain;
+       GroupObject *go = NULL;                 
+       Base *base = NULL;      
+
+       // do flows and fluids
+       if(1)                   
+       {                               
+               Object *otherobj = NULL;                                
+               ModifierData *md = NULL;
+               if(sds->fluid_group) // we use groups since we have 2 domains
+                       go = sds->fluid_group->gobject.first;                           
+               else                                    
+                       base = scene->base.first;
+               while(base || go)
+               {                                       
+                       otherobj = NULL;
+                       if(sds->fluid_group) 
+                       {
+                               if(go->ob)                                                      
+                                       otherobj = go->ob;                                      
+                       }                                       
+                       else                                            
+                               otherobj = base->object;
+                       if(!otherobj)
+                       {
+                               if(sds->fluid_group)
+                                       go = go->next;
+                               else
+                                       base= base->next;
+
+                               continue;
+                       }
+
+                       md = modifiers_findByType(otherobj, eModifierType_Smoke);
+                       
+                       // check for active smoke modifier
+                       if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
+                       {
+                               SmokeModifierData *smd2 = (SmokeModifierData *)md;
+                               
+                               // check for initialized smoke object
+                               if((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow)                                            
+                               {
+                                       // we got nice flow object
+                                       SmokeFlowSettings *sfs = smd2->flow;
+                                       
+                                       if(sfs->psys && sfs->psys->part && sfs->psys->part->type==PART_EMITTER) // is particle system selected
+                                       {
+                                               ParticleSystem *psys = sfs->psys;
+                                               ParticleSettings *part=psys->part;
+                                               ParticleData *pa = NULL;                                                                
+                                               int p = 0;                                                              
+                                               float *density = smoke_get_density(sds->fluid);                                                         
+                                               float *bigdensity = smoke_turbulence_get_density(sds->wt);                                                              
+                                               float *heat = smoke_get_heat(sds->fluid);                                                               
+                                               float *velocity_x = smoke_get_velocity_x(sds->fluid);                                                           
+                                               float *velocity_y = smoke_get_velocity_y(sds->fluid);                                                           
+                                               float *velocity_z = smoke_get_velocity_z(sds->fluid);                                                           
+                                               unsigned char *obstacle = smoke_get_obstacle(sds->fluid);                                                               
+                                               int bigres[3];  
+                                                                                                               
+                                               // mostly copied from particle code                                                             
+                                               for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++)                                                                
+                                               {                                                                       
+                                                       int cell[3];                                                                    
+                                                       size_t i = 0;                                                                   
+                                                       size_t index = 0;                                                                       
+                                                       int badcell = 0;                                                                                                                                                
+                                                       if(pa->alive == PARS_KILLED) continue;                                                                  
+                                                       else if(pa->alive == PARS_UNBORN && (part->flag & PART_UNBORN)==0) continue;                                                                    
+                                                       else if(pa->alive == PARS_DEAD && (part->flag & PART_DIED)==0) continue;                                                                        
+                                                       else if(pa->flag & (PARS_UNEXIST+PARS_NO_DISP)) continue;                                                                                                                                               
+                                                       // VECCOPY(pos, pa->state.co);                                                                  
+                                                       // Mat4MulVecfl (ob->imat, pos);                                                                                                                                                
+                                                       // 1. get corresponding cell                                                                    
+                                                       get_cell(smd->domain->p0, smd->domain->res, smd->domain->dx, pa->state.co, cell, 0);                                                                                                                                    
+                                                       // check if cell is valid (in the domain boundary)                                                                      
+                                                       for(i = 0; i < 3; i++)                                                                  
+                                                       {                                                                               
+                                                               if((cell[i] > sds->res[i] - 1) || (cell[i] < 0))                                                                                
+                                                               {                                                                                       
+                                                                       badcell = 1;                                                                                    
+                                                                       break;                                                                          
+                                                               }                                                                       
+                                                       }                                                                                                                                                       
+                                                       if(badcell)                                                                             
+                                                               continue;                                                                                                                                               
+                                                       // 2. set cell values (heat, density and velocity)                                                                      
+                                                       index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);                                                                                                                                           
+                                                       if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) && !(obstacle[index] & 2)) // this is inflow                                                                      
+                                                       {                                                                               
+                                                               // heat[index] += sfs->temp * 0.1;                                                                              
+                                                               // density[index] += sfs->density * 0.1;
+                                                               heat[index] = sfs->temp;
+                                                               density[index] = sfs->density;
+
+                                                               /*
+                                                               velocity_x[index] = pa->state.vel[0];
+                                                               velocity_y[index] = pa->state.vel[1];
+                                                               velocity_z[index] = pa->state.vel[2];                                                                           
+                                                               */                                                                              
+                                                               
+                                                               // obstacle[index] |= 2;
+                                                               // we need different handling for the high-res feature
+                                                               if(bigdensity)
+                                                               {
+                                                                       // init all surrounding cells according to amplification, too
+                                                                       int i, j, k;
+
+                                                                       smoke_turbulence_get_res(smd->domain->wt, bigres);
+
+                                                                       for(i = 0; i < smd->domain->amplify + 1; i++)
+                                                                               for(j = 0; j < smd->domain->amplify + 1; j++)
+                                                                                       for(k = 0; k < smd->domain->amplify + 1; k++)                                                                                                   
+                                                                                       {                                                                                                               
+                                                                                               index = smoke_get_index((smd->domain->amplify + 1)* cell[0] + i, bigres[0], (smd->domain->amplify + 1)* cell[1] + j, bigres[1], (smd->domain->amplify + 1)* cell[2] + k);                                                                                                               
+                                                                                               bigdensity[index] = sfs->density;                                                                                                       
+                                                                                       }                                                                               
+                                                               }                                                                       
+                                                       }                                                                       
+                                                       else if(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) // outflow                                                                     
+                                                       {                                                                               
+                                                               heat[index] = 0.f;                                                                              
+                                                               density[index] = 0.f;                                                                           
+                                                               velocity_x[index] = 0.f;                                                                                
+                                                               velocity_y[index] = 0.f;                                                                                
+                                                               velocity_z[index] = 0.f;
+                                                               // we need different handling for the high-res feature
+                                                               if(bigdensity)
+                                                               {
+                                                                       // init all surrounding cells according to amplification, too                                                                                   
+                                                                       int i, j, k;
+                                                                       smoke_turbulence_get_res(smd->domain->wt, bigres);
+
+                                                                       for(i = 0; i < smd->domain->amplify + 1; i++)
+                                                                               for(j = 0; j < smd->domain->amplify + 1; j++)
+                                                                                       for(k = 0; k < smd->domain->amplify + 1; k++)
+                                                                                       {                                                                                                               
+                                                                                               index = smoke_get_index((smd->domain->amplify + 1)* cell[0] + i, bigres[0], (smd->domain->amplify + 1)* cell[1] + j, bigres[1], (smd->domain->amplify + 1)* cell[2] + k);                                                                                                               
+                                                                                               bigdensity[index] = 0.f;                                                                                                        
+                                                                                       }                                                                               
+                                                               }                                                                       
+                                                       }       // particles loop                                                       
+                                       }                                                       
+                               }                                                       
+                               else                                                    
+                               {                                                               
+                                       /*                                                              
+                                       for()                                                           
+                                       {                                                                       
+                                               // no psys                                                                      
+                                               BVHTreeNearest nearest;
+                                               nearest.index = -1;
+                                               nearest.dist = FLT_MAX;
+
+                                               BLI_bvhtree_find_nearest(sfs->bvh->tree, pco, &nearest, sfs->bvh->nearest_callback, sfs->bvh);
+                                       }*/                                                     
+                               }                                               
+                       }                                               
+               }
+                       if(sds->fluid_group)
+                               go = go->next;
+                       else
+                               base= base->next;
+               }
+       }
+
+       // do effectors
+       /*
+       if(sds->eff_group)
+       {
+               for(go = sds->eff_group->gobject.first; go; go = go->next) 
+               {
+                       if(go->ob)
+                       {
+                               if(ob->pd)
+                               {
+                                       
+                               }                                       
+                       }
+               }
+       }
+       */
+
+       // do collisions        
+       if(1)
+       {
+               Object *otherobj = NULL;
+               ModifierData *md = NULL;
+
+               if(sds->coll_group) // we use groups since we have 2 domains
+                       go = sds->coll_group->gobject.first;
+               else
+                       base = scene->base.first;
+
+               while(base || go)
+               {
+                       otherobj = NULL;
+                       if(sds->coll_group) 
+                       {                                               
+                               if(go->ob)                                                      
+                                       otherobj = go->ob;                                      
+                       }                                       
+                       else                                            
+                               otherobj = base->object;                                        
+                       if(!otherobj)                                   
+                       {                                               
+                               if(sds->coll_group)                                                     
+                                       go = go->next;                                          
+                               else                                                    
+                                       base= base->next;                                               
+                               continue;                                       
+                       }                       
+                       md = modifiers_findByType(otherobj, eModifierType_Smoke);
+                       
+                       // check for active smoke modifier
+                       if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))                                    
+                       {
+                               SmokeModifierData *smd2 = (SmokeModifierData *)md;
+
+                               if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll)
+                               {
+                                       // we got nice collision object
+                                       SmokeCollSettings *scs = smd2->coll;
+                                       size_t i, j;
+                                       unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid);
+
+                                       for(i = 0; i < scs->numpoints; i++)
+                                       {
+                                               int badcell = 0;
+                                               size_t index = 0;
+                                               int cell[3];
+
+                                               // 1. get corresponding cell
+                                               get_cell(smd->domain->p0, smd->domain->res, smd->domain->dx, &scs->points[3 * i], cell, 0);
+                                       
+                                               // check if cell is valid (in the domain boundary)
+                                               for(j = 0; j < 3; j++)
+                                                       if((cell[j] > sds->res[j] - 1) || (cell[j] < 0))
+                                                       {
+                                                               badcell = 1;
+                                                               break;
+                                                       }
+                                                                                                                               
+                                                       if(badcell)                                                                     
+                                                               continue;
+                                               // 2. set cell values (heat, density and velocity)
+                                               index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
+                                                                                                               
+                                               // printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]);                                                                
+                                               // printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index);                                                                                                                                   
+                                               obstacles[index] = 1;
+                                               // for moving gobstacles                                                                
+                                               /*
+                                               const LbmFloat maxVelVal = 0.1666;
+                                               const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
+
+                                               LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); 
+                                               {                                                               
+                                               const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5;                                                                
+                                               USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz);                                                                 
+                                               if(usqr>maxusqr) {                                                                      
+                                               // cutoff at maxVelVal                                                                  
+                                               for(int jj=0; jj<3; jj++) {                                                                             
+                                               if(objvel[jj]>0.) objvel[jj] =  maxVelVal;                                                                              
+                                               if(objvel[jj]<0.) objvel[jj] = -maxVelVal;                                                                      
+                                               }                                                               
+                                               } 
+                                               }                                                               
+                                               const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) );                                                          
+                                               const LbmVec oldov=objvel; // debug                                                             
+                                               objvel = vec2L((*pNormals)[n]) *dp;                                                             
+                                               */
+                                       }
+                               }
+                       }
 
+                       if(sds->coll_group)
+                               go = go->next;
+                       else
+                               base= base->next;
+               }
+       }
+}
 void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc)
 {      
        if((smd->type & MOD_SMOKE_TYPE_FLOW))
@@ -788,16 +1101,15 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
        }
        else if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
        {
-               PointCache *cache;
+               PointCache *cache = NULL;
                PTCacheID pid;
+               PointCache *cache_wt = NULL;
+               PTCacheID pid_wt;
                float timescale;
-               int cache_result = 0;
-               int startframe, endframe, framenr;
+               int cache_result = 0, cache_result_wt = 0;
+               int startframe, endframe, framenr, badloading = 0;
                SmokeDomainSettings *sds = smd->domain;
-               float light[3] = {0.0,0.0,0.0};
-               int have_lamp = 0;
-
-               // printf("smd->type & MOD_SMOKE_TYPE_DOMAIN\n");
+               float light[3]; 
 
                framenr = scene->r.cfra;
 
@@ -806,51 +1118,26 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
                BKE_ptcache_id_from_smoke(&pid, ob, smd);
                BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, &timescale);
 
+               cache_wt = sds->point_cache[1];
+               BKE_ptcache_id_from_smoke_turbulence(&pid_wt, ob, smd);
+
                /* handle continuous simulation with the play button */
                if(BKE_ptcache_get_continue_physics()) 
                {
-                       cache->flag &= ~PTCACHE_SIMULATION_VALID;
-                       cache->simframe= 0;
-                       cache->last_exact= 0;
-
-                       if(!smokeModifier_init(smd, ob, scene, dm))
-                               return;
-
-                       if(!smd->domain->fluid)
-                               return;
-
-                       smoke_simulate_domain(smd, scene, ob, dm);
-
-                       {
-                               Base *base_tmp = NULL;
-
-                               for(base_tmp = scene->base.first; base_tmp; base_tmp= base_tmp->next) 
-                               {
-                                       if(base_tmp->object->type == OB_LAMP) 
-                                       {
-                                               Lamp *la = (Lamp *)base_tmp->object->data;
-                                               
-                                               if(la->type == LA_LOCAL)
-                                               {
-                                                       VECCOPY(light, base_tmp->object->obmat[3]);
-                                                       have_lamp = 1;
-                                                       break;
-                                               }
-                                       }
-                               }
-                       }
-
-                       smoke_prepare_View(smd, (float)framenr, light, have_lamp);
-
+                       // TODO
                        return;
                }
-               
+
                if(framenr < startframe)
                {
                        cache->flag &= ~PTCACHE_SIMULATION_VALID;
                        cache->simframe= 0;
                        cache->last_exact= 0;
 
+                       cache_wt->flag &= ~PTCACHE_SIMULATION_VALID;
+                       cache_wt->simframe= 0;
+                       cache_wt->last_exact= 0;
+
                        // we got back in time, reset smoke in this case (TODO: use cache later)
                        // smd->time = scene->r.cfra;
                        // smokeModifier_reset(smd);
@@ -860,14 +1147,25 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
                else if(framenr > endframe) 
                {
                        framenr = endframe;
+
+                       // we load last valid frame here 
+                       // and don't update the smd->time variable later
+                       badloading = 1;
                }
 
                if(!(cache->flag & PTCACHE_SIMULATION_VALID))
                {
-                       // printf("reseting\n");
                        BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED);
                }
-       
+               if(!(cache_wt->flag & PTCACHE_SIMULATION_VALID))
+               {
+                       BKE_ptcache_id_reset(scene, &pid_wt, PTCACHE_RESET_OUTDATED);
+               }
+
+               if(smd->time == -1 && framenr!= startframe)
+                       return;
+               
+
                if(!smokeModifier_init(smd, ob, scene, dm))
                        return;
 
@@ -875,7 +1173,7 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
                                return;
 
                /* try to read from cache */
-               cache_result = BKE_ptcache_read_cache(&pid, (float)framenr, scene->r.frs_sec);
+               cache_result =  BKE_ptcache_read_cache(&pid, (float)framenr, scene->r.frs_sec);
                // printf("cache_result: %d\n", cache_result);
 
                if(cache_result == PTCACHE_READ_EXACT) 
@@ -886,417 +1184,148 @@ void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedM
                        cache->simframe= framenr;
                        sds->v3dnum = framenr;
 
+                       if(!badloading)
+                               smd->time = scene->r.cfra;
+
+                       // check for wt cache
+                       if(sds->wt)
+                       {
+                               cache_result_wt = BKE_ptcache_read_cache(&pid_wt, (float)framenr, scene->r.frs_sec);
+                               // printf("cache_result_wt: %d\n", cache_result_wt);
+
+                               // error handling
+                               if(cache_result_wt == PTCACHE_READ_EXACT) 
+                               {
+                                       cache_wt->flag |= PTCACHE_SIMULATION_VALID;
+                                       cache_wt->simframe= framenr;
+                               }
+                               else if(cache_result_wt==PTCACHE_READ_OLD) 
+                               {
+                                       BKE_ptcache_id_reset(scene, &pid_wt, PTCACHE_RESET_FREE);
+                                       cache_wt->flag |= PTCACHE_SIMULATION_VALID;
+                               }
+                               else if(ob->id.lib || (cache_wt->flag & PTCACHE_BAKED)) 
+                               {
+                                       // if baked and nothing in cache, do nothing 
+                                       cache_wt->flag &= ~PTCACHE_SIMULATION_VALID;
+                                       cache_wt->simframe= 0;
+                                       cache_wt->last_exact= 0;
+                               }
+                       }
+
                        // printf("PTCACHE_READ_EXACT\n");
                        return;
                }
                else if(cache_result==PTCACHE_READ_OLD) 
                {
                        BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_FREE);
-
-                       // printf("PTCACHE_READ_OLD\n");
-
                        cache->flag |= PTCACHE_SIMULATION_VALID;
+
+                       BKE_ptcache_id_reset(scene, &pid_wt, PTCACHE_RESET_FREE);
+                       cache_wt->flag |= PTCACHE_SIMULATION_VALID;
                }
                else if(ob->id.lib || (cache->flag & PTCACHE_BAKED)) 
                {
-                       /* if baked and nothing in cache, do nothing */
+                       // if baked and nothing in cache, do nothing 
                        cache->flag &= ~PTCACHE_SIMULATION_VALID;
                        cache->simframe= 0;
                        cache->last_exact= 0;
 
+                       cache_wt->flag &= ~PTCACHE_SIMULATION_VALID;
+                       cache_wt->simframe= 0;
+                       cache_wt->last_exact= 0;
+               
                        // printf("PTCACHE_BAKED\n");
                        return;
                }
-               else if((cache_result==0) && (startframe!=framenr) && !(cache->flag & PTCACHE_SIMULATION_VALID))
+               /*
+               else if((cache_result==0) && ((startframe!=framenr) && !(cache->flag & PTCACHE_SIMULATION_VALID || (framenr == smd->time))))
                {
                        cache->flag &= ~PTCACHE_SIMULATION_VALID;
                        cache->simframe= 0;
                        cache->last_exact= 0;
 
                        return;
-               }
-               
+               }*/
+
+               // printf("framenr: %d, time: %f\n", framenr, smd->time);
+
                /* do simulation */
 
                // low res
                cache->flag |= PTCACHE_SIMULATION_VALID;
                cache->simframe= framenr;
 
-               smoke_simulate_domain(smd, scene, ob, dm);
-
                if(sds->wt)
-                       smoke_turbulence_step(sds->wt, sds->fluid);
-
                {
-                       Base *base_tmp = NULL;
-
-                       for(base_tmp = scene->base.first; base_tmp; base_tmp= base_tmp->next) 
-                       {
-                               if(base_tmp->object->type == OB_LAMP) 
-                               {
-                                       Lamp *la = (Lamp *)base_tmp->object->data;
-                                       
-                                       if(la->type == LA_LOCAL)
-                                       {
-                                               VECCOPY(light, base_tmp->object->obmat[3]);
-                                               have_lamp = 1;
-                                               break;
-                                       }
-                               }
-                       }
+                       cache_wt->flag |= PTCACHE_SIMULATION_VALID;
+                       cache_wt->simframe= framenr;
                }
-
-               smoke_prepare_View(smd, (float)framenr, light, have_lamp);
-
-               BKE_ptcache_write_cache(&pid, framenr);
-
-               // printf("Writing cache_low, %d\n", framenr);
-       
-
-               tend();
-               // printf ( "Frame: %d, Time: %f\n", (int)smd->time, ( float ) tval() );
-       }
-}
-
-void smoke_simulate_domain(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm)
-{
-       GroupObject *go = NULL;
-       Base *base = NULL;
-       SmokeDomainSettings *sds = smd->domain;
-       
-       tstart();
-       
-       if(sds->flags & MOD_SMOKE_DISSOLVE)
-               smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
-
-       // do flows and fluids
-       if(1)
-       {
-               Object *otherobj = NULL;
-               ModifierData *md = NULL;
-
-               if(sds->fluid_group) // we use groups since we have 2 domains
-                       go = sds->fluid_group->gobject.first;
-               else
-                       base = scene->base.first;
-
-               while(base || go)
-               {
-                       otherobj = NULL;
-
-                       if(sds->fluid_group) 
-                       {
-                               if(go->ob)
-                                       otherobj = go->ob;
-                       }
-                       else
-                               otherobj = base->object;
-
-                       if(!otherobj)
-                       {
-                               if(sds->fluid_group)
-                                       go = go->next;
-                               else
-                                       base= base->next;
-
-                               continue;
-                       }
-
-                       md = modifiers_findByType(otherobj, eModifierType_Smoke);
-                       
-                       // check for active smoke modifier
-                       if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
-                       {
-                               SmokeModifierData *smd2 = (SmokeModifierData *)md;
                                
-                               // check for initialized smoke object
-                               if((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow)
-                               {
-                                       // we got nice flow object
-                                       SmokeFlowSettings *sfs = smd2->flow;
-                                       
-                                       if(sfs->psys && sfs->psys->part && sfs->psys->part->type==PART_EMITTER) // is particle system selected
-                                       {
-                                               ParticleSystem *psys = sfs->psys;
-                                               ParticleSettings *part=psys->part;
-                                               ParticleData *pa = NULL;
-                                               int p = 0;
-                                               float *density = smoke_get_density(sds->fluid);
-                                               // float *bigdensity = smoke_turbulence_get_density(sds->wt);
-                                               float *heat = smoke_get_heat(sds->fluid);
-                                               float *velocity_x = smoke_get_velocity_x(sds->fluid);
-                                               float *velocity_y = smoke_get_velocity_y(sds->fluid);
-                                               float *velocity_z = smoke_get_velocity_z(sds->fluid);
-                                               unsigned char *obstacle = smoke_get_obstacle(sds->fluid);
-
-                                               // debug printf("found flow psys\n");
-                                               
-                                               // mostly copied from particle code
-                                               for(p=0, pa=psys->particles; p<psys->totpart; p++, pa++)
-                                               {
-                                                       int cell[3];
-                                                       size_t i = 0;
-                                                       size_t index = 0;
-                                                       int badcell = 0;
-                                                       
-                                                       if(pa->alive == PARS_KILLED) continue;
-                                                       else if(pa->alive == PARS_UNBORN && (part->flag & PART_UNBORN)==0) continue;
-                                                       else if(pa->alive == PARS_DEAD && (part->flag & PART_DIED)==0) continue;
-                                                       else if(pa->flag & (PARS_UNEXIST+PARS_NO_DISP)) continue;
-                                                       
-                                                       // VECCOPY(pos, pa->state.co);
-                                                       // Mat4MulVecfl (ob->imat, pos);
-                                                       
-                                                       // 1. get corresponding cell
-                                                       get_cell(sds->p0, sds->res, sds->dx, pa->state.co, cell, 0);
-                                               
-                                                       // check if cell is valid (in the domain boundary)
-                                                       for(i = 0; i < 3; i++)
-                                                       {
-                                                               if((cell[i] > sds->res[i] - 1) || (cell[i] < 0))
-                                                               {
-                                                                       badcell = 1;
-                                                                       break;
-                                                               }
-                                                       }
-                                                               
-                                                       if(badcell)
-                                                               continue;
-                                                       
-                                                       // 2. set cell values (heat, density and velocity)
-                                                       index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
-                                                       
-                                                       if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) && !(obstacle[index] & 2)) // this is inflow
-                                                       {
-                                                               // heat[index] += sfs->temp * 0.1;
-                                                               // density[index] += sfs->density * 0.1;
-
-                                                               heat[index] = sfs->temp;
-                                                               density[index] = sfs->density;
-
-                                                               /*
-                                                               velocity_x[index] = pa->state.vel[0];
-                                                               velocity_y[index] = pa->state.vel[1];
-                                                               velocity_z[index] = pa->state.vel[2];
-                                                               */
-                                                               obstacle[index] |= 2;
-                                                       }
-                                                       else if(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) // outflow
-                                                       {
-                                                               heat[index] = 0.f;
-                                                               density[index] = 0.f;
-                                                               velocity_x[index] = 0.f;
-                                                               velocity_y[index] = 0.f;
-                                                               velocity_z[index] = 0.f;
-                                                       }
-                                               }
-                                       }
-                                       else
-                                       {
-                                               /*
-                                               for()
-                                               {
-                                                       // no psys
-                                                       BVHTreeNearest nearest;
-
-                                                       nearest.index = -1;
-                                                       nearest.dist = FLT_MAX;
-
-                                                       BLI_bvhtree_find_nearest(sfs->bvh->tree, pco, &nearest, sfs->bvh->nearest_callback, sfs->bvh);
-                                               }*/
-                                       }
-                               }       
-                       }
+               tstart();       
 
-                       if(sds->fluid_group)
-                               go = go->next;
-                       else
-                               base= base->next;
-               }
-       }
-
-       // do effectors
-       /*
-       if(sds->eff_group)
-       {
-               for(go = sds->eff_group->gobject.first; go; go = go->next) 
-               {
-                       if(go->ob)
-                       {
-                               if(ob->pd)
-                               {
-                                       
-                               }
-                       }
+               if(sds->flags & MOD_SMOKE_DISSOLVE)                     
+               {                               
+                       smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);                                                               
                }
-       }
-       */
-
-       // do collisions        
-       if(1)
-       {
-               Object *otherobj = NULL;
-               ModifierData *md = NULL;
 
-               if(sds->coll_group) // we use groups since we have 2 domains
-                       go = sds->coll_group->gobject.first;
-               else
-                       base = scene->base.first;
+               smoke_calc_domain(scene, ob, smd);
+               
+               // set new time
+               smd->time = scene->r.cfra;
 
-               while(base || go)
+               // frame 1 is start, don't simulate anything
+               if(smd->time == 1)
                {
-                       otherobj = NULL;
-
-                       if(sds->coll_group) 
-                       {
-                               if(go->ob)
-                                       otherobj = go->ob;
-                       }
-                       else
-                               otherobj = base->object;
+                       // set new time
+                       smd->time = scene->r.cfra;
 
-                       if(!otherobj)
-                       {
-                               if(sds->coll_group)
-                                       go = go->next;
-                               else
-                                       base= base->next;
+                       BKE_ptcache_write_cache(&pid, framenr);
+                       if(sds->wt)
+                               BKE_ptcache_write_cache(&pid_wt, framenr);
 
-                               continue;
-                       }
-       
-                       md = modifiers_findByType(otherobj, eModifierType_Smoke);
-                       
-                       // check for active smoke modifier
-                       if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
-                       {
-                               SmokeModifierData *smd2 = (SmokeModifierData *)md;
+                       if(get_lamp(scene, light))
+                               smoke_calc_transparency(sds->shadow, smoke_get_density(sds->fluid), sds->p0, sds->p1, sds->res, sds->dx, light, calc_voxel_transp, -7.0*sds->dx);
 
-                               if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll)
-                               {
-                                       // we got nice collision object
-                                       SmokeCollSettings *scs = smd2->coll;
-                                       size_t i, j;
-                                       unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid);
+                       return;
+               }
 
-                                       for(i = 0; i < scs->numpoints; i++)
-                                       {
-                                               int badcell = 0;
-                                               size_t index = 0;
-                                               int cell[3];
+               // simulate the actual smoke (c++ code in intern/smoke)
+               smoke_step(sds->fluid, smd->time);
+               BKE_ptcache_write_cache(&pid, framenr);
 
-                                               // 1. get corresponding cell
-                                               get_cell(sds->p0, sds->res, sds->dx, &scs->points[3 * i], cell, 0);
-                                       
-                                               // check if cell is valid (in the domain boundary)
-                                               for(j = 0; j < 3; j++)
-                                                       if((cell[j] > sds->res[j] - 1) || (cell[j] < 0))
-                                                       {
-                                                               badcell = 1;
-                                                               break;
-                                                       }
-                                                               
-                                               if(badcell)
-                                                       continue;
+               if(sds->wt)
+               {
 
-                                               // 2. set cell values (heat, density and velocity)
-                                               index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
-                                               
-                                               // printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]);
-                                               // printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index);
-                                                       
-                                               obstacles[index] = 1;
+                       if(sds->flags & MOD_SMOKE_DISSOLVE)                                                     
+                               smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
 
-                                               // for moving gobstacles
-                                               /*
-                                               const LbmFloat maxVelVal = 0.1666;
-                                               const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
+                       smoke_turbulence_step(sds->wt, sds->fluid);
+                       BKE_ptcache_write_cache(&pid_wt, framenr);
+               }
 
-                                               LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { 
-                                               const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; 
-                                               USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); 
-                                               if(usqr>maxusqr) { 
-                                                       // cutoff at maxVelVal 
-                                                       for(int jj=0; jj<3; jj++) { 
-                                                               if(objvel[jj]>0.) objvel[jj] =  maxVelVal;  
-                                                               if(objvel[jj]<0.) objvel[jj] = -maxVelVal; 
-                                                       } 
-                                               } } 
-
-                                               const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); 
-                                               const LbmVec oldov=objvel; // debug
-                                               objvel = vec2L((*pNormals)[n]) *dp;
-                                               */
-                                       }
-                               }
-                       }
+               if(get_lamp(scene, light))
+                       smoke_calc_transparency(sds->shadow, smoke_get_density(sds->fluid), sds->p0, sds->p1, sds->res, sds->dx, light, calc_voxel_transp, -7.0*sds->dx);
 
-                       if(sds->coll_group)
-                               go = go->next;
-                       else
-                               base= base->next;
-               }
+               tend();
+               // printf ( "Frame: %d, Time: %f\n", (int)smd->time, ( float ) tval() );
        }
-       
-       // set new time
-       smd->time = scene->r.cfra;
-
-       // simulate the actual smoke (c++ code in intern/smoke)
-       smoke_step(sds->fluid, smd->time);
 }
 
-static int calc_voxel_transp(float *input, int res[3], int *pixel, float *tRay)
+static float calc_voxel_transp(float *result, float *input, int res[3], int *pixel, float *tRay, float correct)
 {
        const size_t index = smoke_get_index(pixel[0], res[0], pixel[1], res[1], pixel[2]);
 
        // T_ray *= T_vox
-       *tRay *= input[index*4];
-
-       return *tRay;
-}
-
-// forward decleration
-void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb);
-
-// update necessary information for 3dview
-void smoke_prepare_View(SmokeModifierData *smd, float framenr, float *light, int have_light)
-{
-       float *density = NULL;
-       int x, y, z;
-       size_t cells, i;
-       SmokeDomainSettings *sds = smd->domain;
-
-       // update 3dview
-       density = smoke_get_density(smd->domain->fluid);
-       for(x = 0; x < smd->domain->res[0]; x++)
-               for(y = 0; y < smd->domain->res[1]; y++)
-                       for(z = 0; z < smd->domain->res[2]; z++)
-                       {
-                               size_t index;
-
-                               index = smoke_get_index(x, smd->domain->res[0], y, smd->domain->res[1], z);
-                               // Transparency computation
-                               // formula taken from "Visual Simulation of Smoke" / Fedkiw et al. pg. 4
-                               // T_vox = exp(-C_ext * h)
-                               // C_ext/sigma_t = density * C_ext
-                               smd->domain->view3d[index * 4] = smd->domain->view3d[index * 4 + 1] = 
-                               smd->domain->view3d[index * 4 + 2] = exp(-density[index] * 7.0 * smd->domain->dx);
-                               smd->domain->view3d[index * 4 + 3] = 1.0 - smd->domain->view3d[index * 4];
-                               
-                       }
-
-       if(have_light)
+       *tRay *= exp(input[index]*correct);
+       
+       if(result[index] < 0.0f)        
        {
-               smoke_calc_transparency(sds->view3d, sds->p0, sds->p1, sds->res, sds->dx, light, calc_voxel_transp);
+#pragma omp critical           
+               result[index] = *tRay;  
+       }       
 
-               cells = smd->domain->res[0]*smd->domain->res[1]*smd->domain->res[2];
-               for(i = 0; i < cells; i++)
-               {
-                       smd->domain->view3d[i * 4] = smd->domain->view3d[i * 4 + 1] = 
-                       smd->domain->view3d[i * 4 + 2] = smd->domain->view3d[i * 4 + 1] * smd->domain->view3d[i * 4 + 0];
-               }
-       }
-       smd->domain->v3dnum = framenr;
+       return *tRay;
 }
 
 long long smoke_get_mem_req(int xres, int yres, int zres, int amplify)
@@ -1317,7 +1346,7 @@ long long smoke_get_mem_req(int xres, int yres, int zres, int amplify)
          return totalMB;
 }
 
-static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, float *tRay, bresenham_callback cb, float *input, int res[3])
+static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, float *tRay, bresenham_callback cb, float *result, float *input, int res[3], float correct)
 {
     int dx, dy, dz, i, l, m, n, x_inc, y_inc, z_inc, err_1, err_2, dx2, dy2, dz2;
     int pixel[3];
@@ -1344,8 +1373,8 @@ static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, f
         err_1 = dy2 - l;
         err_2 = dz2 - l;
         for (i = 0; i < l; i++) {
-               if(cb(input, res, pixel, tRay) < 0.0)
-                       return;
+               if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON)
+                       break;
             if (err_1 > 0) {
                 pixel[1] += y_inc;
                 err_1 -= dx2;
@@ -1362,8 +1391,8 @@ static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, f
         err_1 = dx2 - m;
         err_2 = dz2 - m;
         for (i = 0; i < m; i++) {
-               if(cb(input, res, pixel, tRay) < 0.0f)
-                       return;
+               if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON)
+                       break;
             if (err_1 > 0) {
                 pixel[0] += x_inc;
                 err_1 -= dy2;
@@ -1380,8 +1409,8 @@ static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, f
         err_1 = dy2 - n;
         err_2 = dx2 - n;
         for (i = 0; i < n; i++) {
-               if(cb(input, res, pixel, tRay) < 0.0f)
-                       return;
+               if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON)
+                       break;
             if (err_1 > 0) {
                 pixel[1] += y_inc;
                 err_1 -= dz2;
@@ -1395,7 +1424,7 @@ static void bresenham_linie_3D(int x1, int y1, int z1, int x2, int y2, int z2, f
             pixel[2] += z_inc;
         }
     }
-    cb(input, res, pixel, tRay);
+    cb(result, input, res, pixel, tRay, correct);
 }
 
 static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int correct)
@@ -1419,12 +1448,12 @@ static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int
        }
 }
 
-void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb)
+void smoke_calc_transparency(float *result, float *input, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb, float correct)
 {
        int x, y, z;
        float bv[6];
 
-       // x
+       memset(result, -1, sizeof(float)*res[0]*res[1]*res[2]); // x
        bv[0] = p0[0];
        bv[1] = p1[0];
        // y
@@ -1434,7 +1463,7 @@ void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], fl
        bv[4] = p0[2];
        bv[5] = p1[2];
 
-// #pragma omp parallel for schedule(static) private(y, z)
+#pragma omp parallel for schedule(static) private(y, z)
        for(x = 0; x < res[0]; x++)
                for(y = 0; y < res[1]; y++)
                        for(z = 0; z < res[2]; z++)
@@ -1447,6 +1476,8 @@ void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], fl
 
                                index = smoke_get_index(x, res[0], y, res[1], z);
 
+                               if(result[index] >= 0.0f)                                       
+                                       continue;                                                               
                                voxelCenter[0] = p0[0] + dx *  x + dx * 0.5;
                                voxelCenter[1] = p0[1] + dx *  y + dx * 0.5;
                                voxelCenter[2] = p0[2] + dx *  z + dx * 0.5;
@@ -1463,10 +1494,11 @@ void smoke_calc_transparency(float *result, float *p0, float *p1, int res[3], fl
                                        get_cell(p0, res, dx, light, cell, 1);
                                }
 
-                               bresenham_linie_3D(cell[0], cell[1], cell[2], x, y, z, &tRay, cb, result, res);
+                               bresenham_linie_3D(cell[0], cell[1], cell[2], x, y, z, &tRay, cb, result, input, res, correct);
 
                                // convention -> from a RGBA float array, use G value for tRay
-                               result[index*4 + 1] = tRay;
+// #pragma omp critical
+                               result[index] = tRay;                   
                        }
 }
 
index 618e587d6745f97c116ab140711b7dd66c94377c..b6fc9f7af12e5f57d655f3f520766a134e0b8e13 100644 (file)
@@ -3685,8 +3685,8 @@ static void direct_link_modifiers(FileData *fd, ListBase *lb)
                                smd->domain->smd = smd;
 
                                smd->domain->fluid = NULL;
-                               smd->domain->view3d = NULL;
-                               smd->domain->tex = NULL;
+                               // smd->domain->view3d = NULL;
+                               // smd->domain->tex = NULL;
 
                                direct_link_pointcache_list(fd, &(smd->domain->ptcaches[0]), &(smd->domain->point_cache[0]));
                                direct_link_pointcache_list(fd, &(smd->domain->ptcaches[1]), &(smd->domain->point_cache[1]));
index 20cab3b8aebf66360192b4e1f3ceba3bad65923a..7ed029f3eaf0506a684b443268ccdf4f6628a98c 100644 (file)
@@ -5310,13 +5310,60 @@ void draw_object(Scene *scene, ARegion *ar, View3D *v3d, Base *base, int flag)
        }
 
        /* draw code for smoke */
+       if((md = modifiers_findByType(ob, eModifierType_Smoke)))
        {
-               md = modifiers_findByType(ob, eModifierType_Smoke);
-               if (md) {
-                       SmokeModifierData *smd = (SmokeModifierData *)md;
-                       if(smd->type & MOD_SMOKE_TYPE_DOMAIN && smd->domain && smd->domain->fluid) {
-                               GPU_create_smoke(smd);
-                               draw_volume(scene, ar, v3d, base, smd->domain->tex, smd->domain->res);
+               SmokeModifierData *smd = (SmokeModifierData *)md;
+
+               // draw collision objects
+               if((smd->type & MOD_SMOKE_TYPE_COLL) && smd->coll)
+               {
+                       /*SmokeCollSettings *scs = smd->coll;
+                       if(scs->points)
+                       {
+                               size_t i;
+
+                               wmLoadMatrix(rv3d->viewmat);
+
+                               if(col || (ob->flag & SELECT)) cpack(0xFFFFFF); 
+                               glDepthMask(GL_FALSE);
+                               glEnable(GL_BLEND);
+                               
+
+                               // glPointSize(3.0);
+                               bglBegin(GL_POINTS);
+
+                               for(i = 0; i < scs->numpoints; i++)
+                               {
+                                       bglVertex3fv(&scs->points[3*i]);
+                               }
+
+                               bglEnd();
+                               glPointSize(1.0);
+
+                               wmMultMatrix(ob->obmat);
+                               glDisable(GL_BLEND);
+                               glDepthMask(GL_TRUE);
+                               if(col) cpack(col);
+                               
+                       }
+                       */
+               }
+
+               // only draw domains
+               if(smd->domain && smd->domain->fluid)
+               {
+                       if(!smd->domain->wt || !(smd->domain->viewsettings & MOD_SMOKE_VIEW_SHOWBIG))
+                       {
+                               smd->domain->tex = NULL;
+                               GPU_create_smoke(smd, 0);
+                               draw_volume(scene, ar, v3d, base, smd->domain->tex, smd->domain->res, smd->domain->dx, smd->domain->tex_shadow);
+                               GPU_free_smoke(smd);
+                       }
+                       else if(smd->domain->wt || (smd->domain->viewsettings & MOD_SMOKE_VIEW_SHOWBIG))
+                       {
+                               smd->domain->tex = NULL;
+                               GPU_create_smoke(smd, 1);
+                               draw_volume(scene, ar, v3d, base, smd->domain->tex, smd->domain->res_wt, smd->domain->dx_wt, smd->domain->tex_shadow);
                                GPU_free_smoke(smd);
                        }
                }
index 342acfe92b3f59441cc94f38707de923cc9dcbb7..c8eda10566c60e82414618e49dd3829081975661 100644 (file)
@@ -191,7 +191,7 @@ static int larger_pow2(int n)
        return n*2;
 }
 
-void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture *tex, int res[3])
+void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture *tex, int res[3], float dx, GPUTexture *tex_shadow)
 {
        Object *ob = base->object;
        RegionView3D *rv3d= ar->regiondata;
@@ -204,6 +204,27 @@ void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture
        float cor[3] = {1.,1.,1.};
        int gl_depth = 0, gl_blend = 0;
 
+       /* Fragment program to calculate the 3dview of smoke */
+       /* using 2 textures, density and shadow */
+       const char *text = "!!ARBfp1.0\n"
+                                       "PARAM dx = program.local[0];\n"
+                                       "PARAM darkness = program.local[1];\n"
+                                       "PARAM f = {1.442695041, 1.442695041, 1.442695041, 0.01};\n"
+                                       "TEMP temp, shadow, value;\n"
+                                       "TEX temp, fragment.texcoord[0], texture[0], 3D;\n"
+                                       "TEX shadow, fragment.texcoord[0], texture[1], 3D;\n"
+                                       "MUL value, temp, darkness;\n"
+                                       "MUL value, value, dx;\n"
+                                       "MUL value, value, f;\n"
+                                       "EX2 temp, -value.r;\n"
+                                       "SUB temp.a, 1.0, temp.r;\n"
+                                       "MUL temp.r, temp.r, shadow.r;\n"
+                                       "MUL temp.g, temp.g, shadow.r;\n"
+                                       "MUL temp.b, temp.b, shadow.r;\n"
+                                       "MOV result.color, temp;\n"
+                                       "END\n";
+       unsigned int prog;
+
        glGetBooleanv(GL_BLEND, (GLboolean *)&gl_blend);
        glGetBooleanv(GL_DEPTH_TEST, (GLboolean *)&gl_depth);
 
@@ -234,7 +255,23 @@ void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture
                }
        }
 
+       if(GLEW_ARB_fragment_program)
+       {
+               glGenProgramsARB(1, &prog);
+               glEnable(GL_FRAGMENT_PROGRAM_ARB);
+
+               glBindProgramARB(GL_FRAGMENT_PROGRAM_ARB, prog);
+               glProgramStringARB(GL_FRAGMENT_PROGRAM_ARB, GL_PROGRAM_FORMAT_ASCII_ARB, (GLsizei)strlen(text), text);
+
+               // cell spacing
+               glProgramLocalParameter4fARB (GL_FRAGMENT_PROGRAM_ARB, 0, dx, dx, dx, 1.0);
+               // custom parameter for smoke style (higher = thicker)
+               glProgramLocalParameter4fARB (GL_FRAGMENT_PROGRAM_ARB, 1, 7.0, 7.0, 7.0, 1.0);
+       }
+
        GPU_texture_bind(tex, 0);
+       if(tex_shadow)
+               GPU_texture_bind(tex_shadow, 1);
 
        if (!GLEW_ARB_texture_non_power_of_two) {
                cor[0] = (float)res[0]/(float)larger_pow2(res[0]);
@@ -289,8 +326,17 @@ void draw_volume(Scene *scene, ARegion *ar, View3D *v3d, Base *base, GPUTexture
                n++;
        }
 
+       if(tex_shadow)
+               GPU_texture_unbind(tex_shadow);
        GPU_texture_unbind(tex);
 
+       if(GLEW_ARB_fragment_program)
+       {
+               glDisable(GL_FRAGMENT_PROGRAM_ARB);
+               glDeleteProgramsARB(1, &prog);
+       }
+
+
        MEM_freeN(points);
 
        if(!gl_blend)
index e5e85cf9d16ed19d719f2978b4d438452d021ad7..00b0b5c4fd134ca127be7067728a5f6d91c9dfb5 100644 (file)
@@ -157,7 +157,7 @@ ARegion *view3d_has_buttons_region(ScrArea *sa);
 ARegion *view3d_has_tools_region(ScrArea *sa);
 
 /* draw_volume.c */
-void draw_volume(struct Scene *scene, struct ARegion *ar, struct View3D *v3d, struct Base *base, struct GPUTexture *tex, int res[3]);
+void draw_volume(struct Scene *scene, struct ARegion *ar, struct View3D *v3d, struct Base *base, struct GPUTexture *tex, int res[3], float dx, struct GPUTexture *tex_shadow);
 
 
 #endif /* ED_VIEW3D_INTERN_H */
index 85f4233a251c929dfef574b4bf7a31aa5e8fac8c..279596e5ad737cd9c57aad95d3c6daaad63e6df2 100644 (file)
@@ -28,7 +28,7 @@ FILE(GLOB SRC intern/*.c)
 
 SET(INC 
        . ../blenlib ../blenkernel ../makesdna ../include
-       ../../../extern/glew/include ../../../intern/guardedalloc ../imbuf)
+       ../../../extern/glew/include ../../../intern/guardedalloc ../../../intern/smoke/extern ../imbuf)
 
 BLENDERLIB(bf_gpu "${SRC}" "${INC}")
 
index 00d0e3131e561523dae12e90cdeb8f2d276eda30..fabe1420e83dae4ecdbd7b03dddf9f3176e5f6bd 100644 (file)
@@ -115,7 +115,7 @@ void GPU_free_images(void);
 
 /* smoke drawing functions */
 void GPU_free_smoke(struct SmokeModifierData *smd);
-void GPU_create_smoke(struct SmokeModifierData *smd);
+void GPU_create_smoke(struct SmokeModifierData *smd, int highres);
 
 #ifdef __cplusplus
 }
index a81c7e034558831d64bbedeb801c76d6802f5856..75e8073aafd68df65af3ece544e46c8bc3127ef7 100644 (file)
@@ -64,6 +64,8 @@
 #include "GPU_material.h"
 #include "GPU_draw.h"
 
+#include "smoke_API.h"
+
 /* These are some obscure rendering functions shared between the
  * game engine and the blender, in this module to avoid duplicaten
  * and abstract them away from the rest a bit */
@@ -754,13 +756,21 @@ void GPU_free_smoke(SmokeModifierData *smd)
                if(smd->domain->tex)
                        GPU_texture_free(smd->domain->tex);
                smd->domain->tex = NULL;
+
+               if(smd->domain->tex_shadow)
+                       GPU_texture_free(smd->domain->tex_shadow);
+               smd->domain->tex_shadow = NULL;
        }
 }
 
-void GPU_create_smoke(SmokeModifierData *smd)
+void GPU_create_smoke(SmokeModifierData *smd, int highres)
 {
-       if(smd->type & MOD_SMOKE_TYPE_DOMAIN && smd->domain && !smd->domain->tex)
-               smd->domain->tex = GPU_texture_create_3D(smd->domain->res[0], smd->domain->res[1], smd->domain->res[2], smd->domain->view3d);
+       if(smd->type & MOD_SMOKE_TYPE_DOMAIN && smd->domain && !smd->domain->tex && !highres)
+               smd->domain->tex = GPU_texture_create_3D(smd->domain->res[0], smd->domain->res[1], smd->domain->res[2], smoke_get_density(smd->domain->fluid));
+       else if(smd->type & MOD_SMOKE_TYPE_DOMAIN && smd->domain && !smd->domain->tex && highres)
+               smd->domain->tex = GPU_texture_create_3D(smd->domain->res_wt[0], smd->domain->res_wt[1], smd->domain->res_wt[2], smoke_turbulence_get_density(smd->domain->wt));
+
+       smd->domain->tex_shadow = GPU_texture_create_3D(smd->domain->res[0], smd->domain->res[1], smd->domain->res[2], smd->domain->shadow);
 }
 
 void GPU_free_image(Image *ima)
index 850b46dc28cd612acc3fe63c2bf1b12a4ddc1ead..d7b54e425fd0c780b3982a4899d731fc468b5ce4 100644 (file)
@@ -346,17 +346,17 @@ GPUTexture *GPU_texture_create_3D(int w, int h, int depth, float *fpixels)
        tex->number = 0;
        glBindTexture(tex->target, tex->bindcode);
 
-       type = GL_UNSIGNED_BYTE;
-       format = GL_RGBA;
-       internalformat = GL_RGBA8;
+       type = GL_FLOAT; // GL_UNSIGNED_BYTE
+       format = GL_RED;
+       internalformat = GL_RED;
 
-       if (fpixels)
-               pixels = GPU_texture_convert_pixels(w*h*depth, fpixels);
+       //if (fpixels)
+       //      pixels = GPU_texture_convert_pixels(w*h*depth, fpixels);
 
        glTexImage3D(tex->target, 0, internalformat, tex->w, tex->h, tex->depth, 0, format, type, 0);
 
        if (fpixels) {
-               glTexSubImage3D(tex->target, 0, 0, 0, 0, w, h, depth, format, type, pixels);
+               glTexSubImage3D(tex->target, 0, 0, 0, 0, w, h, depth, format, type, fpixels);
        }
 
        glTexParameterfv(GL_TEXTURE_3D, GL_TEXTURE_BORDER_COLOR, vfBorderColor);
index 7c6c7fab9e493e616fc330a46f9f369b8cc87944..4e4714cdaa1d99047f97cbf099e5da47ccd3fc7d 100644 (file)
@@ -38,9 +38,8 @@
 #define MOD_SMOKE_NOISEWAVE (1<<0)
 #define MOD_SMOKE_NOISEFFT (1<<1)
 #define MOD_SMOKE_NOISECURL (1<<2)
-
 /* viewsettings */
-#define MOD_SMOKE_SHOWHIGHRES (1<<0) /* show high resolution */
+#define MOD_SMOKE_VIEW_SHOWBIG (1<<0)
 
 typedef struct SmokeDomainSettings {
        struct SmokeModifierData *smd; /* for fast RNA access */
@@ -48,35 +47,34 @@ typedef struct SmokeDomainSettings {
        struct Group *fluid_group;
        struct Group *eff_group; // effector group for e.g. wind force
        struct Group *coll_group; // collision objects group
+       struct WTURBULENCE *wt; // WTURBULENCE object, if active
        struct GPUTexture *tex;
-       float *view3d; /* voxel data for display */
-       unsigned int v3dnum; /* number of frame in view3d buffer */
+       struct GPUTexture *tex_wt;
+       struct GPUTexture *tex_shadow;
+       float *shadow;
        float p0[3]; /* start point of BB */
        float p1[3]; /* end point of BB */
        float dx; /* edge length of one cell */
-       float firstframe;
-       float lastframe;
+       float omega; /* smoke color - from 0 to 1 */
        float temp; /* fluid temperature */
        float tempAmb; /* ambient temperature */
        float alpha;
        float beta;
        int res[3]; /* domain resolution */
+       int amplify; /* wavelet amplification */
        int maxres; /* longest axis on the BB gets this resolution assigned */
        int flags; /* show up-res or low res, etc */
+       int pad; 
        int viewsettings;
+       short noise; /* noise type: wave, curl, anisotropic */
        short diss_percent; 
-       short pad;
        int diss_speed;/* in frames */
-       struct PointCache *point_cache[2];      /* definition is in DNA_object_force.h */
-       struct ListBase ptcaches[2];
-       struct WTURBULENCE *wt; // WTURBULENCE object, if active
-       int pad3;
        float strength;
        int res_wt[3];
-       int maxres_wt;
-       short noise; /* noise type: wave, curl, anisotropic */
-       short pad2;
-       int amplify;
+       float dx_wt;
+       int v3dnum;
+       struct PointCache *point_cache[2];      /* definition is in DNA_object_force.h */
+       struct ListBase ptcaches[2];
 } SmokeDomainSettings;
 
 
index b3192b110f4f78d93a6be4925c53b602d7ad85f4..943129c716936bb5dd78999863159f952407fb50 100644 (file)
@@ -46,7 +46,6 @@
 #include "BKE_context.h"
 #include "BKE_depsgraph.h"
 #include "BKE_particle.h"
-#include "BKE_pointcache.h"
 
 #include "ED_object.h"
 
@@ -79,15 +78,6 @@ static void rna_Smoke_reset_dependancy(bContext *C, PointerRNA *ptr)
        rna_Smoke_dependency_update(C, ptr);
 }
 
-#if 0
-static void rna_Smoke_redraw(bContext *C, PointerRNA *ptr)
-{
-       SmokeDomainSettings *settings = (SmokeDomainSettings*)ptr->data;
-
-       settings->flags |= MOD_SMOKE_VIEW_REDRAWNICE;
-}
-#endif
-
 static char *rna_SmokeDomainSettings_path(PointerRNA *ptr)
 {
        SmokeDomainSettings *settings = (SmokeDomainSettings*)ptr->data;
@@ -139,6 +129,29 @@ static void rna_def_smoke_domain_settings(BlenderRNA *brna)
        RNA_def_property_ui_text(prop, "Max Res", "Maximal resolution used in the fluid domain.");
        RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
 
+       prop= RNA_def_property(srna, "amplify", PROP_INT, PROP_NONE);
+       RNA_def_property_int_sdna(prop, NULL, "amplify");
+       RNA_def_property_range(prop, 1, 10);
+       RNA_def_property_ui_range(prop, 1, 10, 1, 0);
+       RNA_def_property_ui_text(prop, "Amplification", "Enhance the resolution of smoke by this factor using noise.");
+       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
+
+       prop= RNA_def_property(srna, "highres", PROP_BOOLEAN, PROP_NONE);
+       RNA_def_property_boolean_sdna(prop, NULL, "flags", MOD_SMOKE_HIGHRES);
+       RNA_def_property_ui_text(prop, "High res", "Enable high resolution (using amplification).");
+       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
+
+       prop= RNA_def_property(srna, "viewhighres", PROP_BOOLEAN, PROP_NONE);
+       RNA_def_property_boolean_sdna(prop, NULL, "viewsettings", MOD_SMOKE_VIEW_SHOWBIG);
+       RNA_def_property_ui_text(prop, "Show High Resolution", "Show high resolution (using amplification).");
+       RNA_def_property_update(prop, NC_OBJECT|ND_DRAW, NULL);
+
+       prop= RNA_def_property(srna, "noise_type", PROP_ENUM, PROP_NONE);
+       RNA_def_property_enum_sdna(prop, NULL, "noise");
+       RNA_def_property_enum_items(prop, prop_noise_type_items);
+       RNA_def_property_ui_text(prop, "Noise Method", "Noise method which is used for creating the high resolution");
+       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
+
        prop= RNA_def_property(srna, "alpha", PROP_FLOAT, PROP_NONE);
        RNA_def_property_float_sdna(prop, NULL, "alpha");
        RNA_def_property_range(prop, -5.0, 5.0);
@@ -174,61 +187,36 @@ static void rna_def_smoke_domain_settings(BlenderRNA *brna)
        RNA_def_property_ui_text(prop, "Effector Group", "Limit effectors to this group.");
        RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset_dependancy");
 
+       prop= RNA_def_property(srna, "strength", PROP_FLOAT, PROP_NONE);
+       RNA_def_property_float_sdna(prop, NULL, "strength");
+       RNA_def_property_range(prop, 1.0, 10.0);
+       RNA_def_property_ui_range(prop, 1.0, 10.0, 1, 2);
+       RNA_def_property_ui_text(prop, "Strength", "Strength of wavelet noise");
+       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
+
        prop= RNA_def_property(srna, "dissolve_speed", PROP_INT, PROP_NONE);
        RNA_def_property_int_sdna(prop, NULL, "diss_speed");
        RNA_def_property_range(prop, 1.0, 100.0);
        RNA_def_property_ui_range(prop, 1.0, 1000.0, 1, 0);
        RNA_def_property_ui_text(prop, "Dissolve Speed", "Dissolve Speed");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
-
-       prop= RNA_def_property(srna, "highres", PROP_BOOLEAN, PROP_NONE);
-       RNA_def_property_boolean_sdna(prop, NULL, "flags", MOD_SMOKE_HIGHRES);
-       RNA_def_property_ui_text(prop, "High Resolution Smoke", "Enable high resolution smoke");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL);
+       RNA_def_property_update(prop, 0, NULL);
 
        prop= RNA_def_property(srna, "dissolve_smoke", PROP_BOOLEAN, PROP_NONE);
        RNA_def_property_boolean_sdna(prop, NULL, "flags", MOD_SMOKE_DISSOLVE);
        RNA_def_property_ui_text(prop, "Dissolve Smoke", "Enable smoke to disappear over time.");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
+       RNA_def_property_update(prop, 0, NULL);
 
        prop= RNA_def_property(srna, "dissolve_smoke_log", PROP_BOOLEAN, PROP_NONE);
        RNA_def_property_boolean_sdna(prop, NULL, "flags", MOD_SMOKE_DISSOLVE_LOG);
        RNA_def_property_ui_text(prop, "Logarithmic dissolve", "Using 1/x ");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
+       RNA_def_property_update(prop, 0, NULL);
 
-       prop= RNA_def_property(srna, "point_cache", PROP_POINTER, PROP_NEVER_NULL);
+       prop= RNA_def_property(srna, "point_cache_low", PROP_POINTER, PROP_NEVER_NULL);
        RNA_def_property_pointer_sdna(prop, NULL, "point_cache[0]");
-       RNA_def_property_struct_type(prop, "PointCache");
        RNA_def_property_ui_text(prop, "Point Cache", "");
 
-       prop= RNA_def_property(srna, "show_highres", PROP_BOOLEAN, PROP_NONE);
-       RNA_def_property_boolean_sdna(prop, NULL, "viewsettings", MOD_SMOKE_SHOWHIGHRES);
-       RNA_def_property_ui_text(prop, "High res", "Show high resolution (using amplification).");
-       RNA_def_property_update(prop, NC_OBJECT|ND_DRAW, NULL);
-
-       prop= RNA_def_property(srna, "noise_type", PROP_ENUM, PROP_NONE);
-       RNA_def_property_enum_sdna(prop, NULL, "noise");
-       RNA_def_property_enum_items(prop, prop_noise_type_items);
-       RNA_def_property_ui_text(prop, "Noise Method", "Noise method which is used for creating the high resolution");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
-
-       prop= RNA_def_property(srna, "amplify", PROP_INT, PROP_NONE);
-       RNA_def_property_int_sdna(prop, NULL, "amplify");
-       RNA_def_property_range(prop, 1, 10);
-       RNA_def_property_ui_range(prop, 1, 10, 1, 0);
-       RNA_def_property_ui_text(prop, "Amplification", "Enhance the resolution of smoke by this factor using noise.");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
-
-       prop= RNA_def_property(srna, "strength", PROP_FLOAT, PROP_NONE);
-       RNA_def_property_float_sdna(prop, NULL, "strength");
-       RNA_def_property_range(prop, 1.0, 10.0);
-       RNA_def_property_ui_range(prop, 1.0, 10.0, 1, 2);
-       RNA_def_property_ui_text(prop, "Strength", "Strength of wavelet noise");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset");
-
-       prop= RNA_def_property(srna, "point_cache_turbulence", PROP_POINTER, PROP_NEVER_NULL);
+       prop= RNA_def_property(srna, "point_cache_high", PROP_POINTER, PROP_NEVER_NULL);
        RNA_def_property_pointer_sdna(prop, NULL, "point_cache[1]");
-       RNA_def_property_struct_type(prop, "PointCache");
        RNA_def_property_ui_text(prop, "Point Cache", "");
 }
 
@@ -247,26 +235,26 @@ static void rna_def_smoke_flow_settings(BlenderRNA *brna)
        RNA_def_property_range(prop, 0.001, 1);
        RNA_def_property_ui_range(prop, 0.001, 1.0, 1.0, 4);
        RNA_def_property_ui_text(prop, "Density", "");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL);
+       RNA_def_property_update(prop, 0, NULL); // NC_OBJECT|ND_MODIFIER
 
        prop= RNA_def_property(srna, "temperature", PROP_FLOAT, PROP_NONE);
        RNA_def_property_float_sdna(prop, NULL, "temp");
        RNA_def_property_range(prop, -10, 10);
        RNA_def_property_ui_range(prop, -10, 10, 1, 1);
        RNA_def_property_ui_text(prop, "Temp. Diff.", "Temperature difference to ambientt temperature.");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL);
+       RNA_def_property_update(prop, 0, NULL);
        
        prop= RNA_def_property(srna, "psys", PROP_POINTER, PROP_NONE);
        RNA_def_property_pointer_sdna(prop, NULL, "psys");
        RNA_def_property_struct_type(prop, "ParticleSystem");
        RNA_def_property_flag(prop, PROP_EDITABLE);
        RNA_def_property_ui_text(prop, "Particle Systems", "Particle systems emitted from the object.");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, "rna_Smoke_reset_dependancy");
+       RNA_def_property_update(prop, 0, "rna_Smoke_reset_dependancy");
 
        prop= RNA_def_property(srna, "outflow", PROP_BOOLEAN, PROP_NONE);
        RNA_def_property_boolean_sdna(prop, NULL, "type", MOD_SMOKE_FLOW_TYPE_OUTFLOW);
        RNA_def_property_ui_text(prop, "Outflow", "Deletes smoke from simulation");
-       RNA_def_property_update(prop, NC_OBJECT|ND_MODIFIER, NULL);
+       RNA_def_property_update(prop, 0, NULL);
 }
 
 static void rna_def_smoke_coll_settings(BlenderRNA *brna)