Soft body engine in C++

Question

I am attempting to make an elementary soft body engine in C++ using SDL2. It works by considering all the vertices of the soft body to be interconnected by springs of same length and rigidity (with the same spring constant k and length natural_length). To make it more realistic, I also introduced a damping constant c.
However, I encountered a frustrating problem. I have been trying to debug it for past 6-7 hours but to no avail. The soft body encounters many strange bugs which I don't understand

• Firstly, the "soft body" is not at all "soft". It becomes a crumpled mess of points every time. I have tried calculating the force of neighbouring points instead only, but it still becomes a crumpled mess.
• The soft body flies up to the top corner (origin) every time even though I haven't implemented any external forces.

Both of these bugs are visible in this image -

The 2 following functions (they are in the same class as with all the variables, hence need not take in any arguments) are the actual simulation part of the code. (I have omitted the rest of the code as it is unnecessary.)
I used a vector of SDL_Points to store every point and a vector of Vector to store their velocities. If you are wondering what Vector is, it is simply a struct I created which simply has 2 float members x and y.
The acceleratePoints() function assigns velocity and position to each point and checkCollision() well... checks for collisions with the walls of the window, which has width scr_w and height scr_h.

    void acceleratePoints()
    {
        vector<SDL_Point> soft_body_copy=soft_body;
        vector<Vector> velocity_copy=velocity;
        for(int i=0;i<soft_body.size();++i)
        {
            for(int j=0;j<soft_body.size();++j)
            {
                if(i!=j)
                {
                    Vector d={(soft_body[j].x-soft_body[i].x)/100.0,(soft_body[j].y-soft_body[i].y)/100.0};
                    float t=atan2(d.y,d.x);
                    float disp=fabs(magnitude(d))-natural_length/100.0;
                    velocity_copy[i].x+=(k*disp*cos(t))/10000.0;
                    velocity_copy[i].y+=(k*disp*sin(t))/10000.0;
                    velocity_copy[i].x-=c*velocity_copy[i].x/100.0;
                    velocity_copy[i].y-=c*velocity_copy[i].y/100.0;
                    soft_body_copy[i].x+=velocity_copy[i].x;
                    soft_body_copy[i].y+=velocity_copy[i].y;
                }
            }
            soft_body=soft_body_copy;
            velocity=velocity_copy;
        }
    }
    void checkCollision()
    {
        for(int k=0;k<soft_body.size();++k)
        {
            if(soft_body[k].x>=scr_w||soft_body[k].x<=0)
            {
                velocity[k].x*=e;
                soft_body[k].x=soft_body[k].x>scr_w/2?scr_w-1:1;
            }
            if(soft_body[k].y>=scr_h||soft_body[k].y<=0)
            {
                velocity[k].y*=e;
                soft_body[k].y=soft_body[k].y>scr_h/2?scr_h-1:1;
            }
        }
    }

The magnitude() function returns the magnitude of a Vector.
The values for coefficient of restitution e, damping constant c and spring constant k, which I used for the image are 0.5, 10 and 100 respectively. Thank you for taking the time to read this! Help would be greatly appreciated.

Edit

Here is the entire code if anyone wants to test it. You'll need SDL and a folder 'img' with a '.bmp' file in 'img/point.bmp'.

Answer 1

Spring-based soft body simulation requires a "rest configuration", wherein without gravity or other external forces, the body will remain internally motionless.

Each spring in the simulation needs its own length. A shared natural_length will cause springs to exert forces within the body where the input separation between connected points is different, which is what appears to be causing the problems you describe.

A typical approach to generate a soft-body (or "point blob") out of a set of points looks like this:

• Generate the Delaunay Triangulation (or Tetrahedralization, in 3D) of the point set.
• Create springs on every edge of the triangulation, using the length of the edge as the spring's rest length.
• Simulate with "gravity" and other external forces!

For simplicity, you can just connect all pairs of points, rather than generating a triangulation.

The simulation itself can be performed in three steps:

• Clear vertex force accumulators
• Add spring forces to their connected vertices
• Update the velocity and position of each vertex

In the spirit of the code you've posted, a Spring might look like:

struct Spring
{
    int point_index[2];
    float rest_length;
};

ASIDE

In your code, you can replace t = atan2(d.y, d.x) and the following sin(t) cos(t) with a unit-length image of d:

float mag_d = magnitude(d); // should not become negative...
float r = (0.0f != mag_d) ? 1.0f/mag_d : 0.0f;
Vector d_hat{d.x*r, d.y*r};
float w = k * (mag_d - spring[i].rest_length);
velocity_copy[i].x += w*d_hat.x;
velocity_copy[i].y += w*d_hat.y;

Not completely optimized, but just a little more stable.

EDIT

I've also noticed you were scaling the rest length (dividing it by 100). Don't do that.

ANOTHER EDIT

Change this:

    }
    soft_body=soft_body_copy;
    velocity=velocity_copy;
}

To this:

    }
}
soft_body = move(soft_body_copy);
velocity = move(velocity_copy);

You were changing vertex positions while calculating forces.