1
$\begingroup$

I am trying to implement a cloth simulation in c++ using verlet integration and I'm getting some odd issues that I can't seem to explain, the points in the cloth seem to fly around and have weird behavior.

to investigate I tried to simply place 2 dots next to each other horizontally with a single constraint and to pin one of the points, this should give the effect of the second point falling until both points are vertical.

what I am seeing is the unpinned point is swinging around the first one faster and faster, the odd thing is that if I am moving the simulation step by step using the keyboard I get the desired result what makes me think it might be a timestamp issue..

here is the code, I used glm for calculations and SFML for display.

#include <iostream>

#include <glm/glm.hpp>
#include <vector>
#include <functional>

#include <SFML/Graphics.hpp>

using vec3 = glm::vec3;

#define NUM_OF_ITERATIONS 1
#define TIMESTAMP 0.01f
#define WIDTH 800
#define HEIGHT 600
#define drag 0.001f

#define logInfo(msg) std::cout << msg << std::endl;
#define logDebug(msg) std::cout << msg << std::endl;
#define logWarning(msg) std::cout << msg << std::endl;
#define logError(msg) std::cout << msg << std::endl;

struct Particle
{
    vec3 m_currentPos;
    vec3 m_previousPos;
    vec3 m_forceAccumulations;
};

struct Constraint
{
    int particleA;
    int particleB;
    float restLength;
};

class ParticleSystem
{
public:
    ParticleSystem();
    void step(float dt);
    void addConstraint(const Constraint& c);
    void addParticle(const Particle& p);
    void draw(std::function<void(const Particle&)> cb) const;
private:
    void verlet(float dt);
    void satisfyContraints();
    void accumulateForces();


private:
    vec3 m_gravity = {0, 9.81f, 0};
    float m_timestamp = TIMESTAMP;

    std::vector<Particle> m_particles;
    std::vector<Constraint> m_constraints;
};

ParticleSystem::ParticleSystem()
{

}

void ParticleSystem::step(float dt)
{
    accumulateForces();
    verlet(dt);
    satisfyContraints();
}

void ParticleSystem::verlet(float dt)
{
    for (auto& p : m_particles)
    {
        vec3& x = p.m_currentPos;
        vec3 temp_x = x;
        vec3& old_x = p.m_previousPos;
        vec3& a = p.m_forceAccumulations;
        x += (x - old_x) + a * (dt * dt);
        old_x = temp_x;
    }
}

void ParticleSystem::satisfyContraints()
{
    for (int i = 0; i < NUM_OF_ITERATIONS; i++)
    {
        //satisfy bounds
        for (auto& p : m_particles)
        {
            vec3& v = p.m_currentPos;
            v = glm::min(glm::max(v, { 0,0,0 }), { WIDTH,HEIGHT,HEIGHT });
        }

        // satisfy contsraints
        for (const auto& c : m_constraints)
        {
            vec3& v1 = m_particles[c.particleA].m_currentPos;
            vec3& v2 = m_particles[c.particleB].m_currentPos;
            vec3 delta = v2 - v1;
            float deltaLength = glm::length(delta);
            float diff = (deltaLength - c.restLength) / deltaLength;
            v1 += delta * .5f * diff;
            v2 -= delta * .5f * diff;
        }

        // Constrain one particle of the cloth to orig
        m_particles[0].m_currentPos = { WIDTH / 2 , HEIGHT / 2, 0 };

    }
}

void ParticleSystem::accumulateForces()
{
    for (auto& p : m_particles)
    {
        p.m_forceAccumulations = m_gravity;
    }
}

void ParticleSystem::addConstraint(const Constraint& c)
{
    m_constraints.push_back(c);
}

void ParticleSystem::addParticle(const Particle& p)
{
    m_particles.push_back(p);
}

void ParticleSystem::draw(std::function<void(const Particle& particle)> cb) const
{
    for (const auto& p : m_particles)
    {
        cb(p);
    }
    
}

int main()
{
    ParticleSystem ps;

    vec3 origin = { WIDTH / 2, HEIGHT / 2, 0 };

    ps.addParticle({ origin, origin, vec3{ 0, 0, 0 } });
    ps.addParticle({ origin + vec3{100,0,0}, origin + vec3{ 100,0,0 }, vec3{ 0, 0, 0 } });

    ps.addConstraint({ 0, 1, 100 });
    
    sf::RenderWindow window(sf::VideoMode({ WIDTH, HEIGHT }), "Cloth simulation", sf::Style::Close);

    while (window.isOpen())
    {
        window.clear();

        // Process events
        for (sf::Event event; window.pollEvent(event);)
        {
            // Close window: exit
            if (event.type == sf::Event::Closed)
                window.close();

            // Escape key: exit
            if ((event.type == sf::Event::KeyPressed) && (event.key.code == sf::Keyboard::Escape))
                window.close();

            // Escape key: exit
            if ((event.type == sf::Event::KeyPressed) && (event.key.code == sf::Keyboard::A))
                ps.step(1.f);
        }
        // WHEN I UNCOMMENT THIS THIS SIMULATION GOES INSANE
        //ps.step(0.01f);

        ps.draw([&](const Particle& particle) {
            sf::Vertex point[1];
            point[0].position = sf::Vector2f(particle.m_currentPos.x, particle.m_currentPos.y);
            point[0].color = sf::Color::Red;
            window.draw(point, 1, sf::Points);
        });

        // Finally, display the rendered frame on screen
        window.display();
    }
}

------------- Update -------------
following @Lutz Lehmann suggestion it appears the framerate was simply too high what caused a rounding error. I tried limiting the frame rate and it did the trick.

$\endgroup$
6
  • 1
    $\begingroup$ You are hard-coding the time step to a somewhat large value, have you tried decreasing the time step? $\endgroup$
    – Kyle Kanos
    Nov 26, 2022 at 18:59
  • $\begingroup$ @KyleKanos makes a good suggestion about reducing the timestep. I only skimmed the code, but it also looks like you're enforcing the constraints by clamping some variables. This could be another source of error -- it's worth looking in the Hairer book on methods for constrained Hamiltonian systems, in particular SHAKE and RATTLE. Writing an energy-conserving method for constrained problems is pretty finicky. I wrote a bit about it here as well. $\endgroup$ Nov 26, 2022 at 19:22
  • $\begingroup$ @KyleKanos I tried decreasing the timestamp to about 0.001f but this causes the simulation to not function due to a rounding error.. $\endgroup$ Nov 26, 2022 at 21:13
  • $\begingroup$ It is possible that in the automatic mode the effective frame rate is so high that you see the unavoidable error accumulation much faster than in manual mode. With float numbers the dt*dt factor is critical in introducing rounding errors, thus losing accuracy. Switching to velocity or leapfrog Verlet distributes this factor as a simple dt factor in position and velocity. $\endgroup$ Nov 27, 2022 at 15:28
  • $\begingroup$ @LutzLehmann You are a life saver! limiting the framerate actually solved the issue, thanks a lot! $\endgroup$ Nov 27, 2022 at 19:12

1 Answer 1

2
$\begingroup$

following @Lutz Lehmann suggestion it appears the framerate was simply too high what caused a rounding error. I tried limiting the frame rate and it did the trick.

$\endgroup$
1
  • $\begingroup$ Just want to say, this help me immediately ! Thank you ! $\endgroup$
    – silverfox
    Jul 30, 2023 at 6:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.