【GAMES101】作业8 显式/半隐式欧拉与Verlet方法模拟绳子运动（含GIF）

ycr的帐号

已于 2022-07-30 11:05:13 修改

阅读量2.4k

点赞数 8

分类专栏： GAMES101-计算机图形学文章标签： c++ 线性代数图形渲染 opencv

于 2022-05-05 09:43:03 首次发布

本文链接：https://blog.csdn.net/ycrsw/article/details/124565133

版权

GAMES101-计算机图形学专栏收录该内容

14 篇文章 89 订阅

订阅专栏

1. 作业描述

1.1 连接绳子的约束

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

每个结点都有质量，称为质点；质点之间的线段是一个弹簧。通过创建一系列的质点和弹簧，你就可以创建一个像弹簧一样运动的物体。

pinned_nodes 设置结点的索引。这些索引对应结点的固定属性 (pinned attribute) 应该设置为真（他们是静止的）。对于每一个结点，你应该构造一个 Mass对象，并在 Mass 对象的构造函数里设置质量和固定属性。（请仔细阅读代码，确定传递给构造函数的参数）。你应该在连续的两个结点之间创建一个弹簧，设置弹簧两端的结点索引和弹簧系数 k，请检查构造函数的签名以确定传入的参数。

运行./ropesim。你应该可以看到屏幕上画出绳子，但它不发生运动。

1.2 显式/半隐式欧拉法

胡克定律表示弹簧连接的两个质点之间的力和他们之间的距离成比例。也就是：
在这里插入图片描述

fb→a = −ks b−a||b−a||(||b − a|| − l) 在 Rope::simulateEuler 中, 首先实现胡克定律。遍历所有的弹簧，对弹簧两端的质点施加正确的弹簧力。保证力的方向是正确的！对每个质点，累加所有的弹簧力。
一旦计算出所有的弹簧力，对每个质点应用物理定律：
在这里插入图片描述

运行./ropesim。仿真应该就开始运行了，但是只有 3 个结点，看起来不够多。在application.cpp 文件的最上方，你应该可以看到欧拉绳子和 Verlet 绳子的定义。
改变两个绳子结点个数（默认为 3 个），比如 16 或者更多。

运行 ./ropesim -s 32 来设置仿真中每帧不同的仿真步数。尝试设置较小的值和较大的值（默认值为 64）。
在这里插入图片描述

1.3 显式 Verlet

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

除此之外，我们可以仿真弹簧系数无限大的弹簧。不用再考虑弹簧力，而是用解约束的方法来更新质点位置：只要简单的移动每个质点的位置使得弹簧的长度保持原长。修正向量应该和两个质点之间的位移成比例，方向为一个质点指向另一质点。每个质点应该移动位移的一半。
在这里插入图片描述

只要对每个弹簧执行这样的操作，我们就可以得到稳定的仿真。为了使运动更加平滑，每一帧可能需要更多的仿真次数。

1.4 阻尼

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

在这里插入图片描述

1.5 你应该修改的函数

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

2. 解

2.1 Rope::rope

这个构造函数内我们要做的就是根据传入的参数定义rope的masses（质点）数组和spring（弹簧）数组，它们的构造函数可以在对应的头文件中找到，masses就是传入起始位置、质点质量以及固定属性，spring就是依次传入两端质点以及弹簧系数k，构造时会自动计算原始弹簧长度rest_length，这里只把第一个点设为固定

    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.
        Vector2D step = (end - start) / (num_nodes - 1);
        for(int i = 0; i < num_nodes; i++){
            masses.push_back(new Mass(start + step * i, node_mass, true));
            if(i){
                springs.push_back(new Spring(masses[i-1], masses[i], k));
                masses[i]->pinned = false;
            }
        }
        // Comment-in this part when you implement the constructor       
    }

编写完这部分后我们运行程序可以得到一条固定的默认3个结点的绳子：
在这里插入图片描述

2.2 Rope::simulateEuler

显示/半隐式欧拉法其实是非常直观的一种模拟方法：
在这里插入图片描述
就是通过控制一定时间Δt内加速度不变来近似模拟物体的下一个位置，显示与半隐式的区别就是下一个时刻位置的计算使用的是本时刻的速度还是下一时刻的速度：

在这里我们是通过两个for循环来实现的，第一个对每段弹簧的遍历里面计算出了所有质点的受力情况，第二个对每个质点的遍历里面将受力加上了重力，最后再用欧拉法计算它们的下一个时刻位置：

    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 length = (s->m2->position - s->m1->position).norm();
            Vector2D dis = s->m2->position - s->m1->position;
            Vector2D force = s->k * (length - s->rest_length) * dis / length;
            s->m1->forces += force;
            s->m2->forces -= force;
            // damping
            Vector2D reve = s->m2->velocity - s->m1->velocity;
            Vector2D force1 = 0.05 * (reve.x * dis.x + reve.y * dis.y) * dis / length;
            s->m1->forces += force1;
            s->m2->forces -= force1;
            // air damping
            s->m1->forces -= 0.005 * s->m1->velocity;
            s->m2->forces -= 0.005 * s->m2->velocity;
        }

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

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

首先如果我们只用assignment里面提到的胡克定律来计算弹簧间质点受力的情况的话，我们会发现质点间是永远都在运动的：
在这里插入图片描述

在这里插入图片描述

所以我这里加上了弹簧内部抑制力dumping（通过两端质点的相对速度计算）
在这里插入图片描述
可见质点间的运动基本趋于平稳：

但是我发现这个运动还是不会停止的，这是因为这里忽略了空气的摩擦力，同样，在质点速度的反方向加上一个air damping，这样绳子的运动最终就会停下来了：
在这里插入图片描述
显式与半隐式欧拉的区别不大，主要是在精度需求方面有些不一样，比如我把每帧步长设置为64时，显式欧拉方法计算得到的结果就出现了很大偏差：

2.3 Rope::simulateVerlet

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)
            float length = (s->m2->position - s->m1->position).norm();
            Vector2D dis = s->m2->position - s->m1->position;
            Vector2D force = s->k * (length - s->rest_length) * dis / length;
            s->m1->forces += force;
            s->m2->forces -= force;
        }

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