Impulsive Behavior
점 P에서 충돌하려는 A, B 오브젝트
점 P의 속도는 각 오브젝트 A와 B의 관점에서 :
Relative Velocity(상대속도)
Relative Normal Velocities
A collision occurs when a point on one body touches a point on another body with a negative relative normal velocity.
충돌은 한 body에 있는 한 점이 다른 body에 있는 한 점과 음의 relative normal velocity를 가지고 닿을 때 발생한다.
It means that
Eq.2 > 0 : 그 점들은 서로 떠나는 중이다 -> 무시
Eq.2 = 0 : 그 점들은 충돌하지도 분리하지도 않는다 -> Contact
Eq.2 < 0 : 서로 부딪치고 있음 -> 관통하는 것을 막아야 한다.
force는 즉시 속도를 바꾸지 않으므로 (적분을 해야해서)
impulse라는 새로운 quantity를 도입한다.
즉, Huge force(infinite force) over a really short period time (infinitesimal time)이
곧 finite impulse가 된다.
Force changes the momentum over time (F = dot(p) = m * dot(v) = ma)
Impulse changes the momentum instantaneously, which in turn changes our velocity.
위의 개념을 모델링 하기 위해 다음의 모델을 쓴다
Newton's Law of Restitution for Instantaneous Collision with No Friction
Instantaneous : The collision process takes no time (very small amount of time)
Restitution : Coefficient of Restitution (e or 𝛜)
이 복구 계수로, 충돌하는 bodies의 복잡한 compression(응축)과 restitution(복구)를 단일의 scalar로 모델링한다. 그리고 이것은 Contact Point의 incoming relative normal velocity와 outcoming relative normal velocity를 연관짓는다.
식으로 표현하자면,
좌변이 Outcoming이고, 우변이 Incoming이다.
e인 coefficient of restitution은 우리에게 그 들어오는 에너지가 충돌동안 얼마나 흩어지는지를 말해준다.
e = 1 -> A totally elastic collision (superball)
e = 0 -> A totally plastic collision (a lump of clay landing on the floor)
이 모델은 마찰이 없도록 가정하여 간단하게 만듬
우리는 impulse를 단일 스칼라 j와 normal의 곱으로 표현할 수 있다즉 우리에게 jn을 준다.
A에 의해 느껴지는 impulse는 jn이고,
B에 의해 느껴지는 impulse는 -jn이다.
- 회전하지 않는 오브젝트에 대한 충돌반응 방정식
CM의 incoming/outcoming velocities는 (현재 미지수의) impulse의 영향하에서 연관지어야 한다.
위 처럼 쓸 수 있는 이유는, impulse가 momentum에서의 변화를 의미하기 때문이다
CM의 속도 -> 개별 bodies에 있는 모든 점들의 속도
따라서, Eq. 1에 있는 v^{AP}는 v^P로 바뀐다.
따라서 이러한 것들을 통해 우리는 j를 풀 수 있다.
CM에 의해서,
그리고 이것을 Eq.3과 Eq.4를 이용하여
j를 구했다.
이 식에서, n이 단위벡터일 필요가 없다. dot(n, n)이 상쇄시켜기 때문이다.
그리고 n이 단위벡터라면, 분모의 dot(n, n) 항을 그냥 없앨 수 있다.
Eq 6은 빌딩이나 땅 같은 고정되는 rigidbody와 움직이는 rigid body 사이의 충돌을 처리할 수 있다. 왜냐하면, 오브젝트 중의 하나의 무게가 증가할 때, impulse의 효과는 줄어든다. 만약 무게가 무한이라면, 그 방정식은 그 오브젝트의 질량을 고려하지 않고, 고정된 오브젝트와의 충돌방정식으로 degeneration이 된다.
A의 mass가 1이고, B의 mass가 무한이고 속도가 0이고, coefficient of restitution이 1이라면, v^A를 normal에 대해반사시키는 것과 같다.
- 회전하고 평행이동하는 오브젝트에 대한 충돌반응 방정식
회전하고 평행이동하는 rigid body의 임의의 점에서의 velocity는
그리고 위에서 처럼 incoming/outcoming velocity를 연관짓는다면,
Eq.8b는 body A에 있는 점 P에 impulse jn을 적용한 결과이다.
이것은 linear force를 torque로 바꾸듯이,
linear impulse를 angular impulse로 변환한것ㅅ이다. (한 점에 perp-dot product를 사용해서, 3D에서는 cross 연산을사용한다).
관성모멘트 I^A로 나눠주는 것은 impulse가 angular momentum을 바꾸기 때문에, 관성모멘트로로 나누어준다. linear 의 케이스와 같이.
따라서 위의 식은 Collision impulse가 어떻게 A의 충돌 전 속도에 영향을 미치는지를 보여준다.
이제 B는 j항에 있는 곳에 -만 붙여준다는 것만 알게 된다면, 이전 케이스처럼 식을 변형하여 풀면 된다.
난 여기서 j를 구하려고 풀 때 방정식이 어떻게 푸는지 이해가 안됐었는데,
CM에 대해 우리는 풀게 되고, j를 구하려고 변환하는 와중에, angular velocity가 남아있는데 Eq.7을 통해서 다시 v^A_1 - v^B_1로 바뀐다는 것을 염두에 두어야 한다.
어쨋든 j에대해 풀자면
이 이후의 것은 물리엔진 관련인데, 그냥 바로 번역하도록 하겠다.
A Little Touch Up
이제 너가 충돌 반응 방정식을 알았으니, 그것들이 우리의 전체 시뮬레이션 루프에서 어떻게 맞아들어가는지 봐보자. Listing1은 샘플 프로그램으로부터의 충돌 탐지와 반응을 지원하는 simulation loop에 대한 슏코드를 보여준다. 나는 슈도코드에 대해 마지막 문제의 step-by-step 알고리즘을 바꾸었다. 왜냐하면 그 loop는collisions을 다루기 위해 확장될 떄 조금 더 복잡해지기 때문이다.
이 새로운 복잡성의 근원은 충돌의 "정확한" 시간을 계산하는 것이다. 우리가 처음에 full time step으로 forward로 integrate하는 것을 알아채라. 그리고 만약 새로운 configuration에서 interpenetration이 있다면, 우리는 그 time interval을 subdivide하고, 다시 시도한다. 그 알고리즘은 충돌시간을 찾는 시간에 대한 binary search를 하는 것으로 종합된다. 이것은 필수적으로 충돌 시간을 찾는 가장 효율적인 방법이 아니다. 왜냐하면 우리는 우리의 모든 이전 integration work를 버리기 때문이다. 그러나 그것은 간단하고 튼튼하다. 그 문제에 대한 다른 솔루션들은 충돌이 언제 발생했는지를 측정하는데 돕기 위해서, 이전의 integration parameters를 사용하는 것을 포함한다. 그리고 충돌이 발생할 시간을 앞서서 예측하려고 하거나, 심지어 interpenetrating coordinates를 사용하고 그것이 이상하게 보이지 않기를 희망한다. 또한 이 discrete collision routine은 "tunneling"을 잡지 않는다, tunneling에서, 빠르게 움직이는 오브젝트들은 완전히 단일 integration step에서 다른 오브젝트들을 뚫고 지나갈 수 있다.
noninterpenetrating configuration이 발견된다면, 우리는 그 충돌을 해결한다 - 만약 존재한다면 - 그리고 그 configuration을 업데이트한다. 그러고나서, 우리는time step을 마무리하기 위해 loop back up하고, 최종적으로 그 오브젝트들을그린다.
나는 이 글에서 몇 가지의 것들을 보여줬다, 그래서 나머지 공간을 Bond-O를 get out하기 위하 나머지 공간을 쓰고, 그 구멍들을 채우자.
충돌 반응 방정식으로부터, 우리가 한 충돌에 대해 4개의 데이터를 알 필요가 있다는 것이 명확해져야 한다: 충돌 시간, 충돌에 참여하는 오브젝트들, 그러한 오브젝트들에 있는 충돌점들, 그리고 충돌 normal. 이러한 파라미터들의 각각은 어떤 subtleties를 가지고, 우리는 차례로 각각 알아볼 것이다.
너는 내가 data의 요구되는 첫 번째의 것을 언급할 때 몇 문장에서 "정확한"이라는 단어를 인용한 것을 눈치챌 것이다: the collision time. 그 이유는 너가 컴퓨터에서 numerically하게 작업할 때 exact collision time과 같은 것은 없기 때문이다. 우리는 충돌 탐지에 대해 tolerance value를 사용하도록 강요되는데, 거기에서, 우리는 우리가 충돌한다고 말하는 것을 동의한다 (관통하거나 닿지않다고 하기보다). 그 샘플코드는 이 기법을 보여준다.
다음 데이터는 - 충돌하는 오브젝트들의 collection - 명백한 것 처럼 보이지만, 우리의 현재 알고리즘이 두 bodies 사이의 단일 충돌만을 다룰 수 있는 것을 주목해라. 유사한 한계가 세 번 째 파라미터인 collision points에 대해서도 유효하다. convex polygons 사이의 2D 충돌에서, 너는 edge/edge collision을 얻을 수 있다는 것을 보는 것은 쉽다. 그리고 그 edge/edge collision은 "collision "manifold"가 더 이상 한 점이 아닌 한 선분을 의미한다 - manifold는 닿고있는 오브젝트들의 부분을 나타내는 공간이다. 너는 이 종류의 충돌에 대해 선분의 정점들을 사용하여 할 수 있지만, 심지어 그것은 우리의 현재 충돌 반응 routine의 powers를 넘어선 것이다. 그것은 오직 단일의 충돌 점만을다룰 수 있다, 다양한 동시의 충돌점들이 아닌. Simulataneous collisions는 더 어렵고, 다른 시간을 기다려야 할 것이다. 3D에서상황은 더 악화된다, 거기에서 너는 point, edge, and face collision을 convex polyhedrons과 얻을 수 있고, 충돌 탐지와 반응은 너가 곡면과 관련될 때 악몽이 된다. 어쨋든, 그 샘플 프로그램의 충돌 탐지는 현재 단일의 충돌점을 반환하고, 비록 우리가 flat-edged bounces를 얻지 못할지라도 (그것은 항상 한 점이 처음 닿는 것처럼 보인다), 그것은 여전히 꽤 좋아보인다.
그 충돌 normal은 애매함이 발생하는 마지막 장소이다. 2D에서, edge/edge or a vertex/edge collision에서, 그 normal vector는 쉽게 edge에 수직인 벡터로서 얻을 수 있는 것이다. 그러나, vertex/vertex collision에서, 너는 collision normal로 사용하기 위해 sensible vector를 고를 필요가 있다. 그 샘플 프로그램은, vertex/vertex 충돌을 vertex/edge collisions로 다루어 이 문제를 피하지만, 비현실적인 행동을 이끌 수 있다.
이 주제에 대한 레퍼런스들은 전체 칼럼의 공간을 채울 것이다, 다시 한 번. 나는 그것들을 내 웹사이트에 올려놓을 것이다. 내가 여기에서 사용한 유도들은 David Baraff의 physically based modeling에 대한 SIGGRAPH tutoriap에서 방정식과 유사하다 (그것은 내 레퍼런스들에 있다). 수학과 물리학에서 대부분의 결과처럼, 그 같은 방정식에 도달하는 많은 방법들이 있다, 에너지와 운동량의 보존 법칙을 기반으로한 유도와, "generalized coordinates"라고 불리는 것에 기반한 유도들을 포함해서. 만약 너가 이 주제를 진지하게 공부한다면, 너는 너가 그것들을 이해하도록 하기 위해 많은 다른 방식들로 그 방정식들을 처리하길 원할 것이다. 너가 더 연습할스록, 너는 더 좋은 수학자이자 물리학자가 될 것이다. ~
Chapter 14 Collision Resolution
나는 이 책의 이 부분의 설명이 그렇게 잘 이해가 되지 않았다. 왜 그러는가 하고 살펴보니, 여러 이론에 대해 순차적인 발달이 없이, 여러 이론을 가져다가 잠깐 설명하고, 나중에 다시 그 설명된 것을 바로 가져와서 쓰니, 이해가 연결성이 안되어 이해하기 어려웠다.
그래서 나는 코드의 실행 순서에 따라 이 챕터를 다시 스스로 설명하려 한다.
Collision Resolution 순서는
prepareContacts()
adjustPositions()
adjustVelocities()
이다. 위의 함수들이 순서대로 호출되어 실행된다.
그런데 책에서는 prepareContacts() <-> adjustVelocities() 부분을 왔다갔다하며 설명하고, 마지막으로 adjsutPositions()를 설명한다.
나는 일단 이 챕터가 길고 복잡하니, prepareContacts() -> adjustVelocities() 에 대한 순서로 먼저 설명하고, adjustPositions()를 설명할 것이다. 왜냐하면, prepareContacts()와 adjustVelocities()에 대한 이해가 완벽해야, adjustPositions()가 하려는 것을 이해할 수 있기 때문이다.
- prepareContacts()
이 함수는 contact->calculateInternals() 를 실행한다.
calculateInternals()의 함수안에는 크게 4가지의 작업이 실행된다. 4가지의 작업에 대해서 간단히 설명한 후, 또 다시 각각 세부적으로 볼 것이다.
1. caculateContactBasis() : contact coordinates system(contact 좌표계) 계산. 이것은 좌표변환을 통해 수학을 좀 더 간단하게 하여 연산을 좀 더 쉽게하려고 하는 것이다.
2. calculate relativeContactPosition : 이것은 함수 호출이 아니라 바로 계산이 되는 것인데, body의 CM으로부터 contactPoint까지 가는 각각의 relativeContactPosition을 구한다. 이것을 하는 이유는, Chris Hecker Article 3에서 나와있듯이, 나중의 angular component를 위해서이다.
3. calculateLocalVelocity() : Chris Hecker Article에 나와있듯이, 접촉점의 속도를 구하려고 하는 것이다. 여기에서 LocalVelocity인 이유는 월드 좌표계가 아닌, contact space에서 구하려고 하기 때문이 LocalVelocity가 된다.
4. calculateDesiredDeltaVelocity() : Chris Hecker Article 3에 따라,
이 j를 구하려고 쓰는 것이다. 위의 Eq.9는 2차원에 대해 j값, 즉 impulse quantity를 구하는 것인데, David Baraff 논문에 따라서 3차원에서 저 식은 다음으로 변한다.
이제 다시 각각의 함수에 대해서, 다시 코드 라인별로 살피고 이론을 살펴보겠다.
1. calculateContactBasis()
contact coordinates system을 구성하기 위해, GPED의 글쓴이는 contactNormal를 X축으로 삼아 그에 따라 좌표계를 구성하려 한다. 그 프로세스는 다음과 같다.
X-axis = contactNormal
World-Yaxis = (0, 1, 0)
Z-axis = Cross(X-axis, World-Yaxis)
Y-axis = Cross(Z-axis, X-axis)
하지만, X-axis가 World-Yaxis와 같은 방향을 가리키고 있다면? 그러면 Z-axis, Y-axis 둘 다 ZeroVector가 되어버려서 올바른 좌표계를 구성할 수가 없다. 그래서 글쓴이는 이를 해결하기 위해, 다음과 같은 loose(느슨한) 추정을 한다.
real_abs(contactNormal.x) > real_abs(contactNormal.y)
이면, contactNormal은 X-axis로 구성되기에 적당하고, World-Yaxis를 가리키고 있지 않다. 왜냐하면 (0, 1, 0) 에서 x component < y component이기 때문이다. 그래서
if(real_abs(contactNormal.x) > real_abs(contactNormal.y))
{
이 경우에는 contactNormal이 contact space에서 X축으로 쓰여진다.
글쓴이는 그런데 여기서 성능 Optimization을 위해서, 다음의 사실을 이용한다.
Z-axis = Cross(X-axis, World-Yaxis)
= (-X-axis.z, 0, X-axis.x)
Z-axis가 저렇게 구성되는데, 좌표계를 구성하는 basis는 표준화되어야 한다는 사실에서, 우리는
const real s = (real)1.0f / real_sqrt(contactNormal.z*contactNormal.z +
contactNormal.x*contactNormal.x);
이것을 계산해서 표준화 하는 것을 좀 더 간단하게 한다. 즉, 정확히 Z-axis는
Z-axis = (-X-axis.z * s, 0, X-axis.x * s);
이렇게 되어 Z-axis도 표준화 된다. 따라서 그 이후의 Y-axis 에 대한 계산은 X-axis와 Z- axis가 이미 표준화 되어있으므로, Y-axis를 표준화 할 필요가 없게 된다.
}
else
{
이 경우에 contactNormal가 World Yaxis(0, 1, 0) 방향을 가리키고 있다고 가정하여,
우리는 이제 다음의 프로세스를 거친다.
Y-axis = contactNormal
World-Xaxis = (1, 0, 0)
Z-axis = Cross(World-Xaxis, Y-axis)
= (0, -Y-axis.z, Y-axis.y)
X-axis = Cross(Y-axis, Z-axis)
그리고 위에도 설명했듯이, 성능 최적화를 위해, const real s를 계산해주어 표준화하는 작업을 줄인다.
}
이제 계산된 contact basis들을 가지고 contactToWorld Matrix를 구성하는데,
contactNormal가 항상 1번째 열 contactTangent[0]가 2번째 열 contactTangent[1]이 3번째 열에 들어가도록 항상 한다.
contactNormal이 Y-axis를 가리키고 있을 상황에도 1번째 열에 넣는 이유는, 우리가 나중에 contact 좌표계에서 x-axis component가 contact normal과 관련된 것을 알고 항상 그 컴포넌트를 쓰게 하기 위해서 이다. 그리고 X-axis, Y-axis, Z-axis가 if문에 따라 상황이 달라져도, 저렇게 똑같이 항상 넣어도 되는 이유는, 어쨋든 3개의 직교하고 표준화된 좌표변환 행렬을 구성하게 되면, 우리는 거기서 어떻게든 작업해도 상관없게 된다. 그래서
// Make a matrix from the three vectors
contactToWorld[0] = contactNormal;
contactToWorld[1] = contactTangent[0];
contactToWorld[2] = contactTangent[1];
이 위의 코드가 항상 제대로 작동하게 된다.
2. calculate relativeContactPosition
이것은 아주 간단한데, 이것을 왜 구하는지 모르겠으면, Chris Hecker Article 3를 반드시 읽고 이해해야 한다. Chris Hecker Article 2를 읽고 이해해야 하는 것은 당연히 선행되어야 한다. 왜냐하면 angular motion을 simulate하기 위해서, 한 body의 CM에서 contact point를 향하는 relativeContactPosition이 필요하기 때문이다. 그래서 코드가 다음 처럼된다.
// Store te relative position of the contact relative to each body
relativeContactPosition[0] = contactPoint - body[0]->getPosition();
if (body[1])
{
relativeContactPosition[1] = contactPoint - body[1]->getPosition();
}
3. calculateLocalVelocity()
아까 위에서 잠깐 설명했듯이, 결국에는 j값, 즉 impulse quantity를 구하기 위해서 이것을 계산한다.
즉, j를 구하는데, 저기 위의 분자의 v_rel을 구하기 위해서 한다. 이 때 dot(p)_a와 dot(p)_b 가 calculateLocalVelocity() 함수가 구하는 부분이다.
그러면 어떤 것의 velocity를 구하느냐인데, Chris Hecker Article 3를 읽으면 알 수 있다. 즉, 각 Body의 기준으로 Contact Point의 속도이다. Contact Point라는 것은 각 Body의 접촉점이기 때문에, 각 body는 그 contactpoint를 알 수 있다. 그 방정식의 유도는 Chris Hecker Articl 2에 자세히 설명되어있다. 어쨋든 그 방정식은
이렇게 된다. 사실 저 식은 2D에서 구성되는데, 3D에서 cross를 이용하면 된다. 저 식을 잘 모르겠으면 Chris Hecker Article 2를 잘 읽고 이해해야 한다. 왜냐하면 저것을 설명하기 위해 글이 엄청 길어지기 때문이다.
어 쨋든 코드로는,
// Work out the velocity of the contact point
glm::vec3 velocity = glm::cross(thisBody->getRotation(),
relativeContactPosition[bodyIndex]);
velocity += thisBody->getVelocity();
이렇게 하면 끝이다. 근데 우리는 계산을 편하게 하기 위해, Contact Coordinates System으로 바꿀려고 한다.
// Turn the velocity into contact-coordinates.
glm::mat3 WorldToContact = glm::transpose(contactToWorld);
glm::vec3 contactVelocity = WorldToContact * velocity;
그래서 이렇게 하여 contact space에서 각 body의 contact point의 local velocity를 구하게 된다.
그래서 첫 번째 오브젝트에 대해서 구했고, 또 두 번째 오브젝트가 존재한다면,
// Find the relative velocity ofthe bodiesat the contact point.
contactVelocity = calculateLocalVelocity(0, duration);
if (body[1])
{
contactVelocity -= calculateLocalVelocity(1, duration);
}
이렇게 구성하여 빼준다. 기억이 난다면 알테지만,
우리는 이렇게 하여 (dot(p)_a - dot(p)_b)를 구한 것이다. 결국에는 j를 구하기 때문에. 여기서 코드에서 내적이 안되어 있는것 처럼 보이는데, 우리는 이미 contact space에 있기 때문에, 아까 transpose한 거에서 저 처리가 이미 되어있다. 그래서 우리는 v_rel을 이미 구하게 된 것이나 다름없다.
4. calculateDesiredDeltaVelocity()
이제 j의 분자를 완성할 시간이다.
desiredDeltaVelocity = -(1 + thisRestitution) * contactVelocity.x
contactVelocity는 contact space안에 있기 때문에, 그리고 우리가 X축이 contact normal축이 되도록 설정했기 때문에 j의 분자가 위의 코드로 완서이 된다.
여기에서 이제 우리는 prepareContacts() 함수가 마무리 된다.
- adjustVelocities()
이제 adjustVelocities() 를 설명할 것이다. 이 함수안을 보면 뭔가 복잡해 보이겠지만, 어쨋든 핵심은 applyVelocityChange() 함수만 제대로 이해하면 된다.
applyVelocityChange()는 Chris Hecker Article 3를 제대로 공부해서 이해했다면 쉬운 코드이다.
Chris Hecker Article 3는 2차원에서 하기 때문에, David Baraff를 이용해서 이해에 도움이 되고자 3차원 식을 가져와보자.
우리가 하고 싶은 것은 위에서 j를 구한 후, 그것을 이용하여 8-12 식에서 처럼, linear velocity와 angular velocity에 변화를 주어서 collision resolution을 하고 싶은 것이다.
이제 구현된 코드를 보자. applyVelocityChange() 함수에서중요한 부분은,
1. calculateFrictionlessimpulse() // 나중에 FrictionImpulse는 다음 챕터에서 설명된다.
2. 변화 적용
위의 두 가지가 된다. 뭔가 복잡해보이지만, 저 수식을 통해 간단히 설명된다.
1. calculateFrictionlessImpulse() : 우리는 prepareContacts()함수에서 j의 분자값을 desiredDeltaVelocity라는 변수에 넣어서 구했다. 그러니까 이제 그 분모를 구해서 j값을 구하기 위해 이 함수를 실행시킨다. 코드를 보면 저 식의 분모식과 똑같이 되어있다는 것을 알 수 있다.
나는 여기에서 GPED 저자가 이것을 설명하는 방식을 이해하려 조금 애썼는데, 저자는 이렇게 설명했다. 이 calculateFirctionlessImpulse에서 분모를 구하는 일을 저자는
Velocity Change Per Unit Impulse 라고 설명했다.
즉, body의
velocity change per unit impulse는 linear component와 angular component에 의해서 구성이 된다. 그래서 linear component는 단순히 inverse mass가 되고,
angular component 다음의 과정으로 구해진다
u = cross(q_rel, contactNormal) : 단위 impulse로 생성된 impulsive torque.
Δdot(θ) = u * inverse_inertia_tensor : 단위impulsive torque에 대한 angular velocity의 변화. 즉 나의 이해대로 하자면, impulsive torque의 힘이 얼마만큼의 변화량을 가져야 하는지.
angular component = glm::cross(Δdot(θ), q_rel) : 이제 그 변화량에 대해, q_rel을 cross하는 것은
이 식에서 linear을뺀 angular velocity즉, 한 점의 angular velocity를 갖게된다. 그래서 이것을 angular component로 보는 것이다.
그래서 우리는 이
linear component + angular component = velocity change per unit impulse
를 얻게 되는 것이고 이것의 j의 분모가 된다.
j의 분자는 우리가 이전에 구했는데, desiredDeltaVelocity의 개념이다.
따라서 우리가 필요한 impulse = desiredDeltaVelocity / velocity change per unit impulse
이렇게 되는 것이다.
그런데 코드에서는 우리가 contact space에서 계산해야하기 때문에, FrictionlessImpulse에서는 x component에만 그 값을 넣어준다.
2. 변화 적용
이제 요구되는 impulse를 계산했으니, 그 contact space에 있는 impulse를 world space로 바꾸고 변화를 적용해야 한다.
glm::vec3 impulse = contactToWorld * impulseContact;
간단히 이렇게해서 월드 좌표계로 바꾼다.
그리고
이 식에 따라 (Chris Hecker Article 3에 왜 식이 이렇게 되어야 하는지 설명한다), applyVelocityChanges()를 마무리 해준다.
// Split in the impulse into linear and rotation components
glm::vec3 impulsiveTorque = glm::cross(relativeContactPosition[0], impulse);
rotationChange[0] = inverseInertiaTensor[0] * impulsiveTorque;
velocityChange[0] = impulse * body[0]->getInverseMass();
// Apply the changes
body[0]->addVelocity(velocityChange[0]);
body[0]->addRotation(rotationChange[0]);
여기에서 우리는 contact space에서 놀면서 이미 impulse에 필요한 n을 곱해주었기 때문에, linear change를 구할 때 그냥 impulse에 inverseMass만 곱해주기만 하면된다.
angular change도 마찬가지이다.
따라서 이렇게 하여 body의 velocity가 바꼈으므로, linear velocity change와 angular velocity change가 변화된 것을 contact space로 바꾸어서 contactVelocity를 업데이트 해준다. 그리고 contactVelocity도 바꼇으므로, 다시 calculateDesiredDeltaVelocity도 다시 구해준다.
이제 모든 Contact Resoluiton이 마쳐진다.
- adjustPositions()
이제 중간 단계인 adjustPositions()를 설명할 차례이다. 여기 설명은 adjustVelocities()와 비슷하다. 그렇기 때문에 이전의 단계를 이해하지 못헀다면, 다시 읽어서 이해할 필요가 있다. adjustVelocities()와 유사하게, 여기도 하나의 주된 함수인 applyPositionChange() 함수를 이해하는 것이 중요하다.
applyPositionChange() 함수는 크게 세 부분으로 이루어져 있는데,
1. Contact Resolution의 j의 분모를 계산하는 부분, 그런데 linear과 angular component부분으로 나눈다.
2. j를 분모값의 각 component를 이용하여, 관통값 penetration을 전체 component에 대한 비율로 linear and angular move 값을 구한다.
2-1 이 때 너무 많이 회전하는 것을 방지하기 위해 special case code 추가
3. 그래서 구해진 linear move와 angular move양을 통해 직접 position과 orientation에 변화를 준다.
이 프로세스를 깊이 들어가기전에, 책에서 설명하는 이론에 대해서 설명하도록 하겠다. interpenetration을 해결하는 네 가지 기법이 있는데 그것들은
Linear Projection
Velocity-Based Resolution
Nonlinear Projection *
Relaxation
위와 같은 것들이 있다. 여기에서 글쓴이는 Nonlinear Projection을 선택했는데, 좀 더 근사가 잘되는 기법이기 때문이다. Nonlinear Projection은 Linear Projection을 기반으로 구현된 것이다.
Linear Projection은 Contact Normal 방향으로 objects가 더 이상 닿지 않을 만큼 가장 작게 이동 시킨다. 충돌 관련 각 오브젝트의 움직이는 야은 inverse mass에 비례한다. 즉 가벼운 오브젝트가 무거운 오브젝트 보다 더 많이 움직인다. 이것은 spherical objects 에게는 매우 빠르고 유용하지만, 나머지 것들은 비현실 적이다.
Nonlinear Projection은 이 linear projection에 angular component를 추가해서, 방향도 바꾸게 하는 것이다. 이론은 같은데, contact normal의 방향으로 더 이상 교차하지 않을 때 까지 오브젝트를 움직이는데, 그 움직임이 linear and angular 둘 다를 가지게 된다.
여기에서 penetration depth가 total amount of movement가 된다. 여기에서 글쓴이는 그 penetration depth를 물리학에서 오브젝트들이 부딪힐 때 발생하는 deformation & restitution의 이론에 기반하여 모델링한다.
penetration depth를 물체가 충돌할 때 deformation 되는 양으로 본다 ->
deformation은 오브젝트를 밀어내는 force를 발생시킨다 ->
그 force는 D'Alembert의 원칙에 따라 linear and angular effects를 발생시킨다 ->
그 각각의 effect들은 각 오브젝트의 inverseMass와 inverse inertia tensor에 의존한다 (즉, contact resolution process의 j값의 분모에서 linear와 angular component를 각각 구한다는 것이다, 이것은 우리가 applyVelocitiesChange() 메소드에서 썼던 것과 같은 수학을 쓰게 해준다) ->
즉, 이것은 접촉점에서 impulse or force가 적용될 때, 그 오브젝트가 회전되는 것을 (angular) 얼마나 거부하냐 또 이동(linear)되는 것을 얼마나 거부하는지가 된다. (즉,
linear inertia는 inverseMass에 의해서 기여되고, angular inertia는
cross(q_rel, normal) (impulsive torque) -> result * InverseInertiaTensor (angular velocity change)-> cross(result, q_rel) (contact position의 only angualr velocity) 의 구해진 값에 contactNormal과 내적하여 구해진다)
위에서 잠깐 설명한 것이 1번에 대한 것인데 코드를 보면 이제 이해되기 시작할 것이다.
이제 코드 기반으로 설명한다.
1. Contact Resolution의 j의 linear and angular component 계산 (즉, inertia)
// We need to work out the inertia of each object in the direction // of the ocntact normal, due to angular inertia only. for (unsigned i = 0; i < 2; ++i) if (body[i]) { glm::mat3 inverseInertiaTensor; body[i]->getInverseInertiaTensorWorld(&inverseInertiaTensor); // Use the same procedure as for calculating frictionless // velocity change to work out the angular inertia. glm::vec3 angularInertiaWorld = glm::cross(relativeContactPosition[i], contactNormal); angularInertiaWorld = inverseInertiaTensor * angularInertiaWorld; angularInertiaWorld = glm::cross(angularInertiaWorld, relativeContactPosition[i]); angularInertia[i] = glm::dot(angularInertiaWorld, contactNormal); // The linear component is simply the inverse mass linearInertia[i] = body[i]->getInverseMass(); // keep track of the total inertia from all components totalInertia += linearInertia[i] + angularInertia[i]; // We break the loop here so that the totalInertia value is // completely calculated (by both iterations) before continuing. }
코드에 보이듯이, 먼저 angular component를 구한다.
glm::cross(angularInertiaWorld, relativeContactPosition[i]); 는
위에서 설명했듯이, impulsive torque에 대해, inverseInertiaTensor를 구하여 angular velocity change per unit impulse를 구하고, 그것을 contact point의 속도를 구하는데 쓴다. 그래서 contact Normal 방향으로 내적을 하여 angularInertia를 구한다.
linear Inertia는 그냥 inverseMass만 더한다. 우리는 나중에 penetration에 대해 movement를 설정하기 위해 totalInertia에 둘 다를 더해준다.
2. j를 분모값의 각 component를 이용하여, 관통값 penetration을 전체 component에 대한 비율로 linear and angular move 값을 구한다.
// Loop through again calculating and applying the changes inverseTotalInertia = real(1) / totalInertia; for (unsigned i = 0; i < 2; ++i) if (body[i]) { // The linear and angular movements required are in proportion to // the two inverse inertias. real sign = (i == 0) ? 1 : -1; angularMove[i] = sign * penetration * (angularInertia[i] * inverseTotalInertia); linearMove[i] = sign * penetration * (linearInertia[i] * inverseTotalInertia);
2-1 이 때 너무 많이 회전하는 것을 방지하기 위해 special case code 추가
/** Chan Comment * Refer to the section "Avoiding the large rotation" */ // To avoid angular projections that are too great (when mass is large // but inertia tensor is small) limit the angular move. glm::vec3 projection = relativeContactPosition[i]; projection += contactNormal * -glm::dot(relativeContactPosition[i], contactNormal); // use the small angle approximation for the sine of the angle (i.e. // the magnitude would be sine(angularLimit) * projection.magnitude // but we approximate sign(angularLimit) to angular Limit). real maxMagnitude = angularLimit * glm::length(projection); if (angularMove[i] < -maxMagnitude) { real totalmove = angularMove[i] + linearMove[i]; angularMove[i] = -maxMagnitude; linearMove[i] = totalmove - angularMove[i]; } else if (angularMove[i] > maxMagnitude) { real totalmove = angularMove[i] + linearMove[i]; angularMove[i] = maxMagnitude; linearMove[i] = totalmove - angularMove[i]; }
이것은 책에도 정확히 설명은 안 되어 있는데, special case code이다. 글쓴이는 회전이 너무 많이 되는 것을 제한하기 위해, const angularLimit = 0.2 를 통해 저렇게 구하는데. 저 수식에 대해서 일단 말해보자.
먼저 projection에 CM에서 contact point로 가는 벡터를 projection에 넣어준다. 그리고,
-glm::dot(relativeContactPosition[i], contactNormal);
이 식을 통해, contactNormal 방향으로 relativeContactPosition을 사영시킨다. 사영시킨다는 의미는, relativeContactPosition이 contactNormal의 방향에 떨어졌을 때 되는 길이이다. 이제 그 길이에다 -1을 곱해, contactNormal의 반대 방향만큼 projection에서 뺀다. 그렇다면 projection은 어떤 값을 가지냐면
(CM에서 contact position을 가는 벡터) - (그 벡터를 contactNormal에 사영시키고, 그 사영시킨 결과로 contact Normal와 반대 방향을 가지고 그 사영된 길이 만큼인 벡터)
따라서 projection은 줄어들게 될 것이다. contactNormal의 반대방향과 관련하여, 그러니까, projection이 더 길어진다고 보면된다.
거기에 그것의 길이를 구하고 angularLimit을 곱해준다. 그래서 최소 최대를 비교해서, angularMove를 제한해준다.
이렇게 하는 이유는 글쓴이에 따르면, angular move를 오브젝트의 크기로 스케일링하려고 한다는 것이다. 왜냐하면 오브젝트의 크기에 따라 angular move가 변경된다. 더 큰 오브젝트들은 더 큰 angular movement를 가지게 된다.
3. 그래서 구해진 linear move와 angular move양을 통해 직접 position과 orientation에 변화를 준다.
// We have the linear amount of movement required by turning // the rigid body (in angularMove[i]). We now need to // calculate the desired rotation to achieve that if (angularMove[i] == 0) { // Easy case - no angular movement means no rotation. angularChange[i] = glm::vec3(real(0)); } else { // Work out the direction we'd like to rotate in glm::vec3 targetAngularDirection = glm::cross(relativeContactPosition[i], contactNormal); glm::mat3 inverseInertiaTensor; body[i]->getInverseInertiaTensorWorld(&inverseInertiaTensor); // Work out the direction we'd need to rotate to achieve that angularChange[i] = inverseInertiaTensor * targetAngularDirection * (angularMove[i] / angularInertia[i]); } // Velocity change is easier - it isjsut the linear movement // along the contact normal. linearChange[i] = contactNormal * linearMove[i]; // Now we can start to apply the values we've calculated // Apply the linear movement glm::vec3 pos; body[i]->getPosition(&pos); pos += linearChange[i]; body[i]->setPosition(pos); // And the change in orientation glm::quat q; body[i]->getOrientation(&q); glm::quat deltaOrien(0, angularChange[i]); deltaOrien *= q; q.w += deltaOrien.w * ((real)0.5); q.x += deltaOrien.x * ((real)0.5); q.y += deltaOrien.y * ((real)0.5); q.z += deltaOrien.z * ((real)0.5); body[i]->setOrientation(q); // We need to calculate the derived data for any body that is // asleep, so that the changes are reflected in the object's // data. Otherwise the resolution will not change the position // of the object, and the next collision detection round will // have the same penetration. if (!body[i]->getAwake()) body[i]->calculateDerivedData();
angular move가 0이라면 angular change도 0이다.
angular move가 0이 아니라면,
// Work out the direction we'd like to rotate in
glm::vec3 targetAngularDirection =
glm::cross(relativeContactPosition[i], contactNormal);
이 코드를 통해 회전할 방향을 구한다. 즉, impulsive torque인 것이다.
inverseInertiaTensor * targetAngularDirection
그리고 이것을 통해 얼마만큼 할지, angular velocity change를 구한다. 이것은 unit change인 것이다. 그래서 (angularMove[i] / angularInertia[i])를 곱해서, 움직여야할 양과, 움직이는 것에 대한 저항양으로 나누어 scale해준다.
// Work out the direction we'd need to rotate to achieve that
angularChange[i] = inverseInertiaTensor *
targetAngularDirection * (angularMove[i] / angularInertia[i]);
linear change는 간단하다
linearChange[i] = contactNormal * linearMove[i];
따라서, 이에 따라 pos 업데이트를 해준다.
angular velocity의 change에 따른 orientation 변화는 쿼터니언 적분을 이용해준다.
그리고 깨어있다면, body의 파생데이터를 다시 계산해준다.
- Collision Resolution Process
위에서 Collision Resolution Process를 거의 다 설명했다. 이 부분은 책을 읽으면서 훑으면서 가도 될듯하다.
여기에서 한 가지 성능 개선을 위해 팁을 준게 있는데, adjustPositions와 adjustVelocities에서 contacts의 모든 loop를 돌아서 최대값을 찾아서 알고리즘을 작동하는데, 더블 링크드리스트로 contacts를 관리하여, 그것을 정렬한 후, 알고리즘을 적용시켜서 변경된 것들 body들을 빠르게 업데이트 시키는 것이 있다. 이러면 더 빠르게 할 수 있다.
물리엔진의 성능 최적화를 하기 위해서 내 생각에는 collision detection의 성능을 올리기 위해 spatial data, 공간 분할 알고리즘을 적용해야 하고, contact resolution의 성능을 올리기 위해서 이 contacts의 자료구조를 나중에 손봐야 할 것이다.
Chapter 15 Resting Contacts and Friction
글쓴이는 Resting Contacts와 Friction을 앞에서 구현한 impulse based resolution에 적용하기 위해 모델링하는 과정을 보여준다. 그런데 내가 보기엔 그런 것들이 직관적으로 그리고 단번에 이해하기에 많이 복잡하고 어려웠다. 어쩃든 조금은 이해하긴 했는데 그거에 대해 이야기 하도록 해야겠다. 모든 내용들을 적기보다는 내가 이해한 것들을 기반으로 이것들의 이론과, 우리의 엔진의 모델에 기반으로 하여, resting contacts와 friction을 어떻게 모델링하고 그것을 어떻게 코드에서 구할지에 대해 설명할 것 같다.
그러니까 각 주제에 대한 보통의 순서는
문제제기 -> 문제에 필요한 이론 -> 구현된 엔진에 입각한 모델링 -> 코드 구현
- Resting Contacts
글쓴이는 contacts에 대한 정의를 다시 한다.
Contact : 오브젝트들이 닿거나 교차하는 어떤 location
Collision : Contact의 한 종류, 오브젝트들이 속도를 가지고 함께 움직이는 contact.
Resting Contact : 관련 오브젝트들이 떨어지거나 함께 움직이지 않는 한 contact
여러가지 설명이 있었지만, 가장 잘 이해가 되었던 것을 평면에 박스가 올려진 그림이였다.
평면에 그냥 박스가 올려져있다고 생각해보자. 중력에 의해 그 오브젝트는 -10의 가속도를 가지고 있다. 그리고 적분을 통해 그 오브젝트는 아래로 내려가는 속도(velocity)를 갖게 될 것이다. 근데, 그 오브젝트가 땅으로 꺼지지 않기 위해, 땅에서 올라가는 힘이 생기게 된다. 뉴턴의 법칙에 따라, 힘에는 작용/반작용이 있다는 것이다. 그래서 같은 크기의 반대로 작용하는 힘이 생기게 된다. (아마 물리 시간에는 수직항력이라고 배운 것을 떠올리면 될 것이다) 그래서 우리도 이렇게 reaction force를 계산해서 적용을 해줘야 한다.
하지만, 우리의 엔진은 현재 impulse기반이고, force를 계산해서 적용하지 않는다. 이러한 reaction force를 하는 것은 simultaneous contact resolution 기법에서 하는 것이므로 Davia Baraff 논문과 Erin Catto 논문을 참고하면 될 것이다.
어쨋든 글쓴이는 이 문제를 해결하기 위해 Micro-Collisions Modeling을 도입한다. 즉, reaction forces를 일련의 impulses로 대체하는 것이다.
즉, 한 impulse는 single moment of time에 적용되는 한 force로 생각
한 force는 일련의 impulses로 생각.
더 간단히 이해되게 표현하자면, Velocity Resolution System은 가속도로 인해 생긴 velocity를 제외하는데 필요한 impulse를 계산한다.
이 impulse는 작은 little impulse인데, micro-collisions라고 불려진다.
글쓴이에 따르면, micro-collisions는 reaction forces를 생성하는 잘 알려진 기법이지만, 불안정한 시뮬레이션을 생산한다는 부당한 평판을 받고있다고 한다.
어쨋든 방법은
1. 이전의 rigidbody update에서 가속도로 축적된 어떤 velocity를 제거한다. 또 추가적으로, 2. 매우 느린 speed일 경우에 복구계수를 줄여서 reaction forces를 반영해준다.
void GPED::Contact::calculateDesiredDeltaVelocity(real duration) { const static real velocityLimit = (real)0.25f; // Calculate the acceleration induced velocity accumulated this frame real velocityFromAcc = 0; if (body[0]->getAwake()) { velocityFromAcc += glm::dot(body[0]->getLastFrameAcceleration() * duration, contactNormal); } if (body[1] && body[1]->getAwake()) { velocityFromAcc -= glm::dot(body[1]->getLastFrameAcceleration() * duration, contactNormal); } // If the velocity is very slow, limit the restitution real thisRestitution = restitution; if (real_abs(contactVelocity.x) < velocityLimit) { thisRestitution = (real)0.0; } // Combine the bounce velocity with the removed acceleration velocity /* 18-11-05 Chanhaeng Lee // desiredDeltaVelocity = -contactVelocity.x - thisRestitution * (contactVelocity.x - velocityFromAcc); I changed the above line into the bottom one to help understand the collision resolution process. the equation of solving the impulse j is, j = -(1 + Restitution) * dot(V^(AB), contactNormal) But you don't need to do the dot product, because the contactVelocity is in the contact coordinates. */ desiredDeltaVelocity = -(1 + thisRestitution) * contactVelocity.x + thisRestitution * velocityFromAcc; }
여기에서 필요한 deltaVelocity를 구할 때, 지난 프레임의 가속도에 duration을 곱해서 속도로 적분을 해주고, 그것을 contact normal 방향 컴포넌트 만을 얻는다. 만약 두 번째 body가 있다면, velocityFromAcc에서 빼준다. 왜냐하면 여기에서 relative Velocity를 얻고 있기 때문이다.
그리고 원래 요구되는 DeltaVelocity는
여기에서
이것으로 바뀌게 된다. 위의 식을 간단하게 정리한다면
이렇게 되어서 여러가지로 쓸 수 있는데, 나는 명료함을 위해서 위의 식의 두 번째를 썼다.
이것으로 1번은 되었고, 2번은 더 간단한데
// If the velocity is very slow, limit the restitution
real thisRestitution = restitution;
if (real_abs(contactVelocity.x) < velocityLimit)
{
thisRestitution = (real)0.0;
}
이 부분이 restitution을 줄이는 것이다, 그러니까 평면위에 박스가 있다했을 때, 그 cotact, 정확히는 resting contacts인데, 그것으로 인해 생기는데 contactVelocity는 높지 않을 것이다. 왜냐하면 프로세스로 생각한다면, 중력가속도에 의해 속도로 적분될 때, delta time이 높지 않을 때, velocity는 큰 값을 갖지 않게 되고, 이에 따라, 평면과 박스가 충돌하고나서 contacts를 생성하고, contact resolution system에 들어와서, 그 contacts점의 속도를 계산한다면, 즉 v = linear Velocity + cross(RotationVelocity, relativePosition)이 그리 크지 않게 된다. 그리고 우리는 contact space에서 계산했기 때문에, abs를 통해서 그 크기를 구할 수가 있고, velocityLimit보다 작다고 된다면, 그것은 resting contact로 판단하여 restitution을 줄여서, 오브젝트가 vibration을 하려는 것을 줄이려한다.
나는 vsync를 끄고 안끄고 했을 때 차이가 극명하게 나는 것을 보았는데, 일정 시간으로 적분을 했을 때, 즉, vsync를 했을 때, object들의 vibration은 그렇게 크지 않고 그나마 visual적으로 나은데, vsync를 끄고 가변 delta time에 들어올 때는 매우 vibration이 컸다.
- Friction
마찰은 크게 두 가지로 나눠진데, static and dynamic friction.
static friction은 간단히 설명하면 오브젝트가 정지해있을 때, 움직이기까지 방해하는 힘이다.
dynamic friction은 간단히 설명하면 오브젝트가 움직이고 있을 때, 움직이는 것을 방해하는 힘이다.
먼저 static friction을 설명할 텐데, 물체가 정지해있을 때 static friction을 받고, 그 static friction을 압도하는 힘이 있을 때 dynamic friction의 세계로 넘어가고, 이제 움직이면서 dynamic friction을 받기 때문이다.
아 그전에 또 friction의 category가 또 있는데, isotropic과 anisotropic friction이 있다.
isotropic friction은 모든 방향에 대해 같은 friction을 갖는 것을 말한다.
anisotropic friction은 다른 방햐에 대해 다른 friction을 갖는 것을 말한다.
글쓴이에 의하면 anisotropic friction을 게임에 적용하는 사례는 거의 없다고 말한다. 그래서 우리가 현재 개발하는 마찰은 isotropic friction이 된다.
static friction의 식들이 있다. 그것은
r은 contact normal 방향의 reaction force이고, f_static은 생성된 friction force, u_static은 static friction 계수이다. 마찰계수는 접촉점에서 모든 물질의 특성을 단일 숫자로 캡슐화 한다. 이것은 두 오브젝트에 의존한다. 그러나 이것은 경험적인 양이고, 물리책을 참고하거나 실험적으로 설정되어질 수 있다.
그리고 나는 이 식을 정확히 이해가 안되는데 어쨋든 써놓는다.
여기에서 f_planar은 contact의 평면에서 오브젝트에 가해지는 total force이다. reaction force인 r과 planar force는 적용되는 total force로 부터 계산된다.
d는 contact normal이고, f는 발휘된 total force이다. 그리고
다시 말해서, 최종 planar force는 contact normal의 방향의 컴포넌트가 제거된 total force이다. total force에 reaction force가 더해져서, contact normal component가 제거된다. 이것은 static friction이 normal reaction force에 의존한다는 것이다. 직관적으로 생각해 본다면 맞는 말이 된다.
Dynamic Friction은 kinetic friction이라고 불려지는데, static friction과 비슷하게 행동하고, 다른 마찰 계수를 갖는다. 식은
u_dynamic은 동적 마찰 계수이다. 마찰의 방향이 바뀌었다는 것을 알아야 한다. planar force의 반대방향이라기 보다, 그것은 그 오브젝트의 속도와 반대방향으로 작동하고 있다. 이것은 어떤 말이냐면, 움직이는 오브젝트에 힘을 가하는 것을 중단한다면, 마찰력은 바로 중단되지 않고, 동적 마찰에 의해 중단될 때 까지 느려질 것이다. 동적 마찰계수도 물리 책이나 실험을 통해서 설정될 수 있다.
그런데 게임 물리엔진에서, static and dynamic friction을 구분하지 않는다. 그냥, 정지해있을 때, 마찰력은 static으로 작용하고, 움직일 떄 dynamic friction에 작용한다. static은 가해지는 힘에 반대로 작용하고, dynamic은 움직임의 방향에 반대로 작용한다. 이 게임 엔진 모델에서는 두 유형의 마찰을 한 가지 값으로 합친다. Rolling Friction이라는 것도 있는데, 게임엔진에서는 구현될 필요가 없다.
- 구현
우리는 impulse-based resolution system 이기 때문에, 마찰력이 impulses와 velocity의 관점에서 무엇을 하는지 이해해야 한다.
static friction은 한 force가 오브젝트에 가해질 때 움직이는 것을 막는다. 그것은 contact plane에서 그 오브젝트의 속도를 0으로 유지하려고 작용한다. 나는 글쓴이의 설명이 잘 이해가 되지 않지만, 우리가 static friction을 구현하기 위해서 오브젝트에 가해지는 힘에 반대방향으로 적용되어야 한다는 것이다. 그 적용되는 것이 static friction이 된다. 그래서 contact space에서
// Find the target velocities to kill
glm::vec3 velKill(desiredDeltaVelocity, -contactVelocity.y, -contactVelocity.z);
이게 우리가 contact space에서 현재 static friction이 된다. static friction이 우리는 impulse가 된다는 것을 염두에 두어서, 이 속도를 가지는데 필요한 impulse를 계산해야 한다. 즉, 최종 사용되는 식이
// Find the impulse to kill target velocities
impulseContact = impulseMatrix * velKill;
이렇게 되는 것이다. 우리는 velKill인 friction impulse를 만들기 위해 impulseMatrix를 형성해야 한다. impulseMatrix는 Chapter 14에 했듯이, 변화되는 단위 속도를 만들기 위해 필요한 linear and angular component가 적용된 것이다. 이것을 구하기 위해 글쓴이는 다음의 순서를 제시한다.
- 우리는 contact 기준의 좌표계에서 작업한다: 이것은 수학을 훨씬 더 쉽게 만든다. 우리는 이러한 새로운 좌표계에서 변환시키는 transform matrix를 만든다.
- 우리는 각 오브젝트 위에 있는 contact point에서의 속도에서 unit impulse당의 변화를 처리한다. impulse가 linear and angular motion을 발생시키기 때문에, 이 값은 두 컴포넌트들을 고려할 필요가 있다.
- 우리는 우리가 보기 원하는 속도 변화를 알 것이다. 그래서 우리는 마지막 단계의 결과를 invert 시켜서, 어떤 주어진 속도 변화를 하는데 필요한 impulse를 찾는다.
- 우리는 contact point에서 그 separating velocity가 무엇이 되어야 하는지를 처리하고, 현재 closing velocity가 무엇인지, 그리고 그 둘 사이의 차이가 무엇인지 처리 한다. 이것이 그 속도에서의 요구되는 변화이다.
- 속도 변화로부터, 우리는 발생되어야 할 impulse를 계산할 수 있다. 우리는 그 impulse를 linear and angular components로 분해하고, 그것들을 각 오브젝트에 제공한다.
이것의 applyVelocities까지의 일반적인 프로세스지만, 여기에서 이제 frictionimpulse를 계산하기 위해 세부사항이 조금 달라진다.
1번은 contactBasis() method를 통해서 처리된다. Step2는 수정이 필요하다. 현재 그것은 contact normal의 방향에서 단위 impulse가 주어진다면 속도 변화를 처리한다. 우리는 이제 모든 contact direction 세 방향을 다룰 필요가 있다. 왜냐하면 마찰력 때문이다.
contactImpulse는 벡터로 표현되는데, x는 contact normal 방향의 impulse를, 그리고 y와 z는 contact의 평면에서의 impulse를 나타낸다.
step2의 결과는
일반적인 root에서
// Build a vector that shows the change in velocity in world space
// for a unit impulse in the direction of the contact normal.
glm::vec3 deltaVelWorld = glm::cross(relativeContactPosition[0], contactNormal);
deltaVelWorld = inverseInertiaTensor[0] * deltaVelWorld;
deltaVelWorld = glm::cross(deltaVelWorld, relativeContactPosition[0]);
// Work out the change in velocity in contact coordinates.
real deltaVelocity = glm::dot(deltaVelWorld, contactNormal);
// Add the linear component of velocity change
deltaVelocity += body[0]->getInverseMass();
// Check if we need to the second body's data
if (body[1])
{
deltaVelWorld = glm::cross(relativeContactPosition[1], contactNormal);
deltaVelWorld = inverseInertiaTensor[1] * deltaVelWorld;
deltaVelWorld = glm::cross(deltaVelWorld, relativeContactPosition[1]);
// Add the change in velocity due to rotation
deltaVelocity += glm::dot(deltaVelWorld, contactNormal);
// Add the change in velocity due to linear motion
deltaVelocity += body[1]->getInverseMass();
}
// Calculate the required size of the impulse
impulseContact.x = desiredDeltaVelocity / deltaVelocity;
impulseContact.y = 0;
impulseContact.z = 0;
이런 식이였다. 그러니까, chapter 14에서 설명했듯이, impulse를 구하기 위해서, 우리는 desired Deltavelocity를 구했었고 미리, 그 unit velocity를 만들어내는데 필요한 impulse를 구하기 위해서, deltaVelocity를 구했는데, 이 detlaVelocity라는 것은 그 적용점에서 오브젝트가 linear and angular 움직임을 보이기 위해 저항되는 숫자라 생각하면 된다. 이 method에서는 너도 보이듯이, contactNormal를 중심으로 한 방향으로 하지만, 우리는 이제 matrix를 만들어서 모든 세 방향을 고려할 수 있게 해야 한다.
그래서 변환된 방식이
// The equivalent of a cross product in matrices is multiplication
// by a skew symmetric matrix - we build the matrix for converting
// between linear and angular quantities.
glm::mat3 impulseToTorque;
getSkewSymmetric(relativeContactPosition[0], impulseToTorque);
// Build the matrix to convert contact impulse to change in velocity in world coordinates
glm::mat3 deltaVelWorld = impulseToTorque;
deltaVelWorld *= inverseInertiaTensor[0] * impulseToTorque;
deltaVelWorld *= -1;
// Check if we need to add body two's data
if (body[1])
{
// SEt the cross product matrix
getSkewSymmetric(relativeContactPosition[1], impulseToTorque);
// Calculate the velocity change matrix
glm::mat3 deltaVelWorld2 = impulseToTorque;
deltaVelWorld2 *= inverseInertiaTensor[1] * impulseToTorque;
deltaVelWorld2 *= -1;
// Add to the total delta velocity
deltaVelWorld += deltaVelWorld2;
// Add to the inverse mass
inverseMass += body[1]->getInverseMass();
}
// Do a change of basis to convert into contact coordinates.
glm::mat3 deltaVelocity = glm::transpose(contactToWorld) * deltaVelWorld * contactToWorld;
// Add in the linear velocity change
deltaVelocity[0][0] += inverseMass;
deltaVelocity[1][1] += inverseMass;
deltaVelocity[2][2] += inverseMass;
이렇게 되는데 간단하게 하기 위해서 body[1]이 없다고 생각하고 밑의 deltaVelocity가 어떤 행렬들이 결합되어있는지 써보자면
deltaVelocity = WorldToContact * relativeContactPosition[0] * inverseInertiaTensor[0] * relativeContactPosition[0] * -1 * contactToWorld
이렇게 되어있다. 밑에 inverseMass는 linear한 부분이니까 일단 빼고 angualr component만 보자.
내가 쓴 위의 식은 모두 다 행렬이라는 것에 유의해라. 우리가 이것을 통해 결국하려는 것은
// Invert to get the impulse needed per unit velocity
glm::mat3 impulseMatrix = glm::inverse(deltaVelocity);
// Find the target velocities to kill
glm::vec3 velKill(desiredDeltaVelocity, -contactVelocity.y, -contactVelocity.z);
// Find the impulse to kill target velocities
impulseContact = impulseMatrix * velKill;
impulseContact를 구하는데 세 방향을 고려한 것이다. 그래서 마찰력이 없을 때
// Calculate the required size of the impulse
impulseContact.x = desiredDeltaVelocity / deltaVelocity;
impulseContact.y = 0;
impulseContact.z = 0;
이렇게 했던 것이 저기 위에서는 행렬을 inverse해서 한 것이다. 그렇다면 deltaVelocity를 역행렬 한다면 역행렬은 다 반대로 된다. 그렇다면 역행렬 성질에 의해
deltaVelocity는 inverse(contactToWorld) * -1 * inverse(relativeContactPosition[0]) * inverse(inverseInertiaTensor[0]) * inverse(relativeContactPosition[0]) * inverse(WorldToContact) 가 된다.
그렇다면 impulseMatrix * velKill를 고려한다면 우리가 아는 chapter 14에서 j를 구하는 식이 된다. 여튼 이 정도로 하면 세 방향에 대한 회전에 대한 component를 찾는게 직관적으로 나마 이해가 될 거라 생각한다. 그리고 좌표변환 행렬이 들어가 있어서 자연스럽게 velKill이 어떻게 바뀌는지 보면 inverse(WorldToContact)는 결국에 ContactToWorld가 된다. 기저행렬의 역행렬은 그것의 전치행렬과 같기 때문이다. 여튼, velKill은 world space로 들어가게 된다. 그리고 relativeContactPosition과 외적의 component, inverseInertiaTensor의 컴포넌트 그리고 다시 그것을 통해 relativeContactPosition[0]과 외적을 하는데 이 때, -1을 곱한 이유가 원래 그 component를 구하는데 곱을 하는 순서는
v = cross(relativeContactPostiion[0], contactNormal)
v = invserInertiatensor * v;
v = cross(v, relativeContactPosition[0]);
이렇게 되는데 반대 방향이여서 -1을 곱해준거다. 그리고 마지막으로 inverse(contactToWorld)인 WorldToContact를 곱해주어서 다시 contact space로 바꿔준다. 따라서,
// Find the impulse to kill target velocities
impulseContact = impulseMatrix * velKill;
contact space에 있는 impulse를 구할 수 있게 된 것이다. 즉 회전에 의해서 발생하는 velocity change component를 구해서 impulse를 구한 것이다. linear motion부분은 간단히 저번에도 헀듯이, inverse mass만 있으면 된다. frictionless impulse를 구할 때 inverseMass를 그냥 더해주는 것을 행렬 상에서 생각한다면 그것은 대각행렬에 inverseMass를 더해주는 것과 같다. 따라서,
// Add in the linear velocity change
deltaVelocity[0][0] += inverseMass;
deltaVelocity[1][1] += inverseMass;
deltaVelocity[2][2] += inverseMass;
// Invert to get the impulse needed per unit velocity
glm::mat3 impulseMatrix = glm::inverse(deltaVelocity);
// Find the target velocities to kill
glm::vec3 velKill(desiredDeltaVelocity, -contactVelocity.y, -contactVelocity.z);
// Find the impulse to kill target velocities
impulseContact = impulseMatrix * velKill;
이렇게 linear component까지 구해줘서 impulseMatrix를 우리는 구할 수 있게 된다.
// Find the impulse to kill target velocities
impulseContact = impulseMatrix * velKill;
// Check for exceeding friction
real planarImpulse = real_sqrt(
impulseContact.y * impulseContact.y
+ impulseContact.z * impulseContact.z);
if (planarImpulse > impulseContact.x * friction)
{
// We need to use dynamic friction
impulseContact.y /= planarImpulse;
impulseContact.z /= planarImpulse;
impulseContact.x = deltaVelocity[0][0] +
deltaVelocity[1][0] * friction * impulseContact.y +
deltaVelocity[2][0] * friction * impulseContact.z;
impulseContact.x = desiredDeltaVelocity / impulseContact.x;
impulseContact.y *= friction * impulseContact.x;
impulseContact.z *= friction * impulseContact.x;
}
return impulseContact;
이제 그 구한 impulse를 통해 밑의 식이 friction의 주요식인데 이제 이것을 이해해보자.
그리고 friction 식들을 모두 모아보자.
글쓴이는 f_planar 를 impulseContact.y와 z 컴포넌트의 제곱의 합으로 구했다. impulseContact.x가 contact 방향의 force라인게 이해가 된다면
if (planarImpulse > impulseContact.x * friction)
impulseContact.x * friction은 == u_static |r|이 된다.
u_static은 마찰계수여서 friction이 그 역할이고 |r|은 contact component의 force그러니까 우리의 엔진에서 impulse가 된다.
여기에서 생각을 해본다면, 글쓴이는 f_static 값이 둘 중 작은 것이 된다고 했는데, 그냥 f_planar이 정의를 생각해본다면, 주어진 force(우리는 impulse)에 contact normal의 방향을 뺸 것이다. 근데 두 번째 식은 f_static의 범위를 만족하기 위해서이다. 어쨋든 f_static이 코드상에서 f_planar로 작용한다는 것을 알면된다. 그래서 f_planar은 contact 좌표계 덕분에 그냥 간단히 impulseContact의 y와 z 컴포넌트로 구할 수 있고 제곱의 합으로 그 크기를 구하게 된다. 그래서 planarImpulse가 더 크다면 이제 dynamic friction의 세계로 들어가게 된다. 그게 아니면 그냥 velKill이 그대로 적용될 것이여서, static friction이 적용될 것이다.
// We need to use dynamic friction
impulseContact.y /= planarImpulse;
impulseContact.z /= planarImpulse;
impulseContact.x = deltaVelocity[0][0] +
deltaVelocity[1][0] * friction * impulseContact.y +
deltaVelocity[2][0] * friction * impulseContact.z;
impulseContact.x = desiredDeltaVelocity / impulseContact.x;
impulseContact.y *= friction * impulseContact.x;
impulseContact.z *= friction * impulseContact.x;
그래서 이제 dynamic friction을 이해하려 해보자.
우선 unit velocity 가 필요하고 그리고 마찰계쑤 그리고 contact normal 방향의 impulse가 필요하다. 그리서 동적 마찰력 방정식에 따라 단위 방향을 갖게하기 위해 planarImpulse로 나눠주고, 나머지는 방정식 식에 맞게 friction == 마찰계수와 |r| == impulseContact.x == contact normal 방향의 impulse(force)를 곱해준다.
impulseContact.x = deltaVelocity[0][0] +
deltaVelocity[1][0] * friction * impulseContact.y +
deltaVelocity[2][0] * friction * impulseContact.z;
근데 솔직히 여기는 왜 해주는지 모르겠다. 여기에 대한 설명도 자세히 되어있지 않다.
impulseContact.x = desiredDeltaVelocity / impulseContact.x;
그리고 이 부분은 뭐 desired DeltaVelocity에 저기 이해안된 식을 나누어서 다시 계산해주는데 내 추측으로는 동적 마찰력에서는 정적 마찰력보다 더 적게 들어가니까 저렇게 friction을 추가해서 이렇게 하는 것 같다.
이 챕터는 내가 보기에 설명도 더 어렵게 해놓고, 추가적인 세부적인 설명도 부족한 챕터였다. 내 생각에는 어쨋든 primer로서 그냥 생각하고, 나중에 force-based simultaneous contact resolution engine을 만들 때 이 마찰력이 정확히 구현되도록 해야겠다. Erin Catto에서도 friction에 대한 부분이 있으니 David Baraff -> Erin Catto 논문으로 넘어가도록 해야 겠다.
Chapter 16 Stability and Optimization
이 파트는 그냥 책보면서 순서대로 정리를 하겠다. 근데 요점정리만.
- Stability
안정성 문제가 발생하는 곳
- 개별적으로 합리적으로 행동하는 소프트웨어들 사이에서 좋지 않은 상호작용
- 사용된 방정식의 부정확성 또는 우리가 한 가정의 불리한 효과들
- 컴퓨터에 의해 수행되는 수학의 내재된 부정확성
글쓴이가 발견한 상대적으로 안정성 문제를 해결하기 쉬운 5가지 문제들
- Transform matrices는 오염될 수 있고, 회전과 평행이동 외에도 skew를 수행할 수 있다.
- 빠르게 움직이는 오브젝트들은 가끔씩 이상하게 충돌에 반응한다 (이것은 실제로 충돌이 탐지되는 것과 독립적이다: 빠르게 움직이는 오브젝트는 충돌이 탐지되지 않고 다른 오브젝트를 지나갈 수 있다.)
- 경사가 있는 평면에 있는 오브젝트들 (또는 경사진 표면을 가진 다른 오브젝트 위에서 resting하는 오브젝트들)은 느리게 미끄러지는 경향이 있다.
- rapid succession에 있는 빠른 속도의 충돌의 한 쌍을 경험하는 오브젝트들은 갑자기 땅과 관통한다.
- 그 시뮬레이션은 크고 작은 양이 섞여있을 때 비현실적으로 보여질 수 있다: 크고 작은 질량, 크고 작은 velocities, 크고 작은 rotations, 등등
1. Quaternion Drift
transformation matrix는 position vector와 orientation quaternion으로 생성되는데, 처리되는 동안 numerical errors를 겪는다. 그래서 body의 integration routine에 글쓴이는 quaternion을 normalization 하는 것을 추가했다.
2. 경사에서의 interpenetration
이것은 interpenetration resolution 과정에서 경사 위에 있는 오브젝트를 원래 위치로 옮기지 않고 조금 아래로 하게 하여 충분히 마찰력이 있음에도 경사 위에 오브젝트가 조금씩 미끄러지게 만드는 상황이다.
이것의 해결방법은 contact의 relative velocity 계산에서 처리를 한다. 즉 contact plane에서 force에 의한 축적을 가지는 어떤 velocity를 제거하는 것이다. 이것은 오브젝트가 contact normal의 방향으로 경사쪽으로 움직이도록 허용하지만, 그것을 따라서 움직이게 하지 않는다. {사실 나는 이게 무슨 말인지 이해가 되지 않는다. 설명이 부족한 것 같다}
glm::vec3 GPED::Contact::calculateLocalVelocity(unsigned bodyIndex, real duration) { RigidBody* thisBody = body[bodyIndex]; // Work out the velocity of the contact point glm::vec3 velocity = glm::cross(thisBody->getRotation(), relativeContactPosition[bodyIndex]); velocity += thisBody->getVelocity(); // Turn the velocity into contact-coordinates. glm::mat3 WorldToContact = glm::transpose(contactToWorld); glm::vec3 contactVelocity = WorldToContact * velocity; // Calculate the amount of velocity that is due to forces without reactions. glm::vec3 accVelocity = thisBody->getLastFrameAcceleration() * duration; // Calculate the velocity in contact-coordinates. accVelocity = WorldToContact * accVelocity; // We ignore any component of acceleration in the contact normal // direction, we are only interested in planar acceleration accVelocity.x = 0; // Add the planar velcities - if there's enough friction they will // be removed during velocity resolution contactVelocity += accVelocity; // And return it. return contactVelocity; }
일단 코드만 보면서 이해하자면, 이전 forces에 의해 쌓인 가속도로부터 velocity를 구하고, 그것을 contact space로 변환해서, contact normal 방향은 제외하고 그것의 plane 방향을 더해준다. 코드를 봐도 모르겠다...
3. Integration Stability
우리가 사용하고 있는 적분 방식은 Newton-Euler integration인데, Newton은 linear component를, Euler는 Angular Component를 담당한다. 그 식은
3. Integration Stability
우리가 사용하고 있는 적분 방식은 Newton-Euler integration인데, Newton은 linear component를, Euler는 Angular Component를 담당한다. 그 식은
이것은 1차 미분에 의존하고, 그러므로 "first-order"라고 불려진다. 전반적인 방법은 완전히 "first-order Newton-Euler" or Newton-Euler 1이라고 불려진다.
Newton-Euler 2는 근사이고, 그 방정식은
이것은 2차 미분에 의존한다. angular updates에 대해서 동일한 식으로, 우리는 second-order Newton-Euler integrator을 얻는다. 그것은 업데이트된 위치를 결정할 때 가속도를 고려한다. 우리가 chapter 3에서 보았듯이, 높은 프레임율에 대해 그 t^2 항이 너무 작기 때문에 우리는 그 가속도항을 무시하는 것이 낫다. 이것은 그러나 가속도가 매우 클 때 사실이 아니다. 이 경우에 그 항은 중요할지도 모르고, Newton-Euler 2로 옮기는 것을 이익이 된다.
Runga-Kutta 4
Newton-Euler integrators 둘 다 가속도가 전체 업데이트 동안 상수로 남을 것이라고 가정한다. 우리가 챕터 6에서 스프링에 대해 볼 때 보았듯이, 가속도가 한 업데이트의 과정동안 변하는 방식은 매우 중요할 수 있다. 사실, 가속도가 변하지 않는다고 가정하여, 우리는 극적인 불안정성과 시뮬레이션의 완전한 고장을 볼 수 있다.
스프링은 가속도를 빠르게 바꾸는 유일한 것이 아니다. 몇 가지 resting contacts들의 패턴들 (특히 simultaneous velocity resolution algorithm이 사용될 떄)은 비슷한 효과를 가질 수 있고, 오브젝트 stacks의 vibration or a dramatic explosion을 이끈다.
두 문제들에 대해, 부분적인 솔루션은 중간 단계에서 필요한 가속도를 처리하는 방법에 있고, 그것이 바로 fourth-order Runga-Kutta algorithm(또는 간단히 RK4) 이다.
4. The Benefit of Pessimistic Collision Detection
예제의 케이스를 해결하기 위해서 pessimistic collision detection을 하는 것인데, 그 말은 충돌 탐지가 가깝지만 실제로 닿고있지 않는 contacts를 반환하도록 하는 것이다. 이것은 한 오브젝트를 둘러싸는 collision geometry를 확장시키고, penetration value에 대한 offset으로 사용하여 얻어질 수 있다. 만약 collision geometry가 그 오브젝트의 시각적 representation보다 한 단위 더 크다면, 그러면 탐지된 충돌의 관통값으로 부터 그 한 단위가 뺄셈된다.
실제로 이것의 어떤 효과를 보는 것은 드물다. 그런데 길고 가벼운 오브젝트들(막대같은)과 ground 사이의 충돌에서는 가장 눈에 띈다. 충돌 탐지를 위해 ground를 조금 더 높게 이동시키고, 생성된 ground collisions로부터 적은 양을 뺄셈하는 것은 trivial change이다.
5. Changing Mathematical Accuracy
float -> double로 바꾸는 것이다. 매크로를 사용해서 함.
- Optimizations
1. Sleep
sleep system의 두 가지 컴포넌트 : 오브젝트들을 sleep 들게 하는 알고리즘과 그것들을 깨우는 알고리즘.
rigidBody class에 다음의 세 개의 데이터 멤버를 추가
- isAwake : 현재 잠들어 있는지 깨워있는지. 깨어있으면 처리될 필요가 있다.
- canSleep : 사용자가 끊임없이 사용하는것은 안재워서 비쥬얼을 향상시켜야 한다.
- motion : 현재 움직임의 속도를 추적한다. 이것은 오브젝트가 잠들어야 하는지를 결정하는데 중요하게 쓰인다.
그리고 motion < sleepEpsilon이면 setAwake(false)로 오브젝트를 잠재운다.
sleep epsilon은 trial-and-error process로 결정된다. 너무 낮게 설정하면 오브젝트들이 결코 잠들지 않을 것이고, 너무 높게 잡으면 오브젝트들이 잠드렉 되고 이상하게 보여질 수 있다. 그래서 저자는 이상한 움직이는 중간에 멈추는것이 명백히지기 전만큼의 가능한한 높게 잡는 경향이 있다.
알고리즘은 간단하지만, motion 값을 계산하는 것에 의존한다. 그 motion value는 linear and angular velocity를 단일 스칼라 값으로 캡슐화 할 필요가 있다. 이것을 하기 위해 오브젝트의 total kinetic energy를 사용한다. 챕터 3에서 한 파티클의 kinetic energy는
m은 body의 질량이고 dot(p)는 그것의 linear velocity이다. 비슷한 식이 rotating rigid body의 kinetic energy에도 유효하다:
i_m은 boy의 회전축에 대한 moment of inertia이다 (즉, 그것은 스칼라 양이다) 그리고 dot(theta)는 그것의 angualr velocity이다.
우리는 kinetic energy를 motion에 대한 값으로 사용할 수 있지만, 그것은 다른 질량에 대해 한 문제를 갖는다 : 두 개의 동일한 오브젝트들이 그것들의 질량에 따라 다른 시간에 잠들 것이다. 이것을 피하기 위해 우리는 식으로 부터 질량을 제거하고 다음을 얻는다
이렇게 계산하고 다음 단계는 이 값이 안정된지를 체크하는 것이다. 우리는 몇 프레임에 대해 현재 motion의 기록을 유지하여 이것을 하고, 그것이 얼마나 변하는지를 본다. 이것은 recency-weighted average(RWA)로 깔끔하게 추적될 수 있는데, 이것은 나의 프로그래밍 레퍼토리에서 가장 많이 사용되는 도구중 하나이다.
recency-weighted average는 다음으로 업데이트 된다
rwa = bias * rwa + (1 - bias) * newValue;
그것은 지난 몇 가지 값들의 평균을 유지하는데, 좀 더 최근의 값을 좀 더 중요하게 만든다. 그 bias parameter는 이전 값에 얼마나 중요성이 갈지를 통제한다. 0의 bias는 RWA가 새로운 값과 같도록 만든다, 업데이트 할 떄마다 (즉, 어떠한 평균화도 없다). 1의 bias는 전적으로 새로운 값을 무시한다.
그 RWA는 input을 smoothing하거나 input이 안정화되었는지를 체크하는 훌륭한 자치이다. 우리의 경우에
motion = bias * motion + (1 - bias) * currentMotion;
만약 currentMotion이 sleep epsilon value아래 였지만, 이전 몇 프레임에서 그 오브젝트가 많이 움직였따면, 그러면 전체 motion value는 여전히 높을 것이다. 한 오브젝트가 잠깐을 썻을 때, 움직이지 않는 것은 RWA가 epsilonv value아래로 떨어지게 만들 것이다.
오브젝트가 매우 빠른 속도로 움직일 수 있기 때문에, 스피드의 짧은 시간의 burst는 RWA가 매우 높아지게 만들고, 적당한 수준으로 돌아오는데 많은 시간이 걸릴 것이다. 이것을 방지하기 위해, 오브젝트들이 더 빠르게 잠들게 하기 위해, 저자는 rwa의 값을 제한하하는 코드를 추가 했다:
if(motion > 10 * sleepEpsilon) motion = 10 * sleepEpsilon;
RWA의 bias는 프레임의 duration에 의존적이여 한다. 더 긴 프레임들은 현재 값이 RWA에 더 짧은 프레임들 보다 더 많이 영향을 미치도록 허용해야 한다. 만약 그렇지 않는다면, 오브젝트들을 더 빠른 프레임율에는 더 빠르게 잠들 것이다. 우리는 우리가 damping에 했던 방식과 똑같이 이것을 할 수 있다:
real bias = real_pow(baseBias, duration);
baseBias는 one-second frames에 대해 우리가 기대하는 bias이고, 나는 보통 0.5 ~ 0.8의 값을 여기에서 사용하지만, 또 다시 몇 가지 실험이 필요하다.
우리는 새로운 충돌에 반응해야 할 때 오브젝트들을 깨울 필요가 있다. 자고 있는 오브젝트들 사이의 충돌들이 생성되고, 자동적으로 collision resolution system에 의해 무시된다.
한 새로운 오브젝트가 (예를들어, 플레이어 또는 발사체) 들어오고 잠자는 오브젝트와 충돌할 때, 우리는 충돌에 의해 영향 받을 수 있는 모든 오브젝트들이 일어나기를 원한다. 어떤 특별한 충돌에 대해, 이것은 만약 관련된 한 body가 잠들어 있고 다른 것이 깨어있으면 그러면 그 잠자는 body가 깨어나질 필요가 있다는 것을 의미한다. 그래서 우리는 Contact class에 이것을 하기 위해 matchAwakeState() 메소드를 추가한다.
그래서 collision resolution 때 이 메소드는 항상 호출된다.
force generator같은 경우도 항상 깨워야 하는데, force generator에 setawake를 하는게 아닌, rigidbody class의 메소드인 addforce, addforceatpoint, addtorque에 추가를 해서 기억하기 어려운 것들을 제외했다.
2. Margins of Error For Penetration and Velocity
또 다른 할 가치가 있는 최적화는 penetration and velocity resolution 알고리즘을 극적으로 가속화하는 것이다.
즉 간단히 말해서, penetration이나 desiredDeltaVelocity 값이 충분히 크지 않아서 그 contact를 resolution하지 않아도 될 정도면 그 resolution process를 빨리 끝내서 가소갛는 것이다.
3. Contact Grouping
우리의 resolution algorithm은 이전의 resolution에 의해 영향을 받은지 보기 위해 모든 가능한 contacts를 체크한다. 또한 해결할 현재 가장 심각한 contact를 찾기 위해 모든 contacts를 훑는다. 이러한 두 연산 모두 더 긴 contacts lists에 대해 더 긴 시간이 걸린다.
따라서, Contact Grouping은 서로 영향 받는 contacts만을 모아서 묶어서 contact resolution에 보낸다.
이 batching은 충돌탐지기의 일이 될 수 있거나, collision resolution system에 의해 전체 contacts의 리스트가 분리되어 batches로 처리될 수 있다.
만약 구현될 수 있다면 첫 번째 접근법이 가장 좋다. 그러나 닿고있지 않는 가까운 오브젝트들의 집합에 대해, 충돌탐지기는 실제로 너무 큰 batches를 반환할 것이다. 만약 그 game level이 매우 그러한 상황을 가진다면, 그것은 충돌 탐지 batching뿐만 아니라 resolver에도 batching을 추가하는 이후의 성능을 개선할 수 있다.
batching process는 contacts의 전체 리스트를 batches로 나눈다. 그것은 한 contact에서 시작하여 contacts와 rigid bodies의 조합을 따라서 한 batch에 있는 모든 contacts를 찾는다. 이것은 그러고나서 processing을 위해 보내진다. 그것은 그러고나서 batch에 아직 있지 않은 또 다른 contact를 찾고, 또 다른 batch를 찾기 위해 contacts를 따라간다. 이것은 처리되지 않은 contacts들이 없는 한 반복된다.
batching processor를 구현하는 것은 각 contact에 관련된 rigid bodies를 빠르게 찾는 것과(우리는 이미 contact data structure에 관련된 rigid bodies를 저장했다), 각 rigid body에 있는 contacts들을 찾을 수 있는 것을 포함한다.
그 두 번째 구현사항은 우리가 chapter 14에서 이야기 한, linked list contacts 구조를 이용한 것을 사용한다면 구현될 수 있다.
4. Code Optimization
코드 최적화를 하기 전에 좋은 profiler를 얻어서, 성능 문제를 일으킨다고 너가 증명할 수 있는 소프트웨어의 부분만을 최적화 해라.
1) contact data Caching
한 contact를 여러 번 resolving할 때, 우리는 현재 deltaVelocity와 다른 값들을 재 계산한다. 이러한 것은 대신에 contact data structure에 저장되어 처음에 필요할 때 만 계산되어질 수 있다. {난 이게 뭔 말인지 모르겠다. 재 계산하는건 맞지만, 처음에 필요할 때 계산되어 저장된 것을 계속 쓰면 다른 contact에 연관된 rigidbody의 deltaVel은 당연히 바뀌어야 하는데, 이것을 어떻게 덮어서 쓰는지 모르겠다.)
2) Vectorizing Mathematics
SIMD연산, Intel Software Optimization Cookbook을 추천하는데, 이 책은 인텔 홈페이지에서 Intel64 and IA-32 Architectures Optimization 이 책 pdf를 구할 수 있다.
3) Twizzling Rigid-Body Data
rigid body 에 대해 한 알고리즘을 rigid body마다 돌리기 보다는 동시에 그것을 실행시킨다. 즉, simultaneous process을 하도록 data를 twizzle 할필요가 있다는 건데, 함꼐 묵어진다는 것이다: 한 배열에 각 오브젝트에 대한 배열, 그리고 모든 velocities 등등. 이것은 한 번에 4개의 rigid bodies를 보유하는 새로운 자료 구조를 사용하여 얻어질 수 있다.
저자는 개인전으로 이러한 물리엔진을 구현한 적이 없다고 하지만, 자기가 참여한 AI 엔진 개발에서 이 구조를 사용하여 4개의 캐릭터들을 한 번에 업데이트 시켰다고 한다.
4) Grouping Data for Areas Of the Level
메모리 관리가 게임 기술을 최적화하는 중요한 부분이다. 물리에 대해 대부분의 게임 머신에서 이용가능한 많은 메모리가 있지만, 그것의 구성은 slow-downs을 일으킨다.
프로세서는 모든 메모리의 부분들에 대해 동일한 access를 가지지 않는다. 그것들은 main memory의 chunks에 있는 data를 가져오고 빠른 속도의 caches에 유지한다. 캐시 데이터를 접근하는 것은 빠르다. 만약 필요한 데이터가 캐시에 없다면, 그러면 그것은 main memory에서 가져와져야 하고, 이것은 매우 시간을 소모할 수 있다. 다른 기계들은 다른 캐시 사이즈를 가지고, 어떤 것은 cache의 몇 가지 층을 가진다.
끊임 없이 새로운 데이터를 cache로 가져오는 것을 피하기 위해서, 데이터는 함께 묵여져야 한다. 작은 게임 levels에 대해, 모든 물리 데이터는 쉽게 함께 유지될 수 있다. 중간 크기의 게임 levels에 대해, 물리 데이터가 간단히 또 다른 자료 구조에 추가되지 않도록 하는 조치가 취해져야만 한다.
예시 코드는 쉽게 연속된 오브젝트들에 대한 rigid-body data가 메모리에 긴 거리로 늘어나게끔 만들 수 있다. 많은 오브젝트들 사이의 collisions를 resolving할 때, 데이터는 꽤 쉽게 대부분의 resolution steps에서 캐시로 가져와질 필요가 있고, 이것은 재앙의 성능 문제를 발생시킨다.
더 좋은 솔루션은 모든 데이터 집합을 별개의 array에서 유지하는 것이다:
큰 게임 levels에 대해, 이것은 여전히 충분하지 않을 것이다. 이 경우에, rigid bodies의 셋을 ordering하는 것은 가치가 있다. 어떻게 하냐면, game level의 다른 구역에 있는 오브젝트들이 함께 유지되도록 하느 ㄴ것이다. 그렇게, contacts들이 처리될 때, 관련된 오브젝트들은 함께 캐시에 나타날 가능성이 높다. 서로 level을 지나는 오브젝트들 사이의 contacts들은 생성되지 않을 것이고, 그래서 그것들은 그 배열에서 분리될 수 있다.
캐시 미스는 악명높게 진단하기 어렵고, 그것들이ㅡ 만연함은 극적으로 너가 디버깅 코드나 겉보기에 관련없는 조절을 추가할 때 바뀌는 경향이 있다. 좋은 프로파일러가 필수적이다.
An introduction to Physicall Based Modeling : Rigid Body Simulation 2 – Nonpenetration Constraints (David Baraff)
댓글 없음:
댓글 쓰기