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 <math.h>
30 #include <stdlib.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_curve_types.h"
38 #include "DNA_listBase.h"
39 #include "DNA_mesh_types.h"
40 #include "DNA_meshdata_types.h"
41 #include "DNA_modifier_types.h"
42 #include "DNA_object_force_types.h"
43 #include "DNA_object_types.h"
44 #include "DNA_particle_types.h"
45 #include "DNA_scene_types.h"
46 #include "DNA_texture_types.h"
47
48 #include "BLI_blenlib.h"
49 #include "BLI_edgehash.h"
50 #include "BLI_kdopbvh.h"
51 #include "BLI_kdtree.h"
52 #include "BLI_linklist.h"
53 #include "BLI_math.h"
54 #include "BLI_rand.h"
55 #include "BLI_string_utils.h"
56 #include "BLI_task.h"
57 #include "BLI_threads.h"
58 #include "BLI_utildefines.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_lib_id.h"
66 #include "BKE_lib_query.h"
67 #include "BKE_particle.h"
68
69 #include "BKE_bvhutils.h"
70 #include "BKE_cloth.h"
71 #include "BKE_collection.h"
72 #include "BKE_lattice.h"
73 #include "BKE_material.h"
74 #include "BKE_mesh.h"
75 #include "BKE_modifier.h"
76 #include "BKE_object.h"
77 #include "BKE_pointcache.h"
78 #include "BKE_scene.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_FLUID
90 # include "DNA_fluid_types.h"
91 # include "manta_fluid_API.h"
92 #endif // WITH_FLUID
93
94 static ThreadRWMutex psys_bvhtree_rwlock = BLI_RWLOCK_INITIALIZER;
95
96 /************************************************/
97 /* Reacting to system events */
98 /************************************************/
99
particles_are_dynamic(ParticleSystem * psys)100 static int particles_are_dynamic(ParticleSystem *psys)
101 {
102 if (psys->pointcache->flag & PTCACHE_BAKED) {
103 return 0;
104 }
105
106 if (psys->part->type == PART_HAIR) {
107 return psys->flag & PSYS_HAIR_DYNAMICS;
108 }
109
110 return ELEM(psys->part->phystype, PART_PHYS_NEWTON, PART_PHYS_BOIDS, PART_PHYS_FLUID);
111 }
112
psys_get_current_display_percentage(ParticleSystem * psys,const bool use_render_params)113 float psys_get_current_display_percentage(ParticleSystem *psys, const bool use_render_params)
114 {
115 ParticleSettings *part = psys->part;
116
117 if ((use_render_params &&
118 !particles_are_dynamic(psys)) || /* non-dynamic particles can be rendered fully */
119 (part->child_nbr && part->childtype) || /* display percentage applies to children */
120 (psys->pointcache->flag & PTCACHE_BAKING)) /* baking is always done with full amount */
121 {
122 return 1.0f;
123 }
124
125 return psys->part->disp / 100.0f;
126 }
127
tot_particles(ParticleSystem * psys,PTCacheID * pid)128 static int tot_particles(ParticleSystem *psys, PTCacheID *pid)
129 {
130 if (pid && psys->pointcache->flag & PTCACHE_EXTERNAL) {
131 return pid->cache->totpoint;
132 }
133 if (psys->part->distr == PART_DISTR_GRID && psys->part->from != PART_FROM_VERT) {
134 return psys->part->grid_res * psys->part->grid_res * psys->part->grid_res - psys->totunexist;
135 }
136
137 return psys->part->totpart - psys->totunexist;
138 }
139
psys_reset(ParticleSystem * psys,int mode)140 void psys_reset(ParticleSystem *psys, int mode)
141 {
142 PARTICLE_P;
143
144 if (ELEM(mode, PSYS_RESET_ALL, PSYS_RESET_DEPSGRAPH)) {
145 if (mode == PSYS_RESET_ALL || !(psys->flag & PSYS_EDITED)) {
146 /* don't free if not absolutely necessary */
147 if (psys->totpart != tot_particles(psys, NULL)) {
148 psys_free_particles(psys);
149 psys->totpart = 0;
150 }
151
152 psys->totkeyed = 0;
153 psys->flag &= ~(PSYS_HAIR_DONE | PSYS_KEYED);
154
155 if (psys->edit && psys->free_edit) {
156 psys->free_edit(psys->edit);
157 psys->edit = NULL;
158 psys->free_edit = NULL;
159 }
160 }
161 }
162 else if (mode == PSYS_RESET_CACHE_MISS) {
163 /* set all particles to be skipped */
164 LOOP_PARTICLES
165 {
166 pa->flag |= PARS_NO_DISP;
167 }
168 }
169
170 /* reset children */
171 if (psys->child) {
172 MEM_freeN(psys->child);
173 psys->child = NULL;
174 }
175
176 psys->totchild = 0;
177
178 /* reset path cache */
179 psys_free_path_cache(psys, psys->edit);
180
181 /* reset point cache */
182 BKE_ptcache_invalidate(psys->pointcache);
183
184 if (psys->fluid_springs) {
185 MEM_freeN(psys->fluid_springs);
186 psys->fluid_springs = NULL;
187 }
188
189 psys->tot_fluidsprings = psys->alloc_fluidsprings = 0;
190 }
191
psys_unique_name(Object * object,ParticleSystem * psys,const char * defname)192 void psys_unique_name(Object *object, ParticleSystem *psys, const char *defname)
193 {
194 BLI_uniquename(&object->particlesystem,
195 psys,
196 defname,
197 '.',
198 offsetof(ParticleSystem, name),
199 sizeof(psys->name));
200 }
201
realloc_particles(ParticleSimulationData * sim,int new_totpart)202 static void realloc_particles(ParticleSimulationData *sim, int new_totpart)
203 {
204 ParticleSystem *psys = sim->psys;
205 ParticleSettings *part = psys->part;
206 ParticleData *newpars = NULL;
207 BoidParticle *newboids = NULL;
208 PARTICLE_P;
209 int totpart, totsaved = 0;
210
211 if (new_totpart < 0) {
212 if ((part->distr == PART_DISTR_GRID) && (part->from != PART_FROM_VERT)) {
213 totpart = part->grid_res;
214 totpart *= totpart * totpart;
215 }
216 else {
217 totpart = part->totpart;
218 }
219 }
220 else {
221 totpart = new_totpart;
222 }
223
224 if (totpart != psys->totpart) {
225 if (psys->edit && psys->free_edit) {
226 psys->free_edit(psys->edit);
227 psys->edit = NULL;
228 psys->free_edit = NULL;
229 }
230
231 if (totpart) {
232 newpars = MEM_callocN(totpart * sizeof(ParticleData), "particles");
233 if (newpars == NULL) {
234 return;
235 }
236
237 if (psys->part->phystype == PART_PHYS_BOIDS) {
238 newboids = MEM_callocN(totpart * sizeof(BoidParticle), "boid particles");
239
240 if (newboids == NULL) {
241 /* allocation error! */
242 if (newpars) {
243 MEM_freeN(newpars);
244 }
245 return;
246 }
247 }
248 }
249
250 if (psys->particles) {
251 totsaved = MIN2(psys->totpart, totpart);
252 /*save old pars*/
253 if (totsaved) {
254 memcpy(newpars, psys->particles, totsaved * sizeof(ParticleData));
255
256 if (psys->particles->boid) {
257 memcpy(newboids, psys->particles->boid, totsaved * sizeof(BoidParticle));
258 }
259 }
260
261 if (psys->particles->keys) {
262 MEM_freeN(psys->particles->keys);
263 }
264
265 if (psys->particles->boid) {
266 MEM_freeN(psys->particles->boid);
267 }
268
269 for (p = 0, pa = newpars; p < totsaved; p++, pa++) {
270 if (pa->keys) {
271 pa->keys = NULL;
272 pa->totkey = 0;
273 }
274 }
275
276 for (p = totsaved, pa = psys->particles + totsaved; p < psys->totpart; p++, pa++) {
277 if (pa->hair) {
278 MEM_freeN(pa->hair);
279 }
280 }
281
282 MEM_freeN(psys->particles);
283 psys_free_pdd(psys);
284 }
285
286 psys->particles = newpars;
287 psys->totpart = totpart;
288
289 if (newboids) {
290 LOOP_PARTICLES
291 {
292 pa->boid = newboids++;
293 }
294 }
295 }
296
297 if (psys->child) {
298 MEM_freeN(psys->child);
299 psys->child = NULL;
300 psys->totchild = 0;
301 }
302 }
303
psys_get_child_number(Scene * scene,ParticleSystem * psys,const bool use_render_params)304 int psys_get_child_number(Scene *scene, ParticleSystem *psys, const bool use_render_params)
305 {
306 int nbr;
307
308 if (!psys->part->childtype) {
309 return 0;
310 }
311
312 if (use_render_params) {
313 nbr = psys->part->ren_child_nbr;
314 }
315 else {
316 nbr = psys->part->child_nbr;
317 }
318
319 return get_render_child_particle_number(&scene->r, nbr, use_render_params);
320 }
321
psys_get_tot_child(Scene * scene,ParticleSystem * psys,const bool use_render_params)322 int psys_get_tot_child(Scene *scene, ParticleSystem *psys, const bool use_render_params)
323 {
324 return psys->totpart * psys_get_child_number(scene, psys, use_render_params);
325 }
326
327 /************************************************/
328 /* Distribution */
329 /************************************************/
330
psys_calc_dmcache(Object * ob,Mesh * mesh_final,Mesh * mesh_original,ParticleSystem * psys)331 void psys_calc_dmcache(Object *ob, Mesh *mesh_final, Mesh *mesh_original, ParticleSystem *psys)
332 {
333 /* use for building derived mesh mapping info:
334 *
335 * node: the allocated links - total derived mesh element count
336 * nodearray: the array of nodes aligned with the base mesh's elements, so
337 * each original elements can reference its derived elements
338 */
339 Mesh *me = (Mesh *)ob->data;
340 bool use_modifier_stack = psys->part->use_modifier_stack;
341 PARTICLE_P;
342
343 /* CACHE LOCATIONS */
344 if (!mesh_final->runtime.deformed_only) {
345 /* Will use later to speed up subsurf/evaluated mesh. */
346 LinkNode *node, *nodedmelem, **nodearray;
347 int totdmelem, totelem, i, *origindex, *origindex_poly = NULL;
348
349 if (psys->part->from == PART_FROM_VERT) {
350 totdmelem = mesh_final->totvert;
351
352 if (use_modifier_stack) {
353 totelem = totdmelem;
354 origindex = NULL;
355 }
356 else {
357 totelem = me->totvert;
358 origindex = CustomData_get_layer(&mesh_final->vdata, CD_ORIGINDEX);
359 }
360 }
361 else { /* FROM_FACE/FROM_VOLUME */
362 totdmelem = mesh_final->totface;
363
364 if (use_modifier_stack) {
365 totelem = totdmelem;
366 origindex = NULL;
367 origindex_poly = NULL;
368 }
369 else {
370 totelem = mesh_original->totface;
371 origindex = CustomData_get_layer(&mesh_final->fdata, CD_ORIGINDEX);
372
373 /* for face lookups we need the poly origindex too */
374 origindex_poly = CustomData_get_layer(&mesh_final->pdata, CD_ORIGINDEX);
375 if (origindex_poly == NULL) {
376 origindex = NULL;
377 }
378 }
379 }
380
381 nodedmelem = MEM_callocN(sizeof(LinkNode) * totdmelem, "psys node elems");
382 nodearray = MEM_callocN(sizeof(LinkNode *) * totelem, "psys node array");
383
384 for (i = 0, node = nodedmelem; i < totdmelem; i++, node++) {
385 int origindex_final;
386 node->link = POINTER_FROM_INT(i);
387
388 /* may be vertex or face origindex */
389 if (use_modifier_stack) {
390 origindex_final = i;
391 }
392 else {
393 origindex_final = origindex ? origindex[i] : ORIGINDEX_NONE;
394
395 /* if we have a poly source, do an index lookup */
396 if (origindex_poly && origindex_final != ORIGINDEX_NONE) {
397 origindex_final = origindex_poly[origindex_final];
398 }
399 }
400
401 if (origindex_final != ORIGINDEX_NONE && origindex_final < totelem) {
402 if (nodearray[origindex_final]) {
403 /* prepend */
404 node->next = nodearray[origindex_final];
405 nodearray[origindex_final] = node;
406 }
407 else {
408 nodearray[origindex_final] = node;
409 }
410 }
411 }
412
413 /* cache the verts/faces! */
414 LOOP_PARTICLES
415 {
416 if (pa->num < 0) {
417 pa->num_dmcache = DMCACHE_NOTFOUND;
418 continue;
419 }
420
421 if (use_modifier_stack) {
422 if (pa->num < totelem) {
423 pa->num_dmcache = DMCACHE_ISCHILD;
424 }
425 else {
426 pa->num_dmcache = DMCACHE_NOTFOUND;
427 }
428 }
429 else {
430 if (psys->part->from == PART_FROM_VERT) {
431 if (pa->num < totelem && nodearray[pa->num]) {
432 pa->num_dmcache = POINTER_AS_INT(nodearray[pa->num]->link);
433 }
434 else {
435 pa->num_dmcache = DMCACHE_NOTFOUND;
436 }
437 }
438 else { /* FROM_FACE/FROM_VOLUME */
439 pa->num_dmcache = psys_particle_dm_face_lookup(
440 mesh_final, mesh_original, pa->num, pa->fuv, nodearray);
441 }
442 }
443 }
444
445 MEM_freeN(nodearray);
446 MEM_freeN(nodedmelem);
447 }
448 else {
449 /* TODO PARTICLE, make the following line unnecessary, each function
450 * should know to use the num or num_dmcache, set the num_dmcache to
451 * an invalid value, just in case */
452
453 LOOP_PARTICLES
454 {
455 pa->num_dmcache = DMCACHE_NOTFOUND;
456 }
457 }
458 }
459
460 /* threaded child particle distribution and path caching */
psys_thread_context_init(ParticleThreadContext * ctx,ParticleSimulationData * sim)461 void psys_thread_context_init(ParticleThreadContext *ctx, ParticleSimulationData *sim)
462 {
463 memset(ctx, 0, sizeof(ParticleThreadContext));
464 ctx->sim = *sim;
465 ctx->mesh = ctx->sim.psmd->mesh_final;
466 ctx->ma = BKE_object_material_get(sim->ob, sim->psys->part->omat);
467 }
468
psys_tasks_create(ParticleThreadContext * ctx,int startpart,int endpart,ParticleTask ** r_tasks,int * r_numtasks)469 void psys_tasks_create(ParticleThreadContext *ctx,
470 int startpart,
471 int endpart,
472 ParticleTask **r_tasks,
473 int *r_numtasks)
474 {
475 ParticleTask *tasks;
476 int numtasks = min_ii(BLI_system_thread_count() * 4, endpart - startpart);
477 float particles_per_task = (float)(endpart - startpart) / (float)numtasks, p, pnext;
478 int i;
479
480 tasks = MEM_callocN(sizeof(ParticleTask) * numtasks, "ParticleThread");
481 *r_numtasks = numtasks;
482 *r_tasks = tasks;
483
484 p = (float)startpart;
485 for (i = 0; i < numtasks; i++, p = pnext) {
486 pnext = p + particles_per_task;
487
488 tasks[i].ctx = ctx;
489 tasks[i].begin = (int)p;
490 tasks[i].end = min_ii((int)pnext, endpart);
491 }
492 }
493
psys_tasks_free(ParticleTask * tasks,int numtasks)494 void psys_tasks_free(ParticleTask *tasks, int numtasks)
495 {
496 int i;
497
498 /* threads */
499 for (i = 0; i < numtasks; i++) {
500 if (tasks[i].rng) {
501 BLI_rng_free(tasks[i].rng);
502 }
503 if (tasks[i].rng_path) {
504 BLI_rng_free(tasks[i].rng_path);
505 }
506 }
507
508 MEM_freeN(tasks);
509 }
510
psys_thread_context_free(ParticleThreadContext * ctx)511 void psys_thread_context_free(ParticleThreadContext *ctx)
512 {
513 /* path caching */
514 if (ctx->vg_length) {
515 MEM_freeN(ctx->vg_length);
516 }
517 if (ctx->vg_clump) {
518 MEM_freeN(ctx->vg_clump);
519 }
520 if (ctx->vg_kink) {
521 MEM_freeN(ctx->vg_kink);
522 }
523 if (ctx->vg_rough1) {
524 MEM_freeN(ctx->vg_rough1);
525 }
526 if (ctx->vg_rough2) {
527 MEM_freeN(ctx->vg_rough2);
528 }
529 if (ctx->vg_roughe) {
530 MEM_freeN(ctx->vg_roughe);
531 }
532 if (ctx->vg_twist) {
533 MEM_freeN(ctx->vg_twist);
534 }
535
536 if (ctx->sim.psys->lattice_deform_data) {
537 BKE_lattice_deform_data_destroy(ctx->sim.psys->lattice_deform_data);
538 ctx->sim.psys->lattice_deform_data = NULL;
539 }
540
541 /* distribution */
542 if (ctx->jit) {
543 MEM_freeN(ctx->jit);
544 }
545 if (ctx->jitoff) {
546 MEM_freeN(ctx->jitoff);
547 }
548 if (ctx->weight) {
549 MEM_freeN(ctx->weight);
550 }
551 if (ctx->index) {
552 MEM_freeN(ctx->index);
553 }
554 if (ctx->seams) {
555 MEM_freeN(ctx->seams);
556 }
557 // if (ctx->vertpart) MEM_freeN(ctx->vertpart);
558 BLI_kdtree_3d_free(ctx->tree);
559
560 if (ctx->clumpcurve != NULL) {
561 BKE_curvemapping_free(ctx->clumpcurve);
562 }
563 if (ctx->roughcurve != NULL) {
564 BKE_curvemapping_free(ctx->roughcurve);
565 }
566 if (ctx->twistcurve != NULL) {
567 BKE_curvemapping_free(ctx->twistcurve);
568 }
569 }
570
init_particle_texture(ParticleSimulationData * sim,ParticleData * pa,int p)571 static void init_particle_texture(ParticleSimulationData *sim, ParticleData *pa, int p)
572 {
573 ParticleSystem *psys = sim->psys;
574 ParticleSettings *part = psys->part;
575 ParticleTexture ptex;
576
577 psys_get_texture(sim, pa, &ptex, PAMAP_INIT, 0.f);
578
579 switch (part->type) {
580 case PART_EMITTER:
581 if (ptex.exist < psys_frand(psys, p + 125)) {
582 pa->flag |= PARS_UNEXIST;
583 }
584 pa->time = part->sta + (part->end - part->sta) * ptex.time;
585 break;
586 case PART_HAIR:
587 if (ptex.exist < psys_frand(psys, p + 125)) {
588 pa->flag |= PARS_UNEXIST;
589 }
590 pa->time = 0.f;
591 break;
592 }
593 }
594
595 /* set particle parameters that don't change during particle's life */
init_particle(ParticleSimulationData * sim,ParticleData * pa)596 void init_particle(ParticleSimulationData *sim, ParticleData *pa)
597 {
598 ParticleSettings *part = sim->psys->part;
599 float birth_time = (float)(pa - sim->psys->particles) / (float)sim->psys->totpart;
600
601 pa->flag &= ~PARS_UNEXIST;
602 pa->time = part->sta + (part->end - part->sta) * birth_time;
603
604 pa->hair_index = 0;
605 /* we can't reset to -1 anymore since we've figured out correct index in distribute_particles */
606 /* usage other than straight after distribute has to handle this index by itself - jahka*/
607 // pa->num_dmcache = DMCACHE_NOTFOUND; /* assume we don't have a derived mesh face */
608 }
609
initialize_all_particles(ParticleSimulationData * sim)610 static void initialize_all_particles(ParticleSimulationData *sim)
611 {
612 ParticleSystem *psys = sim->psys;
613 ParticleSettings *part = psys->part;
614 /* Grid distributionsets UNEXIST flag, need to take care of
615 * it here because later this flag is being reset.
616 *
617 * We can't do it for any distribution, because it'll then
618 * conflict with texture influence, which does not free
619 * unexisting particles and only sets flag.
620 *
621 * It's not so bad, because only grid distribution sets
622 * UNEXIST flag.
623 */
624 const bool emit_from_volume_grid = (part->distr == PART_DISTR_GRID) &&
625 (!ELEM(part->from, PART_FROM_VERT, PART_FROM_CHILD));
626 PARTICLE_P;
627 LOOP_PARTICLES
628 {
629 if (!(emit_from_volume_grid && (pa->flag & PARS_UNEXIST) != 0)) {
630 init_particle(sim, pa);
631 }
632 }
633 }
634
free_unexisting_particles(ParticleSimulationData * sim)635 static void free_unexisting_particles(ParticleSimulationData *sim)
636 {
637 ParticleSystem *psys = sim->psys;
638 PARTICLE_P;
639
640 psys->totunexist = 0;
641
642 LOOP_PARTICLES
643 {
644 if (pa->flag & PARS_UNEXIST) {
645 psys->totunexist++;
646 }
647 }
648
649 if (psys->totpart && psys->totunexist == psys->totpart) {
650 if (psys->particles->boid) {
651 MEM_freeN(psys->particles->boid);
652 }
653
654 MEM_freeN(psys->particles);
655 psys->particles = NULL;
656 psys->totpart = psys->totunexist = 0;
657 }
658
659 if (psys->totunexist) {
660 int newtotpart = psys->totpart - psys->totunexist;
661 ParticleData *npa, *newpars;
662
663 npa = newpars = MEM_callocN(newtotpart * sizeof(ParticleData), "particles");
664
665 for (p = 0, pa = psys->particles; p < newtotpart; p++, pa++, npa++) {
666 while (pa->flag & PARS_UNEXIST) {
667 pa++;
668 }
669
670 memcpy(npa, pa, sizeof(ParticleData));
671 }
672
673 if (psys->particles->boid) {
674 MEM_freeN(psys->particles->boid);
675 }
676 MEM_freeN(psys->particles);
677 psys->particles = newpars;
678 psys->totpart -= psys->totunexist;
679
680 if (psys->particles->boid) {
681 BoidParticle *newboids = MEM_callocN(psys->totpart * sizeof(BoidParticle), "boid particles");
682
683 LOOP_PARTICLES
684 {
685 pa->boid = newboids++;
686 }
687 }
688 }
689 }
690
get_angular_velocity_vector(short avemode,ParticleKey * state,float vec[3])691 static void get_angular_velocity_vector(short avemode, ParticleKey *state, float vec[3])
692 {
693 switch (avemode) {
694 case PART_AVE_VELOCITY:
695 copy_v3_v3(vec, state->vel);
696 break;
697 case PART_AVE_HORIZONTAL: {
698 float zvec[3];
699 zvec[0] = zvec[1] = 0;
700 zvec[2] = 1.f;
701 cross_v3_v3v3(vec, state->vel, zvec);
702 break;
703 }
704 case PART_AVE_VERTICAL: {
705 float zvec[3], temp[3];
706 zvec[0] = zvec[1] = 0;
707 zvec[2] = 1.f;
708 cross_v3_v3v3(temp, state->vel, zvec);
709 cross_v3_v3v3(vec, temp, state->vel);
710 break;
711 }
712 case PART_AVE_GLOBAL_X:
713 vec[0] = 1.f;
714 vec[1] = vec[2] = 0;
715 break;
716 case PART_AVE_GLOBAL_Y:
717 vec[1] = 1.f;
718 vec[0] = vec[2] = 0;
719 break;
720 case PART_AVE_GLOBAL_Z:
721 vec[2] = 1.f;
722 vec[0] = vec[1] = 0;
723 break;
724 }
725 }
726
psys_get_birth_coords(ParticleSimulationData * sim,ParticleData * pa,ParticleKey * state,float dtime,float cfra)727 void psys_get_birth_coords(
728 ParticleSimulationData *sim, ParticleData *pa, ParticleKey *state, float dtime, float cfra)
729 {
730 Object *ob = sim->ob;
731 ParticleSystem *psys = sim->psys;
732 ParticleSettings *part = psys->part;
733 ParticleTexture ptex;
734 float fac, phasefac, nor[3] = {0, 0, 0}, loc[3], vel[3] = {0.0, 0.0, 0.0}, rot[4], q2[4];
735 float r_vel[3], r_ave[3], r_rot[4], vec[3], p_vel[3] = {0.0, 0.0, 0.0};
736 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},
737 rot_vec[3] = {0.0, 0.0, 0.0};
738 float q_phase[4];
739
740 const bool use_boids = ((part->phystype == PART_PHYS_BOIDS) && (pa->boid != NULL));
741 const bool use_tangents = ((use_boids == false) &&
742 ((part->tanfac != 0.0f) || (part->rotmode == PART_ROT_NOR_TAN)));
743
744 int p = pa - psys->particles;
745
746 /* get birth location from object */
747 if (use_tangents) {
748 psys_particle_on_emitter(sim->psmd,
749 part->from,
750 pa->num,
751 pa->num_dmcache,
752 pa->fuv,
753 pa->foffset,
754 loc,
755 nor,
756 utan,
757 vtan,
758 0);
759 }
760 else {
761 psys_particle_on_emitter(
762 sim->psmd, part->from, pa->num, pa->num_dmcache, pa->fuv, pa->foffset, loc, nor, 0, 0, 0);
763 }
764
765 /* get possible textural influence */
766 psys_get_texture(sim, pa, &ptex, PAMAP_IVEL, cfra);
767
768 /* particles live in global space so */
769 /* let's convert: */
770 /* -location */
771 mul_m4_v3(ob->obmat, loc);
772
773 /* -normal */
774 mul_mat3_m4_v3(ob->obmat, nor);
775 normalize_v3(nor);
776
777 /* -tangent */
778 if (use_tangents) {
779 #if 0
780 float phase = vg_rot ?
781 2.0f *
782 (psys_particle_value_from_verts(sim->psmd->dm, part->from, pa, vg_rot) -
783 0.5f) :
784 0.0f;
785 #else
786 float phase = 0.0f;
787 #endif
788 mul_v3_fl(vtan, -cosf((float)M_PI * (part->tanphase + phase)));
789 fac = -sinf((float)M_PI * (part->tanphase + phase));
790 madd_v3_v3fl(vtan, utan, fac);
791
792 mul_mat3_m4_v3(ob->obmat, vtan);
793
794 copy_v3_v3(utan, nor);
795 mul_v3_fl(utan, dot_v3v3(vtan, nor));
796 sub_v3_v3(vtan, utan);
797
798 normalize_v3(vtan);
799 }
800
801 /* -velocity (boids need this even if there's no random velocity) */
802 if (part->randfac != 0.0f || (part->phystype == PART_PHYS_BOIDS && pa->boid)) {
803 r_vel[0] = 2.0f * (psys_frand(psys, p + 10) - 0.5f);
804 r_vel[1] = 2.0f * (psys_frand(psys, p + 11) - 0.5f);
805 r_vel[2] = 2.0f * (psys_frand(psys, p + 12) - 0.5f);
806
807 mul_mat3_m4_v3(ob->obmat, r_vel);
808 normalize_v3(r_vel);
809 }
810
811 /* -angular velocity */
812 if (part->avemode == PART_AVE_RAND) {
813 r_ave[0] = 2.0f * (psys_frand(psys, p + 13) - 0.5f);
814 r_ave[1] = 2.0f * (psys_frand(psys, p + 14) - 0.5f);
815 r_ave[2] = 2.0f * (psys_frand(psys, p + 15) - 0.5f);
816
817 mul_mat3_m4_v3(ob->obmat, r_ave);
818 normalize_v3(r_ave);
819 }
820
821 /* -rotation */
822 if (part->randrotfac != 0.0f) {
823 r_rot[0] = 2.0f * (psys_frand(psys, p + 16) - 0.5f);
824 r_rot[1] = 2.0f * (psys_frand(psys, p + 17) - 0.5f);
825 r_rot[2] = 2.0f * (psys_frand(psys, p + 18) - 0.5f);
826 r_rot[3] = 2.0f * (psys_frand(psys, p + 19) - 0.5f);
827 normalize_qt(r_rot);
828
829 mat4_to_quat(rot, ob->obmat);
830 mul_qt_qtqt(r_rot, r_rot, rot);
831 }
832
833 if (use_boids) {
834 float dvec[3], q[4], mat[3][3];
835
836 copy_v3_v3(state->co, loc);
837
838 /* boids don't get any initial velocity */
839 zero_v3(state->vel);
840
841 /* boids store direction in ave */
842 if (fabsf(nor[2]) == 1.0f) {
843 sub_v3_v3v3(state->ave, loc, ob->obmat[3]);
844 normalize_v3(state->ave);
845 }
846 else {
847 copy_v3_v3(state->ave, nor);
848 }
849
850 /* calculate rotation matrix */
851 project_v3_v3v3(dvec, r_vel, state->ave);
852 sub_v3_v3v3(mat[0], state->ave, dvec);
853 normalize_v3(mat[0]);
854 negate_v3_v3(mat[2], r_vel);
855 normalize_v3(mat[2]);
856 cross_v3_v3v3(mat[1], mat[2], mat[0]);
857
858 /* apply rotation */
859 mat3_to_quat_is_ok(q, mat);
860 copy_qt_qt(state->rot, q);
861 }
862 else {
863 /* conversion done so now we apply new: */
864 /* -velocity from: */
865
866 /* *reactions */
867 if (dtime > 0.f) {
868 sub_v3_v3v3(vel, pa->state.vel, pa->prev_state.vel);
869 }
870
871 /* *emitter velocity */
872 if (dtime != 0.f && part->obfac != 0.f) {
873 sub_v3_v3v3(vel, loc, state->co);
874 mul_v3_fl(vel, part->obfac / dtime);
875 }
876
877 /* *emitter normal */
878 if (part->normfac != 0.f) {
879 madd_v3_v3fl(vel, nor, part->normfac);
880 }
881
882 /* *emitter tangent */
883 if (sim->psmd && part->tanfac != 0.f) {
884 madd_v3_v3fl(vel, vtan, part->tanfac);
885 }
886
887 /* *emitter object orientation */
888 if (part->ob_vel[0] != 0.f) {
889 normalize_v3_v3(vec, ob->obmat[0]);
890 madd_v3_v3fl(vel, vec, part->ob_vel[0]);
891 }
892 if (part->ob_vel[1] != 0.f) {
893 normalize_v3_v3(vec, ob->obmat[1]);
894 madd_v3_v3fl(vel, vec, part->ob_vel[1]);
895 }
896 if (part->ob_vel[2] != 0.f) {
897 normalize_v3_v3(vec, ob->obmat[2]);
898 madd_v3_v3fl(vel, vec, part->ob_vel[2]);
899 }
900
901 /* *texture */
902 /* TODO */
903
904 /* *random */
905 if (part->randfac != 0.f) {
906 madd_v3_v3fl(vel, r_vel, part->randfac);
907 }
908
909 /* *particle */
910 if (part->partfac != 0.f) {
911 madd_v3_v3fl(vel, p_vel, part->partfac);
912 }
913
914 mul_v3_v3fl(state->vel, vel, ptex.ivel);
915
916 /* -location from emitter */
917 copy_v3_v3(state->co, loc);
918
919 /* -rotation */
920 unit_qt(state->rot);
921
922 if (part->rotmode) {
923 bool use_global_space;
924
925 /* create vector into which rotation is aligned */
926 switch (part->rotmode) {
927 case PART_ROT_NOR:
928 case PART_ROT_NOR_TAN:
929 copy_v3_v3(rot_vec, nor);
930 use_global_space = false;
931 break;
932 case PART_ROT_VEL:
933 copy_v3_v3(rot_vec, vel);
934 use_global_space = true;
935 break;
936 case PART_ROT_GLOB_X:
937 case PART_ROT_GLOB_Y:
938 case PART_ROT_GLOB_Z:
939 rot_vec[part->rotmode - PART_ROT_GLOB_X] = 1.0f;
940 use_global_space = true;
941 break;
942 case PART_ROT_OB_X:
943 case PART_ROT_OB_Y:
944 case PART_ROT_OB_Z:
945 copy_v3_v3(rot_vec, ob->obmat[part->rotmode - PART_ROT_OB_X]);
946 use_global_space = false;
947 break;
948 default:
949 use_global_space = true;
950 break;
951 }
952
953 /* create rotation quat */
954
955 if (use_global_space) {
956 negate_v3(rot_vec);
957 vec_to_quat(q2, rot_vec, OB_POSX, OB_POSZ);
958
959 /* randomize rotation quat */
960 if (part->randrotfac != 0.0f) {
961 interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
962 }
963 else {
964 copy_qt_qt(rot, q2);
965 }
966 }
967 else {
968 /* calculate rotation in local-space */
969 float q_obmat[4];
970 float q_imat[4];
971
972 mat4_to_quat(q_obmat, ob->obmat);
973 invert_qt_qt_normalized(q_imat, q_obmat);
974
975 if (part->rotmode != PART_ROT_NOR_TAN) {
976 float rot_vec_local[3];
977
978 /* rot_vec */
979 negate_v3(rot_vec);
980 copy_v3_v3(rot_vec_local, rot_vec);
981 mul_qt_v3(q_imat, rot_vec_local);
982 normalize_v3(rot_vec_local);
983
984 vec_to_quat(q2, rot_vec_local, OB_POSX, OB_POSZ);
985 }
986 else {
987 /* (part->rotmode == PART_ROT_NOR_TAN) */
988 float tmat[3][3];
989
990 /* note: utan_local is not taken from 'utan', we calculate from rot_vec/vtan */
991 /* note: it looks like rotation phase may be applied twice (once with vtan, again below)
992 * however this isn't the case - campbell */
993 float *rot_vec_local = tmat[0];
994 float *vtan_local = tmat[1];
995 float *utan_local = tmat[2];
996
997 /* use tangents */
998 BLI_assert(use_tangents == true);
999
1000 /* rot_vec */
1001 copy_v3_v3(rot_vec_local, rot_vec);
1002 mul_qt_v3(q_imat, rot_vec_local);
1003
1004 /* vtan_local */
1005 copy_v3_v3(vtan_local, vtan); /* flips, cant use */
1006 mul_qt_v3(q_imat, vtan_local);
1007
1008 /* ensure orthogonal matrix (rot_vec aligned) */
1009 cross_v3_v3v3(utan_local, vtan_local, rot_vec_local);
1010 cross_v3_v3v3(vtan_local, utan_local, rot_vec_local);
1011
1012 /* note: no need to normalize */
1013 mat3_to_quat(q2, tmat);
1014 }
1015
1016 /* randomize rotation quat */
1017 if (part->randrotfac != 0.0f) {
1018 mul_qt_qtqt(r_rot, r_rot, q_imat);
1019 interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
1020 }
1021 else {
1022 copy_qt_qt(rot, q2);
1023 }
1024
1025 mul_qt_qtqt(rot, q_obmat, rot);
1026 }
1027
1028 /* rotation phase */
1029 phasefac = part->phasefac;
1030 if (part->randphasefac != 0.0f) {
1031 phasefac += part->randphasefac * psys_frand(psys, p + 20);
1032 }
1033 axis_angle_to_quat(q_phase, x_vec, phasefac * (float)M_PI);
1034
1035 /* combine base rotation & phase */
1036 mul_qt_qtqt(state->rot, rot, q_phase);
1037 }
1038
1039 /* -angular velocity */
1040
1041 zero_v3(state->ave);
1042
1043 if (part->avemode) {
1044 if (part->avemode == PART_AVE_RAND) {
1045 copy_v3_v3(state->ave, r_ave);
1046 }
1047 else {
1048 get_angular_velocity_vector(part->avemode, state, state->ave);
1049 }
1050
1051 normalize_v3(state->ave);
1052 mul_v3_fl(state->ave, part->avefac);
1053 }
1054 }
1055 }
1056
1057 /* recursively evaluate emitter parent anim at cfra */
evaluate_emitter_anim(struct Depsgraph * depsgraph,Scene * scene,Object * ob,float cfra)1058 static void evaluate_emitter_anim(struct Depsgraph *depsgraph,
1059 Scene *scene,
1060 Object *ob,
1061 float cfra)
1062 {
1063 if (ob->parent) {
1064 evaluate_emitter_anim(depsgraph, scene, ob->parent, cfra);
1065 }
1066
1067 BKE_object_where_is_calc_time(depsgraph, scene, ob, cfra);
1068 }
1069
1070 /* sets particle to the emitter surface with initial velocity & rotation */
reset_particle(ParticleSimulationData * sim,ParticleData * pa,float dtime,float cfra)1071 void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, float cfra)
1072 {
1073 ParticleSystem *psys = sim->psys;
1074 ParticleSettings *part;
1075 ParticleTexture ptex;
1076 int p = pa - psys->particles;
1077 part = psys->part;
1078
1079 /* get precise emitter matrix if particle is born */
1080 if (part->type != PART_HAIR && dtime > 0.f && pa->time < cfra && pa->time >= sim->psys->cfra) {
1081 evaluate_emitter_anim(sim->depsgraph, sim->scene, sim->ob, pa->time);
1082
1083 psys->flag |= PSYS_OB_ANIM_RESTORE;
1084 }
1085
1086 psys_get_birth_coords(sim, pa, &pa->state, dtime, cfra);
1087
1088 /* Initialize particle settings which depends on texture.
1089 *
1090 * We could only do it now because we'll need to know coordinate
1091 * before sampling the texture.
1092 */
1093 init_particle_texture(sim, pa, p);
1094
1095 if (part->phystype == PART_PHYS_BOIDS && pa->boid) {
1096 BoidParticle *bpa = pa->boid;
1097
1098 /* and gravity in r_ve */
1099 bpa->gravity[0] = bpa->gravity[1] = 0.0f;
1100 bpa->gravity[2] = -1.0f;
1101 if ((sim->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) &&
1102 (sim->scene->physics_settings.gravity[2] != 0.0f)) {
1103 bpa->gravity[2] = sim->scene->physics_settings.gravity[2];
1104 }
1105
1106 bpa->data.health = part->boids->health;
1107 bpa->data.mode = eBoidMode_InAir;
1108 bpa->data.state_id = ((BoidState *)part->boids->states.first)->id;
1109 bpa->data.acc[0] = bpa->data.acc[1] = bpa->data.acc[2] = 0.0f;
1110 }
1111
1112 if (part->type == PART_HAIR) {
1113 pa->lifetime = 100.0f;
1114 }
1115 else {
1116 /* initialize the lifetime, in case the texture coordinates
1117 * are from Particles/Strands, which would cause undefined values
1118 */
1119 pa->lifetime = part->lifetime * (1.0f - part->randlife * psys_frand(psys, p + 21));
1120 pa->dietime = pa->time + pa->lifetime;
1121
1122 /* get possible textural influence */
1123 psys_get_texture(sim, pa, &ptex, PAMAP_LIFE, cfra);
1124
1125 pa->lifetime = part->lifetime * ptex.life;
1126
1127 if (part->randlife != 0.0f) {
1128 pa->lifetime *= 1.0f - part->randlife * psys_frand(psys, p + 21);
1129 }
1130 }
1131
1132 pa->dietime = pa->time + pa->lifetime;
1133
1134 if ((sim->psys->pointcache) && (sim->psys->pointcache->flag & PTCACHE_BAKED) &&
1135 (sim->psys->pointcache->mem_cache.first)) {
1136 float dietime = psys_get_dietime_from_cache(sim->psys->pointcache, p);
1137 pa->dietime = MIN2(pa->dietime, dietime);
1138 }
1139
1140 if (pa->time > cfra) {
1141 pa->alive = PARS_UNBORN;
1142 }
1143 else if (pa->dietime <= cfra) {
1144 pa->alive = PARS_DEAD;
1145 }
1146 else {
1147 pa->alive = PARS_ALIVE;
1148 }
1149
1150 pa->state.time = cfra;
1151 }
reset_all_particles(ParticleSimulationData * sim,float dtime,float cfra,int from)1152 static void reset_all_particles(ParticleSimulationData *sim, float dtime, float cfra, int from)
1153 {
1154 ParticleData *pa;
1155 int p, totpart = sim->psys->totpart;
1156
1157 for (p = from, pa = sim->psys->particles + from; p < totpart; p++, pa++) {
1158 reset_particle(sim, pa, dtime, cfra);
1159 }
1160 }
1161 /************************************************/
1162 /* Particle targets */
1163 /************************************************/
psys_get_target_system(Object * ob,ParticleTarget * pt)1164 ParticleSystem *psys_get_target_system(Object *ob, ParticleTarget *pt)
1165 {
1166 ParticleSystem *psys = NULL;
1167
1168 if (pt->ob == NULL || pt->ob == ob) {
1169 psys = BLI_findlink(&ob->particlesystem, pt->psys - 1);
1170 }
1171 else {
1172 psys = BLI_findlink(&pt->ob->particlesystem, pt->psys - 1);
1173 }
1174
1175 if (psys) {
1176 pt->flag |= PTARGET_VALID;
1177 }
1178 else {
1179 pt->flag &= ~PTARGET_VALID;
1180 }
1181
1182 return psys;
1183 }
1184 /************************************************/
1185 /* Keyed particles */
1186 /************************************************/
1187 /* Counts valid keyed targets */
psys_count_keyed_targets(ParticleSimulationData * sim)1188 void psys_count_keyed_targets(ParticleSimulationData *sim)
1189 {
1190 ParticleSystem *psys = sim->psys, *kpsys;
1191 ParticleTarget *pt = psys->targets.first;
1192 int keys_valid = 1;
1193 psys->totkeyed = 0;
1194
1195 for (; pt; pt = pt->next) {
1196 kpsys = psys_get_target_system(sim->ob, pt);
1197
1198 if (kpsys && kpsys->totpart) {
1199 psys->totkeyed += keys_valid;
1200 if (psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f) {
1201 psys->totkeyed += 1;
1202 }
1203 }
1204 else {
1205 keys_valid = 0;
1206 }
1207 }
1208
1209 psys->totkeyed *= psys->flag & PSYS_KEYED_TIMING ? 1 : psys->part->keyed_loops;
1210 }
1211
set_keyed_keys(ParticleSimulationData * sim)1212 static void set_keyed_keys(ParticleSimulationData *sim)
1213 {
1214 ParticleSystem *psys = sim->psys;
1215 ParticleSimulationData ksim = {0};
1216 ParticleTarget *pt;
1217 PARTICLE_P;
1218 ParticleKey *key;
1219 int totpart = psys->totpart, k, totkeys = psys->totkeyed;
1220 int keyed_flag = 0;
1221
1222 ksim.depsgraph = sim->depsgraph;
1223 ksim.scene = sim->scene;
1224
1225 /* no proper targets so let's clear and bail out */
1226 if (psys->totkeyed == 0) {
1227 free_keyed_keys(psys);
1228 psys->flag &= ~PSYS_KEYED;
1229 return;
1230 }
1231
1232 if (totpart && psys->particles->totkey != totkeys) {
1233 free_keyed_keys(psys);
1234
1235 key = MEM_callocN(totpart * totkeys * sizeof(ParticleKey), "Keyed keys");
1236
1237 LOOP_PARTICLES
1238 {
1239 pa->keys = key;
1240 pa->totkey = totkeys;
1241 key += totkeys;
1242 }
1243 }
1244
1245 psys->flag &= ~PSYS_KEYED;
1246
1247 pt = psys->targets.first;
1248 for (k = 0; k < totkeys; k++) {
1249 ksim.ob = pt->ob ? pt->ob : sim->ob;
1250 ksim.psys = BLI_findlink(&ksim.ob->particlesystem, pt->psys - 1);
1251 keyed_flag = (ksim.psys->flag & PSYS_KEYED);
1252 ksim.psys->flag &= ~PSYS_KEYED;
1253
1254 LOOP_PARTICLES
1255 {
1256 key = pa->keys + k;
1257 key->time = -1.0; /* use current time */
1258
1259 const int p_ksim = (ksim.psys->totpart) ? p % ksim.psys->totpart : 0;
1260 psys_get_particle_state(&ksim, p_ksim, key, 1);
1261
1262 if (psys->flag & PSYS_KEYED_TIMING) {
1263 key->time = pa->time + pt->time;
1264 if (pt->duration != 0.0f && k + 1 < totkeys) {
1265 copy_particle_key(key + 1, key, 1);
1266 (key + 1)->time = pa->time + pt->time + pt->duration;
1267 }
1268 }
1269 else if (totkeys > 1) {
1270 key->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
1271 }
1272 else {
1273 key->time = pa->time;
1274 }
1275 }
1276
1277 if (psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f) {
1278 k++;
1279 }
1280
1281 ksim.psys->flag |= keyed_flag;
1282
1283 pt = (pt->next && pt->next->flag & PTARGET_VALID) ? pt->next : psys->targets.first;
1284 }
1285
1286 psys->flag |= PSYS_KEYED;
1287 }
1288
1289 /************************************************/
1290 /* Point Cache */
1291 /************************************************/
psys_make_temp_pointcache(Object * ob,ParticleSystem * psys)1292 void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
1293 {
1294 PointCache *cache = psys->pointcache;
1295
1296 if (cache->flag & PTCACHE_DISK_CACHE && BLI_listbase_is_empty(&cache->mem_cache)) {
1297 PTCacheID pid;
1298 BKE_ptcache_id_from_particles(&pid, ob, psys);
1299 cache->flag &= ~PTCACHE_DISK_CACHE;
1300 BKE_ptcache_disk_to_mem(&pid);
1301 cache->flag |= PTCACHE_DISK_CACHE;
1302 }
1303 }
psys_clear_temp_pointcache(ParticleSystem * psys)1304 static void psys_clear_temp_pointcache(ParticleSystem *psys)
1305 {
1306 if (psys->pointcache->flag & PTCACHE_DISK_CACHE) {
1307 BKE_ptcache_free_mem(&psys->pointcache->mem_cache);
1308 }
1309 }
psys_get_pointcache_start_end(Scene * scene,ParticleSystem * psys,int * sfra,int * efra)1310 void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra, int *efra)
1311 {
1312 ParticleSettings *part = psys->part;
1313
1314 *sfra = max_ii(1, (int)part->sta);
1315 *efra = min_ii((int)(part->end + part->lifetime + 1.0f), max_ii(scene->r.pefra, scene->r.efra));
1316 }
1317
1318 /************************************************/
1319 /* Effectors */
1320 /************************************************/
psys_update_particle_bvhtree(ParticleSystem * psys,float cfra)1321 static void psys_update_particle_bvhtree(ParticleSystem *psys, float cfra)
1322 {
1323 if (psys) {
1324 PARTICLE_P;
1325 int totpart = 0;
1326 bool need_rebuild;
1327
1328 BLI_rw_mutex_lock(&psys_bvhtree_rwlock, THREAD_LOCK_READ);
1329 need_rebuild = !psys->bvhtree || psys->bvhtree_frame != cfra;
1330 BLI_rw_mutex_unlock(&psys_bvhtree_rwlock);
1331
1332 if (need_rebuild) {
1333 LOOP_SHOWN_PARTICLES
1334 {
1335 totpart++;
1336 }
1337
1338 BLI_rw_mutex_lock(&psys_bvhtree_rwlock, THREAD_LOCK_WRITE);
1339
1340 BLI_bvhtree_free(psys->bvhtree);
1341 psys->bvhtree = BLI_bvhtree_new(totpart, 0.0, 4, 6);
1342
1343 LOOP_SHOWN_PARTICLES
1344 {
1345 if (pa->alive == PARS_ALIVE) {
1346 if (pa->state.time == cfra) {
1347 BLI_bvhtree_insert(psys->bvhtree, p, pa->prev_state.co, 1);
1348 }
1349 else {
1350 BLI_bvhtree_insert(psys->bvhtree, p, pa->state.co, 1);
1351 }
1352 }
1353 }
1354 BLI_bvhtree_balance(psys->bvhtree);
1355
1356 psys->bvhtree_frame = cfra;
1357
1358 BLI_rw_mutex_unlock(&psys_bvhtree_rwlock);
1359 }
1360 }
1361 }
psys_update_particle_tree(ParticleSystem * psys,float cfra)1362 void psys_update_particle_tree(ParticleSystem *psys, float cfra)
1363 {
1364 if (psys) {
1365 PARTICLE_P;
1366 int totpart = 0;
1367
1368 if (!psys->tree || psys->tree_frame != cfra) {
1369 LOOP_SHOWN_PARTICLES
1370 {
1371 totpart++;
1372 }
1373
1374 BLI_kdtree_3d_free(psys->tree);
1375 psys->tree = BLI_kdtree_3d_new(psys->totpart);
1376
1377 LOOP_SHOWN_PARTICLES
1378 {
1379 if (pa->alive == PARS_ALIVE) {
1380 if (pa->state.time == cfra) {
1381 BLI_kdtree_3d_insert(psys->tree, p, pa->prev_state.co);
1382 }
1383 else {
1384 BLI_kdtree_3d_insert(psys->tree, p, pa->state.co);
1385 }
1386 }
1387 }
1388 BLI_kdtree_3d_balance(psys->tree);
1389
1390 psys->tree_frame = cfra;
1391 }
1392 }
1393 }
1394
psys_update_effectors(ParticleSimulationData * sim)1395 static void psys_update_effectors(ParticleSimulationData *sim)
1396 {
1397 BKE_effectors_free(sim->psys->effectors);
1398 sim->psys->effectors = BKE_effectors_create(
1399 sim->depsgraph, sim->ob, sim->psys, sim->psys->part->effector_weights);
1400 precalc_guides(sim, sim->psys->effectors);
1401 }
1402
integrate_particle(ParticleSettings * part,ParticleData * pa,float dtime,float * external_acceleration,void (* force_func)(void * forcedata,ParticleKey * state,float * force,float * impulse),void * forcedata)1403 static void integrate_particle(
1404 ParticleSettings *part,
1405 ParticleData *pa,
1406 float dtime,
1407 float *external_acceleration,
1408 void (*force_func)(void *forcedata, ParticleKey *state, float *force, float *impulse),
1409 void *forcedata)
1410 {
1411 #define ZERO_F43 \
1412 { \
1413 {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, \
1414 { \
1415 0.0f, 0.0f, 0.0f \
1416 } \
1417 }
1418
1419 ParticleKey states[5];
1420 float force[3], acceleration[3], impulse[3], dx[4][3] = ZERO_F43, dv[4][3] = ZERO_F43, oldpos[3];
1421 float pa_mass = (part->flag & PART_SIZEMASS ? part->mass * pa->size : part->mass);
1422 int i, steps = 1;
1423 int integrator = part->integrator;
1424
1425 #undef ZERO_F43
1426
1427 copy_v3_v3(oldpos, pa->state.co);
1428
1429 /* Verlet integration behaves strangely with moving emitters, so do first step with euler. */
1430 if (pa->prev_state.time < 0.f && integrator == PART_INT_VERLET) {
1431 integrator = PART_INT_EULER;
1432 }
1433
1434 switch (integrator) {
1435 case PART_INT_EULER:
1436 steps = 1;
1437 break;
1438 case PART_INT_MIDPOINT:
1439 steps = 2;
1440 break;
1441 case PART_INT_RK4:
1442 steps = 4;
1443 break;
1444 case PART_INT_VERLET:
1445 steps = 1;
1446 break;
1447 }
1448
1449 for (i = 0; i < steps; i++) {
1450 copy_particle_key(states + i, &pa->state, 1);
1451 }
1452
1453 states->time = 0.f;
1454
1455 for (i = 0; i < steps; i++) {
1456 zero_v3(force);
1457 zero_v3(impulse);
1458
1459 force_func(forcedata, states + i, force, impulse);
1460
1461 /* force to acceleration*/
1462 mul_v3_v3fl(acceleration, force, 1.0f / pa_mass);
1463
1464 if (external_acceleration) {
1465 add_v3_v3(acceleration, external_acceleration);
1466 }
1467
1468 /* calculate next state */
1469 add_v3_v3(states[i].vel, impulse);
1470
1471 switch (integrator) {
1472 case PART_INT_EULER:
1473 madd_v3_v3v3fl(pa->state.co, states->co, states->vel, dtime);
1474 madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
1475 break;
1476 case PART_INT_MIDPOINT:
1477 if (i == 0) {
1478 madd_v3_v3v3fl(states[1].co, states->co, states->vel, dtime * 0.5f);
1479 madd_v3_v3v3fl(states[1].vel, states->vel, acceleration, dtime * 0.5f);
1480 states[1].time = dtime * 0.5f;
1481 /*fra=sim->psys->cfra+0.5f*dfra;*/
1482 }
1483 else {
1484 madd_v3_v3v3fl(pa->state.co, states->co, states[1].vel, dtime);
1485 madd_v3_v3v3fl(pa->state.vel, states->vel, acceleration, dtime);
1486 }
1487 break;
1488 case PART_INT_RK4:
1489 switch (i) {
1490 case 0:
1491 copy_v3_v3(dx[0], states->vel);
1492 mul_v3_fl(dx[0], dtime);
1493 copy_v3_v3(dv[0], acceleration);
1494 mul_v3_fl(dv[0], dtime);
1495
1496 madd_v3_v3v3fl(states[1].co, states->co, dx[0], 0.5f);
1497 madd_v3_v3v3fl(states[1].vel, states->vel, dv[0], 0.5f);
1498 states[1].time = dtime * 0.5f;
1499 /*fra=sim->psys->cfra+0.5f*dfra;*/
1500 break;
1501 case 1:
1502 madd_v3_v3v3fl(dx[1], states->vel, dv[0], 0.5f);
1503 mul_v3_fl(dx[1], dtime);
1504 copy_v3_v3(dv[1], acceleration);
1505 mul_v3_fl(dv[1], dtime);
1506
1507 madd_v3_v3v3fl(states[2].co, states->co, dx[1], 0.5f);
1508 madd_v3_v3v3fl(states[2].vel, states->vel, dv[1], 0.5f);
1509 states[2].time = dtime * 0.5f;
1510 break;
1511 case 2:
1512 madd_v3_v3v3fl(dx[2], states->vel, dv[1], 0.5f);
1513 mul_v3_fl(dx[2], dtime);
1514 copy_v3_v3(dv[2], acceleration);
1515 mul_v3_fl(dv[2], dtime);
1516
1517 add_v3_v3v3(states[3].co, states->co, dx[2]);
1518 add_v3_v3v3(states[3].vel, states->vel, dv[2]);
1519 states[3].time = dtime;
1520 /*fra=cfra;*/
1521 break;
1522 case 3:
1523 add_v3_v3v3(dx[3], states->vel, dv[2]);
1524 mul_v3_fl(dx[3], dtime);
1525 copy_v3_v3(dv[3], acceleration);
1526 mul_v3_fl(dv[3], dtime);
1527
1528 madd_v3_v3v3fl(pa->state.co, states->co, dx[0], 1.0f / 6.0f);
1529 madd_v3_v3fl(pa->state.co, dx[1], 1.0f / 3.0f);
1530 madd_v3_v3fl(pa->state.co, dx[2], 1.0f / 3.0f);
1531 madd_v3_v3fl(pa->state.co, dx[3], 1.0f / 6.0f);
1532
1533 madd_v3_v3v3fl(pa->state.vel, states->vel, dv[0], 1.0f / 6.0f);
1534 madd_v3_v3fl(pa->state.vel, dv[1], 1.0f / 3.0f);
1535 madd_v3_v3fl(pa->state.vel, dv[2], 1.0f / 3.0f);
1536 madd_v3_v3fl(pa->state.vel, dv[3], 1.0f / 6.0f);
1537 }
1538 break;
1539 case PART_INT_VERLET: /* Verlet integration */
1540 madd_v3_v3v3fl(pa->state.vel, pa->prev_state.vel, acceleration, dtime);
1541 madd_v3_v3v3fl(pa->state.co, pa->prev_state.co, pa->state.vel, dtime);
1542
1543 sub_v3_v3v3(pa->state.vel, pa->state.co, oldpos);
1544 mul_v3_fl(pa->state.vel, 1.0f / dtime);
1545 break;
1546 }
1547 }
1548 }
1549
1550 /* -------------------------------------------------------------------- */
1551 /** \name SPH fluid physics
1552 *
1553 * In theory, there could be unlimited implementation of SPH simulators
1554 *
1555 * This code uses in some parts adapted algorithms
1556 * from the pseudo code as outlined in the Research paper:
1557 *
1558 * Titled: Particle-based Viscoelastic Fluid Simulation.
1559 * Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
1560 * Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
1561 *
1562 * Presented at Siggraph, (2005)
1563 *
1564 * \{ */
1565
1566 #define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
sph_spring_add(ParticleSystem * psys,ParticleSpring * spring)1567 static ParticleSpring *sph_spring_add(ParticleSystem *psys, ParticleSpring *spring)
1568 {
1569 /* Are more refs required? */
1570 if (psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
1571 psys->alloc_fluidsprings = PSYS_FLUID_SPRINGS_INITIAL_SIZE;
1572 psys->fluid_springs = (ParticleSpring *)MEM_callocN(
1573 psys->alloc_fluidsprings * sizeof(ParticleSpring), "Particle Fluid Springs");
1574 }
1575 else if (psys->tot_fluidsprings == psys->alloc_fluidsprings) {
1576 /* Double the number of refs allocated */
1577 psys->alloc_fluidsprings *= 2;
1578 psys->fluid_springs = (ParticleSpring *)MEM_reallocN(
1579 psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
1580 }
1581
1582 memcpy(psys->fluid_springs + psys->tot_fluidsprings, spring, sizeof(ParticleSpring));
1583 psys->tot_fluidsprings++;
1584
1585 return psys->fluid_springs + psys->tot_fluidsprings - 1;
1586 }
sph_spring_delete(ParticleSystem * psys,int j)1587 static void sph_spring_delete(ParticleSystem *psys, int j)
1588 {
1589 if (j != psys->tot_fluidsprings - 1) {
1590 psys->fluid_springs[j] = psys->fluid_springs[psys->tot_fluidsprings - 1];
1591 }
1592
1593 psys->tot_fluidsprings--;
1594
1595 if (psys->tot_fluidsprings < psys->alloc_fluidsprings / 2 &&
1596 psys->alloc_fluidsprings > PSYS_FLUID_SPRINGS_INITIAL_SIZE) {
1597 psys->alloc_fluidsprings /= 2;
1598 psys->fluid_springs = (ParticleSpring *)MEM_reallocN(
1599 psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
1600 }
1601 }
sph_springs_modify(ParticleSystem * psys,float dtime)1602 static void sph_springs_modify(ParticleSystem *psys, float dtime)
1603 {
1604 SPHFluidSettings *fluid = psys->part->fluid;
1605 ParticleData *pa1, *pa2;
1606 ParticleSpring *spring = psys->fluid_springs;
1607
1608 float h, d, Rij[3], rij, Lij;
1609 int i;
1610
1611 float yield_ratio = fluid->yield_ratio;
1612 float plasticity = fluid->plasticity_constant;
1613 /* scale things according to dtime */
1614 float timefix = 25.f * dtime;
1615
1616 if ((fluid->flag & SPH_VISCOELASTIC_SPRINGS) == 0 || fluid->spring_k == 0.f) {
1617 return;
1618 }
1619
1620 /* Loop through the springs */
1621 for (i = 0; i < psys->tot_fluidsprings; i++, spring++) {
1622 pa1 = psys->particles + spring->particle_index[0];
1623 pa2 = psys->particles + spring->particle_index[1];
1624
1625 sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
1626 rij = normalize_v3(Rij);
1627
1628 /* adjust rest length */
1629 Lij = spring->rest_length;
1630 d = yield_ratio * timefix * Lij;
1631
1632 if (rij > Lij + d) { // Stretch
1633 spring->rest_length += plasticity * (rij - Lij - d) * timefix;
1634 }
1635 else if (rij < Lij - d) { // Compress
1636 spring->rest_length -= plasticity * (Lij - d - rij) * timefix;
1637 }
1638
1639 h = 4.f * pa1->size;
1640
1641 if (spring->rest_length > h) {
1642 spring->delete_flag = 1;
1643 }
1644 }
1645
1646 /* Loop through springs backwaqrds - for efficient delete function */
1647 for (i = psys->tot_fluidsprings - 1; i >= 0; i--) {
1648 if (psys->fluid_springs[i].delete_flag) {
1649 sph_spring_delete(psys, i);
1650 }
1651 }
1652 }
sph_springhash_build(ParticleSystem * psys)1653 static EdgeHash *sph_springhash_build(ParticleSystem *psys)
1654 {
1655 EdgeHash *springhash = NULL;
1656 ParticleSpring *spring;
1657 int i = 0;
1658
1659 springhash = BLI_edgehash_new_ex(__func__, psys->tot_fluidsprings);
1660
1661 for (i = 0, spring = psys->fluid_springs; i < psys->tot_fluidsprings; i++, spring++) {
1662 BLI_edgehash_insert(
1663 springhash, spring->particle_index[0], spring->particle_index[1], POINTER_FROM_INT(i + 1));
1664 }
1665
1666 return springhash;
1667 }
1668
1669 #define SPH_NEIGHBORS 512
1670 typedef struct SPHNeighbor {
1671 ParticleSystem *psys;
1672 int index;
1673 } SPHNeighbor;
1674
1675 typedef struct SPHRangeData {
1676 SPHNeighbor neighbors[SPH_NEIGHBORS];
1677 int tot_neighbors;
1678
1679 float *data;
1680
1681 ParticleSystem *npsys;
1682 ParticleData *pa;
1683
1684 float h;
1685 float mass;
1686 float massfac;
1687 int use_size;
1688 } SPHRangeData;
1689
sph_evaluate_func(BVHTree * tree,ParticleSystem ** psys,const float co[3],SPHRangeData * pfr,float interaction_radius,BVHTree_RangeQuery callback)1690 static void sph_evaluate_func(BVHTree *tree,
1691 ParticleSystem **psys,
1692 const float co[3],
1693 SPHRangeData *pfr,
1694 float interaction_radius,
1695 BVHTree_RangeQuery callback)
1696 {
1697 int i;
1698
1699 pfr->tot_neighbors = 0;
1700
1701 for (i = 0; i < 10 && psys[i]; i++) {
1702 pfr->npsys = psys[i];
1703 pfr->massfac = psys[i]->part->mass / pfr->mass;
1704 pfr->use_size = psys[i]->part->flag & PART_SIZEMASS;
1705
1706 if (tree) {
1707 BLI_bvhtree_range_query(tree, co, interaction_radius, callback, pfr);
1708 break;
1709 }
1710
1711 BLI_rw_mutex_lock(&psys_bvhtree_rwlock, THREAD_LOCK_READ);
1712
1713 BLI_bvhtree_range_query(psys[i]->bvhtree, co, interaction_radius, callback, pfr);
1714
1715 BLI_rw_mutex_unlock(&psys_bvhtree_rwlock);
1716 }
1717 }
sph_density_accum_cb(void * userdata,int index,const float co[3],float squared_dist)1718 static void sph_density_accum_cb(void *userdata, int index, const float co[3], float squared_dist)
1719 {
1720 SPHRangeData *pfr = (SPHRangeData *)userdata;
1721 ParticleData *npa = pfr->npsys->particles + index;
1722 float q;
1723 float dist;
1724
1725 UNUSED_VARS(co);
1726
1727 if (npa == pfr->pa || squared_dist < FLT_EPSILON) {
1728 return;
1729 }
1730
1731 /* Ugh! One particle has too many neighbors! If some aren't taken into
1732 * account, the forces will be biased by the tree search order. This
1733 * effectively adds energy to the system, and results in a churning motion.
1734 * But, we have to stop somewhere, and it's not the end of the world.
1735 * - jahka and z0r
1736 */
1737 if (pfr->tot_neighbors >= SPH_NEIGHBORS) {
1738 return;
1739 }
1740
1741 pfr->neighbors[pfr->tot_neighbors].index = index;
1742 pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
1743 pfr->tot_neighbors++;
1744
1745 dist = sqrtf(squared_dist);
1746 q = (1.f - dist / pfr->h) * pfr->massfac;
1747
1748 if (pfr->use_size) {
1749 q *= npa->size;
1750 }
1751
1752 pfr->data[0] += q * q;
1753 pfr->data[1] += q * q * q;
1754 }
1755
1756 /*
1757 * Find the Courant number for an SPH particle (used for adaptive time step).
1758 */
sph_particle_courant(SPHData * sphdata,SPHRangeData * pfr)1759 static void sph_particle_courant(SPHData *sphdata, SPHRangeData *pfr)
1760 {
1761 ParticleData *pa, *npa;
1762 int i;
1763 float flow[3], offset[3], dist;
1764
1765 zero_v3(flow);
1766
1767 dist = 0.0f;
1768 if (pfr->tot_neighbors > 0) {
1769 pa = pfr->pa;
1770 for (i = 0; i < pfr->tot_neighbors; i++) {
1771 npa = pfr->neighbors[i].psys->particles + pfr->neighbors[i].index;
1772 sub_v3_v3v3(offset, pa->prev_state.co, npa->prev_state.co);
1773 dist += len_v3(offset);
1774 add_v3_v3(flow, npa->prev_state.vel);
1775 }
1776 dist += sphdata->psys[0]->part->fluid->radius; // TODO: remove this? - z0r
1777 sphdata->element_size = dist / pfr->tot_neighbors;
1778 mul_v3_v3fl(sphdata->flow, flow, 1.0f / pfr->tot_neighbors);
1779 }
1780 else {
1781 sphdata->element_size = FLT_MAX;
1782 copy_v3_v3(sphdata->flow, flow);
1783 }
1784 }
sph_force_cb(void * sphdata_v,ParticleKey * state,float * force,float * UNUSED (impulse))1785 static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse))
1786 {
1787 SPHData *sphdata = (SPHData *)sphdata_v;
1788 ParticleSystem **psys = sphdata->psys;
1789 ParticleData *pa = sphdata->pa;
1790 SPHFluidSettings *fluid = psys[0]->part->fluid;
1791 ParticleSpring *spring = NULL;
1792 SPHRangeData pfr;
1793 SPHNeighbor *pfn;
1794 float *gravity = sphdata->gravity;
1795 EdgeHash *springhash = sphdata->eh;
1796
1797 float q, u, rij, dv[3];
1798 float pressure, near_pressure;
1799
1800 float visc = fluid->viscosity_omega;
1801 float stiff_visc = fluid->viscosity_beta *
1802 (fluid->flag & SPH_FAC_VISCOSITY ? fluid->viscosity_omega : 1.f);
1803
1804 float inv_mass = 1.0f / sphdata->mass;
1805 float spring_constant = fluid->spring_k;
1806
1807 /* 4.0 seems to be a pretty good value */
1808 float interaction_radius = fluid->radius *
1809 (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
1810 float h = interaction_radius * sphdata->hfac;
1811 /* 4.77 is an experimentally determined density factor */
1812 float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f);
1813 float rest_length = fluid->rest_length *
1814 (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f);
1815
1816 float stiffness = fluid->stiffness_k;
1817 float stiffness_near_fac = fluid->stiffness_knear *
1818 (fluid->flag & SPH_FAC_REPULSION ? fluid->stiffness_k : 1.f);
1819
1820 ParticleData *npa;
1821 float vec[3];
1822 float vel[3];
1823 float co[3];
1824 float data[2];
1825 float density, near_density;
1826
1827 int i, spring_index, index = pa - psys[0]->particles;
1828
1829 data[0] = data[1] = 0;
1830 pfr.data = data;
1831 pfr.h = h;
1832 pfr.pa = pa;
1833 pfr.mass = sphdata->mass;
1834
1835 sph_evaluate_func(NULL, psys, state->co, &pfr, interaction_radius, sph_density_accum_cb);
1836
1837 density = data[0];
1838 near_density = data[1];
1839
1840 pressure = stiffness * (density - rest_density);
1841 near_pressure = stiffness_near_fac * near_density;
1842
1843 pfn = pfr.neighbors;
1844 for (i = 0; i < pfr.tot_neighbors; i++, pfn++) {
1845 npa = pfn->psys->particles + pfn->index;
1846
1847 madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
1848
1849 sub_v3_v3v3(vec, co, state->co);
1850 rij = normalize_v3(vec);
1851
1852 q = (1.f - rij / h) * pfn->psys->part->mass * inv_mass;
1853
1854 if (pfn->psys->part->flag & PART_SIZEMASS) {
1855 q *= npa->size;
1856 }
1857
1858 copy_v3_v3(vel, npa->prev_state.vel);
1859
1860 /* Double Density Relaxation */
1861 madd_v3_v3fl(force, vec, -(pressure + near_pressure * q) * q);
1862
1863 /* Viscosity */
1864 if (visc > 0.f || stiff_visc > 0.f) {
1865 sub_v3_v3v3(dv, vel, state->vel);
1866 u = dot_v3v3(vec, dv);
1867
1868 if (u < 0.f && visc > 0.f) {
1869 madd_v3_v3fl(force, vec, 0.5f * q * visc * u);
1870 }
1871
1872 if (u > 0.f && stiff_visc > 0.f) {
1873 madd_v3_v3fl(force, vec, 0.5f * q * stiff_visc * u);
1874 }
1875 }
1876
1877 if (spring_constant > 0.f) {
1878 /* Viscoelastic spring force */
1879 if (pfn->psys == psys[0] && fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash) {
1880 /* BLI_edgehash_lookup appears to be thread-safe. - z0r */
1881 spring_index = POINTER_AS_INT(BLI_edgehash_lookup(springhash, index, pfn->index));
1882
1883 if (spring_index) {
1884 spring = psys[0]->fluid_springs + spring_index - 1;
1885
1886 madd_v3_v3fl(
1887 force, vec, -10.f * spring_constant * (1.f - rij / h) * (spring->rest_length - rij));
1888 }
1889 else if (fluid->spring_frames == 0 ||
1890 (pa->prev_state.time - pa->time) <= fluid->spring_frames) {
1891 ParticleSpring temp_spring;
1892 temp_spring.particle_index[0] = index;
1893 temp_spring.particle_index[1] = pfn->index;
1894 temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : rest_length;
1895 temp_spring.delete_flag = 0;
1896
1897 BLI_buffer_append(&sphdata->new_springs, ParticleSpring, temp_spring);
1898 }
1899 }
1900 else { /* PART_SPRING_HOOKES - Hooke's spring force */
1901 madd_v3_v3fl(force, vec, -10.f * spring_constant * (1.f - rij / h) * (rest_length - rij));
1902 }
1903 }
1904 }
1905
1906 /* Artificial buoyancy force in negative gravity direction */
1907 if (fluid->buoyancy > 0.f && gravity) {
1908 madd_v3_v3fl(force, gravity, fluid->buoyancy * (density - rest_density));
1909 }
1910
1911 if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF) {
1912 sph_particle_courant(sphdata, &pfr);
1913 }
1914 sphdata->pass++;
1915 }
1916
sphclassical_density_accum_cb(void * userdata,int index,const float co[3],float UNUSED (squared_dist))1917 static void sphclassical_density_accum_cb(void *userdata,
1918 int index,
1919 const float co[3],
1920 float UNUSED(squared_dist))
1921 {
1922 SPHRangeData *pfr = (SPHRangeData *)userdata;
1923 ParticleData *npa = pfr->npsys->particles + index;
1924 float q;
1925 float qfac = 21.0f / (256.f * (float)M_PI);
1926 float rij, rij_h;
1927 float vec[3];
1928
1929 /* Exclude particles that are more than 2h away. Can't use squared_dist here
1930 * because it is not accurate enough. Use current state, i.e. the output of
1931 * basic_integrate() - z0r */
1932 sub_v3_v3v3(vec, npa->state.co, co);
1933 rij = len_v3(vec);
1934 rij_h = rij / pfr->h;
1935 if (rij_h > 2.0f) {
1936 return;
1937 }
1938
1939 /* Smoothing factor. Utilize the Wendland kernel. gnuplot:
1940 * q1(x) = (2.0 - x)**4 * ( 1.0 + 2.0 * x)
1941 * plot [0:2] q1(x) */
1942 q = qfac / pow3f(pfr->h) * pow4f(2.0f - rij_h) * (1.0f + 2.0f * rij_h);
1943 q *= pfr->npsys->part->mass;
1944
1945 if (pfr->use_size) {
1946 q *= pfr->pa->size;
1947 }
1948
1949 pfr->data[0] += q;
1950 pfr->data[1] += q / npa->sphdensity;
1951 }
1952
sphclassical_neighbor_accum_cb(void * userdata,int index,const float co[3],float UNUSED (squared_dist))1953 static void sphclassical_neighbor_accum_cb(void *userdata,
1954 int index,
1955 const float co[3],
1956 float UNUSED(squared_dist))
1957 {
1958 SPHRangeData *pfr = (SPHRangeData *)userdata;
1959 ParticleData *npa = pfr->npsys->particles + index;
1960 float rij, rij_h;
1961 float vec[3];
1962
1963 if (pfr->tot_neighbors >= SPH_NEIGHBORS) {
1964 return;
1965 }
1966
1967 /* Exclude particles that are more than 2h away. Can't use squared_dist here
1968 * because it is not accurate enough. Use current state, i.e. the output of
1969 * basic_integrate() - z0r */
1970 sub_v3_v3v3(vec, npa->state.co, co);
1971 rij = len_v3(vec);
1972 rij_h = rij / pfr->h;
1973 if (rij_h > 2.0f) {
1974 return;
1975 }
1976
1977 pfr->neighbors[pfr->tot_neighbors].index = index;
1978 pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys;
1979 pfr->tot_neighbors++;
1980 }
sphclassical_force_cb(void * sphdata_v,ParticleKey * state,float * force,float * UNUSED (impulse))1981 static void sphclassical_force_cb(void *sphdata_v,
1982 ParticleKey *state,
1983 float *force,
1984 float *UNUSED(impulse))
1985 {
1986 SPHData *sphdata = (SPHData *)sphdata_v;
1987 ParticleSystem **psys = sphdata->psys;
1988 ParticleData *pa = sphdata->pa;
1989 SPHFluidSettings *fluid = psys[0]->part->fluid;
1990 SPHRangeData pfr;
1991 SPHNeighbor *pfn;
1992 float *gravity = sphdata->gravity;
1993
1994 float dq, u, rij, dv[3];
1995 float pressure, npressure;
1996
1997 float visc = fluid->viscosity_omega;
1998
1999 float interaction_radius;
2000 float h, hinv;
2001 /* 4.77 is an experimentally determined density factor */
2002 float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.0f);
2003
2004 /* Use speed of sound squared */
2005 float stiffness = pow2f(fluid->stiffness_k);
2006
2007 ParticleData *npa;
2008 float vec[3];
2009 float co[3];
2010 float pressureTerm;
2011
2012 int i;
2013
2014 float qfac2 = 42.0f / (256.0f * (float)M_PI);
2015 float rij_h;
2016
2017 /* 4.0 here is to be consistent with previous formulation/interface */
2018 interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f);
2019 h = interaction_radius * sphdata->hfac;
2020 hinv = 1.0f / h;
2021
2022 pfr.h = h;
2023 pfr.pa = pa;
2024
2025 sph_evaluate_func(
2026 NULL, psys, state->co, &pfr, interaction_radius, sphclassical_neighbor_accum_cb);
2027 pressure = stiffness * (pow7f(pa->sphdensity / rest_density) - 1.0f);
2028
2029 /* multiply by mass so that we return a force, not accel */
2030 qfac2 *= sphdata->mass / pow3f(pfr.h);
2031
2032 pfn = pfr.neighbors;
2033 for (i = 0; i < pfr.tot_neighbors; i++, pfn++) {
2034 npa = pfn->psys->particles + pfn->index;
2035 if (npa == pa) {
2036 /* we do not contribute to ourselves */
2037 continue;
2038 }
2039
2040 /* Find vector to neighbor. Exclude particles that are more than 2h
2041 * away. Can't use current state here because it may have changed on
2042 * another thread - so do own mini integration. Unlike basic_integrate,
2043 * SPH integration depends on neighboring particles. - z0r */
2044 madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time);
2045 sub_v3_v3v3(vec, co, state->co);
2046 rij = normalize_v3(vec);
2047 rij_h = rij / pfr.h;
2048 if (rij_h > 2.0f) {
2049 continue;
2050 }
2051
2052 npressure = stiffness * (pow7f(npa->sphdensity / rest_density) - 1.0f);
2053
2054 /* First derivative of smoothing factor. Utilize the Wendland kernel.
2055 * gnuplot:
2056 * q2(x) = 2.0 * (2.0 - x)**4 - 4.0 * (2.0 - x)**3 * (1.0 + 2.0 * x)
2057 * plot [0:2] q2(x)
2058 * Particles > 2h away are excluded above. */
2059 dq = qfac2 * (2.0f * pow4f(2.0f - rij_h) - 4.0f * pow3f(2.0f - rij_h) * (1.0f + 2.0f * rij_h));
2060
2061 if (pfn->psys->part->flag & PART_SIZEMASS) {
2062 dq *= npa->size;
2063 }
2064
2065 pressureTerm = pressure / pow2f(pa->sphdensity) + npressure / pow2f(npa->sphdensity);
2066
2067 /* Note that 'minus' is removed, because vec = vecBA, not vecAB.
2068 * This applies to the viscosity calculation below, too. */
2069 madd_v3_v3fl(force, vec, pressureTerm * dq);
2070
2071 /* Viscosity */
2072 if (visc > 0.0f) {
2073 sub_v3_v3v3(dv, npa->prev_state.vel, pa->prev_state.vel);
2074 u = dot_v3v3(vec, dv);
2075 /* Apply parameters */
2076 u *= -dq * hinv * visc / (0.5f * npa->sphdensity + 0.5f * pa->sphdensity);
2077 madd_v3_v3fl(force, vec, u);
2078 }
2079 }
2080
2081 /* Artificial buoyancy force in negative gravity direction */
2082 if (fluid->buoyancy > 0.f && gravity) {
2083 madd_v3_v3fl(force, gravity, fluid->buoyancy * (pa->sphdensity - rest_density));
2084 }
2085
2086 if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF) {
2087 sph_particle_courant(sphdata, &pfr);
2088 }
2089 sphdata->pass++;
2090 }
2091
sphclassical_calc_dens(ParticleData * pa,float UNUSED (dfra),SPHData * sphdata)2092 static void sphclassical_calc_dens(ParticleData *pa, float UNUSED(dfra), SPHData *sphdata)
2093 {
2094 ParticleSystem **psys = sphdata->psys;
2095 SPHFluidSettings *fluid = psys[0]->part->fluid;
2096 /* 4.0 seems to be a pretty good value */
2097 float interaction_radius = fluid->radius *
2098 (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
2099 SPHRangeData pfr;
2100 float data[2];
2101
2102 data[0] = 0;
2103 data[1] = 0;
2104 pfr.data = data;
2105 pfr.h = interaction_radius * sphdata->hfac;
2106 pfr.pa = pa;
2107 pfr.mass = sphdata->mass;
2108
2109 sph_evaluate_func(
2110 NULL, psys, pa->state.co, &pfr, interaction_radius, sphclassical_density_accum_cb);
2111 pa->sphdensity = min_ff(max_ff(data[0], fluid->rest_density * 0.9f), fluid->rest_density * 1.1f);
2112 }
2113
psys_sph_init(ParticleSimulationData * sim,SPHData * sphdata)2114 void psys_sph_init(ParticleSimulationData *sim, SPHData *sphdata)
2115 {
2116 ParticleTarget *pt;
2117 int i;
2118
2119 BLI_buffer_field_init(&sphdata->new_springs, ParticleSpring);
2120
2121 /* Add other coupled particle systems. */
2122 sphdata->psys[0] = sim->psys;
2123 for (i = 1, pt = sim->psys->targets.first; i < 10; i++, pt = (pt ? pt->next : NULL)) {
2124 sphdata->psys[i] = pt ? psys_get_target_system(sim->ob, pt) : NULL;
2125 }
2126
2127 if (psys_uses_gravity(sim)) {
2128 sphdata->gravity = sim->scene->physics_settings.gravity;
2129 }
2130 else {
2131 sphdata->gravity = NULL;
2132 }
2133 sphdata->eh = sph_springhash_build(sim->psys);
2134
2135 /* These per-particle values should be overridden later, but just for
2136 * completeness we give them default values now. */
2137 sphdata->pa = NULL;
2138 sphdata->mass = 1.0f;
2139
2140 if (sim->psys->part->fluid->solver == SPH_SOLVER_DDR) {
2141 sphdata->force_cb = sph_force_cb;
2142 sphdata->density_cb = sph_density_accum_cb;
2143 sphdata->hfac = 1.0f;
2144 }
2145 else {
2146 /* SPH_SOLVER_CLASSICAL */
2147 sphdata->force_cb = sphclassical_force_cb;
2148 sphdata->density_cb = sphclassical_density_accum_cb;
2149 sphdata->hfac = 0.5f;
2150 }
2151 }
2152
psys_sph_flush_springs(SPHData * sphdata)2153 static void psys_sph_flush_springs(SPHData *sphdata)
2154 {
2155 for (int i = 0; i < sphdata->new_springs.count; i++) {
2156 /* sph_spring_add is not thread-safe. - z0r */
2157 sph_spring_add(sphdata->psys[0], &BLI_buffer_at(&sphdata->new_springs, ParticleSpring, i));
2158 }
2159
2160 BLI_buffer_field_free(&sphdata->new_springs);
2161 }
2162
psys_sph_finalize(SPHData * sphdata)2163 void psys_sph_finalize(SPHData *sphdata)
2164 {
2165 psys_sph_flush_springs(sphdata);
2166
2167 if (sphdata->eh) {
2168 BLI_edgehash_free(sphdata->eh, NULL);
2169 sphdata->eh = NULL;
2170 }
2171 }
2172
2173 /* Sample the density field at a point in space. */
psys_sph_density(BVHTree * tree,SPHData * sphdata,float co[3],float vars[2])2174 void psys_sph_density(BVHTree *tree, SPHData *sphdata, float co[3], float vars[2])
2175 {
2176 ParticleSystem **psys = sphdata->psys;
2177 SPHFluidSettings *fluid = psys[0]->part->fluid;
2178 /* 4.0 seems to be a pretty good value */
2179 float interaction_radius = fluid->radius *
2180 (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f);
2181 SPHRangeData pfr;
2182 float density[2];
2183
2184 density[0] = density[1] = 0.0f;
2185 pfr.data = density;
2186 pfr.h = interaction_radius * sphdata->hfac;
2187 pfr.mass = sphdata->mass;
2188
2189 sph_evaluate_func(tree, psys, co, &pfr, interaction_radius, sphdata->density_cb);
2190
2191 vars[0] = pfr.data[0];
2192 vars[1] = pfr.data[1];
2193 }
2194
sph_integrate(ParticleSimulationData * sim,ParticleData * pa,float dfra,SPHData * sphdata)2195 static void sph_integrate(ParticleSimulationData *sim,
2196 ParticleData *pa,
2197 float dfra,
2198 SPHData *sphdata)
2199 {
2200 ParticleSettings *part = sim->psys->part;
2201 // float timestep = psys_get_timestep(sim); // UNUSED
2202 float pa_mass = part->mass * (part->flag & PART_SIZEMASS ? pa->size : 1.f);
2203 float dtime = dfra * psys_get_timestep(sim);
2204 // int steps = 1; // UNUSED
2205 float effector_acceleration[3];
2206
2207 sphdata->pa = pa;
2208 sphdata->mass = pa_mass;
2209 sphdata->pass = 0;
2210 // sphdata.element_size and sphdata.flow are set in the callback.
2211
2212 /* restore previous state and treat gravity & effectors as external acceleration*/
2213 sub_v3_v3v3(effector_acceleration, pa->state.vel, pa->prev_state.vel);
2214 mul_v3_fl(effector_acceleration, 1.f / dtime);
2215
2216 copy_particle_key(&pa->state, &pa->prev_state, 0);
2217
2218 integrate_particle(part, pa, dtime, effector_acceleration, sphdata->force_cb, sphdata);
2219 }
2220
2221 /** \} */
2222
2223 /************************************************/
2224 /* Basic physics */
2225 /************************************************/
2226 typedef struct EfData {
2227 ParticleTexture ptex;
2228 ParticleSimulationData *sim;
2229 ParticleData *pa;
2230 } EfData;
basic_force_cb(void * efdata_v,ParticleKey * state,float * force,float * impulse)2231 static void basic_force_cb(void *efdata_v, ParticleKey *state, float *force, float *impulse)
2232 {
2233 EfData *efdata = (EfData *)efdata_v;
2234 ParticleSimulationData *sim = efdata->sim;
2235 ParticleSettings *part = sim->psys->part;
2236 ParticleData *pa = efdata->pa;
2237 EffectedPoint epoint;
2238 RNG *rng = sim->rng;
2239
2240 /* add effectors */
2241 pd_point_from_particle(efdata->sim, efdata->pa, state, &epoint);
2242 if (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR) {
2243 BKE_effectors_apply(sim->psys->effectors,
2244 sim->colliders,
2245 part->effector_weights,
2246 &epoint,
2247 force,
2248 NULL,
2249 impulse);
2250 }
2251
2252 mul_v3_fl(force, efdata->ptex.field);
2253 mul_v3_fl(impulse, efdata->ptex.field);
2254
2255 /* calculate air-particle interaction */
2256 if (part->dragfac != 0.0f) {
2257 madd_v3_v3fl(force, state->vel, -part->dragfac * pa->size * pa->size * len_v3(state->vel));
2258 }
2259
2260 /* brownian force */
2261 if (part->brownfac != 0.0f) {
2262 force[0] += (BLI_rng_get_float(rng) - 0.5f) * part->brownfac;
2263 force[1] += (BLI_rng_get_float(rng) - 0.5f) * part->brownfac;
2264 force[2] += (BLI_rng_get_float(rng) - 0.5f) * part->brownfac;
2265 }
2266
2267 if (part->flag & PART_ROT_DYN && epoint.ave) {
2268 copy_v3_v3(pa->state.ave, epoint.ave);
2269 }
2270 }
2271 /* gathers all forces that effect particles and calculates a new state for the particle */
basic_integrate(ParticleSimulationData * sim,int p,float dfra,float cfra)2272 static void basic_integrate(ParticleSimulationData *sim, int p, float dfra, float cfra)
2273 {
2274 ParticleSettings *part = sim->psys->part;
2275 ParticleData *pa = sim->psys->particles + p;
2276 ParticleKey tkey;
2277 float dtime = dfra * psys_get_timestep(sim), time;
2278 float *gravity = NULL, gr[3];
2279 EfData efdata;
2280
2281 psys_get_texture(sim, pa, &efdata.ptex, PAMAP_PHYSICS, cfra);
2282
2283 efdata.pa = pa;
2284 efdata.sim = sim;
2285
2286 /* add global acceleration (gravitation) */
2287 if (psys_uses_gravity(sim) &&
2288 /* normal gravity is too strong for hair so it's disabled by default */
2289 (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
2290 zero_v3(gr);
2291 madd_v3_v3fl(gr,
2292 sim->scene->physics_settings.gravity,
2293 part->effector_weights->global_gravity * efdata.ptex.gravity);
2294 gravity = gr;
2295 }
2296
2297 /* maintain angular velocity */
2298 copy_v3_v3(pa->state.ave, pa->prev_state.ave);
2299
2300 integrate_particle(part, pa, dtime, gravity, basic_force_cb, &efdata);
2301
2302 /* damp affects final velocity */
2303 if (part->dampfac != 0.f) {
2304 mul_v3_fl(pa->state.vel, 1.f - part->dampfac * efdata.ptex.damp * 25.f * dtime);
2305 }
2306
2307 // copy_v3_v3(pa->state.ave, states->ave);
2308
2309 /* finally we do guides */
2310 time = (cfra - pa->time) / pa->lifetime;
2311 CLAMP(time, 0.0f, 1.0f);
2312
2313 copy_v3_v3(tkey.co, pa->state.co);
2314 copy_v3_v3(tkey.vel, pa->state.vel);
2315 tkey.time = pa->state.time;
2316
2317 if (part->type != PART_HAIR) {
2318 if (do_guides(sim->depsgraph, sim->psys->part, sim->psys->effectors, &tkey, p, time)) {
2319 copy_v3_v3(pa->state.co, tkey.co);
2320 /* guides don't produce valid velocity */
2321 sub_v3_v3v3(pa->state.vel, tkey.co, pa->prev_state.co);
2322 mul_v3_fl(pa->state.vel, 1.0f / dtime);
2323 pa->state.time = tkey.time;
2324 }
2325 }
2326 }
basic_rotate(ParticleSettings * part,ParticleData * pa,float dfra,float timestep)2327 static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
2328 {
2329 float rotfac, rot1[4], rot2[4] = {1.0, 0.0, 0.0, 0.0}, dtime = dfra * timestep, extrotfac;
2330
2331 if ((part->flag & PART_ROTATIONS) == 0) {
2332 unit_qt(pa->state.rot);
2333 return;
2334 }
2335
2336 if (part->flag & PART_ROT_DYN) {
2337 extrotfac = len_v3(pa->state.ave);
2338 }
2339 else {
2340 extrotfac = 0.0f;
2341 }
2342
2343 if ((part->flag & PART_ROT_DYN) &&
2344 ELEM(part->avemode, PART_AVE_VELOCITY, PART_AVE_HORIZONTAL, PART_AVE_VERTICAL)) {
2345 float angle;
2346 float len1 = len_v3(pa->prev_state.vel);
2347 float len2 = len_v3(pa->state.vel);
2348 float vec[3];
2349
2350 if (len1 == 0.0f || len2 == 0.0f) {
2351 zero_v3(pa->state.ave);
2352 }
2353 else {
2354 cross_v3_v3v3(pa->state.ave, pa->prev_state.vel, pa->state.vel);
2355 normalize_v3(pa->state.ave);
2356 angle = dot_v3v3(pa->prev_state.vel, pa->state.vel) / (len1 * len2);
2357 mul_v3_fl(pa->state.ave, saacos(angle) / dtime);
2358 }
2359
2360 get_angular_velocity_vector(part->avemode, &pa->state, vec);
2361 axis_angle_to_quat(rot2, vec, dtime * part->avefac);
2362 }
2363
2364 rotfac = len_v3(pa->state.ave);
2365 if (rotfac == 0.0f || (part->flag & PART_ROT_DYN) == 0 || extrotfac == 0.0f) {
2366 unit_qt(rot1);
2367 }
2368 else {
2369 axis_angle_to_quat(rot1, pa->state.ave, rotfac * dtime);
2370 }
2371 mul_qt_qtqt(pa->state.rot, rot1, pa->prev_state.rot);
2372 mul_qt_qtqt(pa->state.rot, rot2, pa->state.rot);
2373
2374 /* keep rotation quat in good health */
2375 normalize_qt(pa->state.rot);
2376 }
2377
2378 /************************************************
2379 * Collisions
2380 *
2381 * The algorithm is roughly:
2382 * 1. Use a BVH tree to search for faces that a particle may collide with.
2383 * 2. Use Newton's method to find the exact time at which the collision occurs.
2384 * https://en.wikipedia.org/wiki/Newton's_method
2385 *
2386 ************************************************/
2387 #define COLLISION_MIN_RADIUS 0.001f
2388 #define COLLISION_MIN_DISTANCE 0.0001f
2389 #define COLLISION_ZERO 0.00001f
2390 #define COLLISION_INIT_STEP 0.00008f
2391 typedef float (*NRDistanceFunc)(float *p, float radius, ParticleCollisionElement *pce, float *nor);
nr_signed_distance_to_plane(float * p,float radius,ParticleCollisionElement * pce,float * nor)2392 static float nr_signed_distance_to_plane(float *p,
2393 float radius,
2394 ParticleCollisionElement *pce,
2395 float *nor)
2396 {
2397 float p0[3], e1[3], e2[3], d;
2398
2399 sub_v3_v3v3(e1, pce->x1, pce->x0);
2400 sub_v3_v3v3(e2, pce->x2, pce->x0);
2401 sub_v3_v3v3(p0, p, pce->x0);
2402
2403 cross_v3_v3v3(nor, e1, e2);
2404 normalize_v3(nor);
2405
2406 d = dot_v3v3(p0, nor);
2407
2408 if (pce->inv_nor == -1) {
2409 if (d < 0.f) {
2410 pce->inv_nor = 1;
2411 }
2412 else {
2413 pce->inv_nor = 0;
2414 }
2415 }
2416
2417 if (pce->inv_nor == 1) {
2418 negate_v3(nor);
2419 d = -d;
2420 }
2421
2422 return d - radius;
2423 }
nr_distance_to_edge(float * p,float radius,ParticleCollisionElement * pce,float * UNUSED (nor))2424 static float nr_distance_to_edge(float *p,
2425 float radius,
2426 ParticleCollisionElement *pce,
2427 float *UNUSED(nor))
2428 {
2429 float v0[3], v1[3], v2[3], c[3];
2430
2431 sub_v3_v3v3(v0, pce->x1, pce->x0);
2432 sub_v3_v3v3(v1, p, pce->x0);
2433 sub_v3_v3v3(v2, p, pce->x1);
2434
2435 cross_v3_v3v3(c, v1, v2);
2436
2437 return fabsf(len_v3(c) / len_v3(v0)) - radius;
2438 }
nr_distance_to_vert(float * p,float radius,ParticleCollisionElement * pce,float * UNUSED (nor))2439 static float nr_distance_to_vert(float *p,
2440 float radius,
2441 ParticleCollisionElement *pce,
2442 float *UNUSED(nor))
2443 {
2444 return len_v3v3(p, pce->x0) - radius;
2445 }
collision_interpolate_element(ParticleCollisionElement * pce,float t,float fac,ParticleCollision * col)2446 static void collision_interpolate_element(ParticleCollisionElement *pce,
2447 float t,
2448 float fac,
2449 ParticleCollision *col)
2450 {
2451 /* t is the current time for newton rhapson */
2452 /* fac is the starting factor for current collision iteration */
2453 /* The col->fac's are factors for the particle subframe step start
2454 * and end during collision modifier step. */
2455 float f = fac + t * (1.f - fac);
2456 float mul = col->fac1 + f * (col->fac2 - col->fac1);
2457 if (pce->tot > 0) {
2458 madd_v3_v3v3fl(pce->x0, pce->x[0], pce->v[0], mul);
2459
2460 if (pce->tot > 1) {
2461 madd_v3_v3v3fl(pce->x1, pce->x[1], pce->v[1], mul);
2462
2463 if (pce->tot > 2) {
2464 madd_v3_v3v3fl(pce->x2, pce->x[2], pce->v[2], mul);
2465 }
2466 }
2467 }
2468 }
collision_point_velocity(ParticleCollisionElement * pce)2469 static void collision_point_velocity(ParticleCollisionElement *pce)
2470 {
2471 float v[3];
2472
2473 copy_v3_v3(pce->vel, pce->v[0]);
2474
2475 if (pce->tot > 1) {
2476 sub_v3_v3v3(v, pce->v[1], pce->v[0]);
2477 madd_v3_v3fl(pce->vel, v, pce->uv[0]);
2478
2479 if (pce->tot > 2) {
2480 sub_v3_v3v3(v, pce->v[2], pce->v[0]);
2481 madd_v3_v3fl(pce->vel, v, pce->uv[1]);
2482 }
2483 }
2484 }
collision_point_distance_with_normal(float p[3],ParticleCollisionElement * pce,float fac,ParticleCollision * col,float * nor)2485 static float collision_point_distance_with_normal(
2486 float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *nor)
2487 {
2488 if (fac >= 0.f) {
2489 collision_interpolate_element(pce, 0.f, fac, col);
2490 }
2491
2492 switch (pce->tot) {
2493 case 1: {
2494 sub_v3_v3v3(nor, p, pce->x0);
2495 return normalize_v3(nor);
2496 }
2497 case 2: {
2498 float u, e[3], vec[3];
2499 sub_v3_v3v3(e, pce->x1, pce->x0);
2500 sub_v3_v3v3(vec, p, pce->x0);
2501 u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2502
2503 madd_v3_v3v3fl(nor, vec, e, -u);
2504 return normalize_v3(nor);
2505 }
2506 case 3:
2507 return nr_signed_distance_to_plane(p, 0.f, pce, nor);
2508 }
2509 return 0;
2510 }
collision_point_on_surface(const float p[3],ParticleCollisionElement * pce,float fac,ParticleCollision * col,float * co)2511 static void collision_point_on_surface(
2512 const float p[3], ParticleCollisionElement *pce, float fac, ParticleCollision *col, float *co)
2513 {
2514 collision_interpolate_element(pce, 0.f, fac, col);
2515
2516 switch (pce->tot) {
2517 case 1: {
2518 sub_v3_v3v3(co, p, pce->x0);
2519 normalize_v3(co);
2520 madd_v3_v3v3fl(co, pce->x0, co, col->radius);
2521 break;
2522 }
2523 case 2: {
2524 float u, e[3], vec[3], nor[3];
2525 sub_v3_v3v3(e, pce->x1, pce->x0);
2526 sub_v3_v3v3(vec, p, pce->x0);
2527 u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2528
2529 madd_v3_v3v3fl(nor, vec, e, -u);
2530 normalize_v3(nor);
2531
2532 madd_v3_v3v3fl(co, pce->x0, e, pce->uv[0]);
2533 madd_v3_v3fl(co, nor, col->radius);
2534 break;
2535 }
2536 case 3: {
2537 float p0[3], e1[3], e2[3], nor[3];
2538
2539 sub_v3_v3v3(e1, pce->x1, pce->x0);
2540 sub_v3_v3v3(e2, pce->x2, pce->x0);
2541 sub_v3_v3v3(p0, p, pce->x0);
2542
2543 cross_v3_v3v3(nor, e1, e2);
2544 normalize_v3(nor);
2545
2546 if (pce->inv_nor == 1) {
2547 negate_v3(nor);
2548 }
2549
2550 madd_v3_v3v3fl(co, pce->x0, nor, col->radius);
2551 madd_v3_v3fl(co, e1, pce->uv[0]);
2552 madd_v3_v3fl(co, e2, pce->uv[1]);
2553 break;
2554 }
2555 }
2556 }
2557 /* find first root in range [0-1] starting from 0 */
collision_newton_rhapson(ParticleCollision * col,float radius,ParticleCollisionElement * pce,NRDistanceFunc distance_func)2558 static float collision_newton_rhapson(ParticleCollision *col,
2559 float radius,
2560 ParticleCollisionElement *pce,
2561 NRDistanceFunc distance_func)
2562 {
2563 float t0, t1, dt_init, d0, d1, dd, n[3];
2564 int iter;
2565
2566 pce->inv_nor = -1;
2567
2568 if (col->inv_total_time > 0.0f) {
2569 /* Initial step size should be small, but not too small or floating point
2570 * precision errors will appear. - z0r */
2571 dt_init = COLLISION_INIT_STEP * col->inv_total_time;
2572 }
2573 else {
2574 dt_init = 0.001f;
2575 }
2576
2577 /* start from the beginning */
2578 t0 = 0.f;
2579 collision_interpolate_element(pce, t0, col->f, col);
2580 d0 = distance_func(col->co1, radius, pce, n);
2581 t1 = dt_init;
2582 d1 = 0.f;
2583
2584 for (iter = 0; iter < 10; iter++) { //, itersum++) {
2585 /* get current location */
2586 collision_interpolate_element(pce, t1, col->f, col);
2587 interp_v3_v3v3(pce->p, col->co1, col->co2, t1);
2588
2589 d1 = distance_func(pce->p, radius, pce, n);
2590
2591 /* particle already inside face, so report collision */
2592 if (iter == 0 && d0 < 0.f && d0 > -radius) {
2593 copy_v3_v3(pce->p, col->co1);
2594 copy_v3_v3(pce->nor, n);
2595 pce->inside = 1;
2596 return 0.f;
2597 }
2598
2599 /* Zero gradient (no movement relative to element). Can't step from
2600 * here. */
2601 if (d1 == d0) {
2602 /* If first iteration, try from other end where the gradient may be
2603 * greater. Note: code duplicated below. */
2604 if (iter == 0) {
2605 t0 = 1.f;
2606 collision_interpolate_element(pce, t0, col->f, col);
2607 d0 = distance_func(col->co2, radius, pce, n);
2608 t1 = 1.0f - dt_init;
2609 d1 = 0.f;
2610 continue;
2611 }
2612
2613 return -1.f;
2614 }
2615
2616 dd = (t1 - t0) / (d1 - d0);
2617
2618 t0 = t1;
2619 d0 = d1;
2620
2621 t1 -= d1 * dd;
2622
2623 /* Particle moving away from plane could also mean a strangely rotating
2624 * face, so check from end. Note: code duplicated above. */
2625 if (iter == 0 && t1 < 0.f) {
2626 t0 = 1.f;
2627 collision_interpolate_element(pce, t0, col->f, col);
2628 d0 = distance_func(col->co2, radius, pce, n);
2629 t1 = 1.0f - dt_init;
2630 d1 = 0.f;
2631 continue;
2632 }
2633 if (iter == 1 && (t1 < -COLLISION_ZERO || t1 > 1.f)) {
2634 return -1.f;
2635 }
2636
2637 if (d1 <= COLLISION_ZERO && d1 >= -COLLISION_ZERO) {
2638 if (t1 >= -COLLISION_ZERO && t1 <= 1.f) {
2639 if (distance_func == nr_signed_distance_to_plane) {
2640 copy_v3_v3(pce->nor, n);
2641 }
2642
2643 CLAMP(t1, 0.f, 1.f);
2644
2645 return t1;
2646 }
2647
2648 return -1.f;
2649 }
2650 }
2651 return -1.0;
2652 }
collision_sphere_to_tri(ParticleCollision * col,float radius,ParticleCollisionElement * pce,float * t)2653 static int collision_sphere_to_tri(ParticleCollision *col,
2654 float radius,
2655 ParticleCollisionElement *pce,
2656 float *t)
2657 {
2658 ParticleCollisionElement *result = &col->pce;
2659 float ct, u, v;
2660
2661 pce->inv_nor = -1;
2662 pce->inside = 0;
2663
2664 ct = collision_newton_rhapson(col, radius, pce, nr_signed_distance_to_plane);
2665
2666 if (ct >= 0.f && ct < *t && (result->inside == 0 || pce->inside == 1)) {
2667 float e1[3], e2[3], p0[3];
2668 float e1e1, e1e2, e1p0, e2e2, e2p0, inv;
2669
2670 sub_v3_v3v3(e1, pce->x1, pce->x0);
2671 sub_v3_v3v3(e2, pce->x2, pce->x0);
2672 /* XXX: add radius correction here? */
2673 sub_v3_v3v3(p0, pce->p, pce->x0);
2674
2675 e1e1 = dot_v3v3(e1, e1);
2676 e1e2 = dot_v3v3(e1, e2);
2677 e1p0 = dot_v3v3(e1, p0);
2678 e2e2 = dot_v3v3(e2, e2);
2679 e2p0 = dot_v3v3(e2, p0);
2680
2681 inv = 1.f / (e1e1 * e2e2 - e1e2 * e1e2);
2682 u = (e2e2 * e1p0 - e1e2 * e2p0) * inv;
2683 v = (e1e1 * e2p0 - e1e2 * e1p0) * inv;
2684
2685 if (u >= 0.f && u <= 1.f && v >= 0.f && u + v <= 1.f) {
2686 *result = *pce;
2687
2688 /* normal already calculated in pce */
2689
2690 result->uv[0] = u;
2691 result->uv[1] = v;
2692
2693 *t = ct;
2694 return 1;
2695 }
2696 }
2697 return 0;
2698 }
collision_sphere_to_edges(ParticleCollision * col,float radius,ParticleCollisionElement * pce,float * t)2699 static int collision_sphere_to_edges(ParticleCollision *col,
2700 float radius,
2701 ParticleCollisionElement *pce,
2702 float *t)
2703 {
2704 ParticleCollisionElement edge[3], *cur = NULL, *hit = NULL;
2705 ParticleCollisionElement *result = &col->pce;
2706
2707 float ct;
2708 int i;
2709
2710 for (i = 0; i < 3; i++) {
2711 cur = edge + i;
2712 cur->x[0] = pce->x[i];
2713 cur->x[1] = pce->x[(i + 1) % 3];
2714 cur->v[0] = pce->v[i];
2715 cur->v[1] = pce->v[(i + 1) % 3];
2716 cur->tot = 2;
2717 cur->inside = 0;
2718
2719 ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_edge);
2720
2721 if (ct >= 0.f && ct < *t) {
2722 float u, e[3], vec[3];
2723
2724 sub_v3_v3v3(e, cur->x1, cur->x0);
2725 sub_v3_v3v3(vec, cur->p, cur->x0);
2726 u = dot_v3v3(vec, e) / dot_v3v3(e, e);
2727
2728 if (u < 0.f || u > 1.f) {
2729 break;
2730 }
2731
2732 *result = *cur;
2733
2734 madd_v3_v3v3fl(result->nor, vec, e, -u);
2735 normalize_v3(result->nor);
2736
2737 result->uv[0] = u;
2738
2739 hit = cur;
2740 *t = ct;
2741 }
2742 }
2743
2744 return hit != NULL;
2745 }
collision_sphere_to_verts(ParticleCollision * col,float radius,ParticleCollisionElement * pce,float * t)2746 static int collision_sphere_to_verts(ParticleCollision *col,
2747 float radius,
2748 ParticleCollisionElement *pce,
2749 float *t)
2750 {
2751 ParticleCollisionElement vert[3], *cur = NULL, *hit = NULL;
2752 ParticleCollisionElement *result = &col->pce;
2753
2754 float ct;
2755 int i;
2756
2757 for (i = 0; i < 3; i++) {
2758 cur = vert + i;
2759 cur->x[0] = pce->x[i];
2760 cur->v[0] = pce->v[i];
2761 cur->tot = 1;
2762 cur->inside = 0;
2763
2764 ct = collision_newton_rhapson(col, radius, cur, nr_distance_to_vert);
2765
2766 if (ct >= 0.f && ct < *t) {
2767 *result = *cur;
2768
2769 sub_v3_v3v3(result->nor, cur->p, cur->x0);
2770 normalize_v3(result->nor);
2771
2772 hit = cur;
2773 *t = ct;
2774 }
2775 }
2776
2777 return hit != NULL;
2778 }
2779 /* Callback for BVHTree near test */
BKE_psys_collision_neartest_cb(void * userdata,int index,const BVHTreeRay * ray,BVHTreeRayHit * hit)2780 void BKE_psys_collision_neartest_cb(void *userdata,
2781 int index,
2782 const BVHTreeRay *ray,
2783 BVHTreeRayHit *hit)
2784 {
2785 ParticleCollision *col = (ParticleCollision *)userdata;
2786 ParticleCollisionElement pce;
2787 const MVertTri *vt = &col->md->tri[index];
2788 MVert *x = col->md->x;
2789 MVert *v = col->md->current_v;
2790 float t = hit->dist / col->original_ray_length;
2791 int collision = 0;
2792
2793 pce.x[0] = x[vt->tri[0]].co;
2794 pce.x[1] = x[vt->tri[1]].co;
2795 pce.x[2] = x[vt->tri[2]].co;
2796
2797 pce.v[0] = v[vt->tri[0]].co;
2798 pce.v[1] = v[vt->tri[1]].co;
2799 pce.v[2] = v[vt->tri[2]].co;
2800
2801 pce.tot = 3;
2802 pce.inside = 0;
2803 pce.index = index;
2804
2805 collision = collision_sphere_to_tri(col, ray->radius, &pce, &t);
2806 if (col->pce.inside == 0) {
2807 collision += collision_sphere_to_edges(col, ray->radius, &pce, &t);
2808 collision += collision_sphere_to_verts(col, ray->radius, &pce, &t);
2809 }
2810
2811 if (collision) {
2812 hit->dist = col->original_ray_length * t;
2813 hit->index = index;
2814
2815 collision_point_velocity(&col->pce);
2816
2817 col->hit = col->current;
2818 }
2819 }
collision_detect(ParticleData * pa,ParticleCollision * col,BVHTreeRayHit * hit,ListBase * colliders)2820 static int collision_detect(ParticleData *pa,
2821 ParticleCollision *col,
2822 BVHTreeRayHit *hit,
2823 ListBase *colliders)
2824 {
2825 const int raycast_flag = BVH_RAYCAST_DEFAULT & ~(BVH_RAYCAST_WATERTIGHT);
2826 ColliderCache *coll;
2827 float ray_dir[3];
2828
2829 if (BLI_listbase_is_empty(colliders)) {
2830 return 0;
2831 }
2832
2833 sub_v3_v3v3(ray_dir, col->co2, col->co1);
2834 hit->index = -1;
2835 hit->dist = col->original_ray_length = normalize_v3(ray_dir);
2836 col->pce.inside = 0;
2837
2838 /* even if particle is stationary we want to check for moving colliders */
2839 /* if hit.dist is zero the bvhtree_ray_cast will just ignore everything */
2840 if (hit->dist == 0.0f) {
2841 hit->dist = col->original_ray_length = 0.000001f;
2842 }
2843
2844 for (coll = colliders->first; coll; coll = coll->next) {
2845 /* for boids: don't check with current ground object; also skip if permeated */
2846 bool skip = false;
2847
2848 for (int i = 0; i < col->skip_count; i++) {
2849 if (coll->ob == col->skip[i]) {
2850 skip = true;
2851 break;
2852 }
2853 }
2854
2855 if (skip) {
2856 continue;
2857 }
2858
2859 /* particles should not collide with emitter at birth */
2860 if (coll->ob == col->emitter && pa->time < col->cfra && pa->time >= col->old_cfra) {
2861 continue;
2862 }
2863
2864 col->current = coll->ob;
2865 col->md = coll->collmd;
2866 col->fac1 = (col->old_cfra - coll->collmd->time_x) /
2867 (coll->collmd->time_xnew - coll->collmd->time_x);
2868 col->fac2 = (col->cfra - coll->collmd->time_x) /
2869 (coll->collmd->time_xnew - coll->collmd->time_x);
2870
2871 if (col->md && col->md->bvhtree) {
2872 BLI_bvhtree_ray_cast_ex(col->md->bvhtree,
2873 col->co1,
2874 ray_dir,
2875 col->radius,
2876 hit,
2877 BKE_psys_collision_neartest_cb,
2878 col,
2879 raycast_flag);
2880 }
2881 }
2882
2883 return hit->index >= 0;
2884 }
collision_response(ParticleSimulationData * sim,ParticleData * pa,ParticleCollision * col,BVHTreeRayHit * hit,int kill,int dynamic_rotation)2885 static int collision_response(ParticleSimulationData *sim,
2886 ParticleData *pa,
2887 ParticleCollision *col,
2888 BVHTreeRayHit *hit,
2889 int kill,
2890 int dynamic_rotation)
2891 {
2892 ParticleCollisionElement *pce = &col->pce;
2893 PartDeflect *pd = col->hit->pd;
2894 RNG *rng = sim->rng;
2895 /* point of collision */
2896 float co[3];
2897 /* location factor of collision between this iteration */
2898 float x = hit->dist / col->original_ray_length;
2899 /* time factor of collision between timestep */
2900 float f = col->f + x * (1.0f - col->f);
2901 /* time since previous collision (in seconds) */
2902 float dt1 = (f - col->f) * col->total_time;
2903 /* time left after collision (in seconds) */
2904 float dt2 = (1.0f - f) * col->total_time;
2905 /* did particle pass through the collision surface? */
2906 int through = (BLI_rng_get_float(rng) < pd->pdef_perm) ? 1 : 0;
2907
2908 /* calculate exact collision location */
2909 interp_v3_v3v3(co, col->co1, col->co2, x);
2910
2911 /* particle dies in collision */
2912 if (through == 0 && (kill || pd->flag & PDEFLE_KILL_PART)) {
2913 pa->alive = PARS_DYING;
2914 pa->dietime = col->old_cfra + (col->cfra - col->old_cfra) * f;
2915
2916 copy_v3_v3(pa->state.co, co);
2917 interp_v3_v3v3(pa->state.vel, pa->prev_state.vel, pa->state.vel, f);
2918 interp_qt_qtqt(pa->state.rot, pa->prev_state.rot, pa->state.rot, f);
2919 interp_v3_v3v3(pa->state.ave, pa->prev_state.ave, pa->state.ave, f);
2920
2921 /* particle is dead so we don't need to calculate further */
2922 return 0;
2923 }
2924 /* figure out velocity and other data after collision */
2925
2926 /* velocity directly before collision to be modified into velocity directly after collision */
2927 float v0[3];
2928 /* normal component of v0 */
2929 float v0_nor[3];
2930 /* tangential component of v0 */
2931 float v0_tan[3];
2932 /* tangential component of collision surface velocity */
2933 float vc_tan[3];
2934 float v0_dot, vc_dot;
2935 float damp = pd->pdef_damp + pd->pdef_rdamp * 2 * (BLI_rng_get_float(rng) - 0.5f);
2936 float frict = pd->pdef_frict + pd->pdef_rfrict * 2 * (BLI_rng_get_float(rng) - 0.5f);
2937 float distance, nor[3], dot;
2938
2939 CLAMP(damp, 0.0f, 1.0f);
2940 CLAMP(frict, 0.0f, 1.0f);
2941
2942 /* get exact velocity right before collision */
2943 madd_v3_v3v3fl(v0, col->ve1, col->acc, dt1);
2944
2945 /* Convert collider velocity from `1/frame_step` to `1/s` TODO:
2946 * here we assume 1 frame step for collision modifier. */
2947 mul_v3_fl(pce->vel, col->inv_timestep);
2948
2949 /* calculate tangential particle velocity */
2950 v0_dot = dot_v3v3(pce->nor, v0);
2951 madd_v3_v3v3fl(v0_tan, v0, pce->nor, -v0_dot);
2952
2953 /* calculate tangential collider velocity */
2954 vc_dot = dot_v3v3(pce->nor, pce->vel);
2955 madd_v3_v3v3fl(vc_tan, pce->vel, pce->nor, -vc_dot);
2956
2957 /* handle friction effects (tangential and angular velocity) */
2958 if (frict > 0.0f) {
2959 /* angular <-> linear velocity */
2960 if (dynamic_rotation) {
2961 float vr_tan[3], v1_tan[3], ave[3];
2962
2963 /* linear velocity of particle surface */
2964 cross_v3_v3v3(vr_tan, pce->nor, pa->state.ave);
2965 mul_v3_fl(vr_tan, pa->size);
2966
2967 /* change to coordinates that move with the collision plane */
2968 sub_v3_v3v3(v1_tan, v0_tan, vc_tan);
2969
2970 /* The resulting velocity is a weighted average of particle cm & surface
2971 * velocity. This weight (related to particle's moment of inertia) could
2972 * be made a parameter for angular <-> linear conversion.
2973 */
2974 madd_v3_v3fl(v1_tan, vr_tan, -0.4);
2975 mul_v3_fl(v1_tan, 1.0f / 1.4f); /* 1/(1+0.4) */
2976
2977 /* rolling friction is around 0.01 of sliding friction
2978 * (could be made a parameter) */
2979 mul_v3_fl(v1_tan, 1.0f - 0.01f * frict);
2980
2981 /* surface_velocity is opposite to cm velocity */
2982 negate_v3_v3(vr_tan, v1_tan);
2983
2984 /* get back to global coordinates */
2985 add_v3_v3(v1_tan, vc_tan);
2986
2987 /* convert to angular velocity*/
2988 cross_v3_v3v3(ave, vr_tan, pce->nor);
2989 mul_v3_fl(ave, 1.0f / MAX2(pa->size, 0.001f));
2990
2991 /* only friction will cause change in linear & angular velocity */
2992 interp_v3_v3v3(pa->state.ave, pa->state.ave, ave, frict);
2993 interp_v3_v3v3(v0_tan, v0_tan, v1_tan, frict);
2994 }
2995 else {
2996 /* just basic friction (unphysical due to the friction model used in Blender) */
2997 interp_v3_v3v3(v0_tan, v0_tan, vc_tan, frict);
2998 }
2999 }
3000
3001 /* Stickiness was possibly added before,
3002 * so cancel that before calculating new normal velocity.
3003 * Otherwise particles go flying out of the surface
3004 * because of high reversed sticky velocity. */
3005 if (v0_dot < 0.0f) {
3006 v0_dot += pd->pdef_stickness;
3007 if (v0_dot > 0.0f) {
3008 v0_dot = 0.0f;
3009 }
3010 }
3011
3012 /* damping and flipping of velocity around normal */
3013 v0_dot *= 1.0f - damp;
3014 vc_dot *= through ? damp : 1.0f;
3015
3016 /* calculate normal particle velocity */
3017 /* special case for object hitting the particle from behind */
3018 if (through == 0 && ((vc_dot > 0.0f && v0_dot > 0.0f && vc_dot > v0_dot) ||
3019 (vc_dot < 0.0f && v0_dot < 0.0f && vc_dot < v0_dot))) {
3020 mul_v3_v3fl(v0_nor, pce->nor, vc_dot);
3021 }
3022 else if (v0_dot > 0.f) {
3023 mul_v3_v3fl(v0_nor, pce->nor, vc_dot + v0_dot);
3024 }
3025 else {
3026 mul_v3_v3fl(v0_nor, pce->nor, vc_dot + (through ? 1.0f : -1.0f) * v0_dot);
3027 }
3028
3029 /* combine components together again */
3030 add_v3_v3v3(v0, v0_nor, v0_tan);
3031
3032 if (col->boid) {
3033 /* keep boids above ground */
3034 BoidParticle *bpa = pa->boid;
3035 if (bpa->data.mode == eBoidMode_OnLand || co[2] <= col->boid_z) {
3036 co[2] = col->boid_z;
3037 v0[2] = 0.0f;
3038 }
3039 }
3040
3041 /* re-apply acceleration to final location and velocity */
3042 madd_v3_v3v3fl(pa->state.co, co, v0, dt2);
3043 madd_v3_v3fl(pa->state.co, col->acc, 0.5f * dt2 * dt2);
3044 madd_v3_v3v3fl(pa->state.vel, v0, col->acc, dt2);
3045
3046 /* make sure particle stays on the right side of the surface */
3047 if (!through) {
3048 distance = collision_point_distance_with_normal(co, pce, -1.f, col, nor);
3049
3050 if (distance < col->radius + COLLISION_MIN_DISTANCE) {
3051 madd_v3_v3fl(co, nor, col->radius + COLLISION_MIN_DISTANCE - distance);
3052 }
3053
3054 dot = dot_v3v3(nor, v0);
3055 if (dot < 0.f) {
3056 madd_v3_v3fl(v0, nor, -dot);
3057 }
3058
3059 distance = collision_point_distance_with_normal(pa->state.co, pce, 1.f, col, nor);
3060
3061 if (distance < col->radius + COLLISION_MIN_DISTANCE) {
3062 madd_v3_v3fl(pa->state.co, nor, col->radius + COLLISION_MIN_DISTANCE - distance);
3063 }
3064
3065 dot = dot_v3v3(nor, pa->state.vel);
3066 if (dot < 0.f) {
3067 madd_v3_v3fl(pa->state.vel, nor, -dot);
3068 }
3069 }
3070
3071 /* add stickiness to surface */
3072 madd_v3_v3fl(pa->state.vel, pce->nor, -pd->pdef_stickness);
3073
3074 /* set coordinates for next iteration */
3075 copy_v3_v3(col->co1, co);
3076 copy_v3_v3(col->co2, pa->state.co);
3077
3078 copy_v3_v3(col->ve1, v0);
3079 copy_v3_v3(col->ve2, pa->state.vel);
3080
3081 col->f = f;
3082
3083 /* if permeability random roll succeeded, disable collider for this sim step */
3084 if (through) {
3085 col->skip[col->skip_count++] = col->hit;
3086 }
3087
3088 return 1;
3089 }
collision_fail(ParticleData * pa,ParticleCollision * col)3090 static void collision_fail(ParticleData *pa, ParticleCollision *col)
3091 {
3092 /* final chance to prevent total failure, so stick to the surface and hope for the best */
3093 collision_point_on_surface(col->co1, &col->pce, 1.f, col, pa->state.co);
3094
3095 copy_v3_v3(pa->state.vel, col->pce.vel);
3096 mul_v3_fl(pa->state.vel, col->inv_timestep);
3097
3098 /* printf("max iterations\n"); */
3099 }
3100
3101 /* Particle - Mesh collision detection and response
3102 * Features:
3103 * -friction and damping
3104 * -angular momentum <-> linear momentum
3105 * -high accuracy by re-applying particle acceleration after collision
3106 * -handles moving, rotating and deforming meshes
3107 * -uses Newton-Rhapson iteration to find the collisions
3108 * -handles spherical particles and (nearly) point like particles
3109 */
collision_check(ParticleSimulationData * sim,int p,float dfra,float cfra)3110 static void collision_check(ParticleSimulationData *sim, int p, float dfra, float cfra)
3111 {
3112 ParticleSettings *part = sim->psys->part;
3113 ParticleData *pa = sim->psys->particles + p;
3114 ParticleCollision col;
3115 BVHTreeRayHit hit;
3116 int collision_count = 0;
3117
3118 float timestep = psys_get_timestep(sim);
3119
3120 memset(&col, 0, sizeof(ParticleCollision));
3121
3122 col.total_time = timestep * dfra;
3123 col.inv_total_time = 1.0f / col.total_time;
3124 col.inv_timestep = 1.0f / timestep;
3125
3126 col.cfra = cfra;
3127 col.old_cfra = sim->psys->cfra;
3128
3129 /* get acceleration (from gravity, forcefields etc. to be re-applied in collision response) */
3130 sub_v3_v3v3(col.acc, pa->state.vel, pa->prev_state.vel);
3131 mul_v3_fl(col.acc, 1.f / col.total_time);
3132
3133 /* set values for first iteration */
3134 copy_v3_v3(col.co1, pa->prev_state.co);
3135 copy_v3_v3(col.co2, pa->state.co);
3136 copy_v3_v3(col.ve1, pa->prev_state.vel);
3137 copy_v3_v3(col.ve2, pa->state.vel);
3138 col.f = 0.0f;
3139
3140 col.radius = ((part->flag & PART_SIZE_DEFL) || (part->phystype == PART_PHYS_BOIDS)) ?
3141 pa->size :
3142 COLLISION_MIN_RADIUS;
3143
3144 /* override for boids */
3145 if (part->phystype == PART_PHYS_BOIDS && part->boids->options & BOID_ALLOW_LAND) {
3146 col.boid = 1;
3147 col.boid_z = pa->state.co[2];
3148 col.skip[col.skip_count++] = pa->boid->ground;
3149 }
3150
3151 /* 10 iterations to catch multiple collisions */
3152 while (collision_count < PARTICLE_COLLISION_MAX_COLLISIONS) {
3153 if (collision_detect(pa, &col, &hit, sim->colliders)) {
3154
3155 collision_count++;
3156
3157 if (collision_count == PARTICLE_COLLISION_MAX_COLLISIONS) {
3158 collision_fail(pa, &col);
3159 }
3160 else if (collision_response(
3161 sim, pa, &col, &hit, part->flag & PART_DIE_ON_COL, part->flag & PART_ROT_DYN) ==
3162 0) {
3163 return;
3164 }
3165 }
3166 else {
3167 return;
3168 }
3169 }
3170 }
3171 /************************************************/
3172 /* Hair */
3173 /************************************************/
3174 /* check if path cache or children need updating and do it if needed */
psys_update_path_cache(ParticleSimulationData * sim,float cfra,const bool use_render_params)3175 static void psys_update_path_cache(ParticleSimulationData *sim,
3176 float cfra,
3177 const bool use_render_params)
3178 {
3179 ParticleSystem *psys = sim->psys;
3180 ParticleSettings *part = psys->part;
3181 ParticleEditSettings *pset = &sim->scene->toolsettings->particle;
3182 int distr = 0, alloc = 0, skip = 0;
3183
3184 if ((psys->part->childtype &&
3185 psys->totchild != psys_get_tot_child(sim->scene, psys, use_render_params)) ||
3186 psys->recalc & ID_RECALC_PSYS_RESET) {
3187 alloc = 1;
3188 }
3189
3190 if (alloc || psys->recalc & ID_RECALC_PSYS_CHILD ||
3191 (psys->vgroup[PSYS_VG_DENSITY] && (sim->ob && sim->ob->mode & OB_MODE_WEIGHT_PAINT))) {
3192 distr = 1;
3193 }
3194
3195 if (distr) {
3196 if (alloc) {
3197 realloc_particles(sim, sim->psys->totpart);
3198 }
3199
3200 if (psys_get_tot_child(sim->scene, psys, use_render_params)) {
3201 /* don't generate children while computing the hair keys */
3202 if (!(psys->part->type == PART_HAIR) || (psys->flag & PSYS_HAIR_DONE)) {
3203 distribute_particles(sim, PART_FROM_CHILD);
3204
3205 if (part->childtype == PART_CHILD_FACES && part->parents != 0.0f) {
3206 psys_find_parents(sim, use_render_params);
3207 }
3208 }
3209 }
3210 else {
3211 psys_free_children(psys);
3212 }
3213 }
3214
3215 if ((part->type == PART_HAIR || psys->flag & PSYS_KEYED ||
3216 psys->pointcache->flag & PTCACHE_BAKED) == 0) {
3217 skip = 1; /* only hair, keyed and baked stuff can have paths */
3218 }
3219 else if (part->ren_as != PART_DRAW_PATH &&
3220 !(part->type == PART_HAIR && ELEM(part->ren_as, PART_DRAW_OB, PART_DRAW_GR))) {
3221 skip = 1; /* particle visualization must be set as path */
3222 }
3223 else if (DEG_get_mode(sim->depsgraph) != DAG_EVAL_RENDER) {
3224 if (part->draw_as != PART_DRAW_REND) {
3225 skip = 1; /* draw visualization */
3226 }
3227 else if (psys->pointcache->flag & PTCACHE_BAKING) {
3228 skip = 1; /* no need to cache paths while baking dynamics */
3229 }
3230 else if (psys_in_edit_mode(sim->depsgraph, psys)) {
3231 if ((pset->flag & PE_DRAW_PART) == 0) {
3232 skip = 1;
3233 }
3234 else if (part->childtype == 0 &&
3235 (psys->flag & PSYS_HAIR_DYNAMICS && psys->pointcache->flag & PTCACHE_BAKED) == 0) {
3236 skip = 1; /* in edit mode paths are needed for child particles and dynamic hair */
3237 }
3238 }
3239 }
3240
3241 if (!skip) {
3242 psys_cache_paths(sim, cfra, use_render_params);
3243
3244 /* for render, child particle paths are computed on the fly */
3245 if (part->childtype) {
3246 if (!psys->totchild) {
3247 skip = 1;
3248 }
3249 else if (psys->part->type == PART_HAIR && (psys->flag & PSYS_HAIR_DONE) == 0) {
3250 skip = 1;
3251 }
3252
3253 if (!skip) {
3254 psys_cache_child_paths(sim, cfra, 0, use_render_params);
3255 }
3256 }
3257 }
3258 else if (psys->pathcache) {
3259 psys_free_path_cache(psys, NULL);
3260 }
3261 }
3262
psys_hair_use_simulation(ParticleData * pa,float max_length)3263 static bool psys_hair_use_simulation(ParticleData *pa, float max_length)
3264 {
3265 /* Minimum segment length relative to average length.
3266 * Hairs with segments below this length will be excluded from the simulation,
3267 * because otherwise the solver will become unstable.
3268 * The hair system should always make sure the hair segments have reasonable length ratios,
3269 * but this can happen in old files when e.g. cutting hair.
3270 */
3271 const float min_length = 0.1f * max_length;
3272
3273 HairKey *key;
3274 int k;
3275
3276 if (pa->totkey < 2) {
3277 return false;
3278 }
3279
3280 for (k = 1, key = pa->hair + 1; k < pa->totkey; k++, key++) {
3281 float length = len_v3v3(key->co, (key - 1)->co);
3282 if (length < min_length) {
3283 return false;
3284 }
3285 }
3286
3287 return true;
3288 }
3289
hair_set_pinning(MDeformVert * dvert,float weight)3290 static MDeformVert *hair_set_pinning(MDeformVert *dvert, float weight)
3291 {
3292 if (dvert) {
3293 if (!dvert->totweight) {
3294 dvert->dw = MEM_callocN(sizeof(MDeformWeight), "deformWeight");
3295 dvert->totweight = 1;
3296 }
3297
3298 dvert->dw->weight = weight;
3299 dvert++;
3300 }
3301 return dvert;
3302 }
3303
hair_create_input_mesh(ParticleSimulationData * sim,int totpoint,int totedge,Mesh ** r_mesh,ClothHairData ** r_hairdata)3304 static void hair_create_input_mesh(ParticleSimulationData *sim,
3305 int totpoint,
3306 int totedge,
3307 Mesh **r_mesh,
3308 ClothHairData **r_hairdata)
3309 {
3310 ParticleSystem *psys = sim->psys;
3311 ParticleSettings *part = psys->part;
3312 Mesh *mesh;
3313 ClothHairData *hairdata;
3314 MVert *mvert;
3315 MEdge *medge;
3316 MDeformVert *dvert;
3317 HairKey *key;
3318 PARTICLE_P;
3319 int k, hair_index;
3320 float hairmat[4][4];
3321 float max_length;
3322 float hair_radius;
3323
3324 mesh = *r_mesh;
3325 if (!mesh) {
3326 *r_mesh = mesh = BKE_mesh_new_nomain(totpoint, totedge, 0, 0, 0);
3327 CustomData_add_layer(&mesh->vdata, CD_MDEFORMVERT, CD_CALLOC, NULL, mesh->totvert);
3328 BKE_mesh_update_customdata_pointers(mesh, false);
3329 }
3330 mvert = mesh->mvert;
3331 medge = mesh->medge;
3332 dvert = mesh->dvert;
3333
3334 hairdata = *r_hairdata;
3335 if (!hairdata) {
3336 *r_hairdata = hairdata = MEM_mallocN(sizeof(ClothHairData) * totpoint, "hair data");
3337 }
3338
3339 /* calculate maximum segment length */
3340 max_length = 0.0f;
3341 LOOP_PARTICLES
3342 {
3343 if (!(pa->flag & PARS_UNEXIST)) {
3344 for (k = 1, key = pa->hair + 1; k < pa->totkey; k++, key++) {
3345 float length = len_v3v3(key->co, (key - 1)->co);
3346 if (max_length < length) {
3347 max_length = length;
3348 }
3349 }
3350 }
3351 }
3352
3353 psys->clmd->sim_parms->vgroup_mass = 1;
3354
3355 /* XXX placeholder for more flexible future hair settings */
3356 hair_radius = part->size;
3357
3358 /* make vgroup for pin roots etc.. */
3359 hair_index = 1;
3360 LOOP_PARTICLES
3361 {
3362 if (!(pa->flag & PARS_UNEXIST)) {
3363 float root_mat[4][4];
3364 float bending_stiffness;
3365 bool use_hair;
3366
3367 pa->hair_index = hair_index;
3368 use_hair = psys_hair_use_simulation(pa, max_length);
3369
3370 psys_mat_hair_to_object(sim->ob, sim->psmd->mesh_final, psys->part->from, pa, hairmat);
3371 mul_m4_m4m4(root_mat, sim->ob->obmat, hairmat);
3372 normalize_m4(root_mat);
3373
3374 bending_stiffness = CLAMPIS(
3375 1.0f - part->bending_random * psys_frand(psys, p + 666), 0.0f, 1.0f);
3376
3377 for (k = 0, key = pa->hair; k < pa->totkey; k++, key++) {
3378 ClothHairData *hair;
3379 float *co, *co_next;
3380
3381 co = key->co;
3382 co_next = (key + 1)->co;
3383
3384 /* create fake root before actual root to resist bending */
3385 if (k == 0) {
3386 hair = &psys->clmd->hairdata[pa->hair_index - 1];
3387 copy_v3_v3(hair->loc, root_mat[3]);
3388 copy_m3_m4(hair->rot, root_mat);
3389
3390 hair->radius = hair_radius;
3391 hair->bending_stiffness = bending_stiffness;
3392
3393 add_v3_v3v3(mvert->co, co, co);
3394 sub_v3_v3(mvert->co, co_next);
3395 mul_m4_v3(hairmat, mvert->co);
3396
3397 medge->v1 = pa->hair_index - 1;
3398 medge->v2 = pa->hair_index;
3399
3400 dvert = hair_set_pinning(dvert, 1.0f);
3401
3402 mvert++;
3403 medge++;
3404 }
3405
3406 /* store root transform in cloth data */
3407 hair = &psys->clmd->hairdata[pa->hair_index + k];
3408 copy_v3_v3(hair->loc, root_mat[3]);
3409 copy_m3_m4(hair->rot, root_mat);
3410
3411 hair->radius = hair_radius;
3412 hair->bending_stiffness = bending_stiffness;
3413
3414 copy_v3_v3(mvert->co, co);
3415 mul_m4_v3(hairmat, mvert->co);
3416
3417 if (k) {
3418 medge->v1 = pa->hair_index + k - 1;
3419 medge->v2 = pa->hair_index + k;
3420 }
3421
3422 /* roots and disabled hairs should be 1.0, the rest can be anything from 0.0 to 1.0 */
3423 if (use_hair) {
3424 dvert = hair_set_pinning(dvert, key->weight);
3425 }
3426 else {
3427 dvert = hair_set_pinning(dvert, 1.0f);
3428 }
3429
3430 mvert++;
3431 if (k) {
3432 medge++;
3433 }
3434 }
3435
3436 hair_index += pa->totkey + 1;
3437 }
3438 }
3439 }
3440
do_hair_dynamics(ParticleSimulationData * sim)3441 static void do_hair_dynamics(ParticleSimulationData *sim)
3442 {
3443 ParticleSystem *psys = sim->psys;
3444 PARTICLE_P;
3445 EffectorWeights *clmd_effweights;
3446 int totpoint;
3447 int totedge;
3448 float(*deformedVerts)[3];
3449 bool realloc_roots;
3450
3451 if (!psys->clmd) {
3452 psys->clmd = (ClothModifierData *)BKE_modifier_new(eModifierType_Cloth);
3453 psys->clmd->sim_parms->goalspring = 0.0f;
3454 psys->clmd->sim_parms->flags |= CLOTH_SIMSETTINGS_FLAG_RESIST_SPRING_COMPRESS;
3455 psys->clmd->coll_parms->flags &= ~CLOTH_COLLSETTINGS_FLAG_SELF;
3456 }
3457
3458 /* count simulated points */
3459 totpoint = 0;
3460 totedge = 0;
3461 LOOP_PARTICLES
3462 {
3463 if (!(pa->flag & PARS_UNEXIST)) {
3464 /* "out" dm contains all hairs */
3465 totedge += pa->totkey;
3466 totpoint += pa->totkey + 1; /* +1 for virtual root point */
3467 }
3468 }
3469
3470 /* whether hair root info array has to be reallocated */
3471 realloc_roots = false;
3472 if (psys->hair_in_mesh) {
3473 Mesh *mesh = psys->hair_in_mesh;
3474 if (totpoint != mesh->totvert || totedge != mesh->totedge) {
3475 BKE_id_free(NULL, mesh);
3476 psys->hair_in_mesh = NULL;
3477 realloc_roots = true;
3478 }
3479 }
3480
3481 if (!psys->hair_in_mesh || !psys->clmd->hairdata || realloc_roots) {
3482 if (psys->clmd->hairdata) {
3483 MEM_freeN(psys->clmd->hairdata);
3484 psys->clmd->hairdata = NULL;
3485 }
3486 }
3487
3488 hair_create_input_mesh(sim, totpoint, totedge, &psys->hair_in_mesh, &psys->clmd->hairdata);
3489
3490 if (psys->hair_out_mesh) {
3491 BKE_id_free(NULL, psys->hair_out_mesh);
3492 }
3493
3494 psys->clmd->point_cache = psys->pointcache;
3495 /* for hair sim we replace the internal cloth effector weights temporarily
3496 * to use the particle settings
3497 */
3498 clmd_effweights = psys->clmd->sim_parms->effector_weights;
3499 psys->clmd->sim_parms->effector_weights = psys->part->effector_weights;
3500
3501 BKE_id_copy_ex(NULL, &psys->hair_in_mesh->id, (ID **)&psys->hair_out_mesh, LIB_ID_COPY_LOCALIZE);
3502 deformedVerts = BKE_mesh_vert_coords_alloc(psys->hair_out_mesh, NULL);
3503 clothModifier_do(
3504 psys->clmd, sim->depsgraph, sim->scene, sim->ob, psys->hair_in_mesh, deformedVerts);
3505 BKE_mesh_vert_coords_apply(psys->hair_out_mesh, deformedVerts);
3506
3507 MEM_freeN(deformedVerts);
3508
3509 /* restore cloth effector weights */
3510 psys->clmd->sim_parms->effector_weights = clmd_effweights;
3511 }
hair_step(ParticleSimulationData * sim,float cfra,const bool use_render_params)3512 static void hair_step(ParticleSimulationData *sim, float cfra, const bool use_render_params)
3513 {
3514 ParticleSystem *psys = sim->psys;
3515 ParticleSettings *part = psys->part;
3516 PARTICLE_P;
3517 float disp = psys_get_current_display_percentage(psys, use_render_params);
3518
3519 LOOP_PARTICLES
3520 {
3521 pa->size = part->size;
3522 if (part->randsize > 0.0f) {
3523 pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);
3524 }
3525
3526 if (psys_frand(psys, p) > disp) {
3527 pa->flag |= PARS_NO_DISP;
3528 }
3529 else {
3530 pa->flag &= ~PARS_NO_DISP;
3531 }
3532 }
3533
3534 if (psys->recalc & ID_RECALC_PSYS_RESET) {
3535 /* need this for changing subsurf levels */
3536 psys_calc_dmcache(sim->ob, sim->psmd->mesh_final, sim->psmd->mesh_original, psys);
3537
3538 if (psys->clmd) {
3539 cloth_free_modifier(psys->clmd);
3540 }
3541 }
3542
3543 /* dynamics with cloth simulation, psys->particles can be NULL with 0 particles T25519. */
3544 if (psys->part->type == PART_HAIR && psys->flag & PSYS_HAIR_DYNAMICS && psys->particles) {
3545 do_hair_dynamics(sim);
3546 }
3547
3548 /* following lines were removed r29079 but cause bug T22811, see report for details */
3549 psys_update_effectors(sim);
3550 psys_update_path_cache(sim, cfra, use_render_params);
3551
3552 psys->flag |= PSYS_HAIR_UPDATED;
3553 }
3554
save_hair(ParticleSimulationData * sim,float UNUSED (cfra))3555 static void save_hair(ParticleSimulationData *sim, float UNUSED(cfra))
3556 {
3557 Object *ob = sim->ob;
3558 ParticleSystem *psys = sim->psys;
3559 HairKey *key, *root;
3560 PARTICLE_P;
3561
3562 invert_m4_m4(ob->imat, ob->obmat);
3563
3564 psys->lattice_deform_data = psys_create_lattice_deform_data(sim);
3565
3566 if (psys->totpart == 0) {
3567 return;
3568 }
3569
3570 /* save new keys for elements if needed */
3571 LOOP_PARTICLES
3572 {
3573 /* first time alloc */
3574 if (pa->totkey == 0 || pa->hair == NULL) {
3575 pa->hair = MEM_callocN((psys->part->hair_step + 1) * sizeof(HairKey), "HairKeys");
3576 pa->totkey = 0;
3577 }
3578
3579 key = root = pa->hair;
3580 key += pa->totkey;
3581
3582 /* convert from global to geometry space */
3583 copy_v3_v3(key->co, pa->state.co);
3584 mul_m4_v3(ob->imat, key->co);
3585
3586 if (pa->totkey) {
3587 sub_v3_v3(key->co, root->co);
3588 psys_vec_rot_to_face(sim->psmd->mesh_final, pa, key->co);
3589 }
3590
3591 key->time = pa->state.time;
3592
3593 key->weight = 1.0f - key->time / 100.0f;
3594
3595 pa->totkey++;
3596
3597 /* Root is always in the origin of hair space
3598 * so we set it to be so after the last key is saved. */
3599 if (pa->totkey == psys->part->hair_step + 1) {
3600 zero_v3(root->co);
3601 }
3602 }
3603 }
3604
3605 /* Code for an adaptive time step based on the Courant-Friedrichs-Lewy
3606 * condition. */
3607 static const float MIN_TIMESTEP = 1.0f / 101.0f;
3608 /* Tolerance of 1.5 means the last subframe neither favors growing nor
3609 * shrinking (e.g if it were 1.3, the last subframe would tend to be too
3610 * small). */
3611 static const float TIMESTEP_EXPANSION_FACTOR = 0.1f;
3612 static const float TIMESTEP_EXPANSION_TOLERANCE = 1.5f;
3613
3614 /* Calculate the speed of the particle relative to the local scale of the
3615 * simulation. This should be called once per particle during a simulation
3616 * step, after the velocity has been updated. element_size defines the scale of
3617 * the simulation, and is typically the distance to neighboring particles. */
update_courant_num(ParticleSimulationData * sim,ParticleData * pa,float dtime,SPHData * sphdata,SpinLock * spin)3618 static void update_courant_num(
3619 ParticleSimulationData *sim, ParticleData *pa, float dtime, SPHData *sphdata, SpinLock *spin)
3620 {
3621 float relative_vel[3];
3622
3623 sub_v3_v3v3(relative_vel, pa->prev_state.vel, sphdata->flow);
3624
3625 const float courant_num = len_v3(relative_vel) * dtime / sphdata->element_size;
3626 if (sim->courant_num < courant_num) {
3627 BLI_spin_lock(spin);
3628 if (sim->courant_num < courant_num) {
3629 sim->courant_num = courant_num;
3630 }
3631 BLI_spin_unlock(spin);
3632 }
3633 }
get_base_time_step(ParticleSettings * part)3634 static float get_base_time_step(ParticleSettings *part)
3635 {
3636 return 1.0f / (float)(part->subframes + 1);
3637 }
3638 /* Update time step size to suit current conditions. */
update_timestep(ParticleSystem * psys,ParticleSimulationData * sim)3639 static void update_timestep(ParticleSystem *psys, ParticleSimulationData *sim)
3640 {
3641 float dt_target;
3642 if (sim->courant_num == 0.0f) {
3643 dt_target = 1.0f;
3644 }
3645 else {
3646 dt_target = psys->dt_frac * (psys->part->courant_target / sim->courant_num);
3647 }
3648
3649 /* Make sure the time step is reasonable. For some reason, the CLAMP macro
3650 * doesn't work here. The time step becomes too large. - z0r */
3651 if (dt_target < MIN_TIMESTEP) {
3652 dt_target = MIN_TIMESTEP;
3653 }
3654 else if (dt_target > get_base_time_step(psys->part)) {
3655 dt_target = get_base_time_step(psys->part);
3656 }
3657
3658 /* Decrease time step instantly, but increase slowly. */
3659 if (dt_target > psys->dt_frac) {
3660 psys->dt_frac = interpf(dt_target, psys->dt_frac, TIMESTEP_EXPANSION_FACTOR);
3661 }
3662 else {
3663 psys->dt_frac = dt_target;
3664 }
3665 }
3666
sync_timestep(ParticleSystem * psys,float t_frac)3667 static float sync_timestep(ParticleSystem *psys, float t_frac)
3668 {
3669 /* Sync with frame end if it's close. */
3670 if (t_frac == 1.0f) {
3671 return psys->dt_frac;
3672 }
3673 if (t_frac + (psys->dt_frac * TIMESTEP_EXPANSION_TOLERANCE) >= 1.0f) {
3674 return 1.0f - t_frac;
3675 }
3676
3677 return psys->dt_frac;
3678 }
3679
3680 /************************************************/
3681 /* System Core */
3682 /************************************************/
3683
3684 typedef struct DynamicStepSolverTaskData {
3685 ParticleSimulationData *sim;
3686
3687 float cfra;
3688 float timestep;
3689 float dtime;
3690
3691 SpinLock spin;
3692 } DynamicStepSolverTaskData;
3693
dynamics_step_sphdata_reduce(const void * __restrict UNUSED (userdata),void * __restrict UNUSED (join_v),void * __restrict chunk_v)3694 static void dynamics_step_sphdata_reduce(const void *__restrict UNUSED(userdata),
3695 void *__restrict UNUSED(join_v),
3696 void *__restrict chunk_v)
3697 {
3698 SPHData *sphdata = chunk_v;
3699
3700 psys_sph_flush_springs(sphdata);
3701 }
3702
dynamics_step_sph_ddr_task_cb_ex(void * __restrict userdata,const int p,const TaskParallelTLS * __restrict tls)3703 static void dynamics_step_sph_ddr_task_cb_ex(void *__restrict userdata,
3704 const int p,
3705 const TaskParallelTLS *__restrict tls)
3706 {
3707 DynamicStepSolverTaskData *data = userdata;
3708 ParticleSimulationData *sim = data->sim;
3709 ParticleSystem *psys = sim->psys;
3710 ParticleSettings *part = psys->part;
3711
3712 SPHData *sphdata = tls->userdata_chunk;
3713
3714 ParticleData *pa;
3715
3716 if ((pa = psys->particles + p)->state.time <= 0.0f) {
3717 return;
3718 }
3719
3720 /* do global forces & effectors */
3721 basic_integrate(sim, p, pa->state.time, data->cfra);
3722
3723 /* actual fluids calculations */
3724 sph_integrate(sim, pa, pa->state.time, sphdata);
3725
3726 if (sim->colliders) {
3727 collision_check(sim, p, pa->state.time, data->cfra);
3728 }
3729
3730 /* SPH particles are not physical particles, just interpolation
3731 * particles, thus rotation has not a direct sense for them */
3732 basic_rotate(part, pa, pa->state.time, data->timestep);
3733
3734 if (part->time_flag & PART_TIME_AUTOSF) {
3735 update_courant_num(sim, pa, data->dtime, sphdata, &data->spin);
3736 }
3737 }
3738
dynamics_step_sph_classical_basic_integrate_task_cb_ex(void * __restrict userdata,const int p,const TaskParallelTLS * __restrict UNUSED (tls))3739 static void dynamics_step_sph_classical_basic_integrate_task_cb_ex(
3740 void *__restrict userdata, const int p, const TaskParallelTLS *__restrict UNUSED(tls))
3741 {
3742 DynamicStepSolverTaskData *data = userdata;
3743 ParticleSimulationData *sim = data->sim;
3744 ParticleSystem *psys = sim->psys;
3745
3746 ParticleData *pa;
3747
3748 if ((pa = psys->particles + p)->state.time <= 0.0f) {
3749 return;
3750 }
3751
3752 basic_integrate(sim, p, pa->state.time, data->cfra);
3753 }
3754
dynamics_step_sph_classical_calc_density_task_cb_ex(void * __restrict userdata,const int p,const TaskParallelTLS * __restrict tls)3755 static void dynamics_step_sph_classical_calc_density_task_cb_ex(
3756 void *__restrict userdata, const int p, const TaskParallelTLS *__restrict tls)
3757 {
3758 DynamicStepSolverTaskData *data = userdata;
3759 ParticleSimulationData *sim = data->sim;
3760 ParticleSystem *psys = sim->psys;
3761
3762 SPHData *sphdata = tls->userdata_chunk;
3763
3764 ParticleData *pa;
3765
3766 if ((pa = psys->particles + p)->state.time <= 0.0f) {
3767 return;
3768 }
3769
3770 sphclassical_calc_dens(pa, pa->state.time, sphdata);
3771 }
3772
dynamics_step_sph_classical_integrate_task_cb_ex(void * __restrict userdata,const int p,const TaskParallelTLS * __restrict tls)3773 static void dynamics_step_sph_classical_integrate_task_cb_ex(void *__restrict userdata,
3774 const int p,
3775 const TaskParallelTLS *__restrict tls)
3776 {
3777 DynamicStepSolverTaskData *data = userdata;
3778 ParticleSimulationData *sim = data->sim;
3779 ParticleSystem *psys = sim->psys;
3780 ParticleSettings *part = psys->part;
3781
3782 SPHData *sphdata = tls->userdata_chunk;
3783
3784 ParticleData *pa;
3785
3786 if ((pa = psys->particles + p)->state.time <= 0.0f) {
3787 return;
3788 }
3789
3790 /* actual fluids calculations */
3791 sph_integrate(sim, pa, pa->state.time, sphdata);
3792
3793 if (sim->colliders) {
3794 collision_check(sim, p, pa->state.time, data->cfra);
3795 }
3796
3797 /* SPH particles are not physical particles, just interpolation
3798 * particles, thus rotation has not a direct sense for them */
3799 basic_rotate(part, pa, pa->state.time, data->timestep);
3800
3801 if (part->time_flag & PART_TIME_AUTOSF) {
3802 update_courant_num(sim, pa, data->dtime, sphdata, &data->spin);
3803 }
3804 }
3805
3806 /* unbaked particles are calculated dynamically */
dynamics_step(ParticleSimulationData * sim,float cfra)3807 static void dynamics_step(ParticleSimulationData *sim, float cfra)
3808 {
3809 ParticleSystem *psys = sim->psys;
3810 ParticleSettings *part = psys->part;
3811 BoidBrainData bbd;
3812 ParticleTexture ptex;
3813 PARTICLE_P;
3814 float timestep;
3815 /* frame & time changes */
3816 float dfra, dtime;
3817 float birthtime, dietime;
3818
3819 /* where have we gone in time since last time */
3820 dfra = cfra - psys->cfra;
3821
3822 timestep = psys_get_timestep(sim);
3823 dtime = dfra * timestep;
3824
3825 if (dfra < 0.0f) {
3826 LOOP_EXISTING_PARTICLES
3827 {
3828 psys_get_texture(sim, pa, &ptex, PAMAP_SIZE, cfra);
3829 pa->size = part->size * ptex.size;
3830 if (part->randsize > 0.0f) {
3831 pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);
3832 }
3833
3834 reset_particle(sim, pa, dtime, cfra);
3835 }
3836 return;
3837 }
3838
3839 /* for now do both, boids us 'rng' */
3840 sim->rng = BLI_rng_new_srandom(31415926 + (int)cfra + psys->seed);
3841
3842 psys_update_effectors(sim);
3843
3844 if (part->type != PART_HAIR) {
3845 sim->colliders = BKE_collider_cache_create(sim->depsgraph, sim->ob, part->collision_group);
3846 }
3847
3848 /* initialize physics type specific stuff */
3849 switch (part->phystype) {
3850 case PART_PHYS_BOIDS: {
3851 ParticleTarget *pt = psys->targets.first;
3852 bbd.sim = sim;
3853 bbd.part = part;
3854 bbd.cfra = cfra;
3855 bbd.dfra = dfra;
3856 bbd.timestep = timestep;
3857 bbd.rng = sim->rng;
3858
3859 psys_update_particle_tree(psys, cfra);
3860
3861 boids_precalc_rules(part, cfra);
3862
3863 for (; pt; pt = pt->next) {
3864 ParticleSystem *psys_target = psys_get_target_system(sim->ob, pt);
3865 if (psys_target && psys_target != psys) {
3866 psys_update_particle_tree(psys_target, cfra);
3867 }
3868 }
3869 break;
3870 }
3871 case PART_PHYS_FLUID: {
3872 ParticleTarget *pt = psys->targets.first;
3873 psys_update_particle_bvhtree(psys, cfra);
3874
3875 for (; pt;
3876 pt = pt->next) { /* Updating others systems particle tree for fluid-fluid interaction */
3877 if (pt->ob) {
3878 psys_update_particle_bvhtree(BLI_findlink(&pt->ob->particlesystem, pt->psys - 1), cfra);
3879 }
3880 }
3881 break;
3882 }
3883 }
3884 /* initialize all particles for dynamics */
3885 LOOP_SHOWN_PARTICLES
3886 {
3887 copy_particle_key(&pa->prev_state, &pa->state, 1);
3888
3889 psys_get_texture(sim, pa, &ptex, PAMAP_SIZE, cfra);
3890
3891 pa->size = part->size * ptex.size;
3892 if (part->randsize > 0.0f) {
3893 pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);
3894 }
3895
3896 birthtime = pa->time;
3897 dietime = pa->dietime;
3898
3899 /* store this, so we can do multiple loops over particles */
3900 pa->state.time = dfra;
3901
3902 if (dietime <= cfra && psys->cfra < dietime) {
3903 /* particle dies some time between this and last step */
3904 pa->state.time = dietime - ((birthtime > psys->cfra) ? birthtime : psys->cfra);
3905 pa->alive = PARS_DYING;
3906 }
3907 else if (birthtime <= cfra && birthtime >= psys->cfra) {
3908 /* particle is born some time between this and last step*/
3909 reset_particle(sim, pa, dfra * timestep, cfra);
3910 pa->alive = PARS_ALIVE;
3911 pa->state.time = cfra - birthtime;
3912 }
3913 else if (dietime < cfra) {
3914 /* nothing to be done when particle is dead */
3915 }
3916
3917 /* only reset unborn particles if they're shown or if the particle is born soon*/
3918 if (pa->alive == PARS_UNBORN &&
3919 (part->flag & PART_UNBORN || (cfra + psys->pointcache->step > pa->time))) {
3920 reset_particle(sim, pa, dtime, cfra);
3921 }
3922 else if (part->phystype == PART_PHYS_NO) {
3923 reset_particle(sim, pa, dtime, cfra);
3924 }
3925
3926 if (ELEM(pa->alive, PARS_ALIVE, PARS_DYING) == 0 ||
3927 (pa->flag & (PARS_UNEXIST | PARS_NO_DISP))) {
3928 pa->state.time = -1.f;
3929 }
3930 }
3931
3932 switch (part->phystype) {
3933 case PART_PHYS_NEWTON: {
3934 LOOP_DYNAMIC_PARTICLES
3935 {
3936 /* do global forces & effectors */
3937 basic_integrate(sim, p, pa->state.time, cfra);
3938
3939 /* deflection */
3940 if (sim->colliders) {
3941 collision_check(sim, p, pa->state.time, cfra);
3942 }
3943
3944 /* rotations */
3945 basic_rotate(part, pa, pa->state.time, timestep);
3946 }
3947 break;
3948 }
3949 case PART_PHYS_BOIDS: {
3950 LOOP_DYNAMIC_PARTICLES
3951 {
3952 bbd.goal_ob = NULL;
3953
3954 boid_brain(&bbd, p, pa);
3955
3956 if (pa->alive != PARS_DYING) {
3957 boid_body(&bbd, pa);
3958
3959 /* deflection */
3960 if (sim->colliders) {
3961 collision_check(sim, p, pa->state.time, cfra);
3962 }
3963 }
3964 }
3965 break;
3966 }
3967 case PART_PHYS_FLUID: {
3968 SPHData sphdata;
3969 psys_sph_init(sim, &sphdata);
3970
3971 DynamicStepSolverTaskData task_data = {
3972 .sim = sim,
3973 .cfra = cfra,
3974 .timestep = timestep,
3975 .dtime = dtime,
3976 };
3977
3978 BLI_spin_init(&task_data.spin);
3979
3980 if (part->fluid->solver == SPH_SOLVER_DDR) {
3981 /* Apply SPH forces using double-density relaxation algorithm
3982 * (Clavat et. al.) */
3983
3984 TaskParallelSettings settings;
3985 BLI_parallel_range_settings_defaults(&settings);
3986 settings.use_threading = (psys->totpart > 100);
3987 settings.userdata_chunk = &sphdata;
3988 settings.userdata_chunk_size = sizeof(sphdata);
3989 settings.func_reduce = dynamics_step_sphdata_reduce;
3990 BLI_task_parallel_range(
3991 0, psys->totpart, &task_data, dynamics_step_sph_ddr_task_cb_ex, &settings);
3992
3993 sph_springs_modify(psys, timestep);
3994 }
3995 else {
3996 /* SPH_SOLVER_CLASSICAL */
3997 /* Apply SPH forces using classical algorithm (due to Gingold
3998 * and Monaghan). Note that, unlike double-density relaxation,
3999 * this algorithm is separated into distinct loops. */
4000
4001 {
4002 TaskParallelSettings settings;
4003 BLI_parallel_range_settings_defaults(&settings);
4004 settings.use_threading = (psys->totpart > 100);
4005 BLI_task_parallel_range(0,
4006 psys->totpart,
4007 &task_data,
4008 dynamics_step_sph_classical_basic_integrate_task_cb_ex,
4009 &settings);
4010 }
4011
4012 /* calculate summation density */
4013 /* Note that we could avoid copying sphdata for each thread here (it's only read here),
4014 * but doubt this would gain us anything except confusion... */
4015 {
4016 TaskParallelSettings settings;
4017 BLI_parallel_range_settings_defaults(&settings);
4018 settings.use_threading = (psys->totpart > 100);
4019 settings.userdata_chunk = &sphdata;
4020 settings.userdata_chunk_size = sizeof(sphdata);
4021 settings.func_reduce = dynamics_step_sphdata_reduce;
4022 BLI_task_parallel_range(0,
4023 psys->totpart,
4024 &task_data,
4025 dynamics_step_sph_classical_calc_density_task_cb_ex,
4026 &settings);
4027 }
4028
4029 /* do global forces & effectors */
4030 {
4031 TaskParallelSettings settings;
4032 BLI_parallel_range_settings_defaults(&settings);
4033 settings.use_threading = (psys->totpart > 100);
4034 settings.userdata_chunk = &sphdata;
4035 settings.userdata_chunk_size = sizeof(sphdata);
4036 settings.func_reduce = dynamics_step_sphdata_reduce;
4037 BLI_task_parallel_range(0,
4038 psys->totpart,
4039 &task_data,
4040 dynamics_step_sph_classical_integrate_task_cb_ex,
4041 &settings);
4042 }
4043 }
4044
4045 BLI_spin_end(&task_data.spin);
4046
4047 psys_sph_finalize(&sphdata);
4048 break;
4049 }
4050 }
4051
4052 /* finalize particle state and time after dynamics */
4053 LOOP_DYNAMIC_PARTICLES
4054 {
4055 if (pa->alive == PARS_DYING) {
4056 pa->alive = PARS_DEAD;
4057 pa->state.time = pa->dietime;
4058 }
4059 else {
4060 pa->state.time = cfra;
4061 }
4062 }
4063
4064 BKE_collider_cache_free(&sim->colliders);
4065 BLI_rng_free(sim->rng);
4066 sim->rng = NULL;
4067 }
4068
update_children(ParticleSimulationData * sim,const bool use_render_params)4069 static void update_children(ParticleSimulationData *sim, const bool use_render_params)
4070 {
4071 if ((sim->psys->part->type == PART_HAIR) && (sim->psys->flag & PSYS_HAIR_DONE) == 0) {
4072 /* don't generate children while growing hair - waste of time */
4073 psys_free_children(sim->psys);
4074 }
4075 else if (sim->psys->part->childtype) {
4076 if (sim->psys->totchild != psys_get_tot_child(sim->scene, sim->psys, use_render_params)) {
4077 distribute_particles(sim, PART_FROM_CHILD);
4078 }
4079 else {
4080 /* Children are up to date, nothing to do. */
4081 }
4082 }
4083 else {
4084 psys_free_children(sim->psys);
4085 }
4086 }
4087 /* updates cached particles' alive & other flags etc..*/
cached_step(ParticleSimulationData * sim,float cfra,const bool use_render_params)4088 static void cached_step(ParticleSimulationData *sim, float cfra, const bool use_render_params)
4089 {
4090 ParticleSystem *psys = sim->psys;
4091 ParticleSettings *part = psys->part;
4092 ParticleTexture ptex;
4093 PARTICLE_P;
4094 float disp, dietime;
4095
4096 psys_update_effectors(sim);
4097
4098 disp = psys_get_current_display_percentage(psys, use_render_params);
4099
4100 LOOP_PARTICLES
4101 {
4102 psys_get_texture(sim, pa, &ptex, PAMAP_SIZE, cfra);
4103 pa->size = part->size * ptex.size;
4104 if (part->randsize > 0.0f) {
4105 pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);
4106 }
4107
4108 psys->lattice_deform_data = psys_create_lattice_deform_data(sim);
4109
4110 dietime = pa->dietime;
4111
4112 /* update alive status and push events */
4113 if (pa->time > cfra) {
4114 pa->alive = PARS_UNBORN;
4115 if (part->flag & PART_UNBORN && (psys->pointcache->flag & PTCACHE_EXTERNAL) == 0) {
4116 reset_particle(sim, pa, 0.0f, cfra);
4117 }
4118 }
4119 else if (dietime <= cfra) {
4120 pa->alive = PARS_DEAD;
4121 }
4122 else {
4123 pa->alive = PARS_ALIVE;
4124 }
4125
4126 if (psys->lattice_deform_data) {
4127 BKE_lattice_deform_data_destroy(psys->lattice_deform_data);
4128 psys->lattice_deform_data = NULL;
4129 }
4130
4131 if (psys_frand(psys, p) > disp) {
4132 pa->flag |= PARS_NO_DISP;
4133 }
4134 else {
4135 pa->flag &= ~PARS_NO_DISP;
4136 }
4137 }
4138 }
4139
particles_has_flip(short parttype)4140 static bool particles_has_flip(short parttype)
4141 {
4142 return (parttype == PART_FLUID_FLIP);
4143 }
4144
particles_has_tracer(short parttype)4145 static bool particles_has_tracer(short parttype)
4146 {
4147 return (parttype == PART_FLUID_TRACER);
4148 }
4149
particles_has_spray(short parttype)4150 static bool particles_has_spray(short parttype)
4151 {
4152 return ((parttype == PART_FLUID_SPRAY) || (parttype == PART_FLUID_SPRAYFOAM) ||
4153 (parttype == PART_FLUID_SPRAYFOAMBUBBLE));
4154 }
4155
particles_has_bubble(short parttype)4156 static bool particles_has_bubble(short parttype)
4157 {
4158 return ((parttype == PART_FLUID_BUBBLE) || (parttype == PART_FLUID_FOAMBUBBLE) ||
4159 (parttype == PART_FLUID_SPRAYFOAMBUBBLE));
4160 }
4161
particles_has_foam(short parttype)4162 static bool particles_has_foam(short parttype)
4163 {
4164 return ((parttype == PART_FLUID_FOAM) || (parttype == PART_FLUID_SPRAYFOAM) ||
4165 (parttype == PART_FLUID_SPRAYFOAMBUBBLE));
4166 }
4167
particles_fluid_step(ParticleSimulationData * sim,int cfra,const bool use_render_params)4168 static void particles_fluid_step(ParticleSimulationData *sim,
4169 int cfra,
4170 const bool use_render_params)
4171 {
4172 ParticleSystem *psys = sim->psys;
4173 if (psys->particles) {
4174 MEM_freeN(psys->particles);
4175 psys->particles = 0;
4176 psys->totpart = 0;
4177 }
4178
4179 #ifndef WITH_FLUID
4180 UNUSED_VARS(use_render_params, cfra);
4181 #else
4182 {
4183 Object *ob = sim->ob;
4184 FluidModifierData *fmd = (FluidModifierData *)BKE_modifiers_findby_type(ob,
4185 eModifierType_Fluid);
4186
4187 if (fmd && fmd->domain && fmd->domain->fluid) {
4188 FluidDomainSettings *fds = fmd->domain;
4189
4190 ParticleSettings *part = psys->part;
4191 ParticleData *pa = NULL;
4192
4193 int p, totpart = 0, tottypepart = 0;
4194 int flagActivePart, activeParts = 0;
4195 float posX, posY, posZ, velX, velY, velZ;
4196 float resX, resY, resZ;
4197 int upres = 1;
4198 char debugStrBuffer[256];
4199 float tmp[3] = {0}, tmp2[3] = {0};
4200
4201 /* Helper variables for scaling. */
4202 float min[3], max[3], size[3], cell_size_scaled[3], max_size;
4203
4204 /* Sanity check: parts also enabled in fluid domain? */
4205 if ((particles_has_flip(part->type) &&
4206 (fds->particle_type & FLUID_DOMAIN_PARTICLE_FLIP) == 0) ||
4207 (particles_has_spray(part->type) &&
4208 (fds->particle_type & FLUID_DOMAIN_PARTICLE_SPRAY) == 0) ||
4209 (particles_has_bubble(part->type) &&
4210 (fds->particle_type & FLUID_DOMAIN_PARTICLE_BUBBLE) == 0) ||
4211 (particles_has_foam(part->type) &&
4212 (fds->particle_type & FLUID_DOMAIN_PARTICLE_FOAM) == 0) ||
4213 (particles_has_tracer(part->type) &&
4214 (fds->particle_type & FLUID_DOMAIN_PARTICLE_TRACER) == 0)) {
4215 BLI_snprintf(debugStrBuffer,
4216 sizeof(debugStrBuffer),
4217 "particles_fluid_step::error - found particle system that is not enabled in "
4218 "fluid domain\n");
4219 return;
4220 }
4221
4222 /* Count particle amount. tottypepart is only important for snd particles. */
4223 if (part->type == PART_FLUID_FLIP) {
4224 tottypepart = totpart = manta_liquid_get_num_flip_particles(fds->fluid);
4225 }
4226 if (particles_has_spray(part->type) || particles_has_bubble(part->type) ||
4227 particles_has_foam(part->type) || particles_has_tracer(part->type)) {
4228 totpart = manta_liquid_get_num_snd_particles(fds->fluid);
4229
4230 /* tottypepart is the amount of particles of a snd particle type. */
4231 for (p = 0; p < totpart; p++) {
4232 flagActivePart = manta_liquid_get_snd_particle_flag_at(fds->fluid, p);
4233 if (particles_has_spray(part->type) && (flagActivePart & PARTICLE_TYPE_SPRAY)) {
4234 tottypepart++;
4235 }
4236 if (particles_has_bubble(part->type) && (flagActivePart & PARTICLE_TYPE_BUBBLE)) {
4237 tottypepart++;
4238 }
4239 if (particles_has_foam(part->type) && (flagActivePart & PARTICLE_TYPE_FOAM)) {
4240 tottypepart++;
4241 }
4242 if (particles_has_tracer(part->type) && (flagActivePart & PARTICLE_TYPE_TRACER)) {
4243 tottypepart++;
4244 }
4245 }
4246 }
4247 /* Sanity check: no particles present. */
4248 if (!totpart || !tottypepart) {
4249 return;
4250 }
4251
4252 /* How many particles to display? */
4253 tottypepart = (use_render_params) ? tottypepart : (part->disp * tottypepart) / 100;
4254
4255 part->totpart = tottypepart;
4256 part->sta = part->end = 1.0f;
4257 part->lifetime = sim->scene->r.efra + 1;
4258
4259 /* Allocate particles. */
4260 realloc_particles(sim, part->totpart);
4261
4262 /* Set some randomness when choosing which particles to display. */
4263 sim->rng = BLI_rng_new_srandom(31415926 + (int)cfra + psys->seed);
4264 double r, dispProb = (double)part->disp / 100.0;
4265
4266 /* Loop over *all* particles. Will break out of loop before tottypepart amount exceeded. */
4267 for (p = 0, pa = psys->particles; p < totpart; p++) {
4268
4269 /* Apply some randomness and determine which particles to skip. */
4270 r = BLI_rng_get_double(sim->rng);
4271 if (r > dispProb) {
4272 continue;
4273 }
4274
4275 /* flag, res, upres, pos, vel for FLIP and snd particles have different getters. */
4276 if (part->type == PART_FLUID_FLIP) {
4277 flagActivePart = manta_liquid_get_flip_particle_flag_at(fds->fluid, p);
4278
4279 resX = (float)manta_get_res_x(fds->fluid);
4280 resY = (float)manta_get_res_y(fds->fluid);
4281 resZ = (float)manta_get_res_z(fds->fluid);
4282
4283 upres = 1;
4284
4285 posX = manta_liquid_get_flip_particle_position_x_at(fds->fluid, p);
4286 posY = manta_liquid_get_flip_particle_position_y_at(fds->fluid, p);
4287 posZ = manta_liquid_get_flip_particle_position_z_at(fds->fluid, p);
4288
4289 velX = manta_liquid_get_flip_particle_velocity_x_at(fds->fluid, p);
4290 velY = manta_liquid_get_flip_particle_velocity_y_at(fds->fluid, p);
4291 velZ = manta_liquid_get_flip_particle_velocity_z_at(fds->fluid, p);
4292 }
4293 else if (particles_has_spray(part->type) || particles_has_bubble(part->type) ||
4294 particles_has_foam(part->type) || particles_has_tracer(part->type)) {
4295 flagActivePart = manta_liquid_get_snd_particle_flag_at(fds->fluid, p);
4296
4297 resX = (float)manta_liquid_get_particle_res_x(fds->fluid);
4298 resY = (float)manta_liquid_get_particle_res_y(fds->fluid);
4299 resZ = (float)manta_liquid_get_particle_res_z(fds->fluid);
4300
4301 upres = manta_liquid_get_particle_upres(fds->fluid);
4302
4303 posX = manta_liquid_get_snd_particle_position_x_at(fds->fluid, p);
4304 posY = manta_liquid_get_snd_particle_position_y_at(fds->fluid, p);
4305 posZ = manta_liquid_get_snd_particle_position_z_at(fds->fluid, p);
4306
4307 velX = manta_liquid_get_snd_particle_velocity_x_at(fds->fluid, p);
4308 velY = manta_liquid_get_snd_particle_velocity_y_at(fds->fluid, p);
4309 velZ = manta_liquid_get_snd_particle_velocity_z_at(fds->fluid, p);
4310 }
4311 else {
4312 BLI_snprintf(debugStrBuffer,
4313 sizeof(debugStrBuffer),
4314 "particles_fluid_step::error - unknown particle system type\n");
4315 return;
4316 }
4317 # if 0
4318 /* Debugging: Print type of particle system and current particles. */
4319 printf("system type is %d and particle type is %d\n", part->type, flagActivePart);
4320 # endif
4321
4322 /* Type of particle must match current particle system type
4323 * (only important for snd particles). */
4324 if ((flagActivePart & PARTICLE_TYPE_SPRAY) && !particles_has_spray(part->type)) {
4325 continue;
4326 }
4327 if ((flagActivePart & PARTICLE_TYPE_BUBBLE) && !particles_has_bubble(part->type)) {
4328 continue;
4329 }
4330 if ((flagActivePart & PARTICLE_TYPE_FOAM) && !particles_has_foam(part->type)) {
4331 continue;
4332 }
4333 if ((flagActivePart & PARTICLE_TYPE_TRACER) && !particles_has_tracer(part->type)) {
4334 continue;
4335 }
4336 # if 0
4337 /* Debugging: Print type of particle system and current particles. */
4338 printf("system type is %d and particle type is %d\n", part->type, flagActivePart);
4339 # endif
4340 /* Particle system has allocated 'tottypepart' particles - so break early before exceeded.
4341 */
4342 if (activeParts >= tottypepart) {
4343 break;
4344 }
4345
4346 /* Only show active particles, i.e. filter out dead particles that just Mantaflow needs.
4347 * Mantaflow convention: PARTICLE_TYPE_DELETE == inactive particle. */
4348 if ((flagActivePart & PARTICLE_TYPE_DELETE) == 0) {
4349 activeParts++;
4350
4351 /* Use particle system settings for particle size. */
4352 pa->size = part->size;
4353 if (part->randsize > 0.0f) {
4354 pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);
4355 }
4356
4357 /* Get size (dimension) but considering scaling */
4358 copy_v3_v3(cell_size_scaled, fds->cell_size);
4359 mul_v3_v3(cell_size_scaled, ob->scale);
4360 madd_v3fl_v3fl_v3fl_v3i(min, fds->p0, cell_size_scaled, fds->res_min);
4361 madd_v3fl_v3fl_v3fl_v3i(max, fds->p0, cell_size_scaled, fds->res_max);
4362 sub_v3_v3v3(size, max, min);
4363
4364 /* Biggest dimension will be used for up-scaling. */
4365 max_size = MAX3(size[0] / (float)upres, size[1] / (float)upres, size[2] / (float)upres);
4366
4367 /* Set particle position. */
4368 const float posParticle[3] = {posX, posY, posZ};
4369 copy_v3_v3(pa->state.co, posParticle);
4370
4371 /* Normalize to unit cube around 0. */
4372 float resDomain[3] = {resX, resY, resZ};
4373 mul_v3_fl(resDomain, 0.5f);
4374 sub_v3_v3(pa->state.co, resDomain);
4375 mul_v3_fl(pa->state.co, fds->dx);
4376
4377 /* Match domain dimension / size. */
4378 float scaleAbs[3] = {
4379 1. / fabsf(ob->scale[0]), 1. / fabsf(ob->scale[1]), 1. / fabsf(ob->scale[2])};
4380 mul_v3_fl(scaleAbs, max_size);
4381 mul_v3_v3(pa->state.co, scaleAbs);
4382
4383 /* Match domain scale. */
4384 mul_m4_v3(ob->obmat, pa->state.co);
4385
4386 /* Add origin offset to particle position. */
4387 zero_v3(tmp);
4388 zero_v3(tmp2);
4389 sub_v3_v3v3(tmp2, fds->p1, fds->p0);
4390 mul_v3_fl(tmp2, 0.5f);
4391 add_v3_v3v3(tmp, tmp, fds->p1);
4392 sub_v3_v3(tmp, tmp2);
4393 mul_v3_v3(tmp, ob->scale);
4394 add_v3_v3(pa->state.co, tmp);
4395 # if 0
4396 /* Debugging: Print particle coordinates. */
4397 printf("pa->state.co[0]: %f, pa->state.co[1]: %f, pa->state.co[2]: %f\n",
4398 pa->state.co[0], pa->state.co[1], pa->state.co[2]);
4399 # endif
4400 /* Set particle velocity. */
4401 const float velParticle[3] = {velX, velY, velZ};
4402 copy_v3_v3(pa->state.vel, velParticle);
4403 mul_v3_fl(pa->state.vel, fds->dx);
4404 # if 0
4405 /* Debugging: Print particle velocity. */
4406 printf("pa->state.vel[0]: %f, pa->state.vel[1]: %f, pa->state.vel[2]: %f\n",
4407 pa->state.vel[0], pa->state.vel[1], pa->state.vel[2]);
4408 # endif
4409 /* Set default angular velocity and particle rotation. */
4410 zero_v3(pa->state.ave);
4411 unit_qt(pa->state.rot);
4412
4413 pa->time = 1.f;
4414 pa->dietime = sim->scene->r.efra + 1;
4415 pa->lifetime = sim->scene->r.efra;
4416 pa->alive = PARS_ALIVE;
4417
4418 /* Increasing particle settings pointer only for active particles. */
4419 pa++;
4420 }
4421 }
4422 # if 0
4423 /* Debugging: Print number of active particles. */
4424 printf("active parts: %d\n", activeParts);
4425 # endif
4426 totpart = psys->totpart = part->totpart = activeParts;
4427
4428 BLI_rng_free(sim->rng);
4429 sim->rng = NULL;
4430
4431 } /* Fluid sim particles done. */
4432 }
4433 #endif /* WITH_FLUID */
4434 }
4435
emit_particles(ParticleSimulationData * sim,PTCacheID * pid,float UNUSED (cfra))4436 static int emit_particles(ParticleSimulationData *sim, PTCacheID *pid, float UNUSED(cfra))
4437 {
4438 ParticleSystem *psys = sim->psys;
4439 int oldtotpart = psys->totpart;
4440 int totpart = tot_particles(psys, pid);
4441
4442 if (totpart != oldtotpart) {
4443 realloc_particles(sim, totpart);
4444 }
4445
4446 return totpart - oldtotpart;
4447 }
4448
4449 /**
4450 * Calculates the next state for all particles of the system.
4451 * In particles code most 'cfra - ending' are frames,
4452 * 'time - ending' are 'cfra * timestep' (seconds).
4453 *
4454 * 1. Emit particles
4455 * 2. Check cache (if used) and return if frame is cached
4456 * 3. Do dynamics
4457 * 4. Save to cache
4458 */
system_step(ParticleSimulationData * sim,float cfra,const bool use_render_params)4459 static void system_step(ParticleSimulationData *sim, float cfra, const bool use_render_params)
4460 {
4461 ParticleSystem *psys = sim->psys;
4462 ParticleSettings *part = psys->part;
4463 PointCache *cache = psys->pointcache;
4464 PTCacheID ptcacheid, *pid = NULL;
4465 PARTICLE_P;
4466 float disp, cache_cfra = cfra; /*, *vg_vel= 0, *vg_tan= 0, *vg_rot= 0, *vg_size= 0; */
4467 int startframe = 0, endframe = 100, oldtotpart = 0;
4468
4469 /* cache shouldn't be used for hair or "continue physics" */
4470 if (part->type != PART_HAIR) {
4471 psys_clear_temp_pointcache(psys);
4472
4473 /* set suitable cache range automatically */
4474 if ((cache->flag & (PTCACHE_BAKING | PTCACHE_BAKED)) == 0) {
4475 psys_get_pointcache_start_end(sim->scene, psys, &cache->startframe, &cache->endframe);
4476 }
4477
4478 pid = &ptcacheid;
4479 BKE_ptcache_id_from_particles(pid, sim->ob, psys);
4480
4481 BKE_ptcache_id_time(pid, sim->scene, 0.0f, &startframe, &endframe, NULL);
4482
4483 /* clear everything on start frame, or when psys needs full reset! */
4484 if ((cfra == startframe) || (psys->recalc & ID_RECALC_PSYS_RESET)) {
4485 BKE_ptcache_id_reset(sim->scene, pid, PTCACHE_RESET_OUTDATED);
4486 BKE_ptcache_validate(cache, startframe);
4487 cache->flag &= ~PTCACHE_REDO_NEEDED;
4488 }
4489
4490 CLAMP(cache_cfra, startframe, endframe);
4491 }
4492
4493 /* 1. emit particles and redo particles if needed */
4494 oldtotpart = psys->totpart;
4495 if (emit_particles(sim, pid, cfra) || psys->recalc & ID_RECALC_PSYS_RESET) {
4496 distribute_particles(sim, part->from);
4497 initialize_all_particles(sim);
4498 /* reset only just created particles (on startframe all particles are recreated) */
4499 reset_all_particles(sim, 0.0, cfra, oldtotpart);
4500 free_unexisting_particles(sim);
4501
4502 if (psys->fluid_springs) {
4503 MEM_freeN(psys->fluid_springs);
4504 psys->fluid_springs = NULL;
4505 }
4506
4507 psys->tot_fluidsprings = psys->alloc_fluidsprings = 0;
4508
4509 /* flag for possible explode modifiers after this system */
4510 sim->psmd->flag |= eParticleSystemFlag_Pars;
4511
4512 BKE_ptcache_id_clear(pid, PTCACHE_CLEAR_AFTER, cfra);
4513 }
4514
4515 /* 2. try to read from the cache */
4516 if (pid) {
4517 int cache_result = BKE_ptcache_read(pid, cache_cfra, true);
4518
4519 if (ELEM(cache_result, PTCACHE_READ_EXACT, PTCACHE_READ_INTERPOLATED)) {
4520 cached_step(sim, cfra, use_render_params);
4521 update_children(sim, use_render_params);
4522 psys_update_path_cache(sim, cfra, use_render_params);
4523
4524 BKE_ptcache_validate(cache, (int)cache_cfra);
4525
4526 if (cache_result == PTCACHE_READ_INTERPOLATED && cache->flag & PTCACHE_REDO_NEEDED) {
4527 BKE_ptcache_write(pid, (int)cache_cfra);
4528 }
4529
4530 return;
4531 }
4532 /* Cache is supposed to be baked, but no data was found so bail out */
4533 if (cache->flag & PTCACHE_BAKED) {
4534 psys_reset(psys, PSYS_RESET_CACHE_MISS);
4535 return;
4536 }
4537 if (cache_result == PTCACHE_READ_OLD) {
4538 psys->cfra = (float)cache->simframe;
4539 cached_step(sim, psys->cfra, use_render_params);
4540 }
4541
4542 /* if on second frame, write cache for first frame */
4543 if (psys->cfra == startframe && (cache->flag & PTCACHE_OUTDATED || cache->last_exact == 0)) {
4544 BKE_ptcache_write(pid, startframe);
4545 }
4546 }
4547 else {
4548 BKE_ptcache_invalidate(cache);
4549 }
4550
4551 /* 3. do dynamics */
4552 /* set particles to be not calculated TODO: can't work with pointcache */
4553 disp = psys_get_current_display_percentage(psys, use_render_params);
4554
4555 LOOP_PARTICLES
4556 {
4557 if (psys_frand(psys, p) > disp) {
4558 pa->flag |= PARS_NO_DISP;
4559 }
4560 else {
4561 pa->flag &= ~PARS_NO_DISP;
4562 }
4563 }
4564
4565 if (psys->totpart) {
4566 int dframe, totframesback = 0;
4567 float t_frac, dt_frac;
4568
4569 /* handle negative frame start at the first frame by doing
4570 * all the steps before the first frame */
4571 if ((int)cfra == startframe && part->sta < startframe) {
4572 totframesback = (startframe - (int)part->sta);
4573 }
4574
4575 if (!(part->time_flag & PART_TIME_AUTOSF)) {
4576 /* Constant time step */
4577 psys->dt_frac = get_base_time_step(part);
4578 }
4579 else if ((int)cfra == startframe) {
4580 /* Variable time step; initialize to sub-frames. */
4581 psys->dt_frac = get_base_time_step(part);
4582 }
4583 else if (psys->dt_frac < MIN_TIMESTEP) {
4584 /* Variable time step; subsequent frames */
4585 psys->dt_frac = MIN_TIMESTEP;
4586 }
4587
4588 for (dframe = -totframesback; dframe <= 0; dframe++) {
4589 /* simulate each subframe */
4590 dt_frac = psys->dt_frac;
4591 for (t_frac = dt_frac; t_frac <= 1.0f; t_frac += dt_frac) {
4592 sim->courant_num = 0.0f;
4593 dynamics_step(sim, cfra + dframe + t_frac - 1.f);
4594 psys->cfra = cfra + dframe + t_frac - 1.f;
4595
4596 if (part->time_flag & PART_TIME_AUTOSF) {
4597 update_timestep(psys, sim);
4598 }
4599 /* Even without AUTOSF dt_frac may not add up to 1.0 due to float precision. */
4600 dt_frac = sync_timestep(psys, t_frac);
4601 }
4602 }
4603 }
4604
4605 /* 4. only write cache starting from second frame */
4606 if (pid) {
4607 BKE_ptcache_validate(cache, (int)cache_cfra);
4608 if ((int)cache_cfra != startframe) {
4609 BKE_ptcache_write(pid, (int)cache_cfra);
4610 }
4611 }
4612
4613 update_children(sim, use_render_params);
4614
4615 /* cleanup */
4616 if (psys->lattice_deform_data) {
4617 BKE_lattice_deform_data_destroy(psys->lattice_deform_data);
4618 psys->lattice_deform_data = NULL;
4619 }
4620 }
4621
4622 /* system type has changed so set sensible defaults and clear non applicable flags */
psys_changed_type(Object * ob,ParticleSystem * psys)4623 void psys_changed_type(Object *ob, ParticleSystem *psys)
4624 {
4625 ParticleSettings *part = psys->part;
4626 PTCacheID pid;
4627
4628 BKE_ptcache_id_from_particles(&pid, ob, psys);
4629
4630 if (part->phystype != PART_PHYS_KEYED) {
4631 psys->flag &= ~PSYS_KEYED;
4632 }
4633
4634 if (part->type == PART_HAIR) {
4635 if (ELEM(part->ren_as, PART_DRAW_NOT, PART_DRAW_PATH, PART_DRAW_OB, PART_DRAW_GR) == 0) {
4636 part->ren_as = PART_DRAW_PATH;
4637 }
4638
4639 if (part->distr == PART_DISTR_GRID) {
4640 part->distr = PART_DISTR_JIT;
4641 }
4642
4643 if (ELEM(part->draw_as, PART_DRAW_NOT, PART_DRAW_REND, PART_DRAW_PATH) == 0) {
4644 part->draw_as = PART_DRAW_REND;
4645 }
4646
4647 CLAMP(part->path_start, 0.0f, 100.0f);
4648 CLAMP(part->path_end, 0.0f, 100.0f);
4649
4650 BKE_ptcache_id_clear(&pid, PTCACHE_CLEAR_ALL, 0);
4651 }
4652 else {
4653 free_hair(ob, psys, 1);
4654
4655 CLAMP(part->path_start, 0.0f, MAX2(100.0f, part->end + part->lifetime));
4656 CLAMP(part->path_end, 0.0f, MAX2(100.0f, part->end + part->lifetime));
4657 }
4658
4659 psys_reset(psys, PSYS_RESET_ALL);
4660 }
psys_check_boid_data(ParticleSystem * psys)4661 void psys_check_boid_data(ParticleSystem *psys)
4662 {
4663 BoidParticle *bpa;
4664 PARTICLE_P;
4665
4666 pa = psys->particles;
4667
4668 if (!pa) {
4669 return;
4670 }
4671
4672 if (psys->part && psys->part->phystype == PART_PHYS_BOIDS) {
4673 if (!pa->boid) {
4674 bpa = MEM_callocN(psys->totpart * sizeof(BoidParticle), "Boid Data");
4675
4676 LOOP_PARTICLES
4677 {
4678 pa->boid = bpa++;
4679 }
4680 }
4681 }
4682 else if (pa->boid) {
4683 MEM_freeN(pa->boid);
4684 LOOP_PARTICLES
4685 {
4686 pa->boid = NULL;
4687 }
4688 }
4689 }
4690
BKE_particlesettings_fluid_default_settings(ParticleSettings * part)4691 void BKE_particlesettings_fluid_default_settings(ParticleSettings *part)
4692 {
4693 SPHFluidSettings *fluid = part->fluid;
4694
4695 fluid->spring_k = 0.f;
4696 fluid->plasticity_constant = 0.1f;
4697 fluid->yield_ratio = 0.1f;
4698 fluid->rest_length = 1.f;
4699 fluid->viscosity_omega = 2.f;
4700 fluid->viscosity_beta = 0.1f;
4701 fluid->stiffness_k = 1.f;
4702 fluid->stiffness_knear = 1.f;
4703 fluid->rest_density = 1.f;
4704 fluid->buoyancy = 0.f;
4705 fluid->radius = 1.f;
4706 fluid->flag |= SPH_FAC_REPULSION | SPH_FAC_DENSITY | SPH_FAC_RADIUS | SPH_FAC_VISCOSITY |
4707 SPH_FAC_REST_LENGTH;
4708 }
4709
psys_prepare_physics(ParticleSimulationData * sim)4710 static void psys_prepare_physics(ParticleSimulationData *sim)
4711 {
4712 ParticleSettings *part = sim->psys->part;
4713
4714 if (ELEM(part->phystype, PART_PHYS_NO, PART_PHYS_KEYED)) {
4715 PTCacheID pid;
4716 BKE_ptcache_id_from_particles(&pid, sim->ob, sim->psys);
4717 BKE_ptcache_id_clear(&pid, PTCACHE_CLEAR_ALL, 0);
4718 }
4719 else {
4720 free_keyed_keys(sim->psys);
4721 sim->psys->flag &= ~PSYS_KEYED;
4722 }
4723
4724 /* RNA Update must ensure this is true. */
4725 if (part->phystype == PART_PHYS_BOIDS) {
4726 BLI_assert(part->boids != NULL);
4727 }
4728 else if (part->phystype == PART_PHYS_FLUID) {
4729 BLI_assert(part->fluid != NULL);
4730 }
4731
4732 psys_check_boid_data(sim->psys);
4733 }
hair_needs_recalc(ParticleSystem * psys)4734 static int hair_needs_recalc(ParticleSystem *psys)
4735 {
4736 if (!(psys->flag & PSYS_EDITED) && (!psys->edit || !psys->edit->edited) &&
4737 ((psys->flag & PSYS_HAIR_DONE) == 0 || psys->recalc & ID_RECALC_PSYS_RESET ||
4738 (psys->part->flag & PART_HAIR_REGROW && !psys->edit))) {
4739 return 1;
4740 }
4741
4742 return 0;
4743 }
4744
particle_settings_localize(ParticleSettings * particle_settings)4745 static ParticleSettings *particle_settings_localize(ParticleSettings *particle_settings)
4746 {
4747 ParticleSettings *particle_settings_local = (ParticleSettings *)BKE_id_copy_ex(
4748 NULL, (ID *)&particle_settings->id, NULL, LIB_ID_COPY_LOCALIZE);
4749 return particle_settings_local;
4750 }
4751
particle_settings_free_local(ParticleSettings * particle_settings)4752 static void particle_settings_free_local(ParticleSettings *particle_settings)
4753 {
4754 BKE_libblock_free_datablock(&particle_settings->id, 0);
4755 BKE_libblock_free_data(&particle_settings->id, false);
4756 MEM_freeN(particle_settings);
4757 }
4758
4759 /* main particle update call, checks that things are ok on the large scale and
4760 * then advances in to actual particle calculations depending on particle type */
particle_system_update(struct Depsgraph * depsgraph,Scene * scene,Object * ob,ParticleSystem * psys,const bool use_render_params)4761 void particle_system_update(struct Depsgraph *depsgraph,
4762 Scene *scene,
4763 Object *ob,
4764 ParticleSystem *psys,
4765 const bool use_render_params)
4766 {
4767 ParticleSimulationData sim = {0};
4768 ParticleSettings *part = psys->part;
4769 ParticleSystem *psys_orig = psys_orig_get(psys);
4770 float cfra;
4771 ParticleSystemModifierData *psmd = psys_get_modifier(ob, psys);
4772
4773 /* drawdata is outdated after ANY change */
4774 if (psys->pdd) {
4775 psys->pdd->flag &= ~PARTICLE_DRAW_DATA_UPDATED;
4776 }
4777
4778 if (!psys_check_enabled(ob, psys, use_render_params)) {
4779 return;
4780 }
4781
4782 cfra = DEG_get_ctime(depsgraph);
4783
4784 sim.depsgraph = depsgraph;
4785 sim.scene = scene;
4786 sim.ob = ob;
4787 sim.psys = psys;
4788 sim.psmd = psmd;
4789
4790 /* system was already updated from modifier stack */
4791 if (sim.psmd->flag & eParticleSystemFlag_psys_updated) {
4792 sim.psmd->flag &= ~eParticleSystemFlag_psys_updated;
4793 /* make sure it really was updated to cfra */
4794 if (psys->cfra == cfra) {
4795 return;
4796 }
4797 }
4798
4799 if (!sim.psmd->mesh_final) {
4800 return;
4801 }
4802
4803 if (part->from != PART_FROM_VERT) {
4804 BKE_mesh_tessface_ensure(sim.psmd->mesh_final);
4805 }
4806
4807 /* to verify if we need to restore object afterwards */
4808 psys->flag &= ~PSYS_OB_ANIM_RESTORE;
4809
4810 if (psys->recalc & ID_RECALC_PSYS_RESET) {
4811 psys->totunexist = 0;
4812 }
4813
4814 /* setup necessary physics type dependent additional data if it doesn't yet exist */
4815 psys_prepare_physics(&sim);
4816
4817 if (part->type == PART_HAIR) {
4818 /* nothing to do so bail out early */
4819 if (psys->totpart == 0 && part->totpart == 0) {
4820 psys_free_path_cache(psys, NULL);
4821 free_hair(ob, psys, 0);
4822 psys->flag |= PSYS_HAIR_DONE;
4823 }
4824 /* (re-)create hair */
4825 else if (hair_needs_recalc(psys)) {
4826 float hcfra = 0.0f;
4827 int i, recalc = psys->recalc;
4828
4829 free_hair(ob, psys, 0);
4830
4831 if (psys_orig->edit && psys_orig->free_edit) {
4832 psys_orig->free_edit(psys_orig->edit);
4833 psys_orig->edit = NULL;
4834 psys_orig->free_edit = NULL;
4835 }
4836
4837 /* first step is negative so particles get killed and reset */
4838 psys->cfra = 1.0f;
4839
4840 ParticleSettings *part_local = part;
4841 if ((part->flag & PART_HAIR_REGROW) == 0) {
4842 part_local = particle_settings_localize(part);
4843 psys->part = part_local;
4844 }
4845
4846 for (i = 0; i <= part->hair_step; i++) {
4847 hcfra = 100.0f * (float)i / (float)psys->part->hair_step;
4848 if ((part->flag & PART_HAIR_REGROW) == 0) {
4849 const AnimationEvalContext anim_eval_context = BKE_animsys_eval_context_construct(
4850 depsgraph, hcfra);
4851 BKE_animsys_evaluate_animdata(
4852 &part_local->id, part_local->adt, &anim_eval_context, ADT_RECALC_ANIM, false);
4853 }
4854 system_step(&sim, hcfra, use_render_params);
4855 psys->cfra = hcfra;
4856 psys->recalc = 0;
4857 save_hair(&sim, hcfra);
4858 }
4859
4860 if (part_local != part) {
4861 particle_settings_free_local(part_local);
4862 psys->part = part;
4863 }
4864
4865 psys->flag |= PSYS_HAIR_DONE;
4866 psys->recalc = recalc;
4867 }
4868 else if (psys->flag & PSYS_EDITED) {
4869 psys->flag |= PSYS_HAIR_DONE;
4870 }
4871
4872 if (psys->flag & PSYS_HAIR_DONE) {
4873 hair_step(&sim, cfra, use_render_params);
4874 }
4875 }
4876 else if (particles_has_flip(part->type) || particles_has_spray(part->type) ||
4877 particles_has_bubble(part->type) || particles_has_foam(part->type) ||
4878 particles_has_tracer(part->type)) {
4879 particles_fluid_step(&sim, (int)cfra, use_render_params);
4880 }
4881 else {
4882 switch (part->phystype) {
4883 case PART_PHYS_NO:
4884 case PART_PHYS_KEYED: {
4885 PARTICLE_P;
4886 float disp = psys_get_current_display_percentage(psys, use_render_params);
4887 bool free_unexisting = false;
4888
4889 /* Particles without dynamics haven't been reset yet because they don't use pointcache */
4890 if (psys->recalc & ID_RECALC_PSYS_RESET) {
4891 psys_reset(psys, PSYS_RESET_ALL);
4892 }
4893
4894 if (emit_particles(&sim, NULL, cfra) || (psys->recalc & ID_RECALC_PSYS_RESET)) {
4895 free_keyed_keys(psys);
4896 distribute_particles(&sim, part->from);
4897 initialize_all_particles(&sim);
4898 free_unexisting = true;
4899
4900 /* flag for possible explode modifiers after this system */
4901 sim.psmd->flag |= eParticleSystemFlag_Pars;
4902 }
4903
4904 LOOP_EXISTING_PARTICLES
4905 {
4906 pa->size = part->size;
4907 if (part->randsize > 0.0f) {
4908 pa->size *= 1.0f - part->randsize * psys_frand(psys, p + 1);
4909 }
4910
4911 reset_particle(&sim, pa, 0.0, cfra);
4912
4913 if (psys_frand(psys, p) > disp) {
4914 pa->flag |= PARS_NO_DISP;
4915 }
4916 else {
4917 pa->flag &= ~PARS_NO_DISP;
4918 }
4919 }
4920
4921 /* free unexisting after resetting particles */
4922 if (free_unexisting) {
4923 free_unexisting_particles(&sim);
4924 }
4925
4926 if (part->phystype == PART_PHYS_KEYED) {
4927 psys_count_keyed_targets(&sim);
4928 set_keyed_keys(&sim);
4929 psys_update_path_cache(&sim, (int)cfra, use_render_params);
4930 }
4931 break;
4932 }
4933 default: {
4934 /* the main dynamic particle system step */
4935 system_step(&sim, cfra, use_render_params);
4936 break;
4937 }
4938 }
4939 }
4940
4941 /* make sure emitter is left at correct time (particle emission can change this) */
4942 if (psys->flag & PSYS_OB_ANIM_RESTORE) {
4943 evaluate_emitter_anim(depsgraph, scene, ob, cfra);
4944 psys->flag &= ~PSYS_OB_ANIM_RESTORE;
4945 }
4946
4947 if (psys_orig->edit) {
4948 psys_orig->edit->flags |= PT_CACHE_EDIT_UPDATE_PARTICLE_FROM_EVAL;
4949 }
4950
4951 psys->cfra = cfra;
4952 psys->recalc = 0;
4953
4954 if (DEG_is_active(depsgraph)) {
4955 if (psys_orig != psys) {
4956 if (psys_orig->edit != NULL && psys_orig->edit->psys == psys_orig) {
4957 psys_orig->edit->psys_eval = psys;
4958 psys_orig->edit->psmd_eval = psmd;
4959 }
4960 psys_orig->flag = (psys->flag & ~PSYS_SHARED_CACHES);
4961 psys_orig->cfra = psys->cfra;
4962 psys_orig->recalc = psys->recalc;
4963 psys_orig->part->totpart = part->totpart;
4964 }
4965 }
4966
4967 /* Save matrix for duplicators,
4968 * at rendertime the actual dupliobject's matrix is used so don't update! */
4969 invert_m4_m4(psys->imat, ob->obmat);
4970
4971 BKE_particle_batch_cache_dirty_tag(psys, BKE_PARTICLE_BATCH_DIRTY_ALL);
4972 }
4973
4974 /* ID looper */
4975
BKE_particlesystem_id_loop(ParticleSystem * psys,ParticleSystemIDFunc func,void * userdata)4976 void BKE_particlesystem_id_loop(ParticleSystem *psys, ParticleSystemIDFunc func, void *userdata)
4977 {
4978 ParticleTarget *pt;
4979
4980 func(psys, (ID **)&psys->part, userdata, IDWALK_CB_USER | IDWALK_CB_NEVER_NULL);
4981 func(psys, (ID **)&psys->target_ob, userdata, IDWALK_CB_NOP);
4982 func(psys, (ID **)&psys->parent, userdata, IDWALK_CB_NOP);
4983
4984 for (pt = psys->targets.first; pt; pt = pt->next) {
4985 func(psys, (ID **)&pt->ob, userdata, IDWALK_CB_NOP);
4986 }
4987
4988 /* Even though psys->part should never be NULL, this can happen as an exception during deletion.
4989 * See ID_REMAP_SKIP/FORCE/FLAG_NEVER_NULL_USAGE in BKE_library_remap. */
4990 if (psys->part && psys->part->phystype == PART_PHYS_BOIDS) {
4991 ParticleData *pa;
4992 int p;
4993
4994 for (p = 0, pa = psys->particles; p < psys->totpart; p++, pa++) {
4995 func(psys, (ID **)&pa->boid->ground, userdata, IDWALK_CB_NOP);
4996 }
4997 }
4998 }
4999
BKE_particlesystem_reset_all(struct Object * object)5000 void BKE_particlesystem_reset_all(struct Object *object)
5001 {
5002 for (ModifierData *md = object->modifiers.first; md != NULL; md = md->next) {
5003 if (md->type != eModifierType_ParticleSystem) {
5004 continue;
5005 }
5006 ParticleSystemModifierData *psmd = (ParticleSystemModifierData *)md;
5007 ParticleSystem *psys = psmd->psys;
5008 psys->recalc |= ID_RECALC_PSYS_RESET;
5009 }
5010 }
5011
5012 /* **** Depsgraph evaluation **** */
5013
BKE_particle_settings_eval_reset(struct Depsgraph * depsgraph,ParticleSettings * particle_settings)5014 void BKE_particle_settings_eval_reset(struct Depsgraph *depsgraph,
5015 ParticleSettings *particle_settings)
5016 {
5017 DEG_debug_print_eval(depsgraph, __func__, particle_settings->id.name, particle_settings);
5018 particle_settings->id.recalc |= ID_RECALC_PSYS_RESET;
5019 }
5020
BKE_particle_system_eval_init(struct Depsgraph * depsgraph,Object * object)5021 void BKE_particle_system_eval_init(struct Depsgraph *depsgraph, Object *object)
5022 {
5023 DEG_debug_print_eval(depsgraph, __func__, object->id.name, object);
5024 for (ParticleSystem *psys = object->particlesystem.first; psys != NULL; psys = psys->next) {
5025 psys->recalc |= (psys->part->id.recalc & ID_RECALC_PSYS_ALL);
5026 }
5027 }
5028