1 files changed, 352 insertions, 67 deletions

diff --git a/source/blender/blenkernel/intern/particle_system.c b/source/blender/blenkernel/intern/particle_system.c
index 90889d7c09e..090082b333d 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,36 @@ 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];
-
-	/* 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_evaluate_func(BVHTree *tree, ParticleSystem **psys, float co[3], SPHRangeData *pfr, float interaction_radius, BVHTree_RangeQuery callback)
+{
+	int i;
+
+	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 +2449,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 +2503,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 +2517,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; i<pfr.tot_neighbors; i++, pfn++) {
@@ -2593,14 +2597,207 @@ 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));
+		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 UNUSED(squared_dist))
+{
+	SPHRangeData *pfr = (SPHRangeData *)userdata;
+	ParticleData *npa = pfr->npsys->particles + index;
+	float q;
+	float qfac = 21.0f / (256.f * (float)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 UNUSED(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 * (float)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 +2818,46 @@ 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)
 {
@@ -3781,18 +4007,19 @@ static void save_hair(ParticleSimulationData *sim, float UNUSED(cfra))
 
 /* Code for an adaptive time step based on the Courant-Friedrichs-Lewy
  * condition. */
-#define MIN_TIMESTEP 1.0f / 101.0f
+static const float MIN_TIMESTEP = 1.0f / 101.0f;
 /* Tolerance of 1.5 means the last subframe neither favors growing nor
  * shrinking (e.g if it were 1.3, the last subframe would tend to be too
  * small). */
-#define TIMESTEP_EXPANSION_TOLERANCE 1.5f
+static const float TIMESTEP_EXPANSION_FACTOR = 0.1f;
+static const float TIMESTEP_EXPANSION_TOLERANCE = 1.5f;
 
 /* Calculate the speed of the particle relative to the local scale of the
  * simulation. This should be called once per particle during a simulation
  * step, after the velocity has been updated. element_size defines the scale of
  * the simulation, and is typically the distance to neighboring particles. */
 static void update_courant_num(ParticleSimulationData *sim, ParticleData *pa,
-	float dtime, SPHData *sphdata)
+                               float dtime, SPHData *sphdata)
 {
 	float relative_vel[3];
 	float speed;
@@ -3802,14 +4029,31 @@ static void update_courant_num(ParticleSimulationData *sim, ParticleData *pa,
 	if (sim->courant_num < speed * dtime / sphdata->element_size)
 		sim->courant_num = speed * dtime / sphdata->element_size;
 }
+static float get_base_time_step(ParticleSettings *part)
+{
+	return 1.0f / (float) (part->subframes + 1);
+}
 /* Update time step size to suit current conditions. */
 static float update_timestep(ParticleSystem *psys, ParticleSimulationData *sim, float t_frac)
 {
+	float dt_target;
 	if (sim->courant_num == 0.0f)
-		psys->dt_frac = 1.0f;
+		dt_target = 1.0f;
+	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 */
+	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 increase slowly. */
+	if (dt_target > psys->dt_frac)
+		psys->dt_frac = interpf(dt_target, psys->dt_frac, TIMESTEP_EXPANSION_FACTOR);
 	else
-		psys->dt_frac *= (psys->part->courant_target / sim->courant_num);
-	CLAMP(psys->dt_frac, MIN_TIMESTEP, 1.0f);
+		psys->dt_frac = dt_target;
 
 	/* Sync with frame end if it's close. */
 	if (t_frac == 1.0f)
@@ -3973,31 +4217,72 @@ 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);
 
-			#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);
+			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);
+					/* 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);  
+					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);
 
-				#pragma omp critical
-				if (part->time_flag & PART_TIME_AUTOSF)
-					update_courant_num(sim, pa, dtime, &sphdata);
 			}
+			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);
+				}
+
+				/* do global forces & effectors */
+				#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);
+
+					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_springs_modify(psys, timestep);
+					#pragma omp critical
+					if (part->time_flag & PART_TIME_AUTOSF)
+						update_courant_num(sim, pa, dtime, &sphdata);
+				}
+			}
 
-			sph_solver_finalise(&sphdata);
+			psys_sph_finalise(&sphdata);
 			break;
 		}
 	}
@@ -4309,14 +4594,14 @@ static void system_step(ParticleSimulationData *sim, float cfra)
 
 		if (!(part->time_flag & PART_TIME_AUTOSF)) {
 			/* Constant time step */
-			psys->dt_frac = 1.0f / (float) (part->subframes + 1);
+			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. */
-			psys->dt_frac = 1.0;
+			/* 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;
 		}