Cloth force fields fix.

Currently depsgraph code assumes that any dependency cycle is an
error, and it is reasonable to provide no behavior guarantees.

However, in the case of physics collision, cycles are actually
supported in a limited way by the underlying system via using
the saved collision data from the previous frame. It is also
reasonable for the user to want multiple cloth objects to collide
with one another, and it worked that way safely in 2.79.

In order to allow using this behavior, dependency graph should
know which relations are collisions, and prefer deleting them if
possible when a cycle is detected; it should also avoid excessive
message spam.

This adds a new relation flag that marks collision relations.
When resolving a cycle, collision relations going backward in
the alphabetic order of object names are chosen as primary
candidates for removal; this predictably resolves all cycles
between multiple objects all colliding with each other.

The depth first search cycle detection algorithm had to be
enchanced to unwind the stack to the relation being removed;
otherwise it cannot correctly detect all cycles. This would
already been a problem with the godmode flag, but it's not
that obvious there, as that flag is rare.

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

source/blender/depsgraph/DEG_depsgraph_build.h

@@ -182,7 +182,8 @@ void DEG_add_modifier_to_transform_relation(struct DepsNodeHandle *node_handle,
 void DEG_add_object_pointcache_relation(struct DepsNodeHandle *node_handle,
                                         struct Object *object,
                                         eDepsObjectComponentType component,
-                                        const char *description);
+                                        const char *description,
+                                        bool allow_break_cycles);
 
 void DEG_add_special_eval_flag(struct DepsNodeHandle *handle, struct ID *id, uint32_t flag);
 void DEG_add_customdata_mask(struct DepsNodeHandle *handle,

source/blender/depsgraph/intern/builder/deg_builder_cycle.cc

@@ -26,11 +26,13 @@
 // TOO(sergey): Use some wrappers over those?
 #include <cstdio>
 #include <cstdlib>
+#include <cstring>
 
 #include "BLI_utildefines.h"
 #include "BLI_stack.h"
 
 #include "intern/node/deg_node.h"
+#include "intern/node/deg_node_id.h"
 #include "intern/node/deg_node_component.h"
 #include "intern/node/deg_node_operation.h"
 
@@ -47,6 +49,8 @@ enum eCyclicCheckVisitedState {
   NODE_VISITED = 1,
   /* Node has been visited during traversal and is in current stack. */
   NODE_IN_STACK = 2,
+  /* Node processing has been started during traversal, but its incoming edge has been removed. */
+  NODE_STASHED = 3,
 };
 
 struct StackEntry {
@@ -59,6 +63,7 @@ struct CyclesSolverState {
   CyclesSolverState(Depsgraph *graph)
       : graph(graph),
         traversal_stack(BLI_stack_new(sizeof(StackEntry), "DEG detect cycles stack")),
+        stash_stack(BLI_stack_new(sizeof(OperationNode *), "DEG detect cycles stash stack")),
         num_cycles(0)
   {
     /* pass */
@@ -66,12 +71,14 @@ struct CyclesSolverState {
   ~CyclesSolverState()
   {
     BLI_stack_free(traversal_stack);
+    BLI_stack_free(stash_stack);
     if (num_cycles != 0) {
       printf("Detected %d dependency cycles\n", num_cycles);
     }
   }
   Depsgraph *graph;
   BLI_Stack *traversal_stack;
+  BLI_Stack *stash_stack;
   int num_cycles;
 };
 
@@ -130,6 +137,17 @@ void schedule_leaf_nodes(CyclesSolverState *state)
  */
 bool schedule_non_checked_node(CyclesSolverState *state)
 {
+  /* First, try unstashing nodes. */
+  while (!BLI_stack_is_empty(state->stash_stack)) {
+    OperationNode *node;
+    BLI_stack_pop(state->stash_stack, &node);
+
+    if (get_node_visited_state(node) == NODE_STASHED) {
+      schedule_node_to_stack(state, node);
+      return true;
+    }
+  }
+
   for (OperationNode *node : state->graph->operations) {
     if (get_node_visited_state(node) == NODE_NOT_VISITED) {
       schedule_node_to_stack(state, node);
@@ -139,6 +157,18 @@ bool schedule_non_checked_node(CyclesSolverState *state)
   return false;
 }
 
+int get_relation_name_order(Relation *relation)
+{
+  if (relation->from->type == NodeType::OPERATION && relation->to->type == NodeType::OPERATION) {
+    OperationNode *from = (OperationNode *)relation->from;
+    OperationNode *to = (OperationNode *)relation->to;
+
+    return strcmp(from->owner->owner->name.c_str() + 2, to->owner->owner->name.c_str() + 2);
+  }
+
+  return 0;
+}
+
 bool check_relation_can_murder(Relation *relation)
 {
   if (relation->flag & RELATION_FLAG_GODMODE) {
@@ -147,24 +177,58 @@ bool check_relation_can_murder(Relation *relation)
   return true;
 }
 
-Relation *select_relation_to_murder(Relation *relation, StackEntry *cycle_start_entry)
+bool check_relation_can_murder_physics(Relation *relation)
+{
+  /* In physics collision cycles, resolve cycles by sorting objects in their name
+   * order, and making relations in one direction first candidates for removal. */
+  return (relation->flag & RELATION_FLAG_PHYSICS_BREAK_CYCLES) &&
+         check_relation_can_murder(relation) && get_relation_name_order(relation) > 0;
+}
+
+bool check_relation_opposite_physics(Relation *relation)
+{
+  /* For detecting 'safe' cycles caused by two physics relations of opposite direction. */
+  return (relation->flag & RELATION_FLAG_PHYSICS_BREAK_CYCLES) &&
+         get_relation_name_order(relation) < 0;
+}
+
+Relation *select_relation(Relation *relation,
+                          StackEntry *cycle_start_entry,
+                          bool (*check)(Relation *))
 {
   /* More or less Russian roulette solver, which will make sure only
    * specially marked relations are kept alive.
    *
    * TODO(sergey): There might be better strategies here. */
-  if (check_relation_can_murder(relation)) {
+  if (check(relation)) {
     return relation;
   }
   StackEntry *current = cycle_start_entry;
   OperationNode *to_node = (OperationNode *)relation->to;
   while (current->node != to_node) {
-    if (check_relation_can_murder(current->via_relation)) {
+    if (check(current->via_relation)) {
       return current->via_relation;
     }
     current = current->from;
   }
-  return relation;
+  return NULL;
+}
+
+void report_cycle(
+    Relation *rel, StackEntry *entry, OperationNode *to, Relation *mark, const char *message)
+{
+  string cycle_str = "  " + to->full_identifier() + " depends on\n";
+  cycle_str += (rel == mark ? "* " : "  ") + entry->node->full_identifier() + " via '" +
+               rel->name + "'\n";
+  StackEntry *current = entry;
+  while (current->node != to) {
+    BLI_assert(current != NULL);
+    cycle_str += (current->via_relation == mark ? "* " : "  ") +
+                 current->from->node->full_identifier() + " via '" + current->via_relation->name +
+                 "'\n";
+    current = current->from;
+  }
+  printf("%s:\n%s", message, cycle_str.c_str());
 }
 
 /* Solve cycles with all nodes which are scheduled for traversal. */
@@ -178,25 +242,55 @@ void solve_cycles(CyclesSolverState *state)
     const int num_visited = get_node_num_visited_children(node);
     for (int i = num_visited; i < node->outlinks.size(); i++) {
       Relation *rel = node->outlinks[i];
+      /* Ignore already removed relations (can happen due to stashing). */
+      if (rel->flag & RELATION_FLAG_CYCLIC) {
+        continue;
+      }
       if (rel->to->type == NodeType::OPERATION) {
         OperationNode *to = (OperationNode *)rel->to;
         eCyclicCheckVisitedState to_state = get_node_visited_state(to);
         if (to_state == NODE_IN_STACK) {
-          string cycle_str = "  " + to->full_identifier() + " depends on\n  " +
-                             node->full_identifier() + " via '" + rel->name + "'\n";
-          StackEntry *current = entry;
-          while (current->node != to) {
-            BLI_assert(current != NULL);
-            cycle_str += "  " + current->from->node->full_identifier() + " via '" +
-                         current->via_relation->name + "'\n";
-            current = current->from;
+          /* Try removing a physics (collision) relation. */
+          Relation *sacrificial_relation = select_relation(
+              rel, entry, check_relation_can_murder_physics);
+          if (sacrificial_relation) {
+            /* Reduce logging for safe cycles caused by opposite physics dependencies. */
+            if (select_relation(rel, entry, check_relation_opposite_physics)) {
+              printf("Ignoring physics dependency of %s on %s (%s)\n",
+                     ((OperationNode *)sacrificial_relation->to)->owner->owner->name.c_str() + 2,
+                     ((OperationNode *)sacrificial_relation->from)->owner->owner->name.c_str() + 2,
+                     sacrificial_relation->name);
+            }
+            else {
+              report_cycle(
+                  rel, entry, to, sacrificial_relation, "Physics dependency cycle detected");
+            }
+          }
+          /* Remove a random relation. */
+          else {
+            sacrificial_relation = select_relation(rel, entry, check_relation_can_murder);
+            if (!sacrificial_relation) {
+              sacrificial_relation = rel;
+            }
+            report_cycle(rel, entry, to, sacrificial_relation, "Dependency cycle detected");
           }
-          printf("Dependency cycle detected:\n%s", cycle_str.c_str());
-          Relation *sacrificial_relation = select_relation_to_murder(rel, entry);
           sacrificial_relation->flag |= RELATION_FLAG_CYCLIC;
           ++state->num_cycles;
+          /* If removing a relation deep in the stack, it is necessary to unwind. */
+          if (sacrificial_relation != rel) {
+            StackEntry cur_entry;
+            set_node_num_visited_children(node, i);
+            do {
+              BLI_stack_pop(traversal_stack, &cur_entry);
+              set_node_visited_state(cur_entry.node, NODE_STASHED);
+              BLI_stack_push(state->stash_stack, &cur_entry.node);
+            } while (cur_entry.via_relation != sacrificial_relation &&
+                     !BLI_stack_is_empty(traversal_stack));
+            all_child_traversed = false;
+            break;
+          }
         }
-        else if (to_state == NODE_NOT_VISITED) {
+        else if (to_state == NODE_NOT_VISITED || to_state == NODE_STASHED) {
           StackEntry new_entry;
           new_entry.node = to;
           new_entry.from = entry;

source/blender/depsgraph/intern/depsgraph.h

@@ -73,6 +73,8 @@ enum RelationFlag {
   RELATION_FLAG_GODMODE = (1 << 4),
   /* Relation will check existence before being added. */
   RELATION_CHECK_BEFORE_ADD = (1 << 5),
+  /* Physics relation that should be preferred for breaking cycles. */
+  RELATION_FLAG_PHYSICS_BREAK_CYCLES = (1 << 6),
 };
 
 /* B depends on A (A -> B) */

source/blender/depsgraph/intern/depsgraph_build.cc

@@ -131,18 +131,21 @@ void DEG_add_bone_relation(DepsNodeHandle *node_handle,
 void DEG_add_object_pointcache_relation(struct DepsNodeHandle *node_handle,
                                         struct Object *object,
                                         eDepsObjectComponentType component,
-                                        const char *description)
+                                        const char *description,
+                                        bool allow_break_cycles)
 {
   DEG::NodeType type = DEG::nodeTypeFromObjectComponent(component);
   DEG::ComponentKey comp_key(&object->id, type);
   DEG::DepsNodeHandle *deg_node_handle = get_node_handle(node_handle);
   DEG::DepsgraphRelationBuilder *relation_builder = deg_node_handle->builder;
+  int flags = allow_break_cycles ? DEG::RELATION_FLAG_PHYSICS_BREAK_CYCLES : 0;
   /* Add relation from source to the node handle. */
-  relation_builder->add_node_handle_relation(comp_key, deg_node_handle, description);
+  relation_builder->add_node_handle_relation(comp_key, deg_node_handle, description, flags);
   /* Node deduct point cache component and connect source to it. */
   ID *id = DEG_get_id_from_handle(node_handle);
   DEG::ComponentKey point_cache_key(id, DEG::NodeType::POINT_CACHE);
-  DEG::Relation *rel = relation_builder->add_relation(comp_key, point_cache_key, "Point Cache");
+  DEG::Relation *rel = relation_builder->add_relation(
+      comp_key, point_cache_key, "Point Cache", flags);
   if (rel != NULL) {
     rel->flag |= DEG::RELATION_FLAG_FLUSH_USER_EDIT_ONLY;
   }

source/blender/depsgraph/intern/depsgraph_physics.cc

@@ -108,8 +108,8 @@ void DEG_add_collision_relations(DepsNodeHandle *handle,
     }
     if (filter_function == NULL ||
         filter_function(ob1, modifiers_findByType(ob1, (ModifierType)modifier_type))) {
-      DEG_add_object_pointcache_relation(handle, ob1, DEG_OB_COMP_TRANSFORM, name);
-      DEG_add_object_pointcache_relation(handle, ob1, DEG_OB_COMP_GEOMETRY, name);
+      DEG_add_object_pointcache_relation(handle, ob1, DEG_OB_COMP_TRANSFORM, name, false);
+      DEG_add_object_pointcache_relation(handle, ob1, DEG_OB_COMP_GEOMETRY, name, true);
     }
   }
 }
@@ -134,21 +134,21 @@ void DEG_add_forcefield_relations(DepsNodeHandle *handle,
 
     /* Relation to forcefield object, optionally including geometry.
      * Use special point cache relations for automatic cache clearing. */
-    DEG_add_object_pointcache_relation(handle, relation->ob, DEG_OB_COMP_TRANSFORM, name);
+    DEG_add_object_pointcache_relation(handle, relation->ob, DEG_OB_COMP_TRANSFORM, name, false);
 
     if (relation->psys || ELEM(relation->pd->shape, PFIELD_SHAPE_SURFACE, PFIELD_SHAPE_POINTS) ||
         relation->pd->forcefield == PFIELD_GUIDE) {
       /* TODO(sergey): Consider going more granular with more dedicated
        * particle system operation. */
-      DEG_add_object_pointcache_relation(handle, relation->ob, DEG_OB_COMP_GEOMETRY, name);
+      DEG_add_object_pointcache_relation(handle, relation->ob, DEG_OB_COMP_GEOMETRY, name, false);
     }
 
     /* Smoke flow relations. */
     if (relation->pd->forcefield == PFIELD_SMOKEFLOW && relation->pd->f_source != NULL) {
       DEG_add_object_pointcache_relation(
-          handle, relation->pd->f_source, DEG_OB_COMP_TRANSFORM, "Smoke Force Domain");
+          handle, relation->pd->f_source, DEG_OB_COMP_TRANSFORM, "Smoke Force Domain", false);
       DEG_add_object_pointcache_relation(
-          handle, relation->pd->f_source, DEG_OB_COMP_GEOMETRY, "Smoke Force Domain");
+          handle, relation->pd->f_source, DEG_OB_COMP_GEOMETRY, "Smoke Force Domain", false);
     }
 
     /* Absorption forces need collision relation. */

source/blender/physics/intern/BPH_mass_spring.cpp

@@ -746,7 +746,7 @@ static void cloth_calc_force(
 
     for (i = 0; i < cloth->tri_num; i++) {
       const MVertTri *vt = &tri[i];
-      BPH_mass_spring_force_face_wind(data, vt->tri[0], vt->tri[1], vt->tri[2], winvec);
+      BPH_mass_spring_force_face_extern(data, vt->tri[0], vt->tri[1], vt->tri[2], winvec, 0.02f);
     }
 
     /* Hair has only edges */

source/blender/physics/intern/implicit.h

@@ -125,9 +125,13 @@ 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]);
+/* Custom force, acting on the surface of the face (i.e. specified per unit of area) */
+void BPH_mass_spring_force_face_extern(
+    struct Implicit_Data *data, int v1, int v2, int v3, const float (*forcevec)[3], float scale);
 /* Wind force, acting on a face */
 void BPH_mass_spring_force_face_wind(
     struct Implicit_Data *data, int v1, int v2, int v3, const float (*winvec)[3]);
+
 /* Wind force, acting on an edge */
 void BPH_mass_spring_force_edge_wind(struct Implicit_Data *data,
                                      int v1,

source/blender/physics/intern/implicit_blender.c

@@ -1458,6 +1458,40 @@ static float calc_nor_area_tri(float nor[3],
   return normalize_v3(nor) / 2.0f;
 }
 
+/* XXX does not support force jacobians yet, since the effector system does not provide them either
+ */
+void BPH_mass_spring_force_face_extern(
+    Implicit_Data *data, int v1, int v2, int v3, const float (*forcevec)[3], float scale)
+{
+  float force_sum[3], pt_force[3], root_force[3], nor[3], area;
+  int v[3] = {v1, v2, v3};
+  float factor;
+
+  /* Calculate face normal and area. */
+  area = calc_nor_area_tri(nor, data->X[v1], data->X[v2], data->X[v3]);
+
+  /* Integrating a force specified per unit of area over the surface of the
+   * triangle (effectively a pressure vector) yields a formula for forces at
+   * the vertices:
+   *
+   * force[i] = (sum(pressure) + pressure[i]) * area / 12
+   */
+  factor = scale * area / 12.0f;
+
+  zero_v3(force_sum);
+
+  for (int i = 0; i < 3; i++) {
+    add_v3_v3(force_sum, forcevec[v[i]]);
+  }
+
+  for (int i = 0; i < 3; i++) {
+    int vi = v[i];
+    add_v3_v3v3(pt_force, force_sum, forcevec[vi]);
+    world_to_root_v3(data, vi, root_force, pt_force);
+    madd_v3_v3fl(data->F[vi], root_force, factor);
+  }
+}
+
 /* XXX does not support force jacobians yet, since the effector system does not provide them either
  */
 void BPH_mass_spring_force_face_wind(