GAMES101-Assignment8

一、总览

1.1 连接绳子的约束

在rope.cpp 中, 实现Rope 类的构造函数。这个构造函数应该可以创建一个
新的绳子(Rope) 对象，该对象从start 开始，end 结束，包含num_nodes 个节点。也就是如下图所示：

每个结点都有质量，称为质点；质点之间的线段是一个弹簧。通过创建一系列的质点和弹簧，你就可以创建一个像弹簧一样运动的物体。
pinned_nodes 设置结点的索引。这些索引对应结点的固定属性(pinned attribute)应该设置为真（他们是静止的）。对于每一个结点，你应该构造一个Mass对象，并在Mass对象的构造函数里设置质量和固定属性。（请仔细阅读代码，确定传递给构造函数的参数）。你应该在连续的两个结点之间创建一个弹簧，设置弹簧两端的结点索引和弹簧系数k，请检查构造函数的签名以确定传入的参数。
运行./ropesim。你应该可以看到屏幕上画出绳子，但它不发生运动。

1.2 显式/半隐式欧拉法

胡克定律表示弹簧连接的两个质点之间的力和他们之间的距离成比例。也就是：

在Rope::simulateEuler 中, 首先实现胡克定律。遍历所有的弹簧，对弹簧两端的质点施加正确的弹簧力。保证力的方向是正确的！对每个质点，累加所有的弹簧力。
一旦计算出所有的弹簧力，对每个质点应用物理定律：

运行./ropesim。仿真应该就开始运行了，但是只有3 个结点，看起来不够多。在application.cpp 文件的最上方，你应该可以看到欧拉绳子Verlet绳子的定义。改变两个绳子结点个数（默认为3个），比如16或者更多。
运行./ropesim -s 32 来设置仿真中每帧不同的仿真步数。尝试设置较小的
值和较大的值（默认值为64）。

1.3 显式Verlet

Verlet 是另一种精确求解所有约束的方法。这种方法的优点是只处理仿真中
顶点的位置并且保证四阶精度。和欧拉法不同，Verlet 积分按如下的方式来更新下一步位置：

除此之外，我们可以仿真弹簧系数无限大的弹簧。不用再考虑弹簧力，而是用解约束的方法来更新质点位置：只要简单的移动每个质点的位置使得弹簧的长度保持原长。修正向量应该和两个质点之间的位移成比例，方向为一个质点指向另一质点。每个质点应该移动位移的一半。
只要对每个弹簧执行这样的操作，我们就可以得到稳定的仿真。为了使运动更加平滑，每一帧可能需要更多的仿真次数。

1.4 阻尼

向显示Verlet方法积分的胡克定律中加入阻尼。现实中的弹簧不会永远跳动，因为动能会因摩擦而减小。阻尼系数设置为0.00005, 加入阻尼之后质点位置更新如下：

应该修改的函数是:

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

---

二、安装依赖

$ sudo apt install libglu1-mesa-dev freeglut3-dev mesa-common-dev
$ sudo apt install xorg-dev

---

三、参考答案

3.1 rope.cpp中的Rope::rope(…)

实现1.1 连接绳子的约束

    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.
        for(int i = 0; i < num_nodes; i++){// 创建质点
            Vector2D pos = start + (end - start) * (double(i) / double(num_nodes));
            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));
        }
//        Comment-in this part when you implement the constructor
       for (auto &i : pinned_nodes) {
           masses[i]->pinned = true;
       }
    }

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

实现1.2 显式/半隐式欧拉法

• 胡克定律
  

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

• 显式/隐式欧拉法
  

    if (!m->pinned)
    {
        // TODO (Part 2): Add the force due to gravity, then compute the new velocity and position
        // 显式欧拉法：下一个时刻位置用当前速度计算
        auto a = m->forces / m->mass + gravity;
        m->position += m->velocity * delta_t;//根据当前速度更新位置
        m->velocity += a * delta_t;//更新下一时刻速度
        
        // 隐式欧拉法：下一个时刻位置用下一时刻速度计算
        // auto a = m->forces / m->mass + gravity;
        // m->velocity += a * delta_t;//更新速度
        // m->position += m->velocity * delta_t;//根据下一时刻速度更新位置 

        // TODO (Part 2): Add global damping
    }

---

• 完整代码

    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
            auto len = (s->m1->position - s->m2->position).norm();
            s->m1->forces += -s->k * (s->m1->position - s->m2->position) / len * (len - s->rest_length);
            s->m2->forces += -s->k * (s->m2->position - s->m1->position) / len * (len - 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
                // 显式欧拉法：下一个时刻位置用当前速度计算
                auto a = m->forces / m->mass + gravity;
                m->position += m->velocity * delta_t;//根据当前速度更新位置
                m->velocity += a * delta_t;//更新下一时刻速度
                
                // 隐式欧拉法：下一个时刻位置用下一时刻速度计算
                // auto a = m->forces / m->mass + gravity;
                // m->velocity += a * delta_t;//更新速度
                // m->position += m->velocity * delta_t;//根据下一时刻速度更新位置 

                // TODO (Part 2): Add global damping
            }

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

---

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

实现1.3 显式Verlet

• 显式Verlet
  

    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)
            auto len = (s->m1->position - s->m2->position).norm();
            s->m1->forces += -s->k * (s->m1->position - s->m2->position) / len * (len - s->rest_length);
            s->m2->forces += -s->k * (s->m2->position - s->m1->position) / len * (len - 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
                auto 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
            }
            m->forces = Vector2D(0, 0);
        }
    }

3.4 增加阻尼

• rope.cpp 中的void Rope::simulateVerlet(…)
  向显示Verlet方法积分的胡克定律中加入阻尼。现实中的弹簧不会永远跳动，因为动能会因摩擦而减小。阻尼系数设置为0.00005, 加入阻尼之后质点位置更新如下：
  

    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)
            auto len = (s->m1->position - s->m2->position).norm();
            s->m1->forces += -s->k * (s->m1->position - s->m2->position) / len * (len - s->rest_length);
            s->m2->forces += -s->k * (s->m2->position - s->m1->position) / len * (len - 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
                auto 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
                double 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;
            }
            m->forces = Vector2D(0, 0);
        }
    }

• rope.cpp 中的void Rope::simulateEuler(…)
  向欧拉方法中加入阻尼，直接使用?kdv作为阻尼，而不是相对速度

    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
            auto len = (s->m1->position - s->m2->position).norm();
            s->m1->forces += -s->k * (s->m1->position - s->m2->position) / len * (len - s->rest_length);
            s->m2->forces += -s->k * (s->m2->position - s->m1->position) / len * (len - 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
                // 显式欧拉法：下一个时刻位置用当前速度计算
                // auto a = m->forces / m->mass + gravity;
                // m->position += m->velocity * delta_t;//根据当前速度更新位置
                // m->velocity += a * delta_t;//更新下一时刻速度
                
                // 隐式欧拉法：下一个时刻位置用下一时刻速度计算
                // auto 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.005;
                auto 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);
        }
    }

四、编译

$ mkdir build
$ cd build
$ cmake ..
$ make

通过查看rope.cpp中的Rope::rope(…)的引用，可以定位到application.cpp文件，

// Create two ropes 
// 通过修改构造函数的第三个参数指定质点个数(默认三个质点，这里改成了10个)
ropeEuler = new Rope(Vector2D(0, 200), Vector2D(-400, 200), 10, config.mass,
                    config.ks, {0});
ropeVerlet = new Rope(Vector2D(0, 200), Vector2D(-400, 200), 10, config.mass,
                    config.ks, {0});

---

附件

作业7压缩包