【GAMES101】作业8——质点弹簧系统

作业内容：

1.rope():实现绳子的构造函数

2.simulateEuler()：用显示和半隐式两种方式实现欧拉法

3.simulateVerlet()：实现显示Verlet并添加阻尼

作业实现：

1.rope()构造函数

按照rope()的定义rope是由一组Mass(质点)和一组Spring(弹簧)组成并存于vector中，故在构造rope是需要依据起始点的位置、终点的位置以及质点的数量生成Mass并存于vector masses中，然后在质点间生成Spring存于vector springs中。

查看Mass定义在生成质点的时候需要知道其position位置、mass质量、pinned是否被固定

struct Mass {
  Mass(Vector2D position, float mass, bool pinned)
      : start_position(position), position(position), last_position(position),
        mass(mass), pinned(pinned) {}
  //质量
  float mass;
  //是否固定
  bool pinned;
  //起始位置与位置
  Vector2D start_position;
  Vector2D position;
  // explicit Verlet integration Verlet积分法
  //上个时间的位置
  Vector2D last_position;
  // explicit Euler integration 欧拉方法
  //速度velocity与力forces
  Vector2D velocity;
  Vector2D forces;
};

查看Spring定义知道起止质点与弹簧系数即可

struct Spring {
  Spring(Mass *a, Mass *b, float k)
      : m1(a), m2(b), k(k), rest_length((a->position - b->position).norm()) {}
  //弹簧系数
  float k;
  //放松情况长度
  double rest_length;
  //起止质点
  Mass *m1;
  Mass *m2;
};

rope构造函数如下，其中计算每个质点位置的公式为：i*(end-start)/(n-1)+start

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++) {
            //生成质点
            masses.emplace_back(new Mass(start + i * (end - start) / (num_nodes - 1),node_mass,false));
            if (i != 0) {
                //生成弹簧
                springs.emplace_back(new Spring(masses[i-1],masses[i],k));
            }
        }
        //将pinned_nodes里面的质点标记为不可移动
        for (auto &i : pinned_nodes) {
            masses[i]->pinned = true;
        }
    }

2.simulateEuler()欧拉方法

先按照公式计算弹簧spring对于质点的力，然后对于每个质点计算并加上其自身的重力，最后根据力的总和求出其加速度并计算速度与位移。

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
            //胡克定律
            Vector2D b_a=s->m2->position-s->m1->position;
            Vector2D force=-s->k*b_a.unit()*(b_a.norm()-s->rest_length);
            s->m1->forces -= force;
            s->m2->forces += force;
            
        }

        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->mass;//G=mg
                //添加摩擦力
                float k_d=0.1;
                m->forces -= k_d*m->velocity;
                // TODO (Part 2): Add global damping
                //计算加速度方向
                Vector2D a = m->forces/m->mass;
                //先计算velocity为半隐式方法
                m->velocity  += a*delta_t;
                m->position += m->velocity*delta_t;
                //先计算位置position为显示方法
                //m->position += m->velocity*delta_t;
                //m->velocity  += a*delta_t;
            }
            // Reset all forces on each mass  重置force
            m->forces = Vector2D(0, 0);
        }
    }

3.simulateVerlet() Verlet方法

按照公式，依然需要计算弹簧的力、重力来计算加速度，但是在计算位移的时候不用计算速度，难点在于需要记录上一个时间的位置position x(t-1),在Mass中已经定义了last_position用于记录上一个时间的位置，利用temp记录位置，计算完成后再用last_position记录即可。

考虑阻尼的情况下，需要再[x(t)-x(t-1)]上乘以[1-dump]。

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

        for (auto &m : masses)
        {
            if (!m->pinned)
            {
                //添加重力
                m->forces += gravity*m->mass;
                //计算加速度
                Vector2D a = m->forces/m->mass;
                //记录现在的位置
                Vector2D temp_position = m->position;
                // TODO (Part 4): Add global Verlet damping
                //添加阻尼
                double damp=0.00005;
                m->position = m->position+(1-damp)*(m->position-m->last_position)+a*delta_t*delta_t;
                //m->position = m->position+(m->position-m->last_position)+a*delta_t*delta_t;
                //存放x(t-1)
                m->last_position =temp_position;
            }
            m->forces = Vector2D(0,0);
        }
    }

结果：

1.只添加了构造函数，是一根不动的绳子

2.实现了欧拉方法但没有加摩擦力，蓝色的绳子不断抖动停不下来

3.更改质点数目从3到16（application.cpp)中

 ropeEuler = new Rope(Vector2D(0, 200), Vector2D(-400, 200), 16, config.mass,
                       config.ks, {0});
  ropeVerlet = new Rope(Vector2D(0, 200), Vector2D(-400, 200), 16, config.mass,
                        config.ks, {0});