Cloth: allow forces to act parallel to cloth.

Currently all force effectors can only act on cloth when the force is
perpendicular to the surface. This makes sense for wind, but not for
other forces; and the user may want even wind to have some friction.

This changes effector code to output two force vectors - although
you of course can pass the same pointer for both. The force is split
between the two outputs based on a new per-effector setting.

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

3 files changed, 84 insertions, 13 deletions

diff --git a/source/blender/physics/intern/BPH_mass_spring.cpp b/source/blender/physics/intern/BPH_mass_spring.cpp
index 18fab5215a6..3b6e39ef48f 100644
--- a/source/blender/physics/intern/BPH_mass_spring.cpp
+++ b/source/blender/physics/intern/BPH_mass_spring.cpp
@@ -683,23 +683,42 @@ static void cloth_calc_force(
 
   /* handle external forces like wind */
   if (effectors) {
+    bool is_not_hair = (clmd->hairdata == NULL) && (cloth->primitive_num > 0);
+    bool has_wind = false, has_force = false;
+
     /* cache per-vertex forces to avoid redundant calculation */
-    float(*winvec)[3] = (float(*)[3])MEM_callocN(sizeof(float[3]) * mvert_num, "effector forces");
+    float(*winvec)[3] = (float(*)[3])MEM_callocN(sizeof(float[3]) * mvert_num * 2,
+                                                 "effector forces");
+    float(*forcevec)[3] = is_not_hair ? winvec + mvert_num : winvec;
+
     for (i = 0; i < cloth->mvert_num; i++) {
       float x[3], v[3];
       EffectedPoint epoint;
 
       BPH_mass_spring_get_motion_state(data, i, x, v);
       pd_point_from_loc(scene, x, v, i, &epoint);
-      BKE_effectors_apply(
-          effectors, NULL, clmd->sim_parms->effector_weights, &epoint, winvec[i], NULL);
+      BKE_effectors_apply(effectors,
+                          NULL,
+                          clmd->sim_parms->effector_weights,
+                          &epoint,
+                          forcevec[i],
+                          winvec[i],
+                          NULL);
+
+      has_wind = has_wind || !is_zero_v3(winvec[i]);
+      has_force = has_force || !is_zero_v3(forcevec[i]);
     }
 
     /* Hair has only edges. */
-    if ((clmd->hairdata == NULL) && (cloth->primitive_num > 0)) {
+    if (is_not_hair) {
       for (i = 0; i < cloth->primitive_num; i++) {
         const MVertTri *vt = &tri[i];
-        BPH_mass_spring_force_face_wind(data, vt->tri[0], vt->tri[1], vt->tri[2], winvec);
+        if (has_wind) {
+          BPH_mass_spring_force_face_wind(data, vt->tri[0], vt->tri[1], vt->tri[2], winvec);
+        }
+        if (has_force) {
+          BPH_mass_spring_force_face_extern(data, vt->tri[0], vt->tri[1], vt->tri[2], forcevec);
+        }
       }
     }
     else {
diff --git a/source/blender/physics/intern/implicit.h b/source/blender/physics/intern/implicit.h
index 69b50f7fa48..8bc09755180 100644
--- a/source/blender/physics/intern/implicit.h
+++ b/source/blender/physics/intern/implicit.h
@@ -126,9 +126,12 @@ void BPH_mass_spring_force_drag(struct Implicit_Data *data, float drag);
 /* Custom external force */
 void BPH_mass_spring_force_extern(
     struct Implicit_Data *data, int i, const float f[3], float dfdx[3][3], float dfdv[3][3]);
-/* Wind force, acting on a face */
+/* Wind force, acting on a face (only generates pressure from the normal component) */
 void BPH_mass_spring_force_face_wind(
     struct Implicit_Data *data, int v1, int v2, int v3, const float (*winvec)[3]);
+/* Arbitrary per-unit-area vector force field acting on a face. */
+void BPH_mass_spring_force_face_extern(
+    struct Implicit_Data *data, int v1, int v2, int v3, const float (*forcevec)[3]);
 /* Wind force, acting on an edge */
 void BPH_mass_spring_force_edge_wind(struct Implicit_Data *data,
                                      int v1,
diff --git a/source/blender/physics/intern/implicit_blender.c b/source/blender/physics/intern/implicit_blender.c
index 063c224f158..5ec4c750d5d 100644
--- a/source/blender/physics/intern/implicit_blender.c
+++ b/source/blender/physics/intern/implicit_blender.c
@@ -1469,22 +1469,71 @@ void BPH_mass_spring_force_face_wind(
     Implicit_Data *data, int v1, int v2, int v3, const float (*winvec)[3])
 {
   const float effector_scale = 0.02f;
+  int vs[3] = {v1, v2, v3};
   float win[3], 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 = effector_scale * area / 3.0f;
 
-  world_to_root_v3(data, v1, win, winvec[v1]);
-  madd_v3_v3fl(data->F[v1], nor, factor * dot_v3v3(win, nor));
+  /* Calculate wind pressure at each vertex by projecting the wind field on the normal. */
+  for (int i = 0; i < 3; i++) {
+    world_to_root_v3(data, vs[i], win, winvec[vs[i]]);
 
-  world_to_root_v3(data, v2, win, winvec[v2]);
-  madd_v3_v3fl(data->F[v2], nor, factor * dot_v3v3(win, nor));
+    force[i] = dot_v3v3(win, nor);
+  }
 
-  world_to_root_v3(data, v3, win, winvec[v3]);
-  madd_v3_v3fl(data->F[v3], nor, factor * dot_v3v3(win, nor));
+  /* 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.
+   */
+  mul_v3_fl(force, factor / 4.0f);
+
+  base_force = force[0] + force[1] + force[2];
+
+  /* add pressure to each of the face verts */
+  madd_v3_v3fl(data->F[v1], nor, base_force + force[0]);
+  madd_v3_v3fl(data->F[v2], nor, base_force + force[1]);
+  madd_v3_v3fl(data->F[v3], nor, base_force + force[2]);
+}
+
+void BPH_mass_spring_force_face_extern(
+    Implicit_Data *data, int v1, int v2, int v3, const float (*forcevec)[3])
+{
+  const float effector_scale = 0.02f;
+  int vs[3] = {v1, v2, v3};
+  float nor[3], area;
+  float factor, base_force[3];
+  float force[3][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 = effector_scale * area / 3.0f;
+
+  /* Compute common and per-vertex force vectors from the original inputs. */
+  zero_v3(base_force);
+
+  for (int i = 0; i < 3; i++) {
+    world_to_root_v3(data, vs[i], force[i], forcevec[vs[i]]);
+
+    mul_v3_fl(force[i], factor / 4.0f);
+    add_v3_v3(base_force, force[i]);
+  }
+
+  /* Apply the common and vertex components to all vertices. */
+  for (int i = 0; i < 3; i++) {
+    add_v3_v3(force[i], base_force);
+    add_v3_v3(data->F[vs[i]], force[i]);
+  }
 }
 
 float BPH_tri_tetra_volume_signed_6x(Implicit_Data *data, int v1, int v2, int v3)