From 2d5ade4d7f323e5b4cb5ed27af8dbe04cbcb4ad0 Mon Sep 17 00:00:00 2001 From: mattkae Date: Sun, 30 Jan 2022 12:24:55 -0500 Subject: Fixed the damped spring simulation --- 2d/softbody/softbody_1.html | 8 +- 2d/softbody/softbody_1.html.content | 8 +- 2d/softbody/softbody_1/damped.cpp | 148 ++++++-------------------------- 2d/softbody/softbody_1/dist/output.wasm | Bin 69604 -> 58758 bytes 4 files changed, 35 insertions(+), 129 deletions(-) diff --git a/2d/softbody/softbody_1.html b/2d/softbody/softbody_1.html index 4412803..578024e 100644 --- a/2d/softbody/softbody_1.html +++ b/2d/softbody/softbody_1.html @@ -211,13 +211,13 @@ - + - + @@ -252,13 +252,13 @@ - + - + diff --git a/2d/softbody/softbody_1.html.content b/2d/softbody/softbody_1.html.content index c5fb3cc..602d0a8 100644 --- a/2d/softbody/softbody_1.html.content +++ b/2d/softbody/softbody_1.html.content @@ -157,13 +157,13 @@ - + - + @@ -198,13 +198,13 @@ - + - + diff --git a/2d/softbody/softbody_1/damped.cpp b/2d/softbody/softbody_1/damped.cpp index f56f346..0525e25 100644 --- a/2d/softbody/softbody_1/damped.cpp +++ b/2d/softbody/softbody_1/damped.cpp @@ -25,13 +25,6 @@ namespace Damped { void unload(); }; - enum DampedSpringType { - None = 0, - Overdamped = 1, - Underdamped = 2, - Critical = 3 - }; - struct DampedSpring { DampedSpringWeight* weight; @@ -42,24 +35,14 @@ namespace Damped { int32 numVertices = 0; // Initialization variables - float32 initialDisplacement = 0.f; float32 k = 4; // DampedSpring Constant, in N / m (Hooke's Law) float32 c = 1.f; // Viscous damping coefficient (Damping force) - // Discovered during initialization - DampedSpringType type = DampedSpringType::None; - float32 R = 2.f; - float32 gamma = 6.2; - float32 discriminant = 0.f; - float32 omega1 = 0.f; - float32 c1 = 0.f; - float32 c2 = 0.f; - float32 r1 = 0.f; - float32 r2 = 0.f; - // Update variables - float32 displacement = 0.f; - float32 timeElapsed = 0.f; + float32 force = 0.f; + float32 acceleration = 0.f; + float32 velocity = 0.f; + float32 position = 0.f; void load(Renderer2d* renderer, DampedSpringWeight* inWieight, float32 length, float32 loopRadius, float32 initialDisplacement, float32 inK, float32 inC); @@ -122,7 +105,8 @@ namespace Damped { void DampedSpringWeight::load(Renderer2d* renderer, float32 inMass, Vector4 startColor, Vector4 endColor) { mass = inMass; - radius = mass * 16.f; + radius = mass * 8.f; + if (radius > 42.f) radius = 42.f; const int32 numSegments = 96; const float32 radiansPerSegment = (2.f * PI) / static_cast(numSegments); const int32 numVertices = numSegments * 3; @@ -168,8 +152,7 @@ namespace Damped { const float32 epsilon = 0.0001f; - void DampedSpring::load(Renderer2d* renderer, DampedSpringWeight* inWeight, float32 length, float32 loopRadius, float32 inInitialDisplacement, float32 inK, float32 inC) { - initialDisplacement = inInitialDisplacement; + void DampedSpring::load(Renderer2d* renderer, DampedSpringWeight* inWeight, float32 length, float32 loopRadius, float32 initialDisplacement, float32 inK, float32 inC) { weight = inWeight; k = inK; c = inC; @@ -179,92 +162,22 @@ namespace Damped { float32 B = c; float32 C = k; - - discriminant = B * B - (4 * A * C); - timeElapsed = 0; - displacement = 0; + float32 discriminant = B * B - (4 * A * C); + acceleration = 0; + velocity = 0; + position = initialDisplacement; if (discriminant < epsilon && discriminant > -epsilon) { // Real repeated root (~ zero): Critically damped motion printf("Critically damped motion.\n"); - - type = DampedSpringType::Critical; - - r1 = -B / (2 * A); - r2 = r1; - - // Solve for the initial displacement. - c2 = (initialDisplacement - r1) / (r2); - c2 = r2 * c2; - c1 = r1; - printf("Equation: %f e^(%f t) + %f t e^(%f t)\n", c1, r1, c2, r2); } else if (discriminant > 0) { // Two real roots (greater than zero): Overdamped motion printf("Overdamped motion.\n"); - - type = DampedSpringType::Overdamped; - - // Discover the two treal roots. Both of these WILL be negative. - float32 sqrtDis = sqrtf(discriminant); - r1 = (-B + sqrtDis) / (2 * A); - r2 = (-B - sqrtDis) / (2 * A); - - // Now, solve for c1 and c2 given the initial displacement. The equation is: - // x(0) = r1 * c1 * e^(r1 * t) + r2 * c2 * e^(r2 * t) -> - // x(0) = r1 * c1 + r2 * c2 -> - // initialDisplacement = r1 * c1 + r2 * c2; -> if abs(r1) > abs(r2) => c1 = 1, else c2 = 1 - if (r1 < r2) { - c2 = (initialDisplacement - r1) / (r2); - c2 = r2 * c2; - c1 = r1; - } - else { - c1 = (initialDisplacement - r2) / (r1); - c1 = r1 * c1; - c2 = r2; - } - - printf("Equation: %f e^(%f t) + %f e^(%f t)\n", c1, r1, c2, r2); } else { // Complex pair (less than zero): Underdamped motion printf("Underdamped motion.\n"); - - type = DampedSpringType::Underdamped; - float32 sqrtDis = sqrtf(-discriminant); - omega1 = sqrtf(sqrtDis) / (2.f * A); // Note that we negate the discriminant to get the REAL positive part. - - r1 = (-B + sqrtDis) / (2 * A); - r2 = (-B - sqrtDis) / (2 * A); - - // Now, solve for c1 and c2 given the initial displacement. The equation is: - // x(0) = r1 * c1 * e^(r1 * t) + r2 * c2 * e^(r2 * t) -> - // x(0) = r1 * c1 + r2 * c2 -> - // initialDisplacement = r1 * c1 + r2 * c2; -> if abs(r1) > abs(r2) => c1 = 1, else c2 = 1 - if (r1 < r2) { - c2 = (initialDisplacement - r1) / (r2); - c2 = r2 * c2; - c1 = r1; - } - else { - c1 = (initialDisplacement - r2) / (r1); - c1 = r1 * c1; - c2 = r2; - } - - // At this point, we have: - // x(t) = c1 cos(omega1 t) + c2 cos(omega1 t) -> - // x(t) = A cos(omega1 t - Φ) -> - // Now we need to solve for A and Φ for our initial displacement, so: - // initialDisplacement = c1 - R = sqrtf(c1 * c1 + c2 * c2); - if (initialDisplacement < 0) { - R = -R; - } - gamma = atanf(c2 / c1); - - printf("Equation: %f e^(-(%f / %f) t) * cos(%f t - %f)\n", R, c, weight->mass, omega1, gamma); } const int32 verticesPerSegment = 6; @@ -279,13 +192,20 @@ namespace Damped { const float32 offset = 0.25f; int32 vidx = 0; + float32 dx = position; for (int pidx = 0; pidx < numSegments; pidx++) { float32 y1 = lengthIncrement * pidx; float32 x1 = loopWidth * sinf(frequency * y1 + offset); float32 y2 = y1 + lengthIncrement; float32 x2 = loopWidth * sinf(frequency * y2 + offset); - + + float32 y1Offset = dx * (1.f - pidx / static_cast(numSegments)); + float32 y2Offset = dx * (1.f - (pidx + 1) / static_cast(numSegments)); + + y1 += y1Offset; + y2 += y2Offset; + vertices[vidx++].position = Vector2(x1, y1); vertices[vidx++].position = Vector2(x1, y2); vertices[vidx++].position = Vector2(x2, y1); @@ -297,31 +217,17 @@ namespace Damped { shape.load(vertices, numVertices, renderer, GL_DYNAMIC_DRAW); shape.model = Mat4x4().translateByVec2(Vector2(400, 300)); - weight->shape.model = shape.model.translateByVec2(Vector2(0, -weight->radius)); + weight->shape.model = shape.model.translateByVec2(Vector2(0, -weight->radius + dx)); } void DampedSpring::update(float32 dtSeconds) { - timeElapsed += dtSeconds; - - float32 lastDisplacement = displacement; - - switch (type) { - case DampedSpringType::Critical: - displacement = c1 * pow(E, r1 * timeElapsed) + c2 * timeElapsed * pow(E, r2 * timeElapsed); - break; - case DampedSpringType::Overdamped: - displacement = c1 * pow(E, r1 * timeElapsed) + c2 * pow(E, r2 * timeElapsed); - break; - case DampedSpringType::Underdamped: { - float32 exponent = -(c / (2 * weight->mass)) * timeElapsed; - displacement = R * pow(E, exponent) * (cosf(omega1 * timeElapsed - gamma)); - break; - } - default: - break; - } - - float32 dx = displacement - lastDisplacement; + force = -weight->mass * acceleration - c * velocity - k * position; + acceleration = (force * dtSeconds) / weight->mass; + velocity = velocity + acceleration * dtSeconds; + float32 lastPosition = position; + position = position + velocity * dtSeconds; + + float32 dx = position - lastPosition; int32 vidx = 0; for (int pidx = 0; pidx < numSegments; pidx++) { float32 y1Offset = dx * (1.f - pidx / static_cast(numSegments)); diff --git a/2d/softbody/softbody_1/dist/output.wasm b/2d/softbody/softbody_1/dist/output.wasm index b034adc..8f7c621 100755 Binary files a/2d/softbody/softbody_1/dist/output.wasm and b/2d/softbody/softbody_1/dist/output.wasm differ -- cgit v1.2.1