sim: Remove "continue physics" code
[blender.git] / source / blender / blenkernel / intern / particle_system.c
index 941e2990a4b5786e02a03b40e9ba62a5c31bcb62..fda5f6f2ecb8af94f28a35ba2921ae873a1c8a30 100644 (file)
@@ -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 *****
  */
@@ -346,7 +347,7 @@ void psys_calc_dmcache(Object *ob, DerivedMesh *dm, ParticleSystem *psys)
        if (!dm->deformedOnly) {
                /* Will use later to speed up subsurf/derivedmesh */
                LinkNode *node, *nodedmelem, **nodearray;
-               int totdmelem, totelem, i, *origindex;
+               int totdmelem, totelem, i, *origindex, *origindex_poly = NULL;
 
                if (psys->part->from == PART_FROM_VERT) {
                        totdmelem= dm->getNumVerts(dm);
@@ -356,23 +357,38 @@ void psys_calc_dmcache(Object *ob, DerivedMesh *dm, ParticleSystem *psys)
                else { /* FROM_FACE/FROM_VOLUME */
                        totdmelem= dm->getNumTessFaces(dm);
                        totelem= me->totpoly;
-                       origindex= dm->getTessFaceDataArray(dm, CD_ORIGINDEX);
+                       origindex = dm->getTessFaceDataArray(dm, CD_ORIGINDEX);
+                       /* for face lookups we need the poly origindex too */
+                       origindex_poly = dm->getPolyDataArray(dm, CD_ORIGINDEX);
+                       if (origindex_poly == NULL) {
+                               origindex = NULL;
+                       }
                }
        
                nodedmelem= MEM_callocN(sizeof(LinkNode)*totdmelem, "psys node elems");
                nodearray= MEM_callocN(sizeof(LinkNode *)*totelem, "psys node array");
                
-               for (i=0, node=nodedmelem; i<totdmelem; i++, origindex++, node++) {
-                       node->link= SET_INT_IN_POINTER(i);
+               for (i=0, node=nodedmelem; i<totdmelem; i++, node++) {
+                       int origindex_final;
+                       node->link = SET_INT_IN_POINTER(i);
+
+                       /* may be vertex or face origindex */
+                       origindex_final = origindex ? origindex[i] : ORIGINDEX_NONE;
+
+                       /* if we have a poly source, do an index lookup */
+                       if (origindex_poly && origindex_final != ORIGINDEX_NONE) {
+                               origindex_final = origindex_poly[origindex_final];
+                       }
 
-                       if (*origindex != -1) {
-                               if (nodearray[*origindex]) {
+                       if (origindex_final != ORIGINDEX_NONE) {
+                               if (nodearray[origindex_final]) {
                                        /* prepend */
-                                       node->next = nodearray[*origindex];
-                                       nodearray[*origindex]= node;
+                                       node->next = nodearray[origindex_final];
+                                       nodearray[origindex_final] = node;
+                               }
+                               else {
+                                       nodearray[origindex_final] = node;
                                }
-                               else
-                                       nodearray[*origindex]= node;
                        }
                }
                
@@ -452,13 +468,7 @@ static void distribute_grid(DerivedMesh *dm, ParticleSystem *psys)
        mv++;
 
        for (i=1; i<totvert; i++, mv++) {
-               min[0]=MIN2(min[0],mv->co[0]);
-               min[1]=MIN2(min[1],mv->co[1]);
-               min[2]=MIN2(min[2],mv->co[2]);
-
-               max[0]=MAX2(max[0],mv->co[0]);
-               max[1]=MAX2(max[1],mv->co[1]);
-               max[2]=MAX2(max[2],mv->co[2]);
+               minmax_v3v3_v3(min, max, mv->co);
        }
 
        sub_v3_v3v3(delta, max, min);
@@ -512,9 +522,9 @@ static void distribute_grid(DerivedMesh *dm, ParticleSystem *psys)
                        vec[0]/=delta[0];
                        vec[1]/=delta[1];
                        vec[2]/=delta[2];
-                       (pa     + ((int)(vec[0] * (size[0] - 1)) * res +
-                              (int)(vec[1] * (size[1] - 1))) * res +
-                               (int)(vec[2] * (size[2] - 1)))->flag &= ~PARS_UNEXIST;
+                       pa[((int)(vec[0] * (size[0] - 1))  * res +
+                           (int)(vec[1] * (size[1] - 1))) * res +
+                           (int)(vec[2] * (size[2] - 1))].flag &= ~PARS_UNEXIST;
                }
        }
        else if (ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
@@ -556,8 +566,7 @@ static void distribute_grid(DerivedMesh *dm, ParticleSystem *psys)
                                                        else /* store number of intersections */
                                                                (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
                                                }
-                                               
-                                               if (mface->v4) {
+                                               else if (mface->v4) {
                                                        copy_v3_v3(v4, mvert[mface->v4].co);
 
                                                        if (isect_axial_line_tri_v3(a, co1, co2, v4, v1, v3, &lambda)) {
@@ -632,10 +641,10 @@ static void hammersley_create(float *out, int n, int seed, float amount)
        double p, t, offs[2];
        int k, kk;
 
-       rng = rng_new(31415926 + n + seed);
-       offs[0]= rng_getDouble(rng) + (double)amount;
-       offs[1]= rng_getDouble(rng) + (double)amount;
-       rng_free(rng);
+       rng = BLI_rng_new(31415926 + n + seed);
+       offs[0] = BLI_rng_get_double(rng) + (double)amount;
+       offs[1] = BLI_rng_get_double(rng) + (double)amount;
+       BLI_rng_free(rng);
 
        for (k = 0; k < n; k++) {
                t = 0;
@@ -643,8 +652,8 @@ static void hammersley_create(float *out, int n, int seed, float amount)
                        if (kk & 1) /* kk mod 2 = 1 */
                                t += p;
 
-               out[2*k + 0]= fmod((double)k/(double)n + offs[0], 1.0);
-               out[2*k + 1]= fmod(t + offs[1], 1.0);
+               out[2*k + 0] = fmod((double)k/(double)n + offs[0], 1.0);
+               out[2*k + 1] = fmod(t + offs[1], 1.0);
        }
 }
 
@@ -661,13 +670,13 @@ static void init_mv_jit(float *jit, int num, int seed2, float amount)
        rad2= (float)(1.0f/((float)num));
        rad3= (float)sqrt((float)num)/((float)num);
 
-       rng = rng_new(31415926 + num + seed2);
+       rng = BLI_rng_new(31415926 + num + seed2);
        x= 0;
                num2 = 2 * num;
        for (i=0; i<num2; i+=2) {
        
-               jit[i]= x + amount*rad1*(0.5f - rng_getFloat(rng));
-               jit[i+1]= i/(2.0f*num) + amount*rad1*(0.5f - rng_getFloat(rng));
+               jit[i] = x + amount*rad1*(0.5f - BLI_rng_get_float(rng));
+               jit[i+1] = i/(2.0f*num) + amount*rad1*(0.5f - BLI_rng_get_float(rng));
                
                jit[i]-= (float)floor(jit[i]);
                jit[i+1]-= (float)floor(jit[i+1]);
@@ -684,7 +693,7 @@ static void init_mv_jit(float *jit, int num, int seed2, float amount)
                BLI_jitterate2(jit, jit2, num, rad2);
        }
        MEM_freeN(jit2);
-       rng_free(rng);
+       BLI_rng_free(rng);
 }
 
 static void psys_uv_to_w(float u, float v, int quad, float *w)
@@ -698,21 +707,21 @@ static void psys_uv_to_w(float u, float v, int quad, float *w)
                        u= 1.0f-u;
        }
 
-       vert[0][0]= 0.0f; vert[0][1]= 0.0f; vert[0][2]= 0.0f;
-       vert[1][0]= 1.0f; vert[1][1]= 0.0f; vert[1][2]= 0.0f;
-       vert[2][0]= 1.0f; vert[2][1]= 1.0f; vert[2][2]= 0.0f;
+       vert[0][0] = 0.0f; vert[0][1] = 0.0f; vert[0][2] = 0.0f;
+       vert[1][0] = 1.0f; vert[1][1] = 0.0f; vert[1][2] = 0.0f;
+       vert[2][0] = 1.0f; vert[2][1] = 1.0f; vert[2][2] = 0.0f;
 
-       co[0]= u;
-       co[1]= v;
-       co[2]= 0.0f;
+       co[0] = u;
+       co[1] = v;
+       co[2] = 0.0f;
 
        if (quad) {
-               vert[3][0]= 0.0f; vert[3][1]= 1.0f; vert[3][2]= 0.0f;
+               vert[3][0] = 0.0f; vert[3][1] = 1.0f; vert[3][2] = 0.0f;
                interp_weights_poly_v3( w,vert, 4, co);
        }
        else {
                interp_weights_poly_v3( w,vert, 3, co);
-               w[3]= 0.0f;
+               w[3] = 0.0f;
        }
 }
 
@@ -797,13 +806,15 @@ static void distribute_threads_exec(ParticleThread *thread, ParticleData *pa, Ch
                        }
                        else {
                                ctx->jitoff[i] = fmod(ctx->jitoff[i],(float)ctx->jitlevel);
-                               psys_uv_to_w(ctx->jit[2*(int)ctx->jitoff[i]], ctx->jit[2*(int)ctx->jitoff[i]+1], mface->v4, pa->fuv);
-                               ctx->jitoff[i]++;
+                               if (!isnan(ctx->jitoff[i])) {
+                                       psys_uv_to_w(ctx->jit[2*(int)ctx->jitoff[i]], ctx->jit[2*(int)ctx->jitoff[i]+1], mface->v4, pa->fuv);
+                                       ctx->jitoff[i]++;
+                               }
                        }
                        break;
                case PART_DISTR_RAND:
-                       randu= rng_getFloat(thread->rng);
-                       randv= rng_getFloat(thread->rng);
+                       randu= BLI_rng_get_float(thread->rng);
+                       randv= BLI_rng_get_float(thread->rng);
                        rng_skip_tot -= 2;
 
                        psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
@@ -879,8 +890,8 @@ static void distribute_threads_exec(ParticleThread *thread, ParticleData *pa, Ch
 
                mf= dm->getTessFaceData(dm, ctx->index[p], CD_MFACE);
 
-               randu= rng_getFloat(thread->rng);
-               randv= rng_getFloat(thread->rng);
+               randu= BLI_rng_get_float(thread->rng);
+               randv= BLI_rng_get_float(thread->rng);
                rng_skip_tot -= 2;
 
                psys_uv_to_w(randu, randv, mf->v4, cpa->fuv);
@@ -932,7 +943,7 @@ static void distribute_threads_exec(ParticleThread *thread, ParticleData *pa, Ch
        }
 
        if (rng_skip_tot > 0) /* should never be below zero */
-               rng_skip(thread->rng, rng_skip_tot);
+               BLI_rng_skip(thread->rng, rng_skip_tot);
 }
 
 static void *distribute_threads_exec_cb(void *data)
@@ -949,12 +960,12 @@ static void *distribute_threads_exec_cb(void *data)
 
                for (p=0; p<totpart; p++, cpa++) {
                        if (thread->ctx->skip) /* simplification skip */
-                               rng_skip(thread->rng, PSYS_RND_DIST_SKIP * thread->ctx->skip[p]);
+                               BLI_rng_skip(thread->rng, PSYS_RND_DIST_SKIP * thread->ctx->skip[p]);
 
                        if ((p+thread->num) % thread->tot == 0)
                                distribute_threads_exec(thread, NULL, cpa, p);
                        else /* thread skip */
-                               rng_skip(thread->rng, PSYS_RND_DIST_SKIP);
+                               BLI_rng_skip(thread->rng, PSYS_RND_DIST_SKIP);
                }
        }
        else {
@@ -971,8 +982,8 @@ static void *distribute_threads_exec_cb(void *data)
 static int *COMPARE_ORIG_INDEX = NULL;
 static int distribute_compare_orig_index(const void *p1, const void *p2)
 {
-       int index1 = COMPARE_ORIG_INDEX[*(const int*)p1];
-       int index2 = COMPARE_ORIG_INDEX[*(const int*)p2];
+       int index1 = COMPARE_ORIG_INDEX[*(const int *)p1];
+       int index2 = COMPARE_ORIG_INDEX[*(const int *)p2];
 
        if (index1 < index2)
                return -1;
@@ -998,7 +1009,7 @@ static void distribute_invalid(Scene *scene, ParticleSystem *psys, int from)
 
                if (psys->child && totchild) {
                        for (p=0,cpa=psys->child; p<totchild; p++,cpa++) {
-                               cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]= 0.0;
+                               cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3] = 0.0;
                                cpa->foffset= 0.0f;
                                cpa->parent=0;
                                cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
@@ -1009,7 +1020,7 @@ static void distribute_invalid(Scene *scene, ParticleSystem *psys, int from)
        else {
                PARTICLE_P;
                LOOP_PARTICLES {
-                       pa->fuv[0]=pa->fuv[1]=pa->fuv[2]= pa->fuv[3]= 0.0;
+                       pa->fuv[0] = pa->fuv[1] = pa->fuv[2] = pa->fuv[3] = 0.0;
                        pa->foffset= 0.0f;
                        pa->num= -1;
                }
@@ -1111,7 +1122,7 @@ static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, D
 
                if (from == PART_FROM_VERT) {
                        MVert *mv= dm->getVertDataArray(dm, CD_MVERT);
-                       float (*orcodata)[3]= dm->getVertDataArray(dm, CD_ORCO);
+                       float (*orcodata)[3] = dm->getVertDataArray(dm, CD_ORCO);
                        int totvert = dm->getNumVerts(dm);
 
                        tree=BLI_kdtree_new(totvert);
@@ -1243,9 +1254,9 @@ static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, D
        inv_totweight = (totweight > 0.f ? 1.f/totweight : 0.f);
 
        /* Calculate cumulative weights */
-       element_sum[0]= 0.0f;
+       element_sum[0] = 0.0f;
        for (i=0; i<totelem; i++)
-               element_sum[i+1]= element_sum[i] + element_weight[i] * inv_totweight;
+               element_sum[i+1] = element_sum[i] + element_weight[i] * inv_totweight;
        
        /* Finally assign elements to particles */
        if ((part->flag&PART_TRAND) || (part->simplify_flag&PART_SIMPLIFY_ENABLE)) {
@@ -1254,9 +1265,9 @@ static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, D
                for (p=0; p<totpart; p++) {
                        /* In theory element_sum[totelem] should be 1.0, but due to float errors this is not necessarily always true, so scale pos accordingly. */
                        pos= BLI_frand() * element_sum[totelem];
-                       particle_element[p]= distribute_binary_search(element_sum, totelem, pos);
-                       particle_element[p]= MIN2(totelem-1, particle_element[p]);
-                       jitter_offset[particle_element[p]]= pos;
+                       particle_element[p] = distribute_binary_search(element_sum, totelem, pos);
+                       particle_element[p] = MIN2(totelem-1, particle_element[p]);
+                       jitter_offset[particle_element[p]] = pos;
                }
        }
        else {
@@ -1267,16 +1278,16 @@ static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, D
                i= 0;
 
                for (p=0; p<totpart; p++, pos+=step) {
-                       while ((i < totelem) && (pos > element_sum[i+1]))
+                       while ((i < totelem) && (pos > (double)element_sum[i + 1]))
                                i++;
 
-                       particle_element[p]= MIN2(totelem-1, i);
+                       particle_element[p] = MIN2(totelem-1, i);
 
                        /* avoid zero weight face */
                        if (p == totpart-1 && element_weight[particle_element[p]] == 0.0f)
-                               particle_element[p]= particle_element[p-1];
+                               particle_element[p] = particle_element[p-1];
 
-                       jitter_offset[particle_element[p]]= pos;
+                       jitter_offset[particle_element[p]] = pos;
                }
        }
 
@@ -1351,7 +1362,7 @@ static int distribute_threads_init_data(ParticleThread *threads, Scene *scene, D
        
        seed= 31415926 + ctx->sim.psys->seed;
        for (i=0; i<totthread; i++) {
-               threads[i].rng= rng_new(seed);
+               threads[i].rng= BLI_rng_new(seed);
                threads[i].tot= totthread;
        }
 
@@ -1490,9 +1501,9 @@ void psys_threads_free(ParticleThread *threads)
        /* threads */
        for (i=0; i<totthread; i++) {
                if (threads[i].rng)
-                       rng_free(threads[i].rng);
+                       BLI_rng_free(threads[i].rng);
                if (threads[i].rng_path)
-                       rng_free(threads[i].rng_path);
+                       BLI_rng_free(threads[i].rng_path);
        }
 
        MEM_freeN(ctx);
@@ -1576,12 +1587,12 @@ static void initialize_all_particles(ParticleSimulationData *sim)
        }
 }
 
-static void get_angular_velocity_vector(short avemode, ParticleKey *state, float *vec)
+static void get_angular_velocity_vector(short avemode, ParticleKey *state, float vec[3])
 {
        switch (avemode) {
                case PART_AVE_VELOCITY:
                        copy_v3_v3(vec, state->vel);
-                       break;  
+                       break;
                case PART_AVE_HORIZONTAL:
                {
                        float zvec[3];
@@ -1620,9 +1631,9 @@ void psys_get_birth_coordinates(ParticleSimulationData *sim, ParticleData *pa, P
        ParticleSystem *psys = sim->psys;
        ParticleSettings *part;
        ParticleTexture ptex;
-       float fac, phasefac, nor[3]={0,0,0},loc[3],vel[3]={0.0,0.0,0.0},rot[4],q2[4];
-       float r_vel[3],r_ave[3],r_rot[4],vec[3],p_vel[3]={0.0,0.0,0.0};
-       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}, rot_vec[3]={0.0,0.0,0.0};
+       float fac, phasefac, nor[3] = {0,0,0},loc[3],vel[3] = {0.0,0.0,0.0},rot[4],q2[4];
+       float r_vel[3],r_ave[3],r_rot[4],vec[3],p_vel[3] = {0.0,0.0,0.0};
+       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}, rot_vec[3] = {0.0,0.0,0.0};
        float q_phase[4];
        int p = pa - psys->particles;
        part=psys->part;
@@ -1879,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;
        }
@@ -2130,16 +2140,22 @@ static void psys_update_effectors(ParticleSimulationData *sim)
        precalc_guides(sim, sim->psys->effectors);
 }
 
-static void integrate_particle(ParticleSettings *part, ParticleData *pa, float dtime, float *external_acceleration, void (*force_func)(void *forcedata, ParticleKey *state, float *force, float *impulse), void *forcedata)
+static void integrate_particle(ParticleSettings *part, ParticleData *pa, float dtime, float *external_acceleration,
+                               void (*force_func)(void *forcedata, ParticleKey *state, float *force, float *impulse),
+                               void *forcedata)
 {
+#define ZERO_F43 {{0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}, {0.0f, 0.0f, 0.0f}}
+
        ParticleKey states[5];
-       float force[3],acceleration[3],impulse[3],dx[4][3],dv[4][3],oldpos[3];
+       float force[3], acceleration[3], impulse[3], dx[4][3] = ZERO_F43, dv[4][3] = ZERO_F43, oldpos[3];
        float pa_mass= (part->flag & PART_SIZEMASS ? part->mass * pa->size : part->mass);
        int i, steps=1;
        int integrator = part->integrator;
 
+#undef ZERO_F43
+
        copy_v3_v3(oldpos, pa->state.co);
-       
+
        /* Verlet integration behaves strangely with moving emitters, so do first step with euler. */
        if (pa->prev_state.time < 0.f && integrator == PART_INT_VERLET)
                integrator = PART_INT_EULER;
@@ -2159,7 +2175,9 @@ static void integrate_particle(ParticleSettings *part, ParticleData *pa, float d
                        break;
        }
 
-       copy_particle_key(states, &pa->state, 1);
+       for (i=0; i<steps; i++) {
+               copy_particle_key(states + i, &pa->state, 1);
+       }
 
        states->time = 0.f;
 
@@ -2363,43 +2381,46 @@ static EdgeHash *sph_springhash_build(ParticleSystem *psys)
 }
 
 #define SPH_NEIGHBORS 512
-typedef struct SPHNeighbor
-{
+typedef struct SPHNeighbor {
        ParticleSystem *psys;
        int index;
 } SPHNeighbor;
-typedef struct SPHRangeData
-{
+
+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 mass;
        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->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;
@@ -2429,8 +2450,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;
 }
 
 /*
@@ -2442,7 +2463,8 @@ static void sph_particle_courant(SPHData *sphdata, SPHRangeData *pfr)
        int i;
        float flow[3], offset[3], dist;
 
-       flow[0] = flow[1] = flow[2] = 0.0f;
+       zero_v3(flow);
+
        dist = 0.0f;
        if (pfr->tot_neighbors > 0) {
                pa = pfr->pa;
@@ -2470,7 +2492,6 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa
        ParticleSpring *spring = NULL;
        SPHRangeData pfr;
        SPHNeighbor *pfn;
-       float mass = sphdata->mass;
        float *gravity = sphdata->gravity;
        EdgeHash *springhash = sphdata->eh;
 
@@ -2480,10 +2501,12 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa
        float visc = fluid->viscosity_omega;
        float stiff_visc = fluid->viscosity_beta * (fluid->flag & SPH_FAC_VISCOSITY ? fluid->viscosity_omega : 1.f);
 
-       float inv_mass = 1.0f/mass;
+       float inv_mass = 1.0f / sphdata->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);
 
@@ -2494,24 +2517,24 @@ 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;
+       pfr.mass = sphdata->mass;
 
-       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++) {
@@ -2533,7 +2556,7 @@ static void sph_force_cb(void *sphdata_v, ParticleKey *state, float *force, floa
                madd_v3_v3fl(force, vec, -(pressure + near_pressure*q)*q);
 
                /* Viscosity */
-               if (visc > 0.f  || stiff_visc > 0.f) {          
+               if (visc > 0.f  || stiff_visc > 0.f) {
                        sub_v3_v3v3(dv, vel, state->vel);
                        u = dot_v3v3(vec, dv);
 
@@ -2575,14 +2598,209 @@ 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++;
+}
+
+/* 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->npsys->part->mass;
+
+       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);
+
+       // Use speed of sound squared
+       float stiffness = pow2(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 sph_solver_init(ParticleSimulationData *sim, SPHData *sphdata)
+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;
+       pfr.mass = sphdata->mass;
+
+       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;
@@ -2603,17 +2821,47 @@ 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;
+       pfr.mass = sphdata->mass;
+
+       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)
 {
@@ -2641,8 +2889,7 @@ static void sph_integrate(ParticleSimulationData *sim, ParticleData *pa, float d
 /************************************************/
 /*                     Basic physics                                           */
 /************************************************/
-typedef struct EfData
-{
+typedef struct EfData {
        ParticleTexture ptex;
        ParticleSimulationData *sim;
        ParticleData *pa;
@@ -2673,6 +2920,9 @@ static void basic_force_cb(void *efdata_v, ParticleKey *state, float *force, flo
                force[1] += (BLI_frand()-0.5f) * part->brownfac;
                force[2] += (BLI_frand()-0.5f) * part->brownfac;
        }
+
+       if (part->flag & PART_ROT_DYN && epoint.ave)
+               copy_v3_v3(pa->state.ave, epoint.ave);
 }
 /* gathers all forces that effect particles and calculates a new state for the particle */
 static void basic_integrate(ParticleSimulationData *sim, int p, float dfra, float cfra)
@@ -2730,22 +2980,29 @@ static void basic_integrate(ParticleSimulationData *sim, int p, float dfra, floa
 }
 static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
 {
-       float rotfac, rot1[4], rot2[4]={1.0,0.0,0.0,0.0}, dtime=dfra*timestep;
+       float rotfac, rot1[4], rot2[4] = {1.0,0.0,0.0,0.0}, dtime=dfra*timestep, extrotfac;
 
-       if ((part->flag & PART_ROTATIONS)==0) {
-               pa->state.rot[0]=1.0f;
-               pa->state.rot[1]=pa->state.rot[2]=pa->state.rot[3]=0;
+       if ((part->flag & PART_ROTATIONS) == 0) {
+               unit_qt(pa->state.rot);
                return;
        }
 
-       if ((part->flag & PART_ROT_DYN)==0 && ELEM3(part->avemode, PART_AVE_VELOCITY, PART_AVE_HORIZONTAL, PART_AVE_VERTICAL)) {
+       if (part->flag & PART_ROT_DYN) {
+               extrotfac = len_v3(pa->state.ave);
+       }
+       else {
+               extrotfac = 0.0f;
+       }
+
+       if ((part->flag & PART_ROT_DYN) && ELEM3(part->avemode, PART_AVE_VELOCITY, PART_AVE_HORIZONTAL, PART_AVE_VERTICAL)) {
                float angle;
                float len1 = len_v3(pa->prev_state.vel);
                float len2 = len_v3(pa->state.vel);
                float vec[3];
 
-               if (len1==0.0f || len2==0.0f)
-                       pa->state.ave[0] = pa->state.ave[1] = pa->state.ave[2] = 0.0f;
+               if (len1 == 0.0f || len2 == 0.0f) {
+                       zero_v3(pa->state.ave);
+               }
                else {
                        cross_v3_v3v3(pa->state.ave, pa->prev_state.vel, pa->state.vel);
                        normalize_v3(pa->state.ave);
@@ -2758,9 +3015,8 @@ static void basic_rotate(ParticleSettings *part, ParticleData *pa, float dfra, f
        }
 
        rotfac = len_v3(pa->state.ave);
-       if (rotfac == 0.0f) { /* unit_qt(in VecRotToQuat) doesn't give unit quat [1,0,0,0]?? */
-               rot1[0]=1.0f;
-               rot1[1]=rot1[2]=rot1[3]=0;
+       if (rotfac == 0.0f || (part->flag & PART_ROT_DYN)==0 || extrotfac == 0.0f) {
+               unit_qt(rot1);
        }
        else {
                axis_angle_to_quat(rot1,pa->state.ave,rotfac*dtime);
@@ -3154,8 +3410,7 @@ void BKE_psys_collision_neartest_cb(void *userdata, int index, const BVHTreeRay
        if (col->hit == col->current && col->pce.index == index && col->pce.tot == 3)
                return;
 
-       do
-       {       
+       do {
                collision = collision_sphere_to_tri(col, ray->radius, &pce, &t);
                if (col->pce.inside == 0) {
                        collision += collision_sphere_to_edges(col, ray->radius, &pce, &t);
@@ -3317,7 +3572,7 @@ static int collision_response(ParticleData *pa, ParticleCollision *col, BVHTreeR
                        }
                }
 
-               /* stickness was possibly added before, so cancel that before calculating new normal velocity */
+               /* stickiness was possibly added before, so cancel that before calculating new normal velocity */
                /* otherwise particles go flying out of the surface because of high reversed sticky velocity */
                if (v0_dot < 0.0f) {
                        v0_dot += pd->pdef_stickness;
@@ -3376,7 +3631,7 @@ static int collision_response(ParticleData *pa, ParticleCollision *col, BVHTreeR
                                madd_v3_v3fl(pa->state.vel, nor, -dot);
                }
 
-               /* add stickness to surface */
+               /* add stickiness to surface */
                madd_v3_v3fl(pa->state.vel, pce->nor, -pd->pdef_stickness);
 
                /* set coordinates for next iteration */
@@ -3747,25 +4002,28 @@ static void save_hair(ParticleSimulationData *sim, float UNUSED(cfra))
                pa->totkey++;
 
                /* root is always in the origin of hair space so we set it to be so after the last key is saved*/
-               if (pa->totkey == psys->part->hair_step + 1)
-                       root->co[0] = root->co[1] = root->co[2] = 0.0f;
+               if (pa->totkey == psys->part->hair_step + 1) {
+                       zero_v3(root->co);
+               }
+
        }
 }
 
 /* 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 neighbourning particles. */
-void update_courant_num(ParticleSimulationData *sim, ParticleData *pa,
-       float dtime, SPHData *sphdata)
+ * the simulation, and is typically the distance to neighboring particles. */
+static void update_courant_num(ParticleSimulationData *sim, ParticleData *pa,
+                               float dtime, SPHData *sphdata)
 {
        float relative_vel[3];
        float speed;
@@ -3775,15 +4033,31 @@ 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. */
-float update_timestep(ParticleSystem *psys, ParticleSimulationData *sim,
-       float t_frac)
+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
-               psys->dt_frac *= (psys->part->courant_target / sim->courant_num);
-       CLAMP(psys->dt_frac, MIN_TIMESTEP, 1.0f);
+               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 = dt_target;
 
        /* Sync with frame end if it's close. */
        if (t_frac == 1.0f)
@@ -3947,31 +4221,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->solver == 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;
                }
        }
@@ -4063,7 +4378,7 @@ static void particles_fluid_step(ParticleSimulationData *sim, int UNUSED(cfra))
        {
                FluidsimModifierData *fluidmd = (FluidsimModifierData *)modifiers_findByType(sim->ob, eModifierType_Fluidsim);
                
-               if ( fluidmd && fluidmd->fss) { 
+               if ( fluidmd && fluidmd->fss) {
                        FluidsimSettings *fss= fluidmd->fss;
                        ParticleSettings *part = psys->part;
                        ParticleData *pa=NULL;
@@ -4092,7 +4407,7 @@ static void particles_fluid_step(ParticleSimulationData *sim, int UNUSED(cfra))
                        }
        
                        gzread(gzf, &totpart, sizeof(totpart));
-                       totpart = (G.rendering)?totpart:(part->disp*totpart)/100;
+                       totpart = (G.is_rendering)?totpart:(part->disp*totpart) / 100;
                        
                        part->totpart= totpart;
                        part->sta=part->end = 1.0f;
@@ -4108,10 +4423,10 @@ static void particles_fluid_step(ParticleSimulationData *sim, int UNUSED(cfra))
                                int ptype=0;
        
                                gzread(gzf, &ptype, sizeof( ptype )); 
-                               if (ptype&readMask) {
+                               if (ptype & readMask) {
                                        activeParts++;
        
-                                       gzread(gzf, &(pa->size), sizeof( float )); 
+                                       gzread(gzf, &(pa->size), sizeof(float));
        
                                        pa->size /= 10.0f;
        
@@ -4127,15 +4442,14 @@ static void particles_fluid_step(ParticleSimulationData *sim, int UNUSED(cfra))
                                                pa->state.vel[j] = wrf;
                                        }
        
-                                       pa->state.ave[0] = pa->state.ave[1] = pa->state.ave[2] = 0.0f;
-                                       pa->state.rot[0] = 1.0;
-                                       pa->state.rot[1] = pa->state.rot[2] = pa->state.rot[3] = 0.0;
+                                       zero_v3(pa->state.ave);
+                                       unit_qt(pa->state.rot);
        
                                        pa->time = 1.f;
                                        pa->dietime = sim->scene->r.efra + 1;
                                        pa->lifetime = sim->scene->r.efra;
                                        pa->alive = PARS_ALIVE;
-                                       //if (a < 25) fprintf(stderr,"FSPARTICLE debug set %s , a%d = %f,%f,%f , life=%f\n", filename, a, pa->co[0],pa->co[1],pa->co[2], pa->lifetime );
+                                       //if (a < 25) fprintf(stderr,"FSPARTICLE debug set %s, a%d = %f,%f,%f, life=%f\n", filename, a, pa->co[0],pa->co[1],pa->co[2], pa->lifetime );
                                }
                                else {
                                        // skip...
@@ -4186,7 +4500,7 @@ static void system_step(ParticleSimulationData *sim, float cfra)
        int startframe = 0, endframe = 100, oldtotpart = 0;
 
        /* cache shouldn't be used for hair or "continue physics" */
-       if (part->type != PART_HAIR && BKE_ptcache_get_continue_physics() == 0) {
+       if (part->type != PART_HAIR) {
                psys_clear_temp_pointcache(psys);
 
                /* set suitable cache range automatically */
@@ -4284,14 +4598,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;
                }
 
@@ -4444,7 +4758,8 @@ static void psys_prepare_physics(ParticleSimulationData *sim)
 static int hair_needs_recalc(ParticleSystem *psys)
 {
        if (!(psys->flag & PSYS_EDITED) && (!psys->edit || !psys->edit->edited) &&
-               ((psys->flag & PSYS_HAIR_DONE)==0 || psys->recalc & PSYS_RECALC_RESET || (psys->part->flag & PART_HAIR_REGROW && !psys->edit))) {
+           ((psys->flag & PSYS_HAIR_DONE)==0 || psys->recalc & PSYS_RECALC_RESET || (psys->part->flag & PART_HAIR_REGROW && !psys->edit)))
+       {
                return 1;
        }