Posts Tagged ‘całkowanie verleta’

PostHeaderIcon 5. Łączenie ciał

Możliwe jest łączenie wielu ciał sztywnych używając stawów, bolców, zawiasów, itp. Po prostu niech oba ciała sztywne współdzielą cząstkę – w ten sposób będą połączone bolcem. Gdy współdzielą dwie cząstki, będą połączone zawiasem. Patrz rys. 7.

Możliwe jest również, aby połączyć dwa ciała sztywne za pomocą pręta lub jakiegokolwiek innego rodzaju ograniczenia – aby to zrobić, po prostu dodaj odpowiednią poprawkę do kodu do pętli relaksacyjnej.

Takie podejście pozwala na skonstruowanie kompletnego modelu ludzkiego ciała. Aby uzyskać dodatkowy realizm będą także musiały zostać wprowadzone ograniczenia kątowe. Można to osiągnąć różnymi sposobami. Prostym sposobem jest użycie pręta, który wymusza ograniczenia, jeśli tylko odległość pomiędzy dwoma cząstkami spada poniżej pewnego progu (matematycznie mamy jednostronne ograniczenie odległości (nierówność), |x_2 – x_1| > 100). W efekcie, dwie cząstki nigdy nie będą się znajdować zbyt blisko siebie. Patrz rys. 8.

Inną metodą ograniczania kątów jest spełnienie ograniczenia iloczynu skalarnego:

 (x_2-x_0) \cdot (x_1-x_0) < \alpha

Cząsteczki mogą mieć ograniczony zakres ruchu, na przykład w obrębie określonych płaszczyzn. Po raz kolejny, cząstki z pozycjami nie spełniającymi wyżej wymienionych ograniczeń powinny zostać przeniesione – decydowanie w jaki dokładnie sposób jest nieco bardziej skomplikowane, w związku z ograniczeniami pręta.

W Hitmanie korpusy nie składają się ze sztywnych ciał wzorowanych na czworościanach. Są jeszcze prostsze, gdyż składają się z cząsteczek połączonych przez ograniczający pręt w efekcie tworząc prętowe kształty ( :) Stick Figures). Patrz rys. 9. Położenie i orientacja dla każdej kończyny (wektora i macierzy) są następnie uzyskane dla celów renderingu z pozycji cząstek przy użyciu różnych iloczynów wektorowych i normalizacji wektorów (sprawiając, że kolana i łokcie zginają się naturalnie).

Innymi słowy, każda widziana pojedyncza kończyna nie jest ciałem sztywnym ze zwykłymi 6 stopniami swobody. Oznacza to, że fizycznie rotacja wokół osi kończyny nie jest symulowana. Zamiast tego, system animacji szkieletowych użyty do instalacji wielobocznej siatki postaci jest zmuszony do ustalania pozycji nogi, na przykład takiego, że kolana wydają się zginać naturalnie. Ponieważ rotacja nóg i ramion wokół osi nie obejmuje zasadniczego ruchu z upadającego ludzkiego ciała, to działa w porządku i faktycznie w ogromnym stopniu optymalizuje szybkość.

Ograniczenia kątowe są realizowane w celu egzekwowania ograniczeń ludzkiej anatomii. Prosty test sprawdzania kolizji jest spełniany przez strategiczne wprowadzenie ograniczenia nierówności na odległość, jak wspomniano powyżej, na przykład między dwoma kolanami – upewniając się, że nogi nigdy się nie krzyżują.

W przypadku kolizji z otoczeniem, które składa się z trójkątów, każdy pręt jest reprezentowany jako zamknięty cylinder. Gdzieś w systemie kolizji, a ‚subroutine’ obsługują kolizje między cylindrem a trójkątami. Po wykryciu kolizji, głębokość przenikania i punkty są wyodrębniane, a zderzenie jest traktowane jak pręt ( handled for the offending stick), dokładnie tak, jak opisano na początku sekcji 5.

Oczywiście, wiele dodatkowych ulepszeń było koniecznych, aby uzyskać poprawne wyniki.

PostHeaderIcon 4. Ciała sztywne

Równania opisujące ruch ciał sztywnych odkryto na długo przed wynalezieniem nowoczesnych komputerów. Aby móc wyrazić coś przydatnego w tym czasie, matematycy zmuszeni byli do manipulowania wyrażeniami symbolicznymi. W teorii ciał sztywnych, prowadzi to do przydatnych pojęć i narzędzi, takich jak tensory bezwładności, moment pędu, moment obrotowy, kwaterniony reprezentujące kierunki itp. Jednak, z obecną zdolnością do numerycznego przetwarzania ogromnych ilości danych, korzystne byłoby nawet w niektórych przypadkach rozbić obliczenia na prostsze elementy podczas uruchamiania symulacji. W przypadku trójwymiarowych ciał sztywnych, może to oznaczać modelowanie ciała sztywnego przez cztery cząstki i sześć ograniczeń (podając właściwą ilość stopni swobody, 4*3 – 6 = 6). To ułatwia wiele aspektów i jest to dokładnie tym, czym będziemy się dalej zajmować.

Rozważmy czworościan i umieśćmy cząstki w każdym z jego czterech wierzchołków. Ponadto, dla każdej z sześciu krawędzi czworościanu utworzymy ograniczenie na odległość takie jak pręt ograniczeń omówiony w poprzednim rozdziale. Jest to właściwie wystarczające do symulacji ciała sztywnego. Czworościan może zostać uwolniony wcześniej wewnątrz kostki i całkowanie Verleta pozwoli umieścić go prawidłowo. Funkcja SatisfyConstraints() powinna zadbać o dwie rzeczy: 1) cząstki są przechowywane wewnątrz kostki (jak wcześniej) i 2) sześć ograniczeń dotyczących odległości jest spełnionych. Znowu, można to zrobić za pomocą podejścia relaksacyjnego; trzy lub cztery powtórzenia (ew. z aproksymacją funkcji kwadratowej) powinny wystarczyć.

Powiedzmy wyraźnie, sztywne organy w ogóle nie zachowują się jak czworościany kolizyjnie mądre (choć może to zrobić kinetically). Istnieje też inny problem: obecnie wykrywanie kolizji między ciałem sztywnym i światem zewnętrznym bazuje tylko na wierzchołkach, tzn. jeśli wierzchołek znajduje się na zewnątrz świata jest przerzucany ponownie do wewnątrz. Działa to dobrze tak długo, jak długo wnętrze świata jest wypukłe. Jeśli świat nie jest wypukły, wtedy czworościan i świat na zewnątrz rzeczywiście mogą wzajemnie się przenikać, bez żadnego z wierzchołków czworościanu znajdującego się w nielegalnym regionie (patrz rys. 3, gdzie trójkąt reprezentuje w 2D analogię z czworościanu). Ten problem jest podejmowany dalej.

Zajmijmy się najpierw prostszą wersję problemu. Rozważmy wcześniejszy przykład pręta i załóżmy, że świat zewnętrzny przeżywa mały wstrząs. Pręt może teraz przeniknąć na zewnątrz świata, mimo, że żadna z cząstek pręta nie opuszcza świata (patrz rys. 4). Nie będziemy zagłębiać się w tajniki budowy silnika wykrywania kolizji, gdyż jest to nauka sama w sobie. Zamiast tego założymy, że istnieje już dostępny podsystem, który pozwala nam na wykrywanie kolizji. Ponadto zakładamy, że podsystem może ujawnić nam głębokość przenikania i zidentyfikować punkty przenikania w każdym z dwóch kolidujących obiektów. (Jedna z definicji punktów przenikania i głębokości przenikania brzmi tak: odległość przenikania d_p jest najkrótszą odległością, która uniemożliwia dwóm obiektom przenikania, aby jeden z nich przekształcił jeden z obiektów o odległość d_p w pewnym kierunku. Punkty przenikania są punktami na każdym obiekcie, który dokładnie dotyka innego obiektu po translacji, która miała przed chwilą miejsce).

Spójrz ponownie na rys. 4. Tutaj pręt został przeniesiony po kroku Verleta pod wpływem uderzenia. Silnik kolizji określił dwa punkty przenikania, p i q. Na rys. 4a, p jest rzeczywiście identyczny z pozycją cząstki 1, czyli \textbf{x}_1. Na rys. 4b, p leży między \textbf{x}_1oraz \textbf{x}_2 w odległości ¼ długości pręta od \textbf{x}_1. W obu przypadkach punkt p leży na pręcie, a co za tym idzie może być przedstawiony jako kombinacja liniowa \textbf{x}_1 i \textbf{x}_2, \textbf{p}=c_1 \textbf{x}_1 + c_2 \textbf{x}_2 takie, że c_1+c_2=1. W pierwszym, przypadku {c}_1=1, {c}_2=0, natomiast w drugim \textbf{c}_2=0.75, \textbf{c}_1=0.25. Wartości te mówią nam, o ile powinniśmy przemieścić odpowiednie cząstki.

Aby poprawić nieprawidłową konfigurację pręta, powinien on zostać w jakiś sposób przeniesiony do góry.  Naszym celem jest uniknięcie przenikania dzięki przeniesieniu p na tę samą pozycję co q. Robimy to poprzez dostosowanie pozycji dwóch cząsteczek \textbf{x}_1 oraz \textbf{x}_2 w kierunku wektora między p i q, D = qp.

W pierwszym przypadku, po prostu przerzucamy jak wcześniej \textbf{x}_1 z błędnego obszaru (w kierunku q) i to wszystko (\textbf{x}_2 jest nietknięte). W drugim przypadku, p jest nadal najbliżej \textbf{x}_1 i jest to jeden z powodów, dla których \textbf{x}_1 powinno zostać przeniesione bardziej niż \textbf{x}_2. Ponieważ \textbf{p} = 0.75\textbf{x}_1 + 0.25\textbf{x}_2, wybierzemy do przeniesienia \textbf{x}_1 o 0,75 za każdym razem, kiedy przenosimy \textbf{x}_2 o 0.25. Innymi słowy, nowe miejsca cząstek \textbf{x'}_1 oraz \textbf{x'}_2 są dane wyrażeniami:

 {x_1}' = x_1+0.75\lambda \cdot \Delta
 {x_2}'= x_2+0.25\lambda \cdot \Delta

gdzie l jest pewną nieznaną wartością. Nową pozycją p po przeprowadzce obu cząstek jest
\textbf{p'}=c_1 \textbf{x'}_1 + c_2 \textbf{x'}_2.

Przypomnijmy, że chcemy, aby p’ = q, tzn. należy wybrać l dokładnie takie, żeby p’ zbiegło się z q. Ponieważ przenosimy cząstki tylko w kierunku D, również p przesuwa się w kierunku D i w związku z tym rozwiązanie równania p’ = q można znaleźć poprzez obliczenie

 p' \cdot \Delta = q \cdot \Delta

dla l. Rozwijając lewą stronę:

 \begin{tabular}{ r l } \( p' \cdot \Delta \) & \(= (0.75 \cdot {x_1}'+0.25 \cdot {x_2}') \cdot \Delta \) \\ & \(= (0.75 \cdot ( x_1 + 0.75\lambda \cdot \Delta) + 0.25 \cdot ( x_2+0.25\lambda \cdot \Delta)) \cdot \Delta \) \\ & \(= (0.75 \cdot x_1 + 0.25 \cdot x_2) \cdot \Delta + \lambda (0.75^2 +0.25^2) \cdot \Delta^2 \) \\ & \(= p \cdot \Delta + \lambda (0.75^2 + 0.25^2) \cdot \Delta^2 \) \end{tabular}

które z prawą stroną równania (**) daje:

 \lambda = \frac{(\textbf{q}-\textbf{p}) \cdot \Delta} { (0.75^2+0.25^2) \cdot \Delta^2}

Łącząc l z (*) otrzymujemy nowe pozycje cząstek, dla których p’ zbiegnie się z q.

Rysunek 5 obrazuje sytuację po ruszeniu cząsteczek. Nie mamy żadnego przenikającego obiektu, ale teraz długość ograniczającego pręta została naruszona. Aby rozwiązać ten problem, robimy kolejną iterację pętli (lub kilka).

Powyższa strategia działa również dla czworościanu w zupełnie analogiczny sposób. Najpierw znajdowane są punkty przenikające p i q (mogą to być również punkty leżące wewnątrz trójkąta), a **p** jest wyrażone jako liniowa kombinacja czterech cząstek \textbf{p}=c_1 \textbf{x}_1 + c_2 \textbf{x}_2 + c_3 \textbf{x}_3 + c_4 \textbf{x}_4 takich, że c_1 + c_2 + c_3 + c_4 = 1 (wymaga to rozwiązania niewielkiego układu równań liniowych). Po znalezieniu D = qp, obliczamy wartość

 \lambda = \frac{(\textbf{q}-\textbf{p}) \cdot \Delta} { (c_1^2+c_2^2+c_3^2+c_4^2) \cdot \Delta^2}

i nowe pozycje są następnie dane przez
 {x_1}' = {x_1} + c_1 \lambda \cdot \Delta
 {x_2}' = {x_2} + c_2 \lambda \cdot \Delta
 {x_3}' = {x_3}+c_3 \lambda \cdot \Delta
 {x_4}' = {x_4}+c_4 \lambda \cdot \Delta

Tutaj mamy kolizję ciała sztywnego z nieruchomym światem. Powyższa metoda jest ogólną do obsługi kolizji kilku sztywnych ciał. Kolizje są przetwarzane na jedną parę ciał naraz. Zamiast poruszać tylko p, w tym przypadku zarówno p i q są przenoszone do siebie.

Znowu, po skorygowaniu pozycji cząstki w taki sposób, że spełniają one ograniczenia nieprzenikania, należy uwzględnić sześć ograniczeń odległości, które tworzą ciała sztywne, itd. W tej metodzie, czworościan może być nawet zanurzony wewnątrz innego obiektu, który może być stosowany zamiast samego czworościanu do obsługi kolizji. Na rys. 6, czworościan jest osadzony wewnątrz kostki.

Po pierwsze, kostka musi być umocowana do czworościanu w jakiś sposób. Jednym z rozwiązań byłoby wybranie systemu masowej punkcie środkowym 0.25*(x_1 + x_2 + x_3 + x_4) jako pozycję kostki, a następnie dziedziczyć orientacyjną macierz (derive an orientation matrix) poprzez analizę obecnych pozycji cząsteczek. Gdy kolizja/przenikanie zostaną wychwycona, punkt kolizyjny p (który w tym przypadku zostanie położony na kostce) jest traktowana dokładnie tak samo, jak powyżej i pozycje cząsteczek są stosownie aktualizowane. Jako optymalizacja, możliwe jest wcześniejsze przeliczenie wartości c_1 - c_4 dla wszystkich wierzchołków z kostki. Jeśli punktem przenikania p jest wierzchołek, wartości c_1 - c_4 mogą być użyte bezpośrednio. W innym wypadku, p leży na wnętrzu powierzchni trójkąta lub jednego z jego krawędzi i wartości c_1 - c_4 mogą następnie być interpolowane z przeliczonych wcześniej wartości odpowiednich wierzchołków trójkąta.

Zazwyczaj 3 do 4 powtórzeń relaksacji jest wystarczające. Ciała nie będą się zachowywać tak, jakby były całkowicie sztywne, ponieważ iteracje relaksacyjne są zatrzymane przedwcześnie. Jest to przede wszystkim użyteczna funkcja, jako że nie istnieje coś takiego jak doskonale sztywne ciało – zwłaszcza ludzkie. To również sprawia, że system stał się bardziej stabilny.

Przez przemieszczenie cząsteczek, które tworzą czworościan, właściwości fizyczne również mogą być zmienione (matematycznie, bezwładność się zmienia, kiedy pozycje i masy cząstek są zmieniane).

Możliwe są inne ustawienia cząstek niż czworościan, np. umieszczenie cząsteczek w modelu układu współrzędnych, tzn. w (0,0,0), (1,0,0), (0,1,0), (0,0,1). Niech a, b, c będą wektorami od cząstki 1 do cząstek odpowiednio 2, 3 i 4. Ograniczenie pozycji cząstek przez przyjęcie, aby wektory a, b, c miały długość 1, a kąt między każdą z trzech par wektorów wynosił 90 stopni (odpowiadający iloczyn skalarny powinien wynosić zero). (Zauważ, że ten znów daje cztery cząstki i sześć ograniczeń).

PostHeaderIcon 2. Całkowanie Verleta

Sercem symulacji jest system cząsteczek. Zazwyczaj w implementacji systemu cząstek każda cząstka ma dwie główne zmienne: x – pozycja i v – jej prędkość. Następnie w ‚czasokrokowej’ pętli, nowa pozycja x’ i prędkość v’ są obliczane według:

 \textbf{x}' = \textbf{x} + \textbf{v} \cdot \Delta t
 \textbf{v}' = \textbf{v} + \textbf{a} \cdot \Delta t

gdzie dt jest krokiem czasu, a – przyspieszeniem obliczanym z wykorzystaniem prawa Newtona  \textbf{F} = m \cdot \textbf{a} (gdzie F to wypadkowa sił działających na cząstki). Jest to prosta integracja Eulera.

W naszym przypadku jednak wybieramy schemat w mniejszym stopniu oparty na prędkości: zamiast przechowywać pozycję i prędkość każdej cząstki, będziemy przechowywać jej aktualną pozycję x oraz wcześniejszą \textbf{x}^*. Utrzymując stały czasowy krok, pętla aktualizacyjna wygląda wtedy następująco:

 \textbf{x}' = 2\textbf{x} - \textbf{x}^* + \textbf{a} \cdot \Delta t^2
 \textbf{x}^* = \textbf{x}

Nazywamy to całkowaniem Verleta (patrz [((bibcite Verlet))]); jest ono intensywnie wykorzystywane podczas symulowania dynamiki molekularnej. Metoda jest dość stabilna, gdyż prędkość jest domyślnie podana, a co za tym idzie trudniej jest prędkości i pozycji być niezsynchronizowanymi. (Na marginesie warto zauważyć, że podobne podejście wykorzystuje znany efekt tworzenia fal w wodzie). Metoda działa z uwagi na fakt, że  2x - x^* = x + (x - x^*) oraz  x - x^* są przybliżeniem aktualnej prędkości (tak naprawdę jest to droga przebyta od ostatniego kroku). Metoda ta nie zawsze jest wyjątkowo dokładna (energia może np. opuścić system), ale za to szybka i stabilna. Zmieniając aktualizację reguły np. na taką:

 \textbf{x}' = 1,99\textbf{x} - 0,99\textbf{x}^* + \textbf{a} \cdot \Delta t^2

możemy wprowadzić niewielką ilość oporu do systemu.

Pod koniec każdego etapu aktualna pozycja x każdej cząstki zostaje przechowana w odpowiedniej zmiennej \textbf{x}^*. Zauważ, że gdy manipulujemy wieloma cząstkami, przydatną optymalizacją może być po prostu zamiana tablicy wskaźników.

Wynikowy kod będzie wyglądać mniej więcej tak (klasa Vector3 powinna zawierać odpowiednie funkcje i przeciążenia operatorów manipulacji wektorów):

// przykładowy kod symulacji fizyki
class ParticleSystem
{
Vector3   m_x[NUM_PARTICLES];            // aktualne pozycje
Vector3   m_oldx[NUM_PARTICLES];        // poprzednie pozycje
Vector3   m_a[NUM_PARTICLES];            // obliczone siły
Vector3   m_vGravity;                           // grawitacja
float   m_fTimeStep;

public:
void   TimeStep();

private:
void   Verlet();
void   SatisfyConstraints();
void   AccumulateForces();

// konstruktory, inicjalizacje itp. pominięte
};

// krok całkowania Verleta
void ParticleSystem::Verlet()
{
for (int i = 0; i < NUM_PARTICLES; i++)
{
Vector3& x = m_x[i];
Vector3 temp = x;
Vector3& oldx = m_oldx[i];
Vector3& a = m_a[i];
x += x - oldx + a*fTimeStep*fTimeStep;
oldx = temp;
}
}

// funkcja ta powinna obliczać siły dla każdej cząstki
void ParticleSystem::AccumulateForces()
{
// wszystkie cząstki podlegają grawitacji
for (int i = 0; i < NUM_PARTICLES; i++)
m_a[i] = m_vGravity;
}

// tu uwzględniamy ograniczenia
void ParticleSystem::SatisfyConstraints()
{
// ta funkcja jest na razie nieistotna
}

void ParticleSystem::TimeStep()
{
AccumulateForces();
Verlet();
SatisfyConstraints();
}

Powyższy kod został napisany dla przejrzystości, a nie szybkości działania. Jedną z optymalizacji byłoby użycie tablic typu float zamiast Vector3 dla reprezentacji stanów. Mogłoby to również ułatwić implementację na procesorach wektorowych.

Zaprezentowana metoda prawdopodobnie nie wygląda jeszcze bardzo przełomowo. Jednakże korzyści powinny być wyraźne, gdy zaczniemy używać ograniczeń i przełączymy się na ciała sztywne. Zauważymy wtedy, jak powyższa metoda prowadzi do zwiększenia stabilności i zmniejszenia ilości obliczeń w porównaniu do innych metod.

Spróbuj np. ustawić \textbf{a} = (0,0,1) i użyć jako warunku początkowego \textbf{x} = (0,0,0) i \textbf{x}^* = (0,0,0), a następnie wykonać ręcznie kilka powtórzeń, aby zobaczyć, co się dzieje.

PostHeaderIcon Całkowanie Verleta

Spis treści:

  1. Wprowadzenie
  2. Całkowanie Verleta
  3. Kolizje i ich obsługa przez rzutowanie
  4. Ciała sztywne
  5. Łączenie ciał
  6. Komentarze i wnioski