diff --git a/examples/testbed/framework/view.cpp b/examples/testbed/framework/view.cpp index 7173102..7229916 100644 --- a/examples/testbed/framework/view.cpp +++ b/examples/testbed/framework/view.cpp @@ -53,7 +53,7 @@ static inline bool ImGui_GLFW_GL_Init(GLFWwindow* w, bool install_callbacks) // error #endif - + return false; } static inline void ImGui_GLFW_GL_Shutdown() diff --git a/include/bounce/dynamics/cloth/cloth.h b/include/bounce/dynamics/cloth/cloth.h index 0dce055..9897b6e 100644 --- a/include/bounce/dynamics/cloth/cloth.h +++ b/include/bounce/dynamics/cloth/cloth.h @@ -147,6 +147,8 @@ private: // Proxy mesh b3ClothMesh* m_mesh; + + // Mesh density float32 m_density; // Particle pool diff --git a/include/bounce/dynamics/cloth/cloth_mesh.h b/include/bounce/dynamics/cloth/cloth_mesh.h index f25273f..8aea759 100644 --- a/include/bounce/dynamics/cloth/cloth_mesh.h +++ b/include/bounce/dynamics/cloth/cloth_mesh.h @@ -22,6 +22,8 @@ #include #include +class b3Particle; + struct b3ClothMeshTriangle { u32 v1, v2, v3; @@ -41,8 +43,6 @@ struct b3ClothMeshSewingLine u32 v1, v2; }; -class b3Particle; - struct b3ClothMesh { u32 vertexCount; diff --git a/include/bounce/dynamics/cloth/cloth_solver.h b/include/bounce/dynamics/cloth/cloth_solver.h index deb91d0..90007b9 100644 --- a/include/bounce/dynamics/cloth/cloth_solver.h +++ b/include/bounce/dynamics/cloth/cloth_solver.h @@ -24,15 +24,14 @@ class b3StackAllocator; -struct b3DenseVec3; -struct b3DiagMat33; -struct b3SparseSymMat33; -struct b3SymMat33; - class b3Particle; class b3Force; struct b3BodyContact; +struct b3DenseVec3; +struct b3DiagMat33; +struct b3SparseSymMat33; + struct b3ClothSolverDef { b3StackAllocator* stack; @@ -43,11 +42,12 @@ struct b3ClothSolverDef struct b3ClothSolverData { - b3Vec3* x; - b3Vec3* v; - b3Vec3* f; - b3SymMat33* dfdx; - b3SymMat33* dfdv; + b3DenseVec3* x; + b3DenseVec3* v; + b3DenseVec3* f; + b3DenseVec3* y; + b3SparseSymMat33* dfdx; + b3SparseSymMat33* dfdv; float32 dt; float32 invdt; }; @@ -59,45 +59,6 @@ struct b3AccelerationConstraint b3Vec3 p, q, z; }; -struct b3SymMat33 -{ - b3SymMat33(b3StackAllocator* a, u32 m, u32 n) - { - allocator = a; - M = m; - N = n; - values = (b3Mat33*)b3Alloc(M * N * sizeof(b3Mat33)); - } - - ~b3SymMat33() - { - b3Free(values); - } - - b3Mat33& operator()(u32 i, u32 j) - { - return values[i * N + j]; - } - - const b3Mat33& operator()(u32 i, u32 j) const - { - return values[i * N + j]; - } - - void SetZero() - { - for (u32 v = 0; v < M * N; ++v) - { - values[v].SetZero(); - } - } - - u32 M; - u32 N; - b3Mat33* values; - b3StackAllocator* allocator; -}; - class b3ClothSolver { public: @@ -120,7 +81,7 @@ private: void InitializeConstraints(); // Compute A and b in Ax = b - void Compute_A_b(b3SparseSymMat33& A, b3DenseVec3& b, const b3DenseVec3& f, const b3DenseVec3& x, const b3DenseVec3& v, const b3DenseVec3& y) const; + void Compute_A_b(b3SparseSymMat33& A, b3DenseVec3& b) const; // Compute S and z. void Compute_S_z(b3DiagMat33& S, b3DenseVec3& z); diff --git a/include/bounce/dynamics/cloth/sparse_sym_mat33.h b/include/bounce/dynamics/cloth/sparse_sym_mat33.h index ce21e77..a9fa45b 100644 --- a/include/bounce/dynamics/cloth/sparse_sym_mat33.h +++ b/include/bounce/dynamics/cloth/sparse_sym_mat33.h @@ -19,7 +19,6 @@ #ifndef B3_SPARSE_SYM_MAT_33_H #define B3_SPARSE_SYM_MAT_33_H -#include #include #include #include @@ -27,7 +26,7 @@ struct b3SparseSymMat33 { // - b3SparseSymMat33(b3StackAllocator* a, u32 m, u32 n); + b3SparseSymMat33(u32 m, u32 n); // ~b3SparseSymMat33(); @@ -42,8 +41,8 @@ struct b3SparseSymMat33 void Diagonal(b3DiagMat33& out) const; // + void SetZero(); - b3StackAllocator* allocator; u32 M; u32 N; u32* row_ptrs; @@ -53,25 +52,24 @@ struct b3SparseSymMat33 u32* value_columns; }; -inline b3SparseSymMat33::b3SparseSymMat33(b3StackAllocator* a, u32 m, u32 n) +inline b3SparseSymMat33::b3SparseSymMat33(u32 m, u32 n) { B3_ASSERT(m == n); - allocator = a; M = m; N = n; - row_ptrs = (u32*)allocator->Allocate((M + 1) * sizeof(u32)); + row_ptrs = (u32*)b3Alloc((M + 1) * sizeof(u32)); memset(row_ptrs, 0, (M + 1) * sizeof(u32)); value_count = 0; - value_capacity = n * (n + 1) / 2; - values = (b3Mat33*)allocator->Allocate(value_capacity * sizeof(b3Mat33)); - value_columns = (u32*)allocator->Allocate(value_capacity * sizeof(u32)); + value_capacity = M * (M + 1) / 2; + values = (b3Mat33*)b3Alloc(value_capacity * sizeof(b3Mat33)); + value_columns = (u32*)b3Alloc(value_capacity * sizeof(u32)); } inline b3SparseSymMat33::~b3SparseSymMat33() { - allocator->Free(value_columns); - allocator->Free(values); - allocator->Free(row_ptrs); + b3Free(value_columns); + b3Free(values); + b3Free(row_ptrs); } inline const b3Mat33& b3SparseSymMat33::operator()(u32 i, u32 j) const @@ -103,7 +101,7 @@ inline const b3Mat33& b3SparseSymMat33::operator()(u32 i, u32 j) const } inline b3Mat33& b3SparseSymMat33::operator()(u32 i, u32 j) -{ +{ B3_ASSERT(i < M); B3_ASSERT(j < N); @@ -120,57 +118,44 @@ inline b3Mat33& b3SparseSymMat33::operator()(u32 i, u32 j) { u32 row_value_index = row_value_begin + row_value; u32 row_value_column = value_columns[row_value_index]; - + if (row_value_column == j) { return values[row_value_index]; } } - // Insert sorted by column - u32 max_column_row_value = 0; + // Find insert position + u32 row_value_k = 0; for (u32 row_value = 0; row_value < row_value_count; ++row_value) { u32 row_value_index = row_value_begin + row_value; u32 row_value_column = value_columns[row_value_index]; - + if (row_value_column >= j) { - max_column_row_value = row_value; + row_value_k = row_value; break; } } - u32 max_column_row_value_index = row_value_begin + max_column_row_value; + // Shift the values + u32 right_count = value_count - row_value_begin - row_value_k; + memcpy(value_columns + row_value_begin + row_value_k + 1, value_columns + row_value_begin + row_value_k, right_count * sizeof(u32)); + memcpy(values + row_value_begin + row_value_k + 1, values + row_value_begin + row_value_k, right_count * sizeof(b3Mat33)); - // Copy the values to be shifted - u32 shift_count = value_count - max_column_row_value_index; - - b3Mat33* shift_values = (b3Mat33*)allocator->Allocate(shift_count * sizeof(b3Mat33)); - memcpy(shift_values, values + max_column_row_value_index, shift_count * sizeof(b3Mat33)); - - u32* shift_value_columns = (u32*)allocator->Allocate(shift_count * sizeof(u32)); - memcpy(shift_value_columns, value_columns + max_column_row_value_index, shift_count * sizeof(u32)); - - // Insert the new value - B3_ASSERT(value_count < value_capacity); - value_columns[max_column_row_value_index] = j; - ++value_count; - - // Shift the old values - memcpy(values + max_column_row_value_index + 1, shift_values, shift_count * sizeof(b3Mat33)); - memcpy(value_columns + max_column_row_value_index + 1, shift_value_columns, shift_count * sizeof(u32)); - - allocator->Free(shift_value_columns); - allocator->Free(shift_values); - - // Shift the old row pointers as well + // Shift the row pointers for (u32 row_ptr_index = i + 1; row_ptr_index < M + 1; ++row_ptr_index) { ++row_ptrs[row_ptr_index]; } - return values[max_column_row_value_index]; + // Insert the value + value_columns[row_value_begin + row_value_k] = j; + ++value_count; + + // Return the value + return values[row_value_begin + row_value_k]; } inline void b3SparseSymMat33::Diagonal(b3DiagMat33& out) const @@ -198,6 +183,17 @@ inline void b3SparseSymMat33::Diagonal(b3DiagMat33& out) const } } +inline void b3SparseSymMat33::SetZero() +{ + for (u32 i = 0; i < M; ++i) + { + for (u32 j = i; j < M; ++j) + { + (*this)(i, j).SetZero(); + } + } +} + inline void b3Mul(b3DenseVec3& out, const b3SparseSymMat33& A, const b3DenseVec3& v) { B3_ASSERT(A.N == out.n); diff --git a/src/bounce/dynamics/cloth/cloth_solver.cpp b/src/bounce/dynamics/cloth/cloth_solver.cpp index 224b132..328f4c6 100644 --- a/src/bounce/dynamics/cloth/cloth_solver.cpp +++ b/src/bounce/dynamics/cloth/cloth_solver.cpp @@ -25,29 +25,6 @@ #include #include -static B3_FORCE_INLINE void b3Mul(b3DenseVec3& out, const b3SymMat33& A, const b3DenseVec3& v) -{ - B3_ASSERT(A.N == v.n); - B3_ASSERT(out.n == A.M); - - for (u32 row = 0; row < A.M; ++row) - { - out[row].SetZero(); - - for (u32 column = 0; column < A.N; ++column) - { - out[row] += A(row, column) * v[column]; - } - } -} - -static B3_FORCE_INLINE b3DenseVec3 operator*(const b3SymMat33& m, const b3DenseVec3& v) -{ - b3DenseVec3 result(v.n); - b3Mul(result, m, v); - return result; -} - // Here, we solve Ax = b using the Modified Preconditioned Conjugate Gradient (MPCG) algorithm. // described in the paper: // "Large Steps in Cloth Simulation - David Baraff, Andrew Witkin". @@ -168,23 +145,22 @@ void b3ClothSolver::InitializeConstraints() void b3ClothSolver::Solve(float32 dt, const b3Vec3& gravity) { - B3_PROFILE("Integrate"); - b3DenseVec3 sx(m_particleCount); b3DenseVec3 sv(m_particleCount); b3DenseVec3 sf(m_particleCount); b3DenseVec3 sy(m_particleCount); b3DenseVec3 sx0(m_particleCount); - b3SymMat33 dfdx(m_allocator, m_particleCount, m_particleCount); + b3SparseSymMat33 dfdx(m_particleCount, m_particleCount); dfdx.SetZero(); - b3SymMat33 dfdv(m_allocator, m_particleCount, m_particleCount); + b3SparseSymMat33 dfdv(m_particleCount, m_particleCount); dfdv.SetZero(); - m_solverData.x = sx.v; - m_solverData.v = sv.v; - m_solverData.f = sf.v; + m_solverData.x = &sx; + m_solverData.v = &sv; + m_solverData.f = &sf; + m_solverData.y = &sy; m_solverData.dfdx = &dfdx; m_solverData.dfdv = &dfdv; m_solverData.dt = dt; @@ -235,12 +211,12 @@ void b3ClothSolver::Solve(float32 dt, const b3Vec3& gravity) // b = h * (f0 + h * dfdx * v0 + dfdx * y) // A - b3SparseSymMat33 A(m_allocator, m_particleCount, m_particleCount); + b3SparseSymMat33 A(m_particleCount, m_particleCount); // b b3DenseVec3 b(m_particleCount); - Compute_A_b(A, b, sf, sx, sv, sy); + Compute_A_b(A, b); // x b3DenseVec3 x(m_particleCount); @@ -292,11 +268,17 @@ void b3ClothSolver::Solve(float32 dt, const b3Vec3& gravity) } } -void b3ClothSolver::Compute_A_b(b3SparseSymMat33& A, b3DenseVec3& b, const b3DenseVec3& f, const b3DenseVec3& x, const b3DenseVec3& v, const b3DenseVec3& y) const +void b3ClothSolver::Compute_A_b(b3SparseSymMat33& A, b3DenseVec3& b) const { + B3_PROFILE("Compute A, b"); + + b3DenseVec3& x = *m_solverData.x; + b3DenseVec3& v = *m_solverData.v; + b3DenseVec3& f = *m_solverData.f; + b3DenseVec3& y = *m_solverData.y; + b3SparseSymMat33& dfdx = *m_solverData.dfdx; + b3SparseSymMat33& dfdv = *m_solverData.dfdv; float32 h = m_solverData.dt; - b3SymMat33& dfdx = *m_solverData.dfdx; - b3SymMat33& dfdv = *m_solverData.dfdv; // Compute A // A = M - h * dfdv - h * h * dfdx @@ -439,8 +421,6 @@ void b3ClothSolver::Solve(b3DenseVec3& x, u32& iterations, break; } - B3_ASSERT(b3IsValid(delta_new)); - // Convergence check. if (delta_new <= tolerance * tolerance * b_delta) { diff --git a/src/bounce/dynamics/cloth/spring_force.cpp b/src/bounce/dynamics/cloth/spring_force.cpp index bdc91b3..400e164 100644 --- a/src/bounce/dynamics/cloth/spring_force.cpp +++ b/src/bounce/dynamics/cloth/spring_force.cpp @@ -19,6 +19,8 @@ #include #include #include +#include +#include void b3SpringForceDef::Initialize(b3Particle* particle1, b3Particle* particle2, float32 structuralStiffness, float32 dampingStiffness) { @@ -50,14 +52,17 @@ b3SpringForce::~b3SpringForce() void b3SpringForce::Initialize(const b3ClothSolverData* data) { + b3DenseVec3& x = *data->x; + b3DenseVec3& v = *data->v; + u32 i1 = m_p1->m_solverId; u32 i2 = m_p2->m_solverId; - b3Vec3 x1 = data->x[i1]; - b3Vec3 v1 = data->v[i1]; + b3Vec3 x1 = x[i1]; + b3Vec3 v1 = v[i1]; - b3Vec3 x2 = data->x[i2]; - b3Vec3 v2 = data->v[i2]; + b3Vec3 x2 = x[i2]; + b3Vec3 v2 = v[i2]; b3Mat33 I; I.SetIdentity(); @@ -90,9 +95,9 @@ void b3SpringForce::Initialize(const b3ClothSolverData* data) void b3SpringForce::Apply(const b3ClothSolverData* data) { - b3Vec3* f = data->f; - b3SymMat33& dfdx = *data->dfdx; - b3SymMat33& dfdv = *data->dfdv; + b3DenseVec3& f = *data->f; + b3SparseSymMat33& dfdx = *data->dfdx; + b3SparseSymMat33& dfdv = *data->dfdv; u32 i1 = m_p1->m_solverId; u32 i2 = m_p2->m_solverId; @@ -102,21 +107,21 @@ void b3SpringForce::Apply(const b3ClothSolverData* data) b3Mat33 Jx11 = m_Jx; b3Mat33 Jx12 = -Jx11; - b3Mat33 Jx21 = Jx12; + //b3Mat33 Jx21 = Jx12; b3Mat33 Jx22 = Jx11; dfdx(i1, i1) += Jx11; dfdx(i1, i2) += Jx12; - dfdx(i2, i1) += Jx21; + //dfdx(i2, i1) += Jx21; dfdx(i2, i2) += Jx22; b3Mat33 Jv11 = m_Jv; b3Mat33 Jv12 = -Jv11; - b3Mat33 Jv21 = Jv12; + //b3Mat33 Jv21 = Jv12; b3Mat33 Jv22 = Jv11; dfdv(i1, i1) += Jv11; dfdv(i1, i2) += Jv12; - dfdv(i2, i1) += Jv21; + //dfdv(i2, i1) += Jv21; dfdv(i2, i2) += Jv22; } \ No newline at end of file