diff options
Diffstat (limited to 'source/blender')
-rw-r--r-- | source/blender/blenkernel/BKE_cloth.h | 7 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/cloth.c | 1 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/pointcache.c | 62 | ||||
-rw-r--r-- | source/blender/blenloader/intern/writefile.c | 1 | ||||
-rw-r--r-- | source/blender/makesdna/DNA_cloth_types.h | 4 | ||||
-rw-r--r-- | source/blender/makesdna/DNA_pointcache_types.h | 1 | ||||
-rw-r--r-- | source/blender/makesrna/intern/rna_cloth.c | 32 | ||||
-rw-r--r-- | source/blender/physics/intern/BPH_mass_spring.cpp | 220 | ||||
-rw-r--r-- | source/blender/physics/intern/implicit.h | 7 | ||||
-rw-r--r-- | source/blender/physics/intern/implicit_blender.c | 52 |
10 files changed, 290 insertions, 97 deletions
diff --git a/source/blender/blenkernel/BKE_cloth.h b/source/blender/blenkernel/BKE_cloth.h index f5073ff7aa5..f25482570b1 100644 --- a/source/blender/blenkernel/BKE_cloth.h +++ b/source/blender/blenkernel/BKE_cloth.h @@ -93,9 +93,10 @@ typedef struct Cloth { struct Implicit_Data *implicit; /* our implicit solver connects to this pointer */ struct EdgeSet *edgeset; /* used for selfcollisions */ int last_frame; - float initial_mesh_volume; /* Initial volume of the mesh. Used for pressure */ - struct MEdge *edges; /* Used for hair collisions. */ - struct GHash *sew_edge_graph; /* Sewing edges represented using a GHash */ + float initial_mesh_volume; /* Initial volume of the mesh. Used for pressure */ + float average_acceleration[3]; /* Moving average of overall acceleration. */ + struct MEdge *edges; /* Used for hair collisions. */ + struct GHash *sew_edge_graph; /* Sewing edges represented using a GHash */ } Cloth; /** diff --git a/source/blender/blenkernel/intern/cloth.c b/source/blender/blenkernel/intern/cloth.c index 0b9780ac81c..ce6143116a1 100644 --- a/source/blender/blenkernel/intern/cloth.c +++ b/source/blender/blenkernel/intern/cloth.c @@ -146,6 +146,7 @@ void cloth_init(ClothModifierData *clmd) clmd->sim_parms->uniform_pressure_force = 0.0f; clmd->sim_parms->target_volume = 0.0f; clmd->sim_parms->pressure_factor = 1.0f; + clmd->sim_parms->fluid_density = 0.0f; clmd->sim_parms->vgroup_pressure = 0; // also from softbodies diff --git a/source/blender/blenkernel/intern/pointcache.c b/source/blender/blenkernel/intern/pointcache.c index 9cd17310f07..836fe59d7bd 100644 --- a/source/blender/blenkernel/intern/pointcache.c +++ b/source/blender/blenkernel/intern/pointcache.c @@ -138,6 +138,7 @@ static int ptcache_data_size[] = { static int ptcache_extra_datasize[] = { 0, sizeof(ParticleSpring), + sizeof(float) * 3, }; /* forward declarations */ @@ -176,6 +177,23 @@ static int ptcache_basic_header_write(PTCacheFile *pf) return 1; } +static void ptcache_add_extra_data(PTCacheMem *pm, + unsigned int type, + unsigned int count, + void *data) +{ + PTCacheExtra *extra = MEM_callocN(sizeof(PTCacheExtra), "Point cache: extra data descriptor"); + + extra->type = type; + extra->totdata = count; + + size_t size = extra->totdata * ptcache_extra_datasize[extra->type]; + + extra->data = MEM_mallocN(size, "Point cache: extra data"); + memcpy(extra->data, data, size); + + BLI_addtail(&pm->extradata, extra); +} /* Softbody functions */ static int ptcache_softbody_write(int index, void *soft_v, void **data, int UNUSED(cfra)) { @@ -467,21 +485,12 @@ static int ptcache_particle_totwrite(void *psys_v, int cfra) static void ptcache_particle_extra_write(void *psys_v, PTCacheMem *pm, int UNUSED(cfra)) { ParticleSystem *psys = psys_v; - PTCacheExtra *extra = NULL; if (psys->part->phystype == PART_PHYS_FLUID && psys->part->fluid && psys->part->fluid->flag & SPH_VISCOELASTIC_SPRINGS && psys->tot_fluidsprings && psys->fluid_springs) { - extra = MEM_callocN(sizeof(PTCacheExtra), "Point cache: fluid extra data"); - - extra->type = BPHYS_EXTRA_FLUID_SPRINGS; - extra->totdata = psys->tot_fluidsprings; - - extra->data = MEM_callocN(extra->totdata * ptcache_extra_datasize[extra->type], - "Point cache: extra data"); - memcpy(extra->data, psys->fluid_springs, extra->totdata * ptcache_extra_datasize[extra->type]); - - BLI_addtail(&pm->extradata, extra); + ptcache_add_extra_data( + pm, BPHYS_EXTRA_FLUID_SPRINGS, psys->tot_fluidsprings, psys->fluid_springs); } } @@ -575,6 +584,33 @@ static void ptcache_cloth_interpolate( /* should vert->xconst be interpolated somehow too? - jahka */ } +static void ptcache_cloth_extra_write(void *cloth_v, PTCacheMem *pm, int UNUSED(cfra)) +{ + ClothModifierData *clmd = cloth_v; + Cloth *cloth = clmd->clothObject; + + if (!is_zero_v3(cloth->average_acceleration)) { + ptcache_add_extra_data(pm, BPHYS_EXTRA_CLOTH_ACCELERATION, 1, cloth->average_acceleration); + } +} +static void ptcache_cloth_extra_read(void *cloth_v, PTCacheMem *pm, float UNUSED(cfra)) +{ + ClothModifierData *clmd = cloth_v; + Cloth *cloth = clmd->clothObject; + PTCacheExtra *extra = pm->extradata.first; + + zero_v3(cloth->average_acceleration); + + for (; extra; extra = extra->next) { + switch (extra->type) { + case BPHYS_EXTRA_CLOTH_ACCELERATION: { + copy_v3_v3(cloth->average_acceleration, extra->data); + break; + } + } + } +} + static int ptcache_cloth_totpoint(void *cloth_v, int UNUSED(cfra)) { ClothModifierData *clmd = cloth_v; @@ -1681,8 +1717,8 @@ void BKE_ptcache_id_from_cloth(PTCacheID *pid, Object *ob, ClothModifierData *cl pid->write_stream = NULL; pid->read_stream = NULL; - pid->write_extra_data = NULL; - pid->read_extra_data = NULL; + pid->write_extra_data = ptcache_cloth_extra_write; + pid->read_extra_data = ptcache_cloth_extra_read; pid->interpolate_extra_data = NULL; pid->write_header = ptcache_basic_header_write; diff --git a/source/blender/blenloader/intern/writefile.c b/source/blender/blenloader/intern/writefile.c index a3102200411..2b57a8e1f5e 100644 --- a/source/blender/blenloader/intern/writefile.c +++ b/source/blender/blenloader/intern/writefile.c @@ -1376,6 +1376,7 @@ static const char *ptcache_data_struct[] = { static const char *ptcache_extra_struct[] = { "", "ParticleSpring", + "vec3f", }; static void write_pointcaches(BlendWriter *writer, ListBase *ptcaches) { diff --git a/source/blender/makesdna/DNA_cloth_types.h b/source/blender/makesdna/DNA_cloth_types.h index 8f3a26cf9c0..2b66bb771e1 100644 --- a/source/blender/makesdna/DNA_cloth_types.h +++ b/source/blender/makesdna/DNA_cloth_types.h @@ -109,8 +109,10 @@ typedef struct ClothSimSettings { pressure=( (current_volume/target_volume) - 1 + uniform_pressure_force) * pressure_factor */ float pressure_factor; + /* Density of the fluid inside or outside the object for use in the hydrostatic pressure gradient. */ + float fluid_density; short vgroup_pressure; - char _pad7[2]; + char _pad7[6]; /* XXX various hair stuff * should really be separate, this struct is a horrible mess already diff --git a/source/blender/makesdna/DNA_pointcache_types.h b/source/blender/makesdna/DNA_pointcache_types.h index fcdc22ca56d..3c7fc9031de 100644 --- a/source/blender/makesdna/DNA_pointcache_types.h +++ b/source/blender/makesdna/DNA_pointcache_types.h @@ -50,6 +50,7 @@ extern "C" { #define BPHYS_TOT_DATA 8 #define BPHYS_EXTRA_FLUID_SPRINGS 1 +#define BPHYS_EXTRA_CLOTH_ACCELERATION 2 typedef struct PTCacheExtra { struct PTCacheExtra *next, *prev; diff --git a/source/blender/makesrna/intern/rna_cloth.c b/source/blender/makesrna/intern/rna_cloth.c index 70f219259ef..594b77ea1ad 100644 --- a/source/blender/makesrna/intern/rna_cloth.c +++ b/source/blender/makesrna/intern/rna_cloth.c @@ -975,8 +975,10 @@ static void rna_def_cloth_sim_settings(BlenderRNA *brna) prop = RNA_def_property(srna, "use_pressure_volume", PROP_BOOLEAN, PROP_NONE); RNA_def_property_boolean_sdna(prop, NULL, "flags", CLOTH_SIMSETTINGS_FLAG_PRESSURE_VOL); - RNA_def_property_ui_text( - prop, "Use Custom Volume", "Use the Volume parameter as the initial volume"); + RNA_def_property_ui_text(prop, + "Use Custom Volume", + "Use the Target Volume parameter as the initial volume, instead " + "of calculating it from the mesh itself"); RNA_def_property_clear_flag(prop, PROP_ANIMATABLE); RNA_def_property_update(prop, 0, "rna_cloth_update"); @@ -984,10 +986,10 @@ static void rna_def_cloth_sim_settings(BlenderRNA *brna) RNA_def_property_float_sdna(prop, NULL, "uniform_pressure_force"); RNA_def_property_range(prop, -10000.0f, 10000.0f); RNA_def_property_float_default(prop, 0.0f); - RNA_def_property_ui_text( - prop, - "Pressure", - "The uniform pressure that is constantly applied to the mesh. Can be negative"); + RNA_def_property_ui_text(prop, + "Pressure", + "The uniform pressure that is constantly applied to the mesh, in units " + "of Pressure Scale. Can be negative"); RNA_def_property_update(prop, 0, "rna_cloth_update"); prop = RNA_def_property(srna, "target_volume", PROP_FLOAT, PROP_NONE); @@ -997,14 +999,28 @@ static void rna_def_cloth_sim_settings(BlenderRNA *brna) RNA_def_property_ui_text(prop, "Target Volume", "The mesh volume where the inner/outer pressure will be the same. If " - "set to zero the volume will not contribute to the total pressure"); + "set to zero the change in volume will not affect pressure"); RNA_def_property_update(prop, 0, "rna_cloth_update"); prop = RNA_def_property(srna, "pressure_factor", PROP_FLOAT, PROP_NONE); RNA_def_property_float_sdna(prop, NULL, "pressure_factor"); RNA_def_property_range(prop, 0.0f, 10000.0f); RNA_def_property_float_default(prop, 1.0f); - RNA_def_property_ui_text(prop, "Pressure Scale", "Air pressure scaling factor"); + RNA_def_property_ui_text(prop, + "Pressure Scale", + "Ambient pressure (kPa) that balances out between the inside and " + "outside of the object when it has the target volume"); + RNA_def_property_update(prop, 0, "rna_cloth_update"); + + prop = RNA_def_property(srna, "fluid_density", PROP_FLOAT, PROP_NONE); + RNA_def_property_float_sdna(prop, NULL, "fluid_density"); + RNA_def_property_ui_range(prop, -2.0f, 2.0f, 0.05f, 4); + RNA_def_property_ui_text( + prop, + "Fluid Density", + "Density (kg/l) of the fluid contained inside the object, used to create " + "a hydrostatic pressure gradient simulating the weight of the internal fluid, " + "or buoyancy from the surrounding fluid if negative"); RNA_def_property_update(prop, 0, "rna_cloth_update"); prop = RNA_def_property(srna, "vertex_group_pressure", PROP_STRING, PROP_NONE); diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp index 6a951519730..18fab5215a6 100644 --- a/source/blender/physics/intern/BPH_mass_spring.cpp +++ b/source/blender/physics/intern/BPH_mass_spring.cpp @@ -72,12 +72,50 @@ static int cloth_count_nondiag_blocks(Cloth *cloth) return nondiag; } +static bool cloth_get_pressure_weights(ClothModifierData *clmd, + const MVertTri *vt, + float *r_weights) +{ + /* We have custom vertex weights for pressure. */ + if (clmd->sim_parms->vgroup_pressure > 0) { + Cloth *cloth = clmd->clothObject; + ClothVertex *verts = cloth->verts; + + for (unsigned int j = 0; j < 3; j++) { + r_weights[j] = verts[vt->tri[j]].pressure_factor; + + /* Skip the entire triangle if it has a zero weight. */ + if (r_weights[j] == 0.0f) { + return false; + } + } + } + + return true; +} + +static void cloth_calc_pressure_gradient(ClothModifierData *clmd, + const float gradient_vector[3], + float *r_vertex_pressure) +{ + Cloth *cloth = clmd->clothObject; + Implicit_Data *data = cloth->implicit; + unsigned int mvert_num = cloth->mvert_num; + float pt[3]; + + for (unsigned int i = 0; i < mvert_num; i++) { + BPH_mass_spring_get_position(data, i, pt); + r_vertex_pressure[i] = dot_v3v3(pt, gradient_vector); + } +} + static float cloth_calc_volume(ClothModifierData *clmd) { /* Calculate the (closed) cloth volume. */ Cloth *cloth = clmd->clothObject; const MVertTri *tri = cloth->tri; Implicit_Data *data = cloth->implicit; + float weights[3] = {1.0f, 1.0f, 1.0f}; float vol = 0; /* Early exit for hair, as it never has volume. */ @@ -85,38 +123,45 @@ static float cloth_calc_volume(ClothModifierData *clmd) return 0.0f; } - if (clmd->sim_parms->vgroup_pressure > 0) { - for (unsigned int i = 0; i < cloth->primitive_num; i++) { - bool skip_face = false; - /* We have custom vertex weights for pressure. */ - const MVertTri *vt = &tri[i]; - for (unsigned int j = 0; j < 3; j++) { - /* If any weight is zero, don't take this face into account for volume calculation. */ - ClothVertex *verts = clmd->clothObject->verts; - - if (verts[vt->tri[j]].pressure_factor == 0.0f) { - skip_face = true; - } - } - if (skip_face) { - continue; - } + for (unsigned int i = 0; i < cloth->primitive_num; i++) { + const MVertTri *vt = &tri[i]; + if (cloth_get_pressure_weights(clmd, vt, weights)) { vol += BPH_tri_tetra_volume_signed_6x(data, vt->tri[0], vt->tri[1], vt->tri[2]); } } - else { - for (unsigned int i = 0; i < cloth->primitive_num; i++) { - const MVertTri *vt = &tri[i]; - vol += BPH_tri_tetra_volume_signed_6x(data, vt->tri[0], vt->tri[1], vt->tri[2]); - } - } + /* We need to divide by 6 to get the actual volume. */ vol = vol / 6.0f; return vol; } +static float cloth_calc_average_pressure(ClothModifierData *clmd, const float *vertex_pressure) +{ + Cloth *cloth = clmd->clothObject; + const MVertTri *tri = cloth->tri; + Implicit_Data *data = cloth->implicit; + float weights[3] = {1.0f, 1.0f, 1.0f}; + float total_force = 0; + float total_area = 0; + + for (unsigned int i = 0; i < cloth->primitive_num; i++) { + const MVertTri *vt = &tri[i]; + + if (cloth_get_pressure_weights(clmd, vt, weights)) { + float area = BPH_tri_area(data, vt->tri[0], vt->tri[1], vt->tri[2]); + + total_force += (vertex_pressure[vt->tri[0]] + vertex_pressure[vt->tri[1]] + + vertex_pressure[vt->tri[2]]) * + area / 3.0f; + total_area += area; + } + } + + return total_force / total_area; +} + int BPH_cloth_solver_init(Object *UNUSED(ob), ClothModifierData *clmd) { Cloth *cloth = clmd->clothObject; @@ -554,6 +599,7 @@ static void cloth_calc_force( if ((parms->flags & CLOTH_SIMSETTINGS_FLAG_PRESSURE) && (clmd->hairdata == NULL)) { /* The difference in pressure between the inside and outside of the mesh.*/ float pressure_difference = 0.0f; + float volume_factor = 1.0f; float init_vol; if (parms->flags & CLOTH_SIMSETTINGS_FLAG_PRESSURE_VOL) { @@ -568,59 +614,71 @@ static void cloth_calc_force( float f; float vol = cloth_calc_volume(clmd); + /* If the volume is the same don't apply any pressure. */ + volume_factor = init_vol / vol; + pressure_difference = volume_factor - 1; + /* Calculate an artificial maximum value for cloth pressure. */ f = fabs(clmd->sim_parms->uniform_pressure_force) + 200.0f; /* Clamp the cloth pressure to the calculated maximum value. */ - if (vol * f < init_vol) { - pressure_difference = f; - } - else { - /* If the volume is the same don't apply any pressure. */ - pressure_difference = (init_vol / vol) - 1; - } + CLAMP_MAX(pressure_difference, f); } - pressure_difference += clmd->sim_parms->uniform_pressure_force; + pressure_difference += clmd->sim_parms->uniform_pressure_force; pressure_difference *= clmd->sim_parms->pressure_factor; - for (i = 0; i < cloth->primitive_num; i++) { - const MVertTri *vt = &tri[i]; - if (fabs(pressure_difference) > 1E-6f) { - if (clmd->sim_parms->vgroup_pressure > 0) { - /* We have custom vertex weights for pressure. */ - ClothVertex *verts = clmd->clothObject->verts; - int v1, v2, v3; - v1 = vt->tri[0]; - v2 = vt->tri[1]; - v3 = vt->tri[2]; - - float weights[3]; - bool skip_face = false; - - weights[0] = verts[v1].pressure_factor; - weights[1] = verts[v2].pressure_factor; - weights[2] = verts[v3].pressure_factor; - for (unsigned int j = 0; j < 3; j++) { - if (weights[j] == 0.0f) { - /* Exclude faces which has a zero weight vert. */ - skip_face = true; - break; - } - } - if (skip_face) { - continue; - } + /* Compute the hydrostatic pressure gradient if enabled. */ + float fluid_density = clmd->sim_parms->fluid_density * 1000; /* kg/l -> kg/m3 */ + float *hydrostatic_pressure = NULL; - BPH_mass_spring_force_pressure(data, v1, v2, v3, pressure_difference, weights); - } - else { - float weights[3] = {1.0f, 1.0f, 1.0f}; - BPH_mass_spring_force_pressure( - data, vt->tri[0], vt->tri[1], vt->tri[2], pressure_difference, weights); + if (fabs(fluid_density) > 1e-6f) { + float hydrostatic_vector[3]; + copy_v3_v3(hydrostatic_vector, gravity); + + /* When the fluid is inside the object, factor in the acceleration of + * the object into the pressure field, as gravity is indistinguishable + * from acceleration from the inside. */ + if (fluid_density > 0) { + sub_v3_v3(hydrostatic_vector, cloth->average_acceleration); + + /* Preserve the total mass by scaling density to match the change in volume. */ + fluid_density *= volume_factor; + } + + mul_v3_fl(hydrostatic_vector, fluid_density); + + /* Compute an array of per-vertex hydrostatic pressure, and subtract the average. */ + hydrostatic_pressure = (float *)MEM_mallocN(sizeof(float) * mvert_num, + "hydrostatic pressure gradient"); + + cloth_calc_pressure_gradient(clmd, hydrostatic_vector, hydrostatic_pressure); + + pressure_difference -= cloth_calc_average_pressure(clmd, hydrostatic_pressure); + } + + /* Apply pressure. */ + if (hydrostatic_pressure || fabs(pressure_difference) > 1E-6f) { + float weights[3] = {1.0f, 1.0f, 1.0f}; + + for (i = 0; i < cloth->primitive_num; i++) { + const MVertTri *vt = &tri[i]; + + if (cloth_get_pressure_weights(clmd, vt, weights)) { + BPH_mass_spring_force_pressure(data, + vt->tri[0], + vt->tri[1], + vt->tri[2], + pressure_difference, + hydrostatic_pressure, + weights); } } } + + if (hydrostatic_pressure) { + MEM_freeN(hydrostatic_pressure); + } } /* handle external forces like wind */ @@ -1032,6 +1090,30 @@ static void cloth_calc_volume_force(ClothModifierData *clmd) } #endif +static void cloth_calc_average_acceleration(ClothModifierData *clmd, float dt) +{ + Cloth *cloth = clmd->clothObject; + Implicit_Data *data = cloth->implicit; + int i, mvert_num = cloth->mvert_num; + float total[3] = {0.0f, 0.0f, 0.0f}; + + for (i = 0; i < mvert_num; i++) { + float v[3], nv[3]; + + BPH_mass_spring_get_velocity(data, i, v); + BPH_mass_spring_get_new_velocity(data, i, nv); + + sub_v3_v3(nv, v); + add_v3_v3(total, nv); + } + + mul_v3_fl(total, 1.0f / dt / mvert_num); + + /* Smooth the data using a running average to prevent instability. + * This is effectively an abstraction of the wave propagation speed in fluid. */ + interp_v3_v3v3(cloth->average_acceleration, total, cloth->average_acceleration, powf(0.25f, dt)); +} + static void cloth_solve_collisions( Depsgraph *depsgraph, Object *ob, ClothModifierData *clmd, float step, float dt) { @@ -1136,6 +1218,10 @@ int BPH_cloth_solve( float dt = clmd->sim_parms->dt * clmd->sim_parms->timescale; Implicit_Data *id = cloth->implicit; + /* Hydrostatic pressure gradient of the fluid inside the object is affected by acceleration. */ + bool use_acceleration = (clmd->sim_parms->flags & CLOTH_SIMSETTINGS_FLAG_PRESSURE) && + (clmd->sim_parms->fluid_density > 0); + BKE_sim_debug_data_clear_category("collision"); if (!clmd->solver_result) { @@ -1158,6 +1244,10 @@ int BPH_cloth_solve( } } + if (!use_acceleration) { + zero_v3(cloth->average_acceleration); + } + while (step < tf) { ImplicitSolverResult result; @@ -1181,6 +1271,10 @@ int BPH_cloth_solve( cloth_continuum_step(clmd, dt); } + if (use_acceleration) { + cloth_calc_average_acceleration(clmd, dt); + } + BPH_mass_spring_solve_positions(id, dt); BPH_mass_spring_apply_result(id); diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h index 490e727b5f2..69b50f7fa48 100644 --- a/source/blender/physics/intern/implicit.h +++ b/source/blender/physics/intern/implicit.h @@ -79,6 +79,7 @@ void BPH_mass_spring_get_motion_state(struct Implicit_Data *data, float x[3], float v[3]); void BPH_mass_spring_get_position(struct Implicit_Data *data, int index, float x[3]); +void BPH_mass_spring_get_velocity(struct Implicit_Data *data, int index, float v[3]); /* access to modified motion state during solver step */ void BPH_mass_spring_get_new_position(struct Implicit_Data *data, int index, float x[3]); @@ -183,13 +184,15 @@ bool BPH_mass_spring_force_spring_goal(struct Implicit_Data *data, float damping); float BPH_tri_tetra_volume_signed_6x(struct Implicit_Data *data, int v1, int v2, int v3); +float BPH_tri_area(struct Implicit_Data *data, int v1, int v2, int v3); void BPH_mass_spring_force_pressure(struct Implicit_Data *data, int v1, int v2, int v3, - float pressure_difference, - float weights[3]); + float common_pressure, + const float *vertex_pressure, + const float weights[3]); /* ======== Hair Volumetric Forces ======== */ diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c index b1afaa52ef2..063c224f158 100644 --- a/source/blender/physics/intern/implicit_blender.c +++ b/source/blender/physics/intern/implicit_blender.c @@ -1246,6 +1246,11 @@ void BPH_mass_spring_get_position(struct Implicit_Data *data, int index, float x root_to_world_v3(data, index, x, data->X[index]); } +void BPH_mass_spring_get_velocity(struct Implicit_Data *data, int index, float v[3]) +{ + root_to_world_v3(data, index, v, data->V[index]); +} + void BPH_mass_spring_get_new_position(struct Implicit_Data *data, int index, float x[3]) { root_to_world_v3(data, index, x, data->Xnew[index]); @@ -1488,21 +1493,54 @@ float BPH_tri_tetra_volume_signed_6x(Implicit_Data *data, int v1, int v2, int v3 return volume_tri_tetrahedron_signed_v3_6x(data->X[v1], data->X[v2], data->X[v3]); } -void BPH_mass_spring_force_pressure( - Implicit_Data *data, int v1, int v2, int v3, float pressure_difference, float weights[3]) +float BPH_tri_area(struct Implicit_Data *data, int v1, int v2, int v3) +{ + float nor[3]; + + return calc_nor_area_tri(nor, data->X[v1], data->X[v2], data->X[v3]); +} + +void BPH_mass_spring_force_pressure(Implicit_Data *data, + int v1, + int v2, + int v3, + float common_pressure, + const float *vertex_pressure, + const float weights[3]) { float nor[3], area; - float factor; + float factor, base_force; + float force[3]; /* calculate face normal and area */ area = calc_nor_area_tri(nor, data->X[v1], data->X[v2], data->X[v3]); /* The force is calculated and split up evenly for each of the three face verts */ - factor = pressure_difference * area / 3.0f; + factor = area / 3.0f; + base_force = common_pressure * factor; + + /* Compute per-vertex force values from local pressures. + * From integrating the pressure over the triangle and deriving + * equivalent vertex forces, it follows that: + * + * force[idx] = (sum(pressure) + pressure[idx]) * area / 12 + * + * Effectively, 1/4 of the pressure acts just on its vertex, + * while 3/4 is split evenly over all three. + */ + if (vertex_pressure) { + copy_v3_fl3(force, vertex_pressure[v1], vertex_pressure[v2], vertex_pressure[v3]); + mul_v3_fl(force, factor / 4.0f); + + base_force += force[0] + force[1] + force[2]; + } + else { + zero_v3(force); + } /* add pressure to each of the face verts */ - madd_v3_v3fl(data->F[v1], nor, factor * weights[0]); - madd_v3_v3fl(data->F[v2], nor, factor * weights[1]); - madd_v3_v3fl(data->F[v3], nor, factor * weights[2]); + madd_v3_v3fl(data->F[v1], nor, (base_force + force[0]) * weights[0]); + madd_v3_v3fl(data->F[v2], nor, (base_force + force[1]) * weights[1]); + madd_v3_v3fl(data->F[v3], nor, (base_force + force[2]) * weights[2]); } static void edge_wind_vertex(const float dir[3], |