实现要求：

  • 连接绳子约束，正确的构造绳子
  • 半隐式欧拉法
  • 显式欧拉法
  • 显式 Verlet
  • 阻尼

你应该修改的函数是:

  • rope.cpp 中的 Rope::rope(…)
  • rope.cpp 中的 void Rope::simulateEuler(…)
  • rope.cpp 中的 void Rope::simulateVerlet(…)

环境配置

本次作业只能在Linux上面进行，需要在作业路径输入以下指令，将会自动安装依赖。需要安装OpenGL, Freetype，RandR这几样依赖。

1
2
sudo apt install libglu1-mesa-dev freeglut3-dev mesa-common-dev
sudo apt install xorg-dev

构造绳子

根据代码框架，可以得知一个Rope由num_nodes个节点和num_nodes-1条线段组成，我们将每个节点视为带质量的质点，质点之间的线段是一个弹簧，因此，我们需要为每个质点创建它的位置信息以及质量信息，并且为每个线段指定它的头节点和尾节点以及弹簧系数k，然后让每个顶点被视为pinned，即无法移动的状态。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
Rope::Rope(Vector2D start, Vector2D end, int num_nodes, float node_mass, float k, vector<int> pinned_nodes)
{
// TODO (Part 1): Create a rope starting at `start`, ending at `end`, and containing `num_nodes` nodes.

//        Comment-in this part when you implement the constructor
//        for (auto &i : pinned_nodes) {
//            masses[i]->pinned = true;
//        }

    for (int i = 0; i < num_nodes; i++)
    {
        Vector2D pos = start + (end - start) * ((double)i / ((double)num_nodes - 1.0));
        masses.push_back(new Mass(pos, node_mass, false));

    }

    for (int i = 0; i < num_nodes - 1; i++)
    {
        springs.push_back(new Spring(masses[i], masses[i+1], k));
    }

    for (auto &i : pinned_nodes) {
            masses[i]->pinned = true;
    }

}

按照下面命令创建工程：

1
2
3
4
$ mkdir build
$ cd build
$ cmake . .
$ make

之后，可以使用命令./ropesim来运行仿真，我们可以看到我们构造出来了这样一条绳子，但是他现在还不能移动，我们需要完善之后的代码让他动起来。

image-20231102220412878

显式/半隐式欧拉法

使用胡克定律来表示弹簧连接的两个质点之间的力以及他们之间的距离的比例，这里的l是未发生弹性形变的弹簧长度，在这里为s->rest_length。
$$
f_{b\to a} = -K_s \frac{b-a}{\vert \vert b-a \vert \vert}(\vert \vert b-a \vert \vert - l)
$$
下面的代码根据胡克定律为每个绳子上的质点添加作用力：

1
2
3
4
5
6
7
for (auto &s : springs)
{
    // TODO (Part 2): Use Hooke's law to calculate the force on a node
    float l = (s->m1->position - s->m2->position).norm();
    s->m1->forces += -(s->k) * (s->m1->position - s->m2->position) / l * (l - s->rest_length);
    s->m2->forces += -(s->k) * (s->m2->position - s->m1->position) / l * (l - s->rest_length);
}

接着我们要根据加速度来更新每个质点的下一时刻的位置
$$
F=ma \
v(t+1) = v(t) + a(t) * dt \
x(t+1) = x(t) + v(t) * dt \
x(t+1) = x(t) + v(t+1) * dt
$$
下面使用半隐式欧拉法编写代码，首先计算出加速度a，然后使用加速更新速度，最后使用速度更新加速度。

1
2
3
4
5
6
7
8
9
10
11
12
13
for (auto &m : masses)
{
    if (!m->pinned)
    {
        // TODO (Part 2): Add the force due to gravity, then compute the new velocity and position
        Vector2D a = m->forces / m->mass + gravity;
        m->velocity += a * delta_t; 
        m->position += m->velocity * delta_t; 
    }

    // Reset all forces on each mass
    m->forces = Vector2D(0, 0);
}

此外，为了得到更好的视觉效果，可以在application.cpp中改变质点的个数。

在application.cpp这里可以得知，蓝色绳子使用的是Euler方法，绿色绳子使用的是Verlet方法。

1
2
3
4
5
6
7
8
for (int i = 0; i < 2; i++) {
if (i == 0) {
  glColor3f(0.0, 0.0, 1.0);
  rope = ropeEuler;
} else {
  glColor3f(0.0, 1.0, 0.0);
  rope = ropeVerlet;
}

image-20231103013426210

显式 Verlet

Verlet 是另一种精确求解所有约束的方法。这种方法的优点是只处理仿真中顶点的位置并且保证四阶精度。和欧拉法不同，Verlet 积分按如下的方式来更新下一步位置:
$$
x(t+1) = x(t) + [x(t) - x(t-1)] + a(t) * dt * dt
$$
代码编写中，我们先实现每个弹簧上受的作用力，然后对于每个质点，我们记录一个中间时刻的位置，然后计算加速度，使用上面的公式计算下一时刻的位置，最后将中间时刻的位置作为上一时刻的位置供下次更新时使用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
void Rope::simulateVerlet(float delta_t, Vector2D gravity)
{
    for (auto &s : springs)
    {
        // TODO (Part 3): Simulate one timestep of the rope using explicit Verlet （solving constraints)
        float l = (s->m1->position - s->m2->position).norm();
        s->m1->forces += -(s->k) * (s->m1->position - s->m2->position) / l * (l - s->rest_length);
        s->m2->forces += -(s->k) * (s->m2->position - s->m1->position) / l * (l - s->rest_length);

    }

    for (auto &m : masses)
    {
        if (!m->pinned)
        {
            Vector2D temp_position = m->position;
            // TODO (Part 3.1): Set the new position of the rope mass
            Vector2D a = m->forces / m->mass + gravity;
            m->position =  temp_position + (temp_position - m->last_position) + a * delta_t * delta_t;
            m->last_position = temp_position;
        }
        // Reset all forces on each mass
        m->forces = Vector2D(0, 0);
    }
}

阻尼

向显式 Verlet加入阻尼，只要定义一个阻尼系数设置为 0.00005，然后更新质点的方向即可，

加入阻尼系数的公式：
$$
x(t+1) = x(t) + (1 - dampingfactor) * [x(t) - x(t-1)] + a(t) * dt * dt
$$
然后根据公式敲即可，完整代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
void Rope::simulateVerlet(float delta_t, Vector2D gravity)
{
    for (auto &s : springs)
    {
        // TODO (Part 3): Simulate one timestep of the rope using explicit Verlet （solving constraints)
        float l = (s->m1->position - s->m2->position).norm();
        s->m1->forces += -(s->k) * (s->m1->position - s->m2->position) / l * (l - s->rest_length);
        s->m2->forces += -(s->k) * (s->m2->position - s->m1->position) / l * (l - s->rest_length);

    }

    for (auto &m : masses)
    {
        if (!m->pinned)
        {
            Vector2D temp_position = m->position;
            // TODO (Part 3.1): Set the new position of the rope mass
            //Vector2D a = m->forces / m->mass + gravity;
            //m->position =  temp_position + (temp_position - m->last_position) + a * delta_t * delta_t;
            //m->last_position = temp_position;

            // TODO (Part 4): Add global Verlet damping
            Vector2D a = m->forces / m->mass + gravity;
            float damping_factor = 0.00005;
            m->position = temp_position + (1 - damping_factor) * (temp_position - m->last_position) + a * delta_t * delta_t;
            m->last_position = temp_position;

        }
        // Reset all forces on each mass
        m->forces = Vector2D(0, 0);
    }
}

而向半隐式欧拉法中加入阻尼，需要设置一个kd，然后将这个系数用于影响加速度，最终的代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
void Rope::simulateEuler(float delta_t, Vector2D gravity)
{
    for (auto &s : springs)
    {
        // TODO (Part 2): Use Hooke's law to calculate the force on a node
        float l = (s->m1->position - s->m2->position).norm();
        s->m1->forces += -(s->k) * (s->m1->position - s->m2->position) / l * (l - s->rest_length);
        s->m2->forces += -(s->k) * (s->m2->position - s->m1->position) / l * (l - s->rest_length);
    }

    for (auto &m : masses)
    {
        if (!m->pinned)
        {
            // TODO (Part 2): Add the force due to gravity, then compute the new velocity and position
            //Vector2D a = m->forces / m->mass + gravity;
            //m->velocity += a * delta_t; 
            //m->position += m->velocity * delta_t; 


            // TODO (Part 2): Add global damping
            float kd = 0.05;
            Vector2D a = m->forces / m->mass + gravity - kd * m->velocity / m->mass;
            m->velocity += a * delta_t; 
            m->position += m->velocity * delta_t;
        }

        // Reset all forces on each mass
        m->forces = Vector2D(0, 0);
    }
}

加入阻尼之后，欧拉法衰减的速度会比显示Verlet方法衰减的速度要快。

image-20231103013656256