Simulace písečné bouře

České vysoké učení technické v Praze Fakulta elektrotechnická Katedra počítačové grafiky a interakce

Diplomová práce

Simulace písečné bouře Bc. David Eršil

Vedoucí práce: Ing. Jaroslav Sloup

Studijní program: Elektrotechnika a informatika (magisterský), strukturovaný Obor: Výpočetní technika leden 2011

ii

Poděkování Na tomto místě bych rád poděkoval vedoucímu práce, Ing. Jaroslavu Sloupovi, za věcné připomínky. Zároveň děkuji také rodině a přítelkyni, bez jejichž tolerance a podpory by tvorba této práce byla mnohem náročnější. Poděkování patří v neposlední řadě také Bc. Lukáši Oreškovi za jeho příspěvek k výsledné aplikaci v podobě vykreslovacího modulu. iii

iv

Prohlášení Prohlašuji, že jsem svoji diplomovou práci vypracoval samostatně a použil jsem pouze podklady uvedené v přiloženém seznamu. Nemám závažný důvod proti užití tohoto školního díla ve smyslu §60 Zákona č. 121/2000 Sb., o právu autorském, o právech souvisejících s právem autorským a o změně některých zákonů (autorský zákon).

V Lužné dne 2.1. 2011

................................................................... v

vi

Abstract The goal of this work is an implementation of a simulator of evolution and move of still more common natural phenomenon – sand storm – based on fluid model and its physical description by Navier-Stokes equations. The result of this work is a physical module working with adaptive space division, ensuring interaction between storm and surrounding scene and which gives raise to a full simulation application after merging with graphic module.

Abstrakt Cílem této práce je implementace simulátoru vývoje a pohybu stále častějšího přírodního jevu – písečné bouře – na základě fluidního modelu a jeho fyzikálního popisu pomocí Navier-Stokesových rovnic. Výsledkem této práce je fyzikální modul, který pracuje s adaptivním dělením prostoru, zajišťuje interakci bouře s okolní scénou a po sloučení s grafickým modulem vytváří plnohodnotnou simulační aplikaci.

vii

viii

Obsah Kapitola 1 – Úvod...............................................................................1 1.1 Cíl.....................................................................................................................1 1.2 Požadavky.......................................................................................................2 1.3 Písečná bouře.................................................................................................3 1.3.1 Vznik bouře.............................................................................................3 1.3.2 Průběh bouře..........................................................................................7 1.3.3 Mimozemské písečné bouře...................................................................8 Kapitola 2 – Analýza.........................................................................11 2.1 Předchozí práce.............................................................................................11 2.1.1 Vědecké publikace.................................................................................11 2.2 Rozbor požadavků........................................................................................17 2.2.1 Fyzikální vztahy – Navier-Stokesovy rovnice......................................17 2.2.1.1 Vliv externích sil............................................................................19 2.2.1.2 Advekce..........................................................................................19 2.2.1.3 Difuze (rozptyl)..............................................................................21 2.2.1.4 Projekce.........................................................................................23 2.2.1.5 Počáteční a okrajové podmínky....................................................24 2.2.2 Adaptivní dělení...................................................................................25 2.2.3 Simulace pískového podloží................................................................28 2.2.3.1 Vzedmutí............................................................................................29 2.2.3.2 Unášení.........................................................................................29 2.2.3.3 Splaz..............................................................................................29 2.2.4 Statická scéna.......................................................................................30 Kapitola 3 – Návrh řešení.................................................................31 3.1 Volba programovacího jazyka......................................................................31 3.2 Návrh řešení pro hledání sousedů .............................................................32 3.2.1 Sousedé v oktalovém stromu...............................................................32 3.2.2 Sousedé v uniformní mřížce................................................................33 3.2 Přístup k řešení Navier-Stokesových rovnic...............................................35 3.2.1 Návrh řešení advekčního kroku...........................................................35 3.2.2 Návrh řešení projekčního kroku.........................................................36 ix

3.3 Návrh řešení interakce se statickými objekty ve scéně..............................40 3.3.1 Interakce objektů a proudového pole..................................................40 3.3.2 Interakce objektů a unášeného písku.................................................42 3.4 Aktualizace oktalového stromu...................................................................43 Kapitola 4 – Implementace..............................................................45 4.1 Implementace fyzikálního modulu.............................................................45 4.1.1 Struktura modulu – třída Solver..........................................................45 4.1.2 Interakce proudového pole s pískovým podložím...............................52 4.2 Datové struktury simulátoru.......................................................................54 4.2.1 Třída COctree - oktalový strom............................................................55 4.2.2 Struktura SGrid – uniformní mřížka..................................................63 4.2.3 Struktura SVoxel – výpočetní buňka..................................................66 4.3 Plynulost běhu aplikace..............................................................................69 Kapitola 5 – Ladění..........................................................................71 5.1 Podpůrné aplikace........................................................................................71 5.1.1 MemoryValidator..................................................................................71 5.1.2 Very Sleepy............................................................................................74 5.2 Zpracování výsledků testování....................................................................76 Kapitola 6 – Výsledky.......................................................................79 6.1 Adaptivní dělení scény.................................................................................79 6.2 Proudové pole..............................................................................................82 6.3 Bouřkové mračno........................................................................................85 6.4 Podloží.........................................................................................................90 Kapitola 7 – Závěr............................................................................95 Použitá literatura a zdroje...............................................................97 Příloha.............................................................................................99 A – Interpolace v advekci..................................................................................99 B – Uživatelská příručka...................................................................................99 Instalace........................................................................................................99 Spuštění.......................................................................................................100 Ovládání.......................................................................................................102 Odinstalace..................................................................................................103 C – Obsah přiloženého datového média.........................................................103 x

xi

xii

Seznam obrázků Obr.1: Rozdělení úkolů pro jednotlivé části simulátoru....................................................2 Obr.2: Písečná bouře mířící k Chartúmu..........................................................................3 Obr.3: Mapa extrémně vyprahlých, vyprahlých a středně suchých oblastí ve světě.........4 Obr.4: Průběh mezikontinentálního přenosu pouštního písku..........................................6 Obr.5: Mechanismy transportu písku................................................................................7 Obr.6: Přechod písečné bouře přes americkou dálnici I-5...............................................8 Obr.7: Vývoj bouřkového mračna na Marsu v průběhu čtyř týdnů...................................9 Obr.8: Pohled na planetu Mars před a po písečné bouři................................................10 Obr.9: Výsledky simulace horkého vířivého plynu..........................................................12 Obr.10: Výsledky Stamova implicitního fluidního simulátoru........................................13 Obr.11: Porovnání simulované a reálné scény................................................................13 Obr.12: Scéna simulátoru písečné bouře........................................................................14 Obr.13: Výstup fluidního simulátoru s adaptivním prostorový dělením..........................14 Obr.14: Simulace tekutého písku.....................................................................................15 Obr.15: Simulace částicově reprezentované tekutiny za pomoci zjemňovaného vzorkování.......................................................................................................................15 Obr.16: Simulace vzduchové bubliny stoupající vodním objemem s překážkou.............16 Obr.17: Porovnání výsledku simulace kouře s využitím uniformní mřížky a adaptivního dělení...............................................................................................................................17 Obr.18: Schéma operací..................................................................................................19 Obr.19: Grafická znázornění principu zpětného advekčního trasování..........................21 Obr.20: Difuzní schéma...................................................................................................21 Obr.21: Porovnání uniformní mřížky a adaptivního dělení ve 2D (quadtree)................26 Obr.22: Porovnání uniformního a adaptivního dělení v prostoru (octree).....................26 Obr.23: Uzel oktalového stromu a jeho vnitřní struktura...............................................27 Obr.24: Příklad terénu vymodelovaného výškovým polem.............................................28 Obr.25: Schematické znázornění principu unášení.........................................................29 Obr.26: Schéma spolupráce mezi jednotlivými částmi celého výsledného simulátoru...31 Obr.27: Schéma jednotlivých kroků simulátoru (včetně grafické části)..........................31 Obr.28: Indexace potomků v COctree.............................................................................32 xiii

Obr.29: Binární strom pro ilustraci prapředků...............................................................33 Obr.30: Indexace voxelů uvnitř mřížky pro 8 buněk v každé ose....................................34 Obr.31: Znázornění interpolačních koeficientů..............................................................36 Obr.32: Ilustrace předposledních uzlů binárního stromu...............................................37 Obr.33: Ilustrace principu přebírání rychlostí z nižších úrovní......................................37 Obr.34: Ilustrace „podrozdělení“ voxelů doplněná o indexy buněk vyšší úrovně..........38 Obr.35: Rozdílná indexace mřížek a uzlů oktalového stromu.........................................38 Obr.36: Ukázka přizpůsobení adaptivního dělení objektům ve scéně v souvislosti s obtékáním........................................................................................................................41 Obr.37: Reálné obtékání objektu proudovým polem a situace při backtrackingu..........41 Obr.38: Diagram zjednodušených tříd............................................................................54 Obr.39: Úvodní obrazovka aplikace MemoryValidator..................................................72 Obr.40: Výřez okna aplikace MemoryValidator v průběhu testování.............................73 Obr.41: Hlavní okno aplikace Very Sleepy......................................................................74 Obr.42: Výsledky profilování simulátoru písečné bouře.................................................75 Obr.43: Graf závislosti doby běhu fyzikálního modulu na složitosti scény.....................76 Obr.44: Graf vztahu mezi alokovanou pamětí a počtem voxelů......................................77 Obr.45: Přizpůsobení oktalového stromu objektu ve scéně.............................................79 Obr.46: Přizpůsobení adaptivního dělení většímu počtu objektů ve scéně ....................80 Obr.47: Adaptivní dělení složitější scény........................................................................80 Obr.48: Průběh podrozdělování scény v místech výskytu bouřkového mračna..............81 Obr.49: Vliv maximální hloubky dělení stromu na obtékání modelů..............................82 Obr.50: Proudové čáry....................................................................................................83 Obr.51: Chybový stav, při kterém proudové čáry protínají objekty.................................84 Obr.52: Demonstrace zastavení proudového pole v okolí objektu..................................84 Obr.53: Vývoj bouřkového mračna.................................................................................85 Obr.54: Průběh bouře ve scéně s jedním větším modelem..............................................86 Obr.55: Testovací scéna č.1.............................................................................................87 Obr.56: Prostřední část Obr.54 v plné velikosti..............................................................87 Obr.57: Osmnáctý krok simulace první scény v plné velikosti........................................88 Obr.58: Testovací scéna č.2.............................................................................................88 Obr.59: Prostřední část Obr.58 v plné velikosti..............................................................89

xiv

Obr.60: Detailní pohled na obtékání mračna okolo dřevěných klád..............................89 Obr.61: Testovací scéna č.3.............................................................................................89 Obr.62: Prostřední část Obr.61 v plné velikosti..............................................................90 Obr.63: Terén bez vyhlazování a terén zpracovaný metodou creep................................91 Obr.64: Detail terénu a vliv počtu iterací na výsledky vyhlazovací metody creep..........91 Obr.65: Příklady vzedmutí písku v závislosti na proudovém poli...................................92 Obr.66: Výsledná scéna po spojení fyzikálního a grafického modulu.............................93

xv

xvi

KAPITOLA 1 – ÚVOD

1

Kapitola 1 – Úvod Využití počítačové grafiky v oblasti vizuálních simulací nejrůznějších přírodních jevů se s rostoucí výpočetní silou hardwaru stává stále populárnější (např. ve filmovém průmyslu) a tyto simulace zaznamenávají velice dobré výsledky. I přes velké množství simulátorů však na téma písečných bouří zatím mnoho prací vytvořeno nebylo. Existence takového simulátoru je ale poměrně důležitá, neboť písečné bouře mají zásadní vliv na podobu prostředí, v němž se vyskytují, stejně jako na život a zdraví jeho obyvatelů. Jelikož jsem se během magisterské etapy studia na katedře výpočetní techniky ČVUT FEL již jednou podílel na zajímavé semestrální práci, implementující jednodušší simulátor písečné bouře, rozhodl jsem se vrátit se k tomuto tématu ve své diplomové práci a implementovat simulátor nový a lépe zpracovaný. V jednotlivých částech této kapitoly bude proveden základní popis písečných bouří. Ve druhé kapitole budou představeny Navier-Stokesovy rovnice a metoda jejich řešení ve vztahu k simulaci vytvářené v této práci. Třetí kapitola bude zaměřena na návrh řešení simulace a na základní popis použitých postupů. V kapitole čtvrté dojde k podrobnějšímu popisu použitých algoritmů a jejich provázání. Následuje ladění vytvořeného simulátoru a kapitola s výsledky.

1.1 Cíl Tato diplomová práce se zabývá návrhem a implementací prototypu softwarové aplikace, s jejíž pomocí bude možné provádět simulaci písečné bouře, založenou na fyzikálním modelu. S doplňujícím grafickým modulem, který je předmětem souběžně vytvářené implementační diplomové práce Vizualizace písečných bouří Bc. Lukáše Orešky [1], bude následně možné sledovat průběh bouře a interakci s různými podobami scény, v níž působí.

2


Obr.1: Rozdělení úkolů pro jednotlivé části simulátoru

1.2 Požadavky ●

Simulace musí vycházet z odpovídajících fyzikálních vztahů – řešení Navier-Stokesových rovnic.

●

Předmětem simulace má být nejen proudové pole, ale zároveň také pískové podloží.

●

Simulace musí využívat adaptivní dělení prostoru.

●

Simulace bude brát v potaz případné překážky ve scéně. Ta bude složena z písečného podloží a 3D modelů, reprezentovaných trojúhelníkovou sítí.

●

Pohyb a změna tvaru se týká pouze mračna a podloží. Ostatní objekty ve scéně jsou statické.

●

Simulace má být implementována v jazyce C, případně C++ a využívat grafickou knihovnu OpenGL.

●

Simulace by měla probíhat co nejplynuleji, avšak reálný čas není podmínkou.


3

1.3 Písečná bouře Australský duststorm, simúm v saharské oblasti a na Arabském poloostrově, Andhi v Indii či habúb v Súdánu. Různá lokální pojmenování pro jeden a tentýž přírodní jev – písečnou bouři. Husté mračno písečných zrn, rostoucí v extrémních případech až do výše větší, než 1 kilometr nad zemským povrchem, s sebou přes rozsáhlá území rychlostí až 130 kilometrů za hodinu unáší tisíce (ve výjimečných případech až miliony) tun jemného písku.

Obr.2: Písečná bouře mířící k Chartúmu (Súdán), 1. května 2007 (převzato z [O2])

1.3.1 Vznik bouře Písečné bouře, vznikají v důsledku desertifikace, ve vyprahlých písčitých oblastech po celém světě. Takové oblasti pokrývají zhruba třetinu zemské souše, avšak v důsledku lidské činnosti se jejich plocha neustále zvětšuje. V závislosti na poměru spadlých srážek a středního ročního výparu z rostlin a půdy se dělí na tři druhy [2], jejichž výčet je uveden v Tabulce 1.

4

KAPITOLA 1 – ÚVOD VÝSLEDNÝ POMĚR

DRUH OBLASTI

< 0.03

EXTRÉMNĚ VYPRAHLÁ (HYPERARID)

0.03 ~ 0.2

VYPRAHLÁ (ARID)

0.2 ~ 0.5

STŘEDNĚ SUCHÁ (SEMIARID)

Tabulka 1: Rozdělení oblastí podle poměru srážek a výparů

Obr.3: Mapa extrémně vyprahlých, vyprahlých a středně suchých oblastí ve světě. Převzato z [O3]

Studie zemské kůry na dně oceánů dokazují, že písečné bouře se na naší planetě vyskytovaly již v pravěku. Jejich četnost ale strmě roste. Díky záznamům ze staré Číny, sahajících až do období okolo roku 1100 před Kristem, víme, že zatímco tehdy docházelo k bouřím jednou za několik století, během 50. let minulého století se jich na území Číny zformovalo pět. Během osmdesátých let počet výskytů vzrostl na 19 a jen v průběhu roku 2000 jich zde bylo zaznamenáno celkem 12 [3]. Navzdory své bohaté historii se ale písečné bouře staly předmětem intenzivního zkoumání až v 50. letech 20. století. Díky němu a nejrůznějším vědeckým pokusům ve větrných tunelech, jsou v dnešní době


5

fyzikální procesy, odpovědné za vznik písečných bouří, již velmi dobře prozkoumány a byly rozděleny do tří fází [4] – vzedmutí (suspension), unášení (saltation) a splaz (creep), viz Obr.5. Základem písečné bouře je vítr (respektive nestabilní vzdušné proudy) a suchý, sypký povrch, tvořený pískem a prachem. V klidovém stavu jsou zrnka písku udržována na podloží gravitační silou a vzájemnou vazbou. S rostoucí silou větru se však zrnka dávají nejprve do vibračního pohybu, aby se pak následně uvolnila a vznesla do ovzduší. Tento proces je označován pojmem vzedmutí a působí na nejmenší zrnka (s průměrem pod 0,1 mm). Pro vzedmutí těchto částic je zapotřebí vítr o rychlosti zhruba 3,5-4 m/s. Právě obrovské množství mikroskopických a lehkých zrnek je zodpovědné za tvorbu celého neprůhledného bouřkového mraku. Meteorologický úřad Spojených států například vyhlašuje stav písečné bouře při snížení viditelnosti na hodnotu 5 8

míle (0,5 - 1 km). Při poklesu viditelnosti pod

5 16

5 16

až

míle (0,5 km) je již bouře

vyhodnocena jako velmi silná. Tyto částice pak mohou být vyneseny až do výše několika tisíc metrů a v této hladině unášeny přes velké vzdálenosti napříč kontinenty a oceány. Například není výjimkou, že částečky ze saharské písečné bouře překonají za pomoci tropického východního proudění celý Atlantský oceán a spadnou hluboko ve vnitrozemí Jižní Ameriky. Podle výzkumů ([5]) takto každoročně skončí v Amazonské nížině okolo 13 milionů tun saharského písku, což představuje zhruba 190 kg na území o ploše jednoho hektaru. Sekvence radarových snímků Severoamerického úřadu pro letectví a kosmonautiku (NASA) (Obr.4) zachycuje pohyb písečných zrn z pouští Gobi a Sahara v období od 7. do 13. dubna 2001. Takový pohyb má ovšem za následek ovlivnění životního prostředí v cílových oblastech. Ale nemusí jít vždy nutně o vliv negativní. Zatímco poprašek původně saharských písečných zrn na vrcholcích alpských velikánů vyvolává předčasné tání sněhu, do údolí Amazonky přináší fosfor a zúrodňuje zdejší půdu. Písek, který svoji pouť napříč planetou ukončí na hladině oceánu, podporuje tvorbu fytoplanktonu – jednobuněčných organismů, např. sinic a řas – který absorbuje vzdušný oxid uhličitý a v neposlední řadě také

6


slouží jako potrava pro celou škálu mořských živočichů. Ty nejmenší částice, s poloměrem v jednotkách mikrometrů, však mohou zůstat v horních vrstvách atmosféry i několik let.

Obr.4: Průběh mezikontinentálního přenosu pouštního písku (převzato z [O4])


7

Větší částice (s průměrem v rozmezí 0,1 až 0,5 mm), které tvoří 50 až 80% celkového množství transportovaného písku, jsou součástí procesu unášení. Rovněž tato zrna jsou větrným polem zdvihnuta, avšak kvůli své váze zřídkakdy do výšky přesahující 30 centimetrů, aby po několika metrech opět spadla, čímž rozpohybují částice v okolí dopadu. To má za následek doslova řetězovou reakci, při které se dává do pohybu stále větší masa písku. Konečně největší písečná zrna, s průměrem nad 0,5 mm, se účastní třetího mechanismu – splazu. Kvůli své váze se již nedostávají do ovzduší a jsou větrným prouděním unášeny po povrchu. V této fázi dochází k přeskupování písečných dun. Písečné částice o průměru nad 2 milimetry se dávají do pohybu již jen velmi zřídka.

Obr.5: Mechanismy transportu písku (převzato z [O5])

1.3.2 Průběh bouře Podobně jako v případě jiného fenoménu – tornáda – je i v okolí písečného mraku velmi slabý, ba dokonce nulový vítr, avšak tlak vzduchu výrazně klesá, zatímco jeho teplota může prudce vzrůst až o deset stupňů. Narozdíl od tornáda však může písečná bouře probíhat pod zcela jasnou oblohou a trvat obvykle několik hodin. S příchodem bouře do místa pozorovatele se situace výrazně mění. Jako první přichází na řadu extrémně silný vítr, který zdvihne písek a prach. Zároveň dochází k takřka skokovému zvýšení tlaku, doprovázenému naopak rychlým poklesem teploty. Měření, které

8


bylo provedeno 22. dubna 1977 při písečné bouři v Zhangye v čínské provincii Západní Gansu, ukázalo, že zde v průběhu deseti minut vzrostl tlak o 2,8 hPa a teplota vzduchu klesla o 6,8 °C. Směr větru se změnil z východního na západní až severozápadní, přičemž se jeho střední rychlost zvýšila o 20 m/s a ve své nejvyšší hodnotě sahala až ke 30 m/s. Koncentrace písečných zrn a prachových částic v ovzduší více než 40-násobně přesáhla limity pro extrémní znečištění. Hodnota na otevřeném prostranství dosahovala více než 1 gramu na m3, uvnitř budov byla naměřena koncentrace až 80 mg/m3. Bouře, která 20. května 1976 zasáhla letiště v Novém Dillí, s sebou během pouhých dvou minut přinesla snížení viditelnosti ze 4000 metrů na 280, teplota vzduchu klesla při rychlosti 80 km/h z 35 °C na 28 °C a jeho vlhkost se zvýšila z 30 na 70%.

Obr.6: Přechod písečné bouře přes americkou dálnici I-5, spojující Kanadu a Mexiko (převzato z [O6])

1.3.3 Mimozemské písečné bouře Písečné bouře nejsou pouze výsadou naší planety. Mohou se vytvořit všude tam, kde k tomu existují potřebné podmínky – tedy přítomnost atmosféry a rozsáhlých vyprahlých oblastí pokrytých pískem, prachem či jiným sypkým materiálem s částicemi o dostatečně malém průměru. Zároveň jsou často schopny dosáhnout takových rozměrů a intenzit, které by na planetě Zemi měly


9

na svědomí totální zkázu. V červnu 2001 poskytla sonda Mars Global Surveyor společně s Hubbleovým teleskopem lidskému oku možnost spatřit úchvatný jev – písečnou bouři na Marsu, jejíž mračno během čtrnácti dnů pohltila celou planetu na dobu tří měsíců (viz Obr.7).

Obr.7: Vývoj bouřkového mračna (červená barva) na Marsu v průběhu čtyř týdnů (převzato z [O7])

Příčinou této neobvykle silné bouře bylo množství drobného prachu, poletujícího atmosférou [6]. Ten způsoboval lom slunečního záření, což mělo za následek postupné zahřívání vyšších vrstev atmosféry o několik desítek stupňů Celsia nad běžný stav a zároveň zastínění povrchu a tudíž jeho chladnutí. Rozdíl v teplotách a tlacích nejprve vytvořil silné atmosférické proudění nad pánví Hellas Planitia, 9 kilometrů hlubokém a 2300 km širokém kráteru, jehož přítomnost narušuje atmosférickou cirkulaci a stává se tak častým místem vzniku písečných bouří. Proudění na jižní polokouli zajistilo propagaci jevu východním směrem a zároveň také přes rovník, což byl úkaz, který zatím nebyl

10


nikdy v historii zdokumentován. Postupně díky členitosti povrchu docházelo k vytváření dalších ohnisek nových bouří, které se nakonec spojily v jeden obrovský celek, rozprostírající se napříč celou planetou, od jižního pólu k severnímu (Obr.8). Jak je ale možné, že se bouře dostala i nad oba póly, je pro vědce stále hádankou, neboť planetární póly bývají tvořeny atmosférickými kapsami, zabraňujícími průniku teplého proudění z rovníkových oblastí.

Obr.8: Pohled na planetu Mars před a po písečné bouři vlevo – obvyklý pohled na Mars (v JV části patrné víry začínající bouře nad planinou Hellas) vpravo – planeta pohlcená mračnem písku (převzato z [O8])

KAPITOLA 2 – ANALÝZA

11

Kapitola 2 – Analýza V této kapitole uvedu přehled významných prací, zabývajících se fluidní simulací. Následně provedu rozbor Navier-Stokesových rovnic (NSR) a uvedu přístup k jejich řešení v souvislosti s prováděnou simulací. Kapitolu pak zakončí podrobnější rozbor požadavků na aplikaci a přístup k jejich naplnění.

2.1 Předchozí práce V oblasti fluidní dynamiky a její simulace pomocí výpočetní techniky byly vypracovány stovky prací, ať již v podobě vědeckých článků či rovnou celých knih. Avšak prací, orientovaných na fluidní dynamiku ve vztahu k písečným bouřím je o poznání méně. Takové simulátory, používané v meteorologii a fyzikálních experimentech, jsou zaměřeny převážně na zpracování velkého množství reálných dat a snaží se o maximálně přesný výpočet a modelování budoucího vývoje a chování bouří. Jejich komplikovanost však brání jejich použití v praktických grafických aplikacích.

2.1.1 Vědecké publikace Ačkoliv ve fyzice představují Navier-Stokesovy rovnice známý problém, simulace založené na jejich řešení jsou výpočetně velmi náročné a začaly se rozšiřovat až s nástupem výpočetní techniky. První metody pro řešení fluidní simulace byly vyvinuty ve 30. letech 20. století. Při výpočtech toku dvojrozměrného proudového pole v okolí profilu křídla letounu (získaného z kruhu pomocí Žukovského transformace [7]) se využívaly linearizované potenciální rovnice (tj. Navier-Stokesovy rovnice zjednodušené odebráním difuzního a vířivého členu a následně linearizované). Nástup výpočetní techniky znamenal možnost přidání i třetí dimenze. První práce, zabývající se třídimenzionálním

řešením

linearizovaných

potenciálních

rovnic,

byla

publikována roku 1967 [8]. Popis přístupu k výpočetnímu řešení fyzikálních vztahů a následné realizaci vlastní simulace nestlačitelných tekutin v trojrozměrném prostoru shrnul Robert Bridson ve své knize z roku 2008 [9].

12

KAPITOLA 2 – ANALÝZA Většina simulací v počítačové grafice byla navržena s upřednostňováním

vizuální stránky na úkor fyzikálního pozadí. Mezi prvními vědeckými publikacemi, které se zabývaly kvalitním třídimenzionálním grafickým výstupem, navíc podpořeným reálnými fyzikálními vztahy, byla práce Nicka Fostera a Dimitriho Metaxase z roku 1996 [10]. Autoři v ní přiblížili postup na tvorbu systému pro animaci tekutin. Obsahovala přístup k řešení NavierStokesových rovnic pomocí metody konečných diferencí a explicitních výpočtů proudového pole. Nevýhodou tohoto přístupu ale byla numerická nestabilita při nevhodné volbě časového kroku. Roku

1997

Chen,

Lobo,

Hughes

a

Moshell

[11]

simulovali

pohyb

dvoudimenzionální vodní hladiny na základě výpočtu tlaků s použitím NSR. Rovněž tato metoda ale trpěla nestabilitou. V srpnu téhož roku dvojice autorů Foster-Metaxas zveřejnila svoji další práci, ve které popisují tvorbu animací vířivých tekutin [12] (Obr.9).

Obr.9: Výsledky simulace horkého vířivého plynu (převzato z [12])

V roce 1999 představil Jos Stam [13] stabilní implicitní řešení NavierStokesových rovnic (NSR) využívající principu zpětného trasování. Toto řešení se velmi rychle ujalo a posloužilo jako vzor mnoha následujícím fluidním simulátorům. Jeho další významnou vlastností byla možnost interakce uživatele s 3D mřížkou v reálném čase. Jak ale sám autor ve článku přiznal, jeho řešení trpí rychlejším zklidněním proudového pole, než k jakému při stejných podmínkách dochází ve skutečnosti. Výsledky Stamova fluidního simulátoru jsou uvedeny na Obr.10.


13

Obr.10: Výsledky Stamova implicitního fluidního simulátoru (převzato z [13])

Společným znakem výše uvedených metod je použití Eulerovské metody řešení – tedy rozdělení prostoru do uniformní mřížky, neměnné v čase, v níž je díky snadnému přístupu do sousedních buněk umožněno jednodušeji pracovat např. s derivacemi používaných veličin. Protipólem této metody je metoda Lagrangeova, která využívá částicový systém (na modelovanou tekutinu pohlíží jako na shluk nezávislých částic pohybujících se prostorem). Liu, Wang, Gong, Huang a Peng [14] vytvořili roku 2007 simulátor písečné bouře. Jelikož písek reprezentovali neviskózním částicovým systémem, použili při výpočtech dvě proudová pole – jedno pro vzduch a druhé pro pohyb písku. Ten byl navíc doplněný o Reynoldsův tenzor pnutí, umožňující počítat vířivost za

každou

jednotlivou

částicí

unášenou

prouděním.

Vizualizace

byla

implementována pomocí rozptylu světla (scatteringu) a nechyběla ani možnost převedení výpočtů na grafickou jednotku (GPU). Simulaci prováděli nad scénou, jejíž předlohou bylo skutečné prostředí zachycené na fotografii. Na Obr.11 je výsledek simulace porovnán s reálnou scénou.

Obr.11: Porovnání simulované a reálné scény (převzato z [14])

14


Z této práce vycházeli Oreška, Stefanidis a Eršil [15] při implementaci vlastního simulátoru písečné bouře v roce 2008 (viz Obr.12). Narozdíl od své předlohy však tento simulátor místo částicového systému využíval eulerovskou mřížku.

Obr.12: Scéna simulátoru písečné bouře (převzato z [15])

Roku 2004 použili Losasso, Gibou a Fedkiw pro simulaci vody a kouře adaptivní dělení prostoru [16], což s sebou přineslo více detailů například v oblastech s víry či na rozhraní tekutiny a objektů a tím i větší realističnost simulované scény (viz. Obr.13). Informace o tlaku si udržovali uprostřed každé buňky (oktalového stromu), ostatní skalární veličiny pak v jejích rozích. Vektor rychlosti rozložili na jednotlivé komponenty a jejich hodnoty uchovávali v příslušných stěnách buňky. Takovéto rozmístění umožňovalo jednodušší interpolaci hodnot při zjemňování či zhrubnutí stromu. Pro řešení NSR použili diskretizaci, přičemž v případě, kdy ve výpočtech pracovali s buňkami různých úrovní, používali jako vyrovnávací měřítko hodnotu z hrubější buňky.

Obr.13: Výstup fluidního simulátoru s adaptivním prostorový dělením (převzato z [16])


15

Kombinaci Eulerovské a Lagrangeovy metody použili v roce 2005 Y. Zhu a R. Bridson [17] při simulaci písku jako tekutiny. Na uniformní mřížce nejprve vypočítali toky proudového pole a výsledky následně použili pro pohyb písku, reprezentovaného mračnem částic a doplněného o síly vnitřního tření. Výsledkem byl stav připomínající pohyb mokrého písku, zachycený na Obr.14.

Obr.14: Simulace tekutého písku (převzato z [17])

Hong, House a Keyser [18] představili roku 2007 částicový simulátor nestlačitelných tekutin, založený na adaptivním zjemňování částicového vzorkování. Objem každé částice v simulované nestlačitelné tekutině byl úměrný hmotě. V místech ustáleného toku docházelo ke slučování jednotlivých částic, naopak v oblastech podstatných změn toku se mohly částice dále dělit na menší a tím v tekutině vytvářet lepší detaily (tato metoda je ilustrována na Obr.15).

Obr.15: Simulace částicově reprezentované tekutiny za pomoci zjemňovaného vzorkování. Zeleně jsou obarveny částice standardní velikosti (r=1), červeně sloučené částice (r>1) a modře částice vzniklé adaptivním dělením (r<1) (převzato z [18])

16


Kupka a Ďurikovič [19] z Univerzity Komenského v Bratislavě implementovali simulátor tekutin využívající adaptivní prostorové dělení, přičemž pro určení tvaru tekutiny použili hmotu zachovávající metodu VOF (Volume Of Fluid). Tato metoda přiřazuje každé výpočetní buňce ve scéně číslo z intervalu <0;1>, přičemž 0 označuje prázdný prostor, 1 buňku vyplněnou tekutinou a mezilehlá čísla jsou použita na označení objektů (překážek). Použití oktalového stromu umožňuje poskytnout objemu tekutiny větší detaily a dodat jí tak tvar blížící se skutečnosti. Práce se ovšem nezabývala vizualizační částí - výstup simulátoru (včetně pozic modelů) byl ukládán do externích souborů čitelných pro vizualizační software PovRay. Vizualizace těchto výsledků je zachycena na Obr. 16.

Obr.16: Simulace vzduchové bubliny stoupající vodním objemem s překážkou (převzato z [19])

Z roku 2002 pochází publikace [20], která je stěžejním materiálem pro tuto diplomovou práci. Prezentuje metodu pro fyzikálně založenou simulaci kouře s využitím dynamického adaptivního dělení scény pomocí oktalového stromu. Od [16] se ale liší způsobem uchovávání veličin ve stromu. Používá totiž hybridní datovou strukturu, která kombinuje výhody oktalového stromu a uniformní mřížky. Zatímco v [16] jsou veličiny uchovány v každém uzlu stromu, tato práce používá oktalový strom jen na rozdělení scény. Používaný oktalový strom je navíc nekompletní, což znamená, že uzel nemusí mít vždy nutně přesně osm potomků. Uzel totiž dělí pouze v těch oktantech, kde jsou vyžadovány větší detaily. Informace o veličinách jsou uloženy v uniformních mřížkách, nacházejících se v listech a ve všech uzlech s méně než osmi potomky. Každá mřížka má stejný počet buněk a až samotné buňky těchto mřížek jsou pak


17

předmětem vlastních výpočtů. Pro zjemnění a zhrubnutí stromu je zde použito dělení a slučování uzlů, založené na měření maximálních gradientů v mřížkách uvnitř listů. Pokud maximální gradient v listu překročí stanovenou mez, dojde k podrozdělení tohoto uzlu a naopak – je-li maximální gradient všech potomků daného uzlu menší než druhá daná mez, tyto listy se odstraní. Pro řešení Poissonovy tlakové rovnice (viz rovnice (4) v části 2.2.1.4) používá metody vzestupné propagace rychlostí a sestupné propagace tlaků (z listů do kořene, resp. z kořene zpět do listů).

Obr.17: Porovnání výsledku simulace kouře s využitím uniformní mřížky a adaptivního dělení (převzato z [20])

2.2 Rozbor požadavků 2.2.1 Fyzikální vztahy – Navier-Stokesovy rovnice Navier-Stokesovy rovnice popisují fluidní dynamiku a používají se při výpočtech proudění v aerodynamice a hydrodynamice. Anglickým termínem fluid (česky tekutina) se označuje kontinuální, amorfní látka, jejíž molekuly se mohou volně pohybovat, a která má tendenci přebírat tvar svého úložiště; tj. kapaliny a plyny. Autory těchto rovnic jsou francouzský fyzik Claude-Louis Navier (1785-1836) a irský matematik a fyzik George Gabriel Stokes (1819-1903), kteří je odvodili nezávisle na sobě v letech 1827, resp. 1845.

18


Vycházejí z aplikace Newtonova druhého pohybového zákona na newtonské tekutiny (tj. splňující Newtonův zákon viskozity – odpor způsobený vnitřním třením v tekutině je přímo úměrný rychlosti toku). Jeden z možných zápisů NSR je uveden v rovnici (1) [21].

∂ u 1  u⋅∇  u = − ϱ ∇ p ∇ 2 u f ∂t

(1)

Není-li uvedeno jinak, budou v následujícím textu veškeré veličiny uvažovány ve trojrozměrném prostoru. u =u , v , w 

VEKTOR RYCHLOSTI

m/ s

t

ČAS

s

ϱ

HUSTOTA TEKUTINY

kg /m3

p

TLAK

hPa



KINEMATICKÁ VISKOZITA

m /s

f = f X , f Y , f Z 

EXTERNÍ SÍLY

N

2

Levá strana rovnice (1) představuje zrychlení složené z lokální a konvektivní složky v daném směru, na straně pravé pak nalezneme zrychlení způsobené gradientem tlaku (poděleného hustotou tekutiny), zrychlení potřebné k překonání třecích sil a konečně zrychlení způsobené objemovými silami. Lokální zrychlení – změna rychlosti  u s časem t v daném směru

u s prostorem v daném směru Konvektivní zrychlení – změna rychlosti  u

∂ u ∂ u ∂ u v w ∂x ∂y ∂z

Zrychlení pro překonání třecích sil – součin kinematické viskozity  (průtoku dané tekutiny za jednotku času) a druhé parciální derivace rychlosti



∂ 2 u ∂2 u ∂2 u ⋅   ∂ x ∂ y ∂z



Objemová síla – hustota síly působící v objemu tělesa

 f =lim F V 0 V

(velmi často se za tento člen dosazuje tíhové zrychlení g ).


19

Přepíšeme-li rovnici (1) tak, aby na její levé straně zůstala pouze časová změna rychlosti, dostáváme tvar

∂ u 1 2 = − u ⋅∇  u  − ϱ ∇ p ∇ u   f ∂t

(2)

Jednotlivé členy na pravé straně rovnice představují operace, které budeme nad u⋅∇   u ), projekce ( ∇ϱ p ), difuze (  ∇ 2  tekutinou provádět – advekce (   u )a přičtení vlivu externích sil ( f ). Výsledky jednoho kroku budou působit jako vstupní parametry kroku následujícího. Schéma těchto operací, tak, jak jej uvádí [13], je uvedeno na Obr.18.

Obr.18: Schéma operací

2.2.1.1 Vliv externích sil Tato operace je z celého průběhu nejsnadnější. Akceptujeme-li předpoklad, že uvažovaná externí síla se v průběhu časového (výpočetního) kroku výrazně nemění, lze vliv přičtení externích sil v bodě x a čase t vyjádřit rovnicí

u1  x ,t t  = u0  x ,t  f  x ,t ⋅ t tedy budoucí rychlost v daném místě je rovna aktuální rychlosti sečtené s externí silou působící po dobu časového kroku. 2.2.1.2 Advekce Advekcí rozumíme unášení objektů a hustot ve směru toku proudového pole. V tomto kroku ale působí také autoadvekce (samounášení), neboli propagace rychlostí a rušivých změn napříč polem. Dobrým příkladem pro

20


ilustraci advekce může být přítomnost barviva ve vodě. To se postupem času vlivem proudů působících ve vodě rozplyne v objemu. Příkladem pro autoadvekci je situace, která nastane po vhození kamene do stojaté vody - v místě dopadu se vytvoří vlnění a začne se šířit po hladině. Přítomnost advekčního členu způsobuje nelineárnost NSR. Foster a Metaxas [10] tento problém vyřešili použitím metody konečných diferencí a transport částic vyhodnocovali explicitně (Eulerovskou metodou):

r t t =r t   u t ⋅ t , kde r t označuje pozici částice v čase t . Nevýhodou tohoto přístupu je však nutnost uvažování dostatečně malých časových kroků (  t

 , kde   ∣ u∣

označuje krok prostorového dělení). Velikost časového kroku musela být volena s ohledem na velikost prostorového kroku a/nebo rychlostí – uvažovaná částice nesměla za jeden časový krok urazit vzdálenost větší než rozměr buňky (   ). V opačném případě se totiž systém stává numericky nestabilním (simulované fluidum se „rozpadne“). Stamův postup je na rychlostech a velikosti časového kroku nezávislý a systém tak zůstává stabilní po celou dobu trvání simulace. Jeho implicitní přístup spočívá v použití zpětného trasování částic (backtracking). Neurčuje tedy každé částici novou budoucí pozici (novou buňku, do které se daná částice advekcí dostane), ale pro každou buňku zjišťuje (v podstatě převrácením směru jejího aktuálního větru), čí obsah (rychlost, teplota, hustota apod.) do ní bude proudovým polem unesen. Pro zajištění větší přesnosti se po zpětném trasování provádí ještě interpolace z hodnot nejbližších sousedů (na Obr.19 vpravo vyznačeni červenou barvou).


21

Obr.19: Grafická znázornění principu zpětného advekčního trasování. Rychlost v bodě čase

x a

s je získána z rychlosti, kterou měla daná částice na své předchozí pozici (v čase − t ) (převzato z [13] a [24])

Numerická stabilita této metody vychází z faktu, že výsledné nové hodnoty nemohou být nikdy větší než hodnoty v čase

t− t . Rovnicové vyjádření

přenosu veličiny q do bodu x je následující:

q x ,t t =q x−u  x ,t ⋅ t ,t  2.2.1.3 Difuze (rozptyl) Difuzní člen popisuje chování tekutiny v závislosti na její kinematické viskozitě  a proudech v okolních buňkách. Zajišťuje rozptyl veličin (např. rychlosti, hustoty) z buněk s vyšší koncentrací dané veličiny do buněk s koncentrací nižší.

Obr.20: Difuzní schéma (převzato z [O20])

22


Tvar difuzního členu

∂ u3 = ∇ 2 u 2 ∂t

odpovídá tvaru Poissonovy rovnice.

Poissonovy rovnice jsou maticové rovnice ve tvaru A⋅x= b , kde A je matice (v tomto případě reprezentovaná Laplaceovým operátorem =∇ 2 ), x je vektor hledaných řešení a

 b představuje vektor konstant. Častým přístupem

počítačových simulací k řešení Poissonových rovnic je iterativní výpočet. V prvním kroku iterace se použije iniciální odhad řešení, který se pak v průběhu následujících iterací postupně zpřesňuje. Mezi nejjednodušší iterativní postupy patří Jacobiho iterativní metoda a Gauss-Siedelova metoda [22]. Jejich zápisy (níže z prostorových důvodů uvedené ve 2D) jsou si velice podobné, rozdíl mezi nimi spočívá ve využívání hodnot sousedů. Zatímco Jacobiho metoda v iteračním kroku n1 používá výhradně hodnoty z předchozí iterace, metoda Gauss-Siedel bere v úvahu i hodnoty sousedů vypočtené v aktuální iteraci. Z toho vyplývá, že Gauss-Siedelova metoda konverguje k řešení rychleji. Pro dosažení dobré konvergence Jacobiho metody se doporučuje 40 – 80 iteračních cyklů. Rozšíření uvedených metod do třídimenzionálního prostoru je intuitivní – veličina w (zastupující vektor rychlosti u2 ze schématu na Obr.18) získá třetí rozměr k a jmenovatelé zlomků změní svoji hodnotu na 6. Jacob:

wn1 i , j  =

1 n 1  n n 2 ⋅[ wi −1, j wn i1, j w i , j −1w i , j 1 ] ⋅qi , j  h 4 4

Gauss-Siedel:

wn1 i , j  =

1 n1 1  n1 n 2 ⋅[ wi −1, j wn i1, j w i , j −1w i , j 1 ] ⋅qi , j  h 4 4

Členem q i , j  je nahrazen diskrétní tvar divergence ∇⋅w  : −qi , j  = Člen

h2

w −w w −w ∂ wx ∂ w y  ≈ i 1, j 2i −1, j   i , j1 2i , j −1 ∂x ∂y  x   y

vznikl substitucí prostorových kroků (za předpokladu uniformní

mřížky h= x= y ). Rozšířenou praxí počítačových simulací tekutin bývá (pokud to situace umožňuje) považovat je za nestlačitelné a neviskózní. Ačkoliv ve skutečnosti


23

kapaliny i plyny svůj objem mění (jako důkaz poslouží šíření zvukových vln), bývají tyto změny velice malé a těžko proveditelné (je téměř nemožné změnit objem vody či vzduchu, a to i za použití extrémně výkonného čerpadla). Simulace stlačitelných tekutin je navíc mnohem složitější a náročnější. Pracujeme-li s tekutinami při nízkých rychlostech (cca Mach 0.3), lze v případě, že není nutné uvažovat zvukové a rázové vlny, aplikovat nestlačitelnost s dobrými výsledky i na tekutiny stlačitelné. Rozdíl mezi stlačitelnou a nestlačitelnou tekutinou je totiž v takovém okamžiku jak po vizuální, tak po kinematické stránce zanedbatelný. Nestlačitelnost (tedy zachování stejného objemu po celou dobu práce s tekutinou, kdy žádná tekutina ze scény neodchází, ani není doplňována) však ústí v nulovou divergenci rychlosti (zákon zachování hmoty):

∇ u =0

(3)

Nulová divergence ze své definice požaduje, aby všechny rychlosti (tj. v každém jednotlivém

uvažovaném

bodě)

byly

známy

již

na

začátku

výpočtů

(inicializovány v čase t=0 ). Rovněž viskozita tekutin se mění v závislosti na teplotě. Pokud ovšem není předmětem naší simulace tekutina s významnou vazkostí (jako např. med), lze i viskózní člen považovat ze simulačního hlediska za zanedbatelný. Důsledkem takového předpokladu je pak odstranění členu viskózního rozptylu z rovnice (2). Pohyb tekutiny tak bude popsán pouze svojí rychlostí a tlakem. Obě tyto veličiny se samozřejmě mění jak v čase, tak prostoru, a závisejí na počátečních a okrajových podmínkách (tj. hodnotách těchto veličin na počátku simulace, resp. na okrajích modelované scény). 2.2.1.4 Projekce Při aktualizaci rychlostí proudového pole pomocí NSR dle uvedeného řetězového schématu se nevyhneme situaci, kdy je výsledné proudové pole u ≠0 ). Předpoklad nestlačitelnosti, který jsme učinili v divergentní (tj. ∇  předchozí sekci, ale požaduje, aby výpočty v každém časovém kroku končily s nulovou divergencí pole. Naštěstí existuje matematický aparát, který nám

24


umožní situaci napravit – Helmholtz-Hodgeova dekompozice. Tento teorém říká, že každé dostatečně hladké třírozměrné vektorové pole

w 

může být

jednoznačně rozloženo na nedivergentní vektorové pole a gradient skalárového pole [25]:

w u ∇ q  =  Budeme-li za skalární pole uvažovat hodnoty tlaku v jednotlivých buňkách, můžeme nedivergentní pole rychlostí získat z divergentního snadnou úpravou – totiž odečtením tlakového gradientu

u = w  –∇p Pokud na obě strany této rovnice aplikujeme operátor divergence ∇ , dostáváme po úpravě tvar ∇⋅w u ∇ p = ∇⋅ u ∇ 2 p ,  = ∇⋅  ze kterého vlivem rovnice (3) dostáváme Poissonovu tlakovou rovnici 2

∇w  = ∇ p

(4)

Její řešení je stejné jako v případě rovnice viskózního rozptylu (část 2.2.1.3), s tím, že v prvním členu Jacobiho iterace budeme za veličinu w dosazovat tlak a člen q bude představovat divergenci rychlostí. 2.2.1.5 Počáteční a okrajové podmínky Pro úspěšnou simulaci vývoje proudového pole potřebujeme znát stav veličin na jeho počátku, tzn. nadefinovat počáteční podmínky. Rozšířenou praxí je začínat s nulovými počátečními podmínkami, ovšem v obecném hledisku je jejich nastavení záležitostí potřeb každého jednotlivého simulátoru. Při řešení Navier-Stokesových rovnic na diskrétní mřížce vektorů jsou ve výpočtech užívány také hodnoty sousedů aktuální buňky. Dokud se výpočty týkají vnitřních buněk, jež mají ve všech uvažovaných směrech sousedy, je vše v


25

pořádku. Situace se ale změní, když se dostaneme na okraj simulační mřížky. V tom okamžiku je potřebovat zvážit, jakým způsobem se vypořádáme s okolím – tedy stanovit okrajové podmínky. Dva nejčastěji používané druhy okrajových podmínek, určujících chování simulace na hranicích scény, jsou Dirichletova a Neumannova. Dirichletova okrajová podmínka říká, že vektor proudění podél normály stěny, sdílené scénou a okolním prostředím, je nulový (  u⋅n=0 ). To znamená, že do buňky, která není součástí řešení (simulace), nesmí procházet žádné proudění. Scéna se tak stane „uzavřenou“. Oproti tomu Neumannova podmínka požaduje, aby byla nulová změna toku mezi okrajovými a okolními buňkami (tedy derivace podle normály,

∂ u =0 ). ∂ n

Okrajové podmínky se však v tomto případě nebudou aplikovat pouze na okraje scény, ale také na rozhraní tekutina-objekt, čímž bude zajištěno obtékání objektů. Všechny objekty ve scéně budeme považovat za uzavřené a tudíž ve všech částech prostorového dělení, zasahujících do modelů, budeme uvažovat nulovou rychlost proudění a nulovou hustotu písku. Tím pádem backtracking vycházející z těchto buněk nebude přebírat odnikud žádné hodnoty a naopak, budou-li tyto buňky cílem zpětného trasování, žádná veličina se z nich do výchozí buňky nepřenese.

2.2.2 Adaptivní dělení Adaptivní dělení přináší pohled na zacházení se scénou odlišný od uniformní mřížky. Zatímco uniformní dělení provede prosté rozčlenění scény na požadovaný počet stejně velkých krychlí, jejichž množství ani rozměry se během simulace nemění, adaptivní dělení (reprezentované stromem v podstatě libovolné arity n) lépe využívá prostor a přizpůsobuje se obsahu scény (Obr.21 a 22). To znamená, že v okolí modelů a jiných objektů dojde k hustějšímu rozdělení, zatímco volné oblasti mohou být obsazeny třeba i pouhou jednou velkou výpočetní buňkou.

26


Obr.21: Porovnání uniformní mřížky a adaptivního dělení ve 2D (quadtree) (převzato z [O21])

Obr.22: Porovnání uniformního a adaptivního dělení v prostoru (octree) (převzato z [O22])

Kromě již zmíněného hospodaření s prostorem scény (a potažmo pamětí počítače) je rozdíl mezi těmito metodami například v postupu při hledání sousedů. V uniformní mřížce jsme díky indexaci buněk vždy schopni rychle a efektivně najít v požadovaném směru souseda (pokud existuje). V případě adaptivního dělení ovšem musíme použít vhodnou traverzaci (průchod) stromu a vypořádat se s různými úrovněmi hloubek sousedních uzlů. V adaptivním dělení budu vycházet z publikace [20], která se zabývá problematikou fluidní simulace za použití oktalového stromu. Stejně jako zde použiji hybridní datovou strukturou - oktalový strom a uniformní mřížku. Zatímco oktalový strom zajistí adaptivní dělení scény, uniformní mřížka mi umožní diskretizovat výpočty NSR do Eulerovského tvaru. V implementaci však


27

provedu několik drobných změn. Namísto neúplného oktalového stromu použiji úplný - zvýší se sice paměťová náročnost, ale předpokládám, že tak dosáhnu lepších detailů. Veličiny (rychlost a směr proudění, hustota (obsazenost buňky) a tlak) budou uloženy výhradně v listech, avšak kvůli potřebě vzestupných a sestupných propagací veličin zavedu prázdné uniformní mřížky pro mezivýpočty ve všech uzlech stromu. Mřížky budou mít vždy stejný počet voxelů, nicméně jejich rozměry budou určeny umístěním mřížky v různých úrovních stromu. Právě tento mechanismus následně vytvoří ve zobrazovanému mračnu detaily.

Obr.23: Uzel oktalového stromu a jeho vnitřní struktura

Všechny uvažované veličiny budu vztahovat vždy ke středu dané výpočetní buňky (voxelu). Množství písku uvnitř výpočetní buňky bude omezena na interval <0; 63>, všechny vyšší hodnoty budou oříznuty na horní mez tohoto rozsahu. Uvedené omezení je dáno vykreslovacím modulem, který bude při zobrazování bouře používat sadu 64 textur s různým stupněm průhlednosti. U rychlosti proudění je nutné rozmyslet, jestli budu uvažovat pouze rychlost samotného vzduchu, nebo zda zavedu další doplňkovou veličinu (rychlost pohybu pískových zrn), případně použiji dvě proudová pole, stejně jako v [14]. Je zřejmé, že rychlost písku (stejně tak jako směr jeho pohybu) bude velmi podobná rychlosti proudění vzduchu, který jej unáší. Do celé situace je ale nutné zanést také externí síly. V tomto případě je jedinou externí silou gravitační zrychlení, a tak, ačkoliv to není fyzikálně korektní (tíhové zrychlení nemůže ovlivnit nehmotné vzdušné proudění), jelikož by byla hmotnost jednotlivého

28


zrnka natolik malá, že by se gravitace prakticky téměř neprojevila, budu uvnitř buněk uvažovat pro vítr i písečná zrna pouze jednu společnou rychlost. A díky tomu, že nepoužiji částicovou reprezentaci, mi bude stačit pouze jedno proudové pole, kterým bude písečná hustota unášena. Důležitou částí implementace bude hledání sousedů, a to jak pro uzly stromu, tak i pro samotné voxely. Jelikož tuto operaci bude nutné při výpočtech provádět velmi často, je potřeba pro ni vytvořit efektivní algoritmus.

2.2.3 Simulace pískového podloží Efekt písečné bouře by nebyl úplný, kdyby nedocházelo k ovlivňování okolí, obzvláště pískového podloží. Pro zlepšení realističnosti tak součástí simulace budou i procesy, uvedené a popsané v kapitole 1.3.1 – vzedmutí, unášení a splaz. Samotné pískové podloží bude vytvořeno pomocí výškového pole. Jedná se implementačně jednoduchou metodu s dobrými vizuálními výsledky. Výškové pole je reprezentováno trojrozměrným grafem - maticí bodů v rovině XZ, přičemž Y-ová souřadnice každého bodu určuje jeho vertikální pozici (výšku). Při vhodné volbě rozlišení (tj. počtu bodů v každé ose) a vykreslovacího mechanismu vytváří výškové pole velmi realistický efekt členité krajiny (viz Obr. 24). Pro vykreslení lze použít např. síť trojúhelníků, resp. trojúhelníkový pás (triangle strip), případně rozšíření OpenGL nazvané vertex buffer object. Vzhled terénu se může dále vylepšit patřičným obarvením či nanesením vhodné textury.

Obr.24: Příklad terénu vymodelovaného výškovým polem (převzato z [O24])


29

2.2.3.1 Vzedmutí Budeme-li znát pozici konkrétního bodu výškového pole, můžeme provést operaci vzedmutí. Nejprve vyhledáme voxel, do kterého daný bod spadá a přečteme rychlost a směr jeho větru. Pokud splňuje podmínky pro vzedmutí písku (směr a dostatečná síla), odebereme z daného bodu výškového pole určité množství písku a v příslušném voxelu inkrementujeme hustotu. 2.2.3.2 Unášení Unášení písku (saltation) bude, podobně jako advekce, využívat principu zpětného trasování, při kterém se pracuje se dvěma body, výchozím a tzv. backtrackem (bodem nalezeným pomocí zpětného trasování), a dvěma konstantami ( L0 pro vzdálenost a q 0 pro množství přenášeného materiálu). Pro každý bod výškového pole se nejprve zjistí, zda neleží pod objektem a může tak do něho přibýt písek z okolí. Následně se zjistí voxel, v jehož objemu daný bod výškového pole leží, přečte se vítr a ten se, vynásobený konstantou L0 , odečte od pozice aktuálního vrcholu, čímž dostaneme backtrack. Pokud se i písek z tohoto bodu může účastnit unášení, je aktuální hodnota snížena o určité množství písku ( q 0 ), které je přičteno k výšce výchozího bodu výškového pole.

Obr.25: Schematické znázornění principu unášení

2.2.3.3 Splaz Poslední operací pro výškové pole je splaz (creep), který se stará o vyhlazování terénu. Jak bylo uvedeno výše, podloží scény bude implementováno

30


jako výškové pole, tedy prostou číselnou informací o výšce na dané pozici. Jemnost takového pole se dá částečně ovlivňovat jeho rozlišením (počtem buněk v jednotlivých směrech), nicméně ani tak ještě nemusí výsledek vypadat tak, jak bychom si přáli. Ani na velmi jemně naděleném poli totiž nezabráníme vytváření špiček v místech, kde se vůči okolí nahromadí více písku. Vyhlazovací operace bude při své činnosti iterativně procházet celé výškové pole a v každém bodě, který se nenachází pod objektem, zjistí gradient výšky. Pokud vypočítaná hodnota překročí určitou mez strmosti, část materiálu (úměrná gradientu) se z tohoto bodu rozdělí sousedním buňkám. Výsledkem pak bude opticky hladký povrch.

2.2.4 Statická scéna Objekty uvnitř simulovaného prostředí mohou být buď pohyblivé (dynamické) či nepohyblivé (statické). Při použití adaptivního členění by změny pozice objektů (nebo jejich částí) znamenaly zvýšenou zátěž výpočetního aparátu, neboť při každém takovém prostorově významném pohybu by muselo dojít k přestavbě určité oblasti (v mimořádných případech až kompletně celého) dělení. Přípustnost statické scény tak umožňuje výrazné snížení režie na manipulaci.

KAPITOLA 3 – NÁVRH ŘEŠENÍ

31

Kapitola 3 – Návrh řešení Modul pro fyzikální výpočty, který je předmětem této práce, by měl být zcela nezávislý na modulu vykreslovacím (viz [1]), ačkoliv oba pracují nad stejnými daty podle následujícího schématu.

Obr.26: Schéma spolupráce mezi jednotlivými částmi celého výsledného simulátoru

Simulace bude probíhat iterativně, tzn. jednotlivé kroky se budou opakovat ve smyčce po celou dobu simulace. Jejich schéma je znázorněno na Obr.27.

Obr.27: Schéma jednotlivých kroků simulátoru (včetně grafické části, vyznačené šedě)

3.1 Volba programovacího jazyka Jedním z požadavků na tuto diplomovou práci byla její implementace v jazyce C nebo C++. Během téměř celého studia oboru Výpočetní technika na ČVUT FEL mě doprovázel programovací jazyk C++, se kterým jsem se velmi dobře sžil. Rozhodnutí mezi danými možnostmi tak mělo jednoznačného vítěze – jazyk C++.

32


3.2 Návrh řešení pro hledání sousedů V této části budou představeny postupy hledání sousedů pro uzly oktalového stromu a pro výpočetní buňky (voxely) v uniformní mřížce.

3.2.1 Sousedé v oktalovém stromu Při hledání sousedů v octree se budu snažit využít množinu pravidel, vyplývajících z prostorové sousednosti, tak, aby se zrychlil přístup k sousedním uzlům a nebylo nutné traverzovat celý strom vždy znovu od kořene. Nejprve je vhodné uvést, jakým způsobem budou potomci v oktalovém stromu indexováni.

Obr.28: Indexace potomků v COctree (vlevo graf, vpravo prostorové uspořádání)

Pro určité dvojice (index, směr hledání) bude získání souseda triviální: (lichý, -X) (sudý, +X) (<4, -Y) (>3, +Y) ({2,3,6,7}, -Z) ({0,1,4,5}, +Z) V takové situaci je soused sourozencem daného uzlu a stačí tak pouze příslušně upravit index - v ose X zvýšit (+X), resp. snížit (-X) o 1, v ose Y snížit (+Y), případně zvýšit (-Y) o 4 a v ose Z zvýšit (+Z) či snížit (-Z) o 2. Například sousedem pro uzel B1 ve směru +Z je uzel B3. Složitější situace nastane, pokud požadovaná kombinace (index, směr) nebude patřit do výše uvedeného výčtu. Může jít o jakoukoliv jinou dvojici, např. (H2,


33

+Y). Z Obr.28 je patrné, že v takovém případě hledaný soused zcela jistě není sourozencem aktuálního uzlu a musí tak dojít k nalezení jiného rodiče (rodiče jsou na Obr.28 označeni písmeny), jakéhosi společného prapředka. V mezním případě je tímto prapředkem až vlastní kořen, avšak ne vždy k takové situaci dojde. Cestou vzhůru si budeme vytvářet seznam navštívených uzlů. V okamžiku, kdy dorazíme do prapředka jsme vlastně přestoupili z množiny „složitých“ dvojic do množiny dvojic „snadných“, uvedených výše, a proto, neníli prapředkem kořen, nastavíme si pomocný příznak flip (otočení, změna množiny dvojic), který nám následně pomůže při zpětném sestupu. Při něm zpracováváme seznam navštívených uzlů od konce a na jejich indexy aplikujeme pravidla pro „snadné“ dvojice. Pokud jsme ovšem vystoupili až do kořene, hledaný soused ve scéně neexistuje.

Obr.29: Binární strom pro ilustraci prapředků

Pro uzly H a J na Obr.29 je nejbližším společným prapředkem uzel A, zatímco např. pro uzly G a M je nejbližším společným prapředkem kořen celého stromu.

3.2.2 Sousedé v uniformní mřížce Známe-li směr hledání a pozici středu výchozího voxelu, můžeme pomocí kroku dělení mřížky vypočítat pozici středu voxelu cílového. Pokud cílový voxel spadá do stejné mřížky jako výchozí, nahlásíme jej rovnou jako hledaného souseda. Pokud ovšem patří do jiné mřížky, musíme nejprve na základě vypočítaného nového středu nalézt odpovídající uzel oktalového stromu, do kterého voxel patří, a následně tyto souřadnice přepočítat na osové indexy. Z osových indexů budeme poté ještě muset zjistit absolutní index v mřížce, jelikož voxely jsou v mřížkách uloženy v jednorozměrném poli. Tím pádem bude nutné použít mapovací funkci ze 3D na 1D.

34


Numerický index (id) voxelů se bude v každé mřížce pohybovat v rozmezí od 0 do počtu buněk-1. Na každou osu připadne počet buněk řízený konstantou (GRID_SIZE), číslování buněk poroste nejprve ve směru osy X, následně v ose Y a až poté v ose Z, avšak v záporném směru (viz Obr.30).

Obr.30: Indexace voxelů uvnitř mřížky pro 8 buněk v každé ose

Podoba mapovací funkce přepočítávající osové indexy na index absolutní (  x , y , z  id ) bude následující: id = x + y*GRID_SIZE + z*GRID_SIZE*GRID_SIZE;

Například pro hodnoty (x,y,z) = (6,1,2) a GRID_SIZE=8 získáme index id = 6+ 1*8 + 2*64 = 142. Máme-li tento přepočet, bude vhodné najít i jeho protějšek, tedy mapovací funkci z absolutního indexu na osové indexy ( id  x , y , z  ). Její zápis je o něco delší, avšak postup je analogický předchozímu. z y x

= id / (GRID_SIZE*GRID_SIZE); rem = id % (GRID_SIZE*GRID_SIZE); = rem / GRID_SIZE; rem = id % GRID_SIZE; = rem;

Například pro buňku s id = 507 pak dostaneme výsledné osové indexy z=507/64=7, y=59/8=7, x=3, tedy (x,y,z) = (3,7,7).


35

3.2 Přístup k řešení Navier-Stokesových rovnic V této chvíli je potřeba stanovit vztah mezi proudovým polem a unášeným obsahem (pískem). Z popisu písečné bouře v úvodní kapitole vyplývá, že nejtěžší zrna písku jsou větrnými proudy zdvihána do výšky maximálně několika desítek centimetrů a unášena na velice krátkou vzdálenost. Tento proces bude ovládán simulací písečného podloží. Drtivá většina bouřkového mračna

je

tvořena

obrovským

množstvím

velice

drobných,

téměř

mikroskopických prachových a písečných částic. To umožňuje na mrak pohlížet stejně jako na plyn (resp. kouř) a uvažovat o něm tudíž jako o neviskózní a nestlačitelné tekutině. Při implementaci tak bude stačit počítat jedno proudové pole pro vítr a jím nechat drobné částečky písku unášet. Při implementaci NSR budu za jednu z externích sil považovat gravitační

g ). Jako druhá externí síla bude muset vystupovat zdroj větrného zrychlení (  proudění ( us ), který bude do scény v určitém bodě dodávat nový vítr, a bez jehož přítomnosti by se mračno téměř nepohybovalo a brzy by došlo k jeho rozplynutí. Poté, co jsme v kapitole Analýza provedli předpoklad neviskozity, dojde k vypuštění difuzního členu z rovnice (2). Pokud dodržím požadavek počátečních a okrajových podmínek (počáteční bezvětří a použití Dirichletovy okrajové podmínky), rovnice (2) pak bude mít při implementaci následující tvar:

∂ u 1 =  g us − u⋅∇  u− ϱ ∇ p ∂t

(5)

3.2.1 Návrh řešení advekčního kroku Při výpočtech unášení je nutné překonat několik překážek. První z nich je nalezení správné mřížky, ve které je obsažen voxel, z něhož budeme přebírat hodnoty (tzv. backtrack). Pokud je vítr v aktuální buňce (příp. časový krok) dostatečně malý, může se stát, že se backtrack bude nacházet ve stejné mřížce, jako aktuální voxel. Ve většině případů však dojde k situaci, kdy bude nutné nalézt zcela jinou mřížku pomocí traverzace stromu. Poté, co mřížku nalezneme,

36


musíme ze světových souřadnic pro backtrack zjistit, o kterou konkrétní buňku v mřížce se jedná, nalézt sedm jejich nejbližších sousedů (viz Obr.19) a nakonec převzít interpolaci jejich hodnot (v případě, že se některý ze sousedů nachází uvnitř zneplatněného uzlu, přebírají se nulové veličiny). Nové hodnoty musíme navíc následně ukládat do duplikovaného seznamu voxelů, aby nedocházelo ke kolizním přístupům čtení starých a zápisu nových hodnot ve stejné množině voxelů.

Obr.31: Znázornění interpolačních koeficientů (koeficient r2 je pro osu Y)

Interpolace přebírané hodnoty probíhá v závislosti na vzdálenosti bodu dopadu backtracku od okolních buněk (Obr.31). Její rozpis je rozsáhlejší a kvůli přehlednosti je tudíž uveden v příloze A.

3.2.2 Návrh řešení projekčního kroku Při implementaci výpočtu Poissonovy tlakové rovnice (4) se využije vzestupné propagace rychlostí a sestupného dopočítávání tlaků. První krok zajistí přenos rychlostí z listů až do kořene, odkud se pomocí druhého zpracovávají (Jacobiho iterační metodou, viz. 2.2.1.3) na cestě zpět do listů příslušné tlaky tlaky. Při cestě vzhůru se využívá průměrování – zatímco na hladině n je určitý objem reprezentován jedním voxelem, na hladině n+1 jej určuje voxelů osm. Při propagaci rychlostí směrem vzhůru se tedy do voxelu na hladině n uloží aritmetický průměr z rychlostí osmi voxelů hladiny n+1. Pro vzestupnou propagaci je nutné najít tzv. „předposlední“ uzly – tedy uzly, jejichž všichni potomci jsou listy. Na Obr.32 jsou předposledními uzly č. 3 a 4.


37

Obr.32: Ilustrace předposledních uzlů binárního stromu (uzly 3 a 4)

Při jejich hledání začneme od kořene stromu a rekurzivně zpracujeme potomky aktuálního uzlu. Jakmile nalezneme nějaký předposlední uzel, provedeme převzetí a zprůměrování rychlostí z listů a označíme daný uzel za zpracovaný. V okamžiku, kdy je všech osm potomků jednoho rodiče zpracováno, dopočítá se divergence proudového pole, rekurze v daném místě skončí a výpočet se přenese do vyššího patra (což mimochodem vysvětluje přítomnost uniformních mřížek pro mezivýpočty ve všech uzlech). Toto se provádí tak dlouho, dokud nejsou průměrné rychlosti z celé scény doručeny do kořene. Na Obr.33 je tato operace znázorněna schematicky ve čtyřech krocích (předpokládejme, že rekurzivní volání metody na všechny uzly proběhlo v nultém kroku a již jsme nalezli všechny listy). Z prostorových důvodů byl použit pouze binární strom. Zeleně označené jsou zpracovávané uzly – tedy ty, z nichž jsou v daném kroku vybírány a průměrovány rychlosti proudového pole. Šedou barvou jsou zvýrazněny již hotové uzly.

Obr.33: Ilustrace principu přebírání rychlostí z nižších úrovní

38


Po proběhnutí čtvrtého kroku se v kořenu nacházejí zprůměrované rychlosti celé scény a může přijít na řadu iterativní řešení Poissonovy rovnice. Tento výpočet, probíhající rekurzivně shora dolů, vypočte na základě divergence rychlostí tlaky v příslušných buňkách. Gradienty těchto tlaků se odečtou od rychlostí proudového pole (Helmholtz-Hodgeova dekompozice) a dostaneme v každé buňce nové, nedivergující pole rychlostí.

O přenos veličin v obou směrech se stará rodič – tedy buňka na vyšší úrovni. Využívá přitom indexace a prostorové koherence. Existuje-li mřížka na nižší úrovni, pak objem reprezentovaný jedním voxelem na úrovni n je reprezentován osmi voxely na úrovni n+1 (dochází tedy k imaginárnímu podrozdělení voxelu (dělení však ve skutečnosti probíhá pouze nad uzly oktalového stromu)). Poté, co dojde k rozdělení stromového uzlu, bude příslušná oblast, pokrytá aktuálně jednou mřížkou, nově reprezentována osmi novými (jemnějšími) mřížkami. Situace je ilustrována na Obr.34. Stejně obarvené nové voxely reprezentují dohromady totožný objem jako jejich „rodič“.

Obr.34: Ilustrace „podrozdělení“ voxelů doplněná o indexy buněk vyšší úrovně („rodičů“)

Obr.35: Rozdílná indexace mřížek (vlevo) a uzlů oktalového stromu (vpravo)


39

Pro znázornění metody hledání nižších voxelů uvedu příklad, pracující se zelenou buňkou (id=60). Při podrozdělení uzlu vznikne v objemu osm nových oktantů (Obr.34 vlevo), z nichž každý je reprezentován samostatnou mřížkou. To znamená, že na místě původních čtyř buněk (zelená, žlutá, červená a modrá) vznikne v daném směru buněk osm. Zjistíme prostorové indexy původní buňky 60 (x,y,z) = (4,7,0). Z těchto jednotlivých indexů vypočítáme zbytky po dělení čtyřmi (polovičním rozměrem). Dostáváme hodnoty (vx,vy,vz) = (0,3,0), což jsou vlastně prostorové indexy buňky 60 v mřížce 4x4x4. Tyto prostorové indexy použijeme na hledání příslušných voxelů v nové mřížce (Obr.34 vpravo) pomocí mapovací funkce 3D  1D. Jelikož je ale v nové mřížce původní objem reprezentován dvěma buňkami v každé ose, budeme při hledání nových voxelů přeskakovat po dvou. Voxel1 = 2*vx + 2*vy*8 + 2*vz*8*8 = 2*0 + 2*3*8 + 2*0*64 = 48 Voxel2 = 2*vx+1 + 2*vy*8 + 2*vz*8*8 = 1+48+0 = 49 Voxel3 = 2*vx + (2*vy+1)*8 + 2*vz*8*8 = 0+56+0 = 56 Voxel4 = 2*vx+1 + (2*vy+1)*8 + 2*vz*8*8 = 1+56+0 = 57 Voxel5 = 2*vx + 2*vy*8 + (2*vz+1)*8*8 = 0+48+64 = 112 Voxel6 = 2*vx+1 + 2*vy*8 + (2*vz+1)*8*8 = 1+48+64 = 113 Voxel7 = 2*vx + (2*vy+1)*8 + (2*vz+1)*8*8 = 0+56+64 = 120 Voxel8 = 2*vx+1 + (2*vy+1)*8 + (2*vz+1)*8*8 = 1+56+64 = 121 Když známe indexy nových voxelů, musíme ještě znát, ve kterém oktantu (neboli ve které z osmi nových mřížek) je máme hledat. To zjistíme přepočtem osových indexů (x,y,z) na absolutní index v mřížce s rozměrem 2x2x2 (neboť máme 8 oktantů, tedy dva v každé ose). Ze vzorce x/4 + (y/4)*2 + (z/4)*2*2 dostáváme výsledek 4/4 + (7/4)*2 + (0/4)*4 1+2+0 = 3. Voxely se tedy nacházejí v mřížce tři. A ta náleží třetímu potomku daného uzlu stromu (viz Obr.35).

40


3.3 Návrh řešení interakce se statickými objekty ve scéně Tento bod požaduje, aby se statické objekty a písečná bouře vzájemně ovlivňovaly – tedy objekty budou působit jako překážky, kterým se má bouře vyhnout a naopak písek unášený proudovým polem se bude na těchto objektech usazovat. Výsledná scéna tak bude působit mnohem reálněji.

3.3.1 Interakce objektů a proudového pole Statické objekty ve scéně musejí působit jako překážky, kolem kterých bude proudění obtékat. Musí být tedy vytvořen aparát, který toto bude schopen zajistit. K tomu bude potřeba znát pozici objektů ve scéně a ve vztahu k prostorovému rozdělení a oktalový strom v příslušných oblastech co nejvíce podrozdělit. V uzlech stromu (či ještě lépe v samotných voxelech), které budou zasahovat do objektů, nebude možné uvažovat dané veličiny a budou muset být tím pádem také vyřazeny ze všech kroků simulace (především z advekce). Když se bouřkový mrak dostane do malé vzdálenosti od objektů, je potřeba určit, jak se má zachovat. Je jasné, že by neměl procházet skrz objekty. Opticky by se sice jevilo, že mračno objekt těsně obklopuje a pohlcuje, jako ve skutečnosti, ale na druhou stranu při detailním pohledu by výsledky mohly být neuspokojivé (např. kus mraku vycházející z objektu). Mračno by tedy mělo objekt obcházet, a to pokud možno co nejtěsněji. Ovšem „co nejtěsněji“ v důsledku znamená mít v daném místě velmi jemné dělení (tj. velkou hloubku oktalového stromu), což se nevyhnutelně projeví na výpočetní náročnosti. Na Obr.36 je znázorněn princip zjemňování oktalového stromu. Máme-li pouze počáteční uzel stromu (kořen), objekt do něj zasahuje celým svým objemem. Stejná situace je i v další úrovni, do všech listů ještě stále zasahuje nějaká část objektu. Při hloubce 2 se již začíná rýsovat (zatím velmi hrubý) obrys objektu. Ten se ale s rostoucí hloubkou stromu postupně zjemňuje a objekt je tak obklopován se stále větší přesností. S každou kompletní úrovní, o kterou se zvětší hloubka stromu ale vzroste počet listů (a spolu s nimi také výpočetních buněk) 8x.


41

Obr.36: Ukázka přizpůsobení adaptivního dělení objektům ve scéně v souvislosti s obtékáním. Šrafování představuje zneplatněné uzly (s využitím [O36])

Jiný druh obtékání by v navrhované implementaci nebyl dost dobře proveditelný. Například proudové obtékání, tak, jak jej známe z reálného světa (Obr.37 vlevo) naráží na problém správného určení natočení proudu tak, aby sklouzl po povrchu objektu. Tento problém je navíc ještě umocněn používáním zpětného trasování v advekčním kroku (Obr.37 vpravo) – nové hodnoty (vyznačeny tmavě modrou barvou) se vypočítávají ze starších (světle modré), pohybem zpět v čase, bez jakékoliv znalosti prostoru před sebou a jeho geometrie. Mračno se tak v blízkosti zneplatněných buněk pozastaví a pokračuje pouze ve volném okolí modelu. I tak lze ale poměrně dobře simulovat větrný stín za objektem, neboť poté, co se mrak u objektů rozdělí, nějakou chvíli trvá, než se opět spojí dohromady.

Obr.37: Vlevo: Reálné obtékání objektu proudovým polem (převzato z [O37]) Vpravo: Situace při backtrackingu

42


3.3.2 Interakce objektů a unášeného písku Další variantou interakce je kolize objektů s písečnými zrny unášenými bouří. Tato interakce se nejlépe vyhodnocuje při použití částicového systému. Implementace této práce ale počítá namísto částic pouze s hodnotami hustot (obsazení) dané buňky, proto bude potřeba tento údaj dočasně převést na částicovou reprezentaci. Na náhodných pozicích v objemu každé buňky se vygeneruje určitý počet zrn v závislosti na její obsazenosti (hustotě), každé zrno se vyšle ve směru větru příslušné buňky a následně se za použití algoritmu pro test průniku trojúhelník-paprsek vyhodnotí jednotlivé kolize. Pro tuto operaci bude pro snížení zátěže stačit vytvořit v každém iteračním kroku vždy jen jednu instanci zrna a její parametry měnit podle podmínek v dané výpočetní buňce. Po vyhodnocení všech zrn bude tato instance opět smazána, což bude znamenat nižší paměťové nároky, než v případě, kdy by bylo částicově reprezentováno celé mračno (kvůli požadovanému efektu hustého mraku by částic muselo být velmi velké množství). Testování tak velkého počtu vygenerovaných zrn by ovšem znamenalo vysokou výpočetní zátěž a proto bude nezbytné zavést další opatření, která tyto kolizní testy zrychlí. Prvním z nich bude omezení počtu generovaných částic v každé buňce. Druhým mechanismem je obklopení modelů hraničními obálkami (BVH) typu AABB (axis-aligned bounding-box), dělených pomocí kd-stromu. V první fázi tedy provedeme test paprsek-AABB. Pokud částice neprotne kořen kdstromu, víme, že model minula. Když ale ke kolizi dojde, sestoupíme ve stromu do nižší úrovně a provádíme další test, tentokrát vůči oběma polovinám původního objemu, z nichž jednu částice určitě protne. Tento postup opakujeme tak dlouho, dokud po příslušné stromové větvi nedojdeme až do listu. Až teprve v něm je uložena informace o prostorově příslušných trojúhelnících objektu a můžeme přistoupit ke druhé fázi – testu kolize paprsek-trojúhelník. Zavedením BVH se asymptotická časová složitost testování průniku sníží z O(n) na O(log2 n),

kde

n

značí

celkový

počet

objektových

geometrických

primitiv

(trojúhelníků). Vlastní mechanismus testování kolizí paprsek-AABB a paprsektrojúhelník je blíže popsán v práci „Vizualizace písečných bouří“.


43

Rovněž písek přesouvaný po povrchu podloží musí brát v potaz výskyt objektů a zastavovat se o ně. K tomu bude opět nutné znát pozici těchto objektů, tentokrát ve vztahu k podloží a příslušné oblasti podloží označit jako obsazené. V takových oblastech pak nebude docházet k žádným operacím s podložím.

3.4 Aktualizace oktalového stromu Po každém simulačním kroku je nutné aktualizovat prostorové dělení, tzn. spočítat ve všech mřížkách gradienty hustot a tam, kde je to nutné, zajistit podrozdělení a tím přidat detaily, případně naopak již dále nepotřebné uzly opět sloučit a tím snížit aktuální nároky na paměť.

44


KAPITOLA 4 – IMPLEMENTACE

45

Kapitola 4 – Implementace Pro psaní aplikací v jazyce C++ je pravděpodobně nejvhodnějším vývojovým prostředím Microsoft Visual Studio. Při implementaci zde popisovaného simulátoru byla použita verze 2008 Professional Edition. Pozn.: Pro zjednodušení a sjednocení zápisu bude v ukázkách kódů v následujícím textu u všech prováděných přístupů (i v případě ukazatelů) používána tečková notace.

4.1 Implementace fyzikálního modulu V kapitole 2.2.1 jsme si popsali Navier-Stokesovy rovnice a princip jejich řešení. V této části si ukážeme, jakým způsobem bylo toto řešení naimplementováno.

4.1.1 Struktura modulu – třída Solver Třída fyzikálního modulu Solver se skládá z jediné proměnné, avšak o to více obsahuje metod. Jejich výčet společně s následným popisem je uveden níže. class Solver { private: //---ukazatel na kořen oktalového stromu, nad nímž //budou prováděny výpočty COctree *octree; /* hlavní metody pro řešení NSR */ //---přidání externích sil void addForces(); //---advekční krok void advect(); //---gravitace void gravity(); //---projekční krok void project(); /* submetody pro řešení NSR */ //---předání rychlostí z nižší úrovně do vyšší void propagateVelocities(COctree *node); //---výpočet divergence rychlosti ve voxelech void computeDivergence(SGrid *grid); //---výpočet tlaku void computePressure(COctree *node);

46

KAPITOLA 4 – IMPLEMENTACE //---přičtení nových rychlostí void projectNewVelocities(SVoxel *voxel); //---odstranění mezivýpočtů void clearGrid(COctree *node); /* duplikace výsledků pro vykreslovací modul */ //---duplikace voxelů void duplicateVoxels(); /* zpracování písku */ //---splaz void creep(int iterations); //---unášení písku void saltation(); //---vzedmutí písku void suspension(SVoxel *voxel, int x, int z); //---výpočet kolizí písku s objekty void sandShooting(int i); public: //---konstruktor Solver(COctree *tree); //---destruktor ~Solver();

};

//---hlavní metoda modulu void solve();

V konstruktoru se provádí základní inicializace modulu: Solver(tree) { this.octree  tree; initialize random seed; set random wind in voxels close to storm origin; ensure the wind is not set in voxels with vertices; }

Veškerá činnost fyzikálního modulu je volána v metodě solve(), která se ve smyčce opakovaně spouští po celou dobu běhu simulace a volá jednotlivé výpočetní metody v následujícím pořadí:

}

solve() { addForces(); advect(); suspension(); creep(); project(); duplicateVoxels();


47

Podrobnější popis jednotlivých hlavních metod řešení NSR je vhodné uvést ve stejném pořadí, v jakém jsou tyto metody volány při simulaci. Všechny uvedené submetody pro řešení NSR jsou pak součástí projekčního kroku a budou tedy rozebrány až na samém konci této části. Prvním krokem každé simulační iterace je vložení externích sil. Sem patří jednak vliv zdroje proudění, který do scény na stanoveném místě po určitý počet kroků simulace (řízený konstantou SOURCE_LIFE) přidává hustotu mračna a čerstvý vítr a dále také gravitační zrychlení, působící na unášený písek. Jelikož je však písečné mračno tvořeno velmi malými a lehkými částicemi písku, vliv gravitace by se prakticky neprojevil. Proto je použit dočasný přepočet 1 bodu hustoty na určité množství částic (konstanta PARTICLE_PACK = 1000). Hodnotu výsledné působící síla tak spočítáme pomocí vzorce

Fg = ϱvoxel⋅PARTICLEPACK⋅m⋅g Poslední neznámou zůstává váha částice. [23] uvádí rozmezí 10 μg až 0,5 mg. Pro naše účely tedy použijeme údaj z levé hranice intervalu, který ještě následně doplníme o náhodně vygenerovanou hodnotu z intervalu reálných čísel

〈−1 , 1 〉 μg. V závislosti na hustotě obsažené v daném voxelu se tak výsledná gravitační síla může pohybovat v rozmezí 0 až cca 7 Newtonů. Implementaci metody addForces si nyní stručně znázorníme. addForces() { sourceLife--; for each voxel in scene { //---vynechání zakázaných voxelů if (node containing voxel has vertices) skip; //---přidej hustotu a rychlost if (voxel is close to storm origin) { //rychlost větru přidáváme vždy add wind velocity to voxel; //hustotu pouze po dobu činnosti zdroje if (sourceLife > 0) add density to voxel; }

48

KAPITOLA 4 – IMPLEMENTACE //---přidej gravitaci weight  ((100+rand(-10,10)) * 1e-7); gravity  voxel.density*weight*PARTICLE_PACK*GRAV; voxel.addVelocity(0, -gravity, 0); }

}

Dalším krokem po přidání externích sil je vyřešení advekce, neboli unášení tekutiny proudovým polem, v metodě advect. Její stabilitu zajišťuje použití Stamova implicitního přístupu. advect() { newList  allocate new list of voxels; for each voxel in scene { clone  corresponding voxel from newList; //---zjistíme rychlost v daném voxelu velo  voxel.velocity; //---odečteme vítr od pozice středu voxelu _x  voxel.center.x – velo.x * DT; _y  voxel.center.y – velo.y * DT; _z  voxel.center.z – velo.z * DT; correct new position to stay in volume; tmp  find grid containing voxel(_x,_y,_z); if (!tmp) report error; convert coords(_x,_y,_z) to indices(i,j,k); btIndex  GetIDFromIndices(i,j,k); backtrack  tmp.voxel[btIndex]; neighbors[7]  find backtrack's nearest neighboring voxels according to coordinates (_x,_y,_z); if (voxel contains no vertices) { clone.velocity  interpolate velocity from neighbors; clone.density  interpolate velocity from neighbors; } } //---všechny nové hodnoty jsou již napočítány; //údaje z newList nahrajeme do voxelů swap(scene voxels, newList); //---duplikáty smažeme delete newList; }


49

Po dokončení advekčního kroku přichází na řadu krok projekční (metoda project), během něhož dojde ke změně proudového pole z divergentního na nedivergentní. project() { for each voxel in scene voxel.velDiv  compute velocity divergence; //---šíření rychlostí z listů do kořene (pullUp) propagateVelocities(this.octree); //---šíření tlaků z kořene do listů (pushDown) computePressure(this.octree); //---vytvoření nedivergentního proudového pole for each voxel in scene if (voxel's node has no model vertices) projectNewVelocities(voxel); //---odstranění mezivýpočtů z vnitřních uzlů clearGrid(); }

Rekurzivní metoda propagateVelocities se stará o převzetí rychlostí z potomků, jejich zprůměrování a zajištění výpočtu divergence příslušné části proudového pole. propagateVelocities(node) { leaves  0; if (node has children) { for each child {//---rekurzivní sestup do předposledních uzlů if (child has children) propagateVelocities(child); leaves++; } //---pokud jsme prošli všech 8 potomků if (leaves == 8) { //---index uzlu na nižší úrovni int cNode; int ind[8]; for each voxel in this.grid { //---zjisti indexy nižších buněk GetChildVoxelIndices(voxel,cNode,ind); voxel.velocity  0;

50

KAPITOLA 4 – IMPLEMENTACE //---převezmi a zprůměruj rychlosti for (i = 0..7) voxel.velocity += child.voxel[ind[i]].velocity /8;

} }

} this.computeDivergence(node.grid);

}

Pro výpočet divergence je volána samostatná metoda computeDivergence s parametrem (*grid), tedy ukazatelem na mřížku, ve které se má divergence počítat. Vypočtená divergence se uloží do buňky a nalezne uplatnění později, během řešení Poissonovy tlakové rovnice Jacobiho iterační metodou. computeDivergence(grid) { for each voxel in grid { divX  voxel.neighbor(+X).velocity voxel.neighbor(-X).velocity; divY  voxel.neighbor(+Y).velocity voxel.neighbor(-Y).velocity; divZ  voxel.neighbor(+Z).velocity voxel.neighbor(-Z).velocity; voxel.velDiv }



(divX+divY+divZ)/(2*grid.delta);

}

Poté, co jsme doručili všechny rychlosti do kořene stromu a v každé mřížce napočítali divergenci proudového pole, můžeme zahájit iterativní výpočet tlaku v rekurzivní metodě computePressure. Jelikož se tlak bude během pomyslné operace pushDown počítat opět pro všechny uzly na cestách od kořene do listů, přebírá tato metoda v parametru *node ukazatel na uzel, ve kterém má výpočet proběhnout. Pro iteraci byla zvolena a implementována Jacobiho metoda, která sice k řešení konverguje pomaleji, nicméně je implementačně jednodušší, než Gauss-Siedelova. Pro uspokojivé výsledky se doporučuje minimálně 40 iteračních

cyklů.

Tento

počet

je

v

aplikaci

řízen

konstantou

POISSON_ITERATIONS, nastavenou v průběhu vývoje na hodnotu 40.


51

computePressure(node) { for (iter = 0..POISSON_ITERATIONS) for each voxel in node's grid { //---sečti tlaky v sousedech; pokud soused //neexistuje, použij svoji hodnotu p  sum up pressures in voxel's neighbors; k  (node.grid.delta^2 * voxel.velDiv); voxel.pressure  (p-k)/6; } //---rekurzivně vypočítáme tlak v potomcích if (node has children) for each node's child computePressure(child); }

Po dokončení běhu metody computePressure máme ve výpočetních buňkách listů kromě původního (a s velkou pravděpodobností také divergentního) proudového pole napočítané i nové tlaky. V metodě projectNewVelocities provedeme v kontextu Helmholtz-Hodgeovy dekompozice odečtení gradientů těchto tlaků od původních rychlostí, čímž dostaneme nové, nedivergentní proudové pole, a tím dokončíme celý projekční krok simulace. projectNewVelocities(voxel) { //---výpočet tlakového gradientu v ose X gradX  voxel.neighbor(+X).pressure voxel.neighbor(-X).pressure; gradX  gradX / (2*grid.delta); //---výpočet tlakového gradientu v ose Y gradY  voxel.neighbor(+Y).pressure voxel.neighbor(-Y).pressure; gradY  gradY / (2*grid.delta); //---výpočet tlakového gradientu v ose Z gradZ  voxel.neighbor(+Z).pressure voxel.neighbor(-Z).pressure; gradZ  gradZ / (2*grid.delta); //---odečteme gradienty od rychlosti proudění voxel.addVelocity(-gradX*DT,-gradY*DT,-gradZ*DT); }

Po dokončení výpočtů všech simulačních kroků zbývá jediné – zduplikovat všechny výpočetní buňky a předat je vykreslovacímu modulu, běžícímu v samostatném vlákně. Vykreslování výsledků tak bude v ideálním případě probíhat současně s počítáním nového proudového pole. Duplikaci zajišťuje metoda duplicateVoxels. Jelikož ale do duplicitního pole dat budou najednou

52


chtít přistupovat oba dva moduly, je nutné zajistit jejich synchronizaci (jinak by mohlo dojít k tomu, že vykreslovací vlákno obdrží nová data, ale stále ještě nedokončila

vykreslování

předchozích).

Synchronizace

je

zajištěna

mechanismem kritické sekce – tedy chráněné oblasti, do které má v každý okamžik povolen přístup maximálně jeden modul. duplicateVoxels() { enter critical section; create clone space with the size of voxels count; for each voxel in scene { //---zkopíruj všechny hodnoty z voxelu clone  voxel; //---vyšetři kolize písku s objekty shoot voxel's grains towards models; } //---zkopíruj data terénu duplicate terrain; }

leave critical section;

4.1.2 Interakce proudového pole s pískovým podložím Pro ucelení přehledu o fyzikálním modulu nám ještě chybí popsat implementaci tří metod, zajišťujících interakci proudového pole s podložím – suspension (vzedmutí), saltation (unášení) a creep (splaz). Právě jim bude věnována tato část. Jejich princip byl již rozveden v kapitole 3.4.1 a proto se zaměřím především na implementační detaily. Metoda suspension přijímá jako parametry ukazatel na voxel, ze kterého čte vlastnosti větru a případně zapisuje novou hodnotu hustoty, a dále pak číselné hodnoty x a z, což jsou osové indexy výškového pole. suspension(voxel, x, z) { wind  fetch wind from voxel; if (wind points downwards) return; compute the angle between wind and ground; if (angle > threshold) { strength  compute wind strength;


}

53

decrease sand amount in heightfield(x,z); increase density in voxel;

}

Aby se nemuselo výškové pole procházet vícekrát, než je nutné, je metoda suspension volána z metody saltation, která jej rovněž prochází, čímž se jeden průchod ušetří. Metoda saltation využívá stejně jako advekce princip zpětného trasování. Nepřebírá žádné vstupní parametry, nicméně při své činnosti využívá dvě konstanty, které umožňují její škálovatelnost. Konstanta L0 ovlivňuje výslednou sílu větru a tím pádem délku kroku mezi aktuálním bodem a backtrackem. Druhá konstanta, q0, určuje, jak velké množství se má mezi danými dvěma body přenést. V závislosti na nastavení hodnot těchto konstant můžeme ovlivnit především velikost dun, vytvářených při unášení písku a také jejich rozestupy. saltation() { for each x,z in heightfield extents { //---vypočítej absolutní index z osových index  getIndex(x,z); //---pokud na dané buňce výškového pole stojí //model, přeskoč ji if (heightfield(index) lies under model) skip; //---zjisti výšku na dané pozici a najdi, do //kterého uzlu octree výškové pole zasahuje height  heightfield.getHeight(x,z); node  getNodeByPosition(x,height,z); if (node not found) skip; voxel  node.getVoxelByPosition(x,height,z); wind  voxel.windVelocity; if (wind points downwards or is too weak) skip; //---uprav velikost větru konstantou //(škálovatelnost dosahu) wind  wind*L0; //---najdi backtrack, ze kterého se odnese písek backtrack  heightfield(x,z)-wind; //---pokud na backtracku nestojí žádný objekt if (backtrack does not lie under model) { //---odeber z backtracku dané množství decrease backtrack's height by q0; //---přidej písek na aktuální pozici increase current area's height by q0; } suspension(voxel, x, z); } }

54


Poslední metodou fyzikálního modulu v našem výčtu je creep, provádějící odstraňování ostrých přechodů terénu. Je řízena konstantou pro maximální povolený gradient (MAX_GRADIENT) a počtem požadovaných průchodů (vstupní parametr iterations). creep(iterations) { for (i = 0..iterations) for each x,z in heightfield extents { current  heightfield(x,z); for each neighbor if ( (current not under model) and (neighbor not under model) ) { gr  compute gradient; if (gr > MAX_GRADIENT) { amount  gr/MAX_GRADIENT; decrease current's height by amount; increase neighbor's height by amount; } } } }

4.2 Datové struktury simulátoru Hybridní datová struktura popsaná v kapitole 2.2.2 intuitivně poukazuje na užití tří tříd. Jejich základní zjednodušený diagram je uveden na Obr.38.

Obr.38: Diagram zjednodušených tříd


55

Metody obsažené v uzlech stromu poslouží k předávání zpracovaných rychlostí z voxelů do vyšší, resp. nižší úrovně stromu. Tento mechanismus je využit při výpočtu tlakového členu Poissonovy rovnice (4).

4.2.1 Třída COctree - oktalový strom Nejdůležitější součástí simulátoru je třída COctree, reprezentující oktalový strom, použitý pro adaptivní dělení scény. class COctree { private: //---zjištění,zda se mají potomci sloučit bool checkMerging(); //---zjištění,zda se má uzel rozdělit bool checkDivision(); public: //---ukazatel na uzel rodiče COctree *parent;

//---ukazatele na potomky COctree *children[8]; //---identifikátor uzlu (0..7) unsigned short int index; //---hloubka uzlu ve stromu unsigned short int level; //---počet potomků (0 nebo 8) unsigned short int validChildren; //---minimální a maximální souřadnice uzlu tVector3 min,max; //---ukazatel na vnitřní uniformní mřížku SGrid *grid; //---seznam modelových vertexů, zasahujících do uzlu vector vertices; //---číslo iterace poslední opravy stromu int repairedAt; //---konstruktor pro kořen COctree(vector ver); //---konstruktor pro ostatní uzly COctree(COctree *parent, int index, vector ver); //---destruktor ~COctree(); //---vytvoření potomků pro daný uzel void createChildren(); //---vytištění cesty (indexů uzlů na větvi) od //kořene k aktuálnímu uzlu void printPath(bool endline=true); //---funkce pro nalezení sousedního stromového uzlu

56


};

//v požadovaném směru COctree* getNeighbor(int direction); //---nalezení všech předposledních uzlů stromu //(uzlů, jejichž všichni potomci jsou listy) void findSupLeaves(); //---nalezení všech listů stromu void findLeaves(); //---oprava stromu (rozdělení/sjednocení uzlů) void repairTree(); //---nalezení uzlu obsahujícího bod na dané pozici COctree* getNodeByPosition(float posX, float posY, float posZ, short int level = -1); //---test na obsažení daného bodu v aktuálním uzlu bool isInside(float x, float y, float z); //---vyzvednutí požadovaného voxelu z mřížky SVoxel* getVoxel(unsigned int index);

Jelikož je konstrukce kořene odlišná od ostatních uzlů (vnitřních a listů), využívá třída COctree dva konstruktory. //---konstruktor kořene COctree(vertexList) { set coordinates according to scene; create inner uniform grid; this.level  0; this.parent  NULL; for each child pointer this.child  NULL; this.validChildren  0; //---adaptivní dělení podle objektů ve scéně //---zjistíme, kolik vertexů zasahuje do uzlu for each vertex in vertexList { if (this.isInside(vertex)) this.vertices.push(vertex); } //---pokud uzel obsahuje vertexy, podrozdělí se if (this.vertices > 0) this.createChildren(); }

Zatímco do konstruktoru kořene jako parametr vstupoval seznam všech modelových vertexů ve scéně, do konstruktoru potomků je předáván vždy pouze seznam vertexů jeho rodiče.


57

//---konstruktor potomka COctree(*parent, index, vertexList) { compute child's coordinates according to its index; create inner uniform grid; this.parent  parent; this.level  parent.level + 1; this.index  index; for each child pointer this.child  NULL; this.validChildren  0;

}

//---adaptivní dělení podle objektů ve scéně //---zjistíme, kolik vertexů zasahuje do uzlu for each vertex in vertexList { if (this.isInside(vertex)) this.vertices.push(vertex); } //---pokud uzel obsahuje vertexy a jeho hloubka //nepřekračuje limit, uzel se podrozdělí if (this.vertices > 0 and this.level < MAXLEVEL) this.createChildren();

Předávání modelových vertexů má zásadní význam pro obtékání modelů, neboť uzly (společně se svými výpočetními buňkami), obsahující vertexy modelů jsou zneplatněné pro výpočty v advekčním kroku (viz kapitola 4.1.1). Konstantou MAXLEVEL, která určuje maximální povolenou hloubku oktalového stromu, se dá řídit přesnost dělení v okolí modelů (jak těsně bude model obepnut). Tato konstanta ale přímo ovlivňuje výpočetní náročnost simulátoru, a tudíž příliš vysoká hodnota (a tím pádem velmi vysoký počet buněk) bude představovat extrémní nároky na výpočetní čas. Pro úplnost se sluší uvést ještě definici metody createChildren(): createChildren() { this.validChildren = 8; for (index = 0..7) this.child[index]  }

COctree(this,index, this.vertices);

Funkce checkMerging a checkDivision na základě hustoty ve výpočetních buňkách a jejího gradientu určují, zda má ve stromě dojít ke sloučení, resp.

58


podrozdělení uzlů. Zatímco checkDivision pracuje s hustotami ve vlastním uzlu, její protějšek zkoumá hustotu v potomcích. Konečné rozhodnutí se provádí na základě porovnání získaných gradientů s limitní hodnotou T. Tato hranice je v obou funkcích rozdílná a je k ní rovněž přistupováno odlišnými způsoby. Je-li ve funkci checkDivision maximální nalezený gradient větší než mezní hodnota, provede se rozdělení a na výstup je předán příznak true. Nutno podotknout, že hledání gradientu nemusí vždy nutně proběhnout nad celou mřížkou. Pro realizaci rozdělení stačí první překročení mezní hodnoty. Funkce checkMerging provádí sloučení potomků ve chvíli, kdy je maximální gradient ze všech potomků nižší než nastavená mez. Hodnota obou mezí byla zjišťována experimentálně. Uspokojivé výsledky byly získány při Tmax = 10 a Tmin = Tmax/8 = 1,25. Podrozdělování stromu má za následek přítomnost větších detailů v tekutině, zatímco slučování zajistí, že v místě, kde již detaily nejsou potřeba (např. tam, kde již bouře přestala působit), dojde ke snížení rozlišení a tím i urychlení výpočtů a uvolnění paměti. checkDivision() { for each voxel in this.grid { find neighbors in three axes; correct differences in voxel levels; store neighbors' densities; find maximum gradient; if (maximum gradient > Tmax) { createChildren; //---propagace hustot do potomků for each voxel in this.grid { find 8 corresponding child voxels; set these children voxels' density to current voxel's density; } this.repairedAt  iteration; return true; } } return false; }

Uchovávání iterace, ve které došlo k poslední změně daného uzlu má velký význam pro funkci checkMerging. Může totiž dojít k situaci, kdy budou čerstvě


59

vytvořené voxely obsahovat malé rozdíly hustot a tudíž by byly touto funkcí vyhodnoceny jako nepotřebné a okamžitě by došlo k jejich sloučení. Funkce checkMerging a checkDivision by se tak dostaly do kolizní smyčky. Proto checkMerging nesmí slučovat uzly vytvořené v aktuální iteraci. V checkDivision není toto ošetření nutné, neboť se při běhu simulace volá před checkMerging a v dalším iteračním kroku jsou již hustoty ve voxelech přepočítány. checkMerging() { for each child node { childGrad  get maximum density gradient; if (childGrad > maxGrad) maxGrad  childGrad; } if (maxGrad < Tmin and this.repairedAt != iteration) { //---převzetí veličin z potomků for each voxel in this.grid { find 8 corresponding child voxels; accumulate children's density and velocity in current voxel; } mark child tree nodes for removal; this.validChildren  0; return true; } else return false; }

Metoda findSupLeaves zavádí definici tzv. „předposledních“ uzlů – vnitřních uzlů stromu, jejichž všichni potomci jsou listy (nejsou dále dělení). Znalost předposledních uzlů je nezbytná pro to, abychom byli schopni určit, nad potomky kterých uzlů je možné eventuálně provádět slučování za použití funkce checkMerging. Nicméně skutečnost, že je některý uzel předposlední, nám ještě definitivně nezaručuje, že jeho potomci mohou být sloučeni. V potaz je nutné brát také dělení stromu v blízkosti objektů – listy obsahující vertexy modelů rovněž nesmějí být předmětem slučování a ze seznamu uzlů pro sloučení musejí být odebrány.

60

KAPITOLA 4 – IMPLEMENTACE findSupLeaves() { if (this.validChildren == 0) return; //---sečteme vertexy v potomcích a „vnuky“ for each child { vertices += child.vertices; grandChildren += child.validChildren; } //---zkontrolujeme splnění obou podmínek if (vertices==0 and grandChildren==0) add this to supleaves list; else for each child child.findSupLeaves(); }

Ekvivalentní problém vzniká při hledání listů pro podrozdělování v metodě findLeaves. Identifikace listu je snadnější, než v předchozím případě, nicméně i zde je nutné z našich úvah vyloučit listy, do kterých zasahují modely. Takové uzly musejí zůstat neměnné po celou dobu simulace a není na nich tudíž možné provádět slučování ani dělení. findLeaves() { if (this.validChildren == 0) if (this.vertices==0) { add this to leaves list; return; } else for each child child.findLeaves(); }

O zajištění volání všech doposud popsaných metod a funkcí oktalového stromu se stará metoda repairTree. repairTree() { findLeaves(); while (leaves not empty) { leaf  leaves.pop(); leaf.checkDivision(); } findSupLeaves(); while (supLeaves not empty) {


}

61

supLeaf  supLeaves.pop(); supLeaf.checkMerging();

}

Funkce isInside(x,y,z) na základě porovnání zadaných souřadnic bodu s maximálními souřadnicemi uzlu určuje, zda zadaný bod patří do aktuálního uzlu. Výstupem je logická hodnota (true / false). if (

) else

x ϵ (this.min.x; this.max.x) and y ϵ (this.min.y; this.max.y) and z ϵ (this.min.z; this.max.z) return true; return false;

Její činnost je využívána ve funkci getNodeByPosition(posX,posY, posZ,maxLevel), která má za úkol vyhledat a vrátit nejhlubší uzel stromu obsahující daný bod a splňující zadané omezení hloubky. Tato funkce najde uplatnění například při backtrackingu v advekčním kroku. getNodeByPosition(posX,posY,posZ,maxLevel) { if (!this.isInside(posX,posY,posZ)) return NULL; if (this.level == maxLevel) return this; for each child node { if ( child.isInside(posX,posY,posZ) and child.level <= maxLevel ) return child.getNodeByPosition(posX,posY, posZ,level); } return NULL; }

Jednoduchá

funkce

getVoxel(index)

zajišťuje

zkrácený

přístup

do

požadovaného voxelu. Při běžném přístupu je nutné přistoupit přes uzel nejprve do mřížky a z ní teprve do daného voxelu (node.grid.voxel[index]). Tato funkce

provádí

přístup

(node.getVoxel(index)).

do

mřížky

implicitně,

stačí

tedy

volat

62


Funkce získání sousedního uzlu getNeighbor(direction) implementuje postup nastíněný v části 3.2.1. Následující část jejího pseudokódu funkce popisuje uvedený algoritmus pro směr +Y (tj. hledání horního souseda). … if (direction == +Y) { if (this.index >= 4) //---snadná dvojice return this.parent.child[this.index-4]; else { ptr  this; //---cesta vzhůru while (ptr.parent != NULL and flip==false) { if (nodeList not empty) //---porovnáme poslední a aktuální uzel; //pokud mají indexy z opačných množin, //nalezli jsme prapředka; nastavíme příznak if ((ptr.index<4 and nodeList.back>=4) or (ptr.index>=4 and nodeList.back<4)) flip  true; nodeList.push(ptr.index); ptr  ptr.parent; } if (flip) { //---pokud jsme nalezli prapředka (!root) //sestupujeme zpět dolů while (nodeList not empty) {//aplikujeme pravidla pro snadné dvojice if (nodeList.back < 4) { if (ptr.child[nodeList.back+4]!=NULL) ptr  ptr.child[nodeList.back+4]; else return ptr; } else { if (ptr.child[nodeList.back-4]!=NULL) ptr  ptr.child[nodeList.back-4]; else return ptr; } nodeList.pop(); } return ptr; } else //---žádná vyšší osmice neexistuje, //soused nenalezen


}

63

if (nodeList.back-4 < 0) return NULL;

} …

Paměťové nároky struktury: 128 B.

4.2.2 Struktura SGrid – uniformní mřížka Struktura uniformní mřížky je přítomna v každém uzlu oktalového stromu a obsahuje pevně stanovený počet stejně velkých výpočetních buněk (voxelů). V této implementaci je počet voxelů ve všech osách totožný a řízený konstantou GRID_SIZE (=8). Celkový počet buněk v každé mřížce je tedy 83 = 512 (konstanta GRID_SIZE_3). struct SGrid { //---ukazatel na uzel stromu COctree* node; //---seznam buněk vector<SVoxel*> voxels; //---minimální a maximální souřadnice mřížky tVector3 min,max; //---rozměr buněk float delta; //---úroveň mřížky unsigned short int level; //---konstruktory SGrid(); SGrid(tVector3 min, tVector3 max, ushort level); //---destruktor ~SGrid();

};

//---součet hustot všech buněk v mřížce float getDensity(); //---nalezení maximálního gradientu hustoty mezi //všemi buňkami v mřížce float getMaxDensityGradient(); //---rozhodnutí, zda voxel se středem na pozici //(x,y,z) je součástí této mřížky bool isInside(float x, float y, float z); //---nalezení voxelu podle pozice jeho středu SVoxel* getVoxelByPosition(float x,float y,float z);

64


Funkce getMaxDensityGradient() je použita pro nalezení maximálního gradientu hustoty mezi voxely dané mřížky. Znalost této hodnoty je důležitá pro správu adaptivní dělení, jelikož se podle ní určuje, ve kterém okamžiku se má příslušný uzel oktalového stromu rozdělit, resp. sloučit (viz kapitola 4.2.1). Při hledání sousedů může dojít k situaci s rozdílnými úrovněmi mřížek a je proto nutné provést korekci hodnot. Ta využívá předpokladu, že je-li aktuální buňka na hladině n a její soused na hladině n-1, bude rozdíl jejich hustot maximálně osminásobný (hladina n má 8x více buněk než hladina n-1). Tím pádem stačí hustotu v sousedovi vynásobit číslem 8hladina souseda−aktuální hladina . Pro výpočet gradientů v jednotlivých osách se následně použijí diskrétní tvary:

∇ x =  x1, y , z  x−1, y , z −2⋅ x , y , z ∇  y =  x , y1, z  x , y −1, z−2⋅x , y , z ∇  z =  x , y , z1x , y , z−1 −2⋅x , y , z Nakonec se z absolutních hodnot získaných gradientů vybere ta nejvyšší a je-li větší než dosud nalezené maximum, toto maximum se aktualizuje. getMaxDensityGradient() { highestGradient  0; for each voxel v in grid { neighbor[0]  getNeighbor(-AXIS_X, myLevel); neighbor[1]  getNeighbor(+AXIS_X, myLevel); neighbor[2]  getNeighbor(-AXIS_Y, myLevel); neighbor[3]  getNeighbor(+AXIS_Y, myLevel); neighbor[4]  getNeighbor(-AXIS_Z, myLevel); neighbor[5]  getNeighbor(+AXIS_Z, myLevel); for each neighbor n { if (neighbor[n] == NULL) dens[n]  0; else dens[n]  neighbor[n].density * correction(); } gradX  abs(dens[0]+dens[1] - 2*myDensity); gradY  abs(dens[2]+dens[3] – 2*myDensity); gradZ  abs(dens[4]+dens[5] – 2*myDensity); maxGrad  MAX3(gradX,gradY,gradZ); if (maxGrad > highestGradient) highestGradient  maxGrad; } return highestGradient; }


65

Funkce isInside(x,y,z) pracuje na totožném principu jako stejnojmenná funkce ve struktuře COctree (4.2.1). Na základě porovnání zadaných souřadnic středu buňky s maximálními souřadnicemi mřížky určuje, zda voxel s tímto středem patří do aktuální mřížky. getVoxelByPosition(x,y,z) – funkce, která zjistí, zda se v dané mřížce může vyskytovat voxel se středem na pozici (x,y,z). Jakmile je hledaný voxel součástí mřížky, jsou jeho osové indexy přemapovány na absolutní index uvnitř mřížky a je vrácen ukazatel na daný voxel. V opačném případě funkce vrací ukazatel nulový (NULL). Podmínkou úspěšného výpočtu indexu určení, zda se mřížka nachází v kladné či záporné části každé osy (současný výskyt v obou částech se nepředpokládá). getVoxelByPosition(x,y,z) { if (isInside(x,y,z)) { if (grid in positive X) indexX  abs(x-this.min.x)/delta; else indexX  abs(this.min.x-x)/delta; if (grid in positive Y) indexY  abs(y-this.min.y)/delta; else indexY  abs(this.min.y-y)/delta; if (grid in positive Z) indexZ  abs(abs(this.max.z)-z)/delta; else indexZ  abs(abs(z)-abs(this.max.z))/delta; index



getIDFromIndices(x,y,z);

return this.voxels[index]; } else }

return NULL;

Paměťové nároky struktury: 88 B.

66


4.2.3 Struktura SVoxel – výpočetní buňka Hlavním nositelem datové informace je buňka uniformní mřížky, umístěné uvnitř uzlu oktalového stromu. Právě nad těmito jednotlivými buňkami probíhají všechny výpočty fyzikálního modulu. struct SVoxel { //---hustota písku v buňce float density; //---směr a síla větru tVector3 windVelocity; //---tlak uvnitř buňky float pressure; //---divergence rychlosti větru float velDiv; //---minimální souřadnice buňky tVector3 min; //---maximální souřadnice buňky tVector3 max; //---pozice středu buňky tVector3 center; //---numerický identifikátor buňky int id; //---ukazatel na rodičovskou mřížkou SGrid* grid; //---vzdálenost středu buňky od kamery float distance; //---barevné složky buňky float r,g,b; //---seznam sousedních buněk vector<SVoxel*> neighbors; //---konstruktory SVoxel(); SVoxel(tVector3 &min, tVector3 &max); SVoxel(const SVoxel &input); //---destruktor ~SVoxel(); //---operátor přiřazení SVoxel& operator=(const SVoxel& input); //---změna rychlosti proudění void addVelocity(double x,double y,double z); //---získání celkové rychlosti z odpovídajících //buněk mřížky na nižší úrovni tVector3 getVelocity(); //---získání celkové hustoty z odpovídajících buněk //mřížky na nižší úrovni float getDensity(); //---hledání souseda v daném směru s omezením //maximální hloubky SVoxel* getNeighbor(int direction, ushort level); };


67

Poslední tři vyjmenované funkce si zaslouží bližší rozbor. Ve dvou z nich je použito vyzvedávání veličin z voxelů potomků, proto je vhodné si nyní popsat postup získání jejich indexů pomocí globální metody GetChildVoxelIndices, jejíž princip byl nastíněn v kapitole 3.2.2. Používá jeden vstupní a dva návratové parametry. Vstupním parametrem je index voxelu vyšší úrovně, prvním návratovým pak index oktalového uzlu nižší úrovně a druhým pole indexů hledaných voxelů. Rovnež tento pseudokód pracuje pro ilustraci s buňkou id=60 z Obr.34 a je doplněn o komentáře, ve kterých jsou po každé operaci vypsány hodnoty uložené v proměnných. GetChildVoxelIndices(inParentVoxelIndex, &retChildIndex, *retIndices) { halfGrid  GRID_SIZE/2; //Původní mřížka byla zmenšena (oříznuta) na 4x4x4 //buňky -> halfGrid = 8/2 = 4 (parentX,parentY,parentZ)  getIndicesFromID(id); //id=60 -> (parentX,parentY,parentZ) = (4,7,0) childX  parentX % halfGrid; //childX = 4%4 = 0 childY  parentY % halfGrid; //childY = 7%4 = 3 childZ  parentZ % halfGrid; //childZ = 0%4 = 0 //Ve zmenšené mřížce má buňka 60 osové indexy //(0,3,0). //Zjistíme, do které z nových mřížek se dostal objem //buňky 60. Mřížky jsou indexovány stejným způsobem //jako voxely (od 0 do 7, inkrementace nejprve v ose //X, poté v Y a nakonec v záporné části osy Z). childGrid  parentX/halfGrid + parentY/halfGrid * 2 + parentZ/halfGrid * 4; //childGrid = 4/4 + 7/4 *2 + 0/4 * 4 = 1+2+0 = 3 //Původní buňka 60 je součástí mřížky č.3. //Indexace mřížek je ovšem odlišná od indexace uzlů //oktalových stromů (viz Obr.35). Musím tedy ještě //zjistit, ve kterém stromovém uzlu se daná mřížka //nachází. retChildIndex  GridIndexToTreeIndex(childGrid); //retChildIndex = 3 //Nyní zbývá vypočítat osové indexy 8 nových voxelů retIndices[0]  2*childX + 2*childY*GRID_SIZE + 2*childZ*GRID_SIZE*GRID_SIZE; //retIndices[0] = 2*0 + 2*3*8 + 2*0*8*8 = 48 retIndices[1]  2*childX+1 + 2*childY*GRID_SIZE +

68

KAPITOLA 4 – IMPLEMENTACE 2*childZ*GRID_SIZE*GRID_SIZE; //retIndices[1] = 2*0+1 + 2*3*8 + 2*0*8*8 = 49 retIndices[2]  2*childX + (2*childY+1)*GRID_SIZE + 2*childZ*GRID_SIZE*GRID_SIZE; //retIndices[2] = 2*0 + 7*8 + 2*0*8*8 = 56 retIndices[3]  2*childX+1 +(2*childY+1)*GRID_SIZE + 2*childZ*GRID_SIZE*GRID_SIZE; //retIndices[3] = 2*0+1 + (2*3+1)*8 + 2*0*8*8 = 57 retIndices[4]  2*childX + 2*childY*GRID_SIZE + (2*childZ+1)*GRID_SIZE*GRID_SIZE; //retIndices[4] = 2*0 + 2*3*8 + (2*0+1)*8*8 = 112 retIndices[5]  2*childX+1 + 2*childY*GRID_SIZE + (2*childZ+1)*GRID_SIZE*GRID_SIZE; //retIndices[5] = 2*0+1 + 2*3*8 + (2*0+1)*8*8 = 113 retIndices[6]  2*childX + (2*childY+1)*GRID_SIZE +(2*childZ+1)*GRID_SIZE*GRID_SIZE; //retIndices[6]= 2*0 + (2*3+1)*8 + (2*0+1)*8*8 = 120 retIndices[7]  2*childX+1 + (2*childY+1)*GRID_SIZE + (2*childZ+1)*GRID_SIZE*GRID_SIZE; //retIndices[7]= 2*0+1+(2*3+1)*8+(2*0+1)*8*8 = 121 }

getVelocity() - Funkce zjišťující rychlost proudění v daném voxelu (existují-li voxely na dané pozici i v nižší úrovni stromu, funkce pracuje rekurzivně). Pomocí předchozí metody zjišťuje potřebné indexy a z daných voxelů vyzvedne rychlosti, které jsou nasčítány a nakonec zprůměrovány. getVelocity() { int childNode; int index[8]; GetChildVoxelIndices(id, childNode, index); //Poslední krokem je vyzvednutí rychlostí z těchto //potomků a jejich zprůměrování. value  (0,0,0); for (j=0..8) value+= node[childNode].voxel[index[j]].velocity/8; return value; }

Totožný algoritmus pak využívá také funkce getDensity() pro získání hustot. Funkce getNeighbor(direction, level) slouží k vyhledávání sousedního voxelu ve směru direction, reprezentovaném jedním ze tří identifikátorů (AXIS_X, AXIS_Y, AXIS_Z), přičemž negace jednotlivých identifikátorů (-AXIS_X, -AXIS_Y, -AXIS_Z) reprezentují záporné části daných os. Na základě


69

parametricky předaného identifikátoru směru se z pozice středu aktuálního voxelu vypočítá střed hledané buňky (např. při použití identifikátoru AXIS_X dojde k přičtení kroku prostorového dělení k x-ové souřadnici středu). Následně se testuje, zda voxel s takto získanou pozicí středu je součástí aktuální mřížky. Pokud ano, je hledaný soused ihned vrácen. O něco složitější situace nastane v okamžiku, kdy vypočítaný střed spadá do jiné než aktuální mřížky. Tehdy se musí traverzovat oktalový strom a na základě požadované pozice vrátit novou mřížku, ze které se získá požadovaný voxel. Zde přichází na řadu druhý parametr, level, neboli maximální úroveň (hloubka stromu), do které se má traverzace provádět. getNeighbor(direction,level) { newCenter  this.position + direction; if (newCenter lies in this grid) return (voxel with newCenter); else { while (newCenter in octree) traverse octree if (octree is NULL) return NULL; else return (octree.grid.voxel with newCenter); } }

Samozřejmě ani při omezení hloubky není zaručeno, že sousední buňka bude ležet ve stejné hladině jako výchozí. Strom totiž může být v dané větvi rozdělen méně a výsledná hloubka tak může být nižší. V tuto chvíli ale ještě není potřeba činit žádné srovnávání, neboť příslušná opatření pro vyrovnání úrovní se provedou až v okamžiku samotného zpracování sousedů. Paměťové nároky struktury: 160 B.

4.3 Plynulost běhu aplikace Pro efektivnější běh simulátoru je vhodné použít paralelizaci. Problémem je však práce nad stejnými daty – výstup jednoho výpočetního kroku je totiž vstupem pro krok následující, což požaduje dodržení jejich sekvence. Paralelní

70


zpracování tak najde uplatnění pouze mezi modulem pro fyzikální výpočty a modulem vykreslovacím, kdy vykreslovací vlákno pracuje nad daty z iterace i, zatímco fyzikální modul již provádí výpočty v iteraci i+1.

KAPITOLA 5 – LADĚNÍ

71

Kapitola 5 – Ladění Finální

fáze

vývoje

simulátoru

byla

spojena

s

optimalizací

naimplementovaných operací. Cílem bylo co nejvíce snížit nároky nejen na výpočetní čas, ale také na operační paměť.

5.1 Podpůrné aplikace Pro získání potřebného náhledu na tyto parametry byly použity dvě monitorovací aplikace – MemoryValidator pro kontrolu paměťových nároků a Very Sleepy, profiler sledující četnost využívání určitých kusů programového kódu.

5.1.1 MemoryValidator MemoryValidator [26] je produktem britské společnosti Software Verification Limited. Ta toho času nabízí celkem osm převážně komerčních1 testovacích aplikací pro nejpoužívanější programovací jazyky a operační systémy. MemoryValidator představuje komplexní nástroj pro sledování paměti po dobu běhu testované aplikace. Mezi její schopnosti patří kromě analýzy aktuálního množství paměťového prostoru, alokovaného pro aplikaci, také např. detekce přetečení a podtečení zásobníku/haldy, nesprávných dealokací či neinicializovaných dat. V daném místě nabízí také možnost zobrazení odpovídacího zdrojového kódu. Pro nastavení spolupráce testovací a testované aplikace je k dispozici pět různých možností: ●

„CreateProcess“ - vytvoření procesu. Při použití této doporučené volby je hlavní vlákno testované aplikace pozastaveno a její spouštěcí kód je upraven tak, aby umožnil korektní navázání (injection) aplikace MemoryValidator. Po této inicializaci je testované aplikaci opět povolena původní činnost.

1 Ke všem nabízeným komerčním aplikacím je poskytována 30-denní bezplatná zkušební verze. Cena plné verze produktu MemoryValidator 4.95 byla v době psaní této práce (prosinec 2010) $299.

72

KAPITOLA 5 – LADĚNÍ ●

„Normal“ - Testovaná aplikace je spuštěna běžným způsobem a s normální prioritou a MemoryValidator se připojí k jejímu procesu. Tato volba není vhodná pro použití u malých aplikací, jejichž start probíhá velmi rychle a tudíž nemusí dojít ke včasnému svázání a korektnímu sběru všech dat.

●

„Idle“ - Princip této metody je velmi podobný předchozímu. Jediným rozdílem je spuštění testované aplikace s nízkou prioritou, což by mělo zajistit více času pro svázání procesů.

●

„Suspend“ - Hlavní vlákno testované aplikace je potlačeno a dojde k injekci MemoryValidatoru. Korektní funkčnost této metody ovšem není zaručena u všech aplikací.

●

„Paused“ - Při startu testované aplikace dojde k pozastavení všech jejích vláken a následné injekci MemoryValidatoru. Ani tato metoda však nemusí vždy fungovat

Obr.39: Úvodní obrazovka aplikace MemoryValidator


73

Na Obr.40 je uveden výřez jedné z mnoha záložek aplikace MemoryValidator v průběhu sběru dat z připojeného simulátoru písečné bouře. U každého datového objektu (Type) je graficky i číselně uvedena jeho velikost (Size), aktuální (Count), maximální (Max) a celkový (Cumulative) počet alokovaných instancí v průběhu testování. Předposlední sloupec (% Size) udává v daný moment podíl součtu velikostí všech alokovaných instancí daného datového objektu na celkové alokované operační paměti aplikace. Podobně sloupec (% Object) udává poměr počtu instancí daného objektu ku všem aktuálně naalokovaným objektům.

Obr.40: Výřez okna aplikace MemoryValidator v průběhu testování

S pomocí aplikace MemoryValidator bylo odhaleno několik nedostatků, z nichž nejzávažnějším byla nejprve nesprávná (a později neúplná) dealokace výpočetních buněk (struktura SVoxel). Opravou této chyby došlo ke snížení paměťových nároků na přijatelnější hodnoty (1,5 GB vs. 400 MB na dané scéně).

74


5.1.2 Very Sleepy Druhým použitým testovacím nástrojem byla open source aplikace Very Sleepy (verze 0.7) [27]. Jde o profilovací program pro aplikace napsané v jazyce C++, který byl vytvořen modifikací původního profileru Sleepy [28] (oproti svému předchůdci například umožňuje sledovat více vláken současně či exportovat a importovat výsledky). Princip takového nástroje je podobný jako u aplikace MemoryValidator, avšak namísto s operační pamětí operuje s hodnotami registru instrukcí (instruction pointer register).

Obr.41: Hlavní okno aplikace Very Sleepy. Nahoře výpis aktivních systémových procesů, v dolní části pak výpis sledovatelných vláken.

Profiler je spuštěn v samostatném vlákně a ve velmi krátkých časových intervalech (~1ms) přerušuje běh testované aplikace a z jejího kontextu čte hodnoty registru instrukcí. Získané hodnoty adres jsou za pomoci kompilačních údajů převedeny na čísla řádků a názvy příslušných instrukcí. Tyto informace jsou při každém přerušení uchovány a po dokončení testu je pak vytvořen jejich detailní statistický přehled (viz Obr.42). V levé horní části je výpis operací, doplněný o absolutní i relativní množství času, které v nich aplikace při běhu strávila. Při kliknutí na některou z uvedených operací se aktivují zbylé tři části okna. V levém dolním výstupu je vypsán zdrojový kód příslušející dané operaci.


75

V pravé části jsou pak údaje o tom, ze které metody a jak často byla daná operace volána (nahoře), zatímco dolní polovina obsahuje výpis (a podíl) metod volaných.

Obr.42: Výsledky profilování simulátoru písečné bouře

Rovněž použití této aplikace mělo zásadní vliv na finální úpravy zdrojového kódu. Bylo pomocí ní zjištěno, že velká část běhu (cca 30%, viz Obr.42) je strávena ve funkci vykreslovacího modulu pro výpočet fáze osvětlení. Trasování ukázalo, že téměř celý čas, strávený v této funkci, připadá na normalizaci vektoru. To nám názorně předvedlo, jak časově náročné jsou i na dnešním hardwaru době operace odmocňování a dělení. V místech, kde to bylo možné, tak bylo dělení převedeno na násobení (x/2  x*0.5), případně předpočítáno, což mělo zásadní význam například v cyklech, kde se v každé iteraci dělilo konstantou. Podíl pak stačilo jednorázově vypočítat před zahájením cyklu a na dále pracovat přímo s výsledkem. Co se týče odmocňování - jelikož přesná hodnota výsledku nebyla prakticky na žádném místě kriticky důležitá, byla tato operace nahrazena několikanásobně rychlejší logaritmickou aproximací. Tyto dvě úpravy následně přinesly zrychlení výpočtů zhruba o 60%.

76


5.2 Zpracování výsledků testování Vhodným měřítkem pro vyhodnocení simulace je porovnání závislosti doby běhu jednotlivých kroků řešení NSR na složitosti dělení scény.

Obr.43: Graf závislosti doby běhu fyzikálního modulu na složitosti scény

Závislosti všech tří kroků jsou v souladu s očekáváním lineárně závislé na počtu voxelů. Nejdelší čas je věnován projekčnímu kroku (zeleně), ve kterém dochází nejen k dvojnásobnému průchodu celým stromem scény, ale také k iterativnímu řešení

tlakové

rovnice

v

kořenu

stromu.

Vliv

konstanty

POISSON_ITERATIONS (zde nastavené na hodnotu 40) na délku výpočtu bude rovněž lineární – čím více iterací provedeme, tím déle nám to bude trvat. Je ale patrné, že již jen kvůli délce běhu projekčního kroku nemůže běžet simulace v reálném čase. Další velký díl z celkového času běhu si pak odebere ještě vykreslovací část simulátoru.


77

Obr.44: Graf vztahu mezi alokovanou pamětí a počtem voxelů. Počáteční strmý nárůst je způsoben inicializací pole duplikovaných voxelů.

78

KAPITOLA 6 – VÝSLEDKY

79

Kapitola 6 – Výsledky V této kapitole odprezentuji a zhodnotím výsledky různých částí implementovaného simulátoru. Obrázky, ilustrující vývoj bouřkového mračna, byly vytvořeny s pomocí vizualizačního modulu ([1]).

6.1 Adaptivní dělení scény Na Obr.45 jsou uvedeny výsledky adaptivního dělení prostoru v okolí jednoho objektu. Pro lepší ilustraci je použito více pohledů na stejnou scénu. Maximální hloubka stromu je nastavena na hodnotu 4 (nastavení na vyšší hodnotu by výslednou scénu značně znepřehlednilo).

Obr.45: Přizpůsobení oktalového stromu objektu ve scéně.

Z ukázky je zřejmé, že v případě jednoho objektu probíhá adaptivní dělení podle očekávání. Nyní je nutné ověřit, zda stejně kvalitních výsledků dosáhneme i v případě použití většího počtu objektů. Výsledek (opět s maximální hloubkou dělení 4) je uveden na Obr.46 a Obr.47.

80


Obr.46: Přizpůsobení adaptivního dělení většímu počtu objektů ve scéně

Obr.47: Adaptivní dělení složitější scény

Také v tomto případě proběhlo dělení scény v pořádku. Scéna se ale musí správně adaptivně dělit i v oblastech výskytu bouřkového mračna, aby se ukázalo, zda jsou správně naimplementovány funkce pro dělení a slučování uzlů v místech, splňujících podmínku maximálního gradientu. Patrně nejlepší vypovídací hodnotu bude mít tento test v případě, že budeme na scénu pohlížet shora. Výsledek je vyobrazen na Obr.48.


81

Obr.48: Průběh podrozdělování scény v místech výskytu bouřkového mračna (pohled shora)

Z uvedené sekvence je patrné, že se scéna adaptivně dělí správně pouze v oblastech aktuálního výskytu bouře. Oblasti, do kterých ještě nedorazila, jsou, stejně jako oblasti, ve kterých již přestala působit, reprezentovány větším prostorem (méně rozdělenými uzly). Důležitým prvkem simulace je obtékání bouřkového mračno okolo objektů. Těsnost obtékání je závislá na maximální povolené hloubce dělení oktalového stromu. Čím větší bude tato hodnota, tím těsněji bude mračno objekty obtékat, avšak tím více také vzrostou nároky na paměť a výpočetní čas celé simulace.

82


Obr.49: Vliv maximální hloubky dělení stromu na obtékání modelů. Červeně je vyznačena orientační hranice dělení. Vlevo: h=3, uprostřed: h=4, vpravo: h=5

Na Obr.49 je vliv maximální hloubky dělení dobře patrný (hranice dělení je zvýrazněna červenou lomenou čárou). Zatímco při hloubce h=3 je obtékání velmi hrubé, po zvýšení o pouhou jednu další úroveň je již znatelně těsnější. Přidání další úrovně způsobí další zjemnění obtékání a lepší přizpůsobení mračna tvarům objektů. Dobře patrné je uvolnění prostoru mezi hlavou slona a podložím. Jelikož mají výpočetní buňky tvar krychle, dojde k nejlepšímu (nejtěsnějšímu) obtékání u málo členitých, osově zarovnaných objektů.

6.2 Proudové pole Jedním ze způsobů, jak ověřit chování proudového pole a zjistit aktuální směry jeho proudění, je použití tzv. proudových čar (streamlines). Jde o záznam kontinuálního pohybu částice vpuštěné do proudového pole a vzorkování jejího pohybu – celková trasa, kterou částice urazí je ve výsledku znázorněna jako lomená čára. Záznam pohybu několika takových částic je uveden na Obr.50.


83

Obr.50: Proudové čáry; vlevo pohled shora, vpravo pohled z boku

Na scéně bez objektů nevykazují proudové čáry žádné anomálie, proudové pole se pohybuje jedním směrem. Otázkou je, jaký výsledek dostaneme, když do scény přidáme objekty. Tato situace je znázorněna na Obr.51, kde dochází k průchodům proudových čar skrz objekty. Tento výsledek nesplňuje očekávání. Je ale způsoben tím, že proudové čáry jsou tvořeny dopředně, tzn. v aktuální buňce se zjistí směr větru a v tomto směru se částice vyšle, bez ohledu na to, zda cestou cokoliv protne. Naopak advekce, jak již bylo zmíněno výše, využívá zpětné trasování, při kterém se do buňky obsazené modelem nic dostat nemůže. Během simulace se tak proudění v těsné blízkosti zneplatněných buněk zastaví a pokračuje v okolí modelů (viz podkapitola 3.3.1). Ke stejnému problému dochází i na úrovni podloží, kdy proudové pole bohužel nerespektuje případné vysoké duny (neobtéká je).

84


Obr.51: Chybový stav, při kterém proudové čáry protínají objekty.

Průchod proudových čar objekty lze napravit zkrácením časového kroku, používaného při jejich vytváření. V takovém případě dojde v souladu s proudovým polem k jejich pozastavení v okolí objektů (Obr.52).

Obr.52: Demonstrace zastavení proudového pole v okolí objektu


85

6.3 Bouřkové mračno Vývoj bouřkového mračna je vhodné sledovat nejprve na volné (resp. velice málo obsazené) scéně. Jedna taková scéna společně s několika částmi sekvence vývoje bouře je uvedena na Obr.53.

Obr.53: Vývoj bouřkového mračna

Z obrázku je zřejmé, že vývoj bouře na volném prostranství probíhá dle očekávání. Dokud mračno nenarazí na protilehlý okraj scény, má podobu odpovídající reálnému mraku. V okamžiku protnutí okraje simulace správně působí dojmem, že bouře pokračuje i mimo scénu. Průběh vývoje bouře na scéně s jedním větším objektem můžeme pozorovat na Obr.54. Je zde patrné zachování tvaru i při obklopení objektu, stejně jako zmiňované pozastavení části mračna v jeho těsné blízkosti.

86


Obr.54: Průběh bouře ve scéně s jedním větším modelem

Jak již bylo zmíněno, vzhled a průběh obklopení je závislý na členitosti povrchu objektu, potažmo hloubce stromové struktury pro prostorové dělení. Celkový vizuální dojem simulace této jednoduché scény je poměrně uspokojivý. Na řadu přichází test simulace na třech různých scénách. První a nejjednodušší scéna obsahuje tři stromy a dva různě velké modely slona. Výchozí podoba je uvedena v první části Obr.55. Prostřední část reprezentuje stav v osmém kroku simulace a je na ní viditelné ovlivnění mračna prvním objektem (tato část je v plné velikosti vyobrazena na Obr.56). Poslední snímek zachycuje výsledek osmnáctého simulačního kroku, kdy již bouře postoupila ke slonům a začala je obtékat (snímek v původní velikosti je na Obr.57).


87

Obr.55: Testovací scéna č.1

Obr.56: Prostřední část Obr.54 v plné velikosti

88


Obr.57: Osmnáctý krok simulace první scény v plné velikosti

Druhá testovací scéna byla složena z většího počtu různě velkých a různě tvarovaných objektů, rozmístěných v různých částech scény. Stejně jako v předchozím případě je výchozí scéna znázorněna v první části sekvence, zbylé dvě reprezentují čtrnáctý, resp. dvacátý krok simulace.


Původní velikost snímku v prostřední části kompozice je uvedena na Obr.59. Je na něm dobře patrné obtékání budov a zejména dřevěných klád (v detailu na Obr.60).


89


Obr.60: Detailní pohled na obtékání mračna okolo dřevěných klád (3x zvětšeno)

Poslední, třetí scéna představuje část moderního města. Na Obr.61 jsou kromě výchozí podoby scény zachyceny také výsledky dvanáctého a pětadvacátého kroku simulace.


90


Výsledky této scény působí nejvíce realisticky. Na Obr.62, který představuje výsledek dvanáctého simulačního kroku v plné velikosti, je velmi dobře patrný rozdíl mezi obtékáním objektu v podobě osově zarovnaného kvádru (výšková budova) a objektu se složitějším tvarem (Eiffelova věž) při hloubce stromu h=3. Zatímco výšková budova je obklopena velmi těsně, v okolí věže jsou na mraku viditelné „zuby“.


Pozn.: Kompletní výstupní sekvence simulace všech tří testovací scén jsou obsaženy na přiloženém datovém nosiči.

6.4 Podloží Podloží je ovlivněno třemi operacemi – vzedmutí, unášení a splaz, přičemž vliv poslední jmenované je nejlépe viditelný. Bez vyhlazování by se na povrchu tvořily nevzhledné a rušivé špičky. Účinek této operace je znázorněn na Obr.63, který porovnává scénu bez vyhlazování a s vyhlazováním.


91

Obr.63: Vlevo terén bez vyhlazování, vpravo terén zpracovaný metodou creep

Jelikož metoda splazu (creep) byla implementována jako iterativní, je vhodné uvést také vliv počtu iterací na účinek vyhlazování. Tyto výsledky jsou viditelné na Obr.64, který účinnost rostoucí s počtem iterací potvrzuje a tím potvrzuje také správnost implementace této operace.

Obr.64: Detail terénu a vliv počtu iterací na výsledky vyhlazovací metody creep. a) bez iterací, b) 1 iterace, c) 5 iterací, d) 10 iterací, e) 20 iterací, f) 50 iterací

Sama o sobě by ovšem tato metoda neměla v průběhu simulace žádný smysl. Její funkce je totiž závislá na obou zbývajících, tedy vzedmutí a unášení, které

92


způsobují rozdíly ve výškách mezi jednotlivými body výškového pole. Suspense (vzedmutí) je nejvíce patrná ve chvíli, kdy v celém proudovém poli povolíme nenulové rychlosti. To má za následek zdvihání písku i mimo oblast aktuálního působení bouře. Jelikož je však v okolí bouře rychlost větru prakticky nulová, vzedmutí se projeví pouze v těsné blízkosti bouřkového mraku. Oba stavy jsou znázorněny na Obr.65 (z pravé části bylo pro lepší zřetelnost odstraněno podloží).

Obr.65: Příklady vzedmutí písku v závislosti na proudovém poli


Obr.66: Výsledná scéna po spojení fyzikálního a grafického modulu

93

94


KAPITOLA 7 – ZÁVĚR

95

Kapitola 7 – Závěr Cílem této práce měl být simulátor písečné bouře, založený na chování proudového pole fyzikálně popsaném Navier-Stokesovými rovnicemi. Důraz byl kladen především na adaptivní dělení prostoru, a to jak v okolí statických objektů, tak samotné bouře. Pro větší realističnost mělo proudové pole ovlivňovat také podloží scény. Úspěšně jsem provedl implementaci adaptivního dělení a její správnost jsem ověřil na několika různých scénách. Úspěšně byla provedena také implementace řešení Navier-Stokesových rovnic pro účely fyzikálně založené simulace písečné bouře včetně působení bouře na písečné podloží. Jedinou slabinou je obtékání dun a objektů, jehož implementace naráží na problém určení správného směru pohybu na základě povrchů daných objektů. Ve stávajícím provedení je obtékání řešeno pozastavením bouře v blízkosti objektu a jeho vizuální podoba je silně závislá na hloubce stromu prostorového dělení. I přesto je ale modul schopen provádět věrohodnou simulaci pohybu a vývoje mračna písečné bouře a spojením s grafickým modulem tak vznikla finální aplikace schopná modelovat vznik, vývoj a chování bouřkového mračna včetně vizuálně reálného osvětlení celé scény. V samotném závěru diplomové práce byly uvedeny grafické výsledky implementovaného simulátoru, vykreslené s pomocí modulu pro vizualizaci písečných bouří Bc. Lukáše Orešky. Při dalším rozšiřování stávajícího fyzikálního modulu by bylo užitečné vylepšit zejména chování mračna v blízkosti překážek – tedy vhodně upravit obtékání tak, aby i u velmi členitých objektů probíhalo co nejtěsněji a přitom bez zvýšených nároků na výpočetní čas. Obtékání by také mělo brát lépe v potaz tvary a povrchy objektů a přizpůsobit se jim. Námětem na další práci by mohlo být rovněž nalezení a implementace vhodného způsobu, který by dokázal simulaci (včetně následného vykreslování) provádět rychleji a plynuleji.

96

POUŽITÁ LITERATURA A ZDROJE

97

Použitá literatura a zdroje [1] Lukáš Oreška – Vizualizace písečných bouří. Diplomová práce, ČVUT FEL, 2011 [2] Deserts - Causes of deserts, Desert landforms, Deserts in the past http://science.jrank.org/pages/47498/deserts.html , verze z 7.12. 2010 [3] Q. Zhang, X. Zhao, Y. Zhang, L. Li - Preliminary study on sand-dust storm disaster and countermeasures in China, Chinese geographical science, vol. 12, no. 1, pp. 9-13, Science Press, Beijing, China, 2002 [4] Global Alarm: Sand and Dust Storms from the World's Drylands, United Nations Convention to Combat Desertification, 2002, http://www.unccd.int/publicinfo/duststorms/part1-eng.pdf , verze z 10.12.2010 [5] Swap, R. et al., Saharan dust in the Amazon basin. 1992, Tellus, 44B, 2: 133-149 [6] BBC online - http://news.bbc.co.uk/2/hi/sci/tech/1595338.stm , verze z 10.12.2010 [7] L.M. Milne-Thomas – Theoretical Aerodynamics. Dover Publications, 1973. ISBN 048661980X. [8] J. Hess, A.M.O. Smith – Calculation o f Potential Flow About Arbitrary Bodies. Progress in Aeronautics Sciences 8: 1-138, 1967 [9] R. Bridson – Fluid Simulation for Computer Graphics. A K Peters, Ltd., 2008. ISBN 1568813260. [10] N. Foster, D. Metaxas – Realistic Animation of Liquids. Graphical Models and Image Processing, 1996 [11] J.X. Chen, N. da Vittoria Lobo, C.E. Hughes, J.M. Moshell – Real-Time Fluid Simulation in a Dynamic Virtual Environment. IEEE Computer Graphics And Applications, 1997. [12] N. Foster, D. Metaxas – Modeling the Motion of a Hot, Turbulent Gas. Computer Graphics Proceedings, Annual Conference Series, 1997 [13] J. Stam - Stable Fluids. In SIGGRAPH 99 Conference Proceedings, Annual Conference Series, pages 121-128, srpen 1999. [14] S. Liu, Z. Wang, Z. Gong, L. Huang, Q. Peng – Physically Based Animation of Sandstorm. Computer Animation and Virtual Worlds, Volume 18, Number 4-5: 259-269, 2007 [15] L. Oreška, A. Stefanidis, D. Eršil – Physical Implementation of sandstorm based on two-fluid dynamic model. České vysoké učení technické v Praze, 2008 [16] F. Losasso, F. Gibou, R. Fedkiw – Simulating water and smoke with an octree data structure. ACM Transactions on Graphics 23, 3, 2004 [17] Y. Zhu, R. Bridson – Animating Sand as Fluid. ACM Transactions on Graphics. Proceedings of ACM SIGGRAPH 2005 [18] W. Honk, D.H. House, J. Keyser – Adaptive Particles for Incompressible Fluid Simulation. Technical report, Texas A&M Computer Science, 2007 [19] R. Ďurikovič, M. Kupka - VOF Method for Fluids and Solids on Octree Structure, pp.76-81, 2009 Second International Conference in Visualisation, 2009 [20] L. Shi, Y. Yu – Visual Smoke Simulation with Adaptive Octree Refinement, Computer Graphics and Imaging: Seventh IASTED International Conference Proceedings, pages 13-19, 2004 [21] Wikipedie, otevřená encyklopedie – Navier-Stokesova rovnice, http://cs.wikipedia.org/wiki/Navierova-Stokesova_rovnice , verze z 16.12.2010 [22] R. Culver – Numerical Solution of the Poisson Equation. Department of

98

POUŽITÁ LITERATURA A ZDROJE

Mechanical and Aerospace Engineering, University of Carolina, Irvine, 2005 [23] Newton forum - http://www.newton.dep.anl.gov/askasci/gen06/gen06342.htm , verze z 16.12.2010 [24] A. Nealen – Physical Based Simulation and Animation of Gaseous Phenomena in a Periodic Domain. Project report, Department of Computer Science, University of British Columbia, 2002 [25] A.J. Chorin, J.E. Marsden – A Mathematical Introduction to Fluid Mechanics. 3rd ed., Springer, 1993, ISBN 0387979182 [26] Stránky programu MemoryValidator - http://www.softwareverify.com/cpp/memory/ [27] Stránky programu Very Sleepy - http://www.codersnotes.com/sleepy [28] Stránky programu Sleepy - http://sleepy.sourceforge.net/

Zdroje obrázků [O2] Daily Mail Online http://www.dailymail.co.uk/news/article-452050/Return-The-Mummy.html [O3] Science Encyclopedia - http://science.jrank.org/pages/47498/deserts.html [O4] NASA Science News http://science.nasa.gov/science-news/science-at-nasa/2001/ast26jun_1/ [O5] Greg Jasion, Computational Modelling Group, University of Southampton http://cmg.soton.ac.uk/people/gtj104/ [O6] Wikipedia, the free encyclopedia - http://en.wikipedia.org/wiki/Dust_storm [O7] Astronet - Martian Dust Storm - http://www.astronet.ru/db/xware/msg/1169949 [O8] Space Today - Exploring Mars http://www.spacetoday.org/SolSys/Mars/MarsThePlanet/MarsDustStorms.html [O20] Miss Baker's Biology Class Wiki http://missbakersbiologyclasswiki.wikispaces.com/Test+Study+Guide [O21] Pling! Computer Graphics and Visualisation http://www.pling.org.uk/cs/cgv.html [O22] nVidia Developer Zone http://http.developer.nvidia.com/GPUGems3/gpugems3_ch22.html [O24] PowRay Wiki http://wiki.povray.org/content/Documentation:Tutorial_Section_3.1 [O36] Fiddlers Green Paper Models http://www.fiddlersgreen.net/models/aircraft/Mitsubishi-Claude.html [O37] Matlab Central http://www.mathworks.com/matlabcentral/fx_files/4674/1/Immagine.JPG

PŘÍLOHA

99

Příloha A – Interpolace v advekci Pro získání nové hodnoty (ať již hustoty nebo rychlosti) se při zpětném trasování v advekčním kroku používá následující interpolace: newValue = (1 - r1) * ( (1 - r2) * + r2 *

( (1 – r3) * backtrack + r3 * neighbor Z ) ( (1 – r3) * neighbor Y + r3 * neighbor YZ )

(1 - r2) * + r2 *

( (1 - r3) * neighbor X + r3 * neighbor XZ ) ( (1 - r3) * neighbor XY + r3 * neighbor XYZ )

) + r1* (

)

Indexy proměnné neighbor určují, v jakých osách se máme pohybovat při přistupování do souseda. Konkrétní směr (kladná/záporná část) závisí na výsledku výpočtu bližšího souseda. Například, jsme-li blíže k sousedovi vpravo než vlevo, pak proměnná neighbor X představuje souseda v kladném směru osy x, kdežto v opačném případě by znamenala souseda v záporném směru osy x.

B – Uživatelská příručka Instalace Simulační aplikaci je možné velmi snadno nainstalovat pomocí instalačního průvodce programu install.exe, umístěného na přiloženém datovém nosiči. V případě, že cílový počítač není vybaven vývojovým prostředím Microsoft Visual Studio 2008, je nutné po instalaci simulátoru doinstalovat některé

knihovny

tohoto

prostředí

pomocí

!support/vcredist_x86.exe (vytvoří se při instalaci).

programu

100

PŘÍLOHA

Spuštění Simulace se spouští pomocí programu sand-storm.exe. Nastavení scény (tj. světelných zdrojů, objektů a pohledové kamery) je možné provést editací přiloženého souboru settings/config.xml, případně vytvořením nového

XML souboru, předávaného

následně aplikaci

jako

parametr

příkazového řádku. Jeho struktura je následující: <MODEL name path positionX positionY positionZ scale >
= "-150.0" = "300.0" = "250.0" = "0.7" = "0.45" = "0.15" = "1.0" = "1.0" = "1.0" = "0.0" = "0.0" = "0.0"

= = = = = =

"elefant.obj" "models/elefant/" "40.0" "0.0" "20.0" "0.04"

= = = = = =

"200.0" "180.0" "0.0" "0.0" "0.0" "0.0"

PŘÍLOHA

101 upX upY upZ

= "0.0" = "1.0" = "0.0"

>

PROGRAM WIDTH

ŠÍŘKA GRAFICKÉHO OKNA

HEIGHT

VÝŠKA GRAFICKÉHO OKNA

MAXWIDTHHF

ŠÍŘKA VÝŠKOVÉHO POLE (ROZLIŠENÍ V OSE X)

MAXDEPTHHF

HLOUBKA VÝŠKOVÉHO POLE (ROZLIŠENÍ V OSE Z) MODEL

POSITIONX

POZICE MODELU V OSE X

POSITIONY

POZICE MODELU V OSE Y

POSITIONZ

POZICE MODELU V OSE Z

SCALE

MĚŘÍTKO MODELU

NAME

JMÉNO SOUBORU S MODELEM

PATH

CESTA K MODELU LIGHT

POSITIONX

POZICE SVĚTLA V OSE X

POSITIONY

POZICE SVĚTLA V OSE Y

POSITIONZ

POZICE SVĚTLA V OSE Z

AMBIENTR

ČERVENÁ SLOŽKA AMBIENTNÍHO OSVĚTLENÍ

AMBIENTG

ZELENÁ SLOŽKA AMBIENTNÍHO OSVĚTLENÍ

AMBIENTB

MODRÁ SLOŽKA AMBIENTNÍHO OSVĚTLENÍ

DIFFUSER

ČERVENÁ SLOŽKA DIFUZNÍHO OSVĚTLENÍ

DIFFUSEG

ZELENÁ SLOŽKA DIFUZNÍHO OSVĚTLENÍ

DIFFUSEB

MODRÁ SLOŽKA DIFUZNÍHO OSVĚTLENÍ

SPECULARR

ČERVENÁ SLOŽKA ZRCADLOVÉHO (SPEKULÁRNÍHO) OSVĚTLENÍ

SPECULARG

ZELENÁ SLOŽKA ZRCADLOVÉHO (SPEKULÁRNÍHO) OSVĚTLENÍ

SPECULARB

MODRÁ SLOŽKA ZRCADLOVÉHO (SPEKULÁRNÍHO) OSVĚTLENÍ CAMERA

POSITIONX

X SLOŽKA POZICE KAMERY

POSITIONY

Y SLOŽKA POZICE KAMERY

102

PŘÍLOHA POSITIONZ

Z SLOŽKA POZICE KAMERY

VIEWX

X SLOŽKA BODU POHLEDU

VIEWY

Y SLOŽKA BODU POHLEDU

VIEWZ

Z SLOŽKA BODU POHLEDU

UPX

X SLOŽKA UP-VECTORU

UPY

Y SLOŽKA UP-VECTORU

UPZ

Z SLOŽKA UP-VECTORU

Ovládání Samotná fyzikální simulace probíhá bez jakékoliv potřeby interakce s uživatelem. Ovládací prvky tak souvisejí pouze s grafickým modulem. KLÁVESA

AKCE

A

OTOČENÍ DOLEVA

D

OTOČENÍ DOPRAVA

W

POHYB VPŘED

S

POHYB VZAD

Q

PŘEDKLON

E

ZÁKLON

L

PŘEPNUTÍ VYKRESLOVACÍHO MÓDU (DRÁTĚNÝ MODEL / VÝPLŇ)

F1

ZAPNOUT/VYPNOUT VYKRESLOVÁNÍ MRAČNA

F2

ZAPNOUT/VYPNOUT VYKRESLOVÁNÍ MODELŮ

F3

ZAPNOUT/VYPNOUT VYKRESLOVÁNÍ PODLOŽÍ

F4

ZAPNOUT/VYPNOUT VYKRESLENÍ SVĚTELNÝCH ZDROJŮ

F5

ZAPNOUT/VYPNOUT VYKRESLOVÁNÍ VÝPOČETNÍ MŘÍŽKY

F6

ZAPNOUT/VYPNOUT VYKRESLENÍ OBÁLEK MODELŮ (BVH)

F7

ZAPNOUT/VYPNOUT VYKRESLOVÁNÍ RYCHLOSTNÍCH VEKTORŮ

F8

ZAPNOUT/VYPNOUT PROUDOVÉ ČÁRY (STREAMLINES)

F9

ZAPNOUT/VYPNOUT

VYKRESLOVÁNÍ

NORMÁL

MODELOVÝCH

PLOCH

F10

ZAPNOUT/VYPNOUT UKLÁDÁNÍ SNÍMKŮ (VE FORMÁTU BMP)

F11

ZAPNOUT/VYPNOUT OBJECT (VBO)

VYKRESLOVÁNÍ POMOCÍ VERTEX BUFFER

PŘÍLOHA

103

Odinstalace Odebrání simulační aplikace a jejích součástí se provádí spuštěním programu uninstall.exe, který je při instalaci automaticky vytvořen v cílovém adresáři.

C – Obsah přiloženého datového média [DOC]

DOKUMENTACE VE FORMÁTECH HTML A LATEX

[OUTPUT]

GRAFICKÉ VÝSTUPY SIMULÁTORU

[SRC]

ZDROJOVÉ KÓDY A KNIHOVNY

DIP_ERSILDAV.PDF

TENTO TEXT

INSTALL.EXE

INSTALÁTOR SPUSTITELNÉ VERZE SIMULACE

README.TXT

STRUČNÝ POPIS OBSAHU CD

Složka src obsahuje kompletní zdrojové kódy celé aplikace (tj. včetně grafického

modulu).

Předmětem

této

diplomové

práce

byly

soubory

coctree.cpp a solver.cpp a jejich hlavičkové soubory a soubor common.h, nesoucí záznamy o používaných konstantách a makrech.

Simulace písečné bouře

Recommend Documents