Cloth: implement support for a hydrostatic pressure gradient.

When a fluid is put under influence of gravity or acceleration, it
forms an internal pressure gradient, which causes observable effects
like buoyancy. Since now cloth has support for simulating pressure
changes caused by fluid compression or expansion, it makes sense to
also support the effects of gravity.

This is intended for better simulation of objects filled or
surrounded by fluids, especially when constrained by collisions
or pinned vertices, and should result in more realistic shapes.

Obviously, this doesn't actually simulate fluid dynamics; instead
it is assumed that the fluid immediately adapts to changes in the
shape or acceleration of the object without friction or turbulence,
and instantly reaches a new static equilibrium.

Differential Revision: https://developer.blender.org/D6442

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],