Физика «на пальцах»: Position-Based подход

Автор: Александр Санников

В этой статье из цикла «Физика на пальцах» речь пойдёт об одном из подходов, с успехом применяемом в самых разных областях математического моделирования: начиная от динамики жидкости и заканчивая физикой взаимодействия твёрдых тел.

Основные принципы Position-Based подхода
В чём же профит?
Further Improvements
Приятная находка
Ограничения метода

Основные принципы Position-Based подхода

Рассмотрим обычный подход, ещё называемый velocity-based, к моделированию, например, движения материальной точки. Точку в таком подходе принято описывать её мгновенным положением, скоростью и ускорением:

struct ClassicPoint
{
  Vector3 position;
  Vector3 velocity;
  Vector3 acceleration;
  void Move(float dt);
};

И если мы хотим предсказать, как изменится положение точки, когда пройдёт интервал времени dt, мы считаем, что на протяжении этого интервала времени частица движется с постоянным ускорением acceleration. Так как dt мал, мы можем посчитать приближённое изменение скорости и положения, численно проинтегрировав по Эйлеру. Это нетрудно себе представить: скорость изменяется на малую величину acceleration * dt, а положение — на velocity * dt.

+ Для интересующихся мат.частью

В коде это можно представить так:

void ClassicPoint::Move(float dt)
{
  position += velocity * dt;
  velocity += acceleration * dt;
}

здесь * — оператор умножения вектора на скаляр.

Отличие же Position-Based подхода от Velocity-Based состоит в том, что мы не храним скорость и ускорение в явном виде. Вместо этого мы храним положение точки на предыдущем кадре и на текущем. Например, так:

struct PositionBasedPoint
{
  Vector3 position;
  Vector3 prevPosition;
  Vector3 acceleration;
  void Move(float dt);
};

А передвигаем точку из интуитивно понятного принципа: если с предыдущего кадра до текущего прошёл малый интервал времени, и за это время частица сместилась на малое расстояние deltaPos, то, скорее всего, за следующий такой же интервал времени частица сместится на такое же расстояние. За объяснением, откуда взялся член с ускорением, придётся открыть спойлер:

+ Строгое обоснование

В коде метод можно выразить так:

void PositionBasedPoint::Move(float dt)
{
  Vector3 delta = position - prevPosition;
  prevPosition = position;
  position += delta + acceleration * dt * dt;
}

Обратим внимание, что для успешной работы метода dt всегда должен быть постоянным, так как мы считаем, что от предыдущего шага интегрирования до текущего прошло столько же времени, сколько от текущего до следующего.

В чём же профит?

Что же нам даёт этот метод? Рассмотрим, например, взаимодействие большого количества круглых частиц друг с другом:

Сперва расширим класс нашей position-based точки, добавив ей радиус, теперь это будет полноценная частица:

struct PositionBasedParticle
{
  Vector3 position;
  Vector3 prevPosition;
  float radius;
  void Move(float dt);
};

В классическом честном velocity-based подходе нам бы пришлось считать время до ближайшего соударения двух шаров, смещать всю систему на это малое расстояние по времени, высчитывать их углы отскока, заниматься неинтересными и вычислительно дорогими вещами. В нашем же подходе мы делаем следующее: просто расталкиваем на равное расстояние шары, которые на текущей итерации пролезли немного друг в друга:

for(int i = 0; i < particles.size() - 1; i++)
{
  for(int j = i + 1; j < particles.size(); j++)
  {
    Vector3 p1 = particles[i].position;
    Vector3 p2 = particles[j].position;
    if((p1 - p2).SquareLength() < sqr(particles[i].radius + particles[j].radius))
      // расстояние между центрами шаров меньше
      // суммы их радиусов -> они пересекаются
    {
      // шары удобнее всего расталкивать вдоль прямой, проходящей через их центры
      Vector3 penetrationDirection = (p2 - p1).Normalize();
      // глубина проникновения - это сумма радиусов двух шаров
      // минус расстояние между их центрами
      float penetrationDepth = particles[i].radius + particles[j].radius
                - (p2 - p1).Length();

      // расталкиваем шары в противоположенных направлениях
      particles[i].position += -penetrationDirection * penetrationDepth * 0.5f;
      particles[j].position +=  penetrationDirection * penetrationDepth * 0.5f;
    }
  }
}

Данный цикл лучше прогнать несколько раз, так как при выталкивании шара A из шара B, он может ещё глубже залезть в шар C, но растолкнув их попарно несколько раз, в пределе все коллизии разрешатся.

Далее, если мы вызовем для всех шаров функцию ::Move, столкнувшиеся шары «магическим» образом разлетятся в нужном направлении с нужными скоростями!

for(int i = 0; i < particles.size(); i++)
{
  paraticles[i].Move(dt);
}

+ Почему это работает

С виду простой метод возволяет моделировать вполне внушительное количество взаимодействующих объектов:

На этом просчитывается взаимодействие 5000 шариков, расчёт всей физики занимает чуть больше 2мс на одном ядре процессора Core2Duo E6600 2.4GHz.

Но и это ещё не всё. Предположим, мы хотим соединить какие-то пары шариков жёсткими связями. Заведём тип данных, описывающий эту самую связь:

struct ParticleLink
{
  float defLen; //длина пружины в ненапряжённом состоянии
  float stiffness; //параметр, задающий жёсткость связи

  PositionBasedParticle *particle1;
  PositionBasedParticle *particle2;

  void Solve();
}

После каждого шага по времени мы варварски передвигаем частицы на расстояние, соответствующее их смещению с предыдушего кадра к текущему. Разумеется, из-за этого связи могут нарушаться, растягиваясь или сжимаясь. Придерживаясь нашего position-based подхода, всё, что от нас требуется — просто сжать или растянуть связь до её собственной длины. Реализуем это в методе Solve():

void ParticleLink::Solve()
{
  Vector2 point1 = particle1->position;
  Vector2 point2 = particle2->position;

  Vector2 norm = (point2 - point1).Normalize();

  Vector2 goal1 = (point1 + point2) * 0.5f - norm * defLen * 0.5f;
  Vector2 goal2 = (point1 + point2) * 0.5f + norm * defLen * 0.5f;

  particle1->position = particle1->position + (goal1 - particle1->position) * stiffness;
  particle2->position = particle2->position + (goal2 - particle2->position) * stiffness;
}

Обратим внимание на параметр stiffness. Если он равен нулю, то после «решения» связи частицы не сдвинутся с места. Если он равен единице, частицы примут своё положение, как если бы связь мгновенно приняла «удобную» для себя длину. Этот параметр означает что-то вроде жёсткости связи. Его максимальное значение — единица, соостветствует максимальной жёсткости, допустимой численным методом.

Грубо говоря, всё, что от нас требуется для реализации какого-то физического явления — это каждый кадр ставить все частицы в положение, удовлетворяющее этому явлению и более-менее близкое к текущему. Если это затруднительно сделать, как, например, в случае с большим количеством связей не тривиально найти положение шариков, когда не нарушена ни одна связь, мы просто выбираем «менее плохую» конфигурацию, постепенно решая, каждую связь отдельно от других. Всё остальное получается как бы само собой — шарики расталкиваются, конструкция из связей получает жёсткость.

Further Improvements

Одной из неприятных особенностей метода является то, что энергия системы при каждом расталкивании пары шаров гарантированно не увеличивается. Более того, нельзя точно сказать, уменьшается ли она. В описанной реализации отскочут ли шары друг от друга, как мячи, или слипнутся, как пластилиновые, сказать заранее нельзя. Благо, в хаосе из тысяч частиц это не заметно, но если мы реализуем, например, бильярд, на самотёк это дело пускать нельзя.

Нужно найти способ, как явно задавать скорость частиц после соударения. Один из подходов — менять prevPosition, но я считаю его антилогичным. Вместо этого я предлагаю поступить следующим образом:

Вместо того, чтобы хранить пару prevPosition + position, можно хранить пару position + deltaPosition, где deltaPosition — это разница между position и prevPosition. Соотвественно, каждый раз изменяя значение position, нам нужно в точности на такую же величину изменить и deltaPosition, в этом нетрудно убедиться взяв листок бумаги, нарисовав предыдущую и текущую позицию точки и посмотрев, как изменяется вектор, их соединяющий, если двигать текущую позицию.

Соответственно изменится и код класса чистицы:

struct PositionBasedParticle
{
  Vector3 position;
  Vector3 deltaPosition;
  float radius;
  void Move(float dt);
  void Push(Vector3 delta);
};

void PositionBasedParticle::Move(float dt)
{
  deltaPosition += acceleration * dt * dt;
  position += deltaPosition;
}

void PositionBasedParticle::Push(Vector3 delta)
{
  position += delta;
  deltaPosition += delta;
}

Теперь, чтобы корректно обработать отскок пары шаров, мы не только расталкиваем их (обратим внимание, теперь уже не простым присвоением текущей позиции, а методом Push()), но и пересчитываем новые скорости из закона сохранения импульса:

Для интересующихся, откуда был взят код пересчёта скоростей после отскока рекомендую обратиться к статье «Физика «на пальцах»: солверы физических движков» по разрешению контактов в импульсных солверах.

Ничто принципиально не мешает добавить в метод угловые скорости, трение и тела более сложной формы. Но, к сожалению, после таких улучшений подход выглядит уже далеко не так лаконично и красиво.

Немного модифицировав этот алгоритм, чтобы связи образовывались и рвались на определённом расстоянии между частицами, можно добиться интересных эффектов, например, моделировать что-то вроде вязкой жидкости:

Если добавить третье измерение, построить сетку из частиц и выправлять слишком «сложившиеся» её ячейки, получим что-то вроде ткани:

Если мы выделим группу частиц и скажем, что объём, который они ограничивают, должен сохраняться (этакий бурдюк с несжимаемой жидкостью), и будем каждый кадр просто увеличивать его (растягивая вдоль координатных осей) во столько раз, чтобы сохранить объём, получим так называемый baloon:

Наконец, если свяжем частицы не отрезками, а тетраэдрами, мы получим объёмные связи. «Выправляя» их от кадра к кадру можно получить нетривиальные явления, например, таким подходом удобно моделировать деформируемые тела:

или распространение ударного волнового фронта:

Приятная находка

Я ещё давно заметил, что если интегрировать движение тела таким образом:

void Particle::Move(float dt)
{
  velocity += acceleration * dt;
  position += velocity * dt;
}

То при константном шаге по времени точность получается существенно выше, чем, если даже использовать больше членов в разложении Тейлора:

void Particle::Move(float dt)
{
  position += velocity * dt + acceleration * dt * dt / 2.0f;
  velocity += acceleration * dt;
}

Объяснение просто - "волшебная" формула при постоянном шаге является не интегратором по Эйлеру первого порядка, а переходит в position-based третьего порядка, в формуле достаточно сделать замену velocity = deltaPosition / dt и она полностью перейдёт в ту, что описана в Further Improvements.

Ограничения метода

Разумеется, чудес не бывает и не бывает одновременно простого, быстрого и общего метода моделирования физики. За простоту реализации и широкий спектр моделируемых явлений position-based подход заплатил своё:

Любые попытки изменять шаг по времени до добра не доведут. Если оригинальная схема имеет третий порядок точности по времени, то использование нефиксированного шага сбрасывает порядок точности до первого. Если в игре не фиксированный fps, лучше вызывать физику не каждый кадр а, например, тогда, когда с момента предыдущего обновления физики прошло какое-то время, например, 1/50 секунды.
Метод вовсе не дружит с continuous collision detection. Если неаккуратно выбрать размеры или скорости тел, они могут безнаказанно пролетать друг через друга. Если избавиться от быстро летящих частиц нельзя, то можно попытаться обрабатывать их другим, не position-based подходом.

#game dynamics, #физика, #position-based

1 марта 2010 (Обновление: 10 июня 2014)

Комментарии [191]

Тир для пневматического оружия.