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

3 files changed, 207 insertions, 72 deletions

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