ce3e7dbf931d50c3785bedf983ef1f9b399ed37a
[blender.git] / source / blender / blenkernel / intern / effect.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) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): none yet.
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/blenkernel/intern/effect.c
29  *  \ingroup bke
30  */
31
32 #include <stddef.h>
33 #include <stdarg.h>
34
35 #include <math.h>
36 #include <stdlib.h>
37
38 #include "MEM_guardedalloc.h"
39
40 #include "DNA_curve_types.h"
41 #include "DNA_group_types.h"
42 #include "DNA_listBase.h"
43 #include "DNA_meshdata_types.h"
44 #include "DNA_object_types.h"
45 #include "DNA_object_force.h"
46 #include "DNA_particle_types.h"
47 #include "DNA_texture_types.h"
48 #include "DNA_scene_types.h"
49
50 #include "BLI_math.h"
51 #include "BLI_blenlib.h"
52 #include "BLI_noise.h"
53 #include "BLI_rand.h"
54 #include "BLI_utildefines.h"
55 #include "BLI_ghash.h"
56
57 #include "PIL_time.h"
58
59 #include "BKE_anim.h"           /* needed for where_on_path */
60 #include "BKE_collision.h"
61 #include "BKE_curve.h"
62 #include "BKE_displist.h"
63 #include "BKE_DerivedMesh.h"
64 #include "BKE_cdderivedmesh.h"
65 #include "BKE_effect.h"
66 #include "BKE_global.h"
67 #include "BKE_layer.h"
68 #include "BKE_library.h"
69 #include "BKE_modifier.h"
70 #include "BKE_object.h"
71 #include "BKE_particle.h"
72 #include "BKE_scene.h"
73 #include "BKE_smoke.h"
74
75 #include "DEG_depsgraph.h"
76
77 #include "RE_render_ext.h"
78 #include "RE_shader_ext.h"
79
80 /* fluid sim particle import */
81 #ifdef WITH_MOD_FLUID
82 #include "LBM_fluidsim.h"
83 #include <zlib.h>
84 #include <string.h>
85 #endif // WITH_MOD_FLUID
86
87 EffectorWeights *BKE_add_effector_weights(Group *group)
88 {
89         EffectorWeights *weights = MEM_callocN(sizeof(EffectorWeights), "EffectorWeights");
90         int i;
91
92         for (i=0; i<NUM_PFIELD_TYPES; i++)
93                 weights->weight[i] = 1.0f;
94
95         weights->global_gravity = 1.0f;
96
97         weights->group = group;
98
99         return weights;
100 }
101 PartDeflect *object_add_collision_fields(int type)
102 {
103         PartDeflect *pd;
104
105         pd= MEM_callocN(sizeof(PartDeflect), "PartDeflect");
106
107         pd->forcefield = type;
108         pd->pdef_sbdamp = 0.1f;
109         pd->pdef_sbift  = 0.2f;
110         pd->pdef_sboft  = 0.02f;
111         pd->seed = ((unsigned int)(ceil(PIL_check_seconds_timer()))+1) % 128;
112         pd->f_strength = 1.0f;
113         pd->f_damp = 1.0f;
114
115         /* set sensible defaults based on type */
116         switch (type) {
117                 case PFIELD_VORTEX:
118                         pd->shape = PFIELD_SHAPE_PLANE;
119                         break;
120                 case PFIELD_WIND:
121                         pd->shape = PFIELD_SHAPE_PLANE;
122                         pd->f_flow = 1.0f; /* realistic wind behavior */
123                         break;
124                 case PFIELD_TEXTURE:
125                         pd->f_size = 1.0f;
126                         break;
127                 case PFIELD_SMOKEFLOW:
128                         pd->f_flow = 1.0f;
129                         break;
130         }
131         pd->flag = PFIELD_DO_LOCATION|PFIELD_DO_ROTATION;
132
133         return pd;
134 }
135
136 /* ***************** PARTICLES ***************** */
137
138 /* -------------------------- Effectors ------------------ */
139 void free_partdeflect(PartDeflect *pd)
140 {
141         if (!pd)
142                 return;
143
144         if (pd->rng)
145                 BLI_rng_free(pd->rng);
146
147         MEM_freeN(pd);
148 }
149
150 static EffectorCache *new_effector_cache(const struct EvaluationContext *eval_ctx, Scene *scene, Object *ob, ParticleSystem *psys, PartDeflect *pd)
151 {
152         EffectorCache *eff = MEM_callocN(sizeof(EffectorCache), "EffectorCache");
153         eff->eval_ctx = eval_ctx;
154         eff->scene = scene;
155         eff->ob = ob;
156         eff->psys = psys;
157         eff->pd = pd;
158         eff->frame = -1;
159         return eff;
160 }
161 static void add_object_to_effectors(ListBase **effectors, const struct EvaluationContext *eval_ctx, Scene *scene, EffectorWeights *weights, Object *ob, Object *ob_src, bool for_simulation)
162 {
163         EffectorCache *eff = NULL;
164
165         if ( ob == ob_src )
166                 return;
167
168         if (for_simulation) {
169                 if (weights->weight[ob->pd->forcefield] == 0.0f )
170                         return;
171
172                 if (ob->pd->shape == PFIELD_SHAPE_POINTS && !ob->derivedFinal )
173                         return;
174         }
175
176         if (*effectors == NULL)
177                 *effectors = MEM_callocN(sizeof(ListBase), "effectors list");
178
179         eff = new_effector_cache(eval_ctx, scene, ob, NULL, ob->pd);
180
181         /* make sure imat is up to date */
182         invert_m4_m4(ob->imat, ob->obmat);
183
184         BLI_addtail(*effectors, eff);
185 }
186 static void add_particles_to_effectors(ListBase **effectors, const struct EvaluationContext *eval_ctx, Scene *scene, EffectorWeights *weights, Object *ob, ParticleSystem *psys, ParticleSystem *psys_src, bool for_simulation)
187 {
188         ParticleSettings *part= psys->part;
189
190         if ( !psys_check_enabled(ob, psys, G.is_rendering) )
191                 return;
192
193         if ( psys == psys_src && (part->flag & PART_SELF_EFFECT) == 0)
194                 return;
195
196         if ( part->pd && part->pd->forcefield && (!for_simulation || weights->weight[part->pd->forcefield] != 0.0f)) {
197                 if (*effectors == NULL)
198                         *effectors = MEM_callocN(sizeof(ListBase), "effectors list");
199
200                 BLI_addtail(*effectors, new_effector_cache(eval_ctx, scene, ob, psys, part->pd));
201         }
202
203         if (part->pd2 && part->pd2->forcefield && (!for_simulation || weights->weight[part->pd2->forcefield] != 0.0f)) {
204                 if (*effectors == NULL)
205                         *effectors = MEM_callocN(sizeof(ListBase), "effectors list");
206
207                 BLI_addtail(*effectors, new_effector_cache(eval_ctx, scene, ob, psys, part->pd2));
208         }
209 }
210
211 /* returns ListBase handle with objects taking part in the effecting */
212 ListBase *pdInitEffectors(
213         const struct EvaluationContext *eval_ctx, Scene *scene, Object *ob_src, ParticleSystem *psys_src,
214         EffectorWeights *weights, bool for_simulation)
215 {
216         SceneLayer *sl;
217         Base *base;
218         unsigned int layer= ob_src->lay;
219         ListBase *effectors = NULL;
220
221         /* eval_ctx is NULL during deg build */
222         if (eval_ctx) {
223                 sl = eval_ctx->scene_layer;
224         }
225         else {
226                 sl = BKE_scene_layer_context_active_PLACEHOLDER(scene);
227         }
228         
229         if (weights->group) {
230                 GroupObject *go;
231                 
232                 for (go= weights->group->gobject.first; go; go= go->next) {
233                         if ( (go->ob->lay & layer) ) {
234                                 if ( go->ob->pd && go->ob->pd->forcefield )
235                                         add_object_to_effectors(&effectors, eval_ctx, scene, weights, go->ob, ob_src, for_simulation);
236
237                                 if ( go->ob->particlesystem.first ) {
238                                         ParticleSystem *psys= go->ob->particlesystem.first;
239
240                                         for ( ; psys; psys=psys->next )
241                                                 add_particles_to_effectors(&effectors, eval_ctx, scene, weights, go->ob, psys, psys_src, for_simulation);
242                                 }
243                         }
244                 }
245         }
246         else {
247                 for (base = FIRSTBASE(sl); base; base = base->next) {
248                         if ( base->object->pd && base->object->pd->forcefield )
249                                 add_object_to_effectors(&effectors, eval_ctx, scene, weights, base->object, ob_src, for_simulation);
250
251                         if ( base->object->particlesystem.first ) {
252                                 ParticleSystem *psys= base->object->particlesystem.first;
253
254                                 for ( ; psys; psys=psys->next )
255                                         add_particles_to_effectors(&effectors, eval_ctx, scene, weights, base->object, psys, psys_src, for_simulation);
256                         }
257                 }
258         }
259         
260         if (for_simulation)
261                 pdPrecalculateEffectors(eval_ctx, effectors);
262         
263         return effectors;
264 }
265
266 void pdEndEffectors(ListBase **effectors)
267 {
268         if (*effectors) {
269                 EffectorCache *eff = (*effectors)->first;
270
271                 for (; eff; eff=eff->next) {
272                         if (eff->guide_data)
273                                 MEM_freeN(eff->guide_data);
274                 }
275
276                 BLI_freelistN(*effectors);
277                 MEM_freeN(*effectors);
278                 *effectors = NULL;
279         }
280 }
281
282 static void precalculate_effector(const struct EvaluationContext *eval_ctx, EffectorCache *eff)
283 {
284         unsigned int cfra = (unsigned int)(eff->scene->r.cfra >= 0 ? eff->scene->r.cfra : -eff->scene->r.cfra);
285         if (!eff->pd->rng)
286                 eff->pd->rng = BLI_rng_new(eff->pd->seed + cfra);
287         else
288                 BLI_rng_srandom(eff->pd->rng, eff->pd->seed + cfra);
289
290         if (eff->pd->forcefield == PFIELD_GUIDE && eff->ob->type==OB_CURVE) {
291                 Curve *cu= eff->ob->data;
292                 if (cu->flag & CU_PATH) {
293                         if (eff->ob->curve_cache == NULL || eff->ob->curve_cache->path==NULL || eff->ob->curve_cache->path->data==NULL)
294                                 BKE_displist_make_curveTypes(eval_ctx, eff->scene, eff->ob, 0);
295
296                         if (eff->ob->curve_cache->path && eff->ob->curve_cache->path->data) {
297                                 where_on_path(eff->ob, 0.0, eff->guide_loc, eff->guide_dir, NULL, &eff->guide_radius, NULL);
298                                 mul_m4_v3(eff->ob->obmat, eff->guide_loc);
299                                 mul_mat3_m4_v3(eff->ob->obmat, eff->guide_dir);
300                         }
301                 }
302         }
303         else if (eff->pd->shape == PFIELD_SHAPE_SURFACE) {
304                 eff->surmd = (SurfaceModifierData *)modifiers_findByType( eff->ob, eModifierType_Surface );
305                 if (eff->ob->type == OB_CURVE)
306                         eff->flag |= PE_USE_NORMAL_DATA;
307         }
308         else if (eff->psys)
309                 psys_update_particle_tree(eff->psys, eff->scene->r.cfra);
310
311         /* Store object velocity */
312         if (eff->ob) {
313                 float old_vel[3];
314
315                 BKE_object_where_is_calc_time(eval_ctx, eff->scene, eff->ob, cfra - 1.0f);
316                 copy_v3_v3(old_vel, eff->ob->obmat[3]);
317                 BKE_object_where_is_calc_time(eval_ctx, eff->scene, eff->ob, cfra);
318                 sub_v3_v3v3(eff->velocity, eff->ob->obmat[3], old_vel);
319         }
320 }
321
322 void pdPrecalculateEffectors(const struct EvaluationContext *eval_ctx, ListBase *effectors)
323 {
324         if (effectors) {
325                 EffectorCache *eff = effectors->first;
326                 for (; eff; eff=eff->next)
327                         precalculate_effector(eval_ctx, eff);
328         }
329 }
330
331
332 void pd_point_from_particle(ParticleSimulationData *sim, ParticleData *pa, ParticleKey *state, EffectedPoint *point)
333 {
334         ParticleSettings *part = sim->psys->part;
335         point->loc = state->co;
336         point->vel = state->vel;
337         point->index = pa - sim->psys->particles;
338         point->size = pa->size;
339         point->charge = 0.0f;
340         
341         if (part->pd && part->pd->forcefield == PFIELD_CHARGE)
342                 point->charge += part->pd->f_strength;
343
344         if (part->pd2 && part->pd2->forcefield == PFIELD_CHARGE)
345                 point->charge += part->pd2->f_strength;
346
347         point->vel_to_sec = 1.0f;
348         point->vel_to_frame = psys_get_timestep(sim);
349
350         point->flag = 0;
351
352         if (sim->psys->part->flag & PART_ROT_DYN) {
353                 point->ave = state->ave;
354                 point->rot = state->rot;
355         }
356         else
357                 point->ave = point->rot = NULL;
358
359         point->psys = sim->psys;
360 }
361
362 void pd_point_from_loc(Scene *scene, float *loc, float *vel, int index, EffectedPoint *point)
363 {
364         point->loc = loc;
365         point->vel = vel;
366         point->index = index;
367         point->size = 0.0f;
368
369         point->vel_to_sec = (float)scene->r.frs_sec;
370         point->vel_to_frame = 1.0f;
371
372         point->flag = 0;
373
374         point->ave = point->rot = NULL;
375         point->psys = NULL;
376 }
377 void pd_point_from_soft(Scene *scene, float *loc, float *vel, int index, EffectedPoint *point)
378 {
379         point->loc = loc;
380         point->vel = vel;
381         point->index = index;
382         point->size = 0.0f;
383
384         point->vel_to_sec = (float)scene->r.frs_sec;
385         point->vel_to_frame = 1.0f;
386
387         point->flag = PE_WIND_AS_SPEED;
388
389         point->ave = point->rot = NULL;
390
391         point->psys = NULL;
392 }
393 /************************************************/
394 /*                      Effectors               */
395 /************************************************/
396
397 // triangle - ray callback function
398 static void eff_tri_ray_hit(void *UNUSED(userData), int UNUSED(index), const BVHTreeRay *UNUSED(ray), BVHTreeRayHit *hit)
399 {       
400         /* whenever we hit a bounding box, we don't check further */
401         hit->dist = -1;
402         hit->index = 1;
403 }
404
405 // get visibility of a wind ray
406 static float eff_calc_visibility(ListBase *colliders, EffectorCache *eff, EffectorData *efd, EffectedPoint *point)
407 {
408         const int raycast_flag = BVH_RAYCAST_DEFAULT & ~(BVH_RAYCAST_WATERTIGHT);
409         ListBase *colls = colliders;
410         ColliderCache *col;
411         float norm[3], len = 0.0;
412         float visibility = 1.0, absorption = 0.0;
413         
414         if (!(eff->pd->flag & PFIELD_VISIBILITY))
415                 return visibility;
416
417         if (!colls)
418                 colls = get_collider_cache(eff->scene, eff->ob, NULL);
419
420         if (!colls)
421                 return visibility;
422
423         negate_v3_v3(norm, efd->vec_to_point);
424         len = normalize_v3(norm);
425         
426         /* check all collision objects */
427         for (col = colls->first; col; col = col->next) {
428                 CollisionModifierData *collmd = col->collmd;
429
430                 if (col->ob == eff->ob)
431                         continue;
432
433                 if (collmd->bvhtree) {
434                         BVHTreeRayHit hit;
435
436                         hit.index = -1;
437                         hit.dist = len + FLT_EPSILON;
438
439                         /* check if the way is blocked */
440                         if (BLI_bvhtree_ray_cast_ex(
441                                 collmd->bvhtree, point->loc, norm, 0.0f, &hit,
442                                 eff_tri_ray_hit, NULL, raycast_flag) != -1)
443                         {
444                                 absorption= col->ob->pd->absorption;
445
446                                 /* visibility is only between 0 and 1, calculated from 1-absorption */
447                                 visibility *= CLAMPIS(1.0f-absorption, 0.0f, 1.0f);
448                                 
449                                 if (visibility <= 0.0f)
450                                         break;
451                         }
452                 }
453         }
454
455         if (!colliders)
456                 free_collider_cache(&colls);
457         
458         return visibility;
459 }
460
461 // noise function for wind e.g.
462 static float wind_func(struct RNG *rng, float strength)
463 {
464         int random = (BLI_rng_get_int(rng)+1) % 128; // max 2357
465         float force = BLI_rng_get_float(rng) + 1.0f;
466         float ret;
467         float sign = 0;
468         
469         sign = ((float)random > 64.0f) ? 1.0f: -1.0f; // dividing by 2 is not giving equal sign distribution
470         
471         ret = sign*((float)random / force)*strength/128.0f;
472         
473         return ret;
474 }
475
476 /* maxdist: zero effect from this distance outwards (if usemax) */
477 /* mindist: full effect up to this distance (if usemin) */
478 /* power: falloff with formula 1/r^power */
479 static float falloff_func(float fac, int usemin, float mindist, int usemax, float maxdist, float power)
480 {
481         /* first quick checks */
482         if (usemax && fac > maxdist)
483                 return 0.0f;
484
485         if (usemin && fac < mindist)
486                 return 1.0f;
487
488         if (!usemin)
489                 mindist = 0.0;
490
491         return pow((double)(1.0f+fac-mindist), (double)(-power));
492 }
493
494 static float falloff_func_dist(PartDeflect *pd, float fac)
495 {
496         return falloff_func(fac, pd->flag&PFIELD_USEMIN, pd->mindist, pd->flag&PFIELD_USEMAX, pd->maxdist, pd->f_power);
497 }
498
499 static float falloff_func_rad(PartDeflect *pd, float fac)
500 {
501         return falloff_func(fac, pd->flag&PFIELD_USEMINR, pd->minrad, pd->flag&PFIELD_USEMAXR, pd->maxrad, pd->f_power_r);
502 }
503
504 float effector_falloff(EffectorCache *eff, EffectorData *efd, EffectedPoint *UNUSED(point), EffectorWeights *weights)
505 {
506         float temp[3];
507         float falloff = weights ? weights->weight[0] * weights->weight[eff->pd->forcefield] : 1.0f;
508         float fac, r_fac;
509
510         fac = dot_v3v3(efd->nor, efd->vec_to_point2);
511
512         if (eff->pd->zdir == PFIELD_Z_POS && fac < 0.0f)
513                 falloff=0.0f;
514         else if (eff->pd->zdir == PFIELD_Z_NEG && fac > 0.0f)
515                 falloff=0.0f;
516         else {
517                 switch (eff->pd->falloff) {
518                 case PFIELD_FALL_SPHERE:
519                         falloff*= falloff_func_dist(eff->pd, efd->distance);
520                         break;
521
522                 case PFIELD_FALL_TUBE:
523                         falloff*= falloff_func_dist(eff->pd, ABS(fac));
524                         if (falloff == 0.0f)
525                                 break;
526
527                         madd_v3_v3v3fl(temp, efd->vec_to_point2, efd->nor, -fac);
528                         r_fac= len_v3(temp);
529                         falloff*= falloff_func_rad(eff->pd, r_fac);
530                         break;
531                 case PFIELD_FALL_CONE:
532                         falloff*= falloff_func_dist(eff->pd, ABS(fac));
533                         if (falloff == 0.0f)
534                                 break;
535
536                         r_fac= RAD2DEGF(saacos(fac/len_v3(efd->vec_to_point)));
537                         falloff*= falloff_func_rad(eff->pd, r_fac);
538
539                         break;
540                 }
541         }
542
543         return falloff;
544 }
545
546 int closest_point_on_surface(SurfaceModifierData *surmd, const float co[3], float surface_co[3], float surface_nor[3], float surface_vel[3])
547 {
548         BVHTreeNearest nearest;
549
550         nearest.index = -1;
551         nearest.dist_sq = FLT_MAX;
552
553         BLI_bvhtree_find_nearest(surmd->bvhtree->tree, co, &nearest, surmd->bvhtree->nearest_callback, surmd->bvhtree);
554
555         if (nearest.index != -1) {
556                 copy_v3_v3(surface_co, nearest.co);
557
558                 if (surface_nor) {
559                         copy_v3_v3(surface_nor, nearest.no);
560                 }
561
562                 if (surface_vel) {
563                         const MLoop *mloop = surmd->bvhtree->loop;
564                         const MLoopTri *lt = &surmd->bvhtree->looptri[nearest.index];
565                         
566                         copy_v3_v3(surface_vel, surmd->v[mloop[lt->tri[0]].v].co);
567                         add_v3_v3(surface_vel, surmd->v[mloop[lt->tri[1]].v].co);
568                         add_v3_v3(surface_vel, surmd->v[mloop[lt->tri[2]].v].co);
569
570                         mul_v3_fl(surface_vel, (1.0f / 3.0f));
571                 }
572                 return 1;
573         }
574
575         return 0;
576 }
577 int get_effector_data(EffectorCache *eff, EffectorData *efd, EffectedPoint *point, int real_velocity)
578 {
579         float cfra = eff->scene->r.cfra;
580         int ret = 0;
581
582         /* In case surface object is in Edit mode when loading the .blend, surface modifier is never executed
583          * and bvhtree never built, see T48415. */
584         if (eff->pd && eff->pd->shape==PFIELD_SHAPE_SURFACE && eff->surmd && eff->surmd->bvhtree) {
585                 /* closest point in the object surface is an effector */
586                 float vec[3];
587
588                 /* using velocity corrected location allows for easier sliding over effector surface */
589                 copy_v3_v3(vec, point->vel);
590                 mul_v3_fl(vec, point->vel_to_frame);
591                 add_v3_v3(vec, point->loc);
592
593                 ret = closest_point_on_surface(eff->surmd, vec, efd->loc, efd->nor, real_velocity ? efd->vel : NULL);
594
595                 efd->size = 0.0f;
596         }
597         else if (eff->pd && eff->pd->shape==PFIELD_SHAPE_POINTS) {
598
599                 if (eff->ob->derivedFinal) {
600                         DerivedMesh *dm = eff->ob->derivedFinal;
601
602                         dm->getVertCo(dm, *efd->index, efd->loc);
603                         dm->getVertNo(dm, *efd->index, efd->nor);
604
605                         mul_m4_v3(eff->ob->obmat, efd->loc);
606                         mul_mat3_m4_v3(eff->ob->obmat, efd->nor);
607
608                         normalize_v3(efd->nor);
609
610                         efd->size = 0.0f;
611
612                         /**/
613                         ret = 1;
614                 }
615         }
616         else if (eff->psys) {
617                 ParticleData *pa = eff->psys->particles + *efd->index;
618                 ParticleKey state;
619
620                 /* exclude the particle itself for self effecting particles */
621                 if (eff->psys == point->psys && *efd->index == point->index) {
622                         /* pass */
623                 }
624                 else {
625                         ParticleSimulationData sim= {NULL};
626                         sim.eval_ctx = eff->eval_ctx;
627                         sim.scene= eff->scene;
628                         sim.ob= eff->ob;
629                         sim.psys= eff->psys;
630
631                         /* TODO: time from actual previous calculated frame (step might not be 1) */
632                         state.time = cfra - 1.0f;
633                         ret = psys_get_particle_state(&sim, *efd->index, &state, 0);
634
635                         /* TODO */
636                         //if (eff->pd->forcefiled == PFIELD_HARMONIC && ret==0) {
637                         //      if (pa->dietime < eff->psys->cfra)
638                         //              eff->flag |= PE_VELOCITY_TO_IMPULSE;
639                         //}
640
641                         copy_v3_v3(efd->loc, state.co);
642
643                         /* rather than use the velocity use rotated x-axis (defaults to velocity) */
644                         efd->nor[0] = 1.f;
645                         efd->nor[1] = efd->nor[2] = 0.f;
646                         mul_qt_v3(state.rot, efd->nor);
647                 
648                         if (real_velocity)
649                                 copy_v3_v3(efd->vel, state.vel);
650
651                         efd->size = pa->size;
652                 }
653         }
654         else {
655                 /* use center of object for distance calculus */
656                 const Object *ob = eff->ob;
657
658                 /* use z-axis as normal*/
659                 normalize_v3_v3(efd->nor, ob->obmat[2]);
660
661                 if (eff->pd && eff->pd->shape == PFIELD_SHAPE_PLANE) {
662                         float temp[3], translate[3];
663                         sub_v3_v3v3(temp, point->loc, ob->obmat[3]);
664                         project_v3_v3v3(translate, temp, efd->nor);
665
666                         /* for vortex the shape chooses between old / new force */
667                         if (eff->pd->forcefield == PFIELD_VORTEX)
668                                 add_v3_v3v3(efd->loc, ob->obmat[3], translate);
669                         else /* normally efd->loc is closest point on effector xy-plane */
670                                 sub_v3_v3v3(efd->loc, point->loc, translate);
671                 }
672                 else {
673                         copy_v3_v3(efd->loc, ob->obmat[3]);
674                 }
675
676                 if (real_velocity)
677                         copy_v3_v3(efd->vel, eff->velocity);
678
679                 efd->size = 0.0f;
680
681                 ret = 1;
682         }
683
684         if (ret) {
685                 sub_v3_v3v3(efd->vec_to_point, point->loc, efd->loc);
686                 efd->distance = len_v3(efd->vec_to_point);
687
688                 /* rest length for harmonic effector, will have to see later if this could be extended to other effectors */
689                 if (eff->pd && eff->pd->forcefield == PFIELD_HARMONIC && eff->pd->f_size)
690                         mul_v3_fl(efd->vec_to_point, (efd->distance-eff->pd->f_size)/efd->distance);
691
692                 if (eff->flag & PE_USE_NORMAL_DATA) {
693                         copy_v3_v3(efd->vec_to_point2, efd->vec_to_point);
694                         copy_v3_v3(efd->nor2, efd->nor);
695                 }
696                 else {
697                         /* for some effectors we need the object center every time */
698                         sub_v3_v3v3(efd->vec_to_point2, point->loc, eff->ob->obmat[3]);
699                         normalize_v3_v3(efd->nor2, eff->ob->obmat[2]);
700                 }
701         }
702
703         return ret;
704 }
705 static void get_effector_tot(EffectorCache *eff, EffectorData *efd, EffectedPoint *point, int *tot, int *p, int *step)
706 {
707         *p = 0;
708         efd->index = p;
709
710         if (eff->pd->shape == PFIELD_SHAPE_POINTS) {
711                 *tot = eff->ob->derivedFinal ? eff->ob->derivedFinal->numVertData : 1;
712
713                 if (*tot && eff->pd->forcefield == PFIELD_HARMONIC && point->index >= 0) {
714                         *p = point->index % *tot;
715                         *tot = *p+1;
716                 }
717         }
718         else if (eff->psys) {
719                 *tot = eff->psys->totpart;
720                 
721                 if (eff->pd->forcefield == PFIELD_CHARGE) {
722                         /* Only the charge of the effected particle is used for 
723                          * interaction, not fall-offs. If the fall-offs aren't the
724                          * same this will be unphysical, but for animation this
725                          * could be the wanted behavior. If you want physical
726                          * correctness the fall-off should be spherical 2.0 anyways.
727                          */
728                         efd->charge = eff->pd->f_strength;
729                 }
730                 else if (eff->pd->forcefield == PFIELD_HARMONIC && (eff->pd->flag & PFIELD_MULTIPLE_SPRINGS)==0) {
731                         /* every particle is mapped to only one harmonic effector particle */
732                         *p= point->index % eff->psys->totpart;
733                         *tot= *p + 1;
734                 }
735
736                 if (eff->psys->part->effector_amount) {
737                         int totpart = eff->psys->totpart;
738                         int amount = eff->psys->part->effector_amount;
739
740                         *step = (totpart > amount) ? totpart/amount : 1;
741                 }
742         }
743         else {
744                 *tot = 1;
745         }
746 }
747 static void do_texture_effector(EffectorCache *eff, EffectorData *efd, EffectedPoint *point, float *total_force)
748 {
749         TexResult result[4];
750         float tex_co[3], strength, force[3];
751         float nabla = eff->pd->tex_nabla;
752         int hasrgb;
753         short mode = eff->pd->tex_mode;
754         bool scene_color_manage;
755
756         if (!eff->pd->tex)
757                 return;
758
759         result[0].nor = result[1].nor = result[2].nor = result[3].nor = NULL;
760
761         strength= eff->pd->f_strength * efd->falloff;
762
763         copy_v3_v3(tex_co, point->loc);
764
765         if (eff->pd->flag & PFIELD_TEX_OBJECT) {
766                 mul_m4_v3(eff->ob->imat, tex_co);
767
768                 if (eff->pd->flag & PFIELD_TEX_2D)
769                         tex_co[2] = 0.0f;
770         }
771         else if (eff->pd->flag & PFIELD_TEX_2D) {
772                 float fac=-dot_v3v3(tex_co, efd->nor);
773                 madd_v3_v3fl(tex_co, efd->nor, fac);
774         }
775
776         scene_color_manage = BKE_scene_check_color_management_enabled(eff->scene);
777
778         hasrgb = multitex_ext(eff->pd->tex, tex_co, NULL, NULL, 0, result, 0, NULL, scene_color_manage, false);
779
780         if (hasrgb && mode==PFIELD_TEX_RGB) {
781                 force[0] = (0.5f - result->tr) * strength;
782                 force[1] = (0.5f - result->tg) * strength;
783                 force[2] = (0.5f - result->tb) * strength;
784         }
785         else if (nabla != 0) {
786                 strength/=nabla;
787
788                 tex_co[0] += nabla;
789                 multitex_ext(eff->pd->tex, tex_co, NULL, NULL, 0, result+1, 0, NULL, scene_color_manage, false);
790
791                 tex_co[0] -= nabla;
792                 tex_co[1] += nabla;
793                 multitex_ext(eff->pd->tex, tex_co, NULL, NULL, 0, result+2, 0, NULL, scene_color_manage, false);
794
795                 tex_co[1] -= nabla;
796                 tex_co[2] += nabla;
797                 multitex_ext(eff->pd->tex, tex_co, NULL, NULL, 0, result+3, 0, NULL, scene_color_manage, false);
798
799                 if (mode == PFIELD_TEX_GRAD || !hasrgb) { /* if we don't have rgb fall back to grad */
800                         /* generate intensity if texture only has rgb value */
801                         if (hasrgb & TEX_RGB) {
802                                 int i;
803                                 for (i=0; i<4; i++)
804                                         result[i].tin = (1.0f / 3.0f) * (result[i].tr + result[i].tg + result[i].tb);
805                         }
806                         force[0] = (result[0].tin - result[1].tin) * strength;
807                         force[1] = (result[0].tin - result[2].tin) * strength;
808                         force[2] = (result[0].tin - result[3].tin) * strength;
809                 }
810                 else { /*PFIELD_TEX_CURL*/
811                         float dbdy, dgdz, drdz, dbdx, dgdx, drdy;
812
813                         dbdy = result[2].tb - result[0].tb;
814                         dgdz = result[3].tg - result[0].tg;
815                         drdz = result[3].tr - result[0].tr;
816                         dbdx = result[1].tb - result[0].tb;
817                         dgdx = result[1].tg - result[0].tg;
818                         drdy = result[2].tr - result[0].tr;
819
820                         force[0] = (dbdy - dgdz) * strength;
821                         force[1] = (drdz - dbdx) * strength;
822                         force[2] = (dgdx - drdy) * strength;
823                 }
824         }
825         else {
826                 zero_v3(force);
827         }
828
829         if (eff->pd->flag & PFIELD_TEX_2D) {
830                 float fac = -dot_v3v3(force, efd->nor);
831                 madd_v3_v3fl(force, efd->nor, fac);
832         }
833
834         add_v3_v3(total_force, force);
835 }
836 static void do_physical_effector(EffectorCache *eff, EffectorData *efd, EffectedPoint *point, float *total_force)
837 {
838         PartDeflect *pd = eff->pd;
839         RNG *rng = pd->rng;
840         float force[3] = {0, 0, 0};
841         float temp[3];
842         float fac;
843         float strength = pd->f_strength;
844         float damp = pd->f_damp;
845         float noise_factor = pd->f_noise;
846
847         if (noise_factor > 0.0f) {
848                 strength += wind_func(rng, noise_factor);
849
850                 if (ELEM(pd->forcefield, PFIELD_HARMONIC, PFIELD_DRAG))
851                         damp += wind_func(rng, noise_factor);
852         }
853
854         copy_v3_v3(force, efd->vec_to_point);
855
856         switch (pd->forcefield) {
857                 case PFIELD_WIND:
858                         copy_v3_v3(force, efd->nor);
859                         mul_v3_fl(force, strength * efd->falloff);
860                         break;
861                 case PFIELD_FORCE:
862                         normalize_v3(force);
863                         if (pd->flag & PFIELD_GRAVITATION){ /* Option: Multiply by 1/distance^2 */
864                                 if (efd->distance < FLT_EPSILON){
865                                         strength = 0.0f;
866                                 }
867                                 else {
868                                         strength *= powf(efd->distance, -2.0f);
869                                 }
870                         }
871                         mul_v3_fl(force, strength * efd->falloff);
872                         break;
873                 case PFIELD_VORTEX:
874                         /* old vortex force */
875                         if (pd->shape == PFIELD_SHAPE_POINT) {
876                                 cross_v3_v3v3(force, efd->nor, efd->vec_to_point);
877                                 normalize_v3(force);
878                                 mul_v3_fl(force, strength * efd->distance * efd->falloff);
879                         }
880                         else {
881                                 /* new vortex force */
882                                 cross_v3_v3v3(temp, efd->nor2, efd->vec_to_point2);
883                                 mul_v3_fl(temp, strength * efd->falloff);
884                                 
885                                 cross_v3_v3v3(force, efd->nor2, temp);
886                                 mul_v3_fl(force, strength * efd->falloff);
887                                 
888                                 madd_v3_v3fl(temp, point->vel, -point->vel_to_sec);
889                                 add_v3_v3(force, temp);
890                         }
891                         break;
892                 case PFIELD_MAGNET:
893                         if (eff->pd->shape == PFIELD_SHAPE_POINT)
894                                 /* magnetic field of a moving charge */
895                                 cross_v3_v3v3(temp, efd->nor, efd->vec_to_point);
896                         else
897                                 copy_v3_v3(temp, efd->nor);
898
899                         normalize_v3(temp);
900                         mul_v3_fl(temp, strength * efd->falloff);
901                         cross_v3_v3v3(force, point->vel, temp);
902                         mul_v3_fl(force, point->vel_to_sec);
903                         break;
904                 case PFIELD_HARMONIC:
905                         mul_v3_fl(force, -strength * efd->falloff);
906                         copy_v3_v3(temp, point->vel);
907                         mul_v3_fl(temp, -damp * 2.0f * sqrtf(fabsf(strength)) * point->vel_to_sec);
908                         add_v3_v3(force, temp);
909                         break;
910                 case PFIELD_CHARGE:
911                         mul_v3_fl(force, point->charge * strength * efd->falloff);
912                         break;
913                 case PFIELD_LENNARDJ:
914                         fac = pow((efd->size + point->size) / efd->distance, 6.0);
915                         
916                         fac = - fac * (1.0f - fac) / efd->distance;
917
918                         /* limit the repulsive term drastically to avoid huge forces */
919                         fac = ((fac>2.0f) ? 2.0f : fac);
920
921                         mul_v3_fl(force, strength * fac);
922                         break;
923                 case PFIELD_BOID:
924                         /* Boid field is handled completely in boids code. */
925                         return;
926                 case PFIELD_TURBULENCE:
927                         if (pd->flag & PFIELD_GLOBAL_CO) {
928                                 copy_v3_v3(temp, point->loc);
929                         }
930                         else {
931                                 add_v3_v3v3(temp, efd->vec_to_point2, efd->nor2);
932                         }
933                         force[0] = -1.0f + 2.0f * BLI_gTurbulence(pd->f_size, temp[0], temp[1], temp[2], 2, 0, 2);
934                         force[1] = -1.0f + 2.0f * BLI_gTurbulence(pd->f_size, temp[1], temp[2], temp[0], 2, 0, 2);
935                         force[2] = -1.0f + 2.0f * BLI_gTurbulence(pd->f_size, temp[2], temp[0], temp[1], 2, 0, 2);
936                         mul_v3_fl(force, strength * efd->falloff);
937                         break;
938                 case PFIELD_DRAG:
939                         copy_v3_v3(force, point->vel);
940                         fac = normalize_v3(force) * point->vel_to_sec;
941
942                         strength = MIN2(strength, 2.0f);
943                         damp = MIN2(damp, 2.0f);
944
945                         mul_v3_fl(force, -efd->falloff * fac * (strength * fac + damp));
946                         break;
947                 case PFIELD_SMOKEFLOW:
948                         zero_v3(force);
949                         if (pd->f_source) {
950                                 float density;
951                                 if ((density = smoke_get_velocity_at(pd->f_source, point->loc, force)) >= 0.0f) {
952                                         float influence = strength * efd->falloff;
953                                         if (pd->flag & PFIELD_SMOKE_DENSITY)
954                                                 influence *= density;
955                                         mul_v3_fl(force, influence);
956                                         /* apply flow */
957                                         madd_v3_v3fl(total_force, point->vel, -pd->f_flow * influence);
958                                 }
959                         }
960                         break;
961
962         }
963
964         if (pd->flag & PFIELD_DO_LOCATION) {
965                 madd_v3_v3fl(total_force, force, 1.0f/point->vel_to_sec);
966
967                 if (ELEM(pd->forcefield, PFIELD_HARMONIC, PFIELD_DRAG, PFIELD_SMOKEFLOW)==0 && pd->f_flow != 0.0f) {
968                         madd_v3_v3fl(total_force, point->vel, -pd->f_flow * efd->falloff);
969                 }
970         }
971
972         if (point->ave)
973                 zero_v3(point->ave);
974         if (pd->flag & PFIELD_DO_ROTATION && point->ave && point->rot) {
975                 float xvec[3] = {1.0f, 0.0f, 0.0f};
976                 float dave[3];
977                 mul_qt_v3(point->rot, xvec);
978                 cross_v3_v3v3(dave, xvec, force);
979                 if (pd->f_flow != 0.0f) {
980                         madd_v3_v3fl(dave, point->ave, -pd->f_flow * efd->falloff);
981                 }
982                 add_v3_v3(point->ave, dave);
983         }
984 }
985
986 /*  -------- pdDoEffectors() --------
987  * generic force/speed system, now used for particles and softbodies
988  * scene       = scene where it runs in, for time and stuff
989  * lb                   = listbase with objects that take part in effecting
990  * opco         = global coord, as input
991  * force                = force accumulator
992  * speed                = actual current speed which can be altered
993  * cur_time     = "external" time in frames, is constant for static particles
994  * loc_time     = "local" time in frames, range <0-1> for the lifetime of particle
995  * par_layer    = layer the caller is in
996  * flags                = only used for softbody wind now
997  * guide                = old speed of particle
998  */
999 void pdDoEffectors(ListBase *effectors, ListBase *colliders, EffectorWeights *weights, EffectedPoint *point, float *force, float *impulse)
1000 {
1001         /*
1002          * Modifies the force on a particle according to its
1003          * relation with the effector object
1004          * Different kind of effectors include:
1005          *     Forcefields: Gravity-like attractor
1006          *     (force power is related to the inverse of distance to the power of a falloff value)
1007          *     Vortex fields: swirling effectors
1008          *     (particles rotate around Z-axis of the object. otherwise, same relation as)
1009          *     (Forcefields, but this is not done through a force/acceleration)
1010          *     Guide: particles on a path
1011          *     (particles are guided along a curve bezier or old nurbs)
1012          *     (is independent of other effectors)
1013          */
1014         EffectorCache *eff;
1015         EffectorData efd;
1016         int p=0, tot = 1, step = 1;
1017
1018         /* Cycle through collected objects, get total of (1/(gravity_strength * dist^gravity_power)) */
1019         /* Check for min distance here? (yes would be cool to add that, ton) */
1020         
1021         if (effectors) for (eff = effectors->first; eff; eff=eff->next) {
1022                 /* object effectors were fully checked to be OK to evaluate! */
1023
1024                 get_effector_tot(eff, &efd, point, &tot, &p, &step);
1025
1026                 for (; p<tot; p+=step) {
1027                         if (get_effector_data(eff, &efd, point, 0)) {
1028                                 efd.falloff= effector_falloff(eff, &efd, point, weights);
1029                                 
1030                                 if (efd.falloff > 0.0f)
1031                                         efd.falloff *= eff_calc_visibility(colliders, eff, &efd, point);
1032
1033                                 if (efd.falloff <= 0.0f) {
1034                                         /* don't do anything */
1035                                 }
1036                                 else if (eff->pd->forcefield == PFIELD_TEXTURE) {
1037                                         do_texture_effector(eff, &efd, point, force);
1038                                 }
1039                                 else {
1040                                         float temp1[3] = {0, 0, 0}, temp2[3];
1041                                         copy_v3_v3(temp1, force);
1042
1043                                         do_physical_effector(eff, &efd, point, force);
1044                                         
1045                                         /* for softbody backward compatibility */
1046                                         if (point->flag & PE_WIND_AS_SPEED && impulse) {
1047                                                 sub_v3_v3v3(temp2, force, temp1);
1048                                                 sub_v3_v3v3(impulse, impulse, temp2);
1049                                         }
1050                                 }
1051                         }
1052                         else if (eff->flag & PE_VELOCITY_TO_IMPULSE && impulse) {
1053                                 /* special case for harmonic effector */
1054                                 add_v3_v3v3(impulse, impulse, efd.vel);
1055                         }
1056                 }
1057         }
1058 }
1059
1060 /* ======== Simulation Debugging ======== */
1061
1062 SimDebugData *_sim_debug_data = NULL;
1063
1064 unsigned int BKE_sim_debug_data_hash(int i)
1065 {
1066         return BLI_ghashutil_uinthash((unsigned int)i);
1067 }
1068
1069 unsigned int BKE_sim_debug_data_hash_combine(unsigned int kx, unsigned int ky)
1070 {
1071 #define rot(x,k) (((x)<<(k)) | ((x)>>(32-(k))))
1072
1073         unsigned int a, b, c;
1074
1075         a = b = c = 0xdeadbeef + (2 << 2) + 13;
1076         a += kx;
1077         b += ky;
1078
1079         c ^= b; c -= rot(b,14);
1080         a ^= c; a -= rot(c,11);
1081         b ^= a; b -= rot(a,25);
1082         c ^= b; c -= rot(b,16);
1083         a ^= c; a -= rot(c,4);
1084         b ^= a; b -= rot(a,14);
1085         c ^= b; c -= rot(b,24);
1086
1087         return c;
1088
1089 #undef rot
1090 }
1091
1092 static unsigned int debug_element_hash(const void *key)
1093 {
1094         const SimDebugElement *elem = key;
1095         return elem->hash;
1096 }
1097
1098 static bool debug_element_compare(const void *a, const void *b)
1099 {
1100         const SimDebugElement *elem1 = a;
1101         const SimDebugElement *elem2 = b;
1102
1103         if (elem1->hash == elem2->hash) {
1104                 return 0;
1105         }
1106         return 1;
1107 }
1108
1109 static void debug_element_free(void *val)
1110 {
1111         SimDebugElement *elem = val;
1112         MEM_freeN(elem);
1113 }
1114
1115 void BKE_sim_debug_data_set_enabled(bool enable)
1116 {
1117         if (enable) {
1118                 if (!_sim_debug_data) {
1119                         _sim_debug_data = MEM_callocN(sizeof(SimDebugData), "sim debug data");
1120                         _sim_debug_data->gh = BLI_ghash_new(debug_element_hash, debug_element_compare, "sim debug element hash");
1121                 }
1122         }
1123         else {
1124                 BKE_sim_debug_data_free();
1125         }
1126 }
1127
1128 bool BKE_sim_debug_data_get_enabled(void)
1129 {
1130         return _sim_debug_data != NULL;
1131 }
1132
1133 void BKE_sim_debug_data_free(void)
1134 {
1135         if (_sim_debug_data) {
1136                 if (_sim_debug_data->gh)
1137                         BLI_ghash_free(_sim_debug_data->gh, NULL, debug_element_free);
1138                 MEM_freeN(_sim_debug_data);
1139         }
1140 }
1141
1142 static void debug_data_insert(SimDebugData *debug_data, SimDebugElement *elem)
1143 {
1144         SimDebugElement *old_elem = BLI_ghash_lookup(debug_data->gh, elem);
1145         if (old_elem) {
1146                 *old_elem = *elem;
1147                 MEM_freeN(elem);
1148         }
1149         else
1150                 BLI_ghash_insert(debug_data->gh, elem, elem);
1151 }
1152
1153 void BKE_sim_debug_data_add_element(int type, const float v1[3], const float v2[3], const char *str, float r, float g, float b, const char *category, unsigned int hash)
1154 {
1155         unsigned int category_hash = BLI_ghashutil_strhash_p(category);
1156         SimDebugElement *elem;
1157         
1158         if (!_sim_debug_data) {
1159                 if (G.debug & G_DEBUG_SIMDATA)
1160                         BKE_sim_debug_data_set_enabled(true);
1161                 else
1162                         return;
1163         }
1164         
1165         elem = MEM_callocN(sizeof(SimDebugElement), "sim debug data element");
1166         elem->type = type;
1167         elem->category_hash = category_hash;
1168         elem->hash = hash;
1169         elem->color[0] = r;
1170         elem->color[1] = g;
1171         elem->color[2] = b;
1172         if (v1)
1173                 copy_v3_v3(elem->v1, v1);
1174         else
1175                 zero_v3(elem->v1);
1176         if (v2)
1177                 copy_v3_v3(elem->v2, v2);
1178         else
1179                 zero_v3(elem->v2);
1180         if (str)
1181                 BLI_strncpy(elem->str, str, sizeof(elem->str));
1182         else
1183                 elem->str[0] = '\0';
1184         
1185         debug_data_insert(_sim_debug_data, elem);
1186 }
1187
1188 void BKE_sim_debug_data_remove_element(unsigned int hash)
1189 {
1190         SimDebugElement dummy;
1191         if (!_sim_debug_data)
1192                 return;
1193         
1194         dummy.hash = hash;
1195         BLI_ghash_remove(_sim_debug_data->gh, &dummy, NULL, debug_element_free);
1196 }
1197
1198 void BKE_sim_debug_data_clear(void)
1199 {
1200         if (!_sim_debug_data)
1201                 return;
1202         
1203         if (_sim_debug_data->gh)
1204                 BLI_ghash_clear(_sim_debug_data->gh, NULL, debug_element_free);
1205 }
1206
1207 void BKE_sim_debug_data_clear_category(const char *category)
1208 {
1209         int category_hash = (int)BLI_ghashutil_strhash_p(category);
1210         
1211         if (!_sim_debug_data)
1212                 return;
1213         
1214         if (_sim_debug_data->gh) {
1215                 GHashIterator iter;
1216                 BLI_ghashIterator_init(&iter, _sim_debug_data->gh);
1217                 while (!BLI_ghashIterator_done(&iter)) {
1218                         const SimDebugElement *elem = BLI_ghashIterator_getValue(&iter);
1219                         BLI_ghashIterator_step(&iter); /* removing invalidates the current iterator, so step before removing */
1220                         
1221                         if (elem->category_hash == category_hash)
1222                                 BLI_ghash_remove(_sim_debug_data->gh, elem, NULL, debug_element_free);
1223                 }
1224         }
1225 }