Cleanup: style, use braces for blenkernel
[blender.git] / source / blender / blenkernel / intern / particle_system.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) 2007 by Janne Karhu.
17  * All rights reserved.
18  * Adaptive time step
19  * Classical SPH
20  * Copyright 2011-2012 AutoCRC
21  */
22
23 /** \file
24  * \ingroup bke
25  */
26
27 #include <stddef.h>
28
29 #include <stdlib.h>
30 #include <math.h>
31 #include <string.h>
32
33 #include "MEM_guardedalloc.h"
34
35 #include "DNA_anim_types.h"
36 #include "DNA_boid_types.h"
37 #include "DNA_particle_types.h"
38 #include "DNA_mesh_types.h"
39 #include "DNA_meshdata_types.h"
40 #include "DNA_modifier_types.h"
41 #include "DNA_object_force_types.h"
42 #include "DNA_object_types.h"
43 #include "DNA_curve_types.h"
44 #include "DNA_scene_types.h"
45 #include "DNA_texture_types.h"
46 #include "DNA_listBase.h"
47
48 #include "BLI_utildefines.h"
49 #include "BLI_edgehash.h"
50 #include "BLI_rand.h"
51 #include "BLI_math.h"
52 #include "BLI_blenlib.h"
53 #include "BLI_kdtree.h"
54 #include "BLI_kdopbvh.h"
55 #include "BLI_task.h"
56 #include "BLI_threads.h"
57 #include "BLI_linklist.h"
58 #include "BLI_string_utils.h"
59
60 #include "BKE_animsys.h"
61 #include "BKE_boids.h"
62 #include "BKE_collision.h"
63 #include "BKE_colortools.h"
64 #include "BKE_effect.h"
65 #include "BKE_library.h"
66 #include "BKE_library_query.h"
67 #include "BKE_particle.h"
68
69 #include "BKE_collection.h"
70 #include "BKE_object.h"
71 #include "BKE_material.h"
72 #include "BKE_cloth.h"
73 #include "BKE_lattice.h"
74 #include "BKE_pointcache.h"
75 #include "BKE_mesh.h"
76 #include "BKE_modifier.h"
77 #include "BKE_scene.h"
78 #include "BKE_bvhutils.h"
79
80 #include "DEG_depsgraph.h"
81 #include "DEG_depsgraph_physics.h"
82 #include "DEG_depsgraph_query.h"
83
84 #include "PIL_time.h"
85
86 #include "RE_shader_ext.h"
87
88 /* fluid sim particle import */
89 #ifdef WITH_MOD_FLUID
90 #  include "DNA_object_fluidsim_types.h"
91 #  include "LBM_fluidsim.h"
92 #  include <zlib.h>
93 #  include <string.h>
94
95 #endif  // WITH_MOD_FLUID
96
97 static ThreadRWMutex psys_bvhtree_rwlock = BLI_RWLOCK_INITIALIZER;
98
99 /************************************************/
100 /*          Reacting to system events           */
101 /************************************************/
102
103 static int particles_are_dynamic(ParticleSystem *psys)
104 {
105   if (psys->pointcache->flag & PTCACHE_BAKED) {
106     return 0;
107   }
108
109   if (psys->part->type == PART_HAIR) {
110     return psys->flag & PSYS_HAIR_DYNAMICS;
111   }
112   else {
113     return ELEM(psys->part->phystype, PART_PHYS_NEWTON, PART_PHYS_BOIDS, PART_PHYS_FLUID);
114   }
115 }
116
117 float psys_get_current_display_percentage(ParticleSystem *psys, const bool use_render_params)
118 {
119   ParticleSettings *part = psys->part;
120
121   if ((use_render_params &&
122        !particles_are_dynamic(psys)) ||          /* non-dynamic particles can be rendered fully */
123       (part->child_nbr && part->childtype) ||    /* display percentage applies to children */
124       (psys->pointcache->flag & PTCACHE_BAKING)) /* baking is always done with full amount */
125   {
126     return 1.0f;
127   }
128
129   return psys->part->disp / 100.0f;
130 }
131
132 static int tot_particles(ParticleSystem *psys, PTCacheID *pid)
133 {
134   if (pid && psys->pointcache->flag & PTCACHE_EXTERNAL) {
135     return pid->cache->totpoint;
136   }
137   else if (psys->part->distr == PART_DISTR_GRID && psys->part->from != PART_FROM_VERT) {
138     return psys->part->grid_res * psys->part->grid_res * psys->part->grid_res - psys->totunexist;
139   }
140   else {
141     return psys->part->totpart - psys->totunexist;
142   }
143 }
144
145 void psys_reset(ParticleSystem *psys, int mode)
146 {
147   PARTICLE_P;
148
149   if (ELEM(mode, PSYS_RESET_ALL, PSYS_RESET_DEPSGRAPH)) {
150     if (mode == PSYS_RESET_ALL || !(psys->flag & PSYS_EDITED)) {
151       /* don't free if not absolutely necessary */
152       if (psys->totpart != tot_particles(psys, NULL)) {
153         psys_free_particles(psys);
154         psys->totpart = 0;
155       }
156
157       psys->totkeyed = 0;
158       psys->flag &= ~(PSYS_HAIR_DONE | PSYS_KEYED);
159
160       if (psys->edit && psys->free_edit) {
161         psys->free_edit(psys->edit);
162         psys->edit = NULL;
163         psys->free_edit = NULL;
164       }
165     }
166   }
167   else if (mode == PSYS_RESET_CACHE_MISS) {
168     /* set all particles to be skipped */
169     LOOP_PARTICLES
170     {
171       pa->flag |= PARS_NO_DISP;
172     }
173   }
174
175   /* reset children */
176   if (psys->child) {
177     MEM_freeN(psys->child);
178     psys->child = NULL;
179   }
180
181   psys->totchild = 0;
182
183   /* reset path cache */
184   psys_free_path_cache(psys, psys->edit);
185
186   /* reset point cache */
187   BKE_ptcache_invalidate(psys->pointcache);
188
189   if (psys->fluid_springs) {
190     MEM_freeN(psys->fluid_springs);
191     psys->fluid_springs = NULL;
192   }
193
194   psys->tot_fluidsprings = psys->alloc_fluidsprings = 0;
195 }
196
197 void psys_unique_name(Object *object, ParticleSystem *psys, const char *defname)
198 {
199   BLI_uniquename(&object->particlesystem,
200                  psys,
201                  defname,
202                  '.',
203                  offsetof(ParticleSystem, name),
204                  sizeof(psys->name));
205 }
206
207 static void realloc_particles(ParticleSimulationData *sim, int new_totpart)
208 {
209   ParticleSystem *psys = sim->psys;
210   ParticleSettings *part = psys->part;
211   ParticleData *newpars = NULL;
212   BoidParticle *newboids = NULL;
213   PARTICLE_P;
214   int totpart, totsaved = 0;
215
216   if (new_totpart < 0) {
217     if ((part->distr == PART_DISTR_GRID) && (part->from != PART_FROM_VERT)) {
218       totpart = part->grid_res;
219       totpart *= totpart * totpart;
220     }
221     else {
222       totpart = part->totpart;
223     }
224   }
225   else {
226     totpart = new_totpart;
227   }
228
229   if (totpart != psys->totpart) {
230     if (psys->edit && psys->free_edit) {
231       psys->free_edit(psys->edit);
232       psys->edit = NULL;
233       psys->free_edit = NULL;
234     }
235
236     if (totpart) {
237       newpars = MEM_callocN(totpart * sizeof(ParticleData), "particles");
238       if (newpars == NULL) {
239         return;
240       }
241
242       if (psys->part->phystype == PART_PHYS_BOIDS) {
243         newboids = MEM_callocN(totpart * sizeof(BoidParticle), "boid particles");
244
245         if (newboids == NULL) {
246           /* allocation error! */
247           if (newpars) {
248             MEM_freeN(newpars);
249           }
250           return;
251         }
252       }
253     }
254
255     if (psys->particles) {
256       totsaved = MIN2(psys->totpart, totpart);
257       /*save old pars*/
258       if (totsaved) {
259         memcpy(newpars, psys->particles, totsaved * sizeof(ParticleData));
260
261         if (psys->particles->boid) {
262           memcpy(newboids, psys->particles->boid, totsaved * sizeof(BoidParticle));
263         }
264       }
265
266       if (psys->particles->keys) {
267         MEM_freeN(psys->particles->keys);
268       }
269
270       if (psys->particles->boid) {
271         MEM_freeN(psys->particles->boid);
272       }
273
274       for (p = 0, pa = newpars; p < totsaved; p++, pa++) {
275         if (pa->keys) {
276           pa->keys = NULL;
277           pa->totkey = 0;
278         }
279       }
280
281       for (p = totsaved, pa = psys->particles + totsaved; p < psys->totpart; p++, pa++) {
282         if (pa->hair) {
283           MEM_freeN(pa->hair);
284         }
285       }
286
287       MEM_freeN(psys->particles);
288       psys_free_pdd(psys);
289     }
290
291     psys->particles = newpars;
292     psys->totpart = totpart;
293
294     if (newboids) {
295       LOOP_PARTICLES
296       {
297         pa->boid = newboids++;
298       }
299     }
300   }
301
302   if (psys->child) {
303     MEM_freeN(psys->child);
304     psys->child = NULL;
305     psys->totchild = 0;
306   }
307 }
308
309 int psys_get_child_number(Scene *scene, ParticleSystem *psys, const bool use_render_params)
310 {
311   int nbr;
312
313   if (!psys->part->childtype) {
314     return 0;
315   }
316
317   if (use_render_params) {
318     nbr = psys->part->ren_child_nbr;
319   }
320   else {
321     nbr = psys->part->child_nbr;
322   }
323
324   return get_render_child_particle_number(&scene->r, nbr, use_render_params);
325 }
326
327 int psys_get_tot_child(Scene *scene, ParticleSystem *psys, const bool use_render_params)
328 {
329   return psys->totpart * psys_get_child_number(scene, psys, use_render_params);
330 }
331
332 /************************************************/
333 /*          Distribution                        */
334 /************************************************/
335
336 void psys_calc_dmcache(Object *ob, Mesh *mesh_final, Mesh *mesh_original, ParticleSystem *psys)
337 {
338   /* use for building derived mesh mapping info:
339    *
340    * node: the allocated links - total derived mesh element count
341    * nodearray: the array of nodes aligned with the base mesh's elements, so
342    *            each original elements can reference its derived elements
343    */
344   Mesh *me = (Mesh *)ob->data;
345   bool use_modifier_stack = psys->part->use_modifier_stack;
346   PARTICLE_P;
347
348   /* CACHE LOCATIONS */
349   if (!mesh_final->runtime.deformed_only) {
350     /* Will use later to speed up subsurf/evaluated mesh. */
351     LinkNode *node, *nodedmelem, **nodearray;
352     int totdmelem, totelem, i, *origindex, *origindex_poly = NULL;
353
354     if (psys->part->from == PART_FROM_VERT) {
355       totdmelem = mesh_final->totvert;
356
357       if (use_modifier_stack) {
358         totelem = totdmelem;
359         origindex = NULL;
360       }
361       else {
362         totelem = me->totvert;
363         origindex = CustomData_get_layer(&mesh_final->vdata, CD_ORIGINDEX);
364       }
365     }
366     else { /* FROM_FACE/FROM_VOLUME */
367       totdmelem = mesh_final->totface;
368
369       if (use_modifier_stack) {
370         totelem = totdmelem;
371         origindex = NULL;
372         origindex_poly = NULL;
373       }
374       else {
375         totelem = mesh_original->totface;
376         origindex = CustomData_get_layer(&mesh_final->fdata, CD_ORIGINDEX);
377
378         /* for face lookups we need the poly origindex too */
379         origindex_poly = CustomData_get_layer(&mesh_final->pdata, CD_ORIGINDEX);
380         if (origindex_poly == NULL) {
381           origindex = NULL;
382         }
383       }
384     }
385
386     nodedmelem = MEM_callocN(sizeof(LinkNode) * totdmelem, "psys node elems");
387     nodearray = MEM_callocN(sizeof(LinkNode *) * totelem, "psys node array");
388
389     for (i = 0, node = nodedmelem; i < totdmelem; i++, node++) {
390       int origindex_final;
391       node->link = POINTER_FROM_INT(i);
392
393       /* may be vertex or face origindex */
394       if (use_modifier_stack) {
395         origindex_final = i;
396       }
397       else {
398         origindex_final = origindex ? origindex[i] : ORIGINDEX_NONE;
399
400         /* if we have a poly source, do an index lookup */
401         if (origindex_poly && origindex_final != ORIGINDEX_NONE) {
402           origindex_final = origindex_poly[origindex_final];
403         }
404       }
405
406       if (origindex_final != ORIGINDEX_NONE && origindex_final < totelem) {
407         if (nodearray[origindex_final]) {
408           /* prepend */
409           node->next = nodearray[origindex_final];
410           nodearray[origindex_final] = node;
411         }
412         else {
413           nodearray[origindex_final] = node;
414         }
415       }
416     }
417
418     /* cache the verts/faces! */
419     LOOP_PARTICLES
420     {
421       if (pa->num < 0) {
422         pa->num_dmcache = DMCACHE_NOTFOUND;
423         continue;
424       }
425
426       if (use_modifier_stack) {
427         if (pa->num < totelem) {
428           pa->num_dmcache = DMCACHE_ISCHILD;
429         }
430         else {
431           pa->num_dmcache = DMCACHE_NOTFOUND;
432         }
433       }
434       else {
435         if (psys->part->from == PART_FROM_VERT) {
436           if (pa->num < totelem && nodearray[pa->num]) {
437             pa->num_dmcache = POINTER_AS_INT(nodearray[pa->num]->link);
438           }
439           else {
440             pa->num_dmcache = DMCACHE_NOTFOUND;
441           }
442         }
443         else { /* FROM_FACE/FROM_VOLUME */
444           pa->num_dmcache = psys_particle_dm_face_lookup(
445               mesh_final, mesh_original, pa->num, pa->fuv, nodearray);
446         }
447       }
448     }
449
450     MEM_freeN(nodearray);
451     MEM_freeN(nodedmelem);
452   }
453   else {
454     /* TODO PARTICLE, make the following line unnecessary, each function
455      * should know to use the num or num_dmcache, set the num_dmcache to
456      * an invalid value, just in case */
457
458     LOOP_PARTICLES
459     {
460       pa->num_dmcache = DMCACHE_NOTFOUND;
461     }
462   }
463 }
464
465 /* threaded child particle distribution and path caching */
466 void psys_thread_context_init(ParticleThreadContext *ctx, ParticleSimulationData *sim)
467 {
468   memset(ctx, 0, sizeof(ParticleThreadContext));
469   ctx->sim = *sim;
470   ctx->mesh = ctx->sim.psmd->mesh_final;
471   ctx->ma = give_current_material(sim->ob, sim->psys->part->omat);
472 }
473
474 #define MAX_PARTICLES_PER_TASK \
475   256 /* XXX arbitrary - maybe use at least number of points instead for better balancing? */
476
477 BLI_INLINE int ceil_ii(int a, int b)
478 {
479   return (a + b - 1) / b;
480 }
481
482 void psys_tasks_create(ParticleThreadContext *ctx,
483                        int startpart,
484                        int endpart,
485                        ParticleTask **r_tasks,
486                        int *r_numtasks)
487 {
488   ParticleTask *tasks;
489   int numtasks = ceil_ii((endpart - startpart), MAX_PARTICLES_PER_TASK);
490   float particles_per_task = (float)(endpart - startpart) / (float)numtasks, p, pnext;
491   int i;
492
493   tasks = MEM_callocN(sizeof(ParticleTask) * numtasks, "ParticleThread");
494   *r_numtasks = numtasks;
495   *r_tasks = tasks;
496
497   p = (float)startpart;
498   for (i = 0; i < numtasks; i++, p = pnext) {
499     pnext = p + particles_per_task;
500
501     tasks[i].ctx = ctx;
502     tasks[i].begin = (int)p;
503     tasks[i].end = min_ii((int)pnext, endpart);
504   }
505 }
506
507 void psys_tasks_free(ParticleTask *tasks, int numtasks)
508 {
509   int i;
510
511   /* threads */
512   for (i = 0; i < numtasks; ++i) {
513     if (tasks[i].rng) {
514       BLI_rng_free(tasks[i].rng);
515     }
516     if (tasks[i].rng_path) {
517       BLI_rng_free(tasks[i].rng_path);
518     }
519   }
520
521   MEM_freeN(tasks);
522 }
523
524 void psys_thread_context_free(ParticleThreadContext *ctx)
525 {
526   /* path caching */
527   if (ctx->vg_length) {
528     MEM_freeN(ctx->vg_length);
529   }
530   if (ctx->vg_clump) {
531     MEM_freeN(ctx->vg_clump);
532   }
533   if (ctx->vg_kink) {
534     MEM_freeN(ctx->vg_kink);
535   }
536   if (ctx->vg_rough1) {
537     MEM_freeN(ctx->vg_rough1);
538   }
539   if (ctx->vg_rough2) {
540     MEM_freeN(ctx->vg_rough2);
541   }
542   if (ctx->vg_roughe) {
543     MEM_freeN(ctx->vg_roughe);
544   }
545   if (ctx->vg_twist) {
546     MEM_freeN(ctx->vg_twist);
547   }
548
549   if (ctx->sim.psys->lattice_deform_data) {
550     end_latt_deform(ctx->sim.psys->lattice_deform_data);
551     ctx->sim.psys->lattice_deform_data = NULL;
552   }
553
554   /* distribution */
555   if (ctx->jit) {
556     MEM_freeN(ctx->jit);
557   }
558   if (ctx->jitoff) {
559     MEM_freeN(ctx->jitoff);
560   }
561   if (ctx->weight) {
562     MEM_freeN(ctx->weight);
563   }
564   if (ctx->index) {
565     MEM_freeN(ctx->index);
566   }
567   if (ctx->seams) {
568     MEM_freeN(ctx->seams);
569   }
570   //if (ctx->vertpart) MEM_freeN(ctx->vertpart);
571   BLI_kdtree_3d_free(ctx->tree);
572
573   if (ctx->clumpcurve != NULL) {
574     curvemapping_free(ctx->clumpcurve);
575   }
576   if (ctx->roughcurve != NULL) {
577     curvemapping_free(ctx->roughcurve);
578   }
579   if (ctx->twistcurve != NULL) {
580     curvemapping_free(ctx->twistcurve);
581   }
582 }
583
584 static void initialize_particle_texture(ParticleSimulationData *sim, ParticleData *pa, int p)
585 {
586   ParticleSystem *psys = sim->psys;
587   ParticleSettings *part = psys->part;
588   ParticleTexture ptex;
589
590   psys_get_texture(sim, pa, &ptex, PAMAP_INIT, 0.f);
591
592   switch (part->type) {
593     case PART_EMITTER:
594       if (ptex.exist < psys_frand(psys, p + 125)) {
595         pa->flag |= PARS_UNEXIST;
596       }
597       pa->time = part->sta + (part->end - part->sta) * ptex.time;
598       break;
599     case PART_HAIR:
600       if (ptex.exist < psys_frand(psys, p + 125)) {
601         pa->flag |= PARS_UNEXIST;
602       }
603       pa->time = 0.f;
604       break;
605     case PART_FLUID:
606       break;
607   }
608 }
609
610 /* set particle parameters that don't change during particle's life */
611 void initialize_particle(ParticleSimulationData *sim, ParticleData *pa)
612 {
613   ParticleSettings *part = sim->psys->part;
614   float birth_time = (float)(pa - sim->psys->particles) / (float)sim->psys->totpart;
615
616   pa->flag &= ~PARS_UNEXIST;
617   pa->time = part->sta + (part->end - part->sta) * birth_time;
618
619   pa->hair_index = 0;
620   /* we can't reset to -1 anymore since we've figured out correct index in distribute_particles */
621   /* usage other than straight after distribute has to handle this index by itself - jahka*/
622   //pa->num_dmcache = DMCACHE_NOTFOUND; /* assume we don't have a derived mesh face */
623 }
624
625 static void initialize_all_particles(ParticleSimulationData *sim)
626 {
627   ParticleSystem *psys = sim->psys;
628   ParticleSettings *part = psys->part;
629   /* Grid distributionsets UNEXIST flag, need to take care of
630    * it here because later this flag is being reset.
631    *
632    * We can't do it for any distribution, because it'll then
633    * conflict with texture influence, which does not free
634    * unexisting particles and only sets flag.
635    *
636    * It's not so bad, because only grid distribution sets
637    * UNEXIST flag.
638    */
639   const bool emit_from_volume_grid = (part->distr == PART_DISTR_GRID) &&
640                                      (!ELEM(part->from, PART_FROM_VERT, PART_FROM_CHILD));
641   PARTICLE_P;
642   LOOP_PARTICLES
643   {
644     if (!(emit_from_volume_grid && (pa->flag & PARS_UNEXIST) != 0)) {
645       initialize_particle(sim, pa);
646     }
647   }
648 }
649
650 static void free_unexisting_particles(ParticleSimulationData *sim)
651 {
652   ParticleSystem *psys = sim->psys;
653   PARTICLE_P;
654
655   psys->totunexist = 0;
656
657   LOOP_PARTICLES
658   {
659     if (pa->flag & PARS_UNEXIST) {
660       psys->totunexist++;
661     }
662   }
663
664   if (psys->totpart && psys->totunexist == psys->totpart) {
665     if (psys->particles->boid) {
666       MEM_freeN(psys->particles->boid);
667     }
668
669     MEM_freeN(psys->particles);
670     psys->particles = NULL;
671     psys->totpart = psys->totunexist = 0;
672   }
673
674   if (psys->totunexist) {
675     int newtotpart = psys->totpart - psys->totunexist;
676     ParticleData *npa, *newpars;
677
678     npa = newpars = MEM_callocN(newtotpart * sizeof(ParticleData), "particles");
679
680     for (p = 0, pa = psys->particles; p < newtotpart; p++, pa++, npa++) {
681       while (pa->flag & PARS_UNEXIST) {
682         pa++;
683       }
684
685       memcpy(npa, pa, sizeof(ParticleData));
686     }
687
688     if (psys->particles->boid) {
689       MEM_freeN(psys->particles->boid);
690     }
691     MEM_freeN(psys->particles);
692     psys->particles = newpars;
693     psys->totpart -= psys->totunexist;
694
695     if (psys->particles->boid) {
696       BoidParticle *newboids = MEM_callocN(psys->totpart * sizeof(BoidParticle), "boid particles");
697
698       LOOP_PARTICLES
699       {
700         pa->boid = newboids++;
701       }
702     }
703   }
704 }
705
706 static void get_angular_velocity_vector(short avemode, ParticleKey *state, float vec[3])
707 {
708   switch (avemode) {
709     case PART_AVE_VELOCITY:
710       copy_v3_v3(vec, state->vel);
711       break;
712     case PART_AVE_HORIZONTAL: {
713       float zvec[3];
714       zvec[0] = zvec[1] = 0;
715       zvec[2] = 1.f;
716       cross_v3_v3v3(vec, state->vel, zvec);
717       break;
718     }
719     case PART_AVE_VERTICAL: {
720       float zvec[3], temp[3];
721       zvec[0] = zvec[1] = 0;
722       zvec[2] = 1.f;
723       cross_v3_v3v3(temp, state->vel, zvec);
724       cross_v3_v3v3(vec, temp, state->vel);
725       break;
726     }
727     case PART_AVE_GLOBAL_X:
728       vec[0] = 1.f;
729       vec[1] = vec[2] = 0;
730       break;
731     case PART_AVE_GLOBAL_Y:
732       vec[1] = 1.f;
733       vec[0] = vec[2] = 0;
734       break;
735     case PART_AVE_GLOBAL_Z:
736       vec[2] = 1.f;
737       vec[0] = vec[1] = 0;
738       break;
739   }
740 }
741
742 void psys_get_birth_coords(
743     ParticleSimulationData *sim, ParticleData *pa, ParticleKey *state, float dtime, float cfra)
744 {
745   Object *ob = sim->ob;
746   ParticleSystem *psys = sim->psys;
747   ParticleSettings *part = psys->part;
748   ParticleTexture ptex;
749   float fac, phasefac, nor[3] = {0, 0, 0}, loc[3], vel[3] = {0.0, 0.0, 0.0}, rot[4], q2[4];
750   float r_vel[3], r_ave[3], r_rot[4], vec[3], p_vel[3] = {0.0, 0.0, 0.0};
751   float x_vec[3] = {1.0, 0.0, 0.0}, utan[3] = {0.0, 1.0, 0.0}, vtan[3] = {0.0, 0.0, 1.0},
752         rot_vec[3] = {0.0, 0.0, 0.0};
753   float q_phase[4];
754
755   const bool use_boids = ((part->phystype == PART_PHYS_BOIDS) && (pa->boid != NULL));
756   const bool use_tangents = ((use_boids == false) &&
757                              ((part->tanfac != 0.0f) || (part->rotmode == PART_ROT_NOR_TAN)));
758
759   int p = pa - psys->particles;
760
761   /* get birth location from object       */
762   if (use_tangents) {
763     psys_particle_on_emitter(sim->psmd,
764                              part->from,
765                              pa->num,
766                              pa->num_dmcache,
767                              pa->fuv,
768                              pa->foffset,
769                              loc,
770                              nor,
771                              utan,
772                              vtan,
773                              0);
774   }
775   else {
776     psys_particle_on_emitter(
777         sim->psmd, part->from, pa->num, pa->num_dmcache, pa->fuv, pa->foffset, loc, nor, 0, 0, 0);
778   }
779
780   /* get possible textural influence */
781   psys_get_texture(sim, pa, &ptex, PAMAP_IVEL, cfra);
782
783   /* particles live in global space so    */
784   /* let's convert:                       */
785   /* -location                            */
786   mul_m4_v3(ob->obmat, loc);
787
788   /* -normal                              */
789   mul_mat3_m4_v3(ob->obmat, nor);
790   normalize_v3(nor);
791
792   /* -tangent                             */
793   if (use_tangents) {
794     //float phase=vg_rot?2.0f*(psys_particle_value_from_verts(sim->psmd->dm,part->from,pa,vg_rot)-0.5f):0.0f;
795     float phase = 0.0f;
796     mul_v3_fl(vtan, -cosf((float)M_PI * (part->tanphase + phase)));
797     fac = -sinf((float)M_PI * (part->tanphase + phase));
798     madd_v3_v3fl(vtan, utan, fac);
799
800     mul_mat3_m4_v3(ob->obmat, vtan);
801
802     copy_v3_v3(utan, nor);
803     mul_v3_fl(utan, dot_v3v3(vtan, nor));
804     sub_v3_v3(vtan, utan);
805
806     normalize_v3(vtan);
807   }
808
809   /* -velocity (boids need this even if there's no random velocity) */
810   if (part->randfac != 0.0f || (part->phystype == PART_PHYS_BOIDS && pa->boid)) {
811     r_vel[0] = 2.0f * (psys_frand(psys, p + 10) - 0.5f);
812     r_vel[1] = 2.0f * (psys_frand(psys, p + 11) - 0.5f);
813     r_vel[2] = 2.0f * (psys_frand(psys, p + 12) - 0.5f);
814
815     mul_mat3_m4_v3(ob->obmat, r_vel);
816     normalize_v3(r_vel);
817   }
818
819   /* -angular velocity                    */
820   if (part->avemode == PART_AVE_RAND) {
821     r_ave[0] = 2.0f * (psys_frand(psys, p + 13) - 0.5f);
822     r_ave[1] = 2.0f * (psys_frand(psys, p + 14) - 0.5f);
823     r_ave[2] = 2.0f * (psys_frand(psys, p + 15) - 0.5f);
824
825     mul_mat3_m4_v3(ob->obmat, r_ave);
826     normalize_v3(r_ave);
827   }
828
829   /* -rotation                            */
830   if (part->randrotfac != 0.0f) {
831     r_rot[0] = 2.0f * (psys_frand(psys, p + 16) - 0.5f);
832     r_rot[1] = 2.0f * (psys_frand(psys, p + 17) - 0.5f);
833     r_rot[2] = 2.0f * (psys_frand(psys, p + 18) - 0.5f);
834     r_rot[3] = 2.0f * (psys_frand(psys, p + 19) - 0.5f);
835     normalize_qt(r_rot);
836
837     mat4_to_quat(rot, ob->obmat);
838     mul_qt_qtqt(r_rot, r_rot, rot);
839   }
840
841   if (use_boids) {
842     float dvec[3], q[4], mat[3][3];
843
844     copy_v3_v3(state->co, loc);
845
846     /* boids don't get any initial velocity  */
847     zero_v3(state->vel);
848
849     /* boids store direction in ave */
850     if (fabsf(nor[2]) == 1.0f) {
851       sub_v3_v3v3(state->ave, loc, ob->obmat[3]);
852       normalize_v3(state->ave);
853     }
854     else {
855       copy_v3_v3(state->ave, nor);
856     }
857
858     /* calculate rotation matrix */
859     project_v3_v3v3(dvec, r_vel, state->ave);
860     sub_v3_v3v3(mat[0], state->ave, dvec);
861     normalize_v3(mat[0]);
862     negate_v3_v3(mat[2], r_vel);
863     normalize_v3(mat[2]);
864     cross_v3_v3v3(mat[1], mat[2], mat[0]);
865
866     /* apply rotation */
867     mat3_to_quat_is_ok(q, mat);
868     copy_qt_qt(state->rot, q);
869   }
870   else {
871     /* conversion done so now we apply new: */
872     /* -velocity from:                      */
873
874     /*      *reactions                      */
875     if (dtime > 0.f) {
876       sub_v3_v3v3(vel, pa->state.vel, pa->prev_state.vel);
877     }
878
879     /*      *emitter velocity               */
880     if (dtime != 0.f && part->obfac != 0.f) {
881       sub_v3_v3v3(vel, loc, state->co);
882       mul_v3_fl(vel, part->obfac / dtime);
883     }
884
885     /*      *emitter normal                 */
886     if (part->normfac != 0.f) {
887       madd_v3_v3fl(vel, nor, part->normfac);
888     }
889
890     /*      *emitter tangent                */
891     if (sim->psmd && part->tanfac != 0.f) {
892       madd_v3_v3fl(vel, vtan, part->tanfac);
893     }
894
895     /*      *emitter object orientation     */
896     if (part->ob_vel[0] != 0.f) {
897       normalize_v3_v3(vec, ob->obmat[0]);
898       madd_v3_v3fl(vel, vec, part->ob_vel[0]);
899     }
900     if (part->ob_vel[1] != 0.f) {
901       normalize_v3_v3(vec, ob->obmat[1]);
902       madd_v3_v3fl(vel, vec, part->ob_vel[1]);
903     }
904     if (part->ob_vel[2] != 0.f) {
905       normalize_v3_v3(vec, ob->obmat[2]);
906       madd_v3_v3fl(vel, vec, part->ob_vel[2]);
907     }
908
909     /*      *texture                        */
910     /* TODO */
911
912     /*      *random                         */
913     if (part->randfac != 0.f) {
914       madd_v3_v3fl(vel, r_vel, part->randfac);
915     }
916
917     /*      *particle                       */
918     if (part->partfac != 0.f) {
919       madd_v3_v3fl(vel, p_vel, part->partfac);
920     }
921
922     mul_v3_v3fl(state->vel, vel, ptex.ivel);
923
924     /* -location from emitter               */
925     copy_v3_v3(state->co, loc);
926
927     /* -rotation                            */
928     unit_qt(state->rot);
929
930     if (part->rotmode) {
931       bool use_global_space;
932
933       /* create vector into which rotation is aligned */
934       switch (part->rotmode) {
935         case PART_ROT_NOR:
936         case PART_ROT_NOR_TAN:
937           copy_v3_v3(rot_vec, nor);
938           use_global_space = false;
939           break;
940         case PART_ROT_VEL:
941           copy_v3_v3(rot_vec, vel);
942           use_global_space = true;
943           break;
944         case PART_ROT_GLOB_X:
945         case PART_ROT_GLOB_Y:
946         case PART_ROT_GLOB_Z:
947           rot_vec[part->rotmode - PART_ROT_GLOB_X] = 1.0f;
948           use_global_space = true;
949           break;
950         case PART_ROT_OB_X:
951         case PART_ROT_OB_Y:
952         case PART_ROT_OB_Z:
953           copy_v3_v3(rot_vec, ob->obmat[part->rotmode - PART_ROT_OB_X]);
954           use_global_space = false;
955           break;
956         default:
957           use_global_space = true;
958           break;
959       }
960
961       /* create rotation quat */
962
963       if (use_global_space) {
964         negate_v3(rot_vec);
965         vec_to_quat(q2, rot_vec, OB_POSX, OB_POSZ);
966
967         /* randomize rotation quat */
968         if (part->randrotfac != 0.0f) {
969           interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
970         }
971         else {
972           copy_qt_qt(rot, q2);
973         }
974       }
975       else {
976         /* calculate rotation in local-space */
977         float q_obmat[4];
978         float q_imat[4];
979
980         mat4_to_quat(q_obmat, ob->obmat);
981         invert_qt_qt_normalized(q_imat, q_obmat);
982
983         if (part->rotmode != PART_ROT_NOR_TAN) {
984           float rot_vec_local[3];
985
986           /* rot_vec */
987           negate_v3(rot_vec);
988           copy_v3_v3(rot_vec_local, rot_vec);
989           mul_qt_v3(q_imat, rot_vec_local);
990           normalize_v3(rot_vec_local);
991
992           vec_to_quat(q2, rot_vec_local, OB_POSX, OB_POSZ);
993         }
994         else {
995           /* (part->rotmode == PART_ROT_NOR_TAN) */
996           float tmat[3][3];
997
998           /* note: utan_local is not taken from 'utan', we calculate from rot_vec/vtan */
999           /* note: it looks like rotation phase may be applied twice (once with vtan, again below)
1000            * however this isn't the case - campbell */
1001           float *rot_vec_local = tmat[0];
1002           float *vtan_local = tmat[1];
1003           float *utan_local = tmat[2];
1004
1005           /* use tangents */
1006           BLI_assert(use_tangents == true);
1007
1008           /* rot_vec */
1009           copy_v3_v3(rot_vec_local, rot_vec);
1010           mul_qt_v3(q_imat, rot_vec_local);
1011
1012           /* vtan_local */
1013           copy_v3_v3(vtan_local, vtan); /* flips, cant use */
1014           mul_qt_v3(q_imat, vtan_local);
1015
1016           /* ensure orthogonal matrix (rot_vec aligned) */
1017           cross_v3_v3v3(utan_local, vtan_local, rot_vec_local);
1018           cross_v3_v3v3(vtan_local, utan_local, rot_vec_local);
1019
1020           /* note: no need to normalize */
1021           mat3_to_quat(q2, tmat);
1022         }
1023
1024         /* randomize rotation quat */
1025         if (part->randrotfac != 0.0f) {
1026           mul_qt_qtqt(r_rot, r_rot, q_imat);
1027           interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
1028         }
1029         else {
1030           copy_qt_qt(rot, q2);
1031         }
1032
1033         mul_qt_qtqt(rot, q_obmat, rot);
1034       }
1035
1036       /* rotation phase */
1037       phasefac = part->phasefac;
1038       if (part->randphasefac != 0.0f) {
1039         phasefac += part->randphasefac * psys_frand(psys, p + 20);
1040       }
1041       axis_angle_to_quat(q_phase, x_vec, phasefac * (float)M_PI);
1042
1043       /* combine base rotation & phase */
1044       mul_qt_qtqt(state->rot, rot, q_phase);
1045     }
1046
1047     /* -angular velocity                    */
1048
1049     zero_v3(state->ave);
1050
1051     if (part->avemode) {
1052       if (part->avemode == PART_AVE_RAND) {
1053         copy_v3_v3(state->ave, r_ave);
1054       }
1055       else {
1056         get_angular_velocity_vector(part->avemode, state, state->ave);
1057       }
1058
1059       normalize_v3(state->ave);
1060       mul_v3_fl(state->ave, part->avefac);
1061     }
1062   }
1063 }
1064
1065 /* recursively evaluate emitter parent anim at cfra */
1066 static void evaluate_emitter_anim(struct Depsgraph *depsgraph,
1067                                   Scene *scene,
1068                                   Object *ob,
1069                                   float cfra)
1070 {
1071   if (ob->parent) {
1072     evaluate_emitter_anim(depsgraph, scene, ob->parent, cfra);
1073   }
1074
1075   BKE_object_where_is_calc_time(depsgraph, scene, ob, cfra);
1076 }
1077
1078 /* sets particle to the emitter surface with initial velocity & rotation */
1079 void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, float cfra)
1080 {
1081   ParticleSystem *psys = sim->psys;
1082   ParticleSettings *part;
1083   ParticleTexture ptex;
1084   int p = pa - psys->particles;
1085   part = psys->part;
1086
1087   /* get precise emitter matrix if particle is born */
1088   if (part->type != PART_HAIR && dtime > 0.f && pa->time < cfra && pa->time >= sim->psys->cfra) {
1089     evaluate_emitter_anim(sim->depsgraph, sim->scene, sim->ob, pa->time);
1090
1091     psys->flag |= PSYS_OB_ANIM_RESTORE;
1092   }
1093
1094   psys_get_birth_coords(sim, pa, &pa->state, dtime, cfra);
1095
1096   /* Initialize particle settings which depends on texture.
1097    *
1098    * We could only do it now because we'll need to know coordinate
1099    * before sampling the texture.
1100    */
1101   initialize_particle_texture(sim, pa, p);
1102
1103   if (part->phystype == PART_PHYS_BOIDS && pa->boid) {
1104     BoidParticle *bpa = pa->boid;
1105
1106     /* and gravity in r_ve */
1107     bpa->gravity[0] = bpa->gravity[1] = 0.0f;
1108     bpa->gravity[2] = -1.0f;
1109     if ((sim->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) &&
1110         (sim->scene->physics_settings.gravity[2] != 0.0f)) {
1111       bpa->gravity[2] = sim->scene->physics_settings.gravity[2];
1112     }
1113
1114     bpa->data.health = part->boids->health;
1115     bpa->data.mode = eBoidMode_InAir;
1116     bpa->data.state_id = ((BoidState *)part->boids->states.first)->id;
1117     bpa->data.acc[0] = bpa->data.acc[1] = bpa->data.acc[2] = 0.0f;
1118   }
1119
1120   if (part->type == PART_HAIR) {
1121     pa->lifetime = 100.0f;
1122   }
1123   else {
1124     /* initialize the lifetime, in case the texture coordinates
1125      * are from Particles/Strands, which would cause undefined values
1126      */
1127     pa->lifetime = part->lifetime * (1.0f - part->randlife * psys_frand(psys, p + 21));
1128     pa->dietime = pa->time + pa->lifetime;
1129
1130     /* get possible textural influence */
1131     psys_get_texture(sim, pa, &ptex, PAMAP_LIFE, cfra);
1132
1133     pa->lifetime = part->lifetime * ptex.life;
1134
1135     if (part->randlife != 0.0f) {
1136       pa->lifetime *= 1.0f - part->randlife * psys_frand(psys, p + 21);
1137     }
1138   }
1139
1140   pa->dietime = pa->time + pa->lifetime;
1141
1142   if ((sim->psys->pointcache) && (sim->psys->pointcache->flag & PTCACHE_BAKED) &&
1143       (sim->psys->pointcache->mem_cache.first)) {
1144     float dietime = psys_get_dietime_from_cache(sim->psys->pointcache, p);
1145     pa->dietime = MIN2(pa->dietime, dietime);
1146   }
1147
1148   if (pa->time > cfra) {
1149     pa->alive = PARS_UNBORN;
1150   }
1151   else if (pa->dietime <= cfra) {
1152     pa->alive = PARS_DEAD;
1153   }
1154   else {
1155     pa->alive = PARS_ALIVE;
1156   }
1157
1158   pa->state.time = cfra;
1159 }
1160 static void reset_all_particles(ParticleSimulationData *sim, float dtime, float cfra, int from)
1161 {
1162   ParticleData *pa;
1163   int p, totpart = sim->psys->totpart;
1164
1165   for (p = from, pa = sim->psys->particles + from; p < totpart; p++, pa++) {
1166     reset_particle(sim, pa, dtime, cfra);
1167   }
1168 }
1169 /************************************************/
1170 /*          Particle targets                    */
1171 /************************************************/
1172 ParticleSystem *psys_get_target_system(Object *ob, ParticleTarget *pt)
1173 {
1174   ParticleSystem *psys = NULL;
1175
1176   if (pt->ob == NULL || pt->ob == ob) {
1177     psys = BLI_findlink(&ob->particlesystem, pt->psys - 1);
1178   }
1179   else {
1180     psys = BLI_findlink(&pt->ob->particlesystem, pt->psys - 1);
1181   }
1182
1183   if (psys) {
1184     pt->flag |= PTARGET_VALID;
1185   }
1186   else {
1187     pt->flag &= ~PTARGET_VALID;
1188   }
1189
1190   return psys;
1191 }
1192 /************************************************/
1193 /*          Keyed particles                     */
1194 /************************************************/
1195 /* Counts valid keyed targets */
1196 void psys_count_keyed_targets(ParticleSimulationData *sim)
1197 {
1198   ParticleSystem *psys = sim->psys, *kpsys;
1199   ParticleTarget *pt = psys->targets.first;
1200   int keys_valid = 1;
1201   psys->totkeyed = 0;
1202
1203   for (; pt; pt = pt->next) {
1204     kpsys = psys_get_target_system(sim->ob, pt);
1205
1206     if (kpsys && kpsys->totpart) {
1207       psys->totkeyed += keys_valid;
1208       if (psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f) {
1209         psys->totkeyed += 1;
1210       }
1211     }
1212     else {
1213       keys_valid = 0;
1214     }
1215   }
1216
1217   psys->totkeyed *= psys->flag & PSYS_KEYED_TIMING ? 1 : psys->part->keyed_loops;
1218 }
1219
1220 static void set_keyed_keys(ParticleSimulationData *sim)
1221 {
1222   ParticleSystem *psys = sim->psys;
1223   ParticleSimulationData ksim = {0};
1224   ParticleTarget *pt;
1225   PARTICLE_P;
1226   ParticleKey *key;
1227   int totpart = psys->totpart, k, totkeys = psys->totkeyed;
1228   int keyed_flag = 0;
1229
1230   ksim.depsgraph = sim->depsgraph;
1231   ksim.scene = sim->scene;
1232
1233   /* no proper targets so let's clear and bail out */
1234   if (psys->totkeyed == 0) {
1235     free_keyed_keys(psys);
1236     psys->flag &= ~PSYS_KEYED;
1237     return;
1238   }
1239
1240   if (totpart && psys->particles->totkey != totkeys) {
1241     free_keyed_keys(psys);
1242
1243     key = MEM_callocN(totpart * totkeys * sizeof(ParticleKey), "Keyed keys");
1244
1245     LOOP_PARTICLES
1246     {
1247       pa->keys = key;
1248       pa->totkey = totkeys;
1249       key += totkeys;
1250     }
1251   }
1252
1253   psys->flag &= ~PSYS_KEYED;
1254
1255   pt = psys->targets.first;
1256   for (k = 0; k < totkeys; k++) {
1257     ksim.ob = pt->ob ? pt->ob : sim->ob;
1258     ksim.psys = BLI_findlink(&ksim.ob->particlesystem, pt->psys - 1);
1259     keyed_flag = (ksim.psys->flag & PSYS_KEYED);
1260     ksim.psys->flag &= ~PSYS_KEYED;
1261
1262     LOOP_PARTICLES
1263     {
1264       key = pa->keys + k;
1265       key->time = -1.0; /* use current time */
1266
1267       psys_get_particle_state(&ksim, p % ksim.psys->totpart, key, 1);
1268
1269       if (psys->flag & PSYS_KEYED_TIMING) {
1270         key->time = pa->time + pt->time;
1271         if (pt->duration != 0.0f && k + 1 < totkeys) {
1272           copy_particle_key(key + 1, key, 1);
1273           (key + 1)->time = pa->time + pt->time + pt->duration;
1274         }
1275       }
1276       else if (totkeys > 1) {
1277         key->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
1278       }
1279       else {
1280         key->time = pa->time;
1281       }
1282     }
1283
1284     if (psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f) {
1285       k++;
1286     }
1287
1288     ksim.psys->flag |= keyed_flag;
1289
1290     pt = (pt->next && pt->next->flag & PTARGET_VALID) ? pt->next : psys->targets.first;
1291   }
1292
1293   psys->flag |= PSYS_KEYED;
1294 }
1295
1296 /************************************************/
1297 /*          Point Cache                         */
1298 /************************************************/
1299 void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
1300 {
1301   PointCache *cache = psys->pointcache;
1302
1303   if (cache->flag & PTCACHE_DISK_CACHE && BLI_listbase_is_empty(&cache->mem_cache)) {
1304     PTCacheID pid;
1305     BKE_ptcache_id_from_particles(&pid, ob, psys);
1306     cache->flag &= ~PTCACHE_DISK_CACHE;
1307     BKE_ptcache_disk_to_mem(&pid);
1308     cache->flag |= PTCACHE_DISK_CACHE;
1309   }
1310 }
1311 static void psys_clear_temp_pointcache(ParticleSystem *psys)
1312 {
1313   if (psys->pointcache->flag & PTCACHE_DISK_CACHE) {
1314     BKE_ptcache_free_mem(&psys->pointcache->mem_cache);
1315   }
1316 }
1317 void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra, int *efra)
1318 {
1319   ParticleSettings *part = psys->part;
1320
1321   *sfra = max_ii(1, (int)part->sta);
1322   *efra = min_ii((int)(part->end + part->lifetime + 1.0f), max_ii(scene->r.pefra, scene->r.efra));
1323 }
1324
1325 /************************************************/
1326 /*          Effectors                           */
1327 /************************************************/
1328 static void psys_update_particle_bvhtree(ParticleSystem *psys, float cfra)
1329 {
1330   if (psys) {
1331     PARTICLE_P;
1332     int totpart = 0;
1333     bool need_rebuild;
1334
1335     BLI_rw_mutex_lock(&psys_bvhtree_rwlock, THREAD_LOCK_READ);
1336     need_rebuild = !psys->bvhtree || psys->bvhtree_frame != cfra;
1337     BLI_rw_mutex_unlock(&psys_bvhtree_rwlock);
1338
1339     if (need_rebuild) {
1340       LOOP_SHOWN_PARTICLES
1341       {
1342         totpart++;
1343       }
1344
1345       BLI_rw_mutex_lock(&psys_bvhtree_rwlock, THREAD_LOCK_WRITE);
1346
1347       BLI_bvhtree_free(psys->bvhtree);
1348       psys->bvhtree = BLI_bvhtree_new(totpart, 0.0, 4, 6);
1349
1350       LOOP_SHOWN_PARTICLES
1351       {
1352         if (pa->alive == PARS_ALIVE) {
1353           if (pa->state.time == cfra) {
1354             BLI_bvhtree_insert(psys->bvhtree, p, pa->prev_state.co, 1);
1355           }
1356           else {
1357             BLI_bvhtree_insert(psys->bvhtree, p, pa->state.co, 1);
1358           }
1359         }
1360       }
1361       BLI_bvhtree_balance(psys->bvhtree);
1362
1363       psys->bvhtree_frame = cfra;
1364
1365       BLI_rw_mutex_unlock(&psys_bvhtree_rwlock);
1366     }
1367   }
1368 }
1369 void psys_update_particle_tree(ParticleSystem *psys, float cfra)
1370 {
1371   if (psys) {
1372     PARTICLE_P;
1373     int totpart = 0;
1374
1375     if (!psys->tree || psys->tree_frame != cfra) {
1376       LOOP_SHOWN_PARTICLES
1377       {
1378         totpart++;
1379       }
1380
1381       BLI_kdtree_3d_free(psys->tree);
1382       psys->tree = BLI_kdtree_3d_new(psys->totpart);
1383
1384       LOOP_SHOWN_PARTICLES
1385       {
1386         if (pa->alive == PARS_ALIVE) {
1387           if (pa->state.time == cfra) {
1388             BLI_kdtree_3d_insert(psys->tree, p, pa->prev_state.co);
1389           }
1390           else {
1391             BLI_kdtree_3d_insert(psys->tree, p, pa->state.co);
1392           }
1393         }
1394       }
1395       BLI_kdtree_3d_balance(psys->tree);
1396
1397       psys->tree_frame = cfra;
1398     }
1399   }
1400 }
1401
1402 static void psys_update_effectors(ParticleSimulationData *sim)
1403 {
1404   BKE_effectors_free(sim->psys->effectors);
1405   sim->psys->effectors = BKE_effectors_create(
1406       sim->depsgraph, sim->ob, sim->psys, sim->psys->part->effector_weights);
1407   precalc_guides(sim, sim->psys->effectors);
1408 }
1409
1410 static void integrate_particle(
1411     ParticleSettings *part,
1412     ParticleData *pa,
1413     float dtime,
1414     float *external_acceleration,
1415     void (*force_func)(void *forcedata, ParticleKey *state, float *force, float *impulse),
1416     void *forcedata)
1417 {
1418 #define ZERO_F43 \
1419   { \
1420     {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, \
1421     { \
1422       0.0f, 0.0f, 0.0f \
1423     } \
1424   }
1425
1426   ParticleKey states[5];
1427   float force[3], acceleration[3], impulse[3], dx[4][3] = ZERO_F43, dv[4][3] = ZERO_F43, oldpos[3];
1428   float pa_mass = (part->flag & PART_SIZEMASS ? part->mass * pa->size : part->mass);
1429   int i, steps = 1;
1430   int integrator = part->integrator;
1431
1432 #undef ZERO_F43
1433
1434   copy_v3_v3(oldpos, pa->state.co);
1435
1436   /* Verlet integration behaves strangely with moving emitters, so do first step with euler. */
1437   if (pa->prev_state.time < 0.f && integrator == PART_INT_VERLET) {
1438     integrator = PART_INT_EULER;
1439   }
1440
1441   switch (integrator) {
1442     case PART_INT_EULER:
1443       steps = 1;
1444       break;
1445     case PART_INT_MIDPOINT:
1446       steps = 2;
1447       break;
1448     case PART_INT_RK4:
1449       steps = 4;
1450       break;
1451     case PART_INT_VERLET:
1452       steps = 1;
1453       break;
1454   }
1455
1456   for (i = 0; i < steps; i++) {
1457     copy_particle_key(states + i, &pa->state, 1);
1458   }
1459
1460   states->time = 0.f;
1461
1462   for (i = 0; i < steps; i++) {
1463     zero_v3(force);
1464     zero_v3(impulse);
1465
1466     force_func(forcedata, states + i, force, impulse);
1467
1468     /* force to acceleration*/
1469     mul_v3_v3fl(acceleration, force, 1.0f / pa_mass);
1470
1471     if (external_acceleration) {
1472       add_v3_v3(acceleration, external_acceleration);
1473     }
1474
1475     /* calculate next state */
1476     add_v3_v3(states[i].vel, impulse);
1477
1478     switch (integrator) {
1479       case PART_INT_EULER:
1480         madd_v3_v3v3fl(pa->state.co, states->co, states->vel, dtime);
1481         madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
1482         break;
1483       case PART_INT_MIDPOINT:
1484         if (i == 0) {
1485           madd_v3_v3v3fl(states[1].co, states->co, states->vel, dtime * 0.5f);
1486           madd_v3_v3v3fl(states[1].vel, states->vel, acceleration, dtime * 0.5f);
1487           states[1].time = dtime * 0.5f;
1488           /*fra=sim->psys->cfra+0.5f*dfra;*/
1489         }
1490         else {
1491           madd_v3_v3v3fl(pa->state.co, states->co, states[1].vel, dtime);
1492           madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
1493         }
1494         break;
1495       case PART_INT_RK4:
1496         switch (i) {
1497           case 0:
1498             copy_v3_v3(dx[0], states->vel);
1499             mul_v3_fl(dx[0], dtime);
1500             copy_v3_v3(dv[0], acceleration);
1501             mul_v3_fl(dv[0], dtime);
1502
1503             madd_v3_v3v3fl(states[1].co, states->co, dx[0], 0.5f);
1504             madd_v3_v3v3fl(states[1].vel, states->vel, dv[0], 0.5f);
1505             states[1].time = dtime * 0.5f;
1506             /*fra=sim->psys->cfra+0.5f*dfra;*/
1507             break;
1508           case 1:
1509             madd_v3_v3v3fl(dx[1], states->vel, dv[0], 0.5f);
1510             mul_v3_fl(dx[1], dtime);
1511             copy_v3_v3(dv[1], acceleration);
1512             mul_v3_fl(dv[1], dtime);
1513
1514             madd_v3_v3v3fl(states[2].co, states->co, dx[1], 0.5f);
1515             madd_v3_v3v3fl(states[2].vel, states->vel, dv[1], 0.5f);
1516             states[2].time = dtime * 0.5f;
1517             break;
1518           case 2:
1519             madd_v3_v3v3fl(dx[2], states->vel, dv[1], 0.5f);
1520             mul_v3_fl(dx[2], dtime);
1521             copy_v3_v3(dv[2], acceleration);
1522             mul_v3_fl(dv[2], dtime);
1523
1524             add_v3_v3v3(states[3].co, states->co, dx[2]);
1525             add_v3_v3v3(states[3].vel, states->vel, dv[2]);
1526             states[3].time = dtime;
1527             /*fra=cfra;*/
1528             break;
1529           case 3:
1530             add_v3_v3v3(dx[3], states->vel, dv[2]);
1531             mul_v3_fl(dx[3], dtime);
1532             copy_v3_v3(dv[3], acceleration);
1533             mul_v3_fl(dv[3], dtime);
1534
1535             madd_v3_v3v3fl(pa->state.co, states->co, dx[0], 1.0f / 6.0f);
1536             madd_v3_v3fl(pa->state.co, dx[1], 1.0f / 3.0f);
1537             madd_v3_v3fl(pa->state.co, dx[2], 1.0f / 3.0f);
1538             madd_v3_v3fl(pa->state.co, dx[3], 1.0f / 6.0f);
1539
1540             madd_v3_v3v3fl(pa->state.vel, states->vel, dv[0], 1.0f / 6.0f);
1541             madd_v3_v3fl(pa->state.vel, dv[1], 1.0f / 3.0f);
1542             madd_v3_v3fl(pa->state.vel, dv[2], 1.0f / 3.0f);
1543             madd_v3_v3fl(pa->state.vel, dv[3], 1.0f / 6.0f);
1544         }
1545         break;
1546       case PART_INT_VERLET: /* Verlet integration */
1547         madd_v3_v3v3fl(pa->state.vel, pa->prev_state.vel, acceleration, dtime);
1548         madd_v3_v3v3fl(pa->state.co, pa->prev_state.co, pa->state.vel, dtime);
1549
1550         sub_v3_v3v3(pa->state.vel, pa->state.co, oldpos);
1551         mul_v3_fl(pa->state.vel, 1.0f / dtime);
1552         break;
1553     }
1554   }
1555 }
1556
1557 /*********************************************************************************************************
1558  *                    SPH fluid physics
1559  *
1560  * In theory, there could be unlimited implementation of SPH simulators
1561  *
1562  * This code uses in some parts adapted algorithms from the pseudo code as outlined in the Research paper:
1563  *
1564  * Titled: Particle-based Viscoelastic Fluid Simulation.
1565  * Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
1566  * Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
1567  *
1568  * Presented at Siggraph, (2005)
1569  *
1570  * ********************************************************************************************************/
1571 #define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
1572 static ParticleSpring *sph_spring_add(ParticleSystem *psys, ParticleSpring *spring)
1573 {
1574   /* Are more refs required? */
1575   if (psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
1576     psys->alloc_fluidsprings = PSYS_FLUID_SPRINGS_INITIAL_SIZE;
1577     psys->fluid_springs = (ParticleSpring *)MEM_callocN(
1578         psys->alloc_fluidsprings * sizeof(ParticleSpring), "Particle Fluid Springs");
1579   }
1580   else if (psys->tot_fluidsprings == psys->alloc_fluidsprings) {
1581     /* Double the number of refs allocated */
1582     psys->alloc_fluidsprings *= 2;
1583     psys->fluid_springs = (ParticleSpring *)MEM_reallocN(
1584         psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
1585   }
1586
1587   memcpy(psys->fluid_springs + psys->tot_fluidsprings, spring, sizeof(ParticleSpring));
1588   psys->tot_fluidsprings++;
1589
1590   return psys->fluid_springs + psys->tot_fluidsprings - 1;
1591 }
1592 static void sph_spring_delete(ParticleSystem *psys, int j)
1593 {
1594   if (j != psys->tot_fluidsprings - 1) {
1595     psys->fluid_springs[j] = psys->fluid_springs[psys->tot_fluidsprings - 1];
1596   }
1597
1598   psys->tot_fluidsprings--;
1599
1600   if (psys->tot_fluidsprings < psys->alloc_fluidsprings / 2 &&
1601       psys->alloc_fluidsprings > PSYS_FLUID_SPRINGS_INITIAL_SIZE) {
1602     psys->alloc_fluidsprings /= 2;
1603     psys->fluid_springs = (ParticleSpring *)MEM_reallocN(
1604         psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
1605   }
1606 }
1607 static void sph_springs_modify(ParticleSystem *psys, float dtime)
1608 {
1609   SPHFluidSettings *fluid = psys->part->fluid;
1610   ParticleData *pa1, *pa2;
1611   ParticleSpring *spring = psys->fluid_springs;
1612
1613   float h, d, Rij[3], rij, Lij;
1614   int i;
1615
1616   float yield_ratio = fluid->yield_ratio;
1617   float plasticity = fluid->plasticity_constant;
1618   /* scale things according to dtime */
1619   float timefix = 25.f * dtime;
1620
1621   if ((fluid->flag & SPH_VISCOELASTIC_SPRINGS) == 0 || fluid->spring_k == 0.f) {
1622     return;
1623   }
1624
1625   /* Loop through the springs */
1626   for (i = 0; i < psys->tot_fluidsprings; i++, spring++) {
1627     pa1 = psys->particles + spring->particle_index[0];
1628     pa2 = psys->particles + spring->particle_index[1];
1629
1630     sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
1631     rij = normalize_v3(Rij);
1632
1633     /* adjust rest length */
1634     Lij = spring->rest_length;
1635     d = yield_ratio * timefix * Lij;
1636
1637     if (rij > Lij + d) {  // Stretch
1638       spring->rest_length += plasticity * (rij - Lij - d) * timefix;
1639     }
1640     else if (rij < Lij - d) {  // Compress
1641       spring->rest_length -= plasticity * (Lij - d - rij) * timefix;
1642     }
1643
1644     h = 4.f * pa1->size;
1645
1646     if (spring->rest_length > h) {
1647       spring->delete_flag = 1;
1648     }
1649   }
1650
1651   /* Loop through springs backwaqrds - for efficient delete function */
1652   for (i = psys->tot_fluidsprings - 1; i >= 0; i--) {
1653     if (psys->fluid_springs[i].delete_flag) {
1654       sph_spring_delete(psys, i);
1655     }
1656   }
1657 }
1658 static EdgeHash *sph_springhash_build(ParticleSystem *psys)
1659 {
1660   EdgeHash *springhash = NULL;
1661   ParticleSpring *spring;
1662   int i = 0;
1663
1664   springhash = BLI_edgehash_new_ex(__func__, psys->tot_fluidsprings);
1665
1666   for (i = 0, spring = psys->fluid_springs; i < psys->tot_fluidsprings; i++, spring++) {
1667     BLI_edgehash_insert(
1668         springhash, spring->particle_index[0], spring->particle_index[1], POINTER_FROM_INT(i + 1));
1669   }
1670
1671   return springhash;
1672 }
1673
1674 #define SPH_NEIGHBORS 512
1675 typedef struct SPHNeighbor {
1676   ParticleSystem *psys;
1677   int index;
1678 } SPHNeighbor;
1679
1680 typedef struct SPHRangeData {
1681   SPHNeighbor neighbors[SPH_NEIGHBORS];
1682   int tot_neighbors;
1683
1684   float *data;
1685
1686   ParticleSystem *npsys;
1687   ParticleData *pa;
1688
1689   float h;
1690   float mass;
1691   float massfac;
1692   int use_size;
1693 } SPHRangeData;
1694
1695 static void sph_evaluate_func(BVHTree *tree,
1696                               ParticleSystem **psys,
1697                               float co[3],
1698                               SPHRangeData *pfr,
1699                               float interaction_radius,
1700                               BVHTree_RangeQuery callback)
1701 {
1702   int i;
1703
1704   pfr->tot_neighbors = 0;
1705
1706   for (i = 0; i < 10 && psys[i]; i++) {
1707     pfr->npsys = psys[i];
1708     pfr->massfac = psys[i]->part->mass / pfr->mass;
1709     pfr->use_size = psys[i]->part->flag & PART_SIZEMASS;
1710
1711     if (tree) {
1712       BLI_bvhtree_range_query(tree, co, interaction_radius, callback, pfr);
1713       break;
1714     }
1715     else {
1716       BLI_rw_mutex_lock(&psys_bvhtree_rwlock, THREAD_LOCK_READ);
1717
1718       BLI_bvhtree_range_query(psys[i]->bvhtree, co, interaction_radius, callback, pfr);
1719
1720       BLI_rw_mutex_unlock(&psys_bvhtree_rwlock);
1721     }
1722   }
1723 }
1724 static void sph_density_accum_cb(void *userdata, int index, const float co[3], float squared_dist)
1725 {
1726   SPHRangeData *pfr = (SPHRangeData *)userdata;
1727   ParticleData *npa = pfr->npsys->particles + index;
1728   float q;
1729   float dist;
1730
1731   UNUSED_VARS(co);
1732
1733   if (npa == pfr->pa || squared_dist < FLT_EPSILON) {
1734     return;
1735   }
1736
1737   /* Ugh! One particle has too many neighbors! If some aren't taken into
1738    * account, the forces will be biased by the tree search order. This
1739    * effectively adds energy to the system, and results in a churning motion.
1740    * But, we have to stop somewhere, and it's not the end of the world.
1741    * - jahka and z0r
1742    */
1743   if (pfr->tot_neighbors >= SPH_NEIGHBORS) {
1744     return;
1745   }
1746
1747   pfr->neighbors[pfr->tot_neighbors].index = index;
1748   pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
1749   pfr->tot_neighbors++;
1750
1751   dist = sqrtf(squared_dist);
1752   q = (1.f - dist / pfr->h) * pfr->massfac;
1753
1754   if (pfr->use_size) {
1755     q *= npa->size;
1756   }
1757
1758   pfr->data[0] += q * q;
1759   pfr->data[1] += q * q * q;
1760 }
1761
1762 /*
1763  * Find the Courant number for an SPH particle (used for adaptive time step).
1764  */
1765 static void sph_particle_courant(SPHData *sphdata, SPHRangeData *pfr)
1766 {
1767   ParticleData *pa, *npa;
1768   int i;
1769   float flow[3], offset[3], dist;
1770
1771   zero_v3(flow);
1772
1773   dist = 0.0f;
1774   if (pfr->tot_neighbors > 0) {
1775     pa = pfr->pa;
1776     for (i = 0; i < pfr->tot_neighbors; i++) {
1777       npa = pfr->neighbors[i].psys->particles + pfr->neighbors[i].index;
1778       sub_v3_v3v3(offset, pa->prev_state.co, npa->prev_state.co);
1779       dist += len_v3(offset);
1780       add_v3_v3(flow, npa->prev_state.vel);
1781     }
1782     dist += sphdata->psys[0]->part->fluid->radius;  // TODO: remove this? - z0r
1783     sphdata->element_size = dist / pfr->tot_neighbors;
1784     mul_v3_v3fl(sphdata->flow, flow, 1.0f / pfr->tot_neighbors);
1785   }
1786   else {
1787     sphdata->element_size = FLT_MAX;
1788     copy_v3_v3(sphdata->flow, flow);
1789   }
1790 }
1791 static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse))
1792 {
1793   SPHData *sphdata = (SPHData *)sphdata_v;
1794   ParticleSystem **psys = sphdata->psys;
1795   ParticleData *pa = sphdata->pa;
1796   SPHFluidSettings *fluid = psys[0]->part->fluid;
1797   ParticleSpring *spring = NULL;
1798   SPHRangeData pfr;
1799   SPHNeighbor *pfn;
1800   float *gravity = sphdata->gravity;
1801   EdgeHash *springhash = sphdata->eh;
1802
1803   float q, u, rij, dv[3];
1804   float pressure, near_pressure;
1805
1806   float visc = fluid->viscosity_omega;
1807   float stiff_visc = fluid->viscosity_beta *
1808                      (fluid->flag & SPH_FAC_VISCOSITY ? fluid->viscosity_omega : 1.f);
1809
1810   float inv_mass = 1.0f / sphdata->mass;
1811   float spring_constant = fluid->spring_k;
1812
1813   /* 4.0 seems to be a pretty good value */
1814   float interaction_radius = fluid->radius *
1815                              (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
1816   float h = interaction_radius * sphdata->hfac;
1817   /* 4.77 is an experimentally determined density factor */
1818   float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f);
1819   float rest_length = fluid->rest_length *
1820                       (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f);
1821
1822   float stiffness = fluid->stiffness_k;
1823   float stiffness_near_fac = fluid->stiffness_knear *
1824                              (fluid->flag & SPH_FAC_REPULSION ? fluid->stiffness_k : 1.f);
1825
1826   ParticleData *npa;
1827   float vec[3];
1828   float vel[3];
1829   float co[3];
1830   float data[2];
1831   float density, near_density;
1832
1833   int i, spring_index, index = pa - psys[0]->particles;
1834
1835   data[0] = data[1] = 0;
1836   pfr.data = data;
1837   pfr.h = h;
1838   pfr.pa = pa;
1839   pfr.mass = sphdata->mass;
1840
1841   sph_evaluate_func(NULL, psys, state->co, &pfr, interaction_radius, sph_density_accum_cb);
1842
1843   density = data[0];
1844   near_density = data[1];
1845
1846   pressure = stiffness * (density - rest_density);
1847   near_pressure = stiffness_near_fac * near_density;
1848
1849   pfn = pfr.neighbors;
1850   for (i = 0; i < pfr.tot_neighbors; i++, pfn++) {
1851     npa = pfn->psys->particles + pfn->index;
1852
1853     madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
1854
1855     sub_v3_v3v3(vec, co, state->co);
1856     rij = normalize_v3(vec);
1857
1858     q = (1.f - rij / h) * pfn->psys->part->mass * inv_mass;
1859
1860     if (pfn->psys->part->flag & PART_SIZEMASS) {
1861       q *= npa->size;
1862     }
1863
1864     copy_v3_v3(vel, npa->prev_state.vel);
1865
1866     /* Double Density Relaxation */
1867     madd_v3_v3fl(force, vec, -(pressure + near_pressure * q) * q);
1868
1869     /* Viscosity */
1870     if (visc > 0.f || stiff_visc > 0.f) {
1871       sub_v3_v3v3(dv, vel, state->vel);
1872       u = dot_v3v3(vec, dv);
1873
1874       if (u < 0.f && visc > 0.f) {
1875         madd_v3_v3fl(force, vec, 0.5f * q * visc * u);
1876       }
1877
1878       if (u > 0.f && stiff_visc > 0.f) {
1879         madd_v3_v3fl(force, vec, 0.5f * q * stiff_visc * u);
1880       }
1881     }
1882
1883     if (spring_constant > 0.f) {
1884       /* Viscoelastic spring force */
1885       if (pfn->psys == psys[0] && fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash) {
1886         /* BLI_edgehash_lookup appears to be thread-safe. - z0r */
1887         spring_index = POINTER_AS_INT(BLI_edgehash_lookup(springhash, index, pfn->index));
1888
1889         if (spring_index) {
1890           spring = psys[0]->fluid_springs + spring_index - 1;
1891
1892           madd_v3_v3fl(
1893               force, vec, -10.f * spring_constant * (1.f - rij / h) * (spring->rest_length - rij));
1894         }
1895         else if (fluid->spring_frames == 0 ||
1896                  (pa->prev_state.time - pa->time) <= fluid->spring_frames) {
1897           ParticleSpring temp_spring;
1898           temp_spring.particle_index[0] = index;
1899           temp_spring.particle_index[1] = pfn->index;
1900           temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : rest_length;
1901           temp_spring.delete_flag = 0;
1902
1903           /* sph_spring_add is not thread-safe. - z0r */
1904           sph_spring_add(psys[0], &temp_spring);
1905         }
1906       }
1907       else { /* PART_SPRING_HOOKES - Hooke's spring force */
1908         madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij / h) * (rest_length - rij));
1909       }
1910     }
1911   }
1912
1913   /* Artificial buoyancy force in negative gravity direction  */
1914   if (fluid->buoyancy > 0.f && gravity) {
1915     madd_v3_v3fl(force, gravity, fluid->buoyancy * (density - rest_density));
1916   }
1917
1918   if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF) {
1919     sph_particle_courant(sphdata, &pfr);
1920   }
1921   sphdata->pass++;
1922 }
1923
1924 static void sphclassical_density_accum_cb(void *userdata,
1925                                           int index,
1926                                           const float co[3],
1927                                           float UNUSED(squared_dist))
1928 {
1929   SPHRangeData *pfr = (SPHRangeData *)userdata;
1930   ParticleData *npa = pfr->npsys->particles + index;
1931   float q;
1932   float qfac = 21.0f / (256.f * (float)M_PI);
1933   float rij, rij_h;
1934   float vec[3];
1935
1936   /* Exclude particles that are more than 2h away. Can't use squared_dist here
1937    * because it is not accurate enough. Use current state, i.e. the output of
1938    * basic_integrate() - z0r */
1939   sub_v3_v3v3(vec, npa->state.co, co);
1940   rij = len_v3(vec);
1941   rij_h = rij / pfr->h;
1942   if (rij_h > 2.0f) {
1943     return;
1944   }
1945
1946   /* Smoothing factor. Utilise the Wendland kernel. gnuplot:
1947    *     q1(x) = (2.0 - x)**4 * ( 1.0 + 2.0 * x)
1948    *     plot [0:2] q1(x) */
1949   q = qfac / pow3f(pfr->h) * pow4f(2.0f - rij_h) * (1.0f + 2.0f * rij_h);
1950   q *= pfr->npsys->part->mass;
1951
1952   if (pfr->use_size) {
1953     q *= pfr->pa->size;
1954   }
1955
1956   pfr->data[0] += q;
1957   pfr->data[1] += q / npa->sphdensity;
1958 }
1959
1960 static void sphclassical_neighbour_accum_cb(void *userdata,
1961                                             int index,
1962                                             const float co[3],
1963                                             float UNUSED(squared_dist))
1964 {
1965   SPHRangeData *pfr = (SPHRangeData *)userdata;
1966   ParticleData *npa = pfr->npsys->particles + index;
1967   float rij, rij_h;
1968   float vec[3];
1969
1970   if (pfr->tot_neighbors >= SPH_NEIGHBORS) {
1971     return;
1972   }
1973
1974   /* Exclude particles that are more than 2h away. Can't use squared_dist here
1975    * because it is not accurate enough. Use current state, i.e. the output of
1976    * basic_integrate() - z0r */
1977   sub_v3_v3v3(vec, npa->state.co, co);
1978   rij = len_v3(vec);
1979   rij_h = rij / pfr->h;
1980   if (rij_h > 2.0f) {
1981     return;
1982   }
1983
1984   pfr->neighbors[pfr->tot_neighbors].index = index;
1985   pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
1986   pfr->tot_neighbors++;
1987 }
1988 static void sphclassical_force_cb(void *sphdata_v,
1989                                   ParticleKey *state,
1990                                   float *force,
1991                                   float *UNUSED(impulse))
1992 {
1993   SPHData *sphdata = (SPHData *)sphdata_v;
1994   ParticleSystem **psys = sphdata->psys;
1995   ParticleData *pa = sphdata->pa;
1996   SPHFluidSettings *fluid = psys[0]->part->fluid;
1997   SPHRangeData pfr;
1998   SPHNeighbor *pfn;
1999   float *gravity = sphdata->gravity;
2000
2001   float dq, u, rij, dv[3];
2002   float pressure, npressure;
2003
2004   float visc = fluid->viscosity_omega;
2005
2006   float interaction_radius;
2007   float h, hinv;
2008   /* 4.77 is an experimentally determined density factor */
2009   float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.0f);
2010
2011   // Use speed of sound squared
2012   float stiffness = pow2f(fluid->stiffness_k);
2013
2014   ParticleData *npa;
2015   float vec[3];
2016   float co[3];
2017   float pressureTerm;
2018
2019   int i;
2020
2021   float qfac2 = 42.0f / (256.0f * (float)M_PI);
2022   float rij_h;
2023
2024   /* 4.0 here is to be consistent with previous formulation/interface */
2025   interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
2026   h = interaction_radius * sphdata->hfac;
2027   hinv = 1.0f / h;
2028
2029   pfr.h = h;
2030   pfr.pa = pa;
2031
2032   sph_evaluate_func(
2033       NULL, psys, state->co, &pfr, interaction_radius, sphclassical_neighbour_accum_cb);
2034   pressure = stiffness * (pow7f(pa->sphdensity / rest_density) - 1.0f);
2035
2036   /* multiply by mass so that we return a force, not accel */
2037   qfac2 *= sphdata->mass / pow3f(pfr.h);
2038
2039   pfn = pfr.neighbors;
2040   for (i = 0; i < pfr.tot_neighbors; i++, pfn++) {
2041     npa = pfn->psys->particles + pfn->index;
2042     if (npa == pa) {
2043       /* we do not contribute to ourselves */
2044       continue;
2045     }
2046
2047     /* Find vector to neighbor. Exclude particles that are more than 2h
2048      * away. Can't use current state here because it may have changed on
2049      * another thread - so do own mini integration. Unlike basic_integrate,
2050      * SPH integration depends on neighboring particles. - z0r */
2051     madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
2052     sub_v3_v3v3(vec, co, state->co);
2053     rij = normalize_v3(vec);
2054     rij_h = rij / pfr.h;
2055     if (rij_h > 2.0f) {
2056       continue;
2057     }
2058
2059     npressure = stiffness * (pow7f(npa->sphdensity / rest_density) - 1.0f);
2060
2061     /* First derivative of smoothing factor. Utilise the Wendland kernel.
2062      * gnuplot:
2063      *     q2(x) = 2.0 * (2.0 - x)**4 - 4.0 * (2.0 - x)**3 * (1.0 + 2.0 * x)
2064      *     plot [0:2] q2(x)
2065      * Particles > 2h away are excluded above. */
2066     dq = qfac2 * (2.0f * pow4f(2.0f - rij_h) - 4.0f * pow3f(2.0f - rij_h) * (1.0f + 2.0f * rij_h));
2067
2068     if (pfn->psys->part->flag & PART_SIZEMASS) {
2069       dq *= npa->size;
2070     }
2071
2072     pressureTerm = pressure / pow2f(pa->sphdensity) + npressure / pow2f(npa->sphdensity);
2073
2074     /* Note that 'minus' is removed, because vec = vecBA, not vecAB.
2075      * This applies to the viscosity calculation below, too. */
2076     madd_v3_v3fl(force, vec, pressureTerm * dq);
2077
2078     /* Viscosity */
2079     if (visc > 0.0f) {
2080       sub_v3_v3v3(dv, npa->prev_state.vel, pa->prev_state.vel);
2081       u = dot_v3v3(vec, dv);
2082       /* Apply parameters */
2083       u *= -dq * hinv * visc / (0.5f * npa->sphdensity + 0.5f * pa->sphdensity);
2084       madd_v3_v3fl(force, vec, u);
2085     }
2086   }
2087
2088   /* Artificial buoyancy force in negative gravity direction  */
2089   if (fluid->buoyancy > 0.f && gravity) {
2090     madd_v3_v3fl(force, gravity, fluid->buoyancy * (pa->sphdensity - rest_density));
2091   }
2092
2093   if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF) {
2094     sph_particle_courant(sphdata, &pfr);
2095   }
2096   sphdata->pass++;
2097 }
2098
2099 static void sphclassical_calc_dens(ParticleData *pa, float UNUSED(dfra), SPHData *sphdata)
2100 {
2101   ParticleSystem **psys = sphdata->psys;
2102   SPHFluidSettings *fluid = psys[0]->part->fluid;
2103   /* 4.0 seems to be a pretty good value */
2104   float interaction_radius = fluid->radius *
2105                              (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
2106   SPHRangeData pfr;
2107   float data[2];
2108
2109   data[0] = 0;
2110   data[1] = 0;
2111   pfr.data = data;
2112   pfr.h = interaction_radius * sphdata->hfac;
2113   pfr.pa = pa;
2114   pfr.mass = sphdata->mass;
2115
2116   sph_evaluate_func(
2117       NULL, psys, pa->state.co, &pfr, interaction_radius, sphclassical_density_accum_cb);
2118   pa->sphdensity = min_ff(max_ff(data[0], fluid->rest_density * 0.9f), fluid->rest_density * 1.1f);
2119 }
2120
2121 void psys_sph_init(ParticleSimulationData *sim, SPHData *sphdata)
2122 {
2123   ParticleTarget *pt;
2124   int i;
2125
2126   // Add other coupled particle systems.
2127   sphdata->psys[0] = sim->psys;
2128   for (i = 1, pt = sim->psys->targets.first; i < 10; i++, pt = (pt ? pt->next : NULL)) {
2129     sphdata->psys[i] = pt ? psys_get_target_system(sim->ob, pt) : NULL;
2130   }
2131
2132   if (psys_uses_gravity(sim)) {
2133     sphdata->gravity = sim->scene->physics_settings.gravity;
2134   }
2135   else {
2136     sphdata->gravity = NULL;
2137   }
2138   sphdata->eh = sph_springhash_build(sim->psys);
2139
2140   // These per-particle values should be overridden later, but just for
2141   // completeness we give them default values now.
2142   sphdata->pa = NULL;
2143   sphdata->mass = 1.0f;
2144
2145   if (sim->psys->part->fluid->solver == SPH_SOLVER_DDR) {
2146     sphdata->force_cb = sph_force_cb;
2147     sphdata->density_cb = sph_density_accum_cb;
2148     sphdata->hfac = 1.0f;
2149   }
2150   else {
2151     /* SPH_SOLVER_CLASSICAL */
2152     sphdata->force_cb = sphclassical_force_cb;
2153     sphdata->density_cb = sphclassical_density_accum_cb;
2154     sphdata->hfac = 0.5f;
2155   }
2156 }
2157
2158 void psys_sph_finalise(SPHData *sphdata)
2159 {
2160   if (sphdata->eh) {
2161     BLI_edgehash_free(sphdata->eh, NULL);
2162     sphdata->eh = NULL;
2163   }
2164 }
2165 /* Sample the density field at a point in space. */
2166 void psys_sph_density(BVHTree *tree, SPHData *sphdata, float co[3], float vars[2])
2167 {
2168   ParticleSystem **psys = sphdata->psys;
2169   SPHFluidSettings *fluid = psys[0]->part->fluid;
2170   /* 4.0 seems to be a pretty good value */
2171   float interaction_radius = fluid->radius *
2172                              (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
2173   SPHRangeData pfr;
2174   float density[2];
2175
2176   density[0] = density[1] = 0.0f;
2177   pfr.data = density;
2178   pfr.h = interaction_radius * sphdata->hfac;
2179   pfr.mass = sphdata->mass;
2180
2181   sph_evaluate_func(tree, psys, co, &pfr, interaction_radius, sphdata->density_cb);
2182
2183   vars[0] = pfr.data[0];
2184   vars[1] = pfr.data[1];
2185 }
2186
2187 static void sph_integrate(ParticleSimulationData *sim,
2188                           ParticleData *pa,
2189                           float dfra,
2190                           SPHData *sphdata)
2191 {
2192   ParticleSettings *part = sim->psys->part;
2193   // float timestep = psys_get_timestep(sim); // UNUSED
2194   float pa_mass = part->mass * (part->flag & PART_SIZEMASS ? pa->size : 1.f);
2195   float dtime = dfra * psys_get_timestep(sim);
2196   // int steps = 1; // UNUSED
2197   float effector_acceleration[3];
2198
2199   sphdata->pa = pa;
2200   sphdata->mass = pa_mass;
2201   sphdata->pass = 0;
2202   //sphdata.element_size and sphdata.flow are set in the callback.
2203
2204   /* restore previous state and treat gravity & effectors as external acceleration*/
2205   sub_v3_v3v3(effector_acceleration, pa->state.vel, pa->prev_state.vel);
2206   mul_v3_fl(effector_acceleration, 1.f / dtime);
2207
2208   copy_particle_key(&pa->state, &pa->prev_state, 0);
2209
2210   integrate_particle(part, pa, dtime, effector_acceleration, sphdata->force_cb, sphdata);
2211 }
2212
2213 /************************************************/
2214 /*          Basic physics                       */
2215 /************************************************/
2216 typedef struct EfData {
2217   ParticleTexture ptex;
2218   ParticleSimulationData *sim;
2219   ParticleData *pa;
2220 } EfData;
2221 static void basic_force_cb(void *efdata_v, ParticleKey *state, float *force, float *impulse)
2222 {
2223   EfData *efdata = (EfData *)efdata_v;
2224   ParticleSimulationData *sim = efdata->sim;
2225   ParticleSettings *part = sim->psys->part;
2226   ParticleData *pa = efdata->pa;
2227   EffectedPoint epoint;
2228   RNG *rng = sim->rng;
2229
2230   /* add effectors */
2231   pd_point_from_particle(efdata->sim, efdata->pa, state, &epoint);
2232   if (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR) {
2233     BKE_effectors_apply(
2234         sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);
2235   }
2236
2237   mul_v3_fl(force, efdata->ptex.field);
2238   mul_v3_fl(impulse, efdata->ptex.field);
2239
2240   /* calculate air-particle interaction */
2241   if (part->dragfac != 0.0f) {
2242     madd_v3_v3fl(force, state->vel, -part->dragfac * pa->size * pa->size * len_v3(state->vel));
2243   }
2244
2245   /* brownian force */
2246   if (part->brownfac != 0.0f) {
2247     force[0] += (BLI_rng_get_float(rng) - 0.5f) * part->brownfac;
2248     force[1] += (BLI_rng_get_float(rng) - 0.5f) * part->brownfac;
2249     force[2] += (BLI_rng_get_float(rng) - 0.5f) * part->brownfac;
2250   }
2251
2252   if (part->flag & PART_ROT_DYN && epoint.ave) {
2253     copy_v3_v3(pa->state.ave, epoint.ave);
2254   }
2255 }
2256 /* gathers all forces that effect particles and calculates a new state for the particle */
2257 static void basic_integrate(ParticleSimulationData *sim, int p, float dfra, float cfra)
2258 {
2259   ParticleSettings *part = sim->psys->part;
2260   ParticleData *pa = sim->psys->particles + p;
2261   ParticleKey tkey;
2262   float dtime = dfra * psys_get_timestep(sim), time;
2263   float *gravity = NULL, gr[3];
2264   EfData efdata;
2265
2266   psys_get_texture(sim, pa, &efdata.ptex, PAMAP_PHYSICS, cfra);
2267
2268   efdata.pa = pa;
2269   efdata.sim = sim;
2270
2271   /* add global acceleration (gravitation) */
2272   if (psys_uses_gravity(sim) &&
2273       /* normal gravity is too strong for hair so it's disabled by default */
2274       (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
2275     zero_v3(gr);
2276     madd_v3_v3fl(gr,
2277                  sim->scene->physics_settings.gravity,
2278                  part->effector_weights->global_gravity * efdata.ptex.gravity);
2279     gravity = gr;
2280   }
2281
2282   /* maintain angular velocity */
2283   copy_v3_v3(pa->state.ave, pa->prev_state.ave);
2284
2285   integrate_particle(part, pa, dtime, gravity, basic_force_cb, &efdata);
2286
2287   /* damp affects final velocity */
2288   if (part->dampfac != 0.f) {
2289     mul_v3_fl(pa->state.vel, 1.f - part->dampfac * efdata.ptex.damp * 25.f * dtime);
2290   }
2291
2292   //copy_v3_v3(pa->state.ave, states->ave);
2293
2294   /* finally we do guides */
2295   time = (cfra - pa->time) / pa->lifetime;
2296   CLAMP(time, 0.0f, 1.0f);
2297
2298   copy_v3_v3(tkey.co, pa->state.co);
2299   copy_v3_v3(tkey.vel, pa->state.vel);
2300   tkey.time = pa->state.time;
2301
2302   if (part->type != PART_HAIR) {
2303     if (do_guides(sim->depsgraph, sim->psys->part, sim->psys->effectors, &tkey, p, time)) {
2304       copy_v3_v3(pa->state.co, tkey.co);
2305       /* guides don't produce valid velocity */
2306       sub_v3_v3v3(pa->state.vel, tkey.co, pa->prev_state.co);
2307       mul_v3_fl(pa->state.vel, 1.0f / dtime);
2308       pa->state.time = tkey.time;
2309     }
2310   }
2311 }
2312 static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
2313 {
2314   float rotfac, rot1[4], rot2[4] = {1.0, 0.0, 0.0, 0.0}, dtime = dfra * timestep, extrotfac;
2315
2316   if ((part->flag & PART_ROTATIONS) == 0) {
2317     unit_qt(pa->state.rot);
2318     return;
2319   }
2320
2321   if (part->flag & PART_ROT_DYN) {
2322     extrotfac = len_v3(pa->state.ave);
2323   }
2324   else {
2325     extrotfac = 0.0f;
2326   }
2327
2328   if ((part->flag & PART_ROT_DYN) &&
2329       ELEM(part->avemode, PART_AVE_VELOCITY, PART_AVE_HORIZONTAL, PART_AVE_VERTICAL)) {
2330     float angle;
2331     float len1 = len_v3(pa->prev_state.vel);
2332     float len2 = len_v3(pa->state.vel);
2333     float vec[3];
2334
2335     if (len1 == 0.0f || len2 == 0.0f) {
2336       zero_v3(pa->state.ave);
2337     }
2338     else {
2339       cross_v3_v3v3(pa->state.ave, pa->prev_state.vel, pa->state.vel);
2340       normalize_v3(pa->state.ave);
2341       angle = dot_v3v3(pa->prev_state.vel, pa->state.vel) / (len1 * len2);
2342       mul_v3_fl(pa->state.ave, saacos(angle) / dtime);
2343     }
2344
2345     get_angular_velocity_vector(part->avemode, &pa->state, vec);
2346     axis_angle_to_quat(rot2, vec, dtime * part->avefac);
2347   }
2348
2349   rotfac = len_v3(pa->state.ave);
2350   if (rotfac == 0.0f || (part->flag & PART_ROT_DYN) == 0 || extrotfac == 0.0f) {
2351     unit_qt(rot1);
2352   }
2353   else {
2354     axis_angle_to_quat(rot1, pa->state.ave, rotfac * dtime);
2355   }
2356   mul_qt_qtqt(pa->state.rot, rot1, pa->prev_state.rot);
2357   mul_qt_qtqt(pa->state.rot, rot2, pa->state.rot);
2358
2359   /* keep rotation quat in good health */
2360   normalize_qt(pa->state.rot);
2361 }
2362
2363 /************************************************
2364  * Collisions
2365  *
2366  * The algorithm is roughly:
2367  *  1. Use a BVH tree to search for faces that a particle may collide with.
2368  *  2. Use Newton's method to find the exact time at which the collision occurs.
2369  *     https://en.wikipedia.org/wiki/Newton's_method
2370  *
2371  ************************************************/
2372 #define COLLISION_MIN_RADIUS 0.001f
2373 #define COLLISION_MIN_DISTANCE 0.0001f
2374 #define COLLISION_ZERO 0.00001f
2375 #define COLLISION_INIT_STEP 0.00008f
2376 typedef float (*NRDistanceFunc)(float *p, float radius, ParticleCollisionElement *pce, float *nor);
2377 static float nr_signed_distance_to_plane(float *p,
2378                                          float radius,
2379                                          ParticleCollisionElement *pce,
2380                                          float *nor)
2381 {
2382   float p0[3], e1[3], e2[3], d;
2383
2384   sub_v3_v3v3(e1, pce->x1, pce->x0);
2385   sub_v3_v3v3(e2, pce->x2, pce->x0);
2386   sub_v3_v3v3(p0, p, pce->x0);
2387
2388   cross_v3_v3v3(nor, e1, e2);
2389   normalize_v3(nor);
2390
2391   d = dot_v3v3(p0, nor);
2392
2393   if (pce->inv_nor == -1) {
2394     if (d < 0.f) {
2395       pce->inv_nor = 1;
2396     }
2397     else {
2398       pce->inv_nor = 0;
2399     }
2400   }
2401
2402   if (pce->inv_nor == 1) {
2403     negate_v3(nor);
2404     d = -d;
2405   }
2406
2407   return d - radius;
2408 }
2409 static float nr_distance_to_edge(float *p,
2410                                  float radius,
2411                                  ParticleCollisionElement *pce,
2412                                  float *UNUSED(nor))
2413 {
2414   float v0[3], v1[3], v2[3], c[3];
2415
2416   sub_v3_v3v3(v0, pce->x1, pce->x0);
2417   sub_v3_v3v3(v1, p, pce->x0);
2418   sub_v3_v3v3(v2, p, pce->x1);
2419
2420   cross_v3_v3v3(c, v1, v2);
2421
2422   return fabsf(len_v3(c) / len_v3(v0)) - radius;
2423 }
2424 static float nr_distance_to_vert(float *p,
2425                                  float radius,
2426                                  ParticleCollisionElement *pce,
2427                                  float *UNUSED(nor))
2428 {
2429   return len_v3v3(p, pce->x0) - radius;
2430 }
2431 static void collision_interpolate_element(ParticleCollisionElement *pce,
2432                                           float t,
2433                                           float fac,
2434                                           ParticleCollision *col)
2435 {
2436   /* t is the current time for newton rhapson */
2437   /* fac is the starting factor for current collision iteration */
2438   /* the col->fac's are factors for the particle subframe step start and end during collision modifier step */
2439   float f = fac + t * (1.f - fac);
2440   float mul = col->fac1 + f * (col->fac2 - col->fac1);
2441   if (pce->tot > 0) {
2442     madd_v3_v3v3fl(pce->x0, pce->x[0], pce->v[0], mul);
2443
2444     if (pce->tot > 1) {
2445       madd_v3_v3v3fl(pce->x1, pce->x[1], pce->v[1], mul);
2446
2447       if (pce->tot > 2) {
2448         madd_v3_v3v3fl(pce->x2, pce->x[2], pce->v[2], mul);
2449       }
2450     }
2451   }
2452 }
2453 static void collision_point_velocity(ParticleCollisionElement *pce)
2454 {
2455   float v[3];
2456
2457   copy_v3_v3(pce->vel, pce->v[0]);
2458
2459   if (pce->tot > 1) {
2460     sub_v3_v3v3(v, pce->v[1], pce->v[0]);
2461     madd_v3_v3fl(pce->vel, v, pce->uv[0]);
2462
2463     if (pce->tot > 2) {
2464       sub_v3_v3v3(v, pce->v[2], pce->v[0]);
2465       madd_v3_v3fl(pce->vel, v, pce->uv[1]);
2466     }
2467   }
2468 }
2469 static float collision_point_distance_with_normal(
2470     float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *nor)
2471 {
2472   if (fac >= 0.f) {
2473     collision_interpolate_element(pce, 0.f, fac, col);
2474   }
2475
2476   switch (pce->tot) {
2477     case 1: {
2478       sub_v3_v3v3(nor, p, pce->x0);
2479       return normalize_v3(nor);
2480     }
2481     case 2: {
2482       float u, e[3], vec[3];
2483       sub_v3_v3v3(e, pce->x1, pce->x0);
2484       sub_v3_v3v3(vec, p, pce->x0);
2485       u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2486
2487       madd_v3_v3v3fl(nor, vec, e, -u);
2488       return normalize_v3(nor);
2489     }
2490     case 3:
2491       return nr_signed_distance_to_plane(p, 0.f, pce, nor);
2492   }
2493   return 0;
2494 }
2495 static void collision_point_on_surface(
2496     float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *co)
2497 {
2498   collision_interpolate_element(pce, 0.f, fac, col);
2499
2500   switch (pce->tot) {
2501     case 1: {
2502       sub_v3_v3v3(co, p, pce->x0);
2503       normalize_v3(co);
2504       madd_v3_v3v3fl(co, pce->x0, co, col->radius);
2505       break;
2506     }
2507     case 2: {
2508       float u, e[3], vec[3], nor[3];
2509       sub_v3_v3v3(e, pce->x1, pce->x0);
2510       sub_v3_v3v3(vec, p, pce->x0);
2511       u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2512
2513       madd_v3_v3v3fl(nor, vec, e, -u);
2514       normalize_v3(nor);
2515
2516       madd_v3_v3v3fl(co, pce->x0, e, pce->uv[0]);
2517       madd_v3_v3fl(co, nor, col->radius);
2518       break;
2519     }
2520     case 3: {
2521       float p0[3], e1[3], e2[3], nor[3];
2522
2523       sub_v3_v3v3(e1, pce->x1, pce->x0);
2524       sub_v3_v3v3(e2, pce->x2, pce->x0);
2525       sub_v3_v3v3(p0, p, pce->x0);
2526
2527       cross_v3_v3v3(nor, e1, e2);
2528       normalize_v3(nor);
2529
2530       if (pce->inv_nor == 1) {
2531         negate_v3(nor);
2532       }
2533
2534       madd_v3_v3v3fl(co, pce->x0, nor, col->radius);
2535       madd_v3_v3fl(co, e1, pce->uv[0]);
2536       madd_v3_v3fl(co, e2, pce->uv[1]);
2537       break;
2538     }
2539   }
2540 }
2541 /* find first root in range [0-1] starting from 0 */
2542 static float collision_newton_rhapson(ParticleCollision *col,
2543                                       float radius,
2544                                       ParticleCollisionElement *pce,
2545                                       NRDistanceFunc distance_func)
2546 {
2547   float t0, t1, dt_init, d0, d1, dd, n[3];
2548   int iter;
2549
2550   pce->inv_nor = -1;
2551
2552   if (col->inv_total_time > 0.0f) {
2553     /* Initial step size should be small, but not too small or floating point
2554      * precision errors will appear. - z0r */
2555     dt_init = COLLISION_INIT_STEP * col->inv_total_time;
2556   }
2557   else {
2558     dt_init = 0.001f;
2559   }
2560
2561   /* start from the beginning */
2562   t0 = 0.f;
2563   collision_interpolate_element(pce, t0, col->f, col);
2564   d0 = distance_func(col->co1, radius, pce, n);
2565   t1 = dt_init;
2566   d1 = 0.f;
2567
2568   for (iter = 0; iter < 10; iter++) {  //, itersum++) {
2569     /* get current location */
2570     collision_interpolate_element(pce, t1, col->f, col);
2571     interp_v3_v3v3(pce->p, col->co1, col->co2, t1);
2572
2573     d1 = distance_func(pce->p, radius, pce, n);
2574
2575     /* particle already inside face, so report collision */
2576     if (iter == 0 && d0 < 0.f && d0 > -radius) {
2577       copy_v3_v3(pce->p, col->co1);
2578       copy_v3_v3(pce->nor, n);
2579       pce->inside = 1;
2580       return 0.f;
2581     }
2582
2583     /* Zero gradient (no movement relative to element). Can't step from
2584      * here. */
2585     if (d1 == d0) {
2586       /* If first iteration, try from other end where the gradient may be
2587        * greater. Note: code duplicated below. */
2588       if (iter == 0) {
2589         t0 = 1.f;
2590         collision_interpolate_element(pce, t0, col->f, col);
2591         d0 = distance_func(col->co2, radius, pce, n);
2592         t1 = 1.0f - dt_init;
2593         d1 = 0.f;
2594         continue;
2595       }
2596       else {
2597         return -1.f;
2598       }
2599     }
2600
2601     dd = (t1 - t0) / (d1 - d0);
2602
2603     t0 = t1;
2604     d0 = d1;
2605
2606     t1 -= d1 * dd;
2607
2608     /* Particle moving away from plane could also mean a strangely rotating
2609      * face, so check from end. Note: code duplicated above. */
2610     if (iter == 0 && t1 < 0.f) {
2611       t0 = 1.f;
2612       collision_interpolate_element(pce, t0, col->f, col);
2613       d0 = distance_func(col->co2, radius, pce, n);
2614       t1 = 1.0f - dt_init;
2615       d1 = 0.f;
2616       continue;
2617     }
2618     else if (iter == 1 && (t1 < -COLLISION_ZERO || t1 > 1.f)) {
2619       return -1.f;
2620     }
2621
2622     if (d1 <= COLLISION_ZERO && d1 >= -COLLISION_ZERO) {
2623       if (t1 >= -COLLISION_ZERO && t1 <= 1.f) {
2624         if (distance_func == nr_signed_distance_to_plane) {
2625           copy_v3_v3(pce->nor, n);
2626         }
2627
2628         CLAMP(t1, 0.f, 1.f);
2629
2630         return t1;
2631       }
2632       else {
2633         return -1.f;
2634       }
2635     }
2636   }
2637   return -1.0;
2638 }
2639 static int collision_sphere_to_tri(ParticleCollision *col,
2640                                    float radius,
2641                                    ParticleCollisionElement *pce,
2642                                    float *t)
2643 {
2644   ParticleCollisionElement *result = &col->pce;
2645   float ct, u, v;
2646
2647   pce->inv_nor = -1;
2648   pce->inside = 0;
2649
2650   ct = collision_newton_rhapson(col, radius, pce, nr_signed_distance_to_plane);
2651
2652   if (ct >= 0.f && ct < *t && (result->inside == 0 || pce->inside == 1)) {
2653     float e1[3], e2[3], p0[3];
2654     float e1e1, e1e2, e1p0, e2e2, e2p0, inv;
2655
2656     sub_v3_v3v3(e1, pce->x1, pce->x0);
2657     sub_v3_v3v3(e2, pce->x2, pce->x0);
2658     /* XXX: add radius correction here? */
2659     sub_v3_v3v3(p0, pce->p, pce->x0);
2660
2661     e1e1 = dot_v3v3(e1, e1);
2662     e1e2 = dot_v3v3(e1, e2);
2663     e1p0 = dot_v3v3(e1, p0);
2664     e2e2 = dot_v3v3(e2, e2);
2665     e2p0 = dot_v3v3(e2, p0);
2666
2667     inv = 1.f / (e1e1 * e2e2 - e1e2 * e1e2);
2668     u = (e2e2 * e1p0 - e1e2 * e2p0) * inv;
2669     v = (e1e1 * e2p0 - e1e2 * e1p0) * inv;
2670
2671     if (u >= 0.f && u <= 1.f && v >= 0.f && u + v <= 1.f) {
2672       *result = *pce;
2673
2674       /* normal already calculated in pce */
2675
2676       result->uv[0] = u;
2677       result->uv[1] = v;
2678
2679       *t = ct;
2680       return 1;
2681     }
2682   }
2683   return 0;
2684 }
2685 static int collision_sphere_to_edges(ParticleCollision *col,
2686                                      float radius,
2687                                      ParticleCollisionElement *pce,
2688                                      float *t)
2689 {
2690   ParticleCollisionElement edge[3], *cur = NULL, *hit = NULL;
2691   ParticleCollisionElement *result = &col->pce;
2692
2693   float ct;
2694   int i;
2695
2696   for (i = 0; i < 3; i++) {
2697     cur = edge + i;
2698     cur->x[0] = pce->x[i];
2699     cur->x[1] = pce->x[(i + 1) % 3];
2700     cur->v[0] = pce->v[i];
2701     cur->v[1] = pce->v[(i + 1) % 3];
2702     cur->tot = 2;
2703     cur->inside = 0;
2704
2705     ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_edge);
2706
2707     if (ct >= 0.f && ct < *t) {
2708       float u, e[3], vec[3];
2709
2710       sub_v3_v3v3(e, cur->x1, cur->x0);
2711       sub_v3_v3v3(vec, cur->p, cur->x0);
2712       u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2713
2714       if (u < 0.f || u > 1.f) {
2715         break;
2716       }
2717
2718       *result = *cur;
2719
2720       madd_v3_v3v3fl(result->nor, vec, e, -u);
2721       normalize_v3(result->nor);
2722
2723       result->uv[0] = u;
2724
2725       hit = cur;
2726       *t = ct;
2727     }
2728   }
2729
2730   return hit != NULL;
2731 }
2732 static int collision_sphere_to_verts(ParticleCollision *col,
2733                                      float radius,
2734                                      ParticleCollisionElement *pce,
2735                                      float *t)
2736 {
2737   ParticleCollisionElement vert[3], *cur = NULL, *hit = NULL;
2738   ParticleCollisionElement *result = &col->pce;
2739
2740   float ct;
2741   int i;
2742
2743   for (i = 0; i < 3; i++) {
2744     cur = vert + i;
2745     cur->x[0] = pce->x[i];
2746     cur->v[0] = pce->v[i];
2747     cur->tot = 1;
2748     cur->inside = 0;
2749
2750     ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_vert);
2751
2752     if (ct >= 0.f && ct < *t) {
2753       *result = *cur;
2754
2755       sub_v3_v3v3(result->nor, cur->p, cur->x0);
2756       normalize_v3(result->nor);
2757
2758       hit = cur;
2759       *t = ct;
2760     }
2761   }
2762
2763   return hit != NULL;
2764 }
2765 /* Callback for BVHTree near test */
2766 void BKE_psys_collision_neartest_cb(void *userdata,
2767                                     int index,
2768                                     const BVHTreeRay *ray,
2769                                     BVHTreeRayHit *hit)
2770 {
2771   ParticleCollision *col = (ParticleCollision *)userdata;
2772   ParticleCollisionElement pce;
2773   const MVertTri *vt = &col->md->tri[index];
2774   MVert *x = col->md->x;
2775   MVert *v = col->md->current_v;
2776   float t = hit->dist / col->original_ray_length;
2777   int collision = 0;
2778
2779   pce.x[0] = x[vt->tri[0]].co;
2780   pce.x[1] = x[vt->tri[1]].co;
2781   pce.x[2] = x[vt->tri[2]].co;
2782
2783   pce.v[0] = v[vt->tri[0]].co;
2784   pce.v[1] = v[vt->tri[1]].co;
2785   pce.v[2] = v[vt->tri[2]].co;
2786
2787   pce.tot = 3;
2788   pce.inside = 0;
2789   pce.index = index;
2790
2791   collision = collision_sphere_to_tri(col, ray->radius, &pce, &t);
2792   if (col->pce.inside == 0) {
2793     collision += collision_sphere_to_edges(col, ray->radius, &pce, &t);
2794     collision += collision_sphere_to_verts(col, ray->radius, &pce, &t);
2795   }
2796
2797   if (collision) {
2798     hit->dist = col->original_ray_length * t;
2799     hit->index = index;
2800
2801     collision_point_velocity(&col->pce);
2802
2803     col->hit = col->current;
2804   }
2805 }
2806 static int collision_detect(ParticleData *pa,
2807                             ParticleCollision *col,
2808                             BVHTreeRayHit *hit,
2809                             ListBase *colliders)
2810 {
2811   const int raycast_flag = BVH_RAYCAST_DEFAULT & ~(BVH_RAYCAST_WATERTIGHT);
2812   ColliderCache *coll;
2813   float ray_dir[3];
2814
2815   if (BLI_listbase_is_empty(colliders)) {
2816     return 0;
2817   }
2818
2819   sub_v3_v3v3(ray_dir, col->co2, col->co1);
2820   hit->index = -1;
2821   hit->dist = col->original_ray_length = normalize_v3(ray_dir);
2822   col->pce.inside = 0;
2823
2824   /* even if particle is stationary we want to check for moving colliders */
2825   /* if hit.dist is zero the bvhtree_ray_cast will just ignore everything */
2826   if (hit->dist == 0.0f) {
2827     hit->dist = col->original_ray_length = 0.000001f;
2828   }
2829
2830   for (coll = colliders->first; coll; coll = coll->next) {
2831     /* for boids: don't check with current ground object; also skip if permeated */
2832     bool skip = false;
2833
2834     for (int i = 0; i < col->skip_count; i++) {
2835       if (coll->ob == col->skip[i]) {
2836         skip = true;
2837         break;
2838       }
2839     }
2840
2841     if (skip) {
2842       continue;
2843     }
2844
2845     /* particles should not collide with emitter at birth */
2846     if (coll->ob == col->emitter && pa->time < col->cfra && pa->time >= col->old_cfra) {
2847       continue;
2848     }
2849
2850     col->current = coll->ob;
2851     col->md = coll->collmd;
2852     col->fac1 = (col->old_cfra - coll->collmd->time_x) /
2853                 (coll->collmd->time_xnew - coll->collmd->time_x);
2854     col->fac2 = (col->cfra - coll->collmd->time_x) /
2855                 (coll->collmd->time_xnew - coll->collmd->time_x);
2856
2857     if (col->md && col->md->bvhtree) {
2858       BLI_bvhtree_ray_cast_ex(col->md->bvhtree,
2859                               col->co1,
2860                               ray_dir,
2861                               col->radius,
2862                               hit,
2863                               BKE_psys_collision_neartest_cb,
2864                               col,
2865                               raycast_flag);
2866     }
2867   }
2868
2869   return hit->index >= 0;
2870 }
2871 static int collision_response(ParticleSimulationData *sim,
2872                               ParticleData *pa,
2873                               ParticleCollision *col,
2874                               BVHTreeRayHit *hit,
2875                               int kill,
2876                               int dynamic_rotation)
2877 {
2878   ParticleCollisionElement *pce = &col->pce;
2879   PartDeflect *pd = col->hit->pd;
2880   RNG *rng = sim->rng;
2881   /* point of collision */
2882   float co[3];
2883   /* location factor of collision between this iteration */
2884   float x = hit->dist / col->original_ray_length;
2885   /* time factor of collision between timestep */
2886   float f = col->f + x * (1.0f - col->f);
2887   /* time since previous collision (in seconds) */
2888   float dt1 = (f - col->f) * col->total_time;
2889   /* time left after collision (in seconds) */
2890   float dt2 = (1.0f - f) * col->total_time;
2891   /* did particle pass through the collision surface? */
2892   int through = (BLI_rng_get_float(rng) < pd->pdef_perm) ? 1 : 0;
2893
2894   /* calculate exact collision location */
2895   interp_v3_v3v3(co, col->co1, col->co2, x);
2896
2897   /* particle dies in collision */
2898   if (through == 0 && (kill || pd->flag & PDEFLE_KILL_PART)) {
2899     pa->alive = PARS_DYING;
2900     pa->dietime = col->old_cfra + (col->cfra - col->old_cfra) * f;
2901
2902     copy_v3_v3(pa->state.co, co);
2903     interp_v3_v3v3(pa->state.vel, pa->prev_state.vel, pa->state.vel, f);
2904     interp_qt_qtqt(pa->state.rot, pa->prev_state.rot, pa->state.rot, f);
2905     interp_v3_v3v3(pa->state.ave, pa->prev_state.ave, pa->state.ave, f);
2906
2907     /* particle is dead so we don't need to calculate further */
2908     return 0;
2909   }
2910   /* figure out velocity and other data after collision */
2911   else {
2912     /* velocity directly before collision to be modified into velocity directly after collision */
2913     float v0[3];
2914     /* normal component of v0 */
2915     float v0_nor[3];
2916     /* tangential component of v0 */
2917     float v0_tan[3];
2918     /* tangential component of collision surface velocity */
2919     float vc_tan[3];
2920     float v0_dot, vc_dot;
2921     float damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_rng_get_float(rng) - 0.5f);
2922     float frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_rng_get_float(rng) - 0.5f);
2923     float distance, nor[3], dot;
2924
2925     CLAMP(damp, 0.0f, 1.0f);
2926     CLAMP(frict, 0.0f, 1.0f);
2927
2928     /* get exact velocity right before collision */
2929     madd_v3_v3v3fl(v0, col->ve1, col->acc, dt1);
2930
2931     /* convert collider velocity from 1/framestep to 1/s TODO: here we assume 1 frame step for collision modifier */
2932     mul_v3_fl(pce->vel, col->inv_timestep);
2933
2934     /* calculate tangential particle velocity */
2935     v0_dot = dot_v3v3(pce->nor, v0);
2936     madd_v3_v3v3fl(v0_tan, v0, pce->nor, -v0_dot);
2937
2938     /* calculate tangential collider velocity */
2939     vc_dot = dot_v3v3(pce->nor, pce->vel);
2940     madd_v3_v3v3fl(vc_tan, pce->vel, pce->nor, -vc_dot);
2941
2942     /* handle friction effects (tangential and angular velocity) */
2943     if (frict > 0.0f) {
2944       /* angular <-> linear velocity */
2945       if (dynamic_rotation) {
2946         float vr_tan[3], v1_tan[3], ave[3];
2947
2948         /* linear velocity of particle surface */
2949         cross_v3_v3v3(vr_tan, pce->nor, pa->state.ave);
2950         mul_v3_fl(vr_tan, pa->size);
2951
2952         /* change to coordinates that move with the collision plane */
2953         sub_v3_v3v3(v1_tan, v0_tan, vc_tan);
2954
2955         /* The resulting velocity is a weighted average of particle cm & surface
2956          * velocity. This weight (related to particle's moment of inertia) could
2957          * be made a parameter for angular <-> linear conversion.
2958          */
2959         madd_v3_v3fl(v1_tan, vr_tan, -0.4);
2960         mul_v3_fl(v1_tan, 1.0f / 1.4f); /* 1/(1+0.4) */
2961
2962         /* rolling friction is around 0.01 of sliding friction
2963          * (could be made a parameter) */
2964         mul_v3_fl(v1_tan, 1.0f - 0.01f * frict);
2965
2966         /* surface_velocity is opposite to cm velocity */
2967         negate_v3_v3(vr_tan, v1_tan);
2968
2969         /* get back to global coordinates */
2970         add_v3_v3(v1_tan, vc_tan);
2971
2972         /* convert to angular velocity*/
2973         cross_v3_v3v3(ave, vr_tan, pce->nor);
2974         mul_v3_fl(ave, 1.0f / MAX2(pa->size, 0.001f));
2975
2976         /* only friction will cause change in linear & angular velocity */
2977         interp_v3_v3v3(pa->state.ave, pa->state.ave, ave, frict);
2978         interp_v3_v3v3(v0_tan, v0_tan, v1_tan, frict);
2979       }
2980       else {
2981         /* just basic friction (unphysical due to the friction model used in Blender) */
2982         interp_v3_v3v3(v0_tan, v0_tan, vc_tan, frict);
2983       }
2984     }
2985
2986     /* stickiness was possibly added before, so cancel that before calculating new normal velocity */
2987     /* otherwise particles go flying out of the surface because of high reversed sticky velocity */
2988     if (v0_dot < 0.0f) {
2989       v0_dot += pd->pdef_stickness;
2990       if (v0_dot > 0.0f) {
2991         v0_dot = 0.0f;
2992       }
2993     }
2994
2995     /* damping and flipping of velocity around normal */
2996     v0_dot *= 1.0f - damp;
2997     vc_dot *= through ? damp : 1.0f;
2998
2999     /* calculate normal particle velocity */
3000     /* special case for object hitting the particle from behind */
3001     if (through == 0 && ((vc_dot > 0.0f && v0_dot > 0.0f && vc_dot > v0_dot) ||
3002                          (vc_dot < 0.0f && v0_dot < 0.0f && vc_dot < v0_dot))) {
3003       mul_v3_v3fl(v0_nor, pce->nor, vc_dot);
3004     }
3005     else if (v0_dot > 0.f) {
3006       mul_v3_v3fl(v0_nor, pce->nor, vc_dot + v0_dot);
3007     }
3008     else {
3009       mul_v3_v3fl(v0_nor, pce->nor, vc_dot + (through ? 1.0f : -1.0f) * v0_dot);
3010     }
3011
3012     /* combine components together again */
3013     add_v3_v3v3(v0, v0_nor, v0_tan);
3014
3015     if (col->boid) {
3016       /* keep boids above ground */
3017       BoidParticle *bpa = pa->boid;
3018       if (bpa->data.mode == eBoidMode_OnLand || co[2] <= col->boid_z) {
3019         co[2] = col->boid_z;
3020         v0[2] = 0.0f;
3021       }
3022     }
3023
3024     /* re-apply acceleration to final location and velocity */
3025     madd_v3_v3v3fl(pa->state.co, co, v0, dt2);
3026     madd_v3_v3fl(pa->state.co, col->acc, 0.5f * dt2 * dt2);
3027     madd_v3_v3v3fl(pa->state.vel, v0, col->acc, dt2);
3028
3029     /* make sure particle stays on the right side of the surface */
3030     if (!through) {
3031       distance = collision_point_distance_with_normal(co, pce, -1.f, col, nor);
3032
3033       if (distance < col->radius + COLLISION_MIN_DISTANCE) {
3034         madd_v3_v3fl(co, nor, col->radius + COLLISION_MIN_DISTANCE - distance);
3035       }
3036
3037       dot = dot_v3v3(nor, v0);
3038       if (dot < 0.f) {
3039         madd_v3_v3fl(v0, nor, -dot);
3040       }
3041
3042       distance = collision_point_distance_with_normal(pa->state.co, pce, 1.f, col, nor);
3043
3044       if (distance < col->radius + COLLISION_MIN_DISTANCE) {
3045         madd_v3_v3fl(pa->state.co, nor, col->radius + COLLISION_MIN_DISTANCE - distance);
3046       }
3047
3048       dot = dot_v3v3(nor, pa->state.vel);
3049       if (dot < 0.f) {
3050         madd_v3_v3fl(pa->state.vel, nor, -dot);
3051       }
3052     }
3053
3054     /* add stickiness to surface */
3055     madd_v3_v3fl(pa->state.vel, pce->nor, -pd->pdef_stickness);
3056
3057     /* set coordinates for next iteration */
3058     copy_v3_v3(col->co1, co);
3059     copy_v3_v3(col->co2, pa->state.co);
3060
3061     copy_v3_v3(col->ve1, v0);
3062     copy_v3_v3(col->ve2, pa->state.vel);
3063
3064     col->f = f;
3065   }
3066
3067   /* if permeability random roll succeeded, disable collider for this sim step */
3068   if (through) {
3069     col->skip[col->skip_count++] = col->hit;
3070   }
3071
3072   return 1;
3073 }
3074 static void collision_fail(ParticleData *pa, ParticleCollision *col)
3075 {
3076   /* final chance to prevent total failure, so stick to the surface and hope for the best */
3077   collision_point_on_surface(col->co1, &col->pce, 1.f, col, pa->state.co);
3078
3079   copy_v3_v3(pa->state.vel, col->pce.vel);
3080   mul_v3_fl(pa->state.vel, col->inv_timestep);
3081
3082   /* printf("max iterations\n"); */
3083 }
3084
3085 /* Particle - Mesh collision detection and response
3086  * Features:
3087  * -friction and damping
3088  * -angular momentum <-> linear momentum
3089  * -high accuracy by re-applying particle acceleration after collision
3090  * -handles moving, rotating and deforming meshes
3091  * -uses Newton-Rhapson iteration to find the collisions
3092  * -handles spherical particles and (nearly) point like particles
3093  */
3094 static void collision_check(ParticleSimulationData *sim, int p, float dfra, float cfra)
3095 {
3096   ParticleSettings *part = sim->psys->part;
3097   ParticleData *pa = sim->psys->particles + p;
3098   ParticleCollision col;
3099   BVHTreeRayHit hit;
3100   int collision_count = 0;
3101
3102   float timestep = psys_get_timestep(sim);
3103
3104   memset(&col, 0, sizeof(ParticleCollision));
3105
3106   col.total_time = timestep * dfra;
3107   col.inv_total_time = 1.0f / col.total_time;
3108   col.inv_timestep = 1.0f / timestep;
3109
3110   col.cfra = cfra;
3111   col.old_cfra = sim->psys->cfra;
3112
3113   /* get acceleration (from gravity, forcefields etc. to be re-applied in collision response) */
3114   sub_v3_v3v3(col.acc, pa->state.vel, pa->prev_state.vel);
3115   mul_v3_fl(col.acc, 1.f / col.total_time);
3116
3117   /* set values for first iteration */
3118   copy_v3_v3(col.co1, pa->prev_state.co);
3119   copy_v3_v3(col.co2, pa->state.co);
3120   copy_v3_v3(col.ve1, pa->prev_state.vel);
3121   copy_v3_v3(col.ve2, pa->state.vel);
3122   col.f = 0.0f;
3123
3124   col.radius = ((part->flag & PART_SIZE_DEFL) || (part->phystype == PART_PHYS_BOIDS)) ?
3125                    pa->size :
3126                    COLLISION_MIN_RADIUS;
3127
3128   /* override for boids */
3129   if (part->phystype == PART_PHYS_BOIDS && part->boids->options & BOID_ALLOW_LAND) {
3130     col.boid = 1;
3131     col.boid_z = pa->state.co[2];
3132     col.skip[col.skip_count++] = pa->boid->ground;
3133   }
3134
3135   /* 10 iterations to catch multiple collisions */
3136   while (collision_count < PARTICLE_COLLISION_MAX_COLLISIONS) {
3137     if (collision_detect(pa, &col, &hit, sim->colliders)) {
3138
3139       collision_count++;
3140
3141       if (collision_count == PARTICLE_COLLISION_MAX_COLLISIONS) {
3142         collision_fail(pa, &col);
3143       }
3144       else if (collision_response(
3145                    sim, pa, &col, &hit, part->flag & PART_DIE_ON_COL, part->flag & PART_ROT_DYN) ==
3146                0) {
3147         return;
3148       }
3149     }
3150     else {
3151       return;
3152     }
3153   }
3154 }
3155 /************************************************/
3156 /*          Hair                                */
3157 /************************************************/
3158 /* check if path cache or children need updating and do it if needed */
3159 static void psys_update_path_cache(ParticleSimulationData *sim,
3160                                    float cfra,
3161                                    const bool use_render_params)
3162 {
3163   ParticleSystem *psys = sim->psys;
3164   ParticleSettings *part = psys->part;
3165   ParticleEditSettings *pset = &sim->scene->toolsettings->particle;
3166   int distr = 0, alloc = 0, skip = 0;
3167
3168   if ((psys->part->childtype &&
3169        psys->totchild != psys_get_tot_child(sim->scene, psys, use_render_params)) ||
3170       psys->recalc & ID_RECALC_PSYS_RESET) {
3171     alloc = 1;
3172   }
3173
3174   if (alloc || psys->recalc & ID_RECALC_PSYS_CHILD ||
3175       (psys->vgroup[PSYS_VG_DENSITY] && (sim->ob && sim->ob->mode & OB_MODE_WEIGHT_PAINT))) {
3176     distr = 1;
3177   }
3178
3179   if (distr) {
3180     if (alloc) {
3181       realloc_particles(sim, sim->psys->totpart);
3182     }
3183
3184     if (psys_get_tot_child(sim->scene, psys, use_render_params)) {
3185       /* don't generate children while computing the hair keys */
3186       if (!(psys->part->type == PART_HAIR) || (psys->flag & PSYS_HAIR_DONE)) {
3187         distribute_particles(sim, PART_FROM_CHILD);
3188
3189         if (part->childtype == PART_CHILD_FACES && part->parents != 0.0f) {
3190           psys_find_parents(sim, use_render_params);
3191         }
3192       }
3193     }
3194     else {
3195       psys_free_children(psys);
3196     }
3197   }
3198
3199   if ((part->type == PART_HAIR || psys->flag & PSYS_KEYED ||
3200        psys->pointcache->flag & PTCACHE_BAKED) == 0) {
3201     skip = 1; /* only hair, keyed and baked stuff can have paths */
3202   }
3203   else if (part->ren_as != PART_DRAW_PATH &&
3204            !(part->type == PART_HAIR && ELEM(part->ren_as, PART_DRAW_OB, PART_DRAW_GR))) {
3205     skip = 1; /* particle visualization must be set as path */
3206   }
3207   else if (DEG_get_mode(sim->depsgraph) != DAG_EVAL_RENDER) {
3208     if (part->draw_as != PART_DRAW_REND) {
3209       skip = 1; /* draw visualization */
3210     }
3211     else if (psys->pointcache->flag & PTCACHE_BAKING) {
3212       skip = 1; /* no need to cache paths while baking dynamics */
3213     }
3214     else if (psys_in_edit_mode(sim->depsgraph, psys)) {
3215       if ((pset->flag & PE_DRAW_PART) == 0) {
3216         skip = 1;
3217       }
3218       else if (part->childtype == 0 &&
3219                (psys->flag & PSYS_HAIR_DYNAMICS && psys->pointcache->flag & PTCACHE_BAKED) == 0) {
3220         skip = 1; /* in edit mode paths are needed for child particles and dynamic hair */
3221       }
3222     }
3223   }
3224
3225   if (!skip) {
3226     psys_cache_paths(sim, cfra, use_render_params);
3227
3228     /* for render, child particle paths are computed on the fly */
3229     if (part->childtype) {
3230       if (!psys->totchild) {
3231         skip = 1;
3232       }
3233       else if (psys->part->type == PART_HAIR && (psys->flag & PSYS_HAIR_DONE) == 0) {
3234         skip = 1;
3235       }
3236
3237       if (!skip) {
3238         psys_cache_child_paths(sim, cfra, 0, use_render_params);
3239       }
3240     }
3241   }
3242   else if (psys->pathcache) {
3243     psys_free_path_cache(psys, NULL);
3244   }
3245 }
3246
3247 static bool psys_hair_use_simulation(ParticleData *pa, float max_length)
3248 {
3249   /* Minimum segment length relative to average length.
3250    * Hairs with segments below this length will be excluded from the simulation,
3251    * because otherwise the solver will become unstable.
3252    * The hair system should always make sure the hair segments have reasonable length ratios,
3253    * but this can happen in old files when e.g. cutting hair.
3254    */
3255   const float min_length = 0.1f * max_length;
3256
3257   HairKey *key;
3258   int k;
3259
3260   if (pa->totkey < 2) {
3261     return false;
3262   }
3263
3264   for (k = 1, key = pa->hair + 1; k < pa->totkey; k++, key++) {
3265     float length = len_v3v3(key->co, (key - 1)->co);
3266     if (length < min_length) {
3267       return false;
3268     }
3269   }
3270
3271   return true;
3272 }
3273
3274 static MDeformVert *hair_set_pinning(MDeformVert *dvert, float weight)
3275 {
3276   if (dvert) {
3277     if (!dvert->totweight) {
3278       dvert->dw = MEM_callocN(sizeof(MDeformWeight), "deformWeight");
3279       dvert->totweight = 1;
3280     }
3281
3282     dvert->dw->weight = weight;
3283     dvert++;
3284   }
3285   return dvert;
3286 }
3287
3288 static void hair_create_input_mesh(ParticleSimulationData *sim,
3289                                    int totpoint,
3290                                    int totedge,
3291                                    Mesh **r_mesh,
3292                                    ClothHairData **r_hairdata)
3293 {
3294   ParticleSystem *psys = sim->psys;
3295   ParticleSettings *part = psys->part;
3296   Mesh *mesh;
3297   ClothHairData *hairdata;
3298   MVert *mvert;
3299   MEdge *medge;
3300   MDeformVert *dvert;
3301   HairKey *key;
3302   PARTICLE_P;
3303   int k, hair_index;
3304   float hairmat[4][4];
3305   float max_length;
3306   float hair_radius;
3307
3308   mesh = *r_mesh;
3309   if (!mesh) {
3310     *r_mesh = mesh = BKE_mesh_new_nomain(totpoint, totedge, 0, 0, 0);
3311     CustomData_add_layer(&mesh->vdata, CD_MDEFORMVERT, CD_CALLOC, NULL, mesh->totvert);
3312     BKE_mesh_update_customdata_pointers(mesh, false);
3313   }
3314   mvert = mesh->mvert;
3315   medge = mesh->medge;
3316   dvert = mesh->dvert;
3317
3318   hairdata = *r_hairdata;
3319   if (!hairdata) {
3320     *r_hairdata = hairdata = MEM_mallocN(sizeof(ClothHairData) * totpoint, "hair data");
3321   }
3322
3323   /* calculate maximum segment length */
3324   max_length = 0.0f;
3325   LOOP_PARTICLES
3326   {
3327     if (!(pa->flag & PARS_UNEXIST)) {
3328       for (k = 1, key = pa->hair + 1; k < pa->totkey; k++, key++) {
3329         float length = len_v3v3(key->co, (key - 1)->co);
3330         if (max_length < length) {
3331           max_length = length;
3332         }
3333       }
3334     }
3335   }
3336
3337   psys->clmd->sim_parms->vgroup_mass = 1;
3338
3339   /* XXX placeholder for more flexible future hair settings */
3340   hair_radius = part->size;
3341
3342   /* make vgroup for pin roots etc.. */
3343   hair_index = 1;
3344   LOOP_PARTICLES
3345   {
3346     if (!(pa->flag & PARS_UNEXIST)) {
3347       float root_mat[4][4];
3348       float bending_stiffness;
3349       bool use_hair;
3350
3351       pa->hair_index = hair_index;
3352       use_hair = psys_hair_use_simulation(pa, max_length);
3353
3354       psys_mat_hair_to_object(sim->ob, sim->psmd->mesh_final, psys->part->from, pa, hairmat);
3355       mul_m4_m4m4(root_mat, sim->ob->obmat, hairmat);
3356       normalize_m4(root_mat);
3357
3358       bending_stiffness = CLAMPIS(
3359           1.0f - part->bending_random * psys_frand(psys, p + 666), 0.0f, 1.0f);
3360
3361       for (k = 0, key = pa->hair; k < pa->totkey; k++, key++) {
3362         ClothHairData *hair;
3363         float *co, *co_next;
3364
3365         co = key->co;
3366         co_next = (key + 1)->co;
3367
3368         /* create fake root before actual root to resist bending */
3369         if (k == 0) {
3370           hair = &psys->clmd->hairdata[pa->hair_index - 1];
3371           copy_v3_v3(hair->loc, root_mat[3]);
3372           copy_m3_m4(hair->rot, root_mat);
3373
3374           hair->radius = hair_radius;
3375           hair->bending_stiffness = bending_stiffness;
3376
3377           add_v3_v3v3(mvert->co, co, co);
3378           sub_v3_v3(mvert->co, co_next);
3379           mul_m4_v3(hairmat, mvert->co);
3380
3381           medge->v1 = pa->hair_index - 1;
3382           medge->v2 = pa->hair_index;
3383
3384           dvert = hair_set_pinning(dvert, 1.0f);
3385
3386           mvert++;
3387           medge++;
3388         }
3389
3390         /* store root transform in cloth data */
3391         hair = &psys->clmd->hairdata[pa->hair_index + k];
3392         copy_v3_v3(hair->loc, root_mat[3]);
3393         copy_m3_m4(hair->rot, root_mat);
3394
3395         hair->radius = hair_radius;
3396         hair->bending_stiffness = bending_stiffness;
3397
3398         copy_v3_v3(mvert->co, co);
3399         mul_m4_v3(hairmat, mvert->co);
3400
3401         if (k) {
3402           medge->v1 = pa->hair_index + k - 1;
3403           medge->v2 = pa->hair_index + k;
3404         }
3405
3406         /* roots and disabled hairs should be 1.0, the rest can be anything from 0.0 to 1.0 */
3407         if (use_hair) {
3408           dvert = hair_set_pinning(dvert, key->weight);
3409         }
3410         else {
3411           dvert = hair_set_pinning(dvert, 1.0f);
3412         }
3413
3414         mvert++;
3415         if (k) {
3416           medge++;
3417         }
3418       }
3419
3420       hair_index += pa->totkey + 1;
3421     }
3422   }
3423 }
3424
3425 static void do_hair_dynamics(ParticleSimulationData *sim)
3426 {
3427   ParticleSystem *psys = sim->psys;
3428   PARTICLE_P;
3429   EffectorWeights *clmd_effweights;
3430   int totpoint;
3431   int totedge;
3432   float(*deformedVerts)[3];
3433   bool realloc_roots;
3434
3435   if (!psys->clmd) {
3436     psys->clmd = (ClothModifierData *)modifier_new(eModifierType_Cloth);
3437     psys->clmd->sim_parms->goalspring = 0.0f;
3438     psys->clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_RESIST_SPRING_COMPRESS;
3439     psys->clmd->coll_parms->flags &= ~CLOTH_COLLSETTINGS_FLAG_SELF;
3440   }
3441
3442   /* count simulated points */
3443   totpoint = 0;
3444   totedge = 0;
3445   LOOP_PARTICLES
3446   {
3447     if (!(pa->flag & PARS_UNEXIST)) {
3448       /* "out" dm contains all hairs */
3449       totedge += pa->totkey;
3450       totpoint += pa->totkey + 1; /* +1 for virtual root point */
3451     }
3452   }
3453
3454   /* whether hair root info array has to be reallocated */
3455   realloc_roots = false;
3456   if (psys->hair_in_mesh) {
3457     Mesh *mesh = psys->hair_in_mesh;
3458     if (totpoint != mesh->totvert || totedge != mesh->totedge) {
3459       BKE_id_free(NULL, mesh);
3460       psys->hair_in_mesh = NULL;
3461       realloc_roots = true;
3462     }
3463   }
3464
3465   if (!psys->hair_in_mesh || !psys->clmd->hairdata || realloc_roots) {
3466     if (psys->clmd->hairdata) {
3467       MEM_freeN(psys->clmd->hairdata);
3468       psys->clmd->hairdata = NULL;
3469     }
3470   }
3471
3472   hair_create_input_mesh(sim, totpoint, totedge, &psys->hair_in_mesh, &psys->clmd->hairdata);
3473
3474   if (psys->hair_out_mesh) {
3475     BKE_id_free(NULL, psys->hair_out_mesh);
3476   }
3477
3478   psys->clmd->point_cache = psys->pointcache;
3479   /* for hair sim we replace the internal cloth effector weights temporarily
3480    * to use the particle settings
3481    */
3482   clmd_effweights = psys->clmd->sim_parms->effector_weights;
3483   psys->clmd->sim_parms->effector_weights = psys->part->effector_weights;
3484
3485   BKE_id_copy_ex(NULL, &psys->hair_in_mesh->id, (ID **)&psys->hair_out_mesh, LIB_ID_COPY_LOCALIZE);
3486   deformedVerts = BKE_mesh_vertexCos_get(psys->hair_out_mesh, NULL);
3487   clothModifier_do(
3488       psys->clmd, sim->depsgraph, sim->scene, sim->ob, psys->hair_in_mesh, deformedVerts);
3489   BKE_mesh_apply_vert_coords(psys->hair_out_mesh, deformedVerts);
3490
3491   MEM_freeN(deformedVerts);
3492
3493   /* restore cloth effector weights */
3494   psys->clmd->sim_parms->effector_weights = clmd_effweights;
3495 }
3496 static void hair_step(ParticleSimulationData *sim, float cfra, const bool use_render_params)
3497 {
3498   ParticleSystem *psys = sim->psys;
3499   ParticleSettings *part = psys->part;
3500   PARTICLE_P;
3501   float disp = psys_get_current_display_percentage(psys, use_render_params);
3502
3503   LOOP_PARTICLES
3504   {
3505     pa->size = part->size;
3506     if (part->randsize > 0.0f) {
3507       pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);
3508     }
3509
3510     if (psys_frand(psys, p) > disp) {
3511       pa->flag |= PARS_NO_DISP;
3512     }
3513     else {
3514       pa->flag &= ~PARS_NO_DISP;
3515     }
3516   }
3517
3518   if (psys->recalc & ID_RECALC_PSYS_RESET) {
3519     /* need this for changing subsurf levels */
3520     psys_calc_dmcache(sim->ob, sim->psmd->mesh_final, sim->psmd->mesh_original, psys);
3521
3522     if (psys->clmd) {
3523       cloth_free_modifier(psys->clmd);
3524     }
3525   }
3526
3527   /* dynamics with cloth simulation, psys->particles can be NULL with 0 particles [#25519] */
3528   if (psys->part->type == PART_HAIR && psys->flag & PSYS_HAIR_DYNAMICS && psys->particles) {
3529     do_hair_dynamics(sim);
3530   }
3531
3532   /* following lines were removed r29079 but cause bug [#22811], see report for details */
3533   psys_update_effectors(sim);
3534   psys_update_path_cache(sim, cfra, use_render_params);
3535
3536   psys->flag |= PSYS_HAIR_UPDATED;
3537 }
3538
3539 static void save_hair(ParticleSimulationData *sim, float UNUSED(cfra))
3540 {
3541   Object *ob = sim->ob;
3542   ParticleSystem *psys = sim->psys;
3543   HairKey *key, *root;
3544   PARTICLE_P;
3545
3546   invert_m4_m4(ob->imat, ob->obmat);
3547
3548   psys->lattice_deform_data = psys_create_lattice_deform_data(sim);
3549
3550   if (psys->totpart == 0) {
3551     return;
3552   }
3553
3554   /* save new keys for elements if needed */
3555   LOOP_PARTICLES
3556   {
3557     /* first time alloc */
3558     if (pa->totkey == 0 || pa->hair == NULL) {
3559       pa->hair = MEM_callocN((psys->part->hair_step + 1) * sizeof(HairKey), "HairKeys");
3560       pa->totkey = 0;
3561     }
3562
3563     key = root = pa->hair;
3564     key += pa->totkey;
3565
3566     /* convert from global to geometry space */
3567     copy_v3_v3(key->co, pa->state.co);
3568     mul_m4_v3(ob->imat, key->co);
3569
3570     if (pa->totkey) {
3571       sub_v3_v3(key->co, root->co);
3572       psys_vec_rot_to_face(sim->psmd->mesh_final, pa, key->co);
3573     }
3574
3575     key->time = pa->state.time;
3576
3577     key->weight = 1.0f - key->time / 100.0f;
3578
3579     pa->totkey++;
3580
3581     /* root is always in the origin of hair space so we set it to be so after the last key is saved*/
3582     if (pa->totkey == psys->part->hair_step + 1) {
3583       zero_v3(root->co);
3584     }
3585   }
3586 }
3587
3588 /* Code for an adaptive time step based on the Courant-Friedrichs-Lewy
3589  * condition. */
3590 static const float MIN_TIMESTEP = 1.0f / 101.0f;
3591 /* Tolerance of 1.5 means the last subframe neither favors growing nor
3592  * shrinking (e.g if it were 1.3, the last subframe would tend to be too
3593  * small). */
3594 static const float TIMESTEP_EXPANSION_FACTOR = 0.1f;
3595 static const float TIMESTEP_EXPANSION_TOLERANCE = 1.5f;
3596
3597 /* Calculate the speed of the particle relative to the local scale of the
3598  * simulation. This should be called once per particle during a simulation
3599  * step, after the velocity has been updated. element_size defines the scale of
3600  * the simulation, and is typically the distance to neighboring particles. */
3601 static void update_courant_num(
3602     ParticleSimulationData *sim, ParticleData *pa, float dtime, SPHData *sphdata, SpinLock *spin)
3603 {
3604   float relative_vel[3];
3605
3606   sub_v3_v3v3(relative_vel, pa->prev_state.vel, sphdata->flow);