Added unit tests for testing different phases of collision

release/scripts/startup/bl_ui/properties_physics_softbody.py

@@ -37,7 +37,8 @@ class PHYSICS_PT_softbody(PhysicButtonsPanel, Panel):
         md = context.soft_body
         softbody = md.settings
 
-        layout.prop(softbody, "collision_collection")
+        # layout.prop(softbody, "collision_collection")
+
 
 
 class PHYSICS_PT_softbody_object(PhysicButtonsPanel, Panel):
@@ -80,10 +81,24 @@ class PHYSICS_PT_softbody_simulation(PhysicButtonsPanel, Panel):
 
         md = context.soft_body
         softbody = md.settings
+        ob = context.object
 
         layout.enabled = softbody_panel_enabled(md)
+        flow = layout.grid_flow(row_major=True, columns=0, even_columns=True, even_rows=False, align=True)
 
-        layout.prop(softbody, "speed")
+        col = flow.column()
+        col.prop(softbody, "dt")
+
+        col = flow.column()
+        col.prop(softbody, "substep_count")
+
+        col = flow.column()
+        col.prop(softbody, "alpha_vol")
+
+        col = flow.column()
+        col.prop(softbody, "alpha_edge")
+
+        # layout.prop(softbody, "speed")
 
 
 class PHYSICS_PT_softbody_cache(PhysicButtonsPanel, Panel):
@@ -382,20 +397,20 @@ class PHYSICS_PT_softbody_field_weights(PhysicButtonsPanel, Panel):
 
 classes = (
     PHYSICS_PT_softbody,
-    PHYSICS_PT_softbody_object,
+    # PHYSICS_PT_softbody_object,
     PHYSICS_PT_softbody_simulation,
-    PHYSICS_PT_softbody_cache,
-    PHYSICS_PT_softbody_goal,
-    PHYSICS_PT_softbody_goal_settings,
-    PHYSICS_PT_softbody_goal_strengths,
-    PHYSICS_PT_softbody_edge,
-    PHYSICS_PT_softbody_edge_aerodynamics,
-    PHYSICS_PT_softbody_edge_stiffness,
-    PHYSICS_PT_softbody_collision,
-    PHYSICS_PT_softbody_solver,
-    PHYSICS_PT_softbody_solver_diagnostics,
-    PHYSICS_PT_softbody_solver_helpers,
-    PHYSICS_PT_softbody_field_weights,
+    # PHYSICS_PT_softbody_cache,
+    # PHYSICS_PT_softbody_goal,
+    # PHYSICS_PT_softbody_goal_settings,
+    # PHYSICS_PT_softbody_goal_strengths,
+    # PHYSICS_PT_softbody_edge,
+    # PHYSICS_PT_softbody_edge_aerodynamics,
+    # PHYSICS_PT_softbody_edge_stiffness,
+    # PHYSICS_PT_softbody_collision,
+    # PHYSICS_PT_softbody_solver,
+    # PHYSICS_PT_softbody_solver_diagnostics,
+    # PHYSICS_PT_softbody_solver_helpers,
+    # PHYSICS_PT_softbody_field_weights,
 )
 
 

source/blender/blenkernel/BKE_softbody copy.h

@@ -0,0 +1,86 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright Blender Foundation. All rights reserved. */
+#pragma once
+
+/** \file
+ * \ingroup bke
+ */
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct Depsgraph;
+struct Object;
+struct Scene;
+struct SoftBody;
+
+typedef struct BodyPoint {
+  float origS[3], origE[3], origT[3], pos[3], vec[3], force[3];
+  float goal;
+  float prevpos[3], prevvec[3], prevdx[3], prevdv[3]; /* used for Heun integration */
+  float impdv[3], impdx[3];
+  int nofsprings;
+  int *springs;
+  float choke, choke2, frozen;
+  float colball;
+  short loc_flag; /* reserved by locale module specific states */
+  // char octantflag;
+  float mass;
+  float springweight;
+} BodyPoint;
+
+/**
+ * Allocates and initializes general main data.
+ */
+extern struct SoftBody *sbNew(void);
+
+/**
+ * Frees internal data and soft-body itself.
+ */
+extern void sbFree(struct Object *ob);
+
+/**
+ * Frees simulation data to reset simulation.
+ */
+extern void sbFreeSimulation(struct SoftBody *sb);
+
+/**
+ * Do one simulation step, reading and writing vertex locs from given array.
+ * */
+extern void sbObjectStep(struct Depsgraph *depsgraph,
+                         struct Scene *scene,
+                         struct Object *ob,
+                         float cfra,
+                         float (*vertexCos)[3],
+                         int numVerts);
+
+/**
+ * Makes totally fresh start situation, resets time.
+ */
+extern void sbObjectToSoftbody(struct Object *ob);
+
+/**
+ * Soft-body global visible functions.
+ * Links the soft-body module to a 'test for Interrupt' function, pass NULL to clear the callback.
+ */
+extern void sbSetInterruptCallBack(int (*f)(void));
+
+/**
+ * A precise position vector denoting the motion of the center of mass give a rotation/scale matrix
+ * using averaging method, that's why estimate and not calculate see: this is kind of reverse
+ * engineering: having to states of a point cloud and recover what happened our advantage here we
+ * know the identity of the vertex there are others methods giving other results.
+ *
+ * \param ob: Any object that can do soft-body e.g. mesh, lattice, curve.
+ * \param lloc: Output of the calculated location (or NULL).
+ * \param lrot: Output of the calculated rotation (or NULL).
+ * \param lscale: Output for the calculated scale (or NULL).
+ *
+ * For velocity & 2nd order stuff see: #vcloud_estimate_transform_v3.
+ */
+extern void SB_estimate_transform(Object *ob, float lloc[3], float lrot[3][3], float lscale[3][3]);
+
+#ifdef __cplusplus
+}
+#endif

source/blender/blenkernel/BKE_softbody.h

@@ -15,71 +15,54 @@ struct Object;
 struct Scene;
 struct SoftBody;
 
-typedef struct BodyPoint {
-  float origS[3], origE[3], origT[3], pos[3], vec[3], force[3];
-  float goal;
-  float prevpos[3], prevvec[3], prevdx[3], prevdv[3]; /* used for Heun integration */
-  float impdv[3], impdx[3];
-  int nofsprings;
-  int *springs;
-  float choke, choke2, frozen;
-  float colball;
-  short loc_flag; /* reserved by locale module specific states */
-  // char octantflag;
-  float mass;
-  float springweight;
+// typedef struct BodyPoint {
+//   float origS[3], origE[3], origT[3], pos[3], vec[3], force[3];
+//   float goal;
+//   float prevpos[3], prevvec[3], prevdx[3], prevdv[3]; /* used for Heun integration */
+//   float impdv[3], impdx[3];
+//   int nofsprings;
+//   int *springs;
+//   float choke, choke2, frozen;
+//   float colball;
+//   short loc_flag; /* reserved by locale module specific states */
+//   // char octantflag;
+//   float mass;
+//   float springweight;
+// } BodyPoint;
+
+typedef struct BodyPoint{
+  float x[3], v[3], a[3];
+  float x_prev[3], v_prev[3], a_prev[3];
+  float x_ini[3], v_ini[3], a_ini[3];
+
+  float mass_inv; // 1/mass
 } BodyPoint;
 
-/**
- * Allocates and initializes general main data.
- */
-extern struct SoftBody *sbNew(void);
+typedef struct BodyEdge{
+  int u, v;
+} BodyEdge;
 
-/**
- * Frees internal data and soft-body itself.
- */
-extern void sbFree(struct Object *ob);
+typedef struct BodyTet{
+  int verts[4];
+  float initial_volume;
+} BodyTet;
 
-/**
- * Frees simulation data to reset simulation.
- */
-extern void sbFreeSimulation(struct SoftBody *sb);
+struct SoftBody *init_softbody(void);
 
-/**
- * Do one simulation step, reading and writing vertex locs from given array.
- * */
-extern void sbObjectStep(struct Depsgraph *depsgraph,
-                         struct Scene *scene,
-                         struct Object *ob,
-                         float cfra,
-                         float (*vertexCos)[3],
-                         int numVerts);
+void free_softbody_intern(struct SoftBody *sb);
 
-/**
- * Makes totally fresh start situation, resets time.
- */
-extern void sbObjectToSoftbody(struct Object *ob);
+void sbFree(struct Object *ob);
 
-/**
- * Soft-body global visible functions.
- * Links the soft-body module to a 'test for Interrupt' function, pass NULL to clear the callback.
- */
-extern void sbSetInterruptCallBack(int (*f)(void));
+float get_tet_volume(BodyPoint *bpoint, BodyTet *curr_tet);
+
+void mesh_to_softbody(struct Object *ob, float (*vertexCos)[3], int numVerts);
+
+void sb_store_last_frame(struct Depsgraph *depsgraph, struct Object *object, float framenr);
+
+void softbody_to_object(struct Object *ob, float (*vertexCos)[3], int numVerts);
+
+void sbObjectStep(struct Depsgraph *depsgraph, struct Scene *scene, struct Object *ob, float cfra, float (*vertexCos)[3], int numVerts);
 
-/**
- * A precise position vector denoting the motion of the center of mass give a rotation/scale matrix
- * using averaging method, that's why estimate and not calculate see: this is kind of reverse
- * engineering: having to states of a point cloud and recover what happened our advantage here we
- * know the identity of the vertex there are others methods giving other results.
- *
- * \param ob: Any object that can do soft-body e.g. mesh, lattice, curve.
- * \param lloc: Output of the calculated location (or NULL).
- * \param lrot: Output of the calculated rotation (or NULL).
- * \param lscale: Output for the calculated scale (or NULL).
- *
- * For velocity & 2nd order stuff see: #vcloud_estimate_transform_v3.
- */
-extern void SB_estimate_transform(Object *ob, float lloc[3], float lrot[3][3], float lscale[3][3]);
 
 #ifdef __cplusplus
 }

source/blender/blenkernel/CMakeLists.txt

@@ -257,7 +257,7 @@ set(SRC
   intern/shader_fx.c
   intern/shrinkwrap.c
   intern/simulation.cc
-  intern/softbody.c
+  intern/softbody.cc
   intern/sound.c
   intern/speaker.c
   intern/spline_base.cc

source/blender/blenkernel/intern/object.cc

@@ -444,14 +444,14 @@ static void object_foreach_id(ID *id, LibraryForeachIDData *data)
         data, BKE_particlesystem_id_loop(psys, library_foreach_particlesystemsObjectLooper, data));
   }
 
-  if (object->soft) {
-    BKE_LIB_FOREACHID_PROCESS_IDSUPER(data, object->soft->collision_group, IDWALK_CB_NOP);
+  // if (object->soft) {
+  //   BKE_LIB_FOREACHID_PROCESS_IDSUPER(data, object->soft->collision_group, IDWALK_CB_NOP);
 
-    if (object->soft->effector_weights) {
-      BKE_LIB_FOREACHID_PROCESS_IDSUPER(
-          data, object->soft->effector_weights->group, IDWALK_CB_NOP);
-    }
-  }
+  //   if (object->soft->effector_weights) {
+  //     BKE_LIB_FOREACHID_PROCESS_IDSUPER(
+  //         data, object->soft->effector_weights->group, IDWALK_CB_NOP);
+  //   }
+  // }
 }
 
 static void object_foreach_path_pointcache(ListBase *ptcache_list,
@@ -562,12 +562,12 @@ static void object_blend_write(BlendWriter *writer, ID *id, const void *id_addre
   BLO_write_struct(writer, PartDeflect, ob->pd);
   if (ob->soft) {
     /* Set deprecated pointers to prevent crashes of older Blenders */
-    ob->soft->pointcache = ob->soft->shared->pointcache;
-    ob->soft->ptcaches = ob->soft->shared->ptcaches;
+    // ob->soft->pointcache = ob->soft->shared->pointcache;
+    // ob->soft->ptcaches = ob->soft->shared->ptcaches;
     BLO_write_struct(writer, SoftBody, ob->soft);
     BLO_write_struct(writer, SoftBody_Shared, ob->soft->shared);
     BKE_ptcache_blend_write(writer, &(ob->soft->shared->ptcaches));
-    BLO_write_struct(writer, EffectorWeights, ob->soft->effector_weights);
+    // BLO_write_struct(writer, EffectorWeights, ob->soft->effector_weights);
   }
 
   if (ob->rigidbody_object) {
@@ -715,21 +715,21 @@ static void object_blend_read_data(BlendDataReader *reader, ID *id)
     SoftBody *sb = ob->soft;
 
     sb->bpoint = nullptr; /* init pointers so it gets rebuilt nicely */
-    sb->bspring = nullptr;
-    sb->scratch = nullptr;
-    /* although not used anymore */
-    /* still have to be loaded to be compatible with old files */
-    BLO_read_pointer_array(reader, (void **)&sb->keys);
-    if (sb->keys) {
-      for (int a = 0; a < sb->totkey; a++) {
-        BLO_read_data_address(reader, &sb->keys[a]);
-      }
-    }
+    // sb->bspring = nullptr;
+    // sb->scratch = nullptr;
+    // /* although not used anymore */
+    // /* still have to be loaded to be compatible with old files */
+    // BLO_read_pointer_array(reader, (void **)&sb->keys);
+    // if (sb->keys) {
+    //   for (int a = 0; a < sb->totkey; a++) {
+    //     BLO_read_data_address(reader, &sb->keys[a]);
+    //   }
+    // }
 
-    BLO_read_data_address(reader, &sb->effector_weights);
-    if (!sb->effector_weights) {
-      sb->effector_weights = BKE_effector_add_weights(nullptr);
-    }
+    // BLO_read_data_address(reader, &sb->effector_weights);
+    // if (!sb->effector_weights) {
+    //   sb->effector_weights = BKE_effector_add_weights(nullptr);
+    // }
 
     BLO_read_data_address(reader, &sb->shared);
     if (sb->shared == nullptr) {
@@ -737,7 +737,7 @@ static void object_blend_read_data(BlendDataReader *reader, ID *id)
        * We should only do this when sb->shared == nullptr, because those pointers
        * are always set (for compatibility with older Blenders). We mustn't link
        * the same pointcache twice. */
-      BKE_ptcache_blend_read_data(reader, &sb->ptcaches, &sb->pointcache, false);
+      // BKE_ptcache_blend_read_data(reader, &sb->ptcaches, &sb->pointcache, false);
     }
     else {
       /* link caches */
@@ -989,11 +989,11 @@ static void object_blend_read_lib(BlendLibReader *reader, ID *id)
     BKE_particle_partdeflect_blend_read_lib(reader, &ob->id, ob->pd);
   }
 
-  if (ob->soft) {
-    BLO_read_id_address(reader, ob->id.lib, &ob->soft->collision_group);
+  // if (ob->soft) {
+  //   BLO_read_id_address(reader, ob->id.lib, &ob->soft->collision_group);
 
-    BLO_read_id_address(reader, ob->id.lib, &ob->soft->effector_weights->group);
-  }
+  //   BLO_read_id_address(reader, ob->id.lib, &ob->soft->effector_weights->group);
+  // }
 
   BKE_particle_system_blend_read_lib(reader, ob, &ob->id, &ob->particlesystem);
   BKE_modifier_blend_read_lib(reader, ob);
@@ -1108,13 +1108,13 @@ static void object_blend_read_expand(BlendExpander *expander, ID *id)
     BLO_expand(expander, ob->pd->f_source);
   }
 
-  if (ob->soft) {
-    BLO_expand(expander, ob->soft->collision_group);
+  // if (ob->soft) {
+  //   BLO_expand(expander, ob->soft->collision_group);
 
-    if (ob->soft->effector_weights) {
-      BLO_expand(expander, ob->soft->effector_weights->group);
-    }
-  }
+  //   if (ob->soft->effector_weights) {
+  //     BLO_expand(expander, ob->soft->effector_weights->group);
+  //   }
+  // }
 
   if (ob->rigidbody_constraint) {
     BLO_expand(expander, ob->rigidbody_constraint->ob1);
@@ -2332,12 +2332,13 @@ void BKE_object_copy_softbody(Object *ob_dst, const Object *ob_src, const int fl
   SoftBody *sbn = (SoftBody *)MEM_dupallocN(sb);
 
   if ((flag & LIB_ID_COPY_CACHES) == 0) {
-    sbn->totspring = sbn->totpoint = 0;
+    // sbn->totspring = sbn->totpoint = 0;
+    sbn->totpoint = 0;
     sbn->bpoint = nullptr;
-    sbn->bspring = nullptr;
+    // sbn->bspring = nullptr;
   }
   else {
-    sbn->totspring = sb->totspring;
+    // sbn->totspring = sb->totspring;
     sbn->totpoint = sb->totpoint;
 
     if (sbn->bpoint) {
@@ -2345,22 +2346,22 @@ void BKE_object_copy_softbody(Object *ob_dst, const Object *ob_src, const int fl
 
       sbn->bpoint = (BodyPoint *)MEM_dupallocN(sbn->bpoint);
 
-      for (i = 0; i < sbn->totpoint; i++) {
-        if (sbn->bpoint[i].springs) {
-          sbn->bpoint[i].springs = (int *)MEM_dupallocN(sbn->bpoint[i].springs);
-        }
-      }
+      // for (i = 0; i < sbn->totpoint; i++) {
+      //   if (sbn->bpoint[i].springs) {
+      //     sbn->bpoint[i].springs = (int *)MEM_dupallocN(sbn->bpoint[i].springs);
+      //   }
+      // }
     }
 
-    if (sb->bspring) {
-      sbn->bspring = (struct BodySpring *)MEM_dupallocN(sb->bspring);
-    }
+    // if (sb->bspring) {
+    //   sbn->bspring = (struct BodySpring *)MEM_dupallocN(sb->bspring);
+    // }
   }
 
-  sbn->keys = nullptr;
-  sbn->totkey = sbn->totpointkey = 0;
+  // sbn->keys = nullptr;
+  // sbn->totkey = sbn->totpointkey = 0;
 
-  sbn->scratch = nullptr;
+  // sbn->scratch = nullptr;
 
   if (is_orig) {
     sbn->shared = (SoftBody_Shared *)MEM_dupallocN(sb->shared);
@@ -2368,9 +2369,9 @@ void BKE_object_copy_softbody(Object *ob_dst, const Object *ob_src, const int fl
         &sbn->shared->ptcaches, &sb->shared->ptcaches, flag);
   }
 
-  if (sb->effector_weights) {
-    sbn->effector_weights = (EffectorWeights *)MEM_dupallocN(sb->effector_weights);
-  }
+  // if (sb->effector_weights) {
+  //   sbn->effector_weights = (EffectorWeights *)MEM_dupallocN(sb->effector_weights);
+  // }
 
   ob_dst->soft = sbn;
 }

source/blender/blenkernel/intern/pointcache.c

@@ -178,8 +178,11 @@ static int ptcache_softbody_write(int index, void *soft_v, void **data, int UNUS
   SoftBody *soft = soft_v;
   BodyPoint *bp = soft->bpoint + index;
 
-  PTCACHE_DATA_FROM(data, BPHYS_DATA_LOCATION, bp->pos);
-  PTCACHE_DATA_FROM(data, BPHYS_DATA_VELOCITY, bp->vec);
+  // PTCACHE_DATA_FROM(data, BPHYS_DATA_LOCATION, bp->pos);
+  // PTCACHE_DATA_FROM(data, BPHYS_DATA_VELOCITY, bp->vec);
+
+  PTCACHE_DATA_FROM(data, BPHYS_DATA_LOCATION, bp->x);
+  PTCACHE_DATA_FROM(data, BPHYS_DATA_VELOCITY, bp->v);
 
   return 1;
 }
@@ -189,13 +192,22 @@ static void ptcache_softbody_read(
   SoftBody *soft = soft_v;
   BodyPoint *bp = soft->bpoint + index;
 
+  // if (old_data) {
+  //   memcpy(bp->pos, data, sizeof(float[3]));
+  //   memcpy(bp->vec, data + 3, sizeof(float[3]));
+  // }
+  // else {
+  //   PTCACHE_DATA_TO(data, BPHYS_DATA_LOCATION, 0, bp->pos);
+  //   PTCACHE_DATA_TO(data, BPHYS_DATA_VELOCITY, 0, bp->vec);
+  // }
+
   if (old_data) {
-    memcpy(bp->pos, data, sizeof(float[3]));
-    memcpy(bp->vec, data + 3, sizeof(float[3]));
+    memcpy(bp->x, data, sizeof(float[3]));
+    memcpy(bp->v, data + 3, sizeof(float[3]));
   }
   else {
-    PTCACHE_DATA_TO(data, BPHYS_DATA_LOCATION, 0, bp->pos);
-    PTCACHE_DATA_TO(data, BPHYS_DATA_VELOCITY, 0, bp->vec);
+    PTCACHE_DATA_TO(data, BPHYS_DATA_LOCATION, 0, bp->x);
+    PTCACHE_DATA_TO(data, BPHYS_DATA_VELOCITY, 0, bp->v);
   }
 }
 static void ptcache_softbody_interpolate(int index,
@@ -204,7 +216,7 @@ static void ptcache_softbody_interpolate(int index,
                                          float cfra,
                                          float cfra1,
                                          float cfra2,
-                                         const float *old_data)
+                                         const float *old_data) 
 {
   SoftBody *soft = soft_v;
   BodyPoint *bp = soft->bpoint + index;
@@ -215,8 +227,10 @@ static void ptcache_softbody_interpolate(int index,
     return;
   }
 
-  copy_v3_v3(keys[1].co, bp->pos);
-  copy_v3_v3(keys[1].vel, bp->vec);
+  // copy_v3_v3(keys[1].co, bp->pos);
+  // copy_v3_v3(keys[1].vel, bp->vec);
+  copy_v3_v3(keys[1].co, bp->x);
+  copy_v3_v3(keys[1].vel, bp->v);
 
   if (old_data) {
     memcpy(keys[2].co, old_data, sizeof(float[3]));
@@ -235,8 +249,10 @@ static void ptcache_softbody_interpolate(int index,
 
   mul_v3_fl(keys->vel, 1.0f / dfra);
 
-  copy_v3_v3(bp->pos, keys->co);
-  copy_v3_v3(bp->vec, keys->vel);
+  // copy_v3_v3(bp->pos, keys->co);
+  // copy_v3_v3(bp->vec, keys->vel);
+  copy_v3_v3(bp->x, keys->co);
+  copy_v3_v3(bp->v, keys->vel);
 }
 static int ptcache_softbody_totpoint(void *soft_v, int UNUSED(cfra))
 {
@@ -2907,7 +2923,8 @@ int BKE_ptcache_id_reset(Scene *scene, PTCacheID *pid, int mode)
       cloth_free_modifier(pid->calldata);
     }
     else if (pid->type == PTCACHE_TYPE_SOFTBODY) {
-      sbFreeSimulation(pid->calldata);
+      // sbFreeSimulation(pid->calldata);
+      free_softbody_intern(pid->calldata);
     }
     else if (pid->type == PTCACHE_TYPE_PARTICLES) {
       psys_reset(pid->calldata, PSYS_RESET_DEPSGRAPH);

source/blender/blenkernel/intern/softbody copy.c

@@ -3199,7 +3199,7 @@ void sbSetInterruptCallBack(int (*f)(void))
   SB_localInterruptCallBack = f;
 }
 
-static void softbody_update_positions(Object *ob,
+static void   softbody_update_positions(Object *ob,
                                       SoftBody *sb,
                                       float (*vertexCos)[3],
                                       int numVerts)

source/blender/blenkernel/intern/softbody.cc

@@ -0,0 +1,420 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright Blender Foundation. All rights reserved. */
+
+/** \file
+ * \ingroup bke
+ */
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "CLG_log.h"
+
+#include "MEM_guardedalloc.h"
+
+/* types */
+#include "DNA_collection_types.h"
+#include "DNA_curve_types.h"
+#include "DNA_lattice_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_object_force_types.h"
+#include "DNA_object_types.h"
+#include "DNA_scene_types.h"
+
+#include "BLI_ghash.h"
+#include "BLI_listbase.h"
+#include "BLI_math.h"
+#include "BLI_threads.h"
+#include "BLI_utildefines.h"
+
+#include "BLI_vector.hh"
+#include "BLI_map.hh"
+#include "BLI_set.hh"
+#include <bits/stdc++.h>
+
+#include "BKE_collection.h"
+#include "BKE_collision.h"
+#include "BKE_curve.h"
+#include "BKE_deform.h"
+#include "BKE_effect.h"
+#include "BKE_global.h"
+#include "BKE_layer.h"
+#include "BKE_mesh.h"
+#include "BKE_modifier.h"
+#include "BKE_pointcache.h"
+#include "BKE_scene.h"
+#include "BKE_softbody.h"
+
+#include "DEG_depsgraph.h"
+#include "DEG_depsgraph_query.h"
+
+#include "PIL_time.h"
+
+#include "../../simulation/intern/xpbd.h"
+
+using blender::Set;
+// using blender::Map;
+using blender::Vector;
+using namespace std;
+
+struct VectorHasher {
+    int operator()(const Vector<int> &V) const {
+        int hash = V.size();
+        for(auto &i : V) {
+            hash ^= i + 0x9e3779b9 + (hash << 6) + (hash >> 2);
+        }
+        return hash;
+    }
+};
+
+static CLG_LogRef LOG = {"bke.softbody"};
+
+/* callbacks for errors and interrupts and some goo */
+static int (*SB_localInterruptCallBack)(void) = NULL;
+
+int tet_face[4][4] = {{2,1,0,3}, {0,1,3,2}, {1,2,3,0}, {2,0,3,1}}; 
+
+SoftBody *init_softbody()
+{
+  SoftBody *sb;
+
+  sb = (SoftBody *)MEM_callocN(sizeof(SoftBody), "softbody");
+
+  sb->totpoint = 0;
+  sb->tottet = 0;
+  sb->tot_surface_point = 0;
+  sb->tot_surface_tet = 0;
+  sb->bpoint = NULL;
+  sb->btet = NULL;
+  sb->surface_points = NULL;
+  sb->surface_tets = NULL;
+
+  sb->dt = 1.0/30.0;
+  sb->alpha_edge = 0.1;
+  sb->alpha_vol = 0;
+  sb->substep_count = 50;
+  sb->grav = 9.8f;
+
+  sb->shared = (SoftBody_Shared *)MEM_callocN(sizeof(*sb->shared), "SoftBody_Shared");
+  sb->shared->pointcache = BKE_ptcache_add(&sb->shared->ptcaches);
+
+  sb->last_frame = MINFRAME - 1;
+
+  return sb;
+}
+
+void free_softbody_intern(SoftBody *sb){
+  if(sb->bpoint){
+    MEM_freeN(sb->bpoint);
+  }
+  if(sb->bedge){
+    MEM_freeN(sb->bedge);
+  }
+  if(sb->btet){
+    MEM_freeN(sb->btet);
+  }
+
+  sb->bpoint = NULL;
+  sb->bedge = NULL;
+  sb->btet = NULL;
+
+  sb->totpoint = 0;
+  sb->tottet = 0;
+
+  if(sb->surface_points){
+    MEM_freeN(sb->surface_points);
+  }
+  if(sb->surface_tets){
+    MEM_freeN(sb->surface_tets);
+  }
+  sb->surface_points = NULL;
+  sb->surface_tets = NULL;
+  sb->tot_surface_point = 0;
+  sb->tot_surface_tet = 0;
+}
+
+void sbFree(Object *ob)
+{
+  SoftBody *sb = ob->soft;
+  if (sb == NULL) {
+    return;
+  }
+
+  const bool is_orig = (ob->id.tag & LIB_TAG_COPIED_ON_WRITE) == 0;
+
+  free_softbody_intern(sb);
+
+  if (is_orig) {
+    /* Only free shared data on non-CoW copies */
+    BKE_ptcache_free_list(&sb->shared->ptcaches);
+    sb->shared->pointcache = NULL;
+    MEM_freeN(sb->shared);
+  }
+  // if (sb->effector_weights) {
+  //   MEM_freeN(sb->effector_weights);
+  // }
+  MEM_freeN(sb);
+
+  ob->soft = NULL;
+}
+
+void mesh_to_softbody(Object *ob, float (*vertexCos)[3], int numVerts){
+  SoftBody *sb;
+  Mesh *me = (Mesh *)ob->data;
+  MLoop *mloop = me->mloop;
+  BodyPoint *bpoint;
+  BodyTet *btet;
+  int numTets = me->totpoly;
+
+  if(ob->soft == NULL){
+    ob->soft = init_softbody();
+  }
+  else{
+    free_softbody_intern(ob->soft);
+  }
+  sb = ob->soft;
+
+  // Set bodypoints
+  sb->bpoint = (BodyPoint *)MEM_mallocN(numVerts * sizeof(BodyPoint), "bodypoint");
+  for(int i = 0; i<numVerts; i++){
+    
+    float temp[3];
+    copy_v3_v3(temp, vertexCos[i]);
+
+    copy_v3_v3(sb->bpoint[i].x, vertexCos[i]);
+    mul_m4_v3(ob->obmat, sb->bpoint[i].x);
+    copy_v3_v3(sb->bpoint[i].x_ini, sb->bpoint[i].x);
+    sb->bpoint[i].mass_inv = 0.0f;
+
+    copy_v3_fl(sb->bpoint[i].v, 0.0);
+    copy_v3_fl(sb->bpoint[i].v_ini, 0.0);
+    copy_v3_fl(sb->bpoint[i].a, 0.0);
+    copy_v3_fl(sb->bpoint[i].a_ini, 0.0);
+  }
+  sb->totpoint = numVerts;
+  bpoint = sb->bpoint;
+
+  // Set tets
+  sb->btet = (BodyTet *)MEM_mallocN(numTets * sizeof(BodyTet), "bodytet");
+  btet = sb->btet;
+  for(int i = 0; i<numTets; i++){
+    MPoly *face = &me->mpoly[i];
+
+    for(int vert = 0; vert<4; vert++){
+      btet[i].verts[vert] = mloop[face->loopstart + vert].v;
+    }
+
+    BodyTet *curr_tet = &sb->btet[i];
+    curr_tet->initial_volume = get_tet_volume(bpoint, curr_tet);
+
+    for(int vert = 0; vert<4; vert++){
+      bpoint[btet[i].verts[vert]].mass_inv += 4.0/curr_tet->initial_volume; 
+    }
+  }
+  sb->tottet = numTets;
+
+  // Set edges
+  Set<pair<int,int>> edge_set;
+  for(int tetnr = 0; tetnr < sb->tottet; tetnr++){
+    // Brute force, iterating over all vertex pairs for a tet
+    for(int verti = 0; verti<4; verti++){
+        for(int vertj = verti+1; vertj<4; vertj++){
+            int vert1 = btet[tetnr].verts[verti];
+            int vert2 = btet[tetnr].verts[vertj];
+
+            int temp = min(vert1, vert2);
+            vert2 = max(vert1, vert2);
+            vert1 = temp;
+
+            if(edge_set.contains({vert1, vert2})){
+                continue;
+            }
+
+            edge_set.add({vert1, vert2});
+        }
+    }
+  }
+  sb->bedge = (BodyEdge *)MEM_mallocN(edge_set.size() * sizeof(BodyEdge), "bodyedge");
+  int i = 0; 
+  for(auto it : edge_set){
+    sb->bedge[i].u = it.first;
+    sb->bedge[i].v = it.second;
+    i++;
+  }
+  sb->totedge = edge_set.size();
+
+  // Set boundary points and tets
+  unordered_map<Vector<int>, int, VectorHasher> face_tet;
+  for(int tetnr = 0; tetnr < sb->tottet; tetnr++){
+    for(int facenr = 0; facenr<4; facenr++){
+      Vector<int> temp;
+      temp.append(btet[tetnr].verts[tet_face[facenr][0]]);
+      temp.append(btet[tetnr].verts[tet_face[facenr][1]]);
+      temp.append(btet[tetnr].verts[tet_face[facenr][2]]);
+
+      sort(temp.begin(), temp.end());
+
+      if(face_tet.find(temp) != face_tet.end()){
+        face_tet.erase(temp);
+        continue;
+      }
+
+      face_tet[temp] = tetnr;
+    }
+  }
+
+  Set<int> surface_point_set;
+  Set<int> surface_tet_set;
+  // sb->surface_tets = (int *)MEM_mallocN(face_tet.size() * sizeof(int), "surface_tets");
+  // sb->tot_surface_tet = face_tet.size();
+
+  int surface_tet_nr = 0;
+  for(auto it : face_tet){
+    for(int i = 0; i<3; i++){
+      surface_point_set.add(it.first[i]);
+    }
+
+    // sb->surface_tets[surface_tet_nr++] = it.second;
+    surface_tet_set.add(it.second);
+  }
+
+  sb->surface_tets = (int *)MEM_mallocN(face_tet.size() * sizeof(int), "surface_tets");
+  sb->tot_surface_tet = surface_tet_set.size();
+
+  for(auto it : surface_tet_set){
+    sb->surface_tets[surface_tet_nr++] = it;
+  }
+
+  sb->surface_points = (int *)MEM_mallocN(surface_point_set.size() * sizeof(int), "surface points");
+  sb->tot_surface_point = surface_point_set.size();
+  int surface_point_nr = 0;
+  for(auto it : surface_point_set){
+    sb->surface_points[surface_point_nr++] = it;
+  }
+}
+
+void softbody_to_object(Object *ob, float (*vertexCos)[3], int numVerts){
+  invert_m4_m4(ob->imat, ob->obmat);
+  for(int i = 0; i<numVerts; i++){
+    float temp[3];
+    copy_v3_v3(temp, ob->soft->bpoint[i].x);
+
+    copy_v3_v3(vertexCos[i], ob->soft->bpoint[i].x);
+    mul_m4_v3(ob->imat, vertexCos[i]);
+  }
+}
+
+void sb_store_last_frame(struct Depsgraph *depsgraph, Object *object, float framenr)
+{
+  if (!DEG_is_active(depsgraph)) {
+    return;
+  }
+  Object *object_orig = DEG_get_original_object(object);
+  object->soft->last_frame = framenr;
+  object_orig->soft->last_frame = framenr;
+}
+
+void sbObjectStep(struct Depsgraph *depsgraph,
+                  Scene *scene,
+                  Object *ob,
+                  float cfra,
+                  float (*vertexCos)[3],
+                  int numVerts){
+  SoftBody *sb = ob->soft;
+  PointCache *cache;
+  PTCacheID pid;
+  float dtime, timescale;
+  int framedelta, framenr, startframe, endframe;
+  int cache_result;
+  cache = sb->shared->pointcache;
+
+  framenr = (int)cfra;
+  framedelta = framenr - cache->simframe;
+
+  BKE_ptcache_id_from_softbody(&pid, ob, sb);
+  BKE_ptcache_id_time(&pid, scene, framenr, &startframe, &endframe, &timescale);
+
+  /* check for changes in mesh, should only happen in case the mesh
+   * structure changes during an animation */
+  if (sb->bpoint && numVerts != sb->totpoint) {
+    BKE_ptcache_invalidate(cache);
+    return;
+  }
+
+  /* clamp frame ranges */
+  if (framenr < startframe) {
+    BKE_ptcache_invalidate(cache);
+    return;
+  }
+  if (framenr > endframe) {
+    framenr = endframe;
+  }
+
+  /* verify if we need to create the softbody data */
+  if(sb == NULL || sb->bpoint == NULL){
+    mesh_to_softbody(ob, vertexCos, numVerts);
+  }
+
+  /* still no points? go away */
+  if (sb->totpoint == 0) {
+    return;
+  }
+  if(framenr == startframe){
+    BKE_ptcache_id_reset(scene, &pid, PTCACHE_RESET_OUTDATED);
+
+    // shift code that updates sb data and reinitialize it here
+
+    BKE_ptcache_validate(cache, framenr);
+    cache->flag &= ~PTCACHE_REDO_NEEDED;
+
+    sb_store_last_frame(depsgraph, ob, framenr);
+
+    return;
+  }
+
+  /* try to read from cache */
+  bool can_write_cache = DEG_is_active(depsgraph);
+  bool can_simulate = (framenr == sb->last_frame + 1) && !(cache->flag & PTCACHE_BAKED) &&
+                      can_write_cache;
+
+  cache_result = BKE_ptcache_read(&pid, (float)framenr + scene->r.subframe, can_simulate);
+
+  if (cache_result == PTCACHE_READ_EXACT || cache_result == PTCACHE_READ_INTERPOLATED ||
+      (!can_simulate && cache_result == PTCACHE_READ_OLD)) {
+    softbody_to_object(ob, vertexCos, numVerts);
+
+    BKE_ptcache_validate(cache, framenr);
+
+    if (cache_result == PTCACHE_READ_INTERPOLATED && cache->flag & PTCACHE_REDO_NEEDED &&
+        can_write_cache) {
+      BKE_ptcache_write(&pid, framenr);
+    }
+
+    sb_store_last_frame(depsgraph, ob, framenr);
+
+    return;
+  }
+
+  /* if on second frame, write cache for first frame */
+  if (cache->simframe == startframe &&
+      (cache->flag & PTCACHE_OUTDATED || cache->last_exact == 0)) {
+    BKE_ptcache_write(&pid, startframe);
+  }
+
+  for(int i = 0; i<sb->substep_count; i++){
+    xpbd_position_update(sb);
+
+    xpbd_enforce_constraints(sb);
+
+    xpbd_velocity_update(sb);
+  }
+
+  softbody_to_object(ob, vertexCos, numVerts);
+
+  BKE_ptcache_validate(cache, framenr);
+  BKE_ptcache_write(&pid, framenr);
+
+  sb_store_last_frame(depsgraph, ob, framenr);
+}

source/blender/blenloader/intern/versioning_250.c

@@ -934,9 +934,9 @@ void blo_do_versions_250(FileData *fd, Library *lib, Main *bmain)
           }
         }
 
-        if (ob->soft) {
-          ob->soft->effector_weights->global_gravity = ob->soft->grav / 9.81f;
-        }
+        // if (ob->soft) {
+        //   ob->soft->effector_weights->global_gravity = ob->soft->grav / 9.81f;
+        // }
       }
 
       /* Normal wind shape is plane */

source/blender/blenloader/intern/versioning_280.c

@@ -2675,11 +2675,11 @@ void blo_do_versions_280(FileData *fd, Library *UNUSED(lib), Main *bmain)
         }
 
         /* Move shared pointers from deprecated location to current location */
-        sb->shared->pointcache = sb->pointcache;
-        sb->shared->ptcaches = sb->ptcaches;
+        // sb->shared->pointcache = sb->pointcache;
+        // sb->shared->ptcaches = sb->ptcaches;
 
-        sb->pointcache = NULL;
-        BLI_listbase_clear(&sb->ptcaches);
+        // sb->pointcache = NULL;
+        // BLI_listbase_clear(&sb->ptcaches);
       }
     }
 

source/blender/blenloader/intern/versioning_legacy.c

@@ -1259,21 +1259,21 @@ void blo_do_versions_pre250(FileData *fd, Library *lib, Main *bmain)
 
     /* softbody init new vars */
     for (ob = bmain->objects.first; ob; ob = ob->id.next) {
-      if (ob->soft) {
-        if (ob->soft->defgoal == 0.0f) {
-          ob->soft->defgoal = 0.7f;
-        }
-        if (ob->soft->physics_speed == 0.0f) {
-          ob->soft->physics_speed = 1.0f;
-        }
-      }
-      if (ob->soft && ob->soft->vertgroup == 0) {
-        bDeformGroup *locGroup = BKE_object_defgroup_find_name(ob, "SOFTGOAL");
-        if (locGroup) {
-          /* retrieve index for that group */
-          ob->soft->vertgroup = 1 + BLI_findindex(&ob->defbase, locGroup);
-        }
-      }
+      // if (ob->soft) {
+      //   if (ob->soft->defgoal == 0.0f) {
+      //     ob->soft->defgoal = 0.7f;
+      //   }
+      //   if (ob->soft->physics_speed == 0.0f) {
+      //     ob->soft->physics_speed = 1.0f;
+      //   }
+      // }
+      // if (ob->soft && ob->soft->vertgroup == 0) {
+      //   bDeformGroup *locGroup = BKE_object_defgroup_find_name(ob, "SOFTGOAL");
+      //   if (locGroup) {
+      //     /* retrieve index for that group */
+      //     ob->soft->vertgroup = 1 + BLI_findindex(&ob->defbase, locGroup);
+      //   }
+      // }
     }
   }
 
@@ -2018,9 +2018,9 @@ void blo_do_versions_pre250(FileData *fd, Library *lib, Main *bmain)
 
     /* add point caches */
     for (ob = bmain->objects.first; ob; ob = ob->id.next) {
-      if (ob->soft && !ob->soft->pointcache) {
-        ob->soft->pointcache = BKE_ptcache_add(&ob->soft->ptcaches);
-      }
+      // if (ob->soft && !ob->soft->pointcache) {
+      //   ob->soft->pointcache = BKE_ptcache_add(&ob->soft->ptcaches);
+      // }
 
       for (psys = ob->particlesystem.first; psys; psys = psys->next) {
         if (psys->pointcache) {
@@ -2164,21 +2164,21 @@ void blo_do_versions_pre250(FileData *fd, Library *lib, Main *bmain)
       }
       do_version_constraints_245(&ob->constraints);
 
-      if (ob->soft && ob->soft->keys) {
-        SoftBody *sb = ob->soft;
-        int k;
+      // if (ob->soft && ob->soft->keys) {
+      //   SoftBody *sb = ob->soft;
+      //   int k;
 
-        for (k = 0; k < sb->totkey; k++) {
-          if (sb->keys[k]) {
-            MEM_freeN(sb->keys[k]);
-          }
-        }
+      //   for (k = 0; k < sb->totkey; k++) {
+      //     if (sb->keys[k]) {
+      //       MEM_freeN(sb->keys[k]);
+      //     }
+      //   }
 
-        MEM_freeN(sb->keys);
+      //   MEM_freeN(sb->keys);
 
-        sb->keys = NULL;
-        sb->totkey = 0;
-      }
+      //   sb->keys = NULL;
+      //   sb->totkey = 0;
+      // }
     }
   }
 
@@ -2188,21 +2188,21 @@ void blo_do_versions_pre250(FileData *fd, Library *lib, Main *bmain)
     PartEff *paf = NULL;
 
     for (ob = bmain->objects.first; ob; ob = ob->id.next) {
-      if (ob->soft && ob->soft->keys) {
-        SoftBody *sb = ob->soft;
-        int k;
+      // if (ob->soft && ob->soft->keys) {
+      //   SoftBody *sb = ob->soft;
+      //   int k;
 
-        for (k = 0; k < sb->totkey; k++) {
-          if (sb->keys[k]) {
-            MEM_freeN(sb->keys[k]);
-          }
-        }
+      //   for (k = 0; k < sb->totkey; k++) {
+      //     if (sb->keys[k]) {
+      //       MEM_freeN(sb->keys[k]);
+      //     }
+      //   }
 
-        MEM_freeN(sb->keys);
+      //   MEM_freeN(sb->keys);
 
-        sb->keys = NULL;
-        sb->totkey = 0;
-      }
+      //   sb->keys = NULL;
+      //   sb->totkey = 0;
+      // }
 
       /* convert old particles to new system */
       if ((paf = BKE_object_do_version_give_parteff_245(ob))) {
@@ -2385,10 +2385,10 @@ void blo_do_versions_pre250(FileData *fd, Library *lib, Main *bmain)
           strip->scale = 1.0f;
         }
       }
-      if (ob->soft) {
-        ob->soft->inpush = ob->soft->inspring;
-        ob->soft->shearstiff = 1.0f;
-      }
+      // if (ob->soft) {
+      //   ob->soft->inpush = ob->soft->inspring;
+      //   ob->soft->shearstiff = 1.0f;
+      // }
     }
   }
 

source/blender/editors/object/object_modifier.cc

@@ -195,7 +195,7 @@ ModifierData *ED_object_modifier_add(
     /* special cases */
     if (type == eModifierType_Softbody) {
       if (!ob->soft) {
-        ob->soft = sbNew();
+        ob->soft = init_softbody();
         ob->softflag |= OB_SB_GOAL | OB_SB_EDGES;
       }
     }

source/blender/makesdna/DNA_object_force_types copy.h

@@ -0,0 +1,395 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright 2004-2005 Blender Foundation. All rights reserved. */
+
+/** \file
+ * \ingroup DNA
+ */
+
+#pragma once
+
+#include "DNA_defs.h"
+#include "DNA_listBase.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct BodySpring;
+
+/** #PartDeflect.forcefield: Effector Fields types. */
+typedef enum ePFieldType {
+  /** (this is used for general effector weight). */
+  PFIELD_NULL = 0,
+  /** Force away/towards a point depending on force strength. */
+  PFIELD_FORCE = 1,
+  /** Force around the effector normal. */
+  PFIELD_VORTEX = 2,
+  /** Force from the cross product of effector normal and point velocity. */
+  PFIELD_MAGNET = 3,
+  /** Force away and towards a point depending which side of the effector normal the point is. */
+  PFIELD_WIND = 4,
+  /** Force along curve for dynamics, a shaping curve for hair paths. */
+  PFIELD_GUIDE = 5,
+  /** Force based on texture values calculated at point coordinates. */
+  PFIELD_TEXTURE = 6,
+  /** Force of a harmonic (damped) oscillator. */
+  PFIELD_HARMONIC = 7,
+  /** Force away/towards a point depending on point charge. */
+  PFIELD_CHARGE = 8,
+  /** Force due to a Lennard-Jones potential. */
+  PFIELD_LENNARDJ = 9,
+  /** Defines predator / goal for boids. */
+  PFIELD_BOID = 10,
+  /** Force defined by BLI_noise_generic_turbulence. */
+  PFIELD_TURBULENCE = 11,
+  /** Linear & quadratic drag. */
+  PFIELD_DRAG = 12,
+  /** Force based on fluid simulation velocities. */
+  PFIELD_FLUIDFLOW = 13,
+
+  /* Keep last. */
+  NUM_PFIELD_TYPES,
+} ePFieldType;
+
+typedef struct PartDeflect {
+  /** General settings flag. */
+  int flag;
+  /** Deflection flag - does mesh deflect particles. */
+  short deflect;
+  /** Force field type, do the vertices attract / repel particles? */
+  short forcefield;
+  /** Fall-off type. */
+  short falloff;
+  /** Point, plane or surface. */
+  short shape;
+  /** Texture effector. */
+  short tex_mode;
+  /** For curve guide. */
+  short kink, kink_axis;
+  short zdir;
+
+  /* Main effector values */
+  /** The strength of the force (+ or - ). */
+  float f_strength;
+  /** Damping ratio of the harmonic effector. */
+  float f_damp;
+  /**
+   * How much force is converted into "air flow", i.e.
+   * force used as the velocity of surrounding medium. */
+  float f_flow;
+  /** How much force is reduced when acting parallel to a surface, e.g. cloth. */
+  float f_wind_factor;
+
+  char _pad0[4];
+
+  /** Noise size for noise effector, restlength for harmonic effector. */
+  float f_size;
+
+  /* fall-off */
+  /** The power law - real gravitation is 2 (square). */
+  float f_power;
+  /** If indicated, use this maximum. */
+  float maxdist;
+  /** If indicated, use this minimum. */
+  float mindist;
+  /** Radial fall-off power. */
+  float f_power_r;
+  /** Radial versions of above. */
+  float maxrad;
+  float minrad;
+
+  /* particle collisions */
+  /** Damping factor for particle deflection. */
+  float pdef_damp;
+  /** Random element of damping for deflection. */
+  float pdef_rdamp;
+  /** Chance of particle passing through mesh. */
+  float pdef_perm;
+  /** Friction factor for particle deflection. */
+  float pdef_frict;
+  /** Random element of friction for deflection. */
+  float pdef_rfrict;
+  /** Surface particle stickiness. */
+  float pdef_stickness;
+
+  /** Used for forces. */
+  float absorption;
+
+  /* softbody collisions */
+  /** Damping factor for softbody deflection. */
+  float pdef_sbdamp;
+  /** Inner face thickness for softbody deflection. */
+  float pdef_sbift;
+  /** Outer face thickness for softbody deflection. */
+  float pdef_sboft;
+
+  /* guide curve, same as for particle child effects */
+  float clump_fac, clump_pow;
+  float kink_freq, kink_shape, kink_amp, free_end;
+
+  /* texture effector */
+  /** Used for calculating partial derivatives. */
+  float tex_nabla;
+  /** Texture of the texture effector. */
+  struct Tex *tex;
+
+  /* effector noise */
+  /** Random noise generator for e.g. wind. */
+  struct RNG *rng;
+  /** Noise of force. */
+  float f_noise;
+  /** Noise random seed. */
+  int seed;
+
+  /* Display Size */
+  /** Runtime only : start of the curve or draw scale. */
+  float drawvec1[4];
+  /** Runtime only : end of the curve. */
+  float drawvec2[4];
+  /** Runtime only. */
+  float drawvec_falloff_min[3];
+  char _pad1[4];
+  /** Runtime only. */
+  float drawvec_falloff_max[3];
+  char _pad2[4];
+
+  /** Force source object. */
+  struct Object *f_source;
+
+  /** Friction of cloth collisions. */
+  float pdef_cfrict;
+  char _pad[4];
+} PartDeflect;
+
+typedef struct EffectorWeights {
+  /** Only use effectors from this group of objects. */
+  struct Collection *group;
+
+  /** Effector type specific weights. */
+  float weight[14];
+  float global_gravity;
+  short flag;
+  char _pad[2];
+} EffectorWeights;
+
+/* EffectorWeights->flag */
+#define EFF_WEIGHT_DO_HAIR 1
+
+typedef struct SBVertex {
+  float vec[4];
+} SBVertex;
+
+/* Container for data that is shared among CoW copies.
+ *
+ * This is placed in a separate struct so that values can be changed
+ * without having to update all CoW copies. */
+typedef struct SoftBody_Shared {
+  struct PointCache *pointcache;
+  struct ListBase ptcaches;
+} SoftBody_Shared;
+
+typedef struct SoftBody {
+  /* dynamic data */
+  int totpoint, totspring;
+  /** Not saved in file. */
+  struct BodyPoint *bpoint;
+  /** Not saved in file. */
+  struct BodySpring *bspring;
+  char _pad;
+  char msg_lock;
+  short msg_value;
+
+  /* part of UI: */
+
+  /* general options */
+  /** Softbody mass of *vertex*. */
+  float nodemass;
+  /**
+   * Along with it introduce mass painting
+   * starting to fix old bug .. nastiness that VG are indexes
+   * rather find them by name tag to find it -> jow20090613.
+   * MAX_VGROUP_NAME */
+  char namedVG_Mass[64];
+  /** Softbody amount of gravitation to apply. */
+  float grav;
+  /** Friction to env. */
+  float mediafrict;
+  /** Error limit for ODE solver. */
+  float rklimit;
+  /** User control over simulation speed. */
+  float physics_speed;
+
+  /* goal */
+  /** Softbody goal springs. */
+  float goalspring;
+  /** Softbody goal springs friction. */
+  float goalfrict;
+  /** Quick limits for goal. */
+  float mingoal;
+  float maxgoal;
+  /** Default goal for vertices without vgroup. */
+  float defgoal;
+  /** Index starting at 1. */
+  short vertgroup;
+  /**
+   * Starting to fix old bug .. nastiness that VG are indexes
+   * rather find them by name tag to find it -> jow20090613.
+   * MAX_VGROUP_NAME */
+  char namedVG_Softgoal[64];
+
+  short fuzzyness;
+
+  /* springs */
+  /** Softbody inner springs. */
+  float inspring;
+  /** Softbody inner springs friction. */
+  float infrict;
+  /**
+   * Along with it introduce Spring_K painting
+   * starting to fix old bug .. nastiness that VG are indexes
+   * rather find them by name tag to find it -> jow20090613.
+   * MAX_VGROUP_NAME
+   */
+  char namedVG_Spring_K[64];
+
+  /* baking */
+  char _pad1[6];
+  /** Local==1: use local coords for baking. */
+  char local, solverflags;
+
+  /* -- these must be kept for backwards compatibility -- */
+  /** Array of size totpointkey. */
+  SBVertex **keys;
+  /** If totpointkey != totpoint or totkey!- (efra-sfra)/interval -> free keys. */
+  int totpointkey, totkey;
+  /* ---------------------------------------------------- */
+  float secondspring;
+
+  /* Self collision. */
+  /** Fixed collision ball size if > 0. */
+  float colball;
+  /** Cooling down collision response. */
+  float balldamp;
+  /** Pressure the ball is loaded with. */
+  float ballstiff;
+  short sbc_mode;
+  short aeroedge;
+  short minloops;
+  short maxloops;
+  short choke;
+  short solver_ID;
+  short plastic;
+  short springpreload;
+
+  /** Scratchpad/cache on live time not saved in file. */
+  struct SBScratch *scratch;
+  float shearstiff;
+  float inpush;
+
+  struct SoftBody_Shared *shared;
+  /** Moved to SoftBody_Shared. */
+  struct PointCache *pointcache DNA_DEPRECATED;
+  /** Moved to SoftBody_Shared. */
+  struct ListBase ptcaches DNA_DEPRECATED;
+
+  struct Collection *collision_group;
+
+  struct EffectorWeights *effector_weights;
+  /* Reverse estimated object-matrix (run-time data, no need to store in the file). */
+  float lcom[3];
+  float lrot[3][3];
+  float lscale[3][3];
+
+  int last_frame;
+} SoftBody;
+
+/* pd->flag: various settings */
+#define PFIELD_USEMAX (1 << 0)
+// #define PDEFLE_DEFORM         (1 << 1) /* UNUSED */
+/** TODO: do_versions for below */
+#define PFIELD_GUIDE_PATH_ADD (1 << 2)
+/** used for do_versions */
+#define PFIELD_PLANAR (1 << 3)
+#define PDEFLE_KILL_PART (1 << 4)
+/** used for do_versions */
+#define PFIELD_POSZ (1 << 5)
+#define PFIELD_TEX_OBJECT (1 << 6)
+/** used for turbulence */
+#define PFIELD_GLOBAL_CO (1 << 6)
+#define PFIELD_TEX_2D (1 << 7)
+/** used for harmonic force */
+#define PFIELD_MULTIPLE_SPRINGS (1 << 7)
+#define PFIELD_USEMIN (1 << 8)
+#define PFIELD_USEMAXR (1 << 9)
+#define PFIELD_USEMINR (1 << 10)
+#define PFIELD_TEX_ROOTCO (1 << 11)
+/** used for do_versions */
+#define PFIELD_SURFACE (1 << 12)
+#define PFIELD_VISIBILITY (1 << 13)
+#define PFIELD_DO_LOCATION (1 << 14)
+#define PFIELD_DO_ROTATION (1 << 15)
+/** apply curve weights */
+#define PFIELD_GUIDE_PATH_WEIGHT (1 << 16)
+/** multiply smoke force by density */
+#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)
+/** Replace collision direction with collider normal. */
+#define PFIELD_CLOTH_USE_NORMAL (1 << 20)
+
+/* pd->falloff */
+#define PFIELD_FALL_SPHERE 0
+#define PFIELD_FALL_TUBE 1
+#define PFIELD_FALL_CONE 2
+
+/* pd->shape */
+#define PFIELD_SHAPE_POINT 0
+#define PFIELD_SHAPE_PLANE 1
+#define PFIELD_SHAPE_SURFACE 2
+#define PFIELD_SHAPE_POINTS 3
+#define PFIELD_SHAPE_LINE 4
+
+/* pd->tex_mode */
+#define PFIELD_TEX_RGB 0
+#define PFIELD_TEX_GRAD 1
+#define PFIELD_TEX_CURL 2
+
+/* pd->zdir */
+#define PFIELD_Z_BOTH 0
+#define PFIELD_Z_POS 1
+#define PFIELD_Z_NEG 2
+
+/* ob->softflag */
+#define OB_SB_ENABLE 1 /* deprecated, use modifier */
+#define OB_SB_GOAL 2
+#define OB_SB_EDGES 4
+#define OB_SB_QUADS 8
+#define OB_SB_POSTDEF 16
+// #define OB_SB_REDO       32
+// #define OB_SB_BAKESET    64
+// #define OB_SB_BAKEDO 128
+// #define OB_SB_RESET      256
+#define OB_SB_SELF 512
+#define OB_SB_FACECOLL 1024
+#define OB_SB_EDGECOLL 2048
+/* #define OB_SB_COLLFINAL 4096 */ /* deprecated */
+/* #define OB_SB_BIG_UI 8192 */    /* deprecated */
+#define OB_SB_AERO_ANGLE 16384
+
+/* sb->solverflags */
+#define SBSO_MONITOR 1
+#define SBSO_OLDERR 2
+#define SBSO_ESTIMATEIPO 4
+
+/* sb->sbc_mode */
+#define SBC_MODE_MANUAL 0
+#define SBC_MODE_AVG 1
+#define SBC_MODE_MIN 2
+#define SBC_MODE_MAX 3
+#define SBC_MODE_AVGMINMAX 4
+
+#ifdef __cplusplus
+}
+#endif

source/blender/makesdna/DNA_object_force_types.h

@@ -188,121 +188,149 @@ typedef struct SoftBody_Shared {
   struct ListBase ptcaches;
 } SoftBody_Shared;
 
-typedef struct SoftBody {
-  /* dynamic data */
-  int totpoint, totspring;
+typedef struct SoftBody{
+  int totpoint, totedge, tottet, tot_surface_point, tot_surface_tet;
+  char _pad1[4];
   /** Not saved in file. */
   struct BodyPoint *bpoint;
-  /** Not saved in file. */
-  struct BodySpring *bspring;
-  char _pad;
-  char msg_lock;
-  short msg_value;
+  struct BodyEdge *bedge;
+  struct BodyTet *btet;
 
-  /* part of UI: */
+  int *surface_points;
+  int *surface_tets;
 
-  /* general options */
-  /** Softbody mass of *vertex*. */
-  float nodemass;
-  /**
-   * Along with it introduce mass painting
-   * starting to fix old bug .. nastiness that VG are indexes
-   * rather find them by name tag to find it -> jow20090613.
-   * MAX_VGROUP_NAME */
-  char namedVG_Mass[64];
-  /** Softbody amount of gravitation to apply. */
   float grav;
-  /** Friction to env. */
-  float mediafrict;
-  /** Error limit for ODE solver. */
-  float rklimit;
-  /** User control over simulation speed. */
-  float physics_speed;
+  float substep_count;
+  float alpha_vol; //stiffness coefficient
+  float alpha_edge;
+  float dt;
+  
+  int last_frame;
 
-  /* goal */
-  /** Softbody goal springs. */
-  float goalspring;
-  /** Softbody goal springs friction. */
-  float goalfrict;
-  /** Quick limits for goal. */
-  float mingoal;
-  float maxgoal;
-  /** Default goal for vertices without vgroup. */
-  float defgoal;
-  /** Index starting at 1. */
   short vertgroup;
-  /**
-   * Starting to fix old bug .. nastiness that VG are indexes
-   * rather find them by name tag to find it -> jow20090613.
-   * MAX_VGROUP_NAME */
-  char namedVG_Softgoal[64];
-
-  short fuzzyness;
-
-  /* springs */
-  /** Softbody inner springs. */
-  float inspring;
-  /** Softbody inner springs friction. */
-  float infrict;
-  /**
-   * Along with it introduce Spring_K painting
-   * starting to fix old bug .. nastiness that VG are indexes
-   * rather find them by name tag to find it -> jow20090613.
-   * MAX_VGROUP_NAME
-   */
-  char namedVG_Spring_K[64];
-
-  /* baking */
-  char _pad1[6];
-  /** Local==1: use local coords for baking. */
-  char local, solverflags;
-
-  /* -- these must be kept for backwards compatibility -- */
-  /** Array of size totpointkey. */
-  SBVertex **keys;
-  /** If totpointkey != totpoint or totkey!- (efra-sfra)/interval -> free keys. */
-  int totpointkey, totkey;
-  /* ---------------------------------------------------- */
-  float secondspring;
-
-  /* Self collision. */
-  /** Fixed collision ball size if > 0. */
-  float colball;
-  /** Cooling down collision response. */
-  float balldamp;
-  /** Pressure the ball is loaded with. */
-  float ballstiff;
-  short sbc_mode;
-  short aeroedge;
-  short minloops;
-  short maxloops;
-  short choke;
-  short solver_ID;
-  short plastic;
-  short springpreload;
-
-  /** Scratchpad/cache on live time not saved in file. */
-  struct SBScratch *scratch;
-  float shearstiff;
-  float inpush;
+  char _pad[6];
 
   struct SoftBody_Shared *shared;
-  /** Moved to SoftBody_Shared. */
-  struct PointCache *pointcache DNA_DEPRECATED;
-  /** Moved to SoftBody_Shared. */
-  struct ListBase ptcaches DNA_DEPRECATED;
 
-  struct Collection *collision_group;
-
-  struct EffectorWeights *effector_weights;
-  /* Reverse estimated object-matrix (run-time data, no need to store in the file). */
-  float lcom[3];
-  float lrot[3][3];
-  float lscale[3][3];
-
-  int last_frame;
+  // char _pad[4];
 } SoftBody;
 
+
+// typedef struct SoftBody {
+//   /* dynamic data */
+//   int totpoint, tottet;
+//   /** Not saved in file. */
+//   struct BodyPoint *bpoint;
+//   /** Not saved in file. */
+//   struct BodyTet *btet;
+//   char _pad;
+//   char msg_lock;
+//   short msg_value;
+
+//   /* part of UI: */
+
+//   /* general options */
+//   /** Softbody mass of *vertex*. */
+//   float nodemass;
+//   /**
+//    * Along with it introduce mass painting
+//    * starting to fix old bug .. nastiness that VG are indexes
+//    * rather find them by name tag to find it -> jow20090613.
+//    * MAX_VGROUP_NAME */
+//   char namedVG_Mass[64];
+//   /** Softbody amount of gravitation to apply. */
+//   float grav;
+//   /** Friction to env. */
+//   float mediafrict;
+//   /** Error limit for ODE solver. */
+//   float rklimit;
+//   /** User control over simulation speed. */
+//   float physics_speed;
+
+//   /* goal */
+//   /** Softbody goal springs. */
+//   float goalspring;
+//   /** Softbody goal springs friction. */
+//   float goalfrict;
+//   /** Quick limits for goal. */
+//   float mingoal;
+//   float maxgoal;
+//   /** Default goal for vertices without vgroup. */
+//   float defgoal;
+//   /** Index starting at 1. */
+//   short vertgroup;
+//   /**
+//    * Starting to fix old bug .. nastiness that VG are indexes
+//    * rather find them by name tag to find it -> jow20090613.
+//    * MAX_VGROUP_NAME */
+//   char namedVG_Softgoal[64];
+
+//   short fuzzyness;
+
+//   /* springs */
+//   /** Softbody inner springs. */
+//   float inspring;
+//   /** Softbody inner springs friction. */
+//   float infrict;
+//   /**
+//    * Along with it introduce Spring_K painting
+//    * starting to fix old bug .. nastiness that VG are indexes
+//    * rather find them by name tag to find it -> jow20090613.
+//    * MAX_VGROUP_NAME
+//    */
+//   char namedVG_Spring_K[64];
+
+//   /* baking */
+//   char _pad1[6];
+//   /** Local==1: use local coords for baking. */
+//   char local, solverflags;
+
+//   /* -- these must be kept for backwards compatibility -- */
+//   /** Array of size totpointkey. */
+//   SBVertex **keys;
+//   /** If totpointkey != totpoint or totkey!- (efra-sfra)/interval -> free keys. */
+//   int totpointkey, totkey;
+//   /* ---------------------------------------------------- */
+//   float secondspring;
+
+//   /* Self collision. */
+//   /** Fixed collision ball size if > 0. */
+//   float colball;
+//   /** Cooling down collision response. */
+//   float balldamp;
+//   /** Pressure the ball is loaded with. */
+//   float ballstiff;
+//   short sbc_mode;
+//   short aeroedge;
+//   short minloops;
+//   short maxloops;
+//   short choke;
+//   short solver_ID;
+//   short plastic;
+//   short springpreload;
+
+//   /** Scratchpad/cache on live time not saved in file. */
+//   struct SBScratch *scratch;
+//   float shearstiff;
+//   float inpush;
+
+//   struct SoftBody_Shared *shared;
+//   /** Moved to SoftBody_Shared. */
+//   struct PointCache *pointcache DNA_DEPRECATED;
+//   /** Moved to SoftBody_Shared. */
+//   struct ListBase ptcaches DNA_DEPRECATED;
+
+//   struct Collection *collision_group;
+
+//   struct EffectorWeights *effector_weights;
+//   /* Reverse estimated object-matrix (run-time data, no need to store in the file). */
+//   float lcom[3];
+//   float lrot[3][3];
+//   float lscale[3][3];
+
+//   int last_frame;
+// } SoftBody;
+
 /* pd->flag: various settings */
 #define PFIELD_USEMAX (1 << 0)
 // #define PDEFLE_DEFORM         (1 << 1) /* UNUSED */

source/blender/makesrna/intern/rna_object_force.c

@@ -465,159 +465,163 @@ static char *rna_CollisionSettings_path(const PointerRNA *UNUSED(ptr))
 #  endif
 }
 
-static bool rna_SoftBodySettings_use_edges_get(PointerRNA *ptr)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  return (((data->softflag) & OB_SB_EDGES) != 0);
-}
+//----------------------------
+// Soft Body Runtime Functions
+//----------------------------
 
-static void rna_SoftBodySettings_use_edges_set(PointerRNA *ptr, bool value)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  if (value) {
-    data->softflag |= OB_SB_EDGES;
-  }
-  else {
-    data->softflag &= ~OB_SB_EDGES;
-  }
-}
+// static bool rna_SoftBodySettings_use_edges_get(PointerRNA *ptr)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   return (((data->softflag) & OB_SB_EDGES) != 0);
+// }
 
-static bool rna_SoftBodySettings_use_goal_get(PointerRNA *ptr)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  return (((data->softflag) & OB_SB_GOAL) != 0);
-}
+// static void rna_SoftBodySettings_use_edges_set(PointerRNA *ptr, bool value)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   if (value) {
+//     data->softflag |= OB_SB_EDGES;
+//   }
+//   else {
+//     data->softflag &= ~OB_SB_EDGES;
+//   }
+// }
 
-static void rna_SoftBodySettings_use_goal_set(PointerRNA *ptr, bool value)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  if (value) {
-    data->softflag |= OB_SB_GOAL;
-  }
-  else {
-    data->softflag &= ~OB_SB_GOAL;
-  }
-}
+// static bool rna_SoftBodySettings_use_goal_get(PointerRNA *ptr)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   return (((data->softflag) & OB_SB_GOAL) != 0);
+// }
 
-static bool rna_SoftBodySettings_stiff_quads_get(PointerRNA *ptr)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  return (((data->softflag) & OB_SB_QUADS) != 0);
-}
+// static void rna_SoftBodySettings_use_goal_set(PointerRNA *ptr, bool value)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   if (value) {
+//     data->softflag |= OB_SB_GOAL;
+//   }
+//   else {
+//     data->softflag &= ~OB_SB_GOAL;
+//   }
+// }
 
-static void rna_SoftBodySettings_stiff_quads_set(PointerRNA *ptr, bool value)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  if (value) {
-    data->softflag |= OB_SB_QUADS;
-  }
-  else {
-    data->softflag &= ~OB_SB_QUADS;
-  }
-}
+// static bool rna_SoftBodySettings_stiff_quads_get(PointerRNA *ptr)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   return (((data->softflag) & OB_SB_QUADS) != 0);
+// }
 
-static bool rna_SoftBodySettings_self_collision_get(PointerRNA *ptr)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  return (((data->softflag) & OB_SB_SELF) != 0);
-}
+// static void rna_SoftBodySettings_stiff_quads_set(PointerRNA *ptr, bool value)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   if (value) {
+//     data->softflag |= OB_SB_QUADS;
+//   }
+//   else {
+//     data->softflag &= ~OB_SB_QUADS;
+//   }
+// }
 
-static void rna_SoftBodySettings_self_collision_set(PointerRNA *ptr, bool value)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  if (value) {
-    data->softflag |= OB_SB_SELF;
-  }
-  else {
-    data->softflag &= ~OB_SB_SELF;
-  }
-}
+// static bool rna_SoftBodySettings_self_collision_get(PointerRNA *ptr)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   return (((data->softflag) & OB_SB_SELF) != 0);
+// }
 
-static int rna_SoftBodySettings_new_aero_get(PointerRNA *ptr)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  if (data->softflag & OB_SB_AERO_ANGLE) {
-    return 1;
-  }
-  else {
-    return 0;
-  }
-}
+// static void rna_SoftBodySettings_self_collision_set(PointerRNA *ptr, bool value)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   if (value) {
+//     data->softflag |= OB_SB_SELF;
+//   }
+//   else {
+//     data->softflag &= ~OB_SB_SELF;
+//   }
+// }
 
-static void rna_SoftBodySettings_new_aero_set(PointerRNA *ptr, int value)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  if (value == 1) {
-    data->softflag |= OB_SB_AERO_ANGLE;
-  }
-  else { /* value == 0 */
-    data->softflag &= ~OB_SB_AERO_ANGLE;
-  }
-}
+// static int rna_SoftBodySettings_new_aero_get(PointerRNA *ptr)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   if (data->softflag & OB_SB_AERO_ANGLE) {
+//     return 1;
+//   }
+//   else {
+//     return 0;
+//   }
+// }
 
-static bool rna_SoftBodySettings_face_collision_get(PointerRNA *ptr)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  return (((data->softflag) & OB_SB_FACECOLL) != 0);
-}
+// static void rna_SoftBodySettings_new_aero_set(PointerRNA *ptr, int value)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   if (value == 1) {
+//     data->softflag |= OB_SB_AERO_ANGLE;
+//   }
+//   else { /* value == 0 */
+//     data->softflag &= ~OB_SB_AERO_ANGLE;
+//   }
+// }
 
-static void rna_SoftBodySettings_face_collision_set(PointerRNA *ptr, bool value)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  if (value) {
-    data->softflag |= OB_SB_FACECOLL;
-  }
-  else {
-    data->softflag &= ~OB_SB_FACECOLL;
-  }
-}
+// static bool rna_SoftBodySettings_face_collision_get(PointerRNA *ptr)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   return (((data->softflag) & OB_SB_FACECOLL) != 0);
+// }
 
-static bool rna_SoftBodySettings_edge_collision_get(PointerRNA *ptr)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  return (((data->softflag) & OB_SB_EDGECOLL) != 0);
-}
+// static void rna_SoftBodySettings_face_collision_set(PointerRNA *ptr, bool value)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   if (value) {
+//     data->softflag |= OB_SB_FACECOLL;
+//   }
+//   else {
+//     data->softflag &= ~OB_SB_FACECOLL;
+//   }
+// }
 
-static void rna_SoftBodySettings_edge_collision_set(PointerRNA *ptr, bool value)
-{
-  Object *data = (Object *)(ptr->owner_id);
-  if (value) {
-    data->softflag |= OB_SB_EDGECOLL;
-  }
-  else {
-    data->softflag &= ~OB_SB_EDGECOLL;
-  }
-}
+// static bool rna_SoftBodySettings_edge_collision_get(PointerRNA *ptr)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   return (((data->softflag) & OB_SB_EDGECOLL) != 0);
+// }
 
-static void rna_SoftBodySettings_goal_vgroup_get(PointerRNA *ptr, char *value)
-{
-  SoftBody *sb = (SoftBody *)ptr->data;
-  rna_object_vgroup_name_index_get(ptr, value, sb->vertgroup);
-}
+// static void rna_SoftBodySettings_edge_collision_set(PointerRNA *ptr, bool value)
+// {
+//   Object *data = (Object *)(ptr->owner_id);
+//   if (value) {
+//     data->softflag |= OB_SB_EDGECOLL;
+//   }
+//   else {
+//     data->softflag &= ~OB_SB_EDGECOLL;
+//   }
+// }
 
-static int rna_SoftBodySettings_goal_vgroup_length(PointerRNA *ptr)
-{
-  SoftBody *sb = (SoftBody *)ptr->data;
-  return rna_object_vgroup_name_index_length(ptr, sb->vertgroup);
-}
+// static void rna_SoftBodySettings_goal_vgroup_get(PointerRNA *ptr, char *value)
+// {
+//   SoftBody *sb = (SoftBody *)ptr->data;
+//   rna_object_vgroup_name_index_get(ptr, value, sb->vertgroup);
+// }
 
-static void rna_SoftBodySettings_goal_vgroup_set(PointerRNA *ptr, const char *value)
-{
-  SoftBody *sb = (SoftBody *)ptr->data;
-  rna_object_vgroup_name_index_set(ptr, value, &sb->vertgroup);
-}
+// static int rna_SoftBodySettings_goal_vgroup_length(PointerRNA *ptr)
+// {
+//   SoftBody *sb = (SoftBody *)ptr->data;
+//   return rna_object_vgroup_name_index_length(ptr, sb->vertgroup);
+// }
 
-static void rna_SoftBodySettings_mass_vgroup_set(PointerRNA *ptr, const char *value)
-{
-  SoftBody *sb = (SoftBody *)ptr->data;
-  rna_object_vgroup_name_set(ptr, value, sb->namedVG_Mass, sizeof(sb->namedVG_Mass));
-}
+// static void rna_SoftBodySettings_goal_vgroup_set(PointerRNA *ptr, const char *value)
+// {
+//   SoftBody *sb = (SoftBody *)ptr->data;
+//   rna_object_vgroup_name_index_set(ptr, value, &sb->vertgroup);
+// }
 
-static void rna_SoftBodySettings_spring_vgroup_set(PointerRNA *ptr, const char *value)
-{
-  SoftBody *sb = (SoftBody *)ptr->data;
-  rna_object_vgroup_name_set(ptr, value, sb->namedVG_Spring_K, sizeof(sb->namedVG_Spring_K));
-}
+// static void rna_SoftBodySettings_mass_vgroup_set(PointerRNA *ptr, const char *value)
+// {
+//   SoftBody *sb = (SoftBody *)ptr->data;
+//   rna_object_vgroup_name_set(ptr, value, sb->namedVG_Mass, sizeof(sb->namedVG_Mass));
+// }
+
+// static void rna_SoftBodySettings_spring_vgroup_set(PointerRNA *ptr, const char *value)
+// {
+//   SoftBody *sb = (SoftBody *)ptr->data;
+//   rna_object_vgroup_name_set(ptr, value, sb->namedVG_Spring_K, sizeof(sb->namedVG_Spring_K));
+// }
 
 static char *rna_SoftBodySettings_path(const PointerRNA *ptr)
 {
@@ -636,6 +640,10 @@ static int particle_id_check(const PointerRNA *ptr)
   return (GS(id->name) == ID_PA);
 }
 
+//----------------------
+//Field Settings Runtime
+//----------------------
+
 static void rna_FieldSettings_update(Main *UNUSED(bmain), Scene *UNUSED(scene), PointerRNA *ptr)
 {
   if (particle_id_check(ptr)) {
@@ -816,15 +824,15 @@ static char *rna_EffectorWeight_path(const PointerRNA *ptr)
     ModifierData *md;
 
     /* check softbody modifier */
-    md = (ModifierData *)BKE_modifiers_findby_type(ob, eModifierType_Softbody);
-    if (md) {
-      /* no pointer from modifier data to actual softbody storage, would be good to add */
-      if (ob->soft->effector_weights == ew) {
-        char name_esc[sizeof(md->name) * 2];
-        BLI_str_escape(name_esc, md->name, sizeof(name_esc));
-        return BLI_sprintfN("modifiers[\"%s\"].settings.effector_weights", name_esc);
-      }
-    }
+    // md = (ModifierData *)BKE_modifiers_findby_type(ob, eModifierType_Softbody);
+    // if (md) {
+    //   /* no pointer from modifier data to actual softbody storage, would be good to add */
+    //   if (ob->soft->effector_weights == ew) {
+    //     char name_esc[sizeof(md->name) * 2];
+    //     BLI_str_escape(name_esc, md->name, sizeof(name_esc));
+    //     return BLI_sprintfN("modifiers[\"%s\"].settings.effector_weights", name_esc);
+    //   }
+    // }
 
     /* check cloth modifier */
     md = (ModifierData *)BKE_modifiers_findby_type(ob, eModifierType_Cloth);
@@ -1861,318 +1869,344 @@ static void rna_def_softbody(BlenderRNA *brna)
   RNA_def_struct_ui_text(
       srna, "Soft Body Settings", "Soft body simulation settings for an object");
 
-  /* General Settings */
+  /* New Settings */
 
-  prop = RNA_def_property(srna, "friction", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "mediafrict");
-  RNA_def_property_range(prop, 0.0f, 50.0f);
-  RNA_def_property_ui_text(prop, "Friction", "General media friction for point movements");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "mass", PROP_FLOAT, PROP_UNIT_MASS);
-  RNA_def_property_float_sdna(prop, NULL, "nodemass");
-  RNA_def_property_range(prop, 0.0f, 50000.0f);
-  RNA_def_property_ui_text(prop, "Mass", "General Mass value");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "vertex_group_mass", PROP_STRING, PROP_NONE);
-  RNA_def_property_string_sdna(prop, NULL, "namedVG_Mass");
-  RNA_def_property_ui_text(prop, "Mass Vertex Group", "Control point mass values");
-  RNA_def_property_string_funcs(prop, NULL, NULL, "rna_SoftBodySettings_mass_vgroup_set");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  /* no longer used */
-  prop = RNA_def_property(srna, "gravity", PROP_FLOAT, PROP_ACCELERATION);
-  RNA_def_property_float_sdna(prop, NULL, "grav");
-  RNA_def_property_range(prop, -10.0f, 10.0f);
-  RNA_def_property_ui_text(prop, "Gravitation", "Apply gravitation to point movement");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "speed", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "physics_speed");
-  RNA_def_property_range(prop, 0.01f, 100.0f);
-  RNA_def_property_ui_text(
-      prop, "Speed", "Tweak timing for physics to control frequency and speed");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  /* Goal */
-
-  prop = RNA_def_property(srna, "vertex_group_goal", PROP_STRING, PROP_NONE);
-  RNA_def_property_string_sdna(prop, NULL, "vertgroup");
-  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE); /* not impossible .. but not supported yet */
-  RNA_def_property_string_funcs(prop,
-                                "rna_SoftBodySettings_goal_vgroup_get",
-                                "rna_SoftBodySettings_goal_vgroup_length",
-                                "rna_SoftBodySettings_goal_vgroup_set");
-  RNA_def_property_ui_text(prop, "Goal Vertex Group", "Control point weight values");
-
-  prop = RNA_def_property(srna, "goal_min", PROP_FLOAT, PROP_FACTOR);
-  RNA_def_property_float_sdna(prop, NULL, "mingoal");
-  RNA_def_property_range(prop, 0.0f, 1.0f);
-  RNA_def_property_ui_text(
-      prop, "Goal Minimum", "Goal minimum, vertex weights are scaled to match this range");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "goal_max", PROP_FLOAT, PROP_FACTOR);
-  RNA_def_property_float_sdna(prop, NULL, "maxgoal");
-  RNA_def_property_range(prop, 0.0f, 1.0f);
-  RNA_def_property_ui_text(
-      prop, "Goal Maximum", "Goal maximum, vertex weights are scaled to match this range");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "goal_default", PROP_FLOAT, PROP_FACTOR);
-  RNA_def_property_float_sdna(prop, NULL, "defgoal");
-  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
-  RNA_def_property_range(prop, 0.0f, 1.0f);
-  RNA_def_property_ui_text(prop, "Goal Default", "Default Goal (vertex target position) value");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "goal_spring", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "goalspring");
-  RNA_def_property_range(prop, 0.0f, 0.999f);
-  RNA_def_property_ui_text(
-      prop, "Goal Stiffness", "Goal (vertex target position) spring stiffness");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "goal_friction", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "goalfrict");
-  RNA_def_property_range(prop, 0.0f, 50.0f);
-  RNA_def_property_ui_text(prop, "Goal Damping", "Goal (vertex target position) friction");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  /* Edge Spring Settings */
-
-  prop = RNA_def_property(srna, "pull", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "inspring");
-  RNA_def_property_range(prop, 0.0f, 0.999f);
-  RNA_def_property_ui_text(prop, "Pull", "Edge spring stiffness when longer than rest length");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "push", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "inpush");
-  RNA_def_property_range(prop, 0.0f, 0.999f);
-  RNA_def_property_ui_text(prop, "Push", "Edge spring stiffness when shorter than rest length");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "damping", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "infrict");
-  RNA_def_property_range(prop, 0.0f, 50.0f);
-  RNA_def_property_ui_text(prop, "Damp", "Edge spring friction");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "spring_length", PROP_INT, PROP_NONE);
-  RNA_def_property_int_sdna(prop, NULL, "springpreload");
-  RNA_def_property_range(prop, 0.0f, 200.0f);
-  RNA_def_property_ui_text(
-      prop, "View Layer", "Alter spring length to shrink/blow up (unit %) 0 to disable");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "aero", PROP_INT, PROP_NONE);
-  RNA_def_property_int_sdna(prop, NULL, "aeroedge");
-  RNA_def_property_range(prop, 0.0f, 30000.0f);
-  RNA_def_property_ui_text(prop, "Aero", "Make edges 'sail'");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "plastic", PROP_INT, PROP_NONE);
-  RNA_def_property_int_sdna(prop, NULL, "plastic");
-  RNA_def_property_range(prop, 0.0f, 100.0f);
-  RNA_def_property_ui_text(prop, "Plasticity", "Permanent deform");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "bend", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "secondspring");
-  RNA_def_property_range(prop, 0.0f, 10.0f);
-  RNA_def_property_ui_text(prop, "Bending", "Bending Stiffness");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "shear", PROP_FLOAT, PROP_FACTOR);
-  RNA_def_property_float_sdna(prop, NULL, "shearstiff");
-  RNA_def_property_range(prop, 0.0f, 1.0f);
-  RNA_def_property_ui_text(prop, "Shear", "Shear Stiffness");
-
-  prop = RNA_def_property(srna, "vertex_group_spring", PROP_STRING, PROP_NONE);
-  RNA_def_property_string_sdna(prop, NULL, "namedVG_Spring_K");
-  RNA_def_property_ui_text(prop, "Spring Vertex Group", "Control point spring strength values");
-  RNA_def_property_string_funcs(prop, NULL, NULL, "rna_SoftBodySettings_spring_vgroup_set");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  /* Collision */
-
-  prop = RNA_def_property(srna, "collision_type", PROP_ENUM, PROP_NONE);
-  RNA_def_property_enum_sdna(prop, NULL, "sbc_mode");
-  RNA_def_property_enum_items(prop, collision_type_items);
-  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
-  RNA_def_property_ui_text(prop, "Collision Type", "Choose Collision Type");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "ball_size", PROP_FLOAT, PROP_DISTANCE);
-  RNA_def_property_float_sdna(prop, NULL, "colball");
-  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE); /* code is not ready for that yet */
-  RNA_def_property_range(prop, -10.0f, 10.0f);
-  RNA_def_property_ui_text(
-      prop, "Ball Size", "Absolute ball size or factor if not manually adjusted");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "ball_stiff", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "ballstiff");
-  RNA_def_property_range(prop, 0.001f, 100.0f);
-  RNA_def_property_ui_text(prop, "Ball Size", "Ball inflating pressure");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
-
-  prop = RNA_def_property(srna, "ball_damp", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "balldamp");
+  prop = RNA_def_property(srna, "dt", PROP_FLOAT, PROP_NONE);
+  RNA_def_property_float_sdna(prop, NULL, "dt");
   RNA_def_property_range(prop, 0.001f, 1.0f);
-  RNA_def_property_ui_text(prop, "Ball Size", "Blending to inelastic collision");
+  RNA_def_property_ui_text(prop, "dt", "Time duration of a frame");
   RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  /* Solver */
-
-  prop = RNA_def_property(srna, "error_threshold", PROP_FLOAT, PROP_NONE);
-  RNA_def_property_float_sdna(prop, NULL, "rklimit");
-  RNA_def_property_range(prop, 0.001f, 10.0f);
-  RNA_def_property_ui_text(
-      prop,
-      "Error Limit",
-      "The Runge-Kutta ODE solver error limit, low value gives more precision, "
-      "high values speed");
+  prop = RNA_def_property(srna, "substep_count", PROP_FLOAT, PROP_NONE);
+  RNA_def_property_float_sdna(prop, NULL, "substep_count");
+  RNA_def_property_range(prop, 1.0f, 5000.0f);
+  RNA_def_property_ui_text(prop, "Substep Count", "Number of substeps in each step");
   RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "step_min", PROP_INT, PROP_NONE);
-  RNA_def_property_int_sdna(prop, NULL, "minloops");
-  RNA_def_property_range(prop, 0, 30000);
-  RNA_def_property_ui_text(prop, "Min Step", "Minimal # solver steps/frame");
+  prop = RNA_def_property(srna, "alpha_vol", PROP_FLOAT, PROP_NONE);
+  RNA_def_property_float_sdna(prop, NULL, "alpha_vol");
+  RNA_def_property_range(prop, 0.0f, 100.0f);
+  RNA_def_property_ui_text(prop, "Volume Compliance", "Volume stiffness");
   RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "step_max", PROP_INT, PROP_NONE);
-  RNA_def_property_int_sdna(prop, NULL, "maxloops");
-  RNA_def_property_range(prop, 0, 30000);
-  RNA_def_property_ui_text(prop, "Max Step", "Maximal # solver steps/frame");
+  prop = RNA_def_property(srna, "alpha_edge", PROP_FLOAT, PROP_NONE);
+  RNA_def_property_float_sdna(prop, NULL, "alpha_edge");
+  RNA_def_property_range(prop, 0.0f, 100.0f);
+  RNA_def_property_ui_text(prop, "Edge Compliance", "Edge stiffness");
   RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "choke", PROP_INT, PROP_NONE);
-  RNA_def_property_int_sdna(prop, NULL, "choke");
-  RNA_def_property_range(prop, 0, 100);
-  RNA_def_property_ui_text(prop, "Choke", "'Viscosity' inside collision target");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  /* General Settings */
+  
+  // prop = RNA_def_property(srna, "friction", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "mediafrict");
+  // RNA_def_property_range(prop, 0.0f, 50.0f);
+  // RNA_def_property_ui_text(prop, "Friction", "General media friction for point movements");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "fuzzy", PROP_INT, PROP_NONE);
-  RNA_def_property_int_sdna(prop, NULL, "fuzzyness");
-  RNA_def_property_range(prop, 1, 100);
-  RNA_def_property_ui_text(
-      prop,
-      "Fuzzy",
-      "Fuzziness while on collision, high values make collision handling faster "
-      "but less stable");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  // prop = RNA_def_property(srna, "mass", PROP_FLOAT, PROP_UNIT_MASS);
+  // RNA_def_property_float_sdna(prop, NULL, "nodemass");
+  // RNA_def_property_range(prop, 0.0f, 50000.0f);
+  // RNA_def_property_ui_text(prop, "Mass", "General Mass value");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "use_auto_step", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_sdna(prop, NULL, "solverflags", SBSO_OLDERR);
-  RNA_def_property_ui_text(prop, "V", "Use velocities for automagic step sizes");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  // prop = RNA_def_property(srna, "vertex_group_mass", PROP_STRING, PROP_NONE);
+  // RNA_def_property_string_sdna(prop, NULL, "namedVG_Mass");
+  // RNA_def_property_ui_text(prop, "Mass Vertex Group", "Control point mass values");
+  // RNA_def_property_string_funcs(prop, NULL, NULL, "rna_SoftBodySettings_mass_vgroup_set");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "use_diagnose", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_sdna(prop, NULL, "solverflags", SBSO_MONITOR);
-  RNA_def_property_ui_text(
-      prop, "Print Performance to Console", "Turn on SB diagnose console prints");
+  // /* no longer used */
+  // prop = RNA_def_property(srna, "gravity", PROP_FLOAT, PROP_ACCELERATION);
+  // RNA_def_property_float_sdna(prop, NULL, "grav");
+  // RNA_def_property_range(prop, -10.0f, 10.0f);
+  // RNA_def_property_ui_text(prop, "Gravitation", "Apply gravitation to point movement");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "use_estimate_matrix", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_sdna(prop, NULL, "solverflags", SBSO_ESTIMATEIPO);
-  RNA_def_property_ui_text(
-      prop, "Estimate Transforms", "Store the estimated transforms in the soft body settings");
+  // prop = RNA_def_property(srna, "speed", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "physics_speed");
+  // RNA_def_property_range(prop, 0.01f, 100.0f);
+  // RNA_def_property_ui_text(
+  //     prop, "Speed", "Tweak timing for physics to control frequency and speed");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  /***********************************************************************************/
-  /* These are not exactly settings, but reading calculated results
-   * but i did not want to start a new property struct
-   * so rather rename this from SoftBodySettings to SoftBody
-   * translation. */
-  prop = RNA_def_property(srna, "location_mass_center", PROP_FLOAT, PROP_TRANSLATION);
-  RNA_def_property_float_sdna(prop, NULL, "lcom");
-  RNA_def_property_ui_text(prop, "Center of Mass", "Location of center of mass");
-  RNA_def_property_ui_range(prop, -FLT_MAX, FLT_MAX, 1, RNA_TRANSLATION_PREC_DEFAULT);
+  // /* Goal */
 
-  /* matrix */
-  prop = RNA_def_property(srna, "rotation_estimate", PROP_FLOAT, PROP_MATRIX);
-  RNA_def_property_float_sdna(prop, NULL, "lrot");
-  RNA_def_property_multi_array(prop, 2, rna_matrix_dimsize_3x3);
-  RNA_def_property_ui_text(prop, "Rotation Matrix", "Estimated rotation matrix");
+  // prop = RNA_def_property(srna, "vertex_group_goal", PROP_STRING, PROP_NONE);
+  // RNA_def_property_string_sdna(prop, NULL, "vertgroup");
+  // RNA_def_property_clear_flag(prop, PROP_ANIMATABLE); /* not impossible .. but not supported yet */
+  // RNA_def_property_string_funcs(prop,
+  //                               "rna_SoftBodySettings_goal_vgroup_get",
+  //                               "rna_SoftBodySettings_goal_vgroup_length",
+  //                               "rna_SoftBodySettings_goal_vgroup_set");
+  // RNA_def_property_ui_text(prop, "Goal Vertex Group", "Control point weight values");
 
-  prop = RNA_def_property(srna, "scale_estimate", PROP_FLOAT, PROP_MATRIX);
-  RNA_def_property_float_sdna(prop, NULL, "lscale");
-  RNA_def_property_multi_array(prop, 2, rna_matrix_dimsize_3x3);
-  RNA_def_property_ui_text(prop, "Scale Matrix", "Estimated scale matrix");
-  /***********************************************************************************/
+  // prop = RNA_def_property(srna, "goal_min", PROP_FLOAT, PROP_FACTOR);
+  // RNA_def_property_float_sdna(prop, NULL, "mingoal");
+  // RNA_def_property_range(prop, 0.0f, 1.0f);
+  // RNA_def_property_ui_text(
+  //     prop, "Goal Minimum", "Goal minimum, vertex weights are scaled to match this range");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  /* Flags */
+  // prop = RNA_def_property(srna, "goal_max", PROP_FLOAT, PROP_FACTOR);
+  // RNA_def_property_float_sdna(prop, NULL, "maxgoal");
+  // RNA_def_property_range(prop, 0.0f, 1.0f);
+  // RNA_def_property_ui_text(
+  //     prop, "Goal Maximum", "Goal maximum, vertex weights are scaled to match this range");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "use_goal", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_funcs(
-      prop, "rna_SoftBodySettings_use_goal_get", "rna_SoftBodySettings_use_goal_set");
-  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
-  RNA_def_property_ui_text(
-      prop, "Use Goal", "Define forces for vertices to stick to animated position");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  // prop = RNA_def_property(srna, "goal_default", PROP_FLOAT, PROP_FACTOR);
+  // RNA_def_property_float_sdna(prop, NULL, "defgoal");
+  // RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  // RNA_def_property_range(prop, 0.0f, 1.0f);
+  // RNA_def_property_ui_text(prop, "Goal Default", "Default Goal (vertex target position) value");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "use_edges", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_funcs(
-      prop, "rna_SoftBodySettings_use_edges_get", "rna_SoftBodySettings_use_edges_set");
-  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
-  RNA_def_property_ui_text(prop, "Use Edges", "Use Edges as springs");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  // prop = RNA_def_property(srna, "goal_spring", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "goalspring");
+  // RNA_def_property_range(prop, 0.0f, 0.999f);
+  // RNA_def_property_ui_text(
+  //     prop, "Goal Stiffness", "Goal (vertex target position) spring stiffness");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "use_stiff_quads", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_funcs(
-      prop, "rna_SoftBodySettings_stiff_quads_get", "rna_SoftBodySettings_stiff_quads_set");
-  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
-  RNA_def_property_ui_text(prop, "Stiff Quads", "Add diagonal springs on 4-gons");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  // prop = RNA_def_property(srna, "goal_friction", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "goalfrict");
+  // RNA_def_property_range(prop, 0.0f, 50.0f);
+  // RNA_def_property_ui_text(prop, "Goal Damping", "Goal (vertex target position) friction");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "use_edge_collision", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_funcs(
-      prop, "rna_SoftBodySettings_edge_collision_get", "rna_SoftBodySettings_edge_collision_set");
-  RNA_def_property_ui_text(prop, "Edge Collision", "Edges collide too");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  // /* Edge Spring Settings */
 
-  prop = RNA_def_property(srna, "use_face_collision", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_funcs(
-      prop, "rna_SoftBodySettings_face_collision_get", "rna_SoftBodySettings_face_collision_set");
-  RNA_def_property_ui_text(prop, "Face Collision", "Faces collide too, can be very slow");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  // prop = RNA_def_property(srna, "pull", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "inspring");
+  // RNA_def_property_range(prop, 0.0f, 0.999f);
+  // RNA_def_property_ui_text(prop, "Pull", "Edge spring stiffness when longer than rest length");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "aerodynamics_type", PROP_ENUM, PROP_NONE);
-  RNA_def_property_enum_items(prop, aerodynamics_type);
-  RNA_def_property_enum_funcs(
-      prop, "rna_SoftBodySettings_new_aero_get", "rna_SoftBodySettings_new_aero_set", NULL);
-  RNA_def_property_ui_text(
-      prop, "Aerodynamics Type", "Method of calculating aerodynamic interaction");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  // prop = RNA_def_property(srna, "push", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "inpush");
+  // RNA_def_property_range(prop, 0.0f, 0.999f);
+  // RNA_def_property_ui_text(prop, "Push", "Edge spring stiffness when shorter than rest length");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "use_self_collision", PROP_BOOLEAN, PROP_NONE);
-  RNA_def_property_boolean_funcs(
-      prop, "rna_SoftBodySettings_self_collision_get", "rna_SoftBodySettings_self_collision_set");
-  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
-  RNA_def_property_ui_text(prop, "Self Collision", "Enable naive vertex ball self collision");
-  RNA_def_property_update(prop, 0, "rna_softbody_update");
+  // prop = RNA_def_property(srna, "damping", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "infrict");
+  // RNA_def_property_range(prop, 0.0f, 50.0f);
+  // RNA_def_property_ui_text(prop, "Damp", "Edge spring friction");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "collision_collection", PROP_POINTER, PROP_NONE);
-  RNA_def_property_struct_type(prop, "Collection");
-  RNA_def_property_pointer_sdna(prop, NULL, "collision_group");
-  RNA_def_property_flag(prop, PROP_EDITABLE);
-  RNA_def_property_ui_text(prop, "Collision Collection", "Limit colliders to this collection");
-  RNA_def_property_update(prop, 0, "rna_softbody_dependency_update");
+  // prop = RNA_def_property(srna, "spring_length", PROP_INT, PROP_NONE);
+  // RNA_def_property_int_sdna(prop, NULL, "springpreload");
+  // RNA_def_property_range(prop, 0.0f, 200.0f);
+  // RNA_def_property_ui_text(
+  //     prop, "View Layer", "Alter spring length to shrink/blow up (unit %) 0 to disable");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
 
-  prop = RNA_def_property(srna, "effector_weights", PROP_POINTER, PROP_NONE);
-  RNA_def_property_pointer_sdna(prop, NULL, "effector_weights");
-  RNA_def_property_struct_type(prop, "EffectorWeights");
-  RNA_def_property_clear_flag(prop, PROP_EDITABLE);
-  RNA_def_property_override_flag(prop, PROPOVERRIDE_OVERRIDABLE_LIBRARY);
-  RNA_def_property_ui_text(prop, "Effector Weights", "");
+  // prop = RNA_def_property(srna, "aero", PROP_INT, PROP_NONE);
+  // RNA_def_property_int_sdna(prop, NULL, "aeroedge");
+  // RNA_def_property_range(prop, 0.0f, 30000.0f);
+  // RNA_def_property_ui_text(prop, "Aero", "Make edges 'sail'");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "plastic", PROP_INT, PROP_NONE);
+  // RNA_def_property_int_sdna(prop, NULL, "plastic");
+  // RNA_def_property_range(prop, 0.0f, 100.0f);
+  // RNA_def_property_ui_text(prop, "Plasticity", "Permanent deform");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "bend", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "secondspring");
+  // RNA_def_property_range(prop, 0.0f, 10.0f);
+  // RNA_def_property_ui_text(prop, "Bending", "Bending Stiffness");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "shear", PROP_FLOAT, PROP_FACTOR);
+  // RNA_def_property_float_sdna(prop, NULL, "shearstiff");
+  // RNA_def_property_range(prop, 0.0f, 1.0f);
+  // RNA_def_property_ui_text(prop, "Shear", "Shear Stiffness");
+
+  // prop = RNA_def_property(srna, "vertex_group_spring", PROP_STRING, PROP_NONE);
+  // RNA_def_property_string_sdna(prop, NULL, "namedVG_Spring_K");
+  // RNA_def_property_ui_text(prop, "Spring Vertex Group", "Control point spring strength values");
+  // RNA_def_property_string_funcs(prop, NULL, NULL, "rna_SoftBodySettings_spring_vgroup_set");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // /* Collision */
+
+  // prop = RNA_def_property(srna, "collision_type", PROP_ENUM, PROP_NONE);
+  // RNA_def_property_enum_sdna(prop, NULL, "sbc_mode");
+  // RNA_def_property_enum_items(prop, collision_type_items);
+  // RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  // RNA_def_property_ui_text(prop, "Collision Type", "Choose Collision Type");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "ball_size", PROP_FLOAT, PROP_DISTANCE);
+  // RNA_def_property_float_sdna(prop, NULL, "colball");
+  // RNA_def_property_clear_flag(prop, PROP_ANIMATABLE); /* code is not ready for that yet */
+  // RNA_def_property_range(prop, -10.0f, 10.0f);
+  // RNA_def_property_ui_text(
+  //     prop, "Ball Size", "Absolute ball size or factor if not manually adjusted");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "ball_stiff", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "ballstiff");
+  // RNA_def_property_range(prop, 0.001f, 100.0f);
+  // RNA_def_property_ui_text(prop, "Ball Size", "Ball inflating pressure");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "ball_damp", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "balldamp");
+  // RNA_def_property_range(prop, 0.001f, 1.0f);
+  // RNA_def_property_ui_text(prop, "Ball Size", "Blending to inelastic collision");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // /* Solver */
+
+  // prop = RNA_def_property(srna, "error_threshold", PROP_FLOAT, PROP_NONE);
+  // RNA_def_property_float_sdna(prop, NULL, "rklimit");
+  // RNA_def_property_range(prop, 0.001f, 10.0f);
+  // RNA_def_property_ui_text(
+  //     prop,
+  //     "Error Limit",
+  //     "The Runge-Kutta ODE solver error limit, low value gives more precision, "
+  //     "high values speed");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "step_min", PROP_INT, PROP_NONE);
+  // RNA_def_property_int_sdna(prop, NULL, "minloops");
+  // RNA_def_property_range(prop, 0, 30000);
+  // RNA_def_property_ui_text(prop, "Min Step", "Minimal # solver steps/frame");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "step_max", PROP_INT, PROP_NONE);
+  // RNA_def_property_int_sdna(prop, NULL, "maxloops");
+  // RNA_def_property_range(prop, 0, 30000);
+  // RNA_def_property_ui_text(prop, "Max Step", "Maximal # solver steps/frame");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "choke", PROP_INT, PROP_NONE);
+  // RNA_def_property_int_sdna(prop, NULL, "choke");
+  // RNA_def_property_range(prop, 0, 100);
+  // RNA_def_property_ui_text(prop, "Choke", "'Viscosity' inside collision target");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "fuzzy", PROP_INT, PROP_NONE);
+  // RNA_def_property_int_sdna(prop, NULL, "fuzzyness");
+  // RNA_def_property_range(prop, 1, 100);
+  // RNA_def_property_ui_text(
+  //     prop,
+  //     "Fuzzy",
+  //     "Fuzziness while on collision, high values make collision handling faster "
+  //     "but less stable");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "use_auto_step", PROP_BOOLEAN, PROP_NONE);
+  // RNA_def_property_boolean_sdna(prop, NULL, "solverflags", SBSO_OLDERR);
+  // RNA_def_property_ui_text(prop, "V", "Use velocities for automagic step sizes");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "use_diagnose", PROP_BOOLEAN, PROP_NONE);
+  // RNA_def_property_boolean_sdna(prop, NULL, "solverflags", SBSO_MONITOR);
+  // RNA_def_property_ui_text(
+  //     prop, "Print Performance to Console", "Turn on SB diagnose console prints");
+
+  // prop = RNA_def_property(srna, "use_estimate_matrix", PROP_BOOLEAN, PROP_NONE);
+  // RNA_def_property_boolean_sdna(prop, NULL, "solverflags", SBSO_ESTIMATEIPO);
+  // RNA_def_property_ui_text(
+  //     prop, "Estimate Transforms", "Store the estimated transforms in the soft body settings");
+
+  // /***********************************************************************************/
+  // /* These are not exactly settings, but reading calculated results
+  //  * but i did not want to start a new property struct
+  //  * so rather rename this from SoftBodySettings to SoftBody
+  //  * translation. */
+  // prop = RNA_def_property(srna, "location_mass_center", PROP_FLOAT, PROP_TRANSLATION);
+  // RNA_def_property_float_sdna(prop, NULL, "lcom");
+  // RNA_def_property_ui_text(prop, "Center of Mass", "Location of center of mass");
+  // RNA_def_property_ui_range(prop, -FLT_MAX, FLT_MAX, 1, RNA_TRANSLATION_PREC_DEFAULT);
+
+  // /* matrix */
+  // prop = RNA_def_property(srna, "rotation_estimate", PROP_FLOAT, PROP_MATRIX);
+  // RNA_def_property_float_sdna(prop, NULL, "lrot");
+  // RNA_def_property_multi_array(prop, 2, rna_matrix_dimsize_3x3);
+  // RNA_def_property_ui_text(prop, "Rotation Matrix", "Estimated rotation matrix");
+
+  // prop = RNA_def_property(srna, "scale_estimate", PROP_FLOAT, PROP_MATRIX);
+  // RNA_def_property_float_sdna(prop, NULL, "lscale");
+  // RNA_def_property_multi_array(prop, 2, rna_matrix_dimsize_3x3);
+  // RNA_def_property_ui_text(prop, "Scale Matrix", "Estimated scale matrix");
+  // /***********************************************************************************/
+
+  // /* Flags */
+
+  // prop = RNA_def_property(srna, "use_goal", PROP_BOOLEAN, PROP_NONE);
+  // RNA_def_property_boolean_funcs(
+  //     prop, "rna_SoftBodySettings_use_goal_get", "rna_SoftBodySettings_use_goal_set");
+  // RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  // RNA_def_property_ui_text(
+  //     prop, "Use Goal", "Define forces for vertices to stick to animated position");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "use_edges", PROP_BOOLEAN, PROP_NONE);
+  // RNA_def_property_boolean_funcs(
+  //     prop, "rna_SoftBodySettings_use_edges_get", "rna_SoftBodySettings_use_edges_set");
+  // RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  // RNA_def_property_ui_text(prop, "Use Edges", "Use Edges as springs");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "use_stiff_quads", PROP_BOOLEAN, PROP_NONE);
+  // RNA_def_property_boolean_funcs(
+  //     prop, "rna_SoftBodySettings_stiff_quads_get", "rna_SoftBodySettings_stiff_quads_set");
+  // RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  // RNA_def_property_ui_text(prop, "Stiff Quads", "Add diagonal springs on 4-gons");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "use_edge_collision", PROP_BOOLEAN, PROP_NONE);
+  // RNA_def_property_boolean_funcs(
+  //     prop, "rna_SoftBodySettings_edge_collision_get", "rna_SoftBodySettings_edge_collision_set");
+  // RNA_def_property_ui_text(prop, "Edge Collision", "Edges collide too");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "use_face_collision", PROP_BOOLEAN, PROP_NONE);
+  // RNA_def_property_boolean_funcs(
+  //     prop, "rna_SoftBodySettings_face_collision_get", "rna_SoftBodySettings_face_collision_set");
+  // RNA_def_property_ui_text(prop, "Face Collision", "Faces collide too, can be very slow");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "aerodynamics_type", PROP_ENUM, PROP_NONE);
+  // RNA_def_property_enum_items(prop, aerodynamics_type);
+  // RNA_def_property_enum_funcs(
+  //     prop, "rna_SoftBodySettings_new_aero_get", "rna_SoftBodySettings_new_aero_set", NULL);
+  // RNA_def_property_ui_text(
+  //     prop, "Aerodynamics Type", "Method of calculating aerodynamic interaction");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "use_self_collision", PROP_BOOLEAN, PROP_NONE);
+  // RNA_def_property_boolean_funcs(
+  //     prop, "rna_SoftBodySettings_self_collision_get", "rna_SoftBodySettings_self_collision_set");
+  // RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  // RNA_def_property_ui_text(prop, "Self Collision", "Enable naive vertex ball self collision");
+  // RNA_def_property_update(prop, 0, "rna_softbody_update");
+
+  // prop = RNA_def_property(srna, "collision_collection", PROP_POINTER, PROP_NONE);
+  // RNA_def_property_struct_type(prop, "Collection");
+  // RNA_def_property_pointer_sdna(prop, NULL, "collision_group");
+  // RNA_def_property_flag(prop, PROP_EDITABLE);
+  // RNA_def_property_ui_text(prop, "Collision Collection", "Limit colliders to this collection");
+  // RNA_def_property_update(prop, 0, "rna_softbody_dependency_update");
+
+  // prop = RNA_def_property(srna, "effector_weights", PROP_POINTER, PROP_NONE);
+  // RNA_def_property_pointer_sdna(prop, NULL, "effector_weights");
+  // RNA_def_property_struct_type(prop, "EffectorWeights");
+  // RNA_def_property_clear_flag(prop, PROP_EDITABLE);
+  // RNA_def_property_override_flag(prop, PROPOVERRIDE_OVERRIDABLE_LIBRARY);
+  // RNA_def_property_ui_text(prop, "Effector Weights", "");
 }
 
 void RNA_def_object_force(BlenderRNA *brna)
 {
   rna_def_pointcache_active(brna);
   rna_def_collision(brna);
-  rna_def_effector_weight(brna);
+  rna_def_effector_weight(brna); 
   rna_def_field(brna);
   rna_def_softbody(brna);
 }

source/blender/modifiers/intern/MOD_softbody.c

@@ -54,17 +54,17 @@ static bool dependsOnTime(struct Scene *UNUSED(scene), ModifierData *UNUSED(md))
 
 static void updateDepsgraph(ModifierData *UNUSED(md), const ModifierUpdateDepsgraphContext *ctx)
 {
-  if (ctx->object->soft) {
-    /* Actual code uses ccd_build_deflector_hash */
-    DEG_add_collision_relations(ctx->node,
-                                ctx->object,
-                                ctx->object->soft->collision_group,
-                                eModifierType_Collision,
-                                NULL,
-                                "Softbody Collision");
-    DEG_add_forcefield_relations(
-        ctx->node, ctx->object, ctx->object->soft->effector_weights, true, 0, "Softbody Field");
-  }
+  // if (ctx->object->soft) {
+  //   /* Actual code uses ccd_build_deflector_hash */
+  //   DEG_add_collision_relations(ctx->node,
+  //                               ctx->object,
+  //                               ctx->object->soft->collision_group,
+  //                               eModifierType_Collision,
+  //                               NULL,
+  //                               "Softbody Collision");
+  //   DEG_add_forcefield_relations(
+  //       ctx->node, ctx->object, ctx->object->soft->effector_weights, true, 0, "Softbody Field");
+  // }
   /* We need own transformation as well. */
   DEG_add_modifier_to_transform_relation(ctx->node, "SoftBody Modifier");
 }

source/blender/modifiers/intern/MOD_weld copy.cc

@@ -0,0 +1,811 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright 2005 Blender Foundation. All rights reserved. */
+
+/** \file
+ * \ingroup modifiers
+ *
+ * Weld modifier: Remove doubles.
+ */
+
+/* TODOs:
+ * - Review weight and vertex color interpolation.;
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_utildefines.h"
+
+#include "BLI_array.hh"
+#include "BLI_index_range.hh"
+#include "BLI_span.hh"
+
+#include "BLI_math_vector.h"
+#include "BLI_utildefines.h"
+
+#include "BLI_vector.hh"
+#include "BLI_math_vec_types.hh"
+#include "BLI_math_vector.hh"
+
+// #include "bmesh_construct.h"
+#include "bmesh.h"
+#include "bmesh_tools.h"
+
+#include "BLT_translation.h"
+
+#include "DNA_defaults.h"
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_modifier_types.h"
+#include "DNA_screen_types.h"
+
+#include "BKE_bvhutils.h"
+#include "BKE_context.h"
+#include "BKE_deform.h"
+#include "BKE_modifier.h"
+#include "BKE_screen.h"
+#include "BKE_mesh.h"
+
+#include "UI_interface.h"
+#include "UI_resources.h"
+
+#include "RNA_access.h"
+#include "RNA_prototypes.h"
+
+#include "DEG_depsgraph.h"
+
+#include "MOD_modifiertypes.h"
+#include "MOD_ui_common.h"
+
+#include "GEO_mesh_merge_by_distance.hh"
+
+using blender::Array;
+using blender::IndexMask;
+using blender::Span;
+using blender::Vector;
+using namespace blender;
+using namespace blender::math;
+// using namespace std;
+
+
+Vector<Vector<int>> tetFaces({{2,1,0}, {0,1,3}, {1,2,3}, {2,0,3}});
+float directions[6][3] = {{1,0,0}, {-1,0,0}, {0,1,0}, {0,-1,0}, {0,0,1}, {0,0,-1}};
+
+float inf = FLT_MAX;
+float eps = 0.0f;
+
+static float randomEps(){
+  float eps = 0.0001 ;
+  return -eps + 2.0 * (static_cast <float> (rand()) / static_cast <float> (RAND_MAX)) * eps;
+  // return 0.0f;
+}
+
+static bool isInside(float3 vert, BVHTreeFromMesh *treedata){
+  // float eps = 0.0001;
+
+  int count =  0;
+  float min_dist = 0.0f;
+  for(auto dir : directions){
+    float radius = 0.0f;
+    
+    float max_length = FLT_MAX;
+
+    BVHTreeRayHit rayhit = {0};
+    rayhit.index = -1;
+    rayhit.dist = max_length;
+    BLI_bvhtree_ray_cast(treedata->tree, vert, dir, radius, &rayhit, treedata->raycast_callback, treedata);
+
+    if (rayhit.index != -1 && rayhit.dist <= max_length) {
+      if(dot_v3v3(rayhit.no, dir) > 0.0f){
+        count++;
+      }
+
+      if((min_dist > 0.0) && (min_dist - rayhit.dist) > 0.0)
+        return false;
+    }
+  }
+
+  return count > 3;
+}
+
+// Used to populate a tets face normals and planesD
+static void setTetProperties(Vector<float3> &verts, 
+                      Vector<int> &tetVertId,
+                      Vector<float3> &faceNormals, 
+                      Vector<float> &planesD,
+                      int tetNr){
+  for(int i = 0; i<4; i++){ 
+    float3 p0 = verts[tetVertId[4*tetNr + tetFaces[i][0]]];
+    float3 p1 = verts[tetVertId[4*tetNr + tetFaces[i][1]]];
+    float3 p2 = verts[tetVertId[4*tetNr + tetFaces[i][2]]];
+
+    float3 normal = cross(p1 - p0, p2 - p0);
+    normal = normalize(normal);
+
+    faceNormals[4*tetNr + i] = normal;
+    planesD[4*tetNr + i] = dot(p0, normal);
+  }
+}
+
+static float3 getCircumCenter(float3 p0, float3 p1, float3 p2, float3 p3){
+  float3 b = p1 - p0;
+  float3 c = p2 - p0;
+  float3 d = p3 - p0;
+
+  float det = 2.0 * (b.x*(c.y*d.z - c.z*d.y) - b.y*(c.x*d.z - c.z*d.x) + b.z*(c.x*d.y - c.y*d.x));
+  if (det == 0.0f){
+    return p0;
+  }
+  else{
+    float3 v = cross(c, d)*dot(b, b) + cross(d, b)*dot(c, c) + cross(b, c)*dot(d, d);
+    v /= det;
+    return p0 + v;
+  }
+      
+}
+
+// bool isSameSide(float3 p0, float3 p1, float3 p2, float3 p4, float3 vert){
+
+// }
+
+// bool isInsideTet(float3 p0, float3 p1, float3 p2, float3 p3, float3 vert){
+
+// }
+
+static int findContainingTet(Vector<float3> &verts, 
+                    Vector<int> &tetVertId, 
+                    Vector<int> &tetFaceNeighbors, 
+                    Vector<float3> &faceNormals, 
+                    Vector<float> &planesD, 
+                    int tetMarkId, 
+                    Vector<int> &tetMarks,
+                    float3 currVert){
+  
+  /* --------------------
+  Matthias Method
+  ----------------------*/
+
+  // bool found = false;
+  // int tetNr = 0;
+
+  // while(tetNr < tetVertId.size()/4 && tetVertId[4*tetNr]<0)
+  //   tetNr++;
+
+  // float3 center(0.0, 0.0, 0.0);
+  // while(!found){
+  //   if(tetNr < 0 || tetMarks[tetNr] == tetMarkId){
+  //     break;
+  //   }
+    
+  //   tetMarks[tetNr] = tetMarkId;
+
+  //   center = {0.0,0.0,0.0};
+  //   for(int i = 0; i<4; i++){
+  //     center += verts[tetVertId[4*tetNr + i]];
+  //   }
+  //   center*=0.25;
+
+  //   float minT = inf; //
+  //   int minFaceNr = -1;
+
+  //   for(int i = 0; i<4; i++){
+  //     float3 normal = faceNormals[4*tetNr + i];
+  //     float d = planesD[4*tetNr + i];
+
+  //     float hp = dot(normal, currVert) - d;
+  //     float hc = dot(normal, center) - d;
+
+  //     float t = hp - hc;
+  //     if(t <= eps)
+  //       continue;
+
+  //     t = -hc/t;
+
+  //     if((t >= eps) && ((t - minT) > eps)){
+  //       minT = t;
+  //       minFaceNr = i;
+  //     }
+  //   }
+
+  //   if(minT >= 1.0){
+  //     found = true;
+  //   }
+  //   else{
+  //     tetNr = tetFaceNeighbors[4*tetNr + minFaceNr];
+  //   }
+  // }
+
+  // if(found)
+  //   return tetNr;
+  // else
+  //   return -1;
+
+  /* --------------------
+  Brute force finding first violating tet
+  ----------------------*/
+
+  for(int currTet = 0; currTet<tetVertId.size()/4; currTet++){
+    if(tetVertId[4*currTet] == -1)
+      continue;
+    
+    float3 p0 = verts[tetVertId[4*currTet + 0]];
+    float3 p1 = verts[tetVertId[4*currTet + 1]];
+    float3 p2 = verts[tetVertId[4*currTet + 2]];
+    float3 p3 = verts[tetVertId[4*currTet + 3]];
+
+    float3 circumCenter = getCircumCenter(p0, p1, p2, p3);
+    float circumRadius = length(p0 - circumCenter);
+    
+    // if((circumRadius - length(currVert - circumCenter)) >= eps){
+    if(length(currVert - circumCenter) < circumRadius){
+      return currTet;
+    }
+  }
+
+  /* -----------------------
+  My method, checks if point is inside tet, if not finds the face that is connecting to the point and center
+  ------------------------*/
+
+  // bool found = false;
+  // int tetNr = 0;
+  // while(tetNr < tetVertId.size()/4 && tetVertId[4*tetNr] < 0)
+  //   tetNr++;
+
+  // while(!found){
+  //   float3 p0 = verts[tetVertId[4*currTet + 0]];
+  //   float3 p1 = verts[tetVertId[4*currTet + 1]];
+  //   float3 p2 = verts[tetVertId[4*currTet + 2]];
+  //   float3 p3 = verts[tetVertId[4*currTet + 3]];
+    
+  //   if(isInsideTet(p0, p1, p2, p3, currVert))
+  //     return tetNr;
+    
+  //   for(int i = 0; i<4; i++)
+  // }
+
+  return -1;
+}
+
+/*
+The basic assumption employed is that 2 violating tets cannot not be neighbors. 
+Simple BFS approach that checks neighbors of all violating tets and adds them to stack
+
+If a violating tet shares a face with a non-violating tet, that's a boundary face and for each violating tet, list of its boundary faces is returned
+
+Note - BFS can be written better
+*/
+static Vector<std::pair<int, Vector<int>>> getViolatingTets(Vector<float3> &verts, 
+                            Vector<int> &tetVertId, 
+                            Vector<int> &tetFaceNeighbors, 
+                            int tetMarkId, 
+                            Vector<int> &tetMarks,
+                            float3 currVert,
+                            int containingTetNr
+                            ){
+  
+  Vector< std::pair<int,Vector<int>> > violatingTets;
+  Vector<int> stack;
+
+  stack.append(containingTetNr);  
+
+  while(stack.size()){
+    int currTet = stack.last();
+    stack.remove_last();
+
+    if(tetMarks[currTet] == tetMarkId){
+      continue;
+    }
+    tetMarks[currTet] = tetMarkId; 
+    Vector<int> currTetBorderFaces;
+
+    for(int i = 0; i<4; i++){
+      int neighborTet = tetFaceNeighbors[4*currTet + i];
+      
+      if(neighborTet<0){
+        currTetBorderFaces.append(i);
+        continue;
+      }
+      if(tetMarks[neighborTet]==tetMarkId){
+        continue;
+      }
+      
+      float3 p0 = verts[tetVertId[4*neighborTet + 0]];
+      float3 p1 = verts[tetVertId[4*neighborTet + 1]];
+      float3 p2 = verts[tetVertId[4*neighborTet + 2]];
+      float3 p3 = verts[tetVertId[4*neighborTet + 3]];
+
+      float3 circumCenter = getCircumCenter(p0, p1, p2, p3);
+      float circumRadius = length(p0 - circumCenter);
+
+      // if((circumRadius - length(currVert - circumCenter)) >= eps){
+      if(length(currVert - circumCenter) < circumRadius){
+        stack.append(neighborTet);
+        // tetMarks[neighborTet] = tetMarkId;
+      }
+      else{
+        currTetBorderFaces.append(i);
+      }
+    }
+    violatingTets.append({currTet, currTetBorderFaces});
+  }
+
+  return violatingTets;
+}
+
+static Vector<int> createTets(Vector<float3> verts, BVHTreeFromMesh *treedata, float minTetQuality){
+  Vector<int> tetVertId; // Stores indices of vertex that form a tet. Every tet is stored as 4 indices 
+  Vector<int> tetFaceNeighbors; // Stores index of tet that shares the face as per tetFaces order. 4 neighbors per tet, 1 for each face
+
+  int firstFreeTet = -1;
+
+  Vector<float3> faceNormals; // Stores the normal of each of face of the tet
+  Vector<float> planesD;
+
+  int tetMarkId = 0;
+  Vector<int> tetMarks; 
+  /*Used to keep track of visited tets for various processess. In each step where visited tets need to avoided, 
+  on visiting them, a unique number is stored representing that particular step
+  */
+ 
+  int bigTet = verts.size() - 4;  
+
+  for(int i = 0; i<4; i++){
+    tetVertId.append(bigTet + i);
+    tetFaceNeighbors.append(-1);
+    faceNormals.append({0.0, 0.0, 0.0});
+    planesD.append(0.0);
+  }
+  tetMarks.append(0);
+  setTetProperties(verts, tetVertId, faceNormals, planesD, 0);
+
+  for(int vertNr = 0; vertNr<bigTet; vertNr++){
+    float3 currVert = verts[vertNr];
+
+    // std::cout << "Adding vert " << vertNr << std::endl;
+
+    // Find containing tet (need to understand) center can lie outside the tet?
+    tetMarkId += 1;
+    int containingTetNr = -1;
+    containingTetNr = findContainingTet(verts, tetVertId, tetFaceNeighbors, faceNormals, planesD, tetMarkId, tetMarks, currVert);
+    
+    if(containingTetNr == -1){
+      std::cout << "Couldn't add vert " << vertNr << std::endl;
+      continue;
+    }
+
+    // Find Violating Tets - Returns all tets that are violating the Delaunay condition along with a list of face indices that are boundary faces
+    // Boundary faces may be outermost face of the structure or the face shared by a violating tet and a non violating tet
+    tetMarkId += 1;
+    Vector<std::pair<int, Vector<int>>> violatingTets = getViolatingTets(verts, tetVertId, tetFaceNeighbors, tetMarkId, tetMarks, currVert, containingTetNr);
+
+    //Create new tets centered at the new point and including the boundary faces
+    Vector<int> newTets; // Stores tet indices of the new tets formed. Used for making neighbors of the new tets formed
+    for(int violatingTetNr = 0; violatingTetNr<violatingTets.size(); violatingTetNr++){
+      int violatingTet = violatingTets[violatingTetNr].first;
+      Vector<int> boundaryFaces = violatingTets[violatingTetNr].second;
+
+      // Copying old tet information
+      Vector<int> currTetVerts;
+      Vector<int> currTetNeighbors;
+      for(int i = 0; i<4; i++){
+        currTetVerts.append(tetVertId[4*violatingTet + i]);
+        currTetNeighbors.append(tetFaceNeighbors[4*violatingTet + i]);
+      }
+
+      // Deleting old tet. For each deleted tet, index 1 stores the index of next free tet
+      tetVertId[4*violatingTet] = -1;
+      tetVertId[4*violatingTet + 1] = firstFreeTet;
+      firstFreeTet = violatingTet;
+
+      for(int i = 0; i<boundaryFaces.size(); i++){
+        Vector<int> faceVerts;
+        for(int j = 2; j>=0; j--){
+          faceVerts.append(currTetVerts[tetFaces[boundaryFaces[i]][j]]); 
+        }
+
+        // Make new tet
+        int newTetNr = -1;
+        if(firstFreeTet == -1){
+          newTetNr = tetVertId.size()/4;
+          for(int j = 0; j<4; j++){
+            tetVertId.append(j<3 ? faceVerts[j] : vertNr);
+            tetFaceNeighbors.append(-1);
+            faceNormals.append({0.0,0.0,0.0});
+            planesD.append(0.0);
+          }
+
+          tetMarks.append(0);
+        }
+        else{
+          newTetNr = firstFreeTet;
+          firstFreeTet = tetVertId[4*firstFreeTet + 1];
+          for(int j = 0; j<3; j++){
+            tetVertId[4*newTetNr + j] = faceVerts[j];
+            tetFaceNeighbors[4*newTetNr + j] = -1;
+          }
+          tetVertId[4*newTetNr + 3] = vertNr;
+          tetFaceNeighbors[4*newTetNr + 3] = -1;
+          tetMarks[newTetNr] = 0;
+        }
+        newTets.append(newTetNr);
+
+        tetFaceNeighbors[4*newTetNr] = currTetNeighbors[boundaryFaces[i]];
+        // If the boundary face has no neighboring tet
+        if(currTetNeighbors[boundaryFaces[i]] != -1){
+          // Else correcting the neighbors for the shared face
+          // tetFaceNeighbors[4*newTetNr] = currTetNeighbors[boundaryFaces[i]];
+          for(int j = 0; j<4; j++){
+            if(tetFaceNeighbors[4*currTetNeighbors[boundaryFaces[i]] + j] == violatingTet){
+              tetFaceNeighbors[4*currTetNeighbors[boundaryFaces[i]] + j] = newTetNr;
+              break;
+            }
+          }
+        }
+
+        setTetProperties(verts, tetVertId, faceNormals, planesD, newTetNr);
+      }
+    }
+
+    // Setting the neighbors of internal faces of new tets
+    for(int i = 0; i<newTets.size(); i++){
+      for(int j = 0; j<newTets.size(); j++){
+        
+        if(j == i){
+          continue;
+        }
+
+        for(int facei = 0; facei<4; facei++){
+          int vertsI[3];
+          for(int k = 0; k<3; k++){
+            vertsI[k] = tetVertId[4*newTets[i] + tetFaces[facei][k]];
+          }
+
+          int count = 0;
+          for(int k = 0; k<3; k++){
+            for(int l = 0; l<4; l++){
+              if(vertsI[k] == tetVertId[4*newTets[j] + l]){
+                count++;
+                break;
+              }
+            }
+          }
+
+          if(count == 3){
+            tetFaceNeighbors[4*newTets[i] + facei] = newTets[j];
+            // tetFaceNeighbors[4*newTets[i] + facei] = newTets[j];
+          }
+        }
+      }
+    }
+  }
+
+  // Remove empty tets and tets that lie outside the structure
+  int emptyTet = 0;
+  int tetLen = tetVertId.size()/4;
+  for(int tetNr = 0; tetNr<tetLen; tetNr++){
+    
+    int flag = 1;
+    float3 center(0.0f, 0.0f, 0.0f);
+    for(int i = 0; i<4; i++){
+      center += verts[tetVertId[4*tetNr + i]];
+      if(tetVertId[4*tetNr + i]<0 || tetVertId[4*tetNr + i]>=bigTet){
+        flag = 0;
+        break;
+      }
+    }
+    center*=0.25;
+
+    if(!flag || !isInside(center, treedata)){
+      continue;
+    }
+
+    for(int i = 0; i<4; i++){
+      tetVertId[4*emptyTet + i] = tetVertId[4*tetNr + i];
+    }
+    emptyTet++;
+  }
+
+    tetVertId.remove(4*emptyTet, 4*(tetLen - emptyTet));
+
+    return tetVertId;
+}
+
+static Mesh *modifyMesh(ModifierData *md, const ModifierEvalContext *UNUSED(ctx), Mesh *mesh)
+{
+  // Parameters of tetrahedralization, need to be taken as user input, being defined here as placeholders
+  float interiorResolution = 0;
+  float minTetQuality = 0.001; // Exp goes from -4 to 0
+  bool oneFacePerTet = false;
+  float tetScale = 0.8;
+
+  BVHTreeFromMesh treedata = {NULL};
+  BKE_bvhtree_from_mesh_get(&treedata, mesh, BVHTREE_FROM_LOOPTRI, 2);
+
+  Mesh *result;
+  BMesh *bm;
+  bool use_operators = true;
+  
+  BMeshCreateParams bmcParam;
+  bmcParam.use_toolflags = true;
+  bm = BM_mesh_create(&bm_mesh_allocsize_default,
+                      &bmcParam);
+  
+  Vector<float3> tetVerts;
+
+  float3 center(0,0,0);
+  float3 bmin(inf, inf, inf);
+  float3 bmax(-inf, -inf, -inf);
+
+  // Copying surface nodes to tetVerts
+  for(int i = 0; i<mesh->totvert; i++){
+    tetVerts.append({0.0, 0.0, 0.0});
+    tetVerts[i][0] = mesh->mvert[i].co[0] + randomEps();
+    tetVerts[i][1] = mesh->mvert[i].co[1] + randomEps();
+    tetVerts[i][2] = mesh->mvert[i].co[2] + randomEps();
+
+    center += tetVerts[i]; // Can cause overflow?
+    for(int axis = 0; axis<3; axis++){
+      bmin[axis] = min(bmin[axis], tetVerts[i][axis]);
+      bmax[axis] = max(bmax[axis], tetVerts[i][axis]);
+    }
+    
+  } 
+  center /= (float)tetVerts.size();
+
+  // Computing radius of bounding sphere (max dist of node from center)
+  float radius = 0.0;
+  for(int i = 0; i<tetVerts.size(); i++){
+    float dist = length(tetVerts[i] - center);
+    radius = max(dist, radius);
+  }
+  
+  // Interior sampling
+  if(interiorResolution > 0.0){
+    float boundLen[3];
+    sub_v3_v3v3(boundLen, bmax, bmin);
+    float maxBoundLen = max_axis_v3(boundLen);
+    float sampleLen = maxBoundLen/interiorResolution;
+
+    for(int xi = 0; xi<(int)(boundLen[0]/sampleLen); xi++){
+      float x = bmin[0] + xi*sampleLen + randomEps();
+      for(int yi = 0; yi<(int)(boundLen[1]/sampleLen); yi++){
+        float y = bmin[1] + yi*sampleLen + randomEps();
+        for(int zi = 0; zi<(int)(boundLen[2]/sampleLen); zi++){
+          float z = bmin[2] + zi*sampleLen + randomEps(); 
+
+          if(isInside({x,y,z}, &treedata)){
+            tetVerts.append({x,y,z});  
+          }
+        }
+      }
+    }
+  }
+
+  float bigTetSize = radius*5.0;
+
+  tetVerts.append({-bigTetSize, 0.0, -bigTetSize});
+  tetVerts.append({bigTetSize, 0.0, -bigTetSize});
+  tetVerts.append({0.0, bigTetSize, bigTetSize});
+  tetVerts.append({0.0, -bigTetSize, bigTetSize});
+
+  Vector<int> tetVertId = createTets(tetVerts, &treedata, minTetQuality);
+
+  Vector<BMVert *> bmverts;
+  if(oneFacePerTet){
+    for(int i = 0; i<mesh->totvert; i++){
+      bmverts.append(BM_vert_create(bm, mesh->mvert[i].co, NULL, BM_CREATE_NOP));
+    }
+    for(int i = mesh->totvert; i<tetVerts.size(); i++){
+      bmverts.append(BM_vert_create(bm, tetVerts[i], NULL, BM_CREATE_NOP));
+    }
+  }
+  else{
+    // Add vertices multiple times to make the tets distinctly visible
+    for(int tetNr = 0; tetNr<tetVertId.size()/4; tetNr++){
+      float3 center(0,0,0);
+      for(int j = 0; j<4; j++){
+        center += tetVerts[tetVertId[4*tetNr + j]];
+      }
+      center *= 0.25;
+
+      for(int faceNr = 0; faceNr < 4; faceNr++){
+        for(int faceVertNr = 0; faceVertNr < 3; faceVertNr++){
+          float3 vert = tetVerts[tetVertId[4*tetNr + tetFaces[faceNr][faceVertNr]]];
+          vert = center + (vert - center)*tetScale;
+          bmverts.append(BM_vert_create(bm, vert, NULL, BM_CREATE_NOP));
+        }
+      }
+    }
+  }
+
+  int numTets = tetVertId.size()/4;
+  int nr = 0;
+  for(int i = 0; i<numTets; i++){
+    if(oneFacePerTet){
+      if(tetVertId[4*i] < 0)
+        continue;
+      BM_face_create_quad_tri(bm, bmverts[tetVertId[4*i + 0]], bmverts[tetVertId[4*i + 1]], bmverts[tetVertId[4*i + 2]], bmverts[tetVertId[4*i + 3]], NULL, BM_CREATE_NO_DOUBLE);
+    }
+    else{
+      for(int faceNr = 0; faceNr<4; faceNr++){
+        BM_face_create_quad_tri(bm, bmverts[nr], bmverts[nr+1], bmverts[nr+2], NULL, NULL, BM_CREATE_NO_DOUBLE);
+        nr+=3;
+      }
+    }
+  }
+
+  CustomData_MeshMasks cd_mask_extra = {
+    .vmask = CD_MASK_ORIGINDEX, .emask = CD_MASK_ORIGINDEX, .pmask = CD_MASK_ORIGINDEX};
+  
+  result = BKE_mesh_from_bmesh_for_eval_nomain(bm, &cd_mask_extra, mesh);
+  BM_mesh_free(bm);
+
+  return result;    
+}
+
+// --------------------------------------
+// --------------------------------------
+
+// static Span<MDeformVert> get_vertex_group(const Mesh &mesh, const int defgrp_index)
+// {
+//   if (defgrp_index == -1) {
+//     return {};
+//   }
+//   const MDeformVert *vertex_group = static_cast<const MDeformVert *>(
+//       CustomData_get_layer(&mesh.vdata, CD_MDEFORMVERT));
+//   if (!vertex_group) {
+//     return {};
+//   }
+//   return {vertex_group, mesh.totvert};
+// }
+
+// static Vector<int64_t> selected_indices_from_vertex_group(Span<MDeformVert> vertex_group,
+//                                                           const int index,
+//                                                           const bool invert)
+// {
+//   Vector<int64_t> selected_indices;
+//   for (const int i : vertex_group.index_range()) {
+//     const bool found = BKE_defvert_find_weight(&vertex_group[i], index) > 0.0f;
+//     if (found != invert) {
+//       selected_indices.append(i);
+//     }
+//   }
+//   return selected_indices;
+// }
+
+// static Array<bool> selection_array_from_vertex_group(Span<MDeformVert> vertex_group,
+//                                                      const int index,
+//                                                      const bool invert)
+// {
+//   Array<bool> selection(vertex_group.size());
+//   for (const int i : vertex_group.index_range()) {
+//     const bool found = BKE_defvert_find_weight(&vertex_group[i], index) > 0.0f;
+//     selection[i] = (found != invert);
+//   }
+//   return selection;
+// }
+
+// static std::optional<Mesh *> calculate_weld(const Mesh &mesh, const WeldModifierData &wmd)
+// {
+//   const int defgrp_index = BKE_id_defgroup_name_index(&mesh.id, wmd.defgrp_name);
+//   Span<MDeformVert> vertex_group = get_vertex_group(mesh, defgrp_index);
+//   const bool invert = (wmd.flag & MOD_WELD_INVERT_VGROUP) != 0;
+
+//   if (wmd.mode == MOD_WELD_MODE_ALL) {
+//     if (!vertex_group.is_empty()) {
+//       Vector<int64_t> selected_indices = selected_indices_from_vertex_group(
+//           vertex_group, defgrp_index, invert);
+//       return blender::geometry::mesh_merge_by_distance_all(
+//           mesh, IndexMask(selected_indices), wmd.merge_dist);
+//     }
+//     return blender::geometry::mesh_merge_by_distance_all(
+//         mesh, IndexMask(mesh.totvert), wmd.merge_dist);
+//   }
+//   if (wmd.mode == MOD_WELD_MODE_CONNECTED) {
+//     const bool only_loose_edges = (wmd.flag & MOD_WELD_LOOSE_EDGES) != 0;
+//     if (!vertex_group.is_empty()) {
+//       Array<bool> selection = selection_array_from_vertex_group(
+//           vertex_group, defgrp_index, invert);
+//       return blender::geometry::mesh_merge_by_distance_connected(
+//           mesh, selection, wmd.merge_dist, only_loose_edges);
+//     }
+//     Array<bool> selection(mesh.totvert, true);
+//     return blender::geometry::mesh_merge_by_distance_connected(
+//         mesh, selection, wmd.merge_dist, only_loose_edges);
+//   }
+
+//   BLI_assert_unreachable();
+//   return nullptr;
+// }
+
+// static Mesh *modifyMesh(ModifierData *md, const ModifierEvalContext *UNUSED(ctx), Mesh *mesh)
+// {
+// const WeldModifierData &wmd = reinterpret_cast<WeldModifierData &>(*md);
+
+//   std::optional<Mesh *> result = calculate_weld(*mesh, wmd);
+//   if (!result) {
+//     return mesh;
+//   }
+//   return *result;
+// }
+
+static void initData(ModifierData *md)
+{
+  WeldModifierData *wmd = (WeldModifierData *)md;
+
+  BLI_assert(MEMCMP_STRUCT_AFTER_IS_ZERO(wmd, modifier));
+
+  MEMCPY_STRUCT_AFTER(wmd, DNA_struct_default_get(WeldModifierData), modifier);
+}
+
+static void requiredDataMask(Object *UNUSED(ob),
+                             ModifierData *md,
+                             CustomData_MeshMasks *r_cddata_masks)
+{
+  WeldModifierData *wmd = (WeldModifierData *)md;
+
+  /* Ask for vertexgroups if we need them. */
+  if (wmd->defgrp_name[0] != '\0') {
+    r_cddata_masks->vmask |= CD_MASK_MDEFORMVERT;
+  }
+}
+
+static void panel_draw(const bContext *UNUSED(C), Panel *panel)
+{
+  uiLayout *layout = panel->layout;
+
+  PointerRNA ob_ptr;
+  PointerRNA *ptr = modifier_panel_get_property_pointers(panel, &ob_ptr);
+  int weld_mode = RNA_enum_get(ptr, "mode");
+
+  uiLayoutSetPropSep(layout, true);
+
+  uiItemR(layout, ptr, "mode", 0, nullptr, ICON_NONE);
+  uiItemR(layout, ptr, "merge_threshold", 0, IFACE_("Distance"), ICON_NONE);
+  if (weld_mode == MOD_WELD_MODE_CONNECTED) {
+    uiItemR(layout, ptr, "loose_edges", 0, nullptr, ICON_NONE);
+  }
+  modifier_vgroup_ui(layout, ptr, &ob_ptr, "vertex_group", "invert_vertex_group", nullptr);
+
+  modifier_panel_end(layout, ptr);
+}
+
+static void panelRegister(ARegionType *region_type)
+{
+  modifier_panel_register(region_type, eModifierType_Weld, panel_draw);
+}
+
+ModifierTypeInfo modifierType_Weld = {
+    /* name */ "Weld",
+    /* structName */ "WeldModifierData",
+    /* structSize */ sizeof(WeldModifierData),
+    /* srna */ &RNA_WeldModifier,
+    /* type */ eModifierTypeType_Constructive,
+    /* flags */
+    (ModifierTypeFlag)(eModifierTypeFlag_AcceptsMesh | eModifierTypeFlag_SupportsMapping |
+                       eModifierTypeFlag_SupportsEditmode | eModifierTypeFlag_EnableInEditmode |
+                       eModifierTypeFlag_AcceptsCVs),
+    /* icon */ ICON_AUTOMERGE_OFF, /* TODO: Use correct icon. */
+
+    /* copyData */ BKE_modifier_copydata_generic,
+
+    /* deformVerts */ nullptr,
+    /* deformMatrices */ nullptr,
+    /* deformVertsEM */ nullptr,
+    /* deformMatricesEM */ nullptr,
+    /* modifyMesh */ modifyMesh,
+    /* modifyGeometrySet */ nullptr,
+
+    /* initData */ initData,
+    /* requiredDataMask */ requiredDataMask,
+    /* freeData */ nullptr,
+    /* isDisabled */ nullptr,
+    /* updateDepsgraph */ nullptr,
+    /* dependsOnTime */ nullptr,
+    /* dependsOnNormals */ nullptr,
+    /* foreachIDLink */ nullptr,
+    /* foreachTexLink */ nullptr,
+    /* freeRuntimeData */ nullptr,
+    /* panelRegister */ panelRegister,
+    /* blendWrite */ nullptr,
+    /* blendRead */ nullptr,
+};

source/blender/modifiers/intern/MOD_weld.cc

@@ -18,7 +18,16 @@
 #include "BLI_array.hh"
 #include "BLI_index_range.hh"
 #include "BLI_span.hh"
+
+#include "BLI_math_vector.h"
+#include "BLI_utildefines.h"
+
 #include "BLI_vector.hh"
+#include "BLI_math_vec_types.hh"
+#include "BLI_math_vector.hh"
+
+#include "bmesh.h"
+#include "bmesh_tools.h"
 
 #include "BLT_translation.h"
 
@@ -28,14 +37,12 @@
 #include "DNA_modifier_types.h"
 #include "DNA_screen_types.h"
 
-#ifdef USE_BVHTREEKDOP
-#  include "BKE_bvhutils.h"
-#endif
-
+#include "BKE_bvhutils.h"
 #include "BKE_context.h"
 #include "BKE_deform.h"
 #include "BKE_modifier.h"
 #include "BKE_screen.h"
+#include "BKE_mesh.h"
 
 #include "UI_interface.h"
 #include "UI_resources.h"
@@ -54,90 +61,720 @@ using blender::Array;
 using blender::IndexMask;
 using blender::Span;
 using blender::Vector;
+using namespace blender;
+using namespace blender::math;
 
-static Span<MDeformVert> get_vertex_group(const Mesh &mesh, const int defgrp_index)
-{
-  if (defgrp_index == -1) {
-    return {};
-  }
-  const MDeformVert *vertex_group = static_cast<const MDeformVert *>(
-      CustomData_get_layer(&mesh.vdata, CD_MDEFORMVERT));
-  if (!vertex_group) {
-    return {};
-  }
-  return {vertex_group, mesh.totvert};
+int tetFaces[4][3] = {{2,1,0}, {0,1,3}, {1,2,3}, {2,0,3}};
+float directions[6][3] = {{1,0,0}, {-1,0,0}, {0,1,0}, {0,-1,0}, {0,0,1}, {0,0,-1}};
+
+float inf = FLT_MAX;
+float eps = 0.0f;
+
+bool globalFlag = false;
+
+static float randomEps(){
+  float eps = 0.0001 ;
+  return -eps + 2.0 * (static_cast <float> (rand()) / static_cast <float> (RAND_MAX)) * eps;
 }
 
-static Vector<int64_t> selected_indices_from_vertex_group(Span<MDeformVert> vertex_group,
-                                                          const int index,
-                                                          const bool invert)
-{
-  Vector<int64_t> selected_indices;
-  for (const int i : vertex_group.index_range()) {
-    const bool found = BKE_defvert_find_weight(&vertex_group[i], index) > 0.0f;
-    if (found != invert) {
-      selected_indices.append(i);
-    }
-  }
-  return selected_indices;
+static float tetQuality(float3 p0, float3 p1, float3 p2, float3 p3){
+  float3 d0 = p1 - p0;
+  float3 d1 = p2 - p0;
+  float3 d2 = p3 - p0;
+  float3 d3 = p2 - p1;
+  float3 d4 = p3 - p2;
+  float3 d5 = p1 - p3;
+
+  float s0 = length(d0);
+  float s1 = length(d1);
+  float s2 = length(d2);
+  float s3 = length(d3);
+  float s4 = length(d4);
+  float s5 = length(d5);
+
+  float ms = (s0*s0 + s1*s1 + s2*s2 + s3*s3 + s4*s4 + s5*s5) / 6.0;
+  float rms = sqrt(ms);
+
+  float s = 12.0 / sqrt(2.0);
+
+  float vol = dot(d0, cross(d1, d2)) / 6.0;
+  return s * vol / (rms * rms * rms);
 }
 
-static Array<bool> selection_array_from_vertex_group(Span<MDeformVert> vertex_group,
-                                                     const int index,
-                                                     const bool invert)
-{
-  Array<bool> selection(vertex_group.size());
-  for (const int i : vertex_group.index_range()) {
-    const bool found = BKE_defvert_find_weight(&vertex_group[i], index) > 0.0f;
-    selection[i] = (found != invert);
+static bool isInside(float3 vert, BVHTreeFromMesh *treedata){
+
+  int count =  0;
+  float min_dist = 0.0f;
+  for(auto dir : directions){
+    float radius = 0.0f;
+    
+    float max_length = FLT_MAX;
+
+    BVHTreeRayHit rayhit = {0};
+    rayhit.index = -1;
+    rayhit.dist = max_length;
+    BLI_bvhtree_ray_cast(treedata->tree, vert, dir, radius, &rayhit, treedata->raycast_callback, treedata);
+
+    if (rayhit.index != -1 && rayhit.dist <= max_length) {
+      if(dot_v3v3(rayhit.no, dir) > 0.0f){
+        count++;
+      }
+
+      if((min_dist > 0.0) && (min_dist - rayhit.dist) > 0.0)
+        return false;
+    }
   }
-  return selection;
+
+  return count > 3;
 }
 
-static std::optional<Mesh *> calculate_weld(const Mesh &mesh, const WeldModifierData &wmd)
-{
-  const int defgrp_index = BKE_id_defgroup_name_index(&mesh.id, wmd.defgrp_name);
-  Span<MDeformVert> vertex_group = get_vertex_group(mesh, defgrp_index);
-  const bool invert = (wmd.flag & MOD_WELD_INVERT_VGROUP) != 0;
+// Used to populate a tets face normals and planesD
+static void setTetProperties(Vector<float3> &verts, 
+                      Vector<int> &tetVertId,
+                      Vector<float3> &faceNormals, 
+                      Vector<float> &planesD,
+                      int tetNr){
+  for(int i = 0; i<4; i++){ 
+    float3 p0 = verts[tetVertId[4*tetNr + tetFaces[i][0]]];
+    float3 p1 = verts[tetVertId[4*tetNr + tetFaces[i][1]]];
+    float3 p2 = verts[tetVertId[4*tetNr + tetFaces[i][2]]];
 
-  if (wmd.mode == MOD_WELD_MODE_ALL) {
-    if (!vertex_group.is_empty()) {
-      Vector<int64_t> selected_indices = selected_indices_from_vertex_group(
-          vertex_group, defgrp_index, invert);
-      return blender::geometry::mesh_merge_by_distance_all(
-          mesh, IndexMask(selected_indices), wmd.merge_dist);
-    }
-    return blender::geometry::mesh_merge_by_distance_all(
-        mesh, IndexMask(mesh.totvert), wmd.merge_dist);
+    float3 normal = cross(p1 - p0, p2 - p0);
+    normal = normalize(normal);
+
+    faceNormals[4*tetNr + i] = normal;
+    planesD[4*tetNr + i] = dot(p0, normal);
   }
-  if (wmd.mode == MOD_WELD_MODE_CONNECTED) {
-    const bool only_loose_edges = (wmd.flag & MOD_WELD_LOOSE_EDGES) != 0;
-    if (!vertex_group.is_empty()) {
-      Array<bool> selection = selection_array_from_vertex_group(
-          vertex_group, defgrp_index, invert);
-      return blender::geometry::mesh_merge_by_distance_connected(
-          mesh, selection, wmd.merge_dist, only_loose_edges);
+}
+
+static bool edgeCompare(const Vector<int> &e0, const Vector<int> &e1){
+  if((e0[0] < e1[0]) || ((e0[0] == e1[0]) && (e0[1] < e1[1])))
+    return true;
+  else
+    return false;
+}
+
+static float3 getCircumCenter(float3 p0, float3 p1, float3 p2, float3 p3){
+  float eps = 0.000001f;
+
+  float3 b = p1 - p0;
+  float3 c = p2 - p0;
+  float3 d = p3 - p0;
+
+  float det = 2.0 * (b.x*(c.y*d.z - c.z*d.y) - b.y*(c.x*d.z - c.z*d.x) + b.z*(c.x*d.y - c.y*d.x));
+  if (det <= eps && det >= -eps){
+    return p0;
+  }
+  else{
+    float3 v = cross(c, d)*dot(b, b) + cross(d, b)*dot(c, c) + cross(b, c)*dot(d, d);
+    v /= det;
+    return p0 + v;
+  }
+      
+}
+
+static int findContainingTet(Vector<float3> &verts, 
+                    Vector<int> &tetVertId, 
+                    Vector<int> &tetFaceNeighbors, 
+                    Vector<float3> &faceNormals, 
+                    Vector<float> &planesD, 
+                    int tetMarkId, 
+                    Vector<int> &tetMarks,
+                    float3 currVert){
+  
+  /* --------------------
+  Matthias Method
+  ----------------------*/
+
+  // bool found = false;
+  // int tetNr = 0;
+
+  // while(tetNr < tetVertId.size()/4 && tetVertId[4*tetNr]<0)
+  //   tetNr++;
+
+  // float3 center(0.0, 0.0, 0.0);
+  // while(!found){
+  //   if(tetNr < 0 || tetMarks[tetNr] == tetMarkId){
+  //     break;
+  //   }
+    
+  //   tetMarks[tetNr] = tetMarkId;
+
+  //   center = {0.0,0.0,0.0};
+  //   for(int i = 0; i<4; i++){
+  //     center += verts[tetVertId[4*tetNr + i]];
+  //   }
+  //   center*=0.25;
+
+  //   float minT = inf; //
+  //   int minFaceNr = -1;
+
+  //   for(int i = 0; i<4; i++){
+  //     float3 normal = faceNormals[4*tetNr + i];
+  //     float d = planesD[4*tetNr + i];
+
+  //     float hp = dot(normal, currVert) - d;
+  //     float hc = dot(normal, center) - d;
+
+  //     float t = hp - hc;
+  //     if(t <= eps)
+  //       continue;
+
+  //     t = -hc/t;
+
+  //     if((t >= eps) && ((t - minT) > eps)){
+  //       minT = t;
+  //       minFaceNr = i;
+  //     }
+  //   }
+
+  //   if(minT >= 1.0){
+  //     found = true;
+  //   }
+  //   else{
+  //     tetNr = tetFaceNeighbors[4*tetNr + minFaceNr];
+  //   }
+  // }
+
+  // if(found)
+  //   return tetNr;
+  // else
+  //   return -1;
+
+  /* --------------------
+  Brute force finding first violating tet
+  ----------------------*/
+
+  for(int currTet = 0; currTet<tetVertId.size()/4; currTet++){
+    if(tetVertId[4*currTet] == -1)
+      continue;
+    
+    float3 p0 = verts[tetVertId[4*currTet + 0]];
+    float3 p1 = verts[tetVertId[4*currTet + 1]];
+    float3 p2 = verts[tetVertId[4*currTet + 2]];
+    float3 p3 = verts[tetVertId[4*currTet + 3]];
+
+    float3 circumCenter = getCircumCenter(p0, p1, p2, p3);
+    float circumRadius = length(p0 - circumCenter);
+    
+    // if((circumRadius - length(currVert - circumCenter)) >= eps){
+    if(length(currVert - circumCenter) < circumRadius){
+      return currTet;
     }
-    Array<bool> selection(mesh.totvert, true);
-    return blender::geometry::mesh_merge_by_distance_connected(
-        mesh, selection, wmd.merge_dist, only_loose_edges);
   }
 
-  BLI_assert_unreachable();
-  return nullptr;
+  /* -----------------------
+  My method, checks if point is inside tet, if not finds the face that is connecting to the point and center
+  ------------------------*/
+
+  // bool found = false;
+  // int tetNr = 0;
+  // while(tetNr < tetVertId.size()/4 && tetVertId[4*tetNr] < 0)
+  //   tetNr++;
+
+  // while(!found){
+  //   float3 p0 = verts[tetVertId[4*currTet + 0]];
+  //   float3 p1 = verts[tetVertId[4*currTet + 1]];
+  //   float3 p2 = verts[tetVertId[4*currTet + 2]];
+  //   float3 p3 = verts[tetVertId[4*currTet + 3]];
+    
+  //   if(isInsideTet(p0, p1, p2, p3, currVert))
+  //     return tetNr;
+    
+  //   for(int i = 0; i<4; i++)
+  // }
+
+  return -1;
+}
+
+/*
+The basic assumption employed is that 2 violating tets cannot not be neighbors. 
+Simple BFS approach that checks neighbors of all violating tets and adds them to stack
+
+If a violating tet shares a face with a non-violating tet, that's a boundary face and for each violating tet, list of its boundary faces is returned
+*/
+static Vector<int> getViolatingTets(Vector<float3> &verts, 
+                            Vector<int> &tetVertId, 
+                            Vector<int> &tetFaceNeighbors, 
+                            int tetMarkId, 
+                            Vector<int> &tetMarks,
+                            float3 currVert,
+                            int containingTetNr
+                            ){
+  
+  Vector<int> violatingTets;
+  Vector<int> stack;
+
+  stack.append(containingTetNr);  
+  tetMarks[containingTetNr] = tetMarkId; 
+
+  while(stack.size()){
+    int currTet = stack.last();
+    stack.remove_last();
+    violatingTets.append(currTet);
+    
+    for(int i = 0; i<4; i++){
+      int neighborTet = tetFaceNeighbors[4*currTet + i];
+      
+      if(neighborTet<0 || tetMarks[neighborTet]==tetMarkId){
+        continue;
+      }
+      
+      if(tetVertId[4*neighborTet] < 0){
+        globalFlag = true;
+        return {};
+      }
+
+      float3 p0 = verts[tetVertId[4*neighborTet + 0]];
+      float3 p1 = verts[tetVertId[4*neighborTet + 1]];
+      float3 p2 = verts[tetVertId[4*neighborTet + 2]];
+      float3 p3 = verts[tetVertId[4*neighborTet + 3]];
+
+      float3 circumCenter = getCircumCenter(p0, p1, p2, p3);
+      float circumRadius = length(p0 - circumCenter);
+
+      // if((circumRadius - length(currVert - circumCenter)) >= eps){
+      if(length(currVert - circumCenter) < circumRadius){
+        stack.append(neighborTet);
+        tetMarks[neighborTet] = tetMarkId;
+      }
+    }
+  }
+
+  // for(int i = 0; i<violatingTets.size(); i++){
+  //   for(int j = i+1; j<violatingTets.size(); j++)
+  //     if(violatingTets[i] == violatingTets[j])
+  //       std::cout << "Duplicates found in violating tets" << std::endl;
+  // }
+
+  return violatingTets;
+}
+
+static Vector<int> createTets(Vector<float3> verts, BVHTreeFromMesh *treedata, float minTetQuality){
+  Vector<int> tetVertId; // Stores indices of vertex that form a tet. Every tet is stored as 4 indices 
+  Vector<int> tetFaceNeighbors; // Stores index of tet that shares the face as per tetFaces order. 4 neighbors per tet, 1 for each face
+
+  int firstFreeTet = -1;
+
+  Vector<float3> faceNormals; // Stores the normal of each of face of the tet
+  Vector<float> planesD;
+
+  int tetMarkId = 0;
+  Vector<int> tetMarks; 
+  /*Used to keep track of visited tets for various processess. In each step where visited tets need to avoided, 
+  on visiting them, a unique number is stored representing that particular step
+  */
+ 
+  int bigTet = verts.size() - 4;  
+
+  for(int i = 0; i<4; i++){
+    tetVertId.append(bigTet + i);
+    tetFaceNeighbors.append(-1);
+    faceNormals.append({0.0, 0.0, 0.0});
+    planesD.append(0.0);
+  }
+  tetMarks.append(0);
+  setTetProperties(verts, tetVertId, faceNormals, planesD, 0);
+
+  for(int vertNr = 0; vertNr<bigTet; vertNr++){
+    float3 currVert = verts[vertNr];
+
+    // std::cout << "Adding vert " << vertNr << std::endl;
+
+    // Find containing tet (need to understand) center can lie outside the tet?
+    tetMarkId += 1;
+    int containingTetNr = -1;
+    containingTetNr = findContainingTet(verts, tetVertId, tetFaceNeighbors, faceNormals, planesD, tetMarkId, tetMarks, currVert);
+    
+    if(containingTetNr == -1){
+      std::cout << "Couldn't add vert " << vertNr << std::endl;
+      continue;
+    }
+
+    // Find Violating Tets - Returns all tets that are violating the Delaunay condition along with a list of face indices that are boundary faces
+    // Boundary faces may be outermost face of the structure or the face shared by a violating tet and a non violating tet
+    tetMarkId += 1;
+    Vector<int> violatingTets = getViolatingTets(verts, tetVertId, tetFaceNeighbors, tetMarkId, tetMarks, currVert, containingTetNr);
+
+    //Create new tets centered at the new point and including the boundary faces
+    Vector<int> newTets; // Stores tet indices of the new tets formed. Used for making neighbors of the new tets formed
+    Vector<Vector<int>> edges;
+    for(int violatingTetNr = 0; violatingTetNr<violatingTets.size(); violatingTetNr++){
+      int violatingTet = violatingTets[violatingTetNr];
+
+      // Copying old tet information
+      Vector<int> currTetVerts;
+      Vector<int> currTetNeighbors;
+      for(int i = 0; i<4; i++){
+        currTetVerts.append(tetVertId[4*violatingTet + i]);
+        currTetNeighbors.append(tetFaceNeighbors[4*violatingTet + i]);
+      }
+
+      // Deleting old tet. For each deleted tet, index 1 stores the index of next free tet
+      tetVertId[4*violatingTet] = -1;
+      tetVertId[4*violatingTet + 1] = firstFreeTet;
+      firstFreeTet = violatingTet;
+
+      for(int i = 0; i<4; i++){
+        if(currTetNeighbors[i] >= 0 && tetMarks[currTetNeighbors[i]] == tetMarkId){
+          continue;
+        }
+
+        // Make new tet
+        int newTetNr = firstFreeTet;
+        if(firstFreeTet == -1){
+          newTetNr = tetVertId.size()/4;
+          for(int j = 0; j<4; j++){
+            tetVertId.append(-1);
+            tetFaceNeighbors.append(-1);
+            faceNormals.append({0.0,0.0,0.0});
+            planesD.append(0.0);
+          }
+          tetMarks.append(0);
+        }
+        else{
+          firstFreeTet = tetVertId[4*firstFreeTet + 1];
+        }
+
+        int id0 = currTetVerts[tetFaces[i][2]];
+        int id1 = currTetVerts[tetFaces[i][1]];
+        int id2 = currTetVerts[tetFaces[i][0]];
+
+        tetVertId[4 * newTetNr] = id0;
+        tetVertId[4 * newTetNr + 1] = id1;
+        tetVertId[4 * newTetNr + 2] = id2;
+        tetVertId[4 * newTetNr + 3] = vertNr;
+
+        tetFaceNeighbors[4*newTetNr] = currTetNeighbors[i];
+        
+        if(currTetNeighbors[i] >= 0){
+          // Else correcting the neighbors for the shared face
+          for(int j = 0; j<4; j++){
+            if(tetFaceNeighbors[4*currTetNeighbors[i] + j] == violatingTet){
+              tetFaceNeighbors[4*currTetNeighbors[i] + j] = newTetNr;
+              // break;
+            }
+          }
+        }
+
+        for(int j = 1; j<4; j++){
+          tetFaceNeighbors[4 * newTetNr + j] = -1;
+        }
+
+        setTetProperties(verts, tetVertId, faceNormals, planesD, newTetNr);
+
+        edges.append({min(id0, id1), max(id0, id1), newTetNr, 1});
+        edges.append({min(id1, id2), max(id1, id2), newTetNr, 2});
+        edges.append({min(id2, id0), max(id2, id0), newTetNr, 3});
+      }
+    }
+
+    std::sort(edges.begin(), edges.end(), edgeCompare);
+    int nr = 0;
+    int numEdges = edges.size();
+
+    while(nr < numEdges){
+      Vector<int> e0 = edges[nr];
+      nr += 1;
+
+      if((nr < numEdges) && (edges[nr][0] == e0[0]) && (edges[nr][1] == e0[1])){
+        Vector<int> e1 = edges[nr];
+
+        tetFaceNeighbors[4*e0[2] + e0[3]] = e1[2];
+        tetFaceNeighbors[4*e1[2] + e1[3]] = e0[2];
+
+        nr += 1;
+      }
+    }
+  }
+
+  // Remove empty tets and tets that lie outside the structure
+  int emptyTet = 0;
+  int tetLen = tetVertId.size()/4;
+  for(int tetNr = 0; tetNr<tetLen; tetNr++){
+    
+    int flag = 1;
+    float3 center(0.0f, 0.0f, 0.0f);
+    for(int i = 0; i<4; i++){
+      if(tetVertId[4*tetNr + i]<0 || tetVertId[4*tetNr + i]>=bigTet){
+        flag = 0;
+        break;
+      }
+      center += verts[tetVertId[4*tetNr + i]];
+    }
+    center*=0.25;
+
+    if(!flag || !isInside(center, treedata)){
+      continue;
+    }
+
+    float3 p0 = verts[tetVertId[4 * tetNr + 0]];
+    float3 p1 = verts[tetVertId[4 * tetNr + 1]];
+    float3 p2 = verts[tetVertId[4 * tetNr + 2]];
+    float3 p3 = verts[tetVertId[4 * tetNr + 3]];
+
+    if(tetQuality(p0, p1, p2, p3) < minTetQuality){
+      continue;
+    }
+
+    for(int i = 0; i<4; i++){
+      tetVertId[4*emptyTet + i] = tetVertId[4*tetNr + i];
+    }
+    emptyTet++;
+  }
+
+  tetVertId.remove(4*emptyTet, 4*(tetLen - emptyTet));
+  return tetVertId;
 }
 
 static Mesh *modifyMesh(ModifierData *md, const ModifierEvalContext *UNUSED(ctx), Mesh *mesh)
 {
   const WeldModifierData &wmd = reinterpret_cast<WeldModifierData &>(*md);
 
-  std::optional<Mesh *> result = calculate_weld(*mesh, wmd);
-  if (!result) {
-    return mesh;
+  globalFlag = false;
+  // Parameters of tetrahedralization, need to be taken as user input, being defined here as placeholders
+  float interiorResolution = wmd.merge_dist;
+  float minTetQuality = 0.01; // Exp goes from -4 to 0
+  bool oneFacePerTet = true;
+  float tetScale = 0.8;
+
+  BVHTreeFromMesh treedata = {NULL};
+  BKE_bvhtree_from_mesh_get(&treedata, mesh, BVHTREE_FROM_LOOPTRI, 2);
+
+  Mesh *result;
+  BMesh *bm;
+  bool use_operators = true;
+  
+  BMeshCreateParams bmcParam;
+  bmcParam.use_toolflags = true;
+  bm = BM_mesh_create(&bm_mesh_allocsize_default,
+                      &bmcParam);
+  
+  Vector<float3> tetVerts;
+
+  float3 center(0,0,0);
+  float3 bmin(inf, inf, inf);
+  float3 bmax(-inf, -inf, -inf);
+
+  // Copying surface nodes to tetVerts
+  for(int i = 0; i<mesh->totvert; i++){
+    tetVerts.append({0.0, 0.0, 0.0});
+    tetVerts[i][0] = mesh->mvert[i].co[0] + randomEps();
+    tetVerts[i][1] = mesh->mvert[i].co[1] + randomEps();
+    tetVerts[i][2] = mesh->mvert[i].co[2] + randomEps();
+
+    center += tetVerts[i]; // Can cause overflow?
+    for(int axis = 0; axis<3; axis++){
+      bmin[axis] = min(bmin[axis], tetVerts[i][axis]);
+      bmax[axis] = max(bmax[axis], tetVerts[i][axis]);
+    }
+    
+  } 
+  center /= (float)tetVerts.size();
+
+  // Computing radius of bounding sphere (max dist of node from center)
+  float radius = 0.0;
+  for(int i = 0; i<tetVerts.size(); i++){
+    float dist = length(tetVerts[i] - center);
+    radius = max(dist, radius);
   }
-  return *result;
+  
+  // Interior sampling
+  if(interiorResolution > 0.0){
+    float boundLen[3];
+    sub_v3_v3v3(boundLen, bmax, bmin);
+    float maxBoundLen = max_axis_v3(boundLen);
+    float sampleLen = maxBoundLen/interiorResolution;
+
+    for(int xi = 0; xi<(int)(boundLen[0]/sampleLen); xi++){
+      float x = bmin[0] + xi*sampleLen + randomEps();
+      for(int yi = 0; yi<(int)(boundLen[1]/sampleLen); yi++){
+        float y = bmin[1] + yi*sampleLen + randomEps();
+        for(int zi = 0; zi<(int)(boundLen[2]/sampleLen); zi++){
+          float z = bmin[2] + zi*sampleLen + randomEps(); 
+
+          if(isInside({x,y,z}, &treedata)){
+            tetVerts.append({x,y,z});  
+          }
+        }
+      }
+    }
+  }
+
+  float bigTetSize = radius*5.0;
+
+  tetVerts.append({-bigTetSize, 0.0, -bigTetSize});
+  tetVerts.append({bigTetSize, 0.0, -bigTetSize});
+  tetVerts.append({0.0, bigTetSize, bigTetSize});
+  tetVerts.append({0.0, -bigTetSize, bigTetSize});
+
+  Vector<int> tetVertId = createTets(tetVerts, &treedata, minTetQuality);
+
+  // Creates a 2 tets to test collisions
+  bool testing = true;
+
+  if(testing){
+    tetVertId.remove(0, tetVertId.size());
+    tetVerts.remove(0, tetVerts.size());
+
+    tetVerts.append({-1, -1, 0});
+    tetVerts.append({-1, 1, 0});
+    tetVerts.append({1, -1, 0});
+    tetVerts.append({-1, -1, 1});
+
+    tetVerts.append({0, 0, 1});
+    tetVerts.append({0, 0, 2});
+    tetVerts.append({0, 1, 2});
+    tetVerts.append({1, 0, 2});
+
+    for(int i = 0; i<8; i++)
+      tetVertId.append(i);
+  }
+
+  // if(globalFlag){
+  //   return mesh;
+  // }
+
+  Vector<BMVert *> bmverts;
+  if(oneFacePerTet){
+    for(int i = 0; i<mesh->totvert; i++){
+      // bmverts.append(BM_vert_create(bm, mesh->mvert[i].co, NULL, BM_CREATE_NOP));
+      bmverts.append(BM_vert_create(bm, tetVerts[i], NULL, BM_CREATE_NOP));
+    }
+    for(int i = mesh->totvert; i<tetVerts.size()-4; i++){
+      bmverts.append(BM_vert_create(bm, tetVerts[i], NULL, BM_CREATE_NOP));
+    }
+  }
+  else{
+    // Add vertices multiple times to make the tets distinctly visible
+    for(int tetNr = 0; tetNr<tetVertId.size()/4; tetNr++){
+      float3 center(0,0,0);
+      for(int j = 0; j<4; j++){
+        center += tetVerts[tetVertId[4*tetNr + j]];
+      }
+      center *= 0.25;
+
+      for(int faceNr = 0; faceNr < 4; faceNr++){
+        for(int faceVertNr = 0; faceVertNr < 3; faceVertNr++){
+          float3 vert = tetVerts[tetVertId[4*tetNr + tetFaces[faceNr][faceVertNr]]];
+          vert = center + (vert - center)*tetScale;
+          bmverts.append(BM_vert_create(bm, vert, NULL, BM_CREATE_NOP));
+        }
+      }
+    }
+  }
+
+  int numTets = tetVertId.size()/4;
+  int nr = 0;
+  for(int i = 0; i<numTets; i++){
+    if(oneFacePerTet){
+      if(tetVertId[4*i] < 0)
+        continue;
+      BM_face_create_quad_tri(bm, bmverts[tetVertId[4*i + 0]], bmverts[tetVertId[4*i + 1]], bmverts[tetVertId[4*i + 2]], bmverts[tetVertId[4*i + 3]], NULL, BM_CREATE_NO_DOUBLE);
+    }
+    else{
+      for(int faceNr = 0; faceNr<4; faceNr++){
+        BM_face_create_quad_tri(bm, bmverts[nr], bmverts[nr+1], bmverts[nr+2], NULL, NULL, BM_CREATE_NO_DOUBLE);
+        nr+=3;
+      }
+    }
+  }
+
+  CustomData_MeshMasks cd_mask_extra = {
+    .vmask = CD_MASK_ORIGINDEX, .emask = CD_MASK_ORIGINDEX, .pmask = CD_MASK_ORIGINDEX};
+  
+  result = BKE_mesh_from_bmesh_for_eval_nomain(bm, &cd_mask_extra, mesh);
+  BM_mesh_free(bm);
+
+  return result;    
 }
 
+// --------------------------------------
+// --------------------------------------
+
+// static Span<MDeformVert> get_vertex_group(const Mesh &mesh, const int defgrp_index)
+// {
+//   if (defgrp_index == -1) {
+//     return {};
+//   }
+//   const MDeformVert *vertex_group = static_cast<const MDeformVert *>(
+//       CustomData_get_layer(&mesh.vdata, CD_MDEFORMVERT));
+//   if (!vertex_group) {
+//     return {};
+//   }
+//   return {vertex_group, mesh.totvert};
+// }
+
+// static Vector<int64_t> selected_indices_from_vertex_group(Span<MDeformVert> vertex_group,
+//                                                           const int index,
+//                                                           const bool invert)
+// {
+//   Vector<int64_t> selected_indices;
+//   for (const int i : vertex_group.index_range()) {
+//     const bool found = BKE_defvert_find_weight(&vertex_group[i], index) > 0.0f;
+//     if (found != invert) {
+//       selected_indices.append(i);
+//     }
+//   }
+//   return selected_indices;
+// }
+
+// static Array<bool> selection_array_from_vertex_group(Span<MDeformVert> vertex_group,
+//                                                      const int index,
+//                                                      const bool invert)
+// {
+//   Array<bool> selection(vertex_group.size());
+//   for (const int i : vertex_group.index_range()) {
+//     const bool found = BKE_defvert_find_weight(&vertex_group[i], index) > 0.0f;
+//     selection[i] = (found != invert);
+//   }
+//   return selection;
+// }
+
+// static std::optional<Mesh *> calculate_weld(const Mesh &mesh, const WeldModifierData &wmd)
+// {
+//   const int defgrp_index = BKE_id_defgroup_name_index(&mesh.id, wmd.defgrp_name);
+//   Span<MDeformVert> vertex_group = get_vertex_group(mesh, defgrp_index);
+//   const bool invert = (wmd.flag & MOD_WELD_INVERT_VGROUP) != 0;
+
+//   if (wmd.mode == MOD_WELD_MODE_ALL) {
+//     if (!vertex_group.is_empty()) {
+//       Vector<int64_t> selected_indices = selected_indices_from_vertex_group(
+//           vertex_group, defgrp_index, invert);
+//       return blender::geometry::mesh_merge_by_distance_all(
+//           mesh, IndexMask(selected_indices), wmd.merge_dist);
+//     }
+//     return blender::geometry::mesh_merge_by_distance_all(
+//         mesh, IndexMask(mesh.totvert), wmd.merge_dist);
+//   }
+//   if (wmd.mode == MOD_WELD_MODE_CONNECTED) {
+//     const bool only_loose_edges = (wmd.flag & MOD_WELD_LOOSE_EDGES) != 0;
+//     if (!vertex_group.is_empty()) {
+//       Array<bool> selection = selection_array_from_vertex_group(
+//           vertex_group, defgrp_index, invert);
+//       return blender::geometry::mesh_merge_by_distance_connected(
+//           mesh, selection, wmd.merge_dist, only_loose_edges);
+//     }
+//     Array<bool> selection(mesh.totvert, true);
+//     return blender::geometry::mesh_merge_by_distance_connected(
+//         mesh, selection, wmd.merge_dist, only_loose_edges);
+//   }
+
+//   BLI_assert_unreachable();
+//   return nullptr;
+// }
+
+// static Mesh *modifyMesh(ModifierData *md, const ModifierEvalContext *UNUSED(ctx), Mesh *mesh)
+// {
+// const WeldModifierData &wmd = reinterpret_cast<WeldModifierData &>(*md);
+
+//   std::optional<Mesh *> result = calculate_weld(*mesh, wmd);
+//   if (!result) {
+//     return mesh;
+//   }
+//   return *result;
+// }
+
 static void initData(ModifierData *md)
 {
   WeldModifierData *wmd = (WeldModifierData *)md;

source/blender/simulation/CMakeLists.txt

@@ -29,6 +29,9 @@ set(SRC
   intern/eigen_utils.h
   intern/implicit.h
 
+  intern/xpbd.h
+  intern/xpbd.cc
+
   SIM_mass_spring.h
 )
 

source/blender/simulation/intern/xpbd.cc

@@ -0,0 +1,365 @@
+#include <bits/stdc++.h>
+
+#include "BKE_softbody.h"
+#include "DNA_object_force_types.h"
+
+#include "BLI_math.h"   
+
+#include "BLI_set.hh"
+
+#include "xpbd.h"
+
+using blender::Set;
+using namespace std;
+
+int tet_faces[4][4] = {{2,1,0,3}, {0,1,3,2}, {1,2,3,0}, {2,0,3,1}}; 
+// Why does tet_faces variable name not work?
+// Stores indices permutations that make up a face in a cyclic order. The 4th indice of each face is the point not in the face
+
+float get_tet_volume(BodyPoint *bpoint, BodyTet *curr_tet){
+  float diff1[3], diff2[3], diff3[3];
+
+  float vert0[3];
+  float vert1[3];
+  float vert2[3];
+  float vert3[3];
+
+  copy_v3_v3(vert0, bpoint[curr_tet->verts[0]].x);
+  copy_v3_v3(vert1, bpoint[curr_tet->verts[1]].x);
+  copy_v3_v3(vert2, bpoint[curr_tet->verts[2]].x);
+  copy_v3_v3(vert3, bpoint[curr_tet->verts[3]].x);
+
+  sub_v3_v3v3(diff1, bpoint[curr_tet->verts[1]].x, bpoint[curr_tet->verts[0]].x);
+  sub_v3_v3v3(diff2, bpoint[curr_tet->verts[2]].x, bpoint[curr_tet->verts[0]].x);
+  sub_v3_v3v3(diff3, bpoint[curr_tet->verts[3]].x, bpoint[curr_tet->verts[0]].x);
+  
+  float cross[3];
+  cross_v3_v3v3(cross, diff1, diff2);
+
+  return dot_v3v3(cross, diff3);
+}
+
+void xpbd_position_update(SoftBody *sb){
+    BodyPoint *bpoint_arr = sb->bpoint;
+    BodyTet *btet_arr = sb->btet;
+
+    int totpoint = sb->totpoint;
+    int tottet = sb->tottet;
+
+    float sdt = sb->dt/sb->substep_count;
+
+    for(int i = 0; i<totpoint; i++){
+        if(bpoint_arr[i].mass_inv == 0.0f){ 
+            continue;
+        }
+
+        copy_v3_v3(bpoint_arr[i].x_prev, bpoint_arr[i].x);
+        float temp[3];
+        copy_v3_fl3(temp, 0.0, 0.0, -sb->grav);
+        // add_v3_v3(temp, bpoint_arr[i].a);
+        mul_v3_fl(temp, sdt);
+        add_v3_v3(bpoint_arr[i].v, temp);
+
+        copy_v3_v3(temp, bpoint_arr[i].v);
+        mul_v3_fl(temp, sdt);
+        add_v3_v3(bpoint_arr[i].x, temp);
+
+        if(bpoint_arr[i].x[2] < 0.0){
+            // copy_v3_v3(bpoint_arr[i].x, bpoint_arr[i].x_prev);
+            bpoint_arr[i].x[2] = 0.0;
+        }
+    }
+}
+
+void xpbd_velocity_update(SoftBody *sb){
+    int totpoint = sb->totpoint;
+    BodyPoint *bpoint_arr = sb->bpoint;
+
+    float sdt = sb->dt/sb->substep_count;
+
+    for(int i = 0; i<totpoint; i++){
+        if(bpoint_arr[i].mass_inv == 0.0f){
+            continue;
+        }
+
+        sub_v3_v3v3(bpoint_arr[i].v, bpoint_arr[i].x, bpoint_arr[i].x_prev);
+        mul_v3_fl(bpoint_arr[i].v, 1.0/sdt);
+    }
+}
+
+void xpbd_distance_constraint_for_edge(BodyPoint *p1, BodyPoint *p2, float compliance){
+    float ini_length = len_v3v3(p1->x_ini, p2->x_ini);
+    float curr_length = len_v3v3(p1->x, p2->x);
+
+    float lambda = (curr_length - ini_length)/(curr_length * (p1->mass_inv + p2->mass_inv + compliance));
+    float dx_dir[3];
+    sub_v3_v3v3(dx_dir, p1->x, p2->x);
+    mul_v3_fl(dx_dir, lambda);  
+
+    float dx_p1[3] = {1.0f, 1.0f, 1.0f};
+    mul_v3_v3fl(dx_p1, dx_dir, p1->mass_inv);
+    sub_v3_v3(p1->x, dx_p1);
+                    
+    float dx_p2[3] = {1.0f, 1.0f, 1.0f};
+    mul_v3_v3fl(dx_p2, dx_dir, -p2->mass_inv);
+    sub_v3_v3(p2->x, dx_p2); 
+
+    float new_length = len_v3v3(p1->x, p2->x);
+}
+
+void xpbd_enforce_distance_constraint(SoftBody *sb){
+    BodyPoint *bpoint_arr = sb->bpoint;
+    BodyEdge *bedge_arr = sb->bedge;
+    int totedge = sb->totedge;
+
+    float sdt = sb->dt/sb->substep_count;
+
+    for(int i = 0; i<totedge; i++){
+        xpbd_distance_constraint_for_edge(&bpoint_arr[bedge_arr[i].u], &bpoint_arr[bedge_arr[i].v], sb->alpha_edge/(sdt*sdt));
+    }
+}
+
+void xpbd_enforce_volume_constraint(SoftBody *sb){
+    BodyPoint *bpoint_arr = sb->bpoint;
+    BodyTet *btet_arr = sb->btet;
+
+    float sdt = sb->dt/sb->substep_count;
+
+    int totpoint = sb->totpoint;
+    int tottet = sb->tottet;
+
+    float diff10[3], diff20[3], diff30[3], diff21[3], diff31[3], diff32[3];
+    float delC0[3], delC1[3], delC2[3], delC3[3];
+
+    for(int tetnr = 0; tetnr<tottet; tetnr++){
+        float lambda;
+
+        float curr_volume = get_tet_volume(bpoint_arr, &btet_arr[tetnr]);
+        float ini_volume = btet_arr[tetnr].initial_volume;
+
+        int vert0 = btet_arr[tetnr].verts[0];
+        int vert1 = btet_arr[tetnr].verts[1];
+        int vert2 = btet_arr[tetnr].verts[2];
+        int vert3 = btet_arr[tetnr].verts[3];
+
+        sub_v3_v3v3(diff10, bpoint_arr[vert1].x, bpoint_arr[vert0].x);
+        sub_v3_v3v3(diff20, bpoint_arr[vert2].x, bpoint_arr[vert0].x);
+        sub_v3_v3v3(diff30, bpoint_arr[vert3].x, bpoint_arr[vert0].x);
+        sub_v3_v3v3(diff21, bpoint_arr[vert2].x, bpoint_arr[vert1].x);
+        sub_v3_v3v3(diff31, bpoint_arr[vert3].x, bpoint_arr[vert1].x);
+        sub_v3_v3v3(diff32, bpoint_arr[vert3].x, bpoint_arr[vert2].x);
+
+        cross_v3_v3v3(delC0, diff31, diff21);
+        cross_v3_v3v3(delC1, diff20, diff30);
+        cross_v3_v3v3(delC2, diff30, diff10);
+        cross_v3_v3v3(delC3, diff10, diff20);
+
+        lambda = -1*(curr_volume - ini_volume)/(bpoint_arr[vert0].mass_inv * len_squared_v3(delC0) + 
+                                                bpoint_arr[vert1].mass_inv * len_squared_v3(delC1) + 
+                                                bpoint_arr[vert2].mass_inv * len_squared_v3(delC2) +
+                                                bpoint_arr[vert3].mass_inv * len_squared_v3(delC3) +
+                                                sb->alpha_vol/(sdt*sdt));
+
+        mul_v3_fl(delC0, lambda*bpoint_arr[vert0].mass_inv);
+        mul_v3_fl(delC1, lambda*bpoint_arr[vert1].mass_inv);
+        mul_v3_fl(delC2, lambda*bpoint_arr[vert2].mass_inv);
+        mul_v3_fl(delC3, lambda*bpoint_arr[vert3].mass_inv);
+        
+        add_v3_v3(bpoint_arr[vert0].x, delC0);
+        add_v3_v3(bpoint_arr[vert1].x, delC1);
+        add_v3_v3(bpoint_arr[vert2].x, delC2);
+        add_v3_v3(bpoint_arr[vert3].x, delC3);
+    }
+}
+
+void z_axis_collision_constraint(SoftBody *sb){
+    for(int i = 0; i<sb->totpoint; i++){
+        if(sb->bpoint[i].x[2] < 0.0)
+            sb->bpoint[i].x[2] = 0.0; 
+    }
+}
+
+static bool point_part_of_tet(int pointnr, BodyTet *tet){
+    for(int i = 0; i<4; i++){
+        if(tet->verts[i] == pointnr){
+            return true;
+        }
+    }
+
+    return false;
+}
+
+static bool same_side_of_tri(float *p1, float *p2, float *p3, float *p4, float *q){
+    float normal[3];
+    
+    float temp1[3];
+    sub_v3_v3v3(temp1, p2, p1);
+    float temp2[3];
+    sub_v3_v3v3(temp2, p3, p1);
+    cross_v3_v3v3(normal, temp1, temp2);
+
+    float diff_p4p1[3], diff_qp1[3];
+    sub_v3_v3v3(diff_p4p1, p4, p1);
+    sub_v3_v3v3(diff_qp1, q, p1);
+    float dotp4 = dot_v3v3(normal, diff_p4p1);
+    float dotq = dot_v3v3(normal, diff_qp1);
+
+    return (dotq/abs(dotq)) == (dotp4/abs(dotp4));
+    // Could floating point error cause problems here?
+}
+
+static bool point_inside_tet(BodyPoint *bpoint_arr, float *curr_point, BodyTet *btet){
+    for(int facenr = 0; facenr < 4; facenr++){
+        BodyPoint *face_point0 = &bpoint_arr[btet->verts[tet_faces[facenr][0]]];
+        BodyPoint *face_point1 = &bpoint_arr[btet->verts[tet_faces[facenr][1]]];
+        BodyPoint *face_point2 = &bpoint_arr[btet->verts[tet_faces[facenr][2]]];
+
+        BodyPoint *opposite_point = &bpoint_arr[btet->verts[tet_faces[facenr][3]]];
+
+        if(!same_side_of_tri(face_point0->x, face_point1->x, face_point2->x, opposite_point->x, curr_point)){
+            return false;
+        }
+    }
+
+    BodyPoint *tet_point0 = &bpoint_arr[btet->verts[0]];
+    BodyPoint *tet_point1 = &bpoint_arr[btet->verts[1]];
+    BodyPoint *tet_point2 = &bpoint_arr[btet->verts[2]];
+    BodyPoint *tet_point3 = &bpoint_arr[btet->verts[3]];
+
+    return true;
+}
+
+static int find_intersecting_face(BodyPoint *bpoint_arr, BodyPoint *curr_point, BodyTet *curr_tet){
+    float min_dist = FLT_MAX;
+    float min_face = -1;
+
+    float translation_vec[3];
+    sub_v3_v3v3(translation_vec, curr_point->x, curr_point->x_prev);
+    normalize_v3(translation_vec);
+
+    bool ahead_flag = !point_inside_tet(bpoint_arr, curr_point->x_prev, curr_tet);
+
+    for(int facenr = 0; facenr < 4; facenr++){
+        BodyPoint *face_point0 = &bpoint_arr[curr_tet->verts[tet_faces[facenr][0]]];
+        BodyPoint *face_point1 = &bpoint_arr[curr_tet->verts[tet_faces[facenr][1]]];
+        BodyPoint *face_point2 = &bpoint_arr[curr_tet->verts[tet_faces[facenr][2]]];
+        
+        float coeff[4];
+        float diff10[3], diff20[3];
+        sub_v3_v3v3(diff10, face_point1->x, face_point0->x);
+        sub_v3_v3v3(diff20, face_point2->x, face_point0->x);
+        cross_v3_v3v3(coeff, diff10, diff20);
+        coeff[3] = -dot_v3v3(coeff, face_point0->x);
+        
+        float lambda = -(dot_v3v3(coeff, curr_point->x_prev) + coeff[3])/dot_v3v3(coeff, translation_vec);
+        if(!ahead_flag){
+            lambda = -lambda;
+        }
+
+        if(lambda > 0.0f && lambda<min_dist){
+            min_dist = lambda;
+            min_face = facenr;
+        }
+    }
+
+    return min_face;
+
+    // Checking which is the closest face
+    // float min_dist = FLT_MAX;
+    // float min_face = -1;
+
+    // for(int facenr = 0; facenr < 4; facenr++){
+    //     BodyPoint *face_point0 = &bpoint_arr[curr_tet->verts[tet_faces[facenr][0]]];
+    //     BodyPoint *face_point1 = &bpoint_arr[curr_tet->verts[tet_faces[facenr][1]]];
+    //     BodyPoint *face_point2 = &bpoint_arr[curr_tet->verts[tet_faces[facenr][2]]];
+        
+    //     float coeff[4];
+    //     float diff10[3], diff20[3];
+    //     sub_v3_v3v3(diff10, face_point1->x, face_point0->x);
+    //     sub_v3_v3v3(diff20, face_point2->x, face_point0->x);
+    //     cross_v3_v3v3(coeff, diff10, diff20);
+    //     coeff[3] = -dot_v3v3(coeff, face_point0->x);
+    // }
+    
+    // Is the point inside tet getting detected correctly?
+    // Is the intersecting face getting detected correctly?
+    // Is the correction working as expected?
+}
+
+void xpbd_solve_self_collision(SoftBody *sb){
+    BodyPoint *bpoint_arr = sb->bpoint;
+    BodyTet *btet_arr = sb->btet;
+
+    int totpoint = sb->totpoint;
+    int tottet = sb->tottet;
+
+    // for(int pointnr = 0; pointnr<sb->tot_surface_point; pointnr++){
+    for(int pointid = 0; pointid < sb->tot_surface_point; pointid++){
+        int pointnr = sb->surface_points[pointid];
+        BodyPoint *curr_point = &bpoint_arr[pointnr];
+
+        // for(int tetnr = 0; tetnr < tottet; tetnr++){
+        for(int tetid = 0; tetid < sb->tot_surface_tet; tetid++){
+            int tetnr = sb->surface_tets[tetid];
+            BodyTet *curr_tet = &btet_arr[tetnr];
+            if(point_part_of_tet(pointnr, curr_tet)){
+                continue;
+            }
+
+            if(!point_inside_tet(bpoint_arr, curr_point->x, curr_tet)){
+                continue;
+            }
+
+            // Check if x_prev of point is also inside tet?
+
+            int intersecting_face = find_intersecting_face(bpoint_arr, curr_point, curr_tet);
+
+            BodyPoint *p0 = &bpoint_arr[curr_tet->verts[tet_faces[intersecting_face][0]]];
+            BodyPoint *p1 = &bpoint_arr[curr_tet->verts[tet_faces[intersecting_face][1]]];
+            BodyPoint *p2 = &bpoint_arr[curr_tet->verts[tet_faces[intersecting_face][2]]];
+            BodyPoint *p3 = &bpoint_arr[curr_tet->verts[tet_faces[intersecting_face][3]]];
+            BodyPoint *q = curr_point;
+
+            float diff10[3], diff20[3], diffq0[3], diff30[3];
+            sub_v3_v3v3(diff10, p1->x, p0->x);
+            sub_v3_v3v3(diff20, p2->x, p0->x);
+            sub_v3_v3v3(diffq0, q->x, p0->x);
+            
+
+            float delC[3]; // Points outside the tet
+            cross_v3_v3v3(delC, diff10, diff20);
+
+            sub_v3_v3v3(diff30, p0->x, p3->x);
+            if(dot_v3v3(diff30, delC) < 0){
+                mul_v3_fl(delC, -1);
+            }
+            normalize_v3(delC);
+
+            float lambda = dot_v3v3(diffq0, delC)/(p0->mass_inv + p1->mass_inv + p2->mass_inv + q->mass_inv);
+            lambda /= len_squared_v3(delC);
+
+            mul_v3_fl(delC, lambda);
+
+            float dx_q[3];
+            copy_v3_v3(dx_q, delC);
+            mul_v3_fl(dx_q, -q->mass_inv);  
+            add_v3_v3(q->x, dx_q);
+
+            for(int i = 0; i<3; i++){
+                BodyPoint *point = &bpoint_arr[curr_tet->verts[tet_faces[intersecting_face][i]]];
+
+                float dx_point[3];
+                copy_v3_v3(dx_point, delC);
+                mul_v3_fl(dx_point, point->mass_inv);
+                add_v3_v3(point->x, dx_point);
+            }
+        }
+    }
+}
+
+void xpbd_enforce_constraints(SoftBody *sb){
+    // z_axis_collision_constraint(sb);
+    xpbd_enforce_distance_constraint(sb);
+    xpbd_enforce_volume_constraint(sb);
+    xpbd_solve_self_collision(sb);
+}

source/blender/simulation/intern/xpbd.h

@@ -0,0 +1,25 @@
+// #include "BKE_softbody.h"
+// #include "DNA_object_force_types.h"
+
+#pragma once
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct BodyPoint;
+struct BodyTet;
+struct SoftBody;
+
+void xpbd_position_update(struct SoftBody *sb);
+void xpbd_distance_constraint_for_edge(struct BodyPoint *p1, struct BodyPoint *p2, float compliance);
+void xpbd_enforce_distance_constraint(struct SoftBody *sb);
+void xpbd_enforce_volume_constraint(struct SoftBody *sb);
+void z_axis_collision_constraint(struct SoftBody *sb);
+void xpbd_solve_self_collision(struct SoftBody *sb);
+void xpbd_enforce_constraints(struct SoftBody *sb);
+void xpbd_velocity_update(struct SoftBody *sb);
+
+#ifdef __cplusplus
+}
+#endif