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