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