Adjusts EdgexEdge

# Conflicts:
#	source/blender/makesdna/DNA_cloth_types.h

release/scripts/startup/bl_ui/properties_physics_cloth.py

@@ -351,12 +351,8 @@ class PHYSICS_PT_cloth_object_collision(PhysicButtonsPanel, Panel):
         flow = layout.grid_flow(row_major=False, columns=0, even_columns=True, even_rows=False, align=True)
 
         col = flow.column()
+        col.prop(cloth, "primitive")
         col.prop(cloth, "distance_min", slider=True, text="Distance")
-
-        col = flow.column()
-        col.prop(cloth, "impulse_clamp")
-
-        col = flow.column()
         col.prop(cloth, "collection")
 
 
@@ -384,14 +380,10 @@ class PHYSICS_PT_cloth_self_collision(PhysicButtonsPanel, Panel):
         flow = layout.grid_flow(row_major=False, columns=0, even_columns=True, even_rows=False, align=True)
 
         col = flow.column()
+        col.prop(cloth, "self_primitive")
         col.prop(cloth, "self_friction", text="Friction")
-
-        col = flow.column()
         col.prop(cloth, "self_distance_min", slider=True, text="Distance")
 
-        col = flow.column()
-        col.prop(cloth, "self_impulse_clamp")
-
         col = flow.column()
         col.prop_search(cloth, "vertex_group_self_collisions", ob, "vertex_groups", text="Vertex Group")
 

release/scripts/startup/bl_ui/properties_physics_field.py

@@ -418,9 +418,6 @@ class PHYSICS_PT_collision_softbody(PhysicButtonsPanel, Panel):
         col = flow.column()
         col.prop(settings, "cloth_friction")
 
-        col = flow.column()
-        col.prop(settings, "use_culling")
-
         col = flow.column()
         col.prop(settings, "use_normal")
 

source/blender/blenkernel/BKE_cloth.h

@@ -42,10 +42,6 @@ struct Scene;
 /* goal defines */
 #define SOFTGOALSNAP 0.999f
 
-/* This is approximately the smallest number that can be
- * represented by a float, given its precision. */
-#define ALMOST_ZERO FLT_EPSILON
-
 /* Bits to or into the ClothVertex.flags. */
 typedef enum eClothVertexFlag {
   CLOTH_VERT_FLAG_PINNED = (1 << 0),
@@ -78,45 +74,38 @@ typedef struct ClothSolverResult {
  * own connectivity of the mesh based on the actual edges in the mesh.
  */
 typedef struct Cloth {
-  struct ClothVertex *verts;     /* The vertices that represent this cloth. */
-  struct LinkNode *springs;      /* The springs connecting the mesh. */
-  unsigned int numsprings;       /* The count of springs. */
-  unsigned int mvert_num;        /* The number of verts == m * n. */
-  unsigned int primitive_num;    /* Number of triangles for cloth and edges for hair. */
-  unsigned char old_solver_type; /* unused, only 1 solver here */
-  unsigned char pad2;
-  short pad3;
-  struct BVHTree *bvhtree;     /* collision tree for this cloth object */
-  struct BVHTree *bvhselftree; /* collision tree for this cloth object */
-  struct MVertTri *tri;
   struct Implicit_Data *implicit; /* our implicit solver connects to this pointer */
-  struct EdgeSet *edgeset;        /* used for selfcollisions */
+  struct ClothVertex *verts;      /* The vertices that represent this cloth. */
-  int last_frame;
+  struct LinkNode *springs;       /* The springs connecting the mesh. */
-  float initial_mesh_volume;      /* Initial volume of the mesh. Used for pressure */
+  struct MVertTri *tri;
-  float average_acceleration[3];  /* Moving average of overall acceleration. */
   struct MEdge *edges;            /* Used for hair collisions. */
   struct EdgeSet *sew_edge_graph; /* Sewing edges represented using a GHash */
+  struct ClothCollisionData *collision_data;
+  unsigned int numsprings; /* The count of springs. */
+  unsigned int mvert_num;  /* The number of verts == m * n. */
+  unsigned int medge_num;
+  unsigned int mlooptri_num;
+  int last_frame;
+  float initial_mesh_volume;     /* Initial volume of the mesh. Used for pressure */
+  float average_acceleration[3]; /* Moving average of overall acceleration. */
 } Cloth;
 
 /**
  * The definition of a cloth vertex.
  */
 typedef struct ClothVertex {
-  int flags;                  /* General flags per vertex.        */
+  int flags;            /* General flags per vertex.        */
-  float v[3];                 /* The velocity of the point.       */
+  float v[3];           /* The velocity of the point.       */
-  float xconst[3];            /* constrained position         */
+  float xconst[3];      /* constrained position         */
-  float x[3];                 /* The current position of this vertex. */
+  float x[3];           /* The current position of this vertex. */
-  float xold[3];              /* The previous position of this vertex.*/
+  float xold[3];        /* The previous position of this vertex.*/
-  float tx[3];                /* temporary position */
+  float tx[3];          /* temporary position */
-  float txold[3];             /* temporary old position */
+  float txold[3];       /* temporary old position */
-  float tv[3];                /* temporary "velocity", mostly used as tv = tx-txold */
+  float tv[3];          /* temporary "velocity", mostly used as tv = tx-txold */
-  float mass;                 /* mass / weight of the vertex      */
+  float mass;           /* mass / weight of the vertex      */
-  float goal;                 /* goal, from SB            */
+  float goal;           /* goal, from SB            */
-  float impulse[3];           /* used in collision.c */
+  float xrest[3];       /* rest position of the vertex */
-  float xrest[3];             /* rest position of the vertex */
+  float avg_spring_len; /* average length of connected springs */
-  float dcvel[3];             /* delta velocities to be applied by collision response */
-  unsigned int impulse_count; /* same as above */
-  float avg_spring_len;       /* average length of connected springs */
   float struct_stiff;
   float bend_stiff;
   float shear_stiff;
@@ -164,20 +153,6 @@ typedef struct ClothSpring {
     *(v1 + 2) = *(v2 + 2) * aS + *(v3 + 2) * bS; \
   } \
   ((void)0)
-#define VECADDS(v1, v2, v3, bS) \
-  { \
-    *(v1) = *(v2) + *(v3)*bS; \
-    *(v1 + 1) = *(v2 + 1) + *(v3 + 1) * bS; \
-    *(v1 + 2) = *(v2 + 2) + *(v3 + 2) * bS; \
-  } \
-  ((void)0)
-#define VECSUBMUL(v1, v2, aS) \
-  { \
-    *(v1) -= *(v2)*aS; \
-    *(v1 + 1) -= *(v2 + 1) * aS; \
-    *(v1 + 2) -= *(v2 + 2) * aS; \
-  } \
-  ((void)0)
 #define VECSUBS(v1, v2, v3, bS) \
   { \
     *(v1) = *(v2) - *(v3)*bS; \
@@ -224,13 +199,6 @@ typedef struct ColliderContacts {
   int totcollisions;
 } ColliderContacts;
 
-// needed for implicit.c
-int cloth_bvh_collision(struct Depsgraph *depsgraph,
-                        struct Object *ob,
-                        struct ClothModifierData *clmd,
-                        float step,
-                        float dt);
-
 ////////////////////////////////////////////////
 
 /////////////////////////////////////////////////
@@ -249,9 +217,6 @@ void clothModifier_do(struct ClothModifierData *clmd,
 
 int cloth_uses_vgroup(struct ClothModifierData *clmd);
 
-// needed for collision.c
-void bvhtree_update_from_cloth(struct ClothModifierData *clmd, bool moving, bool self);
-
 // needed for button_object.c
 void cloth_clear_cache(struct Object *ob, struct ClothModifierData *clmd, float framenr);
 

source/blender/blenkernel/BKE_cloth_collision.h

@@ -0,0 +1,55 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) Blender Foundation.
+ * All rights reserved.
+ */
+#pragma once
+
+/** \file
+ * \ingroup bke
+ */
+
+#include <float.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+/* types */
+#include "BKE_collision.h"
+#include "DNA_cloth_types.h"
+
+#include "BLI_kdopbvh.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct ClothModifierData;
+struct Depsgraph;
+struct Object;
+
+void BKE_cloth_collision_init(struct ClothModifierData *clmd);
+void BKE_cloth_collision_free(struct ClothModifierData *clmd);
+
+int BKE_cloth_collision_do_step(struct Depsgraph *depsgraph,
+                                struct Object *ob,
+                                struct ClothModifierData *clmd,
+                                float step,
+                                float dt);
+
+#ifdef __cplusplus
+}
+#endif

source/blender/blenkernel/BKE_collision.h

@@ -44,67 +44,6 @@ struct MVert;
 struct MVertTri;
 struct Object;
 
-////////////////////////////////////////
-// used for collisions in collision.c
-////////////////////////////////////////
-
-/* COLLISION FLAGS */
-typedef enum {
-  COLLISION_IN_FUTURE = (1 << 1),
-#ifdef WITH_ELTOPO
-  COLLISION_USE_COLLFACE = (1 << 2),
-  COLLISION_IS_EDGES = (1 << 3),
-#endif
-  COLLISION_INACTIVE = (1 << 4),
-} COLLISION_FLAGS;
-
-////////////////////////////////////////
-// used for collisions in collision.c
-////////////////////////////////////////
-/* used for collisions in collision.c */
-typedef struct CollPair {
-  unsigned int face1; /* cloth face */
-  unsigned int face2; /* object face */
-  float distance;
-  float normal[3];
-  float vector[3];    /* unnormalized collision vector: p2-p1 */
-  float pa[3], pb[3]; /* collision point p1 on face1, p2 on face2 */
-  int flag;
-  float time; /* collision time, from 0 up to 1 */
-
-  /* mesh-mesh collision */
-#ifdef WITH_ELTOPO /*either ap* or bp* can be set, but not both*/
-  float bary[3];
-  int ap1, ap2, ap3, collp, bp1, bp2, bp3;
-  int collface;
-#else
-  int ap1, ap2, ap3, bp1, bp2, bp3;
-#endif
-  int pointsb[4];
-} CollPair;
-
-/* used for collisions in collision.c */
-typedef struct EdgeCollPair {
-  unsigned int p11, p12, p21, p22;
-  float normal[3];
-  float vector[3];
-  float time;
-  int lastsign;
-  float pa[3], pb[3]; /* collision point p1 on face1, p2 on face2 */
-} EdgeCollPair;
-
-/* used for collisions in collision.c */
-typedef struct FaceCollPair {
-  unsigned int p11, p12, p13, p21;
-  float normal[3];
-  float vector[3];
-  float time;
-  int lastsign;
-  float pa[3], pb[3]; /* collision point p1 on face1, p2 on face2 */
-} FaceCollPair;
-
-////////////////////////////////////////
-
 /////////////////////////////////////////////////
 // forward declarations
 /////////////////////////////////////////////////
@@ -113,30 +52,24 @@ typedef struct FaceCollPair {
 // used in modifier.c from collision.c
 /////////////////////////////////////////////////
 
-BVHTree *bvhtree_build_from_mvert(const struct MVert *mvert,
+BVHTree *BKE_collision_bvhtree_tris_build(const struct MVert *mvert,
-                                  const struct MVertTri *tri,
+                                          const struct MVertTri *tri,
-                                  int tri_num,
+                                          int tri_num,
-                                  float epsilon);
+                                          float epsilon);
-void bvhtree_update_from_mvert(BVHTree *bvhtree,
+void BKE_collision_bvhtree_ensure_and_update(struct CollisionModifierData *collmd,
-                               const struct MVert *mvert,
+                                             const bool moving,
-                               const struct MVert *mvert_moving,
+                                             const bool update_bvhtree_tris,
-                               const struct MVertTri *tri,
+                                             const bool update_bvhtree_edges,
-                               int tri_num,
+                                             const bool update_bvhtree_verts);
-                               bool moving);
+void BKE_collision_bvhtree_free_unused(struct Depsgraph *depsgraph, struct Collection *collection);
 
 /////////////////////////////////////////////////
 
 /* move Collision modifier object inter-frame with step = [0,1]
  * defined in collisions.c */
-void collision_move_object(struct CollisionModifierData *collmd,
+void BKE_collision_move_object(struct CollisionModifierData *collmd,
-                           const float step,
+                               const float step,
-                           const float prevstep,
+                               const float prevstep);
-                           const bool moving_bvh);
-
-void collision_get_collider_velocity(float vel_old[3],
-                                     float vel_new[3],
-                                     struct CollisionModifierData *collmd,
-                                     struct CollPair *collpair);
 
 /* Collision relations for dependency graph build. */
 

source/blender/blenkernel/CMakeLists.txt

@@ -93,6 +93,7 @@ set(SRC
   intern/camera.c
   intern/cdderivedmesh.c
   intern/cloth.c
+  intern/cloth_collision.c
   intern/collection.c
   intern/collision.c
   intern/colorband.c
@@ -283,6 +284,7 @@ set(SRC
   BKE_ccg.h
   BKE_cdderivedmesh.h
   BKE_cloth.h
+  BKE_cloth_collision.h
   BKE_collection.h
   BKE_collision.h
   BKE_colorband.h

source/blender/blenkernel/intern/boids.c

@@ -256,8 +256,8 @@ static bool rule_avoid_collision(BoidRule *rule,
       col.current = coll->ob;
       col.md = coll->collmd;
 
-      if (col.md && col.md->bvhtree) {
+      if (col.md && col.md->bvhtree_tris) {
-        BLI_bvhtree_ray_cast_ex(col.md->bvhtree,
+        BLI_bvhtree_ray_cast_ex(col.md->bvhtree_tris,
                                 col.co1,
                                 ray_dir,
                                 radius,
@@ -901,8 +901,8 @@ static Object *boid_find_ground(BoidBrainData *bbd,
     col.md = coll->collmd;
     col.fac1 = col.fac2 = 0.0f;
 
-    if (col.md && col.md->bvhtree) {
+    if (col.md && col.md->bvhtree_tris) {
-      BLI_bvhtree_ray_cast_ex(col.md->bvhtree,
+      BLI_bvhtree_ray_cast_ex(col.md->bvhtree_tris,
                               col.co1,
                               ray_dir,
                               radius,
@@ -933,8 +933,8 @@ static Object *boid_find_ground(BoidBrainData *bbd,
     col.current = coll->ob;
     col.md = coll->collmd;
 
-    if (col.md && col.md->bvhtree) {
+    if (col.md && col.md->bvhtree_tris) {
-      BLI_bvhtree_ray_cast_ex(col.md->bvhtree,
+      BLI_bvhtree_ray_cast_ex(col.md->bvhtree_tris,
                               col.co1,
                               ray_dir,
                               radius,

source/blender/blenkernel/intern/cloth.c

@@ -40,6 +40,7 @@
 
 #include "BKE_bvhutils.h"
 #include "BKE_cloth.h"
+#include "BKE_cloth_collision.h"
 #include "BKE_effect.h"
 #include "BKE_global.h"
 #include "BKE_mesh.h"
@@ -76,140 +77,6 @@ typedef struct BendSpringRef {
  *
  ******************************************************************************/
 
-static BVHTree *bvhtree_build_from_cloth(ClothModifierData *clmd, float epsilon)
-{
-  if (!clmd) {
-    return NULL;
-  }
-
-  Cloth *cloth = clmd->clothObject;
-
-  if (!cloth) {
-    return NULL;
-  }
-
-  ClothVertex *verts = cloth->verts;
-  const MVertTri *vt = cloth->tri;
-
-  /* in the moment, return zero if no faces there */
-  if (!cloth->primitive_num) {
-    return NULL;
-  }
-
-  /* create quadtree with k=26 */
-  BVHTree *bvhtree = BLI_bvhtree_new(cloth->primitive_num, epsilon, 4, 26);
-
-  /* fill tree */
-  if (clmd->hairdata == NULL) {
-    for (int i = 0; i < cloth->primitive_num; i++, vt++) {
-      float co[3][3];
-
-      copy_v3_v3(co[0], verts[vt->tri[0]].xold);
-      copy_v3_v3(co[1], verts[vt->tri[1]].xold);
-      copy_v3_v3(co[2], verts[vt->tri[2]].xold);
-
-      BLI_bvhtree_insert(bvhtree, i, co[0], 3);
-    }
-  }
-  else {
-    MEdge *edges = cloth->edges;
-
-    for (int i = 0; i < cloth->primitive_num; i++) {
-      float co[2][3];
-
-      copy_v3_v3(co[0], verts[edges[i].v1].xold);
-      copy_v3_v3(co[1], verts[edges[i].v2].xold);
-
-      BLI_bvhtree_insert(bvhtree, i, co[0], 2);
-    }
-  }
-
-  /* balance tree */
-  BLI_bvhtree_balance(bvhtree);
-
-  return bvhtree;
-}
-
-void bvhtree_update_from_cloth(ClothModifierData *clmd, bool moving, bool self)
-{
-  unsigned int i = 0;
-  Cloth *cloth = clmd->clothObject;
-  BVHTree *bvhtree;
-  ClothVertex *verts = cloth->verts;
-  const MVertTri *vt;
-
-  BLI_assert(!(clmd->hairdata != NULL && self));
-
-  if (self) {
-    bvhtree = cloth->bvhselftree;
-  }
-  else {
-    bvhtree = cloth->bvhtree;
-  }
-
-  if (!bvhtree) {
-    return;
-  }
-
-  vt = cloth->tri;
-
-  /* update vertex position in bvh tree */
-  if (clmd->hairdata == NULL) {
-    if (verts && vt) {
-      for (i = 0; i < cloth->primitive_num; i++, vt++) {
-        float co[3][3], co_moving[3][3];
-        bool ret;
-
-        /* copy new locations into array */
-        if (moving) {
-          copy_v3_v3(co[0], verts[vt->tri[0]].txold);
-          copy_v3_v3(co[1], verts[vt->tri[1]].txold);
-          copy_v3_v3(co[2], verts[vt->tri[2]].txold);
-
-          /* update moving positions */
-          copy_v3_v3(co_moving[0], verts[vt->tri[0]].tx);
-          copy_v3_v3(co_moving[1], verts[vt->tri[1]].tx);
-          copy_v3_v3(co_moving[2], verts[vt->tri[2]].tx);
-
-          ret = BLI_bvhtree_update_node(bvhtree, i, co[0], co_moving[0], 3);
-        }
-        else {
-          copy_v3_v3(co[0], verts[vt->tri[0]].tx);
-          copy_v3_v3(co[1], verts[vt->tri[1]].tx);
-          copy_v3_v3(co[2], verts[vt->tri[2]].tx);
-
-          ret = BLI_bvhtree_update_node(bvhtree, i, co[0], NULL, 3);
-        }
-
-        /* check if tree is already full */
-        if (ret == false) {
-          break;
-        }
-      }
-
-      BLI_bvhtree_update_tree(bvhtree);
-    }
-  }
-  else {
-    if (verts) {
-      MEdge *edges = cloth->edges;
-
-      for (i = 0; i < cloth->primitive_num; i++) {
-        float co[2][3];
-
-        copy_v3_v3(co[0], verts[edges[i].v1].tx);
-        copy_v3_v3(co[1], verts[edges[i].v2].tx);
-
-        if (!BLI_bvhtree_update_node(bvhtree, i, co[0], NULL, 2)) {
-          break;
-        }
-      }
-
-      BLI_bvhtree_update_tree(bvhtree);
-    }
-  }
-}
-
 void cloth_clear_cache(Object *ob, ClothModifierData *clmd, float framenr)
 {
   PTCacheID pid;
@@ -253,7 +120,6 @@ static bool do_init_cloth(Object *ob, ClothModifierData *clmd, Mesh *result, int
     }
 
     clmd->clothObject->last_frame = MINFRAME - 1;
-    clmd->sim_parms->dt = 1.0f / clmd->sim_parms->stepsPerFrame;
   }
 
   return true;
@@ -467,22 +333,15 @@ void cloth_free_modifier(ClothModifierData *clmd)
     cloth->springs = NULL;
     cloth->numsprings = 0;
 
-    /* free BVH collision tree */
+    BKE_cloth_collision_free(clmd);
-    if (cloth->bvhtree) {
-      BLI_bvhtree_free(cloth->bvhtree);
-    }
-
-    if (cloth->bvhselftree) {
-      BLI_bvhtree_free(cloth->bvhselftree);
-    }
 
     /* we save our faces for collision objects */
     if (cloth->tri) {
       MEM_freeN(cloth->tri);
     }
 
-    if (cloth->edgeset) {
+    if (cloth->edges) {
-      BLI_edgeset_free(cloth->edgeset);
+      MEM_freeN(cloth->edges);
     }
 
     if (cloth->sew_edge_graph) {
@@ -550,21 +409,15 @@ void cloth_free_modifier_extern(ClothModifierData *clmd)
     cloth->numsprings = 0;
 
     /* free BVH collision tree */
-    if (cloth->bvhtree) {
+    BKE_cloth_collision_free(clmd);
-      BLI_bvhtree_free(cloth->bvhtree);
-    }
-
-    if (cloth->bvhselftree) {
-      BLI_bvhtree_free(cloth->bvhselftree);
-    }
 
     /* we save our faces for collision objects */
     if (cloth->tri) {
       MEM_freeN(cloth->tri);
     }
 
-    if (cloth->edgeset) {
+    if (cloth->edges) {
-      BLI_edgeset_free(cloth->edgeset);
+      MEM_freeN(cloth->edges);
     }
 
     if (cloth->sew_edge_graph) {
@@ -725,7 +578,6 @@ static bool cloth_from_object(
   MVert *mvert = NULL;
   ClothVertex *verts = NULL;
   float(*shapekey_rest)[3] = NULL;
-  const float tnull[3] = {0, 0, 0};
 
   /* If we have a clothObject, free it. */
   if (clmd->clothObject != NULL) {
@@ -737,11 +589,7 @@ static bool cloth_from_object(
 
   /* Allocate a new cloth object. */
   clmd->clothObject = MEM_callocN(sizeof(Cloth), "cloth");
-  if (clmd->clothObject) {
+  if (clmd->clothObject == NULL) {
-    clmd->clothObject->old_solver_type = 255;
-    clmd->clothObject->edgeset = NULL;
-  }
-  else {
     BKE_modifier_set_error(ob, &(clmd->modifier), "Out of memory on allocating clmd->clothObject");
     return false;
   }
@@ -786,7 +634,6 @@ static bool cloth_from_object(
 
     /* no GUI interface yet */
     verts->mass = clmd->sim_parms->mass;
-    verts->impulse_count = 0;
 
     if (clmd->sim_parms->vgroup_mass > 0) {
       verts->goal = clmd->sim_parms->defgoal;
@@ -803,9 +650,6 @@ static bool cloth_from_object(
     copy_v3_v3(verts->txold, verts->x);
     copy_v3_v3(verts->tx, verts->x);
     mul_v3_fl(verts->v, 0.0f);
-
-    verts->impulse_count = 0;
-    copy_v3_v3(verts->impulse, tnull);
   }
 
   /* apply / set vertex groups */
@@ -825,8 +669,7 @@ static bool cloth_from_object(
     SIM_cloth_solver_set_positions(clmd);
   }
 
-  clmd->clothObject->bvhtree = bvhtree_build_from_cloth(clmd, clmd->coll_parms->epsilon);
+  BKE_cloth_collision_init(clmd);
-  clmd->clothObject->bvhselftree = bvhtree_build_from_cloth(clmd, clmd->coll_parms->selfepsilon);
 
   return true;
 }
@@ -851,12 +694,9 @@ static void cloth_from_mesh(ClothModifierData *clmd, const Object *ob, Mesh *mes
   }
 
   /* save face information */
-  if (clmd->hairdata == NULL) {
+  clmd->clothObject->mvert_num = mvert_num;
-    clmd->clothObject->primitive_num = looptri_num;
+  clmd->clothObject->medge_num = mesh->totedge;
-  }
+  clmd->clothObject->mlooptri_num = looptri_num;
-  else {
-    clmd->clothObject->primitive_num = mesh->totedge;
-  }
 
   clmd->clothObject->tri = MEM_mallocN(sizeof(MVertTri) * looptri_num, "clothLoopTris");
   if (clmd->clothObject->tri == NULL) {
@@ -868,7 +708,7 @@ static void cloth_from_mesh(ClothModifierData *clmd, const Object *ob, Mesh *mes
   }
   BKE_mesh_runtime_verttri_from_looptri(clmd->clothObject->tri, mloop, looptri, looptri_num);
 
-  clmd->clothObject->edges = mesh->medge;
+  clmd->clothObject->edges = MEM_dupallocN(mesh->medge);
 
   /* Free the springs since they can't be correct if the vertices
    * changed.
@@ -928,11 +768,6 @@ static void cloth_free_errorsprings(Cloth *cloth,
   cloth_free_edgelist(edgelist, cloth->mvert_num);
 
   MEM_SAFE_FREE(spring_ref);
-
-  if (cloth->edgeset) {
-    BLI_edgeset_free(cloth->edgeset);
-    cloth->edgeset = NULL;
-  }
 }
 
 BLI_INLINE void cloth_bend_poly_dir(
@@ -1141,7 +976,7 @@ static void cloth_update_springs(ClothModifierData *clmd)
 
       /* Activate / Deactivate existing springs */
       if ((!(cloth->verts[spring->ij].flags & CLOTH_VERT_FLAG_PINNED)) &&
-          (cloth->verts[spring->ij].goal > ALMOST_ZERO)) {
+          (cloth->verts[spring->ij].goal > FLT_EPSILON)) {
         spring->flags &= ~CLOTH_SPRING_FLAG_DEACTIVATE;
       }
       else {
@@ -1473,7 +1308,6 @@ static bool cloth_build_springs(ClothModifierData *clmd, Mesh *mesh)
   const MLoop *mloop = mesh->mloop;
   int index2 = 0; /* our second vertex index */
   LinkNodePair *edgelist = NULL;
-  EdgeSet *edgeset = NULL;
   LinkNode *search = NULL, *search2 = NULL;
   BendSpringRef *spring_ref = NULL;
 
@@ -1485,10 +1319,8 @@ static bool cloth_build_springs(ClothModifierData *clmd, Mesh *mesh)
   /* NOTE: handling ownership of springs and edgeset is quite sloppy
    * currently they are never initialized but assert just to be sure */
   BLI_assert(cloth->springs == NULL);
-  BLI_assert(cloth->edgeset == NULL);
 
   cloth->springs = NULL;
-  cloth->edgeset = NULL;
 
   if (clmd->sim_parms->bending_model == CLOTH_BENDING_ANGULAR) {
     spring_ref = MEM_callocN(sizeof(*spring_ref) * numedges, "temp bend spring reference");
@@ -1650,9 +1482,6 @@ static bool cloth_build_springs(ClothModifierData *clmd, Mesh *mesh)
     }
   }
 
-  edgeset = BLI_edgeset_new_ex(__func__, numedges);
-  cloth->edgeset = edgeset;
-
   if (numpolys) {
     for (int i = 0; i < numpolys; i++) {
       /* Shear springs. */
@@ -1754,6 +1583,7 @@ static bool cloth_build_springs(ClothModifierData *clmd, Mesh *mesh)
     if (clmd->sim_parms->bending_model == CLOTH_BENDING_LINEAR) {
       search2 = cloth->springs;
 
+      EdgeSet *edgeset = BLI_edgeset_new_ex(__func__, numedges);
       for (int i = struct_springs; i < struct_springs + shear_springs; i++) {
         if (!search2) {
           break;
@@ -1796,6 +1626,7 @@ static bool cloth_build_springs(ClothModifierData *clmd, Mesh *mesh)
 
         search2 = search2->next;
       }
+      BLI_edgeset_free(edgeset);
     }
   }
   else if (struct_springs > 2) {
@@ -1875,20 +1706,6 @@ static bool cloth_build_springs(ClothModifierData *clmd, Mesh *mesh)
     cloth_hair_update_bending_rest_targets(clmd);
   }
 
-  /* note: the edges may already exist so run reinsert */
-
-  /* insert other near springs in edgeset AFTER bending springs are calculated (for selfcolls) */
-  for (int i = 0; i < numedges; i++) { /* struct springs */
-    BLI_edgeset_add(edgeset, medge[i].v1, medge[i].v2);
-  }
-
-  for (int i = 0; i < numpolys; i++) { /* edge springs */
-    if (mpoly[i].totloop == 4) {
-      BLI_edgeset_add(edgeset, mloop[mpoly[i].loopstart + 0].v, mloop[mpoly[i].loopstart + 2].v);
-      BLI_edgeset_add(edgeset, mloop[mpoly[i].loopstart + 1].v, mloop[mpoly[i].loopstart + 3].v);
-    }
-  }
-
   MEM_SAFE_FREE(spring_ref);
 
   cloth->numsprings = struct_springs + shear_springs + bend_springs;

source/blender/blenkernel/intern/effect.c

@@ -93,7 +93,7 @@ PartDeflect *BKE_partdeflect_new(int type)
   pd->pdef_sbdamp = 0.1f;
   pd->pdef_sbift = 0.2f;
   pd->pdef_sboft = 0.02f;
-  pd->pdef_cfrict = 5.0f;
+  pd->pdef_cfrict = 0.05f;
   pd->seed = ((uint)(ceil(PIL_check_seconds_timer())) + 1) % 128;
   pd->f_strength = 1.0f;
   pd->f_damp = 1.0f;
@@ -115,7 +115,7 @@ PartDeflect *BKE_partdeflect_new(int type)
       pd->f_flow = 1.0f;
       break;
   }
-  pd->flag = PFIELD_DO_LOCATION | PFIELD_DO_ROTATION | PFIELD_CLOTH_USE_CULLING;
+  pd->flag = PFIELD_DO_LOCATION | PFIELD_DO_ROTATION;
 
   return pd;
 }
@@ -452,14 +452,14 @@ static float eff_calc_visibility(ListBase *colliders,
     if (col->ob == eff->ob) {
       continue;
     }
-    if (collmd->bvhtree) {
+    if (collmd->bvhtree_tris) {
       BVHTreeRayHit hit;
 
       hit.index = -1;
       hit.dist = len + FLT_EPSILON;
 
       /* check if the way is blocked */
-      if (BLI_bvhtree_ray_cast_ex(collmd->bvhtree,
+      if (BLI_bvhtree_ray_cast_ex(collmd->bvhtree_tris,
                                   point->loc,
                                   norm,
                                   0.0f,

source/blender/blenkernel/intern/particle_system.c

@@ -2870,8 +2870,8 @@ static int collision_detect(ParticleData *pa,
     col->fac2 = (col->cfra - coll->collmd->time_x) /
                 (coll->collmd->time_xnew - coll->collmd->time_x);
 
-    if (col->md && col->md->bvhtree) {
+    if (col->md && col->md->bvhtree_tris) {
-      BLI_bvhtree_ray_cast_ex(col->md->bvhtree,
+      BLI_bvhtree_ray_cast_ex(col->md->bvhtree_tris,
                               col->co1,
                               ray_dir,
                               col->radius,

source/blender/blenkernel/intern/pointcache.c

@@ -2874,6 +2874,7 @@ void BKE_ptcache_id_time(
     }
   }
 }
+
 int BKE_ptcache_id_reset(Scene *scene, PTCacheID *pid, int mode)
 {
   PointCache *cache;

source/blender/blenlib/BLI_math_geom.h

@@ -215,6 +215,8 @@ void closest_to_plane_v3(float r_close[3], const float plane[4], const float pt[
 void closest_to_plane3_normalized_v3(float r_close[3], const float plane[3], const float pt[3]);
 void closest_to_plane3_v3(float r_close[3], const float plane[3], const float pt[3]);
 
+void interp_weights_closest_on_tri_v3(
+    const float p[3], const float v1[3], const float v2[3], const float v3[3], float r_uv[2]);
 /* Set 'r' to the point in triangle (v1, v2, v3) closest to point 'p' */
 void closest_on_tri_to_point_v3(
     float r[3], const float p[3], const float v1[3], const float v2[3], const float v3[3]);
@@ -264,8 +266,8 @@ void isect_seg_seg_v3(const float a0[3],
                       const float a1[3],
                       const float b0[3],
                       const float b1[3],
-                      float r_a[3],
+                      float *r_lambda_a,
-                      float r_b[3]);
+                      float *r_lambda_b);
 
 int isect_seg_seg_v2_int(const int v1[2], const int v2[2], const int v3[2], const int v4[2]);
 int isect_seg_seg_v2_point_ex(const float v0[2],
@@ -413,10 +415,10 @@ bool isect_ray_tri_epsilon_v3(const float ray_origin[3],
                               float *r_lambda,
                               float r_uv[2],
                               const float epsilon);
-bool isect_tri_tri_v3_ex(const float tri_a[3][3],
+bool isect_tri_tri_v3_ex(const float ta[3][3],
-                         const float tri_b[3][3],
+                         const float tb[3][3],
-                         float r_i1[3],
+                         float r_w0[3],
-                         float r_i2[3],
+                         float r_w1[3],
                          int *r_tri_a_edge_isect_count);
 bool isect_tri_tri_v3(const float t_a0[3],
                       const float t_a1[3],
@@ -424,8 +426,8 @@ bool isect_tri_tri_v3(const float t_a0[3],
                       const float t_b0[3],
                       const float t_b1[3],
                       const float t_b2[3],
-                      float r_i1[3],
+                      float r_i0[3],
-                      float r_i2[3]);
+                      float r_i1[3]);
 
 bool isect_tri_tri_v2(const float p1[2],
                       const float q1[2],

source/blender/blenlib/BLI_math_matrix.h

@@ -245,6 +245,7 @@ void transpose_m3_m4(float R[3][3], const float M[4][4]);
 void transpose_m4(float R[4][4]);
 void transpose_m4_m4(float R[4][4], const float M[4][4]);
 
+bool compare_m3m3(const float mat1[3][3], const float mat2[3][3], float limit);
 bool compare_m4m4(const float mat1[4][4], const float mat2[4][4], float limit);
 
 void normalize_m2_ex(float R[2][2], float r_scale[2]) ATTR_NONNULL();

source/blender/blenlib/BLI_math_solvers.h

@@ -68,6 +68,28 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta,
                         const float x_init[3],
                         float result[3]);
 
+bool collision_tri_point_bisect(const float tri_time0[3][3],
+                                const float tri_time1[3][3],
+                                const float point_time0[3],
+                                const float point_time1[3],
+                                const float point_radius,
+                                const int max_bisect_steps,
+                                float *r_time,
+                                float *r_dist_sq,
+                                float r_tri[3],
+                                float r_point[3]);
+
+bool collision_seg_seg_bisect(const float seg_a_time0[2][3],
+                              const float seg_a_time1[2][3],
+                              const float seg_b_time0[2][3],
+                              const float seg_b_time1[2][3],
+                              const float collision_dist,
+                              const int max_bisect_steps,
+                              float *r_time,
+                              float *r_dist_sq,
+                              float r_seg_a[3],
+                              float r_seg_b[3]);
+
 #ifdef BLI_MATH_GCC_WARN_PRAGMA
 #  pragma GCC diagnostic pop
 #endif

source/blender/blenlib/intern/math_geom.c

@@ -1017,6 +1017,81 @@ float dist_squared_to_projected_aabb_simple(const float projmat[4][4],
 }
 /** \} */
 
+/* Adapted from "Real-Time Collision Detection" by Christer Ericson,
+ * published by Morgan Kaufmann Publishers, copyright 2005 Elsevier Inc.
+ *
+ * Set #r_uv to be the baricentric weights of the closest point in the triangle. */
+void interp_weights_closest_on_tri_v3(
+    const float p[3], const float v1[3], const float v2[3], const float v3[3], float r_uv[2])
+{
+  float ab[3], ac[3], ap[3], d1, d2;
+  float bp[3], d3, d4, vc, cp[3], d5, d6, vb, va;
+  float denom, u, v, w;
+
+  /* Check if P in vertex region outside A */
+  sub_v3_v3v3(ab, v2, v1);
+  sub_v3_v3v3(ac, v3, v1);
+  sub_v3_v3v3(ap, p, v1);
+  d1 = dot_v3v3(ab, ap);
+  d2 = dot_v3v3(ac, ap);
+  if (d1 <= 0.0f && d2 <= 0.0f) {
+    /* barycentric coordinates (1,0,0) */
+    copy_v2_fl2(r_uv, 1.0f, 0.0f);
+    return;
+  }
+
+  /* Check if P in vertex region outside B */
+  sub_v3_v3v3(bp, p, v2);
+  d3 = dot_v3v3(ab, bp);
+  d4 = dot_v3v3(ac, bp);
+  if (d3 >= 0.0f && d4 <= d3) {
+    /* barycentric coordinates (0,1,0) */
+    copy_v2_fl2(r_uv, 0.0f, 1.0f);
+    return;
+  }
+  /* Check if P in edge region of AB, if so return projection of P onto AB */
+  vc = d1 * d4 - d3 * d2;
+  if (vc <= 0.0f && d1 >= 0.0f && d3 <= 0.0f) {
+    v = d1 / (d1 - d3);
+    /* barycentric coordinates (1-v,v,0) */
+    copy_v2_fl2(r_uv, 1.0f - v, v);
+    return;
+  }
+  /* Check if P in vertex region outside C */
+  sub_v3_v3v3(cp, p, v3);
+  d5 = dot_v3v3(ab, cp);
+  d6 = dot_v3v3(ac, cp);
+  if (d6 >= 0.0f && d5 <= d6) {
+    /* barycentric coordinates (0,0,1) */
+    copy_v2_fl2(r_uv, 0.0f, 0.0f);
+    return;
+  }
+  /* Check if P in edge region of AC, if so return projection of P onto AC */
+  vb = d5 * d2 - d1 * d6;
+  if (vb <= 0.0f && d2 >= 0.0f && d6 <= 0.0f) {
+    w = d2 / (d2 - d6);
+    /* barycentric coordinates (1-w,0,w) */
+    copy_v2_fl2(r_uv, 1.0f - w, 0.0f);
+    return;
+  }
+  /* Check if P in edge region of BC, if so return projection of P onto BC */
+  va = d3 * d6 - d5 * d4;
+  if (va <= 0.0f && (d4 - d3) >= 0.0f && (d5 - d6) >= 0.0f) {
+    w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
+    /* barycentric coordinates (0,1-w,w) */
+    copy_v2_fl2(r_uv, 0.0f, 1.0f - w);
+    return;
+  }
+
+  denom = 1.0f / (va + vb + vc);
+  u = va * denom;
+  v = vb * denom;
+  // w = vc * denom = 1.0f - u - v ;
+
+  /* barycentric coordinates (u,v,w) */
+  copy_v2_fl2(r_uv, u, v);
+}
+
 /* Adapted from "Real-Time Collision Detection" by Christer Ericson,
  * published by Morgan Kaufmann Publishers, copyright 2005 Elsevier Inc.
  *
@@ -1176,65 +1251,100 @@ void isect_seg_seg_v3(const float a0[3],
                       const float a1[3],
                       const float b0[3],
                       const float b1[3],
-                      float r_a[3],
+                      float *r_lambda_a,
-                      float r_b[3])
+                      float *r_lambda_b)
 {
-  float fac_a, fac_b;
+  float a_dir[3], b_dir[3];
-  float a_dir[3], b_dir[3], a0b0[3], crs_ab[3];
   sub_v3_v3v3(a_dir, a1, a0);
   sub_v3_v3v3(b_dir, b1, b0);
-  sub_v3_v3v3(a0b0, b0, a0);
-  cross_v3_v3v3(crs_ab, b_dir, a_dir);
-  const float nlen = len_squared_v3(crs_ab);
 
-  if (nlen == 0.0f) {
+  float lambda_a, lambda_b;
-    /* Parallel Lines */
+  if (isect_ray_ray_v3(a0, a_dir, b0, b_dir, &lambda_a, &lambda_b)) {
-    /* In this case return any point that
+    bool has_clamp_a = true;
-     * is between the closest segments. */
+    bool has_clamp_b = true;
-    float a0b1[3], a1b0[3], len_a, len_b, fac1, fac2;
+    if (lambda_a < 0.0f) {
-    sub_v3_v3v3(a0b1, b1, a0);
+      lambda_a = 0.0f;
-    sub_v3_v3v3(a1b0, b0, a1);
+    }
-    len_a = len_squared_v3(a_dir);
+    else if (lambda_a > 1.0f) {
-    len_b = len_squared_v3(b_dir);
+      lambda_a = 1.0f;
-
-    if (len_a) {
-      fac1 = dot_v3v3(a0b0, a_dir);
-      fac2 = dot_v3v3(a0b1, a_dir);
-      CLAMP(fac1, 0.0f, len_a);
-      CLAMP(fac2, 0.0f, len_a);
-      fac_a = (fac1 + fac2) / (2 * len_a);
     }
     else {
-      fac_a = 0.0f;
+      has_clamp_a = false;
     }
 
-    if (len_b) {
+    if (lambda_b < 0.0f) {
-      fac1 = -dot_v3v3(a0b0, b_dir);
+      lambda_b = 0.0f;
-      fac2 = -dot_v3v3(a1b0, b_dir);
+    }
-      CLAMP(fac1, 0.0f, len_b);
+    else if (lambda_b > 1.0f) {
-      CLAMP(fac2, 0.0f, len_b);
+      lambda_b = 1.0f;
-      fac_b = (fac1 + fac2) / (2 * len_b);
     }
     else {
-      fac_b = 0.0f;
+      has_clamp_b = false;
+    }
+
+    if (has_clamp_a || has_clamp_b) {
+      float tmp[3];
+      if (has_clamp_a && has_clamp_b) {
+        /* Only one of these clamps is really valid. */
+        lambda_b = closest_to_ray_v3(tmp, lambda_a == 0.0f ? a0 : a1, b0, b_dir);
+        if (IN_RANGE_INCL(lambda_b, 0.0f, 1.0f)) {
+          BLI_assert(IN_RANGE_INCL(lambda_a, 0.0f, 1.0f));
+          *r_lambda_a = lambda_a;
+          *r_lambda_b = lambda_b;
+          return;
+        }
+        CLAMP(lambda_b, 0.0f, 1.0f);
+        has_clamp_a = false;
+      }
+      if (has_clamp_a) {
+        lambda_b = closest_to_ray_v3(tmp, lambda_a == 0.0f ? a0 : a1, b0, b_dir);
+        CLAMP(lambda_b, 0.0f, 1.0f);
+      }
+      else {
+        lambda_a = closest_to_ray_v3(tmp, lambda_b == 0.0f ? b0 : b1, a0, a_dir);
+        CLAMP(lambda_a, 0.0f, 1.0f);
+      }
+      BLI_assert(IN_RANGE_INCL(lambda_a, 0.0f, 1.0f) && IN_RANGE_INCL(lambda_b, 0.0f, 1.0f));
+      *r_lambda_a = lambda_a;
+      *r_lambda_b = lambda_b;
+      return;
     }
   }
   else {
-    float c[3], cray[3];
+    /* Parallel Lines */
-    sub_v3_v3v3(c, crs_ab, a0b0);
+    /* In this case return any point that
+     * is between the closest segments. */
+    float a0b0[3], a0b1[3];
+    sub_v3_v3v3(a0b0, b0, a0);
+    sub_v3_v3v3(a0b1, b1, a0);
+    float len_sq_a = len_squared_v3(a_dir);
+    if (len_sq_a) {
+      float da_a0b0 = dot_v3v3(a_dir, a0b0);
+      float da_a0b1 = dot_v3v3(a_dir, a0b0);
+      CLAMP(da_a0b0, 0.0f, len_sq_a);
+      CLAMP(da_a0b1, 0.0f, len_sq_a);
+      lambda_a = (da_a0b0 + da_a0b1) / (2 * len_sq_a);
+    }
+    else {
+      lambda_a = 0.0f;
+    }
 
-    cross_v3_v3v3(cray, c, b_dir);
+    float len_sq_b = len_squared_v3(b_dir);
-    fac_a = dot_v3v3(cray, crs_ab) / nlen;
+    if (len_sq_b) {
-
+      float db_b0a0 = -dot_v3v3(b_dir, a0b0);
-    cross_v3_v3v3(cray, c, a_dir);
+      float db_b0a1 = -dot_v3v3(b_dir, a0b0);
-    fac_b = dot_v3v3(cray, crs_ab) / nlen;
+      CLAMP(db_b0a0, 0.0f, len_sq_b);
-
+      CLAMP(db_b0a1, 0.0f, len_sq_b);
-    CLAMP(fac_a, 0.0f, 1.0f);
+      lambda_b = (db_b0a0 + db_b0a1) / (2 * len_sq_b);
-    CLAMP(fac_b, 0.0f, 1.0f);
+    }
+    else {
+      lambda_b = 0.0f;
+    }
   }
 
-  madd_v3_v3v3fl(r_a, a0, a_dir, fac_a);
+  BLI_assert(IN_RANGE_INCL(lambda_a, 0.0f, 1.0f) && IN_RANGE_INCL(lambda_b, 0.0f, 1.0f));
-  madd_v3_v3v3fl(r_b, b0, b_dir, fac_b);
+  *r_lambda_a = lambda_a;
+  *r_lambda_b = lambda_b;
 }
 
 /**
@@ -2382,10 +2492,10 @@ bool isect_planes_v3_fn(
  * \note If it exists, \a r_i1 will be a point on the edge of the 1st triangle.
  * \note intersections between coplanar triangles are currently undetected.
  */
-bool isect_tri_tri_v3_ex(const float tri_a[3][3],
+bool isect_tri_tri_v3_ex(const float ta[3][3],
-                         const float tri_b[3][3],
+                         const float tb[3][3],
-                         float r_i1[3],
+                         float r_w1[3],
-                         float r_i2[3],
+                         float r_w2[3],
                          int *r_tri_a_edge_isect_count)
 {
   struct {
@@ -2393,20 +2503,20 @@ bool isect_tri_tri_v3_ex(const float tri_a[3][3],
      * that intersects the planes of the triangles. */
     float min, max;
     /* Intersection point location. */
-    float loc[2][3];
+    float w[2][3];
   } range[2];
 
   float side[2][3];
   double ba[3], bc[3], plane_a[4], plane_b[4];
   *r_tri_a_edge_isect_count = 0;
 
-  sub_v3db_v3fl_v3fl(ba, tri_a[0], tri_a[1]);
+  sub_v3db_v3fl_v3fl(ba, ta[0], ta[1]);
-  sub_v3db_v3fl_v3fl(bc, tri_a[2], tri_a[1]);
+  sub_v3db_v3fl_v3fl(bc, ta[2], ta[1]);
   cross_v3_v3v3_db(plane_a, ba, bc);
-  plane_a[3] = -dot_v3db_v3fl(plane_a, tri_a[1]);
+  plane_a[3] = -dot_v3db_v3fl(plane_a, ta[1]);
-  side[1][0] = (float)(dot_v3db_v3fl(plane_a, tri_b[0]) + plane_a[3]);
+  side[1][0] = (float)(dot_v3db_v3fl(plane_a, tb[0]) + plane_a[3]);
-  side[1][1] = (float)(dot_v3db_v3fl(plane_a, tri_b[1]) + plane_a[3]);
+  side[1][1] = (float)(dot_v3db_v3fl(plane_a, tb[1]) + plane_a[3]);
-  side[1][2] = (float)(dot_v3db_v3fl(plane_a, tri_b[2]) + plane_a[3]);
+  side[1][2] = (float)(dot_v3db_v3fl(plane_a, tb[2]) + plane_a[3]);
 
   if (!side[1][0] && !side[1][1] && !side[1][2]) {
     /* Coplanar case is not supported. */
@@ -2420,13 +2530,13 @@ bool isect_tri_tri_v3_ex(const float tri_a[3][3],
     return false;
   }
 
-  sub_v3db_v3fl_v3fl(ba, tri_b[0], tri_b[1]);
+  sub_v3db_v3fl_v3fl(ba, tb[0], tb[1]);
-  sub_v3db_v3fl_v3fl(bc, tri_b[2], tri_b[1]);
+  sub_v3db_v3fl_v3fl(bc, tb[2], tb[1]);
   cross_v3_v3v3_db(plane_b, ba, bc);
-  plane_b[3] = -dot_v3db_v3fl(plane_b, tri_b[1]);
+  plane_b[3] = -dot_v3db_v3fl(plane_b, tb[1]);
-  side[0][0] = (float)(dot_v3db_v3fl(plane_b, tri_a[0]) + plane_b[3]);
+  side[0][0] = (float)(dot_v3db_v3fl(plane_b, ta[0]) + plane_b[3]);
-  side[0][1] = (float)(dot_v3db_v3fl(plane_b, tri_a[1]) + plane_b[3]);
+  side[0][1] = (float)(dot_v3db_v3fl(plane_b, ta[1]) + plane_b[3]);
-  side[0][2] = (float)(dot_v3db_v3fl(plane_b, tri_a[2]) + plane_b[3]);
+  side[0][2] = (float)(dot_v3db_v3fl(plane_b, ta[2]) + plane_b[3]);
 
   if ((side[0][0] && side[0][1] && side[0][2]) && (side[0][0] < 0.0f) == (side[0][1] < 0.0f) &&
       (side[0][0] < 0.0f) == (side[0][2] < 0.0f)) {
@@ -2439,7 +2549,18 @@ bool isect_tri_tri_v3_ex(const float tri_a[3][3],
   double isect_dir[3];
   cross_v3_v3v3_db(isect_dir, plane_a, plane_b);
   for (int i = 0; i < 2; i++) {
-    const float(*tri)[3] = i == 0 ? tri_a : tri_b;
+    float tri[3][3];
+    if (i == 0) {
+      copy_v3_v3(tri[0], ta[0]);
+      copy_v3_v3(tri[1], ta[1]);
+      copy_v3_v3(tri[2], ta[2]);
+    }
+    else {
+      copy_v3_v3(tri[0], tb[0]);
+      copy_v3_v3(tri[1], tb[1]);
+      copy_v3_v3(tri[2], tb[2]);
+    }
+
     /* Rearrange the triangle so that the vertex that is alone on one side
      * of the plane is located at index 1. */
     int tri_i[3];
@@ -2459,13 +2580,15 @@ bool isect_tri_tri_v3_ex(const float tri_a[3][3],
       tri_i[2] = 2;
     }
 
+    float fac0 = 0.0f;
+    float fac1 = 0.0f;
     double dot_b = dot_v3db_v3fl(isect_dir, tri[tri_i[1]]);
     float sidec = side[i][tri_i[1]];
     if (sidec) {
       double dot_a = dot_v3db_v3fl(isect_dir, tri[tri_i[0]]);
       double dot_c = dot_v3db_v3fl(isect_dir, tri[tri_i[2]]);
-      float fac0 = sidec / (sidec - side[i][tri_i[0]]);
+      fac0 = sidec / (sidec - side[i][tri_i[0]]);
-      float fac1 = sidec / (sidec - side[i][tri_i[2]]);
+      fac1 = sidec / (sidec - side[i][tri_i[2]]);
       double offset0 = fac0 * (dot_a - dot_b);
       double offset1 = fac1 * (dot_c - dot_b);
       if (offset0 > offset1) {
@@ -2477,14 +2600,18 @@ bool isect_tri_tri_v3_ex(const float tri_a[3][3],
 
       range[i].min = (float)(dot_b + offset0);
       range[i].max = (float)(dot_b + offset1);
-      interp_v3_v3v3(range[i].loc[0], tri[tri_i[1]], tri[tri_i[0]], fac0);
-      interp_v3_v3v3(range[i].loc[1], tri[tri_i[1]], tri[tri_i[2]], fac1);
     }
     else {
-      range[i].min = range[i].max = (float)dot_b;
+      range[i].min = range[i].max = (float)(dot_b);
-      copy_v3_v3(range[i].loc[0], tri[tri_i[1]]);
-      copy_v3_v3(range[i].loc[1], tri[tri_i[1]]);
     }
+
+    range[i].w[0][tri_i[0]] = fac0;
+    range[i].w[0][tri_i[1]] = 1.0f - fac0;
+    range[i].w[0][tri_i[2]] = 0.0f;
+
+    range[i].w[1][tri_i[0]] = 0.0f;
+    range[i].w[1][tri_i[1]] = 1.0f - fac1;
+    range[i].w[1][tri_i[2]] = fac1;
   }
 
   if ((range[0].max > range[1].min) && (range[0].min < range[1].max)) {
@@ -2492,25 +2619,25 @@ bool isect_tri_tri_v3_ex(const float tri_a[3][3],
      * Now identify the two points of intersection that are in the middle to get the actual
      * intersection between the triangles. (B--C from A--B--C--D) */
     if (range[0].min >= range[1].min) {
-      copy_v3_v3(r_i1, range[0].loc[0]);
+      copy_v3_v3(r_w1, range[0].w[0]);
       if (range[0].max <= range[1].max) {
-        copy_v3_v3(r_i2, range[0].loc[1]);
+        copy_v3_v3(r_w2, range[0].w[1]);
         *r_tri_a_edge_isect_count = 2;
       }
       else {
-        copy_v3_v3(r_i2, range[1].loc[1]);
+        copy_v3_v3(r_w2, range[1].w[1]);
         *r_tri_a_edge_isect_count = 1;
       }
     }
     else {
       if (range[0].max <= range[1].max) {
-        copy_v3_v3(r_i1, range[0].loc[1]);
+        copy_v3_v3(r_w1, range[0].w[1]);
-        copy_v3_v3(r_i2, range[1].loc[0]);
+        copy_v3_v3(r_w2, range[1].w[0]);
         *r_tri_a_edge_isect_count = 1;
       }
       else {
-        copy_v3_v3(r_i1, range[1].loc[0]);
+        copy_v3_v3(r_w1, range[1].w[0]);
-        copy_v3_v3(r_i2, range[1].loc[1]);
+        copy_v3_v3(r_w2, range[1].w[1]);
       }
     }
     return true;
@@ -2525,18 +2652,34 @@ bool isect_tri_tri_v3(const float t_a0[3],
                       const float t_b0[3],
                       const float t_b1[3],
                       const float t_b2[3],
-                      float r_i1[3],
+                      float r_i0[3],
-                      float r_i2[3])
+                      float r_i1[3])
 {
-  float tri_a[3][3], tri_b[3][3];
+  float w0[3], w1[3];
-  int dummy;
+  int tri_a_edge_isect_count;
-  copy_v3_v3(tri_a[0], t_a0);
+  float tria[3][3], trib[3][3];
-  copy_v3_v3(tri_a[1], t_a1);
+  copy_v3_v3(tria[0], t_a0);
-  copy_v3_v3(tri_a[2], t_a2);
+  copy_v3_v3(tria[1], t_a1);
-  copy_v3_v3(tri_b[0], t_b0);
+  copy_v3_v3(tria[2], t_a2);
-  copy_v3_v3(tri_b[1], t_b1);
+  copy_v3_v3(trib[0], t_b0);
-  copy_v3_v3(tri_b[2], t_b2);
+  copy_v3_v3(trib[1], t_b1);
-  return isect_tri_tri_v3_ex(tri_a, tri_b, r_i1, r_i2, &dummy);
+  copy_v3_v3(trib[2], t_b2);
+  if (isect_tri_tri_v3_ex(tria, trib, w0, w1, &tri_a_edge_isect_count)) {
+    if (tri_a_edge_isect_count == 0) {
+      interp_barycentric_tri_v3(trib, UNPACK2(w0), r_i0);
+      interp_barycentric_tri_v3(trib, UNPACK2(w1), r_i1);
+    }
+    else if (tri_a_edge_isect_count == 1) {
+      interp_barycentric_tri_v3(tria, UNPACK2(w0), r_i0);
+      interp_barycentric_tri_v3(trib, UNPACK2(w1), r_i1);
+    }
+    else {
+      interp_barycentric_tri_v3(tria, UNPACK2(w0), r_i0);
+      interp_barycentric_tri_v3(tria, UNPACK2(w1), r_i1);
+    }
+    return true;
+  }
+  return false;
 }
 
 /* -------------------------------------------------------------------- */
@@ -3186,7 +3329,7 @@ bool isect_ray_ray_epsilon_v3(const float ray_origin_a[3],
   const float nlen = len_squared_v3(n);
 
   /* `nlen` is the square of the area formed by the two vectors. */
-  if (UNLIKELY(nlen < epsilon)) {
+  if (UNLIKELY(nlen <= epsilon)) {
     /* The lines are parallel. */
     return false;
   }

source/blender/blenlib/intern/math_matrix.c

@@ -1403,6 +1403,18 @@ void transpose_m4_m4(float R[4][4], const float M[4][4])
   R[3][3] = M[3][3];
 }
 
+bool compare_m3m3(const float mat1[3][3], const float mat2[3][3], float limit)
+{
+  if (compare_v3v3(mat1[0], mat2[0], limit)) {
+    if (compare_v3v3(mat1[1], mat2[1], limit)) {
+      if (compare_v3v3(mat1[2], mat2[2], limit)) {
+        return true;
+      }
+    }
+  }
+  return false;
+}
+
 bool compare_m4m4(const float mat1[4][4], const float mat2[4][4], float limit)
 {
   if (compare_v4v4(mat1[0], mat2[0], limit)) {

source/blender/blenlib/intern/math_solvers.c

@@ -279,3 +279,447 @@ bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta,
   copy_v3_v3(result, x);
   return success;
 }
+
+bool collision_tri_point_bisect(const float tri_time0[3][3],
+                                const float tri_time1[3][3],
+                                const float point_time0[3],
+                                const float point_time1[3],
+                                const float collision_dist,
+                                int max_bisect_steps,
+                                float *r_time,
+                                float *r_dist_sq,
+                                float r_tri[3],
+                                float r_point[3])
+{
+  /* Compensate movement. */
+  float tri0[3][3], tri1[3][3];
+  sub_v3_v3v3(tri0[0], tri_time0[0], point_time0);
+  sub_v3_v3v3(tri0[1], tri_time0[1], point_time0);
+  sub_v3_v3v3(tri0[2], tri_time0[2], point_time0);
+
+  sub_v3_v3v3(tri1[0], tri_time1[0], point_time1);
+  sub_v3_v3v3(tri1[1], tri_time1[1], point_time1);
+  sub_v3_v3v3(tri1[2], tri_time1[2], point_time1);
+
+  const float zero[3] = {0.0f, 0.0f, 0.0f};
+  const float dist_sq_best = square_f(collision_dist);
+
+  /* Bisect. */
+  struct {
+    float time;
+    float uv[2];
+    float dist_sq;
+  } tmin, tmax, test_1, test_2;
+
+  float vel_rel[3];
+  float time_offset;
+  float move_dist_sq;
+  float dir_curr[3];
+  float tmp_tri[3][3];
+
+  tmin.time = 0.0f;
+  interp_weights_closest_on_tri_v3(zero, UNPACK3(tri0), tmin.uv);
+  interp_barycentric_tri_v3(tri0, UNPACK2(tmin.uv), dir_curr);
+  tmin.dist_sq = len_squared_v3(dir_curr);
+
+  if (tmin.dist_sq <= dist_sq_best) {
+    /* Ditance alread inside margin, so report collision. */
+    copy_v3_v3(r_point, point_time0);
+    add_v3_v3v3(r_tri, r_point, dir_curr);
+    *r_dist_sq = tmin.dist_sq;
+    *r_time = 0.0f;
+    return true;
+  }
+
+  /* Loop. */
+  float tri_vel[3][3];
+  sub_v3_v3v3(tri_vel[0], tri1[0], tri0[0]);
+  sub_v3_v3v3(tri_vel[1], tri1[1], tri0[1]);
+  sub_v3_v3v3(tri_vel[2], tri1[2], tri0[2]);
+
+  /* Test linear motion. */
+  interp_barycentric_tri_v3(tri_vel, UNPACK2(tmin.uv), vel_rel);
+  float vel_rel_len_sq = len_squared_v3(vel_rel);
+  if (vel_rel_len_sq < FLT_EPSILON) {
+    /* Motionless. */
+    tmax = tmin;
+  }
+  else {
+    tmax.time = 1.0f;
+    tmax.dist_sq = FLT_MAX;
+    test_1.time = -dot_v3v3(dir_curr, vel_rel) / vel_rel_len_sq;
+    while (--max_bisect_steps) {
+      float delta_time = tmax.time - tmin.time;
+      float delta_time_sq = square_f(delta_time);
+      if (test_1.time <= tmin.time) {
+        test_1.time = tmin.time + delta_time * (1.0f / 4.0f);
+      }
+      else if (test_1.time >= tmax.time) {
+        test_1.time = tmin.time + delta_time * (3.0f / 4.0f);
+      }
+
+      madd_v3_v3v3fl(tmp_tri[0], tri0[0], tri_vel[0], test_1.time);
+      madd_v3_v3v3fl(tmp_tri[1], tri0[1], tri_vel[1], test_1.time);
+      madd_v3_v3v3fl(tmp_tri[2], tri0[2], tri_vel[2], test_1.time);
+      interp_weights_closest_on_tri_v3(zero, UNPACK3(tmp_tri), test_1.uv);
+      interp_barycentric_tri_v3(tmp_tri, UNPACK2(test_1.uv), dir_curr);
+      test_1.dist_sq = len_squared_v3(dir_curr);
+      if (test_1.dist_sq <= dist_sq_best) {
+        tmax = test_1;
+        break;
+      }
+
+      interp_barycentric_tri_v3(tri_vel, UNPACK2(test_1.uv), vel_rel);
+      vel_rel_len_sq = len_squared_v3(vel_rel);
+      move_dist_sq = vel_rel_len_sq * delta_time_sq;
+      if (move_dist_sq < dist_sq_best) {
+        /* The range is already very small. */
+        break;
+      }
+
+      time_offset = -dot_v3v3(dir_curr, vel_rel) / vel_rel_len_sq;
+      test_2.time = test_1.time + time_offset;
+      if (test_2.time <= tmin.time) {
+        test_2.time = (test_1.time + tmin.time) * 0.5f;
+      }
+      else if (test_2.time >= tmax.time) {
+        test_2.time = (tmax.time + test_1.time) * 0.5f;
+      }
+
+      madd_v3_v3v3fl(tmp_tri[0], tri0[0], tri_vel[0], test_2.time);
+      madd_v3_v3v3fl(tmp_tri[1], tri0[1], tri_vel[1], test_2.time);
+      madd_v3_v3v3fl(tmp_tri[2], tri0[2], tri_vel[2], test_2.time);
+      interp_weights_closest_on_tri_v3(zero, UNPACK3(tmp_tri), test_2.uv);
+      interp_barycentric_tri_v3(tmp_tri, UNPACK2(test_2.uv), dir_curr);
+      test_2.dist_sq = len_squared_v3(dir_curr);
+      if (test_2.dist_sq <= dist_sq_best) {
+        tmax = test_2;
+        break;
+        // if (test_1.time < max.time) {
+        //  min = test_1;
+        //}
+      }
+
+      if (test_1.time <= test_2.time) {
+        if (test_1.dist_sq < test_2.dist_sq) {
+          /* Moving away. */
+          tmax = test_2;
+        }
+        else {
+          /* Aproaxing. */
+          tmin = test_1;
+        }
+      }
+      else {
+        if (test_2.dist_sq < test_1.dist_sq) {
+          /* Moving away. */
+          tmax = test_1;
+        }
+        else {
+          /* Aproaxing. */
+          tmin = test_2;
+        }
+      }
+
+      /* Calc a new time to test. */
+      interp_barycentric_tri_v3(tri_vel, UNPACK2(test_2.uv), vel_rel);
+      vel_rel_len_sq = len_squared_v3(vel_rel);
+      move_dist_sq = vel_rel_len_sq * delta_time_sq;
+      if (move_dist_sq < dist_sq_best) {
+        /* The range is already very small. */
+        break;
+      }
+
+      time_offset = -dot_v3v3(dir_curr, vel_rel) / vel_rel_len_sq;
+      test_1.time = test_2.time + time_offset;
+    }
+  }
+
+  if (tmax.dist_sq > dist_sq_best) {
+    return false;
+  }
+
+  {
+    /* Last refinement. */
+    interp_barycentric_tri_v3(tri_vel, UNPACK2(tmax.uv), vel_rel);
+    BLI_assert(!is_zero_v3(vel_rel));
+
+    float vel_rel_lensq_inv = 1.0f / len_squared_v3(vel_rel);
+    float dir_to_closest[3], p_close[3];
+    if (max_bisect_steps == 0) {
+      madd_v3_v3v3fl(tmp_tri[0], tri0[0], tri_vel[0], tmax.time);
+      madd_v3_v3v3fl(tmp_tri[1], tri0[1], tri_vel[1], tmax.time);
+      madd_v3_v3v3fl(tmp_tri[2], tri0[2], tri_vel[2], tmax.time);
+      interp_barycentric_tri_v3(tmp_tri, UNPACK2(tmax.uv), dir_curr);
+    }
+
+    float mul = -dot_v3v3(dir_curr, vel_rel) * vel_rel_lensq_inv;
+    mul_v3_v3fl(dir_to_closest, vel_rel, mul);
+    add_v3_v3v3(p_close, dir_curr, dir_to_closest);
+
+    float dvec_len_sq = dist_sq_best - len_squared_v3(p_close);
+    float fac_sq = dvec_len_sq * vel_rel_lensq_inv;
+    float fac = sqrtf(fac_sq);
+
+    float offset_dir[3];
+    mul_v3_v3fl(offset_dir, vel_rel, -fac);
+
+    float dvec[3];
+    add_v3_v3v3(dvec, dir_to_closest, offset_dir);
+
+    time_offset = dot_v3v3(dvec, vel_rel) * vel_rel_lensq_inv;
+    test_1.time = tmax.time + time_offset;
+    if (test_1.time < 0.0f) {
+      test_1.time = 0.0f;
+    }
+    madd_v3_v3v3fl(tmp_tri[0], tri0[0], tri_vel[0], test_1.time);
+    madd_v3_v3v3fl(tmp_tri[1], tri0[1], tri_vel[1], test_1.time);
+    madd_v3_v3v3fl(tmp_tri[2], tri0[2], tri_vel[2], test_1.time);
+    interp_weights_closest_on_tri_v3(zero, UNPACK3(tmp_tri), test_1.uv);
+    interp_barycentric_tri_v3(tmp_tri, UNPACK2(test_1.uv), dir_curr);
+    test_1.dist_sq = len_squared_v3(dir_curr);
+  }
+
+  interp_v3_v3v3(r_point, point_time0, point_time1, test_1.time);
+  add_v3_v3v3(r_tri, r_point, dir_curr);
+  *r_dist_sq = test_1.dist_sq;
+  *r_time = test_1.time;
+  return true;
+}
+
+bool collision_seg_seg_bisect(const float seg_a_time0[2][3],
+                              const float seg_a_time1[2][3],
+                              const float seg_b_time0[2][3],
+                              const float seg_b_time1[2][3],
+                              const float collision_dist,
+                              int max_bisect_steps,
+                              float *r_time,
+                              float *r_dist_sq,
+                              float r_seg_a[3],
+                              float r_seg_b[3])
+{
+  /* Compensate movement. */
+  float eb0[2][3], eb1[2][3], ea[2][3];
+  sub_v3_v3v3(eb0[0], seg_b_time0[0], seg_a_time0[0]);
+  sub_v3_v3v3(eb0[1], seg_b_time0[1], seg_a_time0[0]);
+
+  sub_v3_v3v3(eb1[0], seg_b_time1[0], seg_a_time1[0]);
+  sub_v3_v3v3(eb1[1], seg_b_time1[1], seg_a_time1[0]);
+
+  sub_v3_v3v3(ea[0], seg_a_time0[1], seg_a_time0[0]);
+  sub_v3_v3v3(ea[1], seg_a_time1[1], seg_a_time1[0]);
+
+  const float zero[3] = {0.0f, 0.0f, 0.0f};
+  const float dist_sq_best = square_f(collision_dist);
+
+  /* Bisect. */
+  struct {
+    float time;
+    float wa;
+    float wb;
+    float dist_sq;
+  } tmin, tmax, test_1, test_2;
+
+  float vel_a[3];
+  float vel_b[3];
+  float ra[3];
+  float rb[3];
+  float dir_curr[3];
+  float vel_rel[3];
+  float time_offset;
+  float move_dist_sq;
+
+  tmin.time = 0.0f;
+  isect_seg_seg_v3(zero, ea[0], UNPACK2(eb0), &tmin.wa, &tmin.wb);
+  mul_v3_v3fl(ra, ea[0], tmin.wa);
+  interp_v3_v3v3(rb, UNPACK2(eb0), tmin.wb);
+  sub_v3_v3v3(dir_curr, ra, rb);
+  tmin.dist_sq = len_squared_v3(dir_curr);
+
+  if (tmin.dist_sq <= dist_sq_best) {
+    /* Ditance alread inside margin, so report collision. */
+    add_v3_v3v3(r_seg_a, ra, seg_a_time0[0]);
+    add_v3_v3v3(r_seg_b, rb, seg_a_time0[0]);
+    *r_dist_sq = tmin.dist_sq;
+    *r_time = 0.0f;
+    return true;
+  }
+
+  /* Loop. */
+  struct {
+    float ea[3];
+    float eb[2][3];
+  } tmp;
+
+  /* Test linear motion. */
+  float eb_vel[2][3], ea_vel[3];
+  sub_v3_v3v3(eb_vel[0], eb1[0], eb0[0]);
+  sub_v3_v3v3(eb_vel[1], eb1[1], eb0[1]);
+  sub_v3_v3v3(ea_vel, ea[1], ea[0]);
+  mul_v3_v3fl(vel_a, ea_vel, tmin.wa);
+  interp_v3_v3v3(vel_b, UNPACK2(eb_vel), tmin.wb);
+  sub_v3_v3v3(vel_rel, vel_b, vel_a);
+  float vel_rel_len_sq = len_squared_v3(vel_rel);
+  if (vel_rel_len_sq < FLT_EPSILON) {
+    /* Motionless. */
+    tmax = tmin;
+  }
+  else {
+    tmax.time = 1.0f;
+    tmax.dist_sq = FLT_MAX;
+    test_1.time = dot_v3v3(dir_curr, vel_rel) / vel_rel_len_sq;
+    while (--max_bisect_steps) {
+      float delta_time = tmax.time - tmin.time;
+      float delta_time_sq = square_f(delta_time);
+      if (test_1.time <= tmin.time) {
+        test_1.time = tmin.time + delta_time * (1.0f / 4.0f);
+      }
+      else if (test_1.time >= tmax.time) {
+        test_1.time = tmin.time + delta_time * (3.0f / 4.0f);
+      }
+
+      madd_v3_v3v3fl(tmp.ea, ea[0], ea_vel, test_1.time);
+      madd_v3_v3v3fl(tmp.eb[0], eb0[0], eb_vel[0], test_1.time);
+      madd_v3_v3v3fl(tmp.eb[1], eb0[1], eb_vel[1], test_1.time);
+      isect_seg_seg_v3(zero, tmp.ea, UNPACK2(tmp.eb), &test_1.wa, &test_1.wb);
+      mul_v3_v3fl(ra, tmp.ea, test_1.wa);
+      interp_v3_v3v3(rb, UNPACK2(tmp.eb), test_1.wb);
+      sub_v3_v3v3(dir_curr, ra, rb);
+      test_1.dist_sq = len_squared_v3(dir_curr);
+      if (test_1.dist_sq <= dist_sq_best) {
+        tmax = test_1;
+        break;
+      }
+
+      mul_v3_v3fl(vel_a, ea_vel, test_1.wa);
+      interp_v3_v3v3(vel_b, UNPACK2(eb_vel), test_1.wb);
+      sub_v3_v3v3(vel_rel, vel_b, vel_a);
+      vel_rel_len_sq = len_squared_v3(vel_rel);
+      move_dist_sq = vel_rel_len_sq * delta_time_sq;
+      if (move_dist_sq < dist_sq_best) {
+        /* The range is already very small. */
+        break;
+      }
+
+      time_offset = dot_v3v3(dir_curr, vel_rel) / vel_rel_len_sq;
+      test_2.time = test_1.time + time_offset;
+      if (test_2.time <= tmin.time) {
+        test_2.time = (test_1.time + tmin.time) * 0.5f;
+      }
+      else if (test_2.time >= tmax.time) {
+        test_2.time = (tmax.time + test_1.time) * 0.5f;
+      }
+
+      madd_v3_v3v3fl(tmp.ea, ea[0], ea_vel, test_2.time);
+      madd_v3_v3v3fl(tmp.eb[0], eb0[0], eb_vel[0], test_2.time);
+      madd_v3_v3v3fl(tmp.eb[1], eb0[1], eb_vel[1], test_2.time);
+      isect_seg_seg_v3(zero, tmp.ea, UNPACK2(tmp.eb), &test_2.wa, &test_2.wb);
+      mul_v3_v3fl(ra, tmp.ea, test_2.wa);
+      interp_v3_v3v3(rb, UNPACK2(tmp.eb), test_2.wb);
+      sub_v3_v3v3(dir_curr, ra, rb);
+      test_2.dist_sq = len_squared_v3(dir_curr);
+      if (test_2.dist_sq <= dist_sq_best) {
+        tmax = test_2;
+        break;
+        // if (test_1.time < max.time) {
+        //  min = test_1;
+        //}
+      }
+
+      if (test_1.time <= test_2.time) {
+        if (test_1.dist_sq < test_2.dist_sq) {
+          /* Moving away. */
+          tmax = test_2;
+        }
+        else {
+          /* Aproaxing. */
+          tmin = test_1;
+        }
+      }
+      else {
+        if (test_2.dist_sq < test_1.dist_sq) {
+          /* Moving away. */
+          tmax = test_1;
+        }
+        else {
+          /* Aproaxing. */
+          tmin = test_2;
+        }
+      }
+
+      /* Calc a new time to test. */
+      mul_v3_v3fl(vel_a, ea_vel, test_2.wa);
+      interp_v3_v3v3(vel_b, UNPACK2(eb_vel), test_2.wb);
+      sub_v3_v3v3(vel_rel, vel_b, vel_a);
+      vel_rel_len_sq = len_squared_v3(vel_rel);
+      move_dist_sq = vel_rel_len_sq * delta_time_sq;
+      if (move_dist_sq < dist_sq_best) {
+        /* The range is already very small. */
+        break;
+      }
+
+      time_offset = dot_v3v3(dir_curr, vel_rel) / vel_rel_len_sq;
+      test_1.time = test_2.time + time_offset;
+    }
+  }
+
+  if (tmax.dist_sq > dist_sq_best) {
+    return false;
+  }
+
+  {
+    /* Last refinement. */
+    mul_v3_v3fl(vel_a, ea_vel, tmax.wa);
+    interp_v3_v3v3(vel_b, UNPACK2(eb_vel), tmax.wb);
+    sub_v3_v3v3(vel_rel, vel_b, vel_a);
+    BLI_assert(!is_zero_v3(vel_rel));
+
+    float vel_rel_lensq_inv = 1.0f / len_squared_v3(vel_rel);
+    if (max_bisect_steps == 0) {
+      madd_v3_v3v3fl(tmp.ea, ea[0], ea_vel, test_2.time);
+      madd_v3_v3v3fl(tmp.eb[0], eb0[0], eb_vel[0], test_2.time);
+      madd_v3_v3v3fl(tmp.eb[1], eb0[1], eb_vel[1], test_2.time);
+      isect_seg_seg_v3(zero, tmp.ea, UNPACK2(tmp.eb), &test_2.wa, &test_2.wb);
+      mul_v3_v3fl(ra, tmp.ea, test_2.wa);
+      interp_v3_v3v3(rb, UNPACK2(tmp.eb), test_2.wb);
+      sub_v3_v3v3(dir_curr, ra, rb);
+    }
+
+    float dir_to_closest[3], p_close[3];
+    float mul = dot_v3v3(dir_curr, vel_rel) * vel_rel_lensq_inv;
+    mul_v3_v3fl(dir_to_closest, vel_rel, mul);
+    sub_v3_v3v3(p_close, dir_curr, dir_to_closest);
+
+    float dvec_len_sq = dist_sq_best - len_squared_v3(p_close);
+    float fac_sq = dvec_len_sq * vel_rel_lensq_inv;
+    float fac = sqrtf(fac_sq);
+
+    float offset_dir[3];
+    mul_v3_v3fl(offset_dir, vel_rel, -fac);
+
+    add_v3_v3v3(dir_curr, dir_to_closest, offset_dir);
+
+    time_offset = dot_v3v3(dir_curr, vel_rel) * vel_rel_lensq_inv;
+    test_1.time = tmax.time + time_offset;
+    if (test_1.time < 0.0f) {
+      test_1.time = 0.0f;
+    }
+
+    madd_v3_v3v3fl(tmp.ea, ea[0], ea_vel, test_1.time);
+    madd_v3_v3v3fl(tmp.eb[0], eb0[0], eb_vel[0], test_1.time);
+    madd_v3_v3v3fl(tmp.eb[1], eb0[1], eb_vel[1], test_1.time);
+    isect_seg_seg_v3(zero, tmp.ea, UNPACK2(tmp.eb), &test_1.wa, &test_1.wb);
+    mul_v3_v3fl(ra, tmp.ea, test_1.wa);
+    interp_v3_v3v3(rb, UNPACK2(tmp.eb), test_1.wb);
+    sub_v3_v3v3(dir_curr, ra, rb);
+    test_1.dist_sq = len_squared_v3(dir_curr);
+  }
+
+  float r_tmp[3];
+  interp_v3_v3v3(r_tmp, seg_a_time0[0], seg_a_time1[0], test_1.time);
+  add_v3_v3v3(r_seg_a, ra, r_tmp);
+  add_v3_v3v3(r_seg_b, rb, r_tmp);
+  *r_dist_sq = test_1.dist_sq;
+  *r_time = test_1.time;
+  return true;
+}

source/blender/blenloader/intern/versioning_290.c

@@ -29,6 +29,7 @@
 #include "DNA_anim_types.h"
 #include "DNA_brush_types.h"
 #include "DNA_cachefile_types.h"
+#include "DNA_cloth_types.h"
 #include "DNA_constraint_types.h"
 #include "DNA_fluid_types.h"
 #include "DNA_genfile.h"
@@ -38,6 +39,7 @@
 #include "DNA_mesh_types.h"
 #include "DNA_meshdata_types.h"
 #include "DNA_modifier_types.h"
+#include "DNA_object_force_types.h"
 #include "DNA_object_types.h"
 #include "DNA_particle_types.h"
 #include "DNA_pointcloud_types.h"
@@ -840,6 +842,21 @@ void blo_do_versions_290(FileData *fd, Library *UNUSED(lib), Main *bmain)
         }
       }
     }
+
+    /* Change cloth friction value from [0, 100] to [0, 1] ranges. */
+    LISTBASE_FOREACH (Object *, ob, &bmain->objects) {
+      LISTBASE_FOREACH (ModifierData *, md, &ob->modifiers) {
+        if (md->type == eModifierType_Cloth) {
+          ClothModifierData *clmd = (ClothModifierData *)md;
+          clmd->coll_parms->self_friction *= 0.01f;
+        }
+      }
+
+      if (ob->pd) {
+        /* Change collision object friction value from [0, 100] to [0, 1] ranges. */
+        ob->pd->pdef_cfrict *= 0.01f;
+      }
+    }
   }
 
   if (!MAIN_VERSION_ATLEAST(bmain, 291, 2)) {

source/blender/editors/sculpt_paint/sculpt_cloth.c

@@ -607,7 +607,7 @@ static ListBase *cloth_brush_collider_cache_create(Depsgraph *depsgraph)
                              DEG_ITER_OBJECT_FLAG_DUPLI) {
     CollisionModifierData *cmd = (CollisionModifierData *)BKE_modifiers_findby_type(
         ob, eModifierType_Collision);
-    if (cmd && cmd->bvhtree) {
+    if (cmd && cmd->bvhtree_tris) {
       if (cache == NULL) {
         cache = MEM_callocN(sizeof(ListBase), "ColliderCache array");
       }
@@ -615,7 +615,8 @@ static ListBase *cloth_brush_collider_cache_create(Depsgraph *depsgraph)
       ColliderCache *col = MEM_callocN(sizeof(ColliderCache), "ColliderCache");
       col->ob = ob;
       col->collmd = cmd;
-      collision_move_object(cmd, 1.0, 0.0, true);
+      BKE_collision_move_object(cmd, 1.0, 0.0);
+      BKE_collision_bvhtree_ensure_and_update(cmd, true, true, false, false);
       BLI_addtail(cache, col);
     }
   }
@@ -688,7 +689,7 @@ static void cloth_brush_solve_collision(Object *object,
     col.col_data = collmd;
     isect_ray_tri_watertight_v3_precalc(&col.isect_precalc, ray_normal);
 
-    BLI_bvhtree_ray_cast_ex(collmd->bvhtree,
+    BLI_bvhtree_ray_cast_ex(collmd->bvhtree_tris,
                             ray_start,
                             ray_normal,
                             0.3f,

source/blender/makesdna/DNA_cloth_types.h

@@ -52,8 +52,7 @@ typedef struct ClothSimSettings {
   float Cvi;
   /** Gravity/external force vector. */
   float gravity[3];
-  /** This is the duration of our time step, computed..   */
+  float _pad0;
-  float dt;
   /** The mass of the entire cloth. */
   float mass;
   /** Structural spring stiffness. */
@@ -116,7 +115,7 @@ typedef struct ClothSimSettings {
    * gradient. */
   float fluid_density;
   short vgroup_pressure;
-  char _pad7[6];
+  char _pad1[6];
 
   /* XXX various hair stuff
    * should really be separate, this struct is a horrible mess already
@@ -172,12 +171,12 @@ typedef struct ClothSimSettings {
   float internal_spring_max_diversion;
   /** Vertex group for scaling structural stiffness. */
   short vgroup_intern;
-  char _pad1[2];
+  char _pad2[2];
   float internal_tension;
   float internal_compression;
   float max_internal_tension;
   float max_internal_compression;
-  char _pad0[4];
+  char _pad3[4];
 
 } ClothSimSettings;
 
@@ -221,12 +220,12 @@ typedef struct ClothCollSettings {
   struct LinkNode *collision_list;
   /** Min distance for collisions. */
   float epsilon;
-  /** Fiction/damping with self contact. */
+  /** Fiction with self contact. */
   float self_friction;
-  /** Friction/damping applied on contact with other object. */
+  /** Friction applied on contact. */
-  float friction;
+  float friction DNA_DEPRECATED;
-  /** Collision restitution on contact with other object. */
+  /** Collision restitution. */
-  float damping;
+  float damping DNA_DEPRECATED;
   /** For selfcollision. */
   float selfepsilon;
   float repel_force DNA_DEPRECATED;
@@ -237,16 +236,19 @@ typedef struct ClothCollSettings {
   short self_loop_count DNA_DEPRECATED;
   /** How many iterations for the collision loop. */
   short loop_count;
-  char _pad[4];
+  /** Which element will collide, vertex, edge or triangle. */
+  char primitive;
+  char primitive_self;
+  char _pad[2];
   /** Only use colliders from this group of objects. */
   struct Collection *group;
   /** Vgroup to paint which vertices are used for self collisions. */
   short vgroup_selfcol;
   char _pad2[6];
   /** Impulse clamp for object collisions. */
-  float clamp;
+  float clamp DNA_DEPRECATED;
   /** Impulse clamp for self collisions. */
-  float self_clamp;
+  float self_clamp DNA_DEPRECATED;
 } ClothCollSettings;
 
 /* COLLISION FLAGS */
@@ -255,6 +257,12 @@ typedef enum {
   CLOTH_COLLSETTINGS_FLAG_SELF = (1 << 2),    /* enables selfcollisions */
 } CLOTH_COLLISIONSETTINGS_FLAGS;
 
+typedef enum {
+  CLOTH_TRIANGLE,
+  CLOTH_EDGE,
+  CLOTH_VERTEX,
+} eClothPrimitive;
+
 #ifdef __cplusplus
 }
 #endif

source/blender/makesdna/DNA_modifier_defaults.h

@@ -107,7 +107,6 @@
     .mingoal = 0.0f, \
     .Cvi = 1.0f, \
     .gravity = {0.0f, 0.0f, -9.81f}, \
-    .dt = 0.0f, \
     .mass = 0.3f, \
     .shear = 5.0f, \
     .bending = 0.5f, \
@@ -171,12 +170,12 @@
   { \
     .collision_list = NULL, \
     .epsilon = 0.015f, \
-    .self_friction = 5.0f, \
+    .self_friction = 0.05f, \
-    .friction = 5.0f, \
+    .friction = 0.05f, \
-    .damping = 0.0f, \
+    .damping = 1.0f, \
     .selfepsilon = 0.015f, \
     .flags = CLOTH_COLLSETTINGS_FLAG_ENABLED, \
-    .loop_count = 2, \
+    .loop_count = 3, \
     .group = NULL, \
     .vgroup_selfcol = 0, \
     .clamp = 0.0f, \
@@ -207,12 +206,16 @@
     .current_x = NULL, \
     .current_v = NULL, \
     .tri = NULL, \
+    .medge = NULL, \
     .mvert_num = 0, \
+    .medge_num = 0, \
     .tri_num = 0, \
     .time_x = -1000.0f, \
     .time_xnew = -1000.0f, \
     .is_static = false, \
-    .bvhtree = NULL, \
+    .bvhtree_tris = NULL, \
+    .bvhtree_edges = NULL, \
+    .bvhtree_verts = NULL, \
   }
 
 #define _DNA_DEFAULT_CorrectiveSmoothModifierData \

source/blender/makesdna/DNA_modifier_types.h

@@ -839,18 +839,25 @@ typedef struct CollisionModifierData {
   /** (xnew - x) at the actual inter-frame step. */
   struct MVert *current_v;
 
+  struct MEdge *medge;
   struct MVertTri *tri;
 
   unsigned int mvert_num;
+  unsigned int medge_num;
   unsigned int tri_num;
+  unsigned int _pad;
+
   /** Cfra time of modifier. */
   float time_x, time_xnew;
+
   /** Collider doesn't move this frame, i.e. x[].co==xnew[].co. */
   char is_static;
-  char _pad[7];
+  char _pad2[7];
 
-  /** Bounding volume hierarchy for this cloth object. */
+  /** Bounding volume hierarchy for this object. */
-  struct BVHTree *bvhtree;
+  struct BVHTree *bvhtree_tris;
+  struct BVHTree *bvhtree_edges;
+  struct BVHTree *bvhtree_verts;
 } CollisionModifierData;
 
 typedef struct SurfaceModifierData {

source/blender/makesdna/DNA_object_force_types.h

@@ -349,8 +349,7 @@ typedef struct SoftBody {
 #define PFIELD_SMOKE_DENSITY (1 << 17)
 /** used for (simple) force */
 #define PFIELD_GRAVITATION (1 << 18)
-/** Enable cloth collision side detection based on normal. */
+// #define PFIELD_CLOTH_USE_CULLING (1 << 19)  /* deprecated */
-#define PFIELD_CLOTH_USE_CULLING (1 << 19)
 /** Replace collision direction with collider normal. */
 #define PFIELD_CLOTH_USE_NORMAL (1 << 20)
 

source/blender/makesrna/intern/rna_cloth.c

@@ -1071,6 +1071,13 @@ static void rna_def_cloth_collision_settings(BlenderRNA *brna)
   StructRNA *srna;
   PropertyRNA *prop;
 
+  static const EnumPropertyItem prop_primitive_items[] = {
+      {CLOTH_VERTEX, "VERTEX", 0, "Vertex", "Collide only the vertices of the cloth"},
+      {CLOTH_EDGE, "EDGE", 0, "Edge", "Collide only the edges of the cloth"},
+      {CLOTH_TRIANGLE, "TRIANGLE", 0, "Triangle", "Collide only the triangles of the cloth"},
+      {0, NULL, 0, NULL, NULL},
+  };
+
   srna = RNA_def_struct(brna, "ClothCollisionSettings", NULL);
   RNA_def_struct_ui_text(
       srna,
@@ -1095,19 +1102,6 @@ static void rna_def_cloth_collision_settings(BlenderRNA *brna)
       "Minimum distance between collision objects before collision response takes effect");
   RNA_def_property_update(prop, 0, "rna_cloth_update");
 
-  prop = RNA_def_property(srna, "friction", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_range(prop, 0.0f, 80.0f);
-  RNA_def_property_ui_text(
-      prop, "Friction", "Friction force if a collision happened (higher = less movement)");
-  RNA_def_property_update(prop, 0, "rna_cloth_update");
-
-  prop = RNA_def_property(srna, "damping", PROP_FLOAT, PROP_FACTOR);
-  RNA_def_property_float_sdna(prop, NULL, "damping");
-  RNA_def_property_range(prop, 0.0f, 1.0f);
-  RNA_def_property_float_default(prop, 1.0f);
-  RNA_def_property_ui_text(prop, "Restitution", "Amount of velocity lost on collision");
-  RNA_def_property_update(prop, 0, "rna_cloth_update");
-
   prop = RNA_def_property(srna, "collision_quality", PROP_INT, PROP_NONE);
   RNA_def_property_int_sdna(prop, NULL, "loop_count");
   RNA_def_property_range(prop, 1, SHRT_MAX);
@@ -1118,13 +1112,11 @@ static void rna_def_cloth_collision_settings(BlenderRNA *brna)
       "How many collision iterations should be done. (higher is better quality but slower)");
   RNA_def_property_update(prop, 0, "rna_cloth_update");
 
-  prop = RNA_def_property(srna, "impulse_clamp", PROP_FLOAT, PROP_NONE);
+  prop = RNA_def_property(srna, "primitive", PROP_ENUM, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "clamp");
+  RNA_def_property_enum_sdna(prop, NULL, "primitive");
-  RNA_def_property_range(prop, 0.0f, 100.0f);
+  RNA_def_property_enum_items(prop, prop_primitive_items);
   RNA_def_property_ui_text(
-      prop,
+      prop, "Primitive", "Geometry of the cloth that collides with the object");
-      "Impulse Clamping",
-      "Clamp collision impulses to avoid instability (0.0 to disable clamping)");
   RNA_def_property_update(prop, 0, "rna_cloth_update");
 
   /* self collision */
@@ -1143,8 +1135,8 @@ static void rna_def_cloth_collision_settings(BlenderRNA *brna)
       "Minimum distance between cloth faces before collision response takes effect");
   RNA_def_property_update(prop, 0, "rna_cloth_update");
 
-  prop = RNA_def_property(srna, "self_friction", PROP_FLOAT, PROP_NONE);
+  prop = RNA_def_property(srna, "self_friction", PROP_FLOAT, PROP_FACTOR);
-  RNA_def_property_range(prop, 0.0f, 80.0f);
+  RNA_def_property_range(prop, 0.0f, 1.0f);
   RNA_def_property_ui_text(prop, "Self Friction", "Friction with self contact");
   RNA_def_property_update(prop, 0, "rna_cloth_update");
 
@@ -1166,13 +1158,10 @@ static void rna_def_cloth_collision_settings(BlenderRNA *brna)
       "Vertex group to define vertices which are not used during self collisions");
   RNA_def_property_update(prop, 0, "rna_cloth_update");
 
-  prop = RNA_def_property(srna, "self_impulse_clamp", PROP_FLOAT, PROP_NONE);
+  prop = RNA_def_property(srna, "self_primitive", PROP_ENUM, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "self_clamp");
+  RNA_def_property_enum_sdna(prop, NULL, "primitive_self");
-  RNA_def_property_range(prop, 0.0f, 100.0f);
+  RNA_def_property_enum_items(prop, prop_primitive_items);
-  RNA_def_property_ui_text(
+  RNA_def_property_ui_text(prop, "Primitive", "Geometry of the cloth that collided with itself");
-      prop,
-      "Impulse Clamping",
-      "Clamp collision impulses to avoid instability (0.0 to disable clamping)");
   RNA_def_property_update(prop, 0, "rna_cloth_update");
 }
 

source/blender/makesrna/intern/rna_object_force.c

@@ -1198,20 +1198,12 @@ static void rna_def_collision(BlenderRNA *brna)
       "How much of effector force gets lost during collision with this object (in percent)");
   RNA_def_property_update(prop, 0, "rna_CollisionSettings_update");
 
-  prop = RNA_def_property(srna, "cloth_friction", PROP_FLOAT, PROP_NONE);
+  prop = RNA_def_property(srna, "cloth_friction", PROP_FLOAT, PROP_FACTOR);
   RNA_def_property_float_sdna(prop, NULL, "pdef_cfrict");
-  RNA_def_property_range(prop, 0.0f, 80.0f);
+  RNA_def_property_range(prop, 0.0f, 1.0f);
   RNA_def_property_ui_text(prop, "Friction", "Friction for cloth collisions");
   RNA_def_property_update(prop, 0, "rna_CollisionSettings_update");
 
-  prop = RNA_def_property(srna, "use_culling", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_sdna(prop, NULL, "flag", PFIELD_CLOTH_USE_CULLING);
-  RNA_def_property_ui_text(
-      prop,
-      "Single Sided",
-      "Cloth collision acts with respect to the collider normals (improves penetration recovery)");
-  RNA_def_property_update(prop, 0, "rna_CollisionSettings_update");
-
   prop = RNA_def_property(srna, "use_normal", PROP_BOOLEAN, PROP_NONE);
   RNA_def_property_boolean_sdna(prop, NULL, "flag", PFIELD_CLOTH_USE_NORMAL);
   RNA_def_property_ui_text(prop,

source/blender/modifiers/intern/MOD_collision.c

@@ -73,9 +73,17 @@ static void freeData(ModifierData *md)
   CollisionModifierData *collmd = (CollisionModifierData *)md;
 
   if (collmd) { /* Seriously? */
-    if (collmd->bvhtree) {
+    if (collmd->bvhtree_tris) {
-      BLI_bvhtree_free(collmd->bvhtree);
+      BLI_bvhtree_free(collmd->bvhtree_tris);
-      collmd->bvhtree = NULL;
+      collmd->bvhtree_tris = NULL;
+    }
+    if (collmd->bvhtree_edges) {
+      BLI_bvhtree_free(collmd->bvhtree_edges);
+      collmd->bvhtree_edges = NULL;
+    }
+    if (collmd->bvhtree_verts) {
+      BLI_bvhtree_free(collmd->bvhtree_verts);
+      collmd->bvhtree_verts = NULL;
     }
 
     MEM_SAFE_FREE(collmd->x);
@@ -84,10 +92,12 @@ static void freeData(ModifierData *md)
     MEM_SAFE_FREE(collmd->current_xnew);
     MEM_SAFE_FREE(collmd->current_v);
 
+    MEM_SAFE_FREE(collmd->medge);
     MEM_SAFE_FREE(collmd->tri);
 
     collmd->time_x = collmd->time_xnew = -1000;
     collmd->mvert_num = 0;
+    collmd->medge_num = 0;
     collmd->tri_num = 0;
     collmd->is_static = false;
   }
@@ -177,8 +187,14 @@ static void deformVerts(ModifierData *md,
         collmd->tri = tri;
       }
 
+      {
+        /* TODO: This is only required for cloths. */
+        collmd->medge_num = mesh_src->totedge;
+        collmd->medge = MEM_dupallocN(mesh_src->medge);
+      }
+
       /* create bounding box hierarchy */
-      collmd->bvhtree = bvhtree_build_from_mvert(
+      collmd->bvhtree_tris = BKE_collision_bvhtree_tris_build(
           collmd->x, collmd->tri, collmd->tri_num, ob->pd->pdef_sboft);
 
       collmd->time_x = collmd->time_xnew = current_time;
@@ -207,27 +223,21 @@ static void deformVerts(ModifierData *md,
       memcpy(collmd->current_x, collmd->x, mvert_num * sizeof(MVert));
 
       /* check if GUI setting has changed for bvh */
-      if (collmd->bvhtree) {
+      if (collmd->bvhtree_tris) {
-        if (ob->pd->pdef_sboft != BLI_bvhtree_get_epsilon(collmd->bvhtree)) {
+        if (ob->pd->pdef_sboft != BLI_bvhtree_get_epsilon(collmd->bvhtree_tris)) {
-          BLI_bvhtree_free(collmd->bvhtree);
+          BLI_bvhtree_free(collmd->bvhtree_tris);
-          collmd->bvhtree = bvhtree_build_from_mvert(
+          collmd->bvhtree_tris = NULL;
-              collmd->current_x, collmd->tri, collmd->tri_num, ob->pd->pdef_sboft);
         }
       }
 
       /* happens on file load (ONLY when i decomment changes in readfile.c) */
-      if (!collmd->bvhtree) {
+      if (!collmd->bvhtree_tris) {
-        collmd->bvhtree = bvhtree_build_from_mvert(
+        collmd->bvhtree_tris = BKE_collision_bvhtree_tris_build(
             collmd->current_x, collmd->tri, collmd->tri_num, ob->pd->pdef_sboft);
       }
       else if (!collmd->is_static || !is_static) {
         /* recalc static bounding boxes */
-        bvhtree_update_from_mvert(collmd->bvhtree,
+        BKE_collision_bvhtree_ensure_and_update(collmd, true, true, false, false);
-                                  collmd->current_x,
-                                  collmd->current_xnew,
-                                  collmd->tri,
-                                  collmd->tri_num,
-                                  true);
       }
 
       collmd->is_static = is_static;
@@ -286,10 +296,14 @@ static void blendRead(BlendDataReader *UNUSED(reader), ModifierData *md)
   collmd->current_v = NULL;
   collmd->time_x = collmd->time_xnew = -1000;
   collmd->mvert_num = 0;
+  collmd->medge_num = 0;
   collmd->tri_num = 0;
   collmd->is_static = false;
-  collmd->bvhtree = NULL;
+  collmd->bvhtree_tris = NULL;
+  collmd->bvhtree_edges = NULL;
+  collmd->bvhtree_verts = NULL;
   collmd->tri = NULL;
+  collmd->medge = NULL;
 }
 
 ModifierTypeInfo modifierType_Collision = {

source/blender/python/mathutils/mathutils_geometry.c

@@ -1675,6 +1675,199 @@ exit_cdt:
 
 #endif /* MATH_STANDALONE */
 
+/* --------------------------------- COLLISION FUNCTIONS-------------------- */
+
+PyDoc_STRVAR(
+    M_Geometry_collision_tri_point_bisect_doc,
+    ".. function:: collision_tri_point_bisect(tri_t0, tri_t1, p_t0, p_t1, coll_dist, max_steps)\n"
+    "\n"
+    "   Returns a list of points inside all planes given and a list of index values for "
+    "the planes used.\n"
+    "\n"
+    "   :arg planes: List of planes (4D vectors).\n"
+    "   :type planes: list of :class:`mathutils.Vector`\n"
+    "   :return: two lists, once containing the vertices inside the planes, another "
+    "containing the plane indices used\n"
+    "   :rtype: pair of lists\n");
+static PyObject *M_Geometry_collision_tri_point_bisect(PyObject *UNUSED(self), PyObject *args)
+{
+  const char *error_prefix = "collision_tri_point_bisect";
+  PyObject *py_tri_t[2], *py_p_t[2];
+  float tri_t[2][3][3], p_t[2][3], coll_dist;
+  uint max_steps;
+
+  if (!PyArg_ParseTuple(args,
+                        "OOOOfI:collision_tri_point_bisect",
+                        &py_tri_t[0],
+                        &py_tri_t[1],
+                        &py_p_t[0],
+                        &py_p_t[1],
+                        &coll_dist,
+                        &max_steps)) {
+    return NULL;
+  }
+
+  for (int i = 0; i < 2; i++) {
+    /* non list/tuple cases */
+    PyObject *value_fast = NULL;
+    if (!(value_fast = PySequence_Fast(py_tri_t[i], error_prefix))) {
+      /* PySequence_Fast sets the error */
+      return NULL;
+    }
+
+    if (PySequence_Fast_GET_SIZE(value_fast) != 3) {
+      Py_DECREF(value_fast);
+      PyErr_SetString(PyExc_ValueError,
+                      "matrix[begin:end] = []: "
+                      "size mismatch in slice assignment");
+      return NULL;
+    }
+
+    PyObject **value_fast_items = PySequence_Fast_ITEMS(value_fast);
+    int size = (int)PySequence_Fast_GET_SIZE(value_fast);
+    if (size != 3) {
+      PyErr_Format(
+          PyExc_ValueError, "%.200s: sequence size is %d, expected %d", error_prefix, size, 3);
+      Py_DECREF(value_fast);
+      return NULL;
+    }
+
+    /* parse sub items */
+    for (int j = 0; j < 3; j++) {
+      /* parse each sub sequence */
+      PyObject *item = value_fast_items[j];
+
+      if (mathutils_array_parse(tri_t[i][j], 3, 3, item, error_prefix) == -1) {
+        Py_DECREF(value_fast);
+        return NULL;
+      }
+    }
+
+    Py_DECREF(value_fast);
+
+    if (mathutils_array_parse(p_t[i], 3, 3, py_p_t[i], error_prefix) == -1) {
+      return NULL;
+    }
+  }
+
+  PyObject *py_retval = PyTuple_New(4);
+  float r_time = -1.0f, r_dist_sq = -1.0f, r_tri[3] = {0.0f}, r_point[3] = {0.0f};
+  if (collision_tri_point_bisect(UNPACK2(tri_t),
+                                 UNPACK2(p_t),
+                                 coll_dist,
+                                 max_steps,
+                                 &r_time,
+                                 &r_dist_sq,
+                                 r_tri,
+                                 r_point)) {
+    PyTuple_SET_ITEMS(py_retval,
+                      PyFloat_FromDouble(r_time),
+                      PyFloat_FromDouble(r_dist_sq),
+                      Vector_CreatePyObject(r_tri, 3, NULL),
+                      Vector_CreatePyObject(r_point, 3, NULL));
+  }
+  else {
+    PyC_Tuple_Fill(py_retval, Py_None);
+  }
+
+  return py_retval;
+}
+
+PyDoc_STRVAR(
+    M_Geometry_collision_seg_seg_bisect_doc,
+    ".. function:: collision_seg_seg_bisect(sa_t0, sa_t1, sb_t0, sb_t1, coll_dist, max_steps)\n"
+    "\n"
+    "   Returns a list of points inside all planes given and a list of index values for "
+    "the planes used.\n"
+    "\n"
+    "   :arg planes: List of planes (4D vectors).\n"
+    "   :type planes: list of :class:`mathutils.Vector`\n"
+    "   :return: two lists, once containing the vertices inside the planes, another "
+    "containing the plane indices used\n"
+    "   :rtype: pair of lists\n");
+static PyObject *M_Geometry_collision_seg_seg_bisect(PyObject *UNUSED(self), PyObject *args)
+{
+  const char *error_prefix = "collision_seg_seg_bisect";
+  PyObject *py_s_t[4];
+  float s_t[4][2][3], coll_dist;
+  uint max_steps;
+
+  if (!PyArg_ParseTuple(args,
+                        "OOOOfI:collision_seg_seg_bisect",
+                        &py_s_t[0],
+                        &py_s_t[1],
+                        &py_s_t[2],
+                        &py_s_t[3],
+                        &coll_dist,
+                        &max_steps)) {
+    return NULL;
+  }
+
+  for (int i = 0; i < 4; i++) {
+    /* non list/tuple cases */
+    PyObject *value_fast = NULL;
+    if (!(value_fast = PySequence_Fast(py_s_t[i], error_prefix))) {
+      /* PySequence_Fast sets the error */
+      return NULL;
+    }
+
+    if (PySequence_Fast_GET_SIZE(value_fast) != 2) {
+      Py_DECREF(value_fast);
+      PyErr_SetString(PyExc_ValueError,
+                      "matrix[begin:end] = []: "
+                      "size mismatch in slice assignment");
+      return NULL;
+    }
+
+    PyObject **value_fast_items = PySequence_Fast_ITEMS(value_fast);
+    int size = (int)PySequence_Fast_GET_SIZE(value_fast);
+    if (size != 2) {
+      PyErr_Format(
+          PyExc_ValueError, "%.200s: sequence size is %d, expected %d", error_prefix, size, 3);
+      Py_DECREF(value_fast);
+      return NULL;
+    }
+
+    /* parse sub items */
+    for (int j = 0; j < 2; j++) {
+      /* parse each sub sequence */
+      PyObject *item = value_fast_items[j];
+
+      if (mathutils_array_parse(s_t[i][j], 3, 3, item, error_prefix) == -1) {
+        Py_DECREF(value_fast);
+        return NULL;
+      }
+    }
+
+    Py_DECREF(value_fast);
+  }
+
+  PyObject *py_retval = PyTuple_New(4);
+  float r_time = -1.0f, r_dist_sq = -1.0f, r_tri[3] = {0.0f}, r_point[3] = {0.0f};
+  if (collision_seg_seg_bisect(s_t[0],
+                               s_t[1],
+                               s_t[2],
+                               s_t[3],
+                               coll_dist,
+                               max_steps,
+                               &r_time,
+                               &r_dist_sq,
+                               r_tri,
+                               r_point)) {
+
+    PyTuple_SET_ITEMS(py_retval,
+                      PyFloat_FromDouble(r_time),
+                      PyFloat_FromDouble(r_dist_sq),
+                      Vector_CreatePyObject(r_tri, 3, NULL),
+                      Vector_CreatePyObject(r_point, 3, NULL));
+  }
+  else {
+    PyC_Tuple_Fill(py_retval, Py_None);
+  }
+
+  return py_retval;
+}
+
 static PyMethodDef M_Geometry_methods[] = {
     {"intersect_ray_tri",
      (PyCFunction)M_Geometry_intersect_ray_tri,
@@ -1770,6 +1963,14 @@ static PyMethodDef M_Geometry_methods[] = {
     {"box_fit_2d", (PyCFunction)M_Geometry_box_fit_2d, METH_O, M_Geometry_box_fit_2d_doc},
     {"box_pack_2d", (PyCFunction)M_Geometry_box_pack_2d, METH_O, M_Geometry_box_pack_2d_doc},
 #endif
+    {"collision_tri_point_bisect",
+     (PyCFunction)M_Geometry_collision_tri_point_bisect,
+     METH_VARARGS,
+     M_Geometry_collision_tri_point_bisect_doc},
+    {"collision_seg_seg_bisect",
+     (PyCFunction)M_Geometry_collision_seg_seg_bisect,
+     METH_VARARGS,
+     M_Geometry_collision_seg_seg_bisect_doc},
     {NULL, NULL, 0, NULL},
 };
 

source/blender/simulation/CMakeLists.txt

@@ -40,7 +40,6 @@ set(SRC
   intern/SIM_mass_spring.cpp
   intern/hair_volume.cpp
   intern/implicit_blender.c
-  intern/implicit_eigen.cpp
 
   intern/ConstrainedConjugateGradient.h
   intern/eigen_utils.h

source/blender/simulation/intern/SIM_mass_spring.cpp

@@ -35,7 +35,7 @@
 #include "BLI_utildefines.h"
 
 #include "BKE_cloth.h"
-#include "BKE_collision.h"
+#include "BKE_cloth_collision.h"
 #include "BKE_effect.h"
 
 #include "SIM_mass_spring.h"
@@ -123,7 +123,7 @@ static float cloth_calc_volume(ClothModifierData *clmd)
     return 0.0f;
   }
 
-  for (unsigned int i = 0; i < cloth->primitive_num; i++) {
+  for (unsigned int i = 0; i < cloth->mlooptri_num; i++) {
     const MVertTri *vt = &tri[i];
 
     if (cloth_get_pressure_weights(clmd, vt, weights)) {
@@ -151,7 +151,7 @@ static float cloth_calc_rest_volume(ClothModifierData *clmd)
     return 0.0f;
   }
 
-  for (unsigned int i = 0; i < cloth->primitive_num; i++) {
+  for (unsigned int i = 0; i < cloth->mlooptri_num; i++) {
     const MVertTri *vt = &tri[i];
 
     if (cloth_get_pressure_weights(clmd, vt, weights)) {
@@ -175,7 +175,7 @@ static float cloth_calc_average_pressure(ClothModifierData *clmd, const float *v
   float total_force = 0;
   float total_area = 0;
 
-  for (unsigned int i = 0; i < cloth->primitive_num; i++) {
+  for (unsigned int i = 0; i < cloth->mlooptri_num; i++) {
     const MVertTri *vt = &tri[i];
 
     if (cloth_get_pressure_weights(clmd, vt, weights)) {
@@ -272,8 +272,6 @@ static void cloth_setup_constraints(ClothModifierData *clmd)
       /* pinned vertex constraints */
       SIM_mass_spring_add_constraint_ndof0(data, v, ZERO); /* velocity is defined externally */
     }
-
-    verts[v].impulse_count = 0;
   }
 }
 
@@ -689,7 +687,7 @@ static void cloth_calc_force(
     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++) {
+      for (i = 0; i < cloth->mlooptri_num; i++) {
         const MVertTri *vt = &tri[i];
 
         if (cloth_get_pressure_weights(clmd, vt, weights)) {
@@ -711,7 +709,7 @@ static void cloth_calc_force(
 
   /* handle external forces like wind */
   if (effectors) {
-    bool is_not_hair = (clmd->hairdata == nullptr) && (cloth->primitive_num > 0);
+    bool is_not_hair = (clmd->hairdata == nullptr) && (cloth->mlooptri_num > 0);
     bool has_wind = false, has_force = false;
 
     /* cache per-vertex forces to avoid redundant calculation */
@@ -739,7 +737,7 @@ static void cloth_calc_force(
 
     /* Hair has only edges. */
     if (is_not_hair) {
-      for (i = 0; i < cloth->primitive_num; i++) {
+      for (i = 0; i < cloth->mlooptri_num; i++) {
         const MVertTri *vt = &tri[i];
         if (has_wind) {
           SIM_mass_spring_force_face_wind(data, vt->tri[0], vt->tri[1], vt->tri[2], winvec);
@@ -1087,66 +1085,6 @@ static void cloth_continuum_step(ClothModifierData *clmd, float dt)
   }
 }
 
-#if 0
-static void cloth_calc_volume_force(ClothModifierData *clmd)
-{
-  ClothSimSettings *parms = clmd->sim_parms;
-  Cloth *cloth = clmd->clothObject;
-  Implicit_Data *data = cloth->implicit;
-  int mvert_num = cloth->mvert_num;
-  ClothVertex *vert;
-
-  /* 2.0f is an experimental value that seems to give good results */
-  float smoothfac = 2.0f * parms->velocity_smooth;
-  float collfac = 2.0f * parms->collider_friction;
-  float pressfac = parms->pressure;
-  float minpress = parms->pressure_threshold;
-  float gmin[3], gmax[3];
-  int i;
-
-  hair_get_boundbox(clmd, gmin, gmax);
-
-  /* gather velocities & density */
-  if (smoothfac > 0.0f || pressfac > 0.0f) {
-    HairVertexGrid *vertex_grid = SIM_hair_volume_create_vertex_grid(
-        clmd->sim_parms->voxel_res, gmin, gmax);
-
-    vert = cloth->verts;
-    for (i = 0; i < mvert_num; i++, vert++) {
-      float x[3], v[3];
-
-      if (vert->solver_index < 0) {
-        copy_v3_v3(x, vert->x);
-        copy_v3_v3(v, vert->v);
-      }
-      else {
-        SIM_mass_spring_get_motion_state(data, vert->solver_index, x, v);
-      }
-      SIM_hair_volume_add_vertex(vertex_grid, x, v);
-    }
-    SIM_hair_volume_normalize_vertex_grid(vertex_grid);
-
-    vert = cloth->verts;
-    for (i = 0; i < mvert_num; i++, vert++) {
-      float x[3], v[3], f[3], dfdx[3][3], dfdv[3][3];
-
-      if (vert->solver_index < 0) {
-        continue;
-      }
-
-      /* calculate volumetric forces */
-      SIM_mass_spring_get_motion_state(data, vert->solver_index, x, v);
-      SIM_hair_volume_vertex_grid_forces(
-          vertex_grid, x, v, smoothfac, pressfac, minpress, f, dfdx, dfdv);
-      /* apply on hair data */
-      SIM_mass_spring_force_extern(data, vert->solver_index, f, dfdx, dfdv);
-    }
-
-    SIM_hair_volume_free_vertex_grid(vertex_grid);
-  }
-}
-#endif
-
 static void cloth_calc_average_acceleration(ClothModifierData *clmd, float dt)
 {
   Cloth *cloth = clmd->clothObject;
@@ -1178,41 +1116,36 @@ static void cloth_solve_collisions(
   Implicit_Data *id = cloth->implicit;
   ClothVertex *verts = cloth->verts;
   int mvert_num = cloth->mvert_num;
-  const float time_multiplier = 1.0f / (clmd->sim_parms->dt * clmd->sim_parms->timescale);
+  const float time_multiplier = 1.0f / dt;
-  int i;
+  const bool use_pin = (clmd->sim_parms->vgroup_mass > 0);
 
   if (!(clmd->coll_parms->flags &
         (CLOTH_COLLSETTINGS_FLAG_ENABLED | CLOTH_COLLSETTINGS_FLAG_SELF))) {
     return;
   }
 
-  if (!clmd->clothObject->bvhtree) {
+  if (!clmd->clothObject->collision_data) {
     return;
   }
 
-  SIM_mass_spring_solve_positions(id, dt);
-
   /* Update verts to current positions. */
-  for (i = 0; i < mvert_num; i++) {
+  for (int i = 0; i < mvert_num; i++) {
     SIM_mass_spring_get_new_position(id, i, verts[i].tx);
-
+    SIM_mass_spring_get_new_velocity(id, i, verts[i].tv);
-    sub_v3_v3v3(verts[i].tv, verts[i].tx, verts[i].txold);
-    zero_v3(verts[i].dcvel);
   }
 
-  if (cloth_bvh_collision(depsgraph,
+  if (BKE_cloth_collision_do_step(depsgraph,
-                          ob,
+                                  ob,
-                          clmd,
+                                  clmd,
-                          step / clmd->sim_parms->timescale,
+                                  step / clmd->sim_parms->timescale,
-                          dt / clmd->sim_parms->timescale)) {
+                                  dt / clmd->sim_parms->timescale)) {
-    for (i = 0; i < mvert_num; i++) {
+    for (int i = 0; i < mvert_num; i++) {
-      if ((clmd->sim_parms->vgroup_mass > 0) && (verts[i].flags & CLOTH_VERT_FLAG_PINNED)) {
+      if (use_pin && (verts[i].flags & CLOTH_VERT_FLAG_PINNED)) {
         continue;
       }
 
-      SIM_mass_spring_get_new_velocity(id, i, verts[i].tv);
-      madd_v3_v3fl(verts[i].tv, verts[i].dcvel, time_multiplier);
       SIM_mass_spring_set_new_velocity(id, i, verts[i].tv);
+      SIM_mass_spring_set_new_position(id, i, verts[i].tx);
     }
   }
 }
@@ -1272,7 +1205,7 @@ int SIM_cloth_solve(
   Cloth *cloth = clmd->clothObject;
   ClothVertex *verts = cloth->verts /*, *cv*/;
   unsigned int mvert_num = cloth->mvert_num;
-  float dt = clmd->sim_parms->dt * clmd->sim_parms->timescale;
+  const float dt = clmd->sim_parms->timescale / clmd->sim_parms->stepsPerFrame;
   Implicit_Data *id = cloth->implicit;
 
   /* Hydrostatic pressure gradient of the fluid inside the object is affected by acceleration. */
@@ -1321,9 +1254,6 @@ int SIM_cloth_solve(
     SIM_mass_spring_solve_velocities(id, dt, &result);
     cloth_record_result(clmd, &result, dt);
 
-    /* Calculate collision impulses. */
-    cloth_solve_collisions(depsgraph, ob, clmd, step, dt);
-
     if (is_hair) {
       cloth_continuum_step(clmd, dt);
     }
@@ -1333,6 +1263,10 @@ int SIM_cloth_solve(
     }
 
     SIM_mass_spring_solve_positions(id, dt);
+
+    /* Calculate collision impulses. */
+    cloth_solve_collisions(depsgraph, ob, clmd, step, dt);
+
     SIM_mass_spring_apply_result(id);
 
     /* move pinned verts to correct position */

source/blender/simulation/intern/implicit.h

@@ -108,13 +108,6 @@ void SIM_mass_spring_apply_result(struct Implicit_Data *data);
 
 /* Clear the force vector at the beginning of the time step */
 void SIM_mass_spring_clear_forces(struct Implicit_Data *data);
-/* Fictitious forces introduced by moving coordinate systems */
-void SIM_mass_spring_force_reference_frame(struct Implicit_Data *data,
-                                           int index,
-                                           const float acceleration[3],
-                                           const float omega[3],
-                                           const float domega_dt[3],
-                                           float mass);
 /* Simple uniform gravity force */
 void SIM_mass_spring_force_gravity(struct Implicit_Data *data,
                                    int index,
@@ -122,9 +115,6 @@ void SIM_mass_spring_force_gravity(struct Implicit_Data *data,
                                    const float g[3]);
 /* Global drag force (velocity damping) */
 void SIM_mass_spring_force_drag(struct Implicit_Data *data, float drag);
-/* Custom external force */
-void SIM_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 (only generates pressure from the normal component) */
 void SIM_mass_spring_force_face_wind(
     struct Implicit_Data *data, int v1, int v2, int v3, const float (*winvec)[3]);