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