Blender  V2.59
smoke.c
Go to the documentation of this file.
00001 /*
00002  * smoke.c
00003  *
00004  * $Id: smoke.c 38372 2011-07-13 18:40:21Z campbellbarton $
00005  *
00006  * ***** BEGIN GPL LICENSE BLOCK *****
00007  *
00008  * This program is free software; you can redistribute it and/or
00009  * modify it under the terms of the GNU General Public License
00010  * as published by the Free Software Foundation; either version 2
00011  * of the License, or (at your option) any later version.
00012  *
00013  * This program is distributed in the hope that it will be useful,
00014  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00015  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00016  * GNU General Public License for more details.
00017  *
00018  * You should have received a copy of the GNU General Public License
00019  * along with this program; if not, write to the Free Software Foundation,
00020  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
00021  *
00022  * The Original Code is Copyright (C) Blender Foundation.
00023  * All rights reserved.
00024  *
00025  * The Original Code is: all of this file.
00026  *
00027  * Contributor(s): Daniel Genrich
00028  *
00029  * ***** END GPL LICENSE BLOCK *****
00030  */
00031 
00037 /* Part of the code copied from elbeem fluid library, copyright by Nils Thuerey */
00038 
00039 #include <GL/glew.h>
00040 
00041 #include "MEM_guardedalloc.h"
00042 
00043 #include <float.h>
00044 #include <math.h>
00045 #include <stdio.h>
00046 #include <string.h> /* memset */
00047 
00048 #include "BLI_linklist.h"
00049 #include "BLI_rand.h"
00050 #include "BLI_jitter.h"
00051 #include "BLI_blenlib.h"
00052 #include "BLI_math.h"
00053 #include "BLI_edgehash.h"
00054 #include "BLI_kdtree.h"
00055 #include "BLI_kdopbvh.h"
00056 #include "BLI_utildefines.h"
00057 
00058 #include "BKE_bvhutils.h"
00059 #include "BKE_cdderivedmesh.h"
00060 #include "BKE_customdata.h"
00061 #include "BKE_DerivedMesh.h"
00062 #include "BKE_effect.h"
00063 #include "BKE_modifier.h"
00064 #include "BKE_particle.h"
00065 #include "BKE_pointcache.h"
00066 #include "BKE_smoke.h"
00067 
00068 
00069 #include "DNA_customdata_types.h"
00070 #include "DNA_group_types.h"
00071 #include "DNA_lamp_types.h"
00072 #include "DNA_mesh_types.h"
00073 #include "DNA_meshdata_types.h"
00074 #include "DNA_modifier_types.h"
00075 #include "DNA_object_types.h"
00076 #include "DNA_particle_types.h"
00077 #include "DNA_scene_types.h"
00078 #include "DNA_smoke_types.h"
00079 
00080 #include "smoke_API.h"
00081 
00082 #include "BKE_smoke.h"
00083 
00084 #ifdef _WIN32
00085 #include <time.h>
00086 #include <stdio.h>
00087 #include <conio.h>
00088 #include <windows.h>
00089 
00090 static LARGE_INTEGER liFrequency;
00091 static LARGE_INTEGER liStartTime;
00092 static LARGE_INTEGER liCurrentTime;
00093 
00094 static void tstart ( void )
00095 {
00096         QueryPerformanceFrequency ( &liFrequency );
00097         QueryPerformanceCounter ( &liStartTime );
00098 }
00099 static void tend ( void )
00100 {
00101         QueryPerformanceCounter ( &liCurrentTime );
00102 }
00103 static double tval( void )
00104 {
00105         return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart)* (double)1000.0/(double)liFrequency.QuadPart ));
00106 }
00107 #else
00108 #include <sys/time.h>
00109 static struct timeval _tstart, _tend;
00110 static struct timezone tz;
00111 static void tstart ( void )
00112 {
00113         gettimeofday ( &_tstart, &tz );
00114 }
00115 static void tend ( void )
00116 {
00117         gettimeofday ( &_tend,&tz );
00118 }
00119 
00120 #if 0 // unused
00121 static double tval()
00122 {
00123         double t1, t2;
00124         t1 = ( double ) _tstart.tv_sec*1000 + ( double ) _tstart.tv_usec/ ( 1000 );
00125         t2 = ( double ) _tend.tv_sec*1000 + ( double ) _tend.tv_usec/ ( 1000 );
00126         return t2-t1;
00127 }
00128 #endif
00129 #endif
00130 
00131 struct Object;
00132 struct Scene;
00133 struct DerivedMesh;
00134 struct SmokeModifierData;
00135 
00136 // forward declerations
00137 static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int correct);
00138 void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *tris, int numfaces, int numtris, int **tridivs, float cell_len);
00139 static void fill_scs_points(Object *ob, DerivedMesh *dm, SmokeCollSettings *scs);
00140 
00141 #define TRI_UVOFFSET (1./4.)
00142 
00143 /* Stubs to use when smoke is disabled */
00144 #ifndef WITH_SMOKE
00145 struct WTURBULENCE *smoke_turbulence_init(int *UNUSED(res), int UNUSED(amplify), int UNUSED(noisetype)) { return NULL; }
00146 struct FLUID_3D *smoke_init(int *UNUSED(res), float *UNUSED(p0)) { return NULL; }
00147 void smoke_free(struct FLUID_3D *UNUSED(fluid)) {}
00148 void smoke_turbulence_free(struct WTURBULENCE *UNUSED(wt)) {}
00149 void smoke_initWaveletBlenderRNA(struct WTURBULENCE *UNUSED(wt), float *UNUSED(strength)) {}
00150 void smoke_initBlenderRNA(struct FLUID_3D *UNUSED(fluid), float *UNUSED(alpha), float *UNUSED(beta), float *UNUSED(dt_factor), float *UNUSED(vorticity), int *UNUSED(border_colli)) {}
00151 long long smoke_get_mem_req(int UNUSED(xres), int UNUSED(yres), int UNUSED(zres), int UNUSED(amplify)) { return 0; }
00152 void smokeModifier_do(SmokeModifierData *UNUSED(smd), Scene *UNUSED(scene), Object *UNUSED(ob), DerivedMesh *UNUSED(dm)) {}
00153 #endif // WITH_SMOKE
00154 
00155 
00156 static int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, DerivedMesh *dm)
00157 {
00158         if((smd->type & MOD_SMOKE_TYPE_DOMAIN) && smd->domain && !smd->domain->fluid)
00159         {
00160                 size_t i;
00161                 float min[3] = {FLT_MAX, FLT_MAX, FLT_MAX}, max[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
00162                 float size[3];
00163                 MVert *verts = dm->getVertArray(dm);
00164                 float scale = 0.0;
00165                 int res;                
00166 
00167                 res = smd->domain->maxres;
00168 
00169                 // get BB of domain
00170                 for(i = 0; i < dm->getNumVerts(dm); i++)
00171                 {
00172                         float tmp[3];
00173 
00174                         VECCOPY(tmp, verts[i].co);
00175                         mul_m4_v3(ob->obmat, tmp);
00176 
00177                         // min BB
00178                         min[0] = MIN2(min[0], tmp[0]);
00179                         min[1] = MIN2(min[1], tmp[1]);
00180                         min[2] = MIN2(min[2], tmp[2]);
00181 
00182                         // max BB
00183                         max[0] = MAX2(max[0], tmp[0]);
00184                         max[1] = MAX2(max[1], tmp[1]);
00185                         max[2] = MAX2(max[2], tmp[2]);
00186                 }
00187 
00188                 VECCOPY(smd->domain->p0, min);
00189                 VECCOPY(smd->domain->p1, max);
00190 
00191                 // calc other res with max_res provided
00192                 VECSUB(size, max, min);
00193 
00194                 // printf("size: %f, %f, %f\n", size[0], size[1], size[2]);
00195 
00196                 // prevent crash when initializing a plane as domain
00197                 if((size[0] < FLT_EPSILON) || (size[1] < FLT_EPSILON) || (size[2] < FLT_EPSILON))
00198                         return 0;
00199 
00200                 if(size[0] > size[1])
00201                 {
00202                         if(size[0] > size[2])
00203                         {
00204                                 scale = res / size[0];
00205                                 smd->domain->dx = size[0] / res;
00206                                 smd->domain->res[0] = res;
00207                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
00208                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
00209                         }
00210                         else
00211                         {
00212                                 scale = res / size[2];
00213                                 smd->domain->dx = size[2] / res;
00214                                 smd->domain->res[2] = res;
00215                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
00216                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
00217                         }
00218                 }
00219                 else
00220                 {
00221                         if(size[1] > size[2])
00222                         {
00223                                 scale = res / size[1];
00224                                 smd->domain->dx = size[1] / res;
00225                                 smd->domain->res[1] = res;
00226                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
00227                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
00228                         }
00229                         else
00230                         {
00231                                 scale = res / size[2];
00232                                 smd->domain->dx = size[2] / res;
00233                                 smd->domain->res[2] = res;
00234                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
00235                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
00236                         }
00237                 }
00238 
00239                 // printf("smd->domain->dx: %f\n", smd->domain->dx);
00240 
00241                 // TODO: put in failsafe if res<=0 - dg
00242 
00243                 // printf("res[0]: %d, res[1]: %d, res[2]: %d\n", smd->domain->res[0], smd->domain->res[1], smd->domain->res[2]);
00244                 // dt max is 0.1
00245                 smd->domain->fluid = smoke_init(smd->domain->res, smd->domain->p0);
00246                 smd->time = scene->r.cfra;
00247 
00248                 if(smd->domain->flags & MOD_SMOKE_HIGHRES)
00249                 {
00250                         smd->domain->wt = smoke_turbulence_init(smd->domain->res, smd->domain->amplify + 1, smd->domain->noise);
00251                         smd->domain->res_wt[0] = smd->domain->res[0] * (smd->domain->amplify + 1);
00252                         smd->domain->res_wt[1] = smd->domain->res[1] * (smd->domain->amplify + 1);                      
00253                         smd->domain->res_wt[2] = smd->domain->res[2] * (smd->domain->amplify + 1);                      
00254                         smd->domain->dx_wt = smd->domain->dx / (smd->domain->amplify + 1);              
00255                         // printf("smd->domain->amplify: %d\n",  smd->domain->amplify);
00256                         // printf("(smd->domain->flags & MOD_SMOKE_HIGHRES)\n");
00257                 }
00258 
00259                 if(!smd->domain->shadow)
00260                         smd->domain->shadow = MEM_callocN(sizeof(float) * smd->domain->res[0] * smd->domain->res[1] * smd->domain->res[2], "SmokeDomainShadow");
00261 
00262                 smoke_initBlenderRNA(smd->domain->fluid, &(smd->domain->alpha), &(smd->domain->beta), &(smd->domain->time_scale), &(smd->domain->vorticity), &(smd->domain->border_collisions));
00263 
00264                 if(smd->domain->wt)     
00265                 {
00266                         smoke_initWaveletBlenderRNA(smd->domain->wt, &(smd->domain->strength));
00267                         // printf("smoke_initWaveletBlenderRNA\n");
00268                 }
00269                 return 1;
00270         }
00271         else if((smd->type & MOD_SMOKE_TYPE_FLOW) && smd->flow)
00272         {
00273                 // handle flow object here
00274                 // XXX TODO
00275 
00276                 smd->time = scene->r.cfra;
00277 
00278                 // update particle lifetime to be one frame
00279                 // smd->flow->psys->part->lifetime = scene->r.efra + 1;
00280 /*
00281                 if(!smd->flow->bvh)
00282                 {
00283                         // smd->flow->bvh = MEM_callocN(sizeof(BVHTreeFromMesh), "smoke_bvhfromfaces");
00284                         // bvhtree_from_mesh_faces(smd->flow->bvh, dm, 0.0, 2, 6);
00285 
00286                         // copy obmat
00287                         // copy_m4_m4(smd->flow->mat, ob->obmat);
00288                         // copy_m4_m4(smd->flow->mat_old, ob->obmat);
00289                 }
00290 */
00291 
00292                 return 1;
00293         }
00294         else if((smd->type & MOD_SMOKE_TYPE_COLL))
00295         {
00296                 smd->time = scene->r.cfra;
00297 
00298                 // todo: delete this when loading colls work -dg
00299                 if(!smd->coll)
00300                         smokeModifier_createType(smd);
00301 
00302                 if(!smd->coll->points)
00303                 {
00304                         // init collision points
00305                         SmokeCollSettings *scs = smd->coll;
00306 
00307                         // copy obmat
00308                         copy_m4_m4(scs->mat, ob->obmat);
00309                         copy_m4_m4(scs->mat_old, ob->obmat);
00310 
00311                         fill_scs_points(ob, dm, scs);
00312                 }
00313 
00314                 if(!smd->coll->bvhtree)
00315                 {
00316                         smd->coll->bvhtree = NULL; // bvhtree_build_from_smoke ( ob->obmat, dm->getFaceArray(dm), dm->getNumFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm), 0.0 );
00317                 }
00318                 return 1;
00319         }
00320 
00321         return 2;
00322 }
00323 
00324 static void fill_scs_points(Object *ob, DerivedMesh *dm, SmokeCollSettings *scs)
00325 {
00326         MVert *mvert = dm->getVertArray(dm);
00327         MFace *mface = dm->getFaceArray(dm);
00328         int i = 0, divs = 0;
00329         int *tridivs = NULL;
00330         float cell_len = 1.0 / 50.0; // for res = 50
00331         int newdivs = 0;
00332         int quads = 0, facecounter = 0;
00333 
00334         // count quads
00335         for(i = 0; i < dm->getNumFaces(dm); i++)
00336         {
00337                 if(mface[i].v4)
00338                         quads++;
00339         }
00340 
00341         calcTriangleDivs(ob, mvert, dm->getNumVerts(dm), mface,  dm->getNumFaces(dm), dm->getNumFaces(dm) + quads, &tridivs, cell_len);
00342 
00343         // count triangle divisions
00344         for(i = 0; i < dm->getNumFaces(dm) + quads; i++)
00345         {
00346                 divs += (tridivs[3 * i] + 1) * (tridivs[3 * i + 1] + 1) * (tridivs[3 * i + 2] + 1);
00347         }
00348 
00349         // printf("divs: %d\n", divs);
00350 
00351         scs->points = MEM_callocN(sizeof(float) * (dm->getNumVerts(dm) + divs) * 3, "SmokeCollPoints");
00352 
00353         for(i = 0; i < dm->getNumVerts(dm); i++)
00354         {
00355                 float tmpvec[3];
00356                 VECCOPY(tmpvec, mvert[i].co);
00357                 mul_m4_v3(ob->obmat, tmpvec);
00358                 VECCOPY(&scs->points[i * 3], tmpvec);
00359         }
00360         
00361         for(i = 0, facecounter = 0; i < dm->getNumFaces(dm); i++)
00362         {
00363                 int again = 0;
00364                 do
00365                 {
00366                         int j, k;
00367                         int divs1 = tridivs[3 * facecounter + 0];
00368                         int divs2 = tridivs[3 * facecounter + 1];
00369                         //int divs3 = tridivs[3 * facecounter + 2];
00370                         float side1[3], side2[3], trinormorg[3], trinorm[3];
00371                         
00372                         if(again == 1 && mface[i].v4)
00373                         {
00374                                 VECSUB(side1,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
00375                                 VECSUB(side2,  mvert[ mface[i].v4 ].co, mvert[ mface[i].v1 ].co);
00376                         }
00377                         else
00378                         {
00379                                 VECSUB(side1,  mvert[ mface[i].v2 ].co, mvert[ mface[i].v1 ].co);
00380                                 VECSUB(side2,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
00381                         }
00382 
00383                         cross_v3_v3v3(trinormorg, side1, side2);
00384                         normalize_v3(trinormorg);
00385                         VECCOPY(trinorm, trinormorg);
00386                         mul_v3_fl(trinorm, 0.25 * cell_len);
00387 
00388                         for(j = 0; j <= divs1; j++)
00389                         {
00390                                 for(k = 0; k <= divs2; k++)
00391                                 {
00392                                         float p1[3], p2[3], p3[3], p[3]={0,0,0}; 
00393                                         const float uf = (float)(j + TRI_UVOFFSET) / (float)(divs1 + 0.0);
00394                                         const float vf = (float)(k + TRI_UVOFFSET) / (float)(divs2 + 0.0);
00395                                         float tmpvec[3];
00396                                         
00397                                         if(uf+vf > 1.0) 
00398                                         {
00399                                                 // printf("bigger - divs1: %d, divs2: %d\n", divs1, divs2);
00400                                                 continue;
00401                                         }
00402 
00403                                         VECCOPY(p1, mvert[ mface[i].v1 ].co);
00404                                         if(again == 1 && mface[i].v4)
00405                                         {
00406                                                 VECCOPY(p2, mvert[ mface[i].v3 ].co);
00407                                                 VECCOPY(p3, mvert[ mface[i].v4 ].co);
00408                                         }
00409                                         else
00410                                         {
00411                                                 VECCOPY(p2, mvert[ mface[i].v2 ].co);
00412                                                 VECCOPY(p3, mvert[ mface[i].v3 ].co);
00413                                         }
00414 
00415                                         mul_v3_fl(p1, (1.0-uf-vf));
00416                                         mul_v3_fl(p2, uf);
00417                                         mul_v3_fl(p3, vf);
00418                                         
00419                                         VECADD(p, p1, p2);
00420                                         VECADD(p, p, p3);
00421 
00422                                         if(newdivs > divs)
00423                                                 printf("mem problem\n");
00424 
00425                                         // mMovPoints.push_back(p + trinorm);
00426                                         VECCOPY(tmpvec, p);
00427                                         VECADD(tmpvec, tmpvec, trinorm);
00428                                         mul_m4_v3(ob->obmat, tmpvec);
00429                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
00430                                         newdivs++;
00431 
00432                                         if(newdivs > divs)
00433                                                 printf("mem problem\n");
00434 
00435                                         // mMovPoints.push_back(p - trinorm);
00436                                         VECCOPY(tmpvec, p);
00437                                         VECSUB(tmpvec, tmpvec, trinorm);
00438                                         mul_m4_v3(ob->obmat, tmpvec);
00439                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
00440                                         newdivs++;
00441                                 }
00442                         }
00443 
00444                         if(again == 0 && mface[i].v4)
00445                                 again++;
00446                         else
00447                                 again = 0;
00448 
00449                         facecounter++;
00450 
00451                 } while(again!=0);
00452         }
00453 
00454         scs->numpoints = dm->getNumVerts(dm) + newdivs;
00455 
00456         MEM_freeN(tridivs);
00457 }
00458 
00460 void calcTriangleDivs(Object *ob, MVert *verts, int UNUSED(numverts), MFace *faces, int numfaces, int numtris, int **tridivs, float cell_len) 
00461 {
00462         // mTriangleDivs1.resize( faces.size() );
00463         // mTriangleDivs2.resize( faces.size() );
00464         // mTriangleDivs3.resize( faces.size() );
00465 
00466         size_t i = 0, facecounter = 0;
00467         float maxscale[3] = {1,1,1}; // = channelFindMaxVf(mcScale);
00468         float maxpart = ABS(maxscale[0]);
00469         float scaleFac = 0;
00470         float fsTri = 0;
00471         if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
00472         if(ABS(maxscale[2])>maxpart) maxpart = ABS(maxscale[2]);
00473         scaleFac = 1.0 / maxpart;
00474         // featureSize = mLevel[mMaxRefine].nodeSize
00475         fsTri = cell_len * 0.5 * scaleFac;
00476 
00477         if(*tridivs)
00478                 MEM_freeN(*tridivs);
00479 
00480         *tridivs = MEM_callocN(sizeof(int) * numtris * 3, "Smoke_Tridivs");
00481 
00482         for(i = 0, facecounter = 0; i < numfaces; i++) 
00483         {
00484                 float p0[3], p1[3], p2[3];
00485                 float side1[3];
00486                 float side2[3];
00487                 float side3[3];
00488                 int divs1=0, divs2=0, divs3=0;
00489 
00490                 VECCOPY(p0, verts[faces[i].v1].co);
00491                 mul_m4_v3(ob->obmat, p0);
00492                 VECCOPY(p1, verts[faces[i].v2].co);
00493                 mul_m4_v3(ob->obmat, p1);
00494                 VECCOPY(p2, verts[faces[i].v3].co);
00495                 mul_m4_v3(ob->obmat, p2);
00496 
00497                 VECSUB(side1, p1, p0);
00498                 VECSUB(side2, p2, p0);
00499                 VECSUB(side3, p1, p2);
00500 
00501                 if(INPR(side1, side1) > fsTri*fsTri) 
00502                 { 
00503                         float tmp = normalize_v3(side1);
00504                         divs1 = (int)ceil(tmp/fsTri); 
00505                 }
00506                 if(INPR(side2, side2) > fsTri*fsTri) 
00507                 { 
00508                         float tmp = normalize_v3(side2);
00509                         divs2 = (int)ceil(tmp/fsTri); 
00510                         
00511                         /*
00512                         // debug
00513                         if(i==0)
00514                                 printf("b tmp: %f, fsTri: %f, divs2: %d\n", tmp, fsTri, divs2);
00515                         */
00516                 }
00517 
00518                 (*tridivs)[3 * facecounter + 0] = divs1;
00519                 (*tridivs)[3 * facecounter + 1] = divs2;
00520                 (*tridivs)[3 * facecounter + 2] = divs3;
00521 
00522                 // TODO quad case
00523                 if(faces[i].v4)
00524                 {
00525                         divs1=0, divs2=0, divs3=0;
00526 
00527                         facecounter++;
00528                         
00529                         VECCOPY(p0, verts[faces[i].v3].co);
00530                         mul_m4_v3(ob->obmat, p0);
00531                         VECCOPY(p1, verts[faces[i].v4].co);
00532                         mul_m4_v3(ob->obmat, p1);
00533                         VECCOPY(p2, verts[faces[i].v1].co);
00534                         mul_m4_v3(ob->obmat, p2);
00535 
00536                         VECSUB(side1, p1, p0);
00537                         VECSUB(side2, p2, p0);
00538                         VECSUB(side3, p1, p2);
00539 
00540                         if(INPR(side1, side1) > fsTri*fsTri) 
00541                         { 
00542                                 float tmp = normalize_v3(side1);
00543                                 divs1 = (int)ceil(tmp/fsTri); 
00544                         }
00545                         if(INPR(side2, side2) > fsTri*fsTri) 
00546                         { 
00547                                 float tmp = normalize_v3(side2);
00548                                 divs2 = (int)ceil(tmp/fsTri); 
00549                         }
00550 
00551                         (*tridivs)[3 * facecounter + 0] = divs1;
00552                         (*tridivs)[3 * facecounter + 1] = divs2;
00553                         (*tridivs)[3 * facecounter + 2] = divs3;
00554                 }
00555                 facecounter++;
00556         }
00557 }
00558 
00559 static void smokeModifier_freeDomain(SmokeModifierData *smd)
00560 {
00561         if(smd->domain)
00562         {
00563                 if(smd->domain->shadow)
00564                                 MEM_freeN(smd->domain->shadow);
00565                         smd->domain->shadow = NULL;
00566 
00567                 if(smd->domain->fluid)
00568                         smoke_free(smd->domain->fluid);
00569 
00570                 if(smd->domain->wt)
00571                         smoke_turbulence_free(smd->domain->wt);
00572 
00573                 if(smd->domain->effector_weights)
00574                                 MEM_freeN(smd->domain->effector_weights);
00575                 smd->domain->effector_weights = NULL;
00576 
00577                 BKE_ptcache_free_list(&(smd->domain->ptcaches[0]));
00578                 smd->domain->point_cache[0] = NULL;
00579 
00580                 MEM_freeN(smd->domain);
00581                 smd->domain = NULL;
00582         }
00583 }
00584 
00585 static void smokeModifier_freeFlow(SmokeModifierData *smd)
00586 {
00587         if(smd->flow)
00588         {
00589 /*
00590                 if(smd->flow->bvh)
00591                 {
00592                         free_bvhtree_from_mesh(smd->flow->bvh);
00593                         MEM_freeN(smd->flow->bvh);
00594                 }
00595                 smd->flow->bvh = NULL;
00596 */
00597                 MEM_freeN(smd->flow);
00598                 smd->flow = NULL;
00599         }
00600 }
00601 
00602 static void smokeModifier_freeCollision(SmokeModifierData *smd)
00603 {
00604         if(smd->coll)
00605         {
00606                 if(smd->coll->points)
00607                 {
00608                         MEM_freeN(smd->coll->points);
00609                         smd->coll->points = NULL;
00610                 }
00611 
00612                 if(smd->coll->bvhtree)
00613                 {
00614                         BLI_bvhtree_free(smd->coll->bvhtree);
00615                         smd->coll->bvhtree = NULL;
00616                 }
00617 
00618                 if(smd->coll->dm)
00619                         smd->coll->dm->release(smd->coll->dm);
00620                 smd->coll->dm = NULL;
00621 
00622                 MEM_freeN(smd->coll);
00623                 smd->coll = NULL;
00624         }
00625 }
00626 
00627 void smokeModifier_reset_turbulence(struct SmokeModifierData *smd)
00628 {
00629         if(smd && smd->domain && smd->domain->wt)
00630         {
00631                 smoke_turbulence_free(smd->domain->wt);
00632                 smd->domain->wt = NULL;
00633         }
00634 }
00635 
00636 void smokeModifier_reset(struct SmokeModifierData *smd)
00637 {
00638         if(smd)
00639         {
00640                 if(smd->domain)
00641                 {
00642                         if(smd->domain->shadow)
00643                                 MEM_freeN(smd->domain->shadow);
00644                         smd->domain->shadow = NULL;
00645 
00646                         if(smd->domain->fluid)
00647                         {
00648                                 smoke_free(smd->domain->fluid);
00649                                 smd->domain->fluid = NULL;
00650                         }
00651 
00652                         smokeModifier_reset_turbulence(smd);
00653 
00654                         smd->time = -1;
00655 
00656                         // printf("reset domain end\n");
00657                 }
00658                 else if(smd->flow)
00659                 {
00660                         /*
00661                         if(smd->flow->bvh)
00662                         {
00663                                 free_bvhtree_from_mesh(smd->flow->bvh);
00664                                 MEM_freeN(smd->flow->bvh);
00665                         }
00666                         smd->flow->bvh = NULL;
00667                         */
00668                 }
00669                 else if(smd->coll)
00670                 {
00671                         if(smd->coll->points)
00672                         {
00673                                 MEM_freeN(smd->coll->points);
00674                                 smd->coll->points = NULL;
00675                         }
00676 
00677                         if(smd->coll->bvhtree)
00678                         {
00679                                 BLI_bvhtree_free(smd->coll->bvhtree);
00680                                 smd->coll->bvhtree = NULL;
00681                         }
00682 
00683                         if(smd->coll->dm)
00684                                 smd->coll->dm->release(smd->coll->dm);
00685                         smd->coll->dm = NULL;
00686 
00687                 }
00688         }
00689 }
00690 
00691 void smokeModifier_free (SmokeModifierData *smd)
00692 {
00693         if(smd)
00694         {
00695                 smokeModifier_freeDomain(smd);
00696                 smokeModifier_freeFlow(smd);
00697                 smokeModifier_freeCollision(smd);
00698         }
00699 }
00700 
00701 void smokeModifier_createType(struct SmokeModifierData *smd)
00702 {
00703         if(smd)
00704         {
00705                 if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
00706                 {
00707                         if(smd->domain)
00708                                 smokeModifier_freeDomain(smd);
00709 
00710                         smd->domain = MEM_callocN(sizeof(SmokeDomainSettings), "SmokeDomain");
00711 
00712                         smd->domain->smd = smd;
00713 
00714                         smd->domain->point_cache[0] = BKE_ptcache_add(&(smd->domain->ptcaches[0]));
00715                         smd->domain->point_cache[0]->flag |= PTCACHE_DISK_CACHE;
00716                         smd->domain->point_cache[0]->step = 1;
00717 
00718                         /* Deprecated */
00719                         smd->domain->point_cache[1] = NULL;
00720                         smd->domain->ptcaches[1].first = smd->domain->ptcaches[1].last = NULL;
00721                         /* set some standard values */
00722                         smd->domain->fluid = NULL;
00723                         smd->domain->wt = NULL;                 
00724                         smd->domain->eff_group = NULL;
00725                         smd->domain->fluid_group = NULL;
00726                         smd->domain->coll_group = NULL;
00727                         smd->domain->maxres = 32;
00728                         smd->domain->amplify = 1;                       
00729                         smd->domain->omega = 1.0;                       
00730                         smd->domain->alpha = -0.001;
00731                         smd->domain->beta = 0.1;
00732                         smd->domain->time_scale = 1.0;
00733                         smd->domain->vorticity = 2.0;
00734                         smd->domain->border_collisions = 1;             // vertically non-colliding
00735                         smd->domain->flags = MOD_SMOKE_DISSOLVE_LOG | MOD_SMOKE_HIGH_SMOOTH;
00736                         smd->domain->strength = 2.0;
00737                         smd->domain->noise = MOD_SMOKE_NOISEWAVE;
00738                         smd->domain->diss_speed = 5;
00739                         // init 3dview buffer
00740 
00741                         smd->domain->viewsettings = MOD_SMOKE_VIEW_SHOWBIG;
00742                         smd->domain->effector_weights = BKE_add_effector_weights(NULL);
00743                 }
00744                 else if(smd->type & MOD_SMOKE_TYPE_FLOW)
00745                 {
00746                         if(smd->flow)
00747                                 smokeModifier_freeFlow(smd);
00748 
00749                         smd->flow = MEM_callocN(sizeof(SmokeFlowSettings), "SmokeFlow");
00750 
00751                         smd->flow->smd = smd;
00752 
00753                         /* set some standard values */
00754                         smd->flow->density = 1.0;
00755                         smd->flow->temp = 1.0;
00756                         smd->flow->flags = MOD_SMOKE_FLOW_ABSOLUTE;
00757                         smd->flow->vel_multi = 1.0;
00758 
00759                         smd->flow->psys = NULL;
00760 
00761                 }
00762                 else if(smd->type & MOD_SMOKE_TYPE_COLL)
00763                 {
00764                         if(smd->coll)
00765                                 smokeModifier_freeCollision(smd);
00766 
00767                         smd->coll = MEM_callocN(sizeof(SmokeCollSettings), "SmokeColl");
00768 
00769                         smd->coll->smd = smd;
00770                         smd->coll->points = NULL;
00771                         smd->coll->numpoints = 0;
00772                         smd->coll->bvhtree = NULL;
00773                         smd->coll->dm = NULL;
00774                 }
00775         }
00776 }
00777 
00778 void smokeModifier_copy(struct SmokeModifierData *smd, struct SmokeModifierData *tsmd)
00779 {
00780         tsmd->type = smd->type;
00781         tsmd->time = smd->time;
00782         
00783         smokeModifier_createType(tsmd);
00784 
00785         if (tsmd->domain) {
00786                 tsmd->domain->maxres = smd->domain->maxres;
00787                 tsmd->domain->amplify = smd->domain->amplify;
00788                 tsmd->domain->omega = smd->domain->omega;
00789                 tsmd->domain->alpha = smd->domain->alpha;
00790                 tsmd->domain->beta = smd->domain->beta;
00791                 tsmd->domain->flags = smd->domain->flags;
00792                 tsmd->domain->strength = smd->domain->strength;
00793                 tsmd->domain->noise = smd->domain->noise;
00794                 tsmd->domain->diss_speed = smd->domain->diss_speed;
00795                 tsmd->domain->viewsettings = smd->domain->viewsettings;
00796                 tsmd->domain->fluid_group = smd->domain->fluid_group;
00797                 tsmd->domain->coll_group = smd->domain->coll_group;
00798                 tsmd->domain->vorticity = smd->domain->vorticity;
00799                 tsmd->domain->time_scale = smd->domain->time_scale;
00800                 tsmd->domain->border_collisions = smd->domain->border_collisions;
00801                 
00802                 MEM_freeN(tsmd->domain->effector_weights);
00803                 tsmd->domain->effector_weights = MEM_dupallocN(smd->domain->effector_weights);
00804         } else if (tsmd->flow) {
00805                 tsmd->flow->density = smd->flow->density;
00806                 tsmd->flow->temp = smd->flow->temp;
00807                 tsmd->flow->psys = smd->flow->psys;
00808                 tsmd->flow->type = smd->flow->type;
00809                 tsmd->flow->flags = smd->flow->flags;
00810                 tsmd->flow->vel_multi = smd->flow->vel_multi;
00811         } else if (tsmd->coll) {
00812                 ;
00813                 /* leave it as initialised, collision settings is mostly caches */
00814         }
00815 }
00816 
00817 
00818 // forward decleration
00819 static void smoke_calc_transparency(float *result, float *input, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb, float correct);
00820 static float calc_voxel_transp(float *result, float *input, int res[3], int *pixel, float *tRay, float correct);
00821 
00822 #ifdef WITH_SMOKE
00823 
00824 static int get_lamp(Scene *scene, float *light)
00825 {       
00826         Base *base_tmp = NULL;  
00827         int found_lamp = 0;
00828 
00829         // try to find a lamp, preferably local
00830         for(base_tmp = scene->base.first; base_tmp; base_tmp= base_tmp->next) {
00831                 if(base_tmp->object->type == OB_LAMP) {
00832                         Lamp *la = base_tmp->object->data;
00833 
00834                         if(la->type == LA_LOCAL) {
00835                                 copy_v3_v3(light, base_tmp->object->obmat[3]);
00836                                 return 1;
00837                         }
00838                         else if(!found_lamp) {
00839                                 copy_v3_v3(light, base_tmp->object->obmat[3]);
00840                                 found_lamp = 1;
00841                         }
00842                 }
00843         }
00844 
00845         return found_lamp;
00846 }
00847 
00848 static void smoke_calc_domain(Scene *scene, Object *ob, SmokeModifierData *smd)
00849 {
00850         SmokeDomainSettings *sds = smd->domain;
00851         GroupObject *go = NULL;                 
00852         Base *base = NULL;      
00853 
00854         // do collisions, needs to be done before emission, so that smoke isn't emitted inside collision cells
00855         if(1)
00856         {
00857                 Object *otherobj = NULL;
00858                 ModifierData *md = NULL;
00859 
00860                 if(sds->coll_group) // we use groups since we have 2 domains
00861                         go = sds->coll_group->gobject.first;
00862                 else
00863                         base = scene->base.first;
00864 
00865                 while(base || go)
00866                 {
00867                         otherobj = NULL;
00868                         if(sds->coll_group) 
00869                         {                                               
00870                                 if(go->ob)                                                      
00871                                         otherobj = go->ob;                                      
00872                         }                                       
00873                         else                                            
00874                                 otherobj = base->object;                                        
00875                         if(!otherobj)                                   
00876                         {                                               
00877                                 if(sds->coll_group)                                                     
00878                                         go = go->next;                                          
00879                                 else                                                    
00880                                         base= base->next;                                               
00881                                 continue;                                       
00882                         }                       
00883                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
00884                         
00885                         // check for active smoke modifier
00886                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))                                    
00887                         {
00888                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
00889 
00890                                 if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll && smd2->coll->points)
00891                                 {
00892                                         // we got nice collision object
00893                                         SmokeCollSettings *scs = smd2->coll;
00894                                         size_t i, j;
00895                                         unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid);
00896 
00897                                         for(i = 0; i < scs->numpoints; i++)
00898                                         {
00899                                                 int badcell = 0;
00900                                                 size_t index = 0;
00901                                                 int cell[3];
00902 
00903                                                 // 1. get corresponding cell
00904                                                 get_cell(smd->domain->p0, smd->domain->res, smd->domain->dx, &scs->points[3 * i], cell, 0);
00905                                         
00906                                                 // check if cell is valid (in the domain boundary)
00907                                                 for(j = 0; j < 3; j++)
00908                                                         if((cell[j] > sds->res[j] - 1) || (cell[j] < 0))
00909                                                         {
00910                                                                 badcell = 1;
00911                                                                 break;
00912                                                         }
00913                                                                                                                                 
00914                                                         if(badcell)                                                                     
00915                                                                 continue;
00916                                                 // 2. set cell values (heat, density and velocity)
00917                                                 index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
00918                                                                                                                 
00919                                                 // printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]);                                                                
00920                                                 // printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index);                                                                                                                                   
00921                                                 obstacles[index] = 1;
00922                                                 // for moving gobstacles                                                                
00923                                                 /*
00924                                                 const LbmFloat maxVelVal = 0.1666;
00925                                                 const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
00926 
00927                                                 LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); 
00928                                                 {                                                               
00929                                                 const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5;                                                                
00930                                                 USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz);                                                                 
00931                                                 if(usqr>maxusqr) {                                                                      
00932                                                 // cutoff at maxVelVal                                                                  
00933                                                 for(int jj=0; jj<3; jj++) {                                                                             
00934                                                 if(objvel[jj]>0.) objvel[jj] =  maxVelVal;                                                                              
00935                                                 if(objvel[jj]<0.) objvel[jj] = -maxVelVal;                                                                      
00936                                                 }                                                               
00937                                                 } 
00938                                                 }                                                               
00939                                                 const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) );                                                          
00940                                                 const LbmVec oldov=objvel; // debug                                                             
00941                                                 objvel = vec2L((*pNormals)[n]) *dp;                                                             
00942                                                 */
00943                                         }
00944                                 }
00945                         }
00946 
00947                         if(sds->coll_group)
00948                                 go = go->next;
00949                         else
00950                                 base= base->next;
00951                 }
00952         }
00953 
00954         // do flows and fluids
00955         if(1)                   
00956         {
00957                 Object *otherobj = NULL;                                
00958                 ModifierData *md = NULL;
00959                 if(sds->fluid_group) // we use groups since we have 2 domains
00960                         go = sds->fluid_group->gobject.first;                           
00961                 else                                    
00962                         base = scene->base.first;
00963                 while(base || go)
00964                 {                                       
00965                         otherobj = NULL;
00966                         if(sds->fluid_group) 
00967                         {
00968                                 if(go->ob)                                                      
00969                                         otherobj = go->ob;                                      
00970                         }                                       
00971                         else                                            
00972                                 otherobj = base->object;
00973                         if(!otherobj)
00974                         {
00975                                 if(sds->fluid_group)
00976                                         go = go->next;
00977                                 else
00978                                         base= base->next;
00979 
00980                                 continue;
00981                         }
00982 
00983                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
00984                         
00985                         // check for active smoke modifier
00986                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
00987                         {
00988                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
00989                                 
00990                                 // check for initialized smoke object
00991                                 if((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow)                                            
00992                                 {
00993                                         // we got nice flow object
00994                                         SmokeFlowSettings *sfs = smd2->flow;
00995                                         
00996                                         if(sfs && sfs->psys && sfs->psys->part && sfs->psys->part->type==PART_EMITTER) // is particle system selected
00997                                         {
00998                                                 ParticleSimulationData sim;
00999                                                 ParticleSystem *psys = sfs->psys;
01000                                                 int p = 0;                                                              
01001                                                 float *density = smoke_get_density(sds->fluid);                                                         
01002                                                 float *bigdensity = smoke_turbulence_get_density(sds->wt);                                                              
01003                                                 float *heat = smoke_get_heat(sds->fluid);                                                               
01004                                                 float *velocity_x = smoke_get_velocity_x(sds->fluid);                                                           
01005                                                 float *velocity_y = smoke_get_velocity_y(sds->fluid);                                                           
01006                                                 float *velocity_z = smoke_get_velocity_z(sds->fluid);                                                           
01007                                                 unsigned char *obstacle = smoke_get_obstacle(sds->fluid);                                                               
01008                                                 int bigres[3];
01009                                                 short absolute_flow = (sfs->flags & MOD_SMOKE_FLOW_ABSOLUTE);
01010                                                 short high_emission_smoothing = bigdensity ? (smd->domain->flags & MOD_SMOKE_HIGH_SMOOTH) : 0;
01011 
01012                                                 /*
01013                                                 * A temporary volume map used to store whole emissive
01014                                                 * area to be added to smoke density and interpolated
01015                                                 * for high resolution smoke.
01016                                                 */
01017                                                 float *temp_emission_map = NULL;
01018 
01019                                                 sim.scene = scene;
01020                                                 sim.ob = otherobj;
01021                                                 sim.psys = psys;
01022 
01023                                                 // initialize temp emission map
01024                                                 if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW))
01025                                                 {
01026                                                         int i;
01027                                                         temp_emission_map = MEM_callocN(sizeof(float) * sds->res[0]*sds->res[1]*sds->res[2], "SmokeTempEmission");
01028                                                         // set whole volume to 0.0f
01029                                                         for (i=0; i<sds->res[0]*sds->res[1]*sds->res[2]; i++) {
01030                                                                 temp_emission_map[i] = 0.0f;
01031                                                         }
01032                                                 }
01033                                                                                                                 
01034                                                 // mostly copied from particle code                                                             
01035                                                 for(p=0; p<psys->totpart; p++)                                                          
01036                                                 {
01037                                                         int cell[3];
01038                                                         size_t i = 0;
01039                                                         size_t index = 0;
01040                                                         int badcell = 0;
01041                                                         ParticleKey state;
01042 
01043                                                         if(psys->particles[p].flag & (PARS_NO_DISP|PARS_UNEXIST))
01044                                                                 continue;
01045 
01046                                                         state.time = smd->time;
01047 
01048                                                         if(psys_get_particle_state(&sim, p, &state, 0) == 0)
01049                                                                 continue;
01050                                                         
01051                                                         // VECCOPY(pos, pa->state.co);                                                                  
01052                                                         // mul_m4_v3(ob->imat, pos);                                                                                                                                            
01053                                                         // 1. get corresponding cell    
01054                                                         get_cell(smd->domain->p0, smd->domain->res, smd->domain->dx, state.co, cell, 0);                                                                                                                                        
01055                                                         // check if cell is valid (in the domain boundary)                                                                      
01056                                                         for(i = 0; i < 3; i++)                                                                  
01057                                                         {                                                                               
01058                                                                 if((cell[i] > sds->res[i] - 1) || (cell[i] < 0))                                                                                
01059                                                                 {                                                                                       
01060                                                                         badcell = 1;                                                                                    
01061                                                                         break;                                                                          
01062                                                                 }                                                                       
01063                                                         }                                                                                                                                                       
01064                                                         if(badcell)                                                                             
01065                                                                 continue;                                                                                                                                               
01066                                                         // 2. set cell values (heat, density and velocity)                                                                      
01067                                                         index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);                                                                                                                                           
01068                                                         if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) && !(obstacle[index])) // this is inflow
01069                                                         {                                                                               
01070                                                                 // heat[index] += sfs->temp * 0.1;                                                                              
01071                                                                 // density[index] += sfs->density * 0.1;
01072                                                                 heat[index] = sfs->temp;
01073                                                                 
01074                                                                 // Add emitter density to temp emission map
01075                                                                 temp_emission_map[index] = sfs->density;
01076 
01077                                                                 // Uses particle velocity as initial velocity for smoke
01078                                                                 if(sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (psys->part->phystype != PART_PHYS_NO))
01079                                                                 {
01080                                                                         velocity_x[index] = state.vel[0]*sfs->vel_multi;
01081                                                                         velocity_y[index] = state.vel[1]*sfs->vel_multi;
01082                                                                         velocity_z[index] = state.vel[2]*sfs->vel_multi;
01083                                                                 }                                                                               
01084                                                         }                                                                       
01085                                                         else if(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) // outflow                                                                     
01086                                                         {                                                                               
01087                                                                 heat[index] = 0.f;                                                                              
01088                                                                 density[index] = 0.f;                                                                           
01089                                                                 velocity_x[index] = 0.f;                                                                                
01090                                                                 velocity_y[index] = 0.f;                                                                                
01091                                                                 velocity_z[index] = 0.f;
01092                                                                 // we need different handling for the high-res feature
01093                                                                 if(bigdensity)
01094                                                                 {
01095                                                                         // init all surrounding cells according to amplification, too                                                                                   
01096                                                                         int i, j, k;
01097                                                                         smoke_turbulence_get_res(smd->domain->wt, bigres);
01098 
01099                                                                         for(i = 0; i < smd->domain->amplify + 1; i++)
01100                                                                                 for(j = 0; j < smd->domain->amplify + 1; j++)
01101                                                                                         for(k = 0; k < smd->domain->amplify + 1; k++)
01102                                                                                         {                                                                                                               
01103                                                                                                 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);                                                                                                               
01104                                                                                                 bigdensity[index] = 0.f;                                                                                                        
01105                                                                                         }                                                                               
01106                                                                 }
01107                                                         }
01108                                                         }       // particles loop
01109 
01110 
01111                                                         // apply emission values
01112                                                         if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW)) {
01113 
01114                                                                 // initialize variables
01115                                                                 int ii, jj, kk, x, y, z, block_size;
01116                                                                 size_t index, index_big;
01117 
01118                                                                 smoke_turbulence_get_res(smd->domain->wt, bigres);
01119                                                                 block_size = smd->domain->amplify + 1;  // high res block size
01120 
01121 
01122                                                                         // loop through every low res cell
01123                                                                         for(x = 0; x < sds->res[0]; x++)
01124                                                                                 for(y = 0; y < sds->res[1]; y++)
01125                                                                                         for(z = 0; z < sds->res[2]; z++)                                                                                                        
01126                                                                                         {
01127 
01128                                                                                                 // neighbour cell emission densities (for high resolution smoke smooth interpolation)
01129                                                                                                 float c000, c001, c010, c011,  c100, c101, c110, c111;
01130 
01131                                                                                                 c000 = (x>0 && y>0 && z>0) ? temp_emission_map[smoke_get_index(x-1, sds->res[0], y-1, sds->res[1], z-1)] : 0;
01132                                                                                                 c001 = (x>0 && y>0) ? temp_emission_map[smoke_get_index(x-1, sds->res[0], y-1, sds->res[1], z)] : 0;
01133                                                                                                 c010 = (x>0 && z>0) ? temp_emission_map[smoke_get_index(x-1, sds->res[0], y, sds->res[1], z-1)] : 0;
01134                                                                                                 c011 = (x>0) ? temp_emission_map[smoke_get_index(x-1, sds->res[0], y, sds->res[1], z)] : 0;
01135 
01136                                                                                                 c100 = (y>0 && z>0) ? temp_emission_map[smoke_get_index(x, sds->res[0], y-1, sds->res[1], z-1)] : 0;
01137                                                                                                 c101 = (y>0) ? temp_emission_map[smoke_get_index(x, sds->res[0], y-1, sds->res[1], z)] : 0;
01138                                                                                                 c110 = (z>0) ? temp_emission_map[smoke_get_index(x, sds->res[0], y, sds->res[1], z-1)] : 0;
01139                                                                                                 c111 = temp_emission_map[smoke_get_index(x, sds->res[0], y, sds->res[1], z)];                   // this cell
01140 
01141 
01142 
01143                                                                                                 // get cell index
01144                                                                                                 index = smoke_get_index(x, sds->res[0], y, sds->res[1], z);
01145 
01146                                                                                                 // add emission to low resolution density
01147                                                                                                 if (absolute_flow) {if (temp_emission_map[index]>0) density[index] = temp_emission_map[index];}
01148                                                                                                 else {
01149                                                                                                         density[index] += temp_emission_map[index];
01150                                                                                                         if (density[index]>1) density[index]=1.0f;
01151                                                                                                 }
01152 
01153                                                                                                 smoke_turbulence_get_res(smd->domain->wt, bigres);
01154 
01155 
01156 
01157                                                                                                 /*
01158                                                                                                 loop through high res blocks if high res enabled
01159                                                                                                 */
01160                                                                                                 if (bigdensity)
01161                                                                                                 for(ii = 0; ii < block_size; ii++)
01162                                                                                                         for(jj = 0; jj < block_size; jj++)
01163                                                                                                                 for(kk = 0; kk < block_size; kk++)                                                                                                      
01164                                                                                                                 {
01165 
01166                                                                                                                 float fx,fy,fz, interpolated_value;
01167                                                                                                                 int shift_x, shift_y, shift_z;
01168 
01169 
01170                                                                                                                 /*
01171                                                                                                                 * Do volume interpolation if emitter smoothing
01172                                                                                                                 * is enabled
01173                                                                                                                 */
01174                                                                                                                 if (high_emission_smoothing) {
01175                                                                                                                         // convert block position to relative
01176                                                                                                                         // for interpolation smoothing
01177                                                                                                                         fx = (float)ii/block_size + 0.5f/block_size;
01178                                                                                                                         fy = (float)jj/block_size + 0.5f/block_size;
01179                                                                                                                         fz = (float)kk/block_size + 0.5f/block_size;
01180 
01181                                                                                                                         // calculate trilinear interpolation
01182                                                                                                                         interpolated_value = c000 * (1-fx) * (1-fy) * (1-fz) +
01183                                                                                                                         c100 * fx * (1-fy) * (1-fz) +
01184                                                                                                                         c010 * (1-fx) * fy * (1-fz) +
01185                                                                                                                         c001 * (1-fx) * (1-fy) * fz +
01186                                                                                                                         c101 * fx * (1-fy) * fz +
01187                                                                                                                         c011 * (1-fx) * fy * fz +
01188                                                                                                                         c110 * fx * fy * (1-fz) +
01189                                                                                                                         c111 * fx * fy * fz;
01190 
01191 
01192                                                                                                                         // add some contrast / sharpness
01193                                                                                                                         // depending on hi-res block size
01194 
01195                                                                                                                         interpolated_value = (interpolated_value-0.4f*sfs->density)*(block_size/2) + 0.4f*sfs->density;
01196                                                                                                                         if (interpolated_value<0.0f) interpolated_value = 0.0f;
01197                                                                                                                         if (interpolated_value>1.0f) interpolated_value = 1.0f;
01198 
01199                                                                                                                         // shift smoke block index
01200                                                                                                                         // (because pixel center is actually
01201                                                                                                                         // in halfway of the low res block)
01202                                                                                                                         shift_x = (x < 1) ? 0 : block_size/2;
01203                                                                                                                         shift_y = (y < 1) ? 0 : block_size/2;
01204                                                                                                                         shift_z = (z < 1) ? 0 : block_size/2;
01205                                                                                                                 }
01206                                                                                                                 else {
01207                                                                                                                         // without interpolation use same low resolution
01208                                                                                                                         // block value for all hi-res blocks
01209                                                                                                                         interpolated_value = c111;
01210                                                                                                                         shift_x = 0;
01211                                                                                                                         shift_y = 0;
01212                                                                                                                         shift_z = 0;
01213                                                                                                                 }
01214 
01215                                                                                                                 // get shifted index for current high resolution block
01216                                                                                                                 index_big = smoke_get_index(block_size * x + ii - shift_x, bigres[0], block_size * y + jj - shift_y, bigres[1], block_size * z + kk - shift_z);                                                                                                         
01217                                                                                                                 
01218                                                                                                                 // add emission data to high resolution density
01219                                                                                                                 if (absolute_flow) {if (interpolated_value > 0) bigdensity[index_big] = interpolated_value;}
01220                                                                                                                 else {
01221                                                                                                                         bigdensity[index_big] += interpolated_value;
01222                                                                                                                         if (bigdensity[index_big]>1) bigdensity[index_big]=1.0f;
01223                                                                                                                 }
01224 
01225                                                                                                                 } // end of hires loop
01226 
01227                                                                         }       // end of low res loop
01228 
01229                                                                 // free temporary emission map
01230                                                         if (temp_emission_map) MEM_freeN(temp_emission_map);
01231 
01232                                                         }       // end emission
01233 
01234 
01235                                                                                 
01236                                 }                                                       
01237                                 else                                                    
01238                                 {                                                               
01239                                         /*                                                              
01240                                         for()                                                           
01241                                         {                                                                       
01242                                                 // no psys                                                                      
01243                                                 BVHTreeNearest nearest;
01244                                                 nearest.index = -1;
01245                                                 nearest.dist = FLT_MAX;
01246 
01247                                                 BLI_bvhtree_find_nearest(sfs->bvh->tree, pco, &nearest, sfs->bvh->nearest_callback, sfs->bvh);
01248                                         }*/                                                     
01249                                 }
01250                         }                                               
01251                 }
01252                         if(sds->fluid_group)
01253                                 go = go->next;
01254                         else
01255                                 base= base->next;
01256                 }
01257         }
01258 
01259         // do effectors
01260         {
01261                 ListBase *effectors = pdInitEffectors(scene, ob, NULL, sds->effector_weights);
01262 
01263                 if(effectors)
01264                 {
01265                         float *density = smoke_get_density(sds->fluid);
01266                         float *force_x = smoke_get_force_x(sds->fluid);
01267                         float *force_y = smoke_get_force_y(sds->fluid);
01268                         float *force_z = smoke_get_force_z(sds->fluid);
01269                         float *velocity_x = smoke_get_velocity_x(sds->fluid);
01270                         float *velocity_y = smoke_get_velocity_y(sds->fluid);
01271                         float *velocity_z = smoke_get_velocity_z(sds->fluid);
01272                         int x, y, z;
01273 
01274                         // precalculate wind forces
01275                         for(x = 0; x < sds->res[0]; x++)
01276                                 for(y = 0; y < sds->res[1]; y++)
01277                                         for(z = 0; z < sds->res[2]; z++)
01278                         {       
01279                                 EffectedPoint epoint;
01280                                 float voxelCenter[3] = {0,0,0} , vel[3] = {0,0,0} , retvel[3] = {0,0,0};
01281                                 unsigned int index = smoke_get_index(x, sds->res[0], y, sds->res[1], z);
01282 
01283                                 if(density[index] < FLT_EPSILON)                                        
01284                                         continue;       
01285 
01286                                 vel[0] = velocity_x[index];
01287                                 vel[1] = velocity_y[index];
01288                                 vel[2] = velocity_z[index];
01289 
01290                                 voxelCenter[0] = sds->p0[0] + sds->dx *  x + sds->dx * 0.5;
01291                                 voxelCenter[1] = sds->p0[1] + sds->dx *  y + sds->dx * 0.5;
01292                                 voxelCenter[2] = sds->p0[2] + sds->dx *  z + sds->dx * 0.5;
01293 
01294                                 pd_point_from_loc(scene, voxelCenter, vel, index, &epoint);
01295                                 pdDoEffectors(effectors, NULL, sds->effector_weights, &epoint, retvel, NULL);
01296 
01297                                 // TODO dg - do in force!
01298                                 force_x[index] = MIN2(MAX2(-1.0, retvel[0] * 0.2), 1.0); 
01299                                 force_y[index] = MIN2(MAX2(-1.0, retvel[1] * 0.2), 1.0); 
01300                                 force_z[index] = MIN2(MAX2(-1.0, retvel[2] * 0.2), 1.0);
01301                         }
01302                 }
01303 
01304                 pdEndEffectors(&effectors);
01305         }
01306 
01307 }
01308 void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm)
01309 {       
01310         if((smd->type & MOD_SMOKE_TYPE_FLOW))
01311         {
01312                 if(scene->r.cfra >= smd->time)
01313                         smokeModifier_init(smd, ob, scene, dm);
01314 
01315                 if(scene->r.cfra > smd->time)
01316                 {
01317                         // XXX TODO
01318                         smd->time = scene->r.cfra;
01319 
01320                         // rigid movement support
01321                         /*
01322                         copy_m4_m4(smd->flow->mat_old, smd->flow->mat);
01323                         copy_m4_m4(smd->flow->mat, ob->obmat);
01324                         */
01325                 }
01326                 else if(scene->r.cfra < smd->time)
01327                 {
01328                         smd->time = scene->r.cfra;
01329                         smokeModifier_reset(smd);
01330                 }
01331         }
01332         else if(smd->type & MOD_SMOKE_TYPE_COLL)
01333         {
01334                 if(scene->r.cfra >= smd->time)
01335                         smokeModifier_init(smd, ob, scene, dm);
01336 
01337                 if(scene->r.cfra > smd->time)
01338                 {
01339                         // XXX TODO
01340                         smd->time = scene->r.cfra;
01341                         
01342                         if(smd->coll->dm)
01343                                 smd->coll->dm->release(smd->coll->dm);
01344 
01345                         smd->coll->dm = CDDM_copy(dm);
01346 
01347                         // rigid movement support
01348                         copy_m4_m4(smd->coll->mat_old, smd->coll->mat);
01349                         copy_m4_m4(smd->coll->mat, ob->obmat);
01350                 }
01351                 else if(scene->r.cfra < smd->time)
01352                 {
01353                         smd->time = scene->r.cfra;
01354                         smokeModifier_reset(smd);
01355                 }
01356         }
01357         else if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
01358         {
01359                 SmokeDomainSettings *sds = smd->domain;
01360                 float light[3]; 
01361                 PointCache *cache = NULL;
01362                 PTCacheID pid;
01363                 int startframe, endframe, framenr;
01364                 float timescale;
01365 
01366                 framenr = scene->r.cfra;
01367 
01368                 //printf("time: %d\n", scene->r.cfra);
01369 
01370                 cache = sds->point_cache[0];
01371                 BKE_ptcache_id_from_smoke(&pid, ob, smd);
01372                 BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, &timescale);
01373 
01374                 if(!smd->domain->fluid || framenr == startframe)
01375                 {
01376                         BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED);
01377                         BKE_ptcache_validate(cache, framenr);
01378                         cache->flag &= ~PTCACHE_REDO_NEEDED;
01379                 }
01380 
01381                 if(!smd->domain->fluid && (framenr != startframe) && (smd->domain->flags & MOD_SMOKE_FILE_LOAD)==0 && (cache->flag & PTCACHE_BAKED)==0)
01382                         return;
01383 
01384                 smd->domain->flags &= ~MOD_SMOKE_FILE_LOAD;
01385 
01386                 CLAMP(framenr, startframe, endframe);
01387 
01388                 /* If already viewing a pre/after frame, no need to reload */
01389                 if ((smd->time == framenr) && (framenr != scene->r.cfra))
01390                         return;
01391 
01392                 // printf("startframe: %d, framenr: %d\n", startframe, framenr);
01393 
01394                 if(smokeModifier_init(smd, ob, scene, dm)==0)
01395                 {
01396                         printf("bad smokeModifier_init\n");
01397                         return;
01398                 }
01399 
01400                 /* try to read from cache */
01401                 if(BKE_ptcache_read(&pid, (float)framenr) == PTCACHE_READ_EXACT) {
01402                         BKE_ptcache_validate(cache, framenr);
01403                         smd->time = framenr;
01404                         return;
01405                 }
01406                 
01407                 /* only calculate something when we advanced a single frame */
01408                 if(framenr != (int)smd->time+1)
01409                         return;
01410 
01411                 /* don't simulate if viewing start frame, but scene frame is not real start frame */
01412                 if (framenr != scene->r.cfra)
01413                         return;
01414 
01415                 tstart();
01416 
01417                 smoke_calc_domain(scene, ob, smd);
01418 
01419                 /* if on second frame, write cache for first frame */
01420                 if((int)smd->time == startframe && (cache->flag & PTCACHE_OUTDATED || cache->last_exact==0)) {
01421                         // create shadows straight after domain initialization so we get nice shadows for startframe, too
01422                         if(get_lamp(scene, light))
01423                                 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);
01424 
01425                         if(sds->wt)
01426                         {
01427                                 if(sds->flags & MOD_SMOKE_DISSOLVE)
01428                                         smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
01429                                 smoke_turbulence_step(sds->wt, sds->fluid);
01430                         }
01431 
01432                         BKE_ptcache_write(&pid, startframe);
01433                 }
01434                 
01435                 // set new time
01436                 smd->time = scene->r.cfra;
01437 
01438                 /* do simulation */
01439 
01440                 // low res
01441 
01442                 // simulate the actual smoke (c++ code in intern/smoke)
01443                 // DG: interesting commenting this line + deactivating loading of noise files
01444                 if(framenr!=startframe)
01445                 {
01446                         if(sds->flags & MOD_SMOKE_DISSOLVE)
01447                                 smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
01448                         smoke_step(sds->fluid, smd->time, scene->r.frs_sec / scene->r.frs_sec_base);
01449                 }
01450 
01451                 // create shadows before writing cache so they get stored
01452                 if(get_lamp(scene, light))
01453                         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);
01454 
01455                 if(sds->wt)
01456                 {
01457                         if(sds->flags & MOD_SMOKE_DISSOLVE)
01458                                 smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
01459                         smoke_turbulence_step(sds->wt, sds->fluid);
01460                 }
01461         
01462                 BKE_ptcache_validate(cache, framenr);
01463                 if(framenr != startframe)
01464                         BKE_ptcache_write(&pid, framenr);
01465 
01466                 tend();
01467                 //printf ( "Frame: %d, Time: %f\n", (int)smd->time, ( float ) tval() );
01468         }
01469 }
01470 
01471 static float calc_voxel_transp(float *result, float *input, int res[3], int *pixel, float *tRay, float correct)
01472 {
01473         const size_t index = smoke_get_index(pixel[0], res[0], pixel[1], res[1], pixel[2]);
01474 
01475         // T_ray *= T_vox
01476         *tRay *= exp(input[index]*correct);
01477         
01478         if(result[index] < 0.0f)        
01479         {
01480 #pragma omp critical            
01481                 result[index] = *tRay;  
01482         }       
01483 
01484         return *tRay;
01485 }
01486 
01487 long long smoke_get_mem_req(int xres, int yres, int zres, int amplify)
01488 {
01489         int totalCells = xres * yres * zres;
01490         int amplifiedCells = totalCells * amplify * amplify * amplify;
01491 
01492         // print out memory requirements
01493         long long int coarseSize = sizeof(float) * totalCells * 22 +
01494         sizeof(unsigned char) * totalCells;
01495 
01496         long long int fineSize = sizeof(float) * amplifiedCells * 7 + // big grids
01497         sizeof(float) * totalCells * 8 +     // small grids
01498         sizeof(float) * 128 * 128 * 128;     // noise tile
01499 
01500         long long int totalMB = (coarseSize + fineSize) / (1024 * 1024);
01501 
01502         return totalMB;
01503 }
01504 
01505 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)
01506 {
01507         int dx, dy, dz, i, l, m, n, x_inc, y_inc, z_inc, err_1, err_2, dx2, dy2, dz2;
01508         int pixel[3];
01509 
01510         pixel[0] = x1;
01511         pixel[1] = y1;
01512         pixel[2] = z1;
01513 
01514         dx = x2 - x1;
01515         dy = y2 - y1;
01516         dz = z2 - z1;
01517 
01518         x_inc = (dx < 0) ? -1 : 1;
01519         l = abs(dx);
01520         y_inc = (dy < 0) ? -1 : 1;
01521         m = abs(dy);
01522         z_inc = (dz < 0) ? -1 : 1;
01523         n = abs(dz);
01524         dx2 = l << 1;
01525         dy2 = m << 1;
01526         dz2 = n << 1;
01527 
01528         if ((l >= m) && (l >= n)) {
01529                 err_1 = dy2 - l;
01530                 err_2 = dz2 - l;
01531                 for (i = 0; i < l; i++) {
01532                         if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON)
01533                                 break;
01534                         if (err_1 > 0) {
01535                                 pixel[1] += y_inc;
01536                                 err_1 -= dx2;
01537                         }
01538                         if (err_2 > 0) {
01539                                 pixel[2] += z_inc;
01540                                 err_2 -= dx2;
01541                         }
01542                         err_1 += dy2;
01543                         err_2 += dz2;
01544                         pixel[0] += x_inc;
01545                 }
01546         } else if ((m >= l) && (m >= n)) {
01547                 err_1 = dx2 - m;
01548                 err_2 = dz2 - m;
01549                 for (i = 0; i < m; i++) {
01550                         if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON)
01551                                 break;
01552                         if (err_1 > 0) {
01553                                 pixel[0] += x_inc;
01554                                 err_1 -= dy2;
01555                         }
01556                         if (err_2 > 0) {
01557                                 pixel[2] += z_inc;
01558                                 err_2 -= dy2;
01559                         }
01560                         err_1 += dx2;
01561                         err_2 += dz2;
01562                         pixel[1] += y_inc;
01563                 }
01564         } else {
01565                 err_1 = dy2 - n;
01566                 err_2 = dx2 - n;
01567                 for (i = 0; i < n; i++) {
01568                         if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON)
01569                                 break;
01570                         if (err_1 > 0) {
01571                                 pixel[1] += y_inc;
01572                                 err_1 -= dz2;
01573                         }
01574                         if (err_2 > 0) {
01575                                 pixel[0] += x_inc;
01576                                 err_2 -= dz2;
01577                         }
01578                         err_1 += dy2;
01579                         err_2 += dx2;
01580                         pixel[2] += z_inc;
01581                 }
01582         }
01583         cb(result, input, res, pixel, tRay, correct);
01584 }
01585 
01586 static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int correct)
01587 {
01588         float tmp[3];
01589 
01590         VECSUB(tmp, pos, p0);
01591         mul_v3_fl(tmp, 1.0 / dx);
01592 
01593         if(correct)
01594         {
01595                 cell[0] = MIN2(res[0] - 1, MAX2(0, (int)floor(tmp[0])));
01596                 cell[1] = MIN2(res[1] - 1, MAX2(0, (int)floor(tmp[1])));
01597                 cell[2] = MIN2(res[2] - 1, MAX2(0, (int)floor(tmp[2])));
01598         }
01599         else
01600         {
01601                 cell[0] = (int)floor(tmp[0]);
01602                 cell[1] = (int)floor(tmp[1]);
01603                 cell[2] = (int)floor(tmp[2]);
01604         }
01605 }
01606 
01607 static void smoke_calc_transparency(float *result, float *input, float *p0, float *p1, int res[3], float dx, float *light, bresenham_callback cb, float correct)
01608 {
01609         float bv[6];
01610         int a, z, slabsize=res[0]*res[1], size= res[0]*res[1]*res[2];
01611 
01612         for(a=0; a<size; a++)
01613                 result[a]= -1.0f;
01614 
01615         bv[0] = p0[0];
01616         bv[1] = p1[0];
01617         // y
01618         bv[2] = p0[1];
01619         bv[3] = p1[1];
01620         // z
01621         bv[4] = p0[2];
01622         bv[5] = p1[2];
01623 
01624 #pragma omp parallel for schedule(static,1)
01625         for(z = 0; z < res[2]; z++)
01626         {
01627                 size_t index = z*slabsize;
01628                 int x,y;
01629 
01630                 for(y = 0; y < res[1]; y++)
01631                         for(x = 0; x < res[0]; x++, index++)
01632                         {
01633                                 float voxelCenter[3];
01634                                 float pos[3];
01635                                 int cell[3];
01636                                 float tRay = 1.0;
01637 
01638                                 if(result[index] >= 0.0f)                                       
01639                                         continue;                                                               
01640                                 voxelCenter[0] = p0[0] + dx *  x + dx * 0.5;
01641                                 voxelCenter[1] = p0[1] + dx *  y + dx * 0.5;
01642                                 voxelCenter[2] = p0[2] + dx *  z + dx * 0.5;
01643 
01644                                 // get starting position (in voxel coords)
01645                                 if(BLI_bvhtree_bb_raycast(bv, light, voxelCenter, pos) > FLT_EPSILON)
01646                                 {
01647                                         // we're ouside
01648                                         get_cell(p0, res, dx, pos, cell, 1);
01649                                 }
01650                                 else
01651                                 {
01652                                         // we're inside
01653                                         get_cell(p0, res, dx, light, cell, 1);
01654                                 }
01655 
01656                                 bresenham_linie_3D(cell[0], cell[1], cell[2], x, y, z, &tRay, cb, result, input, res, correct);
01657 
01658                                 // convention -> from a RGBA float array, use G value for tRay
01659 // #pragma omp critical
01660                                 result[index] = tRay;                   
01661                         }
01662         }
01663 }
01664 
01665 #endif // WITH_SMOKE