source/blender/blenkernel/intern/particle_system.c

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