source/blender/blenkernel/intern/particle_system.c

   1 /* particle_system.c
   2  *
   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) 2007 by Janne Karhu.
  23  * All rights reserved.
  24  *
  25  * The Original Code is: all of this file.
  26  *
  27  * Contributor(s): Raul Fernandez Hernandez (Farsthary), Stephen Swhitehorn.
  28  *
  29  * ***** END GPL LICENSE BLOCK *****
  30  */
  31
  32 /** \file blender/blenkernel/intern/particle_system.c
  33  *  \ingroup bke
  34  */
  35
  36
  37 #include <stddef.h>
  38 #include "BLI_storage.h" /* _LARGEFILE_SOURCE */
  39
  40 #include <stdlib.h>
  41 #include <math.h>
  42 #include <string.h>
  43
  44 #include "MEM_guardedalloc.h"
  45
  46 #include "DNA_anim_types.h"
  47 #include "DNA_boid_types.h"
  48 #include "DNA_particle_types.h"
  49 #include "DNA_mesh_types.h"
  50 #include "DNA_meshdata_types.h"
  51 #include "DNA_modifier_types.h"
  52 #include "DNA_object_force.h"
  53 #include "DNA_object_types.h"
  54 #include "DNA_material_types.h"
  55 #include "DNA_curve_types.h"
  56 #include "DNA_group_types.h"
  57 #include "DNA_scene_types.h"
  58 #include "DNA_texture_types.h"
  59 #include "DNA_ipo_types.h" // XXX old animation system stuff... to be removed!
  60 #include "DNA_listBase.h"
  61
  62 #include "BLI_rand.h"
  63 #include "BLI_jitter.h"
  64 #include "BLI_math.h"
  65 #include "BLI_blenlib.h"
  66 #include "BLI_kdtree.h"
  67 #include "BLI_kdopbvh.h"
  68 #include "BLI_listbase.h"
  69 #include "BLI_threads.h"
  70 #include "BLI_storage.h" /* For _LARGEFILE64_SOURCE;  zlib needs this on some systems */
  71 #include "BLI_utildefines.h"
  72 #include "BLI_edgehash.h"
  73
  74 #include "BKE_main.h"
  75 #include "BKE_animsys.h"
  76 #include "BKE_boids.h"
  77 #include "BKE_cdderivedmesh.h"
  78 #include "BKE_collision.h"
  79 #include "BKE_displist.h"
  80 #include "BKE_effect.h"
  81 #include "BKE_particle.h"
  82 #include "BKE_global.h"
  83
  84 #include "BKE_DerivedMesh.h"
  85 #include "BKE_object.h"
  86 #include "BKE_material.h"
  87 #include "BKE_cloth.h"
  88 #include "BKE_depsgraph.h"
  89 #include "BKE_lattice.h"
  90 #include "BKE_pointcache.h"
  91 #include "BKE_mesh.h"
  92 #include "BKE_modifier.h"
  93 #include "BKE_scene.h"
  94 #include "BKE_bvhutils.h"
  95
  96 #include "PIL_time.h"
  97
  98 #include "RE_shader_ext.h"
  99
 100 /* fluid sim particle import */
 101 #ifndef DISABLE_ELBEEM
 102 #include "DNA_object_fluidsim.h"
 103 #include "LBM_fluidsim.h"
 104 #include <zlib.h>
 105 #include <string.h>
 106
 107 #ifdef WIN32
 108 #ifndef snprintf
 109 #define snprintf _snprintf
 110 #endif
 111 #endif
 112
 113 #endif // DISABLE_ELBEEM
 114
 115 /************************************************/
 116 /*                      Reacting to system events                       */
 117 /************************************************/
 118
 119 static int particles_are_dynamic(ParticleSystem *psys) {
 120         if(psys->pointcache->flag & PTCACHE_BAKED)
 121                 return 0;
 122
 123         if(psys->part->type == PART_HAIR)
 124                 return psys->flag & PSYS_HAIR_DYNAMICS;
 125         else
 126                 return ELEM3(psys->part->phystype, PART_PHYS_NEWTON, PART_PHYS_BOIDS, PART_PHYS_FLUID);
 127 }
 128
 129 static int psys_get_current_display_percentage(ParticleSystem *psys)
 130 {
 131         ParticleSettings *part=psys->part;
 132
 133         if((psys->renderdata && !particles_are_dynamic(psys)) /* non-dynamic particles can be rendered fully */
 134                 || (part->child_nbr && part->childtype) /* display percentage applies to children */
 135                 || (psys->pointcache->flag & PTCACHE_BAKING)) /* baking is always done with full amount */
 136                 return 100;
 137
 138         return psys->part->disp;
 139 }
 140
 141 static int tot_particles(ParticleSystem *psys, PTCacheID *pid)
 142 {
 143         if(pid && psys->pointcache->flag & PTCACHE_EXTERNAL)
 144                 return pid->cache->totpoint;
 145         else if(psys->part->distr == PART_DISTR_GRID && psys->part->from != PART_FROM_VERT)
 146                 return psys->part->grid_res * psys->part->grid_res * psys->part->grid_res - psys->totunexist;
 147         else
 148                 return psys->part->totpart - psys->totunexist;
 149 }
 150
 151 void psys_reset(ParticleSystem *psys, int mode)
 152 {
 153         PARTICLE_P;
 154
 155         if(ELEM(mode, PSYS_RESET_ALL, PSYS_RESET_DEPSGRAPH)) {
 156                 if(mode == PSYS_RESET_ALL || !(psys->flag & PSYS_EDITED)) {
 157                         /* don't free if not absolutely necessary */
 158                         if(psys->totpart != tot_particles(psys, NULL)) {
 159                                 psys_free_particles(psys);
 160                                 psys->totpart= 0;
 161                         }
 162
 163                         psys->totkeyed= 0;
 164                         psys->flag &= ~(PSYS_HAIR_DONE|PSYS_KEYED);
 165
 166                         if(psys->edit && psys->free_edit) {
 167                                 psys->free_edit(psys->edit);
 168                                 psys->edit = NULL;
 169                                 psys->free_edit = NULL;
 170                         }
 171                 }
 172         }
 173         else if(mode == PSYS_RESET_CACHE_MISS) {
 174                 /* set all particles to be skipped */
 175                 LOOP_PARTICLES
 176                         pa->flag |= PARS_NO_DISP;
 177         }
 178
 179         /* reset children */
 180         if(psys->child) {
 181                 MEM_freeN(psys->child);
 182                 psys->child= NULL;
 183         }
 184
 185         psys->totchild= 0;
 186
 187         /* reset path cache */
 188         psys_free_path_cache(psys, psys->edit);
 189
 190         /* reset point cache */
 191         BKE_ptcache_invalidate(psys->pointcache);
 192
 193         if(psys->fluid_springs) {
 194                 MEM_freeN(psys->fluid_springs);
 195                 psys->fluid_springs = NULL;
 196         }
 197
 198         psys->tot_fluidsprings = psys->alloc_fluidsprings = 0;
 199 }
 200
 201 static void realloc_particles(ParticleSimulationData *sim, int new_totpart)
 202 {
 203         ParticleSystem *psys = sim->psys;
 204         ParticleSettings *part = psys->part;
 205         ParticleData *newpars = NULL;
 206         BoidParticle *newboids = NULL;
 207         PARTICLE_P;
 208         int totpart, totsaved = 0;
 209
 210         if(new_totpart<0) {
 211                 if(part->distr==PART_DISTR_GRID  && part->from != PART_FROM_VERT) {
 212                         totpart= part->grid_res;
 213                         totpart*=totpart*totpart;
 214                 }
 215                 else
 216                         totpart=part->totpart;
 217         }
 218         else
 219                 totpart=new_totpart;
 220
 221         if(totpart != psys->totpart) {
 222                 if(psys->edit && psys->free_edit) {
 223                         psys->free_edit(psys->edit);
 224                         psys->edit = NULL;
 225                         psys->free_edit = NULL;
 226                 }
 227
 228                 if(totpart) {
 229                         newpars= MEM_callocN(totpart*sizeof(ParticleData), "particles");
 230                         if(newpars == NULL)
 231                                 return;
 232
 233                         if(psys->part->phystype == PART_PHYS_BOIDS) {
 234                                 newboids= MEM_callocN(totpart*sizeof(BoidParticle), "boid particles");
 235
 236                                 if(newboids == NULL) {
 237                                          /* allocation error! */
 238                                         if(newpars)
 239                                                 MEM_freeN(newpars);
 240                                         return;
 241                                 }
 242                         }
 243                 }
 244
 245                 if(psys->particles) {
 246                         totsaved=MIN2(psys->totpart,totpart);
 247                         /*save old pars*/
 248                         if(totsaved) {
 249                                 memcpy(newpars,psys->particles,totsaved*sizeof(ParticleData));
 250
 251                                 if(psys->particles->boid)
 252                                         memcpy(newboids, psys->particles->boid, totsaved*sizeof(BoidParticle));
 253                         }
 254
 255                         if(psys->particles->keys)
 256                                 MEM_freeN(psys->particles->keys);
 257
 258                         if(psys->particles->boid)
 259                                 MEM_freeN(psys->particles->boid);
 260
 261                         for(p=0, pa=newpars; p<totsaved; p++, pa++) {
 262                                 if(pa->keys) {
 263                                         pa->keys= NULL;
 264                                         pa->totkey= 0;
 265                                 }
 266                         }
 267
 268                         for(p=totsaved, pa=psys->particles+totsaved; p<psys->totpart; p++, pa++)
 269                                 if(pa->hair) MEM_freeN(pa->hair);
 270
 271                         MEM_freeN(psys->particles);
 272                         psys_free_pdd(psys);
 273                 }
 274
 275                 psys->particles=newpars;
 276                 psys->totpart=totpart;
 277
 278                 if(newboids) {
 279                         LOOP_PARTICLES
 280                                 pa->boid = newboids++;
 281                 }
 282         }
 283
 284         if(psys->child) {
 285                 MEM_freeN(psys->child);
 286                 psys->child=NULL;
 287                 psys->totchild=0;
 288         }
 289 }
 290
 291 static int get_psys_child_number(struct Scene *scene, ParticleSystem *psys)
 292 {
 293         int nbr;
 294
 295         if(!psys->part->childtype)
 296                 return 0;
 297
 298         if(psys->renderdata)
 299                 nbr= psys->part->ren_child_nbr;
 300         else
 301                 nbr= psys->part->child_nbr;
 302
 303         return get_render_child_particle_number(&scene->r, nbr);
 304 }
 305
 306 static int get_psys_tot_child(struct Scene *scene, ParticleSystem *psys)
 307 {
 308         return psys->totpart*get_psys_child_number(scene, psys);
 309 }
 310
 311 static void alloc_child_particles(ParticleSystem *psys, int tot)
 312 {
 313         if(psys->child){
 314                 /* only re-allocate if we have to */
 315                 if(psys->part->childtype && psys->totchild == tot) {
 316                         memset(psys->child, 0, tot*sizeof(ChildParticle));
 317                         return;
 318                 }
 319
 320                 MEM_freeN(psys->child);
 321                 psys->child=NULL;
 322                 psys->totchild=0;
 323         }
 324
 325         if(psys->part->childtype) {
 326                 psys->totchild= tot;
 327                 if(psys->totchild)
 328                         psys->child= MEM_callocN(psys->totchild*sizeof(ChildParticle), "child_particles");
 329         }
 330 }
 331
 332 /************************************************/
 333 /*                      Distribution                                            */
 334 /************************************************/
 335
 336 void psys_calc_dmcache(Object *ob, DerivedMesh *dm, ParticleSystem *psys)
 337 {
 338         /* use for building derived mesh mapping info:
 339
 340            node: the allocated links - total derived mesh element count
 341            nodearray: the array of nodes aligned with the base mesh's elements, so
 342                                   each original elements can reference its derived elements
 343         */
 344         Mesh *me= (Mesh*)ob->data;
 345         PARTICLE_P;
 346
 347         /* CACHE LOCATIONS */
 348         if(!dm->deformedOnly) {
 349                 /* Will use later to speed up subsurf/derivedmesh */
 350                 LinkNode *node, *nodedmelem, **nodearray;
 351                 int totdmelem, totelem, i, *origindex;
 352
 353                 if(psys->part->from == PART_FROM_VERT) {
 354                         totdmelem= dm->getNumVerts(dm);
 355                         totelem= me->totvert;
 356                         origindex= dm->getVertDataArray(dm, CD_ORIGINDEX);
 357                 }
 358                 else { /* FROM_FACE/FROM_VOLUME */
 359                         totdmelem= dm->getNumTessFaces(dm);
 360                         totelem= me->totface;
 361                         origindex= dm->getTessFaceDataArray(dm, CD_ORIGINDEX);
 362                 }
 363
 364                 nodedmelem= MEM_callocN(sizeof(LinkNode)*totdmelem, "psys node elems");
 365                 nodearray= MEM_callocN(sizeof(LinkNode *)*totelem, "psys node array");
 366
 367                 for(i=0, node=nodedmelem; i<totdmelem; i++, origindex++, node++) {
 368                         node->link= SET_INT_IN_POINTER(i);
 369
 370                         if(*origindex != -1) {
 371                                 if(nodearray[*origindex]) {
 372                                         /* prepend */
 373                                         node->next = nodearray[*origindex];
 374                                         nodearray[*origindex]= node;
 375                                 }
 376                                 else
 377                                         nodearray[*origindex]= node;
 378                         }
 379                 }
 380
 381                 /* cache the verts/faces! */
 382                 LOOP_PARTICLES {
 383                         if(pa->num < 0) {
 384                                 pa->num_dmcache = -1;
 385                                 continue;
 386                         }
 387
 388                         if(psys->part->from == PART_FROM_VERT) {
 389                                 if(nodearray[pa->num])
 390                                         pa->num_dmcache= GET_INT_FROM_POINTER(nodearray[pa->num]->link);
 391                         }
 392                         else { /* FROM_FACE/FROM_VOLUME */
 393                                 /* Note that sometimes the pa->num is over the nodearray size, this is bad, maybe there is a better place to fix this,
 394                                  * but for now passing NULL is OK. every face will be searched for the particle so its slower - Campbell */
 395                                 pa->num_dmcache= psys_particle_dm_face_lookup(ob, dm, pa->num, pa->fuv, pa->num < totelem ? nodearray[pa->num] : NULL);
 396                         }
 397                 }
 398
 399                 MEM_freeN(nodearray);
 400                 MEM_freeN(nodedmelem);
 401         }
 402         else {
 403                 /* TODO PARTICLE, make the following line unnecessary, each function
 404                  * should know to use the num or num_dmcache, set the num_dmcache to
 405                  * an invalid value, just incase */
 406
 407                 LOOP_PARTICLES
 408                         pa->num_dmcache = -1;
 409         }
 410 }
 411
 412 static void distribute_simple_children(Scene *scene, Object *ob, DerivedMesh *finaldm, ParticleSystem *psys)
 413 {
 414         ChildParticle *cpa = NULL;
 415         int i, p;
 416         int child_nbr= get_psys_child_number(scene, psys);
 417         int totpart= get_psys_tot_child(scene, psys);
 418
 419         alloc_child_particles(psys, totpart);
 420
 421         cpa = psys->child;
 422         for(i=0; i<child_nbr; i++){
 423                 for(p=0; p<psys->totpart; p++,cpa++){
 424                         float length=2.0;
 425                         cpa->parent=p;
 426
 427                         /* create even spherical distribution inside unit sphere */
 428                         while(length>=1.0f){
 429                                 cpa->fuv[0]=2.0f*BLI_frand()-1.0f;
 430                                 cpa->fuv[1]=2.0f*BLI_frand()-1.0f;
 431                                 cpa->fuv[2]=2.0f*BLI_frand()-1.0f;
 432                                 length=len_v3(cpa->fuv);
 433                         }
 434
 435                         cpa->num=-1;
 436                 }
 437         }
 438         /* dmcache must be updated for parent particles if children from faces is used */
 439         psys_calc_dmcache(ob, finaldm, psys);
 440 }
 441 static void distribute_grid(DerivedMesh *dm, ParticleSystem *psys)
 442 {
 443         ParticleData *pa=NULL;
 444         float min[3], max[3], delta[3], d;
 445         MVert *mv, *mvert = dm->getVertDataArray(dm,0);
 446         int totvert=dm->getNumVerts(dm), from=psys->part->from;
 447         int i, j, k, p, res=psys->part->grid_res, size[3], axis;
 448
 449         mv=mvert;
 450
 451         /* find bounding box of dm */
 452         copy_v3_v3(min, mv->co);
 453         copy_v3_v3(max, mv->co);
 454         mv++;
 455
 456         for(i=1; i<totvert; i++, mv++){
 457                 min[0]=MIN2(min[0],mv->co[0]);
 458                 min[1]=MIN2(min[1],mv->co[1]);
 459                 min[2]=MIN2(min[2],mv->co[2]);
 460
 461                 max[0]=MAX2(max[0],mv->co[0]);
 462                 max[1]=MAX2(max[1],mv->co[1]);
 463                 max[2]=MAX2(max[2],mv->co[2]);
 464         }
 465
 466         VECSUB(delta,max,min);
 467
 468         /* determine major axis */
 469         axis = (delta[0]>=delta[1]) ? 0 : ((delta[1]>=delta[2]) ? 1 : 2);
 470
 471         d = delta[axis]/(float)res;
 472
 473         size[axis] = res;
 474         size[(axis+1)%3] = (int)ceil(delta[(axis+1)%3]/d);
 475         size[(axis+2)%3] = (int)ceil(delta[(axis+2)%3]/d);
 476
 477         /* float errors grrr.. */
 478         size[(axis+1)%3] = MIN2(size[(axis+1)%3],res);
 479         size[(axis+2)%3] = MIN2(size[(axis+2)%3],res);
 480
 481         size[0] = MAX2(size[0], 1);
 482         size[1] = MAX2(size[1], 1);
 483         size[2] = MAX2(size[2], 1);
 484
 485         /* no full offset for flat/thin objects */
 486         min[0]+= d < delta[0] ? d/2.f : delta[0]/2.f;
 487         min[1]+= d < delta[1] ? d/2.f : delta[1]/2.f;
 488         min[2]+= d < delta[2] ? d/2.f : delta[2]/2.f;
 489
 490         for(i=0,p=0,pa=psys->particles; i<res; i++){
 491                 for(j=0; j<res; j++){
 492                         for(k=0; k<res; k++,p++,pa++){
 493                                 pa->fuv[0] = min[0] + (float)i*d;
 494                                 pa->fuv[1] = min[1] + (float)j*d;
 495                                 pa->fuv[2] = min[2] + (float)k*d;
 496                                 pa->flag |= PARS_UNEXIST;
 497                                 pa->hair_index = 0; /* abused in volume calculation */
 498                         }
 499                 }
 500         }
 501
 502         /* enable particles near verts/edges/faces/inside surface */
 503         if(from==PART_FROM_VERT){
 504                 float vec[3];
 505
 506                 pa=psys->particles;
 507
 508                 min[0] -= d/2.0f;
 509                 min[1] -= d/2.0f;
 510                 min[2] -= d/2.0f;
 511
 512                 for(i=0,mv=mvert; i<totvert; i++,mv++){
 513                         sub_v3_v3v3(vec,mv->co,min);
 514                         vec[0]/=delta[0];
 515                         vec[1]/=delta[1];
 516                         vec[2]/=delta[2];
 517                         (pa     +((int)(vec[0]*(size[0]-1))*res
 518                                 +(int)(vec[1]*(size[1]-1)))*res
 519                                 +(int)(vec[2]*(size[2]-1)))->flag &= ~PARS_UNEXIST;
 520                 }
 521         }
 522         else if(ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)){
 523                 float co1[3], co2[3];
 524
 525                 MFace *mface= NULL, *mface_array;
 526                 float v1[3], v2[3], v3[3], v4[4], lambda;
 527                 int a, a1, a2, a0mul, a1mul, a2mul, totface;
 528                 int amax= from==PART_FROM_FACE ? 3 : 1;
 529
 530                 totface=dm->getNumTessFaces(dm);
 531                 mface=mface_array=dm->getTessFaceDataArray(dm,CD_MFACE);
 532
 533                 for(a=0; a<amax; a++){
 534                         if(a==0){ a0mul=res*res; a1mul=res; a2mul=1; }
 535                         else if(a==1){ a0mul=res; a1mul=1; a2mul=res*res; }
 536                         else{ a0mul=1; a1mul=res*res; a2mul=res; }
 537
 538                         for(a1=0; a1<size[(a+1)%3]; a1++){
 539                                 for(a2=0; a2<size[(a+2)%3]; a2++){
 540                                         mface= mface_array;
 541
 542                                         pa = psys->particles + a1*a1mul + a2*a2mul;
 543                                         copy_v3_v3(co1, pa->fuv);
 544                                         co1[a] -= d < delta[a] ? d/2.f : delta[a]/2.f;
 545                                         copy_v3_v3(co2, co1);
 546                                         co2[a] += delta[a] + 0.001f*d;
 547                                         co1[a] -= 0.001f*d;
 548
 549                                         /* lets intersect the faces */
 550                                         for(i=0; i<totface; i++,mface++){
 551                                                 copy_v3_v3(v1, mvert[mface->v1].co);
 552                                                 copy_v3_v3(v2, mvert[mface->v2].co);
 553                                                 copy_v3_v3(v3, mvert[mface->v3].co);
 554
 555                                                 if(isect_axial_line_tri_v3(a, co1, co2, v2, v3, v1, &lambda)){
 556                                                         if(from==PART_FROM_FACE)
 557                                                                 (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
 558                                                         else /* store number of intersections */
 559                                                                 (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
 560                                                 }
 561
 562                                                 if(mface->v4){
 563                                                         copy_v3_v3(v4, mvert[mface->v4].co);
 564
 565                                                         if(isect_axial_line_tri_v3(a, co1, co2, v4, v1, v3, &lambda)){
 566                                                                 if(from==PART_FROM_FACE)
 567                                                                         (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
 568                                                                 else
 569                                                                         (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
 570                                                         }
 571                                                 }
 572                                         }
 573
 574                                         if(from==PART_FROM_VOLUME){
 575                                                 int in=pa->hair_index%2;
 576                                                 if(in) pa->hair_index++;
 577                                                 for(i=0; i<size[0]; i++){
 578                                                         if(in || (pa+i*a0mul)->hair_index%2)
 579                                                                 (pa+i*a0mul)->flag &= ~PARS_UNEXIST;
 580                                                         /* odd intersections == in->out / out->in */
 581                                                         /* even intersections -> in stays same */
 582                                                         in=(in + (pa+i*a0mul)->hair_index) % 2;
 583                                                 }
 584                                         }
 585                                 }
 586                         }
 587                 }
 588         }
 589
 590         if(psys->part->flag & PART_GRID_HEXAGONAL) {
 591                 for(i=0,p=0,pa=psys->particles; i<res; i++){
 592                         for(j=0; j<res; j++){
 593                                 for(k=0; k<res; k++,p++,pa++){
 594                                         if(j%2)
 595                                                 pa->fuv[0] += d/2.f;
 596
 597                                         if(k%2) {
 598                                                 pa->fuv[0] += d/2.f;
 599                                                 pa->fuv[1] += d/2.f;
 600                                         }
 601                                 }
 602                         }
 603                 }
 604         }
 605
 606         if(psys->part->flag & PART_GRID_INVERT){
 607                 for(i=0; i<size[0]; i++){
 608                         for(j=0; j<size[1]; j++){
 609                                 pa=psys->particles + res*(i*res + j);
 610                                 for(k=0; k<size[2]; k++, pa++){
 611                                         pa->flag ^= PARS_UNEXIST;
 612                                 }
 613                         }
 614                 }
 615         }
 616
 617         if(psys->part->grid_rand > 0.f) {
 618                 float rfac = d * psys->part->grid_rand;
 619                 for(p=0,pa=psys->particles; p<psys->totpart; p++,pa++){
 620                         if(pa->flag & PARS_UNEXIST)
 621                                 continue;
 622
 623                         pa->fuv[0] += rfac * (PSYS_FRAND(p + 31) - 0.5f);
 624                         pa->fuv[1] += rfac * (PSYS_FRAND(p + 32) - 0.5f);
 625                         pa->fuv[2] += rfac * (PSYS_FRAND(p + 33) - 0.5f);
 626                 }
 627         }
 628 }
 629
 630 /* modified copy from rayshade.c */
 631 static void hammersley_create(float *out, int n, int seed, float amount)
 632 {
 633         RNG *rng;
 634         double p, t, offs[2];
 635         int k, kk;
 636
 637         rng = rng_new(31415926 + n + seed);
 638         offs[0]= rng_getDouble(rng) + (double)amount;
 639         offs[1]= rng_getDouble(rng) + (double)amount;
 640         rng_free(rng);
 641
 642         for (k = 0; k < n; k++) {
 643                 t = 0;
 644                 for (p = 0.5, kk = k; kk; p *= 0.5, kk >>= 1)
 645                         if (kk & 1) /* kk mod 2 = 1 */
 646                                 t += p;
 647
 648                 out[2*k + 0]= fmod((double)k/(double)n + offs[0], 1.0);
 649                 out[2*k + 1]= fmod(t + offs[1], 1.0);
 650         }
 651 }
 652
 653 /* modified copy from effect.c */
 654 static void init_mv_jit(float *jit, int num, int seed2, float amount)
 655 {
 656         RNG *rng;
 657         float *jit2, x, rad1, rad2, rad3;
 658         int i, num2;
 659
 660         if(num==0) return;
 661
 662         rad1= (float)(1.0f/sqrtf((float)num));
 663         rad2= (float)(1.0f/((float)num));
 664         rad3= (float)sqrt((float)num)/((float)num);
 665
 666         rng = rng_new(31415926 + num + seed2);
 667         x= 0;
 668                 num2 = 2 * num;
 669         for(i=0; i<num2; i+=2) {
 670
 671                 jit[i]= x + amount*rad1*(0.5f - rng_getFloat(rng));
 672                 jit[i+1]= i/(2.0f*num) + amount*rad1*(0.5f - rng_getFloat(rng));
 673
 674                 jit[i]-= (float)floor(jit[i]);
 675                 jit[i+1]-= (float)floor(jit[i+1]);
 676
 677                 x+= rad3;
 678                 x -= (float)floor(x);
 679         }
 680
 681         jit2= MEM_mallocN(12 + 2*sizeof(float)*num, "initjit");
 682
 683         for (i=0 ; i<4 ; i++) {
 684                 BLI_jitterate1(jit, jit2, num, rad1);
 685                 BLI_jitterate1(jit, jit2, num, rad1);
 686                 BLI_jitterate2(jit, jit2, num, rad2);
 687         }
 688         MEM_freeN(jit2);
 689         rng_free(rng);
 690 }
 691
 692 static void psys_uv_to_w(float u, float v, int quad, float *w)
 693 {
 694         float vert[4][3], co[3];
 695
 696         if(!quad) {
 697                 if(u+v > 1.0f)
 698                         v= 1.0f-v;
 699                 else
 700                         u= 1.0f-u;
 701         }
 702
 703         vert[0][0]= 0.0f; vert[0][1]= 0.0f; vert[0][2]= 0.0f;
 704         vert[1][0]= 1.0f; vert[1][1]= 0.0f; vert[1][2]= 0.0f;
 705         vert[2][0]= 1.0f; vert[2][1]= 1.0f; vert[2][2]= 0.0f;
 706
 707         co[0]= u;
 708         co[1]= v;
 709         co[2]= 0.0f;
 710
 711         if(quad) {
 712                 vert[3][0]= 0.0f; vert[3][1]= 1.0f; vert[3][2]= 0.0f;
 713                 interp_weights_poly_v3( w,vert, 4, co);
 714         }
 715         else {
 716                 interp_weights_poly_v3( w,vert, 3, co);
 717                 w[3]= 0.0f;
 718         }
 719 }
 720
 721 /* Find the index in "sum" array before "value" is crossed. */
 722 static int distribute_binary_search(float *sum, int n, float value)
 723 {
 724         int mid, low=0, high=n;
 725
 726         if(value == 0.f)
 727                 return 0;
 728
 729         while(low <= high) {
 730                 mid= (low + high)/2;
 731
 732                 if(sum[mid] < value && value <= sum[mid+1])
 733                         return mid;
 734
 735                 if(sum[mid] >= value)
 736                         high= mid - 1;
 737                 else if(sum[mid] < value)
 738                         low= mid + 1;
 739                 else
 740                         return mid;
 741         }
 742
 743         return low;
 744 }
 745
 746 /* the max number if calls to rng_* funcs within psys_thread_distribute_particle
 747  * be sure to keep up to date if this changes */
 748 #define PSYS_RND_DIST_SKIP 2
 749
 750 /* note: this function must be thread safe, for from == PART_FROM_CHILD */
 751 #define ONLY_WORKING_WITH_PA_VERTS 0
 752 static void distribute_threads_exec(ParticleThread *thread, ParticleData *pa, ChildParticle *cpa, int p)
 753 {
 754         ParticleThreadContext *ctx= thread->ctx;
 755         Object *ob= ctx->sim.ob;
 756         DerivedMesh *dm= ctx->dm;
 757         float *v1, *v2, *v3, *v4, nor[3], orco1[3], co1[3], co2[3], nor1[3];
 758         float cur_d, min_d, randu, randv;
 759         int from= ctx->from;
 760         int cfrom= ctx->cfrom;
 761         int distr= ctx->distr;
 762         int i, intersect, tot;
 763         int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
 764
 765         if(from == PART_FROM_VERT) {
 766                 /* TODO_PARTICLE - use original index */
 767                 pa->num= ctx->index[p];
 768                 pa->fuv[0] = 1.0f;
 769                 pa->fuv[1] = pa->fuv[2] = pa->fuv[3] = 0.0;
 770
 771 #if ONLY_WORKING_WITH_PA_VERTS
 772                 if(ctx->tree){
 773                         KDTreeNearest ptn[3];
 774                         int w, maxw;
 775
 776                         psys_particle_on_dm(ctx->dm,from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co1,0,0,0,orco1,0);
 777                         transform_mesh_orco_verts((Mesh*)ob->data, &orco1, 1, 1);
 778                         maxw = BLI_kdtree_find_n_nearest(ctx->tree,3,orco1,NULL,ptn);
 779
 780                         for(w=0; w<maxw; w++){
 781                                 pa->verts[w]=ptn->num;
 782                         }
 783                 }
 784 #endif
 785         }
 786         else if(from == PART_FROM_FACE || from == PART_FROM_VOLUME) {
 787                 MFace *mface;
 788
 789                 pa->num = i = ctx->index[p];
 790                 mface = dm->getTessFaceData(dm,i,CD_MFACE);
 791
 792                 switch(distr){
 793                 case PART_DISTR_JIT:
 794                         if(ctx->jitlevel == 1) {
 795                                 if(mface->v4)
 796                                         psys_uv_to_w(0.5f, 0.5f, mface->v4, pa->fuv);
 797                                 else
 798                                         psys_uv_to_w(0.33333f, 0.33333f, mface->v4, pa->fuv);
 799                         }
 800                         else {
 801                                 ctx->jitoff[i] = fmod(ctx->jitoff[i],(float)ctx->jitlevel);
 802                                 psys_uv_to_w(ctx->jit[2*(int)ctx->jitoff[i]], ctx->jit[2*(int)ctx->jitoff[i]+1], mface->v4, pa->fuv);
 803                                 ctx->jitoff[i]++;
 804                         }
 805                         break;
 806                 case PART_DISTR_RAND:
 807                         randu= rng_getFloat(thread->rng);
 808                         randv= rng_getFloat(thread->rng);
 809                         rng_skip_tot -= 2;
 810
 811                         psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
 812                         break;
 813                 }
 814                 pa->foffset= 0.0f;
 815
 816                 /* experimental */
 817                 if(from==PART_FROM_VOLUME){
 818                         MVert *mvert=dm->getVertDataArray(dm,CD_MVERT);
 819
 820                         tot=dm->getNumTessFaces(dm);
 821
 822                         psys_interpolate_face(mvert,mface,0,0,pa->fuv,co1,nor,0,0,0,0);
 823
 824                         normalize_v3(nor);
 825                         mul_v3_fl(nor,-100.0);
 826
 827                         VECADD(co2,co1,nor);
 828
 829                         min_d=2.0;
 830                         intersect=0;
 831
 832                         for(i=0,mface=dm->getTessFaceDataArray(dm,CD_MFACE); i<tot; i++,mface++){
 833                                 if(i==pa->num) continue;
 834
 835                                 v1=mvert[mface->v1].co;
 836                                 v2=mvert[mface->v2].co;
 837                                 v3=mvert[mface->v3].co;
 838
 839                                 if(isect_line_tri_v3(co1, co2, v2, v3, v1, &cur_d, 0)){
 840                                         if(cur_d<min_d){
 841                                                 min_d=cur_d;
 842                                                 pa->foffset=cur_d*50.0f; /* to the middle of volume */
 843                                                 intersect=1;
 844                                         }
 845                                 }
 846                                 if(mface->v4){
 847                                         v4=mvert[mface->v4].co;
 848
 849                                         if(isect_line_tri_v3(co1, co2, v4, v1, v3, &cur_d, 0)){
 850                                                 if(cur_d<min_d){
 851                                                         min_d=cur_d;
 852                                                         pa->foffset=cur_d*50.0f; /* to the middle of volume */
 853                                                         intersect=1;
 854                                                 }
 855                                         }
 856                                 }
 857                         }
 858                         if(intersect==0)
 859                                 pa->foffset=0.0;
 860                         else switch(distr){
 861                                 case PART_DISTR_JIT:
 862                                         pa->foffset*= ctx->jit[p%(2*ctx->jitlevel)];
 863                                         break;
 864                                 case PART_DISTR_RAND:
 865                                         pa->foffset*=BLI_frand();
 866                                         break;
 867                         }
 868                 }
 869         }
 870         else if(from == PART_FROM_CHILD) {
 871                 MFace *mf;
 872
 873                 if(ctx->index[p] < 0) {
 874                         cpa->num=0;
 875                         cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]=0.0f;
 876                         cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
 877                         return;
 878                 }
 879
 880                 mf= dm->getTessFaceData(dm, ctx->index[p], CD_MFACE);
 881
 882                 randu= rng_getFloat(thread->rng);
 883                 randv= rng_getFloat(thread->rng);
 884                 rng_skip_tot -= 2;
 885
 886                 psys_uv_to_w(randu, randv, mf->v4, cpa->fuv);
 887
 888                 cpa->num = ctx->index[p];
 889
 890                 if(ctx->tree){
 891                         KDTreeNearest ptn[10];
 892                         int w,maxw;//, do_seams;
 893                         float maxd /*, mind,dd */, totw= 0.0f;
 894                         int parent[10];
 895                         float pweight[10];
 896
 897                         psys_particle_on_dm(dm,cfrom,cpa->num,DMCACHE_ISCHILD,cpa->fuv,cpa->foffset,co1,nor1,NULL,NULL,orco1,NULL);
 898                         transform_mesh_orco_verts((Mesh*)ob->data, &orco1, 1, 1);
 899                         maxw = BLI_kdtree_find_n_nearest(ctx->tree,4,orco1,NULL,ptn);
 900
 901                         maxd=ptn[maxw-1].dist;
 902                         /* mind=ptn[0].dist; */ /* UNUSED */
 903
 904                         /* the weights here could be done better */
 905                         for(w=0; w<maxw; w++){
 906                                 parent[w]=ptn[w].index;
 907                                 pweight[w]=(float)pow(2.0,(double)(-6.0f*ptn[w].dist/maxd));
 908                         }
 909                         for(;w<10; w++){
 910                                 parent[w]=-1;
 911                                 pweight[w]=0.0f;
 912                         }
 913
 914                         for(w=0,i=0; w<maxw && i<4; w++){
 915                                 if(parent[w]>=0){
 916                                         cpa->pa[i]=parent[w];
 917                                         cpa->w[i]=pweight[w];
 918                                         totw+=pweight[w];
 919                                         i++;
 920                                 }
 921                         }
 922                         for(;i<4; i++){
 923                                 cpa->pa[i]=-1;
 924                                 cpa->w[i]=0.0f;
 925                         }
 926
 927                         if(totw>0.0f) for(w=0; w<4; w++)
 928                                 cpa->w[w]/=totw;
 929
 930                         cpa->parent=cpa->pa[0];
 931                 }
 932         }
 933
 934         if(rng_skip_tot > 0) /* should never be below zero */
 935                 rng_skip(thread->rng, rng_skip_tot);
 936 }
 937
 938 static void *distribute_threads_exec_cb(void *data)
 939 {
 940         ParticleThread *thread= (ParticleThread*)data;
 941         ParticleSystem *psys= thread->ctx->sim.psys;
 942         ParticleData *pa;
 943         ChildParticle *cpa;
 944         int p, totpart;
 945
 946         if(thread->ctx->from == PART_FROM_CHILD) {
 947                 totpart= psys->totchild;
 948                 cpa= psys->child;
 949
 950                 for(p=0; p<totpart; p++, cpa++) {
 951                         if(thread->ctx->skip) /* simplification skip */
 952                                 rng_skip(thread->rng, PSYS_RND_DIST_SKIP * thread->ctx->skip[p]);
 953
 954                         if((p+thread->num) % thread->tot == 0)
 955                                 distribute_threads_exec(thread, NULL, cpa, p);
 956                         else /* thread skip */
 957                                 rng_skip(thread->rng, PSYS_RND_DIST_SKIP);
 958                 }
 959         }
 960         else {
 961                 totpart= psys->totpart;
 962                 pa= psys->particles + thread->num;
 963                 for(p=thread->num; p<totpart; p+=thread->tot, pa+=thread->tot)
 964                         distribute_threads_exec(thread, pa, NULL, p);
 965         }
 966
 967         return 0;
 968 }
 969
 970 /* not thread safe, but qsort doesn't take userdata argument */
 971 static int *COMPARE_ORIG_INDEX = NULL;
 972 static int distribute_compare_orig_index(const void *p1, const void *p2)
 973 {
 974         int index1 = COMPARE_ORIG_INDEX[*(const int*)p1];
 975         int index2 = COMPARE_ORIG_INDEX[*(const int*)p2];
 976
 977         if(index1 < index2)
 978                 return -1;
 979         else if(index1 == index2) {
 980                 /* this pointer comparison appears to make qsort stable for glibc,
 981                  * and apparently on solaris too, makes the renders reproducable */
 982                 if(p1 < p2)
 983                         return -1;
 984                 else if(p1 == p2)
 985                         return 0;
 986                 else
 987                         return 1;
 988         }
 989         else
 990                 return 1;
 991 }
 992
 993 static void distribute_invalid(Scene *scene, ParticleSystem *psys, int from)
 994 {
 995         if(from == PART_FROM_CHILD) {
 996                 ChildParticle *cpa;
 997                 int p, totchild = get_psys_tot_child(scene, psys);
 998
 999                 if(psys->child && totchild) {
1000                         for(p=0,cpa=psys->child; p<totchild; p++,cpa++){
1001                                 cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]= 0.0;
1002                                 cpa->foffset= 0.0f;
1003                                 cpa->parent=0;
1004                                 cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
1005                                 cpa->num= -1;
1006                         }
1007                 }
1008         }
1009         else {
1010                 PARTICLE_P;
1011                 LOOP_PARTICLES {
1012                         pa->fuv[0]=pa->fuv[1]=pa->fuv[2]= pa->fuv[3]= 0.0;
1013                         pa->foffset= 0.0f;
1014                         pa->num= -1;
1015                 }
1016         }
1017 }
1018
1019 /* Creates a distribution of coordinates on a DerivedMesh       */
1020 /* This is to denote functionality that does not yet work with mesh - only derived mesh */
1021 static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, DerivedMesh *finaldm, int from)
1022 {
1023         ParticleThreadContext *ctx= threads[0].ctx;
1024         Object *ob= ctx->sim.ob;
1025         ParticleSystem *psys= ctx->sim.psys;
1026         ParticleData *pa=0, *tpars= 0;
1027         ParticleSettings *part;
1028         ParticleSeam *seams= 0;
1029         KDTree *tree=0;
1030         DerivedMesh *dm= NULL;
1031         float *jit= NULL;
1032         int i, seed, p=0, totthread= threads[0].tot;
1033         int cfrom=0;
1034         int totelem=0, totpart, *particle_element=0, children=0, totseam=0;
1035         int jitlevel= 1, distr;
1036         float *element_weight=NULL,*element_sum=NULL,*jitter_offset=NULL, *vweight=NULL;
1037         float cur, maxweight=0.0, tweight, totweight, inv_totweight, co[3], nor[3], orco[3], ornor[3];
1038
1039         if(ELEM3(NULL, ob, psys, psys->part))
1040                 return 0;
1041
1042         part=psys->part;
1043         totpart=psys->totpart;
1044         if(totpart==0)
1045                 return 0;
1046
1047         if (!finaldm->deformedOnly && !finaldm->getTessFaceDataArray(finaldm, CD_ORIGINDEX)) {
1048                 printf("Can't create particles with the current modifier stack, disable destructive modifiers\n");
1049 // XXX          error("Can't paint with the current modifier stack, disable destructive modifiers");
1050                 return 0;
1051         }
1052
1053         /* First handle special cases */
1054         if(from == PART_FROM_CHILD) {
1055                 /* Simple children */
1056                 if(part->childtype != PART_CHILD_FACES) {
1057                         BLI_srandom(31415926 + psys->seed + psys->child_seed);
1058                         distribute_simple_children(scene, ob, finaldm, psys);
1059                         return 0;
1060                 }
1061         }
1062         else {
1063                 /* Grid distribution */
1064                 if(part->distr==PART_DISTR_GRID && from != PART_FROM_VERT){
1065                         BLI_srandom(31415926 + psys->seed);
1066                         dm= CDDM_from_mesh((Mesh*)ob->data, ob);
1067                         distribute_grid(dm,psys);
1068                         dm->release(dm);
1069                         return 0;
1070                 }
1071         }
1072
1073         /* Create trees and original coordinates if needed */
1074         if(from == PART_FROM_CHILD) {
1075                 distr=PART_DISTR_RAND;
1076                 BLI_srandom(31415926 + psys->seed + psys->child_seed);
1077                 dm= finaldm;
1078                 children=1;
1079
1080                 tree=BLI_kdtree_new(totpart);
1081
1082                 for(p=0,pa=psys->particles; p<totpart; p++,pa++){
1083                         psys_particle_on_dm(dm,part->from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co,nor,0,0,orco,ornor);
1084                         transform_mesh_orco_verts((Mesh*)ob->data, &orco, 1, 1);
1085                         BLI_kdtree_insert(tree, p, orco, ornor);
1086                 }
1087
1088                 BLI_kdtree_balance(tree);
1089
1090                 totpart = get_psys_tot_child(scene, psys);
1091                 cfrom = from = PART_FROM_FACE;
1092         }
1093         else {
1094                 distr = part->distr;
1095                 BLI_srandom(31415926 + psys->seed);
1096
1097                 dm= CDDM_from_mesh((Mesh*)ob->data, ob);
1098
1099                 /* we need orco for consistent distributions */
1100                 DM_add_vert_layer(dm, CD_ORCO, CD_ASSIGN, get_mesh_orco_verts(ob));
1101
1102                 if(from == PART_FROM_VERT) {
1103                         MVert *mv= dm->getVertDataArray(dm, CD_MVERT);
1104                         float (*orcodata)[3]= dm->getVertDataArray(dm, CD_ORCO);
1105                         int totvert = dm->getNumVerts(dm);
1106
1107                         tree=BLI_kdtree_new(totvert);
1108
1109                         for(p=0; p<totvert; p++) {
1110                                 if(orcodata) {
1111                                         VECCOPY(co,orcodata[p])
1112                                         transform_mesh_orco_verts((Mesh*)ob->data, &co, 1, 1);
1113                                 }
1114                                 else
1115                                         VECCOPY(co,mv[p].co)
1116                                 BLI_kdtree_insert(tree,p,co,NULL);
1117                         }
1118
1119                         BLI_kdtree_balance(tree);
1120                 }
1121         }
1122
1123         /* Get total number of emission elements and allocate needed arrays */
1124         totelem = (from == PART_FROM_VERT) ? dm->getNumVerts(dm) : dm->getNumTessFaces(dm);
1125
1126         if(totelem == 0){
1127                 distribute_invalid(scene, psys, children ? PART_FROM_CHILD : 0);
1128
1129                 if(G.f & G_DEBUG)
1130                         fprintf(stderr,"Particle distribution error: Nothing to emit from!\n");
1131
1132                 if(dm != finaldm) dm->release(dm);
1133                 return 0;
1134         }
1135
1136         element_weight  = MEM_callocN(sizeof(float)*totelem, "particle_distribution_weights");
1137         particle_element= MEM_callocN(sizeof(int)*totpart, "particle_distribution_indexes");
1138         element_sum             = MEM_callocN(sizeof(float)*(totelem+1), "particle_distribution_sum");
1139         jitter_offset   = MEM_callocN(sizeof(float)*totelem, "particle_distribution_jitoff");
1140
1141         /* Calculate weights from face areas */
1142         if((part->flag&PART_EDISTR || children) && from != PART_FROM_VERT){
1143                 MVert *v1, *v2, *v3, *v4;
1144                 float totarea=0.f, co1[3], co2[3], co3[3], co4[3];
1145                 float (*orcodata)[3];
1146
1147                 orcodata= dm->getVertDataArray(dm, CD_ORCO);
1148
1149                 for(i=0; i<totelem; i++){
1150                         MFace *mf=dm->getTessFaceData(dm,i,CD_MFACE);
1151
1152                         if(orcodata) {
1153                                 VECCOPY(co1, orcodata[mf->v1]);
1154                                 VECCOPY(co2, orcodata[mf->v2]);
1155                                 VECCOPY(co3, orcodata[mf->v3]);
1156                                 transform_mesh_orco_verts((Mesh*)ob->data, &co1, 1, 1);
1157                                 transform_mesh_orco_verts((Mesh*)ob->data, &co2, 1, 1);
1158                                 transform_mesh_orco_verts((Mesh*)ob->data, &co3, 1, 1);
1159                                 if(mf->v4) {
1160                                         VECCOPY(co4, orcodata[mf->v4]);
1161                                         transform_mesh_orco_verts((Mesh*)ob->data, &co4, 1, 1);
1162                                 }
1163                         }
1164                         else {
1165                                 v1= (MVert*)dm->getVertData(dm,mf->v1,CD_MVERT);
1166                                 v2= (MVert*)dm->getVertData(dm,mf->v2,CD_MVERT);
1167                                 v3= (MVert*)dm->getVertData(dm,mf->v3,CD_MVERT);
1168                                 VECCOPY(co1, v1->co);
1169                                 VECCOPY(co2, v2->co);
1170                                 VECCOPY(co3, v3->co);
1171                                 if(mf->v4) {
1172                                         v4= (MVert*)dm->getVertData(dm,mf->v4,CD_MVERT);
1173                                         VECCOPY(co4, v4->co);
1174                                 }
1175                         }
1176
1177                         cur = mf->v4 ? area_quad_v3(co1, co2, co3, co4) : area_tri_v3(co1, co2, co3);
1178
1179                         if(cur > maxweight)
1180                                 maxweight = cur;
1181
1182                         element_weight[i] = cur;
1183                         totarea += cur;
1184                 }
1185
1186                 for(i=0; i<totelem; i++)
1187                         element_weight[i] /= totarea;
1188
1189                 maxweight /= totarea;
1190         }
1191         else{
1192                 float min=1.0f/(float)(MIN2(totelem,totpart));
1193                 for(i=0; i<totelem; i++)
1194                         element_weight[i]=min;
1195                 maxweight=min;
1196         }
1197
1198         /* Calculate weights from vgroup */
1199         vweight = psys_cache_vgroup(dm,psys,PSYS_VG_DENSITY);
1200
1201         if(vweight){
1202                 if(from==PART_FROM_VERT) {
1203                         for(i=0;i<totelem; i++)
1204                                 element_weight[i]*=vweight[i];
1205                 }
1206                 else { /* PART_FROM_FACE / PART_FROM_VOLUME */
1207                         for(i=0;i<totelem; i++){
1208                                 MFace *mf=dm->getTessFaceData(dm,i,CD_MFACE);
1209                                 tweight = vweight[mf->v1] + vweight[mf->v2] + vweight[mf->v3];
1210
1211                                 if(mf->v4) {
1212                                         tweight += vweight[mf->v4];
1213                                         tweight /= 4.0f;
1214                                 }
1215                                 else {
1216                                         tweight /= 3.0f;
1217                                 }
1218
1219                                 element_weight[i]*=tweight;
1220                         }
1221                 }
1222                 MEM_freeN(vweight);
1223         }
1224
1225         /* Calculate total weight of all elements */
1226         totweight= 0.0f;
1227         for(i=0;i<totelem; i++)
1228                 totweight += element_weight[i];
1229
1230         inv_totweight = (totweight > 0.f ? 1.f/totweight : 0.f);
1231
1232         /* Calculate cumulative weights */
1233         element_sum[0]= 0.0f;
1234         for(i=0; i<totelem; i++)
1235                 element_sum[i+1]= element_sum[i] + element_weight[i] * inv_totweight;
1236
1237         /* Finally assign elements to particles */
1238         if((part->flag&PART_TRAND) || (part->simplify_flag&PART_SIMPLIFY_ENABLE)) {
1239                 float pos;
1240
1241                 for(p=0; p<totpart; p++) {
1242                         /* In theory element_sum[totelem] should be 1.0, but due to float errors this is not necessarily always true, so scale pos accordingly. */
1243                         pos= BLI_frand() * element_sum[totelem];
1244                         particle_element[p]= distribute_binary_search(element_sum, totelem, pos);
1245                         particle_element[p]= MIN2(totelem-1, particle_element[p]);
1246                         jitter_offset[particle_element[p]]= pos;
1247                 }
1248         }
1249         else {
1250                 double step, pos;
1251
1252                 step= (totpart < 2) ? 0.5 : 1.0/(double)totpart;
1253                 pos= 1e-6; /* tiny offset to avoid zero weight face */
1254                 i= 0;
1255
1256                 for(p=0; p<totpart; p++, pos+=step) {
1257                         while((i < totelem) && (pos > element_sum[i+1]))
1258                                 i++;
1259
1260                         particle_element[p]= MIN2(totelem-1, i);
1261
1262                         /* avoid zero weight face */
1263                         if(p == totpart-1 && element_weight[particle_element[p]] == 0.0f)
1264                                 particle_element[p]= particle_element[p-1];
1265
1266                         jitter_offset[particle_element[p]]= pos;
1267                 }
1268         }
1269
1270         MEM_freeN(element_sum);
1271
1272         /* For hair, sort by origindex (allows optimizations in rendering), */
1273         /* however with virtual parents the children need to be in random order. */
1274         if(part->type == PART_HAIR && !(part->childtype==PART_CHILD_FACES && part->parents!=0.0f)) {
1275                 COMPARE_ORIG_INDEX = NULL;
1276
1277                 if(from == PART_FROM_VERT) {
1278                         if(dm->numVertData)
1279                                 COMPARE_ORIG_INDEX= dm->getVertDataArray(dm, CD_ORIGINDEX);
1280                 }
1281                 else {
1282                         if(dm->numFaceData)
1283                                 COMPARE_ORIG_INDEX= dm->getTessFaceDataArray(dm, CD_ORIGINDEX);
1284                 }
1285
1286                 if(COMPARE_ORIG_INDEX) {
1287                         qsort(particle_element, totpart, sizeof(int), distribute_compare_orig_index);
1288                         COMPARE_ORIG_INDEX = NULL;
1289                 }
1290         }
1291
1292         /* Create jittering if needed */
1293         if(distr==PART_DISTR_JIT && ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
1294                 jitlevel= part->userjit;
1295
1296                 if(jitlevel == 0) {
1297                         jitlevel= totpart/totelem;
1298                         if(part->flag & PART_EDISTR) jitlevel*= 2;      /* looks better in general, not very scietific */
1299                         if(jitlevel<3) jitlevel= 3;
1300                 }
1301
1302                 jit= MEM_callocN((2+ jitlevel*2)*sizeof(float), "jit");
1303
1304                 /* for small amounts of particles we use regular jitter since it looks
1305                  * a bit better, for larger amounts we switch to hammersley sequence
1306                  * because it is much faster */
1307                 if(jitlevel < 25)
1308                         init_mv_jit(jit, jitlevel, psys->seed, part->jitfac);
1309                 else
1310                         hammersley_create(jit, jitlevel+1, psys->seed, part->jitfac);
1311                 BLI_array_randomize(jit, 2*sizeof(float), jitlevel, psys->seed); /* for custom jit or even distribution */
1312         }
1313
1314         /* Setup things for threaded distribution */
1315         ctx->tree= tree;
1316         ctx->seams= seams;
1317         ctx->totseam= totseam;
1318         ctx->sim.psys= psys;
1319         ctx->index= particle_element;
1320         ctx->jit= jit;
1321         ctx->jitlevel= jitlevel;
1322         ctx->jitoff= jitter_offset;
1323         ctx->weight= element_weight;
1324         ctx->maxweight= maxweight;
1325         ctx->from= (children)? PART_FROM_CHILD: from;
1326         ctx->cfrom= cfrom;
1327         ctx->distr= distr;
1328         ctx->dm= dm;
1329         ctx->tpars= tpars;
1330
1331         if(children) {
1332                 totpart= psys_render_simplify_distribution(ctx, totpart);
1333                 alloc_child_particles(psys, totpart);
1334         }
1335
1336         if(!children || psys->totchild < 10000)
1337                 totthread= 1;
1338
1339         seed= 31415926 + ctx->sim.psys->seed;
1340         for(i=0; i<totthread; i++) {
1341                 threads[i].rng= rng_new(seed);
1342                 threads[i].tot= totthread;
1343         }
1344
1345         return 1;
1346 }
1347
1348 static void distribute_particles_on_dm(ParticleSimulationData *sim, int from)
1349 {
1350         DerivedMesh *finaldm = sim->psmd->dm;
1351         ListBase threads;
1352         ParticleThread *pthreads;
1353         ParticleThreadContext *ctx;
1354         int i, totthread;
1355
1356         pthreads= psys_threads_create(sim);
1357
1358         if(!distribute_threads_init_data(pthreads, sim->scene, finaldm, from)) {
1359                 psys_threads_free(pthreads);
1360                 return;
1361         }
1362
1363         totthread= pthreads[0].tot;
1364         if(totthread > 1) {
1365                 BLI_init_threads(&threads, distribute_threads_exec_cb, totthread);
1366
1367                 for(i=0; i<totthread; i++)
1368                         BLI_insert_thread(&threads, &pthreads[i]);
1369
1370                 BLI_end_threads(&threads);
1371         }
1372         else
1373                 distribute_threads_exec_cb(&pthreads[0]);
1374
1375         psys_calc_dmcache(sim->ob, finaldm, sim->psys);
1376
1377         ctx= pthreads[0].ctx;
1378         if(ctx->dm != finaldm)
1379                 ctx->dm->release(ctx->dm);
1380
1381         psys_threads_free(pthreads);
1382 }
1383
1384 /* ready for future use, to emit particles without geometry */
1385 static void distribute_particles_on_shape(ParticleSimulationData *sim, int UNUSED(from))
1386 {
1387         distribute_invalid(sim->scene, sim->psys, 0);
1388
1389         fprintf(stderr,"Shape emission not yet possible!\n");
1390 }
1391
1392 static void distribute_particles(ParticleSimulationData *sim, int from)
1393 {
1394         PARTICLE_PSMD;
1395         int distr_error=0;
1396
1397         if(psmd){
1398                 if(psmd->dm)
1399                         distribute_particles_on_dm(sim, from);
1400                 else
1401                         distr_error=1;
1402         }
1403         else
1404                 distribute_particles_on_shape(sim, from);
1405
1406         if(distr_error){
1407                 distribute_invalid(sim->scene, sim->psys, from);
1408
1409                 fprintf(stderr,"Particle distribution error!\n");
1410         }
1411 }
1412
1413 /* threaded child particle distribution and path caching */
1414 ParticleThread *psys_threads_create(ParticleSimulationData *sim)
1415 {
1416         ParticleThread *threads;
1417         ParticleThreadContext *ctx;
1418         int i, totthread;
1419
1420         if(sim->scene->r.mode & R_FIXED_THREADS)
1421                 totthread= sim->scene->r.threads;
1422         else
1423                 totthread= BLI_system_thread_count();
1424
1425         threads= MEM_callocN(sizeof(ParticleThread)*totthread, "ParticleThread");
1426         ctx= MEM_callocN(sizeof(ParticleThreadContext), "ParticleThreadContext");
1427
1428         ctx->sim = *sim;
1429         ctx->dm= ctx->sim.psmd->dm;
1430         ctx->ma= give_current_material(sim->ob, sim->psys->part->omat);
1431
1432         memset(threads, 0, sizeof(ParticleThread)*totthread);
1433
1434         for(i=0; i<totthread; i++) {
1435                 threads[i].ctx= ctx;
1436                 threads[i].num= i;
1437                 threads[i].tot= totthread;
1438         }
1439
1440         return threads;
1441 }
1442
1443 void psys_threads_free(ParticleThread *threads)
1444 {
1445         ParticleThreadContext *ctx= threads[0].ctx;
1446         int i, totthread= threads[0].tot;
1447
1448         /* path caching */
1449         if(ctx->vg_length)
1450                 MEM_freeN(ctx->vg_length);
1451         if(ctx->vg_clump)
1452                 MEM_freeN(ctx->vg_clump);
1453         if(ctx->vg_kink)
1454                 MEM_freeN(ctx->vg_kink);
1455         if(ctx->vg_rough1)
1456                 MEM_freeN(ctx->vg_rough1);
1457         if(ctx->vg_rough2)
1458                 MEM_freeN(ctx->vg_rough2);
1459         if(ctx->vg_roughe)
1460                 MEM_freeN(ctx->vg_roughe);
1461
1462         if(ctx->sim.psys->lattice){
1463                 end_latt_deform(ctx->sim.psys->lattice);
1464                 ctx->sim.psys->lattice= NULL;
1465         }
1466
1467         /* distribution */
1468         if(ctx->jit) MEM_freeN(ctx->jit);
1469         if(ctx->jitoff) MEM_freeN(ctx->jitoff);
1470         if(ctx->weight) MEM_freeN(ctx->weight);
1471         if(ctx->index) MEM_freeN(ctx->index);
1472         if(ctx->skip) MEM_freeN(ctx->skip);
1473         if(ctx->seams) MEM_freeN(ctx->seams);
1474         //if(ctx->vertpart) MEM_freeN(ctx->vertpart);
1475         BLI_kdtree_free(ctx->tree);
1476
1477         /* threads */
1478         for(i=0; i<totthread; i++) {
1479                 if(threads[i].rng)
1480                         rng_free(threads[i].rng);
1481                 if(threads[i].rng_path)
1482                         rng_free(threads[i].rng_path);
1483         }
1484
1485         MEM_freeN(ctx);
1486         MEM_freeN(threads);
1487 }
1488
1489 /* set particle parameters that don't change during particle's life */
1490 void initialize_particle(ParticleSimulationData *sim, ParticleData *pa, int p)
1491 {
1492         ParticleSystem *psys = sim->psys;
1493         ParticleSettings *part = psys->part;
1494         ParticleTexture ptex;
1495
1496         pa->flag &= ~PARS_UNEXIST;
1497
1498         if(part->type != PART_FLUID) {
1499                 psys_get_texture(sim, pa, &ptex, PAMAP_INIT, 0.f);
1500
1501                 if(ptex.exist < PSYS_FRAND(p+125))
1502                         pa->flag |= PARS_UNEXIST;
1503
1504                 pa->time = (part->type == PART_HAIR) ? 0.f : part->sta + (part->end - part->sta)*ptex.time;
1505         }
1506
1507         pa->hair_index = 0;
1508         /* we can't reset to -1 anymore since we've figured out correct index in distribute_particles */
1509         /* usage other than straight after distribute has to handle this index by itself - jahka*/
1510         //pa->num_dmcache = DMCACHE_NOTFOUND; /* assume we dont have a derived mesh face */
1511 }
1512 static void initialize_all_particles(ParticleSimulationData *sim)
1513 {
1514         ParticleSystem *psys = sim->psys;
1515         PARTICLE_P;
1516
1517         psys->totunexist = 0;
1518
1519         LOOP_PARTICLES {
1520                 if((pa->flag & PARS_UNEXIST)==0)
1521                         initialize_particle(sim, pa, p);
1522
1523                 if(pa->flag & PARS_UNEXIST)
1524                         psys->totunexist++;
1525         }
1526
1527         /* Free unexisting particles. */
1528         if(psys->totpart && psys->totunexist == psys->totpart) {
1529                 if(psys->particles->boid)
1530                         MEM_freeN(psys->particles->boid);
1531
1532                 MEM_freeN(psys->particles);
1533                 psys->particles = NULL;
1534                 psys->totpart = psys->totunexist = 0;
1535         }
1536
1537         if(psys->totunexist) {
1538                 int newtotpart = psys->totpart - psys->totunexist;
1539                 ParticleData *npa, *newpars;
1540
1541                 npa = newpars = MEM_callocN(newtotpart * sizeof(ParticleData), "particles");
1542
1543                 for(p=0, pa=psys->particles; p<newtotpart; p++, pa++, npa++) {
1544                         while(pa->flag & PARS_UNEXIST)
1545                                 pa++;
1546
1547                         memcpy(npa, pa, sizeof(ParticleData));
1548                 }
1549
1550                 if(psys->particles->boid)
1551                         MEM_freeN(psys->particles->boid);
1552                 MEM_freeN(psys->particles);
1553                 psys->particles = newpars;
1554                 psys->totpart -= psys->totunexist;
1555
1556                 if(psys->particles->boid) {
1557                         BoidParticle *newboids = MEM_callocN(psys->totpart * sizeof(BoidParticle), "boid particles");
1558
1559                         LOOP_PARTICLES
1560                                 pa->boid = newboids++;
1561
1562                 }
1563         }
1564 }
1565 void psys_get_birth_coordinates(ParticleSimulationData *sim, ParticleData *pa, ParticleKey *state, float dtime, float cfra)
1566 {
1567         Object *ob = sim->ob;
1568         ParticleSystem *psys = sim->psys;
1569         ParticleSettings *part;
1570         ParticleTexture ptex;
1571         float fac, phasefac, nor[3]={0,0,0},loc[3],vel[3]={0.0,0.0,0.0},rot[4],q2[4];
1572         float r_vel[3],r_ave[3],r_rot[4],vec[3],p_vel[3]={0.0,0.0,0.0};
1573         float x_vec[3]={1.0,0.0,0.0}, utan[3]={0.0,1.0,0.0}, vtan[3]={0.0,0.0,1.0}, rot_vec[3]={0.0,0.0,0.0};
1574         float q_phase[4];
1575         int p = pa - psys->particles;
1576         part=psys->part;
1577
1578         /* get birth location from object               */
1579         if(part->tanfac != 0.f)
1580                 psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,utan,vtan,0,0);
1581         else
1582                 psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,0,0,0,0);
1583
1584         /* get possible textural influence */
1585         psys_get_texture(sim, pa, &ptex, PAMAP_IVEL, cfra);
1586
1587         /* particles live in global space so    */
1588         /* let's convert:                                               */
1589         /* -location                                                    */
1590         mul_m4_v3(ob->obmat, loc);
1591
1592         /* -normal                                                              */
1593         mul_mat3_m4_v3(ob->obmat, nor);
1594         normalize_v3(nor);
1595
1596         /* -tangent                                                             */
1597         if(part->tanfac!=0.0f){
1598                 //float phase=vg_rot?2.0f*(psys_particle_value_from_verts(sim->psmd->dm,part->from,pa,vg_rot)-0.5f):0.0f;
1599                 float phase=0.0f;
1600                 mul_v3_fl(vtan,-(float)cos((float)M_PI*(part->tanphase+phase)));
1601                 fac=-(float)sin((float)M_PI*(part->tanphase+phase));
1602                 VECADDFAC(vtan,vtan,utan,fac);
1603
1604                 mul_mat3_m4_v3(ob->obmat,vtan);
1605
1606                 VECCOPY(utan,nor);
1607                 mul_v3_fl(utan,dot_v3v3(vtan,nor));
1608                 VECSUB(vtan,vtan,utan);
1609
1610                 normalize_v3(vtan);
1611         }
1612
1613
1614         /* -velocity                                                    */
1615         if(part->randfac != 0.0f){
1616                 r_vel[0] = 2.0f * (PSYS_FRAND(p + 10) - 0.5f);
1617                 r_vel[1] = 2.0f * (PSYS_FRAND(p + 11) - 0.5f);
1618                 r_vel[2] = 2.0f * (PSYS_FRAND(p + 12) - 0.5f);
1619
1620                 mul_mat3_m4_v3(ob->obmat, r_vel);
1621                 normalize_v3(r_vel);
1622         }
1623
1624         /* -angular velocity                                    */
1625         if(part->avemode==PART_AVE_RAND){
1626                 r_ave[0] = 2.0f * (PSYS_FRAND(p + 13) - 0.5f);
1627                 r_ave[1] = 2.0f * (PSYS_FRAND(p + 14) - 0.5f);
1628                 r_ave[2] = 2.0f * (PSYS_FRAND(p + 15) - 0.5f);
1629
1630                 mul_mat3_m4_v3(ob->obmat,r_ave);
1631                 normalize_v3(r_ave);
1632         }
1633
1634         /* -rotation                                                    */
1635         if(part->randrotfac != 0.0f){
1636                 r_rot[0] = 2.0f * (PSYS_FRAND(p + 16) - 0.5f);
1637                 r_rot[1] = 2.0f * (PSYS_FRAND(p + 17) - 0.5f);
1638                 r_rot[2] = 2.0f * (PSYS_FRAND(p + 18) - 0.5f);
1639                 r_rot[3] = 2.0f * (PSYS_FRAND(p + 19) - 0.5f);
1640                 normalize_qt(r_rot);
1641
1642                 mat4_to_quat(rot,ob->obmat);
1643                 mul_qt_qtqt(r_rot,r_rot,rot);
1644         }
1645
1646         if(part->phystype==PART_PHYS_BOIDS && pa->boid) {
1647                 float dvec[3], q[4], mat[3][3];
1648
1649                 copy_v3_v3(state->co,loc);
1650
1651                 /* boids don't get any initial velocity  */
1652                 zero_v3(state->vel);
1653
1654                 /* boids store direction in ave */
1655                 if(fabsf(nor[2])==1.0f) {
1656                         sub_v3_v3v3(state->ave, loc, ob->obmat[3]);
1657                         normalize_v3(state->ave);
1658                 }
1659                 else {
1660                         VECCOPY(state->ave, nor);
1661                 }
1662
1663                 /* calculate rotation matrix */
1664                 project_v3_v3v3(dvec, r_vel, state->ave);
1665                 sub_v3_v3v3(mat[0], state->ave, dvec);
1666                 normalize_v3(mat[0]);
1667                 negate_v3_v3(mat[2], r_vel);
1668                 normalize_v3(mat[2]);
1669                 cross_v3_v3v3(mat[1], mat[2], mat[0]);
1670
1671                 /* apply rotation */
1672                 mat3_to_quat_is_ok( q,mat);
1673                 copy_qt_qt(state->rot, q);
1674         }
1675         else {
1676                 /* conversion done so now we apply new: */
1677                 /* -velocity from:                                              */
1678
1679                 /*              *reactions                                              */
1680                 if(dtime > 0.f){
1681                         sub_v3_v3v3(vel, pa->state.vel, pa->prev_state.vel);
1682                 }
1683
1684                 /*              *emitter velocity                               */
1685                 if(dtime != 0.f && part->obfac != 0.f){
1686                         sub_v3_v3v3(vel, loc, state->co);
1687                         mul_v3_fl(vel, part->obfac/dtime);
1688                 }
1689
1690                 /*              *emitter normal                                 */
1691                 if(part->normfac != 0.f)
1692                         madd_v3_v3fl(vel, nor, part->normfac);
1693
1694                 /*              *emitter tangent                                */
1695                 if(sim->psmd && part->tanfac != 0.f)
1696                         madd_v3_v3fl(vel, vtan, part->tanfac);
1697
1698                 /*              *emitter object orientation             */
1699                 if(part->ob_vel[0] != 0.f) {
1700                         normalize_v3_v3(vec, ob->obmat[0]);
1701                         madd_v3_v3fl(vel, vec, part->ob_vel[0]);
1702                 }
1703                 if(part->ob_vel[1] != 0.f) {
1704                         normalize_v3_v3(vec, ob->obmat[1]);
1705                         madd_v3_v3fl(vel, vec, part->ob_vel[1]);
1706                 }
1707                 if(part->ob_vel[2] != 0.f) {
1708                         normalize_v3_v3(vec, ob->obmat[2]);
1709                         madd_v3_v3fl(vel, vec, part->ob_vel[2]);
1710                 }
1711
1712                 /*              *texture                                                */
1713                 /* TODO */
1714
1715                 /*              *random                                                 */
1716                 if(part->randfac != 0.f)
1717                         madd_v3_v3fl(vel, r_vel, part->randfac);
1718
1719                 /*              *particle                                               */
1720                 if(part->partfac != 0.f)
1721                         madd_v3_v3fl(vel, p_vel, part->partfac);
1722
1723                 mul_v3_v3fl(state->vel, vel, ptex.ivel);
1724
1725                 /* -location from emitter                               */
1726                 copy_v3_v3(state->co,loc);
1727
1728                 /* -rotation                                                    */
1729                 unit_qt(state->rot);
1730
1731                 if(part->rotmode){
1732                         /* create vector into which rotation is aligned */
1733                         switch(part->rotmode){
1734                                 case PART_ROT_NOR:
1735                                         copy_v3_v3(rot_vec, nor);
1736                                         break;
1737                                 case PART_ROT_VEL:
1738                                         copy_v3_v3(rot_vec, vel);
1739                                         break;
1740                                 case PART_ROT_GLOB_X:
1741                                 case PART_ROT_GLOB_Y:
1742                                 case PART_ROT_GLOB_Z:
1743                                         rot_vec[part->rotmode - PART_ROT_GLOB_X] = 1.0f;
1744                                         break;
1745                                 case PART_ROT_OB_X:
1746                                 case PART_ROT_OB_Y:
1747                                 case PART_ROT_OB_Z:
1748                                         copy_v3_v3(rot_vec, ob->obmat[part->rotmode - PART_ROT_OB_X]);
1749                                         break;
1750                         }
1751
1752                         /* create rotation quat */
1753                         negate_v3(rot_vec);
1754                         vec_to_quat( q2,rot_vec, OB_POSX, OB_POSZ);
1755
1756                         /* randomize rotation quat */
1757                         if(part->randrotfac!=0.0f)
1758                                 interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
1759                         else
1760                                 copy_qt_qt(rot,q2);
1761
1762                         /* rotation phase */
1763                         phasefac = part->phasefac;
1764                         if(part->randphasefac != 0.0f)
1765                                 phasefac += part->randphasefac * PSYS_FRAND(p + 20);
1766                         axis_angle_to_quat( q_phase,x_vec, phasefac*(float)M_PI);
1767
1768                         /* combine base rotation & phase */
1769                         mul_qt_qtqt(state->rot, rot, q_phase);
1770                 }
1771
1772                 /* -angular velocity                                    */
1773
1774                 zero_v3(state->ave);
1775
1776                 if(part->avemode){
1777                         switch(part->avemode){
1778                                 case PART_AVE_SPIN:
1779                                         copy_v3_v3(state->ave, vel);
1780                                         break;
1781                                 case PART_AVE_RAND:
1782                                         copy_v3_v3(state->ave, r_ave);
1783                                         break;
1784                         }
1785                         normalize_v3(state->ave);
1786                         mul_v3_fl(state->ave, part->avefac);
1787                 }
1788         }
1789 }
1790 /* sets particle to the emitter surface with initial velocity & rotation */
1791 void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, float cfra)
1792 {
1793         Object *ob = sim->ob;
1794         ParticleSystem *psys = sim->psys;
1795         ParticleSettings *part;
1796         ParticleTexture ptex;
1797         int p = pa - psys->particles;
1798         part=psys->part;
1799
1800         /* get precise emitter matrix if particle is born */
1801         if(part->type!=PART_HAIR && dtime > 0.f && pa->time < cfra && pa->time >= sim->psys->cfra) {
1802                 /* we have to force RECALC_ANIM here since where_is_objec_time only does drivers */
1803                 while(ob) {
1804                         BKE_animsys_evaluate_animdata(&ob->id, ob->adt, pa->time, ADT_RECALC_ANIM);
1805                         ob = ob->parent;
1806                 }
1807                 ob = sim->ob;
1808                 where_is_object_time(sim->scene, ob, pa->time);
1809         }
1810
1811         psys_get_birth_coordinates(sim, pa, &pa->state, dtime, cfra);
1812
1813         if(part->phystype==PART_PHYS_BOIDS && pa->boid) {
1814                 BoidParticle *bpa = pa->boid;
1815
1816                 /* and gravity in r_ve */
1817                 bpa->gravity[0] = bpa->gravity[1] = 0.0f;
1818                 bpa->gravity[2] = -1.0f;
1819                 if((sim->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY)
1820                         && sim->scene->physics_settings.gravity[2]!=0.0f)
1821                         bpa->gravity[2] = sim->scene->physics_settings.gravity[2];
1822
1823                 bpa->data.health = part->boids->health;
1824                 bpa->data.mode = eBoidMode_InAir;
1825                 bpa->data.state_id = ((BoidState*)part->boids->states.first)->id;
1826                 bpa->data.acc[0]=bpa->data.acc[1]=bpa->data.acc[2]=0.0f;
1827         }
1828
1829
1830         if(part->type == PART_HAIR){
1831                 pa->lifetime = 100.0f;
1832         }
1833         else{
1834                 /* get possible textural influence */
1835                 psys_get_texture(sim, pa, &ptex, PAMAP_LIFE, cfra);
1836
1837                 pa->lifetime = part->lifetime * ptex.life;
1838
1839                 if(part->randlife != 0.0f)
1840                         pa->lifetime *= 1.0f - part->randlife * PSYS_FRAND(p + 21);
1841         }
1842
1843         pa->dietime = pa->time + pa->lifetime;
1844
1845         if(sim->psys->pointcache && sim->psys->pointcache->flag & PTCACHE_BAKED &&
1846                 sim->psys->pointcache->mem_cache.first) {
1847                 float dietime = psys_get_dietime_from_cache(sim->psys->pointcache, p);
1848                 pa->dietime = MIN2(pa->dietime, dietime);
1849         }
1850
1851         if(pa->time > cfra)
1852                 pa->alive = PARS_UNBORN;
1853         else if(pa->dietime <= cfra)
1854                 pa->alive = PARS_DEAD;
1855         else
1856                 pa->alive = PARS_ALIVE;
1857
1858         pa->state.time = cfra;
1859 }
1860 static void reset_all_particles(ParticleSimulationData *sim, float dtime, float cfra, int from)
1861 {
1862         ParticleData *pa;
1863         int p, totpart=sim->psys->totpart;
1864
1865         for(p=from, pa=sim->psys->particles+from; p<totpart; p++, pa++)
1866                 reset_particle(sim, pa, dtime, cfra);
1867 }
1868 /************************************************/
1869 /*                      Particle targets                                        */
1870 /************************************************/
1871 ParticleSystem *psys_get_target_system(Object *ob, ParticleTarget *pt)
1872 {
1873         ParticleSystem *psys = NULL;
1874
1875         if(pt->ob == NULL || pt->ob == ob)
1876                 psys = BLI_findlink(&ob->particlesystem, pt->psys-1);
1877         else
1878                 psys = BLI_findlink(&pt->ob->particlesystem, pt->psys-1);
1879
1880         if(psys)
1881                 pt->flag |= PTARGET_VALID;
1882         else
1883                 pt->flag &= ~PTARGET_VALID;
1884
1885         return psys;
1886 }
1887 /************************************************/
1888 /*                      Keyed particles                                         */
1889 /************************************************/
1890 /* Counts valid keyed targets */
1891 void psys_count_keyed_targets(ParticleSimulationData *sim)
1892 {
1893         ParticleSystem *psys = sim->psys, *kpsys;
1894         ParticleTarget *pt = psys->targets.first;
1895         int keys_valid = 1;
1896         psys->totkeyed = 0;
1897
1898         for(; pt; pt=pt->next) {
1899                 kpsys = psys_get_target_system(sim->ob, pt);
1900
1901                 if(kpsys && kpsys->totpart) {
1902                         psys->totkeyed += keys_valid;
1903                         if(psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f)
1904                                 psys->totkeyed += 1;
1905                 }
1906                 else {
1907                         keys_valid = 0;
1908                 }
1909         }
1910
1911         psys->totkeyed *= psys->flag & PSYS_KEYED_TIMING ? 1 : psys->part->keyed_loops;
1912 }
1913
1914 static void set_keyed_keys(ParticleSimulationData *sim)
1915 {
1916         ParticleSystem *psys = sim->psys;
1917         ParticleSimulationData ksim= {0};
1918         ParticleTarget *pt;
1919         PARTICLE_P;
1920         ParticleKey *key;
1921         int totpart = psys->totpart, k, totkeys = psys->totkeyed;
1922         int keyed_flag = 0;
1923
1924         ksim.scene= sim->scene;
1925
1926         /* no proper targets so let's clear and bail out */
1927         if(psys->totkeyed==0) {
1928                 free_keyed_keys(psys);
1929                 psys->flag &= ~PSYS_KEYED;
1930                 return;
1931         }
1932
1933         if(totpart && psys->particles->totkey != totkeys) {
1934                 free_keyed_keys(psys);
1935
1936                 key = MEM_callocN(totpart*totkeys*sizeof(ParticleKey), "Keyed keys");
1937
1938                 LOOP_PARTICLES {
1939                         pa->keys = key;
1940                         pa->totkey = totkeys;
1941                         key += totkeys;
1942                 }
1943         }
1944
1945         psys->flag &= ~PSYS_KEYED;
1946
1947
1948         pt = psys->targets.first;
1949         for(k=0; k<totkeys; k++) {
1950                 ksim.ob = pt->ob ? pt->ob : sim->ob;
1951                 ksim.psys = BLI_findlink(&ksim.ob->particlesystem, pt->psys - 1);
1952                 keyed_flag = (ksim.psys->flag & PSYS_KEYED);
1953                 ksim.psys->flag &= ~PSYS_KEYED;
1954
1955                 LOOP_PARTICLES {
1956                         key = pa->keys + k;
1957                         key->time = -1.0; /* use current time */
1958
1959                         psys_get_particle_state(&ksim, p%ksim.psys->totpart, key, 1);
1960
1961                         if(psys->flag & PSYS_KEYED_TIMING){
1962                                 key->time = pa->time + pt->time;
1963                                 if(pt->duration != 0.0f && k+1 < totkeys) {
1964                                         copy_particle_key(key+1, key, 1);
1965                                         (key+1)->time = pa->time + pt->time + pt->duration;
1966                                 }
1967                         }
1968                         else if(totkeys > 1)
1969                                 key->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
1970                         else
1971                                 key->time = pa->time;
1972                 }
1973
1974                 if(psys->flag & PSYS_KEYED_TIMING && pt->duration!=0.0f)
1975                         k++;
1976
1977                 ksim.psys->flag |= keyed_flag;
1978
1979                 pt = (pt->next && pt->next->flag & PTARGET_VALID)? pt->next : psys->targets.first;
1980         }
1981
1982         psys->flag |= PSYS_KEYED;
1983 }
1984
1985 /************************************************/
1986 /*                      Point Cache                                                     */
1987 /************************************************/
1988 void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
1989 {
1990         PointCache *cache = psys->pointcache;
1991
1992         if(cache->flag & PTCACHE_DISK_CACHE && cache->mem_cache.first == NULL) {
1993                 PTCacheID pid;
1994                 BKE_ptcache_id_from_particles(&pid, ob, psys);
1995                 cache->flag &= ~PTCACHE_DISK_CACHE;
1996                 BKE_ptcache_disk_to_mem(&pid);
1997                 cache->flag |= PTCACHE_DISK_CACHE;
1998         }
1999 }
2000 static void psys_clear_temp_pointcache(ParticleSystem *psys)
2001 {
2002         if(psys->pointcache->flag & PTCACHE_DISK_CACHE)
2003                 BKE_ptcache_free_mem(&psys->pointcache->mem_cache);
2004 }
2005 void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra, int *efra)
2006 {
2007         ParticleSettings *part = psys->part;
2008
2009         *sfra = MAX2(1, (int)part->sta);
2010         *efra = MIN2((int)(part->end + part->lifetime + 1.0f), scene->r.efra);
2011 }
2012
2013 /************************************************/
2014 /*                      Effectors                                                       */
2015 /************************************************/
2016 static void psys_update_particle_bvhtree(ParticleSystem *psys, float cfra)
2017 {
2018         if(psys) {
2019                 PARTICLE_P;
2020                 int totpart = 0;
2021
2022                 if(!psys->bvhtree || psys->bvhtree_frame != cfra) {
2023                         LOOP_SHOWN_PARTICLES {
2024                                 totpart++;
2025                         }
2026
2027                         BLI_bvhtree_free(psys->bvhtree);
2028                         psys->bvhtree = BLI_bvhtree_new(totpart, 0.0, 4, 6);
2029
2030                         LOOP_SHOWN_PARTICLES {
2031                                 if(pa->alive == PARS_ALIVE) {
2032                                         if(pa->state.time == cfra)
2033                                                 BLI_bvhtree_insert(psys->bvhtree, p, pa->prev_state.co, 1);
2034                                         else
2035                                                 BLI_bvhtree_insert(psys->bvhtree, p, pa->state.co, 1);
2036                                 }
2037                         }
2038                         BLI_bvhtree_balance(psys->bvhtree);
2039
2040                         psys->bvhtree_frame = cfra;
2041                 }
2042         }
2043 }
2044 void psys_update_particle_tree(ParticleSystem *psys, float cfra)
2045 {
2046         if(psys) {
2047                 PARTICLE_P;
2048                 int totpart = 0;
2049
2050                 if(!psys->tree || psys->tree_frame != cfra) {
2051                         LOOP_SHOWN_PARTICLES {
2052                                 totpart++;
2053                         }
2054
2055                         BLI_kdtree_free(psys->tree);
2056                         psys->tree = BLI_kdtree_new(psys->totpart);
2057
2058                         LOOP_SHOWN_PARTICLES {
2059                                 if(pa->alive == PARS_ALIVE) {
2060                                         if(pa->state.time == cfra)
2061                                                 BLI_kdtree_insert(psys->tree, p, pa->prev_state.co, NULL);
2062                                         else
2063                                                 BLI_kdtree_insert(psys->tree, p, pa->state.co, NULL);
2064                                 }
2065                         }
2066                         BLI_kdtree_balance(psys->tree);
2067
2068                         psys->tree_frame = cfra;
2069                 }
2070         }
2071 }
2072
2073 static void psys_update_effectors(ParticleSimulationData *sim)
2074 {
2075         pdEndEffectors(&sim->psys->effectors);
2076         sim->psys->effectors = pdInitEffectors(sim->scene, sim->ob, sim->psys, sim->psys->part->effector_weights);
2077         precalc_guides(sim, sim->psys->effectors);
2078 }
2079
2080 static void integrate_particle(ParticleSettings *part, ParticleData *pa, float dtime, float *external_acceleration, void (*force_func)(void *forcedata, ParticleKey *state, float *force, float *impulse), void *forcedata)
2081 {
2082         ParticleKey states[5];
2083         float force[3],acceleration[3],impulse[3],dx[4][3],dv[4][3],oldpos[3];
2084         float pa_mass= (part->flag & PART_SIZEMASS ? part->mass * pa->size : part->mass);
2085         int i, steps=1;
2086         int integrator = part->integrator;
2087
2088         copy_v3_v3(oldpos, pa->state.co);
2089
2090         /* Verlet integration behaves strangely with moving emitters, so do first step with euler. */
2091         if(pa->prev_state.time < 0.f && integrator == PART_INT_VERLET)
2092                 integrator = PART_INT_EULER;
2093
2094         switch(integrator){
2095                 case PART_INT_EULER:
2096                         steps=1;
2097                         break;
2098                 case PART_INT_MIDPOINT:
2099                         steps=2;
2100                         break;
2101                 case PART_INT_RK4:
2102                         steps=4;
2103                         break;
2104                 case PART_INT_VERLET:
2105                         steps=1;
2106                         break;
2107         }
2108
2109         copy_particle_key(states, &pa->state, 1);
2110
2111         states->time = 0.f;
2112
2113         for(i=0; i<steps; i++){
2114                 zero_v3(force);
2115                 zero_v3(impulse);
2116
2117                 force_func(forcedata, states+i, force, impulse);
2118
2119                 /* force to acceleration*/
2120                 mul_v3_v3fl(acceleration, force, 1.0f/pa_mass);
2121
2122                 if(external_acceleration)
2123                         add_v3_v3(acceleration, external_acceleration);
2124
2125                 /* calculate next state */
2126                 add_v3_v3(states[i].vel, impulse);
2127
2128                 switch(integrator){
2129                         case PART_INT_EULER:
2130                                 madd_v3_v3v3fl(pa->state.co, states->co, states->vel, dtime);
2131                                 madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
2132                                 break;
2133                         case PART_INT_MIDPOINT:
2134                                 if(i==0){
2135                                         madd_v3_v3v3fl(states[1].co, states->co, states->vel, dtime*0.5f);
2136                                         madd_v3_v3v3fl(states[1].vel, states->vel, acceleration, dtime*0.5f);
2137                                         states[1].time = dtime*0.5f;
2138                                         /*fra=sim->psys->cfra+0.5f*dfra;*/
2139                                 }
2140                                 else{
2141                                         madd_v3_v3v3fl(pa->state.co, states->co, states[1].vel, dtime);
2142                                         madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
2143                                 }
2144                                 break;
2145                         case PART_INT_RK4:
2146                                 switch(i){
2147                                         case 0:
2148                                                 copy_v3_v3(dx[0], states->vel);
2149                                                 mul_v3_fl(dx[0], dtime);
2150                                                 copy_v3_v3(dv[0], acceleration);
2151                                                 mul_v3_fl(dv[0], dtime);
2152
2153                                                 madd_v3_v3v3fl(states[1].co, states->co, dx[0], 0.5f);
2154                                                 madd_v3_v3v3fl(states[1].vel, states->vel, dv[0], 0.5f);
2155                                                 states[1].time = dtime*0.5f;
2156                                                 /*fra=sim->psys->cfra+0.5f*dfra;*/
2157                                                 break;
2158                                         case 1:
2159                                                 madd_v3_v3v3fl(dx[1], states->vel, dv[0], 0.5f);
2160                                                 mul_v3_fl(dx[1], dtime);
2161                                                 copy_v3_v3(dv[1], acceleration);
2162                                                 mul_v3_fl(dv[1], dtime);
2163
2164                                                 madd_v3_v3v3fl(states[2].co, states->co, dx[1], 0.5f);
2165                                                 madd_v3_v3v3fl(states[2].vel, states->vel, dv[1], 0.5f);
2166                                                 states[2].time = dtime*0.5f;
2167                                                 break;
2168                                         case 2:
2169                                                 madd_v3_v3v3fl(dx[2], states->vel, dv[1], 0.5f);
2170                                                 mul_v3_fl(dx[2], dtime);
2171                                                 copy_v3_v3(dv[2], acceleration);
2172                                                 mul_v3_fl(dv[2], dtime);
2173
2174                                                 add_v3_v3v3(states[3].co, states->co, dx[2]);
2175                                                 add_v3_v3v3(states[3].vel, states->vel, dv[2]);
2176                                                 states[3].time = dtime;
2177                                                 /*fra=cfra;*/
2178                                                 break;
2179                                         case 3:
2180                                                 add_v3_v3v3(dx[3], states->vel, dv[2]);
2181                                                 mul_v3_fl(dx[3], dtime);
2182                                                 copy_v3_v3(dv[3], acceleration);
2183                                                 mul_v3_fl(dv[3], dtime);
2184
2185                                                 madd_v3_v3v3fl(pa->state.co, states->co, dx[0], 1.0f/6.0f);
2186                                                 madd_v3_v3fl(pa->state.co, dx[1], 1.0f/3.0f);
2187                                                 madd_v3_v3fl(pa->state.co, dx[2], 1.0f/3.0f);
2188                                                 madd_v3_v3fl(pa->state.co, dx[3], 1.0f/6.0f);
2189
2190                                                 madd_v3_v3v3fl(pa->state.vel, states->vel, dv[0], 1.0f/6.0f);
2191                                                 madd_v3_v3fl(pa->state.vel, dv[1], 1.0f/3.0f);
2192                                                 madd_v3_v3fl(pa->state.vel, dv[2], 1.0f/3.0f);
2193                                                 madd_v3_v3fl(pa->state.vel, dv[3], 1.0f/6.0f);
2194                                 }
2195                                 break;
2196                         case PART_INT_VERLET:   /* Verlet integration */
2197                                 madd_v3_v3v3fl(pa->state.vel, pa->prev_state.vel, acceleration, dtime);
2198                                 madd_v3_v3v3fl(pa->state.co, pa->prev_state.co, pa->state.vel, dtime);
2199
2200                                 sub_v3_v3v3(pa->state.vel, pa->state.co, oldpos);
2201                                 mul_v3_fl(pa->state.vel, 1.0f/dtime);
2202                                 break;
2203                 }
2204         }
2205 }
2206
2207 /*********************************************************************************************************
2208                     SPH fluid physics
2209
2210  In theory, there could be unlimited implementation of SPH simulators
2211
2212  This code uses in some parts adapted algorithms from the pseudo code as outlined in the Research paper:
2213
2214  Titled: Particle-based Viscoelastic Fluid Simulation.
2215  Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
2216  Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
2217
2218  Presented at Siggraph, (2005)
2219
2220 ***********************************************************************************************************/
2221 #define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
2222 static ParticleSpring *sph_spring_add(ParticleSystem *psys, ParticleSpring *spring)
2223 {
2224         /* Are more refs required? */
2225         if(psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
2226                 psys->alloc_fluidsprings = PSYS_FLUID_SPRINGS_INITIAL_SIZE;
2227                 psys->fluid_springs = (ParticleSpring*)MEM_callocN(psys->alloc_fluidsprings * sizeof(ParticleSpring), "Particle Fluid Springs");
2228         }
2229         else if(psys->tot_fluidsprings == psys->alloc_fluidsprings) {
2230                 /* Double the number of refs allocated */
2231                 psys->alloc_fluidsprings *= 2;
2232                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
2233         }
2234
2235         memcpy(psys->fluid_springs + psys->tot_fluidsprings, spring, sizeof(ParticleSpring));
2236         psys->tot_fluidsprings++;
2237
2238         return psys->fluid_springs + psys->tot_fluidsprings - 1;
2239 }
2240 static void sph_spring_delete(ParticleSystem *psys, int j)
2241 {
2242         if (j != psys->tot_fluidsprings - 1)
2243                 psys->fluid_springs[j] = psys->fluid_springs[psys->tot_fluidsprings - 1];
2244
2245         psys->tot_fluidsprings--;
2246
2247         if (psys->tot_fluidsprings < psys->alloc_fluidsprings/2 && psys->alloc_fluidsprings > PSYS_FLUID_SPRINGS_INITIAL_SIZE){
2248                 psys->alloc_fluidsprings /= 2;
2249                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs,  psys->alloc_fluidsprings * sizeof(ParticleSpring));
2250         }
2251 }
2252 static void sph_springs_modify(ParticleSystem *psys, float dtime){
2253         SPHFluidSettings *fluid = psys->part->fluid;
2254         ParticleData *pa1, *pa2;
2255         ParticleSpring *spring = psys->fluid_springs;
2256
2257         float h, d, Rij[3], rij, Lij;
2258         int i;
2259
2260         float yield_ratio = fluid->yield_ratio;
2261         float plasticity = fluid->plasticity_constant;
2262         /* scale things according to dtime */
2263         float timefix = 25.f * dtime;
2264
2265         if((fluid->flag & SPH_VISCOELASTIC_SPRINGS)==0 || fluid->spring_k == 0.f)
2266                 return;
2267
2268         /* Loop through the springs */
2269         for(i=0; i<psys->tot_fluidsprings; i++, spring++) {
2270                 pa1 = psys->particles + spring->particle_index[0];
2271                 pa2 = psys->particles + spring->particle_index[1];
2272
2273                 sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
2274                 rij = normalize_v3(Rij);
2275
2276                 /* adjust rest length */
2277                 Lij = spring->rest_length;
2278                 d = yield_ratio * timefix * Lij;
2279
2280                 if (rij > Lij + d) // Stretch
2281                         spring->rest_length += plasticity * (rij - Lij - d) * timefix;
2282                 else if(rij < Lij - d) // Compress
2283                         spring->rest_length -= plasticity * (Lij - d - rij) * timefix;
2284
2285                 h = 4.f*pa1->size;
2286
2287                 if(spring->rest_length > h)
2288                         spring->delete_flag = 1;
2289         }
2290
2291         /* Loop through springs backwaqrds - for efficient delete function */
2292         for (i=psys->tot_fluidsprings-1; i >= 0; i--) {
2293                 if(psys->fluid_springs[i].delete_flag)
2294                         sph_spring_delete(psys, i);
2295         }
2296 }
2297 static EdgeHash *sph_springhash_build(ParticleSystem *psys)
2298 {
2299         EdgeHash *springhash = NULL;
2300         ParticleSpring *spring;
2301         int i = 0;
2302
2303         springhash = BLI_edgehash_new();
2304
2305         for(i=0, spring=psys->fluid_springs; i<psys->tot_fluidsprings; i++, spring++)
2306                 BLI_edgehash_insert(springhash, spring->particle_index[0], spring->particle_index[1], SET_INT_IN_POINTER(i+1));
2307
2308         return springhash;
2309 }
2310
2311 typedef struct SPHNeighbor
2312 {
2313         ParticleSystem *psys;
2314         int index;
2315 } SPHNeighbor;
2316 typedef struct SPHRangeData
2317 {
2318         SPHNeighbor neighbors[128];
2319         int tot_neighbors;
2320
2321         float density, near_density;
2322         float h;
2323
2324         ParticleSystem *npsys;
2325         ParticleData *pa;
2326
2327         float massfac;
2328         int use_size;
2329 } SPHRangeData;
2330 typedef struct SPHData {
2331         ParticleSystem *psys[10];
2332         ParticleData *pa;
2333         float mass;
2334         EdgeHash *eh;
2335         float *gravity;
2336 }SPHData;
2337 static void sph_density_accum_cb(void *userdata, int index, float squared_dist)
2338 {
2339         SPHRangeData *pfr = (SPHRangeData *)userdata;
2340         ParticleData *npa = pfr->npsys->particles + index;
2341         float q;
2342
2343         if(npa == pfr->pa || squared_dist < FLT_EPSILON)
2344                 return;
2345
2346         /* Ugh! One particle has over 128 neighbors! Really shouldn't happen,
2347          * but even if it does it shouldn't do any terrible harm if all are
2348          * not taken into account - jahka
2349          */
2350         if(pfr->tot_neighbors >= 128)
2351                 return;
2352
2353         pfr->neighbors[pfr->tot_neighbors].index = index;
2354         pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
2355         pfr->tot_neighbors++;
2356
2357         q = (1.f - sqrtf(squared_dist)/pfr->h) * pfr->massfac;
2358
2359         if(pfr->use_size)
2360                 q *= npa->size;
2361
2362         pfr->density += q*q;
2363         pfr->near_density += q*q*q;
2364 }
2365 static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse))
2366 {
2367         SPHData *sphdata = (SPHData *)sphdata_v;
2368         ParticleSystem **psys = sphdata->psys;
2369         ParticleData *pa = sphdata->pa;
2370         SPHFluidSettings *fluid = psys[0]->part->fluid;
2371         ParticleSpring *spring = NULL;
2372         SPHRangeData pfr;
2373         SPHNeighbor *pfn;
2374         float mass = sphdata->mass;
2375         float *gravity = sphdata->gravity;
2376         EdgeHash *springhash = sphdata->eh;
2377
2378         float q, u, rij, dv[3];
2379         float pressure, near_pressure;
2380
2381         float visc = fluid->viscosity_omega;
2382         float stiff_visc = fluid->viscosity_beta * (fluid->flag & SPH_FAC_VISCOSITY ? fluid->viscosity_omega : 1.f);
2383
2384         float inv_mass = 1.0f/mass;
2385         float spring_constant = fluid->spring_k;
2386
2387         float h = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.f*pa->size : 1.f); /* 4.0 seems to be a pretty good value */
2388         float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f); /* 4.77 is an experimentally determined density factor */
2389         float rest_length = fluid->rest_length * (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f);
2390
2391         float stiffness = fluid->stiffness_k;
2392         float stiffness_near_fac = fluid->stiffness_knear * (fluid->flag & SPH_FAC_REPULSION ? fluid->stiffness_k : 1.f);
2393
2394         ParticleData *npa;
2395         float vec[3];
2396         float vel[3];
2397         float co[3];
2398
2399         int i, spring_index, index = pa - psys[0]->particles;
2400
2401         pfr.tot_neighbors = 0;
2402         pfr.density = pfr.near_density = 0.f;
2403         pfr.h = h;
2404         pfr.pa = pa;
2405
2406         for(i=0; i<10 && psys[i]; i++) {
2407                 pfr.npsys = psys[i];
2408                 pfr.massfac = psys[i]->part->mass*inv_mass;
2409                 pfr.use_size = psys[i]->part->flag & PART_SIZEMASS;
2410
2411                 BLI_bvhtree_range_query(psys[i]->bvhtree, state->co, h, sph_density_accum_cb, &pfr);
2412         }
2413
2414         pressure =  stiffness * (pfr.density - rest_density);
2415         near_pressure = stiffness_near_fac * pfr.near_density;
2416
2417         pfn = pfr.neighbors;
2418         for(i=0; i<pfr.tot_neighbors; i++, pfn++) {
2419                 npa = pfn->psys->particles + pfn->index;
2420
2421                 madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
2422
2423                 sub_v3_v3v3(vec, co, state->co);
2424                 rij = normalize_v3(vec);
2425
2426                 q = (1.f - rij/h) * pfn->psys->part->mass * inv_mass;
2427
2428                 if(pfn->psys->part->flag & PART_SIZEMASS)
2429                         q *= npa->size;
2430
2431                 copy_v3_v3(vel, npa->prev_state.vel);
2432
2433                 /* Double Density Relaxation */
2434                 madd_v3_v3fl(force, vec, -(pressure + near_pressure*q)*q);
2435
2436                 /* Viscosity */
2437                 if(visc > 0.f   || stiff_visc > 0.f) {
2438                         sub_v3_v3v3(dv, vel, state->vel);
2439                         u = dot_v3v3(vec, dv);
2440
2441                         if(u < 0.f && visc > 0.f)
2442                                 madd_v3_v3fl(force, vec, 0.5f * q * visc * u );
2443
2444                         if(u > 0.f && stiff_visc > 0.f)
2445                                 madd_v3_v3fl(force, vec, 0.5f * q * stiff_visc * u );
2446                 }
2447
2448                 if(spring_constant > 0.f) {
2449                         /* Viscoelastic spring force */
2450                         if (pfn->psys == psys[0] && fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash) {
2451                                 spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, pfn->index));
2452
2453                                 if(spring_index) {
2454                                         spring = psys[0]->fluid_springs + spring_index - 1;
2455
2456                                         madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (spring->rest_length - rij));
2457                                 }
2458                                 else if(fluid->spring_frames == 0 || (pa->prev_state.time-pa->time) <= fluid->spring_frames){
2459                                         ParticleSpring temp_spring;
2460                                         temp_spring.particle_index[0] = index;
2461                                         temp_spring.particle_index[1] = pfn->index;
2462                                         temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : rest_length;
2463                                         temp_spring.delete_flag = 0;
2464
2465                                         sph_spring_add(psys[0], &temp_spring);
2466                                 }
2467                         }
2468                         else {/* PART_SPRING_HOOKES - Hooke's spring force */
2469                                 madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (rest_length - rij));
2470                         }
2471                 }
2472         }
2473
2474         /* Artificial buoyancy force in negative gravity direction  */
2475         if (fluid->buoyancy > 0.f && gravity)
2476                 madd_v3_v3fl(force, gravity, fluid->buoyancy * (pfr.density-rest_density));
2477 }
2478
2479 static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, float *gravity, EdgeHash *springhash){
2480         ParticleTarget *pt;
2481         int i;
2482
2483         ParticleSettings *part = sim->psys->part;
2484         // float timestep = psys_get_timestep(sim); // UNUSED
2485         float pa_mass = part->mass * (part->flag & PART_SIZEMASS ? pa->size : 1.f);
2486         float dtime = dfra*psys_get_timestep(sim);
2487         // int steps = 1; // UNUSED
2488         float effector_acceleration[3];
2489         SPHData sphdata;
2490
2491         sphdata.psys[0] = sim->psys;
2492         for(i=1, pt=sim->psys->targets.first; i<10; i++, pt=(pt?pt->next:NULL))
2493                 sphdata.psys[i] = pt ? psys_get_target_system(sim->ob, pt) : NULL;
2494
2495         sphdata.pa = pa;
2496         sphdata.gravity = gravity;
2497         sphdata.mass = pa_mass;
2498         sphdata.eh = springhash;
2499
2500         /* restore previous state and treat gravity & effectors as external acceleration*/
2501         sub_v3_v3v3(effector_acceleration, pa->state.vel, pa->prev_state.vel);
2502         mul_v3_fl(effector_acceleration, 1.f/dtime);
2503
2504         copy_particle_key(&pa->state, &pa->prev_state, 0);
2505
2506         integrate_particle(part, pa, dtime, effector_acceleration, sph_force_cb, &sphdata);
2507 }
2508
2509 /************************************************/
2510 /*                      Basic physics                                           */
2511 /************************************************/
2512 typedef struct EfData
2513 {
2514         ParticleTexture ptex;
2515         ParticleSimulationData *sim;
2516         ParticleData *pa;
2517 } EfData;
2518 static void basic_force_cb(void *efdata_v, ParticleKey *state, float *force, float *impulse)
2519 {
2520         EfData *efdata = (EfData *)efdata_v;
2521         ParticleSimulationData *sim = efdata->sim;
2522         ParticleSettings *part = sim->psys->part;
2523         ParticleData *pa = efdata->pa;
2524         EffectedPoint epoint;
2525
2526         /* add effectors */
2527         pd_point_from_particle(efdata->sim, efdata->pa, state, &epoint);
2528         if(part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)
2529                 pdDoEffectors(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);
2530
2531         mul_v3_fl(force, efdata->ptex.field);
2532         mul_v3_fl(impulse, efdata->ptex.field);
2533
2534         /* calculate air-particle interaction */
2535         if(part->dragfac != 0.0f)
2536                 madd_v3_v3fl(force, state->vel, -part->dragfac * pa->size * pa->size * len_v3(state->vel));
2537
2538         /* brownian force */
2539         if(part->brownfac != 0.0f){
2540                 force[0] += (BLI_frand()-0.5f) * part->brownfac;
2541                 force[1] += (BLI_frand()-0.5f) * part->brownfac;
2542                 force[2] += (BLI_frand()-0.5f) * part->brownfac;
2543         }
2544 }
2545 /* gathers all forces that effect particles and calculates a new state for the particle */
2546 static void basic_integrate(ParticleSimulationData *sim, int p, float dfra, float cfra)
2547 {
2548         ParticleSettings *part = sim->psys->part;
2549         ParticleData *pa = sim->psys->particles + p;
2550         ParticleKey tkey;
2551         float dtime=dfra*psys_get_timestep(sim), time;
2552         float *gravity = NULL, gr[3];
2553         EfData efdata;
2554
2555         psys_get_texture(sim, pa, &efdata.ptex, PAMAP_PHYSICS, cfra);
2556
2557         efdata.pa = pa;
2558         efdata.sim = sim;
2559
2560         /* add global acceleration (gravitation) */
2561         if(psys_uses_gravity(sim)
2562                 /* normal gravity is too strong for hair so it's disabled by default */
2563                 && (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
2564                 zero_v3(gr);
2565                 madd_v3_v3fl(gr, sim->scene->physics_settings.gravity, part->effector_weights->global_gravity * efdata.ptex.gravity);
2566                 gravity = gr;
2567         }
2568
2569         /* maintain angular velocity */
2570         copy_v3_v3(pa->state.ave, pa->prev_state.ave);
2571
2572         integrate_particle(part, pa, dtime, gravity, basic_force_cb, &efdata);
2573
2574         /* damp affects final velocity */
2575         if(part->dampfac != 0.f)
2576                 mul_v3_fl(pa->state.vel, 1.f - part->dampfac * efdata.ptex.damp * 25.f * dtime);
2577
2578         //VECCOPY(pa->state.ave, states->ave);
2579
2580         /* finally we do guides */
2581         time=(cfra-pa->time)/pa->lifetime;
2582         CLAMP(time, 0.0f, 1.0f);
2583
2584         VECCOPY(tkey.co,pa->state.co);
2585         VECCOPY(tkey.vel,pa->state.vel);
2586         tkey.time=pa->state.time;
2587
2588         if(part->type != PART_HAIR) {
2589                 if(do_guides(sim->psys->effectors, &tkey, p, time)) {
2590                         VECCOPY(pa->state.co,tkey.co);
2591                         /* guides don't produce valid velocity */
2592                         VECSUB(pa->state.vel,tkey.co,pa->prev_state.co);
2593                         mul_v3_fl(pa->state.vel,1.0f/dtime);
2594                         pa->state.time=tkey.time;
2595                 }
2596         }
2597 }
2598 static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
2599 {
2600         float rotfac, rot1[4], rot2[4]={1.0,0.0,0.0,0.0}, dtime=dfra*timestep;
2601
2602         if((part->flag & PART_ROT_DYN)==0){
2603                 if(part->avemode==PART_AVE_SPIN){
2604                         float angle;
2605                         float len1 = len_v3(pa->prev_state.vel);
2606                         float len2 = len_v3(pa->state.vel);
2607
2608                         if(len1==0.0f || len2==0.0f)
2609                                 pa->state.ave[0]=pa->state.ave[1]=pa->state.ave[2]=0.0f;
2610                         else{
2611                                 cross_v3_v3v3(pa->state.ave,pa->prev_state.vel,pa->state.vel);
2612                                 normalize_v3(pa->state.ave);
2613                                 angle=dot_v3v3(pa->prev_state.vel,pa->state.vel)/(len1*len2);
2614                                 mul_v3_fl(pa->state.ave,saacos(angle)/dtime);
2615                         }
2616
2617                         axis_angle_to_quat(rot2,pa->state.vel,dtime*part->avefac);
2618                 }
2619         }
2620
2621         rotfac=len_v3(pa->state.ave);
2622         if(rotfac == 0.0f){ /* unit_qt(in VecRotToQuat) doesn't give unit quat [1,0,0,0]?? */
2623                 rot1[0]=1.0f;
2624                 rot1[1]=rot1[2]=rot1[3]=0;
2625         }
2626         else{
2627                 axis_angle_to_quat(rot1,pa->state.ave,rotfac*dtime);
2628         }
2629         mul_qt_qtqt(pa->state.rot,rot1,pa->prev_state.rot);
2630         mul_qt_qtqt(pa->state.rot,rot2,pa->state.rot);
2631
2632         /* keep rotation quat in good health */
2633         normalize_qt(pa->state.rot);
2634 }
2635
2636 /************************************************/
2637 /*                      Collisions                                                      */
2638 /************************************************/
2639 #define COLLISION_MAX_COLLISIONS        10
2640 #define COLLISION_MIN_RADIUS 0.001f
2641 #define COLLISION_MIN_DISTANCE 0.0001f
2642 #define COLLISION_ZERO 0.00001f
2643 typedef float (*NRDistanceFunc)(float *p, float radius, ParticleCollisionElement *pce, float *nor);
2644 static float nr_signed_distance_to_plane(float *p, float radius, ParticleCollisionElement *pce, float *nor)
2645 {
2646         float p0[3], e1[3], e2[3], d;
2647
2648         sub_v3_v3v3(e1, pce->x1, pce->x0);
2649         sub_v3_v3v3(e2, pce->x2, pce->x0);
2650         sub_v3_v3v3(p0, p, pce->x0);
2651
2652         cross_v3_v3v3(nor, e1, e2);
2653         normalize_v3(nor);
2654
2655         d = dot_v3v3(p0, nor);
2656
2657         if(pce->inv_nor == -1) {
2658                 if(d < 0.f)
2659                         pce->inv_nor = 1;
2660                 else
2661                         pce->inv_nor = 0;
2662         }
2663
2664         if(pce->inv_nor == 1) {
2665                 mul_v3_fl(nor, -1.f);
2666                 d = -d;
2667         }
2668
2669         return d - radius;
2670 }
2671 static float nr_distance_to_edge(float *p, float radius, ParticleCollisionElement *pce, float *UNUSED(nor))
2672 {
2673         float v0[3], v1[3], v2[3], c[3];
2674
2675         sub_v3_v3v3(v0, pce->x1, pce->x0);
2676         sub_v3_v3v3(v1, p, pce->x0);
2677         sub_v3_v3v3(v2, p, pce->x1);
2678
2679         cross_v3_v3v3(c, v1, v2);
2680
2681         return fabsf(len_v3(c)/len_v3(v0)) - radius;
2682 }
2683 static float nr_distance_to_vert(float *p, float radius, ParticleCollisionElement *pce, float *UNUSED(nor))
2684 {
2685         return len_v3v3(p, pce->x0) - radius;
2686 }
2687 static void collision_interpolate_element(ParticleCollisionElement *pce, float t, float fac, ParticleCollision *col)
2688 {
2689         /* t is the current time for newton rhapson */
2690         /* fac is the starting factor for current collision iteration */
2691         /* the col->fac's are factors for the particle subframe step start and end during collision modifier step */
2692         float f = fac + t*(1.f-fac);
2693         float mul = col->fac1 + f * (col->fac2-col->fac1);
2694         if(pce->tot > 0) {
2695                 madd_v3_v3v3fl(pce->x0, pce->x[0], pce->v[0], mul);
2696
2697                 if(pce->tot > 1) {
2698                         madd_v3_v3v3fl(pce->x1, pce->x[1], pce->v[1], mul);
2699
2700                         if(pce->tot > 2)
2701                                 madd_v3_v3v3fl(pce->x2, pce->x[2], pce->v[2], mul);
2702                 }
2703         }
2704 }
2705 static void collision_point_velocity(ParticleCollisionElement *pce)
2706 {
2707         float v[3];
2708
2709         copy_v3_v3(pce->vel, pce->v[0]);
2710
2711         if(pce->tot > 1) {
2712                 sub_v3_v3v3(v, pce->v[1], pce->v[0]);
2713                 madd_v3_v3fl(pce->vel, v, pce->uv[0]);
2714
2715                 if(pce->tot > 2) {
2716                         sub_v3_v3v3(v, pce->v[2], pce->v[0]);
2717                         madd_v3_v3fl(pce->vel, v, pce->uv[1]);
2718                 }
2719         }
2720 }
2721 static float collision_point_distance_with_normal(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *nor)
2722 {
2723         if(fac >= 0.f)
2724                 collision_interpolate_element(pce, 0.f, fac, col);
2725
2726         switch(pce->tot) {
2727                 case 1:
2728                 {
2729                         sub_v3_v3v3(nor, p, pce->x0);
2730                         return normalize_v3(nor);
2731                 }
2732                 case 2:
2733                 {
2734                         float u, e[3], vec[3];
2735                         sub_v3_v3v3(e, pce->x1, pce->x0);
2736                         sub_v3_v3v3(vec, p, pce->x0);
2737                         u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2738
2739                         madd_v3_v3v3fl(nor, vec, e, -u);
2740                         return normalize_v3(nor);
2741                 }
2742                 case 3:
2743                         return nr_signed_distance_to_plane(p, 0.f, pce, nor);
2744         }
2745         return 0;
2746 }
2747 static void collision_point_on_surface(float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *co)
2748 {
2749         collision_interpolate_element(pce, 0.f, fac, col);
2750
2751         switch(pce->tot) {
2752                 case 1:
2753                 {
2754                         sub_v3_v3v3(co, p, pce->x0);
2755                         normalize_v3(co);
2756                         madd_v3_v3v3fl(co, pce->x0, co, col->radius);
2757                         break;
2758                 }
2759                 case 2:
2760                 {
2761                         float u, e[3], vec[3], nor[3];
2762                         sub_v3_v3v3(e, pce->x1, pce->x0);
2763                         sub_v3_v3v3(vec, p, pce->x0);
2764                         u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2765
2766                         madd_v3_v3v3fl(nor, vec, e, -u);
2767                         normalize_v3(nor);
2768
2769                         madd_v3_v3v3fl(co, pce->x0, e, pce->uv[0]);
2770                         madd_v3_v3fl(co, nor, col->radius);
2771                         break;
2772                 }
2773                 case 3:
2774                 {
2775                                 float p0[3], e1[3], e2[3], nor[3];
2776
2777                                 sub_v3_v3v3(e1, pce->x1, pce->x0);
2778                                 sub_v3_v3v3(e2, pce->x2, pce->x0);
2779                                 sub_v3_v3v3(p0, p, pce->x0);
2780
2781                                 cross_v3_v3v3(nor, e1, e2);
2782                                 normalize_v3(nor);
2783
2784                                 if(pce->inv_nor == 1)
2785                                         mul_v3_fl(nor, -1.f);
2786
2787                                 madd_v3_v3v3fl(co, pce->x0, nor, col->radius);
2788                                 madd_v3_v3fl(co, e1, pce->uv[0]);
2789                                 madd_v3_v3fl(co, e2, pce->uv[1]);
2790                         break;
2791                 }
2792         }
2793 }
2794 /* find first root in range [0-1] starting from 0 */
2795 static float collision_newton_rhapson(ParticleCollision *col, float radius, ParticleCollisionElement *pce, NRDistanceFunc distance_func)
2796 {
2797         float t0, t1, d0, d1, dd, n[3];
2798         int iter;
2799
2800         pce->inv_nor = -1;
2801
2802         /* start from the beginning */
2803         t0 = 0.f;
2804         collision_interpolate_element(pce, t0, col->f, col);
2805         d0 = distance_func(col->co1, radius, pce, n);
2806         t1 = 0.001f;
2807         d1 = 0.f;
2808
2809         for(iter=0; iter<10; iter++) {//, itersum++) {
2810                 /* get current location */
2811                 collision_interpolate_element(pce, t1, col->f, col);
2812                 interp_v3_v3v3(pce->p, col->co1, col->co2, t1);
2813
2814                 d1 = distance_func(pce->p, radius, pce, n);
2815
2816                 /* no movement, so no collision */
2817                 if(d1 == d0) {
2818                         return -1.f;
2819                 }
2820
2821                 /* particle already inside face, so report collision */
2822                 if(iter == 0 && d0 < 0.f && d0 > -radius) {
2823                         copy_v3_v3(pce->p, col->co1);
2824                         copy_v3_v3(pce->nor, n);
2825                         pce->inside = 1;
2826                         return 0.f;
2827                 }
2828
2829                 dd = (t1-t0)/(d1-d0);
2830
2831                 t0 = t1;
2832                 d0 = d1;
2833
2834                 t1 -= d1*dd;
2835
2836                 /* particle movin away from plane could also mean a strangely rotating face, so check from end */
2837                 if(iter == 0 && t1 < 0.f) {
2838                         t0 = 1.f;
2839                         collision_interpolate_element(pce, t0, col->f, col);
2840                         d0 = distance_func(col->co2, radius, pce, n);
2841                         t1 = 0.999f;
2842                         d1 = 0.f;
2843
2844                         continue;
2845                 }
2846                 else if(iter == 1 && (t1 < -COLLISION_ZERO || t1 > 1.f))
2847                         return -1.f;
2848
2849                 if(d1 <= COLLISION_ZERO && d1 >= -COLLISION_ZERO) {
2850                         if(t1 >= -COLLISION_ZERO && t1 <= 1.f) {
2851                                 if(distance_func == nr_signed_distance_to_plane)
2852                                         copy_v3_v3(pce->nor, n);
2853
2854                                 CLAMP(t1, 0.f, 1.f);
2855
2856                                 return t1;
2857                         }
2858                         else
2859                                 return -1.f;
2860                 }
2861         }
2862         return -1.0;
2863 }
2864 static int collision_sphere_to_tri(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
2865 {
2866         ParticleCollisionElement *result = &col->pce;
2867         float ct, u, v;
2868
2869         pce->inv_nor = -1;
2870         pce->inside = 0;
2871
2872         ct = collision_newton_rhapson(col, radius, pce, nr_signed_distance_to_plane);
2873
2874         if(ct >= 0.f && ct < *t && (result->inside==0 || pce->inside==1) ) {
2875                 float e1[3], e2[3], p0[3];
2876                 float e1e1, e1e2, e1p0, e2e2, e2p0, inv;
2877
2878                 sub_v3_v3v3(e1, pce->x1, pce->x0);
2879                 sub_v3_v3v3(e2, pce->x2, pce->x0);
2880                 /* XXX: add radius correction here? */
2881                 sub_v3_v3v3(p0, pce->p, pce->x0);
2882
2883                 e1e1 = dot_v3v3(e1, e1);
2884                 e1e2 = dot_v3v3(e1, e2);
2885                 e1p0 = dot_v3v3(e1, p0);
2886                 e2e2 = dot_v3v3(e2, e2);
2887                 e2p0 = dot_v3v3(e2, p0);
2888
2889                 inv = 1.f/(e1e1 * e2e2 - e1e2 * e1e2);
2890                 u = (e2e2 * e1p0 - e1e2 * e2p0) * inv;
2891                 v = (e1e1 * e2p0 - e1e2 * e1p0) * inv;
2892
2893                 if(u>=0.f && u<=1.f && v>=0.f && u+v<=1.f) {
2894                         *result = *pce;
2895
2896                         /* normal already calculated in pce */
2897
2898                         result->uv[0] = u;
2899                         result->uv[1] = v;
2900
2901                         *t = ct;
2902                         return 1;
2903                 }
2904         }
2905         return 0;
2906 }
2907 static int collision_sphere_to_edges(ParticleCollision *col, float radius, ParticleCollisionElement *pce, float *t)
2908 {
2909         ParticleCollisionElement edge[3], *cur = NULL, *hit = NULL;
2910         ParticleCollisionElement *result = &col->pce;
2911
2912         float ct;
2913         int i;
2914
2915         for(i=0; i<3; i++) {
2916                 /* in case of a quad, no need to check "edge" that goes through face twice */
2917                 if((pce->x[3] && i==2))
2918                         continue;
2919
2920                 cur = edge+i;
2921                 cur->x[0] = pce->x[i]; cur->x[1] = pce->x[(i+1)%3];
2922                 cur->v[0] = pce->v[i]; cur->v[1] = pce->v[(i+1)%3];
2923                 cur->tot = 2;
2924                 cur->inside = 0;
2925
2926                 ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_edge);
2927
2928                 if(ct >= 0.f && ct < *t) {
2929                         float u, e[3], vec[3];
2930
2931                         sub_v3_v3v3(e, cur->x1, cur->x0);
2932                         sub_v3_v3v3(vec, cur->p, cur->x0);
2933                         u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2934
2935                         if(u < 0.f || u > 1.f)
2936                                 break;
2937
2938                         *result = *cur;
2939
2940                         madd_v3_v3v3fl(result->nor, vec, e, -u);
2941                         normalize_v3(result->nor);
2942
2943                         result->uv[0] = u;
2944
2945
2946                         hit = cur;
2947                         *t = ct;
2948                 }
2949
2950         }
2951
2952         return hit != NULL;