Fix T30941: Add cloth air pressure simulation

This adds some basic simulation of internal air pressure inside of
closed cloth mesh objects.

Reviewed By: Jacques Lucke

Differential Revision: http://developer.blender.org/D5473

4 files changed, 99 insertions, 0 deletions

diff --git a/source/blender/physics/BPH_mass_spring.h b/source/blender/physics/BPH_mass_spring.h
index c0ceff4d8cf..5a8c78812a4 100644
--- a/source/blender/physics/BPH_mass_spring.h
+++ b/source/blender/physics/BPH_mass_spring.h
@@ -53,6 +53,7 @@ int BPH_cloth_solve(struct Depsgraph *depsgraph,
                     struct ClothModifierData *clmd,
                     struct ListBase *effectors);
 void BKE_cloth_solver_set_positions(struct ClothModifierData *clmd);
+void BKE_cloth_solver_set_volume(ClothModifierData *clmd);
 
 #ifdef __cplusplus
 }
diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp
index c057e74b72b..259eed88756 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -74,6 +74,25 @@ static int cloth_count_nondiag_blocks(Cloth *cloth)
   return nondiag;
 }
 
+static float cloth_calc_volume(ClothModifierData *clmd)
+{
+  /* calc the (closed) cloth volume */
+  Cloth *cloth = clmd->clothObject;
+  const MVertTri *tri = cloth->tri;
+  Implicit_Data *data = cloth->implicit;
+  float vol = 0;
+
+  for (unsigned int i = 0; i < cloth->tri_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;
+}
+
 int BPH_cloth_solver_init(Object *UNUSED(ob), ClothModifierData *clmd)
 {
   Cloth *cloth = clmd->clothObject;
@@ -127,6 +146,13 @@ void BKE_cloth_solver_set_positions(ClothModifierData *clmd)
   }
 }
 
+void BKE_cloth_solver_set_volume(ClothModifierData *clmd)
+{
+  Cloth *cloth = clmd->clothObject;
+
+  cloth->initial_mesh_volume = cloth_calc_volume(clmd);
+}
+
 static bool collision_response(ClothModifierData *clmd,
                                CollisionModifierData *collmd,
                                CollPair *collpair,
@@ -526,6 +552,7 @@ static void cloth_calc_force(
 {
   /* Collect forces and derivatives:  F, dFdX, dFdV */
   Cloth *cloth = clmd->clothObject;
+  ClothSimSettings *parms = clmd->sim_parms;
   Implicit_Data *data = cloth->implicit;
   unsigned int i = 0;
   float drag = clmd->sim_parms->Cvi * 0.01f; /* viscosity of air scaled in percent */
@@ -570,6 +597,48 @@ static void cloth_calc_force(
 #ifdef CLOTH_FORCE_DRAG
   BPH_mass_spring_force_drag(data, drag);
 #endif
+  /* handle pressure forces */
+  if (parms->flags & CLOTH_SIMSETTINGS_FLAG_PRESSURE) {
+    /* The difference in pressure between the inside and outside of the mesh.*/
+    float pressure_difference = 0.0f;
+
+    float init_vol;
+    if (parms->flags & CLOTH_SIMSETTINGS_FLAG_PRESSURE_VOL) {
+      init_vol = clmd->sim_parms->target_volume;
+    }
+    else {
+      init_vol = cloth->initial_mesh_volume;
+    }
+
+    /* Check if we need to calculate the volume of the mesh. */
+    if (init_vol > 1E-6f) {
+      float f;
+      float vol = cloth_calc_volume(clmd);
+
+      /* Calculate an artifical 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;
+      }
+    }
+    pressure_difference += clmd->sim_parms->uniform_pressure_force;
+
+    pressure_difference *= clmd->sim_parms->pressure_factor;
+
+    for (i = 0; i < cloth->tri_num; i++) {
+      const MVertTri *vt = &tri[i];
+      if (fabs(pressure_difference) > 1E-6f) {
+        BPH_mass_spring_force_pressure(
+            data, vt->tri[0], vt->tri[1], vt->tri[2], pressure_difference);
+      }
+    }
+  }
 
   /* handle external forces like wind */
   if (effectors) {
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index 32416fa01ab..82a3f61e1e1 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -182,6 +182,11 @@ bool BPH_mass_spring_force_spring_goal(struct Implicit_Data *data,
                                        float stiffness,
                                        float damping);
 
+float BPH_tri_tetra_volume_signed_6x(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);
+
 /* ======== Hair Volumetric Forces ======== */
 
 struct HairGrid;
diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c
index d8b3f647591..fa093f482f7 100644
--- a/source/blender/physics/intern/implicit_blender.c
+++ b/source/blender/physics/intern/implicit_blender.c
@@ -1469,6 +1469,7 @@ void BPH_mass_spring_force_face_wind(
 
   /* 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 = effector_scale * area / 3.0f;
 
   world_to_root_v3(data, v1, win, winvec[v1]);
@@ -1481,6 +1482,29 @@ void BPH_mass_spring_force_face_wind(
   madd_v3_v3fl(data->F[v3], nor, factor * dot_v3v3(win, nor));
 }
 
+float BPH_tri_tetra_volume_signed_6x(Implicit_Data *data, int v1, int v2, int v3)
+{
+  /* The result will be 6x the volume */
+  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 nor[3], area;
+  float factor;
+
+  /* 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;
+
+  /* add pressure to each of the face verts */
+  madd_v3_v3fl(data->F[v1], nor, factor);
+  madd_v3_v3fl(data->F[v2], nor, factor);
+  madd_v3_v3fl(data->F[v3], nor, factor);
+}
+
 static void edge_wind_vertex(const float dir[3],
                              float length,
                              float radius,