How to enforce a constraint between two points

Question

I would like to apply a constraint between two points, whilst applying gravity. The following image I drew demonstrates the start and end positions of point 2, which doesn't include the intermediate time-step positions, and assumes point 1 has a fixed position:

I have a point class defined as follows:

class Point{
  glm::vec3 position;
  glm::vec3 op; // original position
  glm::vec3 velocity;
  float mass;
};

I can define two points and find the original length between two points by using the following:

Point p1;
p1.position = glm::vec3(0, 10, 0);
p1.op = p1.position;
p1.velocity = glm::vec3(0, 0, 0);
p1.mass = 1.0f;

Point p2;
p2.position = glm::vec3(10, 10, 0);
p2.op = p2.position;
p2.velocity = glm::vec3(0, 0, 0);
p2.mass = 1.0f;

float original_length_p1_p2 = glm::length(p2.op- p1.op);

I have an update function inside of the point class which is ran within a certain time-step which should update the point position by applying gravity:

glm::vec3 gravity(0,-9.8,0);
...
void update(float dt){
  velocity += gravity * dt;
  position += velocity * dt;
}

The points are stored inside of a vector and the update function is called like below:

std::vector<Point> myPoints;
...
for(int n = 0; n < myPoints.size(); n++){
  myPoints[n].update(dt);
}

Now I want to be able to apply some spring-like constraint between these two points which would swing like a simple spring-like pendulum. I have tried adding the following to the above for loop:

void applyConstraint(Point &p1, Point &p2, float dt){
    float change = (glm::length(p1.position-p2.position) - glm::length(p1.op-p2.op)) / glm::length(p1.position-p2.position);
    p1.position -= 0.5 * (p1.position-p2.position) * change * dt;
    p2.position += 0.5 * (p1.position-p2.position) * change * dt;
}

But when trying this, p2 falls with no constraint. How could I ensure p2 falls similarly to the image?

Updated applyConstraint:

void Scene::applyConstraint(Point &p1, Point &p2, float dt) {
    float change = (glm::length(p1.position - p2.position) - glm::length(p1.op - p2.op)) / glm::length(p1.position - p2.position);
    glm::vec3 force = 0.5f * (p1.position - p2.position) * change * dt;
    glm::vec3 accel1 = (-force / p1.mass) * dt;
    glm::vec3 accel2 = (force / p2.mass) * dt;
    p1.velocity += accel1 * dt;
    p2.velocity += accel2 * dt;
    p1.position += p1.velocity * dt;
    p2.position += p2.velocity * dt;
}

Answer 1

There is three issues in your code. First, you apply Euler integration for each constraint, but it should be applied only once at the end of each iteration. Second, The point p1 should be fixed. Third, you did not consider the mass in the force calculations.

To fix it, add a force vector in the Point structure and use this code:

// Reset forces
p1.force = glm::vec3(0, 0, 0);
p2.force = glm::vec3(0, 0, 0);

// Add gravity
p1.force += gravity / p1.mass ;
p2.force += gravity / p2.mass ;

// Add spring forces
// To be put in applyConstraint, without dependency on dt
float k = 1 ;
glm::vec3 difference = p1.position - p2.position;
float current_length = glm::length(difference);
float original_length = glm::length(p2.op- p1.op);
float displacement = (current_length - original_length) / current_length;
p1.force -= k * displacement * difference ;
p2.force += k * displacement * difference ;

// Euler integration
p1.velocity += p1.force / p1.mass * dt ;
p2.velocity += p2.force / p2.mass * dt ;
//p1.position += p1.velocity * dt ; // This point is an anchor
p2.position += p2.velocity * dt ;

Change k to adjust the elasticity of the spring. If you know the behavior that you want to have, compute it using the formula given on this website.

You can also add damping to the system using p2.force -= c * p2.velocity, where c is the damping ratio.