diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp index 153993dc9b8..821b08ce334 100644 --- a/source/blender/physics/intern/BPH_mass_spring.cpp +++ b/source/blender/physics/intern/BPH_mass_spring.cpp @@ -453,54 +453,6 @@ static void hair_get_boundbox(ClothModifierData *clmd, float gmin[3], float gmax } } -static void cloth_calc_volume_force(ClothModifierData *clmd) -{ - ClothSimSettings *parms = clmd->sim_parms; - Cloth *cloth = clmd->clothObject; - Implicit_Data *data = cloth->implicit; - int numverts = cloth->numverts; - - /* 2.0f is an experimental value that seems to give good results */ - float smoothfac = 2.0f * parms->velocity_smooth; - // float collfac = 2.0f * parms->collider_friction; - float pressfac = parms->pressure; - float minpress = parms->pressure_threshold; - float gmin[3], gmax[3]; - int i; - - hair_get_boundbox(clmd, gmin, gmax); - - /* gather velocities & density */ - if (smoothfac > 0.0f || pressfac > 0.0f) { - HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax); - - for (i = 0; i < numverts; i++) { - float x[3], v[3]; - - BPH_mass_spring_get_motion_state(data, i, x, v); - BPH_hair_volume_add_vertex(vertex_grid, x, v); - } - BPH_hair_volume_normalize_vertex_grid(vertex_grid); - -#if 0 - /* apply velocity filter */ - BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size); -#endif - - for (i = 0; i < numverts; i++) { - float x[3], v[3], f[3], dfdx[3][3], dfdv[3][3]; - - /* calculate volumetric forces */ - BPH_mass_spring_get_motion_state(data, i, x, v); - BPH_hair_volume_vertex_grid_forces(vertex_grid, x, v, smoothfac, pressfac, minpress, f, dfdx, dfdv); - /* apply on hair data */ - BPH_mass_spring_force_extern(data, i, f, dfdx, dfdv); - } - - BPH_hair_volume_free_vertex_grid(vertex_grid); - } -} - static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListBase *effectors, float time) { /* Collect forces and derivatives: F, dFdX, dFdV */ @@ -525,7 +477,7 @@ static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListB } #endif - cloth_calc_volume_force(clmd); + /* cloth_calc_volume_force(clmd); */ #ifdef CLOTH_FORCE_DRAG BPH_mass_spring_force_drag(data, drag); @@ -570,6 +522,125 @@ static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListB } } +#if 0 +static void cloth_calc_volume_force(ClothModifierData *clmd) +{ + ClothSimSettings *parms = clmd->sim_parms; + Cloth *cloth = clmd->clothObject; + Implicit_Data *data = cloth->implicit; + int numverts = cloth->numverts; + ClothVertex *vert; + + /* 2.0f is an experimental value that seems to give good results */ + float smoothfac = 2.0f * parms->velocity_smooth; + float collfac = 2.0f * parms->collider_friction; + float pressfac = parms->pressure; + float minpress = parms->pressure_threshold; + float gmin[3], gmax[3]; + int i; + + hair_get_boundbox(clmd, gmin, gmax); + + /* gather velocities & density */ + if (smoothfac > 0.0f || pressfac > 0.0f) { + HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax); + + vert = cloth->verts; + for (i = 0; i < numverts; i++, vert++) { + float x[3], v[3]; + + if (vert->solver_index < 0) { + copy_v3_v3(x, vert->x); + copy_v3_v3(v, vert->v); + } + else { + BPH_mass_spring_get_motion_state(data, vert->solver_index, x, v); + } + BPH_hair_volume_add_vertex(vertex_grid, x, v); + } + BPH_hair_volume_normalize_vertex_grid(vertex_grid); + +#if 0 + /* apply velocity filter */ + BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size); +#endif + + vert = cloth->verts; + for (i = 0; i < numverts; i++, vert++) { + float x[3], v[3], f[3], dfdx[3][3], dfdv[3][3]; + + if (vert->solver_index < 0) + continue; + + /* calculate volumetric forces */ + BPH_mass_spring_get_motion_state(data, vert->solver_index, x, v); + BPH_hair_volume_vertex_grid_forces(vertex_grid, x, v, smoothfac, pressfac, minpress, f, dfdx, dfdv); + /* apply on hair data */ + BPH_mass_spring_force_extern(data, vert->solver_index, f, dfdx, dfdv); + } + + BPH_hair_volume_free_vertex_grid(vertex_grid); + } +} +#endif + +static void cloth_continuum_step(ClothModifierData *clmd) +{ + ClothSimSettings *parms = clmd->sim_parms; + Cloth *cloth = clmd->clothObject; + Implicit_Data *data = cloth->implicit; + int numverts = cloth->numverts; + ClothVertex *vert; + + const float fluid_factor = 0.95f; /* blend between PIC and FLIP methods */ + /* 2.0f is an experimental value that seems to give good results */ + float smoothfac = 2.0f * parms->velocity_smooth; + float collfac = 2.0f * parms->collider_friction; + float pressfac = parms->pressure; + float minpress = parms->pressure_threshold; + float gmin[3], gmax[3]; + int i; + + hair_get_boundbox(clmd, gmin, gmax); + + /* gather velocities & density */ + if (smoothfac > 0.0f || pressfac > 0.0f) { + HairVertexGrid *vertex_grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_res, gmin, gmax); + + vert = cloth->verts; + for (i = 0; i < numverts; i++, vert++) { + float x[3], v[3]; + + BPH_mass_spring_get_motion_state(data, i, x, v); + BPH_hair_volume_add_vertex(vertex_grid, x, v); + } + BPH_hair_volume_normalize_vertex_grid(vertex_grid); + +#if 0 + /* apply velocity filter */ + BPH_hair_volume_vertex_grid_filter_box(vertex_grid, clmd->sim_parms->voxel_filter_size); +#endif + + vert = cloth->verts; + for (i = 0; i < numverts; i++, vert++) { + float x[3], v[3], nv[3]; + + /* calculate volumetric velocity influence */ + BPH_mass_spring_get_position(data, i, x); + BPH_mass_spring_get_new_velocity(data, i, v); + + BPH_hair_volume_grid_velocity(vertex_grid, x, v, fluid_factor, nv); + + interp_v3_v3v3(nv, v, nv, smoothfac); + + /* apply on hair data */ + BPH_mass_spring_set_new_velocity(data, i, nv); + } + + BPH_hair_volume_free_vertex_grid(vertex_grid); + } +} + static void cloth_clear_result(ClothModifierData *clmd) { ClothSolverResult *sres = clmd->solver_result; @@ -686,9 +757,13 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase * cloth_calc_force(clmd, frame, effectors, step); // calculate new velocity and position - BPH_mass_spring_solve(id, dt, &result); + BPH_mass_spring_solve_velocities(id, dt, &result); cloth_record_result(clmd, &result, clmd->sim_parms->stepsPerFrame); + cloth_continuum_step(clmd); + + BPH_mass_spring_solve_positions(id, dt); + BPH_mass_spring_apply_result(id); /* move pinned verts to correct position */ diff --git a/source/blender/physics/intern/hair_volume.c b/source/blender/physics/intern/hair_volume.c index e72df493b79..87f202bbb2c 100644 --- a/source/blender/physics/intern/hair_volume.c +++ b/source/blender/physics/intern/hair_volume.c @@ -220,6 +220,18 @@ void BPH_hair_volume_vertex_grid_forces(HairVertexGrid *grid, const float x[3], mul_m3_fl(dfdv, smoothfac); } +void BPH_hair_volume_grid_velocity(HairVertexGrid *grid, const float x[3], const float v[3], + float fluid_factor, + float r_v[3]) +{ + float gdensity, gvelocity[3], ggrad[3], gvelgrad[3][3]; + + hair_grid_interpolate(grid->verts, grid->res, grid->gmin, grid->scale, x, &gdensity, gvelocity, ggrad, gvelgrad); + + /* XXX TODO implement FLIP method and use fluid_factor to blend between FLIP and PIC */ + copy_v3_v3(r_v, gvelocity); +} + BLI_INLINE bool hair_grid_point_valid(const float vec[3], float gmin[3], float gmax[3]) { return !(vec[0] < gmin[0] || vec[1] < gmin[1] || vec[2] < gmin[2] || diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h index 4678b273dcf..012125b4c87 100644 --- a/source/blender/physics/intern/implicit.h +++ b/source/blender/physics/intern/implicit.h @@ -111,20 +111,26 @@ BLI_INLINE int hash_collpair(int type, CollPair *collpair) void BPH_mass_spring_solver_debug_data(struct Implicit_Data *id, struct SimDebugData *debug_data); +void BPH_mass_spring_set_vertex_mass(struct Implicit_Data *data, int index, float mass); void BPH_mass_spring_set_rest_transform(struct Implicit_Data *data, int index, float rot[3][3]); void BPH_mass_spring_set_motion_state(struct Implicit_Data *data, int index, const float x[3], const float v[3]); void BPH_mass_spring_set_position(struct Implicit_Data *data, int index, const float x[3]); void BPH_mass_spring_set_velocity(struct Implicit_Data *data, int index, const float v[3]); void BPH_mass_spring_get_motion_state(struct Implicit_Data *data, int index, float x[3], float v[3]); -void BPH_mass_spring_set_vertex_mass(struct Implicit_Data *data, int index, float mass); +void BPH_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3]); + +/* modified velocity during solver step */ +void BPH_mass_spring_get_new_velocity(struct Implicit_Data *data, int index, float v[3]); +void BPH_mass_spring_set_new_velocity(struct Implicit_Data *data, int index, const float v[3]); void BPH_mass_spring_clear_constraints(struct Implicit_Data *data); void BPH_mass_spring_add_constraint_ndof0(struct Implicit_Data *data, int index, const float dV[3]); void BPH_mass_spring_add_constraint_ndof1(struct Implicit_Data *data, int index, const float c1[3], const float c2[3], const float dV[3]); void BPH_mass_spring_add_constraint_ndof2(struct Implicit_Data *data, int index, const float c1[3], const float dV[3]); -bool BPH_mass_spring_solve(struct Implicit_Data *data, float dt, struct ImplicitSolverResult *result); +bool BPH_mass_spring_solve_velocities(struct Implicit_Data *data, float dt, struct ImplicitSolverResult *result); +bool BPH_mass_spring_solve_positions(struct Implicit_Data *data, float dt); void BPH_mass_spring_apply_result(struct Implicit_Data *data); /* Clear the force vector at the beginning of the time step */ @@ -170,6 +176,20 @@ void BPH_hair_volume_normalize_vertex_grid(struct HairVertexGrid *grid); #if 0 /* XXX weighting is incorrect, disabled for now */ void BPH_hair_volume_vertex_grid_filter_box(struct HairVertexGrid *grid, int kernel_size); #endif + +/* Effect of fluid simulation grid on velocities. + * fluid_factor controls blending between PIC (Particle-in-Cell) + * and FLIP (Fluid-Implicit-Particle) methods (0 = only PIC, 1 = only FLIP) + */ +void BPH_hair_volume_grid_velocity(struct HairVertexGrid *grid, const float x[3], const float v[3], + float fluid_factor, + float r_v[3]); +/* XXX Warning: expressing grid effects on velocity as a force is not very stable, + * due to discontinuities in interpolated values! + * Better use hybrid approaches such as described in + * "Detail Preserving Continuum Simulation of Straight Hair" + * (McAdams, Selle 2009) + */ void BPH_hair_volume_vertex_grid_forces(struct HairVertexGrid *grid, const float x[3], const float v[3], float smoothfac, float pressurefac, float minpressure, float f[3], float dfdx[3][3], float dfdv[3][3]); diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c index 6b82dd3f736..865f100e7e6 100644 --- a/source/blender/physics/intern/implicit_blender.c +++ b/source/blender/physics/intern/implicit_blender.c @@ -1138,7 +1138,7 @@ static int cg_filtered_pre(lfVector *dv, fmatrix3x3 *lA, lfVector *lB, lfVector } #endif -bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult *result) +bool BPH_mass_spring_solve_velocities(Implicit_Data *data, float dt, ImplicitSolverResult *result) { unsigned int numverts = data->dFdV[0].vcount; @@ -1163,14 +1163,22 @@ bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult * // advance velocities add_lfvector_lfvector(data->Vnew, data->V, data->dV, numverts); - // advance positions - add_lfvector_lfvectorS(data->Xnew, data->X, data->Vnew, dt, numverts); del_lfvector(dFdXmV); return result->status == BPH_SOLVER_SUCCESS; } +bool BPH_mass_spring_solve_positions(Implicit_Data *data, float dt) +{ + int numverts = data->M[0].vcount; + + // advance positions + add_lfvector_lfvectorS(data->Xnew, data->X, data->Vnew, dt, numverts); + + return true; +} + void BPH_mass_spring_apply_result(Implicit_Data *data) { int numverts = data->M[0].vcount; @@ -1178,6 +1186,12 @@ void BPH_mass_spring_apply_result(Implicit_Data *data) cp_lfvector(data->V, data->Vnew, numverts); } +void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass) +{ + unit_m3(data->M[index].m); + mul_m3_fl(data->M[index].m, mass); +} + void BPH_mass_spring_set_rest_transform(Implicit_Data *data, int index, float tfm[3][3]) { #ifdef CLOTH_ROOT_FRAME @@ -1210,12 +1224,23 @@ void BPH_mass_spring_get_motion_state(struct Implicit_Data *data, int index, flo if (v) root_to_world_v3(data, index, v, data->V[index]); } -void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass) +void BPH_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3]) { - unit_m3(data->M[index].m); - mul_m3_fl(data->M[index].m, mass); + root_to_world_v3(data, index, x, data->X[index]); } +void BPH_mass_spring_get_new_velocity(struct Implicit_Data *data, int index, float v[3]) +{ + root_to_world_v3(data, index, v, data->Vnew[index]); +} + +void BPH_mass_spring_set_new_velocity(struct Implicit_Data *data, int index, const float v[3]) +{ + world_to_root_v3(data, index, data->Vnew[index], v); +} + +/* -------------------------------- */ + static int BPH_mass_spring_add_block(Implicit_Data *data, int v1, int v2) { int s = data->M[0].vcount + data->num_blocks; /* index from array start */ diff --git a/source/blender/physics/intern/implicit_eigen.cpp b/source/blender/physics/intern/implicit_eigen.cpp index 4543cc0a3d3..ee95e94a843 100644 --- a/source/blender/physics/intern/implicit_eigen.cpp +++ b/source/blender/physics/intern/implicit_eigen.cpp @@ -503,7 +503,7 @@ BLI_INLINE void root_to_world_m3(Implicit_Data *data, int index, float r[3][3], /* ================================ */ -bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult *result) +bool BPH_mass_spring_solve_velocities(Implicit_Data *data, float dt, ImplicitSolverResult *result) { #ifdef USE_EIGEN_CORE typedef ConjugateGradient solver_t; @@ -555,7 +555,6 @@ bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult * #endif data->Vnew = data->V + data->dV; - data->Xnew = data->X + data->Vnew * dt; switch (cg.info()) { case Eigen::Success: result->status = BPH_SOLVER_SUCCESS; break; @@ -567,7 +566,13 @@ bool BPH_mass_spring_solve(Implicit_Data *data, float dt, ImplicitSolverResult * result->iterations = cg.iterations(); result->error = cg.error(); - return cg.info() != Eigen::Success; + return cg.info() == Eigen::Success; +} + +bool BPH_mass_spring_solve_positions(Implicit_Data *data, float dt) +{ + data->Xnew = data->X + data->Vnew * dt; + return true; } /* ================================ */ @@ -578,6 +583,14 @@ void BPH_mass_spring_apply_result(Implicit_Data *data) data->V = data->Vnew; } +void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass) +{ + float m[3][3]; + copy_m3_m3(m, I); + mul_m3_fl(m, mass); + data->iM.add(index, index, m); +} + void BPH_mass_spring_set_rest_transform(Implicit_Data *data, int index, float tfm[3][3]) { #ifdef CLOTH_ROOT_FRAME @@ -610,12 +623,19 @@ void BPH_mass_spring_get_motion_state(struct Implicit_Data *data, int index, flo if (v) root_to_world_v3(data, index, v, data->V.v3(index)); } -void BPH_mass_spring_set_vertex_mass(Implicit_Data *data, int index, float mass) +void BPH_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3]) { - float m[3][3]; - copy_m3_m3(m, I); - mul_m3_fl(m, mass); - data->iM.add(index, index, m); + root_to_world_v3(data, index, x, data->X.v3(index)); +} + +void BPH_mass_spring_get_new_velocity(Implicit_Data *data, int index, float v[3]) +{ + root_to_world_v3(data, index, v, data->V.v3(index)); +} + +void BPH_mass_spring_set_new_velocity(Implicit_Data *data, int index, const float v[3]) +{ + world_to_root_v3(data, index, data->V.v3(index), v); } void BPH_mass_spring_clear_constraints(Implicit_Data *data)