Cleanup: comments (long lines) in blenkernel
[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,
650    * surface modifier is never executed 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,
754      * will have to see later if this could be extended to other effectors. */
755     if (eff->pd && eff->pd->forcefield == PFIELD_HARMONIC && eff->pd->f_size) {
756       mul_v3_fl(efd->vec_to_point, (efd->distance - eff->pd->f_size) / efd->distance);
757     }
758
759     if (eff->flag & PE_USE_NORMAL_DATA) {
760       copy_v3_v3(efd->vec_to_point2, efd->vec_to_point);
761       copy_v3_v3(efd->nor2, efd->nor);
762     }
763     else {
764       /* for some effectors we need the object center every time */
765       sub_v3_v3v3(efd->vec_to_point2, point->loc, eff->ob->obmat[3]);
766       normalize_v3_v3(efd->nor2, eff->ob->obmat[2]);
767     }
768   }
769
770   return ret;
771 }
772 static void get_effector_tot(
773     EffectorCache *eff, EffectorData *efd, EffectedPoint *point, int *tot, int *p, int *step)
774 {
775   *p = 0;
776   efd->index = p;
777
778   if (eff->pd->shape == PFIELD_SHAPE_POINTS) {
779     Mesh *me_eval = eff->ob->runtime.mesh_eval;
780     *tot = me_eval != NULL ? me_eval->totvert : 1;
781
782     if (*tot && eff->pd->forcefield == PFIELD_HARMONIC && point->index >= 0) {
783       *p = point->index % *tot;
784       *tot = *p + 1;
785     }
786   }
787   else if (eff->psys) {
788     *tot = eff->psys->totpart;
789
790     if (eff->pd->forcefield == PFIELD_CHARGE) {
791       /* Only the charge of the effected particle is used for
792        * interaction, not fall-offs. If the fall-offs aren't the
793        * same this will be unphysical, but for animation this
794        * could be the wanted behavior. If you want physical
795        * correctness the fall-off should be spherical 2.0 anyways.
796        */
797       efd->charge = eff->pd->f_strength;
798     }
799     else if (eff->pd->forcefield == PFIELD_HARMONIC &&
800              (eff->pd->flag & PFIELD_MULTIPLE_SPRINGS) == 0) {
801       /* every particle is mapped to only one harmonic effector particle */
802       *p = point->index % eff->psys->totpart;
803       *tot = *p + 1;
804     }
805
806     if (eff->psys->part->effector_amount) {
807       int totpart = eff->psys->totpart;
808       int amount = eff->psys->part->effector_amount;
809
810       *step = (totpart > amount) ? totpart / amount : 1;
811     }
812   }
813   else {
814     *tot = 1;
815   }
816 }
817 static void do_texture_effector(EffectorCache *eff,
818                                 EffectorData *efd,
819                                 EffectedPoint *point,
820                                 float *total_force)
821 {
822   TexResult result[4];
823   float tex_co[3], strength, force[3];
824   float nabla = eff->pd->tex_nabla;
825   int hasrgb;
826   short mode = eff->pd->tex_mode;
827   bool scene_color_manage;
828
829   if (!eff->pd->tex) {
830     return;
831   }
832
833   result[0].nor = result[1].nor = result[2].nor = result[3].nor = NULL;
834
835   strength = eff->pd->f_strength * efd->falloff;
836
837   copy_v3_v3(tex_co, point->loc);
838
839   if (eff->pd->flag & PFIELD_TEX_OBJECT) {
840     mul_m4_v3(eff->ob->imat, tex_co);
841
842     if (eff->pd->flag & PFIELD_TEX_2D) {
843       tex_co[2] = 0.0f;
844     }
845   }
846   else if (eff->pd->flag & PFIELD_TEX_2D) {
847     float fac = -dot_v3v3(tex_co, efd->nor);
848     madd_v3_v3fl(tex_co, efd->nor, fac);
849   }
850
851   scene_color_manage = BKE_scene_check_color_management_enabled(eff->scene);
852
853   hasrgb = multitex_ext(
854       eff->pd->tex, tex_co, NULL, NULL, 0, result, 0, NULL, scene_color_manage, false);
855
856   if (hasrgb && mode == PFIELD_TEX_RGB) {
857     force[0] = (0.5f - result->tr) * strength;
858     force[1] = (0.5f - result->tg) * strength;
859     force[2] = (0.5f - result->tb) * strength;
860   }
861   else if (nabla != 0) {
862     strength /= nabla;
863
864     tex_co[0] += nabla;
865     multitex_ext(
866         eff->pd->tex, tex_co, NULL, NULL, 0, result + 1, 0, NULL, scene_color_manage, false);
867
868     tex_co[0] -= nabla;
869     tex_co[1] += nabla;
870     multitex_ext(
871         eff->pd->tex, tex_co, NULL, NULL, 0, result + 2, 0, NULL, scene_color_manage, false);
872
873     tex_co[1] -= nabla;
874     tex_co[2] += nabla;
875     multitex_ext(
876         eff->pd->tex, tex_co, NULL, NULL, 0, result + 3, 0, NULL, scene_color_manage, false);
877
878     if (mode == PFIELD_TEX_GRAD || !hasrgb) { /* if we don't have rgb fall back to grad */
879       /* generate intensity if texture only has rgb value */
880       if (hasrgb & TEX_RGB) {
881         for (int i = 0; i < 4; i++) {
882           result[i].tin = (1.0f / 3.0f) * (result[i].tr + result[i].tg + result[i].tb);
883         }
884       }
885       force[0] = (result[0].tin - result[1].tin) * strength;
886       force[1] = (result[0].tin - result[2].tin) * strength;
887       force[2] = (result[0].tin - result[3].tin) * strength;
888     }
889     else { /*PFIELD_TEX_CURL*/
890       float dbdy, dgdz, drdz, dbdx, dgdx, drdy;
891
892       dbdy = result[2].tb - result[0].tb;
893       dgdz = result[3].tg - result[0].tg;
894       drdz = result[3].tr - result[0].tr;
895       dbdx = result[1].tb - result[0].tb;
896       dgdx = result[1].tg - result[0].tg;
897       drdy = result[2].tr - result[0].tr;
898
899       force[0] = (dbdy - dgdz) * strength;
900       force[1] = (drdz - dbdx) * strength;
901       force[2] = (dgdx - drdy) * strength;
902     }
903   }
904   else {
905     zero_v3(force);
906   }
907
908   if (eff->pd->flag & PFIELD_TEX_2D) {
909     float fac = -dot_v3v3(force, efd->nor);
910     madd_v3_v3fl(force, efd->nor, fac);
911   }
912
913   add_v3_v3(total_force, force);
914 }
915 static void do_physical_effector(EffectorCache *eff,
916                                  EffectorData *efd,
917                                  EffectedPoint *point,
918                                  float *total_force)
919 {
920   PartDeflect *pd = eff->pd;
921   RNG *rng = pd->rng;
922   float force[3] = {0, 0, 0};
923   float temp[3];
924   float fac;
925   float strength = pd->f_strength;
926   float damp = pd->f_damp;
927   float noise_factor = pd->f_noise;
928
929   if (noise_factor > 0.0f) {
930     strength += wind_func(rng, noise_factor);
931
932     if (ELEM(pd->forcefield, PFIELD_HARMONIC, PFIELD_DRAG)) {
933       damp += wind_func(rng, noise_factor);
934     }
935   }
936
937   copy_v3_v3(force, efd->vec_to_point);
938
939   switch (pd->forcefield) {
940     case PFIELD_WIND:
941       copy_v3_v3(force, efd->nor);
942       mul_v3_fl(force, strength * efd->falloff);
943       break;
944     case PFIELD_FORCE:
945       normalize_v3(force);
946       if (pd->flag & PFIELD_GRAVITATION) { /* Option: Multiply by 1/distance^2 */
947         if (efd->distance < FLT_EPSILON) {
948           strength = 0.0f;
949         }
950         else {
951           strength *= powf(efd->distance, -2.0f);
952         }
953       }
954       mul_v3_fl(force, strength * efd->falloff);
955       break;
956     case PFIELD_VORTEX:
957       /* old vortex force */
958       if (pd->shape == PFIELD_SHAPE_POINT) {
959         cross_v3_v3v3(force, efd->nor, efd->vec_to_point);
960         normalize_v3(force);
961         mul_v3_fl(force, strength * efd->distance * efd->falloff);
962       }
963       else {
964         /* new vortex force */
965         cross_v3_v3v3(temp, efd->nor2, efd->vec_to_point2);
966         mul_v3_fl(temp, strength * efd->falloff);
967
968         cross_v3_v3v3(force, efd->nor2, temp);
969         mul_v3_fl(force, strength * efd->falloff);
970
971         madd_v3_v3fl(temp, point->vel, -point->vel_to_sec);
972         add_v3_v3(force, temp);
973       }
974       break;
975     case PFIELD_MAGNET:
976       if (ELEM(eff->pd->shape, PFIELD_SHAPE_POINT, PFIELD_SHAPE_LINE)) {
977         /* magnetic field of a moving charge */
978         cross_v3_v3v3(temp, efd->nor, efd->vec_to_point);
979       }
980       else {
981         copy_v3_v3(temp, efd->nor);
982       }
983
984       normalize_v3(temp);
985       mul_v3_fl(temp, strength * efd->falloff);
986       cross_v3_v3v3(force, point->vel, temp);
987       mul_v3_fl(force, point->vel_to_sec);
988       break;
989     case PFIELD_HARMONIC:
990       mul_v3_fl(force, -strength * efd->falloff);
991       copy_v3_v3(temp, point->vel);
992       mul_v3_fl(temp, -damp * 2.0f * sqrtf(fabsf(strength)) * point->vel_to_sec);
993       add_v3_v3(force, temp);
994       break;
995     case PFIELD_CHARGE:
996       mul_v3_fl(force, point->charge * strength * efd->falloff);
997       break;
998     case PFIELD_LENNARDJ:
999       fac = pow((efd->size + point->size) / efd->distance, 6.0);
1000
1001       fac = -fac * (1.0f - fac) / efd->distance;
1002
1003       /* limit the repulsive term drastically to avoid huge forces */
1004       fac = ((fac > 2.0f) ? 2.0f : fac);
1005
1006       mul_v3_fl(force, strength * fac);
1007       break;
1008     case PFIELD_BOID:
1009       /* Boid field is handled completely in boids code. */
1010       return;
1011     case PFIELD_TURBULENCE:
1012       if (pd->flag & PFIELD_GLOBAL_CO) {
1013         copy_v3_v3(temp, point->loc);
1014       }
1015       else {
1016         add_v3_v3v3(temp, efd->vec_to_point2, efd->nor2);
1017       }
1018       force[0] = -1.0f + 2.0f * BLI_gTurbulence(pd->f_size, temp[0], temp[1], temp[2], 2, 0, 2);
1019       force[1] = -1.0f + 2.0f * BLI_gTurbulence(pd->f_size, temp[1], temp[2], temp[0], 2, 0, 2);
1020       force[2] = -1.0f + 2.0f * BLI_gTurbulence(pd->f_size, temp[2], temp[0], temp[1], 2, 0, 2);
1021       mul_v3_fl(force, strength * efd->falloff);
1022       break;
1023     case PFIELD_DRAG:
1024       copy_v3_v3(force, point->vel);
1025       fac = normalize_v3(force) * point->vel_to_sec;
1026
1027       strength = MIN2(strength, 2.0f);
1028       damp = MIN2(damp, 2.0f);
1029
1030       mul_v3_fl(force, -efd->falloff * fac * (strength * fac + damp));
1031       break;
1032     case PFIELD_SMOKEFLOW:
1033       zero_v3(force);
1034       if (pd->f_source) {
1035         float density;
1036         if ((density = BKE_smoke_get_velocity_at(pd->f_source, point->loc, force)) >= 0.0f) {
1037           float influence = strength * efd->falloff;
1038           if (pd->flag & PFIELD_SMOKE_DENSITY) {
1039             influence *= density;
1040           }
1041           mul_v3_fl(force, influence);
1042           /* apply flow */
1043           madd_v3_v3fl(total_force, point->vel, -pd->f_flow * influence);
1044         }
1045       }
1046       break;
1047   }
1048
1049   if (pd->flag & PFIELD_DO_LOCATION) {
1050     madd_v3_v3fl(total_force, force, 1.0f / point->vel_to_sec);
1051
1052     if (ELEM(pd->forcefield, PFIELD_HARMONIC, PFIELD_DRAG, PFIELD_SMOKEFLOW) == 0 &&
1053         pd->f_flow != 0.0f) {
1054       madd_v3_v3fl(total_force, point->vel, -pd->f_flow * efd->falloff);
1055     }
1056   }
1057
1058   if (point->ave) {
1059     zero_v3(point->ave);
1060   }
1061   if (pd->flag & PFIELD_DO_ROTATION && point->ave && point->rot) {
1062     float xvec[3] = {1.0f, 0.0f, 0.0f};
1063     float dave[3];
1064     mul_qt_v3(point->rot, xvec);
1065     cross_v3_v3v3(dave, xvec, force);
1066     if (pd->f_flow != 0.0f) {
1067       madd_v3_v3fl(dave, point->ave, -pd->f_flow * efd->falloff);
1068     }
1069     add_v3_v3(point->ave, dave);
1070   }
1071 }
1072
1073 /*  -------- BKE_effectors_apply() --------
1074  * generic force/speed system, now used for particles and softbodies
1075  * scene       = scene where it runs in, for time and stuff
1076  * lb           = listbase with objects that take part in effecting
1077  * opco     = global coord, as input
1078  * force        = force accumulator
1079  * speed        = actual current speed which can be altered
1080  * cur_time = "external" time in frames, is constant for static particles
1081  * loc_time = "local" time in frames, range <0-1> for the lifetime of particle
1082  * par_layer    = layer the caller is in
1083  * flags        = only used for softbody wind now
1084  * guide        = old speed of particle
1085  */
1086 void BKE_effectors_apply(ListBase *effectors,
1087                          ListBase *colliders,
1088                          EffectorWeights *weights,
1089                          EffectedPoint *point,
1090                          float *force,
1091                          float *impulse)
1092 {
1093   /*
1094    * Modifies the force on a particle according to its
1095    * relation with the effector object
1096    * Different kind of effectors include:
1097    *     Forcefields: Gravity-like attractor
1098    *     (force power is related to the inverse of distance to the power of a falloff value)
1099    *     Vortex fields: swirling effectors
1100    *     (particles rotate around Z-axis of the object. otherwise, same relation as)
1101    *     (Forcefields, but this is not done through a force/acceleration)
1102    *     Guide: particles on a path
1103    *     (particles are guided along a curve bezier or old nurbs)
1104    *     (is independent of other effectors)
1105    */
1106   EffectorCache *eff;
1107   EffectorData efd;
1108   int p = 0, tot = 1, step = 1;
1109
1110   /* Cycle through collected objects, get total of (1/(gravity_strength * dist^gravity_power)) */
1111   /* Check for min distance here? (yes would be cool to add that, ton) */
1112
1113   if (effectors) {
1114     for (eff = effectors->first; eff; eff = eff->next) {
1115       /* object effectors were fully checked to be OK to evaluate! */
1116
1117       get_effector_tot(eff, &efd, point, &tot, &p, &step);
1118
1119       for (; p < tot; p += step) {
1120         if (get_effector_data(eff, &efd, point, 0)) {
1121           efd.falloff = effector_falloff(eff, &efd, point, weights);
1122
1123           if (efd.falloff > 0.0f) {
1124             efd.falloff *= eff_calc_visibility(colliders, eff, &efd, point);
1125           }
1126           if (efd.falloff <= 0.0f) {
1127             /* don't do anything */
1128           }
1129           else if (eff->pd->forcefield == PFIELD_TEXTURE) {
1130             do_texture_effector(eff, &efd, point, force);
1131           }
1132           else {
1133             float temp1[3] = {0, 0, 0}, temp2[3];
1134             copy_v3_v3(temp1, force);
1135
1136             do_physical_effector(eff, &efd, point, force);
1137
1138             /* for softbody backward compatibility */
1139             if (point->flag & PE_WIND_AS_SPEED && impulse) {
1140               sub_v3_v3v3(temp2, force, temp1);
1141               sub_v3_v3v3(impulse, impulse, temp2);
1142             }
1143           }
1144         }
1145         else if (eff->flag & PE_VELOCITY_TO_IMPULSE && impulse) {
1146           /* special case for harmonic effector */
1147           add_v3_v3v3(impulse, impulse, efd.vel);
1148         }
1149       }
1150     }
1151   }
1152 }
1153
1154 /* ======== Simulation Debugging ======== */
1155
1156 SimDebugData *_sim_debug_data = NULL;
1157
1158 uint BKE_sim_debug_data_hash(int i)
1159 {
1160   return BLI_ghashutil_uinthash((uint)i);
1161 }
1162
1163 uint BKE_sim_debug_data_hash_combine(uint kx, uint ky)
1164 {
1165 #define rot(x, k) (((x) << (k)) | ((x) >> (32 - (k))))
1166
1167   uint a, b, c;
1168
1169   a = b = c = 0xdeadbeef + (2 << 2) + 13;
1170   a += kx;
1171   b += ky;
1172
1173   c ^= b;
1174   c -= rot(b, 14);
1175   a ^= c;
1176   a -= rot(c, 11);
1177   b ^= a;
1178   b -= rot(a, 25);
1179   c ^= b;
1180   c -= rot(b, 16);
1181   a ^= c;
1182   a -= rot(c, 4);
1183   b ^= a;
1184   b -= rot(a, 14);
1185   c ^= b;
1186   c -= rot(b, 24);
1187
1188   return c;
1189
1190 #undef rot
1191 }
1192
1193 static uint debug_element_hash(const void *key)
1194 {
1195   const SimDebugElement *elem = key;
1196   return elem->hash;
1197 }
1198
1199 static bool debug_element_compare(const void *a, const void *b)
1200 {
1201   const SimDebugElement *elem1 = a;
1202   const SimDebugElement *elem2 = b;
1203
1204   if (elem1->hash == elem2->hash) {
1205     return 0;
1206   }
1207   return 1;
1208 }
1209
1210 static void debug_element_free(void *val)
1211 {
1212   SimDebugElement *elem = val;
1213   MEM_freeN(elem);
1214 }
1215
1216 void BKE_sim_debug_data_set_enabled(bool enable)
1217 {
1218   if (enable) {
1219     if (!_sim_debug_data) {
1220       _sim_debug_data = MEM_callocN(sizeof(SimDebugData), "sim debug data");
1221       _sim_debug_data->gh = BLI_ghash_new(
1222           debug_element_hash, debug_element_compare, "sim debug element hash");
1223     }
1224   }
1225   else {
1226     BKE_sim_debug_data_free();
1227   }
1228 }
1229
1230 bool BKE_sim_debug_data_get_enabled(void)
1231 {
1232   return _sim_debug_data != NULL;
1233 }
1234
1235 void BKE_sim_debug_data_free(void)
1236 {
1237   if (_sim_debug_data) {
1238     if (_sim_debug_data->gh) {
1239       BLI_ghash_free(_sim_debug_data->gh, NULL, debug_element_free);
1240     }
1241     MEM_freeN(_sim_debug_data);
1242   }
1243 }
1244
1245 static void debug_data_insert(SimDebugData *debug_data, SimDebugElement *elem)
1246 {
1247   SimDebugElement *old_elem = BLI_ghash_lookup(debug_data->gh, elem);
1248   if (old_elem) {
1249     *old_elem = *elem;
1250     MEM_freeN(elem);
1251   }
1252   else {
1253     BLI_ghash_insert(debug_data->gh, elem, elem);
1254   }
1255 }
1256
1257 void BKE_sim_debug_data_add_element(int type,
1258                                     const float v1[3],
1259                                     const float v2[3],
1260                                     const char *str,
1261                                     float r,
1262                                     float g,
1263                                     float b,
1264                                     const char *category,
1265                                     uint hash)
1266 {
1267   uint category_hash = BLI_ghashutil_strhash_p(category);
1268   SimDebugElement *elem;
1269
1270   if (!_sim_debug_data) {
1271     if (G.debug & G_DEBUG_SIMDATA) {
1272       BKE_sim_debug_data_set_enabled(true);
1273     }
1274     else {
1275       return;
1276     }
1277   }
1278
1279   elem = MEM_callocN(sizeof(SimDebugElement), "sim debug data element");
1280   elem->type = type;
1281   elem->category_hash = category_hash;
1282   elem->hash = hash;
1283   elem->color[0] = r;
1284   elem->color[1] = g;
1285   elem->color[2] = b;
1286   if (v1) {
1287     copy_v3_v3(elem->v1, v1);
1288   }
1289   else {
1290     zero_v3(elem->v1);
1291   }
1292   if (v2) {
1293     copy_v3_v3(elem->v2, v2);
1294   }
1295   else {
1296     zero_v3(elem->v2);
1297   }
1298   if (str) {
1299     BLI_strncpy(elem->str, str, sizeof(elem->str));
1300   }
1301   else {
1302     elem->str[0] = '\0';
1303   }
1304
1305   debug_data_insert(_sim_debug_data, elem);
1306 }
1307
1308 void BKE_sim_debug_data_remove_element(uint hash)
1309 {
1310   SimDebugElement dummy;
1311   if (!_sim_debug_data) {
1312     return;
1313   }
1314   dummy.hash = hash;
1315   BLI_ghash_remove(_sim_debug_data->gh, &dummy, NULL, debug_element_free);
1316 }
1317
1318 void BKE_sim_debug_data_clear(void)
1319 {
1320   if (!_sim_debug_data) {
1321     return;
1322   }
1323   if (_sim_debug_data->gh) {
1324     BLI_ghash_clear(_sim_debug_data->gh, NULL, debug_element_free);
1325   }
1326 }
1327
1328 void BKE_sim_debug_data_clear_category(const char *category)
1329 {
1330   int category_hash = (int)BLI_ghashutil_strhash_p(category);
1331
1332   if (!_sim_debug_data) {
1333     return;
1334   }
1335
1336   if (_sim_debug_data->gh) {
1337     GHashIterator iter;
1338     BLI_ghashIterator_init(&iter, _sim_debug_data->gh);
1339     while (!BLI_ghashIterator_done(&iter)) {
1340       const SimDebugElement *elem = BLI_ghashIterator_getValue(&iter);
1341       BLI_ghashIterator_step(
1342           &iter); /* removing invalidates the current iterator, so step before removing */
1343
1344       if (elem->category_hash == category_hash) {
1345         BLI_ghash_remove(_sim_debug_data->gh, elem, NULL, debug_element_free);
1346       }
1347     }
1348   }
1349 }