diff options
author | Campbell Barton <ideasman42@gmail.com> | 2012-01-03 06:16:52 +0400 |
---|---|---|
committer | Campbell Barton <ideasman42@gmail.com> | 2012-01-03 06:16:52 +0400 |
commit | c0eec8f379bb0c265e7bf2ad74397c22836bae3c (patch) | |
tree | 6502b4e6eaf44a848bebe6406aacd825738877d9 /source/blender/blenkernel/intern | |
parent | 68c4368001db49bd3daab96352ca3b4fb7d04e1c (diff) | |
parent | 84437bb5e9342e67f448e26adb0c4b0811f3ef5c (diff) |
svn merge ^/trunk/blender -r43062:43085
Diffstat (limited to 'source/blender/blenkernel/intern')
-rw-r--r-- | source/blender/blenkernel/intern/particle_system.c | 157 |
1 files changed, 99 insertions, 58 deletions
diff --git a/source/blender/blenkernel/intern/particle_system.c b/source/blender/blenkernel/intern/particle_system.c index 45b0fb80944..7e69e67fca7 100644 --- a/source/blender/blenkernel/intern/particle_system.c +++ b/source/blender/blenkernel/intern/particle_system.c @@ -39,6 +39,10 @@ #include <math.h> #include <string.h> +#ifdef _OPENMP +#include <omp.h> +#endif + #include "MEM_guardedalloc.h" #include "DNA_anim_types.h" @@ -2302,6 +2306,7 @@ static EdgeHash *sph_springhash_build(ParticleSystem *psys) return springhash; } +#define SPH_NEIGHBORS 512 typedef struct SPHNeighbor { ParticleSystem *psys; @@ -2309,7 +2314,7 @@ typedef struct SPHNeighbor } SPHNeighbor; typedef struct SPHRangeData { - SPHNeighbor neighbors[128]; + SPHNeighbor neighbors[SPH_NEIGHBORS]; int tot_neighbors; float density, near_density; @@ -2320,10 +2325,6 @@ typedef struct SPHRangeData float massfac; int use_size; - - /* Same as SPHData::element_size */ - float element_size; - float flow[3]; } SPHRangeData; typedef struct SPHData { ParticleSystem *psys[10]; @@ -2333,9 +2334,15 @@ typedef struct SPHData { float *gravity; /* Average distance to neighbours (other particles in the support domain), for calculating the Courant number (adaptive time step). */ + int pass; float element_size; float flow[3]; + + /* Integrator callbacks. This allows different SPH implementations. */ + void (*force_cb) (void *sphdata_v, ParticleKey *state, float *force, float *impulse); + void (*density_cb) (void *rangedata_v, int index, float squared_dist); }SPHData; + static void sph_density_accum_cb(void *userdata, int index, float squared_dist) { SPHRangeData *pfr = (SPHRangeData *)userdata; @@ -2346,11 +2353,13 @@ static void sph_density_accum_cb(void *userdata, int index, float squared_dist) if(npa == pfr->pa || squared_dist < FLT_EPSILON) return; - /* Ugh! One particle has over 128 neighbors! Really shouldn't happen, - * but even if it does it shouldn't do any terrible harm if all are - * not taken into account - jahka + /* Ugh! One particle has too many neighbors! If some aren't taken into + * account, the forces will be biased by the tree search order. This + * effectively adds enery to the system, and results in a churning motion. + * But, we have to stop somewhere, and it's not the end of the world. + * - jahka and z0r */ - if(pfr->tot_neighbors >= 128) + if(pfr->tot_neighbors >= SPH_NEIGHBORS) return; pfr->neighbors[pfr->tot_neighbors].index = index; @@ -2360,15 +2369,38 @@ static void sph_density_accum_cb(void *userdata, int index, float squared_dist) dist = sqrtf(squared_dist); q = (1.f - dist/pfr->h) * pfr->massfac; - add_v3_v3(pfr->flow, npa->state.vel); - pfr->element_size += dist; - if(pfr->use_size) q *= npa->size; pfr->density += q*q; pfr->near_density += q*q*q; } +/* + * Find the Courant number for an SPH particle (used for adaptive time step). + */ +static void sph_particle_courant(SPHData *sphdata, SPHRangeData *pfr) { + ParticleData *pa, *npa; + int i; + float flow[3], offset[3], dist; + + flow[0] = flow[1] = flow[2] = 0.0f; + dist = 0.0f; + if (pfr->tot_neighbors > 0) { + pa = pfr->pa; + for (i=0; i < pfr->tot_neighbors; i++) { + npa = pfr->neighbors[i].psys->particles + pfr->neighbors[i].index; + sub_v3_v3v3(offset, pa->prev_state.co, npa->prev_state.co); + dist += len_v3(offset); + add_v3_v3(flow, npa->prev_state.vel); + } + dist += sphdata->psys[0]->part->fluid->radius; // TODO: remove this? - z0r + sphdata->element_size = dist / pfr->tot_neighbors; + mul_v3_v3fl(sphdata->flow, flow, 1.0f / pfr->tot_neighbors); + } else { + sphdata->element_size = MAXFLOAT; + VECCOPY(sphdata->flow, flow); + } +} static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse)) { SPHData *sphdata = (SPHData *)sphdata_v; @@ -2409,24 +2441,14 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa pfr.density = pfr.near_density = 0.f; pfr.h = h; pfr.pa = pa; - pfr.element_size = fluid->radius; - pfr.flow[0] = pfr.flow[1] = pfr.flow[2] = 0.0f; for(i=0; i<10 && psys[i]; i++) { pfr.npsys = psys[i]; pfr.massfac = psys[i]->part->mass*inv_mass; pfr.use_size = psys[i]->part->flag & PART_SIZEMASS; - BLI_bvhtree_range_query(psys[i]->bvhtree, state->co, h, sph_density_accum_cb, &pfr); - } - if (pfr.tot_neighbors > 0) { - pfr.element_size /= pfr.tot_neighbors; - mul_v3_fl(pfr.flow, 1.0f / pfr.tot_neighbors); - } else { - pfr.element_size = MAXFLOAT; + BLI_bvhtree_range_query(psys[i]->bvhtree, state->co, h, sphdata->density_cb, &pfr); } - sphdata->element_size = pfr.element_size; - copy_v3_v3(sphdata->flow, pfr.flow); pressure = stiffness * (pfr.density - rest_density); near_pressure = stiffness_near_fac * pfr.near_density; @@ -2465,6 +2487,7 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa if(spring_constant > 0.f) { /* Viscoelastic spring force */ if (pfn->psys == psys[0] && fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash) { + /* BLI_edgehash_lookup appears to be thread-safe. - z0r */ spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, pfn->index)); if(spring_index) { @@ -2478,7 +2501,9 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa temp_spring.particle_index[1] = pfn->index; temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : rest_length; temp_spring.delete_flag = 0; - + + /* sph_spring_add is not thread-safe. - z0r */ + #pragma omp critical sph_spring_add(psys[0], &temp_spring); } } @@ -2491,29 +2516,52 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa /* Artificial buoyancy force in negative gravity direction */ if (fluid->buoyancy > 0.f && gravity) madd_v3_v3fl(force, gravity, fluid->buoyancy * (pfr.density-rest_density)); + + if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF) + sph_particle_courant(sphdata, &pfr); + sphdata->pass++; } -static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, float *gravity, EdgeHash *springhash, float *element_size, float flow[3]) -{ +static void sph_solver_init(ParticleSimulationData *sim, SPHData *sphdata) { ParticleTarget *pt; int i; + // Add other coupled particle systems. + sphdata->psys[0] = sim->psys; + for(i=1, pt=sim->psys->targets.first; i<10; i++, pt=(pt?pt->next:NULL)) + sphdata->psys[i] = pt ? psys_get_target_system(sim->ob, pt) : NULL; + + if (psys_uses_gravity(sim)) + sphdata->gravity = sim->scene->physics_settings.gravity; + else + sphdata->gravity = NULL; + sphdata->eh = sph_springhash_build(sim->psys); + + // These per-particle values should be overridden later, but just for + // completeness we give them default values now. + sphdata->pa = NULL; + sphdata->mass = 1.0f; + + sphdata->force_cb = sph_force_cb; + sphdata->density_cb = sph_density_accum_cb; +} +static void sph_solver_finalise(SPHData *sphdata) { + if (sphdata->eh) { + BLI_edgehash_free(sphdata->eh, NULL); + sphdata->eh = NULL; + } +} +static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, SPHData *sphdata){ ParticleSettings *part = sim->psys->part; // float timestep = psys_get_timestep(sim); // UNUSED float pa_mass = part->mass * (part->flag & PART_SIZEMASS ? pa->size : 1.f); float dtime = dfra*psys_get_timestep(sim); // int steps = 1; // UNUSED float effector_acceleration[3]; - SPHData sphdata; - sphdata.psys[0] = sim->psys; - for(i=1, pt=sim->psys->targets.first; i<10; i++, pt=(pt?pt->next:NULL)) - sphdata.psys[i] = pt ? psys_get_target_system(sim->ob, pt) : NULL; - - sphdata.pa = pa; - sphdata.gravity = gravity; - sphdata.mass = pa_mass; - sphdata.eh = springhash; + sphdata->pa = pa; + sphdata->mass = pa_mass; + sphdata->pass = 0; //sphdata.element_size and sphdata.flow are set in the callback. /* restore previous state and treat gravity & effectors as external acceleration*/ @@ -2522,9 +2570,7 @@ static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float d copy_particle_key(&pa->state, &pa->prev_state, 0); - integrate_particle(part, pa, dtime, effector_acceleration, sph_force_cb, &sphdata); - *element_size = sphdata.element_size; - copy_v3_v3(flow, sphdata.flow); + integrate_particle(part, pa, dtime, effector_acceleration, sphdata->force_cb, sphdata); } /************************************************/ @@ -3624,15 +3670,15 @@ static void save_hair(ParticleSimulationData *sim, float UNUSED(cfra)) step, after the velocity has been updated. element_size defines the scale of the simulation, and is typically the distance to neighbourning particles. */ void update_courant_num(ParticleSimulationData *sim, ParticleData *pa, - float dtime, float element_size, float flow[3]) + float dtime, SPHData *sphdata) { float relative_vel[3]; float speed; - sub_v3_v3v3(relative_vel, pa->state.vel, flow); + sub_v3_v3v3(relative_vel, pa->prev_state.vel, sphdata->flow); speed = len_v3(relative_vel); - if (sim->courant_num < speed * dtime / element_size) - sim->courant_num = speed * dtime / element_size; + if (sim->courant_num < speed * dtime / sphdata->element_size) + sim->courant_num = speed * dtime / sphdata->element_size; } /* Update time step size to suit current conditions. */ float update_timestep(ParticleSystem *psys, ParticleSimulationData *sim, @@ -3718,11 +3764,11 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra) case PART_PHYS_FLUID: { ParticleTarget *pt = psys->targets.first; - psys_update_particle_bvhtree(psys, psys->cfra); + psys_update_particle_bvhtree(psys, cfra); for(; pt; pt=pt->next) { /* Updating others systems particle tree for fluid-fluid interaction */ if(pt->ob) - psys_update_particle_bvhtree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), psys->cfra); + psys_update_particle_bvhtree(BLI_findlink(&pt->ob->particlesystem, pt->psys-1), cfra); } break; } @@ -3804,37 +3850,32 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra) } case PART_PHYS_FLUID: { - EdgeHash *springhash = sph_springhash_build(psys); - float *gravity = NULL; - float element_size, flow[3]; - - if(psys_uses_gravity(sim)) - gravity = sim->scene->physics_settings.gravity; + SPHData sphdata; + sph_solver_init(sim, &sphdata); + #pragma omp parallel for firstprivate (sphdata) private (pa) schedule(dynamic,5) LOOP_DYNAMIC_PARTICLES { /* do global forces & effectors */ basic_integrate(sim, p, pa->state.time, cfra); /* actual fluids calculations */ - sph_integrate(sim, pa, pa->state.time, gravity, springhash, - &element_size, flow); + sph_integrate(sim, pa, pa->state.time, &sphdata); if(sim->colliders) collision_check(sim, p, pa->state.time, cfra); - /* SPH particles are not physical particles, just interpolation particles, thus rotation has not a direct sense for them */ + /* SPH particles are not physical particles, just interpolation + * particles, thus rotation has not a direct sense for them */ basic_rotate(part, pa, pa->state.time, timestep); + #pragma omp critical if (part->time_flag & PART_TIME_AUTOSF) - update_courant_num(sim, pa, dtime, element_size, flow); + update_courant_num(sim, pa, dtime, &sphdata); } sph_springs_modify(psys, timestep); - if(springhash) { - BLI_edgehash_free(springhash, NULL); - springhash = NULL; - } + sph_solver_finalise(&sphdata); break; } } |