diff options
Diffstat (limited to 'source/blender/physics')
| -rw-r--r-- | source/blender/physics/BPH_mass_spring.h | 5 | ||||
| -rw-r--r-- | source/blender/physics/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | source/blender/physics/intern/BPH_mass_spring.cpp | 223 | ||||
| -rw-r--r-- | source/blender/physics/intern/hair_volume.cpp | 57 | ||||
| -rw-r--r-- | source/blender/physics/intern/implicit.h | 14 | ||||
| -rw-r--r-- | source/blender/physics/intern/implicit_blender.c | 219 |
6 files changed, 279 insertions, 240 deletions
diff --git a/source/blender/physics/BPH_mass_spring.h b/source/blender/physics/BPH_mass_spring.h index c650d83c927..f1eb049dd52 100644 --- a/source/blender/physics/BPH_mass_spring.h +++ b/source/blender/physics/BPH_mass_spring.h @@ -40,6 +40,7 @@ struct Implicit_Data; struct Object; struct ClothModifierData; struct ListBase; +struct Depsgraph; struct VoxelData; typedef enum eMassSpringSolverStatus { @@ -55,11 +56,9 @@ int BPH_mass_spring_solver_numvert(struct Implicit_Data *id); int BPH_cloth_solver_init(struct Object *ob, struct ClothModifierData *clmd); void BPH_cloth_solver_free(struct ClothModifierData *clmd); -int BPH_cloth_solve(struct Object *ob, float frame, struct ClothModifierData *clmd, struct ListBase *effectors); +int BPH_cloth_solve(struct Depsgraph *depsgraph, struct Object *ob, float frame, struct ClothModifierData *clmd, struct ListBase *effectors); void BKE_cloth_solver_set_positions(struct ClothModifierData *clmd); -bool BPH_cloth_solver_get_texture_data(struct Object *ob, struct ClothModifierData *clmd, struct VoxelData *vd); - #ifdef __cplusplus } #endif diff --git a/source/blender/physics/CMakeLists.txt b/source/blender/physics/CMakeLists.txt index 0a4ff3fe0f0..b8663a384a7 100644 --- a/source/blender/physics/CMakeLists.txt +++ b/source/blender/physics/CMakeLists.txt @@ -28,6 +28,7 @@ set(INC intern ../blenlib ../blenkernel + ../depsgraph ../imbuf ../makesdna ../../../intern/guardedalloc diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp index 1b8b05ac752..111d5c5b68f 100644 --- a/source/blender/physics/intern/BPH_mass_spring.cpp +++ b/source/blender/physics/intern/BPH_mass_spring.cpp @@ -51,6 +51,9 @@ extern "C" { #include "BPH_mass_spring.h" #include "implicit.h" +#include "DEG_depsgraph.h" +#include "DEG_depsgraph_query.h" + static float I3[3][3] = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}; /* Number of off-diagonal non-zero matrix blocks. @@ -64,7 +67,7 @@ static int cloth_count_nondiag_blocks(Cloth *cloth) for (link = cloth->springs; link; link = link->next) { ClothSpring *spring = (ClothSpring *)link->link; switch (spring->type) { - case CLOTH_SPRING_TYPE_BENDING_ANG: + case CLOTH_SPRING_TYPE_BENDING_HAIR: /* angular bending combines 3 vertices */ nondiag += 3; break; @@ -338,38 +341,75 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s) Cloth *cloth = clmd->clothObject; ClothSimSettings *parms = clmd->sim_parms; Implicit_Data *data = cloth->implicit; - - bool no_compress = parms->flags & CLOTH_SIMSETTINGS_FLAG_NO_SPRING_COMPRESS; + bool using_angular = parms->bending_model == CLOTH_BENDING_ANGULAR; + bool resist_compress = (parms->flags & CLOTH_SIMSETTINGS_FLAG_RESIST_SPRING_COMPRESS) && !using_angular; s->flags &= ~CLOTH_SPRING_FLAG_NEEDED; - // calculate force of structural + shear springs - if ((s->type & CLOTH_SPRING_TYPE_STRUCTURAL) || (s->type & CLOTH_SPRING_TYPE_SHEAR) || (s->type & CLOTH_SPRING_TYPE_SEWING) ) { -#ifdef CLOTH_FORCE_SPRING_STRUCTURAL + /* Calculate force of bending springs. */ + if ((s->type & CLOTH_SPRING_TYPE_BENDING) && using_angular) { +#ifdef CLOTH_FORCE_SPRING_BEND float k, scaling; s->flags |= CLOTH_SPRING_FLAG_NEEDED; - scaling = parms->structural + s->stiffness * fabsf(parms->max_struct - parms->structural); - k = scaling / (parms->avg_spring_len + FLT_EPSILON); + scaling = parms->bending + s->ang_stiffness * fabsf(parms->max_bend - parms->bending); + k = scaling * s->restlen * 0.1f; /* Multiplying by 0.1, just to scale the forces to more reasonable values. */ + + BPH_mass_spring_force_spring_angular(data, s->ij, s->kl, s->pa, s->pb, s->la, s->lb, + s->restang, k, parms->bending_damping); +#endif + } + + /* Calculate force of structural + shear springs. */ + if (s->type & (CLOTH_SPRING_TYPE_STRUCTURAL | CLOTH_SPRING_TYPE_SEWING)) { +#ifdef CLOTH_FORCE_SPRING_STRUCTURAL + float k_tension, scaling_tension; + + s->flags |= CLOTH_SPRING_FLAG_NEEDED; + + scaling_tension = parms->tension + s->lin_stiffness * fabsf(parms->max_tension - parms->tension); + k_tension = scaling_tension / (parms->avg_spring_len + FLT_EPSILON); if (s->type & CLOTH_SPRING_TYPE_SEWING) { // TODO: verify, half verified (couldn't see error) // sewing springs usually have a large distance at first so clamp the force so we don't get tunnelling through colission objects - BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->Cdis, no_compress, parms->max_sewing); + BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, + k_tension, parms->tension_damp, + 0.0f, 0.0f, false, false, parms->max_sewing); } else { - BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->Cdis, no_compress, 0.0f); + float k_compression, scaling_compression; + scaling_compression = parms->compression + s->lin_stiffness * fabsf(parms->max_compression - parms->compression); + k_compression = scaling_compression / (parms->avg_spring_len + FLT_EPSILON); + + BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, + k_tension, parms->tension_damp, + k_compression, parms->compression_damp, + resist_compress, using_angular, 0.0f); } #endif } + else if (s->type & CLOTH_SPRING_TYPE_SHEAR) { +#ifdef CLOTH_FORCE_SPRING_SHEAR + float k, scaling; + + s->flags |= CLOTH_SPRING_FLAG_NEEDED; + + scaling = parms->shear + s->lin_stiffness * fabsf(parms->max_shear - parms->shear); + k = scaling / (parms->avg_spring_len + FLT_EPSILON); + + BPH_mass_spring_force_spring_linear(data, s->ij, s->kl, s->restlen, k, parms->shear_damp, + 0.0f, 0.0f, resist_compress, false, 0.0f); +#endif + } else if (s->type & CLOTH_SPRING_TYPE_BENDING) { /* calculate force of bending springs */ #ifdef CLOTH_FORCE_SPRING_BEND float kb, cb, scaling; s->flags |= CLOTH_SPRING_FLAG_NEEDED; - scaling = parms->bending + s->stiffness * fabsf(parms->max_bend - parms->bending); + scaling = parms->bending + s->lin_stiffness * fabsf(parms->max_bend - parms->bending); kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON)); // Fix for [#45084] for cloth stiffness must have cb proportional to kb @@ -378,7 +418,7 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s) BPH_mass_spring_force_spring_bending(data, s->ij, s->kl, s->restlen, kb, cb); #endif } - else if (s->type & CLOTH_SPRING_TYPE_BENDING_ANG) { + else if (s->type & CLOTH_SPRING_TYPE_BENDING_HAIR) { #ifdef CLOTH_FORCE_SPRING_BEND float kb, cb, scaling; @@ -388,14 +428,14 @@ BLI_INLINE void cloth_calc_spring_force(ClothModifierData *clmd, ClothSpring *s) * this is crap, but needed due to cloth/hair mixing ... * max_bend factor is not even used for hair, so ... */ - scaling = s->stiffness * parms->bending; + scaling = s->lin_stiffness * parms->bending; kb = scaling / (20.0f * (parms->avg_spring_len + FLT_EPSILON)); // Fix for [#45084] for cloth stiffness must have cb proportional to kb cb = kb * parms->bending_damping; /* XXX assuming same restlen for ij and jk segments here, this can be done correctly for hair later */ - BPH_mass_spring_force_spring_bending_angular(data, s->ij, s->kl, s->mn, s->target, kb, cb); + BPH_mass_spring_force_spring_bending_hair(data, s->ij, s->kl, s->mn, s->target, kb, cb); #if 0 { @@ -433,7 +473,7 @@ static void hair_get_boundbox(ClothModifierData *clmd, float gmin[3], float gmax } } -static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListBase *effectors, float time) +static void cloth_calc_force(Scene *scene, ClothModifierData *clmd, float UNUSED(frame), ListBase *effectors, float time) { /* Collect forces and derivatives: F, dFdX, dFdV */ Cloth *cloth = clmd->clothObject; @@ -447,9 +487,9 @@ static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListB #ifdef CLOTH_FORCE_GRAVITY /* global acceleration (gravitation) */ - if (clmd->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) { + if (scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY) { /* scale gravity force */ - mul_v3_v3fl(gravity, clmd->scene->physics_settings.gravity, 0.001f * clmd->sim_parms->effector_weights->global_gravity); + mul_v3_v3fl(gravity, scene->physics_settings.gravity, 0.001f * clmd->sim_parms->effector_weights->global_gravity); } vert = cloth->verts; @@ -487,8 +527,8 @@ static void cloth_calc_force(ClothModifierData *clmd, float UNUSED(frame), ListB EffectedPoint epoint; BPH_mass_spring_get_motion_state(data, i, x, v); - pd_point_from_loc(clmd->scene, x, v, i, &epoint); - pdDoEffectors(effectors, NULL, clmd->sim_parms->effector_weights, &epoint, winvec[i], NULL); + pd_point_from_loc(scene, x, v, i, &epoint); + BKE_effectors_apply(effectors, NULL, clmd->sim_parms->effector_weights, &epoint, winvec[i], NULL); } for (i = 0; i < cloth->tri_num; i++) { @@ -844,87 +884,41 @@ static void cloth_calc_volume_force(ClothModifierData *clmd) } #endif -/* old collision stuff for cloth, use for continuity - * until a good replacement is ready - */ -static void cloth_collision_solve_extra(Object *ob, ClothModifierData *clmd, ListBase *effectors, float frame, float step, float dt) +static void cloth_solve_collisions(Depsgraph *depsgraph, Object *ob, ClothModifierData *clmd, float step, float dt) { Cloth *cloth = clmd->clothObject; Implicit_Data *id = cloth->implicit; ClothVertex *verts = cloth->verts; int mvert_num = cloth->mvert_num; - const float spf = (float)clmd->sim_parms->stepsPerFrame / clmd->sim_parms->timescale; - - bool do_extra_solve; + const float time_multiplier = 1.0f / (clmd->sim_parms->dt * clmd->sim_parms->timescale); int i; - if (!(clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_ENABLED)) + if (!(clmd->coll_parms->flags & (CLOTH_COLLSETTINGS_FLAG_ENABLED | CLOTH_COLLSETTINGS_FLAG_SELF))) return; + if (!clmd->clothObject->bvhtree) return; - // update verts to current positions + BPH_mass_spring_solve_positions(id, dt); + + /* Update verts to current positions. */ for (i = 0; i < mvert_num; i++) { BPH_mass_spring_get_new_position(id, i, verts[i].tx); sub_v3_v3v3(verts[i].tv, verts[i].tx, verts[i].txold); - copy_v3_v3(verts[i].v, verts[i].tv); + zero_v3(verts[i].dcvel); } -#if 0 /* unused */ - for (i = 0, cv = cloth->verts; i < cloth->mvert_num; i++, cv++) { - copy_v3_v3(initial_cos[i], cv->tx); - } -#endif - - // call collision function - // TODO: check if "step" or "step+dt" is correct - dg - do_extra_solve = cloth_bvh_objcollision(ob, clmd, step / clmd->sim_parms->timescale, dt / clmd->sim_parms->timescale); - - // copy corrected positions back to simulation - for (i = 0; i < mvert_num; i++) { - float curx[3]; - BPH_mass_spring_get_position(id, i, curx); - // correct velocity again, just to be sure we had to change it due to adaptive collisions - sub_v3_v3v3(verts[i].tv, verts[i].tx, curx); - } - - if (do_extra_solve) { -// cloth_calc_helper_forces(ob, clmd, initial_cos, step/clmd->sim_parms->timescale, dt/clmd->sim_parms->timescale); - + if (cloth_bvh_collision(depsgraph, ob, clmd, step / clmd->sim_parms->timescale, dt / clmd->sim_parms->timescale)) { for (i = 0; i < mvert_num; i++) { - - float newv[3]; - - if ((clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) && (verts[i].flags & CLOTH_VERT_FLAG_PINNED)) + if ((clmd->sim_parms->vgroup_mass > 0) && (verts[i].flags & CLOTH_VERT_FLAG_PINNED)) continue; - BPH_mass_spring_set_new_position(id, i, verts[i].tx); - mul_v3_v3fl(newv, verts[i].tv, spf); - BPH_mass_spring_set_new_velocity(id, i, newv); + BPH_mass_spring_get_new_velocity(id, i, verts[i].tv); + madd_v3_v3fl(verts[i].tv, verts[i].dcvel, time_multiplier); + BPH_mass_spring_set_new_velocity(id, i, verts[i].tv); } } - - // X = Xnew; - BPH_mass_spring_apply_result(id); - - if (do_extra_solve) { - ImplicitSolverResult result; - - /* initialize forces to zero */ - BPH_mass_spring_clear_forces(id); - - // calculate forces - cloth_calc_force(clmd, frame, effectors, step); - - // calculate new velocity and position - BPH_mass_spring_solve_velocities(id, dt, &result); -// cloth_record_result(clmd, &result, clmd->sim_parms->stepsPerFrame); - - /* note: positions are advanced only once in the main solver step! */ - - BPH_mass_spring_apply_result(id); - } } static void cloth_clear_result(ClothModifierData *clmd) @@ -937,7 +931,7 @@ static void cloth_clear_result(ClothModifierData *clmd) sres->avg_iterations = 0.0f; } -static void cloth_record_result(ClothModifierData *clmd, ImplicitSolverResult *result, int steps) +static void cloth_record_result(ClothModifierData *clmd, ImplicitSolverResult *result, float dt) { ClothSolverResult *sres = clmd->solver_result; @@ -946,33 +940,34 @@ static void cloth_record_result(ClothModifierData *clmd, ImplicitSolverResult *r if (result->status == BPH_SOLVER_SUCCESS) { sres->min_error = min_ff(sres->min_error, result->error); sres->max_error = max_ff(sres->max_error, result->error); - sres->avg_error += result->error / (float)steps; + sres->avg_error += result->error * dt; } sres->min_iterations = min_ii(sres->min_iterations, result->iterations); sres->max_iterations = max_ii(sres->max_iterations, result->iterations); - sres->avg_iterations += (float)result->iterations / (float)steps; + sres->avg_iterations += (float)result->iterations * dt; } else { /* error only makes sense for successful iterations */ if (result->status == BPH_SOLVER_SUCCESS) { sres->min_error = sres->max_error = result->error; - sres->avg_error += result->error / (float)steps; + sres->avg_error += result->error * dt; } sres->min_iterations = sres->max_iterations = result->iterations; - sres->avg_iterations += (float)result->iterations / (float)steps; + sres->avg_iterations += (float)result->iterations * dt; } sres->status |= result->status; } -int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase *effectors) +int BPH_cloth_solve(Depsgraph *depsgraph, Object *ob, float frame, ClothModifierData *clmd, ListBase *effectors) { /* Hair currently is a cloth sim in disguise ... * Collision detection and volumetrics work differently then. * Bad design, TODO */ + Scene *scene = DEG_get_evaluated_scene(depsgraph); const bool is_hair = (clmd->hairdata != NULL); unsigned int i = 0; @@ -980,7 +975,7 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase * Cloth *cloth = clmd->clothObject; ClothVertex *verts = cloth->verts /*, *cv*/; unsigned int mvert_num = cloth->mvert_num; - float dt = clmd->sim_parms->timescale / clmd->sim_parms->stepsPerFrame; + float dt = clmd->sim_parms->dt * clmd->sim_parms->timescale; Implicit_Data *id = cloth->implicit; ColliderContacts *contacts = NULL; int totcolliders = 0; @@ -991,7 +986,7 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase * clmd->solver_result = (ClothSolverResult *)MEM_callocN(sizeof(ClothSolverResult), "cloth solver result"); cloth_clear_result(clmd); - if (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) { /* do goal stuff */ + if (clmd->sim_parms->vgroup_mass > 0) { /* Do goal stuff. */ for (i = 0; i < mvert_num; i++) { // update velocities with constrained velocities from pinned verts if (verts[i].flags & CLOTH_VERT_FLAG_PINNED) { @@ -1008,16 +1003,10 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase * while (step < tf) { ImplicitSolverResult result; - /* copy velocities for collision */ - for (i = 0; i < mvert_num; i++) { - BPH_mass_spring_get_motion_state(id, i, NULL, verts[i].tv); - copy_v3_v3(verts[i].v, verts[i].tv); - } - if (is_hair) { /* determine contact points */ if (clmd->coll_parms->flags & CLOTH_COLLSETTINGS_FLAG_ENABLED) { - cloth_find_point_contacts(ob, clmd, 0.0f, tf, &contacts, &totcolliders); + cloth_find_point_contacts(depsgraph, ob, clmd, 0.0f, tf, &contacts, &totcolliders); } /* setup vertex constraints for pinned vertices and contacts */ @@ -1031,39 +1020,28 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase * /* initialize forces to zero */ BPH_mass_spring_clear_forces(id); - // damping velocity for artistic reasons - // this is a bad way to do it, should be removed imo - lukas_t - if (clmd->sim_parms->vel_damping != 1.0f) { - for (i = 0; i < mvert_num; i++) { - float v[3]; - BPH_mass_spring_get_motion_state(id, i, NULL, v); - mul_v3_fl(v, clmd->sim_parms->vel_damping); - BPH_mass_spring_set_velocity(id, i, v); - } - } - // calculate forces - cloth_calc_force(clmd, frame, effectors, step); + cloth_calc_force(scene, clmd, frame, effectors, step); // calculate new velocity and position BPH_mass_spring_solve_velocities(id, dt, &result); - cloth_record_result(clmd, &result, clmd->sim_parms->stepsPerFrame); + cloth_record_result(clmd, &result, dt); + + /* Calculate collision impulses. */ + if (!is_hair) { + cloth_solve_collisions(depsgraph, ob, clmd, step, dt); + } if (is_hair) { cloth_continuum_step(clmd, dt); } BPH_mass_spring_solve_positions(id, dt); - - if (!is_hair) { - cloth_collision_solve_extra(ob, clmd, effectors, frame, step, dt); - } - BPH_mass_spring_apply_result(id); /* move pinned verts to correct position */ for (i = 0; i < mvert_num; i++) { - if (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_GOAL) { + if (clmd->sim_parms->vgroup_mass > 0) { if (verts[i].flags & CLOTH_VERT_FLAG_PINNED) { float x[3]; /* divide by time_scale to prevent pinned vertices' delta locations from being multiplied */ @@ -1091,24 +1069,3 @@ int BPH_cloth_solve(Object *ob, float frame, ClothModifierData *clmd, ListBase * return 1; } - -bool BPH_cloth_solver_get_texture_data(Object *UNUSED(ob), ClothModifierData *clmd, VoxelData *vd) -{ - Cloth *cloth = clmd->clothObject; - HairGrid *grid; - float gmin[3], gmax[3]; - - if (!clmd->clothObject || !clmd->clothObject->implicit) - return false; - - hair_get_boundbox(clmd, gmin, gmax); - - grid = BPH_hair_volume_create_vertex_grid(clmd->sim_parms->voxel_cell_size, gmin, gmax); - cloth_continuum_fill_grid(grid, cloth); - - BPH_hair_volume_get_texture_data(grid, vd); - - BPH_hair_volume_free_vertex_grid(grid); - - return true; -} diff --git a/source/blender/physics/intern/hair_volume.cpp b/source/blender/physics/intern/hair_volume.cpp index 8496d95855d..5746bc123ea 100644 --- a/source/blender/physics/intern/hair_volume.cpp +++ b/source/blender/physics/intern/hair_volume.cpp @@ -1064,7 +1064,7 @@ static HairGridVert *hair_volume_create_collision_grid(ClothModifierData *clmd, } /* gather colliders */ - colliders = get_collider_cache(clmd->scene, NULL, NULL); + colliders = BKE_collider_cache_create(depsgraph, NULL, NULL); if (colliders && collfac > 0.0f) { for (col = colliders->first; col; col = col->next) { MVert *loc0 = col->collmd->x; @@ -1097,7 +1097,7 @@ static HairGridVert *hair_volume_create_collision_grid(ClothModifierData *clmd, } } } - free_collider_cache(&colliders); + BKE_collider_cache_free(&colliders); /* divide velocity with density */ for (i = 0; i < size; i++) { @@ -1109,56 +1109,3 @@ static HairGridVert *hair_volume_create_collision_grid(ClothModifierData *clmd, return collgrid; } #endif - -bool BPH_hair_volume_get_texture_data(HairGrid *grid, VoxelData *vd) -{ - int totres, i; - int depth; - - vd->resol[0] = grid->res[0]; - vd->resol[1] = grid->res[1]; - vd->resol[2] = grid->res[2]; - - totres = hair_grid_size(grid->res); - - if (vd->hair_type == TEX_VD_HAIRVELOCITY) { - depth = 4; - vd->data_type = TEX_VD_RGBA_PREMUL; - } - else { - depth = 1; - vd->data_type = TEX_VD_INTENSITY; - } - - if (totres > 0) { - vd->dataset = (float *)MEM_mapallocN(sizeof(float) * depth * (totres), "hair volume texture data"); - - for (i = 0; i < totres; ++i) { - switch (vd->hair_type) { - case TEX_VD_HAIRDENSITY: - vd->dataset[i] = grid->verts[i].density; - break; - - case TEX_VD_HAIRRESTDENSITY: - vd->dataset[i] = 0.0f; // TODO - break; - - case TEX_VD_HAIRVELOCITY: { - vd->dataset[i + 0 * totres] = grid->verts[i].velocity[0]; - vd->dataset[i + 1 * totres] = grid->verts[i].velocity[1]; - vd->dataset[i + 2 * totres] = grid->verts[i].velocity[2]; - vd->dataset[i + 3 * totres] = len_v3(grid->verts[i].velocity); - break; - } - case TEX_VD_HAIRENERGY: - vd->dataset[i] = 0.0f; // TODO - break; - } - } - } - else { - vd->dataset = NULL; - } - - return true; -} diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h index 477bc704aff..11621ac812c 100644 --- a/source/blender/physics/intern/implicit.h +++ b/source/blender/physics/intern/implicit.h @@ -50,6 +50,7 @@ extern "C" { #define CLOTH_FORCE_GRAVITY #define CLOTH_FORCE_DRAG #define CLOTH_FORCE_SPRING_STRUCTURAL +#define CLOTH_FORCE_SPRING_SHEAR #define CLOTH_FORCE_SPRING_BEND #define CLOTH_FORCE_SPRING_GOAL #define CLOTH_FORCE_EFFECTORS @@ -114,12 +115,17 @@ void BPH_mass_spring_force_edge_wind(struct Implicit_Data *data, int v1, int v2, void BPH_mass_spring_force_vertex_wind(struct Implicit_Data *data, int v, float radius, const float(*winvec)[3]); /* Linear spring force between two points */ bool BPH_mass_spring_force_spring_linear(struct Implicit_Data *data, int i, int j, float restlen, - float stiffness, float damping, bool no_compress, float clamp_force); + float stiffness_tension, float damping_tension, + float stiffness_compression, float damping_compression, + bool resist_compress, bool new_compress, float clamp_force); +/* Angular spring force between two polygons */ +bool BPH_mass_spring_force_spring_angular(struct Implicit_Data *data, int i, int j, int *i_a, int *i_b, int len_a, int len_b, + float restang, float stiffness, float damping); /* Bending force, forming a triangle at the base of two structural springs */ bool BPH_mass_spring_force_spring_bending(struct Implicit_Data *data, int i, int j, float restlen, float kb, float cb); /* Angular bending force based on local target vectors */ -bool BPH_mass_spring_force_spring_bending_angular(struct Implicit_Data *data, int i, int j, int k, - const float target[3], float stiffness, float damping); +bool BPH_mass_spring_force_spring_bending_hair(struct Implicit_Data *data, int i, int j, int k, + const float target[3], float stiffness, float damping); /* Global goal spring */ bool BPH_mass_spring_force_spring_goal(struct Implicit_Data *data, int i, const float goal_x[3], const float goal_v[3], float stiffness, float damping); @@ -171,8 +177,6 @@ void BPH_hair_volume_vertex_grid_forces(struct HairGrid *grid, const float x[3], float smoothfac, float pressurefac, float minpressure, float f[3], float dfdx[3][3], float dfdv[3][3]); -bool BPH_hair_volume_get_texture_data(struct HairGrid *grid, struct VoxelData *vd); - #ifdef __cplusplus } #endif diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c index 8ee9513e81b..28546f8ca0d 100644 --- a/source/blender/physics/intern/implicit_blender.c +++ b/source/blender/physics/intern/implicit_blender.c @@ -461,6 +461,13 @@ DO_INLINE void muladd_fmatrix_fvector(float to[3], float matrix[3][3], float fro to[2] += dot_v3v3(matrix[2], from); } +DO_INLINE void muladd_fmatrixT_fvector(float to[3], float matrix[3][3], float from[3]) +{ + to[0] += matrix[0][0] * from[0] + matrix[1][0] * from[1] + matrix[2][0] * from[2]; + to[1] += matrix[0][1] * from[0] + matrix[1][1] * from[1] + matrix[2][1] * from[2]; + to[2] += matrix[0][2] * from[0] + matrix[1][2] * from[1] + matrix[2][2] * from[2]; +} + BLI_INLINE void outerproduct(float r[3][3], const float a[3], const float b[3]) { mul_v3_v3fl(r[0], a, b[0]); @@ -604,7 +611,9 @@ DO_INLINE void mul_bfmatrix_lfvector(float(*to)[3], fmatrix3x3 *from, lfVector * #pragma omp section { for (i = from[0].vcount; i < from[0].vcount + from[0].scount; i++) { - muladd_fmatrix_fvector(to[from[i].c], from[i].m, fLongVector[from[i].r]); + /* This is the lower triangle of the sparse matrix, + * therefore multiplication occurs with transposed submatrices. */ + muladd_fmatrixT_fvector(to[from[i].c], from[i].m, fLongVector[from[i].r]); } } #pragma omp section @@ -617,8 +626,6 @@ DO_INLINE void mul_bfmatrix_lfvector(float(*to)[3], fmatrix3x3 *from, lfVector * add_lfvector_lfvector(to, to, temp, from[0].vcount); del_lfvector(temp); - - } /* SPARSE SYMMETRIC sub big matrix with big matrix*/ @@ -1585,9 +1592,13 @@ BLI_INLINE void apply_spring(Implicit_Data *data, int i, int j, const float f[3] } bool BPH_mass_spring_force_spring_linear(Implicit_Data *data, int i, int j, float restlen, - float stiffness, float damping, bool no_compress, float clamp_force) + float stiffness_tension, float damping_tension, + float stiffness_compression, float damping_compression, + bool resist_compress, bool new_compress, float clamp_force) { float extent[3], length, dir[3], vel[3]; + float f[3], dfdx[3][3], dfdv[3][3]; + float damping = 0; // calculate elonglation spring_length(data, i, j, extent, dir, &length, vel); @@ -1595,29 +1606,41 @@ bool BPH_mass_spring_force_spring_linear(Implicit_Data *data, int i, int j, floa /* This code computes not only the force, but also its derivative. Zero derivative effectively disables the spring for the implicit solver. Thus length > restlen makes cloth unconstrained at the start of simulation. */ - if ((length >= restlen && length > 0) || no_compress) { - float stretch_force, f[3], dfdx[3][3], dfdv[3][3]; + if ((length >= restlen && length > 0) || resist_compress) { + float stretch_force; + + damping = damping_tension; - stretch_force = stiffness * (length - restlen); + stretch_force = stiffness_tension * (length - restlen); if (clamp_force > 0.0f && stretch_force > clamp_force) { stretch_force = clamp_force; } mul_v3_v3fl(f, dir, stretch_force); - // Ascher & Boxman, p.21: Damping only during elonglation - // something wrong with it... - madd_v3_v3fl(f, dir, damping * dot_v3v3(vel, dir)); + dfdx_spring(dfdx, dir, length, restlen, stiffness_tension); + } + else if (new_compress) { + /* This is based on the Choi and Ko bending model, which works surprisingly well for compression. */ + float kb = stiffness_compression; + float cb = kb; /* cb equal to kb seems to work, but a factor can be added if necessary */ - dfdx_spring(dfdx, dir, length, restlen, stiffness); - dfdv_damp(dfdv, dir, damping); + damping = damping_compression; - apply_spring(data, i, j, f, dfdx, dfdv); + mul_v3_v3fl(f, dir, fbstar(length, restlen, kb, cb)); - return true; + outerproduct(dfdx, dir, dir); + mul_m3_fl(dfdx, fbstar_jacobi(length, restlen, kb, cb)); } else { return false; } + + madd_v3_v3fl(f, dir, damping * dot_v3v3(vel, dir)); + dfdv_damp(dfdv, dir, damping); + + apply_spring(data, i, j, f, dfdx, dfdv); + + return true; } /* See "Stable but Responsive Cloth" (Choi, Ko 2005) */ @@ -1648,6 +1671,114 @@ bool BPH_mass_spring_force_spring_bending(Implicit_Data *data, int i, int j, flo } } +BLI_INLINE void poly_avg(lfVector *data, int *inds, int len, float r_avg[3]) +{ + float fact = 1.0f / (float)len; + + zero_v3(r_avg); + + for (int i = 0; i < len; i++) { + madd_v3_v3fl(r_avg, data[inds[i]], fact); + } +} + +BLI_INLINE void poly_norm(lfVector *data, int i, int j, int *inds, int len, float r_dir[3]) +{ + float mid[3]; + + poly_avg(data, inds, len, mid); + + normal_tri_v3(r_dir, data[i], data[j], mid); +} + +BLI_INLINE void edge_avg(lfVector *data, int i, int j, float r_avg[3]) +{ + r_avg[0] = (data[i][0] + data[j][0]) * 0.5f; + r_avg[1] = (data[i][1] + data[j][1]) * 0.5f; + r_avg[2] = (data[i][2] + data[j][2]) * 0.5f; +} + +BLI_INLINE void edge_norm(lfVector *data, int i, int j, float r_dir[3]) +{ + sub_v3_v3v3(r_dir, data[i], data[j]); + normalize_v3(r_dir); +} + +BLI_INLINE float bend_angle(float dir_a[3], float dir_b[3], float dir_e[3]) +{ + float cos, sin; + float tmp[3]; + + cos = dot_v3v3(dir_a, dir_b); + + cross_v3_v3v3(tmp, dir_a, dir_b); + sin = dot_v3v3(tmp, dir_e); + + return atan2f(sin, cos); +} + +BLI_INLINE void spring_angle(Implicit_Data *data, int i, int j, int *i_a, int *i_b, int len_a, int len_b, + float r_dir_a[3], float r_dir_b[3], + float *r_angle, float r_vel_a[3], float r_vel_b[3]) +{ + float dir_e[3], vel_e[3]; + + poly_norm(data->X, j, i, i_a, len_a, r_dir_a); + poly_norm(data->X, i, j, i_b, len_b, r_dir_b); + + edge_norm(data->X, i, j, dir_e); + + *r_angle = bend_angle(r_dir_a, r_dir_b, dir_e); + + poly_avg(data->V, i_a, len_a, r_vel_a); + poly_avg(data->V, i_b, len_b, r_vel_b); + + edge_avg(data->V, i, j, vel_e); + + sub_v3_v3(r_vel_a, vel_e); + sub_v3_v3(r_vel_b, vel_e); +} + +/* Angular springs roughly based on the bending model proposed by Baraff and Witkin in "Large Steps in Cloth Simulation". */ +bool BPH_mass_spring_force_spring_angular(Implicit_Data *data, int i, int j, int *i_a, int *i_b, int len_a, int len_b, + float restang, float stiffness, float damping) +{ + float angle, dir_a[3], dir_b[3], vel_a[3], vel_b[3]; + float f_a[3], f_b[3], f_e[3]; + float force; + int x; + + spring_angle(data, i, j, i_a, i_b, len_a, len_b, + dir_a, dir_b, &angle, vel_a, vel_b); + + /* spring force */ + force = stiffness * (angle - restang); + + /* damping force */ + force += -damping * (dot_v3v3(vel_a, dir_a) + dot_v3v3(vel_b, dir_b)); + + mul_v3_v3fl(f_a, dir_a, force / len_a); + mul_v3_v3fl(f_b, dir_b, force / len_b); + + for (x = 0; x < len_a; x++) { + add_v3_v3(data->F[i_a[x]], f_a); + } + + for (x = 0; x < len_b; x++) { + add_v3_v3(data->F[i_b[x]], f_b); + } + + mul_v3_v3fl(f_a, dir_a, force * 0.5f); + mul_v3_v3fl(f_b, dir_b, force * 0.5f); + + add_v3_v3v3(f_e, f_a, f_b); + + sub_v3_v3(data->F[i], f_e); + sub_v3_v3(data->F[j], f_e); + + return true; +} + /* Jacobian of a direction vector. * Basically the part of the differential orthogonal to the direction, * inversely proportional to the length of the edge. @@ -1671,11 +1802,11 @@ BLI_INLINE void spring_grad_dir(Implicit_Data *data, int i, int j, float edge[3] } } -BLI_INLINE void spring_angbend_forces(Implicit_Data *data, int i, int j, int k, - const float goal[3], - float stiffness, float damping, - int q, const float dx[3], const float dv[3], - float r_f[3]) +BLI_INLINE void spring_hairbend_forces(Implicit_Data *data, int i, int j, int k, + const float goal[3], + float stiffness, float damping, + int q, const float dx[3], const float dv[3], + float r_f[3]) { float edge_ij[3], dir_ij[3]; float edge_jk[3], dir_jk[3]; @@ -1720,10 +1851,10 @@ BLI_INLINE void spring_angbend_forces(Implicit_Data *data, int i, int j, int k, } /* Finite Differences method for estimating the jacobian of the force */ -BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data, int i, int j, int k, - const float goal[3], - float stiffness, float damping, - int q, float dfdx[3][3]) +BLI_INLINE void spring_hairbend_estimate_dfdx(Implicit_Data *data, int i, int j, int k, + const float goal[3], + float stiffness, float damping, + int q, float dfdx[3][3]) { const float delta = 0.00001f; // TODO find a good heuristic for this float dvec_null[3][3], dvec_pos[3][3], dvec_neg[3][3]; @@ -1739,12 +1870,12 @@ BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data, int i, int j, /* XXX TODO offset targets to account for position dependency */ for (a = 0; a < 3; ++a) { - spring_angbend_forces(data, i, j, k, goal, stiffness, damping, - q, dvec_pos[a], dvec_null[a], f); + spring_hairbend_forces(data, i, j, k, goal, stiffness, damping, + q, dvec_pos[a], dvec_null[a], f); copy_v3_v3(dfdx[a], f); - spring_angbend_forces(data, i, j, k, goal, stiffness, damping, - q, dvec_neg[a], dvec_null[a], f); + spring_hairbend_forces(data, i, j, k, goal, stiffness, damping, + q, dvec_neg[a], dvec_null[a], f); sub_v3_v3(dfdx[a], f); for (b = 0; b < 3; ++b) { @@ -1754,10 +1885,10 @@ BLI_INLINE void spring_angbend_estimate_dfdx(Implicit_Data *data, int i, int j, } /* Finite Differences method for estimating the jacobian of the force */ -BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data, int i, int j, int k, - const float goal[3], - float stiffness, float damping, - int q, float dfdv[3][3]) +BLI_INLINE void spring_hairbend_estimate_dfdv(Implicit_Data *data, int i, int j, int k, + const float goal[3], + float stiffness, float damping, + int q, float dfdv[3][3]) { const float delta = 0.00001f; // TODO find a good heuristic for this float dvec_null[3][3], dvec_pos[3][3], dvec_neg[3][3]; @@ -1773,12 +1904,12 @@ BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data, int i, int j, /* XXX TODO offset targets to account for position dependency */ for (a = 0; a < 3; ++a) { - spring_angbend_forces(data, i, j, k, goal, stiffness, damping, - q, dvec_null[a], dvec_pos[a], f); + spring_hairbend_forces(data, i, j, k, goal, stiffness, damping, + q, dvec_null[a], dvec_pos[a], f); copy_v3_v3(dfdv[a], f); - spring_angbend_forces(data, i, j, k, goal, stiffness, damping, - q, dvec_null[a], dvec_neg[a], f); + spring_hairbend_forces(data, i, j, k, goal, stiffness, damping, + q, dvec_null[a], dvec_neg[a], f); sub_v3_v3(dfdv[a], f); for (b = 0; b < 3; ++b) { @@ -1790,8 +1921,8 @@ BLI_INLINE void spring_angbend_estimate_dfdv(Implicit_Data *data, int i, int j, /* Angular spring that pulls the vertex toward the local target * See "Artistic Simulation of Curly Hair" (Pixar technical memo #12-03a) */ -bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, int j, int k, - const float target[3], float stiffness, float damping) +bool BPH_mass_spring_force_spring_bending_hair(Implicit_Data *data, int i, int j, int k, + const float target[3], float stiffness, float damping) { float goal[3]; float fj[3], fk[3]; @@ -1806,18 +1937,18 @@ bool BPH_mass_spring_force_spring_bending_angular(Implicit_Data *data, int i, in world_to_root_v3(data, j, goal, target); - spring_angbend_forces(data, i, j, k, goal, stiffness, damping, k, vecnull, vecnull, fk); + spring_hairbend_forces(data, i, j, k, goal, stiffness, damping, k, vecnull, vecnull, fk); negate_v3_v3(fj, fk); /* counterforce */ - spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, i, dfk_dxi); - spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, j, dfk_dxj); - spring_angbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, k, dfk_dxk); + spring_hairbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, i, dfk_dxi); + spring_hairbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, j, dfk_dxj); + spring_hairbend_estimate_dfdx(data, i, j, k, goal, stiffness, damping, k, dfk_dxk); copy_m3_m3(dfj_dxi, dfk_dxi); negate_m3(dfj_dxi); copy_m3_m3(dfj_dxj, dfk_dxj); negate_m3(dfj_dxj); - spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, i, dfk_dvi); - spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, j, dfk_dvj); - spring_angbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, k, dfk_dvk); + spring_hairbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, i, dfk_dvi); + spring_hairbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, j, dfk_dvj); + spring_hairbend_estimate_dfdv(data, i, j, k, goal, stiffness, damping, k, dfk_dvk); copy_m3_m3(dfj_dvi, dfk_dvi); negate_m3(dfj_dvi); copy_m3_m3(dfj_dvj, dfk_dvj); negate_m3(dfj_dvj); |