source/blender/blenkernel/intern/smoke.c

   1 /*
   2  * smoke.c
   3  *
   4  * $Id$
   5  *
   6  * ***** BEGIN GPL LICENSE BLOCK *****
   7  *
   8  * This program is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU General Public License
  10  * as published by the Free Software Foundation; either version 2
  11  * of the License, or (at your option) any later version.
  12  *
  13  * This program is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  * GNU General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU General Public License
  19  * along with this program; if not, write to the Free Software Foundation,
  20  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  21  *
  22  * The Original Code is Copyright (C) Blender Foundation.
  23  * All rights reserved.
  24  *
  25  * The Original Code is: all of this file.
  26  *
  27  * Contributor(s): Daniel Genrich
  28  *
  29  * ***** END GPL LICENSE BLOCK *****
  30  */
  31
  32 /** \file blender/blenkernel/intern/smoke.c
  33  *  \ingroup bke
  34  */
  35
  36
  37 /* Part of the code copied from elbeem fluid library, copyright by Nils Thuerey */
  38
  39 #include <GL/glew.h>
  40
  41 #include "MEM_guardedalloc.h"
  42
  43 #include <float.h>
  44 #include <math.h>
  45 #include <stdio.h>
  46 #include <string.h> /* memset */
  47
  48 #include "BLI_linklist.h"
  49 #include "BLI_rand.h"
  50 #include "BLI_jitter.h"
  51 #include "BLI_blenlib.h"
  52 #include "BLI_math.h"
  53 #include "BLI_edgehash.h"
  54 #include "BLI_kdtree.h"
  55 #include "BLI_kdopbvh.h"
  56 #include "BLI_utildefines.h"
  57
  58 #include "BKE_bvhutils.h"
  59 #include "BKE_cdderivedmesh.h"
  60 #include "BKE_customdata.h"
  61 #include "BKE_DerivedMesh.h"
  62 #include "BKE_effect.h"
  63 #include "BKE_modifier.h"
  64 #include "BKE_particle.h"
  65 #include "BKE_pointcache.h"
  66 #include "BKE_smoke.h"
  67
  68
  69 #include "DNA_customdata_types.h"
  70 #include "DNA_group_types.h"
  71 #include "DNA_lamp_types.h"
  72 #include "DNA_mesh_types.h"
  73 #include "DNA_meshdata_types.h"
  74 #include "DNA_modifier_types.h"
  75 #include "DNA_object_types.h"
  76 #include "DNA_particle_types.h"
  77 #include "DNA_scene_types.h"
  78 #include "DNA_smoke_types.h"
  79
  80 #include "smoke_API.h"
  81
  82 #include "BKE_smoke.h"
  83
  84 #ifdef _WIN32
  85 #include <time.h>
  86 #include <stdio.h>
  87 #include <conio.h>
  88 #include <windows.h>
  89
  90 static LARGE_INTEGER liFrequency;
  91 static LARGE_INTEGER liStartTime;
  92 static LARGE_INTEGER liCurrentTime;
  93
  94 static void tstart ( void )
  95 {
  96         QueryPerformanceFrequency ( &liFrequency );
  97         QueryPerformanceCounter ( &liStartTime );
  98 }
  99 static void tend ( void )
 100 {
 101         QueryPerformanceCounter ( &liCurrentTime );
 102 }
 103 static double tval( void )
 104 {
 105         return ((double)( (liCurrentTime.QuadPart - liStartTime.QuadPart)* (double)1000.0/(double)liFrequency.QuadPart ));
 106 }
 107 #else
 108 #include <sys/time.h>
 109 static struct timeval _tstart, _tend;
 110 static struct timezone tz;
 111 static void tstart ( void )
 112 {
 113         gettimeofday ( &_tstart, &tz );
 114 }
 115 static void tend ( void )
 116 {
 117         gettimeofday ( &_tend,&tz );
 118 }
 119
 120 #if 0 // unused
 121 static double tval()
 122 {
 123         double t1, t2;
 124         t1 = ( double ) _tstart.tv_sec*1000 + ( double ) _tstart.tv_usec/ ( 1000 );
 125         t2 = ( double ) _tend.tv_sec*1000 + ( double ) _tend.tv_usec/ ( 1000 );
 126         return t2-t1;
 127 }
 128 #endif
 129 #endif
 130
 131 struct Object;
 132 struct Scene;
 133 struct DerivedMesh;
 134 struct SmokeModifierData;
 135
 136 // forward declerations
 137 static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int correct);
 138 void calcTriangleDivs(Object *ob, MVert *verts, int numverts, MFace *tris, int numfaces, int numtris, int **tridivs, float cell_len);
 139 static void fill_scs_points(Object *ob, DerivedMesh *dm, SmokeCollSettings *scs);
 140
 141 #define TRI_UVOFFSET (1./4.)
 142
 143 /* Stubs to use when smoke is disabled */
 144 #ifndef WITH_SMOKE
 145 struct WTURBULENCE *smoke_turbulence_init(int *UNUSED(res), int UNUSED(amplify), int UNUSED(noisetype)) { return NULL; }
 146 struct FLUID_3D *smoke_init(int *UNUSED(res), float *UNUSED(p0)) { return NULL; }
 147 void smoke_free(struct FLUID_3D *UNUSED(fluid)) {}
 148 void smoke_turbulence_free(struct WTURBULENCE *UNUSED(wt)) {}
 149 void smoke_initWaveletBlenderRNA(struct WTURBULENCE *UNUSED(wt), float *UNUSED(strength)) {}
 150 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)) {}
 151 long long smoke_get_mem_req(int UNUSED(xres), int UNUSED(yres), int UNUSED(zres), int UNUSED(amplify)) { return 0; }
 152 void smokeModifier_do(SmokeModifierData *UNUSED(smd), Scene *UNUSED(scene), Object *UNUSED(ob), DerivedMesh *UNUSED(dm)) {}
 153 #endif // WITH_SMOKE
 154
 155
 156 static int smokeModifier_init (SmokeModifierData *smd, Object *ob, Scene *scene, DerivedMesh *dm)
 157 {
 158         if((smd->type & MOD_SMOKE_TYPE_DOMAIN) && smd->domain && !smd->domain->fluid)
 159         {
 160                 size_t i;
 161                 float min[3] = {FLT_MAX, FLT_MAX, FLT_MAX}, max[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
 162                 float size[3];
 163                 MVert *verts = dm->getVertArray(dm);
 164                 float scale = 0.0;
 165                 int res;
 166
 167                 res = smd->domain->maxres;
 168
 169                 // get BB of domain
 170                 for(i = 0; i < dm->getNumVerts(dm); i++)
 171                 {
 172                         float tmp[3];
 173
 174                         VECCOPY(tmp, verts[i].co);
 175                         mul_m4_v3(ob->obmat, tmp);
 176
 177                         // min BB
 178                         min[0] = MIN2(min[0], tmp[0]);
 179                         min[1] = MIN2(min[1], tmp[1]);
 180                         min[2] = MIN2(min[2], tmp[2]);
 181
 182                         // max BB
 183                         max[0] = MAX2(max[0], tmp[0]);
 184                         max[1] = MAX2(max[1], tmp[1]);
 185                         max[2] = MAX2(max[2], tmp[2]);
 186                 }
 187
 188                 VECCOPY(smd->domain->p0, min);
 189                 VECCOPY(smd->domain->p1, max);
 190
 191                 // calc other res with max_res provided
 192                 VECSUB(size, max, min);
 193
 194                 // printf("size: %f, %f, %f\n", size[0], size[1], size[2]);
 195
 196                 // prevent crash when initializing a plane as domain
 197                 if((size[0] < FLT_EPSILON) || (size[1] < FLT_EPSILON) || (size[2] < FLT_EPSILON))
 198                         return 0;
 199
 200                 if(size[0] > size[1])
 201                 {
 202                         if(size[0] > size[2])
 203                         {
 204                                 scale = res / size[0];
 205                                 smd->domain->dx = size[0] / res;
 206                                 smd->domain->res[0] = res;
 207                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
 208                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
 209                         }
 210                         else
 211                         {
 212                                 scale = res / size[2];
 213                                 smd->domain->dx = size[2] / res;
 214                                 smd->domain->res[2] = res;
 215                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
 216                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
 217                         }
 218                 }
 219                 else
 220                 {
 221                         if(size[1] > size[2])
 222                         {
 223                                 scale = res / size[1];
 224                                 smd->domain->dx = size[1] / res;
 225                                 smd->domain->res[1] = res;
 226                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
 227                                 smd->domain->res[2] = (int)(size[2] * scale + 0.5);
 228                         }
 229                         else
 230                         {
 231                                 scale = res / size[2];
 232                                 smd->domain->dx = size[2] / res;
 233                                 smd->domain->res[2] = res;
 234                                 smd->domain->res[0] = (int)(size[0] * scale + 0.5);
 235                                 smd->domain->res[1] = (int)(size[1] * scale + 0.5);
 236                         }
 237                 }
 238
 239                 // printf("smd->domain->dx: %f\n", smd->domain->dx);
 240
 241                 // TODO: put in failsafe if res<=0 - dg
 242
 243                 // printf("res[0]: %d, res[1]: %d, res[2]: %d\n", smd->domain->res[0], smd->domain->res[1], smd->domain->res[2]);
 244                 // dt max is 0.1
 245                 smd->domain->fluid = smoke_init(smd->domain->res, smd->domain->p0);
 246                 smd->time = scene->r.cfra;
 247
 248                 if(smd->domain->flags & MOD_SMOKE_HIGHRES)
 249                 {
 250                         smd->domain->wt = smoke_turbulence_init(smd->domain->res, smd->domain->amplify + 1, smd->domain->noise);
 251                         smd->domain->res_wt[0] = smd->domain->res[0] * (smd->domain->amplify + 1);
 252                         smd->domain->res_wt[1] = smd->domain->res[1] * (smd->domain->amplify + 1);
 253                         smd->domain->res_wt[2] = smd->domain->res[2] * (smd->domain->amplify + 1);
 254                         smd->domain->dx_wt = smd->domain->dx / (smd->domain->amplify + 1);
 255                         // printf("smd->domain->amplify: %d\n",  smd->domain->amplify);
 256                         // printf("(smd->domain->flags & MOD_SMOKE_HIGHRES)\n");
 257                 }
 258
 259                 if(!smd->domain->shadow)
 260                         smd->domain->shadow = MEM_callocN(sizeof(float) * smd->domain->res[0] * smd->domain->res[1] * smd->domain->res[2], "SmokeDomainShadow");
 261
 262                 smoke_initBlenderRNA(smd->domain->fluid, &(smd->domain->alpha), &(smd->domain->beta), &(smd->domain->time_scale), &(smd->domain->vorticity), &(smd->domain->border_collisions));
 263
 264                 if(smd->domain->wt)
 265                 {
 266                         smoke_initWaveletBlenderRNA(smd->domain->wt, &(smd->domain->strength));
 267                         // printf("smoke_initWaveletBlenderRNA\n");
 268                 }
 269                 return 1;
 270         }
 271         else if((smd->type & MOD_SMOKE_TYPE_FLOW) && smd->flow)
 272         {
 273                 // handle flow object here
 274                 // XXX TODO
 275
 276                 smd->time = scene->r.cfra;
 277
 278                 // update particle lifetime to be one frame
 279                 // smd->flow->psys->part->lifetime = scene->r.efra + 1;
 280 /*
 281                 if(!smd->flow->bvh)
 282                 {
 283                         // smd->flow->bvh = MEM_callocN(sizeof(BVHTreeFromMesh), "smoke_bvhfromfaces");
 284                         // bvhtree_from_mesh_faces(smd->flow->bvh, dm, 0.0, 2, 6);
 285
 286                         // copy obmat
 287                         // copy_m4_m4(smd->flow->mat, ob->obmat);
 288                         // copy_m4_m4(smd->flow->mat_old, ob->obmat);
 289                 }
 290 */
 291
 292                 return 1;
 293         }
 294         else if((smd->type & MOD_SMOKE_TYPE_COLL))
 295         {
 296                 smd->time = scene->r.cfra;
 297
 298                 // todo: delete this when loading colls work -dg
 299                 if(!smd->coll)
 300                         smokeModifier_createType(smd);
 301
 302                 if(!smd->coll->points)
 303                 {
 304                         // init collision points
 305                         SmokeCollSettings *scs = smd->coll;
 306
 307                         // copy obmat
 308                         copy_m4_m4(scs->mat, ob->obmat);
 309                         copy_m4_m4(scs->mat_old, ob->obmat);
 310
 311                         fill_scs_points(ob, dm, scs);
 312                 }
 313
 314                 if(!smd->coll->bvhtree)
 315                 {
 316                         smd->coll->bvhtree = NULL; // bvhtree_build_from_smoke ( ob->obmat, dm->getFaceArray(dm), dm->getNumFaces(dm), dm->getVertArray(dm), dm->getNumVerts(dm), 0.0 );
 317                 }
 318                 return 1;
 319         }
 320
 321         return 2;
 322 }
 323
 324 static void fill_scs_points(Object *ob, DerivedMesh *dm, SmokeCollSettings *scs)
 325 {
 326         MVert *mvert = dm->getVertArray(dm);
 327         MFace *mface = dm->getFaceArray(dm);
 328         int i = 0, divs = 0;
 329         int *tridivs = NULL;
 330         float cell_len = 1.0 / 50.0; // for res = 50
 331         int newdivs = 0;
 332         int quads = 0, facecounter = 0;
 333
 334         // count quads
 335         for(i = 0; i < dm->getNumFaces(dm); i++)
 336         {
 337                 if(mface[i].v4)
 338                         quads++;
 339         }
 340
 341         calcTriangleDivs(ob, mvert, dm->getNumVerts(dm), mface,  dm->getNumFaces(dm), dm->getNumFaces(dm) + quads, &tridivs, cell_len);
 342
 343         // count triangle divisions
 344         for(i = 0; i < dm->getNumFaces(dm) + quads; i++)
 345         {
 346                 divs += (tridivs[3 * i] + 1) * (tridivs[3 * i + 1] + 1) * (tridivs[3 * i + 2] + 1);
 347         }
 348
 349         // printf("divs: %d\n", divs);
 350
 351         scs->points = MEM_callocN(sizeof(float) * (dm->getNumVerts(dm) + divs) * 3, "SmokeCollPoints");
 352
 353         for(i = 0; i < dm->getNumVerts(dm); i++)
 354         {
 355                 float tmpvec[3];
 356                 VECCOPY(tmpvec, mvert[i].co);
 357                 mul_m4_v3(ob->obmat, tmpvec);
 358                 VECCOPY(&scs->points[i * 3], tmpvec);
 359         }
 360
 361         for(i = 0, facecounter = 0; i < dm->getNumFaces(dm); i++)
 362         {
 363                 int again = 0;
 364                 do
 365                 {
 366                         int j, k;
 367                         int divs1 = tridivs[3 * facecounter + 0];
 368                         int divs2 = tridivs[3 * facecounter + 1];
 369                         //int divs3 = tridivs[3 * facecounter + 2];
 370                         float side1[3], side2[3], trinormorg[3], trinorm[3];
 371
 372                         if(again == 1 && mface[i].v4)
 373                         {
 374                                 VECSUB(side1,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
 375                                 VECSUB(side2,  mvert[ mface[i].v4 ].co, mvert[ mface[i].v1 ].co);
 376                         }
 377                         else
 378                         {
 379                                 VECSUB(side1,  mvert[ mface[i].v2 ].co, mvert[ mface[i].v1 ].co);
 380                                 VECSUB(side2,  mvert[ mface[i].v3 ].co, mvert[ mface[i].v1 ].co);
 381                         }
 382
 383                         cross_v3_v3v3(trinormorg, side1, side2);
 384                         normalize_v3(trinormorg);
 385                         VECCOPY(trinorm, trinormorg);
 386                         mul_v3_fl(trinorm, 0.25 * cell_len);
 387
 388                         for(j = 0; j <= divs1; j++)
 389                         {
 390                                 for(k = 0; k <= divs2; k++)
 391                                 {
 392                                         float p1[3], p2[3], p3[3], p[3]={0,0,0};
 393                                         const float uf = (float)(j + TRI_UVOFFSET) / (float)(divs1 + 0.0);
 394                                         const float vf = (float)(k + TRI_UVOFFSET) / (float)(divs2 + 0.0);
 395                                         float tmpvec[3];
 396
 397                                         if(uf+vf > 1.0)
 398                                         {
 399                                                 // printf("bigger - divs1: %d, divs2: %d\n", divs1, divs2);
 400                                                 continue;
 401                                         }
 402
 403                                         VECCOPY(p1, mvert[ mface[i].v1 ].co);
 404                                         if(again == 1 && mface[i].v4)
 405                                         {
 406                                                 VECCOPY(p2, mvert[ mface[i].v3 ].co);
 407                                                 VECCOPY(p3, mvert[ mface[i].v4 ].co);
 408                                         }
 409                                         else
 410                                         {
 411                                                 VECCOPY(p2, mvert[ mface[i].v2 ].co);
 412                                                 VECCOPY(p3, mvert[ mface[i].v3 ].co);
 413                                         }
 414
 415                                         mul_v3_fl(p1, (1.0-uf-vf));
 416                                         mul_v3_fl(p2, uf);
 417                                         mul_v3_fl(p3, vf);
 418
 419                                         VECADD(p, p1, p2);
 420                                         VECADD(p, p, p3);
 421
 422                                         if(newdivs > divs)
 423                                                 printf("mem problem\n");
 424
 425                                         // mMovPoints.push_back(p + trinorm);
 426                                         VECCOPY(tmpvec, p);
 427                                         VECADD(tmpvec, tmpvec, trinorm);
 428                                         mul_m4_v3(ob->obmat, tmpvec);
 429                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
 430                                         newdivs++;
 431
 432                                         if(newdivs > divs)
 433                                                 printf("mem problem\n");
 434
 435                                         // mMovPoints.push_back(p - trinorm);
 436                                         VECCOPY(tmpvec, p);
 437                                         VECSUB(tmpvec, tmpvec, trinorm);
 438                                         mul_m4_v3(ob->obmat, tmpvec);
 439                                         VECCOPY(&scs->points[3 * (dm->getNumVerts(dm) + newdivs)], tmpvec);
 440                                         newdivs++;
 441                                 }
 442                         }
 443
 444                         if(again == 0 && mface[i].v4)
 445                                 again++;
 446                         else
 447                                 again = 0;
 448
 449                         facecounter++;
 450
 451                 } while(again!=0);
 452         }
 453
 454         scs->numpoints = dm->getNumVerts(dm) + newdivs;
 455
 456         MEM_freeN(tridivs);
 457 }
 458
 459 /*! init triangle divisions */
 460 void calcTriangleDivs(Object *ob, MVert *verts, int UNUSED(numverts), MFace *faces, int numfaces, int numtris, int **tridivs, float cell_len)
 461 {
 462         // mTriangleDivs1.resize( faces.size() );
 463         // mTriangleDivs2.resize( faces.size() );
 464         // mTriangleDivs3.resize( faces.size() );
 465
 466         size_t i = 0, facecounter = 0;
 467         float maxscale[3] = {1,1,1}; // = channelFindMaxVf(mcScale);
 468         float maxpart = ABS(maxscale[0]);
 469         float scaleFac = 0;
 470         float fsTri = 0;
 471         if(ABS(maxscale[1])>maxpart) maxpart = ABS(maxscale[1]);
 472         if(ABS(maxscale[2])>maxpart) maxpart = ABS(maxscale[2]);
 473         scaleFac = 1.0 / maxpart;
 474         // featureSize = mLevel[mMaxRefine].nodeSize
 475         fsTri = cell_len * 0.5 * scaleFac;
 476
 477         if(*tridivs)
 478                 MEM_freeN(*tridivs);
 479
 480         *tridivs = MEM_callocN(sizeof(int) * numtris * 3, "Smoke_Tridivs");
 481
 482         for(i = 0, facecounter = 0; i < numfaces; i++)
 483         {
 484                 float p0[3], p1[3], p2[3];
 485                 float side1[3];
 486                 float side2[3];
 487                 float side3[3];
 488                 int divs1=0, divs2=0, divs3=0;
 489
 490                 VECCOPY(p0, verts[faces[i].v1].co);
 491                 mul_m4_v3(ob->obmat, p0);
 492                 VECCOPY(p1, verts[faces[i].v2].co);
 493                 mul_m4_v3(ob->obmat, p1);
 494                 VECCOPY(p2, verts[faces[i].v3].co);
 495                 mul_m4_v3(ob->obmat, p2);
 496
 497                 VECSUB(side1, p1, p0);
 498                 VECSUB(side2, p2, p0);
 499                 VECSUB(side3, p1, p2);
 500
 501                 if(INPR(side1, side1) > fsTri*fsTri)
 502                 {
 503                         float tmp = normalize_v3(side1);
 504                         divs1 = (int)ceil(tmp/fsTri);
 505                 }
 506                 if(INPR(side2, side2) > fsTri*fsTri)
 507                 {
 508                         float tmp = normalize_v3(side2);
 509                         divs2 = (int)ceil(tmp/fsTri);
 510
 511                         /*
 512                         // debug
 513                         if(i==0)
 514                                 printf("b tmp: %f, fsTri: %f, divs2: %d\n", tmp, fsTri, divs2);
 515                         */
 516                 }
 517
 518                 (*tridivs)[3 * facecounter + 0] = divs1;
 519                 (*tridivs)[3 * facecounter + 1] = divs2;
 520                 (*tridivs)[3 * facecounter + 2] = divs3;
 521
 522                 // TODO quad case
 523                 if(faces[i].v4)
 524                 {
 525                         divs1=0, divs2=0, divs3=0;
 526
 527                         facecounter++;
 528
 529                         VECCOPY(p0, verts[faces[i].v3].co);
 530                         mul_m4_v3(ob->obmat, p0);
 531                         VECCOPY(p1, verts[faces[i].v4].co);
 532                         mul_m4_v3(ob->obmat, p1);
 533                         VECCOPY(p2, verts[faces[i].v1].co);
 534                         mul_m4_v3(ob->obmat, p2);
 535
 536                         VECSUB(side1, p1, p0);
 537                         VECSUB(side2, p2, p0);
 538                         VECSUB(side3, p1, p2);
 539
 540                         if(INPR(side1, side1) > fsTri*fsTri)
 541                         {
 542                                 float tmp = normalize_v3(side1);
 543                                 divs1 = (int)ceil(tmp/fsTri);
 544                         }
 545                         if(INPR(side2, side2) > fsTri*fsTri)
 546                         {
 547                                 float tmp = normalize_v3(side2);
 548                                 divs2 = (int)ceil(tmp/fsTri);
 549                         }
 550
 551                         (*tridivs)[3 * facecounter + 0] = divs1;
 552                         (*tridivs)[3 * facecounter + 1] = divs2;
 553                         (*tridivs)[3 * facecounter + 2] = divs3;
 554                 }
 555                 facecounter++;
 556         }
 557 }
 558
 559 static void smokeModifier_freeDomain(SmokeModifierData *smd)
 560 {
 561         if(smd->domain)
 562         {
 563                 if(smd->domain->shadow)
 564                                 MEM_freeN(smd->domain->shadow);
 565                         smd->domain->shadow = NULL;
 566
 567                 if(smd->domain->fluid)
 568                         smoke_free(smd->domain->fluid);
 569
 570                 if(smd->domain->wt)
 571                         smoke_turbulence_free(smd->domain->wt);
 572
 573                 if(smd->domain->effector_weights)
 574                                 MEM_freeN(smd->domain->effector_weights);
 575                 smd->domain->effector_weights = NULL;
 576
 577                 BKE_ptcache_free_list(&(smd->domain->ptcaches[0]));
 578                 smd->domain->point_cache[0] = NULL;
 579
 580                 MEM_freeN(smd->domain);
 581                 smd->domain = NULL;
 582         }
 583 }
 584
 585 static void smokeModifier_freeFlow(SmokeModifierData *smd)
 586 {
 587         if(smd->flow)
 588         {
 589 /*
 590                 if(smd->flow->bvh)
 591                 {
 592                         free_bvhtree_from_mesh(smd->flow->bvh);
 593                         MEM_freeN(smd->flow->bvh);
 594                 }
 595                 smd->flow->bvh = NULL;
 596 */
 597                 MEM_freeN(smd->flow);
 598                 smd->flow = NULL;
 599         }
 600 }
 601
 602 static void smokeModifier_freeCollision(SmokeModifierData *smd)
 603 {
 604         if(smd->coll)
 605         {
 606                 if(smd->coll->points)
 607                 {
 608                         MEM_freeN(smd->coll->points);
 609                         smd->coll->points = NULL;
 610                 }
 611
 612                 if(smd->coll->bvhtree)
 613                 {
 614                         BLI_bvhtree_free(smd->coll->bvhtree);
 615                         smd->coll->bvhtree = NULL;
 616                 }
 617
 618                 if(smd->coll->dm)
 619                         smd->coll->dm->release(smd->coll->dm);
 620                 smd->coll->dm = NULL;
 621
 622                 MEM_freeN(smd->coll);
 623                 smd->coll = NULL;
 624         }
 625 }
 626
 627 void smokeModifier_reset_turbulence(struct SmokeModifierData *smd)
 628 {
 629         if(smd && smd->domain && smd->domain->wt)
 630         {
 631                 smoke_turbulence_free(smd->domain->wt);
 632                 smd->domain->wt = NULL;
 633         }
 634 }
 635
 636 void smokeModifier_reset(struct SmokeModifierData *smd)
 637 {
 638         if(smd)
 639         {
 640                 if(smd->domain)
 641                 {
 642                         if(smd->domain->shadow)
 643                                 MEM_freeN(smd->domain->shadow);
 644                         smd->domain->shadow = NULL;
 645
 646                         if(smd->domain->fluid)
 647                         {
 648                                 smoke_free(smd->domain->fluid);
 649                                 smd->domain->fluid = NULL;
 650                         }
 651
 652                         smokeModifier_reset_turbulence(smd);
 653
 654                         smd->time = -1;
 655
 656                         // printf("reset domain end\n");
 657                 }
 658                 else if(smd->flow)
 659                 {
 660                         /*
 661                         if(smd->flow->bvh)
 662                         {
 663                                 free_bvhtree_from_mesh(smd->flow->bvh);
 664                                 MEM_freeN(smd->flow->bvh);
 665                         }
 666                         smd->flow->bvh = NULL;
 667                         */
 668                 }
 669                 else if(smd->coll)
 670                 {
 671                         if(smd->coll->points)
 672                         {
 673                                 MEM_freeN(smd->coll->points);
 674                                 smd->coll->points = NULL;
 675                         }
 676
 677                         if(smd->coll->bvhtree)
 678                         {
 679                                 BLI_bvhtree_free(smd->coll->bvhtree);
 680                                 smd->coll->bvhtree = NULL;
 681                         }
 682
 683                         if(smd->coll->dm)
 684                                 smd->coll->dm->release(smd->coll->dm);
 685                         smd->coll->dm = NULL;
 686
 687                 }
 688         }
 689 }
 690
 691 void smokeModifier_free (SmokeModifierData *smd)
 692 {
 693         if(smd)
 694         {
 695                 smokeModifier_freeDomain(smd);
 696                 smokeModifier_freeFlow(smd);
 697                 smokeModifier_freeCollision(smd);
 698         }
 699 }
 700
 701 void smokeModifier_createType(struct SmokeModifierData *smd)
 702 {
 703         if(smd)
 704         {
 705                 if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
 706                 {
 707                         if(smd->domain)
 708                                 smokeModifier_freeDomain(smd);
 709
 710                         smd->domain = MEM_callocN(sizeof(SmokeDomainSettings), "SmokeDomain");
 711
 712                         smd->domain->smd = smd;
 713
 714                         smd->domain->point_cache[0] = BKE_ptcache_add(&(smd->domain->ptcaches[0]));
 715                         smd->domain->point_cache[0]->flag |= PTCACHE_DISK_CACHE;
 716                         smd->domain->point_cache[0]->step = 1;
 717
 718                         /* Deprecated */
 719                         smd->domain->point_cache[1] = NULL;
 720                         smd->domain->ptcaches[1].first = smd->domain->ptcaches[1].last = NULL;
 721                         /* set some standard values */
 722                         smd->domain->fluid = NULL;
 723                         smd->domain->wt = NULL;
 724                         smd->domain->eff_group = NULL;
 725                         smd->domain->fluid_group = NULL;
 726                         smd->domain->coll_group = NULL;
 727                         smd->domain->maxres = 32;
 728                         smd->domain->amplify = 1;
 729                         smd->domain->omega = 1.0;
 730                         smd->domain->alpha = -0.001;
 731                         smd->domain->beta = 0.1;
 732                         smd->domain->time_scale = 1.0;
 733                         smd->domain->vorticity = 2.0;
 734                         smd->domain->border_collisions = 1;             // vertically non-colliding
 735                         smd->domain->flags = MOD_SMOKE_DISSOLVE_LOG | MOD_SMOKE_HIGH_SMOOTH;
 736                         smd->domain->strength = 2.0;
 737                         smd->domain->noise = MOD_SMOKE_NOISEWAVE;
 738                         smd->domain->diss_speed = 5;
 739                         // init 3dview buffer
 740
 741                         smd->domain->viewsettings = MOD_SMOKE_VIEW_SHOWBIG;
 742                         smd->domain->effector_weights = BKE_add_effector_weights(NULL);
 743                 }
 744                 else if(smd->type & MOD_SMOKE_TYPE_FLOW)
 745                 {
 746                         if(smd->flow)
 747                                 smokeModifier_freeFlow(smd);
 748
 749                         smd->flow = MEM_callocN(sizeof(SmokeFlowSettings), "SmokeFlow");
 750
 751                         smd->flow->smd = smd;
 752
 753                         /* set some standard values */
 754                         smd->flow->density = 1.0;
 755                         smd->flow->temp = 1.0;
 756                         smd->flow->flags = MOD_SMOKE_FLOW_ABSOLUTE;
 757                         smd->flow->vel_multi = 1.0;
 758
 759                         smd->flow->psys = NULL;
 760
 761                 }
 762                 else if(smd->type & MOD_SMOKE_TYPE_COLL)
 763                 {
 764                         if(smd->coll)
 765                                 smokeModifier_freeCollision(smd);
 766
 767                         smd->coll = MEM_callocN(sizeof(SmokeCollSettings), "SmokeColl");
 768
 769                         smd->coll->smd = smd;
 770                         smd->coll->points = NULL;
 771                         smd->coll->numpoints = 0;
 772                         smd->coll->bvhtree = NULL;
 773                         smd->coll->dm = NULL;
 774                 }
 775         }
 776 }
 777
 778 void smokeModifier_copy(struct SmokeModifierData *smd, struct SmokeModifierData *tsmd)
 779 {
 780         tsmd->type = smd->type;
 781         tsmd->time = smd->time;
 782
 783         smokeModifier_createType(tsmd);
 784
 785         if (tsmd->domain) {
 786                 tsmd->domain->maxres = smd->domain->maxres;
 787                 tsmd->domain->amplify = smd->domain->amplify;
 788                 tsmd->domain->omega = smd->domain->omega;
 789                 tsmd->domain->alpha = smd->domain->alpha;
 790                 tsmd->domain->beta = smd->domain->beta;
 791                 tsmd->domain->flags = smd->domain->flags;
 792                 tsmd->domain->strength = smd->domain->strength;
 793                 tsmd->domain->noise = smd->domain->noise;
 794                 tsmd->domain->diss_speed = smd->domain->diss_speed;
 795                 tsmd->domain->viewsettings = smd->domain->viewsettings;
 796                 tsmd->domain->fluid_group = smd->domain->fluid_group;
 797                 tsmd->domain->coll_group = smd->domain->coll_group;
 798                 tsmd->domain->vorticity = smd->domain->vorticity;
 799                 tsmd->domain->time_scale = smd->domain->time_scale;
 800                 tsmd->domain->border_collisions = smd->domain->border_collisions;
 801
 802                 MEM_freeN(tsmd->domain->effector_weights);
 803                 tsmd->domain->effector_weights = MEM_dupallocN(smd->domain->effector_weights);
 804         } else if (tsmd->flow) {
 805                 tsmd->flow->density = smd->flow->density;
 806                 tsmd->flow->temp = smd->flow->temp;
 807                 tsmd->flow->psys = smd->flow->psys;
 808                 tsmd->flow->type = smd->flow->type;
 809                 tsmd->flow->flags = smd->flow->flags;
 810                 tsmd->flow->vel_multi = smd->flow->vel_multi;
 811         } else if (tsmd->coll) {
 812                 ;
 813                 /* leave it as initialised, collision settings is mostly caches */
 814         }
 815 }
 816
 817
 818 // forward decleration
 819 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);
 820 static float calc_voxel_transp(float *result, float *input, int res[3], int *pixel, float *tRay, float correct);
 821
 822 #ifdef WITH_SMOKE
 823
 824 static int get_lamp(Scene *scene, float *light)
 825 {
 826         Base *base_tmp = NULL;
 827         int found_lamp = 0;
 828
 829         // try to find a lamp, preferably local
 830         for(base_tmp = scene->base.first; base_tmp; base_tmp= base_tmp->next) {
 831                 if(base_tmp->object->type == OB_LAMP) {
 832                         Lamp *la = base_tmp->object->data;
 833
 834                         if(la->type == LA_LOCAL) {
 835                                 copy_v3_v3(light, base_tmp->object->obmat[3]);
 836                                 return 1;
 837                         }
 838                         else if(!found_lamp) {
 839                                 copy_v3_v3(light, base_tmp->object->obmat[3]);
 840                                 found_lamp = 1;
 841                         }
 842                 }
 843         }
 844
 845         return found_lamp;
 846 }
 847
 848 static void smoke_calc_domain(Scene *scene, Object *ob, SmokeModifierData *smd)
 849 {
 850         SmokeDomainSettings *sds = smd->domain;
 851         GroupObject *go = NULL;
 852         Base *base = NULL;
 853
 854         // do collisions, needs to be done before emission, so that smoke isn't emitted inside collision cells
 855         if(1)
 856         {
 857                 Object *otherobj = NULL;
 858                 ModifierData *md = NULL;
 859
 860                 if(sds->coll_group) // we use groups since we have 2 domains
 861                         go = sds->coll_group->gobject.first;
 862                 else
 863                         base = scene->base.first;
 864
 865                 while(base || go)
 866                 {
 867                         otherobj = NULL;
 868                         if(sds->coll_group)
 869                         {
 870                                 if(go->ob)
 871                                         otherobj = go->ob;
 872                         }
 873                         else
 874                                 otherobj = base->object;
 875                         if(!otherobj)
 876                         {
 877                                 if(sds->coll_group)
 878                                         go = go->next;
 879                                 else
 880                                         base= base->next;
 881                                 continue;
 882                         }
 883                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
 884
 885                         // check for active smoke modifier
 886                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
 887                         {
 888                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
 889
 890                                 if((smd2->type & MOD_SMOKE_TYPE_COLL) && smd2->coll && smd2->coll->points)
 891                                 {
 892                                         // we got nice collision object
 893                                         SmokeCollSettings *scs = smd2->coll;
 894                                         size_t i, j;
 895                                         unsigned char *obstacles = smoke_get_obstacle(smd->domain->fluid);
 896
 897                                         for(i = 0; i < scs->numpoints; i++)
 898                                         {
 899                                                 int badcell = 0;
 900                                                 size_t index = 0;
 901                                                 int cell[3];
 902
 903                                                 // 1. get corresponding cell
 904                                                 get_cell(smd->domain->p0, smd->domain->res, smd->domain->dx, &scs->points[3 * i], cell, 0);
 905
 906                                                 // check if cell is valid (in the domain boundary)
 907                                                 for(j = 0; j < 3; j++)
 908                                                         if((cell[j] > sds->res[j] - 1) || (cell[j] < 0))
 909                                                         {
 910                                                                 badcell = 1;
 911                                                                 break;
 912                                                         }
 913
 914                                                         if(badcell)
 915                                                                 continue;
 916                                                 // 2. set cell values (heat, density and velocity)
 917                                                 index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
 918
 919                                                 // printf("cell[0]: %d, cell[1]: %d, cell[2]: %d\n", cell[0], cell[1], cell[2]);
 920                                                 // printf("res[0]: %d, res[1]: %d, res[2]: %d, index: %d\n\n", sds->res[0], sds->res[1], sds->res[2], index);
 921                                                 obstacles[index] = 1;
 922                                                 // for moving gobstacles
 923                                                 /*
 924                                                 const LbmFloat maxVelVal = 0.1666;
 925                                                 const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
 926
 927                                                 LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec);
 928                                                 {
 929                                                 const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5;
 930                                                 USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz);
 931                                                 if(usqr>maxusqr) {
 932                                                 // cutoff at maxVelVal
 933                                                 for(int jj=0; jj<3; jj++) {
 934                                                 if(objvel[jj]>0.) objvel[jj] =  maxVelVal;
 935                                                 if(objvel[jj]<0.) objvel[jj] = -maxVelVal;
 936                                                 }
 937                                                 }
 938                                                 }
 939                                                 const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) );
 940                                                 const LbmVec oldov=objvel; // debug
 941                                                 objvel = vec2L((*pNormals)[n]) *dp;
 942                                                 */
 943                                         }
 944                                 }
 945                         }
 946
 947                         if(sds->coll_group)
 948                                 go = go->next;
 949                         else
 950                                 base= base->next;
 951                 }
 952         }
 953
 954         // do flows and fluids
 955         if(1)
 956         {
 957                 Object *otherobj = NULL;
 958                 ModifierData *md = NULL;
 959                 if(sds->fluid_group) // we use groups since we have 2 domains
 960                         go = sds->fluid_group->gobject.first;
 961                 else
 962                         base = scene->base.first;
 963                 while(base || go)
 964                 {
 965                         otherobj = NULL;
 966                         if(sds->fluid_group)
 967                         {
 968                                 if(go->ob)
 969                                         otherobj = go->ob;
 970                         }
 971                         else
 972                                 otherobj = base->object;
 973                         if(!otherobj)
 974                         {
 975                                 if(sds->fluid_group)
 976                                         go = go->next;
 977                                 else
 978                                         base= base->next;
 979
 980                                 continue;
 981                         }
 982
 983                         md = modifiers_findByType(otherobj, eModifierType_Smoke);
 984
 985                         // check for active smoke modifier
 986                         if(md && md->mode & (eModifierMode_Realtime | eModifierMode_Render))
 987                         {
 988                                 SmokeModifierData *smd2 = (SmokeModifierData *)md;
 989
 990                                 // check for initialized smoke object
 991                                 if((smd2->type & MOD_SMOKE_TYPE_FLOW) && smd2->flow)
 992                                 {
 993                                         // we got nice flow object
 994                                         SmokeFlowSettings *sfs = smd2->flow;
 995
 996                                         if(sfs && sfs->psys && sfs->psys->part && sfs->psys->part->type==PART_EMITTER) // is particle system selected
 997                                         {
 998                                                 ParticleSimulationData sim;
 999                                                 ParticleSystem *psys = sfs->psys;
1000                                                 int p = 0;
1001                                                 float *density = smoke_get_density(sds->fluid);
1002                                                 float *bigdensity = smoke_turbulence_get_density(sds->wt);
1003                                                 float *heat = smoke_get_heat(sds->fluid);
1004                                                 float *velocity_x = smoke_get_velocity_x(sds->fluid);
1005                                                 float *velocity_y = smoke_get_velocity_y(sds->fluid);
1006                                                 float *velocity_z = smoke_get_velocity_z(sds->fluid);
1007                                                 unsigned char *obstacle = smoke_get_obstacle(sds->fluid);
1008                                                 int bigres[3];
1009                                                 short absolute_flow = (sfs->flags & MOD_SMOKE_FLOW_ABSOLUTE);
1010                                                 short high_emission_smoothing = bigdensity ? (smd->domain->flags & MOD_SMOKE_HIGH_SMOOTH) : 0;
1011
1012                                                 /*
1013                                                 * A temporary volume map used to store whole emissive
1014                                                 * area to be added to smoke density and interpolated
1015                                                 * for high resolution smoke.
1016                                                 */
1017                                                 float *temp_emission_map = NULL;
1018
1019                                                 sim.scene = scene;
1020                                                 sim.ob = otherobj;
1021                                                 sim.psys = psys;
1022
1023                                                 // initialize temp emission map
1024                                                 if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW))
1025                                                 {
1026                                                         int i;
1027                                                         temp_emission_map = MEM_callocN(sizeof(float) * sds->res[0]*sds->res[1]*sds->res[2], "SmokeTempEmission");
1028                                                         // set whole volume to 0.0f
1029                                                         for (i=0; i<sds->res[0]*sds->res[1]*sds->res[2]; i++) {
1030                                                                 temp_emission_map[i] = 0.0f;
1031                                                         }
1032                                                 }
1033
1034                                                 // mostly copied from particle code
1035                                                 for(p=0; p<psys->totpart; p++)
1036                                                 {
1037                                                         int cell[3];
1038                                                         size_t i = 0;
1039                                                         size_t index = 0;
1040                                                         int badcell = 0;
1041                                                         ParticleKey state;
1042
1043                                                         if(psys->particles[p].flag & (PARS_NO_DISP|PARS_UNEXIST))
1044                                                                 continue;
1045
1046                                                         state.time = smd->time;
1047
1048                                                         if(psys_get_particle_state(&sim, p, &state, 0) == 0)
1049                                                                 continue;
1050
1051                                                         // VECCOPY(pos, pa->state.co);
1052                                                         // mul_m4_v3(ob->imat, pos);
1053                                                         // 1. get corresponding cell
1054                                                         get_cell(smd->domain->p0, smd->domain->res, smd->domain->dx, state.co, cell, 0);
1055                                                         // check if cell is valid (in the domain boundary)
1056                                                         for(i = 0; i < 3; i++)
1057                                                         {
1058                                                                 if((cell[i] > sds->res[i] - 1) || (cell[i] < 0))
1059                                                                 {
1060                                                                         badcell = 1;
1061                                                                         break;
1062                                                                 }
1063                                                         }
1064                                                         if(badcell)
1065                                                                 continue;
1066                                                         // 2. set cell values (heat, density and velocity)
1067                                                         index = smoke_get_index(cell[0], sds->res[0], cell[1], sds->res[1], cell[2]);
1068                                                         if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) && !(obstacle[index])) // this is inflow
1069                                                         {
1070                                                                 // heat[index] += sfs->temp * 0.1;
1071                                                                 // density[index] += sfs->density * 0.1;
1072                                                                 heat[index] = sfs->temp;
1073
1074                                                                 // Add emitter density to temp emission map
1075                                                                 temp_emission_map[index] = sfs->density;
1076
1077                                                                 // Uses particle velocity as initial velocity for smoke
1078                                                                 if(sfs->flags & MOD_SMOKE_FLOW_INITVELOCITY && (psys->part->phystype != PART_PHYS_NO))
1079                                                                 {
1080                                                                         velocity_x[index] = state.vel[0]*sfs->vel_multi;
1081                                                                         velocity_y[index] = state.vel[1]*sfs->vel_multi;
1082                                                                         velocity_z[index] = state.vel[2]*sfs->vel_multi;
1083                                                                 }
1084                                                         }
1085                                                         else if(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW) // outflow
1086                                                         {
1087                                                                 heat[index] = 0.f;
1088                                                                 density[index] = 0.f;
1089                                                                 velocity_x[index] = 0.f;
1090                                                                 velocity_y[index] = 0.f;
1091                                                                 velocity_z[index] = 0.f;
1092                                                                 // we need different handling for the high-res feature
1093                                                                 if(bigdensity)
1094                                                                 {
1095                                                                         // init all surrounding cells according to amplification, too
1096                                                                         int i, j, k;
1097                                                                         smoke_turbulence_get_res(smd->domain->wt, bigres);
1098
1099                                                                         for(i = 0; i < smd->domain->amplify + 1; i++)
1100                                                                                 for(j = 0; j < smd->domain->amplify + 1; j++)
1101                                                                                         for(k = 0; k < smd->domain->amplify + 1; k++)
1102                                                                                         {
1103                                                                                                 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);
1104                                                                                                 bigdensity[index] = 0.f;
1105                                                                                         }
1106                                                                 }
1107                                                         }
1108                                                         }       // particles loop
1109
1110
1111                                                         // apply emission values
1112                                                         if(!(sfs->type & MOD_SMOKE_FLOW_TYPE_OUTFLOW)) {
1113
1114                                                                 // initialize variables
1115                                                                 int ii, jj, kk, x, y, z, block_size;
1116                                                                 size_t index, index_big;
1117
1118                                                                 smoke_turbulence_get_res(smd->domain->wt, bigres);
1119                                                                 block_size = smd->domain->amplify + 1;  // high res block size
1120
1121
1122                                                                         // loop through every low res cell
1123                                                                         for(x = 0; x < sds->res[0]; x++)
1124                                                                                 for(y = 0; y < sds->res[1]; y++)
1125                                                                                         for(z = 0; z < sds->res[2]; z++)
1126                                                                                         {
1127
1128                                                                                                 // neighbour cell emission densities (for high resolution smoke smooth interpolation)
1129                                                                                                 float c000, c001, c010, c011,  c100, c101, c110, c111;
1130
1131                                                                                                 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;
1132                                                                                                 c001 = (x>0 && y>0) ? temp_emission_map[smoke_get_index(x-1, sds->res[0], y-1, sds->res[1], z)] : 0;
1133                                                                                                 c010 = (x>0 && z>0) ? temp_emission_map[smoke_get_index(x-1, sds->res[0], y, sds->res[1], z-1)] : 0;
1134                                                                                                 c011 = (x>0) ? temp_emission_map[smoke_get_index(x-1, sds->res[0], y, sds->res[1], z)] : 0;
1135
1136                                                                                                 c100 = (y>0 && z>0) ? temp_emission_map[smoke_get_index(x, sds->res[0], y-1, sds->res[1], z-1)] : 0;
1137                                                                                                 c101 = (y>0) ? temp_emission_map[smoke_get_index(x, sds->res[0], y-1, sds->res[1], z)] : 0;
1138                                                                                                 c110 = (z>0) ? temp_emission_map[smoke_get_index(x, sds->res[0], y, sds->res[1], z-1)] : 0;
1139                                                                                                 c111 = temp_emission_map[smoke_get_index(x, sds->res[0], y, sds->res[1], z)];                   // this cell
1140
1141
1142
1143                                                                                                 // get cell index
1144                                                                                                 index = smoke_get_index(x, sds->res[0], y, sds->res[1], z);
1145
1146                                                                                                 // add emission to low resolution density
1147                                                                                                 if (absolute_flow) {if (temp_emission_map[index]>0) density[index] = temp_emission_map[index];}
1148                                                                                                 else {
1149                                                                                                         density[index] += temp_emission_map[index];
1150                                                                                                         if (density[index]>1) density[index]=1.0f;
1151                                                                                                 }
1152
1153                                                                                                 smoke_turbulence_get_res(smd->domain->wt, bigres);
1154
1155
1156
1157                                                                                                 /*
1158                                                                                                 loop through high res blocks if high res enabled
1159                                                                                                 */
1160                                                                                                 if (bigdensity)
1161                                                                                                 for(ii = 0; ii < block_size; ii++)
1162                                                                                                         for(jj = 0; jj < block_size; jj++)
1163                                                                                                                 for(kk = 0; kk < block_size; kk++)
1164                                                                                                                 {
1165
1166                                                                                                                 float fx,fy,fz, interpolated_value;
1167                                                                                                                 int shift_x, shift_y, shift_z;
1168
1169
1170                                                                                                                 /*
1171                                                                                                                 * Do volume interpolation if emitter smoothing
1172                                                                                                                 * is enabled
1173                                                                                                                 */
1174                                                                                                                 if (high_emission_smoothing) {
1175                                                                                                                         // convert block position to relative
1176                                                                                                                         // for interpolation smoothing
1177                                                                                                                         fx = (float)ii/block_size + 0.5f/block_size;
1178                                                                                                                         fy = (float)jj/block_size + 0.5f/block_size;
1179                                                                                                                         fz = (float)kk/block_size + 0.5f/block_size;
1180
1181                                                                                                                         // calculate trilinear interpolation
1182                                                                                                                         interpolated_value = c000 * (1-fx) * (1-fy) * (1-fz) +
1183                                                                                                                         c100 * fx * (1-fy) * (1-fz) +
1184                                                                                                                         c010 * (1-fx) * fy * (1-fz) +
1185                                                                                                                         c001 * (1-fx) * (1-fy) * fz +
1186                                                                                                                         c101 * fx * (1-fy) * fz +
1187                                                                                                                         c011 * (1-fx) * fy * fz +
1188                                                                                                                         c110 * fx * fy * (1-fz) +
1189                                                                                                                         c111 * fx * fy * fz;
1190
1191
1192                                                                                                                         // add some contrast / sharpness
1193                                                                                                                         // depending on hi-res block size
1194
1195                                                                                                                         interpolated_value = (interpolated_value-0.4f*sfs->density)*(block_size/2) + 0.4f*sfs->density;
1196                                                                                                                         if (interpolated_value<0.0f) interpolated_value = 0.0f;
1197                                                                                                                         if (interpolated_value>1.0f) interpolated_value = 1.0f;
1198
1199                                                                                                                         // shift smoke block index
1200                                                                                                                         // (because pixel center is actually
1201                                                                                                                         // in halfway of the low res block)
1202                                                                                                                         shift_x = (x < 1) ? 0 : block_size/2;
1203                                                                                                                         shift_y = (y < 1) ? 0 : block_size/2;
1204                                                                                                                         shift_z = (z < 1) ? 0 : block_size/2;
1205                                                                                                                 }
1206                                                                                                                 else {
1207                                                                                                                         // without interpolation use same low resolution
1208                                                                                                                         // block value for all hi-res blocks
1209                                                                                                                         interpolated_value = c111;
1210                                                                                                                         shift_x = 0;
1211                                                                                                                         shift_y = 0;
1212                                                                                                                         shift_z = 0;
1213                                                                                                                 }
1214
1215                                                                                                                 // get shifted index for current high resolution block
1216                                                                                                                 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);
1217
1218                                                                                                                 // add emission data to high resolution density
1219                                                                                                                 if (absolute_flow) {if (interpolated_value > 0) bigdensity[index_big] = interpolated_value;}
1220                                                                                                                 else {
1221                                                                                                                         bigdensity[index_big] += interpolated_value;
1222                                                                                                                         if (bigdensity[index_big]>1) bigdensity[index_big]=1.0f;
1223                                                                                                                 }
1224
1225                                                                                                                 } // end of hires loop
1226
1227                                                                         }       // end of low res loop
1228
1229                                                                 // free temporary emission map
1230                                                         if (temp_emission_map) MEM_freeN(temp_emission_map);
1231
1232                                                         }       // end emission
1233
1234
1235
1236                                 }
1237                                 else
1238                                 {
1239                                         /*
1240                                         for()
1241                                         {
1242                                                 // no psys
1243                                                 BVHTreeNearest nearest;
1244                                                 nearest.index = -1;
1245                                                 nearest.dist = FLT_MAX;
1246
1247                                                 BLI_bvhtree_find_nearest(sfs->bvh->tree, pco, &nearest, sfs->bvh->nearest_callback, sfs->bvh);
1248                                         }*/
1249                                 }
1250                         }
1251                 }
1252                         if(sds->fluid_group)
1253                                 go = go->next;
1254                         else
1255                                 base= base->next;
1256                 }
1257         }
1258
1259         // do effectors
1260         {
1261                 ListBase *effectors = pdInitEffectors(scene, ob, NULL, sds->effector_weights);
1262
1263                 if(effectors)
1264                 {
1265                         float *density = smoke_get_density(sds->fluid);
1266                         float *force_x = smoke_get_force_x(sds->fluid);
1267                         float *force_y = smoke_get_force_y(sds->fluid);
1268                         float *force_z = smoke_get_force_z(sds->fluid);
1269                         float *velocity_x = smoke_get_velocity_x(sds->fluid);
1270                         float *velocity_y = smoke_get_velocity_y(sds->fluid);
1271                         float *velocity_z = smoke_get_velocity_z(sds->fluid);
1272                         int x, y, z;
1273
1274                         // precalculate wind forces
1275                         for(x = 0; x < sds->res[0]; x++)
1276                                 for(y = 0; y < sds->res[1]; y++)
1277                                         for(z = 0; z < sds->res[2]; z++)
1278                         {
1279                                 EffectedPoint epoint;
1280                                 float voxelCenter[3] = {0,0,0} , vel[3] = {0,0,0} , retvel[3] = {0,0,0};
1281                                 unsigned int index = smoke_get_index(x, sds->res[0], y, sds->res[1], z);
1282
1283                                 if(density[index] < FLT_EPSILON)
1284                                         continue;
1285
1286                                 vel[0] = velocity_x[index];
1287                                 vel[1] = velocity_y[index];
1288                                 vel[2] = velocity_z[index];
1289
1290                                 voxelCenter[0] = sds->p0[0] + sds->dx *  x + sds->dx * 0.5;
1291                                 voxelCenter[1] = sds->p0[1] + sds->dx *  y + sds->dx * 0.5;
1292                                 voxelCenter[2] = sds->p0[2] + sds->dx *  z + sds->dx * 0.5;
1293
1294                                 pd_point_from_loc(scene, voxelCenter, vel, index, &epoint);
1295                                 pdDoEffectors(effectors, NULL, sds->effector_weights, &epoint, retvel, NULL);
1296
1297                                 // TODO dg - do in force!
1298                                 force_x[index] = MIN2(MAX2(-1.0, retvel[0] * 0.2), 1.0);
1299                                 force_y[index] = MIN2(MAX2(-1.0, retvel[1] * 0.2), 1.0);
1300                                 force_z[index] = MIN2(MAX2(-1.0, retvel[2] * 0.2), 1.0);
1301                         }
1302                 }
1303
1304                 pdEndEffectors(&effectors);
1305         }
1306
1307 }
1308 void smokeModifier_do(SmokeModifierData *smd, Scene *scene, Object *ob, DerivedMesh *dm)
1309 {
1310         if((smd->type & MOD_SMOKE_TYPE_FLOW))
1311         {
1312                 if(scene->r.cfra >= smd->time)
1313                         smokeModifier_init(smd, ob, scene, dm);
1314
1315                 if(scene->r.cfra > smd->time)
1316                 {
1317                         // XXX TODO
1318                         smd->time = scene->r.cfra;
1319
1320                         // rigid movement support
1321                         /*
1322                         copy_m4_m4(smd->flow->mat_old, smd->flow->mat);
1323                         copy_m4_m4(smd->flow->mat, ob->obmat);
1324                         */
1325                 }
1326                 else if(scene->r.cfra < smd->time)
1327                 {
1328                         smd->time = scene->r.cfra;
1329                         smokeModifier_reset(smd);
1330                 }
1331         }
1332         else if(smd->type & MOD_SMOKE_TYPE_COLL)
1333         {
1334                 if(scene->r.cfra >= smd->time)
1335                         smokeModifier_init(smd, ob, scene, dm);
1336
1337                 if(scene->r.cfra > smd->time)
1338                 {
1339                         // XXX TODO
1340                         smd->time = scene->r.cfra;
1341
1342                         if(smd->coll->dm)
1343                                 smd->coll->dm->release(smd->coll->dm);
1344
1345                         smd->coll->dm = CDDM_copy(dm);
1346
1347                         // rigid movement support
1348                         copy_m4_m4(smd->coll->mat_old, smd->coll->mat);
1349                         copy_m4_m4(smd->coll->mat, ob->obmat);
1350                 }
1351                 else if(scene->r.cfra < smd->time)
1352                 {
1353                         smd->time = scene->r.cfra;
1354                         smokeModifier_reset(smd);
1355                 }
1356         }
1357         else if(smd->type & MOD_SMOKE_TYPE_DOMAIN)
1358         {
1359                 SmokeDomainSettings *sds = smd->domain;
1360                 float light[3];
1361                 PointCache *cache = NULL;
1362                 PTCacheID pid;
1363                 int startframe, endframe, framenr;
1364                 float timescale;
1365
1366                 framenr = scene->r.cfra;
1367
1368                 //printf("time: %d\n", scene->r.cfra);
1369
1370                 cache = sds->point_cache[0];
1371                 BKE_ptcache_id_from_smoke(&pid, ob, smd);
1372                 BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, &timescale);
1373
1374                 if(!smd->domain->fluid || framenr == startframe)
1375                 {
1376                         BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED);
1377                         BKE_ptcache_validate(cache, framenr);
1378                         cache->flag &= ~PTCACHE_REDO_NEEDED;
1379                 }
1380
1381                 if(!smd->domain->fluid && (framenr != startframe) && (smd->domain->flags & MOD_SMOKE_FILE_LOAD)==0 && (cache->flag & PTCACHE_BAKED)==0)
1382                         return;
1383
1384                 smd->domain->flags &= ~MOD_SMOKE_FILE_LOAD;
1385
1386                 CLAMP(framenr, startframe, endframe);
1387
1388                 /* If already viewing a pre/after frame, no need to reload */
1389                 if ((smd->time == framenr) && (framenr != scene->r.cfra))
1390                         return;
1391
1392                 // printf("startframe: %d, framenr: %d\n", startframe, framenr);
1393
1394                 if(smokeModifier_init(smd, ob, scene, dm)==0)
1395                 {
1396                         printf("bad smokeModifier_init\n");
1397                         return;
1398                 }
1399
1400                 /* try to read from cache */
1401                 if(BKE_ptcache_read(&pid, (float)framenr) == PTCACHE_READ_EXACT) {
1402                         BKE_ptcache_validate(cache, framenr);
1403                         smd->time = framenr;
1404                         return;
1405                 }
1406
1407                 /* only calculate something when we advanced a single frame */
1408                 if(framenr != (int)smd->time+1)
1409                         return;
1410
1411                 /* don't simulate if viewing start frame, but scene frame is not real start frame */
1412                 if (framenr != scene->r.cfra)
1413                         return;
1414
1415                 tstart();
1416
1417                 smoke_calc_domain(scene, ob, smd);
1418
1419                 /* if on second frame, write cache for first frame */
1420                 if((int)smd->time == startframe && (cache->flag & PTCACHE_OUTDATED || cache->last_exact==0)) {
1421                         // create shadows straight after domain initialization so we get nice shadows for startframe, too
1422                         if(get_lamp(scene, light))
1423                                 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);
1424
1425                         if(sds->wt)
1426                         {
1427                                 if(sds->flags & MOD_SMOKE_DISSOLVE)
1428                                         smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
1429                                 smoke_turbulence_step(sds->wt, sds->fluid);
1430                         }
1431
1432                         BKE_ptcache_write(&pid, startframe);
1433                 }
1434
1435                 // set new time
1436                 smd->time = scene->r.cfra;
1437
1438                 /* do simulation */
1439
1440                 // low res
1441
1442                 // simulate the actual smoke (c++ code in intern/smoke)
1443                 // DG: interesting commenting this line + deactivating loading of noise files
1444                 if(framenr!=startframe)
1445                 {
1446                         if(sds->flags & MOD_SMOKE_DISSOLVE)
1447                                 smoke_dissolve(sds->fluid, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
1448                         smoke_step(sds->fluid, smd->time, scene->r.frs_sec / scene->r.frs_sec_base);
1449                 }
1450
1451                 // create shadows before writing cache so they get stored
1452                 if(get_lamp(scene, light))
1453                         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);
1454
1455                 if(sds->wt)
1456                 {
1457                         if(sds->flags & MOD_SMOKE_DISSOLVE)
1458                                 smoke_dissolve_wavelet(sds->wt, sds->diss_speed, sds->flags & MOD_SMOKE_DISSOLVE_LOG);
1459                         smoke_turbulence_step(sds->wt, sds->fluid);
1460                 }
1461
1462                 BKE_ptcache_validate(cache, framenr);
1463                 if(framenr != startframe)
1464                         BKE_ptcache_write(&pid, framenr);
1465
1466                 tend();
1467                 //printf ( "Frame: %d, Time: %f\n", (int)smd->time, ( float ) tval() );
1468         }
1469 }
1470
1471 static float calc_voxel_transp(float *result, float *input, int res[3], int *pixel, float *tRay, float correct)
1472 {
1473         const size_t index = smoke_get_index(pixel[0], res[0], pixel[1], res[1], pixel[2]);
1474
1475         // T_ray *= T_vox
1476         *tRay *= exp(input[index]*correct);
1477
1478         if(result[index] < 0.0f)
1479         {
1480 #pragma omp critical
1481                 result[index] = *tRay;
1482         }
1483
1484         return *tRay;
1485 }
1486
1487 long long smoke_get_mem_req(int xres, int yres, int zres, int amplify)
1488 {
1489         int totalCells = xres * yres * zres;
1490         int amplifiedCells = totalCells * amplify * amplify * amplify;
1491
1492         // print out memory requirements
1493         long long int coarseSize = sizeof(float) * totalCells * 22 +
1494         sizeof(unsigned char) * totalCells;
1495
1496         long long int fineSize = sizeof(float) * amplifiedCells * 7 + // big grids
1497         sizeof(float) * totalCells * 8 +     // small grids
1498         sizeof(float) * 128 * 128 * 128;     // noise tile
1499
1500         long long int totalMB = (coarseSize + fineSize) / (1024 * 1024);
1501
1502         return totalMB;
1503 }
1504
1505 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)
1506 {
1507         int dx, dy, dz, i, l, m, n, x_inc, y_inc, z_inc, err_1, err_2, dx2, dy2, dz2;
1508         int pixel[3];
1509
1510         pixel[0] = x1;
1511         pixel[1] = y1;
1512         pixel[2] = z1;
1513
1514         dx = x2 - x1;
1515         dy = y2 - y1;
1516         dz = z2 - z1;
1517
1518         x_inc = (dx < 0) ? -1 : 1;
1519         l = abs(dx);
1520         y_inc = (dy < 0) ? -1 : 1;
1521         m = abs(dy);
1522         z_inc = (dz < 0) ? -1 : 1;
1523         n = abs(dz);
1524         dx2 = l << 1;
1525         dy2 = m << 1;
1526         dz2 = n << 1;
1527
1528         if ((l >= m) && (l >= n)) {
1529                 err_1 = dy2 - l;
1530                 err_2 = dz2 - l;
1531                 for (i = 0; i < l; i++) {
1532                         if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON)
1533                                 break;
1534                         if (err_1 > 0) {
1535                                 pixel[1] += y_inc;
1536                                 err_1 -= dx2;
1537                         }
1538                         if (err_2 > 0) {
1539                                 pixel[2] += z_inc;
1540                                 err_2 -= dx2;
1541                         }
1542                         err_1 += dy2;
1543                         err_2 += dz2;
1544                         pixel[0] += x_inc;
1545                 }
1546         } else if ((m >= l) && (m >= n)) {
1547                 err_1 = dx2 - m;
1548                 err_2 = dz2 - m;
1549                 for (i = 0; i < m; i++) {
1550                         if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON)
1551                                 break;
1552                         if (err_1 > 0) {
1553                                 pixel[0] += x_inc;
1554                                 err_1 -= dy2;
1555                         }
1556                         if (err_2 > 0) {
1557                                 pixel[2] += z_inc;
1558                                 err_2 -= dy2;
1559                         }
1560                         err_1 += dx2;
1561                         err_2 += dz2;
1562                         pixel[1] += y_inc;
1563                 }
1564         } else {
1565                 err_1 = dy2 - n;
1566                 err_2 = dx2 - n;
1567                 for (i = 0; i < n; i++) {
1568                         if(cb(result, input, res, pixel, tRay, correct) <= FLT_EPSILON)
1569                                 break;
1570                         if (err_1 > 0) {
1571                                 pixel[1] += y_inc;
1572                                 err_1 -= dz2;
1573                         }
1574                         if (err_2 > 0) {
1575                                 pixel[0] += x_inc;
1576                                 err_2 -= dz2;
1577                         }
1578                         err_1 += dy2;
1579                         err_2 += dx2;
1580                         pixel[2] += z_inc;
1581                 }
1582         }
1583         cb(result, input, res, pixel, tRay, correct);
1584 }
1585
1586 static void get_cell(float *p0, int res[3], float dx, float *pos, int *cell, int correct)
1587 {
1588         float tmp[3];
1589
1590         VECSUB(tmp, pos, p0);
1591         mul_v3_fl(tmp, 1.0 / dx);
1592
1593         if(correct)
1594         {
1595                 cell[0] = MIN2(res[0] - 1, MAX2(0, (int)floor(tmp[0])));
1596                 cell[1] = MIN2(res[1] - 1, MAX2(0, (int)floor(tmp[1])));
1597                 cell[2] = MIN2(res[2] - 1, MAX2(0, (int)floor(tmp[2])));
1598         }
1599         else
1600         {
1601                 cell[0] = (int)floor(tmp[0]);
1602                 cell[1] = (int)floor(tmp[1]);
1603                 cell[2] = (int)floor(tmp[2]);
1604         }
1605 }
1606
1607 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)
1608 {
1609         float bv[6];
1610         int a, z, slabsize=res[0]*res[1], size= res[0]*res[1]*res[2];
1611
1612         for(a=0; a<size; a++)
1613                 result[a]= -1.0f;
1614
1615         bv[0] = p0[0];
1616         bv[1] = p1[0];
1617         // y
1618         bv[2] = p0[1];
1619         bv[3] = p1[1];
1620         // z
1621         bv[4] = p0[2];
1622         bv[5] = p1[2];
1623
1624 #pragma omp parallel for schedule(static,1)
1625         for(z = 0; z < res[2]; z++)
1626         {
1627                 size_t index = z*slabsize;
1628                 int x,y;
1629
1630                 for(y = 0; y < res[1]; y++)
1631                         for(x = 0; x < res[0]; x++, index++)
1632                         {
1633                                 float voxelCenter[3];
1634                                 float pos[3];
1635                                 int cell[3];
1636                                 float tRay = 1.0;
1637
1638                                 if(result[index] >= 0.0f)
1639                                         continue;
1640                                 voxelCenter[0] = p0[0] + dx *  x + dx * 0.5;
1641                                 voxelCenter[1] = p0[1] + dx *  y + dx * 0.5;
1642                                 voxelCenter[2] = p0[2] + dx *  z + dx * 0.5;
1643
1644                                 // get starting position (in voxel coords)
1645                                 if(BLI_bvhtree_bb_raycast(bv, light, voxelCenter, pos) > FLT_EPSILON)
1646                                 {
1647                                         // we're ouside
1648                                         get_cell(p0, res, dx, pos, cell, 1);
1649                                 }
1650                                 else
1651                                 {
1652                                         // we're inside
1653                                         get_cell(p0, res, dx, light, cell, 1);
1654                                 }
1655
1656                                 bresenham_linie_3D(cell[0], cell[1], cell[2], x, y, z, &tRay, cb, result, input, res, correct);
1657
1658                                 // convention -> from a RGBA float array, use G value for tRay
1659 // #pragma omp critical
1660                                 result[index] = tRay;
1661                         }
1662         }
1663 }
1664
1665 #endif // WITH_SMOKE