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