2 * ***** BEGIN GPL LICENSE BLOCK *****
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.
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.
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.
18 * The Original Code is Copyright (C) 2007 by Janne Karhu.
19 * All rights reserved.
21 * The Original Code is: all of this file.
23 * Contributor(s): Raul Fernandez Hernandez (Farsthary), Stephen Swhitehorn.
27 * Copyright 2011-2012 AutoCRC
29 * ***** END GPL LICENSE BLOCK *****
32 /** \file blender/blenkernel/intern/particle_system.c
47 #include "MEM_guardedalloc.h"
49 #include "DNA_anim_types.h"
50 #include "DNA_boid_types.h"
51 #include "DNA_particle_types.h"
52 #include "DNA_mesh_types.h"
53 #include "DNA_meshdata_types.h"
54 #include "DNA_modifier_types.h"
55 #include "DNA_object_force.h"
56 #include "DNA_object_types.h"
57 #include "DNA_material_types.h"
58 #include "DNA_curve_types.h"
59 #include "DNA_group_types.h"
60 #include "DNA_scene_types.h"
61 #include "DNA_texture_types.h"
62 #include "DNA_ipo_types.h" // XXX old animation system stuff... to be removed!
63 #include "DNA_listBase.h"
65 #include "BLI_edgehash.h"
67 #include "BLI_jitter.h"
69 #include "BLI_blenlib.h"
70 #include "BLI_kdtree.h"
71 #include "BLI_kdopbvh.h"
72 #include "BLI_threads.h"
73 #include "BLI_utildefines.h"
74 #include "BLI_linklist.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"
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"
94 #include "BKE_modifier.h"
95 #include "BKE_scene.h"
96 #include "BKE_bvhutils.h"
100 #include "RE_shader_ext.h"
102 /* fluid sim particle import */
103 #ifdef WITH_MOD_FLUID
104 #include "DNA_object_fluidsim.h"
105 #include "LBM_fluidsim.h"
109 #endif // WITH_MOD_FLUID
111 /************************************************/
112 /* Reacting to system events */
113 /************************************************/
115 static int particles_are_dynamic(ParticleSystem *psys)
117 if (psys->pointcache->flag & PTCACHE_BAKED)
120 if (psys->part->type == PART_HAIR)
121 return psys->flag & PSYS_HAIR_DYNAMICS;
123 return ELEM3(psys->part->phystype, PART_PHYS_NEWTON, PART_PHYS_BOIDS, PART_PHYS_FLUID);
126 static int psys_get_current_display_percentage(ParticleSystem *psys)
128 ParticleSettings *part=psys->part;
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 */
137 return psys->part->disp;
140 static int tot_particles(ParticleSystem *psys, PTCacheID *pid)
142 if (pid && psys->pointcache->flag & PTCACHE_EXTERNAL)
143 return pid->cache->totpoint;
144 else if (psys->part->distr == PART_DISTR_GRID && psys->part->from != PART_FROM_VERT)
145 return psys->part->grid_res * psys->part->grid_res * psys->part->grid_res - psys->totunexist;
147 return psys->part->totpart - psys->totunexist;
150 void psys_reset(ParticleSystem *psys, int mode)
154 if (ELEM(mode, PSYS_RESET_ALL, PSYS_RESET_DEPSGRAPH)) {
155 if (mode == PSYS_RESET_ALL || !(psys->flag & PSYS_EDITED)) {
156 /* don't free if not absolutely necessary */
157 if (psys->totpart != tot_particles(psys, NULL)) {
158 psys_free_particles(psys);
163 psys->flag &= ~(PSYS_HAIR_DONE|PSYS_KEYED);
165 if (psys->edit && psys->free_edit) {
166 psys->free_edit(psys->edit);
168 psys->free_edit = NULL;
172 else if (mode == PSYS_RESET_CACHE_MISS) {
173 /* set all particles to be skipped */
175 pa->flag |= PARS_NO_DISP;
180 MEM_freeN(psys->child);
186 /* reset path cache */
187 psys_free_path_cache(psys, psys->edit);
189 /* reset point cache */
190 BKE_ptcache_invalidate(psys->pointcache);
192 if (psys->fluid_springs) {
193 MEM_freeN(psys->fluid_springs);
194 psys->fluid_springs = NULL;
197 psys->tot_fluidsprings = psys->alloc_fluidsprings = 0;
200 static void realloc_particles(ParticleSimulationData *sim, int new_totpart)
202 ParticleSystem *psys = sim->psys;
203 ParticleSettings *part = psys->part;
204 ParticleData *newpars = NULL;
205 BoidParticle *newboids = NULL;
207 int totpart, totsaved = 0;
210 if ((part->distr == PART_DISTR_GRID) && (part->from != PART_FROM_VERT)) {
211 totpart= part->grid_res;
212 totpart*=totpart*totpart;
215 totpart=part->totpart;
220 if (totpart != psys->totpart) {
221 if (psys->edit && psys->free_edit) {
222 psys->free_edit(psys->edit);
224 psys->free_edit = NULL;
228 newpars= MEM_callocN(totpart*sizeof(ParticleData), "particles");
232 if (psys->part->phystype == PART_PHYS_BOIDS) {
233 newboids= MEM_callocN(totpart*sizeof(BoidParticle), "boid particles");
235 if (newboids == NULL) {
236 /* allocation error! */
244 if (psys->particles) {
245 totsaved=MIN2(psys->totpart,totpart);
248 memcpy(newpars,psys->particles,totsaved*sizeof(ParticleData));
250 if (psys->particles->boid)
251 memcpy(newboids, psys->particles->boid, totsaved*sizeof(BoidParticle));
254 if (psys->particles->keys)
255 MEM_freeN(psys->particles->keys);
257 if (psys->particles->boid)
258 MEM_freeN(psys->particles->boid);
260 for (p=0, pa=newpars; p<totsaved; p++, pa++) {
267 for (p=totsaved, pa=psys->particles+totsaved; p<psys->totpart; p++, pa++)
268 if (pa->hair) MEM_freeN(pa->hair);
270 MEM_freeN(psys->particles);
274 psys->particles=newpars;
275 psys->totpart=totpart;
279 pa->boid = newboids++;
284 MEM_freeN(psys->child);
290 static int get_psys_child_number(struct Scene *scene, ParticleSystem *psys)
294 if (!psys->part->childtype)
297 if (psys->renderdata)
298 nbr= psys->part->ren_child_nbr;
300 nbr= psys->part->child_nbr;
302 return get_render_child_particle_number(&scene->r, nbr);
305 static int get_psys_tot_child(struct Scene *scene, ParticleSystem *psys)
307 return psys->totpart*get_psys_child_number(scene, psys);
310 static void alloc_child_particles(ParticleSystem *psys, int tot)
313 /* only re-allocate if we have to */
314 if (psys->part->childtype && psys->totchild == tot) {
315 memset(psys->child, 0, tot*sizeof(ChildParticle));
319 MEM_freeN(psys->child);
324 if (psys->part->childtype) {
327 psys->child= MEM_callocN(psys->totchild*sizeof(ChildParticle), "child_particles");
331 /************************************************/
333 /************************************************/
335 void psys_calc_dmcache(Object *ob, DerivedMesh *dm, ParticleSystem *psys)
337 /* use for building derived mesh mapping info:
339 * node: the allocated links - total derived mesh element count
340 * nodearray: the array of nodes aligned with the base mesh's elements, so
341 * each original elements can reference its derived elements
343 Mesh *me= (Mesh*)ob->data;
346 /* CACHE LOCATIONS */
347 if (!dm->deformedOnly) {
348 /* Will use later to speed up subsurf/derivedmesh */
349 LinkNode *node, *nodedmelem, **nodearray;
350 int totdmelem, totelem, i, *origindex, *origindex_poly = NULL;
352 if (psys->part->from == PART_FROM_VERT) {
353 totdmelem= dm->getNumVerts(dm);
354 totelem= me->totvert;
355 origindex= dm->getVertDataArray(dm, CD_ORIGINDEX);
357 else { /* FROM_FACE/FROM_VOLUME */
358 totdmelem= dm->getNumTessFaces(dm);
359 totelem= me->totpoly;
360 origindex = dm->getTessFaceDataArray(dm, CD_ORIGINDEX);
361 /* for face lookups we need the poly origindex too */
362 origindex_poly = dm->getPolyDataArray(dm, CD_ORIGINDEX);
363 if (origindex_poly == NULL) {
368 nodedmelem= MEM_callocN(sizeof(LinkNode)*totdmelem, "psys node elems");
369 nodearray= MEM_callocN(sizeof(LinkNode *)*totelem, "psys node array");
371 for (i=0, node=nodedmelem; i<totdmelem; i++, node++) {
373 node->link = SET_INT_IN_POINTER(i);
375 /* may be vertex or face origindex */
376 origindex_final = origindex ? origindex[i] : ORIGINDEX_NONE;
378 /* if we have a poly source, do an index lookup */
379 if (origindex_poly && origindex_final != ORIGINDEX_NONE) {
380 origindex_final = origindex_poly[origindex_final];
383 if (origindex_final != ORIGINDEX_NONE) {
384 if (nodearray[origindex_final]) {
386 node->next = nodearray[origindex_final];
387 nodearray[origindex_final] = node;
390 nodearray[origindex_final] = node;
395 /* cache the verts/faces! */
398 pa->num_dmcache = -1;
402 if (psys->part->from == PART_FROM_VERT) {
403 if (nodearray[pa->num])
404 pa->num_dmcache= GET_INT_FROM_POINTER(nodearray[pa->num]->link);
406 else { /* FROM_FACE/FROM_VOLUME */
407 /* Note that sometimes the pa->num is over the nodearray size, this is bad, maybe there is a better place to fix this,
408 * but for now passing NULL is OK. every face will be searched for the particle so its slower - Campbell */
409 pa->num_dmcache= psys_particle_dm_face_lookup(ob, dm, pa->num, pa->fuv, pa->num < totelem ? nodearray[pa->num] : NULL);
413 MEM_freeN(nodearray);
414 MEM_freeN(nodedmelem);
417 /* TODO PARTICLE, make the following line unnecessary, each function
418 * should know to use the num or num_dmcache, set the num_dmcache to
419 * an invalid value, just in case */
422 pa->num_dmcache = -1;
426 static void distribute_simple_children(Scene *scene, Object *ob, DerivedMesh *finaldm, ParticleSystem *psys)
428 ChildParticle *cpa = NULL;
430 int child_nbr= get_psys_child_number(scene, psys);
431 int totpart= get_psys_tot_child(scene, psys);
433 alloc_child_particles(psys, totpart);
436 for (i=0; i<child_nbr; i++) {
437 for (p=0; p<psys->totpart; p++,cpa++) {
441 /* create even spherical distribution inside unit sphere */
442 while (length>=1.0f) {
443 cpa->fuv[0]=2.0f*BLI_frand()-1.0f;
444 cpa->fuv[1]=2.0f*BLI_frand()-1.0f;
445 cpa->fuv[2]=2.0f*BLI_frand()-1.0f;
446 length=len_v3(cpa->fuv);
452 /* dmcache must be updated for parent particles if children from faces is used */
453 psys_calc_dmcache(ob, finaldm, psys);
455 static void distribute_grid(DerivedMesh *dm, ParticleSystem *psys)
457 ParticleData *pa=NULL;
458 float min[3], max[3], delta[3], d;
459 MVert *mv, *mvert = dm->getVertDataArray(dm,0);
460 int totvert=dm->getNumVerts(dm), from=psys->part->from;
461 int i, j, k, p, res=psys->part->grid_res, size[3], axis;
465 /* find bounding box of dm */
466 copy_v3_v3(min, mv->co);
467 copy_v3_v3(max, mv->co);
470 for (i=1; i<totvert; i++, mv++) {
471 minmax_v3v3_v3(min, max, mv->co);
474 sub_v3_v3v3(delta, max, min);
476 /* determine major axis */
477 axis = (delta[0]>=delta[1]) ? 0 : ((delta[1]>=delta[2]) ? 1 : 2);
479 d = delta[axis]/(float)res;
482 size[(axis+1)%3] = (int)ceil(delta[(axis+1)%3]/d);
483 size[(axis+2)%3] = (int)ceil(delta[(axis+2)%3]/d);
485 /* float errors grrr.. */
486 size[(axis+1)%3] = MIN2(size[(axis+1)%3],res);
487 size[(axis+2)%3] = MIN2(size[(axis+2)%3],res);
489 size[0] = MAX2(size[0], 1);
490 size[1] = MAX2(size[1], 1);
491 size[2] = MAX2(size[2], 1);
493 /* no full offset for flat/thin objects */
494 min[0]+= d < delta[0] ? d/2.f : delta[0]/2.f;
495 min[1]+= d < delta[1] ? d/2.f : delta[1]/2.f;
496 min[2]+= d < delta[2] ? d/2.f : delta[2]/2.f;
498 for (i=0,p=0,pa=psys->particles; i<res; i++) {
499 for (j=0; j<res; j++) {
500 for (k=0; k<res; k++,p++,pa++) {
501 pa->fuv[0] = min[0] + (float)i*d;
502 pa->fuv[1] = min[1] + (float)j*d;
503 pa->fuv[2] = min[2] + (float)k*d;
504 pa->flag |= PARS_UNEXIST;
505 pa->hair_index = 0; /* abused in volume calculation */
510 /* enable particles near verts/edges/faces/inside surface */
511 if (from==PART_FROM_VERT) {
520 for (i=0,mv=mvert; i<totvert; i++,mv++) {
521 sub_v3_v3v3(vec,mv->co,min);
525 (pa + ((int)(vec[0] * (size[0] - 1)) * res +
526 (int)(vec[1] * (size[1] - 1))) * res +
527 (int)(vec[2] * (size[2] - 1)))->flag &= ~PARS_UNEXIST;
530 else if (ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
531 float co1[3], co2[3];
533 MFace *mface= NULL, *mface_array;
534 float v1[3], v2[3], v3[3], v4[4], lambda;
535 int a, a1, a2, a0mul, a1mul, a2mul, totface;
536 int amax= from==PART_FROM_FACE ? 3 : 1;
538 totface=dm->getNumTessFaces(dm);
539 mface=mface_array=dm->getTessFaceDataArray(dm,CD_MFACE);
541 for (a=0; a<amax; a++) {
542 if (a==0) { a0mul=res*res; a1mul=res; a2mul=1; }
543 else if (a==1) { a0mul=res; a1mul=1; a2mul=res*res; }
544 else { a0mul=1; a1mul=res*res; a2mul=res; }
546 for (a1=0; a1<size[(a+1)%3]; a1++) {
547 for (a2=0; a2<size[(a+2)%3]; a2++) {
550 pa = psys->particles + a1*a1mul + a2*a2mul;
551 copy_v3_v3(co1, pa->fuv);
552 co1[a] -= d < delta[a] ? d/2.f : delta[a]/2.f;
553 copy_v3_v3(co2, co1);
554 co2[a] += delta[a] + 0.001f*d;
557 /* lets intersect the faces */
558 for (i=0; i<totface; i++,mface++) {
559 copy_v3_v3(v1, mvert[mface->v1].co);
560 copy_v3_v3(v2, mvert[mface->v2].co);
561 copy_v3_v3(v3, mvert[mface->v3].co);
563 if (isect_axial_line_tri_v3(a, co1, co2, v2, v3, v1, &lambda)) {
564 if (from==PART_FROM_FACE)
565 (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
566 else /* store number of intersections */
567 (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
569 else if (mface->v4) {
570 copy_v3_v3(v4, mvert[mface->v4].co);
572 if (isect_axial_line_tri_v3(a, co1, co2, v4, v1, v3, &lambda)) {
573 if (from==PART_FROM_FACE)
574 (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
576 (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
581 if (from==PART_FROM_VOLUME) {
582 int in=pa->hair_index%2;
583 if (in) pa->hair_index++;
584 for (i=0; i<size[0]; i++) {
585 if (in || (pa+i*a0mul)->hair_index%2)
586 (pa+i*a0mul)->flag &= ~PARS_UNEXIST;
587 /* odd intersections == in->out / out->in */
588 /* even intersections -> in stays same */
589 in=(in + (pa+i*a0mul)->hair_index) % 2;
597 if (psys->part->flag & PART_GRID_HEXAGONAL) {
598 for (i=0,p=0,pa=psys->particles; i<res; i++) {
599 for (j=0; j<res; j++) {
600 for (k=0; k<res; k++,p++,pa++) {
613 if (psys->part->flag & PART_GRID_INVERT) {
614 for (i=0; i<size[0]; i++) {
615 for (j=0; j<size[1]; j++) {
616 pa=psys->particles + res*(i*res + j);
617 for (k=0; k<size[2]; k++, pa++) {
618 pa->flag ^= PARS_UNEXIST;
624 if (psys->part->grid_rand > 0.f) {
625 float rfac = d * psys->part->grid_rand;
626 for (p=0,pa=psys->particles; p<psys->totpart; p++,pa++) {
627 if (pa->flag & PARS_UNEXIST)
630 pa->fuv[0] += rfac * (PSYS_FRAND(p + 31) - 0.5f);
631 pa->fuv[1] += rfac * (PSYS_FRAND(p + 32) - 0.5f);
632 pa->fuv[2] += rfac * (PSYS_FRAND(p + 33) - 0.5f);
637 /* modified copy from rayshade.c */
638 static void hammersley_create(float *out, int n, int seed, float amount)
641 double p, t, offs[2];
644 rng = BLI_rng_new(31415926 + n + seed);
645 offs[0] = BLI_rng_get_double(rng) + (double)amount;
646 offs[1] = BLI_rng_get_double(rng) + (double)amount;
649 for (k = 0; k < n; k++) {
651 for (p = 0.5, kk = k; kk; p *= 0.5, kk >>= 1)
652 if (kk & 1) /* kk mod 2 = 1 */
655 out[2*k + 0] = fmod((double)k/(double)n + offs[0], 1.0);
656 out[2*k + 1] = fmod(t + offs[1], 1.0);
660 /* modified copy from effect.c */
661 static void init_mv_jit(float *jit, int num, int seed2, float amount)
664 float *jit2, x, rad1, rad2, rad3;
669 rad1= (float)(1.0f/sqrtf((float)num));
670 rad2= (float)(1.0f/((float)num));
671 rad3= (float)sqrt((float)num)/((float)num);
673 rng = BLI_rng_new(31415926 + num + seed2);
676 for (i=0; i<num2; i+=2) {
678 jit[i] = x + amount*rad1*(0.5f - BLI_rng_get_float(rng));
679 jit[i+1] = i/(2.0f*num) + amount*rad1*(0.5f - BLI_rng_get_float(rng));
681 jit[i]-= (float)floor(jit[i]);
682 jit[i+1]-= (float)floor(jit[i+1]);
685 x -= (float)floor(x);
688 jit2= MEM_mallocN(12 + 2*sizeof(float)*num, "initjit");
690 for (i=0 ; i<4 ; i++) {
691 BLI_jitterate1(jit, jit2, num, rad1);
692 BLI_jitterate1(jit, jit2, num, rad1);
693 BLI_jitterate2(jit, jit2, num, rad2);
699 static void psys_uv_to_w(float u, float v, int quad, float *w)
701 float vert[4][3], co[3];
710 vert[0][0] = 0.0f; vert[0][1] = 0.0f; vert[0][2] = 0.0f;
711 vert[1][0] = 1.0f; vert[1][1] = 0.0f; vert[1][2] = 0.0f;
712 vert[2][0] = 1.0f; vert[2][1] = 1.0f; vert[2][2] = 0.0f;
719 vert[3][0] = 0.0f; vert[3][1] = 1.0f; vert[3][2] = 0.0f;
720 interp_weights_poly_v3( w,vert, 4, co);
723 interp_weights_poly_v3( w,vert, 3, co);
728 /* Find the index in "sum" array before "value" is crossed. */
729 static int distribute_binary_search(float *sum, int n, float value)
731 int mid, low=0, high=n;
736 while (low <= high) {
739 if (sum[mid] < value && value <= sum[mid+1])
742 if (sum[mid] >= value)
744 else if (sum[mid] < value)
753 /* the max number if calls to rng_* funcs within psys_thread_distribute_particle
754 * be sure to keep up to date if this changes */
755 #define PSYS_RND_DIST_SKIP 2
757 /* note: this function must be thread safe, for from == PART_FROM_CHILD */
758 #define ONLY_WORKING_WITH_PA_VERTS 0
759 static void distribute_threads_exec(ParticleThread *thread, ParticleData *pa, ChildParticle *cpa, int p)
761 ParticleThreadContext *ctx= thread->ctx;
762 Object *ob= ctx->sim.ob;
763 DerivedMesh *dm= ctx->dm;
764 float *v1, *v2, *v3, *v4, nor[3], orco1[3], co1[3], co2[3], nor1[3];
765 float cur_d, min_d, randu, randv;
767 int cfrom= ctx->cfrom;
768 int distr= ctx->distr;
769 int i, intersect, tot;
770 int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
772 if (from == PART_FROM_VERT) {
773 /* TODO_PARTICLE - use original index */
774 pa->num= ctx->index[p];
776 pa->fuv[1] = pa->fuv[2] = pa->fuv[3] = 0.0;
778 #if ONLY_WORKING_WITH_PA_VERTS
780 KDTreeNearest ptn[3];
783 psys_particle_on_dm(ctx->dm,from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co1,0,0,0,orco1,0);
784 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &orco1, 1, 1);
785 maxw = BLI_kdtree_find_n_nearest(ctx->tree,3,orco1,NULL,ptn);
787 for (w=0; w<maxw; w++) {
788 pa->verts[w]=ptn->num;
793 else if (from == PART_FROM_FACE || from == PART_FROM_VOLUME) {
796 pa->num = i = ctx->index[p];
797 mface = dm->getTessFaceData(dm,i,CD_MFACE);
801 if (ctx->jitlevel == 1) {
803 psys_uv_to_w(0.5f, 0.5f, mface->v4, pa->fuv);
805 psys_uv_to_w(0.33333f, 0.33333f, mface->v4, pa->fuv);
808 ctx->jitoff[i] = fmod(ctx->jitoff[i],(float)ctx->jitlevel);
809 if (!isnan(ctx->jitoff[i])) {
810 psys_uv_to_w(ctx->jit[2*(int)ctx->jitoff[i]], ctx->jit[2*(int)ctx->jitoff[i]+1], mface->v4, pa->fuv);
815 case PART_DISTR_RAND:
816 randu= BLI_rng_get_float(thread->rng);
817 randv= BLI_rng_get_float(thread->rng);
820 psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
826 if (from==PART_FROM_VOLUME) {
827 MVert *mvert=dm->getVertDataArray(dm,CD_MVERT);
829 tot=dm->getNumTessFaces(dm);
831 psys_interpolate_face(mvert,mface,0,0,pa->fuv,co1,nor,0,0,0,0);
834 mul_v3_fl(nor,-100.0);
836 add_v3_v3v3(co2,co1,nor);
841 for (i=0,mface=dm->getTessFaceDataArray(dm,CD_MFACE); i<tot; i++,mface++) {
842 if (i==pa->num) continue;
844 v1=mvert[mface->v1].co;
845 v2=mvert[mface->v2].co;
846 v3=mvert[mface->v3].co;
848 if (isect_line_tri_v3(co1, co2, v2, v3, v1, &cur_d, 0)) {
851 pa->foffset=cur_d*50.0f; /* to the middle of volume */
856 v4=mvert[mface->v4].co;
858 if (isect_line_tri_v3(co1, co2, v4, v1, v3, &cur_d, 0)) {
861 pa->foffset=cur_d*50.0f; /* to the middle of volume */
872 pa->foffset *= ctx->jit[p % (2 * ctx->jitlevel)];
874 case PART_DISTR_RAND:
875 pa->foffset *= BLI_frand();
881 else if (from == PART_FROM_CHILD) {
884 if (ctx->index[p] < 0) {
886 cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]=0.0f;
887 cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
891 mf= dm->getTessFaceData(dm, ctx->index[p], CD_MFACE);
893 randu= BLI_rng_get_float(thread->rng);
894 randv= BLI_rng_get_float(thread->rng);
897 psys_uv_to_w(randu, randv, mf->v4, cpa->fuv);
899 cpa->num = ctx->index[p];
902 KDTreeNearest ptn[10];
903 int w,maxw;//, do_seams;
904 float maxd /*, mind,dd */, totw= 0.0f;
908 psys_particle_on_dm(dm,cfrom,cpa->num,DMCACHE_ISCHILD,cpa->fuv,cpa->foffset,co1,nor1,NULL,NULL,orco1,NULL);
909 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &orco1, 1, 1);
910 maxw = BLI_kdtree_find_n_nearest(ctx->tree,4,orco1,NULL,ptn);
912 maxd=ptn[maxw-1].dist;
913 /* mind=ptn[0].dist; */ /* UNUSED */
915 /* the weights here could be done better */
916 for (w=0; w<maxw; w++) {
917 parent[w]=ptn[w].index;
918 pweight[w]=(float)pow(2.0,(double)(-6.0f*ptn[w].dist/maxd));
925 for (w=0,i=0; w<maxw && i<4; w++) {
927 cpa->pa[i]=parent[w];
928 cpa->w[i]=pweight[w];
938 if (totw>0.0f) for (w=0; w<4; w++)
941 cpa->parent=cpa->pa[0];
945 if (rng_skip_tot > 0) /* should never be below zero */
946 BLI_rng_skip(thread->rng, rng_skip_tot);
949 static void *distribute_threads_exec_cb(void *data)
951 ParticleThread *thread= (ParticleThread*)data;
952 ParticleSystem *psys= thread->ctx->sim.psys;
957 if (thread->ctx->from == PART_FROM_CHILD) {
958 totpart= psys->totchild;
961 for (p=0; p<totpart; p++, cpa++) {
962 if (thread->ctx->skip) /* simplification skip */
963 BLI_rng_skip(thread->rng, PSYS_RND_DIST_SKIP * thread->ctx->skip[p]);
965 if ((p+thread->num) % thread->tot == 0)
966 distribute_threads_exec(thread, NULL, cpa, p);
967 else /* thread skip */
968 BLI_rng_skip(thread->rng, PSYS_RND_DIST_SKIP);
972 totpart= psys->totpart;
973 pa= psys->particles + thread->num;
974 for (p=thread->num; p<totpart; p+=thread->tot, pa+=thread->tot)
975 distribute_threads_exec(thread, pa, NULL, p);
981 /* not thread safe, but qsort doesn't take userdata argument */
982 static int *COMPARE_ORIG_INDEX = NULL;
983 static int distribute_compare_orig_index(const void *p1, const void *p2)
985 int index1 = COMPARE_ORIG_INDEX[*(const int *)p1];
986 int index2 = COMPARE_ORIG_INDEX[*(const int *)p2];
990 else if (index1 == index2) {
991 /* this pointer comparison appears to make qsort stable for glibc,
992 * and apparently on solaris too, makes the renders reproducible */
1004 static void distribute_invalid(Scene *scene, ParticleSystem *psys, int from)
1006 if (from == PART_FROM_CHILD) {
1008 int p, totchild = get_psys_tot_child(scene, psys);
1010 if (psys->child && totchild) {
1011 for (p=0,cpa=psys->child; p<totchild; p++,cpa++) {
1012 cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3] = 0.0;
1015 cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
1023 pa->fuv[0] = pa->fuv[1] = pa->fuv[2] = pa->fuv[3] = 0.0;
1030 /* Creates a distribution of coordinates on a DerivedMesh */
1031 /* This is to denote functionality that does not yet work with mesh - only derived mesh */
1032 static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, DerivedMesh *finaldm, int from)
1034 ParticleThreadContext *ctx= threads[0].ctx;
1035 Object *ob= ctx->sim.ob;
1036 ParticleSystem *psys= ctx->sim.psys;
1037 ParticleData *pa=0, *tpars= 0;
1038 ParticleSettings *part;
1039 ParticleSeam *seams= 0;
1041 DerivedMesh *dm= NULL;
1043 int i, seed, p=0, totthread= threads[0].tot;
1045 int totelem=0, totpart, *particle_element=0, children=0, totseam=0;
1046 int jitlevel= 1, distr;
1047 float *element_weight=NULL,*element_sum=NULL,*jitter_offset=NULL, *vweight=NULL;
1048 float cur, maxweight=0.0, tweight, totweight, inv_totweight, co[3], nor[3], orco[3], ornor[3];
1050 if (ELEM3(NULL, ob, psys, psys->part))
1054 totpart=psys->totpart;
1058 if (!finaldm->deformedOnly && !finaldm->getTessFaceDataArray(finaldm, CD_ORIGINDEX)) {
1059 printf("Can't create particles with the current modifier stack, disable destructive modifiers\n");
1060 // XXX error("Can't paint with the current modifier stack, disable destructive modifiers");
1064 /* First handle special cases */
1065 if (from == PART_FROM_CHILD) {
1066 /* Simple children */
1067 if (part->childtype != PART_CHILD_FACES) {
1068 BLI_srandom(31415926 + psys->seed + psys->child_seed);
1069 distribute_simple_children(scene, ob, finaldm, psys);
1074 /* Grid distribution */
1075 if (part->distr==PART_DISTR_GRID && from != PART_FROM_VERT) {
1076 BLI_srandom(31415926 + psys->seed);
1077 dm= CDDM_from_mesh((Mesh*)ob->data, ob);
1078 DM_ensure_tessface(dm);
1079 distribute_grid(dm,psys);
1085 /* Create trees and original coordinates if needed */
1086 if (from == PART_FROM_CHILD) {
1087 distr=PART_DISTR_RAND;
1088 BLI_srandom(31415926 + psys->seed + psys->child_seed);
1092 DM_ensure_tessface(dm);
1096 tree=BLI_kdtree_new(totpart);
1098 for (p=0,pa=psys->particles; p<totpart; p++,pa++) {
1099 psys_particle_on_dm(dm,part->from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co,nor,0,0,orco,ornor);
1100 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &orco, 1, 1);
1101 BLI_kdtree_insert(tree, p, orco, ornor);
1104 BLI_kdtree_balance(tree);
1106 totpart = get_psys_tot_child(scene, psys);
1107 cfrom = from = PART_FROM_FACE;
1110 distr = part->distr;
1111 BLI_srandom(31415926 + psys->seed);
1113 dm= CDDM_from_mesh((Mesh*)ob->data, ob);
1115 /* BMESH ONLY, for verts we don't care about tessfaces */
1116 if (from != PART_FROM_VERT) {
1117 DM_ensure_tessface(dm);
1120 /* we need orco for consistent distributions */
1121 DM_add_vert_layer(dm, CD_ORCO, CD_ASSIGN, BKE_mesh_orco_verts_get(ob));
1123 if (from == PART_FROM_VERT) {
1124 MVert *mv= dm->getVertDataArray(dm, CD_MVERT);
1125 float (*orcodata)[3] = dm->getVertDataArray(dm, CD_ORCO);
1126 int totvert = dm->getNumVerts(dm);
1128 tree=BLI_kdtree_new(totvert);
1130 for (p=0; p<totvert; p++) {
1132 copy_v3_v3(co,orcodata[p]);
1133 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co, 1, 1);
1136 copy_v3_v3(co,mv[p].co);
1137 BLI_kdtree_insert(tree,p,co,NULL);
1140 BLI_kdtree_balance(tree);
1144 /* Get total number of emission elements and allocate needed arrays */
1145 totelem = (from == PART_FROM_VERT) ? dm->getNumVerts(dm) : dm->getNumTessFaces(dm);
1148 distribute_invalid(scene, psys, children ? PART_FROM_CHILD : 0);
1150 if (G.debug & G_DEBUG)
1151 fprintf(stderr,"Particle distribution error: Nothing to emit from!\n");
1153 if (dm != finaldm) dm->release(dm);
1155 BLI_kdtree_free(tree);
1160 element_weight = MEM_callocN(sizeof(float)*totelem, "particle_distribution_weights");
1161 particle_element= MEM_callocN(sizeof(int)*totpart, "particle_distribution_indexes");
1162 element_sum = MEM_callocN(sizeof(float)*(totelem+1), "particle_distribution_sum");
1163 jitter_offset = MEM_callocN(sizeof(float)*totelem, "particle_distribution_jitoff");
1165 /* Calculate weights from face areas */
1166 if ((part->flag&PART_EDISTR || children) && from != PART_FROM_VERT) {
1167 MVert *v1, *v2, *v3, *v4;
1168 float totarea=0.f, co1[3], co2[3], co3[3], co4[3];
1169 float (*orcodata)[3];
1171 orcodata= dm->getVertDataArray(dm, CD_ORCO);
1173 for (i=0; i<totelem; i++) {
1174 MFace *mf=dm->getTessFaceData(dm,i,CD_MFACE);
1177 copy_v3_v3(co1, orcodata[mf->v1]);
1178 copy_v3_v3(co2, orcodata[mf->v2]);
1179 copy_v3_v3(co3, orcodata[mf->v3]);
1180 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co1, 1, 1);
1181 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co2, 1, 1);
1182 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co3, 1, 1);
1184 copy_v3_v3(co4, orcodata[mf->v4]);
1185 BKE_mesh_orco_verts_transform((Mesh*)ob->data, &co4, 1, 1);
1189 v1= (MVert*)dm->getVertData(dm,mf->v1,CD_MVERT);
1190 v2= (MVert*)dm->getVertData(dm,mf->v2,CD_MVERT);
1191 v3= (MVert*)dm->getVertData(dm,mf->v3,CD_MVERT);
1192 copy_v3_v3(co1, v1->co);
1193 copy_v3_v3(co2, v2->co);
1194 copy_v3_v3(co3, v3->co);
1196 v4= (MVert*)dm->getVertData(dm,mf->v4,CD_MVERT);
1197 copy_v3_v3(co4, v4->co);
1201 cur = mf->v4 ? area_quad_v3(co1, co2, co3, co4) : area_tri_v3(co1, co2, co3);
1203 if (cur > maxweight)
1206 element_weight[i] = cur;
1210 for (i=0; i<totelem; i++)
1211 element_weight[i] /= totarea;
1213 maxweight /= totarea;
1216 float min=1.0f/(float)(MIN2(totelem,totpart));
1217 for (i=0; i<totelem; i++)
1218 element_weight[i]=min;
1222 /* Calculate weights from vgroup */
1223 vweight = psys_cache_vgroup(dm,psys,PSYS_VG_DENSITY);
1226 if (from==PART_FROM_VERT) {
1227 for (i=0;i<totelem; i++)
1228 element_weight[i]*=vweight[i];
1230 else { /* PART_FROM_FACE / PART_FROM_VOLUME */
1231 for (i=0;i<totelem; i++) {
1232 MFace *mf=dm->getTessFaceData(dm,i,CD_MFACE);
1233 tweight = vweight[mf->v1] + vweight[mf->v2] + vweight[mf->v3];
1236 tweight += vweight[mf->v4];
1243 element_weight[i]*=tweight;
1249 /* Calculate total weight of all elements */
1251 for (i=0;i<totelem; i++)
1252 totweight += element_weight[i];
1254 inv_totweight = (totweight > 0.f ? 1.f/totweight : 0.f);
1256 /* Calculate cumulative weights */
1257 element_sum[0] = 0.0f;
1258 for (i=0; i<totelem; i++)
1259 element_sum[i+1] = element_sum[i] + element_weight[i] * inv_totweight;
1261 /* Finally assign elements to particles */
1262 if ((part->flag&PART_TRAND) || (part->simplify_flag&PART_SIMPLIFY_ENABLE)) {
1265 for (p=0; p<totpart; p++) {
1266 /* In theory element_sum[totelem] should be 1.0, but due to float errors this is not necessarily always true, so scale pos accordingly. */
1267 pos= BLI_frand() * element_sum[totelem];
1268 particle_element[p] = distribute_binary_search(element_sum, totelem, pos);
1269 particle_element[p] = MIN2(totelem-1, particle_element[p]);
1270 jitter_offset[particle_element[p]] = pos;
1276 step= (totpart < 2) ? 0.5 : 1.0/(double)totpart;
1277 pos= 1e-6; /* tiny offset to avoid zero weight face */
1280 for (p=0; p<totpart; p++, pos+=step) {
1281 while ((i < totelem) && (pos > (double)element_sum[i + 1]))
1284 particle_element[p] = MIN2(totelem-1, i);
1286 /* avoid zero weight face */
1287 if (p == totpart-1 && element_weight[particle_element[p]] == 0.0f)
1288 particle_element[p] = particle_element[p-1];
1290 jitter_offset[particle_element[p]] = pos;
1294 MEM_freeN(element_sum);
1296 /* For hair, sort by origindex (allows optimization's in rendering), */
1297 /* however with virtual parents the children need to be in random order. */
1298 if (part->type == PART_HAIR && !(part->childtype==PART_CHILD_FACES && part->parents!=0.0f)) {
1299 COMPARE_ORIG_INDEX = NULL;
1301 if (from == PART_FROM_VERT) {
1302 if (dm->numVertData)
1303 COMPARE_ORIG_INDEX= dm->getVertDataArray(dm, CD_ORIGINDEX);
1306 if (dm->numTessFaceData)
1307 COMPARE_ORIG_INDEX= dm->getTessFaceDataArray(dm, CD_ORIGINDEX);
1310 if (COMPARE_ORIG_INDEX) {
1311 qsort(particle_element, totpart, sizeof(int), distribute_compare_orig_index);
1312 COMPARE_ORIG_INDEX = NULL;
1316 /* Create jittering if needed */
1317 if (distr==PART_DISTR_JIT && ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
1318 jitlevel= part->userjit;
1320 if (jitlevel == 0) {
1321 jitlevel= totpart/totelem;
1322 if (part->flag & PART_EDISTR) jitlevel*= 2; /* looks better in general, not very scietific */
1323 if (jitlevel<3) jitlevel= 3;
1326 jit= MEM_callocN((2+ jitlevel*2)*sizeof(float), "jit");
1328 /* for small amounts of particles we use regular jitter since it looks
1329 * a bit better, for larger amounts we switch to hammersley sequence
1330 * because it is much faster */
1332 init_mv_jit(jit, jitlevel, psys->seed, part->jitfac);
1334 hammersley_create(jit, jitlevel+1, psys->seed, part->jitfac);
1335 BLI_array_randomize(jit, 2*sizeof(float), jitlevel, psys->seed); /* for custom jit or even distribution */
1338 /* Setup things for threaded distribution */
1341 ctx->totseam= totseam;
1342 ctx->sim.psys= psys;
1343 ctx->index= particle_element;
1345 ctx->jitlevel= jitlevel;
1346 ctx->jitoff= jitter_offset;
1347 ctx->weight= element_weight;
1348 ctx->maxweight= maxweight;
1349 ctx->from= (children)? PART_FROM_CHILD: from;
1356 totpart= psys_render_simplify_distribution(ctx, totpart);
1357 alloc_child_particles(psys, totpart);
1360 if (!children || psys->totchild < 10000)
1363 seed= 31415926 + ctx->sim.psys->seed;
1364 for (i=0; i<totthread; i++) {
1365 threads[i].rng= BLI_rng_new(seed);
1366 threads[i].tot= totthread;
1372 static void distribute_particles_on_dm(ParticleSimulationData *sim, int from)
1374 DerivedMesh *finaldm = sim->psmd->dm;
1376 ParticleThread *pthreads;
1377 ParticleThreadContext *ctx;
1380 pthreads= psys_threads_create(sim);
1382 if (!distribute_threads_init_data(pthreads, sim->scene, finaldm, from)) {
1383 psys_threads_free(pthreads);
1387 totthread= pthreads[0].tot;
1388 if (totthread > 1) {
1389 BLI_init_threads(&threads, distribute_threads_exec_cb, totthread);
1391 for (i=0; i<totthread; i++)
1392 BLI_insert_thread(&threads, &pthreads[i]);
1394 BLI_end_threads(&threads);
1397 distribute_threads_exec_cb(&pthreads[0]);
1399 psys_calc_dmcache(sim->ob, finaldm, sim->psys);
1401 ctx= pthreads[0].ctx;
1402 if (ctx->dm != finaldm)
1403 ctx->dm->release(ctx->dm);
1405 psys_threads_free(pthreads);
1408 /* ready for future use, to emit particles without geometry */
1409 static void distribute_particles_on_shape(ParticleSimulationData *sim, int UNUSED(from))
1411 distribute_invalid(sim->scene, sim->psys, 0);
1413 fprintf(stderr,"Shape emission not yet possible!\n");
1416 static void distribute_particles(ParticleSimulationData *sim, int from)
1423 distribute_particles_on_dm(sim, from);
1428 distribute_particles_on_shape(sim, from);
1431 distribute_invalid(sim->scene, sim->psys, from);
1433 fprintf(stderr,"Particle distribution error!\n");
1437 /* threaded child particle distribution and path caching */
1438 ParticleThread *psys_threads_create(ParticleSimulationData *sim)
1440 ParticleThread *threads;
1441 ParticleThreadContext *ctx;
1444 if (sim->scene->r.mode & R_FIXED_THREADS)
1445 totthread= sim->scene->r.threads;
1447 totthread= BLI_system_thread_count();
1449 threads= MEM_callocN(sizeof(ParticleThread)*totthread, "ParticleThread");
1450 ctx= MEM_callocN(sizeof(ParticleThreadContext), "ParticleThreadContext");
1453 ctx->dm= ctx->sim.psmd->dm;
1454 ctx->ma= give_current_material(sim->ob, sim->psys->part->omat);
1456 memset(threads, 0, sizeof(ParticleThread)*totthread);
1458 for (i=0; i<totthread; i++) {
1459 threads[i].ctx= ctx;
1461 threads[i].tot= totthread;
1467 void psys_threads_free(ParticleThread *threads)
1469 ParticleThreadContext *ctx= threads[0].ctx;
1470 int i, totthread= threads[0].tot;
1474 MEM_freeN(ctx->vg_length);
1476 MEM_freeN(ctx->vg_clump);
1478 MEM_freeN(ctx->vg_kink);
1480 MEM_freeN(ctx->vg_rough1);
1482 MEM_freeN(ctx->vg_rough2);
1484 MEM_freeN(ctx->vg_roughe);
1486 if (ctx->sim.psys->lattice) {
1487 end_latt_deform(ctx->sim.psys->lattice);
1488 ctx->sim.psys->lattice= NULL;
1492 if (ctx->jit) MEM_freeN(ctx->jit);
1493 if (ctx->jitoff) MEM_freeN(ctx->jitoff);
1494 if (ctx->weight) MEM_freeN(ctx->weight);
1495 if (ctx->index) MEM_freeN(ctx->index);
1496 if (ctx->skip) MEM_freeN(ctx->skip);
1497 if (ctx->seams) MEM_freeN(ctx->seams);
1498 //if (ctx->vertpart) MEM_freeN(ctx->vertpart);
1499 BLI_kdtree_free(ctx->tree);
1502 for (i=0; i<totthread; i++) {
1504 BLI_rng_free(threads[i].rng);
1505 if (threads[i].rng_path)
1506 BLI_rng_free(threads[i].rng_path);
1513 /* set particle parameters that don't change during particle's life */
1514 void initialize_particle(ParticleSimulationData *sim, ParticleData *pa, int p)
1516 ParticleSystem *psys = sim->psys;
1517 ParticleSettings *part = psys->part;
1518 ParticleTexture ptex;
1520 pa->flag &= ~PARS_UNEXIST;
1522 if (part->type != PART_FLUID) {
1523 psys_get_texture(sim, pa, &ptex, PAMAP_INIT, 0.f);
1525 if (ptex.exist < PSYS_FRAND(p+125))
1526 pa->flag |= PARS_UNEXIST;
1528 pa->time = (part->type == PART_HAIR) ? 0.f : part->sta + (part->end - part->sta)*ptex.time;
1532 /* we can't reset to -1 anymore since we've figured out correct index in distribute_particles */
1533 /* usage other than straight after distribute has to handle this index by itself - jahka*/
1534 //pa->num_dmcache = DMCACHE_NOTFOUND; /* assume we don't have a derived mesh face */
1536 static void initialize_all_particles(ParticleSimulationData *sim)
1538 ParticleSystem *psys = sim->psys;
1541 psys->totunexist = 0;
1544 if ((pa->flag & PARS_UNEXIST)==0)
1545 initialize_particle(sim, pa, p);
1547 if (pa->flag & PARS_UNEXIST)
1551 /* Free unexisting particles. */
1552 if (psys->totpart && psys->totunexist == psys->totpart) {
1553 if (psys->particles->boid)
1554 MEM_freeN(psys->particles->boid);
1556 MEM_freeN(psys->particles);
1557 psys->particles = NULL;
1558 psys->totpart = psys->totunexist = 0;
1561 if (psys->totunexist) {
1562 int newtotpart = psys->totpart - psys->totunexist;
1563 ParticleData *npa, *newpars;
1565 npa = newpars = MEM_callocN(newtotpart * sizeof(ParticleData), "particles");
1567 for (p=0, pa=psys->particles; p<newtotpart; p++, pa++, npa++) {
1568 while (pa->flag & PARS_UNEXIST)
1571 memcpy(npa, pa, sizeof(ParticleData));
1574 if (psys->particles->boid)
1575 MEM_freeN(psys->particles->boid);
1576 MEM_freeN(psys->particles);
1577 psys->particles = newpars;
1578 psys->totpart -= psys->totunexist;
1580 if (psys->particles->boid) {
1581 BoidParticle *newboids = MEM_callocN(psys->totpart * sizeof(BoidParticle), "boid particles");
1584 pa->boid = newboids++;
1590 static void get_angular_velocity_vector(short avemode, ParticleKey *state, float vec[3])
1593 case PART_AVE_VELOCITY:
1594 copy_v3_v3(vec, state->vel);
1596 case PART_AVE_HORIZONTAL:
1599 zvec[0] = zvec[1] = 0;
1601 cross_v3_v3v3(vec, state->vel, zvec);
1604 case PART_AVE_VERTICAL:
1606 float zvec[3], temp[3];
1607 zvec[0] = zvec[1] = 0;
1609 cross_v3_v3v3(temp, state->vel, zvec);
1610 cross_v3_v3v3(vec, temp, state->vel);
1613 case PART_AVE_GLOBAL_X:
1615 vec[1] = vec[2] = 0;
1617 case PART_AVE_GLOBAL_Y:
1619 vec[0] = vec[2] = 0;
1621 case PART_AVE_GLOBAL_Z:
1623 vec[0] = vec[1] = 0;
1628 void psys_get_birth_coordinates(ParticleSimulationData *sim, ParticleData *pa, ParticleKey *state, float dtime, float cfra)
1630 Object *ob = sim->ob;
1631 ParticleSystem *psys = sim->psys;
1632 ParticleSettings *part;
1633 ParticleTexture ptex;
1634 float fac, phasefac, nor[3] = {0,0,0},loc[3],vel[3] = {0.0,0.0,0.0},rot[4],q2[4];
1635 float r_vel[3],r_ave[3],r_rot[4],vec[3],p_vel[3] = {0.0,0.0,0.0};
1636 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};
1638 int p = pa - psys->particles;
1641 /* get birth location from object */
1642 if (part->tanfac != 0.f)
1643 psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,utan,vtan,0,0);
1645 psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,0,0,0,0);
1647 /* get possible textural influence */
1648 psys_get_texture(sim, pa, &ptex, PAMAP_IVEL, cfra);
1650 /* particles live in global space so */
1651 /* let's convert: */
1653 mul_m4_v3(ob->obmat, loc);
1656 mul_mat3_m4_v3(ob->obmat, nor);
1660 if (part->tanfac!=0.0f) {
1661 //float phase=vg_rot?2.0f*(psys_particle_value_from_verts(sim->psmd->dm,part->from,pa,vg_rot)-0.5f):0.0f;
1663 mul_v3_fl(vtan,-cosf((float)M_PI*(part->tanphase+phase)));
1664 fac= -sinf((float)M_PI*(part->tanphase+phase));
1665 madd_v3_v3fl(vtan, utan, fac);
1667 mul_mat3_m4_v3(ob->obmat,vtan);
1669 copy_v3_v3(utan, nor);
1670 mul_v3_fl(utan,dot_v3v3(vtan,nor));
1671 sub_v3_v3(vtan, utan);
1677 /* -velocity (boids need this even if there's no random velocity) */
1678 if (part->randfac != 0.0f || (part->phystype==PART_PHYS_BOIDS && pa->boid)) {
1679 r_vel[0] = 2.0f * (PSYS_FRAND(p + 10) - 0.5f);
1680 r_vel[1] = 2.0f * (PSYS_FRAND(p + 11) - 0.5f);
1681 r_vel[2] = 2.0f * (PSYS_FRAND(p + 12) - 0.5f);
1683 mul_mat3_m4_v3(ob->obmat, r_vel);
1684 normalize_v3(r_vel);
1687 /* -angular velocity */
1688 if (part->avemode==PART_AVE_RAND) {
1689 r_ave[0] = 2.0f * (PSYS_FRAND(p + 13) - 0.5f);
1690 r_ave[1] = 2.0f * (PSYS_FRAND(p + 14) - 0.5f);
1691 r_ave[2] = 2.0f * (PSYS_FRAND(p + 15) - 0.5f);
1693 mul_mat3_m4_v3(ob->obmat,r_ave);
1694 normalize_v3(r_ave);
1698 if (part->randrotfac != 0.0f) {
1699 r_rot[0] = 2.0f * (PSYS_FRAND(p + 16) - 0.5f);
1700 r_rot[1] = 2.0f * (PSYS_FRAND(p + 17) - 0.5f);
1701 r_rot[2] = 2.0f * (PSYS_FRAND(p + 18) - 0.5f);
1702 r_rot[3] = 2.0f * (PSYS_FRAND(p + 19) - 0.5f);
1703 normalize_qt(r_rot);
1705 mat4_to_quat(rot,ob->obmat);
1706 mul_qt_qtqt(r_rot,r_rot,rot);
1709 if (part->phystype==PART_PHYS_BOIDS && pa->boid) {
1710 float dvec[3], q[4], mat[3][3];
1712 copy_v3_v3(state->co,loc);
1714 /* boids don't get any initial velocity */
1715 zero_v3(state->vel);
1717 /* boids store direction in ave */
1718 if (fabsf(nor[2])==1.0f) {
1719 sub_v3_v3v3(state->ave, loc, ob->obmat[3]);
1720 normalize_v3(state->ave);
1723 copy_v3_v3(state->ave, nor);
1726 /* calculate rotation matrix */
1727 project_v3_v3v3(dvec, r_vel, state->ave);
1728 sub_v3_v3v3(mat[0], state->ave, dvec);
1729 normalize_v3(mat[0]);
1730 negate_v3_v3(mat[2], r_vel);
1731 normalize_v3(mat[2]);
1732 cross_v3_v3v3(mat[1], mat[2], mat[0]);
1734 /* apply rotation */
1735 mat3_to_quat_is_ok( q,mat);
1736 copy_qt_qt(state->rot, q);
1739 /* conversion done so now we apply new: */
1740 /* -velocity from: */
1744 sub_v3_v3v3(vel, pa->state.vel, pa->prev_state.vel);
1747 /* *emitter velocity */
1748 if (dtime != 0.f && part->obfac != 0.f) {
1749 sub_v3_v3v3(vel, loc, state->co);
1750 mul_v3_fl(vel, part->obfac/dtime);
1753 /* *emitter normal */
1754 if (part->normfac != 0.f)
1755 madd_v3_v3fl(vel, nor, part->normfac);
1757 /* *emitter tangent */
1758 if (sim->psmd && part->tanfac != 0.f)
1759 madd_v3_v3fl(vel, vtan, part->tanfac);
1761 /* *emitter object orientation */
1762 if (part->ob_vel[0] != 0.f) {
1763 normalize_v3_v3(vec, ob->obmat[0]);
1764 madd_v3_v3fl(vel, vec, part->ob_vel[0]);
1766 if (part->ob_vel[1] != 0.f) {
1767 normalize_v3_v3(vec, ob->obmat[1]);
1768 madd_v3_v3fl(vel, vec, part->ob_vel[1]);
1770 if (part->ob_vel[2] != 0.f) {
1771 normalize_v3_v3(vec, ob->obmat[2]);
1772 madd_v3_v3fl(vel, vec, part->ob_vel[2]);
1779 if (part->randfac != 0.f)
1780 madd_v3_v3fl(vel, r_vel, part->randfac);
1783 if (part->partfac != 0.f)
1784 madd_v3_v3fl(vel, p_vel, part->partfac);
1786 mul_v3_v3fl(state->vel, vel, ptex.ivel);
1788 /* -location from emitter */
1789 copy_v3_v3(state->co,loc);
1792 unit_qt(state->rot);
1794 if (part->rotmode) {
1795 /* create vector into which rotation is aligned */
1796 switch (part->rotmode) {
1798 copy_v3_v3(rot_vec, nor);
1801 copy_v3_v3(rot_vec, vel);
1803 case PART_ROT_GLOB_X:
1804 case PART_ROT_GLOB_Y:
1805 case PART_ROT_GLOB_Z:
1806 rot_vec[part->rotmode - PART_ROT_GLOB_X] = 1.0f;
1811 copy_v3_v3(rot_vec, ob->obmat[part->rotmode - PART_ROT_OB_X]);
1815 /* create rotation quat */
1817 vec_to_quat( q2,rot_vec, OB_POSX, OB_POSZ);
1819 /* randomize rotation quat */
1820 if (part->randrotfac!=0.0f)
1821 interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
1825 /* rotation phase */
1826 phasefac = part->phasefac;
1827 if (part->randphasefac != 0.0f)
1828 phasefac += part->randphasefac * PSYS_FRAND(p + 20);
1829 axis_angle_to_quat( q_phase,x_vec, phasefac*(float)M_PI);
1831 /* combine base rotation & phase */
1832 mul_qt_qtqt(state->rot, rot, q_phase);
1835 /* -angular velocity */
1837 zero_v3(state->ave);
1839 if (part->avemode) {
1840 if (part->avemode == PART_AVE_RAND)
1841 copy_v3_v3(state->ave, r_ave);
1843 get_angular_velocity_vector(part->avemode, state, state->ave);
1845 normalize_v3(state->ave);
1846 mul_v3_fl(state->ave, part->avefac);
1850 /* sets particle to the emitter surface with initial velocity & rotation */
1851 void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, float cfra)
1853 Object *ob = sim->ob;
1854 ParticleSystem *psys = sim->psys;
1855 ParticleSettings *part;
1856 ParticleTexture ptex;
1857 int p = pa - psys->particles;
1860 /* get precise emitter matrix if particle is born */
1861 if (part->type!=PART_HAIR && dtime > 0.f && pa->time < cfra && pa->time >= sim->psys->cfra) {
1862 /* we have to force RECALC_ANIM here since where_is_objec_time only does drivers */
1864 BKE_animsys_evaluate_animdata(sim->scene, &ob->id, ob->adt, pa->time, ADT_RECALC_ANIM);
1868 BKE_object_where_is_calc_time(sim->scene, ob, pa->time);
1870 psys->flag |= PSYS_OB_ANIM_RESTORE;
1873 psys_get_birth_coordinates(sim, pa, &pa->state, dtime, cfra);
1875 if (part->phystype==PART_PHYS_BOIDS && pa->boid) {
1876 BoidParticle *bpa = pa->boid;
1878 /* and gravity in r_ve */
1879 bpa->gravity[0] = bpa->gravity[1] = 0.0f;
1880 bpa->gravity[2] = -1.0f;
1881 if ((sim->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) &&
1882 (sim->scene->physics_settings.gravity[2] != 0.0f))
1884 bpa->gravity[2] = sim->scene->physics_settings.gravity[2];
1887 bpa->data.health = part->boids->health;
1888 bpa->data.mode = eBoidMode_InAir;
1889 bpa->data.state_id = ((BoidState*)part->boids->states.first)->id;
1890 bpa->data.acc[0]=bpa->data.acc[1]=bpa->data.acc[2]=0.0f;
1893 if (part->type == PART_HAIR) {
1894 pa->lifetime = 100.0f;
1897 /* get possible textural influence */
1898 psys_get_texture(sim, pa, &ptex, PAMAP_LIFE, cfra);
1900 pa->lifetime = part->lifetime * ptex.life;
1902 if (part->randlife != 0.0f)
1903 pa->lifetime *= 1.0f - part->randlife * PSYS_FRAND(p + 21);
1906 pa->dietime = pa->time + pa->lifetime;
1908 if (sim->psys->pointcache && sim->psys->pointcache->flag & PTCACHE_BAKED &&
1909 sim->psys->pointcache->mem_cache.first) {
1910 float dietime = psys_get_dietime_from_cache(sim->psys->pointcache, p);
1911 pa->dietime = MIN2(pa->dietime, dietime);
1914 if (pa->time > cfra)
1915 pa->alive = PARS_UNBORN;
1916 else if (pa->dietime <= cfra)
1917 pa->alive = PARS_DEAD;
1919 pa->alive = PARS_ALIVE;
1921 pa->state.time = cfra;
1923 static void reset_all_particles(ParticleSimulationData *sim, float dtime, float cfra, int from)
1926 int p, totpart=sim->psys->totpart;
1928 for (p=from, pa=sim->psys->particles+from; p<totpart; p++, pa++)
1929 reset_particle(sim, pa, dtime, cfra);
1931 /************************************************/
1932 /* Particle targets */
1933 /************************************************/
1934 ParticleSystem *psys_get_target_system(Object *ob, ParticleTarget *pt)
1936 ParticleSystem *psys = NULL;
1938 if (pt->ob == NULL || pt->ob == ob)
1939 psys = BLI_findlink(&ob->particlesystem, pt->psys-1);
1941 psys = BLI_findlink(&pt->ob->particlesystem, pt->psys-1);
1944 pt->flag |= PTARGET_VALID;
1946 pt->flag &= ~PTARGET_VALID;
1950 /************************************************/
1951 /* Keyed particles */
1952 /************************************************/
1953 /* Counts valid keyed targets */
1954 void psys_count_keyed_targets(ParticleSimulationData *sim)
1956 ParticleSystem *psys = sim->psys, *kpsys;
1957 ParticleTarget *pt = psys->targets.first;
1961 for (; pt; pt=pt->next) {
1962 kpsys = psys_get_target_system(sim->ob, pt);
1964 if (kpsys && kpsys->totpart) {
1965 psys->totkeyed += keys_valid;
1966 if (psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f)
1967 psys->totkeyed += 1;
1974 psys->totkeyed *= psys->flag & PSYS_KEYED_TIMING ? 1 : psys->part->keyed_loops;
1977 static void set_keyed_keys(ParticleSimulationData *sim)
1979 ParticleSystem *psys = sim->psys;
1980 ParticleSimulationData ksim= {0};
1984 int totpart = psys->totpart, k, totkeys = psys->totkeyed;
1987 ksim.scene= sim->scene;
1989 /* no proper targets so let's clear and bail out */
1990 if (psys->totkeyed==0) {
1991 free_keyed_keys(psys);
1992 psys->flag &= ~PSYS_KEYED;
1996 if (totpart && psys->particles->totkey != totkeys) {
1997 free_keyed_keys(psys);
1999 key = MEM_callocN(totpart*totkeys*sizeof(ParticleKey), "Keyed keys");
2003 pa->totkey = totkeys;
2008 psys->flag &= ~PSYS_KEYED;
2011 pt = psys->targets.first;
2012 for (k=0; k<totkeys; k++) {
2013 ksim.ob = pt->ob ? pt->ob : sim->ob;
2014 ksim.psys = BLI_findlink(&ksim.ob->particlesystem, pt->psys - 1);
2015 keyed_flag = (ksim.psys->flag & PSYS_KEYED);
2016 ksim.psys->flag &= ~PSYS_KEYED;
2020 key->time = -1.0; /* use current time */
2022 psys_get_particle_state(&ksim, p%ksim.psys->totpart, key, 1);
2024 if (psys->flag & PSYS_KEYED_TIMING) {
2025 key->time = pa->time + pt->time;
2026 if (pt->duration != 0.0f && k+1 < totkeys) {
2027 copy_particle_key(key+1, key, 1);
2028 (key+1)->time = pa->time + pt->time + pt->duration;
2031 else if (totkeys > 1)
2032 key->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
2034 key->time = pa->time;
2037 if (psys->flag & PSYS_KEYED_TIMING && pt->duration!=0.0f)
2040 ksim.psys->flag |= keyed_flag;
2042 pt = (pt->next && pt->next->flag & PTARGET_VALID)? pt->next : psys->targets.first;
2045 psys->flag |= PSYS_KEYED;
2048 /************************************************/
2050 /************************************************/
2051 void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
2053 PointCache *cache = psys->pointcache;
2055 if (cache->flag & PTCACHE_DISK_CACHE && cache->mem_cache.first == NULL) {
2057 BKE_ptcache_id_from_particles(&pid, ob, psys);
2058 cache->flag &= ~PTCACHE_DISK_CACHE;
2059 BKE_ptcache_disk_to_mem(&pid);
2060 cache->flag |= PTCACHE_DISK_CACHE;
2063 static void psys_clear_temp_pointcache(ParticleSystem *psys)
2065 if (psys->pointcache->flag & PTCACHE_DISK_CACHE)
2066 BKE_ptcache_free_mem(&psys->pointcache->mem_cache);
2068 void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra, int *efra)
2070 ParticleSettings *part = psys->part;
2072 *sfra = MAX2(1, (int)part->sta);
2073 *efra = MIN2((int)(part->end + part->lifetime + 1.0f), scene->r.efra);
2076 /************************************************/
2078 /************************************************/
2079 static void psys_update_particle_bvhtree(ParticleSystem *psys, float cfra)
2085 if (!psys->bvhtree || psys->bvhtree_frame != cfra) {
2086 LOOP_SHOWN_PARTICLES {
2090 BLI_bvhtree_free(psys->bvhtree);
2091 psys->bvhtree = BLI_bvhtree_new(totpart, 0.0, 4, 6);
2093 LOOP_SHOWN_PARTICLES {
2094 if (pa->alive == PARS_ALIVE) {
2095 if (pa->state.time == cfra)
2096 BLI_bvhtree_insert(psys->bvhtree, p, pa->prev_state.co, 1);
2098 BLI_bvhtree_insert(psys->bvhtree, p, pa->state.co, 1);
2101 BLI_bvhtree_balance(psys->bvhtree);
2103 psys->bvhtree_frame = cfra;
2107 void psys_update_particle_tree(ParticleSystem *psys, float cfra)
2113 if (!psys->tree || psys->tree_frame != cfra) {
2114 LOOP_SHOWN_PARTICLES {
2118 BLI_kdtree_free(psys->tree);
2119 psys->tree = BLI_kdtree_new(psys->totpart);
2121 LOOP_SHOWN_PARTICLES {
2122 if (pa->alive == PARS_ALIVE) {
2123 if (pa->state.time == cfra)
2124 BLI_kdtree_insert(psys->tree, p, pa->prev_state.co, NULL);
2126 BLI_kdtree_insert(psys->tree, p, pa->state.co, NULL);
2129 BLI_kdtree_balance(psys->tree);
2131 psys->tree_frame = cfra;
2136 static void psys_update_effectors(ParticleSimulationData *sim)
2138 pdEndEffectors(&sim->psys->effectors);
2139 sim->psys->effectors = pdInitEffectors(sim->scene, sim->ob, sim->psys, sim->psys->part->effector_weights);
2140 precalc_guides(sim, sim->psys->effectors);
2143 static void integrate_particle(ParticleSettings *part, ParticleData *pa, float dtime, float *external_acceleration,
2144 void (*force_func)(void *forcedata, ParticleKey *state, float *force, float *impulse),
2147 #define ZERO_F43 {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}
2149 ParticleKey states[5];
2150 float force[3], acceleration[3], impulse[3], dx[4][3] = ZERO_F43, dv[4][3] = ZERO_F43, oldpos[3];
2151 float pa_mass= (part->flag & PART_SIZEMASS ? part->mass * pa->size : part->mass);
2153 int integrator = part->integrator;
2157 copy_v3_v3(oldpos, pa->state.co);
2159 /* Verlet integration behaves strangely with moving emitters, so do first step with euler. */
2160 if (pa->prev_state.time < 0.f && integrator == PART_INT_VERLET)
2161 integrator = PART_INT_EULER;
2163 switch (integrator) {
2164 case PART_INT_EULER:
2167 case PART_INT_MIDPOINT:
2173 case PART_INT_VERLET:
2178 for (i=0; i<steps; i++) {
2179 copy_particle_key(states + i, &pa->state, 1);
2184 for (i=0; i<steps; i++) {
2188 force_func(forcedata, states+i, force, impulse);
2190 /* force to acceleration*/
2191 mul_v3_v3fl(acceleration, force, 1.0f/pa_mass);
2193 if (external_acceleration)
2194 add_v3_v3(acceleration, external_acceleration);
2196 /* calculate next state */
2197 add_v3_v3(states[i].vel, impulse);
2199 switch (integrator) {
2200 case PART_INT_EULER:
2201 madd_v3_v3v3fl(pa->state.co, states->co, states->vel, dtime);
2202 madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
2204 case PART_INT_MIDPOINT:
2206 madd_v3_v3v3fl(states[1].co, states->co, states->vel, dtime*0.5f);
2207 madd_v3_v3v3fl(states[1].vel, states->vel, acceleration, dtime*0.5f);
2208 states[1].time = dtime*0.5f;
2209 /*fra=sim->psys->cfra+0.5f*dfra;*/
2212 madd_v3_v3v3fl(pa->state.co, states->co, states[1].vel, dtime);
2213 madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
2219 copy_v3_v3(dx[0], states->vel);
2220 mul_v3_fl(dx[0], dtime);
2221 copy_v3_v3(dv[0], acceleration);
2222 mul_v3_fl(dv[0], dtime);
2224 madd_v3_v3v3fl(states[1].co, states->co, dx[0], 0.5f);
2225 madd_v3_v3v3fl(states[1].vel, states->vel, dv[0], 0.5f);
2226 states[1].time = dtime*0.5f;
2227 /*fra=sim->psys->cfra+0.5f*dfra;*/
2230 madd_v3_v3v3fl(dx[1], states->vel, dv[0], 0.5f);
2231 mul_v3_fl(dx[1], dtime);
2232 copy_v3_v3(dv[1], acceleration);
2233 mul_v3_fl(dv[1], dtime);
2235 madd_v3_v3v3fl(states[2].co, states->co, dx[1], 0.5f);
2236 madd_v3_v3v3fl(states[2].vel, states->vel, dv[1], 0.5f);
2237 states[2].time = dtime*0.5f;
2240 madd_v3_v3v3fl(dx[2], states->vel, dv[1], 0.5f);
2241 mul_v3_fl(dx[2], dtime);
2242 copy_v3_v3(dv[2], acceleration);
2243 mul_v3_fl(dv[2], dtime);
2245 add_v3_v3v3(states[3].co, states->co, dx[2]);
2246 add_v3_v3v3(states[3].vel, states->vel, dv[2]);
2247 states[3].time = dtime;
2251 add_v3_v3v3(dx[3], states->vel, dv[2]);
2252 mul_v3_fl(dx[3], dtime);
2253 copy_v3_v3(dv[3], acceleration);
2254 mul_v3_fl(dv[3], dtime);
2256 madd_v3_v3v3fl(pa->state.co, states->co, dx[0], 1.0f/6.0f);
2257 madd_v3_v3fl(pa->state.co, dx[1], 1.0f/3.0f);
2258 madd_v3_v3fl(pa->state.co, dx[2], 1.0f/3.0f);
2259 madd_v3_v3fl(pa->state.co, dx[3], 1.0f/6.0f);
2261 madd_v3_v3v3fl(pa->state.vel, states->vel, dv[0], 1.0f/6.0f);
2262 madd_v3_v3fl(pa->state.vel, dv[1], 1.0f/3.0f);
2263 madd_v3_v3fl(pa->state.vel, dv[2], 1.0f/3.0f);
2264 madd_v3_v3fl(pa->state.vel, dv[3], 1.0f/6.0f);
2267 case PART_INT_VERLET: /* Verlet integration */
2268 madd_v3_v3v3fl(pa->state.vel, pa->prev_state.vel, acceleration, dtime);
2269 madd_v3_v3v3fl(pa->state.co, pa->prev_state.co, pa->state.vel, dtime);
2271 sub_v3_v3v3(pa->state.vel, pa->state.co, oldpos);
2272 mul_v3_fl(pa->state.vel, 1.0f/dtime);
2278 /*********************************************************************************************************
2281 * In theory, there could be unlimited implementation of SPH simulators
2283 * This code uses in some parts adapted algorithms from the pseudo code as outlined in the Research paper:
2285 * Titled: Particle-based Viscoelastic Fluid Simulation.
2286 * Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
2287 * Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
2289 * Presented at Siggraph, (2005)
2291 * ********************************************************************************************************/
2292 #define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
2293 static ParticleSpring *sph_spring_add(ParticleSystem *psys, ParticleSpring *spring)
2295 /* Are more refs required? */
2296 if (psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
2297 psys->alloc_fluidsprings = PSYS_FLUID_SPRINGS_INITIAL_SIZE;
2298 psys->fluid_springs = (ParticleSpring*)MEM_callocN(psys->alloc_fluidsprings * sizeof(ParticleSpring), "Particle Fluid Springs");
2300 else if (psys->tot_fluidsprings == psys->alloc_fluidsprings) {
2301 /* Double the number of refs allocated */
2302 psys->alloc_fluidsprings *= 2;
2303 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
2306 memcpy(psys->fluid_springs + psys->tot_fluidsprings, spring, sizeof(ParticleSpring));
2307 psys->tot_fluidsprings++;
2309 return psys->fluid_springs + psys->tot_fluidsprings - 1;
2311 static void sph_spring_delete(ParticleSystem *psys, int j)
2313 if (j != psys->tot_fluidsprings - 1)
2314 psys->fluid_springs[j] = psys->fluid_springs[psys->tot_fluidsprings - 1];
2316 psys->tot_fluidsprings--;
2318 if (psys->tot_fluidsprings < psys->alloc_fluidsprings/2 && psys->alloc_fluidsprings > PSYS_FLUID_SPRINGS_INITIAL_SIZE) {
2319 psys->alloc_fluidsprings /= 2;
2320 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
2323 static void sph_springs_modify(ParticleSystem *psys, float dtime)
2325 SPHFluidSettings *fluid = psys->part->fluid;
2326 ParticleData *pa1, *pa2;
2327 ParticleSpring *spring = psys->fluid_springs;
2329 float h, d, Rij[3], rij, Lij;
2332 float yield_ratio = fluid->yield_ratio;
2333 float plasticity = fluid->plasticity_constant;
2334 /* scale things according to dtime */
2335 float timefix = 25.f * dtime;
2337 if ((fluid->flag & SPH_VISCOELASTIC_SPRINGS)==0 || fluid->spring_k == 0.f)
2340 /* Loop through the springs */
2341 for (i=0; i<psys->tot_fluidsprings; i++, spring++) {
2342 pa1 = psys->particles + spring->particle_index[0];
2343 pa2 = psys->particles + spring->particle_index[1];
2345 sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
2346 rij = normalize_v3(Rij);
2348 /* adjust rest length */
2349 Lij = spring->rest_length;
2350 d = yield_ratio * timefix * Lij;
2352 if (rij > Lij + d) // Stretch
2353 spring->rest_length += plasticity * (rij - Lij - d) * timefix;
2354 else if (rij < Lij - d) // Compress
2355 spring->rest_length -= plasticity * (Lij - d - rij) * timefix;
2359 if (spring->rest_length > h)
2360 spring->delete_flag = 1;
2363 /* Loop through springs backwaqrds - for efficient delete function */
2364 for (i=psys->tot_fluidsprings-1; i >= 0; i--) {
2365 if (psys->fluid_springs[i].delete_flag)
2366 sph_spring_delete(psys, i);
2369 static EdgeHash *sph_springhash_build(ParticleSystem *psys)
2371 EdgeHash *springhash = NULL;
2372 ParticleSpring *spring;
2375 springhash = BLI_edgehash_new();
2377 for (i=0, spring=psys->fluid_springs; i<psys->tot_fluidsprings; i++, spring++)
2378 BLI_edgehash_insert(springhash, spring->particle_index[0], spring->particle_index[1], SET_INT_IN_POINTER(i+1));
2383 #define SPH_NEIGHBORS 512
2384 typedef struct SPHNeighbor {
2385 ParticleSystem *psys;
2389 typedef struct SPHRangeData {
2390 SPHNeighbor neighbors[SPH_NEIGHBORS];
2395 ParticleSystem *npsys;
2403 static void sph_evaluate_func(BVHTree *tree, ParticleSystem **psys, float co[3], SPHRangeData *pfr, float interaction_radius, BVHTree_RangeQuery callback) {
2406 pfr->tot_neighbors = 0;
2408 for (i=0; i < 10 && psys[i]; i++) {
2409 pfr->npsys = psys[i];
2410 pfr->massfac = psys[i]->part->mass;
2411 pfr->use_size = psys[i]->part->flag & PART_SIZEMASS;
2414 BLI_bvhtree_range_query(tree, co, interaction_radius, callback, pfr);
2417 BLI_bvhtree_range_query(psys[i]->bvhtree, co, interaction_radius, callback, pfr);
2421 static void sph_density_accum_cb(void *userdata, int index, float squared_dist)
2423 SPHRangeData *pfr = (SPHRangeData *)userdata;
2424 ParticleData *npa = pfr->npsys->particles + index;
2428 if (npa == pfr->pa || squared_dist < FLT_EPSILON)
2431 /* Ugh! One particle has too many neighbors! If some aren't taken into
2432 * account, the forces will be biased by the tree search order. This
2433 * effectively adds enery to the system, and results in a churning motion.
2434 * But, we have to stop somewhere, and it's not the end of the world.
2437 if (pfr->tot_neighbors >= SPH_NEIGHBORS)
2440 pfr->neighbors[pfr->tot_neighbors].index = index;
2441 pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
2442 pfr->tot_neighbors++;
2444 dist = sqrtf(squared_dist);
2445 q = (1.f - dist/pfr->h) * pfr->massfac;
2450 pfr->data[0] += q*q;
2451 pfr->data[1] += q*q*q;
2455 * Find the Courant number for an SPH particle (used for adaptive time step).
2457 static void sph_particle_courant(SPHData *sphdata, SPHRangeData *pfr)
2459 ParticleData *pa, *npa;
2461 float flow[3], offset[3], dist;
2466 if (pfr->tot_neighbors > 0) {
2468 for (i=0; i < pfr->tot_neighbors; i++) {
2469 npa = pfr->neighbors[i].psys->particles + pfr->neighbors[i].index;
2470 sub_v3_v3v3(offset, pa->prev_state.co, npa->prev_state.co);
2471 dist += len_v3(offset);
2472 add_v3_v3(flow, npa->prev_state.vel);
2474 dist += sphdata->psys[0]->part->fluid->radius; // TODO: remove this? - z0r
2475 sphdata->element_size = dist / pfr->tot_neighbors;
2476 mul_v3_v3fl(sphdata->flow, flow, 1.0f / pfr->tot_neighbors);
2479 sphdata->element_size = MAXFLOAT;
2480 copy_v3_v3(sphdata->flow, flow);
2483 static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse))
2485 SPHData *sphdata = (SPHData *)sphdata_v;
2486 ParticleSystem **psys = sphdata->psys;
2487 ParticleData *pa = sphdata->pa;
2488 SPHFluidSettings *fluid = psys[0]->part->fluid;
2489 ParticleSpring *spring = NULL;
2492 float mass = sphdata->mass;
2493 float *gravity = sphdata->gravity;
2494 EdgeHash *springhash = sphdata->eh;
2496 float q, u, rij, dv[3];
2497 float pressure, near_pressure;
2499 float visc = fluid->viscosity_omega;
2500 float stiff_visc = fluid->viscosity_beta * (fluid->flag & SPH_FAC_VISCOSITY ? fluid->viscosity_omega : 1.f);
2502 float inv_mass = 1.0f/mass;
2503 float spring_constant = fluid->spring_k;
2505 /* 4.0 seems to be a pretty good value */
2506 float interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
2507 float h = interaction_radius * sphdata->hfac;
2508 float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f); /* 4.77 is an experimentally determined density factor */
2509 float rest_length = fluid->rest_length * (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f);
2511 float stiffness = fluid->stiffness_k;
2512 float stiffness_near_fac = fluid->stiffness_knear * (fluid->flag & SPH_FAC_REPULSION ? fluid->stiffness_k : 1.f);
2519 float density, near_density;
2521 int i, spring_index, index = pa - psys[0]->particles;
2523 data[0] = data[1] = 0;
2528 sph_evaluate_func( NULL, psys, state->co, &pfr, interaction_radius, sph_density_accum_cb);
2531 near_density = data[1];
2533 pressure = stiffness * (density - rest_density);
2534 near_pressure = stiffness_near_fac * near_density;
2536 pfn = pfr.neighbors;
2537 for (i=0; i<pfr.tot_neighbors; i++, pfn++) {
2538 npa = pfn->psys->particles + pfn->index;
2540 madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
2542 sub_v3_v3v3(vec, co, state->co);
2543 rij = normalize_v3(vec);
2545 q = (1.f - rij/h) * pfn->psys->part->mass * inv_mass;
2547 if (pfn->psys->part->flag & PART_SIZEMASS)
2550 copy_v3_v3(vel, npa->prev_state.vel);
2552 /* Double Density Relaxation */
2553 madd_v3_v3fl(force, vec, -(pressure + near_pressure*q)*q);
2556 if (visc > 0.f || stiff_visc > 0.f) {
2557 sub_v3_v3v3(dv, vel, state->vel);
2558 u = dot_v3v3(vec, dv);
2560 if (u < 0.f && visc > 0.f)
2561 madd_v3_v3fl(force, vec, 0.5f * q * visc * u );
2563 if (u > 0.f && stiff_visc > 0.f)
2564 madd_v3_v3fl(force, vec, 0.5f * q * stiff_visc * u );
2567 if (spring_constant > 0.f) {
2568 /* Viscoelastic spring force */
2569 if (pfn->psys == psys[0] && fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash) {
2570 /* BLI_edgehash_lookup appears to be thread-safe. - z0r */
2571 spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, pfn->index));
2574 spring = psys[0]->fluid_springs + spring_index - 1;
2576 madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (spring->rest_length - rij));
2578 else if (fluid->spring_frames == 0 || (pa->prev_state.time-pa->time) <= fluid->spring_frames) {
2579 ParticleSpring temp_spring;
2580 temp_spring.particle_index[0] = index;
2581 temp_spring.particle_index[1] = pfn->index;
2582 temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : rest_length;
2583 temp_spring.delete_flag = 0;
2585 /* sph_spring_add is not thread-safe. - z0r */
2586 #pragma omp critical
2587 sph_spring_add(psys[0], &temp_spring);
2590 else {/* PART_SPRING_HOOKES - Hooke's spring force */
2591 madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij/h) * (rest_length - rij));
2596 /* Artificial buoyancy force in negative gravity direction */
2597 if (fluid->buoyancy > 0.f && gravity)
2598 madd_v3_v3fl(force, gravity, fluid->buoyancy * (density-rest_density));
2600 if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF)
2601 sph_particle_courant(sphdata, &pfr);
2605 /* powf is really slow for raising to integer powers. */
2606 MINLINE float pow2(float x) {
2609 MINLINE float pow3(float x) {
2612 MINLINE float pow4(float x) {
2613 return pow2(pow2(x));
2615 MINLINE float pow7(float x) {
2616 return pow2(pow3(x)) * x;
2619 static void sphclassical_density_accum_cb(void *userdata, int index, float squared_dist)
2621 SPHRangeData *pfr = (SPHRangeData *)userdata;
2622 ParticleData *npa = pfr->npsys->particles + index;
2624 float qfac = 21.0f / (256.f * M_PI);
2628 /* Exclude particles that are more than 2h away. Can't use squared_dist here
2629 * because it is not accurate enough. Use current state, i.e. the output of
2630 * basic_integrate() - z0r */
2631 sub_v3_v3v3(vec, npa->state.co, pfr->pa->state.co);
2633 rij_h = rij / pfr->h;
2637 /* Smoothing factor. Utilise the Wendland kernel. gnuplot:
2638 * q1(x) = (2.0 - x)**4 * ( 1.0 + 2.0 * x)
2639 * plot [0:2] q1(x) */
2640 q = qfac / pow3(pfr->h) * pow4(2.0f - rij_h) * ( 1.0f + 2.0f * rij_h);
2647 pfr->data[1] += q / npa->sphdensity;
2650 static void sphclassical_neighbour_accum_cb(void *userdata, int index, float squared_dist)
2652 SPHRangeData *pfr = (SPHRangeData *)userdata;
2653 ParticleData *npa = pfr->npsys->particles + index;
2657 if (pfr->tot_neighbors >= SPH_NEIGHBORS)
2660 /* Exclude particles that are more than 2h away. Can't use squared_dist here
2661 * because it is not accurate enough. Use current state, i.e. the output of
2662 * basic_integrate() - z0r */
2663 sub_v3_v3v3(vec, npa->state.co, pfr->pa->state.co);
2665 rij_h = rij / pfr->h;
2669 pfr->neighbors[pfr->tot_neighbors].index = index;
2670 pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
2671 pfr->tot_neighbors++;
2673 static void sphclassical_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse))
2675 SPHData *sphdata = (SPHData *)sphdata_v;
2676 ParticleSystem **psys = sphdata->psys;
2677 ParticleData *pa = sphdata->pa;
2678 SPHFluidSettings *fluid = psys[0]->part->fluid;
2681 float *gravity = sphdata->gravity;
2683 float dq, u, rij, dv[3];
2684 float pressure, npressure;
2686 float visc = fluid->viscosity_omega;
2688 float interaction_radius;
2690 /* 4.77 is an experimentally determined density factor */
2691 float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.0f);
2693 float stiffness = fluid->stiffness_k;
2702 float qfac2 = 42.0f / (256.0f * M_PI);
2705 /* 4.0 here is to be consistent with previous formulation/interface */
2706 interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
2707 h = interaction_radius * sphdata->hfac;
2713 sph_evaluate_func(NULL, psys, state->co, &pfr, interaction_radius, sphclassical_neighbour_accum_cb);
2714 pressure = stiffness * (pow7(pa->sphdensity / rest_density) - 1.0f);
2716 /* multiply by mass so that we return a force, not accel */
2717 qfac2 *= sphdata->mass / pow3(pfr.h);
2719 pfn = pfr.neighbors;
2720 for (i = 0; i < pfr.tot_neighbors; i++, pfn++) {
2721 npa = pfn->psys->particles + pfn->index;
2723 /* we do not contribute to ourselves */
2727 /* Find vector to neighbour. Exclude particles that are more than 2h
2728 * away. Can't use current state here because it may have changed on
2729 * another thread - so do own mini integration. Unlike basic_integrate,
2730 * SPH integration depends on neighbouring particles. - z0r */
2731 madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
2732 sub_v3_v3v3(vec, co, state->co);
2733 rij = normalize_v3(vec);
2734 rij_h = rij / pfr.h;
2738 npressure = stiffness * (pow7(npa->sphdensity / rest_density) - 1.0f);
2740 /* First derivative of smoothing factor. Utilise the Wendland kernel.
2742 * q2(x) = 2.0 * (2.0 - x)**4 - 4.0 * (2.0 - x)**3 * (1.0 + 2.0 * x)
2744 * Particles > 2h away are excluded above. */
2745 dq = qfac2 * (2.0f * pow4(2.0f - rij_h) - 4.0f * pow3(2.0f - rij_h) * (1.0f + 2.0f * rij_h) );
2747 if (pfn->psys->part->flag & PART_SIZEMASS)
2750 pressureTerm = pressure / pow2(pa->sphdensity) + npressure / pow2(npa->sphdensity);
2752 /* Note that 'minus' is removed, because vec = vecBA, not vecAB.
2753 * This applies to the viscosity calculation below, too. */
2754 madd_v3_v3fl(force, vec, pressureTerm * dq);
2758 sub_v3_v3v3(dv, npa->prev_state.vel, pa->prev_state.vel);
2759 u = dot_v3v3(vec, dv);
2760 /* Apply parameters */
2761 u *= -dq * hinv * visc / (0.5f * npa->sphdensity + 0.5f * pa->sphdensity);
2762 madd_v3_v3fl(force, vec, u);
2766 /* Artificial buoyancy force in negative gravity direction */
2767 if (fluid->buoyancy > 0.f && gravity)
2768 madd_v3_v3fl(force, gravity, fluid->buoyancy * (pa->sphdensity - rest_density));
2770 if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF)
2771 sph_particle_courant(sphdata, &pfr);
2775 static void sphclassical_calc_dens(ParticleData *pa, float UNUSED(dfra), SPHData *sphdata){
2776 ParticleSystem **psys = sphdata->psys;
2777 SPHFluidSettings *fluid = psys[0]->part->fluid;
2778 /* 4.0 seems to be a pretty good value */
2779 float interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
2786 pfr.h = interaction_radius * sphdata->hfac;
2789 sph_evaluate_func( NULL, psys, pa->state.co, &pfr, interaction_radius, sphclassical_density_accum_cb);
2790 pa->sphdensity = MIN2(MAX2(data[0], fluid->rest_density * 0.9f), fluid->rest_density * 1.1f);
2793 void psys_sph_init(ParticleSimulationData *sim, SPHData *sphdata)
2798 // Add other coupled particle systems.
2799 sphdata->psys[0] = sim->psys;
2800 for (i=1, pt=sim->psys->targets.first; i<10; i++, pt=(pt?pt->next:NULL))
2801 sphdata->psys[i] = pt ? psys_get_target_system(sim->ob, pt) : NULL;
2803 if (psys_uses_gravity(sim))
2804 sphdata->gravity = sim->scene->physics_settings.gravity;
2806 sphdata->gravity = NULL;
2807 sphdata->eh = sph_springhash_build(sim->psys);
2809 // These per-particle values should be overridden later, but just for
2810 // completeness we give them default values now.
2812 sphdata->mass = 1.0f;
2814 if (sim->psys->part->fluid->solver == SPH_SOLVER_DDR) {
2815 sphdata->force_cb = sph_force_cb;
2816 sphdata->density_cb = sph_density_accum_cb;
2817 sphdata->hfac = 1.0f;
2819 /* SPH_SOLVER_CLASSICAL */
2820 sphdata->force_cb = sphclassical_force_cb;
2821 sphdata->density_cb = sphclassical_density_accum_cb;
2822 sphdata->hfac = 0.5f;
2827 void psys_sph_finalise(SPHData *sphdata)
2830 BLI_edgehash_free(sphdata->eh, NULL);
2834 /* Sample the density field at a point in space. */
2835 void psys_sph_density(BVHTree *tree, SPHData *sphdata, float co[3], float vars[2]) {
2836 ParticleSystem **psys = sphdata->psys;
2837 SPHFluidSettings *fluid = psys[0]->part->fluid;
2838 /* 4.0 seems to be a pretty good value */
2839 float interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
2843 density[0] = density[1] = 0.0f;
2845 pfr.h = interaction_radius*sphdata->hfac;
2847 sph_evaluate_func(tree, psys, co, &pfr, interaction_radius, sphdata->density_cb);
2849 vars[0] = pfr.data[0];
2850 vars[1] = pfr.data[1];
2853 static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, SPHData *sphdata)
2855 ParticleSettings *part = sim->psys->part;
2856 // float timestep = psys_get_timestep(sim); // UNUSED
2857 float pa_mass = part->mass * (part->flag & PART_SIZEMASS ? pa->size : 1.f);
2858 float dtime = dfra*psys_get_timestep(sim);
2859 // int steps = 1; // UNUSED
2860 float effector_acceleration[3];
2863 sphdata->mass = pa_mass;
2865 //sphdata.element_size and sphdata.flow are set in the callback.
2867 /* restore previous state and treat gravity & effectors as external acceleration*/
2868 sub_v3_v3v3(effector_acceleration, pa->state.vel, pa->prev_state.vel);
2869 mul_v3_fl(effector_acceleration, 1.f/dtime);
2871 copy_particle_key(&pa->state, &pa->prev_state, 0);
2873 integrate_particle(part, pa, dtime, effector_acceleration, sphdata->force_cb, sphdata);
2876 /************************************************/
2878 /************************************************/
2879 typedef struct EfData {
2880 ParticleTexture ptex;
2881 ParticleSimulationData *sim;
2884 static void basic_force_cb(void *efdata_v, ParticleKey *state, float *force, float *impulse)
2886 EfData *efdata = (EfData *)efdata_v;
2887 ParticleSimulationData *sim = efdata->sim;
2888 ParticleSettings *part = sim->psys->part;
2889 ParticleData *pa = efdata->pa;
2890 EffectedPoint epoint;
2893 pd_point_from_particle(efdata->sim, efdata->pa, state, &epoint);
2894 if (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)
2895 pdDoEffectors(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);
2897 mul_v3_fl(force, efdata->ptex.field);
2898 mul_v3_fl(impulse, efdata->ptex.field);
2900 /* calculate air-particle interaction */
2901 if (part->dragfac != 0.0f)
2902 madd_v3_v3fl(force, state->vel, -part->dragfac * pa->size * pa->size * len_v3(state->vel));
2904 /* brownian force */
2905 if (part->brownfac != 0.0f) {
2906 force[0] += (BLI_frand()-0.5f) * part->brownfac;
2907 force[1] += (BLI_frand()-0.5f) * part->brownfac;
2908 force[2] += (BLI_frand()-0.5f) * part->brownfac;
2911 if (part->flag & PART_ROT_DYN && epoint.ave)
2912 copy_v3_v3(pa->state.ave, epoint.ave);
2914 /* gathers all forces that effect particles and calculates a new state for the particle */
2915 static void basic_integrate(ParticleSimulationData *sim, int p, float dfra, float cfra)
2917 ParticleSettings *part = sim->psys->part;
2918 ParticleData *pa = sim->psys->particles + p;
2920 float dtime=dfra*psys_get_timestep(sim), time;
2921 float *gravity = NULL, gr[3];
2924 psys_get_texture(sim, pa, &efdata.ptex, PAMAP_PHYSICS, cfra);
2929 /* add global acceleration (gravitation) */
2930 if (psys_uses_gravity(sim) &&
2931 /* normal gravity is too strong for hair so it's disabled by default */
2932 (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR))
2935 madd_v3_v3fl(gr, sim->scene->physics_settings.gravity, part->effector_weights->global_gravity * efdata.ptex.gravity);
2939 /* maintain angular velocity */
2940 copy_v3_v3(pa->state.ave, pa->prev_state.ave);
2942 integrate_particle(part, pa, dtime, gravity, basic_force_cb, &efdata);
2944 /* damp affects final velocity */
2945 if (part->dampfac != 0.f)
2946 mul_v3_fl(pa->state.vel, 1.f - part->dampfac * efdata.ptex.damp * 25.f * dtime);
2948 //copy_v3_v3(pa->state.ave, states->ave);
2950 /* finally we do guides */
2951 time=(cfra-pa->time)/pa->lifetime;
2952 CLAMP(time, 0.0f, 1.0f);
2954 copy_v3_v3(tkey.co,pa->state.co);
2955 copy_v3_v3(tkey.vel,pa->state.vel);
2956 tkey.time=pa->state.time;
2958 if (part->type != PART_HAIR) {
2959 if (do_guides(sim->psys->effectors, &tkey, p, time)) {
2960 copy_v3_v3(pa->state.co,tkey.co);
2961 /* guides don't produce valid velocity */
2962 sub_v3_v3v3(pa->state.vel, tkey.co, pa->prev_state.co);
2963 mul_v3_fl(pa->state.vel,1.0f/dtime);
2964 pa->state.time=tkey.time;
2968 static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
2970 float rotfac, rot1[4], rot2[4] = {1.0,0.0,0.0,0.0}, dtime=dfra*timestep, extrotfac;
2972 if ((part->flag & PART_ROTATIONS) == 0) {
2973 unit_qt(pa->state.rot);
2977 if (part->flag & PART_ROT_DYN) {
2978 extrotfac = len_v3(pa->state.ave);
2984 if ((part->flag & PART_ROT_DYN) && ELEM3(part->avemode, PART_AVE_VELOCITY, PART_AVE_HORIZONTAL, PART_AVE_VERTICAL)) {
2986 float len1 = len_v3(pa->prev_state.vel);
2987 float len2 = len_v3(pa->state.vel);
2990 if (len1 == 0.0f || len2 == 0.0f) {
2991 zero_v3(pa->state.ave);
2994 cross_v3_v3v3(pa->state.ave, pa->prev_state.vel, pa->state.vel);
2995 normalize_v3(pa->state.ave);
2996 angle = dot_v3v3(pa->prev_state.vel, pa->state.vel) / (len1 * len2);
2997 mul_v3_fl(pa->state.ave, saacos(angle) / dtime);
3000 get_angular_velocity_vector(part->avemode, &pa->state, vec);
3001 axis_angle_to_quat(rot2, vec, dtime*part->avefac);
3004 rotfac = len_v3(pa->state.ave);
3005 if (rotfac == 0.0f || (part->flag & PART_ROT_DYN)==0 || extrotfac == 0.0f) {
3009 axis_angle_to_quat(rot1,pa->state.ave,rotfac*dtime);
3011 mul_qt_qtqt(pa->state.rot,rot1,pa->prev_state.rot);
3012 mul_qt_qtqt(pa->state.rot,rot2,pa->state.rot);
3014 /* keep rotation quat in good health */
3015 normalize_qt(pa->state.rot);
3018 /************************************************/
3020 /************************************************/
3021 #define COLLISION_MAX_COLLISIONS 10
3022 #define COLLISION_MIN_RADIUS 0.001f