diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodInitConstraints.comp b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodInitConstraints.comp new file mode 100644 index 0000000..e04fbbf --- /dev/null +++ b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodInitConstraints.comp @@ -0,0 +1,73 @@ +#extension GL_EXT_shader_explicit_arithmetic_types_int8 : require +#extension GL_KHR_shader_subgroup_basic : require +#extension GL_KHR_shader_subgroup_ballot : require +#extension GL_KHR_shader_subgroup_vote : require +#extension GL_ARB_shading_language_include : enable +#extension GL_EXT_control_flow_attributes : require + +#define DYNAMIC_STRANDS_SET 0 +#include "DynamicStrands.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); + +uint baseID = gl_WorkGroupID.x * EXT_INVOCATIONS_PER_TASK; +uint laneID = gl_LocalInvocationID.x; + +layout(local_size_x = WORKGROUP_SIZE, local_size_y = 1, local_size_z = 1) in; + +#define BARRIER() \ +memoryBarrierShared(); \ +barrier(); + +layout(push_constant) uniform PUSH_CONSTANTS { + uint constraint_size; + float time_step; + float inv_time_step; +}; + +#define DYNAMIC_STRANDS_PHYSICS_SET 1 +#include "DynamicStrandsPhysics.glsl" + +void init_constraint(in uint rod_constraint_handle); +//Line 52 +void main(){ + + [[unroll]] + for (uint i = 0; i < ITERATIONS_PER_TASK; i++) + { + uint local_index = laneID + i * WORKGROUP_SIZE; + uint constraint_index = baseID + local_index; + + if(constraint_index >= constraint_size) break; + init_constraint(constraint_index); + } +} + +void init_constraint(in uint rod_constraint_handle){ + RodConstraint rod_constraint = rod_constraints[rod_constraint_handle]; + + vec3 stretch_compliance; + vec3 bending_and_torsion_compliance; + vec6 lambda_sum; + + init_stretch_bending_twisting_constraint( + rod_constraint.stiffness_coefficient_k_segment0_index.xyz, + inv_time_step, + rod_constraint.stretch_compliance_average_segment_length.w, + stretch_compliance, + bending_and_torsion_compliance, + lambda_sum); + + rod_constraints[rod_constraint_handle].stretch_compliance_average_segment_length.xyz = stretch_compliance; + rod_constraints[rod_constraint_handle].bending_and_torsion_compliance_next_constraint_handle.xyz = bending_and_torsion_compliance; + + rod_constraints[rod_constraint_handle].lambda_sum0 = lambda_sum.v[0]; + rod_constraints[rod_constraint_handle].lambda_sum1 = lambda_sum.v[1]; + rod_constraints[rod_constraint_handle].lambda_sum2 = lambda_sum.v[2]; + rod_constraints[rod_constraint_handle].lambda_sum3 = lambda_sum.v[3]; + rod_constraints[rod_constraint_handle].lambda_sum4 = lambda_sum.v[4]; + rod_constraints[rod_constraint_handle].lambda_sum5 = lambda_sum.v[5]; + +} \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodProjectConstraints.comp b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodProjectConstraints.comp new file mode 100644 index 0000000..a04c6f5 --- /dev/null +++ b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodProjectConstraints.comp @@ -0,0 +1,102 @@ +#extension GL_EXT_shader_explicit_arithmetic_types_int8 : require +#extension GL_KHR_shader_subgroup_basic : require +#extension GL_KHR_shader_subgroup_ballot : require +#extension GL_KHR_shader_subgroup_vote : require +#extension GL_ARB_shading_language_include : enable +#extension GL_EXT_control_flow_attributes : require +#define DYNAMIC_STRANDS_SET 0 +#include "DynamicStrands.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); + +uint baseID = gl_WorkGroupID.x * EXT_INVOCATIONS_PER_TASK; +uint laneID = gl_LocalInvocationID.x; + +layout(local_size_x = WORKGROUP_SIZE, local_size_y = 1, local_size_z = 1) in; + +#define BARRIER() \ +memoryBarrierShared(); \ +barrier(); + +layout(push_constant) uniform PUSH_CONSTANTS { + uint strand_size; +}; + +#define DYNAMIC_STRANDS_PHYSICS_SET 1 +#include "DynamicStrandsPhysics.glsl" + +void project_constraint(in int rod_constraint_handle); + +void main(){ + + [[unroll]] + for (uint i = 0; i < ITERATIONS_PER_TASK; i++) + { + uint local_index = laneID + i * WORKGROUP_SIZE; + uint strand_index = baseID + local_index; + + if(strand_index >= strand_size) break; + + PerStrandData strand_data = per_strand_data_list[strand_index]; + int rod_constraint_handle = strand_data.begin_rod_constraint_handle; + while(rod_constraint_handle != strand_data.end_rod_constraint_handle){ + project_constraint(rod_constraint_handle); + rod_constraint_handle = floatBitsToInt(rod_constraints[rod_constraint_handle].bending_and_torsion_compliance_next_constraint_handle.w); + } + project_constraint(rod_constraint_handle); + } +} + +void project_constraint(in int rod_constraint_handle){ + RodConstraint rod_constraint = rod_constraints[rod_constraint_handle]; + + vec3 x0_correction, x1_correction; + vec4 q0_correction, q1_correction; + vec6 lambda_sum; + lambda_sum.v[0] = rod_constraint.lambda_sum0; + lambda_sum.v[1] = rod_constraint.lambda_sum1; + lambda_sum.v[2] = rod_constraint.lambda_sum2; + lambda_sum.v[3] = rod_constraint.lambda_sum3; + lambda_sum.v[4] = rod_constraint.lambda_sum4; + lambda_sum.v[5] = rod_constraint.lambda_sum5; + + int segment0_handle = floatBitsToInt(rod_constraint.stiffness_coefficient_k_segment0_index.w); + int segment1_handle = floatBitsToInt(rod_constraint.rest_darboux_vector_segment1_index.w); + StrandSegment segment0 = strand_segments[segment0_handle]; + StrandSegment segment1 = strand_segments[segment1_handle]; + + mat3 inverse_inertia_tensor_w0, inverse_inertia_tensor_w1; + 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, + segment1.inv_inertia_tensor_inv_mass.w, segment1.x.xyz, inverse_inertia_tensor_w1, segment1.q, + rod_constraint.rest_darboux_vector_segment1_index.xyz, + rod_constraint.stretch_compliance_average_segment_length.w, + rod_constraint.stretch_compliance_average_segment_length.xyz, + rod_constraint.bending_and_torsion_compliance_next_constraint_handle.xyz, + rod_constraint.constraint_info, + x0_correction, q0_correction, + x1_correction, q1_correction, + lambda_sum); + + if(segment0.inv_inertia_tensor_inv_mass.w != 0.0){ + strand_segments[segment0_handle].x.xyz = segment0.x.xyz + x0_correction; + strand_segments[segment0_handle].q = normalize(segment0.q + q0_correction); + } + + if(segment1.inv_inertia_tensor_inv_mass.w != 0.0){ + strand_segments[segment1_handle].x.xyz = segment1.x.xyz + x1_correction; + strand_segments[segment1_handle].q = normalize(segment1.q + q1_correction); + } + + rod_constraints[rod_constraint_handle].lambda_sum0 = lambda_sum.v[0]; + rod_constraints[rod_constraint_handle].lambda_sum1 = lambda_sum.v[1]; + rod_constraints[rod_constraint_handle].lambda_sum2 = lambda_sum.v[2]; + rod_constraints[rod_constraint_handle].lambda_sum3 = lambda_sum.v[3]; + rod_constraints[rod_constraint_handle].lambda_sum4 = lambda_sum.v[4]; + rod_constraints[rod_constraint_handle].lambda_sum5 = lambda_sum.v[5]; +} \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodUpdateConstraints.comp b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodUpdateConstraints.comp new file mode 100644 index 0000000..d4cb359 --- /dev/null +++ b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Constraints/StiffRodUpdateConstraints.comp @@ -0,0 +1,55 @@ +#extension GL_EXT_shader_explicit_arithmetic_types_int8 : require +#extension GL_KHR_shader_subgroup_basic : require +#extension GL_KHR_shader_subgroup_ballot : require +#extension GL_KHR_shader_subgroup_vote : require +#extension GL_ARB_shading_language_include : enable +#extension GL_EXT_control_flow_attributes : require + +#define DYNAMIC_STRANDS_SET 0 +#include "DynamicStrands.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); + +uint baseID = gl_WorkGroupID.x * EXT_INVOCATIONS_PER_TASK; +uint laneID = gl_LocalInvocationID.x; + +layout(local_size_x = WORKGROUP_SIZE, local_size_y = 1, local_size_z = 1) in; + +#define BARRIER() \ +memoryBarrierShared(); \ +barrier(); + +layout(push_constant) uniform PUSH_CONSTANTS { + uint constraint_size; +}; + +#define DYNAMIC_STRANDS_PHYSICS_SET 1 +#include "DynamicStrandsPhysics.glsl" + +void update_constraint(in uint rod_constraint_handle); + +void main(){ + + [[unroll]] + for (uint i = 0; i < ITERATIONS_PER_TASK; i++) + { + uint local_index = laneID + i * WORKGROUP_SIZE; + uint constraint_index = baseID + local_index; + + if(constraint_index >= constraint_size) break; + update_constraint(constraint_index); + } +} + +void update_constraint(in uint rod_constraint_handle){ + RodConstraint rod_constraint = rod_constraints[rod_constraint_handle]; + mat4 constraint_info = rod_constraint.constraint_info; + + StrandSegment strand_segment0 = strand_segments[floatBitsToInt(rod_constraint.stiffness_coefficient_k_segment0_index.w)]; + StrandSegment strand_segment1 = strand_segments[floatBitsToInt(rod_constraint.rest_darboux_vector_segment1_index.w)]; + + update_stretch_bending_twisting_constraint(strand_segment0.x.xyz, strand_segment0.q, strand_segment1.x.xyz, strand_segment1.q, constraint_info); + rod_constraints[rod_constraint_handle].constraint_info = constraint_info; +} \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Operators/ExternalForce.comp b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Operators/ExternalForce.comp new file mode 100644 index 0000000..bcc60d5 --- /dev/null +++ b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Operators/ExternalForce.comp @@ -0,0 +1,52 @@ +#extension GL_EXT_shader_explicit_arithmetic_types_int8 : require +#extension GL_KHR_shader_subgroup_basic : require +#extension GL_KHR_shader_subgroup_ballot : require +#extension GL_KHR_shader_subgroup_vote : require +#extension GL_ARB_shading_language_include : enable + +#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); + +uint baseID = gl_WorkGroupID.x * EXT_INVOCATIONS_PER_TASK; +uint laneID = gl_LocalInvocationID.x; + +layout(local_size_x = WORKGROUP_SIZE, local_size_y = 1, local_size_z = 1) in; + +#define BARRIER() \ +memoryBarrierShared(); \ +barrier(); + +struct ExternalForce{ + vec4 force_index; +}; + +layout(std430, set = 1, binding = 0) readonly buffer EXTERNAL_FORCE_BLOCK { + ExternalForce external_forces[]; +}; + +layout(push_constant) uniform PUSH_CONSTANTS { + uint command_size; + float delta_time; +}; + +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; + + ExternalForce external_force = external_forces[global_index]; + 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; + } +} \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Operators/PositionUpdate.comp b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Operators/PositionUpdate.comp new file mode 100644 index 0000000..b20fdac --- /dev/null +++ b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/Operators/PositionUpdate.comp @@ -0,0 +1,56 @@ +#extension GL_EXT_shader_explicit_arithmetic_types_int8 : require +#extension GL_KHR_shader_subgroup_basic : require +#extension GL_KHR_shader_subgroup_ballot : require +#extension GL_KHR_shader_subgroup_vote : require +#extension GL_ARB_shading_language_include : enable + +#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); + +uint baseID = gl_WorkGroupID.x * EXT_INVOCATIONS_PER_TASK; +uint laneID = gl_LocalInvocationID.x; + +layout(local_size_x = WORKGROUP_SIZE, local_size_y = 1, local_size_z = 1) in; + +#define BARRIER() \ +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; +}; + +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); + + vec3 new_position = position_update.new_position_index.xyz; + vec3 old_position = strand_segments[segment_index].x.xyz; + vec3 delta_position = new_position - old_position; + + strand_segments[segment_index].x.xyz = new_position; + strand_segments[segment_index].last_x.xyz = strand_segments[segment_index].last_x.xyz + delta_position; + } +} \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrandsPrecompute.comp b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/PostStep.comp similarity index 63% rename from EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrandsPrecompute.comp rename to EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/PostStep.comp index bdd0e4d..7a9e602 100644 --- a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrandsPrecompute.comp +++ b/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,28 +21,18 @@ layout(local_size_x = WORKGROUP_SIZE, local_size_y = 1, local_size_z = 1) in; memoryBarrierShared(); \ barrier(); -struct PhysicsCommand{ - int strand_segment_handle; - vec3 position_offset; -}; - -//#define PHYSICS_COMMANDS_SET 0 -//layout(std430, set = PHYSICS_COMMANDS_SET, binding = 0) readonly buffer PHYSICS_COMMAND_BLOCK { -// PhysicsCommand physics_commands[]; -//}; -layout(push_constant) uniform STRANDS_PRECOMPUTE_CONSTANTS { - uint strands_size; - uint strand_segments_size; - uint command_size; +layout(push_constant) uniform PUSH_CONSTANTS { + 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; - //strand_segments[global_index].end_position_thickness.y = 0.0; + uint segment_index = baseID + local_index; + + if(segment_index >= segment_size) break; + StrandSegment segment = strand_segments[segment_index]; } } \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/PreStep.comp b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/PreStep.comp new file mode 100644 index 0000000..ff07082 --- /dev/null +++ b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Compute/DynamicStrands/PreStep.comp @@ -0,0 +1,106 @@ +#extension GL_EXT_shader_explicit_arithmetic_types_int8 : require +#extension GL_KHR_shader_subgroup_basic : require +#extension GL_KHR_shader_subgroup_ballot : require +#extension GL_KHR_shader_subgroup_vote : require +#extension GL_ARB_shading_language_include : enable + +#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); + +uint baseID = gl_WorkGroupID.x * EXT_INVOCATIONS_PER_TASK; +uint laneID = gl_LocalInvocationID.x; + +layout(local_size_x = WORKGROUP_SIZE, local_size_y = 1, local_size_z = 1) in; + +#define BARRIER() \ +memoryBarrierShared(); \ +barrier(); + +layout(push_constant) uniform PUSH_CONSTANTS { + uint segment_size; + float time_step; + float inv_time_step; +}; + +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); + + 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 segment_index = baseID + local_index; + + if(segment_index >= segment_size) break; + + StrandSegment segment = strand_segments[segment_index]; + update_inertia_w(segment); + + //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)); + } + + //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; + } +} \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Graphics/Mesh/DynamicStrandsRendering.mesh b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Graphics/Mesh/DynamicStrandsRendering.mesh index 4347bdd..4ac2908 100644 --- a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Graphics/Mesh/DynamicStrandsRendering.mesh +++ b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Graphics/Mesh/DynamicStrandsRendering.mesh @@ -4,6 +4,8 @@ #extension GL_ARB_shading_language_include : enable #extension GL_EXT_shader_explicit_arithmetic_types_int8 : require +#include "Math.glsl" + #include "PerFrame.glsl" #define DYNAMIC_STRANDS_SET 1 @@ -76,66 +78,20 @@ taskPayloadSharedEXT Task ts_in; uint strand_segment_id = ts_in.baseID + ts_in.deltaIDs[gl_WorkGroupID.x]; uint laneID = gl_LocalInvocationID.x; -mat4 translate(in vec3 d){ - return mat4(1, 0, 0, 0, - 0, 1, 0, 0, - 0, 0, 1, 0, - d.x, d.y, d.z, 1); -} - -mat4 scale(in vec3 s){ - return mat4(s.x, 0, 0, 0, - 0, s.y, 0, 0, - 0, 0, s.z, 0, - 0, 0, 0, 1); -} - -mat4 mat4_cast(in vec4 q){ - float qxx = q.x * q.x; - float qyy = q.y * q.y; - float qzz = q.z * q.z; - float qxz = q.x * q.z; - float qxy = q.x * q.y; - float qyz = q.y * q.z; - float qwx = q.w * q.x; - float qwy = q.w * q.y; - float qwz = q.w * q.z; - mat4 ret_val = mat4(0, 0, 0, 0, - 0, 0, 0, 0, - 0, 0, 0, 0, - 0, 0, 0, 1); - ret_val[0][0] = 1.0 - 2.0 * (qyy + qzz); - ret_val[0][1] = 2.0 * (qxy + qwz); - ret_val[0][2] = 2.0 * (qxz - qwy); - - ret_val[1][0] = 2.0 * (qxy - qwz); - ret_val[1][1] = 1.0 - 2.0 * (qxx + qzz); - ret_val[1][2] = 2.0 * (qyz + qwx); - - ret_val[2][0] = 2.0 * (qxz + qwy); - ret_val[2][1] = 2.0 * (qyz - qwx); - ret_val[2][2] = 1.0 - 2.0 * (qxx + qyy); - - return ret_val; -} void main(){ + if(strand_segment_id >= strand_segments_size) return; SetMeshOutputsEXT(SEGMENT_VERTICES_SIZE, SEGMENT_TRIANGLE_SIZE); StrandSegment strand_segment = strand_segments[strand_segment_id]; - vec3 start_position = strand_segment_particles[strand_segment.start_particle_index].position_thickness.xyz; - float start_thickness = strand_segment_particles[strand_segment.start_particle_index].position_thickness.w; - vec4 start_color = strand_segment_particles[strand_segment.start_particle_index].color; - vec3 end_position = strand_segment_particles[strand_segment.end_particle_index].position_thickness.xyz; - float end_thickness = strand_segment_particles[strand_segment.end_particle_index].position_thickness.w; - vec4 end_color = strand_segment_particles[strand_segment.end_particle_index].color; - - vec3 center_position = (start_position + end_position) * 0.5; - float center_thickness = (start_thickness + end_thickness) * 0.5; - vec4 rotation = strand_segment.rotation; - mat4 instance_matrix = translate(end_position) * mat4_cast(rotation) * scale(vec3(center_thickness, center_thickness, distance(start_position, end_position))); + vec3 center_position = strand_segment.x.xyz; + float center_thickness = strand_segment.color_radius.w * 2.0; + + vec4 rotation = strand_segment.q; + + mat4 instance_matrix = translate(center_position) * mat4_cast(rotation) * scale(vec3(center_thickness, center_thickness, strand_segment.x0_length.w)); mat4 transform = EE_CAMERAS[camera_index].projection_view * instance_matrix; @@ -149,7 +105,7 @@ void main(){ ms_v_out[vertex_index].frag_pos = vec4(instance_matrix * vec4(vertex_position, 1.0)).xyz; ms_v_out[vertex_index].normal = vec3(0, 1, 0); ms_v_out[vertex_index].tangent = vec3(0, 0, 1); - ms_v_out[vertex_index].color = end_color; + ms_v_out[vertex_index].color = vec4(strand_segment.color_radius.xyz, 1.0); gl_MeshVerticesEXT[vertex_index].gl_Position = transform * vec4(vertex_position, 1.0); } diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Includes/DynamicStrands.glsl b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Includes/DynamicStrands.glsl index 4da2231..4aa0bcb 100644 --- a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Includes/DynamicStrands.glsl +++ b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Includes/DynamicStrands.glsl @@ -4,46 +4,43 @@ struct StrandSegment { int next_handle; int strand_handle; int index; + + vec4 color_radius; - vec4 rotation; + vec4 q0; + vec4 q; + vec4 last_q; + vec4 old_q; + vec4 torque; - int start_particle_index; - int end_particle_index; + vec4 x0_length; + vec4 x; + vec4 last_x; + vec4 old_x; + vec4 acceleration; - int neighbors[10]; -}; + vec4 inertia_tensor_mass; -struct Strand { - int first_strand_segment_handle; - int last_strand_segment_handle; -}; + vec4 inv_inertia_tensor_inv_mass; -struct StrandSegmentParticle { - vec4 position_thickness; - vec4 color; -}; + mat4 inertia_w; + mat4 inv_inertia_w; -layout(std430, set = DYNAMIC_STRANDS_SET, binding = 0) readonly buffer REF_STRAND_SEGMENTS_BLOCK { - StrandSegment ref_strand_segments[]; + int neighbors[8]; }; -layout(std430, set = DYNAMIC_STRANDS_SET, binding = 1) readonly buffer REF_STRANDS_BLOCK { - Strand ref_strands[]; -}; +struct Strand { + int begin_segment_handle; + int end_segment_handle; -layout(std430, set = DYNAMIC_STRANDS_SET, binding = 2) readonly buffer REF_STRAND_SEGMENT_PARTICLES_BLOCK { - StrandSegmentParticle ref_strand_segment_particles[]; + int segment_size; + int padding0; }; -layout(std430, set = DYNAMIC_STRANDS_SET, binding = 3) buffer STRAND_SEGMENTS_BLOCK { +layout(std430, set = DYNAMIC_STRANDS_SET, binding = 0) buffer STRAND_SEGMENTS_BLOCK { StrandSegment strand_segments[]; }; -layout(std430, set = DYNAMIC_STRANDS_SET, binding = 4) buffer STRANDS_BLOCK { +layout(std430, set = DYNAMIC_STRANDS_SET, binding = 1) buffer STRANDS_BLOCK { Strand strands[]; }; - -layout(std430, set = DYNAMIC_STRANDS_SET, binding = 5) buffer STRAND_SEGMENT_PARTICLES_BLOCK { - StrandSegmentParticle strand_segment_particles[]; -}; - diff --git a/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Includes/DynamicStrandsPhysics.glsl b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Includes/DynamicStrandsPhysics.glsl new file mode 100644 index 0000000..731afd8 --- /dev/null +++ b/EvoEngine_Plugins/EcoSysLab/Internals/EcoSysLabResources/Shaders/Includes/DynamicStrandsPhysics.glsl @@ -0,0 +1,432 @@ + +#include "Math.glsl" +struct RodProperties { + float density; + float i1; + float j; + float e; + + float g; + float padding0; + float padding1; + float padding2; +}; + +layout(set = DYNAMIC_STRANDS_PHYSICS_SET, binding = 0) readonly uniform ROD_PROPERTIES_BLOCK { + RodProperties rod_properties; +}; + +struct PerStrandData { + int begin_rod_constraint_handle; + int end_rod_constraint_handle; + int segment_size; + int padding1; +}; + +layout(std430, set = DYNAMIC_STRANDS_PHYSICS_SET, binding = 1) readonly buffer PER_STRAND_DATA_LIST_BLOCK { + PerStrandData per_strand_data_list[]; +}; + +struct RodConstraint { + vec4 stiffness_coefficient_k_segment0_index; + vec4 rest_darboux_vector_segment1_index; + vec4 stretch_compliance_average_segment_length; + vec4 bending_and_torsion_compliance_next_constraint_handle; + /** + * \brief + * constraintInfo contains + * 0: connector in segment 0 (local) + * 1: connector in segment 1 (local) + * 2: connector in segment 0 (global) + * 3: connector in segment 1 (global) + */ + mat4 constraint_info; + + float lambda_sum0; + float lambda_sum1; + float lambda_sum2; + float lambda_sum3; + + float lambda_sum4; + float lambda_sum5; + float padding0; + float padding1; +}; + +layout(std430, set = DYNAMIC_STRANDS_PHYSICS_SET, binding = 2) buffer ROD_CONSTRAINT_BLOCK { + RodConstraint rod_constraints[]; +}; + +struct PerSegmentData { + vec4 lambda_ss; +}; + +layout(std430, set = DYNAMIC_STRANDS_PHYSICS_SET, binding = 3) buffer PER_SEGMENT_DATA_LIST_BLOCK { + PerSegmentData per_segment_data_list[]; +}; + +struct vec6 { + float v[6]; +}; + +struct mat6 { + vec6 v[6]; +}; + +void compute_bending_and_torsion_jacobians(in vec4 q0, in vec4 q1, float segment_length, out mat4 j_omega_0, + out mat4 j_omega_1); +void compute_matrix_g(in vec4 q, out mat4 g); +void compute_matrix_k(in vec3 connector, in float inv_mass, in vec3 x, in mat3 inverse_inertia, out mat3 k); +mat3 cross_product_matrix(in vec3 v); +void solve_ldlt(in mat6 a, out vec6 x, in vec6 b); + +mat3 compute_inverse_inertia_tensor_w(in StrandSegment strand_segment); + +void init_stretch_bending_twisting_constraint(in vec3 stiffness_coefficient_k, in float inverse_time_step_size, + in float average_segment_length, out vec3 stretch_compliance, + out vec3 bending_and_torsion_compliance, out vec6 lambda_sum); + +void update_stretch_bending_twisting_constraint(in vec3 x0, in vec4 q0, in vec3 x1, in vec4 q1, + inout mat4 constraint_info); + +void solve_stretch_bending_twisting_constraints(in float inv_mass0, in vec3 x0, in mat3 inverse_inertia0, in vec4 q0, + in float inv_mass1, in vec3 x1, in mat3 inverse_inertia1, in vec4 q1, + in vec3 rest_darboux_vector, in float average_segment_length, + in vec3 stretch_compliance, in vec3 bending_and_torsion_compliance, + in mat4 constraint_info, out vec3 x0_correction, out vec4 q0_correction, + out vec3 x1_correction, out vec4 q1_correction, inout vec6 lambda_sum); + +//------------------- +//| Implementations | +//------------------- + +void compute_bending_and_torsion_jacobians(in vec4 q0, in vec4 q1, float segment_length, out mat4 j_omega_0, + out mat4 j_omega_1) { + j_omega_0 = mat4(-q1[3], -q1[2], q1[1], q1[0], q1[2], -q1[3], -q1[0], q1[1], -q1[1], q1[0], -q1[3], q1[2], 0.0, 0.0, + 0.0, 0.0); + j_omega_1 = mat4(q0[3], q0[2], -q0[1], -q0[0], -q0[2], q0[3], q0[0], -q0[1], q0[1], -q0[0], q0[3], -q0[2], 0.0, 0.0, + 0.0, 0.0); + j_omega_0 = j_omega_0 * 2.0 / segment_length; + j_omega_1 = j_omega_1 * 2.0 / segment_length; +} + +void compute_matrix_g(in vec4 q, out mat4 g) { + g = mat4(0.5 * q[3], 0.5 * q[2], -0.5 * q[1], 0.0, -0.5 * q[2], 0.5 * q[3], 0.5 * q[0], 0.0, 0.5 * q[1], -0.5 * q[0], + 0.5 * q[3], 0.0, -0.5 * q[0], -0.5 * q[1], -0.5 * q[2], 0.0); +} + +void compute_matrix_k(in vec3 connector, in float inv_mass, in vec3 x, in mat3 inverse_inertia, out mat3 k) { + if (inv_mass != 0.0) { + vec3 v = connector - x; + float a = v[0]; + float b = v[1]; + float c = v[2]; + + float j11 = inverse_inertia[0][0]; + float j12 = inverse_inertia[0][1]; + float j13 = inverse_inertia[0][2]; + float j22 = inverse_inertia[1][1]; + float j23 = inverse_inertia[1][2]; + float j33 = inverse_inertia[2][2]; + + k[0][0] = c * c * j22 - b * c * (j23 + j23) + b * b * j33 + inv_mass; + k[0][1] = -(c * c * j12) + a * c * j23 + b * c * j13 - a * b * j33; + k[0][2] = b * c * j12 - a * c * j22 - b * b * j13 + a * b * j23; + k[1][0] = k[0][1]; + k[1][1] = c * c * j11 - a * c * (j13 + j13) + a * a * j33 + inv_mass; + k[1][2] = -(b * c * j11) + a * c * j12 + a * b * j13 - a * a * j23; + k[2][0] = k[0][2]; + k[2][1] = k[1][2]; + k[2][2] = b * b * j11 - a * b * (j12 + j12) + a * a * j22 + inv_mass; + } else { + k = mat3(0, 0, 0, 0, 0, 0, 0, 0, 0); + } +} + +mat3 cross_product_matrix(in vec3 v) { + return mat3(0, -v[2], v[1], v[2], 0, -v[0], -v[1], v[0], 0); +} + +void solve_ldlt(in mat6 A, out vec6 x, in vec6 b) { + mat6 L; + L.v[0].v[0] = 1.0; + L.v[0].v[1] = 1.0; + L.v[0].v[2] = 1.0; + L.v[0].v[3] = 1.0; + L.v[0].v[4] = 1.0; + L.v[0].v[5] = 1.0; + + L.v[1].v[0] = 1.0; + L.v[1].v[1] = 1.0; + L.v[1].v[2] = 1.0; + L.v[1].v[3] = 1.0; + L.v[1].v[4] = 1.0; + L.v[1].v[5] = 1.0; + + L.v[2].v[0] = 1.0; + L.v[2].v[1] = 1.0; + L.v[2].v[2] = 1.0; + L.v[2].v[3] = 1.0; + L.v[2].v[4] = 1.0; + L.v[2].v[5] = 1.0; + + L.v[3].v[0] = 1.0; + L.v[3].v[1] = 1.0; + L.v[3].v[2] = 1.0; + L.v[3].v[3] = 1.0; + L.v[3].v[4] = 1.0; + L.v[3].v[5] = 1.0; + + L.v[4].v[0] = 1.0; + L.v[4].v[1] = 1.0; + L.v[4].v[2] = 1.0; + L.v[4].v[3] = 1.0; + L.v[4].v[4] = 1.0; + L.v[4].v[5] = 1.0; + + L.v[5].v[0] = 1.0; + L.v[5].v[1] = 1.0; + L.v[5].v[2] = 1.0; + L.v[5].v[3] = 1.0; + L.v[5].v[4] = 1.0; + L.v[5].v[5] = 1.0; + + vec6 D; + D.v[0] = 0.0; + D.v[1] = 0.0; + D.v[2] = 0.0; + D.v[3] = 0.0; + D.v[4] = 0.0; + D.v[5] = 0.0; + + // Decomposition + for (int j = 0; j < 6; j++) { + // Compute diagonal elements of D + D.v[j] = A.v[j].v[j]; + for (int k = 0; k < j; k++) { + D.v[j] -= L.v[j].v[k] * L.v[j].v[k] * D.v[k]; + } + // Compute elements of L (below the diagonal) + for (int i = j + 1; i < 6; i++) { + L.v[i].v[j] = A.v[i].v[j]; + for (int k = 0; k < j; k++) { + L.v[i].v[j] -= L.v[i].v[k] * L.v[j].v[k] * D.v[k]; + } + L.v[i].v[j] /= D.v[j]; + } + } + + // Forward substitution + vec6 y; + for (int i = 0; i < 6; i++) { + y.v[i] = b.v[i]; + for (int j = 0; j < i; j++) { + y.v[i] -= L.v[i].v[j] * y.v[j]; + } + } + + // Diagonal solve + vec6 z; + for (int i = 0; i < 6; i++) { + z.v[i] = y.v[i] / D.v[i]; + } + + // Backward substitution + for (int i = 5; i >= 0; i--) { + x.v[i] = z.v[i]; + for (int j = i + 1; j < 6; j++) { + x.v[i] -= L.v[j].v[i] * x.v[j]; + } + } +} + +mat3 compute_inverse_inertia_tensor_w(in StrandSegment strand_segment) { + mat3 rot = mat3_cast(strand_segment.q); + mat3 diag = + mat3( + strand_segment.inv_inertia_tensor_inv_mass.x, 0.0, 0.0, + 0.0, strand_segment.inv_inertia_tensor_inv_mass.y, 0.0, + 0.0, 0.0, strand_segment.inv_inertia_tensor_inv_mass.z); + return rot * diag * transpose(rot); +} + +void init_stretch_bending_twisting_constraint(in vec3 stiffness_coefficient_k, in float inverse_time_step_size, + in float average_segment_length, out vec3 stretch_compliance, + out vec3 bending_and_torsion_compliance, out vec6 lambda_sum) { + float inverse_time_step_quadratic = inverse_time_step_size * inverse_time_step_size; + float stretch_regularization_parameter = 0.1; + + stretch_compliance = vec3(stretch_regularization_parameter * inverse_time_step_quadratic, + stretch_regularization_parameter * inverse_time_step_quadratic, + stretch_regularization_parameter * inverse_time_step_quadratic); + + bending_and_torsion_compliance = vec3(inverse_time_step_quadratic / stiffness_coefficient_k[0], + inverse_time_step_quadratic / stiffness_coefficient_k[1], + inverse_time_step_quadratic / stiffness_coefficient_k[2]); + + lambda_sum.v[0] = 0.0; + lambda_sum.v[1] = 0.0; + lambda_sum.v[2] = 0.0; + lambda_sum.v[3] = 0.0; + lambda_sum.v[4] = 0.0; + lambda_sum.v[5] = 0.0; +} + +void update_stretch_bending_twisting_constraint(in vec3 x0, in vec4 q0, in vec3 x1, in vec4 q1, + inout mat4 constraint_info) { + // constraintInfo contains + // 0: connector in segment 0 (local) + // 1: connector in segment 1 (local) + // 2: connector in segment 0 (global) + // 3: connector in segment 1 (global) + + // compute world space positions of connectors + mat3 rot0 = mat3_cast(q0); + mat3 rot1 = mat3_cast(q1); + + constraint_info[2].xyz = rot0 * constraint_info[0].xyz + x0; + constraint_info[3].xyz = rot1 * constraint_info[1].xyz + x1; +} + +void solve_stretch_bending_twisting_constraints(in float inv_mass0, in vec3 x0, in mat3 inverse_inertia0, in vec4 q0, + in float inv_mass1, in vec3 x1, in mat3 inverse_inertia1, in vec4 q1, + in vec3 rest_darboux_vector, in float average_segment_length, + in vec3 stretch_compliance, in vec3 bending_and_torsion_compliance, + in mat4 constraint_info, out vec3 x0_correction, out vec4 q0_correction, + out vec3 x1_correction, out vec4 q1_correction, inout vec6 lambda_sum) { + vec3 omega = compute_darboux_vector(q0, q1, average_segment_length); + mat4 j_omega0, j_omega1; + compute_bending_and_torsion_jacobians(q0, q1, average_segment_length, j_omega0, j_omega1); + mat4 g0, g1; + compute_matrix_g(q0, g0); + compute_matrix_g(q1, g1); + + mat3 j_omega_g0 = mat3(j_omega0 * g0); + mat3 j_omega_g1 = mat3(j_omega1 * g1); + + vec3 connector0 = constraint_info[2].xyz; + vec3 connector1 = constraint_info[3].xyz; + + vec3 stretch_violation = connector0 - connector1; + vec3 bending_and_torsion_violation = omega - rest_darboux_vector; + + stretch_violation = stretch_violation + stretch_compliance * vec3(lambda_sum.v[0], lambda_sum.v[1], lambda_sum.v[2]); + bending_and_torsion_violation = + bending_and_torsion_violation + + bending_and_torsion_compliance * vec3(lambda_sum.v[3], lambda_sum.v[4], lambda_sum.v[5]); + + vec6 rhs; + rhs.v[0] = -stretch_violation.x; + rhs.v[1] = -stretch_violation.y; + rhs.v[2] = -stretch_violation.z; + rhs.v[3] = -bending_and_torsion_violation.x; + rhs.v[4] = -bending_and_torsion_violation.y; + rhs.v[5] = -bending_and_torsion_violation.z; + + mat3 k0, k1; + compute_matrix_k(connector0, inv_mass0, x0, inverse_inertia0, k0); + compute_matrix_k(connector1, inv_mass1, x1, inverse_inertia1, k1); + mat3 k = k0 + k1; + mat6 jmjt; + jmjt.v[0].v[0] = k[0][0]; + jmjt.v[0].v[1] = k[0][1]; + jmjt.v[0].v[2] = k[0][2]; + jmjt.v[1].v[0] = k[1][0]; + jmjt.v[1].v[1] = k[1][1]; + jmjt.v[1].v[2] = k[1][2]; + jmjt.v[2].v[0] = k[2][0]; + jmjt.v[2].v[1] = k[2][1]; + jmjt.v[2].v[2] = k[2][2]; + + vec3 r0 = connector0 - x0; + vec3 r1 = connector1 - x1; + + mat3 r0_cross_t = cross_product_matrix(-r0); + mat3 r1_cross_t = cross_product_matrix(-r1); + + mat3 off_diag = mat3(0, 0, 0, 0, 0, 0, 0, 0, 0); + if (inv_mass0 != 0.0) { + off_diag = j_omega_g0 * inverse_inertia0 * r0_cross_t * -1.0; + } + if (inv_mass1 != 0.0) { + off_diag = off_diag + j_omega_g1 * inverse_inertia1 * r1_cross_t * -1.0; + } + jmjt.v[3].v[0] = off_diag[0][0]; + jmjt.v[3].v[1] = off_diag[0][1]; + jmjt.v[3].v[2] = off_diag[0][2]; + jmjt.v[4].v[0] = off_diag[1][0]; + jmjt.v[4].v[1] = off_diag[1][1]; + jmjt.v[4].v[2] = off_diag[1][2]; + jmjt.v[5].v[0] = off_diag[2][0]; + jmjt.v[5].v[1] = off_diag[2][1]; + jmjt.v[5].v[2] = off_diag[2][2]; + + mat3 off_diag_t = transpose(off_diag); + jmjt.v[0].v[3] = off_diag_t[0][0]; + jmjt.v[0].v[4] = off_diag_t[0][1]; + jmjt.v[0].v[5] = off_diag_t[0][2]; + jmjt.v[1].v[3] = off_diag_t[1][0]; + jmjt.v[1].v[4] = off_diag_t[1][1]; + jmjt.v[1].v[5] = off_diag_t[1][2]; + jmjt.v[2].v[3] = off_diag_t[2][0]; + jmjt.v[2].v[4] = off_diag_t[2][1]; + jmjt.v[2].v[5] = off_diag_t[2][2]; + + mat3 m_inv_jt0 = inverse_inertia0 * transpose(j_omega_g0); + mat3 m_inv_jt1 = inverse_inertia1 * transpose(j_omega_g1); + + mat3 jmjt_omega; + if (inv_mass0 != 0.0) { + jmjt_omega = j_omega_g0 * m_inv_jt0; + } + if (inv_mass1 != 0.0) { + jmjt_omega = jmjt_omega + j_omega_g1 * m_inv_jt1; + } + + jmjt.v[3].v[3] = jmjt_omega[0][0]; + jmjt.v[3].v[4] = jmjt_omega[0][1]; + jmjt.v[3].v[5] = jmjt_omega[0][2]; + jmjt.v[4].v[3] = jmjt_omega[1][0]; + jmjt.v[4].v[4] = jmjt_omega[1][1]; + jmjt.v[4].v[5] = jmjt_omega[1][2]; + jmjt.v[5].v[3] = jmjt_omega[2][0]; + jmjt.v[5].v[4] = jmjt_omega[2][1]; + jmjt.v[5].v[5] = jmjt_omega[2][2]; + + vec6 delta_lambda; + solve_ldlt(jmjt, delta_lambda, rhs); + + lambda_sum.v[0] += delta_lambda.v[0]; + lambda_sum.v[1] += delta_lambda.v[1]; + lambda_sum.v[2] += delta_lambda.v[2]; + lambda_sum.v[3] += delta_lambda.v[3]; + lambda_sum.v[4] += delta_lambda.v[4]; + lambda_sum.v[5] += delta_lambda.v[5]; + + vec3 delta_lambda_stretch = vec3(delta_lambda.v[0], delta_lambda.v[1], delta_lambda.v[2]); + vec3 delta_lambda_bending_and_torsion = vec3(delta_lambda.v[3], delta_lambda.v[4], delta_lambda.v[5]); + + x0_correction = vec3(0, 0, 0); + q0_correction = vec4(0, 0, 0, 0); + x1_correction = vec3(0, 0, 0); + q1_correction = vec4(0, 0, 0, 0); + + if (inv_mass0 != 0.0) { + x0_correction = inv_mass0 * delta_lambda_stretch; + vec3 v = + inverse_inertia0 * r0_cross_t * (-1.0 * delta_lambda_stretch) + m_inv_jt0 * delta_lambda_bending_and_torsion; + q0_correction[0] = g0[0][0] * v[0] + g0[0][1] * v[1] + g0[0][2] * v[2]; + q0_correction[1] = g0[1][0] * v[0] + g0[1][1] * v[1] + g0[1][2] * v[2]; + q0_correction[2] = g0[2][0] * v[0] + g0[2][1] * v[1] + g0[2][2] * v[2]; + q0_correction[3] = g0[3][0] * v[0] + g0[3][1] * v[1] + g0[3][2] * v[2]; + } + + if (inv_mass1 != 0.0) { + x1_correction = -inv_mass1 * delta_lambda_stretch; + vec3 v = inverse_inertia1 * r1_cross_t * delta_lambda_stretch + m_inv_jt1 * delta_lambda_bending_and_torsion; + q1_correction[0] = g1[0][0] * v[0] + g1[0][1] * v[1] + g1[0][2] * v[2]; + q1_correction[1] = g1[1][0] * v[0] + g1[1][1] * v[1] + g1[1][2] * v[2]; + q1_correction[2] = g1[2][0] * v[0] + g1[2][1] * v[1] + g1[2][2] * v[2]; + q1_correction[3] = g1[3][0] * v[0] + g1[3][1] * v[1] + g1[3][2] * v[2]; + } +} \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/include/StrandModel/BrokenBranch.hpp b/EvoEngine_Plugins/EcoSysLab/include/StrandModel/BrokenBranch.hpp index ca33652..0df2758 100644 --- a/EvoEngine_Plugins/EcoSysLab/include/StrandModel/BrokenBranch.hpp +++ b/EvoEngine_Plugins/EcoSysLab/include/StrandModel/BrokenBranch.hpp @@ -9,22 +9,24 @@ class BrokenBranch : public IPrivateComponent { bool enable_simulation = true; - PrivateComponentRef tree_ref; - StrandModelStrandGroup strand_group; - StrandModelStrandGroup subdivided_strand_group; + PrivateComponentRef tree_ref{}; + StrandModelStrandGroup strand_group{}; + StrandModelStrandGroup subdivided_strand_group{}; - DynamicStrands dynamic_strands; - void UploadStrands(); + DynamicStrands::InitializeParameters initialize_parameters{}; + DynamicStrands::StepParameters step_parameters{}; + DynamicStrands dynamic_strands{}; + void UpdateDynamicStrands(); void Serialize(YAML::Emitter& out) const override; void Deserialize(const YAML::Node& in) override; bool OnInspect(const std::shared_ptr& editor_layer) override; - void Update() override; void LateUpdate() override; + void FixedUpdate() override; void OnCreate() override; void OnDestroy() override; void CollectAssetRef(std::vector& list) override; void ExperimentSetup(); - static void Subdivide(float segment_length, const StrandModelStrandGroup& src, StrandModelStrandGroup& dst); + void Subdivide(float segment_length, const StrandModelStrandGroup& src); void InitializeStrandParticles(const StrandModelStrandGroup& strand_group) const; void ClearStrandParticles() const; diff --git a/EvoEngine_Plugins/EcoSysLab/include/StrandModel/DynamicStrands.hpp b/EvoEngine_Plugins/EcoSysLab/include/StrandModel/DynamicStrands.hpp index e4a702e..78c29ed 100644 --- a/EvoEngine_Plugins/EcoSysLab/include/StrandModel/DynamicStrands.hpp +++ b/EvoEngine_Plugins/EcoSysLab/include/StrandModel/DynamicStrands.hpp @@ -4,72 +4,195 @@ using namespace evo_engine; namespace eco_sys_lab_plugin { +class DynamicStrandsPreStep; +class DynamicStrandsPostStep; +class IDynamicStrandsOperator; +class IDynamicStrandsConstraint; class DynamicStrands { - void Upload() const; - bool current_left = true; - uint32_t current_version = 0; - uint32_t device_buffer_version[2]; - std::shared_ptr device_ref_strand_segments_buffer[2]; - std::shared_ptr device_ref_strands_buffer[2]; - std::shared_ptr device_ref_strand_segment_particles_buffer[2]; - - std::shared_ptr device_strand_segments_buffer[2]; - std::shared_ptr device_strands_buffer[2]; - std::shared_ptr device_strand_segment_particles_buffer[2]; - public: + DynamicStrands(); +#pragma region Initialization struct InitializeParameters { + bool static_root = true; + float wood_density = 1.f; + float youngs_modulus = 1.f; + float torsion_modulus = 1.f; GlobalTransform root_transform{}; }; + template + void InitializeStrandsGroup(const InitializeParameters& initialize_parameters, + const StrandGroup& strand_group); +#pragma endregion +#pragma region Step + struct PhysicsParameters { + float time_step = 0.01f; + uint32_t sub_step = 1; + uint32_t max_iteration; + }; + + struct RenderParameters { + std::shared_ptr target_camera{}; + }; struct StepParameters { - bool physics = true; - struct PhysicsParameters { - } physics_parameters{}; + bool physics = false; + PhysicsParameters physics_parameters{}; + + bool render = false; + RenderParameters render_parameters; + }; - bool render = true; - struct RenderParameters { - std::shared_ptr target_camera{}; - } render_parameters; - } step_parameters{}; + + std::vector> operators; + std::shared_ptr pre_step; + std::vector> constraints; + std::shared_ptr post_step; + void Step(const StepParameters& target_step_parameters); +#pragma endregion +#pragma region Shared Data struct GpuStrand { - StrandSegmentHandle first_strand_segment_handle = -1; - StrandSegmentHandle last_strand_segment_handle = -1; + int begin_segment_handle = -1; + int end_segment_handle = -1; + + int segment_size; + int padding0; }; + struct GpuStrandSegment { int prev_handle = -1; int next_handle = -1; int strand_handle = -1; int index = -1; - glm::quat rotation; + glm::vec3 color; + float radius; - int start_particle_index; - int end_particle_index; - int neighbors[10]; - }; - struct GpuStrandSegmentParticle { - glm::vec3 position; - float thickness; - glm::vec4 color; + // Initial rotation + glm::quat q0; + // Current rotation + glm::quat q; + // 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; + float length; + // Current position + glm::vec4 x; + // Last frame position + 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]; }; - std::vector ref_strand_segments; - std::vector ref_strand_segment_particles; - std::vector ref_strands; + inline static std::shared_ptr strands_layout{}; - DynamicStrands(); - void UpdateData(const InitializeParameters& initialize_parameters, const std::vector& target_strands, - const std::vector& target_strand_segments); - void Step(const StepParameters& target_step_parameters); - void Step(); + + std::shared_ptr device_strand_segments_buffer; + std::shared_ptr device_strands_buffer; + + std::vector strands; + std::vector strand_segments; +#pragma endregion + + void Upload() const; + void Download(); + void Clear(); + std::vector> strands_descriptor_sets; + private: - void BindStrandsDescriptorSet(const std::shared_ptr& target_descriptor_set) const; - void Render(const StepParameters::RenderParameters& render_parameters) const; - void Physics(const StepParameters::PhysicsParameters& physics_parameters) const; + 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; + void Physics(const PhysicsParameters& physics_parameters, + const std::vector>& target_operators, + const std::shared_ptr& target_pre_step, + const std::vector>& target_constraints, + const std::shared_ptr& target_post_step); }; + +template +void DynamicStrands::InitializeStrandsGroup(const InitializeParameters& initialize_parameters, + const StrandGroup& strand_group) { + Clear(); + assert(initialize_parameters.root_transform.GetScale() == glm::vec3(1.0f)); + const auto& target_strands = strand_group.PeekStrands(); + const auto& target_strand_segments = strand_group.PeekStrandSegments(); + strands.resize(target_strands.size()); + + Jobs::RunParallelFor(target_strands.size(), [&](const size_t i) { + auto& strand = strands[i]; + const auto& target_strand = target_strands[i]; + const auto& handles = target_strand.PeekStrandSegmentHandles(); + strand.segment_size = handles.size(); + if (!handles.empty()) { + strand.begin_segment_handle = handles.front(); + strand.end_segment_handle = handles.back(); + }else { + strand.begin_segment_handle = -1; + strand.end_segment_handle = -1; + } + }); + + strand_segments.resize(target_strand_segments.size()); + Jobs::RunParallelFor(target_strand_segments.size(), [&](const size_t i) { + auto& segment = strand_segments[i]; + const auto& target_strand_segment = target_strand_segments[i]; + segment.prev_handle = target_strand_segment.GetPrevHandle(); + segment.next_handle = target_strand_segment.GetNextHandle(); + segment.strand_handle = target_strand_segment.GetStrandHandle(); + segment.index = target_strand_segment.GetIndex(); + + segment.color = target_strand_segment.end_color; + segment.radius = target_strand_segment.end_thickness * .5f; + + 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(i))), 0.0); + + segment.x0 = segment.x; + segment.acceleration = glm::vec3(0.f); + + 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() * initialize_parameters.wood_density; + segment.inv_mass = 1.f / segment.mass; + segment.inertia_tensor = ComputeInertiaTensorRod(segment.mass, segment.radius, segment.length); + segment.inv_inertia_tensor = 1.f / segment.inertia_tensor; + } + + segment.length = strand_group.GetStrandSegmentLength(static_cast(i)); + }); + Upload(); + + InitializeOperators(initialize_parameters); + InitializeConstraints(initialize_parameters); +} } // namespace eco_sys_lab_plugin \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/include/StrandModel/DynamicStrandsPhysics.hpp b/EvoEngine_Plugins/EcoSysLab/include/StrandModel/DynamicStrandsPhysics.hpp new file mode 100644 index 0000000..940fed9 --- /dev/null +++ b/EvoEngine_Plugins/EcoSysLab/include/StrandModel/DynamicStrandsPhysics.hpp @@ -0,0 +1,197 @@ +#pragma once +#include "DynamicStrands.hpp" +using namespace evo_engine; + +namespace eco_sys_lab_plugin { +class DynamicStrandsPreStep { + public: + DynamicStrandsPreStep(); + + struct PreStepPushConstant { + uint32_t segment_size = 0; + float time_step = 0.01f; + float inv_time_step = 100.f; + }; + + inline static std::shared_ptr pre_step_pipeline; + void Execute(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands); +}; + +class DynamicStrandsPostStep { + public: + DynamicStrandsPostStep(); + + struct PostStepPushConstant { + uint32_t segment_size = 0; + }; + + inline static std::shared_ptr post_step_pipeline; + void Execute(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands); +}; +class IDynamicStrandsOperator { + public: + virtual ~IDynamicStrandsOperator() = default; + virtual void InitializeData(const DynamicStrands::InitializeParameters& initialize_parameters, + const DynamicStrands& target_dynamic_strands) { + } + virtual void Execute(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) = 0; + virtual void DownloadData() { + } +}; +class IDynamicStrandsConstraint { + public: + virtual void InitializeData(const DynamicStrands::InitializeParameters& initialize_parameters, + const DynamicStrands& target_dynamic_strands){}; + virtual void Project(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) = 0; + + virtual void DownloadData() { + } +}; + +#pragma region Operators +class DsPositionUpdate final : public IDynamicStrandsOperator { + public: + DsPositionUpdate(); + + struct PositionUpdate { + glm::vec3 new_position = glm::vec3(0.f); + uint32_t strand_segment_index = 0; + }; + + inline static std::shared_ptr layout{}; + + std::vector commands; + std::vector> commands_buffer; + + struct PositionUpdatePushConstant { + uint32_t commands_size = 0; + }; + + inline static std::shared_ptr external_force_pipeline; + std::vector> commands_descriptor_sets; + + void Execute(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) override; +}; +class DsExternalForce final : public IDynamicStrandsOperator { + public: + DsExternalForce(); + + struct ExternalForce { + glm::vec3 force = glm::vec3(0.f); + uint32_t strand_segment_index = 0; + }; + + inline static std::shared_ptr layout{}; + + std::vector commands; + std::vector> commands_buffer; + + struct ExternalForcePushConstant { + uint32_t commands_size = 0; + }; + + inline static std::shared_ptr external_force_pipeline{}; + std::vector> commands_descriptor_sets{}; + + void Execute(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) override; +}; + +#pragma endregion + +#pragma region Constraints +class DsStiffRod final : public IDynamicStrandsConstraint { + public: + DsStiffRod(); + inline static std::shared_ptr layout{}; + + struct RodProperties { + float density; + float i1; + float j; + float e; + + float g; + float padding0; + float padding1; + float padding2; + }; + + struct GpuRodConstraint { + glm::vec3 stiffness_coefficient_k; + int segment0_index; + + glm::vec3 rest_darboux_vector; + int segment1_index; + + glm::vec3 stretch_compliance; + float average_segment_length; + + glm::vec3 bending_and_torsion_compliance; + int next_handle = -1; + + /** + * \brief + * constraintInfo contains + * 0: connector in segment 0 (local) + * 1: connector in segment 1 (local) + * 2: connector in segment 0 (global) + * 3: connector in segment 1 (global) + */ + glm::mat4 constraint_info; + + float lambda_sum[8]; + }; + + struct PerStrandData { + int begin_rod_constraint_handle = -1; + int end_rod_constraint_handle = -1; + int padding0 = -1; + int padding1 = -1; + }; + + RodProperties rod_properties; + + std::vector per_strand_data_list; + std::vector rod_constraints; + + std::vector> rod_properties_buffer; + + std::shared_ptr per_strand_data_list_buffer; + std::shared_ptr rod_constraints_buffer; + + struct StiffRodInitConstraintConstant { + uint32_t constraint_size = 0; + float time_step; + float inv_time_step; + }; + + struct StiffRodUpdateConstraintConstant { + uint32_t constraint_size = 0; + }; + + struct StiffRodProjectConstraintConstant { + uint32_t strand_size = 0; + }; + + inline static std::shared_ptr init_constraint_pipeline{}; + inline static std::shared_ptr update_constraint_pipeline{}; + inline static std::shared_ptr project_constraint_pipeline{}; + + std::vector> strands_physics_descriptor_sets{}; + + void InitializeData(const DynamicStrands::InitializeParameters& initialize_parameters, + const DynamicStrands& target_dynamic_strands) override; + + void Project(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) override; + void DownloadData() override; + static glm::vec3 ComputeDarbouxVector(const glm::quat& q0, const glm::quat& q1, float average_segment_length); +}; +#pragma endregion +} // namespace eco_sys_lab_plugin \ No newline at end of file diff --git a/EvoEngine_Plugins/EcoSysLab/include/Structures/StrandGroup.hpp b/EvoEngine_Plugins/EcoSysLab/include/Structures/StrandGroup.hpp index f5ad490..4af7f82 100644 --- a/EvoEngine_Plugins/EcoSysLab/include/Structures/StrandGroup.hpp +++ b/EvoEngine_Plugins/EcoSysLab/include/Structures/StrandGroup.hpp @@ -440,16 +440,15 @@ void StrandGroup::CalculateRotat for (const auto& strand : strands_) { if (strand.strand_segment_handles_.empty()) continue; glm::vec3 prev_position = strand.start_position; - const glm::vec3 initial_front = - glm::normalize(strand_segments_[strand.strand_segment_handles_[0]].end_position - prev_position); - auto prev_up = glm::vec3(initial_front.y, initial_front.z, initial_front.x); - strand_segments_[strand.strand_segment_handles_[0]].rotation = glm::quatLookAt(initial_front, prev_up); + const auto& end_position = strand_segments_[strand.strand_segment_handles_[0]].end_position; + glm::vec3 front = glm::normalize(end_position - prev_position); + strand_segments_[strand.strand_segment_handles_[0]].rotation = glm::quatLookAt(front, glm::vec3(front.y, front.z, front.x)); + + prev_position = end_position; for (uint32_t i = 0; i < strand.strand_segment_handles_.size() - 1; i++) { auto& segment = strand_segments_[strand.strand_segment_handles_[i + 1]]; - - auto front = glm::normalize(segment.end_position - prev_position); - prev_up = glm::normalize(glm::cross(glm::cross(front, prev_up), front)); - segment.rotation = glm::quatLookAt(front, prev_up); + front = glm::normalize(segment.end_position - prev_position); + segment.rotation = glm::quatLookAt(front, glm::vec3(front.y, front.z, front.x)); prev_position = segment.end_position; } } diff --git a/EvoEngine_Plugins/EcoSysLab/src/BrokenBranch.cpp b/EvoEngine_Plugins/EcoSysLab/src/BrokenBranch.cpp index cf26759..ae52732 100644 --- a/EvoEngine_Plugins/EcoSysLab/src/BrokenBranch.cpp +++ b/EvoEngine_Plugins/EcoSysLab/src/BrokenBranch.cpp @@ -1,17 +1,35 @@ #include "BrokenBranch.hpp" - #include "Delaunay.hpp" +#include "DynamicStrandsPhysics.hpp" #include "Tree.hpp" using namespace eco_sys_lab_plugin; -void BrokenBranch::UploadStrands() { - DynamicStrands::InitializeParameters initialize_parameters{}; +void BrokenBranch::UpdateDynamicStrands() { + if (step_parameters.physics || step_parameters.render) { + EVOENGINE_WARNING("Auto turned physics and render off."); + } + step_parameters.physics = false; + step_parameters.render = false; const auto owner = GetOwner(); const auto scene = GetScene(); initialize_parameters.root_transform = scene->GetDataComponent(owner); - dynamic_strands.UpdateData(initialize_parameters, subdivided_strand_group.PeekStrands(), - subdivided_strand_group.PeekStrandSegments()); + + + dynamic_strands.InitializeStrandsGroup(initialize_parameters, subdivided_strand_group); + const auto& external_force = std::dynamic_pointer_cast(dynamic_strands.operators.front()); + external_force->commands.resize(dynamic_strands.strand_segments.size()); + Jobs::RunParallelFor(external_force->commands.size(), [&](const size_t i) { + external_force->commands[i].strand_segment_index = static_cast(i); + external_force->commands[i].force = glm::vec3(0, -1, 0) * dynamic_strands.strand_segments[i].mass; + }); + const auto& position_update = std::dynamic_pointer_cast(dynamic_strands.operators.back()); + position_update->commands.resize(dynamic_strands.strands.size()); + Jobs::RunParallelFor(position_update->commands.size(), [&](const size_t i) { + position_update->commands[i].strand_segment_index = dynamic_strands.strands[i].begin_segment_handle; + position_update->commands[i].new_position = + dynamic_strands.strand_segments[position_update->commands[i].strand_segment_index].x0; + }); } void BrokenBranch::Serialize(YAML::Emitter& out) const { @@ -22,117 +40,148 @@ void BrokenBranch::Deserialize(const YAML::Node& in) { bool BrokenBranch::OnInspect(const std::shared_ptr& editor_layer) { static bool auto_subdivide = true; - static float segment_length = 0.01f; + static float segment_length = 0.05f; ImGui::DragFloat("Subdivision length", &segment_length, .0001f, 0.0001f, 1.f, "%.5f"); ImGui::Checkbox("Auto subdivide", &auto_subdivide); if (EditorLayer::DragAndDropButton(tree_ref, "Download Strands from Tree...")) { if (const auto tree = tree_ref.Get()) { - if (tree->strand_model.strand_model_skeleton.data.strand_group.PeekStrands().empty()) { - tree->BuildStrandModel(); - } + tree->BuildStrandModel(); strand_group = tree->strand_model.strand_model_skeleton.data.strand_group; if (auto_subdivide) { - Subdivide(segment_length, strand_group, subdivided_strand_group); - dynamic_strands.UpdateData({}, subdivided_strand_group.PeekStrands(), - subdivided_strand_group.PeekStrandSegments()); + Subdivide(segment_length, strand_group); } tree_ref.Clear(); } } if (ImGui::Button("Subdivide strands")) { - Subdivide(segment_length, strand_group, subdivided_strand_group); - dynamic_strands.UpdateData({}, subdivided_strand_group.PeekStrands(), subdivided_strand_group.PeekStrandSegments()); + Subdivide(segment_length, strand_group); } if (ImGui::Button("Experiment")) { ExperimentSetup(); } - - ImGui::Text((std::string("Original strand count: ") + std::to_string(strand_group.PeekStrands().size())).c_str()); - ImGui::Text( - (std::string("Original strand segment count: ") + std::to_string(strand_group.PeekStrandSegments().size())) - .c_str()); - ImGui::Text((std::string("Subdivided strand count: ") + std::to_string(subdivided_strand_group.PeekStrands().size())) - .c_str()); - ImGui::Text((std::string("Subdivided strand segment count: ") + - std::to_string(subdivided_strand_group.PeekStrandSegments().size())) - .c_str()); - - static float min_trunk_length = 0.3f; - static float max_trunk_length = 0.4f; - ImGui::DragFloat("Min trunk length", &min_trunk_length, .01f, 0.01f, max_trunk_length, "%.5f"); - ImGui::DragFloat("Max trunk length", &max_trunk_length, .01f, min_trunk_length, 1.f, "%.5f"); - - static float noise_frequency = 100.f; - ImGui::DragFloat("Noise frequency", &noise_frequency, .1f, 0.1, 100.f, "%.1f"); - static bool keep_upper = true; - static glm::vec3 upper_offset = glm::vec3(0, 0.1f, 0.0); - ImGui::Checkbox("Keep upper", &keep_upper); - if (keep_upper) { - ImGui::DragFloat3("Upper offset", &upper_offset.x, 0.01, 0.0f, 1.0f); + if (ImGui::TreeNode("Stats")) { + ImGui::Text((std::string("Original strand count: ") + std::to_string(strand_group.PeekStrands().size())).c_str()); + ImGui::Text( + (std::string("Original strand segment count: ") + std::to_string(strand_group.PeekStrandSegments().size())) + .c_str()); + ImGui::Text( + (std::string("Subdivided strand count: ") + std::to_string(subdivided_strand_group.PeekStrands().size())) + .c_str()); + ImGui::Text((std::string("Subdivided strand segment count: ") + + std::to_string(subdivided_strand_group.PeekStrandSegments().size())) + .c_str()); + ImGui::TreePop(); } - if (ImGui::Button("Cut trunk")) { - auto temp_strand_group = strand_group; - const auto size = temp_strand_group.PeekStrands().size(); - for (StrandHandle strand_handle = 0; strand_handle < temp_strand_group.PeekStrands().size(); strand_handle++) { - const auto& strand = temp_strand_group.PeekStrand(strand_handle); - StrandSegmentHandle segment_handle; - float t; - temp_strand_group.FindStrandT(strand_handle, segment_handle, t, - min_trunk_length + glm::perlin(strand.start_position * noise_frequency) * - (max_trunk_length - min_trunk_length)); - if (t <= 0.f) - continue; - const auto new_strand_handle = temp_strand_group.Cut(segment_handle, t); - if (new_strand_handle == -1) - continue; - auto& new_strand = temp_strand_group.RefStrand(new_strand_handle); - new_strand.start_position += upper_offset; - for (const auto& i : new_strand.PeekStrandSegmentHandles()) { - temp_strand_group.RefStrandSegment(i).end_position += upper_offset; - } - if (!keep_upper && new_strand_handle >= size) { - temp_strand_group.RemoveStrand(new_strand_handle); + if (ImGui::TreeNode("Trunk cutting")) { + static float min_trunk_length = 0.3f; + static float max_trunk_length = 0.4f; + ImGui::DragFloat("Min trunk length", &min_trunk_length, .01f, 0.01f, max_trunk_length, "%.5f"); + ImGui::DragFloat("Max trunk length", &max_trunk_length, .01f, min_trunk_length, 1.f, "%.5f"); + + static float noise_frequency = 100.f; + ImGui::DragFloat("Noise frequency", &noise_frequency, .1f, 0.1, 100.f, "%.1f"); + static bool keep_upper = true; + static glm::vec3 upper_offset = glm::vec3(0, 0.1f, 0.0); + ImGui::Checkbox("Keep upper", &keep_upper); + if (keep_upper) { + ImGui::DragFloat3("Upper offset", &upper_offset.x, 0.01, 0.0f, 1.0f); + } + if (ImGui::Button("Cut trunk")) { + auto temp_strand_group = strand_group; + const auto size = temp_strand_group.PeekStrands().size(); + for (StrandHandle strand_handle = 0; strand_handle < temp_strand_group.PeekStrands().size(); strand_handle++) { + const auto& strand = temp_strand_group.PeekStrand(strand_handle); + StrandSegmentHandle segment_handle; + float t; + temp_strand_group.FindStrandT(strand_handle, segment_handle, t, + min_trunk_length + glm::perlin(strand.start_position * noise_frequency) * + (max_trunk_length - min_trunk_length)); + if (t <= 0.f) + continue; + const auto new_strand_handle = temp_strand_group.Cut(segment_handle, t); + if (new_strand_handle == -1) + continue; + auto& new_strand = temp_strand_group.RefStrand(new_strand_handle); + new_strand.start_position += upper_offset; + for (const auto& i : new_strand.PeekStrandSegmentHandles()) { + temp_strand_group.RefStrandSegment(i).end_position += upper_offset; + } + if (!keep_upper && new_strand_handle >= size) { + temp_strand_group.RemoveStrand(new_strand_handle); + } } + Subdivide(segment_length, temp_strand_group); } - Subdivide(segment_length, temp_strand_group, subdivided_strand_group); - - } - if (ImGui::Button("Build particles for original")) { - InitializeStrandParticles(strand_group); - } - ImGui::SameLine(); - if (ImGui::Button("Build particles for subdivided")) { - InitializeStrandParticles(subdivided_strand_group); + ImGui::TreePop(); } - ImGui::SameLine(); - if (ImGui::Button("Clear particles")) { - ClearStrandParticles(); + if (ImGui::TreeNode("Particles")) { + if (ImGui::Button("Build particles for original")) { + InitializeStrandParticles(strand_group); + } + + if (ImGui::Button("Build particles for subdivided")) { + InitializeStrandParticles(subdivided_strand_group); + } + if (ImGui::Button("Clear particles")) { + ClearStrandParticles(); + } + ImGui::TreePop(); } - static float limit = 0.02f; - ImGui::DragFloat("Limit", &limit); - if (ImGui::Button("Build concave mesh")) { - InitializeStrandConcaveMesh(subdivided_strand_group, limit); + ImGui::Checkbox("Physics", &step_parameters.physics); + ImGui::Checkbox("Render", &step_parameters.render); + if (!step_parameters.physics) { + if (ImGui::Button("Simulate 1 step")) { + step_parameters.physics = true; + const bool resume_render = step_parameters.render; + step_parameters.render = false; + dynamic_strands.Step(step_parameters); + step_parameters.physics = false; + step_parameters.render = resume_render; + } } - if (ImGui::Button("Clear concave mesh")) { - ClearStrandConcaveMesh(); + if (!step_parameters.physics && !step_parameters.render) { + if (ImGui::Button("Download strands")) { + dynamic_strands.Download(); + EVOENGINE_LOG("Downloaded data from GPU") + } } - ImGui::Checkbox("Render strands", &enable_simulation); - return false; } - -void BrokenBranch::Update() { -} - void BrokenBranch::LateUpdate() { - if (enable_simulation) { + if (enable_simulation && !dynamic_strands.strand_segments.empty()) { + const auto owner = GetOwner(); + const auto scene = GetScene(); + const auto current_root_transform = scene->GetDataComponent(owner); + GlobalTransform original_inverse_root_transform; + original_inverse_root_transform.value = glm::inverse(initialize_parameters.root_transform.value); + const auto& position_update = std::dynamic_pointer_cast(dynamic_strands.operators.back()); + Jobs::RunParallelFor(position_update->commands.size(), [&](const size_t i) { + position_update->commands[i].strand_segment_index = dynamic_strands.strands[i].begin_segment_handle; + position_update->commands[i] + .new_position = current_root_transform.TransformPoint(original_inverse_root_transform.TransformPoint( + dynamic_strands.strand_segments[position_update->commands[i].strand_segment_index].x0)); + }); const auto editor_layer = Application::GetLayer(); - dynamic_strands.step_parameters.render_parameters.target_camera = editor_layer->GetSceneCamera(); - dynamic_strands.Step(); + step_parameters.render_parameters.target_camera = editor_layer->GetSceneCamera(); + dynamic_strands.Step(step_parameters); } } +void BrokenBranch::FixedUpdate() { + +} + void BrokenBranch::OnCreate() { + dynamic_strands.operators.clear(); + dynamic_strands.constraints.clear(); + + dynamic_strands.operators.emplace_back(std::make_shared()); + dynamic_strands.operators.emplace_back(std::make_shared()); + + dynamic_strands.pre_step = std::make_shared(); + dynamic_strands.constraints.emplace_back(std::make_shared()); + + dynamic_strands.post_step = std::make_shared(); } void BrokenBranch::OnDestroy() { @@ -153,10 +202,10 @@ void BrokenBranch::ExperimentSetup() { strand_group.CalculateRotations(); } -void BrokenBranch::Subdivide(const float segment_length, const StrandModelStrandGroup& src, - StrandModelStrandGroup& dst) { - src.UniformlySubdivide(dst, segment_length, segment_length * .01f); - dst.RandomAssignColor(); +void BrokenBranch::Subdivide(const float segment_length, const StrandModelStrandGroup& src) { + src.UniformlySubdivide(subdivided_strand_group, segment_length, segment_length * .01f); + subdivided_strand_group.RandomAssignColor(); + UpdateDynamicStrands(); } void BrokenBranch::InitializeStrandParticles(const StrandModelStrandGroup& target_strand_group) const { diff --git a/EvoEngine_Plugins/EcoSysLab/src/DynamicStrands.cpp b/EvoEngine_Plugins/EcoSysLab/src/DynamicStrands.cpp index f893ee8..fb8ed5e 100644 --- a/EvoEngine_Plugins/EcoSysLab/src/DynamicStrands.cpp +++ b/EvoEngine_Plugins/EcoSysLab/src/DynamicStrands.cpp @@ -1,10 +1,20 @@ #include "DynamicStrands.hpp" + +#include "DynamicStrandsPhysics.hpp" #include "Shader.hpp" +#include "glm/gtc/matrix_access.hpp" +#include "glm/gtx/quaternion.hpp" using namespace eco_sys_lab_plugin; DynamicStrands::DynamicStrands() { - const auto max_frame_count = Platform::GetMaxFramesInFlight(); - assert(max_frame_count == 2); + if (!strands_layout) { + strands_layout = std::make_shared(); + strands_layout->PushDescriptorBinding(0, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_ALL, 0); + strands_layout->PushDescriptorBinding(1, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_ALL, 0); + + strands_layout->Initialize(); + } + VkBufferCreateInfo buffer_create_info{}; buffer_create_info.sType = VK_STRUCTURE_TYPE_BUFFER_CREATE_INFO; buffer_create_info.usage = VK_BUFFER_USAGE_TRANSFER_DST_BIT | VK_BUFFER_USAGE_STORAGE_BUFFER_BIT; @@ -12,158 +22,78 @@ DynamicStrands::DynamicStrands() { buffer_create_info.size = 1; VmaAllocationCreateInfo buffer_vma_allocation_create_info{}; buffer_vma_allocation_create_info.usage = VMA_MEMORY_USAGE_AUTO_PREFER_DEVICE; - current_version = 0; - for (uint32_t i = 0; i < max_frame_count; i++) { - device_ref_strand_segments_buffer[i] = - std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); - device_ref_strands_buffer[i] = std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); - device_ref_strand_segment_particles_buffer[i] = - std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); - - device_strand_segments_buffer[i] = std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); - device_strands_buffer[i] = std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); - device_strand_segment_particles_buffer[i] = - std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); - - device_buffer_version[i] = current_version; - } - if (!strands_layout) { - strands_layout = std::make_shared(); - strands_layout->PushDescriptorBinding(0, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_ALL, 0); - strands_layout->PushDescriptorBinding(1, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_ALL, 0); - strands_layout->PushDescriptorBinding(2, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_ALL, 0); - strands_layout->PushDescriptorBinding(3, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_ALL, 0); - strands_layout->PushDescriptorBinding(4, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_ALL, 0); - strands_layout->PushDescriptorBinding(5, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_ALL, 0); - strands_layout->Initialize(); - } -} + device_strand_segments_buffer = std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); + device_strands_buffer = std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); -void DynamicStrands::Step() { - Step(step_parameters); + const auto max_frame_in_flight = Platform::GetMaxFramesInFlight(); + strands_descriptor_sets.resize(max_frame_in_flight); + for (auto& i : strands_descriptor_sets) { + i = std::make_shared(strands_layout); + } } void DynamicStrands::Step(const StepParameters& target_step_parameters) { - current_left = !current_left; - - const auto current_frame_index = current_left ? 0 : 1; - if (device_buffer_version[current_frame_index] != current_version) { - Upload(); - device_buffer_version[current_frame_index] = current_version; - } + if (strand_segments.empty()) + return; + const auto current_frame_index = Platform::GetCurrentFrameIndex(); + + strands_descriptor_sets[current_frame_index]->UpdateBufferDescriptorBinding(0, device_strand_segments_buffer, 0); + strands_descriptor_sets[current_frame_index]->UpdateBufferDescriptorBinding(1, device_strands_buffer, 0); if (target_step_parameters.physics) { - Physics(target_step_parameters.physics_parameters); + Physics(target_step_parameters.physics_parameters, operators, pre_step, constraints, post_step); } + if (target_step_parameters.render) { Render(target_step_parameters.render_parameters); } } -void DynamicStrands::UpdateData(const InitializeParameters& initialize_parameters, - const std::vector& target_strands, - const std::vector& target_strand_segments) { - Clear(); - assert(initialize_parameters.root_transform.GetScale() == glm::vec3(1.0f)); - ref_strands.resize(target_strands.size()); - ref_strand_segments.resize(target_strand_segments.size()); - ref_strand_segment_particles.resize(ref_strands.size() + ref_strand_segments.size()); - int particle_index = 0; - for (uint32_t strand_index = 0; strand_index < target_strands.size(); strand_index++) { - auto& target_strand = target_strands[strand_index]; - auto& first_particle = ref_strand_segment_particles[particle_index]; - first_particle.position = target_strand.start_position; - first_particle.thickness = target_strand.start_thickness; - first_particle.color = target_strand.start_color; - particle_index++; - const auto& handles = target_strand.PeekStrandSegmentHandles(); - for (const auto& i : handles) { - ref_strand_segments[i].start_particle_index = particle_index - 1; - ref_strand_segments[i].end_particle_index = particle_index; - auto& particle = ref_strand_segment_particles[particle_index]; - particle.position = target_strand_segments[i].end_position; - particle.thickness = target_strand_segments[i].end_thickness; - particle.color = target_strand_segments[i].end_color; - particle_index++; - } - ref_strands[strand_index].first_strand_segment_handle = handles.front(); - ref_strands[strand_index].last_strand_segment_handle = handles.back(); - } - Jobs::RunParallelFor(target_strand_segments.size(), [&](const size_t i) { - ref_strand_segments[i].prev_handle = target_strand_segments[i].GetPrevHandle(); - ref_strand_segments[i].next_handle = target_strand_segments[i].GetNextHandle(); - ref_strand_segments[i].strand_handle = target_strand_segments[i].GetStrandHandle(); - ref_strand_segments[i].index = target_strand_segments[i].GetIndex(); - ref_strand_segments[i].rotation = - initialize_parameters.root_transform.GetRotation() * target_strand_segments[i].rotation; - }); - current_version++; +void DynamicStrands::Upload() const { + device_strands_buffer->UploadVector(strands); + device_strand_segments_buffer->UploadVector(strand_segments); } -/* -void DynamicStrands::Simulate(const SimulateParameters& simulate_parameters) { - //1. Update velocity. - Jobs::RunParallelFor(ref_strand_segments.size(), [&](const auto& i) { - auto& segment = ref_strand_segments[i]; - segment.end_velocity += glm::vec3(0, -9.81f, 0) * simulate_parameters.d_t * segment.inv_end_mass; - }); - - //2. Update position. - Jobs::RunParallelFor(ref_strand_segments.size(), [&](const auto& i) { - auto& segment = ref_strand_segments[i]; - segment.end_position += segment.end_velocity * simulate_parameters.d_t; - }); - - //3. Update angular velocity. - Jobs::RunParallelFor(ref_strand_segments.size(), [&](const auto& i) { - auto& segment = ref_strand_segments[i]; - }); - - // 4. Update rotation. - Jobs::RunParallelFor(ref_strand_segments.size(), [&](const auto& i) { - auto& segment = ref_strand_segments[i]; - auto dq = glm::quat(0.f, segment.angular_velocity * simulate_parameters.d_t); - segment.rotation += dq * segment.rotation * .5f; - segment.rotation = glm::normalize(segment.rotation); - }); - - for (uint32_t iteration = 0; iteration < simulate_parameters.solver_iterations; iteration++) { - for (uint32_t i = 0; i < ref_strand_segments.size(); i++) { - if (i % 2 == 0) { - - } - } +void DynamicStrands::Download() { + device_strands_buffer->DownloadVector(strands, strands.size()); + device_strand_segments_buffer->DownloadVector(strand_segments, strand_segments.size()); + for (const auto& op : operators) { + op->DownloadData(); + } + for (const auto& c : constraints) { + c->DownloadData(); } -}*/ - -void DynamicStrands::Upload() const { - const auto current_frame_index = current_left ? 0 : 1; - device_ref_strand_segments_buffer[current_frame_index]->UploadVector(ref_strand_segments); - device_ref_strands_buffer[current_frame_index]->UploadVector(ref_strands); - device_ref_strand_segment_particles_buffer[current_frame_index]->UploadVector(ref_strand_segment_particles); - - device_strand_segments_buffer[current_frame_index]->UploadVector(ref_strand_segments); - device_strands_buffer[current_frame_index]->UploadVector(ref_strands); - device_strand_segment_particles_buffer[current_frame_index]->UploadVector(ref_strand_segment_particles); } void DynamicStrands::Clear() { - ref_strand_segments.clear(); - ref_strands.clear(); - ref_strand_segment_particles.clear(); + strands.clear(); + strand_segments.clear(); +} + +glm::vec3 DynamicStrands::ComputeInertiaTensorBox(const float mass, const float width, const float height, const float depth) { + return { + mass / 12.f * (height * height + depth * depth), + mass / 12.f * (width * width + depth * depth), + mass / 12.f * (width * width + height * height), + }; +} + +glm::vec3 DynamicStrands::ComputeInertiaTensorRod(const float mass, const float radius, const float length) { + return { + mass / 12.f * (radius * radius + length * length), + mass / 12.f * (radius * radius + length * length), + mass / 4.f * (radius * radius + radius * radius), + }; +} + +void DynamicStrands::InitializeOperators(const InitializeParameters& initialize_parameters) const { + for (const auto& i : operators) + i->InitializeData(initialize_parameters, *this); } -void DynamicStrands::BindStrandsDescriptorSet(const std::shared_ptr& target_descriptor_set) const { - const auto current_frame_index = current_left ? 0 : 1; - target_descriptor_set->UpdateBufferDescriptorBinding(0, device_ref_strand_segments_buffer[current_frame_index], 0); - target_descriptor_set->UpdateBufferDescriptorBinding(1, device_ref_strands_buffer[current_frame_index], 0); - target_descriptor_set->UpdateBufferDescriptorBinding( - 2, device_ref_strand_segment_particles_buffer[current_frame_index], 0); - - target_descriptor_set->UpdateBufferDescriptorBinding(3, device_strand_segments_buffer[current_frame_index], 0); - target_descriptor_set->UpdateBufferDescriptorBinding(4, device_strands_buffer[current_frame_index], 0); - target_descriptor_set->UpdateBufferDescriptorBinding(5, device_strand_segment_particles_buffer[current_frame_index], - 0); +void DynamicStrands::InitializeConstraints(const InitializeParameters& initialize_parameters) const { + for (const auto& i : constraints) + i->InitializeData(initialize_parameters, *this); } diff --git a/EvoEngine_Plugins/EcoSysLab/src/DynamicStrandsPhysics.cpp b/EvoEngine_Plugins/EcoSysLab/src/DynamicStrandsPhysics.cpp index f7c0188..8c5bfe1 100644 --- a/EvoEngine_Plugins/EcoSysLab/src/DynamicStrandsPhysics.cpp +++ b/EvoEngine_Plugins/EcoSysLab/src/DynamicStrandsPhysics.cpp @@ -1,66 +1,514 @@ -#include "DynamicStrands.hpp" +#include "DynamicStrandsPhysics.hpp" #include "Shader.hpp" using namespace eco_sys_lab_plugin; -void DynamicStrands::Physics(const StepParameters::PhysicsParameters& physics_parameters) const { - const auto current_frame_index = current_left ? 0 : 1; - static std::shared_ptr precompute_pipeline{}; - static std::vector> strands_descriptor_sets{}; - - struct PrecomputePushConstant { - uint32_t strand_size = 0; - uint32_t strand_segment_size = 0; - }; - - if (!precompute_pipeline) { - const auto max_frame_in_flight = Platform::GetMaxFramesInFlight(); - strands_descriptor_sets.resize(max_frame_in_flight); - for (auto& i : strands_descriptor_sets) { - i = std::make_shared(strands_layout); - } - static std::shared_ptr precompute_shader{}; - precompute_shader = std::make_shared(); - precompute_shader->Set( +void DynamicStrands::Physics(const PhysicsParameters& physics_parameters, + const std::vector>& target_operators, + const std::shared_ptr& target_pre_step, + const std::vector>& target_constraints, + const std::shared_ptr& target_post_step) { + for (const auto& op : target_operators) { + op->Execute(physics_parameters, *this); + } + if (target_pre_step) + target_pre_step->Execute(physics_parameters, *this); + for (const auto& c : target_constraints) { + c->Project(physics_parameters, *this); + } + if (target_post_step) + target_post_step->Execute(physics_parameters, *this); +} + +DynamicStrandsPreStep::DynamicStrandsPreStep() { + if (!pre_step_pipeline) { + static std::shared_ptr shader{}; + shader = std::make_shared(); + shader->Set( ShaderType::Compute, Platform::Constants::shader_global_defines, - std::filesystem::path("./EcoSysLabResources") / "Shaders/Compute/DynamicStrandsPrecompute.comp"); + std::filesystem::path("./EcoSysLabResources") / "Shaders/Compute/DynamicStrands/PreStep.comp"); + + pre_step_pipeline = std::make_shared(); + pre_step_pipeline->compute_shader = shader; + pre_step_pipeline->descriptor_set_layouts.emplace_back(DynamicStrands::strands_layout); + + auto& push_constant_range = pre_step_pipeline->push_constant_ranges.emplace_back(); + push_constant_range.size = sizeof(PreStepPushConstant); + push_constant_range.offset = 0; + push_constant_range.stageFlags = VK_SHADER_STAGE_COMPUTE_BIT; + + pre_step_pipeline->Initialize(); + } +} + +void DynamicStrandsPreStep::Execute(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) { + const auto current_frame_index = Platform::GetCurrentFrameIndex(); + + PreStepPushConstant push_constant; + push_constant.segment_size = target_dynamic_strands.strand_segments.size(); + push_constant.time_step = physics_parameters.time_step; + push_constant.inv_time_step = 1.f / push_constant.time_step; + const uint32_t task_work_group_invocations = + Platform::GetSelectedPhysicalDevice()->mesh_shader_properties_ext.maxPreferredTaskWorkGroupInvocations; + + Platform::RecordCommandsMainQueue([&](const VkCommandBuffer vk_command_buffer) { + pre_step_pipeline->Bind(vk_command_buffer); + pre_step_pipeline->BindDescriptorSet( + vk_command_buffer, 0, + target_dynamic_strands.strands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + + pre_step_pipeline->PushConstant(vk_command_buffer, 0, push_constant); + vkCmdDispatch(vk_command_buffer, Platform::DivUp(push_constant.segment_size, task_work_group_invocations), 1, 1); + Platform::EverythingBarrier(vk_command_buffer); + }); +} + +DynamicStrandsPostStep::DynamicStrandsPostStep() { + if (!post_step_pipeline) { + static std::shared_ptr shader{}; + shader = std::make_shared(); + shader->Set(ShaderType::Compute, Platform::Constants::shader_global_defines, + std::filesystem::path("./EcoSysLabResources") / "Shaders/Compute/DynamicStrands/PostStep.comp"); + + post_step_pipeline = std::make_shared(); + post_step_pipeline->compute_shader = shader; + post_step_pipeline->descriptor_set_layouts.emplace_back(DynamicStrands::strands_layout); + + auto& push_constant_range = post_step_pipeline->push_constant_ranges.emplace_back(); + push_constant_range.size = sizeof(PostStepPushConstant); + push_constant_range.offset = 0; + push_constant_range.stageFlags = VK_SHADER_STAGE_COMPUTE_BIT; + + post_step_pipeline->Initialize(); + } +} + +void DynamicStrandsPostStep::Execute(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) { + const auto current_frame_index = Platform::GetCurrentFrameIndex(); + + PostStepPushConstant push_constant; + push_constant.segment_size = target_dynamic_strands.strand_segments.size(); + const uint32_t task_work_group_invocations = + Platform::GetSelectedPhysicalDevice()->mesh_shader_properties_ext.maxPreferredTaskWorkGroupInvocations; + + Platform::RecordCommandsMainQueue([&](const VkCommandBuffer vk_command_buffer) { + post_step_pipeline->Bind(vk_command_buffer); + post_step_pipeline->BindDescriptorSet( + vk_command_buffer, 0, + target_dynamic_strands.strands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + + post_step_pipeline->PushConstant(vk_command_buffer, 0, push_constant); + vkCmdDispatch(vk_command_buffer, Platform::DivUp(push_constant.segment_size, task_work_group_invocations), 1, 1); + Platform::EverythingBarrier(vk_command_buffer); + }); +} + +DsPositionUpdate::DsPositionUpdate() { + if (!layout) { + layout = std::make_shared(); + layout->PushDescriptorBinding(0, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_COMPUTE_BIT, 0); + layout->Initialize(); + } + VkBufferCreateInfo buffer_create_info{}; + buffer_create_info.sType = VK_STRUCTURE_TYPE_BUFFER_CREATE_INFO; + buffer_create_info.usage = VK_BUFFER_USAGE_TRANSFER_DST_BIT | VK_BUFFER_USAGE_STORAGE_BUFFER_BIT; + buffer_create_info.sharingMode = VK_SHARING_MODE_EXCLUSIVE; + buffer_create_info.size = 1; + VmaAllocationCreateInfo buffer_vma_allocation_create_info{}; + buffer_vma_allocation_create_info.usage = VMA_MEMORY_USAGE_AUTO_PREFER_DEVICE; + const auto max_frame_in_flight = Platform::GetMaxFramesInFlight(); + + commands_buffer.resize(max_frame_in_flight); + + for (int frame_index = 0; frame_index < max_frame_in_flight; frame_index++) { + commands_buffer[frame_index] = std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); + } + + if (!external_force_pipeline) { + static std::shared_ptr shader{}; + shader = std::make_shared(); + shader->Set( + ShaderType::Compute, Platform::Constants::shader_global_defines, + std::filesystem::path("./EcoSysLabResources") / "Shaders/Compute/DynamicStrands/Operators/PositionUpdate.comp"); + + external_force_pipeline = std::make_shared(); + external_force_pipeline->compute_shader = shader; + + external_force_pipeline->descriptor_set_layouts.emplace_back(DynamicStrands::strands_layout); + external_force_pipeline->descriptor_set_layouts.emplace_back(layout); + + auto& push_constant_range = external_force_pipeline->push_constant_ranges.emplace_back(); + push_constant_range.size = sizeof(PositionUpdatePushConstant); + push_constant_range.offset = 0; + push_constant_range.stageFlags = VK_SHADER_STAGE_COMPUTE_BIT; + + external_force_pipeline->Initialize(); + } + + commands_descriptor_sets.resize(max_frame_in_flight); + for (auto& i : commands_descriptor_sets) { + i = std::make_shared(layout); + } +} + +void DsPositionUpdate::Execute(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) { + if (commands.empty()) + return; + + const auto current_frame_index = Platform::GetCurrentFrameIndex(); + commands_buffer[current_frame_index]->UploadVector(commands); + commands_descriptor_sets[current_frame_index]->UpdateBufferDescriptorBinding(0, commands_buffer[current_frame_index]); + + PositionUpdatePushConstant push_constant; + push_constant.commands_size = commands.size(); + const uint32_t task_work_group_invocations = + Platform::GetSelectedPhysicalDevice()->mesh_shader_properties_ext.maxPreferredTaskWorkGroupInvocations; + + Platform::RecordCommandsMainQueue([&](const VkCommandBuffer vk_command_buffer) { + external_force_pipeline->Bind(vk_command_buffer); + external_force_pipeline->BindDescriptorSet( + vk_command_buffer, 0, + target_dynamic_strands.strands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + external_force_pipeline->BindDescriptorSet(vk_command_buffer, 1, + commands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + + external_force_pipeline->PushConstant(vk_command_buffer, 0, push_constant); + + vkCmdDispatch(vk_command_buffer, Platform::DivUp(push_constant.commands_size, task_work_group_invocations), 1, 1); + Platform::EverythingBarrier(vk_command_buffer); + }); +} + +DsExternalForce::DsExternalForce() { + if (!layout) { + layout = std::make_shared(); + layout->PushDescriptorBinding(0, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_COMPUTE_BIT, 0); + + layout->Initialize(); + } + + VkBufferCreateInfo buffer_create_info{}; + buffer_create_info.sType = VK_STRUCTURE_TYPE_BUFFER_CREATE_INFO; + buffer_create_info.usage = VK_BUFFER_USAGE_TRANSFER_DST_BIT | VK_BUFFER_USAGE_STORAGE_BUFFER_BIT; + buffer_create_info.sharingMode = VK_SHARING_MODE_EXCLUSIVE; + buffer_create_info.size = 1; + VmaAllocationCreateInfo buffer_vma_allocation_create_info{}; + buffer_vma_allocation_create_info.usage = VMA_MEMORY_USAGE_AUTO_PREFER_DEVICE; + const auto max_frame_in_flight = Platform::GetMaxFramesInFlight(); - precompute_pipeline = std::make_shared(); - precompute_pipeline->compute_shader = precompute_shader; + commands_buffer.resize(max_frame_in_flight); - auto per_frame_layout = Platform::GetDescriptorSetLayout("PER_FRAME_LAYOUT"); - precompute_pipeline->descriptor_set_layouts.emplace_back(strands_layout); + for (int frame_index = 0; frame_index < max_frame_in_flight; frame_index++) { + commands_buffer[frame_index] = std::make_shared(buffer_create_info, buffer_vma_allocation_create_info); + } + + if (!external_force_pipeline) { + static std::shared_ptr shader{}; + shader = std::make_shared(); + shader->Set( + ShaderType::Compute, Platform::Constants::shader_global_defines, + std::filesystem::path("./EcoSysLabResources") / "Shaders/Compute/DynamicStrands/Operators/ExternalForce.comp"); + external_force_pipeline = std::make_shared(); + external_force_pipeline->compute_shader = shader; - auto& push_constant_range = precompute_pipeline->push_constant_ranges.emplace_back(); - push_constant_range.size = sizeof(PrecomputePushConstant); + external_force_pipeline->descriptor_set_layouts.emplace_back(DynamicStrands::strands_layout); + external_force_pipeline->descriptor_set_layouts.emplace_back(layout); + + auto& push_constant_range = external_force_pipeline->push_constant_ranges.emplace_back(); + push_constant_range.size = sizeof(ExternalForcePushConstant); push_constant_range.offset = 0; - push_constant_range.stageFlags = VK_SHADER_STAGE_ALL; + push_constant_range.stageFlags = VK_SHADER_STAGE_COMPUTE_BIT; - precompute_pipeline->Initialize(); + external_force_pipeline->Initialize(); + } + commands_descriptor_sets.resize(max_frame_in_flight); + for (auto& i : commands_descriptor_sets) { + i = std::make_shared(layout); } - if (ref_strand_segments.empty()) +} + +void DsExternalForce::Execute(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) { + if (commands.empty()) return; - BindStrandsDescriptorSet(strands_descriptor_sets[current_frame_index]); - PrecomputePushConstant push_constant{}; - push_constant.strand_segment_size = ref_strand_segments.size(); - push_constant.strand_size = ref_strands.size(); + const auto current_frame_index = Platform::GetCurrentFrameIndex(); + commands_buffer[current_frame_index]->UploadVector(commands); + + commands_descriptor_sets[current_frame_index]->UpdateBufferDescriptorBinding(0, commands_buffer[current_frame_index]); + + ExternalForcePushConstant push_constant; + push_constant.commands_size = commands.size(); const uint32_t task_work_group_invocations = Platform::GetSelectedPhysicalDevice()->mesh_shader_properties_ext.maxPreferredTaskWorkGroupInvocations; Platform::RecordCommandsMainQueue([&](const VkCommandBuffer vk_command_buffer) { - precompute_pipeline->Bind(vk_command_buffer); - precompute_pipeline->BindDescriptorSet(vk_command_buffer, 0, - strands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); - precompute_pipeline->PushConstant(vk_command_buffer, 0, push_constant); - /** - * Dispatch! - */ - vkCmdDispatch(vk_command_buffer, Platform::DivUp(ref_strand_segments.size(), task_work_group_invocations), 1, 1); - /** - * Remember, many of vulkan commands are executed without ordering. So we have this Platform::EverythingBarrier() to - * make sure that the above commands finishes before moving on. This is syncronization on GPU, not between GPU and - * CPU. - */ + external_force_pipeline->Bind(vk_command_buffer); + external_force_pipeline->BindDescriptorSet( + vk_command_buffer, 0, + target_dynamic_strands.strands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + external_force_pipeline->BindDescriptorSet(vk_command_buffer, 1, + commands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + + external_force_pipeline->PushConstant(vk_command_buffer, 0, push_constant); + + vkCmdDispatch(vk_command_buffer, Platform::DivUp(push_constant.commands_size, task_work_group_invocations), 1, 1); Platform::EverythingBarrier(vk_command_buffer); }); } + +DsStiffRod::DsStiffRod() { + if (!layout) { + layout = std::make_shared(); + layout->PushDescriptorBinding(0, VK_DESCRIPTOR_TYPE_UNIFORM_BUFFER, VK_SHADER_STAGE_COMPUTE_BIT, 0); + layout->PushDescriptorBinding(1, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_COMPUTE_BIT, 0); + layout->PushDescriptorBinding(2, VK_DESCRIPTOR_TYPE_STORAGE_BUFFER, VK_SHADER_STAGE_COMPUTE_BIT, 0); + + layout->Initialize(); + } + + VkBufferCreateInfo uniform_buffer_create_info{}; + uniform_buffer_create_info.sType = VK_STRUCTURE_TYPE_BUFFER_CREATE_INFO; + uniform_buffer_create_info.usage = VK_BUFFER_USAGE_TRANSFER_DST_BIT | VK_BUFFER_USAGE_UNIFORM_BUFFER_BIT; + uniform_buffer_create_info.sharingMode = VK_SHARING_MODE_EXCLUSIVE; + uniform_buffer_create_info.size = 1; + VkBufferCreateInfo storage_buffer_create_info{}; + storage_buffer_create_info.sType = VK_STRUCTURE_TYPE_BUFFER_CREATE_INFO; + storage_buffer_create_info.usage = VK_BUFFER_USAGE_TRANSFER_DST_BIT | VK_BUFFER_USAGE_STORAGE_BUFFER_BIT; + storage_buffer_create_info.sharingMode = VK_SHARING_MODE_EXCLUSIVE; + storage_buffer_create_info.size = 1; + VmaAllocationCreateInfo buffer_vma_allocation_create_info{}; + buffer_vma_allocation_create_info.usage = VMA_MEMORY_USAGE_AUTO_PREFER_DEVICE; + const auto max_frame_in_flight = Platform::GetMaxFramesInFlight(); + + rod_properties_buffer.resize(max_frame_in_flight); + + for (int frame_index = 0; frame_index < max_frame_in_flight; frame_index++) { + rod_properties_buffer[frame_index] = + std::make_shared(uniform_buffer_create_info, buffer_vma_allocation_create_info); + } + + per_strand_data_list_buffer = std::make_shared(storage_buffer_create_info, buffer_vma_allocation_create_info); + rod_constraints_buffer = std::make_shared(storage_buffer_create_info, buffer_vma_allocation_create_info); + + if (!init_constraint_pipeline) { + static std::shared_ptr init_constraint_shader{}; + init_constraint_shader = std::make_shared(); + init_constraint_shader->Set(ShaderType::Compute, Platform::Constants::shader_global_defines, + std::filesystem::path("./EcoSysLabResources") / + "Shaders/Compute/DynamicStrands/Constraints/StiffRodInitConstraints.comp"); + + init_constraint_pipeline = std::make_shared(); + init_constraint_pipeline->compute_shader = init_constraint_shader; + + init_constraint_pipeline->descriptor_set_layouts.emplace_back(DynamicStrands::strands_layout); + init_constraint_pipeline->descriptor_set_layouts.emplace_back(layout); + + auto& init_push_constant_range = init_constraint_pipeline->push_constant_ranges.emplace_back(); + init_push_constant_range.size = sizeof(StiffRodInitConstraintConstant); + init_push_constant_range.offset = 0; + init_push_constant_range.stageFlags = VK_SHADER_STAGE_COMPUTE_BIT; + + init_constraint_pipeline->Initialize(); + + static std::shared_ptr update_constraint_shader{}; + update_constraint_shader = std::make_shared(); + update_constraint_shader->Set(ShaderType::Compute, Platform::Constants::shader_global_defines, + std::filesystem::path("./EcoSysLabResources") / + "Shaders/Compute/DynamicStrands/Constraints/StiffRodUpdateConstraints.comp"); + + update_constraint_pipeline = std::make_shared(); + update_constraint_pipeline->compute_shader = update_constraint_shader; + + update_constraint_pipeline->descriptor_set_layouts.emplace_back(DynamicStrands::strands_layout); + update_constraint_pipeline->descriptor_set_layouts.emplace_back(layout); + + auto& update_push_constant_range = update_constraint_pipeline->push_constant_ranges.emplace_back(); + update_push_constant_range.size = sizeof(StiffRodUpdateConstraintConstant); + update_push_constant_range.offset = 0; + update_push_constant_range.stageFlags = VK_SHADER_STAGE_COMPUTE_BIT; + + update_constraint_pipeline->Initialize(); + + static std::shared_ptr project_constraint_shader{}; + project_constraint_shader = std::make_shared(); + project_constraint_shader->Set(ShaderType::Compute, Platform::Constants::shader_global_defines, + std::filesystem::path("./EcoSysLabResources") / + "Shaders/Compute/DynamicStrands/Constraints/StiffRodProjectConstraints.comp"); + + project_constraint_pipeline = std::make_shared(); + project_constraint_pipeline->compute_shader = project_constraint_shader; + + project_constraint_pipeline->descriptor_set_layouts.emplace_back(DynamicStrands::strands_layout); + project_constraint_pipeline->descriptor_set_layouts.emplace_back(layout); + + auto& project_push_constant_range = project_constraint_pipeline->push_constant_ranges.emplace_back(); + project_push_constant_range.size = sizeof(StiffRodProjectConstraintConstant); + project_push_constant_range.offset = 0; + project_push_constant_range.stageFlags = VK_SHADER_STAGE_COMPUTE_BIT; + + project_constraint_pipeline->Initialize(); + } + + if (strands_physics_descriptor_sets.empty()) { + strands_physics_descriptor_sets.resize(max_frame_in_flight); + for (auto& i : strands_physics_descriptor_sets) { + i = std::make_shared(layout); + } + } +} + +void DsStiffRod::InitializeData(const DynamicStrands::InitializeParameters& initialize_parameters, + const DynamicStrands& target_dynamic_strands) { + per_strand_data_list.resize(target_dynamic_strands.strands.size()); + rod_constraints.clear(); + rod_constraints.reserve(glm::max(0, static_cast(target_dynamic_strands.strand_segments.size()) - + static_cast(target_dynamic_strands.strands.size()))); + for (uint32_t strand_index = 0; strand_index < target_dynamic_strands.strands.size(); strand_index++) { + const auto& gpu_strand = target_dynamic_strands.strands[strand_index]; + auto& stiff_rod_per_strand_data = per_strand_data_list[strand_index]; + if (gpu_strand.begin_segment_handle == -1) { + stiff_rod_per_strand_data.begin_rod_constraint_handle = -1; + stiff_rod_per_strand_data.end_rod_constraint_handle = -1; + continue; + } + if (gpu_strand.begin_segment_handle == gpu_strand.end_segment_handle) { + stiff_rod_per_strand_data.begin_rod_constraint_handle = -1; + stiff_rod_per_strand_data.end_rod_constraint_handle = -1; + continue; + } + stiff_rod_per_strand_data.begin_rod_constraint_handle = static_cast(rod_constraints.size()); + int segment_handle = gpu_strand.begin_segment_handle; + int next_segment_handle = target_dynamic_strands.strand_segments[segment_handle].next_handle; + while (next_segment_handle != -1) { + rod_constraints.emplace_back(); + auto& rod_constraint = rod_constraints.back(); + rod_constraint.next_handle = rod_constraints.size(); + rod_constraint.segment0_index = segment_handle; + rod_constraint.segment1_index = next_segment_handle; + const auto& segment0 = target_dynamic_strands.strand_segments[segment_handle]; + const auto& segment1 = target_dynamic_strands.strand_segments[next_segment_handle]; + const auto& q0 = segment0.q0; + const auto& q1 = segment1.q0; + const auto& l0 = segment0.length; + const auto& l1 = segment1.length; + const float segment_length = (l0 + l1) * .5f; + rod_constraint.rest_darboux_vector = ComputeDarbouxVector(q0, q1, segment_length); + rod_constraint.average_segment_length = segment_length; + const auto rotation0_transpose = glm::transpose(glm::mat3_cast(q0)); + const auto rotation1_transpose = glm::transpose(glm::mat3_cast(q1)); + const auto& p0 = segment0.x0; + const auto& p1 = segment1.x0; + const auto constraint_position0 = p0 + q0 * glm::vec3(0, 0, -1) * segment0.length * .5f; + const auto constraint_position1 = p1 + q1 * glm::vec3(0, 0, 1) * segment1.length * .5f; + assert(glm::distance(constraint_position0, constraint_position1) < 0.0001f); + const auto constraint_position = (constraint_position0 + constraint_position1) * .5f; + rod_constraint.constraint_info = + glm::mat4(glm::vec4(rotation0_transpose * (constraint_position - p0), 0.f), + glm::vec4(rotation1_transpose * (constraint_position - p1), 0.f), + glm::vec4(constraint_position, 0.f), glm::vec4(constraint_position, 0.f)); + + // compute bending and torsion stiffness of the K matrix diagonal; assumption: the rod axis follows the y-axis of + // the local frame as with Blender's armatures +#define M_PI_4 0.785398163397448309616 // pi/4 + const auto second_moment_of_area(static_cast(M_PI_4) * + std::pow((segment0.radius + segment1.radius) * .5f, 4.f)); + const auto bending_stiffness(initialize_parameters.youngs_modulus * second_moment_of_area); + const auto torsion_stiffness(2.f * initialize_parameters.torsion_modulus * second_moment_of_area); + rod_constraint.stiffness_coefficient_k = glm::vec3(bending_stiffness, torsion_stiffness, bending_stiffness); + + segment_handle = next_segment_handle; + next_segment_handle = target_dynamic_strands.strand_segments[segment_handle].next_handle; + } + rod_constraints.back().next_handle = -1; + stiff_rod_per_strand_data.end_rod_constraint_handle = static_cast(rod_constraints.size()) - 1; + } + + per_strand_data_list_buffer->UploadVector(per_strand_data_list); + rod_constraints_buffer->UploadVector(rod_constraints); +} + +void DsStiffRod::Project(const DynamicStrands::PhysicsParameters& physics_parameters, + const DynamicStrands& target_dynamic_strands) { + const auto current_frame_index = Platform::GetCurrentFrameIndex(); + rod_properties_buffer[current_frame_index]->Upload(rod_properties); + + strands_physics_descriptor_sets[current_frame_index]->UpdateBufferDescriptorBinding( + 0, rod_properties_buffer[current_frame_index]); + strands_physics_descriptor_sets[current_frame_index]->UpdateBufferDescriptorBinding(1, per_strand_data_list_buffer); + strands_physics_descriptor_sets[current_frame_index]->UpdateBufferDescriptorBinding(2, rod_constraints_buffer); + + int num_sub_step = 1; + int num_iteration = 1; + + StiffRodInitConstraintConstant init_push_constant; + init_push_constant.constraint_size = rod_constraints.size(); + init_push_constant.time_step = physics_parameters.time_step / num_sub_step; + init_push_constant.inv_time_step = 1.f / init_push_constant.time_step; + + StiffRodUpdateConstraintConstant update_constraint_constant; + update_constraint_constant.constraint_size = rod_constraints.size(); + + StiffRodProjectConstraintConstant project_constraint_constant; + project_constraint_constant.strand_size = per_strand_data_list.size(); + + + + const uint32_t task_work_group_invocations = + Platform::GetSelectedPhysicalDevice()->mesh_shader_properties_ext.maxPreferredTaskWorkGroupInvocations; + + Platform::RecordCommandsMainQueue([&](const VkCommandBuffer vk_command_buffer) { + Platform::EverythingBarrier(vk_command_buffer); + for (int sub_step_i = 0; sub_step_i < num_sub_step; sub_step_i++) { + init_constraint_pipeline->Bind(vk_command_buffer); + init_constraint_pipeline->BindDescriptorSet( + vk_command_buffer, 0, + target_dynamic_strands.strands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + init_constraint_pipeline->BindDescriptorSet( + vk_command_buffer, 1, strands_physics_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + + init_constraint_pipeline->PushConstant(vk_command_buffer, 0, init_push_constant); + vkCmdDispatch(vk_command_buffer, Platform::DivUp(init_push_constant.constraint_size, task_work_group_invocations), + 1, 1); + Platform::EverythingBarrier(vk_command_buffer); + for (int iteration_i = 0; iteration_i < num_iteration; iteration_i++) { + update_constraint_pipeline->Bind(vk_command_buffer); + update_constraint_pipeline->BindDescriptorSet( + vk_command_buffer, 0, + target_dynamic_strands.strands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + update_constraint_pipeline->BindDescriptorSet( + vk_command_buffer, 1, strands_physics_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + + update_constraint_pipeline->PushConstant(vk_command_buffer, 0, update_constraint_constant); + vkCmdDispatch(vk_command_buffer, + Platform::DivUp(update_constraint_constant.constraint_size, task_work_group_invocations), 1, 1); + Platform::EverythingBarrier(vk_command_buffer); + + project_constraint_pipeline->Bind(vk_command_buffer); + project_constraint_pipeline->BindDescriptorSet( + vk_command_buffer, 0, + target_dynamic_strands.strands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + project_constraint_pipeline->BindDescriptorSet( + vk_command_buffer, 1, strands_physics_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); + + project_constraint_pipeline->PushConstant(vk_command_buffer, 0, project_constraint_constant); + vkCmdDispatch(vk_command_buffer, + Platform::DivUp(project_constraint_constant.strand_size, task_work_group_invocations), 1, 1); + Platform::EverythingBarrier(vk_command_buffer); + } + } + }); +} + +void DsStiffRod::DownloadData() { + rod_constraints_buffer->DownloadVector(rod_constraints, rod_constraints.size()); + per_strand_data_list_buffer->DownloadVector(per_strand_data_list, per_strand_data_list.size()); +} + +glm::vec3 DsStiffRod::ComputeDarbouxVector(const glm::quat& q0, const glm::quat& q1, + const float average_segment_length) { + const auto relative_rotation = glm::conjugate(q0) * q1; + return 2.f / average_segment_length * glm::vec3(relative_rotation.x, relative_rotation.y, relative_rotation.z); +} diff --git a/EvoEngine_Plugins/EcoSysLab/src/DynamicStrandsRender.cpp b/EvoEngine_Plugins/EcoSysLab/src/DynamicStrandsRender.cpp index 28ee5ab..e96efbe 100644 --- a/EvoEngine_Plugins/EcoSysLab/src/DynamicStrandsRender.cpp +++ b/EvoEngine_Plugins/EcoSysLab/src/DynamicStrandsRender.cpp @@ -2,7 +2,7 @@ #include "Shader.hpp" using namespace eco_sys_lab_plugin; -void DynamicStrands::Render(const StepParameters::RenderParameters& render_parameters) const { +void DynamicStrands::Render(const RenderParameters& render_parameters) const { if (!Platform::Constants::support_mesh_shader) { EVOENGINE_LOG("Failed to render! Mesh shader unsupported!") return; @@ -12,10 +12,9 @@ void DynamicStrands::Render(const StepParameters::RenderParameters& render_param EVOENGINE_LOG("Failed to render! RenderLayer not present!") return; } - const auto current_frame_index = current_left ? 0 : 1; + const auto current_frame_index = Platform::GetCurrentFrameIndex(); static std::shared_ptr render_pipeline{}; - static std::vector> strands_descriptor_sets{}; struct RenderPushConstant { uint32_t camera_index = 0; uint32_t strand_size; @@ -26,12 +25,6 @@ void DynamicStrands::Render(const StepParameters::RenderParameters& render_param static std::shared_ptr task_shader{}; static std::shared_ptr mesh_shader{}; static std::shared_ptr frag_shader{}; - const auto max_frame_in_flight = Platform::GetMaxFramesInFlight(); - strands_descriptor_sets.resize(max_frame_in_flight); - for (auto& i : strands_descriptor_sets) { - i = std::make_shared(strands_layout); - } - // Load shader task_shader = std::make_shared(); task_shader->Set( @@ -71,17 +64,12 @@ void DynamicStrands::Render(const StepParameters::RenderParameters& render_param render_pipeline->Initialize(); } - if (ref_strand_segments.empty()) - return; - - BindStrandsDescriptorSet(strands_descriptor_sets[current_frame_index]); - const uint32_t task_work_group_invocations = Platform::GetSelectedPhysicalDevice()->mesh_shader_properties_ext.maxPreferredTaskWorkGroupInvocations; RenderPushConstant push_constant; push_constant.camera_index = render_layer->GetCameraIndex(render_parameters.target_camera->GetHandle()); - push_constant.strand_size = ref_strands.size(); - push_constant.strand_segment_size = ref_strand_segments.size(); + push_constant.strand_size = strands.size(); + push_constant.strand_segment_size = strand_segments.size(); Platform::RecordCommandsMainQueue([&](const VkCommandBuffer vk_command_buffer) { #pragma region Viewport and scissor @@ -111,7 +99,7 @@ void DynamicStrands::Render(const StepParameters::RenderParameters& render_param render_pipeline->BindDescriptorSet(vk_command_buffer, 1, strands_descriptor_sets[current_frame_index]->GetVkDescriptorSet()); render_pipeline->PushConstant(vk_command_buffer, 0, push_constant); - const uint32_t count = Platform::DivUp(ref_strand_segments.size(), task_work_group_invocations); + const uint32_t count = Platform::DivUp(strand_segments.size(), task_work_group_invocations); vkCmdDrawMeshTasksEXT(vk_command_buffer, count, 1, 1); }); }); diff --git a/EvoEngine_SDK/Internals/DefaultResources/Shaders/Includes/Math.glsl b/EvoEngine_SDK/Internals/DefaultResources/Shaders/Includes/Math.glsl new file mode 100644 index 0000000..affbb9c --- /dev/null +++ b/EvoEngine_SDK/Internals/DefaultResources/Shaders/Includes/Math.glsl @@ -0,0 +1,126 @@ + + +mat4 translate(in vec3 d) { + return mat4(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, d.x, d.y, d.z, 1); +} + +mat4 scale(in vec3 s) { + return mat4(s.x, 0, 0, 0, 0, s.y, 0, 0, 0, 0, s.z, 0, 0, 0, 0, 1); +} + +mat4 mat4_cast(in vec4 q) { + float qxx = q.x * q.x; + float qyy = q.y * q.y; + float qzz = q.z * q.z; + float qxz = q.x * q.z; + float qxy = q.x * q.y; + float qyz = q.y * q.z; + float qwx = q.w * q.x; + float qwy = q.w * q.y; + float qwz = q.w * q.z; + mat4 ret_val = mat4(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1); + ret_val[0][0] = 1.0 - 2.0 * (qyy + qzz); + ret_val[0][1] = 2.0 * (qxy + qwz); + ret_val[0][2] = 2.0 * (qxz - qwy); + + ret_val[1][0] = 2.0 * (qxy - qwz); + ret_val[1][1] = 1.0 - 2.0 * (qxx + qzz); + ret_val[1][2] = 2.0 * (qyz + qwx); + + ret_val[2][0] = 2.0 * (qxz + qwy); + ret_val[2][1] = 2.0 * (qyz - qwx); + ret_val[2][2] = 1.0 - 2.0 * (qxx + qyy); + + return ret_val; +} + +mat3 mat3_cast(in vec4 q) { + float qxx = q.x * q.x; + float qyy = q.y * q.y; + float qzz = q.z * q.z; + float qxz = q.x * q.z; + float qxy = q.x * q.y; + float qyz = q.y * q.z; + float qwx = q.w * q.x; + float qwy = q.w * q.y; + float qwz = q.w * q.z; + mat3 ret_val; + ret_val[0][0] = 1.0 - 2.0 * (qyy + qzz); + ret_val[0][1] = 2.0 * (qxy + qwz); + ret_val[0][2] = 2.0 * (qxz - qwy); + + ret_val[1][0] = 2.0 * (qxy - qwz); + ret_val[1][1] = 1.0 - 2.0 * (qxx + qzz); + ret_val[1][2] = 2.0 * (qyz + qwx); + + ret_val[2][0] = 2.0 * (qxz + qwy); + ret_val[2][1] = 2.0 * (qyz - qwx); + ret_val[2][2] = 1.0 - 2.0 * (qxx + qyy); + + return ret_val; +} + +vec4 quat_cast(in mat3 m) { + float fourXSquaredMinus1 = m[0][0] - m[1][1] - m[2][2]; + float fourYSquaredMinus1 = m[1][1] - m[0][0] - m[2][2]; + float fourZSquaredMinus1 = m[2][2] - m[0][0] - m[1][1]; + float fourWSquaredMinus1 = m[0][0] + m[1][1] + m[2][2]; + + int biggestIndex = 0; + float fourBiggestSquaredMinus1 = fourWSquaredMinus1; + if (fourXSquaredMinus1 > fourBiggestSquaredMinus1) { + fourBiggestSquaredMinus1 = fourXSquaredMinus1; + biggestIndex = 1; + } + if (fourYSquaredMinus1 > fourBiggestSquaredMinus1) { + fourBiggestSquaredMinus1 = fourYSquaredMinus1; + biggestIndex = 2; + } + if (fourZSquaredMinus1 > fourBiggestSquaredMinus1) { + fourBiggestSquaredMinus1 = fourZSquaredMinus1; + biggestIndex = 3; + } + + float biggestVal = sqrt(fourBiggestSquaredMinus1 + 1.0) * 0.5; + float mult = 0.25 / biggestVal; + + switch (biggestIndex) { + case 0: + return vec4((m[1][2] - m[2][1]) * mult, (m[2][0] - m[0][2]) * mult, (m[0][1] - m[1][0]) * mult, biggestVal); + case 1: + return vec4(biggestVal, (m[0][1] + m[1][0]) * mult, (m[2][0] + m[0][2]) * mult, (m[1][2] - m[2][1]) * mult); + case 2: + return vec4((m[0][1] + m[1][0]) * mult, biggestVal, (m[1][2] + m[2][1]) * mult, (m[2][0] - m[0][2]) * mult); + case 3: + return vec4((m[2][0] + m[0][2]) * mult, (m[1][2] + m[2][1]) * mult, biggestVal, (m[0][1] - m[1][0]) * mult); + default: + return vec4(0, 0, 0, 1); + } +} + +vec4 quatLookAt(in vec3 direction, in vec3 up) { + mat3 Result; + + Result[2] = -direction; + vec3 Right = cross(up, Result[2]); + Result[0] = Right * inversesqrt(max(0.00001, dot(Right, Right))); + Result[1] = cross(Result[2], Result[0]); + + return quat_cast(Result); +} + +vec4 conjugate(in vec4 q) { + return vec4(-q.x, -q.y, -q.z, q.w); +} + +vec4 quat_mul(in vec4 a, in vec4 b) { + return vec4(a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y, + a.w * b.y + a.y * b.w + a.z * b.x - a.x * b.z, + a.w * b.z + a.z * b.w + a.x * b.y - a.y * b.x, + a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z); +} + +vec3 compute_darboux_vector(in vec4 q0, in vec4 q1, float length) { + return (quat_mul(conjugate(q0), q1) * 2.0 / length).xyz; // Eq. 27.5 +} + diff --git a/EvoEngine_SDK/include/Rendering/Platform/Platform.hpp b/EvoEngine_SDK/include/Rendering/Platform/Platform.hpp index 2d96f47..63cecb6 100644 --- a/EvoEngine_SDK/include/Rendering/Platform/Platform.hpp +++ b/EvoEngine_SDK/include/Rendering/Platform/Platform.hpp @@ -149,9 +149,7 @@ class Platform final { void RecreateSwapChain(); void OnDestroy(); - void SwapChainSwapImage(); void SubmitPresent(); - void WaitForCommandsComplete() const; void Submit(); void ResetCommandBuffers(); @@ -217,7 +215,6 @@ class Platform final { inline static uint32_t max_spot_light_size = 16; }; - class Constants { public: inline static bool support_mesh_shader = true; diff --git a/EvoEngine_SDK/src/Platform.cpp b/EvoEngine_SDK/src/Platform.cpp index beed025..b2c27aa 100644 --- a/EvoEngine_SDK/src/Platform.cpp +++ b/EvoEngine_SDK/src/Platform.cpp @@ -12,8 +12,8 @@ #include "WindowLayer.hpp" #include "vk_mem_alloc.h" -#ifndef NDEBUG -#define GRAPHICS_VALIDATION true +#ifndef NDEBUG +# define GRAPHICS_VALIDATION true #endif using namespace evo_engine; @@ -1453,33 +1453,6 @@ void Platform::OnDestroy() { #pragma endregion } -void Platform::SwapChainSwapImage() { - if (const auto& window_layer = Application::GetLayer(); - window_layer->window_size_.x == 0 || window_layer->window_size_.y == 0) - return; - const auto just_now = Times::Now(); - vkDeviceWaitIdle(vk_device_); - GeometryStorage::DeviceSync(); - TextureStorage::DeviceSync(); - const VkFence in_flight_fences[] = {in_flight_fences_[current_frame_index_]->GetVkFence()}; - vkWaitForFences(vk_device_, 1, in_flight_fences, VK_TRUE, UINT64_MAX); - cpu_wait_time = Times::Now() - just_now; - auto result = vkAcquireNextImageKHR(vk_device_, swapchain_->GetVkSwapchain(), UINT64_MAX, - image_available_semaphores_[current_frame_index_]->GetVkSemaphore(), - VK_NULL_HANDLE, &next_image_index_); - while (result == VK_ERROR_OUT_OF_DATE_KHR || result == VK_SUBOPTIMAL_KHR || recreate_swap_chain_) { - RecreateSwapChain(); - result = vkAcquireNextImageKHR(vk_device_, swapchain_->GetVkSwapchain(), UINT64_MAX, - image_available_semaphores_[current_frame_index_]->GetVkSemaphore(), VK_NULL_HANDLE, - &next_image_index_); - recreate_swap_chain_ = false; - } - if (result != VK_SUCCESS) { - throw std::runtime_error("failed to acquire swap chain image!"); - } - vkResetFences(vk_device_, 1, in_flight_fences); -} - void Platform::SubmitPresent() { if (const auto& window_layer = Application::GetLayer(); window_layer->window_size_.x == 0 || window_layer->window_size_.y == 0) @@ -1502,15 +1475,6 @@ void Platform::SubmitPresent() { current_frame_index_ = (current_frame_index_ + 1) % max_frame_in_flight_; } -void Platform::WaitForCommandsComplete() const { - vkDeviceWaitIdle(vk_device_); - GeometryStorage::DeviceSync(); - TextureStorage::DeviceSync(); - const VkFence in_flight_fences[] = {in_flight_fences_[current_frame_index_]->GetVkFence()}; - vkWaitForFences(vk_device_, 1, in_flight_fences, VK_TRUE, UINT64_MAX); - vkResetFences(vk_device_, 1, in_flight_fences); -} - void Platform::Submit() { main_queue_->Submit(command_buffer_pool_[current_frame_index_], 0, used_command_buffer_size_, {}, {}, in_flight_fences_[current_frame_index_]); @@ -1527,7 +1491,6 @@ void Platform::ResetCommandBuffers() { #pragma endregion - void Platform::PostResourceLoadingInitialization() { const auto& graphics = GetInstance(); PrepareDescriptorSetLayouts(); @@ -1543,15 +1506,46 @@ void Platform::PreUpdate() { auto& graphics = GetInstance(); const auto window_layer = Application::GetLayer(); const auto render_layer = Application::GetLayer(); + + const auto vulkan_update = [&](const std::function& action) { + vkDeviceWaitIdle(graphics.vk_device_); + GeometryStorage::DeviceSync(); + TextureStorage::DeviceSync(); + const VkFence in_flight_fences[] = {graphics.in_flight_fences_[graphics.current_frame_index_]->GetVkFence()}; + vkWaitForFences(graphics.vk_device_, 1, in_flight_fences, VK_TRUE, UINT64_MAX); + action(); + vkResetFences(graphics.vk_device_, 1, in_flight_fences); + }; + if (window_layer) { if (glfwWindowShouldClose(window_layer->window_)) { Application::End(); } if (render_layer || Application::GetLayer()) { - graphics.SwapChainSwapImage(); + if (window_layer->window_size_.x != 0 || window_layer->window_size_.y != 0) { + const auto just_now = Times::Now(); + vulkan_update([&]() { + graphics.cpu_wait_time = Times::Now() - just_now; + auto result = vkAcquireNextImageKHR( + graphics.vk_device_, graphics.swapchain_->GetVkSwapchain(), UINT64_MAX, + graphics.image_available_semaphores_[graphics.current_frame_index_]->GetVkSemaphore(), VK_NULL_HANDLE, + &graphics.next_image_index_); + while (result == VK_ERROR_OUT_OF_DATE_KHR || result == VK_SUBOPTIMAL_KHR || graphics.recreate_swap_chain_) { + graphics.RecreateSwapChain(); + result = vkAcquireNextImageKHR( + graphics.vk_device_, graphics.swapchain_->GetVkSwapchain(), UINT64_MAX, + graphics.image_available_semaphores_[graphics.current_frame_index_]->GetVkSemaphore(), VK_NULL_HANDLE, + &graphics.next_image_index_); + graphics.recreate_swap_chain_ = false; + } + if (result != VK_SUCCESS) { + throw std::runtime_error("failed to acquire swap chain image!"); + } + }); + } } } else { - graphics.WaitForCommandsComplete(); + vulkan_update({}); } graphics.ResetCommandBuffers();