[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [52989] trunk/blender: Adding a new SPH solver that is more physically accurate.
Alex Fraser
alex at phatcore.com
Fri Dec 14 05:57:32 CET 2012
Revision: 52989
http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=52989
Author: z0r
Date: 2012-12-14 04:57:26 +0000 (Fri, 14 Dec 2012)
Log Message:
-----------
Adding a new SPH solver that is more physically accurate. See patch #29681
http://projects.blender.org/tracker/index.php?func=detail&aid=29681&group_id=9&atid=127
The solver was mostly implemented by John Mansour at VPAC, with help from me and with funding from the AutoCRC. The SPH formulation is due to Gingold and Monaghan, and the smoothing kernel is due to Wendland.
This solver does not replace the old one; it is available as an option. Note that the new solver uses different units than the old one. The patch page has a couple of attachments that can be used to test the new solver, particularly sphclassical_dam_s0.01_grav.blend (ignore the earlier tests). The simulation in that file compares well with a physical experimental dam break; details in a paper by Changhong Hu and Makoto Sueyoshi, also referred to on that page.
Modified Paths:
--------------
trunk/blender/release/scripts/startup/bl_ui/properties_particle.py
trunk/blender/source/blender/blenkernel/BKE_particle.h
trunk/blender/source/blender/blenkernel/intern/particle_system.c
trunk/blender/source/blender/makesdna/DNA_particle_types.h
trunk/blender/source/blender/makesrna/intern/rna_particle.c
Modified: trunk/blender/release/scripts/startup/bl_ui/properties_particle.py
===================================================================
--- trunk/blender/release/scripts/startup/bl_ui/properties_particle.py 2012-12-14 04:38:52 UTC (rev 52988)
+++ trunk/blender/release/scripts/startup/bl_ui/properties_particle.py 2012-12-14 04:57:26 UTC (rev 52989)
@@ -506,7 +506,11 @@
fluid = part.fluid
split = layout.split()
+ sub = split.row()
+ sub.prop(fluid, "solver", expand=True)
+ split = layout.split()
+
col = split.column()
col.label(text="Fluid properties:")
col.prop(fluid, "stiffness", text="Stiffness")
@@ -516,13 +520,14 @@
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,28 +537,38 @@
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()
sub = split.column()
Modified: trunk/blender/source/blender/blenkernel/BKE_particle.h
===================================================================
--- trunk/blender/source/blender/blenkernel/BKE_particle.h 2012-12-14 04:38:52 UTC (rev 52988)
+++ trunk/blender/source/blender/blenkernel/BKE_particle.h 2012-12-14 04:57:26 UTC (rev 52989)
@@ -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 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 @@
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 @@
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,
Modified: trunk/blender/source/blender/blenkernel/intern/particle_system.c
===================================================================
--- trunk/blender/source/blender/blenkernel/intern/particle_system.c 2012-12-14 04:38:52 UTC (rev 52988)
+++ trunk/blender/source/blender/blenkernel/intern/particle_system.c 2012-12-14 04:57:26 UTC (rev 52989)
@@ -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 @@
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 @@
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 @@
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 @@
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 @@
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];
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list