diff --git a/release/scripts/startup/bl_ui/properties_particle.py b/release/scripts/startup/bl_ui/properties_particle.py index 75154550bcb..2c2ced5db0c 100644 --- a/release/scripts/startup/bl_ui/properties_particle.py +++ b/release/scripts/startup/bl_ui/properties_particle.py @@ -505,6 +505,10 @@ class PARTICLE_PT_physics(ParticleButtonsPanel, Panel): if part.physics_type == 'FLUID': fluid = part.fluid + split = layout.split() + sub = split.row() + sub.prop(fluid, "solver", expand=True) + split = layout.split() col = split.column() @@ -516,13 +520,14 @@ class PARTICLE_PT_physics(ParticleButtonsPanel, Panel): col = split.column() col.label(text="Advanced:") - sub = col.row() - sub.prop(fluid, "repulsion", slider=fluid.factor_repulsion) - sub.prop(fluid, "factor_repulsion", text="") + if fluid.solver == 'DDR': + sub = col.row() + sub.prop(fluid, "repulsion", slider=fluid.factor_repulsion) + sub.prop(fluid, "factor_repulsion", text="") - sub = col.row() - sub.prop(fluid, "stiff_viscosity", slider=fluid.factor_stiff_viscosity) - sub.prop(fluid, "factor_stiff_viscosity", text="") + sub = col.row() + sub.prop(fluid, "stiff_viscosity", slider=fluid.factor_stiff_viscosity) + sub.prop(fluid, "factor_stiff_viscosity", text="") sub = col.row() sub.prop(fluid, "fluid_radius", slider=fluid.factor_radius) @@ -532,27 +537,37 @@ class PARTICLE_PT_physics(ParticleButtonsPanel, Panel): sub.prop(fluid, "rest_density", slider=fluid.factor_density) sub.prop(fluid, "factor_density", text="") - split = layout.split() + if fluid.solver == 'CLASSICAL': + # With the classical solver, it is possible to calculate the + # spacing between particles when the fluid is at rest. This + # makes it easier to set stable initial conditions. + particle_volume = part.mass / fluid.rest_density + spacing = pow(particle_volume, 1/3) + sub = col.row() + sub.label(text="Spacing: %g" % spacing) - col = split.column() - col.label(text="Springs:") - col.prop(fluid, "spring_force", text="Force") - col.prop(fluid, "use_viscoelastic_springs") - sub = col.column(align=True) - sub.active = fluid.use_viscoelastic_springs - sub.prop(fluid, "yield_ratio", slider=True) - sub.prop(fluid, "plasticity", slider=True) + elif fluid.solver == 'DDR': + split = layout.split() - col = split.column() - col.label(text="Advanced:") - sub = col.row() - sub.prop(fluid, "rest_length", slider=fluid.factor_rest_length) - sub.prop(fluid, "factor_rest_length", text="") - col.label(text="") - sub = col.column() - sub.active = fluid.use_viscoelastic_springs - sub.prop(fluid, "use_initial_rest_length") - sub.prop(fluid, "spring_frames", text="Frames") + col = split.column() + col.label(text="Springs:") + col.prop(fluid, "spring_force", text="Force") + col.prop(fluid, "use_viscoelastic_springs") + sub = col.column(align=True) + sub.active = fluid.use_viscoelastic_springs + sub.prop(fluid, "yield_ratio", slider=True) + sub.prop(fluid, "plasticity", slider=True) + + col = split.column() + col.label(text="Advanced:") + sub = col.row() + sub.prop(fluid, "rest_length", slider=fluid.factor_rest_length) + sub.prop(fluid, "factor_rest_length", text="") + col.label(text="") + sub = col.column() + sub.active = fluid.use_viscoelastic_springs + sub.prop(fluid, "use_initial_rest_length") + sub.prop(fluid, "spring_frames", text="Frames") elif part.physics_type == 'KEYED': split = layout.split() diff --git a/source/blender/blenkernel/BKE_particle.h b/source/blender/blenkernel/BKE_particle.h index bdaffef6818..f15ad296e4a 100644 --- a/source/blender/blenkernel/BKE_particle.h +++ b/source/blender/blenkernel/BKE_particle.h @@ -20,7 +20,9 @@ * * The Original Code is: all of this file. * - * Contributor(s): none yet. + * Adaptive time step + * Classical SPH + * Copyright 2011-2012 AutoCRC * * ***** END GPL LICENSE BLOCK ***** */ @@ -58,6 +60,7 @@ struct RNG; struct SurfaceModifierData; struct BVHTreeRay; struct BVHTreeRayHit; +struct EdgeHash; #define PARTICLE_P ParticleData * pa; int p #define LOOP_PARTICLES for (p = 0, pa = psys->particles; p < psys->totpart; p++, pa++) @@ -85,6 +88,24 @@ typedef struct ParticleSimulationData { float courant_num; } ParticleSimulationData; +typedef struct SPHData { + ParticleSystem *psys[10]; + ParticleData *pa; + float mass; + struct EdgeHash *eh; + float *gravity; + float hfac; + /* 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; + typedef struct ParticleTexture { float ivel; /* used in reset */ float time, life, exist, size; /* used in init */ @@ -283,6 +304,10 @@ float psys_get_child_size(struct ParticleSystem *psys, struct ChildParticle *cpa void psys_get_particle_on_path(struct ParticleSimulationData *sim, int pa_num, struct ParticleKey *state, int vel); int psys_get_particle_state(struct ParticleSimulationData *sim, int p, struct ParticleKey *state, int always); +void psys_sph_init(struct ParticleSimulationData *sim, struct SPHData *sphdata); +void psys_sph_finalise(struct SPHData *sphdata); +void psys_sph_density(struct BVHTree *tree, struct SPHData* data, float co[3], float vars[2]); + /* for anim.c */ void psys_get_dupli_texture(struct ParticleSystem *psys, struct ParticleSettings *part, struct ParticleSystemModifierData *psmd, struct ParticleData *pa, struct ChildParticle *cpa, diff --git a/source/blender/blenkernel/intern/particle_system.c b/source/blender/blenkernel/intern/particle_system.c index d5fabf60079..7c95265b27b 100644 --- a/source/blender/blenkernel/intern/particle_system.c +++ b/source/blender/blenkernel/intern/particle_system.c @@ -23,7 +23,8 @@ * Contributor(s): Raul Fernandez Hernandez (Farsthary), Stephen Swhitehorn. * * Adaptive time step - * Copyright 2011 AutoCRC + * Classical SPH + * Copyright 2011-2012 AutoCRC * * ***** END GPL LICENSE BLOCK ***** */ @@ -1889,7 +1890,6 @@ void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, bpa->data.acc[0]=bpa->data.acc[1]=bpa->data.acc[2]=0.0f; } - if (part->type == PART_HAIR) { pa->lifetime = 100.0f; } @@ -2390,33 +2390,34 @@ typedef struct SPHRangeData { SPHNeighbor neighbors[SPH_NEIGHBORS]; int tot_neighbors; - float density, near_density; - float h; + float* data; ParticleSystem *npsys; ParticleData *pa; + float h; float massfac; int use_size; } SPHRangeData; -typedef struct SPHData { - ParticleSystem *psys[10]; - ParticleData *pa; - float mass; - EdgeHash *eh; - float *gravity; - /* Average distance to neighbors (other particles in the support domain), - * for calculating the Courant number (adaptive time step). */ - int pass; - float element_size; - float flow[3]; +static void sph_evaluate_func(BVHTree *tree, ParticleSystem **psys, float co[3], SPHRangeData *pfr, float interaction_radius, BVHTree_RangeQuery callback) { + int i; - /* 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; + pfr->tot_neighbors = 0; + for (i=0; i < 10 && psys[i]; i++) { + pfr->npsys = psys[i]; + pfr->massfac = psys[i]->part->mass; + pfr->use_size = psys[i]->part->flag & PART_SIZEMASS; + + if (tree) { + BLI_bvhtree_range_query(tree, co, interaction_radius, callback, pfr); + break; + } else { + BLI_bvhtree_range_query(psys[i]->bvhtree, co, interaction_radius, callback, pfr); + } + } +} static void sph_density_accum_cb(void *userdata, int index, float squared_dist) { SPHRangeData *pfr = (SPHRangeData *)userdata; @@ -2446,8 +2447,8 @@ static void sph_density_accum_cb(void *userdata, int index, float squared_dist) if (pfr->use_size) q *= npa->size; - pfr->density += q*q; - pfr->near_density += q*q*q; + pfr->data[0] += q*q; + pfr->data[1] += q*q*q; } /* @@ -2500,8 +2501,10 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa float inv_mass = 1.0f/mass; float spring_constant = fluid->spring_k; - - float h = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.f*pa->size : 1.f); /* 4.0 seems to be a pretty good value */ + + /* 4.0 seems to be a pretty good value */ + float interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f); + float h = interaction_radius * sphdata->hfac; float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.f); /* 4.77 is an experimentally determined density factor */ float rest_length = fluid->rest_length * (fluid->flag & SPH_FAC_REST_LENGTH ? 2.588f * pa->size : 1.f); @@ -2512,24 +2515,23 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa float vec[3]; float vel[3]; float co[3]; + float data[2]; + float density, near_density; int i, spring_index, index = pa - psys[0]->particles; - pfr.tot_neighbors = 0; - pfr.density = pfr.near_density = 0.f; + data[0] = data[1] = 0; + pfr.data = data; pfr.h = h; pfr.pa = pa; - 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; + sph_evaluate_func( NULL, psys, state->co, &pfr, interaction_radius, sph_density_accum_cb); - BLI_bvhtree_range_query(psys[i]->bvhtree, state->co, h, sphdata->density_cb, &pfr); - } + density = data[0]; + near_density = data[1]; - pressure = stiffness * (pfr.density - rest_density); - near_pressure = stiffness_near_fac * pfr.near_density; + pressure = stiffness * (density - rest_density); + near_pressure = stiffness_near_fac * near_density; pfn = pfr.neighbors; for (i=0; ibuoyancy > 0.f && gravity) - madd_v3_v3fl(force, gravity, fluid->buoyancy * (pfr.density-rest_density)); + madd_v3_v3fl(force, gravity, fluid->buoyancy * (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_solver_init(ParticleSimulationData *sim, SPHData *sphdata) +/* powf is really slow for raising to integer powers. */ +MINLINE float pow2(float x) { + return x * x; +} +MINLINE float pow3(float x) { + return pow2(x) * x; +} +MINLINE float pow4(float x) { + return pow2(pow2(x)); +} +MINLINE float pow7(float x) { + return pow2(pow3(x)) * x; +} + +static void sphclassical_density_accum_cb(void *userdata, int index, float squared_dist) +{ + SPHRangeData *pfr = (SPHRangeData *)userdata; + ParticleData *npa = pfr->npsys->particles + index; + float q; + float qfac = 21.0f / (256.f * M_PI); + float rij, rij_h; + float vec[3]; + + /* Exclude particles that are more than 2h away. Can't use squared_dist here + * because it is not accurate enough. Use current state, i.e. the output of + * basic_integrate() - z0r */ + sub_v3_v3v3(vec, npa->state.co, pfr->pa->state.co); + rij = len_v3(vec); + rij_h = rij / pfr->h; + if (rij_h > 2.0f) + return; + + /* Smoothing factor. Utilise the Wendland kernel. gnuplot: + * q1(x) = (2.0 - x)**4 * ( 1.0 + 2.0 * x) + * plot [0:2] q1(x) */ + q = qfac / pow3(pfr->h) * pow4(2.0f - rij_h) * ( 1.0f + 2.0f * rij_h); + q *= pfr->massfac; + + if (pfr->use_size) + q *= pfr->pa->size; + + pfr->data[0] += q; + pfr->data[1] += q / npa->sphdensity; +} + +static void sphclassical_neighbour_accum_cb(void *userdata, int index, float squared_dist) +{ + SPHRangeData *pfr = (SPHRangeData *)userdata; + ParticleData *npa = pfr->npsys->particles + index; + float rij, rij_h; + float vec[3]; + + if (pfr->tot_neighbors >= SPH_NEIGHBORS) + return; + + /* Exclude particles that are more than 2h away. Can't use squared_dist here + * because it is not accurate enough. Use current state, i.e. the output of + * basic_integrate() - z0r */ + sub_v3_v3v3(vec, npa->state.co, pfr->pa->state.co); + rij = len_v3(vec); + rij_h = rij / pfr->h; + if (rij_h > 2.0f) + return; + + pfr->neighbors[pfr->tot_neighbors].index = index; + pfr->neighbors[pfr->tot_neighbors].psys = pfr->npsys; + pfr->tot_neighbors++; +} +static void sphclassical_force_cb(void *sphdata_v, ParticleKey *state, float *force, float *UNUSED(impulse)) +{ + SPHData *sphdata = (SPHData *)sphdata_v; + ParticleSystem **psys = sphdata->psys; + ParticleData *pa = sphdata->pa; + SPHFluidSettings *fluid = psys[0]->part->fluid; + SPHRangeData pfr; + SPHNeighbor *pfn; + float *gravity = sphdata->gravity; + + float dq, u, rij, dv[3]; + float pressure, npressure; + + float visc = fluid->viscosity_omega; + + float interaction_radius; + float h, hinv; + /* 4.77 is an experimentally determined density factor */ + float rest_density = fluid->rest_density * (fluid->flag & SPH_FAC_DENSITY ? 4.77f : 1.0f); + + float stiffness = fluid->stiffness_k; + + ParticleData *npa; + float vec[3]; + float co[3]; + float pressureTerm; + + int i; + + float qfac2 = 42.0f / (256.0f * M_PI); + float rij_h; + + /* 4.0 here is to be consistent with previous formulation/interface */ + interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * pa->size : 1.0f); + h = interaction_radius * sphdata->hfac; + hinv = 1.0f / h; + + pfr.h = h; + pfr.pa = pa; + + sph_evaluate_func(NULL, psys, state->co, &pfr, interaction_radius, sphclassical_neighbour_accum_cb); + pressure = stiffness * (pow7(pa->sphdensity / rest_density) - 1.0f); + + /* multiply by mass so that we return a force, not accel */ + qfac2 *= sphdata->mass / pow3(pfr.h); + + pfn = pfr.neighbors; + for (i = 0; i < pfr.tot_neighbors; i++, pfn++) { + npa = pfn->psys->particles + pfn->index; + if (npa == pa) { + /* we do not contribute to ourselves */ + continue; + } + + /* Find vector to neighbour. Exclude particles that are more than 2h + * away. Can't use current state here because it may have changed on + * another thread - so do own mini integration. Unlike basic_integrate, + * SPH integration depends on neighbouring particles. - z0r */ + madd_v3_v3v3fl(co, npa->prev_state.co, npa->prev_state.vel, state->time); + sub_v3_v3v3(vec, co, state->co); + rij = normalize_v3(vec); + rij_h = rij / pfr.h; + if (rij_h > 2.0f) + continue; + + npressure = stiffness * (pow7(npa->sphdensity / rest_density) - 1.0f); + + /* First derivative of smoothing factor. Utilise the Wendland kernel. + * gnuplot: + * q2(x) = 2.0 * (2.0 - x)**4 - 4.0 * (2.0 - x)**3 * (1.0 + 2.0 * x) + * plot [0:2] q2(x) + * Particles > 2h away are excluded above. */ + dq = qfac2 * (2.0f * pow4(2.0f - rij_h) - 4.0f * pow3(2.0f - rij_h) * (1.0f + 2.0f * rij_h) ); + + if (pfn->psys->part->flag & PART_SIZEMASS) + dq *= npa->size; + + pressureTerm = pressure / pow2(pa->sphdensity) + npressure / pow2(npa->sphdensity); + + /* Note that 'minus' is removed, because vec = vecBA, not vecAB. + * This applies to the viscosity calculation below, too. */ + madd_v3_v3fl(force, vec, pressureTerm * dq); + + /* Viscosity */ + if (visc > 0.0f) { + sub_v3_v3v3(dv, npa->prev_state.vel, pa->prev_state.vel); + u = dot_v3v3(vec, dv); + /* Apply parameters */ + u *= -dq * hinv * visc / (0.5f * npa->sphdensity + 0.5f * pa->sphdensity); + madd_v3_v3fl(force, vec, u); + } + } + + /* Artificial buoyancy force in negative gravity direction */ + if (fluid->buoyancy > 0.f && gravity) + madd_v3_v3fl(force, gravity, fluid->buoyancy * (pa->sphdensity - rest_density)); + + if (sphdata->pass == 0 && psys[0]->part->time_flag & PART_TIME_AUTOSF) + sph_particle_courant(sphdata, &pfr); + sphdata->pass++; +} + +static void sphclassical_calc_dens(ParticleData *pa, float UNUSED(dfra), SPHData *sphdata){ + ParticleSystem **psys = sphdata->psys; + SPHFluidSettings *fluid = psys[0]->part->fluid; + /* 4.0 seems to be a pretty good value */ + float interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f); + SPHRangeData pfr; + float data[2]; + + data[0] = 0; + data[1] = 0; + pfr.data = data; + pfr.h = interaction_radius * sphdata->hfac; + pfr.pa = pa; + + sph_evaluate_func( NULL, psys, pa->state.co, &pfr, interaction_radius, sphclassical_density_accum_cb); + pa->sphdensity = MIN2(MAX2(data[0], fluid->rest_density * 0.9f), fluid->rest_density * 1.1f); +} + +void psys_sph_init(ParticleSimulationData *sim, SPHData *sphdata) { ParticleTarget *pt; int i; @@ -2621,17 +2811,44 @@ static void sph_solver_init(ParticleSimulationData *sim, SPHData *sphdata) sphdata->pa = NULL; sphdata->mass = 1.0f; - sphdata->force_cb = sph_force_cb; - sphdata->density_cb = sph_density_accum_cb; + if (sim->psys->part->fluid->solver == SPH_SOLVER_DDR) { + sphdata->force_cb = sph_force_cb; + sphdata->density_cb = sph_density_accum_cb; + sphdata->hfac = 1.0f; + } else { + /* SPH_SOLVER_CLASSICAL */ + sphdata->force_cb = sphclassical_force_cb; + sphdata->density_cb = sphclassical_density_accum_cb; + sphdata->hfac = 0.5f; + } + } -static void sph_solver_finalise(SPHData *sphdata) +void psys_sph_finalise(SPHData *sphdata) { if (sphdata->eh) { BLI_edgehash_free(sphdata->eh, NULL); sphdata->eh = NULL; } } +/* Sample the density field at a point in space. */ +void psys_sph_density(BVHTree *tree, SPHData *sphdata, float co[3], float vars[2]) { + ParticleSystem **psys = sphdata->psys; + SPHFluidSettings *fluid = psys[0]->part->fluid; + /* 4.0 seems to be a pretty good value */ + float interaction_radius = fluid->radius * (fluid->flag & SPH_FAC_RADIUS ? 4.0f * psys[0]->part->size : 1.0f); + SPHRangeData pfr; + float density[2]; + + density[0] = density[1] = 0.0f; + pfr.data = density; + pfr.h = interaction_radius*sphdata->hfac; + + sph_evaluate_func(tree, psys, co, &pfr, interaction_radius, sphdata->density_cb); + + vars[0] = pfr.data[0]; + vars[1] = pfr.data[1]; +} static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float dfra, SPHData *sphdata) { @@ -3815,15 +4032,14 @@ static float update_timestep(ParticleSystem *psys, ParticleSimulationData *sim, else dt_target = psys->dt_frac * (psys->part->courant_target / sim->courant_num); - /* Make sure the time step is reasonable. - * For some reason, the CLAMP macro doesn't work here. The time step becomes - * too large. - z0r */ + /* Make sure the time step is reasonable. For some reason, the CLAMP macro + * doesn't work here. The time step becomes too large. - z0r */ if (dt_target < MIN_TIMESTEP) dt_target = MIN_TIMESTEP; else if (dt_target > get_base_time_step(psys->part)) dt_target = get_base_time_step(psys->part); - /* Decrease time step instantly, but expand slowly. */ + /* Decrease time step instantly, but increase slowly. */ if (dt_target > psys->dt_frac) psys->dt_frac = interpf(dt_target, psys->dt_frac, TIMESTEP_EXPANSION_FACTOR); else @@ -3991,31 +4207,71 @@ static void dynamics_step(ParticleSimulationData *sim, float cfra) case PART_PHYS_FLUID: { SPHData sphdata; - sph_solver_init(sim, &sphdata); + ParticleSettings *part = sim->psys->part; + psys_sph_init(sim, &sphdata); + + if (part->fluid->flag & SPH_SOLVER_DDR) { + /* Apply SPH forces using double-density relaxation algorithm + * (Clavat et. al.) */ + #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, &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 */ + basic_rotate(part, pa, pa->state.time, timestep); + + #pragma omp critical + if (part->time_flag & PART_TIME_AUTOSF) + update_courant_num(sim, pa, dtime, &sphdata); + } + + sph_springs_modify(psys, timestep); + + } else { + /* SPH_SOLVER_CLASSICAL */ + /* Apply SPH forces using classical algorithm (due to Gingold + * and Monaghan). Note that, unlike double-density relaxation, + * this algorthim is separated into distinct loops. */ + + #pragma omp parallel for firstprivate (sphdata) private (pa) schedule(dynamic,5) + LOOP_DYNAMIC_PARTICLES { + basic_integrate(sim, p, pa->state.time, cfra); + } + + /* calculate summation density */ + #pragma omp parallel for firstprivate (sphdata) private (pa) schedule(dynamic,5) + LOOP_DYNAMIC_PARTICLES { + sphclassical_calc_dens(pa, pa->state.time, &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); + #pragma omp parallel for firstprivate (sphdata) private (pa) schedule(dynamic,5) + LOOP_DYNAMIC_PARTICLES { + /* actual fluids calculations */ + sph_integrate(sim, pa, pa->state.time, &sphdata); - /* actual fluids calculations */ - sph_integrate(sim, pa, pa->state.time, &sphdata); - - if (sim->colliders) - collision_check(sim, p, pa->state.time, cfra); + 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 */ - basic_rotate(part, pa, pa->state.time, timestep); + /* 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, &sphdata); + #pragma omp critical + if (part->time_flag & PART_TIME_AUTOSF) + update_courant_num(sim, pa, dtime, &sphdata); + } } - sph_springs_modify(psys, timestep); - - sph_solver_finalise(&sphdata); + psys_sph_finalise(&sphdata); break; } } @@ -4325,18 +4581,16 @@ static void system_step(ParticleSimulationData *sim, float cfra) if ((int)cfra == startframe && part->sta < startframe) totframesback = (startframe - (int)part->sta); - /* Initialise time step. This may change if automatic subframes are - * enabled. */ if (!(part->time_flag & PART_TIME_AUTOSF)) { /* Constant time step */ psys->dt_frac = get_base_time_step(part); } else if ((int)cfra == startframe) { - /* Variable time step; use a very conservative value at the start. - * If it doesn't need to be so small, it will quickly grow. */ + /* Variable time step; initialise to subframes */ psys->dt_frac = get_base_time_step(part); } else if (psys->dt_frac < MIN_TIMESTEP) { + /* Variable time step; subsequent frames */ psys->dt_frac = MIN_TIMESTEP; } diff --git a/source/blender/makesdna/DNA_particle_types.h b/source/blender/makesdna/DNA_particle_types.h index 5952aa8afb0..1cf16fb52b3 100644 --- a/source/blender/makesdna/DNA_particle_types.h +++ b/source/blender/makesdna/DNA_particle_types.h @@ -116,6 +116,9 @@ typedef struct ParticleData { float size; /* size and multiplier so that we can update size when ever */ + float sphdensity; /* density of sph particle */ + int pad; + int hair_index; short flag; short alive; /* the life state of a particle */ @@ -130,6 +133,8 @@ typedef struct SPHFluidSettings { float stiffness_k, stiffness_knear, rest_density; float buoyancy; int flag, spring_frames; + short solver; + short pad[3]; } SPHFluidSettings; /* fluid->flag */ @@ -141,6 +146,10 @@ typedef struct SPHFluidSettings { #define SPH_FAC_VISCOSITY 32 #define SPH_FAC_REST_LENGTH 64 +/* fluid->solver (numerical ID field, not bitfield) */ +#define SPH_SOLVER_DDR 0 +#define SPH_SOLVER_CLASSICAL 1 + typedef struct ParticleSettings { ID id; struct AnimData *adt; diff --git a/source/blender/makesrna/intern/rna_particle.c b/source/blender/makesrna/intern/rna_particle.c index 96b6197cf69..93f940b1aa3 100644 --- a/source/blender/makesrna/intern/rna_particle.c +++ b/source/blender/makesrna/intern/rna_particle.c @@ -1112,12 +1112,25 @@ static void rna_def_fluid_settings(BlenderRNA *brna) { StructRNA *srna; PropertyRNA *prop; - + + static EnumPropertyItem sph_solver_items[] = { + {SPH_SOLVER_DDR, "DDR", 0, "Double-Density", "An artistic solver with strong surface tension effects (original)"}, + {SPH_SOLVER_CLASSICAL, "CLASSICAL", 0, "Classical", "A more physically-accurate solver"}, + {0, NULL, 0, NULL, NULL} + }; + srna = RNA_def_struct(brna, "SPHFluidSettings", NULL); RNA_def_struct_path_func(srna, "rna_SPHFluidSettings_path"); RNA_def_struct_ui_text(srna, "SPH Fluid Settings", "Settings for particle fluids physics"); - + /* Fluid settings */ + prop = RNA_def_property(srna, "solver", PROP_ENUM, PROP_NONE); + RNA_def_property_enum_sdna(prop, NULL, "solver"); + RNA_def_property_clear_flag(prop, PROP_ANIMATABLE); + RNA_def_property_enum_items(prop, sph_solver_items); + RNA_def_property_ui_text(prop, "SPH Solver", "The code used to calculate internal forces on particles"); + RNA_def_property_update(prop, 0, "rna_Particle_reset"); + prop = RNA_def_property(srna, "spring_force", PROP_FLOAT, PROP_NONE); RNA_def_property_float_sdna(prop, NULL, "spring_k"); RNA_def_property_range(prop, 0.0f, 100.0f); @@ -1186,7 +1199,7 @@ static void rna_def_fluid_settings(BlenderRNA *brna) /* Double density relaxation */ prop = RNA_def_property(srna, "stiffness", PROP_FLOAT, PROP_NONE); RNA_def_property_float_sdna(prop, NULL, "stiffness_k"); - RNA_def_property_range(prop, 0.0f, 100.0f); + RNA_def_property_range(prop, 0.0f, 100000.0f); RNA_def_property_ui_range(prop, 0.0f, 10.0f, 1, 3); RNA_def_property_ui_text(prop, "Stiffness", "How incompressible the fluid is"); RNA_def_property_update(prop, 0, "rna_Particle_reset"); @@ -1201,7 +1214,7 @@ static void rna_def_fluid_settings(BlenderRNA *brna) prop = RNA_def_property(srna, "rest_density", PROP_FLOAT, PROP_NONE); RNA_def_property_float_sdna(prop, NULL, "rest_density"); - RNA_def_property_range(prop, 0.0f, 100.0f); + RNA_def_property_range(prop, 0.0f, 10000.0f); RNA_def_property_ui_range(prop, 0.0f, 2.0f, 1, 3); RNA_def_property_ui_text(prop, "Rest Density", "Fluid rest density"); RNA_def_property_update(prop, 0, "rna_Particle_reset"); @@ -2281,7 +2294,7 @@ static void rna_def_particle_settings(BlenderRNA *brna) /* physical properties */ prop = RNA_def_property(srna, "mass", PROP_FLOAT, PROP_NONE); - RNA_def_property_range(prop, 0.001f, 100000.0f); + RNA_def_property_range(prop, 0.00000001f, 100000.0f); RNA_def_property_ui_range(prop, 0.01, 100, 1, 3); RNA_def_property_ui_text(prop, "Mass", "Mass of the particles"); RNA_def_property_update(prop, 0, "rna_Particle_reset");