diff --git a/CMakeLists.txt b/CMakeLists.txt index 1b0e26b..26ba37c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,10 +1,17 @@ cmake_minimum_required(VERSION 3.10.0) +if(POLICY CMP0072) + cmake_policy(SET CMP0072 NEW) +endif() + project(fast-mass-spring) +set(OpenGL_GL_PREFERENCE GLVND) + set(Sources ClothApp/app.cpp ClothApp/MassSpringSolver.cpp + ClothApp/PBDSolver.cpp ClothApp/Mesh.cpp ClothApp/Renderer.cpp ClothApp/Shader.cpp @@ -27,6 +34,7 @@ FetchContent_Declare( FetchContent_Declare( eigen GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git + GIT_TAG 3.4.0 ) FetchContent_Declare( @@ -43,9 +51,11 @@ endif() FetchContent_GetProperties(eigen) if(NOT eigen_POPULATED) FetchContent_Populate(eigen) - add_subdirectory(${eigen_SOURCE_DIR} ${eigen_BINARY_DIR}) endif() +add_library(eigen INTERFACE) +target_include_directories(eigen INTERFACE ${eigen_SOURCE_DIR}) + FetchContent_GetProperties(glm) if(NOT glm_POPULATED) @@ -58,9 +68,15 @@ if(WIN32) add_definitions(-D_USE_MATH_DEFINES) endif() -# copy shaders to binary directory -file(INSTALL ClothApp/shaders/ DESTINATION shaders/) - # create executable add_executable(fast-mass-spring ${Sources}) target_link_libraries(fast-mass-spring ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES} ${GLEW_LIBRARIES} OpenMeshCore eigen glm) + +# copy shaders to binary directory on every build. +add_custom_target(copy_shaders ALL + COMMAND ${CMAKE_COMMAND} -E copy_directory + ${CMAKE_CURRENT_SOURCE_DIR}/ClothApp/shaders + ${CMAKE_CURRENT_BINARY_DIR}/shaders +) + +add_dependencies(fast-mass-spring copy_shaders) diff --git a/ClothApp/MassSpringSolver.cpp b/ClothApp/MassSpringSolver.cpp index 83e314a..d847484 100644 --- a/ClothApp/MassSpringSolver.cpp +++ b/ClothApp/MassSpringSolver.cpp @@ -1,5 +1,6 @@ #include "MassSpringSolver.h" #include +#include // S Y S T E M ////////////////////////////////////////////////////////////////////////////////////// mass_spring_system::mass_spring_system( diff --git a/ClothApp/MassSpringSolver.h b/ClothApp/MassSpringSolver.h index d01a604..aee51d0 100644 --- a/ClothApp/MassSpringSolver.h +++ b/ClothApp/MassSpringSolver.h @@ -112,6 +112,13 @@ class MassSpringBuilder { // Constraint Graph class CgNodeVisitor; // Constraint graph node visitor +class FixedPointController { +public: + virtual ~FixedPointController() = default; + virtual void fixPoint(unsigned int i) = 0; + virtual void releasePoint(unsigned int i) = 0; +}; + // Constraint graph node class CgNode { protected: @@ -158,7 +165,7 @@ class CgRootNode : public CgSpringNode { virtual bool accept(CgNodeVisitor& visitor); }; -class CgPointFixNode : public CgPointNode { +class CgPointFixNode : public CgPointNode, public FixedPointController { protected: typedef Eigen::Vector3f Vector3f; std::unordered_map fix_map; diff --git a/ClothApp/Mesh.cpp b/ClothApp/Mesh.cpp index 1417623..9abef74 100644 --- a/ClothApp/Mesh.cpp +++ b/ClothApp/Mesh.cpp @@ -1,10 +1,10 @@ #include "Mesh.h" // M E S H ///////////////////////////////////////////////////////////////////////////////////// -float* Mesh::vbuff() { return VERTEX_DATA(this); } -float* Mesh::nbuff() { return NORMAL_DATA(this); } -float* Mesh::tbuff() { return TEXTURE_DATA(this); } -unsigned int* Mesh::ibuff() { return &_ibuff[0]; } +float* Mesh::vbuff() { return VERTEX_DATA(this); } // vertex buffer +float* Mesh::nbuff() { return NORMAL_DATA(this); } // normal buffer +float* Mesh::tbuff() { return TEXTURE_DATA(this); } // texture coordinate buffer +unsigned int* Mesh::ibuff() { return &_ibuff[0]; } // index buffer void Mesh::useIBuff(std::vector& _ibuff) { this->_ibuff = _ibuff; } unsigned int Mesh::vbuffLen() { return (unsigned int)n_vertices() * 3; } @@ -35,7 +35,7 @@ void MeshBuilder::uniformGrid(float w, int n) { for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { - handle_table[j + i * n] = result->add_vertex(o + d*j*ux + d*i*uy); // add vertex + handle_table[j + i * n] = result->add_vertex(o + d*j*ux + d*i*uy); // add vertex, vertexId(i, j) = i * n + j result->set_texcoord2D(handle_table[j + i * n], OpenMesh::Vec2f(ud*j, ud*i)); // add texture coordinates //add connectivity diff --git a/ClothApp/PBDSolver.cpp b/ClothApp/PBDSolver.cpp new file mode 100644 index 0000000..a7c6e90 --- /dev/null +++ b/ClothApp/PBDSolver.cpp @@ -0,0 +1,2148 @@ +#include "PBDSolver.h" +#include +#include +#include +#include +#include +#include +#include + +namespace PBDDefaultParam { + static const unsigned int solverIterations = 15; + static const float dampingFactor = 0.07f; + static const float collisionEps = 1e-4f; + static const float collisionStiffness = 1.0f; + static const float contactFriction = 0.15f; + static const float selfCollisionStiffness = 0.3f; + static const unsigned int maxSelfCollisionContactsPerVertex = 4u; // 12u + static const float velocitySleepThreshold = 5e-3f; + static const Eigen::Vector3f gravity(0.0f, 0.0f, -9.81f); + static const float twoPi = 6.28318530718f; +} + +namespace { +// Apply the iteration-corrected stiffness +// k' = 1 - (1 - k)^(1 / n_s) +// so that the effective stiffness stays consistent when the number of solver +// iterations n_s changes. +float correctedStiffness(float stiffness, unsigned int solverIterations) { + if (stiffness <= 0.0f) return 0.0f; + if (stiffness >= 1.0f) return 1.0f; + if (solverIterations == 0u) return stiffness; + return 1.0f - std::pow(1.0f - stiffness, 1.0f / static_cast(solverIterations)); +} + +struct SharedEdgeKey { + unsigned int a; + unsigned int b; + + bool operator==(const SharedEdgeKey& other) const { + return a == other.a && b == other.b; + } +}; + +struct SharedEdgeKeyHash { + std::size_t operator()(const SharedEdgeKey& key) const { + return (static_cast(key.a) << 32) ^ static_cast(key.b); + } +}; + +struct PendingTriangleEdge { + unsigned int edge0; + unsigned int edge1; + unsigned int opposite; +}; + +struct SpatialHashKey { + int x; + int y; + int z; + + bool operator==(const SpatialHashKey& other) const { + return x == other.x && y == other.y && z == other.z; + } +}; + +struct SpatialHashKeyHash { + std::size_t operator()(const SpatialHashKey& key) const { + const std::size_t hx = std::hash{}(key.x); + const std::size_t hy = std::hash{}(key.y); + const std::size_t hz = std::hash{}(key.z); + return hx ^ (hy << 1) ^ (hz << 2); + } +}; + +SpatialHashKey hashPosition(const Eigen::Vector3f& position, float cellSize) { + assert(cellSize > 0.0f); + return SpatialHashKey{ + static_cast(std::floor(position.x() / cellSize)), + static_cast(std::floor(position.y() / cellSize)), + static_cast(std::floor(position.z() / cellSize)) + }; +} + +void insertSweptVertexIntoHash( + const Eigen::Vector3f& previousPosition, + const Eigen::Vector3f& predictedPosition, + float padding, + float cellSize, + unsigned int vertexIndex, + std::unordered_map, SpatialHashKeyHash>& verticesByCell +) { + const Eigen::Vector3f minCorner = previousPosition.cwiseMin(predictedPosition) + - Eigen::Vector3f::Constant(padding); + const Eigen::Vector3f maxCorner = previousPosition.cwiseMax(predictedPosition) + + Eigen::Vector3f::Constant(padding); + + const SpatialHashKey minCell = hashPosition(minCorner, cellSize); + const SpatialHashKey maxCell = hashPosition(maxCorner, cellSize); + for (int cellX = minCell.x; cellX <= maxCell.x; ++cellX) { + for (int cellY = minCell.y; cellY <= maxCell.y; ++cellY) { + for (int cellZ = minCell.z; cellZ <= maxCell.z; ++cellZ) { + verticesByCell[SpatialHashKey{ cellX, cellY, cellZ }].push_back(vertexIndex); + } + } + } +} + +bool barycentricCoordinates( + const Eigen::Vector3f& point, + const Eigen::Vector3f& a, + const Eigen::Vector3f& b, + const Eigen::Vector3f& c, + Eigen::Vector3f& barycentric +) { + const Eigen::Vector3f v0 = b - a; + const Eigen::Vector3f v1 = c - a; + const Eigen::Vector3f v2 = point - a; + + const float d00 = v0.dot(v0); + const float d01 = v0.dot(v1); + const float d11 = v1.dot(v1); + const float d20 = v2.dot(v0); + const float d21 = v2.dot(v1); + const float denominator = d00 * d11 - d01 * d01; + if (std::abs(denominator) <= 1e-8f) { + return false; + } + + const float w1 = (d11 * d20 - d01 * d21) / denominator; + const float w2 = (d00 * d21 - d01 * d20) / denominator; + const float w0 = 1.0f - w1 - w2; + barycentric = Eigen::Vector3f(w0, w1, w2); + return true; +} + +bool pointProjectsInsideTriangle( + const Eigen::Vector3f& point, + const Eigen::Vector3f& a, + const Eigen::Vector3f& b, + const Eigen::Vector3f& c, + Eigen::Vector3f& barycentric +) { + if (!barycentricCoordinates(point, a, b, c, barycentric)) { + return false; + } + + const float tolerance = -1e-4f; + return barycentric.x() >= tolerance + && barycentric.y() >= tolerance + && barycentric.z() >= tolerance; +} + +bool triangleUnitNormal( + const std::vector& positions, + unsigned int p1Index, + unsigned int p2Index, + unsigned int p3Index, + Eigen::Vector3f& outNormal, + float* outDoubleArea = nullptr +) { + if (p1Index >= positions.size() + || p2Index >= positions.size() + || p3Index >= positions.size()) { + return false; + } + + const Eigen::Vector3f& p1 = positions[p1Index]; + const Eigen::Vector3f& p2 = positions[p2Index]; + const Eigen::Vector3f& p3 = positions[p3Index]; + outNormal = (p2 - p1).cross(p3 - p1); + const float normalLength = outNormal.norm(); + if (normalLength <= 1e-8f) { + return false; + } + + if (outDoubleArea != nullptr) { + *outDoubleArea = normalLength; + } + outNormal /= normalLength; + return true; +} + +bool triangleNormalAndSignedDistance( + const std::vector& positions, + unsigned int vertexIndex, + unsigned int p1Index, + unsigned int p2Index, + unsigned int p3Index, + Eigen::Vector3f& outNormal, + float& outSignedDistance +) { + if (vertexIndex >= positions.size() + || p1Index >= positions.size() + || p2Index >= positions.size() + || p3Index >= positions.size()) { + return false; + } + + const Eigen::Vector3f& p1 = positions[p1Index]; + if (!triangleUnitNormal(positions, p1Index, p2Index, p3Index, outNormal)) { + return false; + } + outSignedDistance = (positions[vertexIndex] - p1).dot(outNormal); + return true; +} + +bool previousFrameSideAndNormal( + const std::vector& previousPositions, + const std::vector& referencePositions, + unsigned int vertexIndex, + unsigned int p1Index, + unsigned int p2Index, + unsigned int p3Index, + bool& flipNormal, + float& outSignedDistance +) { + Eigen::Vector3f normal; + + // Use the previous frame to decide which side the vertex came from so the + // contact direction follows the paper's above/below logic. The undeformed + // pose remains a fallback if the previous triangle is degenerate or nearly + // coplanar with the vertex. + const bool previousValid = triangleNormalAndSignedDistance( + previousPositions, + vertexIndex, + p1Index, + p2Index, + p3Index, + normal, + outSignedDistance + ); + if (!previousValid || std::abs(outSignedDistance) <= 1e-6f) { + if (!triangleNormalAndSignedDistance( + referencePositions, + vertexIndex, + p1Index, + p2Index, + p3Index, + normal, + outSignedDistance + )) { + return false; + } + } + + flipNormal = outSignedDistance < 0.0f; + if (flipNormal) outSignedDistance = -outSignedDistance; + return true; +} + +float averageEdgeLength(const Eigen::VectorXf& restLengths) { + float sum = 0.0f; + unsigned int count = 0u; + for (int i = 0; i < restLengths.size(); ++i) { + const float restLength = restLengths[i]; + if (restLength <= 1e-6f) continue; + sum += restLength; + ++count; + } + if (count == 0u) return 0.0f; + return sum / static_cast(count); +} + +Eigen::Vector3f closestPointOnPlane( + const Eigen::Vector3f& point, + const Eigen::Vector3f& planePoint, + const Eigen::Vector3f& planeNormal, + float* outSignedDistance = nullptr +) { + const float signedDistance = (point - planePoint).dot(planeNormal); + if (outSignedDistance != nullptr) { + *outSignedDistance = signedDistance; + } + return point - signedDistance * planeNormal; +} + +bool segmentPlaneContactPoint( + const Eigen::Vector3f& previousPosition, + const Eigen::Vector3f& predictedPosition, + const Eigen::Vector3f& planePoint, + const Eigen::Vector3f& planeNormal, + float previousSignedDistance, + float predictedSignedDistance, + Eigen::Vector3f& outContactPoint +); + +bool pointInsideBox( + const Eigen::Vector3f& point, + const Eigen::Vector3f& center, + const Eigen::Vector3f& halfExtents +) { + const Eigen::Vector3f local = point - center; + const Eigen::Vector3f absLocal = local.cwiseAbs(); + return absLocal.x() <= halfExtents.x() + && absLocal.y() <= halfExtents.y() + && absLocal.z() <= halfExtents.z(); +} + +struct BoxFaceContact { + Eigen::Vector3f point; + Eigen::Vector3f normal; + float signedDistance; + float penetrationDepth; +}; + +std::vector boxFaceContacts( + const Eigen::Vector3f& point, + const Eigen::Vector3f& center, + const Eigen::Vector3f& halfExtents, + float edgeCornerTolerance +) { + std::vector contacts; + if (!pointInsideBox(point, center, halfExtents)) { + return contacts; + } + + const Eigen::Vector3f normals[6] = { + Eigen::Vector3f(1.0f, 0.0f, 0.0f), + Eigen::Vector3f(-1.0f, 0.0f, 0.0f), + Eigen::Vector3f(0.0f, 1.0f, 0.0f), + Eigen::Vector3f(0.0f, -1.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, 1.0f), + Eigen::Vector3f(0.0f, 0.0f, -1.0f) + }; + const Eigen::Vector3f planePoints[6] = { + center + Eigen::Vector3f(halfExtents.x(), 0.0f, 0.0f), + center - Eigen::Vector3f(halfExtents.x(), 0.0f, 0.0f), + center + Eigen::Vector3f(0.0f, halfExtents.y(), 0.0f), + center - Eigen::Vector3f(0.0f, halfExtents.y(), 0.0f), + center + Eigen::Vector3f(0.0f, 0.0f, halfExtents.z()), + center - Eigen::Vector3f(0.0f, 0.0f, halfExtents.z()) + }; + + std::array candidates; + float minimumPenetration = std::numeric_limits::infinity(); + for (int faceIndex = 0; faceIndex < 6; ++faceIndex) { + const float signedDistance = (point - planePoints[faceIndex]).dot(normals[faceIndex]); + const float penetrationDepth = -signedDistance; + candidates[faceIndex] = BoxFaceContact{ + planePoints[faceIndex], + normals[faceIndex], + signedDistance, + penetrationDepth + }; + minimumPenetration = std::min(minimumPenetration, penetrationDepth); + } + + for (const BoxFaceContact& candidate : candidates) { + if (candidate.penetrationDepth <= minimumPenetration + edgeCornerTolerance) { + contacts.push_back(candidate); + } + } + + return contacts; +} + +Eigen::Vector3f closestPointOnBox( + const Eigen::Vector3f& point, + const Eigen::Vector3f& center, + const Eigen::Vector3f& halfExtents +) { + const Eigen::Vector3f boxMin = center - halfExtents; + const Eigen::Vector3f boxMax = center + halfExtents; + return point.cwiseMax(boxMin).cwiseMin(boxMax); +} + +std::vector closestPointBoxContacts( + const Eigen::Vector3f& point, + const Eigen::Vector3f& center, + const Eigen::Vector3f& halfExtents, + float proximityTolerance, + float collisionEps +) { + std::vector contacts; + if (pointInsideBox(point, center, halfExtents)) { + return contacts; + } + + const Eigen::Vector3f closestPoint = closestPointOnBox(point, center, halfExtents); + Eigen::Vector3f normal = point - closestPoint; + float distance = normal.norm(); + if (distance > proximityTolerance) { + return contacts; + } + + if (distance <= 1e-8f) { + const Eigen::Vector3f local = point - center; + int dominantAxis = 0; + if (std::abs(local.y()) > std::abs(local[dominantAxis])) dominantAxis = 1; + if (std::abs(local.z()) > std::abs(local[dominantAxis])) dominantAxis = 2; + normal = Eigen::Vector3f::Zero(); + normal[dominantAxis] = (local[dominantAxis] >= 0.0f) ? 1.0f : -1.0f; + distance = 0.0f; + } + else { + normal /= distance; + } + + const Eigen::Vector3f contactPlanePoint = closestPoint + proximityTolerance * normal; + const float signedDistance = (point - contactPlanePoint).dot(normal); + if (signedDistance > collisionEps) { + return contacts; + } + + contacts.push_back(BoxFaceContact{ + contactPlanePoint, + normal, + signedDistance, + -std::min(signedDistance, 0.0f) + }); + return contacts; +} + +std::vector sweptBoxFaceContacts( + const Eigen::Vector3f& previousPosition, + const Eigen::Vector3f& predictedPosition, + const Eigen::Vector3f& center, + const Eigen::Vector3f& halfExtents, + float faceBoundsTolerance, + float collisionEps +) { + std::vector contacts; + const Eigen::Vector3f normals[6] = { + Eigen::Vector3f(1.0f, 0.0f, 0.0f), + Eigen::Vector3f(-1.0f, 0.0f, 0.0f), + Eigen::Vector3f(0.0f, 1.0f, 0.0f), + Eigen::Vector3f(0.0f, -1.0f, 0.0f), + Eigen::Vector3f(0.0f, 0.0f, 1.0f), + Eigen::Vector3f(0.0f, 0.0f, -1.0f) + }; + const Eigen::Vector3f planePoints[6] = { + center + Eigen::Vector3f(halfExtents.x(), 0.0f, 0.0f), + center - Eigen::Vector3f(halfExtents.x(), 0.0f, 0.0f), + center + Eigen::Vector3f(0.0f, halfExtents.y(), 0.0f), + center - Eigen::Vector3f(0.0f, halfExtents.y(), 0.0f), + center + Eigen::Vector3f(0.0f, 0.0f, halfExtents.z()), + center - Eigen::Vector3f(0.0f, 0.0f, halfExtents.z()) + }; + + for (int faceIndex = 0; faceIndex < 6; ++faceIndex) { + const Eigen::Vector3f& planePoint = planePoints[faceIndex]; + const Eigen::Vector3f& planeNormal = normals[faceIndex]; + const float previousSignedDistance = (previousPosition - planePoint).dot(planeNormal); + const float predictedSignedDistance = (predictedPosition - planePoint).dot(planeNormal); + const bool crossedFace = previousSignedDistance > collisionEps && predictedSignedDistance < 0.0f; + if (!crossedFace) continue; + + Eigen::Vector3f contactPoint = Eigen::Vector3f::Zero(); + if (!segmentPlaneContactPoint( + previousPosition, + predictedPosition, + planePoint, + planeNormal, + previousSignedDistance, + predictedSignedDistance, + contactPoint + )) { + continue; + } + + const Eigen::Vector3f local = contactPoint - center; + bool insideFaceBounds = true; + for (int axis = 0; axis < 3; ++axis) { + if (std::abs(planeNormal[axis]) > 0.5f) continue; + if (std::abs(local[axis]) > halfExtents[axis] + faceBoundsTolerance) { + insideFaceBounds = false; + break; + } + } + if (!insideFaceBounds) continue; + + contacts.push_back(BoxFaceContact{ + contactPoint, + planeNormal, + predictedSignedDistance, + -std::min(predictedSignedDistance, 0.0f) + }); + } + + return contacts; +} + +bool segmentPlaneContactPoint( + const Eigen::Vector3f& previousPosition, + const Eigen::Vector3f& predictedPosition, + const Eigen::Vector3f& planePoint, + const Eigen::Vector3f& planeNormal, + float previousSignedDistance, + float predictedSignedDistance, + Eigen::Vector3f& outContactPoint +) { + const float denominator = previousSignedDistance - predictedSignedDistance; + if (std::abs(denominator) <= 1e-8f) { + return false; + } + + const float alpha = std::max(0.0f, std::min(1.0f, previousSignedDistance / denominator)); + outContactPoint = previousPosition + alpha * (predictedPosition - previousPosition); + float residual = 0.0f; + outContactPoint = closestPointOnPlane(outContactPoint, planePoint, planeNormal, &residual); + return true; +} + +void applySphereContactDamping( + std::vector& positions, + const std::vector& previousPositions, + const std::vector& invMass, + const std::vector& sphereColliders, + float collisionEps, + float collisionThickness +) { + if (sphereColliders.empty()) return; + + const float friction = std::max(0.0f, std::min(1.0f, PBDDefaultParam::contactFriction)); + if (friction <= 0.0f) return; + + const float contactBand = std::max(collisionThickness, collisionEps) + collisionEps; + for (const SphereCollider& collider : sphereColliders) { + for (unsigned int i = 0; i < positions.size() && i < previousPositions.size() && i < invMass.size(); ++i) { + if (invMass[i] == 0.0f) continue; + + Eigen::Vector3f radial = positions[i] - collider.center; + const float radialLength = radial.norm(); + if (radialLength > collider.radius + contactBand) continue; + + if (radialLength <= 1e-8f) { + radial = Eigen::Vector3f(0.0f, 0.0f, 1.0f); + } + else { + radial /= radialLength; + } + + const Eigen::Vector3f displacement = positions[i] - previousPositions[i]; + const Eigen::Vector3f normalDisplacement = displacement.dot(radial) * radial; + const Eigen::Vector3f tangentialDisplacement = displacement - normalDisplacement; + positions[i] = previousPositions[i] + normalDisplacement + (1.0f - friction) * tangentialDisplacement; + } + } +} + +float dihedralAngleFromPositions( + /* + C_bend(p1, p2, p3, p4) = + acos( + ( (p2 - p1) x (p3 - p1) ) . ( (p2 - p1) x (p4 - p1) ) + / + ( |(p2 - p1) x (p3 - p1)| * |(p2 - p1) x (p4 - p1)| ) + ) + */ + const std::vector& positions, + unsigned int edge0, + unsigned int edge1, + unsigned int opposite0, + unsigned int opposite1, + bool* valid = nullptr +) { + if (valid != nullptr) *valid = false; + if (edge0 >= positions.size() || edge1 >= positions.size() + || opposite0 >= positions.size() || opposite1 >= positions.size()) { + return 0.0f; + } + + const Eigen::Vector3f& p1 = positions[edge0]; + const Eigen::Vector3f& p2 = positions[edge1]; + const Eigen::Vector3f& p3 = positions[opposite0]; + const Eigen::Vector3f& p4 = positions[opposite1]; + + // Formulate bending on the two triangles adjacent to a shared edge. + // (p1, p2) is the shared edge, and p3/p4 are the opposite vertices. + const Eigen::Vector3f normal1 = (p2 - p1).cross(p3 - p1); + const Eigen::Vector3f normal2 = (p2 - p1).cross(p4 - p1); + const float area1 = normal1.norm(); + const float area2 = normal2.norm(); + + // Degenerate triangles have undefined normals and therefore no meaningful dihedral angle. + // In this case we skip evaluation/projection for robustness. + if (area1 <= 1e-8f || area2 <= 1e-8f) { + return 0.0f; + } + + // The bend constraint uses theta = acos(n1 . n2). Clamp the dot product to + // stay in the valid acos domain under floating-point roundoff. + const float normalizedDot = normal1.dot(normal2) / (area1 * area2); + const float dotValue = std::max(-1.0f, std::min(1.0f, normalizedDot)); + if (valid != nullptr) *valid = true; + return std::acos(dotValue); +} +} + +PBDConstraint::PBDConstraint( + std::vector particleIndices, + float stiffnessValue, + PBDConstraintType constraintType +) + : particleIndices(std::move(particleIndices)), + stiffnessValue(stiffnessValue), + constraintType(constraintType) {} + +std::size_t PBDConstraint::cardinality() const { + return particleIndices.size(); +} + +const std::vector& PBDConstraint::indices() const { + return particleIndices; +} + +float PBDConstraint::stiffness() const { + return stiffnessValue; +} + +void PBDConstraint::setStiffness(float stiffness) { + stiffnessValue = std::max(0.0f, std::min(1.0f, stiffness)); +} + +PBDConstraintType PBDConstraint::type() const { + return constraintType; +} + +bool PBDConstraint::isViolated(const std::vector& positions, float epsilon) const { + const float value = evaluate(positions); + + if (constraintType == PBDConstraintType::Equality) { + return std::abs(value) > epsilon; + } + + // Inequality constraints are only projected when violated. + // For C(p) >= 0, violation means C(p) < 0. + return value < 0.0f; +} + +void PBDConstraint::project( + std::vector& positions, + const std::vector& invMass, + unsigned int solverIterations, + float epsilon +) const { + if (!isViolated(positions, epsilon)) return; + + std::vector gradientValues; + gradients(positions, gradientValues); + if (gradientValues.size() != particleIndices.size()) return; + + /* + Generic PBD projection from the paper: + q_j = grad_{p_j} C(p) + s = C(positions) / ( sum_j ( w_j * | grad_{p_j} C(positions) |^2 ) ) + Delta p_i = -k' * s * w_i * grad_{p_i} C(positions) + Every concrete constraint in the solver only needs to provide C(.) and its gradients q_j; + the correction rule itself stays identical for: + distance, collision, fixed-point, and dihedral bending constraints. + */ + + float denominator = 0.0f; + for (std::size_t j = 0; j < particleIndices.size(); ++j) { + const unsigned int particle = particleIndices[j]; + if (particle >= positions.size() || particle >= invMass.size()) return; + denominator += invMass[particle] * gradientValues[j].squaredNorm(); + } + + if (denominator <= epsilon) return; + + const float constraintValue = evaluate(positions); + const float scale = -correctedStiffness(stiffnessValue, solverIterations) + * (constraintValue / denominator); + + for (std::size_t j = 0; j < particleIndices.size(); ++j) { + const unsigned int particle = particleIndices[j]; + positions[particle] += scale * invMass[particle] * gradientValues[j]; + } +} + +DistanceConstraint::DistanceConstraint(unsigned int i, unsigned int j, float restLength, float stiffness) + : PBDConstraint(std::vector{ i, j }, stiffness, PBDConstraintType::Equality), + restLength(restLength) {} + +float DistanceConstraint::evaluate(const std::vector& positions) const { + // Distance constraint: + // C(p1, p2) = |p1 - p2| - d + const Vector3f delta = positions[particleIndices[0]] - positions[particleIndices[1]]; + return delta.norm() - restLength; +} + +void DistanceConstraint::gradients( + const std::vector& positions, + std::vector& outGradients +) const { + outGradients.assign(2, Vector3f::Zero()); + + const Vector3f delta = positions[particleIndices[0]] - positions[particleIndices[1]]; + const float length = delta.norm(); + if (length <= 1e-8f) return; + + // grad_{p1} C(positions) = (p1 - p2) / |p1 - p2| + // grad_{p2} C(positions) = -grad_{p1} C(positions) + const Vector3f direction = delta / length; + outGradients[0] = direction; + outGradients[1] = -direction; +} + +DihedralBendConstraint::DihedralBendConstraint( + unsigned int edge0, + unsigned int edge1, + unsigned int opposite0, + unsigned int opposite1, + float restAngle, + float stiffness +) + : PBDConstraint( + std::vector{ edge0, edge1, opposite0, opposite1 }, + stiffness, + PBDConstraintType::Equality + ), + restAngle(restAngle) {} + +float DihedralBendConstraint::angle(const std::vector& positions, bool* valid) const { + return dihedralAngleFromPositions( + positions, + particleIndices[0], + particleIndices[1], + particleIndices[2], + particleIndices[3], + valid + ); +} + +float DihedralBendConstraint::evaluate(const std::vector& positions) const { + bool valid = false; + const float currentAngle = angle(positions, &valid); + if (!valid) return 0.0f; + + // C_bend(p1, p2, p3, p4) = theta - theta_0, + // where theta is the current dihedral angle between the two incident triangles, + // and theta_0 is the rest angle measured in the reference pose. + // This formulation measure the deviation from the rest angle, + // and therefore projects toward the rest angle instead of toward a flat fold. + return currentAngle - restAngle; +} + +void DihedralBendConstraint::gradients( + const std::vector& positions, + std::vector& outGradients +) const { + outGradients.assign(4, Vector3f::Zero()); + if (particleIndices.size() != 4) return; + + const unsigned int iShared0 = particleIndices[0]; + const unsigned int iShared1 = particleIndices[1]; + const unsigned int iOpp0 = particleIndices[2]; + const unsigned int iOpp1 = particleIndices[3]; + if (iShared0 >= positions.size() || iShared1 >= positions.size() + || iOpp0 >= positions.size() || iOpp1 >= positions.size()) { + return; + } + + // Local particle order in this solver is: + // [0] shared edge vertex p1 + // [1] shared edge vertex p2 + // [2] opposite vertex p3 of the first triangle + // [3] opposite vertex p4 of the second triangle + // + // For the analytic derivative it is convenient to switch to the standard + // dihedral notation used in many PBD derivations: + // p0 = opposite vertex of triangle 1 + // p1 = opposite vertex of triangle 2 + // p2 = shared edge vertex 1 + // p3 = shared edge vertex 2 + const Vector3f& p0 = positions[iOpp0]; + const Vector3f& p1 = positions[iOpp1]; + const Vector3f& p2 = positions[iShared0]; + const Vector3f& p3 = positions[iShared1]; + + // The two normalized face normals encode the fold between adjacent + // triangles. Their dot product is cos(theta), where theta is the dihedral + // angle across the shared edge. + Vector3f n1 = (p2 - p0).cross(p3 - p0); + Vector3f n2 = (p3 - p1).cross(p2 - p1); + const float n1Len = n1.norm(); + const float n2Len = n2.norm(); + + const float eps = 1e-6f; + if (n1Len < eps || n2Len < eps) return; + + n1 /= n1Len; + n2 /= n2Len; + const float area1 = 0.5f * n1Len; + const float area2 = 0.5f * n2Len; + + float d = n1.dot(n2); + d = std::max(-1.0f, std::min(1.0f, d)); + + // C = acos(d) - theta_0, so by the chain rule: + // grad C = -grad d / sqrt(1 - d^2) + // The denominator becomes singular near perfectly flat or fully folded states, + // so clamp it away from zero for robustness (avoid division by zero when d = +1 or d = -1). + const float sinTheta = std::sqrt(std::max(1.0f - d * d, eps)); + + // These q-vectors are the analytic gradients of d = cos(theta) with respect + // to the four dihedral vertices. They replace the previous finite-difference + // approximation while leaving the generic PBD projection rule unchanged. + // The Muller-style formulas are expressed in terms of triangle areas. Since + // |(a x b)| is twice the triangle area, use 0.5 * |cross| here rather than + // the raw cross-product norm to avoid an overall factor-of-two error. + const Vector3f q0 = + ((p2 - p1).cross(n1) + n2.cross(p2 - p1) * d) / area2; + const Vector3f q1 = + -((p3 - p1).cross(n1) + n2.cross(p3 - p1) * d) / area2; + const Vector3f q2 = + -((p3 - p0).cross(n2) + n1.cross(p3 - p0) * d) / area1; + const Vector3f q3 = + ((p2 - p0).cross(n2) + n1.cross(p2 - p0) * d) / area1; + + const Vector3f g0 = -q0 / sinTheta; + const Vector3f g1 = -q1 / sinTheta; + const Vector3f g2 = -q2 / sinTheta; + const Vector3f g3 = -q3 / sinTheta; + + // Map the paper-style order back to this solver's local particle order. + outGradients[0] = g2; + outGradients[1] = g3; + outGradients[2] = g0; + outGradients[3] = g1; +} + +FixedPointConstraint::FixedPointConstraint(unsigned int i, const Eigen::Vector3f& fixedPosition, float stiffness) + : PBDConstraint(std::vector{ i }, stiffness, PBDConstraintType::Equality), + fixedPosition(fixedPosition) {} + +void FixedPointConstraint::setFixedPosition(const Eigen::Vector3f& position) { + fixedPosition = position; +} + +const Eigen::Vector3f& FixedPointConstraint::position() const { + return fixedPosition; +} + +float FixedPointConstraint::evaluate(const std::vector& positions) const { + // For logging/debugging, measure the distance to the target anchor. + return (positions[particleIndices[0]] - fixedPosition).norm(); +} + +void FixedPointConstraint::gradients( + const std::vector& positions, + std::vector& outGradients +) const { + outGradients.assign(1, Vector3f::Zero()); + const Vector3f delta = positions[particleIndices[0]] - fixedPosition; + const float length = delta.norm(); + if (length <= 1e-8f) return; + outGradients[0] = delta / length; +} + +void FixedPointConstraint::project( + std::vector& positions, + const std::vector& invMass, + unsigned int solverIterations, + float epsilon +) const { + const unsigned int particle = particleIndices[0]; + if (particle >= positions.size() || particle >= invMass.size()) return; + + // A fixed point is treated as a moving positional anchor. + // Using k' here keeps the formulation consistent with the paper: + // p_i <- p_i + k' * (p_fixed - p_i) + const float kPrime = correctedStiffness(stiffnessValue, solverIterations); + if ((positions[particle] - fixedPosition).norm() <= epsilon) return; + positions[particle] += kPrime * (fixedPosition - positions[particle]); +} + +CollisionConstraint::CollisionConstraint( + unsigned int i, + const Eigen::Vector3f& center, + float radius, + float stiffness +) + : PBDConstraint(std::vector{ i }, stiffness, PBDConstraintType::Inequality), + collisionKind(CollisionKind::Sphere), + center(center), + radius(radius), + offset(0.0f), + planeNormal(Eigen::Vector3f::Zero()), + selfCollisionNormal(Eigen::Vector3f::Zero()), + selfCollisionBarycentric(Eigen::Vector3f::Zero()) {} + +CollisionConstraint::CollisionConstraint( + unsigned int i, + const Eigen::Vector3f& planePoint, + const Eigen::Vector3f& planeNormal, + float stiffness +) + : PBDConstraint(std::vector{ i }, stiffness, PBDConstraintType::Inequality), + collisionKind(CollisionKind::Plane), + center(planePoint), + radius(0.0f), + offset(0.0f), + planeNormal(planeNormal.normalized()), + edgeSampleWeights(Eigen::Vector2f::Zero()), + selfCollisionNormal(Eigen::Vector3f::Zero()), + selfCollisionBarycentric(Eigen::Vector3f::Zero()) {} + +CollisionConstraint::CollisionConstraint( + unsigned int edge0, + unsigned int edge1, + const Eigen::Vector2f& weights, + const Eigen::Vector3f& planePoint, + const Eigen::Vector3f& planeNormal, + float stiffness +) + : PBDConstraint(std::vector{ edge0, edge1 }, stiffness, PBDConstraintType::Inequality), + collisionKind(CollisionKind::EdgePlane), + center(planePoint), + radius(0.0f), + offset(0.0f), + planeNormal(planeNormal.normalized()), + edgeSampleWeights(weights), + selfCollisionNormal(Eigen::Vector3f::Zero()), + selfCollisionBarycentric(Eigen::Vector3f::Zero()) {} + +CollisionConstraint::CollisionConstraint( + unsigned int vertex, + unsigned int p1, + unsigned int p2, + unsigned int p3, + float thickness, + const Eigen::Vector3f& normal, + const Eigen::Vector3f& barycentric, + float stiffness +) + : PBDConstraint( + std::vector{ vertex, p1, p2, p3 }, + stiffness, + PBDConstraintType::Inequality + ), + collisionKind(CollisionKind::SelfVertexTriangle), + center(Eigen::Vector3f::Zero()), + radius(0.0f), + offset(thickness), + planeNormal(Eigen::Vector3f::Zero()), + selfCollisionNormal(normal), + selfCollisionBarycentric(barycentric) {} + +float CollisionConstraint::evaluate(const std::vector& positions) const { + if (collisionKind == CollisionKind::Sphere) { + // Collision inequality: + // C(p_i) = |p_i - c| - r >= 0 + return (positions[particleIndices[0]] - center).norm() - radius; + } + + if (collisionKind == CollisionKind::Plane) { + return (positions[particleIndices[0]] - center).dot(planeNormal); + } + + if (collisionKind == CollisionKind::EdgePlane) { + if (particleIndices.size() != 2) return 0.0f; + const Vector3f samplePoint = edgeSampleWeights[0] * positions[particleIndices[0]] + + edgeSampleWeights[1] * positions[particleIndices[1]]; + return (samplePoint - center).dot(planeNormal); + } + + if (particleIndices.size() != 4) return 0.0f; + + const Vector3f& q = positions[particleIndices[0]]; + const Vector3f& p1 = positions[particleIndices[1]]; + + // Self-collision inequality from Muller et al. 2007, Sec. 4.3: + // C(q, p1, p2, p3) = (q - p1) . n - h >= 0 + return (q - p1).dot(selfCollisionNormal) - offset; +} + +void CollisionConstraint::gradients( + const std::vector& positions, + std::vector& outGradients +) const { + if (collisionKind == CollisionKind::Sphere) { + outGradients.assign(1, Vector3f::Zero()); + Vector3f delta = positions[particleIndices[0]] - center; + const float length = delta.norm(); + + if (length <= 1e-8f) { + outGradients[0] = Vector3f(0.0f, 0.0f, 1.0f); + return; + } + + // grad_{p_i} C = (p_i - c) / |p_i - c| + outGradients[0] = delta / length; + return; + } + + if (collisionKind == CollisionKind::Plane) { + outGradients.assign(1, Vector3f::Zero()); + outGradients[0] = planeNormal; + return; + } + + if (collisionKind == CollisionKind::EdgePlane) { + outGradients.assign(2, Vector3f::Zero()); + outGradients[0] = edgeSampleWeights[0] * planeNormal; + outGradients[1] = edgeSampleWeights[1] * planeNormal; + return; + } + + outGradients.assign(4, Vector3f::Zero()); + if (particleIndices.size() != 4) return; + + // Using barycentric weights for the triangle vertices reproduces the paper's + // point-triangle correction split inside the generic PBD projection rule. + outGradients[0] = selfCollisionNormal; + outGradients[1] = -selfCollisionBarycentric[0] * selfCollisionNormal; + outGradients[2] = -selfCollisionBarycentric[1] * selfCollisionNormal; + outGradients[3] = -selfCollisionBarycentric[2] * selfCollisionNormal; +} + +// ----------------------------- +// 1. Constructor / Initialization +// ----------------------------- + +PBDSolver::PBDSolver(pbd_system* system, float* vbuff) + : system(system), vbuff(vbuff), + solverIterations(PBDDefaultParam::solverIterations), + dampingFactor(PBDDefaultParam::dampingFactor), + collisionEps(PBDDefaultParam::collisionEps), + structuralStiffness(1.0f), + shearStiffness(0.7f), + bendStiffness(0.03f), + planeFriction(PBDDefaultParam::contactFriction), + selfCollisionThickness(PBDDefaultParam::collisionEps), + selfCollisionStiffness(PBDDefaultParam::selfCollisionStiffness), + selfCollisionCellSize(PBDDefaultParam::collisionEps), + maxSelfCollisionContactsPerVertex(PBDDefaultParam::maxSelfCollisionContactsPerVertex), + velocitySleepThreshold(PBDDefaultParam::velocitySleepThreshold), + gravity(PBDDefaultParam::gravity), + windConfig(), + simulationTime(0.0f), + windSpeedState(0.0f) { + assert(system != nullptr); + assert(vbuff != nullptr); + + initializeState(); + meshAdjacency.resize(system->n_points); + for (const Edge& edge : system->spring_list) { + if (edge.first >= meshAdjacency.size() || edge.second >= meshAdjacency.size()) continue; + meshAdjacency[edge.first].insert(edge.second); + meshAdjacency[edge.second].insert(edge.first); + } + + float minRestLength = std::numeric_limits::max(); + for (int i = 0; i < system->rest_lengths.size(); ++i) { + const float restLength = system->rest_lengths[i]; + if (restLength > 1e-6f) { + minRestLength = std::min(minRestLength, restLength); + } + } + const float meanRestLength = averageEdgeLength(system->rest_lengths); + + if (minRestLength < std::numeric_limits::max()) { + // Self-collision should model a thin cloth thickness, not half an edge + // length. Large thickness inflates folded cloth and causes false + // repulsion between nearby layers. + selfCollisionThickness = std::max(0.05f * minRestLength, collisionEps); + } + if (meanRestLength > 0.0f) { + // Broad-phase hashing is more reliable when cells match the cloth's edge + // scale instead of the much smaller thickness band. + selfCollisionCellSize = std::max(meanRestLength, collisionEps); + } +} + +void PBDSolver::initializeState() { + const unsigned int n = system->n_points; + // Pseudo-code (1)-(3): initialize x, p, v, w. + restX.resize(n); + x.resize(n); + p.resize(n); + v.resize(n); + invMass.resize(n); + planeContactPoints.assign(n, Vector3f::Zero()); + planeContactNormals.assign(n, Vector3f(0.0f, 0.0f, 1.0f)); + planeContactSignedDistances.assign(n, std::numeric_limits::infinity()); + + for (unsigned int i = 0; i < n; ++i) { + // Read initial position from render vertex buffer. + Vector3f pos( + vbuff[3 * i + 0], + vbuff[3 * i + 1], + vbuff[3 * i + 2] + ); + + // Preserve the initial pose as the reference side for future + // vertex-triangle self-collision tests. + restX[i] = pos; + x[i] = pos; + p[i] = pos; + + // Initial velocity = 0. + v[i] = Vector3f(0.0f, 0.0f, 0.0f); + + // Inverse mass w_i = 1 / m_i. + const float mass = system->masses[i]; + invMass[i] = (mass > 0.0f) ? (1.0f / mass) : 0.0f; + } +} + +// ----------------------------- +// 2. Main Simulation Entry / Solver Loop +// ----------------------------- + +void PBDSolver::step(float dt) { + // Paper-aligned solver structure: + // 1. apply external forces + // 2. damp velocities + // 3. predict positions + // 4. generate collision constraints + // 5. iterate projections + // 6. update velocities from projected positions + // 7. apply post-collision velocity manipulation (friction) + // 8. commit positions + applyExternalForces(dt); + dampVelocities(); + predictPositions(dt); + generateCollisionConstraints(); + + for (unsigned int iteration = 0; iteration < solverIterations; ++iteration) { + projectConstraints(persistentConstraints); + projectConstraints(generatedCollisionConstraints); + } + applySphereContactDamping( + p, + x, + invMass, + sphereColliders, + collisionEps, + selfCollisionThickness + ); + updateSelfCollisionDebugStats(p, true); + + updateVelocities(dt); + applyPlaneContactVelocityDamping(); + commitPositions(); +} + +void PBDSolver::solve(unsigned int n) { + // Wrapper used by the app to override the number of constraint-projection + // iterations for a single simulation time step. + // + // Important: n is NOT the number of outer time steps. This function still + // advances the simulation by exactly one dt via step(system->time_step). + // Instead, n temporarily replaces solverIterations, which controls the + // inner PBD loop inside step(): + // for each solver iteration: + // project persistent constraints + // project generated collision constraints + // + // After that one step is finished, the previous default iteration count is + // restored. This lets the caller choose the quality/cost of one frame update + // without permanently changing solver configuration. + const unsigned int previousIterations = solverIterations; + solverIterations = n; + step(system->time_step); + solverIterations = previousIterations; +} + +// ----------------------------- +// 3. Core PBD Pipeline Stages +// ----------------------------- + +void PBDSolver::applyExternalForces(float dt) { + const unsigned int n = system->n_points; + // Pseudo-code (5): + // v_i <- v_i + dt * f_ext_i / m_i + // Gravity is stored directly as an acceleration vector. + for (unsigned int i = 0; i < n; ++i) { + if (invMass[i] == 0.0f) continue; // fixed particles do not accelerate + v[i] += dt * gravity; + } + + // Wind is treated as an external aerodynamic force instead of a positional + // constraint so the solver remains lightweight and localized to one force + // accumulation stage. + applyAerodynamicForces(dt); + simulationTime += std::max(dt, 0.0f); +} + +bool PBDSolver::isWindEnabled() const { + return windConfig.inputMode != PBDWindInputMode::Disabled; +} + +float PBDSolver::triangleFaceArea(unsigned int i0, unsigned int i1, unsigned int i2) const { + if (i0 >= x.size() || i1 >= x.size() || i2 >= x.size()) return 0.0f; + + // Triangle area used by the aerodynamic model: + // A = 0.5 * || (x_1 - x_0) x (x_2 - x_0) || + // Larger faces catch more air and therefore produce larger drag forces. + Vector3f faceNormal = Vector3f::Zero(); + float doubleArea = 0.0f; + if (!triangleUnitNormal(x, i0, i1, i2, faceNormal, &doubleArea)) return 0.0f; + return 0.5f * doubleArea; +} + +PBDSolver::Vector3f PBDSolver::triangleFaceNormal(unsigned int i0, unsigned int i1, unsigned int i2) const { + // Reuse the shared triangle-normal helper already used by collision code: + // n = ((x_1 - x_0) x (x_2 - x_0)) / || (x_1 - x_0) x (x_2 - x_0) || + Vector3f normal = Vector3f::Zero(); + if (!triangleUnitNormal(x, i0, i1, i2, normal)) return Vector3f::Zero(); + return normal; +} + +PBDSolver::Vector3f PBDSolver::averageTriangleVelocity(unsigned int i0, unsigned int i1, unsigned int i2) const { + if (i0 >= v.size() || i1 >= v.size() || i2 >= v.size()) return Vector3f::Zero(); + + // Paper-style face velocity approximation: + // v_object = (v_0 + v_1 + v_2) / 3 + // The drag model uses the triangle-average velocity instead of a per-vertex + // air speed so each face receives one aerodynamic force sample. + return (v[i0] + v[i1] + v[i2]) / 3.0f; +} + +float PBDSolver::proceduralWindNoise(const Vector3f& triangleCenter) const { + // Lightweight flutter term layered on top of the deterministic gust so the + // cloth does not move as if every face saw the exact same wind history. + const float phase0 = 1.73f * simulationTime + 0.91f * triangleCenter.x() + 0.37f * triangleCenter.z(); + const float phase1 = 2.41f * simulationTime + 0.63f * triangleCenter.y() - 0.52f * triangleCenter.x(); + return 0.5f * std::sin(phase0) + 0.5f * std::cos(phase1); +} + +float PBDSolver::currentWindSpeed(const Vector3f& triangleCenter, float dt) { + if (!isWindEnabled()) return 0.0f; + + // User input selects either a direct speed or an acceleration integrated in + // time. The gust model follows the requested paper-style modulation: + // u(t) = u_base + A_gust * sin(2 * pi * f_gust * t) + noise(t, x) + // The final speed is clamped so gust/noise cannot create extreme forces. + const float maxWindSpeed = std::max(0.0f, windConfig.maxWindSpeed); + if (windConfig.inputMode == PBDWindInputMode::Speed) { + windSpeedState = std::max(0.0f, windConfig.baseSpeed); + } + else if (windConfig.inputMode == PBDWindInputMode::Acceleration) { + windSpeedState = std::max(0.0f, std::min(maxWindSpeed, windSpeedState + dt * windConfig.baseAcceleration)); + } + + const float gust = windConfig.gustAmplitude + * std::sin(PBDDefaultParam::twoPi * windConfig.gustFrequency * simulationTime); + const float noise = windConfig.noiseStrength * proceduralWindNoise(triangleCenter); + return std::max(0.0f, std::min(maxWindSpeed, windSpeedState + gust + noise)); +} + +void PBDSolver::applyAerodynamicForces(float dt) { + if (!isWindEnabled() || dt <= 0.0f) return; + if (system->triangle_indices.size() < 3 || system->triangle_indices.size() % 3 != 0) return; + + for (std::size_t triangle = 0; triangle + 2 < system->triangle_indices.size(); triangle += 3) { + const unsigned int i0 = system->triangle_indices[triangle + 0]; + const unsigned int i1 = system->triangle_indices[triangle + 1]; + const unsigned int i2 = system->triangle_indices[triangle + 2]; + if (i0 >= x.size() || i1 >= x.size() || i2 >= x.size()) continue; + + const float faceArea = triangleFaceArea(i0, i1, i2); + if (faceArea <= 1e-8f) continue; + + // Face area controls how much air the triangle catches, while the face + // normal controls whether wind strikes the broad face or mostly slides + // along the edge. + const Vector3f faceNormal = triangleFaceNormal(i0, i1, i2); + if (faceNormal.squaredNorm() <= 1e-12f) continue; + + const Vector3f triangleCenter = (x[i0] + x[i1] + x[i2]) / 3.0f; + const float modulatedWindSpeed = currentWindSpeed(triangleCenter, dt); + if (modulatedWindSpeed <= 1e-6f) continue; + + const Vector3f objectVelocity = averageTriangleVelocity(i0, i1, i2); + // Wind velocity for the current face sample: + // u_wind = d_wind * currentWindSpeed + const Vector3f windVelocity = windConfig.windDirection * modulatedWindSpeed; + // Relative air speed from the paper-style drag model: + // v_rel = v_object - u_wind + // Drag depends on the relative velocity between the moving cloth face and + // the surrounding air, using the face-averaged triangle velocity. + const Vector3f relativeVelocity = objectVelocity - windVelocity; + const float relativeSpeed = relativeVelocity.norm(); + if (relativeSpeed <= 1e-6f) continue; + const Vector3f vRelHat = relativeVelocity / relativeSpeed; + + const Vector3f incomingFlowDirection = -vRelHat; + const float orientationFactor = std::abs(faceNormal.dot(incomingFlowDirection)); + if (orientationFactor <= 1e-4f) continue; + + // Drag-only aerodynamic force: + // F_drag = -0.5 * rho * C_D * A * |v_rel|^2 * exposure * v_rel_hat + // where exposure is approximated with |n . (-v_rel_hat)| so broad faces + // catch more air than edge-on faces. + const float dragMagnitude = 0.5f + * windConfig.airDensity + * windConfig.dragCoefficient + * faceArea + * relativeSpeed + * relativeSpeed + * orientationFactor; + // Drag acts opposite the relative air velocity. + const Vector3f dragForce = -dragMagnitude * vRelHat; + + Vector3f liftForce = Vector3f::Zero(); + if (windConfig.liftCoefficient > 1e-6f) { + // Use a two-sided cloth normal so lift reacts consistently whether the + // wind hits the front or back side of the triangle. + const Vector3f adjustedNormal = (faceNormal.dot(relativeVelocity) > 0.0f) + ? faceNormal + : -faceNormal; + + // Paper-style lift direction: + // liftDir = (n_adj x v_rel_hat) x v_rel_hat + // Lift acts perpendicular to relative velocity and lies in the plane + // spanned by the relative velocity and the cloth normal. + Vector3f liftDirection = adjustedNormal.cross(vRelHat).cross(vRelHat); + const float liftDirectionLength = liftDirection.norm(); + if (liftDirectionLength > 1e-6f && std::isfinite(liftDirectionLength)) { + liftDirection /= liftDirectionLength; + + // Stable lift exposure approximation based on the angle between the + // adjusted two-sided normal and the relative airflow direction. + const float liftOrientationTerm = std::abs(adjustedNormal.dot(vRelHat)); + if (liftOrientationTerm > 1e-4f) { + // Optional lift force: + // F_lift = 0.5 * rho * C_L * A * |v_rel|^2 * orientationTerm * liftDir + // Lift is optional and should be tuned carefully because large C_L + // values can destabilize thin cloth quickly. + const float liftMagnitude = 0.5f + * windConfig.airDensity + * windConfig.liftCoefficient + * faceArea + * relativeSpeed + * relativeSpeed + * liftOrientationTerm; + liftForce = liftMagnitude * liftDirection; + } + } + } + + // Gust and noise are small time-varying speed modulations that create a + // lightweight flag-like flutter without introducing a full wind-field + // solver. + const Vector3f totalForce = dragForce + liftForce; + const Vector3f perVertexForce = totalForce / 3.0f; + const unsigned int triangleVertices[3] = { i0, i1, i2 }; + for (unsigned int localVertex = 0; localVertex < 3; ++localVertex) { + const unsigned int vertex = triangleVertices[localVertex]; + if (vertex >= invMass.size() || invMass[vertex] == 0.0f) continue; + // Convert force to acceleration with inverse mass and advance velocity: + // a_i = F_i / m_i = F_i * w_i + // v_i <- v_i + dt * a_i + const Vector3f acceleration = perVertexForce * invMass[vertex]; + v[vertex] += dt * acceleration; + } + } +} + +void PBDSolver::dampVelocities() { + const unsigned int n = system->n_points; + if (n == 0) return; + + // Use strictly dissipative damping so residual cloth motion decays instead + // of preserving rigid-body modes indefinitely. + const float factor = std::max(0.0f, 1.0f - dampingFactor); + for (unsigned int i = 0; i < n; ++i) { + if (invMass[i] == 0.0f) continue; + v[i] *= factor; + } +} + +void PBDSolver::predictPositions(float dt) { + const unsigned int n = system->n_points; + for (unsigned int i = 0; i < n; ++i) { + if (invMass[i] == 0.0f) { + p[i] = x[i]; + continue; + } + + // Pseudo-code (7): predict positions + // p_i = x_i + dt * v_i + p[i] = x[i] + dt * v[i]; + } +} + +void PBDSolver::projectConstraints(const ConstraintList& constraints) { + for (const ConstraintPtr& constraint : constraints) { + constraint->project(p, invMass, solverIterations, collisionEps); + } +} + +void PBDSolver::updateVelocities(float dt) { + // Pseudo-code (16): + // v_i = (p_i - x_i) / dt + if (dt <= 0.0f) return; + + const unsigned int n = system->n_points; + for (unsigned int i = 0; i < n; ++i) { + if (invMass[i] == 0.0f) { + v[i] = Vector3f(0.0f, 0.0f, 0.0f); + continue; + } + + v[i] = (p[i] - x[i]) / dt; + if (v[i].norm() < velocitySleepThreshold) { + v[i] = Vector3f::Zero(); + } + } +} + +// ----------------------------- +// 4. Persistent Constraint Setup / Constraint Construction +// ----------------------------- + +void PBDSolver::pinPoint(unsigned int i) { + fixPoint(i); +} + +void PBDSolver::fixPoint(unsigned int i) { + if (i >= x.size()) return; + + const Vector3f fixedPosition( + vbuff[3 * i + 0], + vbuff[3 * i + 1], + vbuff[3 * i + 2] + ); + + // Mouse dragging acts like a moving fixed point. If the point is already + // fixed, only update its target position. Otherwise add a new persistent + // fixed-point constraint. + auto existing = fixedPointConstraints.find(i); + if (existing != fixedPointConstraints.end()) { + existing->second->setFixedPosition(fixedPosition); + } + else { + auto constraint = std::make_unique(i, fixedPosition, 1.0f); + fixedPointConstraints[i] = constraint.get(); + persistentConstraints.push_back(std::move(constraint)); + } + + // Mark the point as immovable for all other constraints. + x[i] = fixedPosition; + p[i] = fixedPosition; + v[i] = Vector3f(0.0f, 0.0f, 0.0f); + invMass[i] = 0.0f; +} + +void PBDSolver::releasePoint(unsigned int i) { + if (i >= x.size()) return; + + auto existing = fixedPointConstraints.find(i); + if (existing == fixedPointConstraints.end()) return; + + FixedPointConstraint* target = existing->second; + fixedPointConstraints.erase(existing); + + persistentConstraints.erase( + std::remove_if( + persistentConstraints.begin(), + persistentConstraints.end(), + [target](const ConstraintPtr& constraint) { return constraint.get() == target; } + ), + persistentConstraints.end() + ); + + const float mass = system->masses[i]; + invMass[i] = (mass > 0.0f) ? (1.0f / mass) : 0.0f; +} + +void PBDSolver::addDistanceConstraints( + const std::vector& indices, + float stiffness, + std::vector* constraintGroup +) { + for (unsigned int index : indices) { + if (index >= system->spring_list.size()) continue; + + const Edge& edge = system->spring_list[index]; + auto constraint = std::make_unique( + edge.first, + edge.second, + system->rest_lengths[index], + stiffness + ); + if (constraintGroup != nullptr) { + constraintGroup->push_back(constraint.get()); + } + persistentConstraints.push_back(std::move(constraint)); + } +} + +void PBDSolver::addDihedralBendConstraints(float stiffness) { + if (system->triangle_indices.size() < 6 || system->triangle_indices.size() % 3 != 0) return; + + std::unordered_map pendingEdges; + pendingEdges.reserve(system->triangle_indices.size()); + + // Build one bending constraint per interior mesh edge. Each triangle is read + // from the render mesh index buffer, and when the same undirected edge is + // seen a second time we have found the adjacent triangle pair required by the + // paper's 4-particle dihedral constraint. + for (std::size_t triangle = 0; triangle < system->triangle_indices.size(); triangle += 3) { + const std::array vertices = { + system->triangle_indices[triangle + 0], + system->triangle_indices[triangle + 1], + system->triangle_indices[triangle + 2] + }; + + for (int edgeIndex = 0; edgeIndex < 3; ++edgeIndex) { + const unsigned int edge0 = vertices[edgeIndex]; + const unsigned int edge1 = vertices[(edgeIndex + 1) % 3]; + const unsigned int opposite = vertices[(edgeIndex + 2) % 3]; + + SharedEdgeKey key{ std::min(edge0, edge1), std::max(edge0, edge1) }; + auto existing = pendingEdges.find(key); + if (existing == pendingEdges.end()) { + pendingEdges.emplace(key, PendingTriangleEdge{ edge0, edge1, opposite }); + continue; + } + + const PendingTriangleEdge firstTriangle = existing->second; + pendingEdges.erase(existing); + + bool valid = false; + // The rest angle theta_0 is taken from the initial cloth configuration, + // so the solver preserves the reference fold across this shared edge. + const float restAngle = dihedralAngleFromPositions( + x, + firstTriangle.edge0, + firstTriangle.edge1, + firstTriangle.opposite, + opposite, + &valid + ); + if (!valid) continue; + + auto constraint = std::make_unique( + firstTriangle.edge0, + firstTriangle.edge1, + firstTriangle.opposite, + opposite, + restAngle, + stiffness + ); + bendConstraints.push_back(constraint.get()); + persistentConstraints.push_back(std::move(constraint)); + } + } +} + +void PBDSolver::addStructuralConstraints(const std::vector& indices, float stiffness) { + structuralStiffness = std::max(0.0f, std::min(1.0f, stiffness)); + addDistanceConstraints(indices, structuralStiffness, &structuralConstraints); +} + +void PBDSolver::addShearConstraints(const std::vector& indices, float stiffness) { + shearStiffness = std::max(0.0f, std::min(1.0f, stiffness)); + addDistanceConstraints(indices, shearStiffness, &shearConstraints); +} + +void PBDSolver::addBendConstraints(const std::vector&, float stiffness) { + // Bending is generated from adjacent triangle pairs rather than from longer + // spring edges. This measures the cloth's fold angle directly and therefore + // remains meaningful even when in-plane stretching changes edge lengths. + bendStiffness = std::max(0.0f, std::min(1.0f, stiffness)); + addDihedralBendConstraints(bendStiffness); + + + // addDistanceConstraints(indices, stiffness); for distance-based bending --- IGNORE --- +} + +// ----------------------------- +// 5. Collision Handling +// ----------------------------- + +void PBDSolver::generateCollisionConstraints() { + // The paper separates collision detection from constraint projection. + // Persistent colliders stay in sphereColliders, while actual contact + // constraints are generated per step from predicted positions p. + generatedCollisionConstraints.clear(); + generatedSelfCollisionConstraints.clear(); + selfCollisionDebugStats = SelfCollisionDebugStats{}; + planeContactPoints.assign(system->n_points, Vector3f::Zero()); + planeContactNormals.assign(system->n_points, Vector3f(0.0f, 0.0f, 1.0f)); + planeContactSignedDistances.assign(system->n_points, std::numeric_limits::infinity()); + + for (const SphereCollider& collider : sphereColliders) { + for (unsigned int i = 0; i < system->n_points; ++i) { + if (invMass[i] == 0.0f) continue; + + const float contactRadius = collider.radius + std::max(selfCollisionThickness, collisionEps); + const Vector3f delta = p[i] - collider.center; + if (delta.norm() >= contactRadius) continue; + + generatedCollisionConstraints.push_back( + std::make_unique( + i, + collider.center, + contactRadius, + PBDDefaultParam::collisionStiffness + ) + ); + } + } + + for (const BoxCollider& collider : boxColliders) { + const Vector3f expandedHalfExtents = collider.halfExtents + + Vector3f::Constant(std::max(selfCollisionThickness, collisionEps)); + const float edgeCornerTolerance = std::max(collisionEps, 0.5f * selfCollisionThickness); + const float sweptFaceBoundsTolerance = std::max(collisionEps, selfCollisionThickness); + const float closestPointTolerance = std::max(collisionEps, selfCollisionThickness); + + for (unsigned int i = 0; i < system->n_points; ++i) { + if (invMass[i] == 0.0f) continue; + + std::vector contacts = boxFaceContacts( + p[i], + collider.center, + expandedHalfExtents, + edgeCornerTolerance + ); + if (contacts.empty()) { + contacts = sweptBoxFaceContacts( + x[i], + p[i], + collider.center, + expandedHalfExtents, + sweptFaceBoundsTolerance, + collisionEps + ); + } + if (contacts.empty()) { + contacts = closestPointBoxContacts( + p[i], + collider.center, + expandedHalfExtents, + closestPointTolerance, + collisionEps + ); + } + if (contacts.empty()) { + continue; + } + + for (const BoxFaceContact& contact : contacts) { + if (contact.signedDistance < planeContactSignedDistances[i]) { + planeContactSignedDistances[i] = contact.signedDistance; + planeContactPoints[i] = contact.point; + planeContactNormals[i] = contact.normal; + } + + generatedCollisionConstraints.push_back( + std::make_unique( + i, + contact.point, + contact.normal, + PBDDefaultParam::collisionStiffness + ) + ); + } + } + + const Eigen::Vector2f midpointWeights(0.5f, 0.5f); + for (const Edge& edge : system->spring_list) { + if (edge.first >= p.size() || edge.second >= p.size()) continue; + if (edge.first >= x.size() || edge.second >= x.size()) continue; + const float edgeInvMass = invMass[edge.first] + invMass[edge.second]; + if (edgeInvMass <= 0.0f) continue; + + const Vector3f predictedMidpoint = 0.5f * (p[edge.first] + p[edge.second]); + const Vector3f previousMidpoint = 0.5f * (x[edge.first] + x[edge.second]); + + std::vector contacts = boxFaceContacts( + predictedMidpoint, + collider.center, + expandedHalfExtents, + edgeCornerTolerance + ); + if (contacts.empty()) { + contacts = sweptBoxFaceContacts( + previousMidpoint, + predictedMidpoint, + collider.center, + expandedHalfExtents, + sweptFaceBoundsTolerance, + collisionEps + ); + } + if (contacts.empty()) { + contacts = closestPointBoxContacts( + predictedMidpoint, + collider.center, + expandedHalfExtents, + closestPointTolerance, + collisionEps + ); + } + if (contacts.empty()) continue; + + for (const BoxFaceContact& contact : contacts) { + if (contact.signedDistance < planeContactSignedDistances[edge.first]) { + planeContactSignedDistances[edge.first] = contact.signedDistance; + planeContactPoints[edge.first] = contact.point; + planeContactNormals[edge.first] = contact.normal; + } + if (contact.signedDistance < planeContactSignedDistances[edge.second]) { + planeContactSignedDistances[edge.second] = contact.signedDistance; + planeContactPoints[edge.second] = contact.point; + planeContactNormals[edge.second] = contact.normal; + } + + generatedCollisionConstraints.push_back( + std::make_unique( + edge.first, + edge.second, + midpointWeights, + contact.point, + contact.normal, + PBDDefaultParam::collisionStiffness + ) + ); + } + } + } + + for (const PlaneCollider& collider : planeColliders) { + for (unsigned int i = 0; i < system->n_points; ++i) { + if (invMass[i] == 0.0f) continue; + + const float previousSignedDistance = (x[i] - collider.point).dot(collider.normal); + const float predictedSignedDistance = (p[i] - collider.point).dot(collider.normal); + const bool crossedPlane = previousSignedDistance > collisionEps && predictedSignedDistance < 0.0f; + const bool insideOrNearPlane = predictedSignedDistance < collisionEps; + if (!crossedPlane && !insideOrNearPlane) continue; + + Vector3f contactPoint = Vector3f::Zero(); + if (crossedPlane) { + if (!segmentPlaneContactPoint( + x[i], + p[i], + collider.point, + collider.normal, + previousSignedDistance, + predictedSignedDistance, + contactPoint + )) { + contactPoint = closestPointOnPlane(p[i], collider.point, collider.normal); + } + } + else { + // Static fallback: when the predicted point is already on or below the + // floor, use the plane closest point as the contact anchor. + contactPoint = closestPointOnPlane(p[i], collider.point, collider.normal); + } + + if (predictedSignedDistance < planeContactSignedDistances[i]) { + planeContactSignedDistances[i] = predictedSignedDistance; + planeContactPoints[i] = contactPoint; + planeContactNormals[i] = collider.normal; + } + + generatedCollisionConstraints.push_back( + std::make_unique( + i, + contactPoint, + collider.normal, + PBDDefaultParam::collisionStiffness + ) + ); + } + } + + generateSelfCollisionConstraints(); + updateSelfCollisionDebugStats(p, false); +} + +void PBDSolver::addSphereCollider(const Vector3f& center, float radius) { + // Store persistent collision geometry. Actual contact constraints are + // generated each step from predicted positions. + sphereColliders.push_back(SphereCollider{ center, radius }); +} + +void PBDSolver::addBoxCollider(const Vector3f& center, const Vector3f& halfExtents) { + if (halfExtents.x() <= 0.0f || halfExtents.y() <= 0.0f || halfExtents.z() <= 0.0f) return; + boxColliders.push_back(BoxCollider{ center, halfExtents }); +} + +void PBDSolver::addPlaneCollider(const Vector3f& point, const Vector3f& normal) { + const float normalLength = normal.norm(); + if (normalLength <= 1e-8f) return; + planeColliders.push_back(PlaneCollider{ point, normal / normalLength }); +} + +void PBDSolver::applyPlaneContactVelocityDamping() { + if (planeContactSignedDistances.empty()) return; + + const float friction = std::max(0.0f, std::min(1.0f, planeFriction)); + const float tangentialSleepSpeed = velocitySleepThreshold; + for (unsigned int i = 0; i < v.size() && i < invMass.size() + && i < planeContactSignedDistances.size() + && i < planeContactNormals.size(); ++i) { + if (invMass[i] == 0.0f) continue; + if (!std::isfinite(planeContactSignedDistances[i])) continue; + + const Vector3f& contactNormal = planeContactNormals[i]; + const float normalSpeed = v[i].dot(contactNormal); + const float separatingSpeed = std::max(0.0f, normalSpeed); + const Vector3f normalVelocity = separatingSpeed * contactNormal; + Vector3f tangentialVelocity = v[i] - normalSpeed * contactNormal; + if (friction > 0.0f) { + tangentialVelocity *= (1.0f - friction); + } + if (tangentialVelocity.norm() <= tangentialSleepSpeed) { + tangentialVelocity = Vector3f::Zero(); + } + + v[i] = normalVelocity + tangentialVelocity; + } +} + +// ----------------------------- +// 6. Self-Collision Utilities / Helpers +// ----------------------------- + +void PBDSolver::generateSelfCollisionConstraints() { + if (system->triangle_indices.size() < 3 || selfCollisionCellSize <= 0.0f) return; + + std::unordered_map, SpatialHashKeyHash> verticesByCell; + verticesByCell.reserve(p.size()); + std::vector contactsPerVertex(system->n_points, 0u); + + for (unsigned int vertex = 0; vertex < p.size(); ++vertex) { + insertSweptVertexIntoHash( + x[vertex], + p[vertex], + selfCollisionThickness, + selfCollisionCellSize, + vertex, + verticesByCell + ); + } + + for (std::size_t triangle = 0; triangle + 2 < system->triangle_indices.size(); triangle += 3) { + const unsigned int p1Index = system->triangle_indices[triangle + 0]; + const unsigned int p2Index = system->triangle_indices[triangle + 1]; + const unsigned int p3Index = system->triangle_indices[triangle + 2]; + + const Vector3f& p1 = p[p1Index]; + const Vector3f& p2 = p[p2Index]; + const Vector3f& p3 = p[p3Index]; + + Vector3f currentNormal = (p2 - p1).cross(p3 - p1); + const float currentNormalLength = currentNormal.norm(); + if (currentNormalLength <= 1e-8f) continue; + currentNormal /= currentNormalLength; + + const Vector3f& x1 = x[p1Index]; + const Vector3f& x2 = x[p2Index]; + const Vector3f& x3 = x[p3Index]; + const Vector3f minCorner = x1.cwiseMin(x2).cwiseMin(x3) + .cwiseMin(p1).cwiseMin(p2).cwiseMin(p3) + - Vector3f::Constant(selfCollisionThickness); + const Vector3f maxCorner = x1.cwiseMax(x2).cwiseMax(x3) + .cwiseMax(p1).cwiseMax(p2).cwiseMax(p3) + + Vector3f::Constant(selfCollisionThickness); + + const SpatialHashKey minCell = hashPosition(minCorner, selfCollisionCellSize); + const SpatialHashKey maxCell = hashPosition(maxCorner, selfCollisionCellSize); + std::unordered_set processedVertices; + + for (int cellX = minCell.x; cellX <= maxCell.x; ++cellX) { + for (int cellY = minCell.y; cellY <= maxCell.y; ++cellY) { + for (int cellZ = minCell.z; cellZ <= maxCell.z; ++cellZ) { + const SpatialHashKey cellKey{ cellX, cellY, cellZ }; + auto cellVertices = verticesByCell.find(cellKey); + if (cellVertices == verticesByCell.end()) continue; + + for (unsigned int vertexIndex : cellVertices->second) { + if (!processedVertices.insert(vertexIndex).second) continue; + if (contactsPerVertex[vertexIndex] >= maxSelfCollisionContactsPerVertex) continue; + if (vertexIndex == p1Index || vertexIndex == p2Index || vertexIndex == p3Index) continue; + if (meshAdjacency[vertexIndex].count(p1Index) != 0 + || meshAdjacency[vertexIndex].count(p2Index) != 0 + || meshAdjacency[vertexIndex].count(p3Index) != 0) { + continue; + } + + const Vector3f& q = p[vertexIndex]; + const float unsignedDistance = std::abs((q - p1).dot(currentNormal)); + if (unsignedDistance >= selfCollisionThickness) continue; + + bool flipNormal = false; + float previousSignedDistance = 0.0f; + if (!previousFrameSideAndNormal( + x, + restX, + vertexIndex, + p1Index, + p2Index, + p3Index, + flipNormal, + previousSignedDistance + )) { + continue; + } + + // Reorient the contact to the side the vertex occupied in the + // previous frame before building the current-step inequality. + const unsigned int orientedP2 = flipNormal ? p3Index : p2Index; + const unsigned int orientedP3 = flipNormal ? p2Index : p3Index; + const Vector3f& orientedP1 = p[p1Index]; + const Vector3f& orientedP2Pos = p[orientedP2]; + const Vector3f& orientedP3Pos = p[orientedP3]; + + Vector3f orientedNormal = (orientedP2Pos - orientedP1).cross(orientedP3Pos - orientedP1); + const float orientedNormalLength = orientedNormal.norm(); + if (orientedNormalLength <= 1e-8f) continue; + orientedNormal /= orientedNormalLength; + + const float signedDistance = (q - orientedP1).dot(orientedNormal); + if (signedDistance >= selfCollisionThickness) continue; + + // Generate a contact when the vertex crosses the plane over the + // step, enters the thickness band, or is still overlapping from a + // previously missed/self-persistent contact. + const bool crossedPlane = previousSignedDistance > collisionEps && signedDistance < 0.0f; + const bool enteredThicknessBand = previousSignedDistance >= selfCollisionThickness + && signedDistance < selfCollisionThickness; + const bool persistentOverlap = previousSignedDistance < selfCollisionThickness; + if (!crossedPlane && !enteredThicknessBand && !persistentOverlap) continue; + + const Vector3f projectedPoint = q - signedDistance * orientedNormal; + Eigen::Vector3f barycentric; + if (!pointProjectsInsideTriangle( + projectedPoint, + orientedP1, + orientedP2Pos, + orientedP3Pos, + barycentric + )) { + continue; + } + + auto constraint = std::make_unique( + vertexIndex, + p1Index, + orientedP2, + orientedP3, + selfCollisionThickness, + orientedNormal, + barycentric, + selfCollisionStiffness + ); + generatedSelfCollisionConstraints.push_back(constraint.get()); + generatedCollisionConstraints.push_back(std::move(constraint)); + // Limiting contacts per vertex reduces conflicting constraints in + // dense folds, which helps suppress pinching and spike artifacts. + ++contactsPerVertex[vertexIndex]; + } + } + } + } + } +} + +// ----------------------------- +// 7. Parameter / Tuning Interface +// ----------------------------- + +void PBDSolver::setConstraintGroupStiffness(const std::vector& constraints, float stiffness) { + const float clamped = std::max(0.0f, std::min(1.0f, stiffness)); + for (PBDConstraint* constraint : constraints) { + if (constraint != nullptr) { + constraint->setStiffness(clamped); + } + } +} + +void PBDSolver::setGravity(float gravityMagnitude) { + gravity = Vector3f(0.0f, 0.0f, -std::max(0.0f, gravityMagnitude)); +} + +float PBDSolver::getGravity() const { + return -gravity.z(); +} + +void PBDSolver::setDampingFactor(float damping) { + dampingFactor = std::max(0.0f, std::min(1.0f, damping)); +} + +float PBDSolver::getDampingFactor() const { + return dampingFactor; +} + +void PBDSolver::setPlaneFriction(float friction) { + planeFriction = std::max(0.0f, std::min(1.0f, friction)); +} + +float PBDSolver::getPlaneFriction() const { + return planeFriction; +} + +void PBDSolver::setSolverIterations(unsigned int iterations) { + solverIterations = std::max(1u, iterations); +} + +unsigned int PBDSolver::getSolverIterations() const { + return solverIterations; +} + +void PBDSolver::setVelocitySleepThreshold(float threshold) { + velocitySleepThreshold = std::max(0.0f, threshold); +} + +float PBDSolver::getVelocitySleepThreshold() const { + return velocitySleepThreshold; +} + +void PBDSolver::setStructuralStiffness(float stiffness) { + structuralStiffness = std::max(0.0f, std::min(1.0f, stiffness)); + setConstraintGroupStiffness(structuralConstraints, structuralStiffness); +} + +float PBDSolver::getStructuralStiffness() const { + return structuralStiffness; +} + +void PBDSolver::setShearStiffness(float stiffness) { + shearStiffness = std::max(0.0f, std::min(1.0f, stiffness)); + setConstraintGroupStiffness(shearConstraints, shearStiffness); +} + +float PBDSolver::getShearStiffness() const { + return shearStiffness; +} + +void PBDSolver::setBendStiffness(float stiffness) { + bendStiffness = std::max(0.0f, std::min(1.0f, stiffness)); + setConstraintGroupStiffness(bendConstraints, bendStiffness); +} + +float PBDSolver::getBendStiffness() const { + return bendStiffness; +} + +void PBDSolver::setSelfCollisionStiffness(float stiffness) { + selfCollisionStiffness = std::max(0.0f, std::min(1.0f, stiffness)); +} + +float PBDSolver::getSelfCollisionStiffness() const { + return selfCollisionStiffness; +} + +float PBDSolver::getSelfCollisionThickness() const { + return selfCollisionThickness; +} + +void PBDSolver::setMaxSelfCollisionContactsPerVertex(unsigned int maxContacts) { + maxSelfCollisionContactsPerVertex = std::max(1u, maxContacts); +} + +unsigned int PBDSolver::getMaxSelfCollisionContactsPerVertex() const { + return maxSelfCollisionContactsPerVertex; +} + +void PBDSolver::setWindConfig(const PBDWindConfig& config) { + windConfig = config; + if (windConfig.windDirection.squaredNorm() <= 1e-12f) { + windConfig.windDirection = Vector3f(1.0f, 0.0f, 0.0f); + } + else { + windConfig.windDirection.normalize(); + } + + windConfig.baseSpeed = std::max(0.0f, windConfig.baseSpeed); + windConfig.baseAcceleration = std::max(0.0f, windConfig.baseAcceleration); + windConfig.gustAmplitude = std::max(0.0f, windConfig.gustAmplitude); + windConfig.gustFrequency = std::max(0.0f, windConfig.gustFrequency); + windConfig.noiseStrength = std::max(0.0f, windConfig.noiseStrength); + windConfig.dragCoefficient = std::max(0.0f, windConfig.dragCoefficient); + windConfig.liftCoefficient = std::max(0.0f, std::min(2.0f, windConfig.liftCoefficient)); + windConfig.airDensity = std::max(0.0f, windConfig.airDensity); + windConfig.maxWindSpeed = std::max(0.0f, windConfig.maxWindSpeed); + + if (windConfig.inputMode == PBDWindInputMode::Disabled) { + windSpeedState = 0.0f; + } + else if (windConfig.inputMode == PBDWindInputMode::Speed) { + windSpeedState = std::min(windConfig.baseSpeed, windConfig.maxWindSpeed); + } + else { + windSpeedState = 0.0f; + } +} + +const PBDWindConfig& PBDSolver::getWindConfig() const { + return windConfig; +} + +// ----------------------------- +// 8. Debug / Diagnostics +// ----------------------------- + +void PBDSolver::updateSelfCollisionDebugStats(const std::vector& positions, bool afterProjection) { + if (!afterProjection) { + selfCollisionDebugStats.generatedContacts = static_cast(generatedSelfCollisionConstraints.size()); + selfCollisionDebugStats.initiallyViolatedContacts = 0u; + selfCollisionDebugStats.maxInitialPenetration = 0.0f; + } + else { + selfCollisionDebugStats.remainingViolatedContacts = 0u; + selfCollisionDebugStats.maxRemainingPenetration = 0.0f; + } + + for (const CollisionConstraint* constraint : generatedSelfCollisionConstraints) { + if (constraint == nullptr) continue; + const float value = constraint->evaluate(positions); + if (value >= 0.0f) continue; + + const float penetration = -value; + if (!afterProjection) { + ++selfCollisionDebugStats.initiallyViolatedContacts; + selfCollisionDebugStats.maxInitialPenetration = std::max(selfCollisionDebugStats.maxInitialPenetration, penetration); + } + else { + ++selfCollisionDebugStats.remainingViolatedContacts; + selfCollisionDebugStats.maxRemainingPenetration = std::max(selfCollisionDebugStats.maxRemainingPenetration, penetration); + } + } +} + +// ----------------------------- +// 9. Rendering / Buffer Synchronization +// ----------------------------- + +void PBDSolver::writeBackToVBuff() { + const unsigned int n = system->n_points; + + for (unsigned int i = 0; i < n; ++i) { + vbuff[3 * i + 0] = x[i][0]; // x + vbuff[3 * i + 1] = x[i][1]; // y + vbuff[3 * i + 2] = x[i][2]; // z + } +} + +void PBDSolver::commitPositions() { + const unsigned int n = system->n_points; + for (unsigned int i = 0; i < n; ++i) { + x[i] = p[i]; + } + + writeBackToVBuff(); +} \ No newline at end of file diff --git a/ClothApp/PBDSolver.h b/ClothApp/PBDSolver.h new file mode 100644 index 0000000..3163314 --- /dev/null +++ b/ClothApp/PBDSolver.h @@ -0,0 +1,453 @@ +#pragma once +#include "MassSpringSolver.h" +#include +#include +#include +#include +#include +#include + +// ----------------------------- +// Generic PBD constraint system +// ----------------------------- +// In Muller et al., each constraint j is described by: +// - a cardinality n_j +// - a scalar constraint function C_j(.) +// - particle indices i_1 ... i_n +// - a stiffness k in [0, 1] +// - an equality or inequality type +// +// For a course project, a polymorphic base class is a good fit because each +// constraint naturally owns its own C(.) and gradient code. That keeps the +// solver loop close to the paper and avoids a large switch statement. +enum class PBDConstraintType { + Equality, + Inequality +}; + +class PBDConstraint { +protected: + typedef Eigen::Vector3f Vector3f; + + std::vector particleIndices; // involved particles {i_1, ..., i_n} + float stiffnessValue; // k in [0, 1] + PBDConstraintType constraintType; // equality or inequality + +public: + PBDConstraint( + std::vector particleIndices, + float stiffnessValue, + PBDConstraintType constraintType + ); + virtual ~PBDConstraint() = default; + + std::size_t cardinality() const; + const std::vector& indices() const; + float stiffness() const; + void setStiffness(float stiffness); + PBDConstraintType type() const; + + // Evaluate the scalar constraint C(p). + virtual float evaluate(const std::vector& positions) const = 0; + + // Return gradients [grad_{p_i1} C, ..., grad_{p_in} C]. + virtual void gradients( + const std::vector& positions, + std::vector& outGradients + ) const = 0; + + // Equality constraints project toward C = 0. + // Inequality constraints project only when violated, i.e. C < 0. + virtual bool isViolated(const std::vector& positions, float epsilon) const; + + // Generic PBD projection step: + // s = C / sum_j w_j * |grad_{p_j} C|^2 + // Delta p_i = -k' * s * w_i * grad_{p_i} C + virtual void project( + std::vector& positions, + const std::vector& invMass, + unsigned int solverIterations, + float epsilon + ) const; +}; + +// Distance constraint: +// C(p1, p2) = |p1 - p2| - d +class DistanceConstraint : public PBDConstraint { +private: + float restLength; + +public: + DistanceConstraint(unsigned int i, unsigned int j, float restLength, float stiffness); + + virtual float evaluate(const std::vector& positions) const override; + virtual void gradients( + const std::vector& positions, + std::vector& outGradients + ) const override; +}; + +// Dihedral bending constraint: +// C_bend = theta - theta_0 +// where theta is the dihedral angle between the two triangles incident to a +// shared cloth edge. Unlike a distance-based bend proxy, this directly +// measures folding and is therefore much less coupled to stretching. +class DihedralBendConstraint : public PBDConstraint { +private: + float restAngle; + +public: + DihedralBendConstraint( + unsigned int edge0, + unsigned int edge1, + unsigned int opposite0, + unsigned int opposite1, + float restAngle, + float stiffness + ); + + float angle(const std::vector& positions, bool* valid = nullptr) const; + + virtual float evaluate(const std::vector& positions) const override; + virtual void gradients( + const std::vector& positions, + std::vector& outGradients + ) const override; +}; + +// Fixed point constraint: +// C(p_i) = p_i - p_fixed +// This is vector-valued in full generality, so for a practical cloth solver we +// keep one object and project it directly to the target position. +class FixedPointConstraint : public PBDConstraint { +private: + Eigen::Vector3f fixedPosition; + +public: + FixedPointConstraint(unsigned int i, const Eigen::Vector3f& fixedPosition, float stiffness = 1.0f); + + void setFixedPosition(const Eigen::Vector3f& position); + const Eigen::Vector3f& position() const; + + virtual float evaluate(const std::vector& positions) const override; + virtual void gradients( + const std::vector& positions, + std::vector& outGradients + ) const override; + virtual void project( + std::vector& positions, + const std::vector& invMass, + unsigned int solverIterations, + float epsilon + ) const override; +}; + +struct SphereCollider { + Eigen::Vector3f center; + float radius; +}; + +struct PlaneCollider { + Eigen::Vector3f point; + Eigen::Vector3f normal; +}; + +struct BoxCollider { + Eigen::Vector3f center; + Eigen::Vector3f halfExtents; +}; + +// Generated collision constraint: +// - sphere contact: C(p_i) = |p_i - c| - r +// - plane/static contact: C(p_i) = (p_i - q_c) . n_c +// - self vertex-triangle: C(q, p1, p2, p3) = (q - p1) . n - h +// Both are inequality constraints and are satisfied when C >= 0. +class CollisionConstraint : public PBDConstraint { +private: + enum class CollisionKind { + Sphere, + Plane, + EdgePlane, + SelfVertexTriangle + }; + + CollisionKind collisionKind; + Eigen::Vector3f center; + float radius; + float offset; + Eigen::Vector3f planeNormal; + Eigen::Vector2f edgeSampleWeights; + Eigen::Vector3f selfCollisionNormal; + Eigen::Vector3f selfCollisionBarycentric; + + bool selfCollisionGeometry( + const std::vector& positions, + Eigen::Vector3f& outNormal, + Eigen::Vector3f& outBarycentric + ) const; + +public: + CollisionConstraint(unsigned int i, const Eigen::Vector3f& center, float radius, float stiffness = 1.0f); + CollisionConstraint( + unsigned int i, + const Eigen::Vector3f& planePoint, + const Eigen::Vector3f& planeNormal, + float stiffness = 1.0f + ); + CollisionConstraint( + unsigned int edge0, + unsigned int edge1, + const Eigen::Vector2f& weights, + const Eigen::Vector3f& planePoint, + const Eigen::Vector3f& planeNormal, + float stiffness = 1.0f + ); + CollisionConstraint( + unsigned int vertex, + unsigned int p1, + unsigned int p2, + unsigned int p3, + float thickness, + const Eigen::Vector3f& normal, + const Eigen::Vector3f& barycentric, + float stiffness = 1.0f + ); + + virtual float evaluate(const std::vector& positions) const override; + virtual void gradients( + const std::vector& positions, + std::vector& outGradients + ) const override; +}; + +// ----------------------------- +// PBD system struct +// ----------------------------- +struct pbd_system { + typedef std::pair Edge; + typedef std::vector EdgeList; + + unsigned int n_points; + unsigned int n_constraints; + float time_step; + + EdgeList spring_list; + Eigen::VectorXf rest_lengths; + Eigen::VectorXf masses; + std::vector triangle_indices; +}; + +struct SelfCollisionDebugStats { + unsigned int generatedContacts = 0u; + unsigned int initiallyViolatedContacts = 0u; + unsigned int remainingViolatedContacts = 0u; + float maxInitialPenetration = 0.0f; + float maxRemainingPenetration = 0.0f; +}; + +enum class PBDWindInputMode { + Disabled, + Speed, + Acceleration +}; + +struct PBDWindConfig { + PBDWindInputMode inputMode = PBDWindInputMode::Disabled; + Eigen::Vector3f windDirection = Eigen::Vector3f(1.0f, 0.0f, 0.0f); + float baseSpeed = 0.0f; + float baseAcceleration = 0.0f; + float gustAmplitude = 0.0f; + float gustFrequency = 1.35f; + float noiseStrength = 0.0f; + float dragCoefficient = 1.15f; + float liftCoefficient = 0.35f; + float airDensity = 1.225f; + float maxWindSpeed = 15.0f; +}; + +// ----------------------------- +// PBD Solver +// ----------------------------- +class PBDSolver : public FixedPointController { +private: + typedef Eigen::Vector3f Vector3f; // 3D vector type + typedef std::pair Edge; + typedef std::unique_ptr ConstraintPtr; + typedef std::vector ConstraintList; + + // system / render buffer + pbd_system* system; // pointer to PBD system + float* vbuff; + + // particle state + std::vector restX; // reference positions used to preserve the original cloth side in self-collision + std::vector x; // current positions x_i + std::vector p; // predicted positions p_i + std::vector v; // velocities v_i + std::vector invMass; // inverse mass w_i = 1 / m_i + + // Persistent constraints are part of the cloth model and exist every frame. + ConstraintList persistentConstraints; + std::unordered_map fixedPointConstraints; + std::vector structuralConstraints; + std::vector shearConstraints; + std::vector bendConstraints; + + // Collision primitives persist, but actual collision constraints are generated + // fresh each step from the predicted positions x -> p. + std::vector sphereColliders; + std::vector boxColliders; + std::vector planeColliders; + ConstraintList generatedCollisionConstraints; + std::vector generatedSelfCollisionConstraints; + std::vector planeContactPoints; + std::vector planeContactNormals; + std::vector planeContactSignedDistances; + std::vector> meshAdjacency; + SelfCollisionDebugStats selfCollisionDebugStats; + + // simulation parameters + unsigned int solverIterations; + float dampingFactor; + float collisionEps; + float structuralStiffness; + float shearStiffness; + float bendStiffness; + float planeFriction; + float selfCollisionThickness; + float selfCollisionStiffness; + float selfCollisionCellSize; + unsigned int maxSelfCollisionContactsPerVertex; + float velocitySleepThreshold; + Vector3f gravity; + PBDWindConfig windConfig; + float simulationTime; + float windSpeedState; + + // ----------------------------- + // 1. Constructor / Initialization + // ----------------------------- + void initializeState(); + + // ----------------------------- + // 3. Core PBD Pipeline Stages + // ----------------------------- + void applyExternalForces(float dt); + void applyAerodynamicForces(float dt); + void dampVelocities(); + void predictPositions(float dt); + void projectConstraints(const ConstraintList& constraints); + void updateVelocities(float dt); + bool isWindEnabled() const; + float triangleFaceArea(unsigned int i0, unsigned int i1, unsigned int i2) const; + Vector3f triangleFaceNormal(unsigned int i0, unsigned int i1, unsigned int i2) const; + Vector3f averageTriangleVelocity(unsigned int i0, unsigned int i1, unsigned int i2) const; + float currentWindSpeed(const Vector3f& triangleCenter, float dt); + float proceduralWindNoise(const Vector3f& triangleCenter) const; + + // ----------------------------- + // 5. Collision Handling + // ----------------------------- + void generateCollisionConstraints(); + void applyPlaneContactVelocityDamping(); + + // ----------------------------- + // 6. Self-Collision Utilities / Helpers + // ----------------------------- + void generateSelfCollisionConstraints(); + + // ----------------------------- + // 4. Persistent Constraint Setup / Constraint Construction + // ----------------------------- + void addDistanceConstraints( + const std::vector& indices, + float stiffness, + std::vector* constraintGroup = nullptr + ); + void addDihedralBendConstraints(float stiffness); + void setConstraintGroupStiffness(const std::vector& constraints, float stiffness); + + // ----------------------------- + // 8. Debug / Diagnostics + // ----------------------------- + void updateSelfCollisionDebugStats(const std::vector& positions, bool afterProjection); + + // ----------------------------- + // 9. Rendering / Buffer Synchronization + // ----------------------------- + void writeBackToVBuff(); + void commitPositions(); + +public: + // ----------------------------- + // 1. Constructor / Initialization + // ----------------------------- + PBDSolver(pbd_system* system, float* vbuff); + + // ----------------------------- + // 2. Main Simulation Entry / Solver Loop + // ----------------------------- + void step(float dt); + void solve(unsigned int n); + + // ----------------------------- + // 4. Persistent Constraint Setup / Constraint Construction + // ----------------------------- + void pinPoint(unsigned int i); + virtual void fixPoint(unsigned int i) override; + virtual void releasePoint(unsigned int i) override; + void addStructuralConstraints(const std::vector& indices, float stiffness = 1.0f); + void addShearConstraints(const std::vector& indices, float stiffness = 1.0f); + void addBendConstraints(const std::vector& indices, float stiffness = 1.0f); + + // ----------------------------- + // 5. Collision Handling + // ----------------------------- + void addSphereCollider(const Vector3f& center, float radius); + void addBoxCollider(const Vector3f& center, const Vector3f& halfExtents); + void addPlaneCollider(const Vector3f& point, const Vector3f& normal); + + // ----------------------------- + // 7. Parameter / Tuning Interface + // ----------------------------- + void setGravity(float gravityMagnitude); + float getGravity() const; + void setDampingFactor(float damping); + float getDampingFactor() const; + void setPlaneFriction(float friction); + float getPlaneFriction() const; + void setSolverIterations(unsigned int iterations); + unsigned int getSolverIterations() const; + void setVelocitySleepThreshold(float threshold); + float getVelocitySleepThreshold() const; + void setSelfCollisionThickness(float thickness) { + if (thickness <= 0.0f) return; + selfCollisionThickness = std::max(thickness, collisionEps); + } + void setStructuralStiffness(float stiffness); + float getStructuralStiffness() const; + void setShearStiffness(float stiffness); + float getShearStiffness() const; + void setBendStiffness(float stiffness); + float getBendStiffness() const; + void setSelfCollisionStiffness(float stiffness); + float getSelfCollisionStiffness() const; + float getSelfCollisionThickness() const; + void setMaxSelfCollisionContactsPerVertex(unsigned int maxContacts); + unsigned int getMaxSelfCollisionContactsPerVertex() const; + void setWindConfig(const PBDWindConfig& config); + const PBDWindConfig& getWindConfig() const; + + // ----------------------------- + // 8. Debug / Diagnostics + // ----------------------------- + const SelfCollisionDebugStats& getSelfCollisionDebugStats() const { return selfCollisionDebugStats; } + + // ----------------------------- + // 9. Rendering / Buffer Synchronization + // ----------------------------- + // accessors + std::vector& getPositions() { return x; } + std::vector& getPredictedPositions() { return p; } + std::vector& getVelocities() { return v; } +}; \ No newline at end of file diff --git a/ClothApp/Shader.cpp b/ClothApp/Shader.cpp index 2eb3722..7cdcf77 100644 --- a/ClothApp/Shader.cpp +++ b/ClothApp/Shader.cpp @@ -16,8 +16,9 @@ GLShader::~GLShader() { void GLShader::compile(const char* source) { GLint compiled = 0; // Compiled flag + const GLint sourceLen = static_cast(std::strlen(source)); const char *ptrs[] = { source }; - const GLint lens[] = { std::strlen(source) }; + const GLint lens[] = { sourceLen }; glShaderSource(handle, 1, ptrs, lens); glCompileShader(handle); glGetShaderiv(handle, GL_COMPILE_STATUS, &compiled); @@ -32,13 +33,23 @@ void GLShader::compile(const char* source) { } void GLShader::compile(std::ifstream& source) { + if (!source) { + throw std::runtime_error("Failed to read shader source."); + } + std::vector text; source.seekg(0, std::ios_base::end); std::streampos fileSize = source.tellg(); - text.resize(fileSize); + if (fileSize <= 0) { + throw std::runtime_error("Shader source is empty."); + } + text.resize(static_cast(fileSize) + 1, '\0'); source.seekg(0, std::ios_base::beg); source.read(&text[0], fileSize); + if (!source) { + throw std::runtime_error("Failed to read shader source."); + } compile(&text[0]); } @@ -154,6 +165,9 @@ void PhongShader::postLink() { uAlbedo = glGetUniformLocation(handle, "uAlbedo"); uAmbient = glGetUniformLocation(handle, "uAmbient"); uLight = glGetUniformLocation(handle, "uLight"); + uSpecularStrength = glGetUniformLocation(handle, "uSpecularStrength"); + uShininess = glGetUniformLocation(handle, "uShininess"); + uUseFlagPattern = glGetUniformLocation(handle, "uUseFlagPattern"); } void PhongShader::setAlbedo(const glm::vec3& albedo) { assert(uAlbedo >= 0); @@ -173,6 +187,38 @@ void PhongShader::setLight(const glm::vec3& light) { glUniform3f(uLight, light[0], light[1], light[2]); glUseProgram(0); } +void PhongShader::setSpecularStrength(float strength) { + assert(uSpecularStrength >= 0); + glUseProgram(*this); + glUniform1f(uSpecularStrength, strength); + glUseProgram(0); +} +void PhongShader::setShininess(float shininess) { + assert(uShininess >= 0); + glUseProgram(*this); + glUniform1f(uShininess, shininess); + glUseProgram(0); +} + +void PhongShader::setUseFlagPattern(bool useFlagPattern) { + assert(uUseFlagPattern >= 0); + glUseProgram(*this); + glUniform1i(uUseFlagPattern, useFlagPattern ? 1 : 0); + glUseProgram(0); +} + +ShadowShader::ShadowShader() : GLProgram() {} + +void ShadowShader::postLink() { + uShadowColor = glGetUniformLocation(handle, "uShadowColor"); +} + +void ShadowShader::setShadowColor(const glm::vec4& color) { + assert(uShadowColor >= 0); + glUseProgram(*this); + glUniform4f(uShadowColor, color[0], color[1], color[2], color[3]); + glUseProgram(0); +} PickShader::PickShader() : GLProgram() {} diff --git a/ClothApp/Shader.h b/ClothApp/Shader.h index 652eb7f..2cc85b3 100644 --- a/ClothApp/Shader.h +++ b/ClothApp/Shader.h @@ -64,8 +64,8 @@ class ProgramInput : public NonCopyable { class PhongShader : public GLProgram { private: - // Albedo | Ambient Light | Light Direction - GLuint uAlbedo, uAmbient, uLight; + // Albedo | Ambient Light | Light Direction | Specular Strength | Shininess + GLuint uAlbedo, uAmbient, uLight, uSpecularStrength, uShininess, uUseFlagPattern; public: PhongShader(); @@ -73,6 +73,19 @@ class PhongShader : public GLProgram { void setAlbedo(const glm::vec3& albedo); void setAmbient(const glm::vec3& ambient); void setLight(const glm::vec3& light); + void setSpecularStrength(float strength); + void setShininess(float shininess); + void setUseFlagPattern(bool useFlagPattern); +}; + +class ShadowShader : public GLProgram { +private: + GLuint uShadowColor; + +public: + ShadowShader(); + virtual void postLink(); + void setShadowColor(const glm::vec4& color); }; diff --git a/ClothApp/UserInteraction.cpp b/ClothApp/UserInteraction.cpp index c66e241..fe3260d 100644 --- a/ClothApp/UserInteraction.cpp +++ b/ClothApp/UserInteraction.cpp @@ -3,7 +3,7 @@ #include #include -UserInteraction::UserInteraction(Renderer* renderer, CgPointFixNode* fixer, float* vbuff) +UserInteraction::UserInteraction(Renderer* renderer, FixedPointController* fixer, float* vbuff) : renderer(renderer), vbuff(vbuff), fixer(fixer), i(-1) {} void UserInteraction::setModelview(const glm::mat4& mv) { renderer->setModelview(mv); } @@ -32,12 +32,18 @@ void UserInteraction::releasePoint() { if (i == -1) return; fixer->releasePoint( void UserInteraction::movePoint(vec3 v) { if (i == -1) return; fixer->releasePoint(i); - for(int j = 0; j < 3; j++) - vbuff[3 * i + j] += v[j]; + vec3 target( + vbuff[3 * i + 0] + v[0], + vbuff[3 * i + 1] + v[1], + vbuff[3 * i + 2] + v[2] + ); + for(int j = 0; j < 3; j++) { + vbuff[3 * i + j] = target[j]; + } fixer->fixPoint(i); } -GridMeshUI::GridMeshUI(Renderer* renderer, CgPointFixNode* fixer, float* vbuff, unsigned int n) +GridMeshUI::GridMeshUI(Renderer* renderer, FixedPointController* fixer, float* vbuff, unsigned int n) : UserInteraction(renderer, fixer, vbuff), n(n) {} int GridMeshUI::colorToIndex(color c) const { diff --git a/ClothApp/UserInteraction.h b/ClothApp/UserInteraction.h index 93ae39e..ef570b1 100644 --- a/ClothApp/UserInteraction.h +++ b/ClothApp/UserInteraction.h @@ -11,12 +11,12 @@ class UserInteraction { int i; // index of fixed point float* vbuff; // vertex buffer - CgPointFixNode* fixer; // point fixer + FixedPointController* fixer; // point fixer Renderer* renderer; // pick shader renderer virtual int colorToIndex(color c) const = 0; public: - UserInteraction(Renderer* renderer, CgPointFixNode* fixer, float* vbuff); + UserInteraction(Renderer* renderer, FixedPointController* fixer, float* vbuff); void setModelview(const glm::mat4& mv); void setProjection(const glm::mat4& p); @@ -31,5 +31,5 @@ class GridMeshUI : public UserInteraction { const unsigned int n; // grid width virtual int colorToIndex(color c) const; public: - GridMeshUI(Renderer* renderer, CgPointFixNode* fixer, float* vbuff, unsigned int n); + GridMeshUI(Renderer* renderer, FixedPointController* fixer, float* vbuff, unsigned int n); }; \ No newline at end of file diff --git a/ClothApp/app.cpp b/ClothApp/app.cpp index c842879..c93a9bc 100644 --- a/ClothApp/app.cpp +++ b/ClothApp/app.cpp @@ -1,60 +1,127 @@ +// ============================================================ +// 1. Includes and Global Constants +// ============================================================ + #include #include #include #include +#include +#include +#include #include #include #include +#include #include +#include #include "Shader.h" #include "Mesh.h" #include "Renderer.h" #include "MassSpringSolver.h" +#include "PBDSolver.h" #include "UserInteraction.h" -// G L O B A L S /////////////////////////////////////////////////////////////////// +// ============================================================ +// 2. Global Runtime State +// ============================================================ -// Window -static int g_windowWidth = 640, g_windowHeight = 640; +// Window and pointer input state. +static int g_windowWidth = 1280, g_windowHeight = 1280; static bool g_mouseClickDown = false; static bool g_mouseLClickButton, g_mouseRClickButton, g_mouseMClickButton; static int g_mouseClickX; static int g_mouseClickY; -// User Interaction +// Shared UI picking and reusable collider visuals. static UserInteraction* UI; static Renderer* g_pickRenderer; +static ProgramInput* g_floor_target; +static ProgramInput* g_sphere_target; +static unsigned int g_sphere_index_count = 0u; +static ProgramInput* g_cube_target; +static unsigned int g_cube_index_count = 0u; -// Constants +// Rendering constants and shared scene colors. static const float PI = glm::pi(); - -// Shader Handles +static const glm::vec3 g_floor_albedo(0.55f, 0.55f, 0.58f); +static const glm::vec3 g_floor_ambient(0.04f, 0.04f, 0.04f); +static const glm::vec3 g_sphere_albedo(0.36f, 0.38f, 0.40f); +static const glm::vec3 g_sphere_ambient(0.05f, 0.05f, 0.05f); +static const glm::vec3 g_cube_albedo(0.36f, 0.38f, 0.40f); +static const glm::vec3 g_cube_ambient(0.05f, 0.05f, 0.05f); +static const float g_specular_strength = 0.28f; +static const float g_shininess = 32.0f; +static const glm::vec4 g_shadow_color(0.0f, 0.0f, 0.0f, 0.32f); +static const float g_floor_collision_height = -1.75f; +static const float g_floor_render_offset = -0.002f; +static const float g_floor_extent = 3.5f; + +// Shader handles. static PhongShader* g_phongShader; // linked phong shader +static ShadowShader* g_shadowShader; // linked shadow shader static PickShader* g_pickShader; // linked pick shader -// Shader parameters +// Shader material/light parameters. static const glm::vec3 g_albedo(0.0f, 0.3f, 0.7f); static const glm::vec3 g_ambient(0.01f, 0.01f, 0.01f); static const glm::vec3 g_light(1.0f, 1.0f, -1.0f); -// Mesh +// Mesh and render buffer state. static Mesh* g_clothMesh; // halfedge data structure -// Render Target static ProgramInput* g_render_target; // vertex, normal, texutre, index -// Animation +// Animation timing for the legacy mass-spring path. static const int g_fps = 60; // frames per second | 60 static const int g_iter = 5; // iterations per time step | 10 static const int g_frame_time = 15; // approximate time for frame calculations | 15 static const int g_animation_timer = (int) ((1.0f / g_fps) * 1000 - g_frame_time); -// Mass Spring System +// Mass-spring runtime objects. static mass_spring_system* g_system; static MassSpringSolver* g_solver; -// System parameters +// PBD runtime objects and shared debug flags. +static pbd_system* g_pbdSystem; +static PBDSolver* g_pbdSolver; +static float g_selfCollisionThicknessOverride = -1.0f; +static unsigned int g_pbdFrameCounter = 0u; +static bool g_enableDebugDiagnostics = false; + +enum class WindCliMode { + None, + Speed +}; + +// CLI overrides and startup argument state. +static WindCliMode g_windCliMode = WindCliMode::None; +static float g_windSpeed = 0.0f; +static bool g_windDirectionCliProvided = false; +static Eigen::Vector3f g_windDirection(1.0f, 0.0f, 0.0f); +static unsigned int g_pbdHangIterationsOverride = 0u; +static unsigned int g_pbdDropIterationsOverride = 0u; +static unsigned int g_pbdFloorIterationsOverride = 0u; +static unsigned int g_pbdDualIterationsOverride = 0u; +static float g_pbdFloorTimestepOverride = -1.0f; + +// Constraint graph for the legacy mass-spring path. +static CgRootNode* g_cgRootNode; + +// Camera / scene framing parameters. +static const float g_camera_distance = 4.2f; +static const float g_flag_camera_distance = 7.2f; + +// Scene matrices. +static glm::mat4 g_ModelViewMatrix; +static glm::mat4 g_ProjectionMatrix; + +// ============================================================ +// 3. Demo Parameter Namespaces +// ============================================================ + +// System parameters for fast-mass-spring namespace SystemParam { static const int n = 33; // must be odd, n * n = n_vertices | 61 static const float w = 2.0f; // width | 2.0f @@ -66,37 +133,450 @@ namespace SystemParam { static const float g = 9.8f * m; // gravitational force | 9.8f } -// Constraint Graph -static CgRootNode* g_cgRootNode; +// System parameters for PBD +namespace PBDSystemParam { + static const int n = 33; // must be odd, n * n = n_vertices + static const float w = 2.0f; // cloth width + static const float h = 0.008f; // time step + static const float r = w / (n - 1); // rest length + static const float m = 0.25f / (n * n); // point mass + static const float g = 9.8f; // gravitational acceleration + + static const int n_iter = 20; // solver iterations | 15 + static const float a = 0.02f; // damping factor + static const float eps = 1e-4f; // collision epsilon + static const float k_stretch = 0.9f; // stretch stiffness | 0.9f + static const float k_shear = 0.8f; // shear stiffness | 0.8f + static const float k_bend = 0.01f; // bend stiffness | 0.01f + static const float sphere_radius = 0.64f; // radius of sphere collider in drop demo | 0.64f +} -// Scene parameters -static const float g_camera_distance = 4.2f; +namespace PBDFloorDemoParam { + static const float h = 0.003f; // Smaller time step for floor-contact stability. + static const int n_iter = 28; // More solver iterations to resolve floor/self contacts robustly. + static const float selfCollisionStiffness = 0.45f; // Extra stiffness for self-collision corrections in floor demos. + static const unsigned int maxSelfCollisionContactsPerVertex = 4u; // Cap on generated self-collision contacts per vertex. +} -// Scene matrices -static glm::mat4 g_ModelViewMatrix; -static glm::mat4 g_ProjectionMatrix; +namespace PBDDebugParam { + static const unsigned int debugPrintPeriod = 20u; +} + +namespace PBDHangCliParam { + static const unsigned int minIterations = 1u; + static const unsigned int maxIterations = 80u; +} + +namespace PBDHangControlParam { + static const float stretchMin = 0.0f; + static const float stretchMax = 1.0f; + static const float shearMin = 0.0f; + static const float shearMax = 1.0f; + static const float bendMin = 0.0f; + static const float bendMax = 0.10f; + static const float dampingMin = 0.0f; + static const float dampingMax = 0.20f; + static const float stretchStep = 0.05f; + static const float shearStep = 0.05f; + static const float bendStep = 0.005f; + static const float dampingStep = 0.01f; + static const float defaultDamping = 0.07f; +} + +namespace PBDDropCliParam { + static const float minRadius = 0.1f; + static const float maxRadius = 1.5f; + static const unsigned int minIterations = 1u; + static const unsigned int maxIterations = 80u; +} + +namespace PBDDropControlParam { + static const float stretchMin = 0.0f; + static const float stretchMax = 1.0f; + static const float shearMin = 0.0f; + static const float shearMax = 1.0f; + static const float bendMin = 0.0f; + static const float bendMax = 0.10f; + static const float dampingMin = 0.0f; + static const float dampingMax = 0.20f; + static const float stretchStep = 0.05f; + static const float shearStep = 0.05f; + static const float bendStep = 0.005f; + static const float dampingStep = 0.01f; + static const float defaultDamping = 0.07f; + // Keep the drop demo on odd grid resolutions because the cloth builder and + // spring layout are authored for odd n values. Adjusting by 2 preserves that. + static const unsigned int minMeshResolution = 17u; + static const unsigned int maxMeshResolution = 65u; + static const unsigned int meshResolutionStep = 2u; + static const unsigned int defaultMeshResolution = 33u; +} + +namespace PBDFloorCliParam { + static const unsigned int minIterations = 1u; + static const unsigned int maxIterations = 80u; + static const float minTimestep = 0.001f; + static const float maxTimestep = 0.01f; +} + +namespace PBDFloorControlParam { + static const float selfCollisionStiffnessMin = 0.0f; + static const float selfCollisionStiffnessMax = 1.0f; + static const float selfCollisionStiffnessStep = 0.05f; + static const unsigned int minSelfCollisionContacts = 1u; + static const unsigned int maxSelfCollisionContacts = 12u; + static const float selfCollisionThicknessMin = 0.001f; + static const float selfCollisionThicknessMax = 0.03f; + static const float selfCollisionThicknessStep = 0.001f; + static const float floorFrictionMin = 0.0f; + static const float floorFrictionMax = 0.8f; + static const float floorFrictionStep = 0.05f; + static const float bendMin = 0.0f; + static const float bendMax = 0.05f; + static const float bendStep = 0.0025f; + static const float dampingMin = 0.0f; + static const float dampingMax = 0.20f; + static const float dampingStep = 0.01f; + static const std::array meshChoices = { 17u, 25u, 33u, 41u, 55u, 65u }; + static const unsigned int defaultMeshResolution = 33u; + static const std::array speedPresets = { 0.003f, 0.0045f, 0.006f }; + static const std::array speedPresetLabels = { "Stable", "Balanced", "Fast" }; +}; + +namespace PBDDualCliParam { + static const unsigned int minIterations = 1u; + static const unsigned int maxIterations = 80u; +} + +namespace PBDDualControlParam { + static const float stretchMin = 0.0f; + static const float stretchMax = 1.0f; + static const float shearMin = 0.0f; + static const float shearMax = 1.0f; + static const float bendMin = 0.0f; + static const float bendMax = 0.10f; + static const float stretchStep = 0.05f; + static const float shearStep = 0.05f; + static const float bendStep = 0.005f; + static const std::array meshChoices = { 17u, 25u, 33u, 41u, 55u, 65u }; + static const unsigned int defaultMeshResolution = 33u; + static const float defaultSphereRadius = 0.30f; + static const float defaultCubeSize = 0.44f; + static const float sphereRadiusMin = 0.15f; + static const float sphereRadiusMax = 0.80f; + static const float cubeSizeMin = 0.15f; + static const float cubeSizeMax = 0.80f; + static const float objectSizeStep = 0.05f; +} + +namespace PBDWindCliParam { + static const float minValue = 0.0f; // Safe lower bound for user-provided --wind-speed. + static const float maxValue = 15.0f; // Safe upper bound for user-provided --wind-speed. + static const glm::vec3 windDirection(1.0f, 0.0f, 0.0f); // wind direction + static const float minDirectionNorm = 1e-4f; + static const float gustFraction = 0.30f; // Gust amplitude scale, used in A_gust = gustFraction * userWindValue. + static const float gustFrequency = 1.35f; // Gust frequency in u(t) = u_base + A_gust * sin(2 * pi * gustFrequency * t) + noise(t, x). + static const float noiseFraction = 0.08f; // Small procedural flutter scale added to u(t) after the sinusoidal gust term. + static const float dragCoefficient = 1.15f; // Drag coefficient C_D in F_drag = 0.5 * rho * C_D * A * |v_rel|^2 * exposure. + static const float liftCoefficient = 0.35f; // Lift coefficient C_L in F_lift = 0.5 * rho * C_L * A * |v_rel|^2 * orientationTerm. + static const float airDensity = 1.225f; // Air density rho used by the drag-only aerodynamic force. + static const float dampingFactor = 0.06f; // Velocity damping for the wind demo to keep flutter stable. + static const float bendStiffness = 0.008f; // Softer bending stiffness so the flag can ripple under gusts. + static const float flagHeightOffset = 1.15f; // Vertical lift applied when rotating the cloth into the flag pose. +} + +namespace PBDWindControlParam { + static const float speedMin = 0.0f; + static const float speedMax = 15.0f; + static const float speedStep = 0.5f; + static const float directionMin = -1.0f; + static const float directionMax = 1.0f; + static const float directionStep = 0.1f; + static const float verticalWarningThreshold = 0.35f; +} -// F U N C T I O N S ////////////////////////////////////////////////////////////// -// state initialization +// ============================================================ +// 4. Runtime State Structs +// ============================================================ + +struct PBDHangRuntimeState { + bool initialized = false; + bool paused = false; + float defaultStretch = PBDSystemParam::k_stretch; + float defaultShear = PBDSystemParam::k_shear; + float defaultBend = PBDSystemParam::k_bend; + float defaultDamping = PBDHangControlParam::defaultDamping; + std::vector initialPositions; +}; + +static PBDHangRuntimeState g_pbdHangRuntime; + +struct PBDWindRuntimeState { + bool initialized = false; + bool paused = false; + float currentWindSpeed = 0.0f; + float startupWindSpeed = 0.0f; + Eigen::Vector3f currentDirectionInput = Eigen::Vector3f(1.0f, 0.0f, 0.0f); + Eigen::Vector3f appliedWindDirection = Eigen::Vector3f(1.0f, 0.0f, 0.0f); +}; + +static PBDWindRuntimeState g_pbdWindRuntime; + +struct PBDDropRuntimeState { + bool initialized = false; + bool paused = false; + float defaultStretch = PBDSystemParam::k_stretch; + float defaultShear = PBDSystemParam::k_shear; + float defaultBend = PBDSystemParam::k_bend; + float defaultDamping = PBDDropControlParam::defaultDamping; + unsigned int currentMeshResolution = PBDDropControlParam::defaultMeshResolution; + unsigned int pendingMeshResolution = PBDDropControlParam::defaultMeshResolution; + float sphereRadius = PBDSystemParam::sphere_radius; + unsigned int solverIterations = static_cast(PBDSystemParam::n_iter); +}; + +static PBDDropRuntimeState g_pbdDropRuntime; +static float g_pbdDropStartupSphereRadius = PBDSystemParam::sphere_radius; + +struct PBDFloorRuntimeState { + bool initialized = false; + bool paused = false; + float defaultSelfCollisionStiffness = PBDFloorDemoParam::selfCollisionStiffness; + unsigned int defaultMaxSelfCollisionContacts = PBDFloorDemoParam::maxSelfCollisionContactsPerVertex; + float defaultSelfCollisionThickness = 0.0f; + float defaultFloorFriction = 0.15f; + float defaultBend = PBDSystemParam::k_bend; + float defaultDamping = PBDSystemParam::a; + unsigned int solverIterations = static_cast(PBDFloorDemoParam::n_iter); + float currentTimestep = PBDFloorDemoParam::h; + int speedPresetIndex = 0; + bool customTimestep = false; + unsigned int currentMeshResolution = PBDFloorControlParam::defaultMeshResolution; + unsigned int pendingMeshResolution = PBDFloorControlParam::defaultMeshResolution; + bool debugEnabled = false; +}; + +static PBDFloorRuntimeState g_pbdFloorRuntime; + +struct PBDDualRuntimeState { + bool initialized = false; + bool paused = false; + float defaultStretch = PBDSystemParam::k_stretch; + float defaultShear = PBDSystemParam::k_shear; + float defaultBend = PBDSystemParam::k_bend; + float currentTimestep = PBDFloorDemoParam::h; + int speedPresetIndex = 0; + bool customTimestep = false; + unsigned int currentMeshResolution = PBDDualControlParam::defaultMeshResolution; + unsigned int pendingMeshResolution = PBDDualControlParam::defaultMeshResolution; + float currentSphereRadius = PBDDualControlParam::defaultSphereRadius; + float pendingSphereRadius = PBDDualControlParam::defaultSphereRadius; + float currentCubeSize = PBDDualControlParam::defaultCubeSize; + float pendingCubeSize = PBDDualControlParam::defaultCubeSize; + unsigned int solverIterations = static_cast(PBDFloorDemoParam::n_iter); +}; + +static PBDDualRuntimeState g_pbdDualRuntime; + +// ============================================================ +// 5. Forward Declarations +// ============================================================ + +// CLI parsing and validation. static void initGlutState(int, char**); static void initGLState(); - +static void parseSimMode(int, char**); +static void parseOptionalArgs(int argc, char** argv, int startIndex); +static void validateParsedOptions(); +static const char* usageMessage(); +static bool tryParseSingleTokenMode(int argc, char** argv); +static bool tryParseTwoTokenMode(int argc, char** argv); +static bool parseSharedOption(int argc, char** argv, int& index); +static bool parseWindOption(int argc, char** argv, int& index); +static bool parseIterationOption(int argc, char** argv, int& index); +static bool parseFloorOption(int argc, char** argv, int& index); +static bool parseDropOption(int argc, char** argv, int& index); + +// OpenGL / GLUT initialization and asset setup. static void initShaders(); // Read, compile and link shaders static void initCloth(); // Generate cloth mesh +static void initFloor(); // Generate floor mesh +static void initSphereColliderVisual(float radius, const glm::vec3& center); // Generate sphere collider mesh +static void initCubeColliderVisual(const glm::vec3& center, const glm::vec3& halfExtents); // Generate cube collider mesh static void initScene(); // Generate scene matrices +static void initMouseInteraction(FixedPointController*, unsigned int); +static void rebuildClothMesh(unsigned int resolution); +static void clearMouseButtonState(); + +// Wind helpers. +static Eigen::Vector3f normalizedWindDirectionOrFallback(const Eigen::Vector3f& direction, const Eigen::Vector3f& fallback, bool* usedFallback = nullptr); +static PBDWindConfig makeWindConfig(float windSpeed, const Eigen::Vector3f& direction); +static void configurePBDWindSolver(float windSpeed, const Eigen::Vector3f& directionInput, bool captureDefaults); +static void applyPBDWindSettings(float windSpeed, const Eigen::Vector3f& directionInput); +static void resetPBDWindDemo(bool resetParameters); +static void logPBDWindControlState(const std::string& reason); + +// Hang helpers. +static void logPBDHangControlState(const std::string& reason); + +// Drop helpers. +static unsigned int previousPBDDropResolution(unsigned int resolution); +static unsigned int nextPBDDropResolution(unsigned int resolution); +static void configurePBDDropSolver(float stretch, float shear, float bend, float damping, unsigned int iterations, float sphereRadius, unsigned int meshResolution, bool captureDefaults); +static void resetPBDDropDemo(bool resetParameters); +static void logPBDDropControlState(const std::string& reason); + +// Drop-floor helpers. +static unsigned int previousPBDFloorResolution(unsigned int resolution); +static unsigned int nextPBDFloorResolution(unsigned int resolution); +static int closestPBDFloorSpeedPreset(float dt); +static void setPBDFloorRuntimeTimestep(float dt, int presetIndex, bool custom); +static void configurePBDFloorSolver( + float selfCollisionStiffness, + unsigned int maxContactsPerVertex, + float selfCollisionThickness, + float floorFriction, + float bend, + float damping, + unsigned int iterations, + float timestep, + unsigned int meshResolution, + bool captureDefaults +); +static void resetPBDFloorDemo(bool resetParameters); +static void logPBDFloorControlState(const std::string& reason); + +// Drop-floor-dual helpers. +static unsigned int previousPBDDualResolution(unsigned int resolution); +static unsigned int nextPBDDualResolution(unsigned int resolution); +static void setPBDDualRuntimeTimestep(float dt, int presetIndex, bool custom); +static void configurePBDDualSolver( + float stretch, + float shear, + float bend, + unsigned int iterations, + float timestep, + unsigned int meshResolution, + float sphereRadius, + float cubeSize, + bool captureDefaults +); +static void resetPBDDualDemo(bool resetParameters); +static void logPBDDualControlState(const std::string& reason); + +// Overlay / text helpers. +static void drawBitmapText(float x, float y, const std::string& text); +static void drawWindDirectionIndicator(const Eigen::Vector3f& normalizedDirection); +static void drawPBDHangOverlay(); +static void drawPBDWindOverlay(); +static void drawPBDDropOverlay(); +static void drawPBDFloorOverlay(); +static void drawPBDDualOverlay(); +static glm::mat4 floorShadowMatrix(float planeHeight, const glm::vec3& lightDirection); +static void orientClothForFloorDrop(); +static void orientClothFlatForDualFloorDrop(float sphereRadius, float cubeSize); +static void orientClothForWindFlag(); +static void logPBDSelfCollisionDiagnostics(); +static bool isPBDMode(); +static bool hasSphereColliderVisual(); +static bool hasCubeColliderVisual(); +static unsigned int activeGridSize(); +static float activeClothWidth(); +static pbd_system* buildPBDSystem(const mass_spring_system& system); +static bool handlePBDWindKey(unsigned char key); +static bool handlePBDHangKey(unsigned char key); +static bool handlePBDDropKey(unsigned char key); +static bool handlePBDFloorKey(unsigned char key); +static bool handlePBDDualKey(unsigned char key); + +struct AnalyticBoxDefinition { + glm::vec3 center; + glm::vec3 halfExtents; +}; // demos -static void demo_hang(); // curtain hanging from top corners -static void demo_drop(); // curtain dropping on sphere -static void(*g_demo)() = demo_drop; +namespace PBDDualObstacleDemoParam { + static const glm::vec3 sphereCenter(0.42f, 0.0f, g_floor_collision_height + 0.30f); + static const float sphereRadius = PBDDualControlParam::defaultSphereRadius; + static const AnalyticBoxDefinition cube = { + glm::vec3(-0.42f, 0.0f, g_floor_collision_height + 0.44f), + glm::vec3(PBDDualControlParam::defaultCubeSize, PBDDualControlParam::defaultCubeSize, PBDDualControlParam::defaultCubeSize) + }; +} + +enum class SimMode { + MassSpringHang, + MassSpringDrop, + PBDHang, + PBDHangWind, + PBDDrop, + PBDDropFloor, + PBDDropFloorDual +}; + +static SimMode g_mode = SimMode::MassSpringHang; // default to mass-spring hanging demo, switch to other demos later +// demos +static void demo_hang(); +static void demo_drop(); +static void demo_pbd_hang(); +static void demo_pbd_hang_wind(); +static void demo_pbd_drop(); +static void demo_pbd_drop_floor(); +static void demo_pbd_drop_floor_dual(); +static void(*g_demo)() = demo_hang; + +static bool isFloorDemo() { + return g_mode == SimMode::PBDDropFloor || g_mode == SimMode::PBDDropFloorDual; +} + +static bool hasSceneFloor() { + return isFloorDemo() || g_mode == SimMode::PBDHangWind; +} + +static bool hasSphereColliderVisual() { + return g_mode == SimMode::MassSpringDrop || g_mode == SimMode::PBDDrop || g_mode == SimMode::PBDDropFloorDual; +} + +static bool hasCubeColliderVisual() { + return g_mode == SimMode::PBDDropFloorDual; +} + +static void selectDemo() { + switch (g_mode) { + case SimMode::PBDDropFloorDual: + g_demo = demo_pbd_drop_floor_dual; + break; + case SimMode::MassSpringHang: + g_demo = demo_hang; + break; + case SimMode::MassSpringDrop: + g_demo = demo_drop; + break; + case SimMode::PBDHang: + g_demo = demo_pbd_hang; + break; + case SimMode::PBDHangWind: + g_demo = demo_pbd_hang_wind; + break; + case SimMode::PBDDrop: + g_demo = demo_pbd_drop; + break; + case SimMode::PBDDropFloor: + g_demo = demo_pbd_drop_floor; + break; + } +} // glut callbacks static void display(); static void reshape(int, int); +static void keyboard(unsigned char, int, int); static void mouse(int, int, int, int); static void motion(int, int); // draw cloth function +static void drawFloor(); +static void drawFloorShadows(); static void drawCloth(); static void animateCloth(int value); @@ -112,15 +592,30 @@ void checkGlErrors(); -// M A I N ////////////////////////////////////////////////////////////////////////// +// ============================================================ +// 6. CLI Parsing and Validation +// ============================================================ + +static const char* usageMessage() { + return "Usage: ./fast-mass-spring [mass-spring|ms] [hang|drop] [--self-thickness value] [--debug], ./fast-mass-spring pbd hang [--iters value] [--self-thickness value] [--debug], ./fast-mass-spring pbd hang-wind [--wind-speed value] [--wind-dir x y z] [--self-thickness value] [--debug], ./fast-mass-spring pbd drop [--radius value] [--iters value] [--self-thickness value] [--debug], ./fast-mass-spring pbd drop-floor [--iters value] [--dt value] [--self-thickness value] [--debug], ./fast-mass-spring pbd drop-floor-dual [--iters value] [--self-thickness value] [--debug], or ./fast-mass-spring [ms-hang|ms-drop|pbd-hang|pbd-hang-wind|pbd-drop|pbd-drop-floor|pbd-drop-floor-dual] [--self-thickness value] [--debug] [--wind-speed value] [--wind-dir x y z] [--iters value]"; +} + +// ============================================================ +// 7. OpenGL / GLUT Initialization +// ============================================================ + int main(int argc, char** argv) { try { + parseSimMode(argc, argv); + validateParsedOptions(); initGlutState(argc, argv); glewInit(); initGLState(); + selectDemo(); initShaders(); initCloth(); + initFloor(); initScene(); glutTimerFunc(g_animation_timer, animateCloth, 0); @@ -135,8 +630,297 @@ int main(int argc, char** argv) { } } +static bool tryParseSingleTokenMode(int argc, char** argv) { + if (argc <= 1) return false; + + const std::string arg1(argv[1]); + if (arg1 == "mass-spring-hang" || arg1 == "ms-hang") { + g_mode = SimMode::MassSpringHang; + parseOptionalArgs(argc, argv, 2); + return true; + } + if (arg1 == "mass-spring-drop" || arg1 == "ms-drop") { + g_mode = SimMode::MassSpringDrop; + parseOptionalArgs(argc, argv, 2); + return true; + } + if (arg1 == "pbd-hang") { + g_mode = SimMode::PBDHang; + parseOptionalArgs(argc, argv, 2); + return true; + } + if (arg1 == "pbd-hang_wind" || arg1 == "pbd-hang-wind") { + g_mode = SimMode::PBDHangWind; + parseOptionalArgs(argc, argv, 2); + return true; + } + if (arg1 == "pbd-drop") { + g_mode = SimMode::PBDDrop; + parseOptionalArgs(argc, argv, 2); + return true; + } + if (arg1 == "pbd-drop-floor") { + g_mode = SimMode::PBDDropFloor; + parseOptionalArgs(argc, argv, 2); + return true; + } + if (arg1 == "pbd-drop-floor-dual") { + g_mode = SimMode::PBDDropFloorDual; + parseOptionalArgs(argc, argv, 2); + return true; + } + + return false; +} + +static bool tryParseTwoTokenMode(int argc, char** argv) { + if (argc < 3) return false; + + const std::string solver(argv[1]); + const std::string scene(argv[2]); + if ((solver == "mass-spring" || solver == "ms") && scene == "hang") { + g_mode = SimMode::MassSpringHang; + parseOptionalArgs(argc, argv, 3); + return true; + } + if ((solver == "mass-spring" || solver == "ms") && scene == "drop") { + g_mode = SimMode::MassSpringDrop; + parseOptionalArgs(argc, argv, 3); + return true; + } + if (solver == "pbd" && scene == "hang") { + g_mode = SimMode::PBDHang; + parseOptionalArgs(argc, argv, 3); + return true; + } + if (solver == "pbd" && (scene == "hang_wind" || scene == "hang-wind")) { + g_mode = SimMode::PBDHangWind; + parseOptionalArgs(argc, argv, 3); + return true; + } + if (solver == "pbd" && scene == "drop") { + g_mode = SimMode::PBDDrop; + parseOptionalArgs(argc, argv, 3); + return true; + } + if (solver == "pbd" && scene == "drop-floor") { + g_mode = SimMode::PBDDropFloor; + parseOptionalArgs(argc, argv, 3); + return true; + } + if (solver == "pbd" && scene == "drop-floor-dual") { + g_mode = SimMode::PBDDropFloorDual; + parseOptionalArgs(argc, argv, 3); + return true; + } + + return false; +} + +static void parseSimMode(int argc, char** argv) { + if (argc <= 1) return; + if (tryParseSingleTokenMode(argc, argv)) return; + if (tryParseTwoTokenMode(argc, argv)) return; + throw std::runtime_error(usageMessage()); +} + +static bool parseSharedOption(int argc, char** argv, int& index) { + const std::string arg(argv[index]); + if (arg == "--debug") { + g_enableDebugDiagnostics = true; + return true; + } + + if (arg != "--self-thickness") { + return false; + } + if (index + 1 >= argc) { + throw std::runtime_error("Missing value after --self-thickness"); + } + + std::stringstream valueStream(argv[++index]); + float thickness = -1.0f; + valueStream >> thickness; + if (!valueStream || !valueStream.eof() || thickness <= 0.0f) { + throw std::runtime_error("--self-thickness expects a positive float value"); + } + + g_selfCollisionThicknessOverride = thickness; + return true; +} + +static bool parseWindOption(int argc, char** argv, int& index) { + const std::string arg(argv[index]); + if (arg == "--wind-speed") { + if (g_mode != SimMode::PBDHangWind) { + throw std::runtime_error("--wind-speed is only valid for the pbd hang-wind demo"); + } + if (index + 1 >= argc) { + throw std::runtime_error("Missing value after --wind-speed"); + } + if (g_windCliMode != WindCliMode::None) { + throw std::runtime_error("Specify --wind-speed only once"); + } + + std::stringstream valueStream(argv[++index]); + float windSpeed = 0.0f; + valueStream >> windSpeed; + if (!valueStream || !valueStream.eof() + || windSpeed < PBDWindCliParam::minValue + || windSpeed > PBDWindCliParam::maxValue) { + throw std::runtime_error("--wind-speed expects a float value in [0, 15]"); + } + + g_windCliMode = WindCliMode::Speed; + g_windSpeed = windSpeed; + return true; + } + + if (arg != "--wind-dir") { + return false; + } + if (g_mode != SimMode::PBDHangWind) { + throw std::runtime_error("--wind-dir is only valid for the pbd hang-wind demo"); + } + if (index + 3 >= argc) { + throw std::runtime_error("--wind-dir expects three float values: x y z"); + } + if (g_windDirectionCliProvided) { + throw std::runtime_error("Specify --wind-dir only once"); + } + + std::stringstream xStream(argv[++index]); + std::stringstream yStream(argv[++index]); + std::stringstream zStream(argv[++index]); + float x = 0.0f; + float y = 0.0f; + float z = 0.0f; + xStream >> x; + yStream >> y; + zStream >> z; + if (!xStream || !xStream.eof() || !yStream || !yStream.eof() || !zStream || !zStream.eof()) { + throw std::runtime_error("--wind-dir expects three float values: x y z"); + } + + const Eigen::Vector3f direction(x, y, z); + if (direction.norm() < PBDWindCliParam::minDirectionNorm) { + throw std::runtime_error("--wind-dir expects a non-zero vector"); + } + + g_windDirection = direction; + g_windDirectionCliProvided = true; + return true; +} + +static bool parseIterationOption(int argc, char** argv, int& index) { + const std::string arg(argv[index]); + if (arg != "--iters") { + return false; + } + if (g_mode != SimMode::PBDHang && g_mode != SimMode::PBDDrop && g_mode != SimMode::PBDDropFloor && g_mode != SimMode::PBDDropFloorDual) { + throw std::runtime_error("--iters is only valid for the pbd hang, pbd drop, pbd drop-floor, or pbd drop-floor-dual demo"); + } + if (index + 1 >= argc) { + throw std::runtime_error("Missing value after --iters"); + } + if ((g_mode == SimMode::PBDHang && g_pbdHangIterationsOverride != 0u) + || (g_mode == SimMode::PBDDrop && g_pbdDropIterationsOverride != 0u) + || (g_mode == SimMode::PBDDropFloor && g_pbdFloorIterationsOverride != 0u) + || (g_mode == SimMode::PBDDropFloorDual && g_pbdDualIterationsOverride != 0u)) { + throw std::runtime_error("Specify --iters only once"); + } + + std::stringstream valueStream(argv[++index]); + int iterations = 0; + valueStream >> iterations; + if (!valueStream || !valueStream.eof() + || iterations < static_cast(PBDHangCliParam::minIterations) + || iterations > static_cast(PBDHangCliParam::maxIterations)) { + throw std::runtime_error("--iters expects an integer value in [1, 80]"); + } + + if (g_mode == SimMode::PBDHang) { + g_pbdHangIterationsOverride = static_cast(iterations); + } + else if (g_mode == SimMode::PBDDrop) { + g_pbdDropIterationsOverride = static_cast(iterations); + } + else if (g_mode == SimMode::PBDDropFloor) { + g_pbdFloorIterationsOverride = static_cast(iterations); + } + else { + g_pbdDualIterationsOverride = static_cast(iterations); + } + return true; +} + +static bool parseFloorOption(int argc, char** argv, int& index) { + const std::string arg(argv[index]); + if (arg != "--dt") { + return false; + } + if (g_mode != SimMode::PBDDropFloor) { + throw std::runtime_error("--dt is only valid for the pbd drop-floor demo"); + } + if (index + 1 >= argc) { + throw std::runtime_error("Missing value after --dt"); + } + + std::stringstream valueStream(argv[++index]); + float dt = 0.0f; + valueStream >> dt; + if (!valueStream || !valueStream.eof() + || dt < PBDFloorCliParam::minTimestep + || dt > PBDFloorCliParam::maxTimestep) { + throw std::runtime_error("--dt expects a float value in [0.001, 0.01]"); + } + + g_pbdFloorTimestepOverride = dt; + return true; +} + +static bool parseDropOption(int argc, char** argv, int& index) { + const std::string arg(argv[index]); + if (arg != "--radius") { + return false; + } + if (g_mode != SimMode::PBDDrop) { + throw std::runtime_error("--radius is only valid for the pbd drop demo"); + } + if (index + 1 >= argc) { + throw std::runtime_error("Missing value after --radius"); + } + + std::stringstream valueStream(argv[++index]); + float radius = 0.0f; + valueStream >> radius; + if (!valueStream || !valueStream.eof() + || radius < PBDDropCliParam::minRadius + || radius > PBDDropCliParam::maxRadius) { + throw std::runtime_error("--radius expects a float value in [0.1, 1.5]"); + } + + g_pbdDropStartupSphereRadius = radius; + return true; +} + +static void parseOptionalArgs(int argc, char** argv, int startIndex) { + for (int i = startIndex; i < argc; ++i) { + if (parseSharedOption(argc, argv, i)) continue; + if (parseWindOption(argc, argv, i)) continue; + if (parseIterationOption(argc, argv, i)) continue; + if (parseFloorOption(argc, argv, i)) continue; + if (parseDropOption(argc, argv, i)) continue; + throw std::runtime_error("Unknown argument: " + std::string(argv[i])); + } +} + +static void validateParsedOptions() { + if (g_mode == SimMode::PBDHangWind && g_windCliMode == WindCliMode::None) { + throw std::runtime_error("The pbd hang-wind demo requires --wind-speed with a value in [0, 15]"); + } +} -// S T A T E I N I T I A L I Z A T O N ///////////////////////////////////////////// static void initGlutState(int argc, char** argv) { glutInit(&argc, argv); glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH); @@ -145,6 +929,7 @@ static void initGlutState(int argc, char** argv) { glutDisplayFunc(display); glutReshapeFunc(reshape); + glutKeyboardFunc(keyboard); glutMouseFunc(mouse); glutMotionFunc(motion); } @@ -162,39 +947,76 @@ static void initGLState() { checkGlErrors(); } +static bool isPBDMode() { + return g_mode == SimMode::PBDHang + || g_mode == SimMode::PBDHangWind + || g_mode == SimMode::PBDDrop + || g_mode == SimMode::PBDDropFloor + || g_mode == SimMode::PBDDropFloorDual; +} + +static unsigned int activeGridSize() { + return isPBDMode() ? PBDSystemParam::n : SystemParam::n; +} + +static float activeClothWidth() { + return isPBDMode() ? PBDSystemParam::w : SystemParam::w; +} + +static void rebuildClothMesh(unsigned int resolution) { + delete g_clothMesh; + g_clothMesh = nullptr; + delete g_render_target; + g_render_target = nullptr; + + MeshBuilder meshBuilder; + meshBuilder.uniformGrid(PBDSystemParam::w, static_cast(resolution)); + g_clothMesh = meshBuilder.getResult(); + + g_render_target = new ProgramInput; + g_render_target->setPositionData(g_clothMesh->vbuff(), g_clothMesh->vbuffLen()); + g_render_target->setNormalData(g_clothMesh->nbuff(), g_clothMesh->nbuffLen()); + g_render_target->setTextureData(g_clothMesh->tbuff(), g_clothMesh->tbuffLen()); + g_render_target->setIndexData(g_clothMesh->ibuff(), g_clothMesh->ibuffLen()); +} + static void initShaders() { GLShader basic_vert(GL_VERTEX_SHADER); GLShader phong_frag(GL_FRAGMENT_SHADER); + GLShader shadow_frag(GL_FRAGMENT_SHADER); GLShader pick_frag(GL_FRAGMENT_SHADER); auto ibasic = std::ifstream("./shaders/basic.vshader"); auto iphong = std::ifstream("./shaders/phong.fshader"); + auto ishadow = std::ifstream("./shaders/shadow.fshader"); auto ifrag = std::ifstream("./shaders/pick.fshader"); basic_vert.compile(ibasic); phong_frag.compile(iphong); + shadow_frag.compile(ishadow); pick_frag.compile(ifrag); g_phongShader = new PhongShader; + g_shadowShader = new ShadowShader; g_pickShader = new PickShader; g_phongShader->link(basic_vert, phong_frag); + g_shadowShader->link(basic_vert, shadow_frag); g_pickShader->link(basic_vert, pick_frag); checkGlErrors(); } static void initCloth() { - // short hand - const int n = SystemParam::n; - const float w = SystemParam::w; + const unsigned int n = activeGridSize(); + const float w = activeClothWidth(); // generate mesh MeshBuilder meshBuilder; - meshBuilder.uniformGrid(w, n); - g_clothMesh = meshBuilder.getResult(); + meshBuilder.uniformGrid(w, n); // generate uniform grid mesh with width w and n vertices per side + g_clothMesh = meshBuilder.getResult(); // halfedge data structure // fill program input - g_render_target = new ProgramInput; + g_render_target = new ProgramInput; // vertex, normal, texutre, index check Shader.h for details g_render_target->setPositionData(g_clothMesh->vbuff(), g_clothMesh->vbuffLen()); g_render_target->setNormalData(g_clothMesh->nbuff(), g_clothMesh->nbuffLen()); g_render_target->setTextureData(g_clothMesh->tbuff(), g_clothMesh->tbuffLen()); @@ -207,15 +1029,1154 @@ static void initCloth() { g_demo(); } +static void initFloor() { + const float floorRenderHeight = g_floor_collision_height + g_floor_render_offset; + const float floorPositions[] = { + -g_floor_extent, -g_floor_extent, floorRenderHeight, + g_floor_extent, -g_floor_extent, floorRenderHeight, + g_floor_extent, g_floor_extent, floorRenderHeight, + -g_floor_extent, g_floor_extent, floorRenderHeight + }; + const float floorNormals[] = { + 0.0f, 0.0f, 1.0f, + 0.0f, 0.0f, 1.0f, + 0.0f, 0.0f, 1.0f, + 0.0f, 0.0f, 1.0f + }; + const float floorTexcoords[] = { + 0.0f, 0.0f, + 1.0f, 0.0f, + 1.0f, 1.0f, + 0.0f, 1.0f + }; + unsigned int floorIndices[] = { + 0, 1, 2, + 0, 2, 3 + }; + + g_floor_target = new ProgramInput; + g_floor_target->setPositionData(const_cast(floorPositions), 12); + g_floor_target->setNormalData(const_cast(floorNormals), 12); + g_floor_target->setTextureData(const_cast(floorTexcoords), 8); + g_floor_target->setIndexData(floorIndices, 6); +} + +static void initSphereColliderVisual(float radius, const glm::vec3& center) { + const unsigned int stacks = 24u; + const unsigned int slices = 48u; + std::vector positions; + std::vector normals; + std::vector texcoords; + std::vector indices; + + positions.reserve((stacks + 1u) * (slices + 1u) * 3u); + normals.reserve((stacks + 1u) * (slices + 1u) * 3u); + texcoords.reserve((stacks + 1u) * (slices + 1u) * 2u); + indices.reserve(stacks * slices * 6u); + + for (unsigned int stack = 0; stack <= stacks; ++stack) { + const float v = static_cast(stack) / static_cast(stacks); + const float phi = PI * v; + const float sinPhi = std::sin(phi); + const float cosPhi = std::cos(phi); + + for (unsigned int slice = 0; slice <= slices; ++slice) { + const float u = static_cast(slice) / static_cast(slices); + const float theta = 2.0f * PI * u; + const float sinTheta = std::sin(theta); + const float cosTheta = std::cos(theta); + const glm::vec3 normal(sinPhi * cosTheta, sinPhi * sinTheta, cosPhi); + const glm::vec3 position = center + radius * normal; + + positions.push_back(position.x); + positions.push_back(position.y); + positions.push_back(position.z); + normals.push_back(normal.x); + normals.push_back(normal.y); + normals.push_back(normal.z); + texcoords.push_back(u); + texcoords.push_back(v); + } + } + + for (unsigned int stack = 0; stack < stacks; ++stack) { + for (unsigned int slice = 0; slice < slices; ++slice) { + const unsigned int rowStart = stack * (slices + 1u); + const unsigned int nextRowStart = (stack + 1u) * (slices + 1u); + const unsigned int topLeft = rowStart + slice; + const unsigned int topRight = topLeft + 1u; + const unsigned int bottomLeft = nextRowStart + slice; + const unsigned int bottomRight = bottomLeft + 1u; + + indices.push_back(topLeft); + indices.push_back(bottomLeft); + indices.push_back(topRight); + indices.push_back(topRight); + indices.push_back(bottomLeft); + indices.push_back(bottomRight); + } + } + + delete g_sphere_target; + g_sphere_target = new ProgramInput; + g_sphere_target->setPositionData(positions.data(), static_cast(positions.size())); + g_sphere_target->setNormalData(normals.data(), static_cast(normals.size())); + g_sphere_target->setTextureData(texcoords.data(), static_cast(texcoords.size())); + g_sphere_target->setIndexData(indices.data(), static_cast(indices.size())); + g_sphere_index_count = static_cast(indices.size()); +} + +static void initCubeColliderVisual(const glm::vec3& center, const glm::vec3& halfExtents) { + const glm::vec3 corners[8] = { + center + glm::vec3(-halfExtents.x, -halfExtents.y, -halfExtents.z), + center + glm::vec3( halfExtents.x, -halfExtents.y, -halfExtents.z), + center + glm::vec3( halfExtents.x, halfExtents.y, -halfExtents.z), + center + glm::vec3(-halfExtents.x, halfExtents.y, -halfExtents.z), + center + glm::vec3(-halfExtents.x, -halfExtents.y, halfExtents.z), + center + glm::vec3( halfExtents.x, -halfExtents.y, halfExtents.z), + center + glm::vec3( halfExtents.x, halfExtents.y, halfExtents.z), + center + glm::vec3(-halfExtents.x, halfExtents.y, halfExtents.z) + }; + const unsigned int faceCorners[6][4] = { + { 0, 3, 2, 1 }, + { 4, 5, 6, 7 }, + { 0, 4, 7, 3 }, + { 1, 2, 6, 5 }, + { 0, 1, 5, 4 }, + { 3, 7, 6, 2 } + }; + const glm::vec3 faceNormals[6] = { + glm::vec3(0.0f, 0.0f, -1.0f), + glm::vec3(0.0f, 0.0f, 1.0f), + glm::vec3(-1.0f, 0.0f, 0.0f), + glm::vec3(1.0f, 0.0f, 0.0f), + glm::vec3(0.0f, -1.0f, 0.0f), + glm::vec3(0.0f, 1.0f, 0.0f) + }; + const glm::vec2 faceTexcoords[4] = { + glm::vec2(0.0f, 0.0f), + glm::vec2(1.0f, 0.0f), + glm::vec2(1.0f, 1.0f), + glm::vec2(0.0f, 1.0f) + }; + + std::vector positions; + std::vector normals; + std::vector texcoords; + std::vector indices; + positions.reserve(24u * 3u); + normals.reserve(24u * 3u); + texcoords.reserve(24u * 2u); + indices.reserve(36u); + + for (unsigned int face = 0; face < 6u; ++face) { + const unsigned int baseIndex = static_cast(positions.size() / 3u); + for (unsigned int corner = 0; corner < 4u; ++corner) { + const glm::vec3& position = corners[faceCorners[face][corner]]; + positions.push_back(position.x); + positions.push_back(position.y); + positions.push_back(position.z); + normals.push_back(faceNormals[face].x); + normals.push_back(faceNormals[face].y); + normals.push_back(faceNormals[face].z); + texcoords.push_back(faceTexcoords[corner].x); + texcoords.push_back(faceTexcoords[corner].y); + } + + indices.push_back(baseIndex + 0u); + indices.push_back(baseIndex + 1u); + indices.push_back(baseIndex + 2u); + indices.push_back(baseIndex + 0u); + indices.push_back(baseIndex + 2u); + indices.push_back(baseIndex + 3u); + } + + delete g_cube_target; + g_cube_target = new ProgramInput; + g_cube_target->setPositionData(positions.data(), static_cast(positions.size())); + g_cube_target->setNormalData(normals.data(), static_cast(normals.size())); + g_cube_target->setTextureData(texcoords.data(), static_cast(texcoords.size())); + g_cube_target->setIndexData(indices.data(), static_cast(indices.size())); + g_cube_index_count = static_cast(indices.size()); +} + +static glm::mat4 floorShadowMatrix(float planeHeight, const glm::vec3& lightDirection) { + const glm::vec4 plane(0.0f, 0.0f, 1.0f, -planeHeight); + const glm::vec4 light(lightDirection.x, lightDirection.y, lightDirection.z, 0.0f); + const float dot = glm::dot(plane, light); + if (std::abs(dot) <= 1e-6f) { + return glm::mat4(1.0f); + } + + glm::mat4 shadow(0.0f); + for (int row = 0; row < 4; ++row) { + for (int col = 0; col < 4; ++col) { + shadow[col][row] = ((row == col) ? dot : 0.0f) - light[row] * plane[col]; + } + } + return shadow; +} + static void initScene() { + const float cameraDistance = (g_mode == SimMode::PBDHangWind) + ? g_flag_camera_distance + : g_camera_distance; + const glm::vec3 focusPoint = (g_mode == SimMode::PBDHangWind) + ? glm::vec3(0.0f, 0.0f, PBDWindCliParam::flagHeightOffset + activeClothWidth() / 4.0f) + : glm::vec3(0.0f, 0.0f, -1.0f); g_ModelViewMatrix = glm::lookAt( - glm::vec3(0.618, -0.786, 0.3f) * g_camera_distance, - glm::vec3(0.0f, 0.0f, -1.0f), + glm::vec3(0.618, -0.786, 0.3f) * cameraDistance, + focusPoint, glm::vec3(0.0f, 0.0f, 1.0f) - ) * glm::translate(glm::mat4(1), glm::vec3(0.0f, 0.0f, SystemParam::w / 4)); + ) * glm::translate(glm::mat4(1), glm::vec3(0.0f, 0.0f, activeClothWidth() / 4)); updateProjection(); } +static void orientClothForFloorDrop() { + if (g_clothMesh == nullptr || g_render_target == nullptr) return; + + const float width = activeClothWidth(); + const float halfWidth = 0.5f * width; + const float lift = g_floor_collision_height + 0.5f * width + 0.35f; + float* const positions = g_clothMesh->vbuff(); + const unsigned int vertexCount = g_clothMesh->n_vertices(); + + for (unsigned int i = 0; i < vertexCount; ++i) { + const float x = positions[3 * i + 0]; + const float y = positions[3 * i + 1]; + const float xNormalized = (halfWidth > 1e-6f) ? (x / halfWidth) : 0.0f; + const float heightNormalized = (halfWidth > 1e-6f) ? (-y / halfWidth) : 0.0f; + const float lateralOffset = 0.035f * xNormalized + 0.0125f * heightNormalized * heightNormalized; + + positions[3 * i + 0] = x; + positions[3 * i + 1] = lateralOffset; + positions[3 * i + 2] = -y + lift; + } + + g_clothMesh->request_face_normals(); + g_clothMesh->update_normals(); + g_clothMesh->release_face_normals(); + updateRenderTarget(); +} + +static void orientClothFlatForDualFloorDrop(float sphereRadius, float cubeSize) { + if (g_clothMesh == nullptr || g_render_target == nullptr) return; + + const float sphereTop = PBDDualObstacleDemoParam::sphereCenter.z + sphereRadius; + const float cubeTop = PBDDualObstacleDemoParam::cube.center.z + cubeSize; + const float targetHeight = std::max(sphereTop, cubeTop) + 0.8f; + float* const positions = g_clothMesh->vbuff(); + const unsigned int vertexCount = g_clothMesh->n_vertices(); + + for (unsigned int i = 0; i < vertexCount; ++i) { + positions[3 * i + 2] = targetHeight; + } + + g_clothMesh->request_face_normals(); + g_clothMesh->update_normals(); + g_clothMesh->release_face_normals(); + updateRenderTarget(); +} + +static void orientClothForWindFlag() { + if (g_clothMesh == nullptr || g_render_target == nullptr) return; + + float* const positions = g_clothMesh->vbuff(); + const unsigned int vertexCount = g_clothMesh->n_vertices(); + for (unsigned int i = 0; i < vertexCount; ++i) { + const float x = positions[3 * i + 0]; + const float y = positions[3 * i + 1]; + + positions[3 * i + 0] = 0.0f; + positions[3 * i + 1] = x; + positions[3 * i + 2] = y + PBDWindCliParam::flagHeightOffset; + } + + g_clothMesh->request_face_normals(); + g_clothMesh->update_normals(); + g_clothMesh->release_face_normals(); + updateRenderTarget(); +} + +static void initMouseInteraction(FixedPointController* mouseFixer, unsigned int n) { + if (UI != nullptr) { + UI->releasePoint(); + delete UI; + UI = nullptr; + } + delete g_pickRenderer; + g_pickRenderer = nullptr; + + g_pickRenderer = new Renderer(); + g_pickRenderer->setProgram(g_pickShader); + g_pickRenderer->setProgramInput(g_render_target); + g_pickRenderer->setElementCount(g_clothMesh->ibuffLen()); + g_pickShader->setTessFact(n); + UI = new GridMeshUI(g_pickRenderer, mouseFixer, g_clothMesh->vbuff(), n); + +} + +static void configurePBDHangSolver( + float stretch, + float shear, + float bend, + float damping, + unsigned int iterations, + bool captureDefaults +) { + const unsigned int n = PBDSystemParam::n; + + delete g_pbdSolver; + g_pbdSolver = nullptr; + delete g_pbdSystem; + g_pbdSystem = nullptr; + + MassSpringBuilder builder; + builder.uniformGrid( + PBDSystemParam::n, + PBDSystemParam::h, + PBDSystemParam::r, + 1.0f, + PBDSystemParam::m, + PBDSystemParam::a, + PBDSystemParam::g + ); + + mass_spring_system* temp = builder.getResult(); + g_pbdSystem = buildPBDSystem(*temp); + delete temp; + g_pbdSolver = new PBDSolver(g_pbdSystem, g_clothMesh->vbuff()); + if (g_selfCollisionThicknessOverride > 0.0f) { + g_pbdSolver->setSelfCollisionThickness(g_selfCollisionThicknessOverride); + } + + // Stretch and shear are distance-constraint stiffness controls. + g_pbdSolver->addStructuralConstraints(builder.getStructIndex(), stretch); + g_pbdSolver->addShearConstraints(builder.getShearIndex(), shear); + // Bend controls dihedral bending stiffness. + g_pbdSolver->addBendConstraints(builder.getBendIndex(), bend); + // Damping controls velocity energy decay without changing gravity. + g_pbdSolver->setDampingFactor(damping); + g_pbdSolver->setSolverIterations(iterations); + g_pbdSolver->pinPoint(0); + g_pbdSolver->pinPoint(n - 1); + g_pbdFrameCounter = 0u; + initMouseInteraction(g_pbdSolver, n); + + if (captureDefaults || !g_pbdHangRuntime.initialized) { + g_pbdHangRuntime.defaultStretch = g_pbdSolver->getStructuralStiffness(); + g_pbdHangRuntime.defaultShear = g_pbdSolver->getShearStiffness(); + g_pbdHangRuntime.defaultBend = g_pbdSolver->getBendStiffness(); + g_pbdHangRuntime.defaultDamping = g_pbdSolver->getDampingFactor(); + g_pbdHangRuntime.initialPositions.assign( + g_clothMesh->vbuff(), + g_clothMesh->vbuff() + g_clothMesh->vbuffLen() + ); + g_pbdHangRuntime.initialized = true; + } +} + +static void resetPBDHangDemo(bool resetParameters) { + if (g_mode != SimMode::PBDHang || g_clothMesh == nullptr || !g_pbdHangRuntime.initialized) return; + + const float stretch = resetParameters + ? g_pbdHangRuntime.defaultStretch + : g_pbdSolver->getStructuralStiffness(); + const float shear = resetParameters + ? g_pbdHangRuntime.defaultShear + : g_pbdSolver->getShearStiffness(); + const float bend = resetParameters + ? g_pbdHangRuntime.defaultBend + : g_pbdSolver->getBendStiffness(); + const float damping = resetParameters + ? g_pbdHangRuntime.defaultDamping + : g_pbdSolver->getDampingFactor(); + const unsigned int iterations = (g_pbdSolver != nullptr) + ? g_pbdSolver->getSolverIterations() + : ((g_pbdHangIterationsOverride > 0u) ? g_pbdHangIterationsOverride : static_cast(PBDSystemParam::n_iter)); + + if (UI != nullptr) { + UI->releasePoint(); + } + clearMouseButtonState(); + + std::copy( + g_pbdHangRuntime.initialPositions.begin(), + g_pbdHangRuntime.initialPositions.end(), + g_clothMesh->vbuff() + ); + g_clothMesh->request_face_normals(); + g_clothMesh->update_normals(); + g_clothMesh->release_face_normals(); + updateRenderTarget(); + + configurePBDHangSolver(stretch, shear, bend, damping, iterations, false); + logPBDHangControlState(resetParameters ? "reset-parameters" : "reset-cloth"); + glutPostRedisplay(); +} + +static void logPBDHangControlState(const std::string& reason) { + if (g_mode != SimMode::PBDHang || g_pbdSolver == nullptr) return; + std::cout + << "[pbd hang controls] " << reason + << " stretch=" << g_pbdSolver->getStructuralStiffness() + << ", shear=" << g_pbdSolver->getShearStiffness() + << ", bend=" << g_pbdSolver->getBendStiffness() + << ", damping=" << g_pbdSolver->getDampingFactor() + << ", iters=" << g_pbdSolver->getSolverIterations() + << (g_pbdHangRuntime.paused ? ", paused" : ", running") + << std::endl; +} + +static Eigen::Vector3f normalizedWindDirectionOrFallback(const Eigen::Vector3f& direction, const Eigen::Vector3f& fallback, bool* usedFallback) { + const float norm = direction.norm(); + if (norm < PBDWindCliParam::minDirectionNorm) { + if (usedFallback != nullptr) { + *usedFallback = true; + } + const float fallbackNorm = fallback.norm(); + if (fallbackNorm < PBDWindCliParam::minDirectionNorm) { + return Eigen::Vector3f(1.0f, 0.0f, 0.0f); + } + return fallback / fallbackNorm; + } + if (usedFallback != nullptr) { + *usedFallback = false; + } + return direction / norm; +} + +static PBDWindConfig makeWindConfig(float windSpeed, const Eigen::Vector3f& direction) { + PBDWindConfig windConfig; + windConfig.inputMode = PBDWindInputMode::Speed; + // Wind direction controls the direction of the applied aerodynamic force and + // is normalized before it is given to the solver. + windConfig.windDirection = normalizedWindDirectionOrFallback(direction, Eigen::Vector3f(1.0f, 0.0f, 0.0f), nullptr); + // Wind speed controls the magnitude of the external wind velocity and is + // clamped to [0, 15] for stability. + windConfig.baseSpeed = std::max(PBDWindControlParam::speedMin, std::min(PBDWindControlParam::speedMax, windSpeed)); + windConfig.baseAcceleration = 0.0f; + const float gustReference = windConfig.baseSpeed; + // Drag, lift, gust, and noise remain internal defaults to keep the demo simple. + windConfig.gustAmplitude = PBDWindCliParam::gustFraction * gustReference; + windConfig.gustFrequency = PBDWindCliParam::gustFrequency; + windConfig.noiseStrength = PBDWindCliParam::noiseFraction * gustReference; + windConfig.dragCoefficient = PBDWindCliParam::dragCoefficient; + windConfig.liftCoefficient = PBDWindCliParam::liftCoefficient; + windConfig.airDensity = PBDWindCliParam::airDensity; + windConfig.maxWindSpeed = PBDWindCliParam::maxValue; + return windConfig; +} + +static void applyPBDWindSettings(float windSpeed, const Eigen::Vector3f& directionInput) { + if (g_pbdSolver == nullptr) return; + + const float clampedSpeed = std::max(PBDWindControlParam::speedMin, std::min(PBDWindControlParam::speedMax, windSpeed)); + const Eigen::Vector3f normalizedDirection = normalizedWindDirectionOrFallback( + directionInput, + g_pbdWindRuntime.appliedWindDirection, + nullptr + ); + g_pbdSolver->setWindConfig(makeWindConfig(clampedSpeed, directionInput)); + g_pbdWindRuntime.currentWindSpeed = clampedSpeed; + g_pbdWindRuntime.currentDirectionInput = directionInput; + g_pbdWindRuntime.appliedWindDirection = normalizedDirection; +} + +static void configurePBDWindSolver(float windSpeed, const Eigen::Vector3f& directionInput, bool captureDefaults) { + const unsigned int n = PBDSystemParam::n; + const Eigen::Vector3f floorPoint(0.0f, 0.0f, g_floor_collision_height); + const Eigen::Vector3f floorNormal(0.0f, 0.0f, 1.0f); + + if (UI != nullptr) { + UI->releasePoint(); + } + + delete g_pbdSolver; + g_pbdSolver = nullptr; + delete g_pbdSystem; + g_pbdSystem = nullptr; + + orientClothForWindFlag(); + + MassSpringBuilder builder; + builder.uniformGrid( + PBDSystemParam::n, + PBDSystemParam::h, + PBDSystemParam::r, + 1.0f, + PBDSystemParam::m, + PBDSystemParam::a, + PBDSystemParam::g + ); + + mass_spring_system* temp = builder.getResult(); + g_pbdSystem = buildPBDSystem(*temp); + delete temp; + g_pbdSolver = new PBDSolver(g_pbdSystem, g_clothMesh->vbuff()); + if (g_selfCollisionThicknessOverride > 0.0f) { + g_pbdSolver->setSelfCollisionThickness(g_selfCollisionThicknessOverride); + } + + applyPBDWindSettings(windSpeed, directionInput); + g_pbdSolver->setDampingFactor(PBDWindCliParam::dampingFactor); + g_pbdSolver->addStructuralConstraints(builder.getStructIndex(), PBDSystemParam::k_stretch); + g_pbdSolver->addShearConstraints(builder.getShearIndex(), PBDSystemParam::k_shear); + g_pbdSolver->addBendConstraints(builder.getBendIndex(), PBDWindCliParam::bendStiffness); + g_pbdSolver->addPlaneCollider(floorPoint, floorNormal); + for (unsigned int row = 0; row < n; ++row) { + g_pbdSolver->pinPoint(row * n); + } + + g_pbdFrameCounter = 0u; + initMouseInteraction(g_pbdSolver, n); + + if (captureDefaults || !g_pbdWindRuntime.initialized) { + g_pbdWindRuntime.startupWindSpeed = g_pbdWindRuntime.currentWindSpeed; + g_pbdWindRuntime.initialized = true; + } +} + +static void resetPBDWindDemo(bool resetParameters) { + if (g_mode != SimMode::PBDHangWind || g_pbdSolver == nullptr || !g_pbdWindRuntime.initialized) return; + + if (resetParameters) { + applyPBDWindSettings( + g_pbdWindRuntime.startupWindSpeed, + Eigen::Vector3f(PBDWindCliParam::windDirection.x, PBDWindCliParam::windDirection.y, PBDWindCliParam::windDirection.z) + ); + logPBDWindControlState("reset-parameters"); + glutPostRedisplay(); + return; + } + + clearMouseButtonState(); + + configurePBDWindSolver(g_pbdWindRuntime.currentWindSpeed, g_pbdWindRuntime.currentDirectionInput, false); + logPBDWindControlState("reset-cloth"); + glutPostRedisplay(); +} + +static void logPBDWindControlState(const std::string& reason) { + if (g_mode != SimMode::PBDHangWind || g_pbdSolver == nullptr) return; + const Eigen::Vector3f& rawDirection = g_pbdWindRuntime.currentDirectionInput; + const Eigen::Vector3f& normalizedDirection = g_pbdWindRuntime.appliedWindDirection; + std::cout + << "[pbd wind controls] " << reason + << " speed=" << g_pbdWindRuntime.currentWindSpeed + << " range=[0, 15]" + << ", direction=(" << rawDirection.x() << ", " << rawDirection.y() << ", " << rawDirection.z() << ")" + << ", normalized=(" << normalizedDirection.x() << ", " << normalizedDirection.y() << ", " << normalizedDirection.z() << ")" + << (g_pbdWindRuntime.paused ? ", paused" : ", running") + << std::endl; +} + +static unsigned int previousPBDDropResolution(unsigned int resolution) { + if (resolution <= PBDDropControlParam::minMeshResolution) { + return PBDDropControlParam::minMeshResolution; + } + return std::max( + PBDDropControlParam::minMeshResolution, + resolution - PBDDropControlParam::meshResolutionStep + ); +} + +static unsigned int nextPBDDropResolution(unsigned int resolution) { + if (resolution >= PBDDropControlParam::maxMeshResolution) { + return PBDDropControlParam::maxMeshResolution; + } + return std::min( + PBDDropControlParam::maxMeshResolution, + resolution + PBDDropControlParam::meshResolutionStep + ); +} + +static void configurePBDDropSolver( + float stretch, + float shear, + float bend, + float damping, + unsigned int iterations, + float sphereRadius, + unsigned int meshResolution, + bool captureDefaults +) { + if (UI != nullptr) { + UI->releasePoint(); + } + + delete g_pbdSolver; + g_pbdSolver = nullptr; + delete g_pbdSystem; + g_pbdSystem = nullptr; + + // Mesh resolution changes particle and constraint count, so applying a new + // resolution requires rebuilding the cloth mesh, render buffers, and PBD system. + rebuildClothMesh(meshResolution); + const float totalClothMass = PBDSystemParam::m * static_cast(PBDSystemParam::n * PBDSystemParam::n); + const float pointMass = totalClothMass / static_cast(meshResolution * meshResolution); + + MassSpringBuilder builder; + builder.uniformGrid( + meshResolution, + PBDSystemParam::h, + PBDSystemParam::w / static_cast(meshResolution - 1u), + 1.0f, + pointMass, + PBDSystemParam::a, + PBDSystemParam::g + ); + + mass_spring_system* temp = builder.getResult(); + g_pbdSystem = buildPBDSystem(*temp); + delete temp; + g_pbdSolver = new PBDSolver(g_pbdSystem, g_clothMesh->vbuff()); + if (g_selfCollisionThicknessOverride > 0.0f) { + g_pbdSolver->setSelfCollisionThickness(g_selfCollisionThicknessOverride); + } + + // Stretch and shear are distance-constraint stiffness controls. + g_pbdSolver->addStructuralConstraints(builder.getStructIndex(), stretch); + g_pbdSolver->addShearConstraints(builder.getShearIndex(), shear); + // Bend controls dihedral bending stiffness. + g_pbdSolver->addBendConstraints(builder.getBendIndex(), bend); + // Damping controls velocity energy decay. + g_pbdSolver->setDampingFactor(damping); + // Iterations control the number of PBD projection passes per timestep. + g_pbdSolver->setSolverIterations(iterations); + + const glm::vec3 sphereCenter(0.0f, 0.0f, -1.0f); + // Sphere radius changes the collision equation C(p) = |p - c| - r >= 0. + g_pbdSolver->addSphereCollider(Eigen::Vector3f(sphereCenter.x, sphereCenter.y, sphereCenter.z), sphereRadius); + initSphereColliderVisual(sphereRadius, sphereCenter); + + g_pbdFrameCounter = 0u; + initMouseInteraction(g_pbdSolver, meshResolution); + + g_pbdDropRuntime.currentMeshResolution = meshResolution; + g_pbdDropRuntime.pendingMeshResolution = meshResolution; + g_pbdDropRuntime.sphereRadius = sphereRadius; + g_pbdDropRuntime.solverIterations = iterations; + + if (captureDefaults || !g_pbdDropRuntime.initialized) { + g_pbdDropRuntime.defaultStretch = g_pbdSolver->getStructuralStiffness(); + g_pbdDropRuntime.defaultShear = g_pbdSolver->getShearStiffness(); + g_pbdDropRuntime.defaultBend = g_pbdSolver->getBendStiffness(); + g_pbdDropRuntime.defaultDamping = g_pbdSolver->getDampingFactor(); + g_pbdDropRuntime.initialized = true; + } +} + +static void resetPBDDropDemo(bool resetParameters) { + if (g_mode != SimMode::PBDDrop || g_pbdSolver == nullptr || !g_pbdDropRuntime.initialized) return; + + if (resetParameters) { + g_pbdSolver->setStructuralStiffness(g_pbdDropRuntime.defaultStretch); + g_pbdSolver->setShearStiffness(g_pbdDropRuntime.defaultShear); + g_pbdSolver->setBendStiffness(g_pbdDropRuntime.defaultBend); + g_pbdSolver->setDampingFactor(g_pbdDropRuntime.defaultDamping); + g_pbdDropRuntime.pendingMeshResolution = PBDDropControlParam::defaultMeshResolution; + logPBDDropControlState("reset-parameters"); + glutPostRedisplay(); + return; + } + + const float stretch = g_pbdSolver->getStructuralStiffness(); + const float shear = g_pbdSolver->getShearStiffness(); + const float bend = g_pbdSolver->getBendStiffness(); + const float damping = g_pbdSolver->getDampingFactor(); + + clearMouseButtonState(); + + configurePBDDropSolver( + stretch, + shear, + bend, + damping, + g_pbdDropRuntime.solverIterations, + g_pbdDropRuntime.sphereRadius, + g_pbdDropRuntime.pendingMeshResolution, + false + ); + logPBDDropControlState("reset-cloth"); + glutPostRedisplay(); +} + +static void logPBDDropControlState(const std::string& reason) { + if (g_mode != SimMode::PBDDrop || g_pbdSolver == nullptr) return; + std::cout + << "[pbd drop controls] " << reason + << " radius=" << g_pbdDropRuntime.sphereRadius + << ", iters=" << g_pbdSolver->getSolverIterations() + << ", stretch=" << g_pbdSolver->getStructuralStiffness() + << ", shear=" << g_pbdSolver->getShearStiffness() + << ", bend=" << g_pbdSolver->getBendStiffness() + << ", damping=" << g_pbdSolver->getDampingFactor() + << ", mesh=" << g_pbdDropRuntime.currentMeshResolution << "x" << g_pbdDropRuntime.currentMeshResolution; + if (g_pbdDropRuntime.pendingMeshResolution != g_pbdDropRuntime.currentMeshResolution) { + std::cout << ", pending-mesh=" << g_pbdDropRuntime.pendingMeshResolution << "x" << g_pbdDropRuntime.pendingMeshResolution; + } + std::cout << (g_pbdDropRuntime.paused ? ", paused" : ", running") << std::endl; +} + +static unsigned int previousPBDFloorResolution(unsigned int resolution) { + for (std::size_t i = 0; i < PBDFloorControlParam::meshChoices.size(); ++i) { + if (PBDFloorControlParam::meshChoices[i] == resolution) { + return (i == 0u) + ? PBDFloorControlParam::meshChoices.front() + : PBDFloorControlParam::meshChoices[i - 1u]; + } + } + return PBDFloorControlParam::meshChoices.front(); +} + +static unsigned int nextPBDFloorResolution(unsigned int resolution) { + for (std::size_t i = 0; i < PBDFloorControlParam::meshChoices.size(); ++i) { + if (PBDFloorControlParam::meshChoices[i] == resolution) { + return (i + 1u >= PBDFloorControlParam::meshChoices.size()) + ? PBDFloorControlParam::meshChoices.back() + : PBDFloorControlParam::meshChoices[i + 1u]; + } + } + return PBDFloorControlParam::meshChoices.back(); +} + +static int closestPBDFloorSpeedPreset(float dt) { + int bestIndex = 0; + float bestDistance = std::abs(dt - PBDFloorControlParam::speedPresets[0]); + for (std::size_t i = 1; i < PBDFloorControlParam::speedPresets.size(); ++i) { + const float distance = std::abs(dt - PBDFloorControlParam::speedPresets[i]); + if (distance < bestDistance) { + bestDistance = distance; + bestIndex = static_cast(i); + } + } + return bestIndex; +} + +static void setPBDFloorRuntimeTimestep(float dt, int presetIndex, bool custom) { + g_pbdFloorRuntime.currentTimestep = std::max(PBDFloorCliParam::minTimestep, std::min(PBDFloorCliParam::maxTimestep, dt)); + g_pbdFloorRuntime.speedPresetIndex = presetIndex; + g_pbdFloorRuntime.customTimestep = custom; + if (g_pbdSystem != nullptr) { + g_pbdSystem->time_step = g_pbdFloorRuntime.currentTimestep; + } +} + +static void configurePBDFloorSolver( + float selfCollisionStiffness, + unsigned int maxContactsPerVertex, + float selfCollisionThickness, + float floorFriction, + float bend, + float damping, + unsigned int iterations, + float timestep, + unsigned int meshResolution, + bool captureDefaults +) { + if (UI != nullptr) { + UI->releasePoint(); + } + + delete g_pbdSolver; + g_pbdSolver = nullptr; + delete g_pbdSystem; + g_pbdSystem = nullptr; + + // Mesh resolution changes particle, constraint, and contact counts, so it + // must be applied by rebuilding the cloth system instead of changing live. + rebuildClothMesh(meshResolution); + orientClothForFloorDrop(); + + const float totalClothMass = PBDSystemParam::m * static_cast(PBDSystemParam::n * PBDSystemParam::n); + const float pointMass = totalClothMass / static_cast(meshResolution * meshResolution); + + MassSpringBuilder builder; + builder.uniformGrid( + meshResolution, + PBDSystemParam::h, + PBDSystemParam::w / static_cast(meshResolution - 1u), + 1.0f, + pointMass, + PBDSystemParam::a, + PBDSystemParam::g + ); + + mass_spring_system* temp = builder.getResult(); + g_pbdSystem = buildPBDSystem(*temp); + delete temp; + g_pbdSystem->time_step = timestep; + g_pbdSolver = new PBDSolver(g_pbdSystem, g_clothMesh->vbuff()); + g_pbdSolver->setSolverIterations(iterations); + + // Self-collision stiffness controls how strongly self-collision inequality + // constraints are projected each iteration. + g_pbdSolver->setSelfCollisionStiffness(std::max( + PBDFloorControlParam::selfCollisionStiffnessMin, + std::min(PBDFloorControlParam::selfCollisionStiffnessMax, selfCollisionStiffness) + )); + // Max contacts per vertex prevents over-constraining/conflicting contacts. + g_pbdSolver->setMaxSelfCollisionContactsPerVertex(std::max( + PBDFloorControlParam::minSelfCollisionContacts, + std::min(PBDFloorControlParam::maxSelfCollisionContacts, maxContactsPerVertex) + )); + + if (g_selfCollisionThicknessOverride > 0.0f) { + g_pbdSolver->setSelfCollisionThickness(g_selfCollisionThicknessOverride); + } + else { + // Self-collision thickness is the h term in + // C(q,p1,p2,p3) = (q - p1) . n - h >= 0. + g_pbdSolver->setSelfCollisionThickness(std::max( + PBDFloorControlParam::selfCollisionThicknessMin, + std::min(PBDFloorControlParam::selfCollisionThicknessMax, selfCollisionThickness) + )); + } + + g_pbdSolver->addStructuralConstraints(builder.getStructIndex(), PBDSystemParam::k_stretch); + g_pbdSolver->addShearConstraints(builder.getShearIndex(), PBDSystemParam::k_shear); + g_pbdSolver->addBendConstraints(builder.getBendIndex(), bend); + + const Eigen::Vector3f floorPoint(0.0f, 0.0f, g_floor_collision_height); + const Eigen::Vector3f floorNormal(0.0f, 0.0f, 1.0f); + g_pbdSolver->addPlaneCollider(floorPoint, floorNormal); + // Floor friction is velocity-level tangential damping after contact. + g_pbdSolver->setPlaneFriction(std::max( + PBDFloorControlParam::floorFrictionMin, + std::min(PBDFloorControlParam::floorFrictionMax, floorFriction) + )); + // Bend stiffness controls resistance to folding/unrolling. + g_pbdSolver->setBendStiffness(std::max( + PBDFloorControlParam::bendMin, + std::min(PBDFloorControlParam::bendMax, bend) + )); + // Damping controls velocity energy decay. + g_pbdSolver->setDampingFactor(std::max( + PBDFloorControlParam::dampingMin, + std::min(PBDFloorControlParam::dampingMax, damping) + )); + + g_pbdFrameCounter = 0u; + initMouseInteraction(g_pbdSolver, meshResolution); + + g_pbdFloorRuntime.currentMeshResolution = meshResolution; + g_pbdFloorRuntime.pendingMeshResolution = meshResolution; + g_pbdFloorRuntime.solverIterations = iterations; + g_pbdFloorRuntime.debugEnabled = g_enableDebugDiagnostics; + + const int closestPreset = closestPBDFloorSpeedPreset(timestep); + const bool customTimestep = std::abs(timestep - PBDFloorControlParam::speedPresets[closestPreset]) > 1e-5f; + setPBDFloorRuntimeTimestep(timestep, closestPreset, customTimestep); + + if (captureDefaults || !g_pbdFloorRuntime.initialized) { + g_pbdFloorRuntime.defaultSelfCollisionStiffness = g_pbdSolver->getSelfCollisionStiffness(); + g_pbdFloorRuntime.defaultMaxSelfCollisionContacts = g_pbdSolver->getMaxSelfCollisionContactsPerVertex(); + g_pbdFloorRuntime.defaultSelfCollisionThickness = g_pbdSolver->getSelfCollisionThickness(); + g_pbdFloorRuntime.defaultFloorFriction = g_pbdSolver->getPlaneFriction(); + g_pbdFloorRuntime.defaultBend = g_pbdSolver->getBendStiffness(); + g_pbdFloorRuntime.defaultDamping = g_pbdSolver->getDampingFactor(); + g_pbdFloorRuntime.initialized = true; + } +} + +static void resetPBDFloorDemo(bool resetParameters) { + if (g_mode != SimMode::PBDDropFloor || g_pbdSolver == nullptr || !g_pbdFloorRuntime.initialized) return; + + if (resetParameters) { + g_pbdSolver->setSelfCollisionStiffness(g_pbdFloorRuntime.defaultSelfCollisionStiffness); + g_pbdSolver->setMaxSelfCollisionContactsPerVertex(g_pbdFloorRuntime.defaultMaxSelfCollisionContacts); + if (g_selfCollisionThicknessOverride > 0.0f) { + g_pbdSolver->setSelfCollisionThickness(g_selfCollisionThicknessOverride); + } + else { + g_pbdSolver->setSelfCollisionThickness(g_pbdFloorRuntime.defaultSelfCollisionThickness); + } + g_pbdSolver->setPlaneFriction(g_pbdFloorRuntime.defaultFloorFriction); + g_pbdSolver->setBendStiffness(g_pbdFloorRuntime.defaultBend); + g_pbdSolver->setDampingFactor(g_pbdFloorRuntime.defaultDamping); + g_pbdFloorRuntime.pendingMeshResolution = PBDFloorControlParam::defaultMeshResolution; + setPBDFloorRuntimeTimestep(PBDFloorControlParam::speedPresets[0], 0, false); + g_enableDebugDiagnostics = g_pbdFloorRuntime.debugEnabled; + logPBDFloorControlState("reset-parameters"); + glutPostRedisplay(); + return; + } + + clearMouseButtonState(); + + configurePBDFloorSolver( + g_pbdSolver->getSelfCollisionStiffness(), + g_pbdSolver->getMaxSelfCollisionContactsPerVertex(), + g_pbdSolver->getSelfCollisionThickness(), + g_pbdSolver->getPlaneFriction(), + g_pbdSolver->getBendStiffness(), + g_pbdSolver->getDampingFactor(), + g_pbdFloorRuntime.solverIterations, + g_pbdFloorRuntime.currentTimestep, + g_pbdFloorRuntime.pendingMeshResolution, + false + ); + logPBDFloorControlState("reset-cloth"); + glutPostRedisplay(); +} + +static void logPBDFloorControlState(const std::string& reason) { + if (g_mode != SimMode::PBDDropFloor || g_pbdSolver == nullptr) return; + std::cout + << "[pbd drop-floor controls] " << reason + << " dt=" << g_pbdFloorRuntime.currentTimestep + << ", iters=" << g_pbdFloorRuntime.solverIterations + << ", self-k=" << g_pbdSolver->getSelfCollisionStiffness() + << ", max-contacts=" << g_pbdSolver->getMaxSelfCollisionContactsPerVertex() + << ", self-thickness=" << g_pbdSolver->getSelfCollisionThickness() + << ", floor-friction=" << g_pbdSolver->getPlaneFriction() + << ", bend=" << g_pbdSolver->getBendStiffness() + << ", damping=" << g_pbdSolver->getDampingFactor() + << ", mesh=" << g_pbdFloorRuntime.currentMeshResolution << "x" << g_pbdFloorRuntime.currentMeshResolution; + if (g_pbdFloorRuntime.pendingMeshResolution != g_pbdFloorRuntime.currentMeshResolution) { + std::cout << ", pending-mesh=" << g_pbdFloorRuntime.pendingMeshResolution << "x" << g_pbdFloorRuntime.pendingMeshResolution; + } + std::cout + << ", debug=" << (g_enableDebugDiagnostics ? "on" : "off") + << (g_pbdFloorRuntime.paused ? ", paused" : ", running") + << std::endl; + if (g_pbdFloorRuntime.speedPresetIndex > 0) { + std::cout << "[pbd drop-floor controls] faster timestep may increase penetration, jitter, or missed self-collision" << std::endl; + } +} + +static unsigned int previousPBDDualResolution(unsigned int resolution) { + for (std::size_t i = 0; i < PBDDualControlParam::meshChoices.size(); ++i) { + if (PBDDualControlParam::meshChoices[i] == resolution) { + return (i == 0u) + ? PBDDualControlParam::meshChoices.front() + : PBDDualControlParam::meshChoices[i - 1u]; + } + } + return PBDDualControlParam::meshChoices.front(); +} + +static unsigned int nextPBDDualResolution(unsigned int resolution) { + for (std::size_t i = 0; i < PBDDualControlParam::meshChoices.size(); ++i) { + if (PBDDualControlParam::meshChoices[i] == resolution) { + return (i + 1u >= PBDDualControlParam::meshChoices.size()) + ? PBDDualControlParam::meshChoices.back() + : PBDDualControlParam::meshChoices[i + 1u]; + } + } + return PBDDualControlParam::meshChoices.back(); +} + +static void setPBDDualRuntimeTimestep(float dt, int presetIndex, bool custom) { + g_pbdDualRuntime.currentTimestep = std::max(PBDFloorCliParam::minTimestep, std::min(PBDFloorCliParam::maxTimestep, dt)); + g_pbdDualRuntime.speedPresetIndex = presetIndex; + g_pbdDualRuntime.customTimestep = custom; + if (g_pbdSystem != nullptr) { + g_pbdSystem->time_step = g_pbdDualRuntime.currentTimestep; + } +} + +static void configurePBDDualSolver( + float stretch, + float shear, + float bend, + unsigned int iterations, + float timestep, + unsigned int meshResolution, + float sphereRadius, + float cubeSize, + bool captureDefaults +) { + if (UI != nullptr) { + UI->releasePoint(); + } + + delete g_pbdSolver; + g_pbdSolver = nullptr; + delete g_pbdSystem; + g_pbdSystem = nullptr; + + // Mesh resolution changes particle, constraint, and contact count, so it + // must be applied through reset by rebuilding the cloth and solver state. + rebuildClothMesh(meshResolution); + + const float clampedSphereRadius = std::max( + PBDDualControlParam::sphereRadiusMin, + std::min(PBDDualControlParam::sphereRadiusMax, sphereRadius) + ); + const float clampedCubeSize = std::max( + PBDDualControlParam::cubeSizeMin, + std::min(PBDDualControlParam::cubeSizeMax, cubeSize) + ); + orientClothFlatForDualFloorDrop(clampedSphereRadius, clampedCubeSize); + + const float totalClothMass = PBDSystemParam::m * static_cast(PBDSystemParam::n * PBDSystemParam::n); + const float pointMass = totalClothMass / static_cast(meshResolution * meshResolution); + + MassSpringBuilder builder; + builder.uniformGrid( + meshResolution, + PBDSystemParam::h, + PBDSystemParam::w / static_cast(meshResolution - 1u), + 1.0f, + pointMass, + PBDSystemParam::a, + PBDSystemParam::g + ); + + mass_spring_system* temp = builder.getResult(); + g_pbdSystem = buildPBDSystem(*temp); + delete temp; + g_pbdSystem->time_step = std::max(PBDFloorCliParam::minTimestep, std::min(PBDFloorCliParam::maxTimestep, timestep)); + g_pbdSolver = new PBDSolver(g_pbdSystem, g_clothMesh->vbuff()); + g_pbdSolver->setSolverIterations(iterations); + g_pbdSolver->setSelfCollisionStiffness(PBDFloorDemoParam::selfCollisionStiffness); + g_pbdSolver->setMaxSelfCollisionContactsPerVertex(PBDFloorDemoParam::maxSelfCollisionContactsPerVertex); + if (g_selfCollisionThicknessOverride > 0.0f) { + g_pbdSolver->setSelfCollisionThickness(g_selfCollisionThicknessOverride); + } + + // Stretch and shear are distance-constraint stiffness controls. + g_pbdSolver->addStructuralConstraints( + builder.getStructIndex(), + std::max(PBDDualControlParam::stretchMin, std::min(PBDDualControlParam::stretchMax, stretch)) + ); + g_pbdSolver->addShearConstraints( + builder.getShearIndex(), + std::max(PBDDualControlParam::shearMin, std::min(PBDDualControlParam::shearMax, shear)) + ); + // Bend controls dihedral bending stiffness. + g_pbdSolver->addBendConstraints( + builder.getBendIndex(), + std::max(PBDDualControlParam::bendMin, std::min(PBDDualControlParam::bendMax, bend)) + ); + + const Eigen::Vector3f floorPoint(0.0f, 0.0f, g_floor_collision_height); + const Eigen::Vector3f floorNormal(0.0f, 0.0f, 1.0f); + g_pbdSolver->addPlaneCollider(floorPoint, floorNormal); + + const glm::vec3 cubeHalfExtents(clampedCubeSize, clampedCubeSize, clampedCubeSize); + // Sphere radius changes the smooth obstacle collision size. + g_pbdSolver->addSphereCollider( + Eigen::Vector3f( + PBDDualObstacleDemoParam::sphereCenter.x, + PBDDualObstacleDemoParam::sphereCenter.y, + PBDDualObstacleDemoParam::sphereCenter.z + ), + clampedSphereRadius + ); + // Cube size changes the sharp analytic box collider and visual cube together. + g_pbdSolver->addBoxCollider( + Eigen::Vector3f( + PBDDualObstacleDemoParam::cube.center.x, + PBDDualObstacleDemoParam::cube.center.y, + PBDDualObstacleDemoParam::cube.center.z + ), + Eigen::Vector3f(cubeHalfExtents.x, cubeHalfExtents.y, cubeHalfExtents.z) + ); + // Pending object sizes are applied only on reset so visuals and colliders + // stay synchronized after the rebuild. + initSphereColliderVisual(clampedSphereRadius, PBDDualObstacleDemoParam::sphereCenter); + initCubeColliderVisual(PBDDualObstacleDemoParam::cube.center, cubeHalfExtents); + + g_pbdFrameCounter = 0u; + initMouseInteraction(g_pbdSolver, meshResolution); + + g_pbdDualRuntime.currentMeshResolution = meshResolution; + g_pbdDualRuntime.pendingMeshResolution = meshResolution; + g_pbdDualRuntime.currentSphereRadius = clampedSphereRadius; + g_pbdDualRuntime.pendingSphereRadius = clampedSphereRadius; + g_pbdDualRuntime.currentCubeSize = clampedCubeSize; + g_pbdDualRuntime.pendingCubeSize = clampedCubeSize; + g_pbdDualRuntime.solverIterations = iterations; + + const int closestPreset = closestPBDFloorSpeedPreset(g_pbdSystem->time_step); + const bool customTimestep = std::abs(g_pbdSystem->time_step - PBDFloorControlParam::speedPresets[closestPreset]) > 1e-5f; + setPBDDualRuntimeTimestep(g_pbdSystem->time_step, closestPreset, customTimestep); + + if (captureDefaults || !g_pbdDualRuntime.initialized) { + g_pbdDualRuntime.defaultStretch = g_pbdSolver->getStructuralStiffness(); + g_pbdDualRuntime.defaultShear = g_pbdSolver->getShearStiffness(); + g_pbdDualRuntime.defaultBend = g_pbdSolver->getBendStiffness(); + g_pbdDualRuntime.initialized = true; + } +} + +static void resetPBDDualDemo(bool resetParameters) { + if (g_mode != SimMode::PBDDropFloorDual || g_pbdSolver == nullptr || !g_pbdDualRuntime.initialized) return; + + if (resetParameters) { + g_pbdSolver->setStructuralStiffness(g_pbdDualRuntime.defaultStretch); + g_pbdSolver->setShearStiffness(g_pbdDualRuntime.defaultShear); + g_pbdSolver->setBendStiffness(g_pbdDualRuntime.defaultBend); + setPBDDualRuntimeTimestep(PBDFloorControlParam::speedPresets[0], 0, false); + g_pbdDualRuntime.pendingMeshResolution = PBDDualControlParam::defaultMeshResolution; + g_pbdDualRuntime.pendingSphereRadius = PBDDualControlParam::defaultSphereRadius; + g_pbdDualRuntime.pendingCubeSize = PBDDualControlParam::defaultCubeSize; + logPBDDualControlState("reset-parameters"); + glutPostRedisplay(); + return; + } + + clearMouseButtonState(); + + configurePBDDualSolver( + g_pbdSolver->getStructuralStiffness(), + g_pbdSolver->getShearStiffness(), + g_pbdSolver->getBendStiffness(), + g_pbdDualRuntime.solverIterations, + g_pbdDualRuntime.currentTimestep, + g_pbdDualRuntime.pendingMeshResolution, + g_pbdDualRuntime.pendingSphereRadius, + g_pbdDualRuntime.pendingCubeSize, + false + ); + logPBDDualControlState("reset-cloth"); + glutPostRedisplay(); +} + +static void logPBDDualControlState(const std::string& reason) { + if (g_mode != SimMode::PBDDropFloorDual || g_pbdSolver == nullptr) return; + std::cout + << "[pbd dual-obstacle controls] " << reason + << " iters=" << g_pbdDualRuntime.solverIterations + << ", dt=" << g_pbdDualRuntime.currentTimestep + << ", stretch=" << g_pbdSolver->getStructuralStiffness() + << ", shear=" << g_pbdSolver->getShearStiffness() + << ", bend=" << g_pbdSolver->getBendStiffness() + << ", mesh=" << g_pbdDualRuntime.currentMeshResolution << "x" << g_pbdDualRuntime.currentMeshResolution + << ", sphere=" << g_pbdDualRuntime.currentSphereRadius + << ", cube=" << g_pbdDualRuntime.currentCubeSize; + if (g_pbdDualRuntime.pendingMeshResolution != g_pbdDualRuntime.currentMeshResolution) { + std::cout << ", pending-mesh=" << g_pbdDualRuntime.pendingMeshResolution << "x" << g_pbdDualRuntime.pendingMeshResolution; + } + if (std::abs(g_pbdDualRuntime.pendingSphereRadius - g_pbdDualRuntime.currentSphereRadius) > 1e-5f) { + std::cout << ", pending-sphere=" << g_pbdDualRuntime.pendingSphereRadius; + } + if (std::abs(g_pbdDualRuntime.pendingCubeSize - g_pbdDualRuntime.currentCubeSize) > 1e-5f) { + std::cout << ", pending-cube=" << g_pbdDualRuntime.pendingCubeSize; + } + std::cout << (g_pbdDualRuntime.paused ? ", paused" : ", running") << std::endl; +} + +static pbd_system* buildPBDSystem(const mass_spring_system& system) { + pbd_system* pbdSystem = new pbd_system; + pbdSystem->n_points = system.n_points; + pbdSystem->n_constraints = system.n_springs; + pbdSystem->time_step = system.time_step; + pbdSystem->spring_list = system.spring_list; + pbdSystem->rest_lengths = system.rest_lengths; + pbdSystem->masses = system.masses; + if (g_clothMesh != nullptr) { + pbdSystem->triangle_indices.assign( + g_clothMesh->ibuff(), + g_clothMesh->ibuff() + g_clothMesh->ibuffLen() + ); + } + return pbdSystem; +} + static void demo_hang() { // short hand const int n = SystemParam::n; @@ -253,13 +2214,8 @@ static void demo_hang() { cornerFixer->fixPoint(n - 1); // initialize user interaction - g_pickRenderer = new Renderer(); - g_pickRenderer->setProgram(g_pickShader); - g_pickRenderer->setProgramInput(g_render_target); - g_pickRenderer->setElementCount(g_clothMesh->ibuffLen()); - g_pickShader->setTessFact(SystemParam::n); CgPointFixNode* mouseFixer = new CgPointFixNode(g_system, g_clothMesh->vbuff()); - UI = new GridMeshUI(g_pickRenderer, mouseFixer, g_clothMesh->vbuff(), n); + initMouseInteraction(mouseFixer, n); // build constraint graph g_cgRootNode = new CgRootNode(g_system, g_clothMesh->vbuff()); @@ -304,6 +2260,7 @@ static void demo_drop() { // sphere collision constraint CgSphereCollisionNode* sphereCollisionNode = new CgSphereCollisionNode(g_system, g_clothMesh->vbuff(), radius, center); + initSphereColliderVisual(radius, glm::vec3(center[0], center[1], center[2])); // spring deformation constraint CgSpringDeformationNode* deformationNode = @@ -312,13 +2269,8 @@ static void demo_drop() { deformationNode->addSprings(massSpringBuilder.getStructIndex()); // initialize user interaction - g_pickRenderer = new Renderer(); - g_pickRenderer->setProgram(g_pickShader); - g_pickRenderer->setProgramInput(g_render_target); - g_pickRenderer->setElementCount(g_clothMesh->ibuffLen()); - g_pickShader->setTessFact(SystemParam::n); CgPointFixNode* mouseFixer = new CgPointFixNode(g_system, g_clothMesh->vbuff()); - UI = new GridMeshUI(g_pickRenderer, mouseFixer, g_clothMesh->vbuff(), n); + initMouseInteraction(mouseFixer, n); // build constraint graph g_cgRootNode = new CgRootNode(g_system, g_clothMesh->vbuff()); @@ -331,10 +2283,125 @@ static void demo_drop() { deformationNode->addChild(mouseFixer); } +static void demo_pbd_hang() { + const unsigned int iterations = (g_pbdHangIterationsOverride > 0u) + ? g_pbdHangIterationsOverride + : static_cast(PBDSystemParam::n_iter); + configurePBDHangSolver( + PBDSystemParam::k_stretch, + PBDSystemParam::k_shear, + PBDSystemParam::k_bend, + PBDHangControlParam::defaultDamping, + iterations, + true + ); + logPBDHangControlState("startup"); +} + +static void demo_pbd_hang_wind() { + configurePBDWindSolver(g_windSpeed, g_windDirection, true); + logPBDWindControlState("startup"); +} + +static void demo_pbd_drop() { + const unsigned int iterations = (g_pbdDropIterationsOverride > 0u) + ? g_pbdDropIterationsOverride + : static_cast(PBDSystemParam::n_iter); + configurePBDDropSolver( + PBDSystemParam::k_stretch, + PBDSystemParam::k_shear, + PBDSystemParam::k_bend, + PBDDropControlParam::defaultDamping, + iterations, + g_pbdDropStartupSphereRadius, + PBDDropControlParam::defaultMeshResolution, + true + ); + logPBDDropControlState("startup"); +} + +static void demo_pbd_drop_floor() { + const unsigned int iterations = (g_pbdFloorIterationsOverride > 0u) + ? g_pbdFloorIterationsOverride + : static_cast(PBDFloorDemoParam::n_iter); + const float timestep = (g_pbdFloorTimestepOverride > 0.0f) + ? g_pbdFloorTimestepOverride + : PBDFloorDemoParam::h; + configurePBDFloorSolver( + PBDFloorDemoParam::selfCollisionStiffness, + PBDFloorDemoParam::maxSelfCollisionContactsPerVertex, + (g_selfCollisionThicknessOverride > 0.0f) ? g_selfCollisionThicknessOverride : PBDFloorControlParam::selfCollisionThicknessMin, + 0.15f, + PBDSystemParam::k_bend, + 0.07f, + iterations, + timestep, + PBDFloorControlParam::defaultMeshResolution, + true + ); + logPBDFloorControlState("startup"); +} + +static void demo_pbd_drop_floor_dual() { + const unsigned int iterations = (g_pbdDualIterationsOverride > 0u) + ? g_pbdDualIterationsOverride + : static_cast(PBDFloorDemoParam::n_iter); + configurePBDDualSolver( + PBDSystemParam::k_stretch, + PBDSystemParam::k_shear, + PBDSystemParam::k_bend, + iterations, + PBDFloorDemoParam::h, + PBDDualControlParam::defaultMeshResolution, + PBDDualControlParam::defaultSphereRadius, + PBDDualControlParam::defaultCubeSize, + true + ); + logPBDDualControlState("startup"); +} // G L U T C A L L B A C K S ////////////////////////////////////////////////////// static void display() { glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + if (hasSceneFloor()) { + drawFloor(); + drawFloorShadows(); + } + if (hasSphereColliderVisual() && g_sphere_target != nullptr) { + Renderer renderer; + renderer.setProgram(g_phongShader); + renderer.setModelview(g_ModelViewMatrix); + renderer.setProjection(g_ProjectionMatrix); + g_phongShader->setUseFlagPattern(false); + g_phongShader->setAlbedo(g_sphere_albedo); + g_phongShader->setAmbient(g_sphere_ambient); + g_phongShader->setLight(g_light); + g_phongShader->setSpecularStrength(g_specular_strength); + g_phongShader->setShininess(g_shininess); + renderer.setProgramInput(g_sphere_target); + renderer.setElementCount(g_sphere_index_count); + renderer.draw(); + } + if (hasCubeColliderVisual() && g_cube_target != nullptr) { + Renderer renderer; + renderer.setProgram(g_phongShader); + renderer.setModelview(g_ModelViewMatrix); + renderer.setProjection(g_ProjectionMatrix); + g_phongShader->setUseFlagPattern(false); + g_phongShader->setAlbedo(g_sphere_albedo); + g_phongShader->setAmbient(g_sphere_ambient); + g_phongShader->setLight(g_light); + g_phongShader->setSpecularStrength(g_specular_strength); + g_phongShader->setShininess(g_shininess); + renderer.setProgramInput(g_cube_target); + renderer.setElementCount(g_cube_index_count); + renderer.draw(); + } drawCloth(); + drawPBDHangOverlay(); + drawPBDWindOverlay(); + drawPBDDropOverlay(); + drawPBDFloorOverlay(); + drawPBDDualOverlay(); glutSwapBuffers(); checkGlErrors(); @@ -348,6 +2415,544 @@ static void reshape(int w, int h) { glutPostRedisplay(); } +static void clearMouseButtonState() { + g_mouseClickDown = false; + g_mouseLClickButton = false; + g_mouseRClickButton = false; + g_mouseMClickButton = false; +} + +static bool handlePBDWindKey(unsigned char key) { + switch (key) { + case '1': + applyPBDWindSettings( + g_pbdWindRuntime.currentWindSpeed - PBDWindControlParam::speedStep, + g_pbdWindRuntime.currentDirectionInput + ); + logPBDWindControlState("speed-"); + return true; + case '2': + applyPBDWindSettings( + g_pbdWindRuntime.currentWindSpeed + PBDWindControlParam::speedStep, + g_pbdWindRuntime.currentDirectionInput + ); + logPBDWindControlState("speed+"); + return true; + case '3': { + Eigen::Vector3f direction = g_pbdWindRuntime.currentDirectionInput; + direction.x() = std::max(PBDWindControlParam::directionMin, direction.x() - PBDWindControlParam::directionStep); + if (direction.norm() >= PBDWindCliParam::minDirectionNorm) { + applyPBDWindSettings(g_pbdWindRuntime.currentWindSpeed, direction); + } + logPBDWindControlState("dir-x-"); + return true; + } + case '4': { + Eigen::Vector3f direction = g_pbdWindRuntime.currentDirectionInput; + direction.x() = std::min(PBDWindControlParam::directionMax, direction.x() + PBDWindControlParam::directionStep); + if (direction.norm() >= PBDWindCliParam::minDirectionNorm) { + applyPBDWindSettings(g_pbdWindRuntime.currentWindSpeed, direction); + } + logPBDWindControlState("dir-x+"); + return true; + } + case '5': { + Eigen::Vector3f direction = g_pbdWindRuntime.currentDirectionInput; + direction.y() = std::max(PBDWindControlParam::directionMin, direction.y() - PBDWindControlParam::directionStep); + if (direction.norm() >= PBDWindCliParam::minDirectionNorm) { + applyPBDWindSettings(g_pbdWindRuntime.currentWindSpeed, direction); + } + logPBDWindControlState("dir-y-"); + return true; + } + case '6': { + Eigen::Vector3f direction = g_pbdWindRuntime.currentDirectionInput; + direction.y() = std::min(PBDWindControlParam::directionMax, direction.y() + PBDWindControlParam::directionStep); + if (direction.norm() >= PBDWindCliParam::minDirectionNorm) { + applyPBDWindSettings(g_pbdWindRuntime.currentWindSpeed, direction); + } + logPBDWindControlState("dir-y+"); + return true; + } + case '7': { + Eigen::Vector3f direction = g_pbdWindRuntime.currentDirectionInput; + direction.z() = std::max(PBDWindControlParam::directionMin, direction.z() - PBDWindControlParam::directionStep); + if (direction.norm() >= PBDWindCliParam::minDirectionNorm) { + applyPBDWindSettings(g_pbdWindRuntime.currentWindSpeed, direction); + } + logPBDWindControlState("dir-z-"); + return true; + } + case '8': { + Eigen::Vector3f direction = g_pbdWindRuntime.currentDirectionInput; + direction.z() = std::min(PBDWindControlParam::directionMax, direction.z() + PBDWindControlParam::directionStep); + if (direction.norm() >= PBDWindCliParam::minDirectionNorm) { + applyPBDWindSettings(g_pbdWindRuntime.currentWindSpeed, direction); + } + logPBDWindControlState("dir-z+"); + return true; + } + case 'r': + case 'R': + resetPBDWindDemo(false); + return true; + case 't': + case 'T': + resetPBDWindDemo(true); + return true; + case 'p': + case 'P': + g_pbdWindRuntime.paused = !g_pbdWindRuntime.paused; + logPBDWindControlState(g_pbdWindRuntime.paused ? "pause" : "resume"); + return true; + default: + return false; + } +} + +static bool handlePBDHangKey(unsigned char key) { + switch (key) { + case '1': + g_pbdSolver->setStructuralStiffness(std::max( + PBDHangControlParam::stretchMin, + g_pbdSolver->getStructuralStiffness() - PBDHangControlParam::stretchStep + )); + logPBDHangControlState("stretch-"); + return true; + case '2': + g_pbdSolver->setStructuralStiffness(std::min( + PBDHangControlParam::stretchMax, + g_pbdSolver->getStructuralStiffness() + PBDHangControlParam::stretchStep + )); + logPBDHangControlState("stretch+"); + return true; + case '3': + g_pbdSolver->setShearStiffness(std::max( + PBDHangControlParam::shearMin, + g_pbdSolver->getShearStiffness() - PBDHangControlParam::shearStep + )); + logPBDHangControlState("shear-"); + return true; + case '4': + g_pbdSolver->setShearStiffness(std::min( + PBDHangControlParam::shearMax, + g_pbdSolver->getShearStiffness() + PBDHangControlParam::shearStep + )); + logPBDHangControlState("shear+"); + return true; + case '5': + g_pbdSolver->setBendStiffness(std::max( + PBDHangControlParam::bendMin, + g_pbdSolver->getBendStiffness() - PBDHangControlParam::bendStep + )); + logPBDHangControlState("bend-"); + return true; + case '6': + g_pbdSolver->setBendStiffness(std::min( + PBDHangControlParam::bendMax, + g_pbdSolver->getBendStiffness() + PBDHangControlParam::bendStep + )); + logPBDHangControlState("bend+"); + return true; + case '7': + g_pbdSolver->setDampingFactor(std::max( + PBDHangControlParam::dampingMin, + g_pbdSolver->getDampingFactor() - PBDHangControlParam::dampingStep + )); + logPBDHangControlState("damping-"); + return true; + case '8': + g_pbdSolver->setDampingFactor(std::min( + PBDHangControlParam::dampingMax, + g_pbdSolver->getDampingFactor() + PBDHangControlParam::dampingStep + )); + logPBDHangControlState("damping+"); + return true; + case 'r': + case 'R': + resetPBDHangDemo(false); + return true; + case 't': + case 'T': + resetPBDHangDemo(true); + return true; + case 'p': + case 'P': + g_pbdHangRuntime.paused = !g_pbdHangRuntime.paused; + logPBDHangControlState(g_pbdHangRuntime.paused ? "pause" : "resume"); + return true; + default: + return false; + } +} + +static bool handlePBDDropKey(unsigned char key) { + switch (key) { + case '1': + g_pbdSolver->setStructuralStiffness(std::max( + PBDDropControlParam::stretchMin, + g_pbdSolver->getStructuralStiffness() - PBDDropControlParam::stretchStep + )); + logPBDDropControlState("stretch-"); + return true; + case '2': + g_pbdSolver->setStructuralStiffness(std::min( + PBDDropControlParam::stretchMax, + g_pbdSolver->getStructuralStiffness() + PBDDropControlParam::stretchStep + )); + logPBDDropControlState("stretch+"); + return true; + case '3': + g_pbdSolver->setShearStiffness(std::max( + PBDDropControlParam::shearMin, + g_pbdSolver->getShearStiffness() - PBDDropControlParam::shearStep + )); + logPBDDropControlState("shear-"); + return true; + case '4': + g_pbdSolver->setShearStiffness(std::min( + PBDDropControlParam::shearMax, + g_pbdSolver->getShearStiffness() + PBDDropControlParam::shearStep + )); + logPBDDropControlState("shear+"); + return true; + case '5': + g_pbdSolver->setBendStiffness(std::max( + PBDDropControlParam::bendMin, + g_pbdSolver->getBendStiffness() - PBDDropControlParam::bendStep + )); + logPBDDropControlState("bend-"); + return true; + case '6': + g_pbdSolver->setBendStiffness(std::min( + PBDDropControlParam::bendMax, + g_pbdSolver->getBendStiffness() + PBDDropControlParam::bendStep + )); + logPBDDropControlState("bend+"); + return true; + case '7': + g_pbdSolver->setDampingFactor(std::max( + PBDDropControlParam::dampingMin, + g_pbdSolver->getDampingFactor() - PBDDropControlParam::dampingStep + )); + logPBDDropControlState("damping-"); + return true; + case '8': + g_pbdSolver->setDampingFactor(std::min( + PBDDropControlParam::dampingMax, + g_pbdSolver->getDampingFactor() + PBDDropControlParam::dampingStep + )); + logPBDDropControlState("damping+"); + return true; + case '[': + g_pbdDropRuntime.pendingMeshResolution = previousPBDDropResolution(g_pbdDropRuntime.pendingMeshResolution); + logPBDDropControlState("pending-mesh-"); + return true; + case ']': + g_pbdDropRuntime.pendingMeshResolution = nextPBDDropResolution(g_pbdDropRuntime.pendingMeshResolution); + logPBDDropControlState("pending-mesh+"); + return true; + case 'r': + case 'R': + resetPBDDropDemo(false); + return true; + case 't': + case 'T': + resetPBDDropDemo(true); + return true; + case 'p': + case 'P': + g_pbdDropRuntime.paused = !g_pbdDropRuntime.paused; + logPBDDropControlState(g_pbdDropRuntime.paused ? "pause" : "resume"); + return true; + default: + return false; + } +} + +static bool handlePBDFloorKey(unsigned char key) { + switch (key) { + case '1': + g_pbdSolver->setSelfCollisionStiffness(std::max( + PBDFloorControlParam::selfCollisionStiffnessMin, + g_pbdSolver->getSelfCollisionStiffness() - PBDFloorControlParam::selfCollisionStiffnessStep + )); + logPBDFloorControlState("self-k-"); + return true; + case '2': + g_pbdSolver->setSelfCollisionStiffness(std::min( + PBDFloorControlParam::selfCollisionStiffnessMax, + g_pbdSolver->getSelfCollisionStiffness() + PBDFloorControlParam::selfCollisionStiffnessStep + )); + logPBDFloorControlState("self-k+"); + return true; + case '3': + g_pbdSolver->setMaxSelfCollisionContactsPerVertex(std::max( + PBDFloorControlParam::minSelfCollisionContacts, + g_pbdSolver->getMaxSelfCollisionContactsPerVertex() - 1u + )); + logPBDFloorControlState("max-contacts-"); + return true; + case '4': + g_pbdSolver->setMaxSelfCollisionContactsPerVertex(std::min( + PBDFloorControlParam::maxSelfCollisionContacts, + g_pbdSolver->getMaxSelfCollisionContactsPerVertex() + 1u + )); + logPBDFloorControlState("max-contacts+"); + return true; + case '5': + g_pbdSolver->setSelfCollisionThickness(std::max( + PBDFloorControlParam::selfCollisionThicknessMin, + g_pbdSolver->getSelfCollisionThickness() - PBDFloorControlParam::selfCollisionThicknessStep + )); + logPBDFloorControlState("self-thickness-"); + return true; + case '6': + g_pbdSolver->setSelfCollisionThickness(std::min( + PBDFloorControlParam::selfCollisionThicknessMax, + g_pbdSolver->getSelfCollisionThickness() + PBDFloorControlParam::selfCollisionThicknessStep + )); + logPBDFloorControlState("self-thickness+"); + return true; + case '7': + g_pbdSolver->setPlaneFriction(std::max( + PBDFloorControlParam::floorFrictionMin, + g_pbdSolver->getPlaneFriction() - PBDFloorControlParam::floorFrictionStep + )); + logPBDFloorControlState("floor-friction-"); + return true; + case '8': + g_pbdSolver->setPlaneFriction(std::min( + PBDFloorControlParam::floorFrictionMax, + g_pbdSolver->getPlaneFriction() + PBDFloorControlParam::floorFrictionStep + )); + logPBDFloorControlState("floor-friction+"); + return true; + case 'q': + case 'Q': + g_pbdSolver->setBendStiffness(std::max( + PBDFloorControlParam::bendMin, + g_pbdSolver->getBendStiffness() - PBDFloorControlParam::bendStep + )); + logPBDFloorControlState("bend-"); + return true; + case 'w': + case 'W': + g_pbdSolver->setBendStiffness(std::min( + PBDFloorControlParam::bendMax, + g_pbdSolver->getBendStiffness() + PBDFloorControlParam::bendStep + )); + logPBDFloorControlState("bend+"); + return true; + case 'a': + case 'A': + g_pbdSolver->setDampingFactor(std::max( + PBDFloorControlParam::dampingMin, + g_pbdSolver->getDampingFactor() - PBDFloorControlParam::dampingStep + )); + logPBDFloorControlState("damping-"); + return true; + case 's': + case 'S': + g_pbdSolver->setDampingFactor(std::min( + PBDFloorControlParam::dampingMax, + g_pbdSolver->getDampingFactor() + PBDFloorControlParam::dampingStep + )); + logPBDFloorControlState("damping+"); + return true; + case '[': + g_pbdFloorRuntime.pendingMeshResolution = previousPBDFloorResolution(g_pbdFloorRuntime.pendingMeshResolution); + logPBDFloorControlState("pending-mesh-"); + return true; + case ']': + g_pbdFloorRuntime.pendingMeshResolution = nextPBDFloorResolution(g_pbdFloorRuntime.pendingMeshResolution); + logPBDFloorControlState("pending-mesh+"); + return true; + case '9': { + const int nextIndex = std::max(0, g_pbdFloorRuntime.speedPresetIndex - 1); + setPBDFloorRuntimeTimestep(PBDFloorControlParam::speedPresets[nextIndex], nextIndex, false); + logPBDFloorControlState("speed-"); + return true; + } + case '0': { + const int nextIndex = std::min( + static_cast(PBDFloorControlParam::speedPresets.size()) - 1, + g_pbdFloorRuntime.speedPresetIndex + 1 + ); + setPBDFloorRuntimeTimestep(PBDFloorControlParam::speedPresets[nextIndex], nextIndex, false); + logPBDFloorControlState("speed+"); + return true; + } + case 'r': + case 'R': + resetPBDFloorDemo(false); + return true; + case 't': + case 'T': + resetPBDFloorDemo(true); + return true; + case 'p': + case 'P': + g_pbdFloorRuntime.paused = !g_pbdFloorRuntime.paused; + logPBDFloorControlState(g_pbdFloorRuntime.paused ? "pause" : "resume"); + return true; + case 'd': + case 'D': + g_enableDebugDiagnostics = !g_enableDebugDiagnostics; + g_pbdFloorRuntime.debugEnabled = g_enableDebugDiagnostics; + logPBDFloorControlState(g_enableDebugDiagnostics ? "debug-on" : "debug-off"); + return true; + default: + return false; + } +} + +static bool handlePBDDualKey(unsigned char key) { + switch (key) { + case '1': + g_pbdSolver->setStructuralStiffness(std::max( + PBDDualControlParam::stretchMin, + g_pbdSolver->getStructuralStiffness() - PBDDualControlParam::stretchStep + )); + logPBDDualControlState("stretch-"); + return true; + case '2': + g_pbdSolver->setStructuralStiffness(std::min( + PBDDualControlParam::stretchMax, + g_pbdSolver->getStructuralStiffness() + PBDDualControlParam::stretchStep + )); + logPBDDualControlState("stretch+"); + return true; + case '3': + g_pbdSolver->setShearStiffness(std::max( + PBDDualControlParam::shearMin, + g_pbdSolver->getShearStiffness() - PBDDualControlParam::shearStep + )); + logPBDDualControlState("shear-"); + return true; + case '4': + g_pbdSolver->setShearStiffness(std::min( + PBDDualControlParam::shearMax, + g_pbdSolver->getShearStiffness() + PBDDualControlParam::shearStep + )); + logPBDDualControlState("shear+"); + return true; + case '5': + g_pbdSolver->setBendStiffness(std::max( + PBDDualControlParam::bendMin, + g_pbdSolver->getBendStiffness() - PBDDualControlParam::bendStep + )); + logPBDDualControlState("bend-"); + return true; + case '6': + g_pbdSolver->setBendStiffness(std::min( + PBDDualControlParam::bendMax, + g_pbdSolver->getBendStiffness() + PBDDualControlParam::bendStep + )); + logPBDDualControlState("bend+"); + return true; + case '[': + g_pbdDualRuntime.pendingMeshResolution = previousPBDDualResolution(g_pbdDualRuntime.pendingMeshResolution); + logPBDDualControlState("pending-mesh-"); + return true; + case ']': + g_pbdDualRuntime.pendingMeshResolution = nextPBDDualResolution(g_pbdDualRuntime.pendingMeshResolution); + logPBDDualControlState("pending-mesh+"); + return true; + case 'q': + case 'Q': + g_pbdDualRuntime.pendingSphereRadius = std::max( + PBDDualControlParam::sphereRadiusMin, + g_pbdDualRuntime.pendingSphereRadius - PBDDualControlParam::objectSizeStep + ); + logPBDDualControlState("pending-sphere-"); + return true; + case 'w': + case 'W': + g_pbdDualRuntime.pendingSphereRadius = std::min( + PBDDualControlParam::sphereRadiusMax, + g_pbdDualRuntime.pendingSphereRadius + PBDDualControlParam::objectSizeStep + ); + logPBDDualControlState("pending-sphere+"); + return true; + case 'a': + case 'A': + g_pbdDualRuntime.pendingCubeSize = std::max( + PBDDualControlParam::cubeSizeMin, + g_pbdDualRuntime.pendingCubeSize - PBDDualControlParam::objectSizeStep + ); + logPBDDualControlState("pending-cube-"); + return true; + case 's': + case 'S': + g_pbdDualRuntime.pendingCubeSize = std::min( + PBDDualControlParam::cubeSizeMax, + g_pbdDualRuntime.pendingCubeSize + PBDDualControlParam::objectSizeStep + ); + logPBDDualControlState("pending-cube+"); + return true; + case '9': { + const int nextIndex = std::max(0, g_pbdDualRuntime.speedPresetIndex - 1); + setPBDDualRuntimeTimestep(PBDFloorControlParam::speedPresets[nextIndex], nextIndex, false); + logPBDDualControlState("speed-"); + return true; + } + case '0': { + const int nextIndex = std::min( + static_cast(PBDFloorControlParam::speedPresets.size()) - 1, + g_pbdDualRuntime.speedPresetIndex + 1 + ); + setPBDDualRuntimeTimestep(PBDFloorControlParam::speedPresets[nextIndex], nextIndex, false); + logPBDDualControlState("speed+"); + return true; + } + case 'r': + case 'R': + resetPBDDualDemo(false); + return true; + case 't': + case 'T': + resetPBDDualDemo(true); + return true; + case 'p': + case 'P': + g_pbdDualRuntime.paused = !g_pbdDualRuntime.paused; + logPBDDualControlState(g_pbdDualRuntime.paused ? "pause" : "resume"); + return true; + default: + return false; + } +} + +static void keyboard(unsigned char key, int, int) { + if (g_pbdSolver == nullptr) return; + + bool handled = false; + switch (g_mode) { + case SimMode::PBDHangWind: + handled = handlePBDWindKey(key); + break; + case SimMode::PBDHang: + handled = handlePBDHangKey(key); + break; + case SimMode::PBDDrop: + handled = handlePBDDropKey(key); + break; + case SimMode::PBDDropFloor: + handled = handlePBDFloorKey(key); + break; + case SimMode::PBDDropFloorDual: + handled = handlePBDDualKey(key); + break; + default: + return; + } + + if (handled) { + glutPostRedisplay(); + } +} + static void mouse(const int button, const int state, const int x, const int y) { g_mouseClickX = x; g_mouseClickY = g_windowHeight - y - 1; @@ -364,19 +2969,19 @@ static void mouse(const int button, const int state, const int x, const int y) { // TODO: move to UserInteraction class: add renderer member variable // pick point - if (g_mouseLClickButton) { + if (g_mouseLClickButton && UI != nullptr) { UI->setModelview(g_ModelViewMatrix); UI->setProjection(g_ProjectionMatrix); UI->grabPoint(g_mouseClickX, g_mouseClickY); } - else UI->releasePoint(); + else if (UI != nullptr) UI->releasePoint(); } static void motion(const int x, const int y) { const float dx = float(x - g_mouseClickX); const float dy = float (-(g_windowHeight - y - 1 - g_mouseClickY)); - if (g_mouseLClickButton) { + if (g_mouseLClickButton && UI != nullptr) { //glm::vec3 ux(g_ModelViewMatrix * glm::vec4(1, 0, 0, 0)); //glm::vec3 uy(g_ModelViewMatrix * glm::vec4(0, 1, 0, 0)); glm::vec3 ux(0, 1, 0); @@ -389,28 +2994,470 @@ static void motion(const int x, const int y) { } // C L O T H /////////////////////////////////////////////////////////////////////// +static void drawFloor() { + Renderer renderer; + renderer.setProgram(g_phongShader); + renderer.setModelview(g_ModelViewMatrix); + renderer.setProjection(g_ProjectionMatrix); + g_phongShader->setUseFlagPattern(false); + g_phongShader->setAlbedo(g_floor_albedo); + g_phongShader->setAmbient(g_floor_ambient); + g_phongShader->setLight(g_light); + g_phongShader->setSpecularStrength(0.12f); + g_phongShader->setShininess(18.0f); + renderer.setProgramInput(g_floor_target); + renderer.setElementCount(6); + renderer.draw(); +} + +static void drawFloorShadows() { + if (!hasSceneFloor() || g_shadowShader == nullptr) return; + + const float shadowPlaneHeight = g_floor_collision_height + g_floor_render_offset + 1e-3f; + const glm::mat4 shadowModelView = g_ModelViewMatrix * floorShadowMatrix(shadowPlaneHeight, glm::normalize(g_light)); + + glEnable(GL_BLEND); + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); + glDepthMask(GL_FALSE); + + Renderer renderer; + renderer.setProgram(g_shadowShader); + renderer.setProjection(g_ProjectionMatrix); + renderer.setModelview(shadowModelView); + g_shadowShader->setShadowColor(g_shadow_color); + + renderer.setProgramInput(g_render_target); + renderer.setElementCount(g_clothMesh->ibuffLen()); + renderer.draw(); + + if (hasSphereColliderVisual() && g_sphere_target != nullptr) { + renderer.setProgramInput(g_sphere_target); + renderer.setElementCount(g_sphere_index_count); + renderer.draw(); + } + if (hasCubeColliderVisual() && g_cube_target != nullptr) { + renderer.setProgramInput(g_cube_target); + renderer.setElementCount(g_cube_index_count); + renderer.draw(); + } + + glDepthMask(GL_TRUE); + glDisable(GL_BLEND); +} + static void drawCloth() { Renderer renderer; renderer.setProgram(g_phongShader); renderer.setModelview(g_ModelViewMatrix); renderer.setProjection(g_ProjectionMatrix); + g_phongShader->setUseFlagPattern(g_mode == SimMode::PBDHangWind); g_phongShader->setAlbedo(g_albedo); g_phongShader->setAmbient(g_ambient); g_phongShader->setLight(g_light); + g_phongShader->setSpecularStrength(g_specular_strength); + g_phongShader->setShininess(g_shininess); renderer.setProgramInput(g_render_target); renderer.setElementCount(g_clothMesh->ibuffLen()); renderer.draw(); } -static void animateCloth(int value) { +static void drawBitmapText(float x, float y, const std::string& text) { + glRasterPos2f(x, y); + for (char c : text) { + glutBitmapCharacter(GLUT_BITMAP_8_BY_13, c); + } +} + +static void drawWindDirectionIndicator(const Eigen::Vector3f& normalizedDirection) { + glm::vec2 horizontalDirection(normalizedDirection.x(), normalizedDirection.y()); + const float horizontalMagnitude = glm::length(horizontalDirection); + + const glm::vec2 boxMin(18.0f, 18.0f); + const glm::vec2 boxMax(130.0f, 130.0f); + const glm::vec2 center = 0.5f * (boxMin + boxMax); + const float shaftLength = 34.0f; + const float headLength = 12.0f; + const float headWidth = 7.0f; + if (horizontalMagnitude > 1e-5f) { + horizontalDirection /= horizontalMagnitude; + } + const glm::vec2 tip = center + shaftLength * horizontalDirection; + const glm::vec2 headBase = tip - headLength * horizontalDirection; + const glm::vec2 normal(-horizontalDirection.y, horizontalDirection.x); + + const glm::vec2 verticalAnchor(boxMax.x + 26.0f, center.y); + const float zMagnitude = std::min(1.0f, std::abs(normalizedDirection.z())); + const float verticalShaftLength = 18.0f + 18.0f * zMagnitude; + const float verticalHeadLength = 10.0f; + const float verticalHeadWidth = 6.0f; + const float verticalSign = (normalizedDirection.z() >= 0.0f) ? 1.0f : -1.0f; + const glm::vec2 verticalDirection(0.0f, verticalSign); + const glm::vec2 verticalTip = verticalAnchor + verticalShaftLength * verticalDirection; + const glm::vec2 verticalHeadBase = verticalTip - verticalHeadLength * verticalDirection; + + glColor3f(0.72f, 0.78f, 0.88f); + glLineWidth(1.5f); + glBegin(GL_LINE_LOOP); + glVertex2f(boxMin.x, boxMin.y); + glVertex2f(boxMax.x, boxMin.y); + glVertex2f(boxMax.x, boxMax.y); + glVertex2f(boxMin.x, boxMax.y); + glEnd(); + + glBegin(GL_LINES); + glVertex2f(center.x - 5.0f, center.y); + glVertex2f(center.x + 5.0f, center.y); + glVertex2f(center.x, center.y - 5.0f); + glVertex2f(center.x, center.y + 5.0f); + glEnd(); + + glColor3f(0.95f, 0.97f, 1.0f); + if (horizontalMagnitude > 1e-5f) { + glLineWidth(2.5f); + glBegin(GL_LINES); + glVertex2f(center.x, center.y); + glVertex2f(tip.x, tip.y); + glEnd(); + + glBegin(GL_TRIANGLES); + glVertex2f(tip.x, tip.y); + glVertex2f(headBase.x + headWidth * normal.x, headBase.y + headWidth * normal.y); + glVertex2f(headBase.x - headWidth * normal.x, headBase.y - headWidth * normal.y); + glEnd(); + } + + glLineWidth(1.5f); + glBegin(GL_LINES); + glVertex2f(verticalAnchor.x, boxMin.y); + glVertex2f(verticalAnchor.x, boxMax.y); + glEnd(); + + if (zMagnitude > 1e-5f) { + glLineWidth(2.5f); + glBegin(GL_LINES); + glVertex2f(verticalAnchor.x, verticalAnchor.y); + glVertex2f(verticalTip.x, verticalTip.y); + glEnd(); + + glBegin(GL_TRIANGLES); + glVertex2f(verticalTip.x, verticalTip.y); + glVertex2f(verticalHeadBase.x + verticalHeadWidth, verticalHeadBase.y); + glVertex2f(verticalHeadBase.x - verticalHeadWidth, verticalHeadBase.y); + glEnd(); + } + + glLineWidth(1.0f); + drawBitmapText(boxMin.x + 30.0f, boxMin.y - 14.0f, "Wind"); + drawBitmapText(verticalAnchor.x - 7.0f, boxMin.y - 14.0f, "Z"); +} + +static void drawPBDHangOverlay() { + if (g_mode != SimMode::PBDHang || g_pbdSolver == nullptr) return; + + const GLboolean depthEnabled = glIsEnabled(GL_DEPTH_TEST); + glDisable(GL_DEPTH_TEST); - // solve two time-steps - g_solver->solve(g_iter); - g_solver->solve(g_iter); + glMatrixMode(GL_PROJECTION); + glPushMatrix(); + glLoadIdentity(); + gluOrtho2D(0.0, static_cast(g_windowWidth), 0.0, static_cast(g_windowHeight)); + + glMatrixMode(GL_MODELVIEW); + glPushMatrix(); + glLoadIdentity(); + glColor3f(1.0f, 1.0f, 1.0f); + + std::ostringstream valueLine; + valueLine << std::fixed << std::setprecision(3) + << "stretch=" << g_pbdSolver->getStructuralStiffness() + << " shear=" << g_pbdSolver->getShearStiffness() + << " bend=" << g_pbdSolver->getBendStiffness() + << " damping=" << g_pbdSolver->getDampingFactor(); + + std::ostringstream iterLine; + iterLine << "iters=" << g_pbdSolver->getSolverIterations() + << " state=" << (g_pbdHangRuntime.paused ? "paused" : "running"); + + drawBitmapText(16.0f, g_windowHeight - 22.0f, "PBD hang tuning: 1/2 stretch 3/4 shear 5/6 bend 7/8 damping R reset cloth T reset params P pause"); + drawBitmapText(16.0f, g_windowHeight - 40.0f, valueLine.str()); + drawBitmapText(16.0f, g_windowHeight - 58.0f, iterLine.str()); + + glPopMatrix(); + glMatrixMode(GL_PROJECTION); + glPopMatrix(); + glMatrixMode(GL_MODELVIEW); + + if (depthEnabled) { + glEnable(GL_DEPTH_TEST); + } +} + +static void drawPBDWindOverlay() { + if (g_mode != SimMode::PBDHangWind || g_pbdSolver == nullptr) return; + + const GLboolean depthEnabled = glIsEnabled(GL_DEPTH_TEST); + glDisable(GL_DEPTH_TEST); + + glMatrixMode(GL_PROJECTION); + glPushMatrix(); + glLoadIdentity(); + gluOrtho2D(0.0, static_cast(g_windowWidth), 0.0, static_cast(g_windowHeight)); + + glMatrixMode(GL_MODELVIEW); + glPushMatrix(); + glLoadIdentity(); + glColor3f(1.0f, 1.0f, 1.0f); + + const Eigen::Vector3f& rawDirection = g_pbdWindRuntime.currentDirectionInput; + const Eigen::Vector3f& normalizedDirection = g_pbdWindRuntime.appliedWindDirection; + std::ostringstream line1; + line1 << std::fixed << std::setprecision(3) + << "PBD Wind Demo" + << " speed=" << g_pbdWindRuntime.currentWindSpeed + << " range=[0, 15]" + << " state=" << (g_pbdWindRuntime.paused ? "paused" : "running"); + + std::ostringstream line2; + line2 << std::fixed << std::setprecision(3) + << "dir=(" << rawDirection.x() << ", " << rawDirection.y() << ", " << rawDirection.z() << ")" + << " normalized=(" << normalizedDirection.x() << ", " << normalizedDirection.y() << ", " << normalizedDirection.z() << ")"; + + drawBitmapText(16.0f, g_windowHeight - 22.0f, "Controls: 1/2 speed 3/4 dir-x 5/6 dir-y 7/8 dir-z R reset cloth T reset wind P pause"); + drawBitmapText(16.0f, g_windowHeight - 40.0f, line1.str()); + drawBitmapText(16.0f, g_windowHeight - 58.0f, line2.str()); + drawWindDirectionIndicator(normalizedDirection); + if (std::abs(normalizedDirection.z()) > PBDWindControlParam::verticalWarningThreshold) { + drawBitmapText(16.0f, g_windowHeight - 76.0f, "Warning: wind has a vertical component. For flag-style motion, use z near 0."); + } + + glPopMatrix(); + glMatrixMode(GL_PROJECTION); + glPopMatrix(); + glMatrixMode(GL_MODELVIEW); + + if (depthEnabled) { + glEnable(GL_DEPTH_TEST); + } +} + +static void drawPBDDropOverlay() { + if (g_mode != SimMode::PBDDrop || g_pbdSolver == nullptr) return; + + const GLboolean depthEnabled = glIsEnabled(GL_DEPTH_TEST); + glDisable(GL_DEPTH_TEST); + + glMatrixMode(GL_PROJECTION); + glPushMatrix(); + glLoadIdentity(); + gluOrtho2D(0.0, static_cast(g_windowWidth), 0.0, static_cast(g_windowHeight)); + + glMatrixMode(GL_MODELVIEW); + glPushMatrix(); + glLoadIdentity(); + glColor3f(1.0f, 1.0f, 1.0f); + + std::ostringstream valueLine; + valueLine << std::fixed << std::setprecision(3) + << "radius=" << g_pbdDropRuntime.sphereRadius + << " iters=" << g_pbdSolver->getSolverIterations() + << " stretch=" << g_pbdSolver->getStructuralStiffness() + << " shear=" << g_pbdSolver->getShearStiffness(); + + std::ostringstream valueLine2; + valueLine2 << std::fixed << std::setprecision(3) + << "bend=" << g_pbdSolver->getBendStiffness() + << " damping=" << g_pbdSolver->getDampingFactor() + << " mesh=" << g_pbdDropRuntime.currentMeshResolution << "x" << g_pbdDropRuntime.currentMeshResolution + << " state=" << (g_pbdDropRuntime.paused ? "paused" : "running"); + + drawBitmapText(16.0f, g_windowHeight - 22.0f, "PBD drop tuning: 1/2 stretch 3/4 shear 5/6 bend 7/8 damping [/ ] pending mesh R apply reset T reset params P pause"); + drawBitmapText(16.0f, g_windowHeight - 40.0f, valueLine.str()); + drawBitmapText(16.0f, g_windowHeight - 58.0f, valueLine2.str()); + if (g_pbdDropRuntime.pendingMeshResolution != g_pbdDropRuntime.currentMeshResolution) { + std::ostringstream pendingLine; + pendingLine << "Pending mesh: " + << g_pbdDropRuntime.pendingMeshResolution << "x" << g_pbdDropRuntime.pendingMeshResolution + << ", press R to apply"; + drawBitmapText(16.0f, g_windowHeight - 76.0f, pendingLine.str()); + } - // fix points - CgSatisfyVisitor visitor; - visitor.satisfy(*g_cgRootNode); + glPopMatrix(); + glMatrixMode(GL_PROJECTION); + glPopMatrix(); + glMatrixMode(GL_MODELVIEW); + + if (depthEnabled) { + glEnable(GL_DEPTH_TEST); + } +} + +static void drawPBDFloorOverlay() { + if (g_mode != SimMode::PBDDropFloor || g_pbdSolver == nullptr) return; + + const GLboolean depthEnabled = glIsEnabled(GL_DEPTH_TEST); + glDisable(GL_DEPTH_TEST); + + glMatrixMode(GL_PROJECTION); + glPushMatrix(); + glLoadIdentity(); + gluOrtho2D(0.0, static_cast(g_windowWidth), 0.0, static_cast(g_windowHeight)); + + glMatrixMode(GL_MODELVIEW); + glPushMatrix(); + glLoadIdentity(); + glColor3f(1.0f, 1.0f, 1.0f); + + std::ostringstream line1; + line1 << std::fixed << std::setprecision(4) + << "PBD Drop-Floor Self-Collision Demo" + << " iters=" << g_pbdFloorRuntime.solverIterations + << " dt=" << g_pbdFloorRuntime.currentTimestep; + if (!g_pbdFloorRuntime.customTimestep) { + line1 << " (" << PBDFloorControlParam::speedPresetLabels[g_pbdFloorRuntime.speedPresetIndex] << ")"; + } + + std::ostringstream line2; + line2 << std::fixed << std::setprecision(3) + << "self-k=" << g_pbdSolver->getSelfCollisionStiffness() + << " max-contacts=" << g_pbdSolver->getMaxSelfCollisionContactsPerVertex() + << " self-thickness=" << g_pbdSolver->getSelfCollisionThickness() + << " floor-friction=" << g_pbdSolver->getPlaneFriction(); + + std::ostringstream line3; + line3 << std::fixed << std::setprecision(3) + << "bend=" << g_pbdSolver->getBendStiffness() + << " damping=" << g_pbdSolver->getDampingFactor() + << " mesh=" << g_pbdFloorRuntime.currentMeshResolution << "x" << g_pbdFloorRuntime.currentMeshResolution + << " debug=" << (g_enableDebugDiagnostics ? "on" : "off") + << " state=" << (g_pbdFloorRuntime.paused ? "paused" : "running"); + + drawBitmapText(16.0f, g_windowHeight - 22.0f, line1.str()); + drawBitmapText(16.0f, g_windowHeight - 40.0f, "Controls: 1/2 self-k 3/4 max contacts 5/6 thickness 7/8 friction Q/W bend A/S damping [/ ] mesh 9/0 speed R reset T defaults P pause D debug"); + drawBitmapText(16.0f, g_windowHeight - 58.0f, line2.str()); + drawBitmapText(16.0f, g_windowHeight - 76.0f, line3.str()); + if (g_pbdFloorRuntime.pendingMeshResolution != g_pbdFloorRuntime.currentMeshResolution) { + std::ostringstream pendingLine; + pendingLine << "Pending mesh: " << g_pbdFloorRuntime.pendingMeshResolution << "x" << g_pbdFloorRuntime.pendingMeshResolution << ", press R to apply"; + drawBitmapText(16.0f, g_windowHeight - 94.0f, pendingLine.str()); + } + if (g_pbdFloorRuntime.speedPresetIndex > 0) { + drawBitmapText(16.0f, g_windowHeight - 112.0f, "Warning: faster timestep may increase penetration, jitter, or missed self-collision"); + } + + glPopMatrix(); + glMatrixMode(GL_PROJECTION); + glPopMatrix(); + glMatrixMode(GL_MODELVIEW); + + if (depthEnabled) { + glEnable(GL_DEPTH_TEST); + } +} + +static void drawPBDDualOverlay() { + if (g_mode != SimMode::PBDDropFloorDual || g_pbdSolver == nullptr) return; + + const GLboolean depthEnabled = glIsEnabled(GL_DEPTH_TEST); + glDisable(GL_DEPTH_TEST); + + glMatrixMode(GL_PROJECTION); + glPushMatrix(); + glLoadIdentity(); + gluOrtho2D(0.0, static_cast(g_windowWidth), 0.0, static_cast(g_windowHeight)); + + glMatrixMode(GL_MODELVIEW); + glPushMatrix(); + glLoadIdentity(); + glColor3f(1.0f, 1.0f, 1.0f); + + std::ostringstream line1; + line1 << std::fixed << std::setprecision(4) + << "PBD Dual-Obstacle Demo" + << " iters=" << g_pbdDualRuntime.solverIterations + << " dt=" << g_pbdDualRuntime.currentTimestep + << " stretch=" << g_pbdSolver->getStructuralStiffness() + << " shear=" << g_pbdSolver->getShearStiffness(); + if (!g_pbdDualRuntime.customTimestep) { + line1 << " (" << PBDFloorControlParam::speedPresetLabels[g_pbdDualRuntime.speedPresetIndex] << ")"; + } + + std::ostringstream line2; + line2 << std::fixed << std::setprecision(3) + << "bend=" << g_pbdSolver->getBendStiffness() + << " mesh=" << g_pbdDualRuntime.currentMeshResolution << "x" << g_pbdDualRuntime.currentMeshResolution + << " sphere=" << g_pbdDualRuntime.currentSphereRadius + << " cube=" << g_pbdDualRuntime.currentCubeSize + << " state=" << (g_pbdDualRuntime.paused ? "paused" : "running"); + + drawBitmapText(16.0f, g_windowHeight - 22.0f, "Controls: 1/2 stretch 3/4 shear 5/6 bend [/ ] mesh Q/W sphere A/S cube 9/0 speed R reset cloth T reset params P pause"); + drawBitmapText(16.0f, g_windowHeight - 40.0f, line1.str()); + drawBitmapText(16.0f, g_windowHeight - 58.0f, line2.str()); + + const bool hasPendingMesh = g_pbdDualRuntime.pendingMeshResolution != g_pbdDualRuntime.currentMeshResolution; + const bool hasPendingSphere = std::abs(g_pbdDualRuntime.pendingSphereRadius - g_pbdDualRuntime.currentSphereRadius) > 1e-5f; + const bool hasPendingCube = std::abs(g_pbdDualRuntime.pendingCubeSize - g_pbdDualRuntime.currentCubeSize) > 1e-5f; + float y = g_windowHeight - 76.0f; + if (hasPendingMesh) { + std::ostringstream pendingLine; + pendingLine << "Pending mesh: " << g_pbdDualRuntime.pendingMeshResolution << "x" << g_pbdDualRuntime.pendingMeshResolution; + drawBitmapText(16.0f, y, pendingLine.str()); + y -= 18.0f; + } + if (hasPendingSphere) { + std::ostringstream pendingLine; + pendingLine << std::fixed << std::setprecision(3) << "Pending sphere radius: " << g_pbdDualRuntime.pendingSphereRadius; + drawBitmapText(16.0f, y, pendingLine.str()); + y -= 18.0f; + } + if (hasPendingCube) { + std::ostringstream pendingLine; + pendingLine << std::fixed << std::setprecision(3) << "Pending cube size: " << g_pbdDualRuntime.pendingCubeSize; + drawBitmapText(16.0f, y, pendingLine.str()); + y -= 18.0f; + } + if (hasPendingMesh || hasPendingSphere || hasPendingCube) { + drawBitmapText(16.0f, y, "Pending changes: press R to apply"); + y -= 18.0f; + } + if (g_pbdDualRuntime.speedPresetIndex > 0) { + drawBitmapText(16.0f, y, "Warning: faster timestep may increase penetration, jitter, or missed self-collision"); + } + + glPopMatrix(); + glMatrixMode(GL_PROJECTION); + glPopMatrix(); + glMatrixMode(GL_MODELVIEW); + + if (depthEnabled) { + glEnable(GL_DEPTH_TEST); + } +} + +static void animateCloth(int value) { + if (isPBDMode()) { + const bool paused = (g_mode == SimMode::PBDHang && g_pbdHangRuntime.paused) + || (g_mode == SimMode::PBDHangWind && g_pbdWindRuntime.paused) + || (g_mode == SimMode::PBDDrop && g_pbdDropRuntime.paused) + || (g_mode == SimMode::PBDDropFloor && g_pbdFloorRuntime.paused) + || (g_mode == SimMode::PBDDropFloorDual && g_pbdDualRuntime.paused); + if (!paused) { + const unsigned int iterationCount = (isFloorDemo() || g_mode == SimMode::PBDHang || g_mode == SimMode::PBDDrop) + ? g_pbdSolver->getSolverIterations() + : static_cast(PBDSystemParam::n_iter); + g_pbdSolver->solve(iterationCount); + ++g_pbdFrameCounter; + if (g_enableDebugDiagnostics) { + logPBDSelfCollisionDiagnostics(); + } + } + } + else { + g_solver->solve(g_iter); + g_solver->solve(g_iter); + + CgSatisfyVisitor visitor; + visitor.satisfy(*g_cgRootNode); + } // update normals g_clothMesh->request_face_normals(); @@ -427,6 +3474,28 @@ static void animateCloth(int value) { glutTimerFunc(g_animation_timer, animateCloth, 0); } +static void logPBDSelfCollisionDiagnostics() { + if (g_pbdSolver == nullptr) return; + if (PBDDebugParam::debugPrintPeriod == 0u) return; + if ((g_pbdFrameCounter % PBDDebugParam::debugPrintPeriod) != 0u) return; + + const SelfCollisionDebugStats& stats = g_pbdSolver->getSelfCollisionDebugStats(); + const char* modeLabel = + (g_mode == SimMode::PBDDropFloorDual) + ? "pbd drop-floor-dual" + : (isFloorDemo() + ? "pbd drop-floor" + : (g_mode == SimMode::PBDDrop ? "pbd drop" : (g_mode == SimMode::PBDHangWind ? "pbd hang_wind" : "pbd hang"))); + std::cout + << "[" << modeLabel << " self-collision] frame=" << g_pbdFrameCounter + << " generated=" << stats.generatedContacts + << " initial-violations=" << stats.initiallyViolatedContacts + << " remaining-violations=" << stats.remainingViolatedContacts + << " max-pen-before=" << stats.maxInitialPenetration + << " max-pen-after=" << stats.maxRemainingPenetration + << std::endl; +} + // S C E N E U P D A T E /////////////////////////////////////////////////////////// static void updateProjection() { g_ProjectionMatrix = glm::perspective(PI / 4.0f, @@ -453,10 +3522,16 @@ static void cleanUp() { // delete render target delete g_render_target; + delete g_floor_target; + delete g_sphere_target; + delete g_cube_target; // delete mass-spring system delete g_system; delete g_solver; + delete g_pbdSystem; + delete g_pbdSolver; + delete g_shadowShader; // delete constraint graph // TODO diff --git a/ClothApp/shaders/basic.vshader b/ClothApp/shaders/basic.vshader index b6128bf..a28073c 100644 --- a/ClothApp/shaders/basic.vshader +++ b/ClothApp/shaders/basic.vshader @@ -1,19 +1,24 @@ #version 430 uniform mat4 uModelViewMatrix; -// uniform mat4 uNormalMatrix; uniform mat4 uProjectionMatrix; +uniform vec3 uLight; layout(location=0) in vec3 aPosition; layout(location=1) in vec3 aNormal; layout(location=2) in vec2 aTexCoord; out vec3 vNormal; +out vec3 vLightDir; +out vec3 vViewPosition; out vec2 vTexCoord; void main(){ - vNormal = aNormal; + mat3 normalMatrix = transpose(inverse(mat3(uModelViewMatrix))); + vec4 viewPosition = uModelViewMatrix * vec4(aPosition, 1.0); + vNormal = normalize(normalMatrix * aNormal); + vLightDir = normalize(mat3(uModelViewMatrix) * (-uLight)); + vViewPosition = viewPosition.xyz; vTexCoord = aTexCoord; - vec4 position = vec4(aPosition, 1.0); - gl_Position = uProjectionMatrix * uModelViewMatrix * position; + gl_Position = uProjectionMatrix * viewPosition; } \ No newline at end of file diff --git a/ClothApp/shaders/phong.fshader b/ClothApp/shaders/phong.fshader index c473e87..d0ca3b6 100644 --- a/ClothApp/shaders/phong.fshader +++ b/ClothApp/shaders/phong.fshader @@ -1,22 +1,80 @@ #version 430 in vec3 vNormal; +in vec3 vLightDir; +in vec3 vViewPosition; +in vec2 vTexCoord; out vec4 fragColor; uniform vec3 uAlbedo; uniform vec3 uAmbient; -uniform vec3 uLight; +uniform float uSpecularStrength; +uniform float uShininess; +uniform bool uUseFlagPattern; + +float starMask(vec2 uv, vec2 center, float radiusOuter, float radiusInner) { + vec2 delta = uv - center; + float angle = atan(delta.y, delta.x); + float distanceValue = length(delta); + float sector = floor((angle + 3.14159265) * 5.0 / 3.14159265); + float localAngle = angle + sector * 3.14159265 / 5.0; + float targetRadius = mod(sector, 2.0) < 0.5 ? radiusOuter : radiusInner; + return step(distanceValue, targetRadius * max(0.35, abs(cos(localAngle)))); +} + +vec3 americanFlagColor(vec2 uv) { + const vec3 red = vec3(0.698, 0.132, 0.203); + const vec3 white = vec3(0.96, 0.96, 0.94); + const vec3 blue = vec3(0.235, 0.233, 0.431); + + vec2 clampedUv = clamp(uv, 0.0, 1.0); + float stripeIndex = floor((1.0 - clampedUv.y) * 13.0); + vec3 color = mod(stripeIndex, 2.0) < 0.5 ? red : white; + + bool inUnion = clampedUv.x <= 0.4 && clampedUv.y >= (1.0 - 7.0 / 13.0); + if (inUnion) { + color = blue; + vec2 unionUv = vec2(clampedUv.x / 0.4, (clampedUv.y - (1.0 - 7.0 / 13.0)) / (7.0 / 13.0)); + float stars = 0.0; + for (int row = 0; row < 9; ++row) { + float starsInRow = (row % 2 == 0) ? 6.0 : 5.0; + float xOffset = (row % 2 == 0) ? 1.0 / 12.0 : 1.0 / 6.0; + for (int col = 0; col < int(starsInRow); ++col) { + vec2 center = vec2(xOffset + float(col) / starsInRow, 1.0 - (float(row) + 0.5) / 9.0); + stars = max(stars, starMask(unionUv, center, 0.045, 0.018)); + } + } + color = mix(color, white, stars); + } + + return color; +} + void main(){ vec3 normal = normalize(vNormal); vec3 albedo = uAlbedo; + vec2 texCoord = vec2(clamp(vTexCoord.x, 0.0, 1.0), 1.0 - clamp(vTexCoord.y, 0.0, 1.0)); + if(!gl_FrontFacing) { - normal = -normal; - albedo = 1.0 - albedo; + normal = -normal; + if (!uUseFlagPattern) { + albedo = 1.0 - albedo; + } + } + + if (uUseFlagPattern) { + albedo = americanFlagColor(texCoord); } - vec3 toLight = normalize(-uLight); + vec3 toLight = normalize(vLightDir); + vec3 toView = normalize(-vViewPosition); float diffuse = max(0, dot(toLight, normal)); + vec3 halfVector = normalize(toLight + toView); + float specular = 0.0; + if (diffuse > 0.0) { + specular = uSpecularStrength * pow(max(0.0, dot(normal, halfVector)), uShininess); + } - vec3 color = diffuse * albedo + uAmbient * albedo; + vec3 color = diffuse * albedo + uAmbient * albedo + specular * vec3(1.0); fragColor = vec4(color, 1.0); } \ No newline at end of file diff --git a/ClothApp/shaders/shadow.fshader b/ClothApp/shaders/shadow.fshader new file mode 100644 index 0000000..d99781d --- /dev/null +++ b/ClothApp/shaders/shadow.fshader @@ -0,0 +1,9 @@ +#version 430 + +out vec4 fragColor; + +uniform vec4 uShadowColor; + +void main() { + fragColor = uShadowColor; +} \ No newline at end of file diff --git a/Demos/Drop_floor_trim.gif b/Demos/Drop_floor_trim.gif new file mode 100644 index 0000000..6f3a87f Binary files /dev/null and b/Demos/Drop_floor_trim.gif differ diff --git a/Demos/Drop_sphere_trim.gif b/Demos/Drop_sphere_trim.gif new file mode 100644 index 0000000..f64ccb8 Binary files /dev/null and b/Demos/Drop_sphere_trim.gif differ diff --git a/Demos/Dual_trim.gif b/Demos/Dual_trim.gif new file mode 100644 index 0000000..e61c10a Binary files /dev/null and b/Demos/Dual_trim.gif differ diff --git a/Demos/Flag_trim.gif b/Demos/Flag_trim.gif new file mode 100644 index 0000000..26cea2c Binary files /dev/null and b/Demos/Flag_trim.gif differ diff --git a/Demos/Hang_trim.gif b/Demos/Hang_trim.gif new file mode 100644 index 0000000..b7d98dd Binary files /dev/null and b/Demos/Hang_trim.gif differ diff --git a/readme.md b/readme.md index e19933b..6cbfc09 100644 --- a/readme.md +++ b/readme.md @@ -1,43 +1,514 @@ -### Fast Mass-Spring System Simulator +# PBD Cloth Simulation -A C++ implementation of *Fast Simulation of Mass-Spring Systems* [1], rendered with OpenGL. -The dynamic inverse procedure described in [2] was implemented to constrain spring deformations -and prevent the "super-elastic" effect when using large time-steps. +**Name:** Vincent Chen +**SBU ID:** 115598737 +**Course:** Stony Brook University — CSE 328: Computer Graphics -### Dependencies +This project was originally forked from [sam007961/FastMassSpring](https://github.com/sam007961/FastMassSpring). + +--- + +## Project Summary + +This project is a C++ OpenGL cloth-simulation application built on top of the original **Fast Simulation of Mass-Spring Systems** codebase [1] and extended with a separate **Position-Based Dynamics (PBD)** solver based on Müller et al. [3]. It keeps the original mass-spring implementation for baseline comparison while adding a broader set of cloth behaviors and interaction-focused PBD demos. + +The current program includes both the original mass-spring hanging and sphere-drop scenes and a set of PBD demos for hanging cloth, sphere drop, floor drop with self-collision, dual-obstacle draping over a sphere and cube, and a flag-style wind scene. The PBD path supports fixed-point, stretch, shear, and bending constraints; sphere, floor, and analytic box collision; vertex-triangle self-collision with optional diagnostics; and interactive runtime controls for parameters such as stiffness, damping, timestep, mesh resolution, obstacle size, and wind direction/speed. -* **OpenGL, freeGLUT, GLEW, GLM** for rendering. -* **OpenMesh** for computing normals. -* **Eigen** for sparse matrix algebra. +The project uses the shared rendering, mesh, shader, and picking infrastructure from the original application, so both solvers can be launched from the same executable and compared through a common visual interface. -### Building +--- -You need to install OpenGL, GLEW and GLUT on your system to build. The rest of the dependencies -will be automatically fetched by cmake during configuration. +## Key Features + +- Original **fast mass-spring** solver preserved for side-by-side comparison +- Separate **Position-Based Dynamics** cloth solver integrated into the same executable +- Fixed-point, stretch, shear, and dihedral bending constraints +- Iteration-corrected stiffness handling for more consistent tuning across solver iteration counts +- Sphere, floor, and analytic box collision handling +- Vertex-triangle self-collision with optional debug diagnostics +- Dedicated PBD hang, sphere-drop, floor-drop, dual-obstacle, and wind demos +- Runtime controls for stiffness, damping, timestep, mesh resolution, obstacle size, and wind settings +- Shared rendering, shader, mesh, and picking infrastructure across both solver modes + +--- + +## Build Instructions + +### Dependencies -To build, run the following commands from the root directory of the project: -``` bash +- **OpenGL, freeGLUT, GLEW, GLM** for rendering +- **OpenMesh** for computing normals +- **Eigen** for vector/matrix operations + +Install the required OpenGL/GLUT/GLEW dependencies, then build with CMake: + +```bash mkdir build cd build cmake .. cmake --build . ``` -On Windows, you will likely need to specify the directories containing GLUT and GLEW in CMAKE_PREFIX_PATH so that cmake can find them. +Notes: + +- Run the executable from the `build` directory so shader paths resolve correctly. +- Shader files are copied from `ClothApp/shaders` into `build/shaders` during each build. +- The project fetches **OpenMesh**, **Eigen**, and **GLM** through CMake. +- Eigen is used as a header-only dependency in the current CMake setup. + +On Windows, you may need to specify the directories containing GLUT and GLEW in `CMAKE_PREFIX_PATH`: + +```bash +cmake .. -DCMAKE_PREFIX_PATH:PATH=/path/to/libs +``` + +You may also need to copy the required DLLs into the build directory if they are not available globally. + +--- + +## Running + +Run the executable from the `build` directory: + +```bash +cd build +./fast-mass-spring +``` + +If no arguments are provided, the program defaults to the original mass-spring hanging cloth demo. + +### Command Syntax + +Long-form mode selection: + +```bash +./fast-mass-spring [mass-spring|ms] [hang|drop] [--self-thickness value] [--debug] +./fast-mass-spring pbd [hang|hang-wind|drop|drop-floor|drop-floor-dual] [--self-thickness value] [--debug] [--wind-speed value] [--wind-dir x y z] [--iters value] [--radius value] [--dt value] +``` + +Short-form mode selection: -``` bash -cmake .. -DCMAKE_PERFIX_PATH:PATH=/path/to/libs +```bash +./fast-mass-spring [ms-hang|ms-drop|pbd-hang|pbd-hang-wind|pbd-drop|pbd-drop-floor|pbd-drop-floor-dual] [--self-thickness value] [--debug] [--wind-speed value] [--wind-dir x y z] [--iters value] [--radius value] [--dt value] ``` -You will also need to copy the DLLs to the build directory if they are not available globally. +### Demo Modes + +| Demo | Command | Purpose | +|---|---|---| +| Mass-spring hang | `./fast-mass-spring ms hang` | Original baseline hanging cloth | +| Mass-spring drop | `./fast-mass-spring ms drop` | Original baseline sphere-drop demo | +| PBD hang | `./fast-mass-spring pbd hang --iters 20` | Interactive constraint tuning demo | +| PBD drop | `./fast-mass-spring pbd drop` | PBD cloth dropping onto a sphere | +| PBD drop-floor | `./fast-mass-spring pbd drop-floor --iters 28 --dt 0.003` | Floor contact and self-collision stress test | +| PBD drop-floor-dual | `./fast-mass-spring pbd drop-floor-dual --iters 28` | Cloth interaction with sphere and box obstacles | +| PBD hang-wind | `./fast-mass-spring pbd hang-wind --wind-speed 5 --wind-dir 1 0 0` | Flag-style wind demo | + +### Example Commands + +```bash +./fast-mass-spring mass-spring hang +./fast-mass-spring ms drop +./fast-mass-spring pbd hang --iters 20 +./fast-mass-spring pbd drop +./fast-mass-spring pbd drop-floor +./fast-mass-spring pbd drop-floor --iters 28 --dt 0.003 +./fast-mass-spring pbd drop-floor --debug +./fast-mass-spring pbd drop-floor --self-thickness 0.02 +./fast-mass-spring pbd drop-floor-dual +./fast-mass-spring pbd drop-floor-dual --iters 28 +./fast-mass-spring pbd hang-wind --wind-speed 5 --wind-dir 1 0 0 +./fast-mass-spring pbd-hang-wind --wind-speed 8 +``` + +### Demo Gallery + +#### Hang Demo + +![Hang demo](Demos/Hang_trim.gif) + +#### Wind Flag Demo + +![Wind flag demo](Demos/Flag_trim.gif) + +#### Sphere Drop Demo + +![Sphere drop demo](Demos/Drop_sphere_trim.gif) + +#### Floor Drop Demo + +![Floor drop demo](Demos/Drop_floor_trim.gif) + +#### Dual-Obstacle Demo + +![Dual-obstacle demo](Demos/Dual_trim.gif) + +--- + +## Controls and Runtime Interface + +### PBD Hang Controls + +The `pbd hang` demo includes an interactive tuning interface for studying how different constraints affect cloth behavior. + +Startup option: + +- `--iters value`: sets the number of PBD projection passes per timestep for the hang demo; accepts an integer in `[1, 80]` + +Keyboard controls: + +| Key | Action | +|---|---| +| `1 / 2` | Decrease / increase stretch stiffness | +| `3 / 4` | Decrease / increase shear stiffness | +| `5 / 6` | Decrease / increase bend stiffness | +| `7 / 8` | Decrease / increase damping factor | +| `R` | Reset cloth to the initial hanging pose while keeping current tuning values | +| `T` | Reset hang-demo tuning values to defaults | +| `P` | Pause / resume simulation | + +Parameter meaning: + +- **Stretch stiffness** controls structural distance constraints. +- **Shear stiffness** controls diagonal/shear distance constraints. +- **Bend stiffness** controls dihedral bending constraints. +- **Damping** controls velocity energy decay. +- **Iterations** control the number of PBD projection passes per timestep. + +The current hang-demo values are displayed in the on-screen overlay while the demo is running. + +### PBD Drop Controls + +The `pbd drop` demo includes an interactive tuning interface for sphere collision, material tuning, damping, and mesh-resolution comparison. + +Startup options: + +- `--radius value`: sets the sphere collider radius for the drop demo; accepts a float in `[0.1, 1.5]` +- `--iters value`: sets the number of PBD projection passes per timestep for the drop demo; accepts an integer in `[1, 80]` + +Keyboard controls: + +| Key | Action | +|---|---| +| `1 / 2` | Decrease / increase stretch stiffness | +| `3 / 4` | Decrease / increase shear stiffness | +| `5 / 6` | Decrease / increase bend stiffness | +| `7 / 8` | Decrease / increase damping factor | +| `[` / `]` | Decrease / increase the pending mesh resolution by `2` | +| `R` | Reset cloth and apply the pending mesh resolution | +| `T` | Reset drop-demo tuning values to defaults and reset pending mesh resolution to `33` | +| `P` | Pause / resume simulation | + +Parameter meaning: + +- **Sphere radius** changes the collision equation $C(p) = \lVert p - c \rVert - r \ge 0$. +- **Stretch stiffness** controls structural distance constraints. +- **Shear stiffness** controls diagonal/shear distance constraints. +- **Bend stiffness** controls dihedral bending constraints. +- **Damping** controls velocity energy decay. +- **Iterations** control the number of PBD projection passes per timestep. +- **Mesh resolution** changes particle and constraint count and is applied only when the cloth system is rebuilt on reset. + +The current and pending mesh resolutions are shown in the on-screen overlay. If the pending mesh differs from the current mesh, the overlay prints `Pending mesh: NxN, press R to apply`. + +### PBD Drop-Floor Controls + +The `pbd drop-floor` demo includes an interactive tuning interface for self-collision robustness, floor response, timestep speed, and mesh-resolution comparison. + +Startup options: + +- `--iters value`: sets the number of PBD projection passes per timestep for the drop-floor demo; accepts an integer in `[1, 80]` +- `--dt value`: sets the simulation timestep for the drop-floor demo; accepts a float in `[0.001, 0.01]` + +Keyboard controls: + +| Key | Action | +|---|---| +| `1 / 2` | Decrease / increase self-collision stiffness | +| `3 / 4` | Decrease / increase max self-collision contacts per vertex | +| `5 / 6` | Decrease / increase self-collision thickness | +| `7 / 8` | Decrease / increase floor friction | +| `Q / W` | Decrease / increase bend stiffness | +| `A / S` | Decrease / increase damping factor | +| `[` / `]` | Decrease / increase the pending mesh resolution | +| `9 / 0` | Step to a slower / faster timestep preset | +| `R` | Reset cloth and apply the pending mesh resolution | +| `T` | Reset drop-floor tuning values to stable defaults | +| `P` | Pause / resume simulation | +| `D` | Toggle debug diagnostics | + +Parameter meaning: + +- **Self-collision stiffness** controls how strongly vertex self-collision constraints are projected apart. +- **Max self-collision contacts** limits how many self-collision contacts are processed per vertex each solver step. +- **Self-collision thickness** is the minimum separation band enforced between cloth layers. +- **Floor friction** damps tangential sliding after plane contact. +- **Bend stiffness** controls dihedral bending constraints. +- **Damping** controls velocity energy decay. +- **Iterations** control the number of PBD projection passes per timestep. +- **Timestep** changes simulation speed and contact robustness; faster presets are less stable and are labeled in the overlay. +- **Mesh resolution** changes particle and constraint count and is applied only when the cloth system is rebuilt on reset. + +The overlay shows the current self-collision, floor, timestep, and mesh settings. If the pending mesh differs from the current mesh, the overlay prints `Pending mesh: NxN, press R to apply`. + +### PBD Drop-Floor-Dual Controls + +The `pbd drop-floor-dual` demo includes a simple interactive interface for comparing how the cloth drapes over a smooth sphere versus a sharp cube. -### Demonstration +Startup option: -![curtain_hang](https://user-images.githubusercontent.com/24758349/79005907-97ad1100-7b60-11ea-9e27-90375461beaf.gif) +- `--iters value`: sets the number of PBD projection passes per timestep for the dual-obstacle demo; accepts an integer in `[1, 80]` + +Keyboard controls: + +| Key | Action | +|---|---| +| `1 / 2` | Decrease / increase stretch stiffness | +| `3 / 4` | Decrease / increase shear stiffness | +| `5 / 6` | Decrease / increase bend stiffness | +| `[` / `]` | Decrease / increase the pending mesh resolution | +| `Q / W` | Decrease / increase the pending sphere radius | +| `A / S` | Decrease / increase the pending cube size | +| `9 / 0` | Step to a slower / faster timestep preset | +| `R` | Reset cloth and apply the pending mesh and obstacle sizes | +| `T` | Reset material values and pending mesh/object sizes to defaults | +| `P` | Pause / resume simulation | + +Parameter meaning: + +- **Stretch stiffness** controls structural distance constraints. +- **Shear stiffness** controls diagonal/shear distance constraints. +- **Bend stiffness** controls dihedral bending constraints. +- **Mesh resolution** changes particle, constraint, and contact count and is applied only when the cloth system is rebuilt on reset. +- **Sphere radius** changes the smooth obstacle collision size. +- **Cube size** changes the sharp analytic box collider and visual cube together. +- **Iterations** control the number of PBD projection passes per timestep. +- **Timestep** changes simulation speed and contact robustness; faster presets are less stable and are labeled in the overlay. + +The overlay shows the current material values, timestep, and currently applied obstacle sizes. If the pending mesh, sphere radius, or cube size differs from the current state, the overlay prints `Pending changes: press R to apply`. + +### Wind Demo Controls + +The `pbd hang-wind` demo configures the cloth like a flag: + +- one side edge is pinned like cloth attached to a pole +- aerodynamic forcing uses triangle-based drag and lift +- drag, lift, gust, and noise stay as internal defaults so the interface stays focused on wind input only + +Startup options: + +- `--wind-speed value`: optional startup wind speed for `pbd hang-wind`; accepts a float in `[0, 15]` +- `--wind-dir x y z`: optional startup wind direction; defaults to `(1, 0, 0)` and is normalized internally + +Keyboard controls: + +| Key | Action | +|---|---| +| `1 / 2` | Decrease / increase wind speed | +| `3 / 4` | Decrease / increase wind direction x | +| `5 / 6` | Decrease / increase wind direction y | +| `7 / 8` | Decrease / increase wind direction z | +| `R` | Reset the cloth while keeping the current wind values | +| `T` | Reset wind values to defaults | +| `P` | Pause / resume simulation | + +Parameter meaning: + +- **Wind speed** controls the magnitude of the external wind velocity and is clamped to `[0, 15]`. +- **Wind direction** controls the direction of the applied aerodynamic force and is normalized before use. + +The overlay shows the current wind speed, the valid speed range, the current direction components, the normalized applied direction, and pause state. It also includes a framed lower-left wind widget: the boxed arrow shows the normalized horizontal `x/y` direction, and a separate `Z` arrow on the right shows whether the normalized vertical component points upward or downward. If the normalized wind direction has a large vertical component, the overlay warns that `z` should stay near `0` for flag-style motion. + +### Debug Flags + +- `--debug`: enables PBD diagnostic output and prints self-collision counters to the terminal +- `--self-thickness value`: overrides the PBD self-collision thickness with a positive float value + +--- + +## Implementation Details + +### PBD Solver Pipeline + +The PBD solver follows the standard position-based simulation pipeline: + +1. Apply external forces +2. Apply damping +3. Predict positions +4. Generate collision constraints +5. Iteratively project persistent and generated constraints +6. Update velocities from projected positions +7. Apply post-collision velocity damping/friction +8. Commit positions to the render buffer + +### Constraint Formulations + +The solver currently supports: + +- **Fixed-point constraints** for pinned and interactively dragged particles +- **Structural distance constraints** for edge-length preservation +- **Shear distance constraints** for diagonal deformation control +- **Dihedral bending constraints** for fold-angle preservation +- **Iteration-corrected stiffness** so the effective stiffness is more consistent when changing the number of solver iterations + +Stretch and shear constraints are distance-based. Bending is handled with a dihedral-angle constraint over adjacent triangle pairs rather than only using longer distance springs. + +### Collision Handling + +The solver includes several collision types: + +- **Sphere collision** as an inequality constraint +- **Plane/floor collision** as an inequality constraint +- **Analytic box collision** using face contacts and closest-point style handling +- **Edge-midpoint sampling** for improving box obstacle collision near sharp edges +- **Vertex-triangle self-collision** for cloth self-intersection handling + +Collision constraints are generated from predicted positions each timestep, then solved together with the persistent cloth constraints during projection. + +### Self-Collision Diagnostics + +When `--debug` is enabled, the solver prints self-collision statistics such as: + +- generated self-collision contacts +- initially violated contacts +- remaining violated contacts after projection +- maximum penetration before projection +- maximum penetration after projection + +These diagnostics are useful for determining whether a self-collision problem is caused by missed contact generation or insufficient projection convergence. + +### Wind Model + +The wind demo implements wind as an external aerodynamic force rather than a full wind-field solver. + +For each cloth triangle, the solver computes: + +- triangle face normal +- triangle face area +- average triangle velocity +- wind velocity +- relative velocity between the cloth face and the air + +The aerodynamic model includes: + +- **drag**, which acts opposite the relative velocity +- **lift**, which acts perpendicular to relative velocity in the plane formed by relative velocity and the face normal +- **gust modulation**, which varies the base wind speed over time +- **procedural noise**, which adds small local flutter variation + +This is a lightweight approximation inspired by the aerodynamic force model in Keckeisen et al. [4]. It does not implement a full Navier-Stokes or particle-tracing wind field. + +--- + +## Demo Notes + +The demo gallery above shows the current PBD scenes included in this project. For baseline reference, the original forked mass-spring implementation is still available and remains useful for qualitative comparison against the newer PBD demos. + +### Fast Mass-Spring Reference + +The original solver from the forked repository is based on Liu et al. [1]. It is preserved here as the baseline solver mode. + +![curtain_hang](https://user-images.githubusercontent.com/24758349/79005907-97ad1100-7b60-11ea-9e27-90375461beaf.gif) ![curtain_ball](https://user-images.githubusercontent.com/24758349/79005924-9d0a5b80-7b60-11ea-8ce4-d9fc683441d7.gif) -### References +--- + +## Project Status + +### Setup / Integration + +- [x] Build and understand existing framework +- [x] Add `PBDSolver` to project / CMake +- [x] Hook PBDSolver into app with solver mode switch +- [x] Preserve original mass-spring solver for comparison + +### Core Solver Foundation + +- [x] Implement `PBDSolver` skeleton +- [x] Initialize particle state (`x`, `p`, `v`, `invMass`) +- [x] Implement external force / gravity update +- [x] Implement damping +- [x] Implement position prediction / velocity reconstruction + +### Constraint System + +- [x] Implement fixed-point constraints +- [x] Implement structural distance constraints +- [x] Implement shear constraints +- [x] Implement dihedral bend constraints +- [x] Add stiffness / iteration-corrected stiffness handling +- [x] Add runtime parameter controls across the PBD demo set + +### Collision Handling + +- [x] Implement plane collision +- [x] Implement sphere collision +- [x] Implement analytic box collision +- [x] Add explicit collision constraint generation stage +- [x] Add self-collision detection and response +- [x] Add self-collision debug diagnostics +- [x] Tune collision robustness / epsilon + +### Rendering / Output + +- [x] Write updated positions to render buffer +- [x] Recompute normals +- [x] Verify shading / mesh integrity after simulation +- [x] Add visual support for floor, sphere, cube, and simple shadows + +### Testing / Validation + +- [x] Validate single-particle / two-particle cases +- [x] Validate hanging cloth behavior +- [x] Validate cloth drop / draping demo +- [x] Test timestep / iteration / resolution stability +- [x] Test floor-contact and self-collision diagnostics +- [ ] Compare quantitatively against original solver + +### Advanced / Stretch Goals + +- [x] Add wind / extra aerodynamic forces +- [x] Add runtime parameter controls +- [x] Add dual-obstacle collision demo +- [ ] Implement full mesh collision +- [ ] Implement full edge-edge cloth self-collision +- [ ] Benchmark performance +- [ ] Collect additional comparison metrics + +### Final Deliverables + +- [x] Stable PBD cloth demos +- [x] Updated README / documentation +- [x] Demo gallery / GIFs +- [ ] Final report / slides / citations + +--- + +## Known Limitations and Future Work + +- Self-collision is discrete, not full continuous collision detection. +- Self-collision currently uses vertex-triangle constraints; full edge-edge self-collision is not implemented. +- Box collision is analytic/proxy-based, not full mesh collision. +- Edge-midpoint sampling improves obstacle collision near sharp edges but is still an approximation. +- Wind is modeled as an external aerodynamic force, not a full Navier-Stokes or particle-tracing wind field. +- Quantitative performance benchmarking is not yet complete. + +--- + +## Credits and License + +This project is based on the original **FastMassSpring** repository by Samer Itani, licensed under the MIT License. + +Modifications and extensions, including the PBD solver, are developed as part of the Stony Brook University CSE 328 final project. + +--- + +## References + +[1] Liu, T., Bargteil, A. W., O'Brien, J. F., & Kavan, L. (2013). Fast simulation of mass-spring systems. *ACM Transactions on Graphics, 32*(6), 1–7. https://doi.org/10.1145/2508363.2508406 + +[2] Provot, X. (1995). Deformation constraints in a mass-spring model to describe rigid cloth behavior. In *Graphics Interface* 1995, 147–154. -[1] Liu, T., Bargteil, A. W., Obrien, J. F., & Kavan, L. (2013). Fast simulation of mass-spring systems. *ACM Transactions on Graphics,32*(6), 1-7. doi:10.1145/2508363.2508406 +[3] Müller, M., Heidelberger, B., Hennix, M., & Ratcliff, J. (2007). *Position Based Dynamics.* In C. Mendoza & I. Navazo (Eds.), Proceedings of the 3rd Workshop in Virtual Reality Interactions and Physical Simulation (VRIPHYS 2006). -[2] Provot, X. (1995). Deformation constraints in a mass-spring modelto describe rigid cloth behavior. *InGraphics Interface* 1995,147-154. +[4] Keckeisen, M., Kimmerle, S., Thomaszewski, B., & Wacker, M. (2004). *Modelling Effects of Wind Fields in Cloth Animations.* Journal of WSCG, 12(1–3). \ No newline at end of file