Checkpoint.

EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodInitConstraints.comp

@@ -8,8 +8,6 @@
 #define DYNAMIC_STRANDS_SET 0
 #include "DynamicStrands.glsl"
 
-#include "Math.glsl"
-
 const uint WORKGROUP_SIZE = EXT_TASK_SUBGROUP_COUNT * EXT_TASK_SUBGROUP_SIZE;
 
 const uint ITERATIONS_PER_TASK = ((EXT_INVOCATIONS_PER_TASK + WORKGROUP_SIZE - 1) / WORKGROUP_SIZE);

EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodProjectConstraints.comp

@@ -7,8 +7,6 @@
 #define DYNAMIC_STRANDS_SET 0
 #include "DynamicStrands.glsl"
 
-#include "Math.glsl"
-
 const uint WORKGROUP_SIZE = EXT_TASK_SUBGROUP_COUNT * EXT_TASK_SUBGROUP_SIZE;
 
 const uint ITERATIONS_PER_TASK = ((EXT_INVOCATIONS_PER_TASK + WORKGROUP_SIZE - 1) / WORKGROUP_SIZE);
@@ -70,8 +68,8 @@ void project_constraint(in int rod_constraint_handle){
 	StrandSegment segment1 = strand_segments[segment1_handle];
 
 	mat3 inverse_inertia_tensor_w0, inverse_inertia_tensor_w1;
-	inverse_inertia_tensor_w0 = compute_inverse_inertia_tensor_w(segment0);
-	inverse_inertia_tensor_w1 = compute_inverse_inertia_tensor_w(segment1);
+	inverse_inertia_tensor_w0 = mat3(segment0.inv_inertia_w);
+	inverse_inertia_tensor_w1 = mat3(segment1.inv_inertia_w);
 
 	solve_stretch_bending_twisting_constraints(
 		segment0.inv_inertia_tensor_inv_mass.w, segment0.x.xyz, inverse_inertia_tensor_w0, segment0.q,

EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodUpdateConstraints.comp

@@ -8,8 +8,6 @@
 #define DYNAMIC_STRANDS_SET 0
 #include "DynamicStrands.glsl"
 
-#include "Math.glsl"
-
 const uint WORKGROUP_SIZE = EXT_TASK_SUBGROUP_COUNT * EXT_TASK_SUBGROUP_SIZE;
 
 const uint ITERATIONS_PER_TASK = ((EXT_INVOCATIONS_PER_TASK + WORKGROUP_SIZE - 1) / WORKGROUP_SIZE);

EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Operators/ExternalForce.comp

@@ -7,8 +7,7 @@
 #define DYNAMIC_STRANDS_SET 0
 #include "DynamicStrands.glsl"
 
-#define DYNAMIC_STRANDS_PHYSICS_SET 1
-#include "DynamicStrandsPhysics.glsl"
+#include "Math.glsl"
 
 const uint WORKGROUP_SIZE = EXT_TASK_SUBGROUP_COUNT * EXT_TASK_SUBGROUP_SIZE;
 
@@ -24,10 +23,10 @@ memoryBarrierShared(); \
 barrier();
 
 struct ExternalForce{
- vec4 index_force;
+ vec4 force_index;
 };
 
-layout(std430, set = 2, binding = 0) readonly buffer EXTERNAL_FORCE_BLOCK {
+layout(std430, set = 1, binding = 0) readonly buffer EXTERNAL_FORCE_BLOCK {
   ExternalForce external_forces[];
 };
 
@@ -44,17 +43,10 @@ void main(){
 
 		if(global_index >= command_size) break;
 
-		//strand_segments[global_index].end_position_thickness.y = 0.0;
 		ExternalForce external_force = external_forces[global_index];
-		uint strand_particle_index = floatBitsToUint(external_force.index_force.x);
-		vec3 force = external_force.index_force.yzw;
-
-		StrandSegmentParticlePhysics particle_physics = strand_segment_particle_physics_list[strand_particle_index];
-		float inv_mass = particle_physics.last_position_inv_mass.w;
-		vec3 acceleration = force * inv_mass;
-		vec4 position_thickness = strand_segment_particles[strand_particle_index].position_thickness;
-
-		vec3 velocity = (position_thickness.xyz - particle_physics.last_position_inv_mass.xyz) + acceleration * delta_time * .5;
-		position_thickness.xyz = position_thickness.xyz + velocity * delta_time;
+		uint segment_index = floatBitsToUint(external_force.force_index.w);
+		vec3 force = external_force.force_index.xyz;
+		StrandSegment segment = strand_segments[segment_index];
+		strand_segments[segment_index].acceleration.xyz = force * segment.inv_inertia_tensor_inv_mass.w;
 	}
 }

EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Operators/PositionUpdate.comp

@@ -7,6 +7,8 @@
 #define DYNAMIC_STRANDS_SET 0
 #include "DynamicStrands.glsl"
 
+#include "Math.glsl"
+
 const uint WORKGROUP_SIZE = EXT_TASK_SUBGROUP_COUNT * EXT_TASK_SUBGROUP_SIZE;
 
 const uint ITERATIONS_PER_TASK = ((EXT_INVOCATIONS_PER_TASK + WORKGROUP_SIZE - 1) / WORKGROUP_SIZE);

EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/PostStep.comp

@@ -6,6 +6,7 @@
 
 #define DYNAMIC_STRANDS_SET 0
 #include "DynamicStrands.glsl"
+#include "Math.glsl"
 
 const uint WORKGROUP_SIZE = EXT_TASK_SUBGROUP_COUNT * EXT_TASK_SUBGROUP_SIZE;
 
@@ -20,32 +21,18 @@ layout(local_size_x = WORKGROUP_SIZE, local_size_y = 1, local_size_z = 1) in;
 memoryBarrierShared(); \
 barrier();
 
-struct PositionUpdate{
-	vec4 new_position_index;
-};
-
-layout(std430, set = 1, binding = 0) readonly buffer POSITION_UPDATE_BLOCK {
-  PositionUpdate position_updates[];
-};
 
 layout(push_constant) uniform PUSH_CONSTANTS {
-	uint command_size;
+	uint segment_size;
 };
 
 void main(){
 	for (uint i = 0; i < ITERATIONS_PER_TASK; i++)
 	{
 		uint local_index  = laneID + i * WORKGROUP_SIZE;
-		uint global_index = baseID + local_index;
-
-		if(global_index >= command_size) break;
-
-		PositionUpdate position_update = position_updates[global_index];
-
-		uint segment_index = floatBitsToUint(position_update.new_position_index.w);
-
-		StrandSegment strand_segment = strand_segments[segment_index];
-
+		uint segment_index = baseID + local_index;
 
+		if(segment_index >= segment_size) break;
+		StrandSegment segment = strand_segments[segment_index];
 	}
 }

EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/PreStep.comp

@@ -6,6 +6,7 @@
 
 #define DYNAMIC_STRANDS_SET 0
 #include "DynamicStrands.glsl"
+#include "Math.glsl"
 
 const uint WORKGROUP_SIZE = EXT_TASK_SUBGROUP_COUNT * EXT_TASK_SUBGROUP_SIZE;
 
@@ -20,32 +21,86 @@ layout(local_size_x = WORKGROUP_SIZE, local_size_y = 1, local_size_z = 1) in;
 memoryBarrierShared(); \
 barrier();
 
-struct PositionUpdate{
-	vec4 new_position_index;
+layout(push_constant) uniform PUSH_CONSTANTS {
+	uint segment_size;
+	float time_step;
+	float inv_time_step;
 };
 
-layout(std430, set = 1, binding = 0) readonly buffer POSITION_UPDATE_BLOCK {
-  PositionUpdate position_updates[];
-};
+void update_inertia_w(inout StrandSegment segment){
+	//Update w
+	mat3 rot = mat3_cast(segment.q);
+	mat3 inertia_tensor_diag =
+		mat3(
+			segment.inertia_tensor_mass.x, 0.0, 0.0, 
+			0.0, segment.inertia_tensor_mass.y, 0.0, 
+			0.0, 0.0, segment.inertia_tensor_mass.z);
+	mat3 inertia_w = rot * inertia_tensor_diag * transpose(rot);
 
-layout(push_constant) uniform PUSH_CONSTANTS {
-	uint command_size;
-};
+	mat3 inverse_inertia_tensor_diag =
+		mat3(
+			segment.inv_inertia_tensor_inv_mass.x, 0.0, 0.0, 
+			0.0, segment.inv_inertia_tensor_inv_mass.y, 0.0, 
+			0.0, 0.0, segment.inv_inertia_tensor_inv_mass.z);
+	mat3 inverse_inertia_w = rot * inverse_inertia_tensor_diag * transpose(rot);
+
+	segment.inertia_w = mat4(
+		inertia_w[0][0], inertia_w[0][1], inertia_w[0][2], 0.0,
+		inertia_w[1][0], inertia_w[1][1], inertia_w[1][2], 0.0,
+		inertia_w[2][0], inertia_w[2][1], inertia_w[2][2], 0.0,
+		0.0, 0.0, 0.0, 0.0
+		);
+	segment.inv_inertia_w = mat4(
+		inverse_inertia_w[0][0], inverse_inertia_w[0][1], inverse_inertia_w[0][2], 0.0,
+		inverse_inertia_w[1][0], inverse_inertia_w[1][1], inverse_inertia_w[1][2], 0.0,
+		inverse_inertia_w[2][0], inverse_inertia_w[2][1], inverse_inertia_w[2][2], 0.0,
+		0.0, 0.0, 0.0, 0.0
+		);
+}
 
 void main(){
 	for (uint i = 0; i < ITERATIONS_PER_TASK; i++)
 	{
 		uint local_index  = laneID + i * WORKGROUP_SIZE;
-		uint global_index = baseID + local_index;
+		uint segment_index = baseID + local_index;
 
-		if(global_index >= command_size) break;
+		if(segment_index >= segment_size) break;
 
-		PositionUpdate position_update = position_updates[global_index];
+		StrandSegment segment = strand_segments[segment_index];
+		update_inertia_w(segment);
 
-		uint segment_index = floatBitsToUint(position_update.new_position_index.w);
+		//Calculate velocity and apply acceleration.
+		vec3 velocity = vec3(0, 0, 0);
+		if(segment.inv_inertia_tensor_inv_mass.w != 0.0){
+			velocity = inv_time_step * (1.5 * segment.x.xyz - 2.0 * segment.old_x.xyz + 0.5 * segment.last_x.xyz);
+			velocity += segment.acceleration.xyz * time_step;
+		}
+		//Calculate angular velocity and apply torque.
+		vec3 angular_velocity = vec3(0, 0, 0);
+		if(segment.inv_inertia_tensor_inv_mass.w != 0.0){
+			vec4 rot = quat_mul(segment.q, conjugate(segment.old_q));
+			angular_velocity = rot.xyz * 2.0 * inv_time_step;
+			angular_velocity += time_step * mat3(segment.inv_inertia_w) * (segment.torque.xyz - cross(angular_velocity, mat3(segment.inertia_w) * angular_velocity));
+		}
 
-		StrandSegment strand_segment = strand_segments[segment_index];
+		//Shift position values
+		segment.last_x.xyz = segment.old_x.xyz;
+		segment.old_x.xyz = segment.x.xyz;
+		//Apply velocity
+		segment.x.xyz += time_step * velocity;
+
+
+		//Shift rotation values
+		segment.last_q = segment.old_q;
+		segment.old_q = segment.q;
+		//Apply angular velocity
+		vec4 angular_velocity_q = vec4(angular_velocity.x, angular_velocity.y, angular_velocity.z, 0.0);
 
+		segment.q += 0.5 * time_step * quat_mul(angular_velocity_q, segment.q);
+		segment.q = normalize(segment.q);
 
+		//Update w
+		update_inertia_w(segment);
+		strand_segments[segment_index] = segment;
 	}
 }

EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Includes/DynamicStrands.glsl

@@ -11,16 +11,21 @@ struct StrandSegment {
   vec4 q;
   vec4 last_q;
   vec4 old_q;
+  vec4 torque;
 
   vec4 x0_length;
   vec4 x;
   vec4 last_x;
   vec4 old_x;
+  vec4 acceleration;
 
   vec4 inertia_tensor_mass;
 
   vec4 inv_inertia_tensor_inv_mass;
 
+  mat4 inertia_w;
+  mat4 inv_inertia_w;
+
   int neighbors[8];
 };
 

EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Includes/DynamicStrandsPhysics.glsl

@@ -1,4 +1,5 @@
 
+#include "Math.glsl"
 struct RodProperties {
   float density;
   float i1;

EvoEngine_Plugins/EcoSysLab/include/StrandModel/BrokenBranch.hpp

@@ -21,6 +21,7 @@ class BrokenBranch : public IPrivateComponent {
   void Deserialize(const YAML::Node& in) override;
   bool OnInspect(const std::shared_ptr<EditorLayer>& editor_layer) override;
   void LateUpdate() override;
+  void FixedUpdate() override;
   void OnCreate() override;
   void OnDestroy() override;
   void CollectAssetRef(std::vector<AssetRef>& list) override;

EvoEngine_Plugins/EcoSysLab/include/StrandModel/DynamicStrands.hpp

@@ -75,6 +75,8 @@ class DynamicStrands {
     // Last frame rotation
     glm::quat last_q;
     glm::quat old_q;
+    glm::vec3 torque = glm::vec3(0.f);
+    float padding0;
 
     // Initial position
     glm::vec3 x0;
@@ -85,12 +87,18 @@ class DynamicStrands {
     glm::vec4 last_x;
     glm::vec4 old_x;
 
+    glm::vec3 acceleration = glm::vec3(0.f);
+    float padding1;
+
     glm::vec3 inertia_tensor;
     float mass;
 
     glm::vec3 inv_inertia_tensor;
     float inv_mass;
 
+    glm::mat4 inertia_w;
+    glm::mat4 inverse_inertia_w;
+
     int neighbors[8];
   };
 
@@ -111,6 +119,9 @@ class DynamicStrands {
   std::vector<std::shared_ptr<DescriptorSet>> strands_descriptor_sets;
 
  private:
+
+  static glm::vec3 ComputeInertiaTensorBox(float mass, float width, float height, float depth);
+  static glm::vec3 ComputeInertiaTensorRod(float mass, float radius, float length);
   void InitializeOperators(const InitializeParameters& initialize_parameters) const;
   void InitializeConstraints(const InitializeParameters& initialize_parameters) const;
   void Render(const RenderParameters& render_parameters) const;
@@ -158,20 +169,21 @@ void DynamicStrands::InitializeStrandsGroup(const InitializeParameters& initiali
 
     segment.q0 = segment.q = segment.last_q = segment.old_q =
         initialize_parameters.root_transform.GetRotation() * target_strand_segment.rotation;
-
+    segment.torque = glm::vec3(0.f);
+
     segment.x = segment.last_x = segment.old_x =
         glm::vec4(initialize_parameters.root_transform.TransformPoint(strand_group.GetStrandSegmentCenter(static_cast<int>(i))), 0.0);
 
     segment.x0 = segment.x;
+    segment.acceleration = glm::vec3(0.f);
 
-    if (initialize_parameters.static_root) {
+    if (segment.prev_handle == -1 && initialize_parameters.static_root) {
       segment.mass = segment.inv_mass = 0.0f;
       segment.inertia_tensor = segment.inv_inertia_tensor = glm::vec3(0.0f);
     }else {
       segment.mass = segment.radius * segment.radius * glm::pi<float>() * initialize_parameters.wood_density;
       segment.inv_mass = 1.f / segment.mass;
-      const float inertia = .2f * segment.mass * segment.radius * segment.radius;
-      segment.inertia_tensor = glm::vec3(inertia, inertia, 2.f * inertia);
+      segment.inertia_tensor = ComputeInertiaTensorRod(segment.mass, segment.radius, segment.length);
       segment.inv_inertia_tensor = 1.f / segment.inertia_tensor;
     }
 

EvoEngine_Plugins/EcoSysLab/include/StrandModel/DynamicStrandsPhysics.hpp

@@ -9,6 +9,8 @@ class DynamicStrandsPreStep {
 
   struct PreStepPushConstant {
     uint32_t segment_size = 0;
+    float time_step = 0.01f;
+    float inv_time_step = 100.f;
   };
 
   inline static std::shared_ptr<ComputePipeline> pre_step_pipeline;
@@ -73,8 +75,8 @@ class DsExternalForce final : public IDynamicStrandsOperator {
   DsExternalForce();
 
   struct ExternalForce {
-    uint32_t strand_particle_index = 0;
     glm::vec3 force = glm::vec3(0.f);
+    uint32_t strand_segment_index = 0;
   };
 
   inline static std::shared_ptr<DescriptorSetLayout> layout{};