[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [26250] trunk/blender: Smoke: The well known Miika H?\195?\164m?\195?\164l?\195?\164inen (aka MiikaH) patch (http://blenderartists.org/forum/showthread.php?t=158317&page=42)
Daniel Genrich
daniel.genrich at gmx.net
Mon Jan 25 16:10:14 CET 2010
Revision: 26250
http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=26250
Author: genscher
Date: 2010-01-25 16:10:14 +0100 (Mon, 25 Jan 2010)
Log Message:
-----------
Smoke: The well known Miika H?\195?\164m?\195?\164l?\195?\164inen (aka MiikaH) patch (http://blenderartists.org/forum/showthread.php?t=158317&page=42)
* Better (and windows enabled) OpenMP handling (> 2x-5x speed)
* More Volumetric Texture mapping options (heat, etc) <-- Matt if that's not to your liking, just revert that part, it's separate anyway
* Initial velocity taken from particle settings (no more slow starting)
* Option to select compression method (there seem to be a bug in my high compression usage, at least it's been reported to result in exploding smoke - better use low compression for the time being)
It's been tested since a while but as usual please report any (new!) bugs. ;-)
Modified Paths:
--------------
trunk/blender/intern/smoke/intern/EIGENVALUE_HELPER.h
trunk/blender/intern/smoke/intern/FLUID_3D.cpp
trunk/blender/intern/smoke/intern/FLUID_3D.h
trunk/blender/intern/smoke/intern/FLUID_3D_SOLVERS.cpp
trunk/blender/intern/smoke/intern/FLUID_3D_STATIC.cpp
trunk/blender/intern/smoke/intern/LU_HELPER.h
trunk/blender/intern/smoke/intern/WTURBULENCE.cpp
trunk/blender/intern/smoke/intern/smoke_API.cpp
trunk/blender/release/scripts/ui/properties_physics_smoke.py
trunk/blender/release/scripts/ui/properties_texture.py
trunk/blender/source/blender/blenkernel/intern/pointcache.c
trunk/blender/source/blender/blenkernel/intern/smoke.c
trunk/blender/source/blender/makesdna/DNA_smoke_types.h
trunk/blender/source/blender/makesdna/DNA_texture_types.h
trunk/blender/source/blender/makesrna/intern/rna_smoke.c
trunk/blender/source/blender/makesrna/intern/rna_texture.c
trunk/blender/source/blender/render/intern/source/voxeldata.c
Added Paths:
-----------
trunk/blender/intern/smoke/intern/EIGENVALUE_HELPER.cpp
trunk/blender/intern/smoke/intern/LU_HELPER.cpp
Added: trunk/blender/intern/smoke/intern/EIGENVALUE_HELPER.cpp
===================================================================
--- trunk/blender/intern/smoke/intern/EIGENVALUE_HELPER.cpp (rev 0)
+++ trunk/blender/intern/smoke/intern/EIGENVALUE_HELPER.cpp 2010-01-25 15:10:14 UTC (rev 26250)
@@ -0,0 +1,885 @@
+
+#include "EIGENVALUE_HELPER.h"
+
+
+void Eigentred2(sEigenvalue& eval) {
+
+ // This is derived from the Algol procedures tred2 by
+ // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+
+ int n=eval.n;
+
+ for (int j = 0; j < n; j++) {
+ eval.d[j] = eval.V[n-1][j];
+ }
+
+ // Householder reduction to tridiagonal form.
+
+ for (int i = n-1; i > 0; i--) {
+
+ // Scale to avoid under/overflow.
+
+ float scale = 0.0;
+ float h = 0.0;
+ for (int k = 0; k < i; k++) {
+ scale = scale + fabs(eval.d[k]);
+ }
+ if (scale == 0.0) {
+ eval.e[i] = eval.d[i-1];
+ for (int j = 0; j < i; j++) {
+ eval.d[j] = eval.V[i-1][j];
+ eval.V[i][j] = 0.0;
+ eval.V[j][i] = 0.0;
+ }
+ } else {
+
+ // Generate Householder vector.
+
+ for (int k = 0; k < i; k++) {
+ eval.d[k] /= scale;
+ h += eval.d[k] * eval.d[k];
+ }
+ float f = eval.d[i-1];
+ float g = sqrt(h);
+ if (f > 0) {
+ g = -g;
+ }
+ eval.e[i] = scale * g;
+ h = h - f * g;
+ eval.d[i-1] = f - g;
+ for (int j = 0; j < i; j++) {
+ eval.e[j] = 0.0;
+ }
+
+ // Apply similarity transformation to remaining columns.
+
+ for (int j = 0; j < i; j++) {
+ f = eval.d[j];
+ eval.V[j][i] = f;
+ g = eval.e[j] + eval.V[j][j] * f;
+ for (int k = j+1; k <= i-1; k++) {
+ g += eval.V[k][j] * eval.d[k];
+ eval.e[k] += eval.V[k][j] * f;
+ }
+ eval.e[j] = g;
+ }
+ f = 0.0;
+ for (int j = 0; j < i; j++) {
+ eval.e[j] /= h;
+ f += eval.e[j] * eval.d[j];
+ }
+ float hh = f / (h + h);
+ for (int j = 0; j < i; j++) {
+ eval.e[j] -= hh * eval.d[j];
+ }
+ for (int j = 0; j < i; j++) {
+ f = eval.d[j];
+ g = eval.e[j];
+ for (int k = j; k <= i-1; k++) {
+ eval.V[k][j] -= (f * eval.e[k] + g * eval.d[k]);
+ }
+ eval.d[j] = eval.V[i-1][j];
+ eval.V[i][j] = 0.0;
+ }
+ }
+ eval.d[i] = h;
+ }
+
+ // Accumulate transformations.
+
+ for (int i = 0; i < n-1; i++) {
+ eval.V[n-1][i] = eval.V[i][i];
+ eval.V[i][i] = 1.0;
+ float h = eval.d[i+1];
+ if (h != 0.0) {
+ for (int k = 0; k <= i; k++) {
+ eval.d[k] = eval.V[k][i+1] / h;
+ }
+ for (int j = 0; j <= i; j++) {
+ float g = 0.0;
+ for (int k = 0; k <= i; k++) {
+ g += eval.V[k][i+1] * eval.V[k][j];
+ }
+ for (int k = 0; k <= i; k++) {
+ eval.V[k][j] -= g * eval.d[k];
+ }
+ }
+ }
+ for (int k = 0; k <= i; k++) {
+ eval.V[k][i+1] = 0.0;
+ }
+ }
+ for (int j = 0; j < n; j++) {
+ eval.d[j] = eval.V[n-1][j];
+ eval.V[n-1][j] = 0.0;
+ }
+ eval.V[n-1][n-1] = 1.0;
+ eval.e[0] = 0.0;
+}
+
+void Eigencdiv(sEigenvalue& eval, float xr, float xi, float yr, float yi) {
+ float r,d;
+ if (fabs(yr) > fabs(yi)) {
+ r = yi/yr;
+ d = yr + r*yi;
+ eval.cdivr = (xr + r*xi)/d;
+ eval.cdivi = (xi - r*xr)/d;
+ } else {
+ r = yr/yi;
+ d = yi + r*yr;
+ eval.cdivr = (r*xr + xi)/d;
+ eval.cdivi = (r*xi - xr)/d;
+ }
+ }
+
+void Eigentql2 (sEigenvalue& eval) {
+
+ // This is derived from the Algol procedures tql2, by
+ // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
+ // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutine in EISPACK.
+
+ int n=eval.n;
+
+ for (int i = 1; i < n; i++) {
+ eval.e[i-1] = eval.e[i];
+ }
+ eval.e[n-1] = 0.0;
+
+ float f = 0.0;
+ float tst1 = 0.0;
+ float eps = pow(2.0,-52.0);
+ for (int l = 0; l < n; l++) {
+
+ // Find small subdiagonal element
+
+ tst1 = max(tst1,fabs(eval.d[l]) + fabs(eval.e[l]));
+ int m = l;
+
+ // Original while-loop from Java code
+ while (m < n) {
+ if (fabs(eval.e[m]) <= eps*tst1) {
+ break;
+ }
+ m++;
+ }
+
+
+ // If m == l, d[l] is an eigenvalue,
+ // otherwise, iterate.
+
+ if (m > l) {
+ int iter = 0;
+ do {
+ iter = iter + 1; // (Could check iteration count here.)
+
+ // Compute implicit shift
+
+ float g = eval.d[l];
+ float p = (eval.d[l+1] - g) / (2.0 * eval.e[l]);
+ float r = hypot(p,1.0);
+ if (p < 0) {
+ r = -r;
+ }
+ eval.d[l] = eval.e[l] / (p + r);
+ eval.d[l+1] = eval.e[l] * (p + r);
+ float dl1 = eval.d[l+1];
+ float h = g - eval.d[l];
+ for (int i = l+2; i < n; i++) {
+ eval.d[i] -= h;
+ }
+ f = f + h;
+
+ // Implicit QL transformation.
+
+ p = eval.d[m];
+ float c = 1.0;
+ float c2 = c;
+ float c3 = c;
+ float el1 = eval.e[l+1];
+ float s = 0.0;
+ float s2 = 0.0;
+ for (int i = m-1; i >= l; i--) {
+ c3 = c2;
+ c2 = c;
+ s2 = s;
+ g = c * eval.e[i];
+ h = c * p;
+ r = hypot(p,eval.e[i]);
+ eval.e[i+1] = s * r;
+ s = eval.e[i] / r;
+ c = p / r;
+ p = c * eval.d[i] - s * g;
+ eval.d[i+1] = h + s * (c * g + s * eval.d[i]);
+
+ // Accumulate transformation.
+
+ for (int k = 0; k < n; k++) {
+ h = eval.V[k][i+1];
+ eval.V[k][i+1] = s * eval.V[k][i] + c * h;
+ eval.V[k][i] = c * eval.V[k][i] - s * h;
+ }
+ }
+ p = -s * s2 * c3 * el1 * eval.e[l] / dl1;
+ eval.e[l] = s * p;
+ eval.d[l] = c * p;
+
+ // Check for convergence.
+
+ } while (fabs(eval.e[l]) > eps*tst1);
+ }
+ eval.d[l] = eval.d[l] + f;
+ eval.e[l] = 0.0;
+ }
+
+ // Sort eigenvalues and corresponding vectors.
+
+ for (int i = 0; i < n-1; i++) {
+ int k = i;
+ float p = eval.d[i];
+ for (int j = i+1; j < n; j++) {
+ if (eval.d[j] < p) {
+ k = j;
+ p = eval.d[j];
+ }
+ }
+ if (k != i) {
+ eval.d[k] = eval.d[i];
+ eval.d[i] = p;
+ for (int j = 0; j < n; j++) {
+ p = eval.V[j][i];
+ eval.V[j][i] = eval.V[j][k];
+ eval.V[j][k] = p;
+ }
+ }
+ }
+}
+
+void Eigenorthes (sEigenvalue& eval) {
+
+ // This is derived from the Algol procedures orthes and ortran,
+ // by Martin and Wilkinson, Handbook for Auto. Comp.,
+ // Vol.ii-Linear Algebra, and the corresponding
+ // Fortran subroutines in EISPACK.
+
+ int n=eval.n;
+
+ int low = 0;
+ int high = n-1;
+
+ for (int m = low+1; m <= high-1; m++) {
+
+ // Scale column.
+
+ float scale = 0.0;
+ for (int i = m; i <= high; i++) {
+ scale = scale + fabs(eval.H[i][m-1]);
+ }
+ if (scale != 0.0) {
+
+ // Compute Householder transformation.
+
+ float h = 0.0;
+ for (int i = high; i >= m; i--) {
+ eval.ort[i] = eval.H[i][m-1]/scale;
+ h += eval.ort[i] * eval.ort[i];
+ }
+ float g = sqrt(h);
+ if (eval.ort[m] > 0) {
+ g = -g;
+ }
+ h = h - eval.ort[m] * g;
+ eval.ort[m] = eval.ort[m] - g;
+
+ // Apply Householder similarity transformation
+ // H = (I-u*u'/h)*H*(I-u*u')/h)
+
+ for (int j = m; j < n; j++) {
+ float f = 0.0;
+ for (int i = high; i >= m; i--) {
+ f += eval.ort[i]*eval.H[i][j];
+ }
+ f = f/h;
+ for (int i = m; i <= high; i++) {
+ eval.H[i][j] -= f*eval.ort[i];
+ }
+ }
+
+ for (int i = 0; i <= high; i++) {
+ float f = 0.0;
+ for (int j = high; j >= m; j--) {
+ f += eval.ort[j]*eval.H[i][j];
+ }
+ f = f/h;
+ for (int j = m; j <= high; j++) {
+ eval.H[i][j] -= f*eval.ort[j];
+ }
+ }
+ eval.ort[m] = scale*eval.ort[m];
+ eval.H[m][m-1] = scale*g;
+ }
+ }
+
+ // Accumulate transformations (Algol's ortran).
+
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ eval.V[i][j] = (i == j ? 1.0 : 0.0);
+ }
+ }
+
+ for (int m = high-1; m >= low+1; m--) {
+ if (eval.H[m][m-1] != 0.0) {
+ for (int i = m+1; i <= high; i++) {
+ eval.ort[i] = eval.H[i][m-1];
+ }
+ for (int j = m; j <= high; j++) {
+ float g = 0.0;
+ for (int i = m; i <= high; i++) {
+ g += eval.ort[i] * eval.V[i][j];
+ }
+ // Double division avoids possible underflow
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list